CTBTO FLEXPART WO4 (2015-10-15)
 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 
57  implicit none
58 
59  !***********************************************************************
60  ! Subroutine Parameters: *
61  ! input *
62  ! indj indicates number of the wind field to be read in *
63  ! n temporal index for meteorological fields (1 to 3)*
64  ! uuh,vvh, wwh wind components in ECMWF model levels *
65  integer :: indj,n
66  real(kind=4) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
67  real(kind=4) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
68  real(kind=4) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
69  !***********************************************************************
70 
71  !***********************************************************************
72  ! Local variables: *
73  ! *
74  ! HSO variables for grib_api: *
75  ! ifile grib file to be opened and read in *
76  ! iret is a return code for successful or not open *
77  ! igrib grib edition number (whether GRIB1 or GRIB2) *
78  ! gribVer where info on igrib is kept, 1 for GRIB1 and *
79  ! 2 for GRIB2 *
80  ! parCat parameter category ( = number) , how FLEXPART *
81  ! identifies fields *
82  ! parNum parameter number by product discipline and *
83  ! parameter category *
84  ! typSurf type of first fixed surface *
85  ! valSurf level *
86  ! discipl discipline of processed data contained in a *
87  ! GRIB message *
88  integer :: ifile
89  integer :: iret
90  integer :: igrib
91  integer :: gribver,parcat,parnum,typsurf,valsurf,discipl
92  ! *
93  ! NCEP specific variables: *
94  ! numpt number of pressure levels for temperature *
95  ! numpu number of pressure levels for U velocity *
96  ! numpv number of pressure levels for V velocity *
97  ! numpw number of pressure levels for W velocity *
98  ! numprh number of pressure levels for relative humidity *
99  ! help temporary variable to store fields *
100  ! temp temporary variable to store tth and tt2 variables*
101  ! ew function to calculate Goff-Gratch formula, *
102  ! the saturation water vapor pressure *
103  ! elev elevation in meters *
104  ! ulev1 U velocity at the lowest sigma level *
105  ! vlev1 V velocity at the lowest sigma level *
106  ! tlev1 temperature at the lowest sigma level *
107  ! qvh2 2m dew point temperature *
108  ! i179 number of pints in x *
109  ! i180 i179 +1, x direction +1 *
110  ! i181 i180 +1 *
111  integer :: numpt,numpu,numpv,numpw,numprh
112  real :: help, temp, ew
113  real :: elev
114  real :: ulev1(0:nxmax-1,0:nymax-1),vlev1(0:nxmax-1,0:nymax-1)
115  real :: tlev1(0:nxmax-1,0:nymax-1)
116  real :: qvh2(0:nxmax-1,0:nymax-1)
117  integer :: i179,i180,i181
118  ! *
119  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING *
120  !HSO kept isec1, isec2 and zsec4 for consistency with gribex GRIB input*
121  !
122  integer :: isec1(8),isec2(3)
123  real(kind=4) :: zsec4(jpunp)
124  real(kind=4) :: xaux,yaux,xaux0,yaux0
125  real(kind=8) :: xauxin,yauxin
126  real(kind=4) :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
127  real :: plev1,hlev1,ff10m,fflev1
128  logical :: hflswitch,strswitch
129  ! *
130  ! Other variables: *
131  ! i,j,k loop control indices in each direction *
132  ! levdiff2 number of model levels - number of model levels *
133  ! for the staggered vertical wind +1 *
134  ! ifield index to control field read in *
135  ! iumax *
136  ! iwmax *
137  integer :: ii,i,j,k,levdiff2,ifield,iumax,iwmax
138  !***********************************************************************
139 
140  !***********************************************************************
141  ! Local constants *
142  real,parameter :: eps=1.e-4
143  !HSO grib api error messages:
144  character(len=24) :: griberrormsg = 'Error reading grib file'
145  character(len=20) :: gribfunction = 'readwind_gfs'
146  !***********************************************************************
147 
148  !***********************************************************************
149  ! Global variables *
150  ! from par_mod and com_mod *
151  ! wfname File name of data to be read in *
152  ! nx,ny,nuvz,nwz expected field dimensions *
153  ! nxfield same as nx but for limited area fields *
154  ! nxmin1,nymin1 nx-1, ny-1, respectively *
155  ! nxmax,nymax maximum dimension of wind fields in x and y *
156  ! direction, respectively *
157  ! nuvzmax,nwzmax maximum dimension of (u,v) and (w) wind fields in z
158  ! direction (for fields on eta levels) *
159  ! nlev_ec number of vertical levels ecmwf model *
160  ! u10,v10 wind fields at 10 m *
161  ! tt2,td2 temperature and dew point temperature at 2 m *
162  ! tth,qvh temperature and specific humidity on original *
163  ! eta levels *
164  ! ps surface pressure *
165  ! sd snow depth (but not used afterwards!) *
166  ! msl mean sea level pressure *
167  ! ttc total cloud cover *
168  ! lsprec large scale precipitation *
169  ! convprec convective precipitation *
170  ! sshf sensible heat flux *
171  ! ssr solar radiation *
172  ! surfstr surface stress *
173  ! oro orography *
174  ! excessoro standard deviation of orography *
175  ! lsm land sea mask *
176  ! hmix mixing height *
177  !***********************************************************************
178 
179 !-----------------------------------------------------------------------------
180 
181  hflswitch=.false.
182  strswitch=.false.
183  levdiff2=nlev_ec-nwz+1
184  iumax=0
185  iwmax=0
186 
187 
188  ! OPENING OF DATA FILE (GRIB CODE)
189 
190  !HSO
191 5 call grib_open_file(ifile,path(3)(1:length(3)) &
192  //trim(wfname(indj)),'r',iret)
193  if (iret.ne.grib_success) then
194  goto 888 ! ERROR DETECTED
195  endif
196  !turn on support for multi fields messages
197  call grib_multi_support_on
198 
199  numpt=0
200  numpu=0
201  numpv=0
202  numpw=0
203  numprh=0
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 50 ! EOF DETECTED
212  elseif (iret.ne.grib_success) then
213  goto 888 ! ERROR DETECTED
214  endif
215 
216  !first see if we read GRIB1 or GRIB2
217  call grib_get_int(igrib,'editionNumber',gribver,iret)
218  call grib_check(iret,gribfunction,griberrormsg)
219 
220  if (gribver.eq.1) then ! GRIB Edition 1
221 
222  !read the grib1 identifiers
223  call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
224  call grib_check(iret,gribfunction,griberrormsg)
225  call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
226  call grib_check(iret,gribfunction,griberrormsg)
227  call grib_get_int(igrib,'level',isec1(8),iret)
228  call grib_check(iret,gribfunction,griberrormsg)
229 
230  else ! GRIB Edition 2
231 
232  !read the grib2 identifiers
233  call grib_get_int(igrib,'discipline',discipl,iret)
234  call grib_check(iret,gribfunction,griberrormsg)
235  call grib_get_int(igrib,'parameterCategory',parcat,iret)
236  call grib_check(iret,gribfunction,griberrormsg)
237  call grib_get_int(igrib,'parameterNumber',parnum,iret)
238  call grib_check(iret,gribfunction,griberrormsg)
239  call grib_get_int(igrib,'typeOfFirstFixedSurface',typsurf,iret)
240  call grib_check(iret,gribfunction,griberrormsg)
241  call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
242  valsurf,iret)
243  call grib_check(iret,gribfunction,griberrormsg)
244 
245  !convert to grib1 identifiers
246  isec1(6)=-1
247  isec1(7)=-1
248  isec1(8)=-1
249  if ((parcat.eq.0).and.(parnum.eq.0).and.(typsurf.eq.100)) then ! T
250  isec1(6)=11 ! indicatorOfParameter
251  isec1(7)=100 ! indicatorOfTypeOfLevel
252  isec1(8)=valsurf/100 ! level, convert to hPa
253  elseif ((parcat.eq.2).and.(parnum.eq.2).and.(typsurf.eq.100)) then ! U
254  isec1(6)=33 ! indicatorOfParameter
255  isec1(7)=100 ! indicatorOfTypeOfLevel
256  isec1(8)=valsurf/100 ! level, convert to hPa
257  elseif ((parcat.eq.2).and.(parnum.eq.3).and.(typsurf.eq.100)) then ! V
258  isec1(6)=34 ! indicatorOfParameter
259  isec1(7)=100 ! indicatorOfTypeOfLevel
260  isec1(8)=valsurf/100 ! level, convert to hPa
261  elseif ((parcat.eq.2).and.(parnum.eq.8).and.(typsurf.eq.100)) then ! W
262  isec1(6)=39 ! indicatorOfParameter
263  isec1(7)=100 ! indicatorOfTypeOfLevel
264  isec1(8)=valsurf/100 ! level, convert to hPa
265  elseif ((parcat.eq.1).and.(parnum.eq.1).and.(typsurf.eq.100)) then ! RH
266  isec1(6)=52 ! indicatorOfParameter
267  isec1(7)=100 ! indicatorOfTypeOfLevel
268  isec1(8)=valsurf/100 ! level, convert to hPa
269  elseif ((parcat.eq.1).and.(parnum.eq.1).and.(typsurf.eq.103)) then ! RH2
270  isec1(6)=52 ! indicatorOfParameter
271  isec1(7)=105 ! indicatorOfTypeOfLevel
272  isec1(8)=2
273  elseif ((parcat.eq.0).and.(parnum.eq.0).and.(typsurf.eq.103)) then ! T2
274  isec1(6)=11 ! indicatorOfParameter
275  isec1(7)=105 ! indicatorOfTypeOfLevel
276  isec1(8)=2
277  elseif ((parcat.eq.2).and.(parnum.eq.2).and.(typsurf.eq.103)) then ! U10
278  isec1(6)=33 ! indicatorOfParameter
279  isec1(7)=105 ! indicatorOfTypeOfLevel
280  isec1(8)=10
281  elseif ((parcat.eq.2).and.(parnum.eq.3).and.(typsurf.eq.103)) then ! V10
282  isec1(6)=34 ! indicatorOfParameter
283  isec1(7)=105 ! indicatorOfTypeOfLevel
284  isec1(8)=10
285  elseif ((parcat.eq.3).and.(parnum.eq.1).and.(typsurf.eq.101)) then ! SLP
286  isec1(6)=2 ! indicatorOfParameter
287  isec1(7)=102 ! indicatorOfTypeOfLevel
288  isec1(8)=0
289  elseif ((parcat.eq.3).and.(parnum.eq.0).and.(typsurf.eq.1)) then ! SP
290  isec1(6)=1 ! indicatorOfParameter
291  isec1(7)=1 ! indicatorOfTypeOfLevel
292  isec1(8)=0
293  elseif ((parcat.eq.1).and.(parnum.eq.13).and.(typsurf.eq.1)) then ! SNOW
294  isec1(6)=66 ! indicatorOfParameter
295  isec1(7)=1 ! indicatorOfTypeOfLevel
296  isec1(8)=0
297  elseif ((parcat.eq.0).and.(parnum.eq.0).and.(typsurf.eq.104)) then ! T sigma 0
298  isec1(6)=11 ! indicatorOfParameter
299  isec1(7)=107 ! indicatorOfTypeOfLevel
300  isec1(8)=0.995 ! lowest sigma level
301  elseif ((parcat.eq.2).and.(parnum.eq.2).and.(typsurf.eq.104)) then ! U sigma 0
302  isec1(6)=33 ! indicatorOfParameter
303  isec1(7)=107 ! indicatorOfTypeOfLevel
304  isec1(8)=0.995 ! lowest sigma level
305  elseif ((parcat.eq.2).and.(parnum.eq.3).and.(typsurf.eq.104)) then ! V sigma 0
306  isec1(6)=34 ! indicatorOfParameter
307  isec1(7)=107 ! indicatorOfTypeOfLevel
308  isec1(8)=0.995 ! lowest sigma level
309  elseif ((parcat.eq.3).and.(parnum.eq.5).and.(typsurf.eq.1)) then ! TOPO
310  isec1(6)=7 ! indicatorOfParameter
311  isec1(7)=1 ! indicatorOfTypeOfLevel
312  isec1(8)=0
313  elseif ((parcat.eq.0).and.(parnum.eq.0).and.(typsurf.eq.1) &
314  .and.(discipl.eq.2)) then ! LSM
315  isec1(6)=81 ! indicatorOfParameter
316  isec1(7)=1 ! indicatorOfTypeOfLevel
317  isec1(8)=0
318  elseif ((parcat.eq.3).and.(parnum.eq.196).and.(typsurf.eq.1)) then ! BLH
319  isec1(6)=221 ! indicatorOfParameter
320  isec1(7)=1 ! indicatorOfTypeOfLevel
321  isec1(8)=0
322  elseif ((parcat.eq.1).and.(parnum.eq.7).and.(typsurf.eq.1)) then ! LSP/TP
323  isec1(6)=62 ! indicatorOfParameter
324  isec1(7)=1 ! indicatorOfTypeOfLevel
325  isec1(8)=0
326  elseif ((parcat.eq.1).and.(parnum.eq.196).and.(typsurf.eq.1)) then ! CP
327  isec1(6)=63 ! indicatorOfParameter
328  isec1(7)=1 ! indicatorOfTypeOfLevel
329  isec1(8)=0
330  endif
331 
332  endif ! gribVer
333 
334  if (isec1(6).ne.-1) then
335  ! get the size and data of the values array
336  call grib_get_real4_array(igrib,'values',zsec4,iret)
337  call grib_check(iret,gribfunction,griberrormsg)
338  endif
339 
340  if(ifield.eq.1) then
341 
342  !get the required fields from section 2
343  !store compatible to gribex input
344  call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
345  isec2(2),iret)
346  call grib_check(iret,gribfunction,griberrormsg)
347  call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
348  isec2(3),iret)
349  call grib_check(iret,gribfunction,griberrormsg)
350  call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
351  xauxin,iret)
352  call grib_check(iret,gribfunction,griberrormsg)
353  call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
354  yauxin,iret)
355  call grib_check(iret,gribfunction,griberrormsg)
356  xaux=xauxin+real(nxshift)*dx
357  yaux=yauxin
358 
359  ! CHECK GRID SPECIFICATIONS
360 
361  if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
362  if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
363  if(xaux.eq.0.) xaux=-179.0 ! NCEP DATA
364  xaux0=xlon0
365  yaux0=ylat0
366  if(xaux.lt.0.) xaux=xaux+360.
367  if(yaux.lt.0.) yaux=yaux+360.
368  if(xaux0.lt.0.) xaux0=xaux0+360.
369  if(yaux0.lt.0.) yaux0=yaux0+360.
370  if(abs(xaux-xaux0).gt.eps) &
371  stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
372  if(abs(yaux-yaux0).gt.eps) &
373  stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
374  endif
375  !HSO end of edits
376 
377  i179=nint(179./dx)
378  if (dx.lt.0.7) then
379  i180=nint(180./dx)+1 ! 0.5 deg data
380  else
381  i180=nint(179./dx)+1 ! 1 deg data
382  endif
383  i181=i180+1
384 
385  if (isec1(6).ne.-1) then
386 
387  do j=0,nymin1
388  do i=0,nxfield-1
389  if((isec1(6).eq.011).and.(isec1(7).eq.100)) then
390  ! TEMPERATURE
391  if((i.eq.0).and.(j.eq.0)) then
392  do ii=1,nuvz
393  if ((isec1(8)*100.0).eq.akz(ii)) numpt=ii
394  end do
395  endif
396  help=zsec4(nxfield*(ny-j-1)+i+1)
397  if(i.le.i180) then
398  tth(i179+i,j,numpt,n)=help
399  else
400  tth(i-i181,j,numpt,n)=help
401  endif
402  endif
403  if((isec1(6).eq.033).and.(isec1(7).eq.100)) then
404  ! U VELOCITY
405  if((i.eq.0).and.(j.eq.0)) then
406  do ii=1,nuvz
407  if ((isec1(8)*100.0).eq.akz(ii)) numpu=ii
408  end do
409  endif
410  help=zsec4(nxfield*(ny-j-1)+i+1)
411  if(i.le.i180) then
412  uuh(i179+i,j,numpu)=help
413  else
414  uuh(i-i181,j,numpu)=help
415  endif
416  endif
417  if((isec1(6).eq.034).and.(isec1(7).eq.100)) then
418  ! V VELOCITY
419  if((i.eq.0).and.(j.eq.0)) then
420  do ii=1,nuvz
421  if ((isec1(8)*100.0).eq.akz(ii)) numpv=ii
422  end do
423  endif
424  help=zsec4(nxfield*(ny-j-1)+i+1)
425  if(i.le.i180) then
426  vvh(i179+i,j,numpv)=help
427  else
428  vvh(i-i181,j,numpv)=help
429  endif
430  endif
431  if((isec1(6).eq.052).and.(isec1(7).eq.100)) then
432  ! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER
433  if((i.eq.0).and.(j.eq.0)) then
434  do ii=1,nuvz
435  if ((isec1(8)*100.0).eq.akz(ii)) numprh=ii
436  end do
437  endif
438  help=zsec4(nxfield*(ny-j-1)+i+1)
439  if(i.le.i180) then
440  qvh(i179+i,j,numprh,n)=help
441  else
442  qvh(i-i181,j,numprh,n)=help
443  endif
444  endif
445  if((isec1(6).eq.001).and.(isec1(7).eq.001)) then
446  ! SURFACE PRESSURE
447  help=zsec4(nxfield*(ny-j-1)+i+1)
448  if(i.le.i180) then
449  ps(i179+i,j,1,n)=help
450  else
451  ps(i-i181,j,1,n)=help
452  endif
453  endif
454  if((isec1(6).eq.039).and.(isec1(7).eq.100)) then
455  ! W VELOCITY
456  if((i.eq.0).and.(j.eq.0)) then
457  do ii=1,nuvz
458  if ((isec1(8)*100.0).eq.akz(ii)) numpw=ii
459  end do
460  endif
461  help=zsec4(nxfield*(ny-j-1)+i+1)
462  if(i.le.i180) then
463  wwh(i179+i,j,numpw)=help
464  else
465  wwh(i-i181,j,numpw)=help
466  endif
467  endif
468  if((isec1(6).eq.066).and.(isec1(7).eq.001)) then
469  ! SNOW DEPTH
470  help=zsec4(nxfield*(ny-j-1)+i+1)
471  if(i.le.i180) then
472  sd(i179+i,j,1,n)=help
473  else
474  sd(i-i181,j,1,n)=help
475  endif
476  endif
477  if((isec1(6).eq.002).and.(isec1(7).eq.102)) then
478  ! MEAN SEA LEVEL PRESSURE
479  help=zsec4(nxfield*(ny-j-1)+i+1)
480  if(i.le.i180) then
481  msl(i179+i,j,1,n)=help
482  else
483  msl(i-i181,j,1,n)=help
484  endif
485  endif
486  if((isec1(6).eq.071).and.(isec1(7).eq.244)) then
487  ! TOTAL CLOUD COVER
488  help=zsec4(nxfield*(ny-j-1)+i+1)
489  if(i.le.i180) then
490  tcc(i179+i,j,1,n)=help
491  else
492  tcc(i-i181,j,1,n)=help
493  endif
494  endif
495  if((isec1(6).eq.033).and.(isec1(7).eq.105).and. &
496  (isec1(8).eq.10)) then
497  ! 10 M U VELOCITY
498  help=zsec4(nxfield*(ny-j-1)+i+1)
499  if(i.le.i180) then
500  u10(i179+i,j,1,n)=help
501  else
502  u10(i-i181,j,1,n)=help
503  endif
504  endif
505  if((isec1(6).eq.034).and.(isec1(7).eq.105).and. &
506  (isec1(8).eq.10)) then
507  ! 10 M V VELOCITY
508  help=zsec4(nxfield*(ny-j-1)+i+1)
509  if(i.le.i180) then
510  v10(i179+i,j,1,n)=help
511  else
512  v10(i-i181,j,1,n)=help
513  endif
514  endif
515  if((isec1(6).eq.011).and.(isec1(7).eq.105).and. &
516  (isec1(8).eq.02)) then
517  ! 2 M TEMPERATURE
518  help=zsec4(nxfield*(ny-j-1)+i+1)
519  if(i.le.i180) then
520  tt2(i179+i,j,1,n)=help
521  else
522  tt2(i-i181,j,1,n)=help
523  endif
524  endif
525  if((isec1(6).eq.017).and.(isec1(7).eq.105).and. &
526  (isec1(8).eq.02)) then
527  ! 2 M DEW POINT TEMPERATURE
528  help=zsec4(nxfield*(ny-j-1)+i+1)
529  if(i.le.i180) then
530  td2(i179+i,j,1,n)=help
531  else
532  td2(i-i181,j,1,n)=help
533  endif
534  endif
535  if((isec1(6).eq.062).and.(isec1(7).eq.001)) then
536  ! LARGE SCALE PREC.
537  help=zsec4(nxfield*(ny-j-1)+i+1)
538  if(i.le.i180) then
539  lsprec(i179+i,j,1,n)=help
540  else
541  lsprec(i-i181,j,1,n)=help
542  endif
543  endif
544  if((isec1(6).eq.063).and.(isec1(7).eq.001)) then
545  ! CONVECTIVE PREC.
546  help=zsec4(nxfield*(ny-j-1)+i+1)
547  if(i.le.i180) then
548  convprec(i179+i,j,1,n)=help
549  else
550  convprec(i-i181,j,1,n)=help
551  endif
552  endif
553  if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
554  ! TOPOGRAPHY
555  help=zsec4(nxfield*(ny-j-1)+i+1)
556  if(i.le.i180) then
557  oro(i179+i,j)=help
558  excessoro(i179+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
559  else
560  oro(i-i181,j)=help
561  excessoro(i-i181,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
562  endif
563  endif
564  if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
565  ! LAND SEA MASK
566  help=zsec4(nxfield*(ny-j-1)+i+1)
567  if(i.le.i180) then
568  lsm(i179+i,j)=help
569  else
570  lsm(i-i181,j)=help
571  endif
572  endif
573  if((isec1(6).eq.221).and.(isec1(7).eq.001)) then
574  ! MIXING HEIGHT
575  help=zsec4(nxfield*(ny-j-1)+i+1)
576  if(i.le.i180) then
577  hmix(i179+i,j,1,n)=help
578  else
579  hmix(i-i181,j,1,n)=help
580  endif
581  endif
582  if((isec1(6).eq.052).and.(isec1(7).eq.105).and. &
583  (isec1(8).eq.02)) then
584  ! 2 M RELATIVE HUMIDITY
585  help=zsec4(nxfield*(ny-j-1)+i+1)
586  if(i.le.i180) then
587  qvh2(i179+i,j)=help
588  else
589  qvh2(i-i181,j)=help
590  endif
591  endif
592  if((isec1(6).eq.011).and.(isec1(7).eq.107)) then
593  ! TEMPERATURE LOWEST SIGMA LEVEL
594  help=zsec4(nxfield*(ny-j-1)+i+1)
595  if(i.le.i180) then
596  tlev1(i179+i,j)=help
597  else
598  tlev1(i-i181,j)=help
599  endif
600  endif
601  if((isec1(6).eq.033).and.(isec1(7).eq.107)) then
602  ! U VELOCITY LOWEST SIGMA LEVEL
603  help=zsec4(nxfield*(ny-j-1)+i+1)
604  if(i.le.i180) then
605  ulev1(i179+i,j)=help
606  else
607  ulev1(i-i181,j)=help
608  endif
609  endif
610  if((isec1(6).eq.034).and.(isec1(7).eq.107)) then
611  ! V VELOCITY LOWEST SIGMA LEVEL
612  help=zsec4(nxfield*(ny-j-1)+i+1)
613  if(i.le.i180) then
614  vlev1(i179+i,j)=help
615  else
616  vlev1(i-i181,j)=help
617  endif
618  endif
619 
620  end do
621  end do
622 
623  endif
624 
625  if((isec1(6).eq.33).and.(isec1(7).eq.100)) then
626  ! NCEP ISOBARIC LEVELS
627  iumax=iumax+1
628  endif
629 
630  call grib_release(igrib)
631  goto 10 !! READ NEXT LEVEL OR PARAMETER
632  !
633  ! CLOSING OF INPUT DATA FILE
634  !
635 
636  !HSO close grib file
637 50 continue
638  call grib_close_file(ifile)
639 
640  ! SENS. HEAT FLUX
641  sshf(:,:,1,n)=0.0 ! not available from gfs.tccz.pgrbfxx files
642  hflswitch=.false. ! Heat flux not available
643  ! SOLAR RADIATIVE FLUXES
644  ssr(:,:,1,n)=0.0 ! not available from gfs.tccz.pgrbfxx files
645  ! EW SURFACE STRESS
646  ewss=0.0 ! not available from gfs.tccz.pgrbfxx files
647  ! NS SURFACE STRESS
648  nsss=0.0 ! not available from gfs.tccz.pgrbfxx files
649  strswitch=.false. ! stress not available
650 
651  ! CONVERT TP TO LSP (GRIB2 only)
652  if (gribver.eq.2) then
653  do j=0,nymin1
654  do i=0,nxfield-1
655  if(i.le.i180) then
656  if (convprec(i179+i,j,1,n).lt.lsprec(i179+i,j,1,n)) then ! neg precip would occur
657  lsprec(i179+i,j,1,n)= &
658  lsprec(i179+i,j,1,n)-convprec(i179+i,j,1,n)
659  else
660  lsprec(i179+i,j,1,n)=0
661  endif
662  else
663  if (convprec(i-i181,j,1,n).lt.lsprec(i-i181,j,1,n)) then
664  lsprec(i-i181,j,1,n)= &
665  lsprec(i-i181,j,1,n)-convprec(i-i181,j,1,n)
666  else
667  lsprec(i-i181,j,1,n)=0
668  endif
669  endif
670  enddo
671  enddo
672  endif
673  !HSO end edits
674 
675 
676  ! TRANSFORM RH TO SPECIFIC HUMIDITY
677 
678  do j=0,ny-1
679  do i=0,nxfield-1
680  do k=1,nuvz
681  help=qvh(i,j,k,n)
682  temp=tth(i,j,k,n)
683  plev1=akm(k)+bkm(k)*ps(i,j,1,n)
684  elev=ew(temp)*help/100.0
685  qvh(i,j,k,n)=xmwml*(elev/(plev1-((1.0-xmwml)*elev)))
686  end do
687  end do
688  end do
689 
690  ! CALCULATE 2 M DEW POINT FROM 2 M RELATIVE HUMIDITY
691  ! USING BOLTON'S (1980) FORMULA
692  ! BECAUSE td2 IS NOT AVAILABLE FROM NCEP GFS DATA
693 
694  do j=0,ny-1
695  do i=0,nxfield-1
696  help=qvh2(i,j)
697  temp=tt2(i,j,1,n)
698  elev=ew(temp)/100.*help/100. !vapour pressure in hPa
699  td2(i,j,1,n)=243.5/(17.67/log(elev/6.112)-1)+273.
700  if (help.le.0.) td2(i,j,1,n)=tt2(i,j,1,n)
701  end do
702  end do
703 
704  if(levdiff2.eq.0) then
705  iwmax=nlev_ec+1
706  do i=0,nxmin1
707  do j=0,nymin1
708  wwh(i,j,nlev_ec+1)=0.
709  end do
710  end do
711  endif
712 
713 
714  ! For global fields, assign the leftmost data column also to the rightmost
715  ! data column; if required, shift whole grid by nxshift grid points
716  !*************************************************************************
717 
718  if (xglobal) then
719  call shift_field_0(ewss,nxfield,ny)
720  call shift_field_0(nsss,nxfield,ny)
721  call shift_field_0(oro,nxfield,ny)
722  call shift_field_0(excessoro,nxfield,ny)
723  call shift_field_0(lsm,nxfield,ny)
724  call shift_field_0(ulev1,nxfield,ny)
725  call shift_field_0(vlev1,nxfield,ny)
726  call shift_field_0(tlev1,nxfield,ny)
727  call shift_field_0(qvh2,nxfield,ny)
728  call shift_field(ps,nxfield,ny,1,1,2,n)
729  call shift_field(sd,nxfield,ny,1,1,2,n)
730  call shift_field(msl,nxfield,ny,1,1,2,n)
731  call shift_field(tcc,nxfield,ny,1,1,2,n)
732  call shift_field(u10,nxfield,ny,1,1,2,n)
733  call shift_field(v10,nxfield,ny,1,1,2,n)
734  call shift_field(tt2,nxfield,ny,1,1,2,n)
735  call shift_field(td2,nxfield,ny,1,1,2,n)
736  call shift_field(lsprec,nxfield,ny,1,1,2,n)
737  call shift_field(convprec,nxfield,ny,1,1,2,n)
738  call shift_field(sshf,nxfield,ny,1,1,2,n)
739  call shift_field(ssr,nxfield,ny,1,1,2,n)
740  call shift_field(hmix,nxfield,ny,1,1,2,n)
741  call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n)
742  call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n)
743  call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
744  call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
745  call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
746  endif
747 
748  do i=0,nxmin1
749  do j=0,nymin1
750  ! Convert precip. from mm/s -> mm/hour
751  convprec(i,j,1,n)=convprec(i,j,1,n)*3600.
752  lsprec(i,j,1,n)=lsprec(i,j,1,n)*3600.
753  surfstr(i,j,1,n)=sqrt(ewss(i,j)**2+nsss(i,j)**2)
754  end do
755  end do
756 
757  if ((.not.hflswitch).or.(.not.strswitch)) then
758  ! write(*,*) 'WARNING: No flux data contained in GRIB file ',
759  ! + wfname(indj)
760 
761  ! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
762  !***************************************************************************
763 
764  do i=0,nxmin1
765  do j=0,nymin1
766  hlev1=30.0 ! HEIGHT OF FIRST MODEL SIGMA LAYER
767  ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2)
768  fflev1=sqrt(ulev1(i,j)**2+vlev1(i,j)**2)
769  call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
770  tt2(i,j,1,n),tlev1(i,j),ff10m,fflev1, &
771  surfstr(i,j,1,n),sshf(i,j,1,n))
772  if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
773  if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
774  end do
775  end do
776  endif
777 
778 
779  if(iumax.ne.nuvz) stop 'READWIND: NUVZ NOT CONSISTENT'
780  if(iumax.ne.nwz) stop 'READWIND: NWZ NOT CONSISTENT'
781 
782  return
783 888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD #### '
784  write(*,*) ' #### ',wfname(indj),' #### '
785  write(*,*) ' #### IS NOT GRIB FORMAT !!! #### '
786  stop 'Execution terminated'
787 999 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD #### '
788  write(*,*) ' #### ',wfname(indj),' #### '
789  write(*,*) ' #### CANNOT BE OPENED !!! #### '
790  stop 'Execution terminated'
791 
792 end subroutine readwind_gfs
ew
real function ew(x)
Definition: ew.f90:22
shift_field
subroutine shift_field(field, nxf, nyf, nzfmax, nzf, nmax, n)
Definition: shift_field.f90:22
par_mod
Definition: par_mod.F90:34
readwind_gfs
subroutine readwind_gfs(indj, n, uuh, vvh, wwh)
Definition: readwind_gfs.f90:22
pbl_profile
subroutine pbl_profile(ps, td2m, zml1, t2m, tml1, u10m, uml1, stress, hf)
Definition: pbl_profile.f90:22
shift_field_0
subroutine shift_field_0(field, nxf, nyf)
Definition: shift_field_0.f90:22
com_mod
Definition: com_mod.f90:13
Generated on Thu Oct 15 2015 06:20:58 for CTBTO FLEXPART WO4 (2015-10-15) by   doxygen 1.8.6