source: flexpart.git/src/readwind_gfs.f90 @ 38b7917

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 38b7917 was 8a65cb0, checked in by Espen Sollum ATMOS <espen@…>, 9 years ago

Added code, makefile for dev branch

  • Property mode set to 100644
File size: 24.8 KB
Line 
1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22subroutine readwind(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  !* INPUT:                                                              *
49  !* indj               indicates number of the wind field to be read in *
50  !* n                  temporal index for meteorological fields (1 to 3)*
51  !*                                                                     *
52  !* IMPORTANT VARIABLES FROM COMMON BLOCK:                              *
53  !*                                                                     *
54  !* wfname             File name of data to be read in                  *
55  !* nx,ny,nuvz,nwz     expected field dimensions                        *
56  !* nlev_ec            number of vertical levels ecmwf model            *
57  !* uu,vv,ww           wind fields                                      *
58  !* tt,qv              temperature and specific humidity                *
59  !* ps                 surface pressure                                 *
60  !*                                                                     *
61  !***********************************************************************
62
63  use grib_api
64  use par_mod
65  use com_mod
66
67  implicit none
68
69  !HSO  new parameters for grib_api
70  integer :: ifile
71  integer :: iret
72  integer :: igrib
73  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
74  !HSO end edits
75  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
76  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
77  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
78  integer :: ii,indj,i,j,k,n,levdiff2,ifield,iumax,iwmax
79
80  ! NCEP
81  integer :: numpt,numpu,numpv,numpw,numprh
82  real :: help, temp, ew
83  real :: elev
84  real :: ulev1(0:nxmax-1,0:nymax-1),vlev1(0:nxmax-1,0:nymax-1)
85  real :: tlev1(0:nxmax-1,0:nymax-1)
86  real :: qvh2(0:nxmax-1,0:nymax-1)
87
88  integer :: i179,i180,i181
89
90  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
91  !HSO kept isec1, isec2 and zsec4 for consistency with gribex GRIB input
92
93  integer :: isec1(8),isec2(3)
94  real(kind=4) :: zsec4(jpunp)
95  real(kind=4) :: xaux,yaux,xaux0,yaux0
96  real(kind=8) :: xauxin,yauxin
97  real,parameter :: eps=1.e-4
98  real(kind=4) :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
99  real :: plev1,hlev1,ff10m,fflev1
100
101  logical :: hflswitch,strswitch
102
103  !HSO  for grib api error messages
104  character(len=24) :: gribErrorMsg = 'Error reading grib file'
105  character(len=20) :: gribFunction = 'readwind_gfs'
106
107
108  hflswitch=.false.
109  strswitch=.false.
110  levdiff2=nlev_ec-nwz+1
111  iumax=0
112  iwmax=0
113
114
115  ! OPENING OF DATA FILE (GRIB CODE)
116
117  !HSO
1185   call grib_open_file(ifile,path(3)(1:length(3)) &
119         //trim(wfname(indj)),'r',iret)
120  if (iret.ne.GRIB_SUCCESS) then
121    goto 888   ! ERROR DETECTED
122  endif
123  !turn on support for multi fields messages
124  call grib_multi_support_on
125
126  numpt=0
127  numpu=0
128  numpv=0
129  numpw=0
130  numprh=0
131  ifield=0
13210   ifield=ifield+1
133  !
134  ! GET NEXT FIELDS
135  !
136  call grib_new_from_file(ifile,igrib,iret)
137  if (iret.eq.GRIB_END_OF_FILE)  then
138    goto 50    ! EOF DETECTED
139  elseif (iret.ne.GRIB_SUCCESS) then
140    goto 888   ! ERROR DETECTED
141  endif
142
143  !first see if we read GRIB1 or GRIB2
144  call grib_get_int(igrib,'editionNumber',gribVer,iret)
145!  call grib_check(iret,gribFunction,gribErrorMsg)
146
147  if (gribVer.eq.1) then ! GRIB Edition 1
148
149  !read the grib1 identifiers
150  call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
151!  call grib_check(iret,gribFunction,gribErrorMsg)
152  call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
153!  call grib_check(iret,gribFunction,gribErrorMsg)
154  call grib_get_int(igrib,'level',isec1(8),iret)
155!  call grib_check(iret,gribFunction,gribErrorMsg)
156
157  else ! GRIB Edition 2
158
159  !read the grib2 identifiers
160  call grib_get_int(igrib,'discipline',discipl,iret)
161!  call grib_check(iret,gribFunction,gribErrorMsg)
162  call grib_get_int(igrib,'parameterCategory',parCat,iret)
163!  call grib_check(iret,gribFunction,gribErrorMsg)
164  call grib_get_int(igrib,'parameterNumber',parNum,iret)
165!  call grib_check(iret,gribFunction,gribErrorMsg)
166  call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
167!  call grib_check(iret,gribFunction,gribErrorMsg)
168  call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
169       valSurf,iret)
170!  call grib_check(iret,gribFunction,gribErrorMsg)
171
172  !convert to grib1 identifiers
173  isec1(6)=-1
174  isec1(7)=-1
175  isec1(8)=-1
176  if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T
177    isec1(6)=11          ! indicatorOfParameter
178    isec1(7)=100         ! indicatorOfTypeOfLevel
179    isec1(8)=valSurf/100 ! level, convert to hPa
180  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
181    isec1(6)=33          ! indicatorOfParameter
182    isec1(7)=100         ! indicatorOfTypeOfLevel
183    isec1(8)=valSurf/100 ! level, convert to hPa
184  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.100)) then ! V
185    isec1(6)=34          ! indicatorOfParameter
186    isec1(7)=100         ! indicatorOfTypeOfLevel
187    isec1(8)=valSurf/100 ! level, convert to hPa
188  elseif ((parCat.eq.2).and.(parNum.eq.8).and.(typSurf.eq.100)) then ! W
189    isec1(6)=39          ! indicatorOfParameter
190    isec1(7)=100         ! indicatorOfTypeOfLevel
191    isec1(8)=valSurf/100 ! level, convert to hPa
192  elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.100)) then ! RH
193    isec1(6)=52          ! indicatorOfParameter
194    isec1(7)=100         ! indicatorOfTypeOfLevel
195    isec1(8)=valSurf/100 ! level, convert to hPa
196  elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.103)) then ! RH2
197    isec1(6)=52          ! indicatorOfParameter
198    isec1(7)=105         ! indicatorOfTypeOfLevel
199    isec1(8)=2
200  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! T2
201    isec1(6)=11          ! indicatorOfParameter
202    isec1(7)=105         ! indicatorOfTypeOfLevel
203    isec1(8)=2
204  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! U10
205    isec1(6)=33          ! indicatorOfParameter
206    isec1(7)=105         ! indicatorOfTypeOfLevel
207    isec1(8)=10
208  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! V10
209    isec1(6)=34          ! indicatorOfParameter
210    isec1(7)=105         ! indicatorOfTypeOfLevel
211    isec1(8)=10
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    isec1(8)=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    isec1(8)=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    isec1(8)=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    isec1(8)=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    isec1(8)=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    isec1(8)=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    isec1(8)=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    isec1(8)=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    isec1(8)=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    isec1(8)=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    isec1(8)=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=-179.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) &
298         stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
299    if(abs(yaux-yaux0).gt.eps) &
300         stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
301  endif
302  !HSO end of edits
303
304  i179=nint(179./dx)
305  if (dx.lt.0.7) then
306    i180=nint(180./dx)+1    ! 0.5 deg data
307  else
308    i180=nint(179./dx)+1    ! 1 deg data
309  endif
310  i181=i180+1
311
312  if (isec1(6).ne.-1) then
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
320              if ((isec1(8)*100.0).eq.akz(ii)) numpt=ii
321            end do
322        endif
323        help=zsec4(nxfield*(ny-j-1)+i+1)
324        if(i.le.i180) then
325          tth(i179+i,j,numpt,n)=help
326        else
327          tth(i-i181,j,numpt,n)=help
328        endif
329      endif
330      if((isec1(6).eq.033).and.(isec1(7).eq.100)) then
331  ! U VELOCITY
332         if((i.eq.0).and.(j.eq.0)) then
333            do ii=1,nuvz
334              if ((isec1(8)*100.0).eq.akz(ii)) numpu=ii
335            end do
336        endif
337        help=zsec4(nxfield*(ny-j-1)+i+1)
338        if(i.le.i180) then
339          uuh(i179+i,j,numpu)=help
340        else
341          uuh(i-i181,j,numpu)=help
342        endif
343      endif
344      if((isec1(6).eq.034).and.(isec1(7).eq.100)) then
345  ! V VELOCITY
346         if((i.eq.0).and.(j.eq.0)) then
347            do ii=1,nuvz
348              if ((isec1(8)*100.0).eq.akz(ii)) numpv=ii
349            end do
350        endif
351        help=zsec4(nxfield*(ny-j-1)+i+1)
352        if(i.le.i180) then
353          vvh(i179+i,j,numpv)=help
354        else
355          vvh(i-i181,j,numpv)=help
356        endif
357      endif
358      if((isec1(6).eq.052).and.(isec1(7).eq.100)) then
359  ! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER
360         if((i.eq.0).and.(j.eq.0)) then
361            do ii=1,nuvz
362              if ((isec1(8)*100.0).eq.akz(ii)) numprh=ii
363            end do
364        endif
365        help=zsec4(nxfield*(ny-j-1)+i+1)
366        if(i.le.i180) then
367          qvh(i179+i,j,numprh,n)=help
368        else
369          qvh(i-i181,j,numprh,n)=help
370        endif
371      endif
372      if((isec1(6).eq.001).and.(isec1(7).eq.001)) then
373  ! SURFACE PRESSURE
374        help=zsec4(nxfield*(ny-j-1)+i+1)
375        if(i.le.i180) then
376          ps(i179+i,j,1,n)=help
377        else
378          ps(i-i181,j,1,n)=help
379        endif
380      endif
381      if((isec1(6).eq.039).and.(isec1(7).eq.100)) then
382  ! W VELOCITY
383         if((i.eq.0).and.(j.eq.0)) then
384            do ii=1,nuvz
385              if ((isec1(8)*100.0).eq.akz(ii)) numpw=ii
386            end do
387        endif
388        help=zsec4(nxfield*(ny-j-1)+i+1)
389        if(i.le.i180) then
390          wwh(i179+i,j,numpw)=help
391        else
392          wwh(i-i181,j,numpw)=help
393        endif
394      endif
395      if((isec1(6).eq.066).and.(isec1(7).eq.001)) then
396  ! SNOW DEPTH
397        help=zsec4(nxfield*(ny-j-1)+i+1)
398        if(i.le.i180) then
399          sd(i179+i,j,1,n)=help
400        else
401          sd(i-i181,j,1,n)=help
402        endif
403      endif
404      if((isec1(6).eq.002).and.(isec1(7).eq.102)) then
405  ! MEAN SEA LEVEL PRESSURE
406        help=zsec4(nxfield*(ny-j-1)+i+1)
407        if(i.le.i180) then
408          msl(i179+i,j,1,n)=help
409        else
410          msl(i-i181,j,1,n)=help
411        endif
412      endif
413      if((isec1(6).eq.071).and.(isec1(7).eq.244)) then
414  ! TOTAL CLOUD COVER
415        help=zsec4(nxfield*(ny-j-1)+i+1)
416        if(i.le.i180) then
417          tcc(i179+i,j,1,n)=help
418        else
419          tcc(i-i181,j,1,n)=help
420        endif
421      endif
422      if((isec1(6).eq.033).and.(isec1(7).eq.105).and. &
423           (isec1(8).eq.10)) then
424  ! 10 M U VELOCITY
425        help=zsec4(nxfield*(ny-j-1)+i+1)
426        if(i.le.i180) then
427          u10(i179+i,j,1,n)=help
428        else
429          u10(i-i181,j,1,n)=help
430        endif
431      endif
432      if((isec1(6).eq.034).and.(isec1(7).eq.105).and. &
433           (isec1(8).eq.10)) then
434  ! 10 M V VELOCITY
435        help=zsec4(nxfield*(ny-j-1)+i+1)
436        if(i.le.i180) then
437          v10(i179+i,j,1,n)=help
438        else
439          v10(i-i181,j,1,n)=help
440        endif
441      endif
442      if((isec1(6).eq.011).and.(isec1(7).eq.105).and. &
443           (isec1(8).eq.02)) then
444  ! 2 M TEMPERATURE
445        help=zsec4(nxfield*(ny-j-1)+i+1)
446        if(i.le.i180) then
447          tt2(i179+i,j,1,n)=help
448        else
449          tt2(i-i181,j,1,n)=help
450        endif
451      endif
452      if((isec1(6).eq.017).and.(isec1(7).eq.105).and. &
453           (isec1(8).eq.02)) then
454  ! 2 M DEW POINT TEMPERATURE
455        help=zsec4(nxfield*(ny-j-1)+i+1)
456        if(i.le.i180) then
457          td2(i179+i,j,1,n)=help
458        else
459          td2(i-i181,j,1,n)=help
460        endif
461      endif
462      if((isec1(6).eq.062).and.(isec1(7).eq.001)) then
463  ! LARGE SCALE PREC.
464        help=zsec4(nxfield*(ny-j-1)+i+1)
465        if(i.le.i180) then
466          lsprec(i179+i,j,1,n)=help
467        else
468          lsprec(i-i181,j,1,n)=help
469        endif
470      endif
471      if((isec1(6).eq.063).and.(isec1(7).eq.001)) then
472  ! CONVECTIVE PREC.
473        help=zsec4(nxfield*(ny-j-1)+i+1)
474        if(i.le.i180) then
475          convprec(i179+i,j,1,n)=help
476        else
477          convprec(i-i181,j,1,n)=help
478        endif
479      endif
480      if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
481  ! TOPOGRAPHY
482        help=zsec4(nxfield*(ny-j-1)+i+1)
483        if(i.le.i180) then
484          oro(i179+i,j)=help
485          excessoro(i179+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
486        else
487          oro(i-i181,j)=help
488          excessoro(i-i181,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
489        endif
490      endif
491      if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
492  ! LAND SEA MASK
493        help=zsec4(nxfield*(ny-j-1)+i+1)
494        if(i.le.i180) then
495          lsm(i179+i,j)=help
496        else
497          lsm(i-i181,j)=help
498        endif
499      endif
500      if((isec1(6).eq.221).and.(isec1(7).eq.001)) then
501  ! MIXING HEIGHT
502        help=zsec4(nxfield*(ny-j-1)+i+1)
503        if(i.le.i180) then
504          hmix(i179+i,j,1,n)=help
505        else
506          hmix(i-i181,j,1,n)=help
507        endif
508      endif
509      if((isec1(6).eq.052).and.(isec1(7).eq.105).and. &
510           (isec1(8).eq.02)) then
511  ! 2 M RELATIVE HUMIDITY
512        help=zsec4(nxfield*(ny-j-1)+i+1)
513        if(i.le.i180) then
514          qvh2(i179+i,j)=help
515        else
516          qvh2(i-i181,j)=help
517        endif
518      endif
519      if((isec1(6).eq.011).and.(isec1(7).eq.107)) then
520  ! TEMPERATURE LOWEST SIGMA LEVEL
521        help=zsec4(nxfield*(ny-j-1)+i+1)
522        if(i.le.i180) then
523          tlev1(i179+i,j)=help
524        else
525          tlev1(i-i181,j)=help
526        endif
527      endif
528      if((isec1(6).eq.033).and.(isec1(7).eq.107)) then
529  ! U VELOCITY LOWEST SIGMA LEVEL
530        help=zsec4(nxfield*(ny-j-1)+i+1)
531        if(i.le.i180) then
532          ulev1(i179+i,j)=help
533        else
534          ulev1(i-i181,j)=help
535        endif
536      endif
537      if((isec1(6).eq.034).and.(isec1(7).eq.107)) then
538  ! V VELOCITY LOWEST SIGMA LEVEL
539        help=zsec4(nxfield*(ny-j-1)+i+1)
540        if(i.le.i180) then
541          vlev1(i179+i,j)=help
542        else
543          vlev1(i-i181,j)=help
544        endif
545      endif
546
547    end do
548  end do
549
550  endif
551
552  if((isec1(6).eq.33).and.(isec1(7).eq.100)) then
553  ! NCEP ISOBARIC LEVELS
554    iumax=iumax+1
555  endif
556
557  call grib_release(igrib)
558  goto 10                      !! READ NEXT LEVEL OR PARAMETER
559  !
560  ! CLOSING OF INPUT DATA FILE
561  !
562
563  !HSO close grib file
56450   continue
565  call grib_close_file(ifile)
566
567  ! SENS. HEAT FLUX
568  sshf(:,:,1,n)=0.0     ! not available from gfs.tccz.pgrbfxx files
569  hflswitch=.false.    ! Heat flux not available
570  ! SOLAR RADIATIVE FLUXES
571  ssr(:,:,1,n)=0.0      ! not available from gfs.tccz.pgrbfxx files
572  ! EW SURFACE STRESS
573  ewss=0.0         ! not available from gfs.tccz.pgrbfxx files
574  ! NS SURFACE STRESS
575  nsss=0.0         ! not available from gfs.tccz.pgrbfxx files
576  strswitch=.false.    ! stress not available
577
578  ! CONVERT TP TO LSP (GRIB2 only)
579  if (gribVer.eq.2) then
580    do j=0,nymin1
581    do i=0,nxfield-1
582     if(i.le.i180) then
583     if (convprec(i179+i,j,1,n).lt.lsprec(i179+i,j,1,n)) then ! neg precip would occur
584         lsprec(i179+i,j,1,n)= &
585              lsprec(i179+i,j,1,n)-convprec(i179+i,j,1,n)
586     else
587         lsprec(i179+i,j,1,n)=0
588     endif
589     else
590     if (convprec(i-i181,j,1,n).lt.lsprec(i-i181,j,1,n)) then
591          lsprec(i-i181,j,1,n)= &
592               lsprec(i-i181,j,1,n)-convprec(i-i181,j,1,n)
593     else
594          lsprec(i-i181,j,1,n)=0
595     endif
596     endif
597    enddo
598    enddo
599  endif
600  !HSO end edits
601
602
603  ! TRANSFORM RH TO SPECIFIC HUMIDITY
604
605  do j=0,ny-1
606    do i=0,nxfield-1
607      do k=1,nuvz
608        help=qvh(i,j,k,n)
609        temp=tth(i,j,k,n)
610        plev1=akm(k)+bkm(k)*ps(i,j,1,n)
611        elev=ew(temp)*help/100.0
612        qvh(i,j,k,n)=xmwml*(elev/(plev1-((1.0-xmwml)*elev)))
613      end do
614    end do
615  end do
616
617  ! CALCULATE 2 M DEW POINT FROM 2 M RELATIVE HUMIDITY
618  ! USING BOLTON'S (1980) FORMULA
619  ! BECAUSE td2 IS NOT AVAILABLE FROM NCEP GFS DATA
620
621  do j=0,ny-1
622    do i=0,nxfield-1
623        help=qvh2(i,j)
624        temp=tt2(i,j,1,n)
625        elev=ew(temp)/100.*help/100.   !vapour pressure in hPa
626        td2(i,j,1,n)=243.5/(17.67/log(elev/6.112)-1)+273.
627        if (help.le.0.) td2(i,j,1,n)=tt2(i,j,1,n)
628    end do
629  end do
630
631  if(levdiff2.eq.0) then
632    iwmax=nlev_ec+1
633    do i=0,nxmin1
634      do j=0,nymin1
635        wwh(i,j,nlev_ec+1)=0.
636      end do
637    end do
638  endif
639
640
641  ! For global fields, assign the leftmost data column also to the rightmost
642  ! data column; if required, shift whole grid by nxshift grid points
643  !*************************************************************************
644
645  if (xglobal) then
646    call shift_field_0(ewss,nxfield,ny)
647    call shift_field_0(nsss,nxfield,ny)
648    call shift_field_0(oro,nxfield,ny)
649    call shift_field_0(excessoro,nxfield,ny)
650    call shift_field_0(lsm,nxfield,ny)
651    call shift_field_0(ulev1,nxfield,ny)
652    call shift_field_0(vlev1,nxfield,ny)
653    call shift_field_0(tlev1,nxfield,ny)
654    call shift_field_0(qvh2,nxfield,ny)
655    call shift_field(ps,nxfield,ny,1,1,2,n)
656    call shift_field(sd,nxfield,ny,1,1,2,n)
657    call shift_field(msl,nxfield,ny,1,1,2,n)
658    call shift_field(tcc,nxfield,ny,1,1,2,n)
659    call shift_field(u10,nxfield,ny,1,1,2,n)
660    call shift_field(v10,nxfield,ny,1,1,2,n)
661    call shift_field(tt2,nxfield,ny,1,1,2,n)
662    call shift_field(td2,nxfield,ny,1,1,2,n)
663    call shift_field(lsprec,nxfield,ny,1,1,2,n)
664    call shift_field(convprec,nxfield,ny,1,1,2,n)
665    call shift_field(sshf,nxfield,ny,1,1,2,n)
666    call shift_field(ssr,nxfield,ny,1,1,2,n)
667    call shift_field(hmix,nxfield,ny,1,1,2,n)
668    call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n)
669    call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n)
670    call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
671    call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
672    call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
673  endif
674
675  do i=0,nxmin1
676    do j=0,nymin1
677  ! Convert precip. from mm/s -> mm/hour
678      convprec(i,j,1,n)=convprec(i,j,1,n)*3600.
679      lsprec(i,j,1,n)=lsprec(i,j,1,n)*3600.
680      surfstr(i,j,1,n)=sqrt(ewss(i,j)**2+nsss(i,j)**2)
681    end do
682  end do
683
684  if ((.not.hflswitch).or.(.not.strswitch)) then
685  !  write(*,*) 'WARNING: No flux data contained in GRIB file ',
686  !    +  wfname(indj)
687
688  ! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
689  !***************************************************************************
690
691    do i=0,nxmin1
692      do j=0,nymin1
693        hlev1=30.0                     ! HEIGHT OF FIRST MODEL SIGMA LAYER
694        ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2)
695        fflev1=sqrt(ulev1(i,j)**2+vlev1(i,j)**2)
696        call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
697             tt2(i,j,1,n),tlev1(i,j),ff10m,fflev1, &
698             surfstr(i,j,1,n),sshf(i,j,1,n))
699        if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
700        if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
701      end do
702    end do
703  endif
704
705  if(iumax.ne.nuvz) stop 'READWIND: NUVZ NOT CONSISTENT'
706  if(iumax.ne.nwz)    stop 'READWIND: NWZ NOT CONSISTENT'
707
708  return
709888   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
710  write(*,*) ' #### ',wfname(indj),'                    #### '
711  write(*,*) ' #### IS NOT GRIB FORMAT !!!                  #### '
712  stop 'Execution terminated'
713999   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
714  write(*,*) ' #### ',wfname(indj),'                    #### '
715  write(*,*) ' #### CANNOT BE OPENED !!!                    #### '
716  stop 'Execution terminated'
717
718end subroutine readwind
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG