readwind_gfs.f90 25.1 KB
Newer Older
1
2
! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
! SPDX-License-Identifier: GPL-3.0-or-later
Matthias Langer's avatar
 
Matthias Langer committed
3

4
subroutine readwind_gfs(indj,n,uuh,vvh,wwh)
Matthias Langer's avatar
 
Matthias Langer committed
5

Espen Sollum's avatar
Espen Sollum committed
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
!***********************************************************************
!*                                                                     *
!*             TRAJECTORY MODEL SUBROUTINE READWIND                    *
!*                                                                     *
!***********************************************************************
!*                                                                     *
!*             AUTHOR:      G. WOTAWA                                  *
!*             DATE:        1997-08-05                                 *
!*             LAST UPDATE: 2000-10-17, Andreas Stohl                  *
!*             CHANGE: 01/02/2001, Bernd C. Krueger, Variables tth and *
!*                     qvh (on eta coordinates) in common block        *
!*             CHANGE: 16/11/2005, Caroline Forster, GFS data          *
!*             CHANGE: 11/01/2008, Harald Sodemann, Input of GRIB1/2   *
!*                     data with the ECMWF grib_api library            *
!*             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
!*                                 ECMWF grib_api                      *
!                                                                      *
!   Unified ECMWF and GFS builds                                       *
!   Marian Harustak, 12.5.2017                                         *
!     - Renamed routine from readwind to readwind_gfs                  *
!*                                                                     *
!***********************************************************************
!*                                                                     *
!* DESCRIPTION:                                                        *
!*                                                                     *
!* READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE   *
!* INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE          *
!*                                                                     *
!* INPUT:                                                              *
!* indj               indicates number of the wind field to be read in *
!* n                  temporal index for meteorological fields (1 to 3)*
!*                                                                     *
!* IMPORTANT VARIABLES FROM COMMON BLOCK:                              *
!*                                                                     *
!* wfname             File name of data to be read in                  *
!* nx,ny,nuvz,nwz     expected field dimensions                        *
!* nlev_ec            number of vertical levels ecmwf model            *
!* uu,vv,ww           wind fields                                      *
!* tt,qv              temperature and specific humidity                *
!* ps                 surface pressure                                 *
!*                                                                     *
!***********************************************************************
Matthias Langer's avatar
 
Matthias Langer committed
48

49
  use eccodes
Matthias Langer's avatar
 
Matthias Langer committed
50
51
52
53
54
  use par_mod
  use com_mod

  implicit none

Espen Sollum's avatar
Espen Sollum committed
55
!HSO  new parameters for grib_api
Matthias Langer's avatar
 
Matthias Langer committed
56
57
58
  integer :: ifile
  integer :: iret
  integer :: igrib
Espen Sollum's avatar
Espen Sollum committed
59
60
  integer :: gribVer,parCat,parNum,typSurf,discipl,valSurf
!HSO end edits
Matthias Langer's avatar
 
Matthias Langer committed
61
62
63
64
65
  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
  integer :: ii,indj,i,j,k,n,levdiff2,ifield,iumax,iwmax

Espen Sollum's avatar
Espen Sollum committed
66
! NCEP
67
  integer :: numpt,numpu,numpv,numpw,numprh,numpclwch
Matthias Langer's avatar
 
Matthias Langer committed
68
69
70
71
72
73
  real :: help, temp, ew
  real :: elev
  real :: ulev1(0:nxmax-1,0:nymax-1),vlev1(0:nxmax-1,0:nymax-1)
  real :: tlev1(0:nxmax-1,0:nymax-1)
  real :: qvh2(0:nxmax-1,0:nymax-1)

74
  integer :: i180
Matthias Langer's avatar
 
Matthias Langer committed
75

Espen Sollum's avatar
Espen Sollum committed
76
77
! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
!HSO kept isec1, isec2 and zsec4 for consistency with gribex GRIB input
Matthias Langer's avatar
 
Matthias Langer committed
78
79

  integer :: isec1(8),isec2(3)
80
  real           :: xsec18
Matthias Langer's avatar
 
Matthias Langer committed
81
82
  real(kind=4) :: zsec4(jpunp)
  real(kind=4) :: xaux,yaux,xaux0,yaux0
83
  real,parameter :: eps=spacing(2.0_4*360.0_4)
Matthias Langer's avatar
 
Matthias Langer committed
84
85
86
87
88
89
  real(kind=8) :: xauxin,yauxin
  real(kind=4) :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
  real :: plev1,hlev1,ff10m,fflev1

  logical :: hflswitch,strswitch

Espen Sollum's avatar
Espen Sollum committed
90
!HSO  for grib api error messages
Matthias Langer's avatar
 
Matthias Langer committed
91
92
  character(len=24) :: gribErrorMsg = 'Error reading grib file'
  character(len=20) :: gribFunction = 'readwind_gfs'
Sabine's avatar
Sabine committed
93
  character(len=20) :: shortname
Matthias Langer's avatar
 
Matthias Langer committed
94
95
96
97
98
99
100
101
102


  hflswitch=.false.
  strswitch=.false.
  levdiff2=nlev_ec-nwz+1
  iumax=0
  iwmax=0


Espen Sollum's avatar
Espen Sollum committed
103
! OPENING OF DATA FILE (GRIB CODE)
Matthias Langer's avatar
 
Matthias Langer committed
104

Espen Sollum's avatar
Espen Sollum committed
105
!HSO
Sabine's avatar
Sabine committed
106
  call grib_open_file(ifile,path(3)(1:length(3)) &
Espen Sollum's avatar
Espen Sollum committed
107
       //trim(wfname(indj)),'r',iret)
Matthias Langer's avatar
 
Matthias Langer committed
108
109
110
  if (iret.ne.GRIB_SUCCESS) then
    goto 888   ! ERROR DETECTED
  endif
Espen Sollum's avatar
Espen Sollum committed
111
!turn on support for multi fields messages
Matthias Langer's avatar
 
Matthias Langer committed
112
113
114
115
116
117
118
  call grib_multi_support_on

  numpt=0
  numpu=0
  numpv=0
  numpw=0
  numprh=0
119
  numpclwch=0
Matthias Langer's avatar
 
Matthias Langer committed
120
  ifield=0
Espen Sollum's avatar
Espen Sollum committed
121
122
123
124
10 ifield=ifield+1
!
! GET NEXT FIELDS
!
Matthias Langer's avatar
 
Matthias Langer committed
125
126
127
128
129
130
131
  call grib_new_from_file(ifile,igrib,iret)
  if (iret.eq.GRIB_END_OF_FILE)  then
    goto 50    ! EOF DETECTED
  elseif (iret.ne.GRIB_SUCCESS) then
    goto 888   ! ERROR DETECTED
  endif

Espen Sollum's avatar
Espen Sollum committed
132
!first see if we read GRIB1 or GRIB2
Matthias Langer's avatar
 
Matthias Langer committed
133
  call grib_get_int(igrib,'editionNumber',gribVer,iret)
134
!  call grib_check(iret,gribFunction,gribErrorMsg)
Matthias Langer's avatar
 
Matthias Langer committed
135
136
137

  if (gribVer.eq.1) then ! GRIB Edition 1

Espen Sollum's avatar
Espen Sollum committed
138
139
140
141
142
143
144
145
146
147
148
!read the grib1 identifiers

    call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
    call grib_check(iret,gribFunction,gribErrorMsg)
    call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
    call grib_check(iret,gribFunction,gribErrorMsg)
    call grib_get_int(igrib,'level',isec1(8),iret)
    call grib_check(iret,gribFunction,gribErrorMsg)
!JMA / SH: isec1(8) not evaluated any more below
!b/c with GRIB 2 this may be a real variable
    xsec18 = real(isec1(8))
Matthias Langer's avatar
 
Matthias Langer committed
149
150
151

  else ! GRIB Edition 2

Espen Sollum's avatar
Espen Sollum committed
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
!read the grib2 identifiers
    
    call grib_get_string(igrib,'shortName',shortname,iret)

    call grib_get_int(igrib,'discipline',discipl,iret)
    call grib_check(iret,gribFunction,gribErrorMsg)
    call grib_get_int(igrib,'parameterCategory',parCat,iret)
    call grib_check(iret,gribFunction,gribErrorMsg)
    call grib_get_int(igrib,'parameterNumber',parNum,iret)
    call grib_check(iret,gribFunction,gribErrorMsg)
    call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
    call grib_check(iret,gribFunction,gribErrorMsg)
!
    call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
         valSurf,iret)
    call grib_check(iret,gribFunction,gribErrorMsg)

!    write(*,*) 'Field: ',ifield,parCat,parNum,typSurf,shortname
!convert to grib1 identifiers
    isec1(6:8)=-1
    xsec18  =-1.0
!JMA / SH: isec1(8) not evaluated any more below
!b/c with GRIB 2 this may be a real variable
    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T
      isec1(6)=11          ! indicatorOfParameter
      isec1(7)=100         ! indicatorOfTypeOfLevel
      xsec18=valSurf/100.0 ! level, convert to hPa
    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
      isec1(6)=33          ! indicatorOfParameter
      isec1(7)=100         ! indicatorOfTypeOfLevel
      xsec18=valSurf/100.0 ! level, convert to hPa
    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.100)) then ! V
      isec1(6)=34          ! indicatorOfParameter
      isec1(7)=100         ! indicatorOfTypeOfLevel
      xsec18=valSurf/100.0 ! level, convert to hPa
    elseif ((parCat.eq.2).and.(parNum.eq.8).and.(typSurf.eq.100)) then ! W
      isec1(6)=39          ! indicatorOfParameter
      isec1(7)=100         ! indicatorOfTypeOfLevel
      xsec18=valSurf/100.0 ! level, convert to hPa
    elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.100)) then ! RH
      isec1(6)=52          ! indicatorOfParameter
      isec1(7)=100         ! indicatorOfTypeOfLevel
      xsec18=valSurf/100.0 ! level, convert to hPa
    elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.103)) then ! RH2
      isec1(6)=52          ! indicatorOfParameter
      isec1(7)=105         ! indicatorOfTypeOfLevel
      xsec18=real(2)
    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! T2
      isec1(6)=11          ! indicatorOfParameter
      isec1(7)=105         ! indicatorOfTypeOfLevel
      xsec18=real(2)
    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! U10
      isec1(6)=33          ! indicatorOfParameter
      isec1(7)=105         ! indicatorOfTypeOfLevel
      xsec18=real(10)
    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! V10
      isec1(6)=34          ! indicatorOfParameter
      isec1(7)=105         ! indicatorOfTypeOfLevel
      xsec18=real(10)
    elseif ((parCat.eq.1).and.(parNum.eq.22).and.(typSurf.eq.100)) then ! CLWMR Cloud Mixing Ratio [kg/kg]:
      isec1(6)=153         ! indicatorOfParameter
      isec1(7)=100         ! indicatorOfTypeOfLevel
      xsec18=valSurf/100.0 ! level, convert to hPa
    elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP
      isec1(6)=2           ! indicatorOfParameter
      isec1(7)=102         ! indicatorOfTypeOfLevel
      xsec18=real(0)
Espen Sollum's avatar
Espen Sollum committed
219
220
   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1).and. &
        (discipl.eq.0)) then ! SP
Espen Sollum's avatar
Espen Sollum committed
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
      isec1(6)=1           ! indicatorOfParameter
      isec1(7)=1           ! indicatorOfTypeOfLevel
      xsec18=real(0)
    elseif ((parCat.eq.1).and.(parNum.eq.13).and.(typSurf.eq.1)) then ! SNOW
      isec1(6)=66          ! indicatorOfParameter
      isec1(7)=1           ! indicatorOfTypeOfLevel
      xsec18=real(0)
    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.104)) then ! T sigma 0
      isec1(6)=11          ! indicatorOfParameter
      isec1(7)=107         ! indicatorOfTypeOfLevel
      xsec18=0.995         ! lowest sigma level
    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.104)) then ! U sigma 0
      isec1(6)=33          ! indicatorOfParameter
      isec1(7)=107         ! indicatorOfTypeOfLevel
      xsec18=0.995         ! lowest sigma level
    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.104)) then ! V sigma 0
      isec1(6)=34          ! indicatorOfParameter
      isec1(7)=107         ! indicatorOfTypeOfLevel
      xsec18=0.995         ! lowest sigma level
    elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
      isec1(6)=7           ! indicatorOfParameter
      isec1(7)=1           ! indicatorOfTypeOfLevel
      xsec18=real(0)
    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
         .and.(discipl.eq.2)) then ! LSM
      isec1(6)=81          ! indicatorOfParameter
      isec1(7)=1           ! indicatorOfTypeOfLevel
      xsec18=real(0)
    elseif ((parCat.eq.3).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! BLH
      isec1(6)=221         ! indicatorOfParameter
      isec1(7)=1           ! indicatorOfTypeOfLevel
      xsec18=real(0)
    elseif ((parCat.eq.1).and.(parNum.eq.7).and.(typSurf.eq.1)) then ! LSP/TP
      isec1(6)=62          ! indicatorOfParameter
      isec1(7)=1           ! indicatorOfTypeOfLevel
      xsec18=real(0)
    elseif ((parCat.eq.1).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! CP
      isec1(6)=63          ! indicatorOfParameter
      isec1(7)=1           ! indicatorOfTypeOfLevel
      xsec18=real(0)
    endif
Matthias Langer's avatar
 
Matthias Langer committed
262
263
264
265

  endif ! gribVer

  if (isec1(6).ne.-1) then
Espen Sollum's avatar
Espen Sollum committed
266
!  get the size and data of the values array
Matthias Langer's avatar
 
Matthias Langer committed
267
    call grib_get_real4_array(igrib,'values',zsec4,iret)
268
    call grib_check(iret,gribFunction,gribErrorMsg)
Matthias Langer's avatar
 
Matthias Langer committed
269
270
271
272
  endif

  if(ifield.eq.1) then

Espen Sollum's avatar
Espen Sollum committed
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
!get the required fields from section 2
!store compatible to gribex input
    call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
         isec2(2),iret)
    call grib_check(iret,gribFunction,gribErrorMsg)
    call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
         isec2(3),iret)
    call grib_check(iret,gribFunction,gribErrorMsg)
    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
         xauxin,iret)
    call grib_check(iret,gribFunction,gribErrorMsg)
    call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
         yauxin,iret)
    call grib_check(iret,gribFunction,gribErrorMsg)
    xaux=xauxin+real(nxshift)*dx
    yaux=yauxin

! CHECK GRID SPECIFICATIONS
Matthias Langer's avatar
 
Matthias Langer committed
291
292
293

    if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
    if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
294
    if(xaux.eq.0.) xaux=-180.0     ! NCEP DATA
Matthias Langer's avatar
 
Matthias Langer committed
295
296
297
298
299
300
    xaux0=xlon0
    yaux0=ylat0
    if(xaux.lt.0.) xaux=xaux+360.
    if(yaux.lt.0.) yaux=yaux+360.
    if(xaux0.lt.0.) xaux0=xaux0+360.
    if(yaux0.lt.0.) yaux0=yaux0+360.
301
    if(abs(xaux-xaux0).gt.eps) then
Espen Sollum's avatar
Espen Sollum committed
302
303
      write (*, *) xaux, xaux0
      stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
304
305
    endif
    if(abs(yaux-yaux0).gt.eps) then
Espen Sollum's avatar
Espen Sollum committed
306
307
      write (*, *) yaux, yaux0
      stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
308
    end if
Matthias Langer's avatar
 
Matthias Langer committed
309
  endif
Espen Sollum's avatar
Espen Sollum committed
310
!HSO end of edits
Matthias Langer's avatar
 
Matthias Langer committed
311

312
  i180=nint(180./dx)
Matthias Langer's avatar
 
Matthias Langer committed
313
314
315

  if (isec1(6).ne.-1) then

316
317
! write (*, *) 'nxfield: ', nxfield, i180

Espen Sollum's avatar
Espen Sollum committed
318
319
320
    do j=0,nymin1
      do i=0,nxfield-1
        if((isec1(6).eq.011).and.(isec1(7).eq.100)) then
321
! TEMPERATURE
Espen Sollum's avatar
Espen Sollum committed
322
323
324
325
326
327
328
329
330
331
332
333
334
          if((i.eq.0).and.(j.eq.0)) then
            numpt=minloc(abs(xsec18*100.0-akz),dim=1)
          endif
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if (help.eq.0) then
            write (*, *) 'i, j: ', i, j
            stop 'help == 0.0'
          endif
          if(i.lt.i180) then
            tth(i180+i,j,numpt,n)=help
          else
            tth(i-i180,j,numpt,n)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
335
        endif
Espen Sollum's avatar
Espen Sollum committed
336
337
338
339
340
341
342
343
344
345
346
        if((isec1(6).eq.033).and.(isec1(7).eq.100)) then
! U VELOCITY
          if((i.eq.0).and.(j.eq.0)) then
            numpu=minloc(abs(xsec18*100.0-akz),dim=1)
          endif
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            uuh(i180+i,j,numpu)=help
          else
            uuh(i-i180,j,numpu)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
347
        endif
Espen Sollum's avatar
Espen Sollum committed
348
349
350
351
352
353
354
355
356
357
358
        if((isec1(6).eq.034).and.(isec1(7).eq.100)) then
! V VELOCITY
          if((i.eq.0).and.(j.eq.0)) then
            numpv=minloc(abs(xsec18*100.0-akz),dim=1)
          endif
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            vvh(i180+i,j,numpv)=help
          else
            vvh(i-i180,j,numpv)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
359
        endif
Espen Sollum's avatar
Espen Sollum committed
360
361
362
363
364
365
366
367
368
369
370
        if((isec1(6).eq.052).and.(isec1(7).eq.100)) then
! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER
          if((i.eq.0).and.(j.eq.0)) then
            numprh=minloc(abs(xsec18*100.0-akz),dim=1)
          endif
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            qvh(i180+i,j,numprh,n)=help
          else
            qvh(i-i180,j,numprh,n)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
371
        endif
Espen Sollum's avatar
Espen Sollum committed
372
373
374
375
376
377
378
379
        if((isec1(6).eq.001).and.(isec1(7).eq.001)) then
! SURFACE PRESSURE
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            ps(i180+i,j,1,n)=help
          else
            ps(i-i180,j,1,n)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
380
        endif
Espen Sollum's avatar
Espen Sollum committed
381
382
383
384
385
386
387
388
389
        if((isec1(6).eq.039).and.(isec1(7).eq.100)) then
! W VELOCITY
          if((i.eq.0).and.(j.eq.0)) numpw=minloc(abs(xsec18*100.0-akz),dim=1)
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            wwh(i180+i,j,numpw)=help
          else
            wwh(i-i180,j,numpw)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
390
        endif
Espen Sollum's avatar
Espen Sollum committed
391
392
393
394
395
396
397
398
        if((isec1(6).eq.066).and.(isec1(7).eq.001)) then
! SNOW DEPTH
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            sd(i180+i,j,1,n)=help
          else
            sd(i-i180,j,1,n)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
399
        endif
Espen Sollum's avatar
Espen Sollum committed
400
401
402
403
404
405
406
407
        if((isec1(6).eq.002).and.(isec1(7).eq.102)) then
! MEAN SEA LEVEL PRESSURE
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            msl(i180+i,j,1,n)=help
          else
            msl(i-i180,j,1,n)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
408
        endif
Espen Sollum's avatar
Espen Sollum committed
409
410
411
412
413
414
415
416
        if((isec1(6).eq.071).and.(isec1(7).eq.244)) then
! TOTAL CLOUD COVER
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            tcc(i180+i,j,1,n)=help
          else
            tcc(i-i180,j,1,n)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
417
        endif
Espen Sollum's avatar
Espen Sollum committed
418
419
420
421
422
423
424
425
426
        if((isec1(6).eq.033).and.(isec1(7).eq.105).and. &
             (nint(xsec18).eq.10)) then
! 10 M U VELOCITY
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            u10(i180+i,j,1,n)=help
          else
            u10(i-i180,j,1,n)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
427
        endif
Espen Sollum's avatar
Espen Sollum committed
428
429
430
431
432
433
434
435
436
        if((isec1(6).eq.034).and.(isec1(7).eq.105).and. &
             (nint(xsec18).eq.10)) then
! 10 M V VELOCITY
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            v10(i180+i,j,1,n)=help
          else
            v10(i-i180,j,1,n)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
437
        endif
Espen Sollum's avatar
Espen Sollum committed
438
439
440
441
442
443
444
445
446
        if((isec1(6).eq.011).and.(isec1(7).eq.105).and. &
             (nint(xsec18).eq.2)) then
! 2 M TEMPERATURE
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            tt2(i180+i,j,1,n)=help
          else
            tt2(i-i180,j,1,n)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
447
        endif
Espen Sollum's avatar
Espen Sollum committed
448
449
450
451
452
453
454
455
456
        if((isec1(6).eq.017).and.(isec1(7).eq.105).and. &
             (nint(xsec18).eq.2)) then
! 2 M DEW POINT TEMPERATURE
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            td2(i180+i,j,1,n)=help
          else
            td2(i-i180,j,1,n)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
457
        endif
Espen Sollum's avatar
Espen Sollum committed
458
459
460
461
462
463
464
465
        if((isec1(6).eq.062).and.(isec1(7).eq.001)) then
! LARGE SCALE PREC.
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            lsprec(i180+i,j,1,n)=help
          else
            lsprec(i-i180,j,1,n)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
466
        endif
Espen Sollum's avatar
Espen Sollum committed
467
468
469
470
471
472
473
474
        if((isec1(6).eq.063).and.(isec1(7).eq.001)) then
! CONVECTIVE PREC.
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            convprec(i180+i,j,1,n)=help
          else
            convprec(i-i180,j,1,n)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
475
        endif
Espen Sollum's avatar
Espen Sollum committed
476
477
478
479
480
481
482
483
484
485
        if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
! TOPOGRAPHY
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            oro(i180+i,j)=help
            excessoro(i180+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
          else
            oro(i-i180,j)=help
            excessoro(i-i180,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
          endif
Matthias Langer's avatar
 
Matthias Langer committed
486
        endif
Espen Sollum's avatar
Espen Sollum committed
487
488
489
490
491
492
493
494
        if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
! LAND SEA MASK
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            lsm(i180+i,j)=help
          else
            lsm(i-i180,j)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
495
        endif
Espen Sollum's avatar
Espen Sollum committed
496
497
498
499
500
501
502
503
        if((isec1(6).eq.221).and.(isec1(7).eq.001)) then
! MIXING HEIGHT
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            hmix(i180+i,j,1,n)=help
          else
            hmix(i-i180,j,1,n)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
504
        endif
Espen Sollum's avatar
Espen Sollum committed
505
506
507
508
509
510
511
512
513
        if((isec1(6).eq.052).and.(isec1(7).eq.105).and. &
             (nint(xsec18).eq.02)) then
! 2 M RELATIVE HUMIDITY
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            qvh2(i180+i,j)=help
          else
            qvh2(i-i180,j)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
514
        endif
Espen Sollum's avatar
Espen Sollum committed
515
516
517
518
519
520
521
522
        if((isec1(6).eq.011).and.(isec1(7).eq.107)) then
! TEMPERATURE LOWEST SIGMA LEVEL
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            tlev1(i180+i,j)=help
          else
            tlev1(i-i180,j)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
523
        endif
Espen Sollum's avatar
Espen Sollum committed
524
525
526
527
528
529
530
531
        if((isec1(6).eq.033).and.(isec1(7).eq.107)) then
! U VELOCITY LOWEST SIGMA LEVEL
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            ulev1(i180+i,j)=help
          else
            ulev1(i-i180,j)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
532
        endif
Espen Sollum's avatar
Espen Sollum committed
533
534
535
536
537
538
539
540
        if((isec1(6).eq.034).and.(isec1(7).eq.107)) then
! V VELOCITY LOWEST SIGMA LEVEL
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            vlev1(i180+i,j)=help
          else
            vlev1(i-i180,j)=help
          endif
Matthias Langer's avatar
 
Matthias Langer committed
541
        endif
Sabine's avatar
Sabine committed
542
! SEC & IP 12/2018 read GFS clouds
Espen Sollum's avatar
Espen Sollum committed
543
544
545
546
547
548
549
550
551
552
553
554
555
        if((isec1(6).eq.153).and.(isec1(7).eq.100)) then  !! CLWCR  Cloud liquid water content [kg/kg] 
          if((i.eq.0).and.(j.eq.0)) then
            numpclwch=minloc(abs(xsec18*100.0-akz),dim=1)
          endif
          help=zsec4(nxfield*(ny-j-1)+i+1)
          if(i.lt.i180) then
            clwch(i180+i,j,numpclwch,n)=help
          else
            clwch(i-i180,j,numpclwch,n)=help
          endif
          readclouds=.true.
          sumclouds=.true.
        endif
Sabine's avatar
Sabine committed
556

Matthias Langer's avatar
 
Matthias Langer committed
557

Espen Sollum's avatar
Espen Sollum committed
558
      end do
Matthias Langer's avatar
 
Matthias Langer committed
559
560
561
562
563
    end do

  endif

  if((isec1(6).eq.33).and.(isec1(7).eq.100)) then
Espen Sollum's avatar
Espen Sollum committed
564
! NCEP ISOBARIC LEVELS
Matthias Langer's avatar
 
Matthias Langer committed
565
566
567
568
569
    iumax=iumax+1
  endif

  call grib_release(igrib)
  goto 10                      !! READ NEXT LEVEL OR PARAMETER
Espen Sollum's avatar
Espen Sollum committed
570
571
572
!
! CLOSING OF INPUT DATA FILE
!
Matthias Langer's avatar
 
Matthias Langer committed
573

Espen Sollum's avatar
Espen Sollum committed
574
575
!HSO close grib file
50 continue
Matthias Langer's avatar
 
Matthias Langer committed
576
577
  call grib_close_file(ifile)

Espen Sollum's avatar
Espen Sollum committed
578
! SENS. HEAT FLUX
Matthias Langer's avatar
 
Matthias Langer committed
579
580
  sshf(:,:,1,n)=0.0     ! not available from gfs.tccz.pgrbfxx files
  hflswitch=.false.    ! Heat flux not available
Espen Sollum's avatar
Espen Sollum committed
581
! SOLAR RADIATIVE FLUXES
Matthias Langer's avatar
 
Matthias Langer committed
582
  ssr(:,:,1,n)=0.0      ! not available from gfs.tccz.pgrbfxx files
Espen Sollum's avatar
Espen Sollum committed
583
! EW SURFACE STRESS
Matthias Langer's avatar
 
Matthias Langer committed
584
  ewss=0.0         ! not available from gfs.tccz.pgrbfxx files
Espen Sollum's avatar
Espen Sollum committed
585
! NS SURFACE STRESS
Matthias Langer's avatar
 
Matthias Langer committed
586
587
588
  nsss=0.0         ! not available from gfs.tccz.pgrbfxx files
  strswitch=.false.    ! stress not available

Espen Sollum's avatar
Espen Sollum committed
589
! CONVERT TP TO LSP (GRIB2 only)
Matthias Langer's avatar
 
Matthias Langer committed
590
  if (gribVer.eq.2) then
591
    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) )
Matthias Langer's avatar
 
Matthias Langer committed
592
  endif
Espen Sollum's avatar
Espen Sollum committed
593
!HSO end edits
Matthias Langer's avatar
 
Matthias Langer committed
594
595


Espen Sollum's avatar
Espen Sollum committed
596
! TRANSFORM RH TO SPECIFIC HUMIDITY
Matthias Langer's avatar
 
Matthias Langer committed
597
598
599
600
601
602

  do j=0,ny-1
    do i=0,nxfield-1
      do k=1,nuvz
        help=qvh(i,j,k,n)
        temp=tth(i,j,k,n)
603
        if (temp .eq. 0.0) then 
Espen Sollum's avatar
Espen Sollum committed
604
          write (*, *) i, j, k, n
Espen Sollum's avatar
Espen Sollum committed
605
606
!          temp = 273.0
          stop
607
        endif
Matthias Langer's avatar
 
Matthias Langer committed
608
609
610
611
612
613
614
        plev1=akm(k)+bkm(k)*ps(i,j,1,n)
        elev=ew(temp)*help/100.0
        qvh(i,j,k,n)=xmwml*(elev/(plev1-((1.0-xmwml)*elev)))
      end do
    end do
  end do

Espen Sollum's avatar
Espen Sollum committed
615
616
617
! CALCULATE 2 M DEW POINT FROM 2 M RELATIVE HUMIDITY
! USING BOLTON'S (1980) FORMULA
! BECAUSE td2 IS NOT AVAILABLE FROM NCEP GFS DATA
Matthias Langer's avatar
 
Matthias Langer committed
618
619
620

  do j=0,ny-1
    do i=0,nxfield-1
Espen Sollum's avatar
Espen Sollum committed
621
622
623
624
625
626
627
628
629
      help=qvh2(i,j)
      temp=tt2(i,j,1,n)
      if (temp .eq. 0.0) then 
        write (*, *) i, j, n
        temp = 273.0
      endif
      elev=ew(temp)/100.*help/100.   !vapour pressure in hPa
      td2(i,j,1,n)=243.5/(17.67/log(elev/6.112)-1)+273.
      if (help.le.0.) td2(i,j,1,n)=tt2(i,j,1,n)
Matthias Langer's avatar
 
Matthias Langer committed
630
631
632
633
634
635
636
637
638
639
640
641
642
    end do
  end do

  if(levdiff2.eq.0) then
    iwmax=nlev_ec+1
    do i=0,nxmin1
      do j=0,nymin1
        wwh(i,j,nlev_ec+1)=0.
      end do
    end do
  endif


Espen Sollum's avatar
Espen Sollum committed
643
644
645
! For global fields, assign the leftmost data column also to the rightmost
! data column; if required, shift whole grid by nxshift grid points
!*************************************************************************
Matthias Langer's avatar
 
Matthias Langer committed
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674

  if (xglobal) then
    call shift_field_0(ewss,nxfield,ny)
    call shift_field_0(nsss,nxfield,ny)
    call shift_field_0(oro,nxfield,ny)
    call shift_field_0(excessoro,nxfield,ny)
    call shift_field_0(lsm,nxfield,ny)
    call shift_field_0(ulev1,nxfield,ny)
    call shift_field_0(vlev1,nxfield,ny)
    call shift_field_0(tlev1,nxfield,ny)
    call shift_field_0(qvh2,nxfield,ny)
    call shift_field(ps,nxfield,ny,1,1,2,n)
    call shift_field(sd,nxfield,ny,1,1,2,n)
    call shift_field(msl,nxfield,ny,1,1,2,n)
    call shift_field(tcc,nxfield,ny,1,1,2,n)
    call shift_field(u10,nxfield,ny,1,1,2,n)
    call shift_field(v10,nxfield,ny,1,1,2,n)
    call shift_field(tt2,nxfield,ny,1,1,2,n)
    call shift_field(td2,nxfield,ny,1,1,2,n)
    call shift_field(lsprec,nxfield,ny,1,1,2,n)
    call shift_field(convprec,nxfield,ny,1,1,2,n)
    call shift_field(sshf,nxfield,ny,1,1,2,n)
    call shift_field(ssr,nxfield,ny,1,1,2,n)
    call shift_field(hmix,nxfield,ny,1,1,2,n)
    call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n)
    call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n)
    call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
    call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
    call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
Sabine's avatar
Sabine committed
675
676
! IP & SEC adding GFS Clouds 20181205
    call shift_field(clwch,nxfield,ny,nuvzmax,nuvz,2,n)
Matthias Langer's avatar
 
Matthias Langer committed
677
678
679
680
  endif

  do i=0,nxmin1
    do j=0,nymin1
Espen Sollum's avatar
Espen Sollum committed
681
! Convert precip. from mm/s -> mm/hour
Matthias Langer's avatar
 
Matthias Langer committed
682
683
684
685
686
687
688
      convprec(i,j,1,n)=convprec(i,j,1,n)*3600.
      lsprec(i,j,1,n)=lsprec(i,j,1,n)*3600.
      surfstr(i,j,1,n)=sqrt(ewss(i,j)**2+nsss(i,j)**2)
    end do
  end do

  if ((.not.hflswitch).or.(.not.strswitch)) then
Espen Sollum's avatar
Espen Sollum committed
689
690
!  write(*,*) 'WARNING: No flux data contained in GRIB file ',
!    +  wfname(indj)
Matthias Langer's avatar
 
Matthias Langer committed
691

Espen Sollum's avatar
Espen Sollum committed
692
693
! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
!***************************************************************************
Matthias Langer's avatar
 
Matthias Langer committed
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712

    do i=0,nxmin1
      do j=0,nymin1
        hlev1=30.0                     ! HEIGHT OF FIRST MODEL SIGMA LAYER
        ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2)
        fflev1=sqrt(ulev1(i,j)**2+vlev1(i,j)**2)
        call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
             tt2(i,j,1,n),tlev1(i,j),ff10m,fflev1, &
             surfstr(i,j,1,n),sshf(i,j,1,n))
        if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
        if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
      end do
    end do
  endif

  if(iumax.ne.nuvz) stop 'READWIND: NUVZ NOT CONSISTENT'
  if(iumax.ne.nwz)    stop 'READWIND: NWZ NOT CONSISTENT'

  return
Espen Sollum's avatar
Espen Sollum committed
713
888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
Matthias Langer's avatar
 
Matthias Langer committed
714
715
716
  write(*,*) ' #### ',wfname(indj),'                    #### '
  write(*,*) ' #### IS NOT GRIB FORMAT !!!                  #### '
  stop 'Execution terminated'
Espen Sollum's avatar
Espen Sollum committed
717
999 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
Matthias Langer's avatar
 
Matthias Langer committed
718
719
720
721
  write(*,*) ' #### ',wfname(indj),'                    #### '
  write(*,*) ' #### CANNOT BE OPENED !!!                    #### '
  stop 'Execution terminated'

722
end subroutine readwind_gfs