FLEXPART CTBTO WO8
 All Classes Files Functions Variables
readwind_gfs.F90
Go to the documentation of this file.
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
ew
real function ew(x)
Definition: ew.f90:22
class_vtable::vtable
Definition: class_vtable_mod.F90:53
readwind_gfs
subroutine readwind_gfs(indj, n, uuh, vvh, wwh)
Definition: readwind_gfs.F90:22
class_vtable::vtable_get_fpname
character(len=15) function, public vtable_get_fpname(igrib, vtable_object)
Definition: class_vtable_mod.F90:369
shift_field
subroutine shift_field(field, nxf, nyf, nzfmax, nzf, nmax, n)
Definition: shift_field.f90:22
par_mod
Definition: par_mod.F90:34
pbl_profile
subroutine pbl_profile(ps, td2m, zml1, t2m, tml1, u10m, uml1, stress, hf)
Definition: pbl_profile.f90:22
class_vtable::vtable_load_by_name
subroutine, public vtable_load_by_name(vtable_name, the_vtable_data)
Definition: class_vtable_mod.F90:135
shift_field_0
subroutine shift_field_0(field, nxf, nyf)
Definition: shift_field_0.f90:22
class_vtable::vtable_detect_gribfile_type
integer function, public vtable_detect_gribfile_type(gribfilename)
Definition: class_vtable_mod.F90:65
class_vtable
Definition: class_vtable_mod.F90:2
com_mod
Definition: com_mod.f90:13
Generated on Sun Mar 27 2016 01:00:33 for FLEXPART CTBTO WO8 by   doxygen 1.8.6