Context Navigation

gridcheck.f90 @ fdc0f03

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

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

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

Download in other formats:

Original Format