Context Navigation

readwind_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: 30.7 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 readwind_gfs(indj,n,uuh,vvh,wwh)
23
24	!***********************************************************************
25	!* *
26	!* TRAJECTORY MODEL SUBROUTINE READWIND *
27	!* *
28	!***********************************************************************
29	!* *
30	!* AUTHOR: G. WOTAWA *
31	!* DATE: 1997-08-05 *
32	!* LAST UPDATE: 2000-10-17, Andreas Stohl *
33	!* CHANGE: 01/02/2001, Bernd C. Krueger, Variables tth and *
34	!* qvh (on eta coordinates) in common block *
35	!* CHANGE: 16/11/2005, Caroline Forster, GFS data *
36	!* CHANGE: 11/01/2008, Harald Sodemann, Input of GRIB1/2 *
37	!* data with the ECMWF grib_api library *
38	!* CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
39	!* ECMWF grib_api *
40	!* *
41	!***********************************************************************
42	!* *
43	!* DESCRIPTION: *
44	!* *
45	!* READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE *
46	!* INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE *
47	!* *
48	!***********************************************************************
49	! Changes Arnold, D. and Morton, D. (2015): *
50	! -- description of local and common variables *
51	!***********************************************************************
52
53	use grib_api
54	use par_mod
55	use com_mod
56	use class_vtable
57
58	implicit none
59
60	!***********************************************************************
61	! Subroutine Parameters: *
62	! input *
63	! indj indicates number of the wind field to be read in *
64	! n temporal index for meteorological fields (1 to 3)*
65	! uuh,vvh, wwh wind components in ECMWF model levels *
66	integer :: indj,n
67	real(kind=4) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
68	real(kind=4) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
69	real(kind=4) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
70	!***********************************************************************
71
72	!***********************************************************************
73	! Local variables: *
74	! *
75	! HSO variables for grib_api: *
76	! ifile grib file to be opened and read in *
77	! iret is a return code for successful or not open *
78	! igrib grib edition number (whether GRIB1 or GRIB2) *
79	! gribVer where info on igrib is kept, 1 for GRIB1 and *
80	! 2 for GRIB2 *
81	! parCat parameter category ( = number) , how FLEXPART *
82	! identifies fields *
83	! parNum parameter number by product discipline and *
84	! parameter category *
85	! typSurf type of first fixed surface *
86	! valSurf level *
87	! discipl discipline of processed data contained in a *
88	! GRIB message *
89	integer :: ifile
90	integer :: iret
91	integer :: igrib
92	integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
93	! *
94	! NCEP specific variables: *
95	! numpt number of pressure levels for temperature *
96	! numpu number of pressure levels for U velocity *
97	! numpv number of pressure levels for V velocity *
98	! numpw number of pressure levels for W velocity *
99	! numprh number of pressure levels for relative humidity *
100	! help temporary variable to store fields *
101	! temp temporary variable to store tth and tt2 variables*
102	! ew function to calculate Goff-Gratch formula, *
103	! the saturation water vapor pressure *
104	! elev elevation in meters *
105	! ulev1 U velocity at the lowest sigma level *
106	! vlev1 V velocity at the lowest sigma level *
107	! tlev1 temperature at the lowest sigma level *
108	! qvh2 2m dew point temperature *
109	! i179 number of pints in x *
110	! i180 i179 +1, x direction +1 *
111	! i181 i180 +1 *
112	integer :: numpt,numpu,numpv,numpw,numprh
113	real :: help, temp, ew
114	real :: elev
115	real :: ulev1(0:nxmax-1,0:nymax-1),vlev1(0:nxmax-1,0:nymax-1)
116	real :: tlev1(0:nxmax-1,0:nymax-1)
117	real :: qvh2(0:nxmax-1,0:nymax-1)
118	integer :: i179,i180,i181
119
120
121	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
122	!!!! Vtable related variables
123
124	! Paths to Vtables (current implementation will assume they are in the cwd
125	! This is a crappy place for these parameters. Need to move them.
126	character(LEN=255), parameter :: VTABLE_NCEP_GRIB1_PATH = &
127	"Vtable_ncep_grib1", &
128	VTABLE_NCEP_GRIB2_PATH = &
129	"Vtable_ncep_grib2"
130
131
132
133	integer :: gribfile_type
134	integer :: current_grib_level ! This "was" isec1(8) in previous version
135	character(len=255) :: gribfile_name
136	character(len=255) :: vtable_path
137	character(len=15) :: fpname
138
139	type(Vtable) :: my_vtable ! unallocated
140	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
141
142
143
144
145
146
147
148
149
150	! *
151	! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING *
152	!HSO kept isec1, isec2 and zsec4 for consistency with gribex GRIB input*
153	!
154	integer :: isec2(3)
155	real(kind=4) :: zsec4(jpunp)
156	real(kind=4) :: xaux,yaux,xaux0,yaux0
157	real(kind=8) :: xauxin,yauxin
158	real(kind=4) :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
159	real :: plev1,hlev1,ff10m,fflev1
160	logical :: hflswitch,strswitch
161	! *
162	! Other variables: *
163	! i,j,k loop control indices in each direction *
164	! levdiff2 number of model levels - number of model levels *
165	! for the staggered vertical wind +1 *
166	! ifield index to control field read in *
167	! iumax *
168	! iwmax *
169	integer :: ii,i,j,k,levdiff2,ifield,iumax,iwmax
170	!***********************************************************************
171
172	!***********************************************************************
173	! Local constants *
174	real,parameter :: eps=1.e-4
175	!HSO grib api error messages:
176	character(len=24) :: gribErrorMsg = 'Error reading grib file'
177	character(len=20) :: gribFunction = 'readwind_gfs'
178	!***********************************************************************
179
180	!***********************************************************************
181	! Global variables *
182	! from par_mod and com_mod *
183	! wfname File name of data to be read in *
184	! nx,ny,nuvz,nwz expected field dimensions *
185	! nxfield same as nx but for limited area fields *
186	! nxmin1,nymin1 nx-1, ny-1, respectively *
187	! nxmax,nymax maximum dimension of wind fields in x and y *
188	! direction, respectively *
189	! nuvzmax,nwzmax maximum dimension of (u,v) and (w) wind fields in z
190	! direction (for fields on eta levels) *
191	! nlev_ec number of vertical levels ecmwf model *
192	! u10,v10 wind fields at 10 m *
193	! tt2,td2 temperature and dew point temperature at 2 m *
194	! tth,qvh temperature and specific humidity on original *
195	! eta levels *
196	! ps surface pressure *
197	! sd snow depth (but not used afterwards!) *
198	! msl mean sea level pressure *
199	! ttc total cloud cover *
200	! lsprec large scale precipitation *
201	! convprec convective precipitation *
202	! sshf sensible heat flux *
203	! ssr solar radiation *
204	! surfstr surface stress *
205	! oro orography *
206	! excessoro standard deviation of orography *
207	! lsm land sea mask *
208	! hmix mixing height *
209	!***********************************************************************
210
211	!-----------------------------------------------------------------------------
212
213
214
215
216
217
218
219
220	hflswitch=.false.
221	strswitch=.false.
222	levdiff2=nlev_ec-nwz+1
223	iumax=0
224	iwmax=0
225
226
227
228	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
229	!!!!!!!!!!!!!!!!!!! VTABLE code
230	!!!!!!! Vtable choice
231	gribfile_name = path(3)(1:length(3))//trim(wfname(indj))
232	print *, 'gribfile_name: ', gribfile_name
233
234	gribfile_type = vtable_detect_gribfile_type( gribfile_name )
235
236	print *, 'gribfile_type: ', gribfile_type
237
238	if (gribfile_type .eq. VTABLE_GRIBFILE_TYPE_NCEP_GRIB1) then
239	vtable_path = VTABLE_NCEP_GRIB1_PATH
240	else if (gribfile_type .eq. VTABLE_GRIBFILE_TYPE_NCEP_GRIB2) then
241	vtable_path = VTABLE_NCEP_GRIB2_PATH
242	else
243	print *, 'Unsupported gribfile_type: ', gribfile_type
244	stop
245	endif
246
247
248	! Load the Vtable into 'my_vtable'
249	print *, 'Loading Vtable: ', vtable_path
250	call vtable_load_by_name(vtable_path, my_vtable)
251	print *, 'Vtable Initialized: ', my_vtable%initialized
252	print *, 'Vtable num_entries: ', my_vtable%num_entries
253	!!!!!!!!!!!!!!!!!!! VTABLE code
254	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
255
256
257
258	! OPENING OF DATA FILE (GRIB CODE)
259
260	!HSO
261	5 call grib_open_file(ifile,path(3)(1:length(3)) &
262	//trim(wfname(indj)),'r',iret)
263	if (iret.ne.GRIB_SUCCESS) then
264	goto 888 ! ERROR DETECTED
265	endif
266	!turn on support for multi fields messages
267	call grib_multi_support_on
268
269	numpt=0
270	numpu=0
271	numpv=0
272	numpw=0
273	numprh=0
274	ifield=0
275	10 ifield=ifield+1
276	!
277	! GET NEXT FIELDS
278	!
279	call grib_new_from_file(ifile,igrib,iret)
280	if (iret.eq.GRIB_END_OF_FILE) then
281	goto 50 ! EOF DETECTED
282	elseif (iret.ne.GRIB_SUCCESS) then
283	goto 888 ! ERROR DETECTED
284	endif
285
286
287	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
288	!!!!!!!!!!!!!!!!!!! VTABLE code
289	! Get the fpname
290	fpname = vtable_get_fpname(igrib, my_vtable)
291	!print *, 'fpname: ', trim(fpname)
292
293
294	!!!!!!!!!!!!!!!!!!! VTABLE code
295	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
296
297
298
299	!first see if we read GRIB1 or GRIB2
300	call grib_get_int(igrib,'editionNumber',gribVer,iret)
301	call grib_check(iret,gribFunction,gribErrorMsg)
302
303	if (gribVer.eq.1) then ! GRIB Edition 1
304
305	call grib_get_int(igrib,'level', current_grib_level, iret)
306	call grib_check(iret,gribFunction,gribErrorMsg)
307
308	!!! Added by DJM 2016-03-02 - if this is GRIB1 we assume that
309	!!! level units are hPa and need to be multiplied by 100 for Pa
310	current_grib_level = current_grib_level*100.0
311
312	else ! GRIB Edition 2
313
314	call grib_get_int(igrib, 'scaledValueOfFirstFixedSurface', &
315	current_grib_level, iret)
316	call grib_check(iret,gribFunction,gribErrorMsg)
317
318	endif ! gribVer
319
320	if (trim(fpname) .ne. 'None') then
321	! get the size and data of the values array
322	call grib_get_real4_array(igrib,'values',zsec4,iret)
323	call grib_check(iret,gribFunction,gribErrorMsg)
324	endif
325
326	if(ifield.eq.1) then
327
328	!get the required fields from section 2
329	!store compatible to gribex input
330	call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
331	isec2(2),iret)
332	call grib_check(iret,gribFunction,gribErrorMsg)
333	call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
334	isec2(3),iret)
335	call grib_check(iret,gribFunction,gribErrorMsg)
336	call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
337	xauxin,iret)
338	call grib_check(iret,gribFunction,gribErrorMsg)
339	call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
340	yauxin,iret)
341	call grib_check(iret,gribFunction,gribErrorMsg)
342	xaux=xauxin+real(nxshift)*dx
343	!print *, 'xauxin: ', xauxin
344	!print *, 'nxshift: ', nxshift
345	!print *, 'dx: ', dx
346	!print *, 'xlon0: ', xlon0
347	yaux=yauxin
348
349	! CHECK GRID SPECIFICATIONS
350
351	if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
352	if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
353
354	if(xaux.eq.0.) xaux=-179.0 ! NCEP GRIB DATA
355	xaux0=xlon0
356	yaux0=ylat0
357	if(xaux.lt.0.) xaux=xaux+360.
358	if(yaux.lt.0.) yaux=yaux+360.
359	if(xaux0.lt.0.) xaux0=xaux0+360.
360	if(yaux0.lt.0.) yaux0=yaux0+360.
361	if(abs(xaux-xaux0).gt.eps) &
362	stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
363	if(abs(yaux-yaux0).gt.eps) &
364	stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
365	endif
366	!HSO end of edits
367
368	i179=nint(179./dx)
369	if (dx.lt.0.7) then
370	i180=nint(180./dx)+1 ! 0.5 deg data
371	else
372	i180=nint(179./dx)+1 ! 1 deg data
373	endif
374	i181=i180+1
375
376	!!!! if (trim(fpname) .ne. 'None') then
377
378
379	! TEMPERATURE
380	if( trim(fpname) .eq. 'TT' ) then
381	do j=0,nymin1
382	do i=0,nxfield-1
383	if((i.eq.0).and.(j.eq.0)) then
384	do ii=1,nuvz
385	if (current_grib_level .eq. akz(ii)) numpt=ii
386	end do
387	endif
388	help=zsec4(nxfield*(ny-j-1)+i+1)
389	if(i.le.i180) then
390	tth(i179+i,j,numpt,n)=help
391	else
392	tth(i-i181,j,numpt,n)=help
393	endif
394	end do
395	end do
396	endif
397
398
399	! U VELOCITY
400	if( trim(fpname) .eq. 'UU' ) then
401	do j=0,nymin1
402	do i=0,nxfield-1
403	if((i.eq.0).and.(j.eq.0)) then
404	do ii=1,nuvz
405	if (current_grib_level .eq. akz(ii)) numpu=ii
406	end do
407	endif
408	help=zsec4(nxfield*(ny-j-1)+i+1)
409	if(i.le.i180) then
410	uuh(i179+i,j,numpu)=help
411	else
412	uuh(i-i181,j,numpu)=help
413	endif
414	end do
415	end do
416	endif
417
418	! V VELOCITY
419	if( trim(fpname) .eq. 'VV' ) then
420	do j=0,nymin1
421	do i=0,nxfield-1
422	if((i.eq.0).and.(j.eq.0)) then
423	do ii=1,nuvz
424	if (current_grib_level .eq. akz(ii)) numpv=ii
425	end do
426	endif
427	help=zsec4(nxfield*(ny-j-1)+i+1)
428	if(i.le.i180) then
429	vvh(i179+i,j,numpv)=help
430	else
431	vvh(i-i181,j,numpv)=help
432	endif
433	end do
434	end do
435	endif
436
437
438	! W VELOCITY
439	if( trim(fpname) .eq. 'WW' ) then
440	do j=0,nymin1
441	do i=0,nxfield-1
442	if((i.eq.0).and.(j.eq.0)) then
443	do ii=1,nuvz
444	if (current_grib_level .eq. akz(ii)) numpw=ii
445	end do
446	endif
447	help=zsec4(nxfield*(ny-j-1)+i+1)
448	if(i.le.i180) then
449	wwh(i179+i,j,numpw)=help
450	else
451	wwh(i-i181,j,numpw)=help
452	endif
453	end do
454	end do
455	endif
456
457
458	! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER
459	if( trim(fpname) .eq. 'RH' ) then
460	do j=0,nymin1
461	do i=0,nxfield-1
462	if((i.eq.0).and.(j.eq.0)) then
463	do ii=1,nuvz
464	if (current_grib_level .eq. akz(ii)) numprh=ii
465	end do
466	endif
467	help=zsec4(nxfield*(ny-j-1)+i+1)
468	if(i.le.i180) then
469	qvh(i179+i,j,numprh,n)=help
470	else
471	qvh(i-i181,j,numprh,n)=help
472	endif
473	end do
474	end do
475	endif
476
477
478	! SURFACE PRESSURE
479	if( trim(fpname) .eq. 'PS' ) then
480	do j=0,nymin1
481	do i=0,nxfield-1
482	help=zsec4(nxfield*(ny-j-1)+i+1)
483	if(i.le.i180) then
484	ps(i179+i,j,1,n)=help
485	else
486	ps(i-i181,j,1,n)=help
487	endif
488	end do
489	end do
490	endif
491
492
493	! SNOW DEPTH
494	if( trim(fpname) .eq. 'SD' ) then
495	do j=0,nymin1
496	do i=0,nxfield-1
497	help=zsec4(nxfield*(ny-j-1)+i+1)
498	if(i.le.i180) then
499	sd(i179+i,j,1,n)=help
500	else
501	sd(i-i181,j,1,n)=help
502	endif
503	end do
504	end do
505	endif
506
507
508	! MEAN SEA LEVEL PRESSURE
509	if( trim(fpname) .eq. 'SLP' ) then
510	do j=0,nymin1
511	do i=0,nxfield-1
512	help=zsec4(nxfield*(ny-j-1)+i+1)
513	if(i.le.i180) then
514	msl(i179+i,j,1,n)=help
515	else
516	msl(i-i181,j,1,n)=help
517	endif
518	end do
519	end do
520	endif
521
522
523	! TOTAL CLOUD COVER
524	if( trim(fpname) .eq. 'TCC' ) then
525	do j=0,nymin1
526	do i=0,nxfield-1
527	help=zsec4(nxfield*(ny-j-1)+i+1)
528	if(i.le.i180) then
529	tcc(i179+i,j,1,n)=help
530	else
531	tcc(i-i181,j,1,n)=help
532	endif
533	end do
534	end do
535	endif
536
537	! 10 M U VELOCITY
538	if( trim(fpname) .eq. 'U10' ) then
539	do j=0,nymin1
540	do i=0,nxfield-1
541	help=zsec4(nxfield*(ny-j-1)+i+1)
542	if(i.le.i180) then
543	u10(i179+i,j,1,n)=help
544	else
545	u10(i-i181,j,1,n)=help
546	endif
547	end do
548	end do
549	endif
550
551	! 10 M V VELOCITY
552	if( trim(fpname) .eq. 'V10' ) then
553	do j=0,nymin1
554	do i=0,nxfield-1
555	help=zsec4(nxfield*(ny-j-1)+i+1)
556	if(i.le.i180) then
557	v10(i179+i,j,1,n)=help
558	else
559	v10(i-i181,j,1,n)=help
560	endif
561	end do
562	end do
563	endif
564
565
566	! 2 M TEMPERATURE
567	if( trim(fpname) .eq. 'T2' ) then
568	do j=0,nymin1
569	do i=0,nxfield-1
570	help=zsec4(nxfield*(ny-j-1)+i+1)
571	if(i.le.i180) then
572	tt2(i179+i,j,1,n)=help
573	else
574	tt2(i-i181,j,1,n)=help
575	endif
576	end do
577	end do
578	endif
579
580	! 2 M DEW POINT TEMPERATURE
581	if( trim(fpname) .eq. 'TD2' ) then
582	do j=0,nymin1
583	do i=0,nxfield-1
584	help=zsec4(nxfield*(ny-j-1)+i+1)
585	if(i.le.i180) then
586	td2(i179+i,j,1,n)=help
587	else
588	td2(i-i181,j,1,n)=help
589	endif
590	end do
591	end do
592	endif
593
594
595	! LARGE SCALE PREC.
596	if( trim(fpname) .eq. 'LSPREC' ) then
597	do j=0,nymin1
598	do i=0,nxfield-1
599	help=zsec4(nxfield*(ny-j-1)+i+1)
600	if(i.le.i180) then
601	lsprec(i179+i,j,1,n)=help
602	else
603	lsprec(i-i181,j,1,n)=help
604	endif
605	end do
606	end do
607	endif
608
609
610	! CONVECTIVE PREC.
611	if( trim(fpname) .eq. 'CONVPREC' ) then
612	do j=0,nymin1
613	do i=0,nxfield-1
614	help=zsec4(nxfield*(ny-j-1)+i+1)
615	if(i.le.i180) then
616	convprec(i179+i,j,1,n)=help
617	else
618	convprec(i-i181,j,1,n)=help
619	endif
620	end do
621	end do
622	endif
623
624
625	! TOPOGRAPHY
626	if( trim(fpname) .eq. 'ORO' ) then
627	do j=0,nymin1
628	do i=0,nxfield-1
629	help=zsec4(nxfield*(ny-j-1)+i+1)
630	if(i.le.i180) then
631	oro(i179+i,j)=help
632	excessoro(i179+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID
633	! TERRAIN DISREGARDED
634	else
635	oro(i-i181,j)=help
636	excessoro(i-i181,j)=0.0 ! ISOBARIC SURFACES: SUBGRID
637	! TERRAIN DISREGARDED
638	endif
639	end do
640	end do
641	endif
642
643	! LAND SEA MASK
644	if( trim(fpname) .eq. 'LSM' ) then
645	do j=0,nymin1
646	do i=0,nxfield-1
647	help=zsec4(nxfield*(ny-j-1)+i+1)
648	if(i.le.i180) then
649	lsm(i179+i,j)=help
650	else
651	lsm(i-i181,j)=help
652	endif
653	end do
654	end do
655	endif
656
657
658	! MIXING HEIGHT
659	if( trim(fpname) .eq. 'HMIX' ) then
660	do j=0,nymin1
661	do i=0,nxfield-1
662	help=zsec4(nxfield*(ny-j-1)+i+1)
663	if(i.le.i180) then
664	hmix(i179+i,j,1,n)=help
665	else
666	hmix(i-i181,j,1,n)=help
667	endif
668	end do
669	end do
670	endif
671
672
673	! 2 M RELATIVE HUMIDITY
674	if( trim(fpname) .eq. 'RH2' ) then
675	do j=0,nymin1
676	do i=0,nxfield-1
677	help=zsec4(nxfield*(ny-j-1)+i+1)
678	if(i.le.i180) then
679	qvh2(i179+i,j)=help
680	else
681	qvh2(i-i181,j)=help
682	endif
683	end do
684	end do
685	endif
686
687
688	! TEMPERATURE LOWEST SIGMA LEVEL
689	if( trim(fpname) .eq. 'TSIG1' ) then
690	do j=0,nymin1
691	do i=0,nxfield-1
692	help=zsec4(nxfield*(ny-j-1)+i+1)
693	if(i.le.i180) then
694	tlev1(i179+i,j)=help
695	else
696	tlev1(i-i181,j)=help
697	endif
698	end do
699	end do
700	endif
701
702
703	! U VELOCITY LOWEST SIGMA LEVEL
704	if( trim(fpname) .eq. 'USIG1' ) then
705	do j=0,nymin1
706	do i=0,nxfield-1
707	help=zsec4(nxfield*(ny-j-1)+i+1)
708	if(i.le.i180) then
709	ulev1(i179+i,j)=help
710	else
711	ulev1(i-i181,j)=help
712	endif
713	end do
714	end do
715	endif
716
717	! V VELOCITY LOWEST SIGMA LEVEL
718	if( trim(fpname) .eq. 'VSIG1' ) then
719	do j=0,nymin1
720	do i=0,nxfield-1
721	help=zsec4(nxfield*(ny-j-1)+i+1)
722	if(i.le.i180) then
723	vlev1(i179+i,j)=help
724	else
725	vlev1(i-i181,j)=help
726	endif
727	end do
728	end do
729	endif
730
731
732
733
734
735	!!!! endif
736
737	if( trim(fpname) .eq. 'UU') then
738	! NCEP ISOBARIC LEVELS - this counts the number of pressure levels
739	! by counting the number of occurrences of UU, which us u-wind on pressure levels
740	iumax=iumax+1
741	endif
742
743	call grib_release(igrib)
744	goto 10 !! READ NEXT LEVEL OR PARAMETER
745	!
746	! CLOSING OF INPUT DATA FILE
747	!
748
749	!HSO close grib file
750	50 continue
751	call grib_close_file(ifile)
752
753	! SENS. HEAT FLUX
754	sshf(:,:,1,n)=0.0 ! not available from gfs.tccz.pgrbfxx files
755	hflswitch=.false. ! Heat flux not available
756	! SOLAR RADIATIVE FLUXES
757	ssr(:,:,1,n)=0.0 ! not available from gfs.tccz.pgrbfxx files
758	! EW SURFACE STRESS
759	ewss=0.0 ! not available from gfs.tccz.pgrbfxx files
760	! NS SURFACE STRESS
761	nsss=0.0 ! not available from gfs.tccz.pgrbfxx files
762	strswitch=.false. ! stress not available
763
764	! CONVERT TP TO LSP (GRIB2 only)
765	if (gribVer.eq.2) then
766	do j=0,nymin1
767	do i=0,nxfield-1
768	if(i.le.i180) then
769	if (convprec(i179+i,j,1,n).lt.lsprec(i179+i,j,1,n)) then ! neg precip would occur
770	lsprec(i179+i,j,1,n)= &
771	lsprec(i179+i,j,1,n)-convprec(i179+i,j,1,n)
772	else
773	lsprec(i179+i,j,1,n)=0
774	endif
775	else
776	if (convprec(i-i181,j,1,n).lt.lsprec(i-i181,j,1,n)) then
777	lsprec(i-i181,j,1,n)= &
778	lsprec(i-i181,j,1,n)-convprec(i-i181,j,1,n)
779	else
780	lsprec(i-i181,j,1,n)=0
781	endif
782	endif
783	enddo
784	enddo
785	endif
786	!HSO end edits
787
788
789	! TRANSFORM RH TO SPECIFIC HUMIDITY
790
791	do j=0,ny-1
792	do i=0,nxfield-1
793	do k=1,nuvz
794	help=qvh(i,j,k,n)
795	temp=tth(i,j,k,n)
796	plev1=akm(k)+bkm(k)*ps(i,j,1,n)
797	elev=ew(temp)*help/100.0
798	qvh(i,j,k,n)=xmwml(elev/(plev1-((1.0-xmwml)elev)))
799	end do
800	end do
801	end do
802
803	! CALCULATE 2 M DEW POINT FROM 2 M RELATIVE HUMIDITY
804	! USING BOLTON'S (1980) FORMULA
805	! BECAUSE td2 IS NOT AVAILABLE FROM NCEP GFS DATA
806
807	do j=0,ny-1
808	do i=0,nxfield-1
809	help=qvh2(i,j)
810	temp=tt2(i,j,1,n)
811	elev=ew(temp)/100.*help/100. !vapour pressure in hPa
812	td2(i,j,1,n)=243.5/(17.67/log(elev/6.112)-1)+273.
813	if (help.le.0.) td2(i,j,1,n)=tt2(i,j,1,n)
814	end do
815	end do
816
817	if(levdiff2.eq.0) then
818	iwmax=nlev_ec+1
819	do i=0,nxmin1
820	do j=0,nymin1
821	wwh(i,j,nlev_ec+1)=0.
822	end do
823	end do
824	endif
825
826
827	! For global fields, assign the leftmost data column also to the rightmost
828	! data column; if required, shift whole grid by nxshift grid points
829	!*************************************************************************
830
831	if (xglobal) then
832	call shift_field_0(ewss,nxfield,ny)
833	call shift_field_0(nsss,nxfield,ny)
834	call shift_field_0(oro,nxfield,ny)
835	call shift_field_0(excessoro,nxfield,ny)
836	call shift_field_0(lsm,nxfield,ny)
837	call shift_field_0(ulev1,nxfield,ny)
838	call shift_field_0(vlev1,nxfield,ny)
839	call shift_field_0(tlev1,nxfield,ny)
840	call shift_field_0(qvh2,nxfield,ny)
841	call shift_field(ps,nxfield,ny,1,1,2,n)
842	call shift_field(sd,nxfield,ny,1,1,2,n)
843	call shift_field(msl,nxfield,ny,1,1,2,n)
844	call shift_field(tcc,nxfield,ny,1,1,2,n)
845	call shift_field(u10,nxfield,ny,1,1,2,n)
846	call shift_field(v10,nxfield,ny,1,1,2,n)
847	call shift_field(tt2,nxfield,ny,1,1,2,n)
848	call shift_field(td2,nxfield,ny,1,1,2,n)
849	call shift_field(lsprec,nxfield,ny,1,1,2,n)
850	call shift_field(convprec,nxfield,ny,1,1,2,n)
851	call shift_field(sshf,nxfield,ny,1,1,2,n)
852	call shift_field(ssr,nxfield,ny,1,1,2,n)
853	call shift_field(hmix,nxfield,ny,1,1,2,n)
854	call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n)
855	call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n)
856	call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
857	call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
858	call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
859	endif
860
861	do i=0,nxmin1
862	do j=0,nymin1
863	! Convert precip. from mm/s -> mm/hour
864	convprec(i,j,1,n)=convprec(i,j,1,n)*3600.
865	lsprec(i,j,1,n)=lsprec(i,j,1,n)*3600.
866	surfstr(i,j,1,n)=sqrt(ewss(i,j)2+nsss(i,j)2)
867	end do
868	end do
869
870	if ((.not.hflswitch).or.(.not.strswitch)) then
871	! write(,) 'WARNING: No flux data contained in GRIB file ',
872	! + wfname(indj)
873
874	! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
875	!***************************************************************************
876
877	do i=0,nxmin1
878	do j=0,nymin1
879	hlev1=30.0 ! HEIGHT OF FIRST MODEL SIGMA LAYER
880	ff10m= sqrt(u10(i,j,1,n)2+v10(i,j,1,n)2)
881	fflev1=sqrt(ulev1(i,j)2+vlev1(i,j)2)
882	call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
883	tt2(i,j,1,n),tlev1(i,j),ff10m,fflev1, &
884	surfstr(i,j,1,n),sshf(i,j,1,n))
885	if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
886	if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
887	end do
888	end do
889	endif
890
891
892	if(iumax.ne.nuvz) stop 'READWIND: NUVZ NOT CONSISTENT'
893	if(iumax.ne.nwz) stop 'READWIND: NWZ NOT CONSISTENT'
894
895	return
896	888 write(,) ' #### FLEXPART MODEL ERROR! WINDFIELD #### '
897	write(,) ' #### ',wfname(indj),' #### '
898	write(,) ' #### IS NOT GRIB FORMAT !!! #### '
899	stop 'Execution terminated'
900	999 write(,) ' #### FLEXPART MODEL ERROR! WINDFIELD #### '
901	write(,) ' #### ',wfname(indj),' #### '
902	write(,) ' #### CANNOT BE OPENED !!! #### '
903	stop 'Execution terminated'
904
905	end subroutine readwind_gfs

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

Download in other formats:

Original Format