Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

gridcheck.f90 @ 13

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

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

Download in other formats:

Original Format