Context Navigation

gridcheck.f90 @ 8a65cb0

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie

Last change on this file since 8a65cb0 was 8a65cb0, checked in by Espen Sollum ATMOS <espen@…>, 9 years ago
Added code, makefile for dev branch
Property mode set to `100755`
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	! 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	!NIK 16.02.2015
108	tot_blc_count=0 !count for total number of occurences of below cloud scavenging
109	tot_inc_count=0 !count for total number of occurences of in cloud scavenging
110
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
166
167	!print*,'GRiB Edition 2'
168	!read the grib2 identifiers
169	call grib_get_int(igrib,'discipline',discipl,iret)
170	call grib_check(iret,gribFunction,gribErrorMsg)
171	call grib_get_int(igrib,'parameterCategory',parCat,iret)
172	call grib_check(iret,gribFunction,gribErrorMsg)
173	call grib_get_int(igrib,'parameterNumber',parNum,iret)
174	call grib_check(iret,gribFunction,gribErrorMsg)
175	call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
176	call grib_check(iret,gribFunction,gribErrorMsg)
177	call grib_get_int(igrib,'level',valSurf,iret)
178	call grib_check(iret,gribFunction,gribErrorMsg)
179	call grib_get_int(igrib,'paramId',parId,iret)
180	call grib_check(iret,gribFunction,gribErrorMsg)
181
182	!print*,discipl,parCat,parNum,typSurf,valSurf
183
184	!convert to grib1 identifiers
185	isec1(6)=-1
186	isec1(7)=-1
187	isec1(8)=-1
188	isec1(8)=valSurf ! level
189	if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
190	isec1(6)=130 ! indicatorOfParameter
191	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
192	isec1(6)=131 ! indicatorOfParameter
193	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
194	isec1(6)=132 ! indicatorOfParameter
195	elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
196	isec1(6)=133 ! indicatorOfParameter
197	!hg
198	elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
199	isec1(6)=246 ! indicatorOfParameter
200	elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
201	isec1(6)=247 ! indicatorOfParameter
202	!hg end
203	elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
204	isec1(6)=134 ! indicatorOfParameter
205	elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
206	isec1(6)=135 ! indicatorOfParameter
207	elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
208	isec1(6)=135 ! indicatorOfParameter
209	elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
210	isec1(6)=151 ! indicatorOfParameter
211	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
212	isec1(6)=165 ! indicatorOfParameter
213	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
214	isec1(6)=166 ! indicatorOfParameter
215	elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
216	isec1(6)=167 ! indicatorOfParameter
217	elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
218	isec1(6)=168 ! indicatorOfParameter
219	elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
220	isec1(6)=141 ! indicatorOfParameter
221	elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
222	isec1(6)=164 ! indicatorOfParameter
223	elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
224	isec1(6)=142 ! indicatorOfParameter
225	elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
226	isec1(6)=143 ! indicatorOfParameter
227	elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
228	isec1(6)=146 ! indicatorOfParameter
229	elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
230	isec1(6)=176 ! indicatorOfParameter
231	elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
232	isec1(6)=180 ! indicatorOfParameter
233	elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
234	isec1(6)=181 ! indicatorOfParameter
235	elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
236	isec1(6)=129 ! indicatorOfParameter
237	elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
238	isec1(6)=160 ! indicatorOfParameter
239	elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
240	(typSurf.eq.1)) then ! LSM
241	isec1(6)=172 ! indicatorOfParameter
242	else
243	print,'**ERROR: undefined GRiB2 message found!',discipl, &
244	parCat,parNum,typSurf
245	endif
246	if(parId .ne. isec1(6) .and. parId .ne. 77) then
247	write(,) 'parId',parId, 'isec1(6)',isec1(6)
248	! stop
249	endif
250
251	endif
252
253	!get the size and data of the values array
254	if (isec1(6).ne.-1) then
255	call grib_get_real4_array(igrib,'values',zsec4,iret)
256	call grib_check(iret,gribFunction,gribErrorMsg)
257	endif
258
259	if (ifield.eq.1) then
260
261	!HSO get the required fields from section 2 in a gribex compatible manner
262	call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
263	isec2(2),iret)
264	call grib_check(iret,gribFunction,gribErrorMsg)
265	call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
266	isec2(3),iret)
267	call grib_check(iret,gribFunction,gribErrorMsg)
268	call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
269	xaux1in,iret)
270	call grib_check(iret,gribFunction,gribErrorMsg)
271	call grib_get_int(igrib,'numberOfVerticalCoordinateValues', &
272	isec2(12),iret)
273	call grib_check(iret,gribFunction,gribErrorMsg)
274
275	! get the size and data of the vertical coordinate array
276	call grib_get_real4_array(igrib,'pv',zsec2,iret)
277	call grib_check(iret,gribFunction,gribErrorMsg)
278
279	nxfield=isec2(2)
280	ny=isec2(3)
281	nlev_ec=isec2(12)/2-1
282	endif
283
284	!HSO get the second part of the grid dimensions only from GRiB1 messages
285	if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
286	call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
287	xaux2in,iret)
288	call grib_check(iret,gribFunction,gribErrorMsg)
289	call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
290	yaux1in,iret)
291	call grib_check(iret,gribFunction,gribErrorMsg)
292	call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees', &
293	yaux2in,iret)
294	call grib_check(iret,gribFunction,gribErrorMsg)
295	xaux1=xaux1in
296	xaux2=xaux2in
297	yaux1=yaux1in
298	yaux2=yaux2in
299	if (xaux1.gt.180.) xaux1=xaux1-360.0
300	if (xaux2.gt.180.) xaux2=xaux2-360.0
301	if (xaux1.lt.-180.) xaux1=xaux1+360.0
302	if (xaux2.lt.-180.) xaux2=xaux2+360.0
303	if (xaux2.lt.xaux1) xaux2=xaux2+360.0
304	xlon0=xaux1
305	ylat0=yaux1
306	dx=(xaux2-xaux1)/real(nxfield-1)
307	dy=(yaux2-yaux1)/real(ny-1)
308	dxconst=180./(dxr_earthpi)
309	dyconst=180./(dyr_earthpi)
310	gotGrid=1
311	! Check whether fields are global
312	! If they contain the poles, specify polar stereographic map
313	! projections using the stlmbr- and stcm2p-calls
314	!***********************************************************
315
316	xauxa=abs(xaux2+dx-360.-xaux1)
317	if (xauxa.lt.0.001) then
318	nx=nxfield+1 ! field is cyclic
319	xglobal=.true.
320	if (abs(nxshift).ge.nx) &
321	stop 'nxshift in file par_mod is too large'
322	xlon0=xlon0+real(nxshift)*dx
323	else
324	nx=nxfield
325	xglobal=.false.
326	if (nxshift.ne.0) &
327	stop 'nxshift (par_mod) must be zero for non-global domain'
328	endif
329	nxmin1=nx-1
330	nymin1=ny-1
331	if (xlon0.gt.180.) xlon0=xlon0-360.
332	xauxa=abs(yaux1+90.)
333	if (xglobal.and.xauxa.lt.0.001) then
334	sglobal=.true. ! field contains south pole
335	! Enhance the map scale by factor 3 (*2=6) compared to north-south
336	! map scale
337	sizesouth=6.*(switchsouth+90.)/dy
338	call stlmbr(southpolemap,-90.,0.)
339	call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
340	sizesouth,switchsouth,180.)
341	switchsouthg=(switchsouth-ylat0)/dy
342	else
343	sglobal=.false.
344	switchsouthg=999999.
345	endif
346	xauxa=abs(yaux2-90.)
347	if (xglobal.and.xauxa.lt.0.001) then
348	nglobal=.true. ! field contains north pole
349	! Enhance the map scale by factor 3 (*2=6) compared to north-south
350	! map scale
351	sizenorth=6.*(90.-switchnorth)/dy
352	call stlmbr(northpolemap,90.,0.)
353	call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
354	sizenorth,switchnorth,180.)
355	switchnorthg=(switchnorth-ylat0)/dy
356	else
357	nglobal=.false.
358	switchnorthg=999999.
359	endif
360	if (nxshift.lt.0) &
361	stop 'nxshift (par_mod) must not be negative'
362	if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large'
363	endif ! gotGrid
364
365	k=isec1(8)
366	if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
367	if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
368
369	if(isec1(6).eq.129) then
370	do jy=0,ny-1
371	do ix=0,nxfield-1
372	oro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)/ga
373	end do
374	end do
375	endif
376	if(isec1(6).eq.172) then
377	do jy=0,ny-1
378	do ix=0,nxfield-1
379	lsm(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
380	end do
381	end do
382	endif
383	if(isec1(6).eq.160) then
384	do jy=0,ny-1
385	do ix=0,nxfield-1
386	excessoro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
387	end do
388	end do
389	endif
390
391	call grib_release(igrib)
392	goto 10 !! READ NEXT LEVEL OR PARAMETER
393	!
394	! CLOSING OF INPUT DATA FILE
395	!
396
397	30 call grib_close_file(ifile)
398
399	!error message if no fields found with correct first longitude in it
400	if (gotGrid.eq.0) then
401	print,'**ERROR: input file needs to contain GRiB1 formatted'// &
402	'messages'
403	stop
404	endif
405
406	nuvz=iumax
407	nwz =iwmax
408	if(nuvz.eq.nlev_ec) nwz=nlev_ec+1
409
410	if (nx.gt.nxmax) then
411	write(,) 'FLEXPART error: Too many grid points in x direction.'
412	write(,) 'Reduce resolution of wind fields.'
413	write(,) 'Or change parameter settings in file par_mod.'
414	write(,) nx,nxmax
415	stop
416	endif
417
418	if (ny.gt.nymax) then
419	write(,) 'FLEXPART error: Too many grid points in y direction.'
420	write(,) 'Reduce resolution of wind fields.'
421	write(,) 'Or change parameter settings in file par_mod.'
422	write(,) ny,nymax
423	stop
424	endif
425
426	if (nuvz+1.gt.nuvzmax) then
427	write(,) 'FLEXPART error: Too many u,v grid points in z '// &
428	'direction.'
429	write(,) 'Reduce resolution of wind fields.'
430	write(,) 'Or change parameter settings in file par_mod.'
431	write(,) nuvz+1,nuvzmax
432	stop
433	endif
434
435	if (nwz.gt.nwzmax) then
436	write(,) 'FLEXPART error: Too many w grid points in z '// &
437	'direction.'
438	write(,) 'Reduce resolution of wind fields.'
439	write(,) 'Or change parameter settings in file par_mod.'
440	write(,) nwz,nwzmax
441	stop
442	endif
443
444	! If desired, shift all grids by nxshift grid cells
445	!**************************************************
446
447	if (xglobal) then
448	call shift_field_0(oro,nxfield,ny)
449	call shift_field_0(lsm,nxfield,ny)
450	call shift_field_0(excessoro,nxfield,ny)
451	endif
452
453	! Output of grid info
454	!********************
455
456	if (lroot) then
457	write(*,'(a,2i7)') ' Vertical levels in ECMWF data: ', &
458	nuvz+1,nwz
459	write(,)
460	write(*,'(a)') ' Mother domain:'
461	write(*,'(a,f10.5,a,f10.5,a,f10.5)') ' Longitude range: ', &
462	xlon0,' to ',xlon0+(nx-1)*dx,' Grid distance: ',dx
463	write(*,'(a,f10.5,a,f10.5,a,f10.5)') ' Latitude range : ', &
464	ylat0,' to ',ylat0+(ny-1)*dy,' Grid distance: ',dy
465	write(,)
466	end if
467
468	! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL
469	! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
470
471	numskip=nlev_ec-nuvz ! number of ecmwf model layers not used
472	! by trajectory model
473	!do 8940 i=1,244
474	! write (,) 'zsec2:',i,ifield,zsec2(i),numskip
475	!940 continue
476	! stop
477	! SEC SEC SEC
478	! for unknown reason zsec 1 to 10 is filled in this version
479	! compared to the old one
480	! SEC SEC SE
481	do i=1,nwz
482	j=numskip+i
483	k=nlev_ec+1+numskip+i
484	akm(nwz-i+1)=zsec2(j)
485	! write (,) 'ifield:',ifield,k,j,zsec2(10+j)
486	bkm(nwz-i+1)=zsec2(k)
487	end do
488
489	!
490	! CALCULATION OF AKZ, BKZ
491	! AKZ,BKZ: model discretization parameters at the center of each model
492	! layer
493	!
494	! Assign the 10 m winds to an artificial model level with akz=0 and bkz=1.0,
495	! i.e. ground level
496	!*****************************************************************************
497
498	akz(1)=0.
499	bkz(1)=1.0
500	do i=1,nuvz
501	akz(i+1)=0.5*(akm(i+1)+akm(i))
502	bkz(i+1)=0.5*(bkm(i+1)+bkm(i))
503	end do
504	nuvz=nuvz+1
505
506	! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
507	! upon the transformation to z levels. In order to save computer memory, this is
508	! not done anymore in the standard version. However, this option can still be
509	! switched on by replacing the following lines with those below, that are
510	! currently commented out. For this, similar changes are necessary in
511	! verttransform.f and verttranform_nests.f
512	!*****************************************************************************
513
514	nz=nuvz
515	if (nz.gt.nzmax) stop 'nzmax too small'
516	do i=1,nuvz
517	aknew(i)=akz(i)
518	bknew(i)=bkz(i)
519	end do
520
521	! Switch on following lines to use doubled vertical resolution
522	!*************************************************************
523	!nz=nuvz+nwz-1
524	!if (nz.gt.nzmax) stop 'nzmax too small'
525	!do 100 i=1,nwz
526	! aknew(2*(i-1)+1)=akm(i)
527	!00 bknew(2*(i-1)+1)=bkm(i)
528	!do 110 i=2,nuvz
529	! aknew(2*(i-1))=akz(i)
530	!10 bknew(2*(i-1))=bkz(i)
531	! End doubled vertical resolution
532
533
534	! Determine the uppermost level for which the convection scheme shall be applied
535	! by assuming that there is no convection above 50 hPa (for standard SLP)
536	!*****************************************************************************
537
538	do i=1,nuvz-2
539	pint=akz(i)+bkz(i)*101325.
540	if (pint.lt.5000.) goto 96
541	end do
542	96 nconvlev=i
543	if (nconvlev.gt.nconvlevmax-1) then
544	nconvlev=nconvlevmax-1
545	write(,) 'Attention, convection only calculated up to ', &
546	akz(nconvlev)+bkz(nconvlev)*1013.25,' hPa'
547	endif
548
549	return
550
551	999 write(,)
552	write(,) ' ###########################################'// &
553	'###### '
554	write(,) ' TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
555	write(,) ' CAN NOT OPEN INPUT DATA FILE '//wfname(ifn)
556	write(,) ' ###########################################'// &
557	'###### '
558	write(,)
559	write(*,'(a)') '!!! PLEASE INSERT A NEW CD-ROM AND !!!'
560	write(*,'(a)') '!!! PRESS ANY KEY TO CONTINUE... !!!'
561	write(*,'(a)') '!!! ...OR TERMINATE FLEXPART PRESSING!!!'
562	write(*,'(a)') '!!! THE "X" KEY... !!!'
563	write(,)
564	read(*,'(a)') opt
565	if(opt.eq.'X') then
566	stop
567	else
568	goto 5
569	endif
570
571	end subroutine gridcheck
572

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

Download in other formats:

Original Format