Context Navigation

gridcheck.f90 @ 37

Last change on this file since 37 was 37, checked in by pesei, 9 years ago
Wet dep quick fix and other small changes. Wet depo quick fix not final yet.
Property svn:executable set to ``*
File size: 21.1 KB

Line
1	!**********************************************************************
2	! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
3	! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
4	! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
5	! *
6	! This file is part of FLEXPART. *
7	! *
8	! FLEXPART is free software: you can redistribute it and/or modify *
9	! it under the terms of the GNU General Public License as published by*
10	! the Free Software Foundation, either version 3 of the License, or *
11	! (at your option) any later version. *
12	! *
13	! FLEXPART is distributed in the hope that it will be useful, *
14	! but WITHOUT ANY WARRANTY; without even the implied warranty of *
15	! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
16	! GNU General Public License for more details. *
17	! *
18	! You should have received a copy of the GNU General Public License *
19	! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
20	!**********************************************************************
21
22	subroutine gridcheck
23
24	!**********************************************************************
25	! *
26	! FLEXPART MODEL SUBROUTINE GRIDCHECK *
27	! *
28	!**********************************************************************
29	! *
30	! AUTHOR: G. WOTAWA *
31	! DATE: 1997-08-06 *
32	! LAST UPDATE: 1997-10-10 *
33	! *
34	!**********************************************************************
35	! *
36	! Update: 1999-02-08, global fields allowed, A. Stohl*
37	! CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
38	! ECMWF grib_api *
39	! CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
40	! ECMWF grib_api *
41	!
42	! 3/2015, PS: add GRIB2 coding for deta/dt dp/deta used by ZAMG
43	!
44	!
45	!**********************************************************************
46	! *
47	! DESCRIPTION: *
48	! *
49	! THIS SUBROUTINE DETERMINES THE GRID SPECIFICATIONS (LOWER LEFT *
50	! LONGITUDE, LOWER LEFT LATITUDE, NUMBER OF GRID POINTS, GRID DIST- *
51	! ANCE AND VERTICAL DISCRETIZATION OF THE ECMWF MODEL) FROM THE *
52	! GRIB HEADER OF THE FIRST INPUT FILE. THE CONSISTANCY (NO CHANGES *
53	! WITHIN ONE FLEXPART RUN) IS CHECKED IN THE ROUTINE "READWIND" AT *
54	! ANY CALL. *
55	! *
56	! XLON0 geographical longitude of lower left gridpoint *
57	! YLAT0 geographical latitude of lower left gridpoint *
58	! NX number of grid points x-direction *
59	! NY number of grid points y-direction *
60	! DX grid distance x-direction *
61	! DY grid distance y-direction *
62	! NUVZ number of grid points for horizontal wind *
63	! components in z direction *
64	! NWZ number of grid points for vertical wind *
65	! component in z direction *
66	! sizesouth, sizenorth give the map scale (i.e. number of virtual grid*
67	! points of the polar stereographic grid): *
68	! used to check the CFL criterion *
69	! UVHEIGHT(1)- heights of gridpoints where u and v are *
70	! UVHEIGHT(NUVZ) given *
71	! WHEIGHT(1)- heights of gridpoints where w is given *
72	! WHEIGHT(NWZ) *
73	! *
74	!**********************************************************************
75
76	use grib_api
77	use par_mod
78	use com_mod
79	use conv_mod
80	use cmapf_mod, only: stlmbr,stcm2p
81
82	implicit none
83
84	!HSO parameters for grib_api
85	integer :: ifile
86	integer :: iret
87	integer :: igrib
88	integer :: gotGrid
89	real(kind=4) :: xaux1,xaux2,yaux1,yaux2
90	real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
91	integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parID
92	!HSO end
93	integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
94	real :: sizesouth,sizenorth,xauxa,pint,conversion_factor
95
96	! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
97
98	! dimension of isec2 at least (22+n), where n is the number of parallels or
99	! meridians in a quasi-regular (reduced) Gaussian or lat/long grid
100
101	! dimension of zsec2 at least (10+nn), where nn is the number of vertical
102	! coordinate parameters
103
104	integer :: isec1(56),isec2(22+nxmax+nymax)
105	real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
106	character(len=1) :: opt
107
108	!HSO grib api error messages
109	character(len=24) :: gribErrorMsg = 'Error reading grib file'
110	character(len=20) :: gribFunction = 'gridcheck'
111
112	iumax=0
113	iwmax=0
114
115	if(ideltas.gt.0) then
116	ifn=1
117	else
118	ifn=numbwf
119	endif
120	!
121	! OPENING OF DATA FILE (GRIB CODE)
122	!
123	5 call grib_open_file(ifile,path(3)(1:length(3)) &
124	//trim(wfname(ifn)),'r',iret)
125	if (iret.ne.GRIB_SUCCESS) then
126	goto 999 ! ERROR DETECTED
127	endif
128	!turn on support for multi fields messages
129	!call grib_multi_support_on
130
131	gotGrid=0
132	ifield=0
133	10 ifield=ifield+1
134
135	!
136	! GET NEXT FIELDS
137	!
138	call grib_new_from_file(ifile,igrib,iret)
139	if (iret.eq.GRIB_END_OF_FILE ) then
140	goto 30 ! EOF DETECTED
141	elseif (iret.ne.GRIB_SUCCESS) then
142	goto 999 ! ERROR DETECTED
143	endif
144
145	!first see if we read GRIB1 or GRIB2
146	call grib_get_int(igrib,'editionNumber',gribVer,iret)
147	call grib_check(iret,gribFunction,gribErrorMsg)
148
149	if (gribVer.eq.1) then ! GRIB Edition 1
150
151	!print*,'GRiB Edition 1'
152	!read the grib2 identifiers
153	call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
154	call grib_check(iret,gribFunction,gribErrorMsg)
155	call grib_get_int(igrib,'level',isec1(8),iret)
156	call grib_check(iret,gribFunction,gribErrorMsg)
157
158	!change code for etadot to code for omega
159	if (isec1(6).eq.77) then
160	isec1(6)=135
161	endif
162
163	!print*,isec1(6),isec1(8)
164
165	else ! GRIB 2
166
167	!print*,'GRiB Edition 2'
168	! read the grib2 identifiers
169
170	call grib_get_int(igrib,'discipline',discipl,iret)
171	call grib_check(iret,gribFunction,gribErrorMsg)
172	call grib_get_int(igrib,'parameterCategory',parCat,iret)
173	call grib_check(iret,gribFunction,gribErrorMsg)
174	call grib_get_int(igrib,'parameterNumber',parNum,iret)
175	call grib_check(iret,gribFunction,gribErrorMsg)
176	call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
177	call grib_check(iret,gribFunction,gribErrorMsg)
178	call grib_get_int(igrib,'level',valSurf,iret)
179	call grib_check(iret,gribFunction,gribErrorMsg)
180	call grib_get_int(igrib,'paramId',parId,iret)
181	call grib_check(iret,gribFunction,gribErrorMsg)
182
183	!print*,discipl,parCat,parNum,typSurf,valSurf
184
185	! convert to grib1 identifiers
186
187	isec1(6)=-1
188	isec1(7)=-1
189	isec1(8)=-1
190	isec1(8)=valSurf ! level
191	if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
192	isec1(6)=130 ! indicatorOfParameter
193	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
194	isec1(6)=131 ! indicatorOfParameter
195	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
196	isec1(6)=132 ! indicatorOfParameter
197	elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
198	isec1(6)=133 ! indicatorOfParameter
199	elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
200	isec1(6)=134 ! indicatorOfParameter
201	elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually detadtdpdeta
202	isec1(6)=135 ! indicatorOfParameter
203	elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually detadtdpdeta
204	isec1(6)=135 ! indicatorOfParameter
205	elseif ((parCat.eq.2).and.(parNum.eq.8)) then ! Omega, actually detadtdpdeta
206	isec1(6)=135 ! indicatorOfParameter
207	elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
208	isec1(6)=151 ! indicatorOfParameter
209	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
210	isec1(6)=165 ! indicatorOfParameter
211	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
212	isec1(6)=166 ! indicatorOfParameter
213	elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
214	isec1(6)=167 ! indicatorOfParameter
215	elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
216	isec1(6)=168 ! indicatorOfParameter
217	elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
218	isec1(6)=141 ! indicatorOfParameter
219	elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
220	isec1(6)=164 ! indicatorOfParameter
221	elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
222	isec1(6)=142 ! indicatorOfParameter
223	elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
224	isec1(6)=143 ! indicatorOfParameter
225	elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
226	isec1(6)=146 ! indicatorOfParameter
227	elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
228	isec1(6)=176 ! indicatorOfParameter
229	elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
230	isec1(6)=180 ! indicatorOfParameter
231	elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
232	isec1(6)=181 ! indicatorOfParameter
233	elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
234	isec1(6)=129 ! indicatorOfParameter
235	elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
236	isec1(6)=160 ! indicatorOfParameter
237	elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
238	(typSurf.eq.1)) then ! LSM
239	isec1(6)=172 ! indicatorOfParameter
240	else
241	print,'**ERROR: undefined GRiB2 message found!',discipl, &
242	parCat,parNum,typSurf
243	endif
244	if(parId .ne. isec1(6) .and. parId .ne. 77) then
245	write(,) 'parId',parId, 'isec1(6)',isec1(6)
246	! stop
247	endif
248
249	endif ! GRIB 1 / GRIB 2 cases
250
251	!get the size and data of the values array
252	if (isec1(6).ne.-1) then
253	call grib_get_real4_array(igrib,'values',zsec4,iret)
254	call grib_check(iret,gribFunction,gribErrorMsg)
255	endif
256
257	if (ifield.eq.1) then
258
259	!HSO get the required fields from section 2 in a gribex compatible manner
260	call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
261	isec2(2),iret)
262	call grib_check(iret,gribFunction,gribErrorMsg)
263	call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
264	isec2(3),iret)
265	call grib_check(iret,gribFunction,gribErrorMsg)
266	call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
267	xaux1in,iret)
268	call grib_check(iret,gribFunction,gribErrorMsg)
269	call grib_get_int(igrib,'numberOfVerticalCoordinateValues', &
270	isec2(12),iret)
271	call grib_check(iret,gribFunction,gribErrorMsg)
272
273	! get the size and data of the vertical coordinate array
274	call grib_get_real4_array(igrib,'pv',zsec2,iret)
275	call grib_check(iret,gribFunction,gribErrorMsg)
276
277	nxfield=isec2(2)
278	ny=isec2(3)
279	nlev_ec=isec2(12)/2-1
280	endif
281
282	!HSO get the second part of the grid dimensions only from GRiB1 messages
283	if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
284	call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
285	xaux2in,iret)
286	call grib_check(iret,gribFunction,gribErrorMsg)
287	call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
288	yaux1in,iret)
289	call grib_check(iret,gribFunction,gribErrorMsg)
290	call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees', &
291	yaux2in,iret)
292	call grib_check(iret,gribFunction,gribErrorMsg)
293	xaux1=xaux1in
294	xaux2=xaux2in
295	yaux1=yaux1in
296	yaux2=yaux2in
297	if (xaux1.gt.180.) xaux1=xaux1-360.0
298	if (xaux2.gt.180.) xaux2=xaux2-360.0
299	if (xaux1.lt.-180.) xaux1=xaux1+360.0
300	if (xaux2.lt.-180.) xaux2=xaux2+360.0
301	if (xaux2.lt.xaux1) xaux2=xaux2+360.0
302	xlon0=xaux1
303	ylat0=yaux1
304	dx=(xaux2-xaux1)/real(nxfield-1)
305	dy=(yaux2-yaux1)/real(ny-1)
306	dxconst=180./(dxr_earthpi)
307	dyconst=180./(dyr_earthpi)
308	gotGrid=1
309	! Check whether fields are global
310	! If they contain the poles, specify polar stereographic map
311	! projections using the stlmbr- and stcm2p-calls
312	!***********************************************************
313
314	xauxa=abs(xaux2+dx-360.-xaux1)
315	if (xauxa.lt.0.001) then
316	nx=nxfield+1 ! field is cyclic
317	xglobal=.true.
318	if (abs(nxshift).ge.nx) &
319	stop 'nxshift in file par_mod is too large'
320	xlon0=xlon0+real(nxshift)*dx
321	else
322	nx=nxfield
323	xglobal=.false.
324	if (nxshift.ne.0) &
325	stop 'nxshift (par_mod) must be zero for non-global domain'
326	endif
327	nxmin1=nx-1
328	nymin1=ny-1
329	if (xlon0.gt.180.) xlon0=xlon0-360.
330	xauxa=abs(yaux1+90.)
331	if (xglobal.and.xauxa.lt.0.001) then
332	sglobal=.true. ! field contains south pole
333	! Enhance the map scale by factor 3 (*2=6) compared to north-south
334	! map scale
335	sizesouth=6.*(switchsouth+90.)/dy
336	call stlmbr(southpolemap,-90.,0.)
337	call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
338	sizesouth,switchsouth,180.)
339	switchsouthg=(switchsouth-ylat0)/dy
340	else
341	sglobal=.false.
342	switchsouthg=999999.
343	endif
344	xauxa=abs(yaux2-90.)
345	if (xglobal.and.xauxa.lt.0.001) then
346	nglobal=.true. ! field contains north pole
347	! Enhance the map scale by factor 3 (*2=6) compared to north-south
348	! map scale
349	sizenorth=6.*(90.-switchnorth)/dy
350	call stlmbr(northpolemap,90.,0.)
351	call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
352	sizenorth,switchnorth,180.)
353	switchnorthg=(switchnorth-ylat0)/dy
354	else
355	nglobal=.false.
356	switchnorthg=999999.
357	endif
358	if (nxshift.lt.0) &
359	stop 'nxshift (par_mod) must not be negative'
360	if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large'
361	endif ! gotGrid
362
363	k=isec1(8)
364	if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
365	if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
366
367	if(isec1(6).eq.129) then
368	do jy=0,ny-1
369	do ix=0,nxfield-1
370	oro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)/ga
371	end do
372	end do
373	endif
374	if(isec1(6).eq.172) then
375	do jy=0,ny-1
376	do ix=0,nxfield-1
377	lsm(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
378	end do
379	end do
380	endif
381	if(isec1(6).eq.160) then
382	do jy=0,ny-1
383	do ix=0,nxfield-1
384	excessoro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
385	end do
386	end do
387	endif
388
389	call grib_release(igrib)
390	goto 10 !! READ NEXT LEVEL OR PARAMETER
391	!
392	! CLOSING OF INPUT DATA FILE
393	!
394
395	30 call grib_close_file(ifile)
396
397	!error message if no fields found with correct first longitude in it
398	if (gotGrid.eq.0) then
399	print,'**ERROR: input file needs to contain GRiB1 formatted'// &
400	'messages'
401	stop
402	endif
403
404	nuvz=iumax
405	nwz =iwmax
406	if(nuvz.eq.nlev_ec) nwz=nlev_ec+1
407
408	if (nx.gt.nxmax) then
409	write(,) 'FLEXPART error: Too many grid points in x direction.'
410	write(,) 'Reduce resolution of wind fields.'
411	write(,) 'Or change parameter settings in file par_mod.'
412	write(,) nx,nxmax
413	stop
414	endif
415
416	if (ny.gt.nymax) then
417	write(,) 'FLEXPART error: Too many grid points in y direction.'
418	write(,) 'Reduce resolution of wind fields.'
419	write(,) 'Or change parameter settings in file par_mod.'
420	write(,) ny,nymax
421	stop
422	endif
423
424	if (nuvz+1.gt.nuvzmax) then
425	write(,) 'FLEXPART error: Too many u,v grid points in z '// &
426	'direction.'
427	write(,) 'Reduce resolution of wind fields.'
428	write(,) 'Or change parameter settings in file par_mod.'
429	write(,) nuvz+1,nuvzmax
430	stop
431	endif
432
433	if (nwz.gt.nwzmax) then
434	write(,) 'FLEXPART error: Too many w grid points in z '// &
435	'direction.'
436	write(,) 'Reduce resolution of wind fields.'
437	write(,) 'Or change parameter settings in file par_mod.'
438	write(,) nwz,nwzmax
439	stop
440	endif
441
442	! If desired, shift all grids by nxshift grid cells
443	!**************************************************
444
445	if (xglobal) then
446	call shift_field_0(oro,nxfield,ny)
447	call shift_field_0(lsm,nxfield,ny)
448	call shift_field_0(excessoro,nxfield,ny)
449	endif
450
451	! Output of grid info
452	!********************
453
454	write(*,'(a,2i7)') ' Vertical levels in ECMWF data: ', &
455	nuvz+1,nwz
456	write(,)
457	write(*,'(a)') ' Mother domain:'
458	write(*,'(a,f10.5,a,f10.5,a,f10.5)') ' Longitude range: ', &
459	xlon0,' to ',xlon0+(nx-1)*dx,' Grid distance: ',dx
460	write(*,'(a,f10.5,a,f10.5,a,f10.5)') ' Latitude range : ', &
461	ylat0,' to ',ylat0+(ny-1)*dy,' Grid distance: ',dy
462	write(,)
463
464
465	! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL
466	! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
467
468	numskip=nlev_ec-nuvz ! number of ecmwf model layers not used
469	! by trajectory model
470	!do 8940 i=1,244
471	! write (,) 'zsec2:',i,ifield,zsec2(i),numskip
472	!940 continue
473	! stop
474	! SEC SEC SEC
475	! for unknown reason zsec 1 to 10 is filled in this version
476	! compared to the old one
477	! SEC SEC SE
478	do i=1,nwz
479	j=numskip+i
480	k=nlev_ec+1+numskip+i
481	akm(nwz-i+1)=zsec2(j)
482	! write (,) 'ifield:',ifield,k,j,zsec2(10+j)
483	bkm(nwz-i+1)=zsec2(k)
484	end do
485
486	!
487	! CALCULATION OF AKZ, BKZ
488	! AKZ,BKZ: model discretization parameters at the center of each model
489	! layer
490	!
491	! Assign the 10 m winds to an artificial model level with akz=0 and bkz=1.0,
492	! i.e. ground level
493	!*****************************************************************************
494
495	akz(1)=0.
496	bkz(1)=1.0
497	do i=1,nuvz
498	akz(i+1)=0.5*(akm(i+1)+akm(i))
499	bkz(i+1)=0.5*(bkm(i+1)+bkm(i))
500	end do
501	nuvz=nuvz+1
502
503	! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
504	! upon the transformation to z levels. In order to save computer memory, this is
505	! not done anymore in the standard version. However, this option can still be
506	! switched on by replacing the following lines with those below, that are
507	! currently commented out. For this, similar changes are necessary in
508	! verttransform.f and verttranform_nests.f
509	!*****************************************************************************
510
511	nz=nuvz
512	if (nz.gt.nzmax) stop 'nzmax too small'
513	do i=1,nuvz
514	aknew(i)=akz(i)
515	bknew(i)=bkz(i)
516	end do
517
518	! Switch on following lines to use doubled vertical resolution
519	!*************************************************************
520	!nz=nuvz+nwz-1
521	!if (nz.gt.nzmax) stop 'nzmax too small'
522	!do 100 i=1,nwz
523	! aknew(2*(i-1)+1)=akm(i)
524	!00 bknew(2*(i-1)+1)=bkm(i)
525	!do 110 i=2,nuvz
526	! aknew(2*(i-1))=akz(i)
527	!10 bknew(2*(i-1))=bkz(i)
528	! End doubled vertical resolution
529
530
531	! Determine the uppermost level for which the convection scheme shall be applied
532	! by assuming that there is no convection above 50 hPa (for standard SLP)
533	!*****************************************************************************
534
535	do i=1,nuvz-2
536	pint=akz(i)+bkz(i)*101325.
537	if (pint.lt.5000.) goto 96
538	end do
539	96 nconvlev=i
540	if (nconvlev.gt.nconvlevmax-1) then
541	nconvlev=nconvlevmax-1
542	write(,) 'Attention, convection only calculated up to ', &
543	akz(nconvlev)+bkz(nconvlev)*1013.25,' hPa'
544	endif
545
546	return
547
548	999 write(,)
549	write(,) ' ###########################################'// &
550	'###### '
551	write(,) ' TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
552	write(,) ' CAN NOT OPEN INPUT DATA FILE '//wfname(ifn)
553	write(,) ' ###########################################'// &
554	'###### '
555	write(,)
556	write(*,'(a)') '!!! PLEASE INSERT A NEW CD-ROM AND !!!'
557	write(*,'(a)') '!!! PRESS ANY KEY TO CONTINUE... !!!'
558	write(*,'(a)') '!!! ...OR TERMINATE FLEXPART PRESSING!!!'
559	write(*,'(a)') '!!! THE "X" KEY... !!!'
560	write(,)
561	read(*,'(a)') opt
562	if(opt.eq.'X') then
563	stop
564	else
565	goto 5
566	endif
567
568	end subroutine gridcheck
569

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format