source: flexpart.git/src/readwind_gfs.f90 @ 467460a

GFS_025dev
Last change on this file since 467460a was 467460a, checked in by Espen Sollum ATMOS <eso@…>, 3 years ago

First commit of files from Hachinger

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