Context Navigation

readwind_gfs.f90

bugfixes+enhancements

Last change on this file was dba4221, checked in by Pirmin Kaufmann <pirmin.kaufmann@…>, 18 months ago

Bugfixes:

options/SPECIES/SPECIES_009: corrected wrong number format, replaced comma by decimal point
options/SPECIES/SPECIES_028: corrected wrong number format, moved sign of exponent to after the E
options/SPECIES/specoverview.f90: added namelist parameters that appear in SPECIES files but were missing here
src/FLEXPART.f90: replaced compiler-specific command line argument routines by standard Fortran intrinsic routines
src/FLEXPART_MPI.f90: ditto
src/gridcheck_ecmwf.f90: corrected handling of vertical levels when input files do not contain uppermost layers
src/gridcheck_nests.f90: ditto
src/readwind_ecmwf.f90: corrected handling of vertical levels when input files do not contain uppermost layers
readwind_ecmwf_mpi.f90: ditto

Code enhancements:

options/OUTGRID: added comments describing contents
options/SPECIES/SPECIES_*: aligned comments
options/SPECIES/specoverview.f90: removed commented lines, rectified lines indenting
src/FLEXPART.f90: rectified lines indenting, updated date in version string
src/FLEXPART_MPI.f90: ditto, and realigned code with src/FLEXPART.f90
src/gridcheck_*.f90: added code to write out name of file before it is opened (helps a lot when an input file causes troubles)
src/par_mod.f90: added comment explaining relevance of nuvzmax for GRIB input
src/readreleases.f90: write out warning if too few particles are used to randomize release
src/readspecies.f90: write out name of SPECIES file before it is read
src/readwind_*.f90: write out name of input file before opening it
src/writeheader_txt.f90: removed wrong comment

Property mode set to 100644

File size: 24.9 KB

Line
1	! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
2	! SPDX-License-Identifier: GPL-3.0-or-later
3
4	subroutine readwind_gfs(indj,n,uuh,vvh,wwh)
5
6	!***********************************************************************
7	!* *
8	!* TRAJECTORY MODEL SUBROUTINE READWIND *
9	!* *
10	!***********************************************************************
11	!* *
12	!* AUTHOR: G. WOTAWA *
13	!* DATE: 1997-08-05 *
14	!* LAST UPDATE: 2000-10-17, Andreas Stohl *
15	!* CHANGE: 01/02/2001, Bernd C. Krueger, Variables tth and *
16	!* qvh (on eta coordinates) in common block *
17	!* CHANGE: 16/11/2005, Caroline Forster, GFS data *
18	!* CHANGE: 11/01/2008, Harald Sodemann, Input of GRIB1/2 *
19	!* data with the ECMWF grib_api library *
20	!* CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
21	!* ECMWF grib_api *
22	! *
23	! Unified ECMWF and GFS builds *
24	! Marian Harustak, 12.5.2017 *
25	! - Renamed routine from readwind to readwind_gfs *
26	!* *
27	!***********************************************************************
28	!* *
29	!* DESCRIPTION: *
30	!* *
31	!* READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE *
32	!* INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE *
33	!* *
34	!* INPUT: *
35	!* indj indicates number of the wind field to be read in *
36	!* n temporal index for meteorological fields (1 to 3)*
37	!* *
38	!* IMPORTANT VARIABLES FROM COMMON BLOCK: *
39	!* *
40	!* wfname File name of data to be read in *
41	!* nx,ny,nuvz,nwz expected field dimensions *
42	!* nlev_ec number of vertical levels ecmwf model *
43	!* uu,vv,ww wind fields *
44	!* tt,qv temperature and specific humidity *
45	!* ps surface pressure *
46	!* *
47	!***********************************************************************
48
49	use grib_api
50	use par_mod
51	use com_mod
52
53	implicit none
54
55	!HSO new parameters for grib_api
56	integer :: ifile
57	integer :: iret
58	integer :: igrib
59	integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
60	!HSO end edits
61	real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
62	real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
63	real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
64	integer :: ii,indj,i,j,k,n,levdiff2,ifield,iumax,iwmax
65
66	! NCEP
67	integer :: numpt,numpu,numpv,numpw,numprh,numpclwch
68	real :: help, temp, ew
69	real :: elev
70	real :: ulev1(0:nxmax-1,0:nymax-1),vlev1(0:nxmax-1,0:nymax-1)
71	real :: tlev1(0:nxmax-1,0:nymax-1)
72	real :: qvh2(0:nxmax-1,0:nymax-1)
73
74	integer :: i179,i180,i181
75
76	! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
77	!HSO kept isec1, isec2 and zsec4 for consistency with gribex GRIB input
78
79	integer :: isec1(8),isec2(3)
80	real(kind=4) :: zsec4(jpunp)
81	real(kind=4) :: xaux,yaux,xaux0,yaux0
82	real(kind=8) :: xauxin,yauxin
83	real,parameter :: eps=1.e-4
84	real(kind=4) :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
85	real :: plev1,hlev1,ff10m,fflev1
86
87	logical :: hflswitch,strswitch
88
89	!HSO for grib api error messages
90	character(len=24) :: gribErrorMsg = 'Error reading grib file'
91	character(len=20) :: gribFunction = 'readwind_gfs'
92	character(len=20) :: shortname
93
94
95	hflswitch=.false.
96	strswitch=.false.
97	levdiff2=nlev_ec-nwz+1
98	iumax=0
99	iwmax=0
100
101
102	! OPENING OF DATA FILE (GRIB CODE)
103
104	!HSO
105	write(,) 'Reading: '//path(3)(1:length(3)) &
106	//trim(wfname(indj))
107	call grib_open_file(ifile,path(3)(1:length(3)) &
108	//trim(wfname(indj)),'r',iret)
109	if (iret.ne.GRIB_SUCCESS) then
110	goto 888 ! ERROR DETECTED
111	endif
112	!turn on support for multi fields messages
113	call grib_multi_support_on
114
115	numpt=0
116	numpu=0
117	numpv=0
118	numpw=0
119	numprh=0
120	numpclwch=0
121	ifield=0
122	10 ifield=ifield+1
123	!
124	! GET NEXT FIELDS
125	!
126	call grib_new_from_file(ifile,igrib,iret)
127	if (iret.eq.GRIB_END_OF_FILE) then
128	goto 50 ! EOF DETECTED
129	elseif (iret.ne.GRIB_SUCCESS) then
130	goto 888 ! ERROR DETECTED
131	endif
132
133	!first see if we read GRIB1 or GRIB2
134	call grib_get_int(igrib,'editionNumber',gribVer,iret)
135	! call grib_check(iret,gribFunction,gribErrorMsg)
136
137	if (gribVer.eq.1) then ! GRIB Edition 1
138
139	!read the grib1 identifiers
140	call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
141	! call grib_check(iret,gribFunction,gribErrorMsg)
142	call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
143	! call grib_check(iret,gribFunction,gribErrorMsg)
144	call grib_get_int(igrib,'level',isec1(8),iret)
145	! call grib_check(iret,gribFunction,gribErrorMsg)
146
147	else ! GRIB Edition 2
148
149	!read the grib2 identifiers
150	call grib_get_string(igrib,'shortName',shortname,iret)
151
152	call grib_get_int(igrib,'discipline',discipl,iret)
153	! call grib_check(iret,gribFunction,gribErrorMsg)
154	call grib_get_int(igrib,'parameterCategory',parCat,iret)
155	! call grib_check(iret,gribFunction,gribErrorMsg)
156	call grib_get_int(igrib,'parameterNumber',parNum,iret)
157	! call grib_check(iret,gribFunction,gribErrorMsg)
158	call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
159	! call grib_check(iret,gribFunction,gribErrorMsg)
160	call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
161	valSurf,iret)
162	! call grib_check(iret,gribFunction,gribErrorMsg)
163
164	! write(,) 'Field: ',ifield,parCat,parNum,typSurf,shortname
165	!convert to grib1 identifiers
166	isec1(6)=-1
167	isec1(7)=-1
168	isec1(8)=-1
169	if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T
170	isec1(6)=11 ! indicatorOfParameter
171	isec1(7)=100 ! indicatorOfTypeOfLevel
172	isec1(8)=valSurf/100 ! level, convert to hPa
173	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
174	isec1(6)=33 ! indicatorOfParameter
175	isec1(7)=100 ! indicatorOfTypeOfLevel
176	isec1(8)=valSurf/100 ! level, convert to hPa
177	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.100)) then ! V
178	isec1(6)=34 ! indicatorOfParameter
179	isec1(7)=100 ! indicatorOfTypeOfLevel
180	isec1(8)=valSurf/100 ! level, convert to hPa
181	elseif ((parCat.eq.2).and.(parNum.eq.8).and.(typSurf.eq.100)) then ! W
182	isec1(6)=39 ! indicatorOfParameter
183	isec1(7)=100 ! indicatorOfTypeOfLevel
184	isec1(8)=valSurf/100 ! level, convert to hPa
185	elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.100)) then ! RH
186	isec1(6)=52 ! indicatorOfParameter
187	isec1(7)=100 ! indicatorOfTypeOfLevel
188	isec1(8)=valSurf/100 ! level, convert to hPa
189	elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.103)) then ! RH2
190	isec1(6)=52 ! indicatorOfParameter
191	isec1(7)=105 ! indicatorOfTypeOfLevel
192	isec1(8)=2
193	elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! T2
194	isec1(6)=11 ! indicatorOfParameter
195	isec1(7)=105 ! indicatorOfTypeOfLevel
196	isec1(8)=2
197	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! U10
198	isec1(6)=33 ! indicatorOfParameter
199	isec1(7)=105 ! indicatorOfTypeOfLevel
200	isec1(8)=10
201	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! V10
202	isec1(6)=34 ! indicatorOfParameter
203	isec1(7)=105 ! indicatorOfTypeOfLevel
204	isec1(8)=10
205	elseif ((parCat.eq.1).and.(parNum.eq.22).and.(typSurf.eq.100)) then ! CLWMR Cloud Mixing Ratio [kg/kg]:
206	isec1(6)=153 ! indicatorOfParameter
207	isec1(7)=100 ! indicatorOfTypeOfLevel
208	isec1(8)=valSurf/100 ! level, convert to hPa
209	elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP
210	isec1(6)=2 ! indicatorOfParameter
211	isec1(7)=102 ! indicatorOfTypeOfLevel
212	isec1(8)=0
213	elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then ! SP
214	isec1(6)=1 ! indicatorOfParameter
215	isec1(7)=1 ! indicatorOfTypeOfLevel
216	isec1(8)=0
217	elseif ((parCat.eq.1).and.(parNum.eq.13).and.(typSurf.eq.1)) then ! SNOW
218	isec1(6)=66 ! indicatorOfParameter
219	isec1(7)=1 ! indicatorOfTypeOfLevel
220	isec1(8)=0
221	elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.104)) then ! T sigma 0
222	isec1(6)=11 ! indicatorOfParameter
223	isec1(7)=107 ! indicatorOfTypeOfLevel
224	isec1(8)=0.995 ! lowest sigma level
225	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.104)) then ! U sigma 0
226	isec1(6)=33 ! indicatorOfParameter
227	isec1(7)=107 ! indicatorOfTypeOfLevel
228	isec1(8)=0.995 ! lowest sigma level
229	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.104)) then ! V sigma 0
230	isec1(6)=34 ! indicatorOfParameter
231	isec1(7)=107 ! indicatorOfTypeOfLevel
232	isec1(8)=0.995 ! lowest sigma level
233	elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
234	isec1(6)=7 ! indicatorOfParameter
235	isec1(7)=1 ! indicatorOfTypeOfLevel
236	isec1(8)=0
237	elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
238	.and.(discipl.eq.2)) then ! LSM
239	isec1(6)=81 ! indicatorOfParameter
240	isec1(7)=1 ! indicatorOfTypeOfLevel
241	isec1(8)=0
242	elseif ((parCat.eq.3).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! BLH
243	isec1(6)=221 ! indicatorOfParameter
244	isec1(7)=1 ! indicatorOfTypeOfLevel
245	isec1(8)=0
246	elseif ((parCat.eq.1).and.(parNum.eq.7).and.(typSurf.eq.1)) then ! LSP/TP
247	isec1(6)=62 ! indicatorOfParameter
248	isec1(7)=1 ! indicatorOfTypeOfLevel
249	isec1(8)=0
250	elseif ((parCat.eq.1).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! CP
251	isec1(6)=63 ! indicatorOfParameter
252	isec1(7)=1 ! indicatorOfTypeOfLevel
253	isec1(8)=0
254	endif
255
256	endif ! gribVer
257
258	if (isec1(6).ne.-1) then
259	! get the size and data of the values array
260	call grib_get_real4_array(igrib,'values',zsec4,iret)
261	! call grib_check(iret,gribFunction,gribErrorMsg)
262	endif
263
264	if(ifield.eq.1) then
265
266	!get the required fields from section 2
267	!store compatible to gribex input
268	call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
269	isec2(2),iret)
270	! call grib_check(iret,gribFunction,gribErrorMsg)
271	call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
272	isec2(3),iret)
273	! call grib_check(iret,gribFunction,gribErrorMsg)
274	call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
275	xauxin,iret)
276	! call grib_check(iret,gribFunction,gribErrorMsg)
277	call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
278	yauxin,iret)
279	! call grib_check(iret,gribFunction,gribErrorMsg)
280	xaux=xauxin+real(nxshift)*dx
281	yaux=yauxin
282
283	! CHECK GRID SPECIFICATIONS
284
285	if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
286	if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
287	if(xaux.eq.0.) xaux=-179.0 ! NCEP DATA
288	xaux0=xlon0
289	yaux0=ylat0
290	if(xaux.lt.0.) xaux=xaux+360.
291	if(yaux.lt.0.) yaux=yaux+360.
292	if(xaux0.lt.0.) xaux0=xaux0+360.
293	if(yaux0.lt.0.) yaux0=yaux0+360.
294	if(abs(xaux-xaux0).gt.eps) &
295	stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
296	if(abs(yaux-yaux0).gt.eps) &
297	stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
298	endif
299	!HSO end of edits
300
301	i179=nint(179./dx)
302	if (dx.lt.0.7) then
303	i180=nint(180./dx)+1 ! 0.5 deg data
304	else
305	i180=nint(179./dx)+1 ! 1 deg data
306	endif
307	i181=i180+1
308
309	if (isec1(6).ne.-1) then
310
311	do j=0,nymin1
312	do i=0,nxfield-1
313	if((isec1(6).eq.011).and.(isec1(7).eq.100)) then
314	! TEMPERATURE
315	if((i.eq.0).and.(j.eq.0)) then
316	do ii=1,nuvz
317	if ((isec1(8)*100.0).eq.akz(ii)) numpt=ii
318	end do
319	endif
320	help=zsec4(nxfield*(ny-j-1)+i+1)
321	if(i.le.i180) then
322	tth(i179+i,j,numpt,n)=help
323	else
324	tth(i-i181,j,numpt,n)=help
325	endif
326	endif
327	if((isec1(6).eq.033).and.(isec1(7).eq.100)) then
328	! U VELOCITY
329	if((i.eq.0).and.(j.eq.0)) then
330	do ii=1,nuvz
331	if ((isec1(8)*100.0).eq.akz(ii)) numpu=ii
332	end do
333	endif
334	help=zsec4(nxfield*(ny-j-1)+i+1)
335	if(i.le.i180) then
336	uuh(i179+i,j,numpu)=help
337	else
338	uuh(i-i181,j,numpu)=help
339	endif
340	endif
341	if((isec1(6).eq.034).and.(isec1(7).eq.100)) then
342	! V VELOCITY
343	if((i.eq.0).and.(j.eq.0)) then
344	do ii=1,nuvz
345	if ((isec1(8)*100.0).eq.akz(ii)) numpv=ii
346	end do
347	endif
348	help=zsec4(nxfield*(ny-j-1)+i+1)
349	if(i.le.i180) then
350	vvh(i179+i,j,numpv)=help
351	else
352	vvh(i-i181,j,numpv)=help
353	endif
354	endif
355	if((isec1(6).eq.052).and.(isec1(7).eq.100)) then
356	! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER
357	if((i.eq.0).and.(j.eq.0)) then
358	do ii=1,nuvz
359	if ((isec1(8)*100.0).eq.akz(ii)) numprh=ii
360	end do
361	endif
362	help=zsec4(nxfield*(ny-j-1)+i+1)
363	if(i.le.i180) then
364	qvh(i179+i,j,numprh,n)=help
365	else
366	qvh(i-i181,j,numprh,n)=help
367	endif
368	endif
369	if((isec1(6).eq.001).and.(isec1(7).eq.001)) then
370	! SURFACE PRESSURE
371	help=zsec4(nxfield*(ny-j-1)+i+1)
372	if(i.le.i180) then
373	ps(i179+i,j,1,n)=help
374	else
375	ps(i-i181,j,1,n)=help
376	endif
377	endif
378	if((isec1(6).eq.039).and.(isec1(7).eq.100)) then
379	! W VELOCITY
380	if((i.eq.0).and.(j.eq.0)) then
381	do ii=1,nuvz
382	if ((isec1(8)*100.0).eq.akz(ii)) numpw=ii
383	end do
384	endif
385	help=zsec4(nxfield*(ny-j-1)+i+1)
386	if(i.le.i180) then
387	wwh(i179+i,j,numpw)=help
388	else
389	wwh(i-i181,j,numpw)=help
390	endif
391	endif
392	if((isec1(6).eq.066).and.(isec1(7).eq.001)) then
393	! SNOW DEPTH
394	help=zsec4(nxfield*(ny-j-1)+i+1)
395	if(i.le.i180) then
396	sd(i179+i,j,1,n)=help
397	else
398	sd(i-i181,j,1,n)=help
399	endif
400	endif
401	if((isec1(6).eq.002).and.(isec1(7).eq.102)) then
402	! MEAN SEA LEVEL PRESSURE
403	help=zsec4(nxfield*(ny-j-1)+i+1)
404	if(i.le.i180) then
405	msl(i179+i,j,1,n)=help
406	else
407	msl(i-i181,j,1,n)=help
408	endif
409	endif
410	if((isec1(6).eq.071).and.(isec1(7).eq.244)) then
411	! TOTAL CLOUD COVER
412	help=zsec4(nxfield*(ny-j-1)+i+1)
413	if(i.le.i180) then
414	tcc(i179+i,j,1,n)=help
415	else
416	tcc(i-i181,j,1,n)=help
417	endif
418	endif
419	if((isec1(6).eq.033).and.(isec1(7).eq.105).and. &
420	(isec1(8).eq.10)) then
421	! 10 M U VELOCITY
422	help=zsec4(nxfield*(ny-j-1)+i+1)
423	if(i.le.i180) then
424	u10(i179+i,j,1,n)=help
425	else
426	u10(i-i181,j,1,n)=help
427	endif
428	endif
429	if((isec1(6).eq.034).and.(isec1(7).eq.105).and. &
430	(isec1(8).eq.10)) then
431	! 10 M V VELOCITY
432	help=zsec4(nxfield*(ny-j-1)+i+1)
433	if(i.le.i180) then
434	v10(i179+i,j,1,n)=help
435	else
436	v10(i-i181,j,1,n)=help
437	endif
438	endif
439	if((isec1(6).eq.011).and.(isec1(7).eq.105).and. &
440	(isec1(8).eq.02)) then
441	! 2 M TEMPERATURE
442	help=zsec4(nxfield*(ny-j-1)+i+1)
443	if(i.le.i180) then
444	tt2(i179+i,j,1,n)=help
445	else
446	tt2(i-i181,j,1,n)=help
447	endif
448	endif
449	if((isec1(6).eq.017).and.(isec1(7).eq.105).and. &
450	(isec1(8).eq.02)) then
451	! 2 M DEW POINT TEMPERATURE
452	help=zsec4(nxfield*(ny-j-1)+i+1)
453	if(i.le.i180) then
454	td2(i179+i,j,1,n)=help
455	else
456	td2(i-i181,j,1,n)=help
457	endif
458	endif
459	if((isec1(6).eq.062).and.(isec1(7).eq.001)) then
460	! LARGE SCALE PREC.
461	help=zsec4(nxfield*(ny-j-1)+i+1)
462	if(i.le.i180) then
463	lsprec(i179+i,j,1,n)=help
464	else
465	lsprec(i-i181,j,1,n)=help
466	endif
467	endif
468	if((isec1(6).eq.063).and.(isec1(7).eq.001)) then
469	! CONVECTIVE PREC.
470	help=zsec4(nxfield*(ny-j-1)+i+1)
471	if(i.le.i180) then
472	convprec(i179+i,j,1,n)=help
473	else
474	convprec(i-i181,j,1,n)=help
475	endif
476	endif
477	if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
478	! TOPOGRAPHY
479	help=zsec4(nxfield*(ny-j-1)+i+1)
480	if(i.le.i180) then
481	oro(i179+i,j)=help
482	excessoro(i179+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
483	else
484	oro(i-i181,j)=help
485	excessoro(i-i181,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
486	endif
487	endif
488	if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
489	! LAND SEA MASK
490	help=zsec4(nxfield*(ny-j-1)+i+1)
491	if(i.le.i180) then
492	lsm(i179+i,j)=help
493	else
494	lsm(i-i181,j)=help
495	endif
496	endif
497	if((isec1(6).eq.221).and.(isec1(7).eq.001)) then
498	! MIXING HEIGHT
499	help=zsec4(nxfield*(ny-j-1)+i+1)
500	if(i.le.i180) then
501	hmix(i179+i,j,1,n)=help
502	else
503	hmix(i-i181,j,1,n)=help
504	endif
505	endif
506	if((isec1(6).eq.052).and.(isec1(7).eq.105).and. &
507	(isec1(8).eq.02)) then
508	! 2 M RELATIVE HUMIDITY
509	help=zsec4(nxfield*(ny-j-1)+i+1)
510	if(i.le.i180) then
511	qvh2(i179+i,j)=help
512	else
513	qvh2(i-i181,j)=help
514	endif
515	endif
516	if((isec1(6).eq.011).and.(isec1(7).eq.107)) then
517	! TEMPERATURE LOWEST SIGMA LEVEL
518	help=zsec4(nxfield*(ny-j-1)+i+1)
519	if(i.le.i180) then
520	tlev1(i179+i,j)=help
521	else
522	tlev1(i-i181,j)=help
523	endif
524	endif
525	if((isec1(6).eq.033).and.(isec1(7).eq.107)) then
526	! U VELOCITY LOWEST SIGMA LEVEL
527	help=zsec4(nxfield*(ny-j-1)+i+1)
528	if(i.le.i180) then
529	ulev1(i179+i,j)=help
530	else
531	ulev1(i-i181,j)=help
532	endif
533	endif
534	if((isec1(6).eq.034).and.(isec1(7).eq.107)) then
535	! V VELOCITY LOWEST SIGMA LEVEL
536	help=zsec4(nxfield*(ny-j-1)+i+1)
537	if(i.le.i180) then
538	vlev1(i179+i,j)=help
539	else
540	vlev1(i-i181,j)=help
541	endif
542	endif
543	! SEC & IP 12/2018 read GFS clouds
544	if((isec1(6).eq.153).and.(isec1(7).eq.100)) then !! CLWCR Cloud liquid water content [kg/kg]
545	if((i.eq.0).and.(j.eq.0)) then
546	do ii=1,nuvz
547	if ((isec1(8)*100.0).eq.akz(ii)) numpclwch=ii
548	end do
549	endif
550	help=zsec4(nxfield*(ny-j-1)+i+1)
551	if(i.le.i180) then
552	clwch(i179+i,j,numpclwch,n)=help
553	else
554	clwch(i-i181,j,numpclwch,n)=help
555	endif
556	readclouds=.true.
557	sumclouds=.true.
558	! readclouds=.false.
559	! sumclouds=.false.
560	endif
561
562
563	end do
564	end do
565
566	endif
567
568	if((isec1(6).eq.33).and.(isec1(7).eq.100)) then
569	! NCEP ISOBARIC LEVELS
570	iumax=iumax+1
571	endif
572
573	call grib_release(igrib)
574	goto 10 !! READ NEXT LEVEL OR PARAMETER
575	!
576	! CLOSING OF INPUT DATA FILE
577	!
578
579	!HSO close grib file
580	50 continue
581	call grib_close_file(ifile)
582
583	! SENS. HEAT FLUX
584	sshf(:,:,1,n)=0.0 ! not available from gfs.tccz.pgrbfxx files
585	hflswitch=.false. ! Heat flux not available
586	! SOLAR RADIATIVE FLUXES
587	ssr(:,:,1,n)=0.0 ! not available from gfs.tccz.pgrbfxx files
588	! EW SURFACE STRESS
589	ewss=0.0 ! not available from gfs.tccz.pgrbfxx files
590	! NS SURFACE STRESS
591	nsss=0.0 ! not available from gfs.tccz.pgrbfxx files
592	strswitch=.false. ! stress not available
593
594	! CONVERT TP TO LSP (GRIB2 only)
595	if (gribVer.eq.2) then
596	do j=0,nymin1
597	do i=0,nxfield-1
598	if(i.le.i180) then
599	if (convprec(i179+i,j,1,n).lt.lsprec(i179+i,j,1,n)) then ! neg precip would occur
600	lsprec(i179+i,j,1,n)= &
601	lsprec(i179+i,j,1,n)-convprec(i179+i,j,1,n)
602	else
603	lsprec(i179+i,j,1,n)=0
604	endif
605	else
606	if (convprec(i-i181,j,1,n).lt.lsprec(i-i181,j,1,n)) then
607	lsprec(i-i181,j,1,n)= &
608	lsprec(i-i181,j,1,n)-convprec(i-i181,j,1,n)
609	else
610	lsprec(i-i181,j,1,n)=0
611	endif
612	endif
613	enddo
614	enddo
615	endif
616	!HSO end edits
617
618
619	! TRANSFORM RH TO SPECIFIC HUMIDITY
620
621	do j=0,ny-1
622	do i=0,nxfield-1
623	do k=1,nuvz
624	help=qvh(i,j,k,n)
625	temp=tth(i,j,k,n)
626	plev1=akm(k)+bkm(k)*ps(i,j,1,n)
627	elev=ew(temp)*help/100.0
628	qvh(i,j,k,n)=xmwml(elev/(plev1-((1.0-xmwml)elev)))
629	end do
630	end do
631	end do
632
633	! CALCULATE 2 M DEW POINT FROM 2 M RELATIVE HUMIDITY
634	! USING BOLTON'S (1980) FORMULA
635	! BECAUSE td2 IS NOT AVAILABLE FROM NCEP GFS DATA
636
637	do j=0,ny-1
638	do i=0,nxfield-1
639	help=qvh2(i,j)
640	temp=tt2(i,j,1,n)
641	elev=ew(temp)/100.*help/100. !vapour pressure in hPa
642	td2(i,j,1,n)=243.5/(17.67/log(elev/6.112)-1)+273.
643	if (help.le.0.) td2(i,j,1,n)=tt2(i,j,1,n)
644	end do
645	end do
646
647	if(levdiff2.eq.0) then
648	iwmax=nlev_ec+1
649	do i=0,nxmin1
650	do j=0,nymin1
651	wwh(i,j,nlev_ec+1)=0.
652	end do
653	end do
654	endif
655
656
657	! For global fields, assign the leftmost data column also to the rightmost
658	! data column; if required, shift whole grid by nxshift grid points
659	!*************************************************************************
660
661	if (xglobal) then
662	call shift_field_0(ewss,nxfield,ny)
663	call shift_field_0(nsss,nxfield,ny)
664	call shift_field_0(oro,nxfield,ny)
665	call shift_field_0(excessoro,nxfield,ny)
666	call shift_field_0(lsm,nxfield,ny)
667	call shift_field_0(ulev1,nxfield,ny)
668	call shift_field_0(vlev1,nxfield,ny)
669	call shift_field_0(tlev1,nxfield,ny)
670	call shift_field_0(qvh2,nxfield,ny)
671	call shift_field(ps,nxfield,ny,1,1,2,n)
672	call shift_field(sd,nxfield,ny,1,1,2,n)
673	call shift_field(msl,nxfield,ny,1,1,2,n)
674	call shift_field(tcc,nxfield,ny,1,1,2,n)
675	call shift_field(u10,nxfield,ny,1,1,2,n)
676	call shift_field(v10,nxfield,ny,1,1,2,n)
677	call shift_field(tt2,nxfield,ny,1,1,2,n)
678	call shift_field(td2,nxfield,ny,1,1,2,n)
679	call shift_field(lsprec,nxfield,ny,1,1,2,n)
680	call shift_field(convprec,nxfield,ny,1,1,2,n)
681	call shift_field(sshf,nxfield,ny,1,1,2,n)
682	call shift_field(ssr,nxfield,ny,1,1,2,n)
683	call shift_field(hmix,nxfield,ny,1,1,2,n)
684	call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n)
685	call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n)
686	call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
687	call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
688	call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
689	! IP & SEC adding GFS Clouds 20181205
690	call shift_field(clwch,nxfield,ny,nuvzmax,nuvz,2,n)
691	endif
692
693	do i=0,nxmin1
694	do j=0,nymin1
695	! Convert precip. from mm/s -> mm/hour
696	convprec(i,j,1,n)=convprec(i,j,1,n)*3600.
697	lsprec(i,j,1,n)=lsprec(i,j,1,n)*3600.
698	surfstr(i,j,1,n)=sqrt(ewss(i,j)2+nsss(i,j)2)
699	end do
700	end do
701
702	if ((.not.hflswitch).or.(.not.strswitch)) then
703	! write(,) 'WARNING: No flux data contained in GRIB file ',
704	! + wfname(indj)
705
706	! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
707	!***************************************************************************
708
709	do i=0,nxmin1
710	do j=0,nymin1
711	hlev1=30.0 ! HEIGHT OF FIRST MODEL SIGMA LAYER
712	ff10m= sqrt(u10(i,j,1,n)2+v10(i,j,1,n)2)
713	fflev1=sqrt(ulev1(i,j)2+vlev1(i,j)2)
714	call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
715	tt2(i,j,1,n),tlev1(i,j),ff10m,fflev1, &
716	surfstr(i,j,1,n),sshf(i,j,1,n))
717	if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
718	if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
719	end do
720	end do
721	endif
722
723	if(iumax.ne.nuvz) stop 'READWIND: NUVZ NOT CONSISTENT'
724	if(iumax.ne.nwz) stop 'READWIND: NWZ NOT CONSISTENT'
725
726	return
727	888 write(,) ' #### FLEXPART MODEL ERROR! WINDFIELD #### '
728	write(,) ' #### ',wfname(indj),' #### '
729	write(,) ' #### IS NOT GRIB FORMAT !!! #### '
730	stop 'Execution terminated'
731	999 write(,) ' #### FLEXPART MODEL ERROR! WINDFIELD #### '
732	write(,) ' #### ',wfname(indj),' #### '
733	write(,) ' #### CANNOT BE OPENED !!! #### '
734	stop 'Execution terminated'
735
736	end subroutine readwind_gfs

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

Download in other formats:

Original Format