source: branches/jerome/src/readwind_gfs.f90 @ 15

Last change on this file since 15 was 15, checked in by jebri, 9 years ago

fixes for flexpart91

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