source: trunk/src/readwind_gfs.f90 @ 28

Last change on this file since 28 was 20, checked in by igpis, 10 years ago

move version 9.1.8 form branches to trunk. Contributions from HSO, saeck, pesei, NIK, RT, XKF, IP and others

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