Context Navigation

gridcheck_gfs.F90 @ b398fb6

FPv9.3.1FPv9.3.1b_testingFPv9.3.2fp9.3.1-20161214-nc4grib2nc4_repair

Last change on this file since b398fb6 was 496c607, checked in by Don Morton <Don.Morton@…>, 8 years ago

Initial commit of FPv9.3.1

Currently, this is a clone of snapshot FPv9.3.0

Property mode set to 100644

File size: 19.8 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_gfs
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: 17/11/2005, Caroline Forster, GFS data *
36	! CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
37	! ECMWF grib_api *
38	! CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
39	! ECMWF grib_api *
40	! *
41	!**********************************************************************
42	! *
43	! DESCRIPTION: *
44	! *
45	! THIS SUBROUTINE DETERMINES THE GRID SPECIFICATIONS (LOWER LEFT *
46	! LONGITUDE, LOWER LEFT LATITUDE, NUMBER OF GRID POINTS, GRID DIST- *
47	! ANCE AND VERTICAL DISCRETIZATION OF THE ECMWF MODEL) FROM THE *
48	! GRIB HEADER OF THE FIRST INPUT FILE. THE CONSISTANCY (NO CHANGES *
49	! WITHIN ONE FLEXPART RUN) IS CHECKED IN THE ROUTINE "READWIND" AT *
50	! ANY CALL. *
51	! *
52	! XLON0 geographical longitude of lower left gridpoint *
53	! YLAT0 geographical latitude of lower left gridpoint *
54	! NX number of grid points x-direction *
55	! NY number of grid points y-direction *
56	! DX grid distance x-direction *
57	! DY grid distance y-direction *
58	! NUVZ number of grid points for horizontal wind *
59	! components in z direction *
60	! NWZ number of grid points for vertical wind *
61	! component in z direction *
62	! sizesouth, sizenorth give the map scale (i.e. number of virtual grid*
63	! points of the polar stereographic grid): *
64	! used to check the CFL criterion *
65	! UVHEIGHT(1)- heights of gridpoints where u and v are *
66	! UVHEIGHT(NUVZ) given *
67	! WHEIGHT(1)- heights of gridpoints where w is given *
68	! WHEIGHT(NWZ) *
69	! *
70	!**********************************************************************
71
72	use grib_api
73	use par_mod
74	use com_mod
75	use conv_mod
76	use cmapf_mod, only: stlmbr,stcm2p
77	use class_vtable
78
79	implicit none
80
81	!HSO parameters for grib_api
82	integer :: ifile
83	integer :: iret
84	integer :: igrib
85	real(kind=4) :: xaux1,xaux2,yaux1,yaux2
86	real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
87	integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
88	!HSO end
89	integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
90	real :: sizesouth,sizenorth,xauxa,pint
91	real :: akm_usort(nwzmax)
92	real,parameter :: eps=0.0001
93
94	! NCEP GFS
95	real :: pres(nwzmax), help
96
97	integer :: i179,i180,i181
98
99	! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
100
101	integer :: isec2(3)
102	real(kind=4) :: zsec4(jpunp)
103	character(len=1) :: opt
104
105	!HSO grib api error messages
106	character(len=24) :: gribErrorMsg = 'Error reading grib file'
107	character(len=20) :: gribFunction = 'gridcheckwind_gfs'
108	!
109
110
111
112	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
113	!!!! Vtable related variables
114
115	! Paths to Vtables (current implementation will assume they are in the cwd
116	! This is a crappy place for these parameters. Need to move them.
117	character(LEN=255), parameter :: VTABLE_NCEP_GRIB1_PATH = &
118	"Vtable_ncep_grib1", &
119	VTABLE_NCEP_GRIB2_PATH = &
120	"Vtable_ncep_grib2"
121
122
123
124	integer :: gribfile_type
125	integer :: current_grib_level ! This "was" isec1(8) in previous version
126	character(len=255) :: gribfile_name
127	character(len=255) :: vtable_path
128	character(len=15) :: fpname
129
130	type(Vtable) :: my_vtable ! unallocated
131	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
132
133
134
135
136
137
138
139
140	if (numbnests.ge.1) then
141	write(,) ' ###########################################'
142	write(,) ' FLEXPART ERROR SUBROUTINE GRIDCHECK:'
143	write(,) ' NO NESTED WINDFIELDAS ALLOWED FOR GFS! '
144	write(,) ' ###########################################'
145	stop
146	endif
147
148	iumax=0
149	iwmax=0
150
151	if(ideltas.gt.0) then
152	ifn=1
153	else
154	ifn=numbwf
155	endif
156
157
158
159	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
160	!!!!!!!!!!!!!!!!!!! VTABLE code
161	!!!!!!! Vtable choice
162	gribfile_name = path(3)(1:length(3))//trim(wfname(ifn))
163	print *, 'gribfile_name: ', gribfile_name
164
165	gribfile_type = vtable_detect_gribfile_type( gribfile_name )
166
167	print *, 'gribfile_type: ', gribfile_type
168
169	if (gribfile_type .eq. VTABLE_GRIBFILE_TYPE_NCEP_GRIB1) then
170	vtable_path = VTABLE_NCEP_GRIB1_PATH
171	else if (gribfile_type .eq. VTABLE_GRIBFILE_TYPE_NCEP_GRIB2) then
172	vtable_path = VTABLE_NCEP_GRIB2_PATH
173	else
174	print *, 'Unsupported gribfile_type: ', gribfile_type
175	stop
176	endif
177
178
179	! Load the Vtable into 'my_vtable'
180	print *, 'Loading Vtable: ', vtable_path
181	call vtable_load_by_name(vtable_path, my_vtable)
182	print *, 'Vtable Initialized: ', my_vtable%initialized
183	print *, 'Vtable num_entries: ', my_vtable%num_entries
184	!!!!!!!!!!!!!!!!!!! VTABLE code
185	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
186
187
188
189
190
191
192
193	!
194	! OPENING OF DATA FILE (GRIB CODE)
195	!
196	5 call grib_open_file(ifile,path(3)(1:length(3)) &
197	//trim(wfname(ifn)),'r',iret)
198	if (iret.ne.GRIB_SUCCESS) then
199	goto 999 ! ERROR DETECTED
200	endif
201	!turn on support for multi fields messages
202	call grib_multi_support_on
203
204	ifield=0
205	10 ifield=ifield+1
206	!
207	! GET NEXT FIELDS
208	!
209	call grib_new_from_file(ifile,igrib,iret)
210	if (iret.eq.GRIB_END_OF_FILE ) then
211	goto 30 ! EOF DETECTED
212	elseif (iret.ne.GRIB_SUCCESS) then
213	goto 999 ! ERROR DETECTED
214	endif
215
216
217	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
218	!!!!!!!!!!!!!!!!!!! VTABLE code
219	! Get the fpname
220	fpname = vtable_get_fpname(igrib, my_vtable)
221	!print *, 'fpname: ', trim(fpname)
222
223
224	!!!!!!!!!!!!!!!!!!! VTABLE code
225	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
226
227
228
229	!first see if we read GRIB1 or GRIB2
230	call grib_get_int(igrib,'editionNumber',gribVer,iret)
231	call grib_check(iret,gribFunction,gribErrorMsg)
232
233	if (gribVer.eq.1) then ! GRIB Edition 1
234
235	call grib_get_int(igrib,'level',current_grib_level,iret)
236	call grib_check(iret,gribFunction,gribErrorMsg)
237
238	!!! Added by DJM 2016-03-02 - if this is GRIB1 we assume that
239	!!! level units are hPa and need to be multiplied by 100 for Pa
240	current_grib_level = current_grib_level*100.0
241
242	else ! GRIB Edition 2
243
244	call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
245	current_grib_level,iret)
246	call grib_check(iret,gribFunction,gribErrorMsg)
247
248	!!! Added by DJM 2016-03-02 - if this is GRIB2 we assume that
249	!!! level units are Pa and don't need to be modified
250
251	endif ! gribVer
252
253	if ( trim(fpname) .ne. 'None' ) then
254	! get the size and data of the values array
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	!get the required fields from section 2
262	!store compatible to gribex input
263	call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
264	isec2(2),iret)
265	call grib_check(iret,gribFunction,gribErrorMsg)
266	call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
267	isec2(3),iret)
268	call grib_check(iret,gribFunction,gribErrorMsg)
269	call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
270	xaux1in,iret)
271	call grib_check(iret,gribFunction,gribErrorMsg)
272	call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
273	xaux2in,iret)
274	call grib_check(iret,gribFunction,gribErrorMsg)
275	call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
276	yaux1in,iret)
277	call grib_check(iret,gribFunction,gribErrorMsg)
278	call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees', &
279	yaux2in,iret)
280	call grib_check(iret,gribFunction,gribErrorMsg)
281
282
283	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
284	!!!!!!!! DJM ARTIFICIAL CHANGE FOR NCEP GRIB1 - change value from -1 to 359
285	!!!!!!!! See flexpart.eu CTBTO project Ticket #112
286	if (xaux2in .lt. 0) xaux2in = 359.0
287	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
288
289
290
291	xaux1=xaux1in
292	xaux2=xaux2in
293	yaux1=yaux1in
294	yaux2=yaux2in
295
296	nxfield=isec2(2)
297	ny=isec2(3)
298	if((abs(xaux1).lt.eps).and.(xaux2.ge.359)) then ! NCEP DATA FROM 0 TO
299	xaux1=-179.0 ! 359 DEG EAST ->
300	xaux2=-179.0+360.-360./real(nxfield) ! TRANSFORMED TO -179
301	endif ! TO 180 DEG EAST
302	if (xaux1.gt.180) xaux1=xaux1-360.0
303	if (xaux2.gt.180) xaux2=xaux2-360.0
304	if (xaux1.lt.-180) xaux1=xaux1+360.0
305	if (xaux2.lt.-180) xaux2=xaux2+360.0
306	if (xaux2.lt.xaux1) xaux2=xaux2+360.
307	xlon0=xaux1
308	ylat0=yaux1
309	dx=(xaux2-xaux1)/real(nxfield-1)
310	dy=(yaux2-yaux1)/real(ny-1)
311	dxconst=180./(dxr_earthpi)
312	dyconst=180./(dyr_earthpi)
313	!HSO end edits
314
315
316	! Check whether fields are global
317	! If they contain the poles, specify polar stereographic map
318	! projections using the stlmbr- and stcm2p-calls
319	!***********************************************************
320
321	xauxa=abs(xaux2+dx-360.-xaux1)
322	if (xauxa.lt.0.001) then
323	nx=nxfield+1 ! field is cyclic
324	xglobal=.true.
325	if (abs(nxshift).ge.nx) &
326	stop 'nxshift in file par_mod is too large'
327	xlon0=xlon0+real(nxshift)*dx
328	else
329	nx=nxfield
330	xglobal=.false.
331	if (nxshift.ne.0) &
332	stop 'nxshift (par_mod) must be zero for non-global domain'
333	endif
334	nxmin1=nx-1
335	nymin1=ny-1
336	if (xlon0.gt.180.) xlon0=xlon0-360.
337	xauxa=abs(yaux1+90.)
338	if (xglobal.and.xauxa.lt.0.001) then
339	sglobal=.true. ! field contains south pole
340	! Enhance the map scale by factor 3 (*2=6) compared to north-south
341	! map scale
342	sizesouth=6.*(switchsouth+90.)/dy
343	call stlmbr(southpolemap,-90.,0.)
344	call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
345	sizesouth,switchsouth,180.)
346	switchsouthg=(switchsouth-ylat0)/dy
347	else
348	sglobal=.false.
349	switchsouthg=999999.
350	endif
351	xauxa=abs(yaux2-90.)
352	if (xglobal.and.xauxa.lt.0.001) then
353	nglobal=.true. ! field contains north pole
354	! Enhance the map scale by factor 3 (*2=6) compared to north-south
355	! map scale
356	sizenorth=6.*(90.-switchnorth)/dy
357	call stlmbr(northpolemap,90.,0.)
358	call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
359	sizenorth,switchnorth,180.)
360	switchnorthg=(switchnorth-ylat0)/dy
361	else
362	nglobal=.false.
363	switchnorthg=999999.
364	endif
365	endif ! ifield.eq.1
366
367	if (nxshift.lt.0) stop 'nxshift (par_mod) must not be negative'
368	if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large'
369
370	! NCEP ISOBARIC LEVELS
371	!*********************
372
373	if( trim(fpname) .eq. 'UU' ) then ! check for U wind
374	iumax=iumax+1
375	pres(iumax) = real(current_grib_level)
376	endif
377
378
379	i179=nint(179./dx)
380	if (dx.lt.0.7) then
381	i180=nint(180./dx)+1 ! 0.5 deg data
382	else
383	i180=nint(179./dx)+1 ! 1 deg data
384	endif
385	i181=i180+1
386
387
388	! NCEP TERRAIN
389	!*************
390
391	if( trim(fpname) .eq. 'ORO' ) then
392	do jy=0,ny-1
393	do ix=0,nxfield-1
394	help=zsec4(nxfield*(ny-jy-1)+ix+1)
395	if(ix.le.i180) then
396	oro(i179+ix,jy)=help
397	excessoro(i179+ix,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
398	else
399	oro(ix-i181,jy)=help
400	excessoro(ix-i181,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
401	endif
402	end do
403	end do
404	endif
405
406	! NCEP LAND SEA MASK
407	!*******************
408
409	if( trim(fpname) .eq. 'LSM' ) then
410	do jy=0,ny-1
411	do ix=0,nxfield-1
412	help=zsec4(nxfield*(ny-jy-1)+ix+1)
413	if(ix.le.i180) then
414	lsm(i179+ix,jy)=help
415	else
416	lsm(ix-i181,jy)=help
417	endif
418	end do
419	end do
420	endif
421
422	call grib_release(igrib)
423
424	goto 10 !! READ NEXT LEVEL OR PARAMETER
425	!
426	! CLOSING OF INPUT DATA FILE
427	!
428
429	! HSO
430	30 continue
431	call grib_close_file(ifile)
432	! HSO end edits
433
434	nuvz=iumax
435	nwz =iumax
436	nlev_ec=iumax
437
438	if (nx.gt.nxmax) then
439	write(,) 'FLEXPART error: Too many grid points in x direction.'
440	write(,) 'Reduce resolution of wind fields.'
441	write(,) 'Or change parameter settings in file par_mod.'
442	write(,) nx,nxmax
443	stop
444	endif
445
446	if (ny.gt.nymax) then
447	write(,) 'FLEXPART error: Too many grid points in y direction.'
448	write(,) 'Reduce resolution of wind fields.'
449	write(,) 'Or change parameter settings in file par_mod.'
450	write(,) ny,nymax
451	stop
452	endif
453
454	if (nuvz.gt.nuvzmax) then
455	write(,) 'FLEXPART error: Too many u,v grid points in z '// &
456	'direction.'
457	write(,) 'Reduce resolution of wind fields.'
458	write(,) 'Or change parameter settings in file par_mod.'
459	write(,) nuvz,nuvzmax
460	stop
461	endif
462
463	if (nwz.gt.nwzmax) then
464	write(,) 'FLEXPART error: Too many w grid points in z '// &
465	'direction.'
466	write(,) 'Reduce resolution of wind fields.'
467	write(,) 'Or change parameter settings in file par_mod.'
468	write(,) nwz,nwzmax
469	stop
470	endif
471
472	! If desired, shift all grids by nxshift grid cells
473	!**************************************************
474
475	if (xglobal) then
476	call shift_field_0(oro,nxfield,ny)
477	call shift_field_0(lsm,nxfield,ny)
478	call shift_field_0(excessoro,nxfield,ny)
479	endif
480
481	! Output of grid info
482	!********************
483
484	write(,)
485	write(,)
486	write(*,'(a,2i7)') '# of vertical levels in NCEP data: ', &
487	nuvz,nwz
488	write(,)
489	write(*,'(a)') 'Mother domain:'
490	write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Longitude range: ', &
491	xlon0,' to ',xlon0+(nx-1)*dx,' Grid distance: ',dx
492	write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Latitude range: ', &
493	ylat0,' to ',ylat0+(ny-1)*dy,' Grid distance: ',dy
494	write(,)
495
496
497	! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL
498	! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
499
500	numskip=nlev_ec-nuvz ! number of ecmwf model layers not used
501	! by trajectory model
502	do i=1,nwz
503	j=numskip+i
504	k=nlev_ec+1+numskip+i
505	akm_usort(nwz-i+1)=pres(nwz-i+1)
506	bkm(nwz-i+1)=0.0
507	end do
508
509	!******************************
510	! change Sabine Eckhardt: akm should always be in descending order ... readwind adapted!
511	!******************************
512	do i=1,nwz
513	if (akm_usort(1).gt.akm_usort(2)) then
514	akm(i)=akm_usort(i)
515	else
516	akm(i)=akm_usort(nwz-i+1)
517	endif
518	end do
519
520	!
521	! CALCULATION OF AKZ, BKZ
522	! AKZ,BKZ: model discretization parameters at the center of each model
523	! layer
524	!
525	! Assign the 10 m winds to an artificial model level with akz=0 and bkz=1.0,
526	! i.e. ground level
527	!*****************************************************************************
528
529	do i=1,nuvz
530	akz(i)=akm(i)
531	bkz(i)=bkm(i)
532	end do
533
534	! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
535	! upon the transformation to z levels. In order to save computer memory, this is
536	! not done anymore in the standard version. However, this option can still be
537	! switched on by replacing the following lines with those below, that are
538	! currently commented out. For this, similar changes are necessary in
539	! verttransform.f and verttranform_nests.f
540	!*****************************************************************************
541
542	nz=nuvz
543	if (nz.gt.nzmax) stop 'nzmax too small'
544	do i=1,nuvz
545	aknew(i)=akz(i)
546	bknew(i)=bkz(i)
547	end do
548
549	! Switch on following lines to use doubled vertical resolution
550	!*************************************************************
551	!nz=nuvz+nwz-1
552	!if (nz.gt.nzmax) stop 'nzmax too small'
553	!do 100 i=1,nwz
554	! aknew(2*(i-1)+1)=akm(i)
555	!00 bknew(2*(i-1)+1)=bkm(i)
556	!do 110 i=2,nuvz
557	! aknew(2*(i-1))=akz(i)
558	!10 bknew(2*(i-1))=bkz(i)
559	! End doubled vertical resolution
560
561
562	! Determine the uppermost level for which the convection scheme shall be applied
563	! by assuming that there is no convection above 50 hPa (for standard SLP)
564	!*****************************************************************************
565
566	do i=1,nuvz-2
567	pint=akz(i)+bkz(i)*101325.
568	if (pint.lt.5000.) goto 96
569	end do
570	96 nconvlev=i
571	if (nconvlev.gt.nconvlevmax-1) then
572	nconvlev=nconvlevmax-1
573	write(,) 'Attention, convection only calculated up to ', &
574	akz(nconvlev)+bkz(nconvlev)*1013.25,' hPa'
575	endif
576
577	return
578
579	999 write(,)
580	write(,) ' ###########################################'// &
581	'###### '
582	write(,) ' TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
583	write(,) ' CAN NOT OPEN INPUT DATA FILE '//wfname(ifn)
584	write(,) ' ###########################################'// &
585	'###### '
586	write(,)
587	write(*,'(a)') '!!! PLEASE INSERT A NEW CD-ROM AND !!!'
588	write(*,'(a)') '!!! PRESS ANY KEY TO CONTINUE... !!!'
589	write(*,'(a)') '!!! ...OR TERMINATE FLEXPART PRESSING!!!'
590	write(*,'(a)') '!!! THE "X" KEY... !!!'
591	write(,)
592	read(*,'(a)') opt
593	if(opt.eq.'X') then
594	stop
595	else
596	goto 5
597	endif
598
599	end subroutine gridcheck_gfs

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

Download in other formats:

Original Format