Context Navigation

gridcheck.f90 @ 24

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

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

Download in other formats:

Original Format