source: flexpart.git/src/readwind_ecmwf.f90 @ 02095e3

devrelease-10univie
Last change on this file since 02095e3 was 61e07ba, checked in by Espen Sollum ATMOS <eso@…>, 2 years ago

Adding files from CTBTO project wo_17

  • Property mode set to 100644
File size: 20.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_ecmwf(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: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
34!                                 ECMWF grib_api                      *
35!             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
36!                                 ECMWF grib_api                      *
37!                                                                     *
38!**********************************************************************
39!  Changes, Bernd C. Krueger, Feb. 2001:
40!   Variables tth and qvh (on eta coordinates) in common block
41!
42!   Unified ECMWF and GFS builds                                     
43!   Marian Harustak, 12.5.2017                                       
44!     - Renamed from readwind to readwind_ecmwf                     
45!
46!**********************************************************************
47!                                                                     *
48! DESCRIPTION:                                                        *
49!                                                                     *
50! READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE   *
51! INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE          *
52!                                                                     *
53! INPUT:                                                              *
54! indj               indicates number of the wind field to be read in *
55! n                  temporal index for meteorological fields (1 to 3)*
56!                                                                     *
57! IMPORTANT VARIABLES FROM COMMON BLOCK:                              *
58!                                                                     *
59! wfname             File name of data to be read in                  *
60! nx,ny,nuvz,nwz     expected field dimensions                        *
61! nlev_ec            number of vertical levels ecmwf model            *
62! uu,vv,ww           wind fields                                      *
63! tt,qv              temperature and specific humidity                *
64! ps                 surface pressure                                 *
65!                                                                     *
66!**********************************************************************
67
68  use grib_api
69  use par_mod
70  use com_mod
71
72  implicit none
73
74!  include 'grib_api.h'
75
76!HSO  parameters for grib_api
77  integer :: ifile
78  integer :: iret
79  integer :: igrib
80  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parId
81  integer :: gotGrid
82!HSO  end
83
84  real(sp) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
85  real(sp) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
86  real(sp) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
87  integer :: indj,i,j,k,n,levdiff2,ifield,iumax,iwmax
88
89! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
90
91! dimension of isec2 at least (22+n), where n is the number of parallels or
92! meridians in a quasi-regular (reduced) Gaussian or lat/long grid
93
94! dimension of zsec2 at least (10+nn), where nn is the number of vertical
95! coordinate parameters
96
97  integer :: isec1(56),isec2(22+nxmax+nymax)
98  real(sp) :: zsec4(jpunp)
99  real(sp) :: xaux,yaux
100  real(dp) :: xauxin,yauxin
101  real(sp),parameter :: eps=1.e-4
102  real(sp) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
103  real(sp) :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
104
105  logical :: hflswitch,strswitch!,readcloud
106
107!HSO  grib api error messages
108  character(len=24) :: gribErrorMsg = 'Error reading grib file'
109  character(len=20) :: gribFunction = 'readwind'
110
111  hflswitch=.false.
112  strswitch=.false.
113!ZHG test the grib fields that have lcwc without using them
114!  readcloud=.false.
115
116  levdiff2=nlev_ec-nwz+1
117  iumax=0
118  iwmax=0
119
120!
121! OPENING OF DATA FILE (GRIB CODE)
122!
1235 call grib_open_file(ifile,path(3)(1:length(3)) &
124       //trim(wfname(indj)),'r',iret)
125  if (iret.ne.GRIB_SUCCESS) then
126    goto 888   ! ERROR DETECTED
127  endif
128!turn on support for multi fields messages */
129!call grib_multi_support_on
130
131  gotGrid=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!print*,'GRiB Edition 1'
151!read the grib2 identifiers
152    call grib_get_int(igrib,'indicatorOfParameter',isec1(6),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!change code for etadot to code for omega
158    if (isec1(6).eq.77) then
159      isec1(6)=135
160    endif
161
162    conversion_factor=1.
163
164  else
165
166!print*,'GRiB Edition 2'
167!read the grib2 identifiers
168    call grib_get_int(igrib,'discipline',discipl,iret)
169    call grib_check(iret,gribFunction,gribErrorMsg)
170    call grib_get_int(igrib,'parameterCategory',parCat,iret)
171    call grib_check(iret,gribFunction,gribErrorMsg)
172    call grib_get_int(igrib,'parameterNumber',parNum,iret)
173    call grib_check(iret,gribFunction,gribErrorMsg)
174    call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
175    call grib_check(iret,gribFunction,gribErrorMsg)
176    call grib_get_int(igrib,'level',valSurf,iret)
177    call grib_check(iret,gribFunction,gribErrorMsg)
178    call grib_get_int(igrib,'paramId',parId,iret)
179    call grib_check(iret,gribFunction,gribErrorMsg)
180
181!print*,discipl,parCat,parNum,typSurf,valSurf
182
183!convert to grib1 identifiers
184    isec1(6)=-1
185    isec1(7)=-1
186    isec1(8)=-1
187    isec1(8)=valSurf     ! level
188    conversion_factor=1.
189    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
190      isec1(6)=130         ! indicatorOfParameter
191    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
192      isec1(6)=131         ! indicatorOfParameter
193    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
194      isec1(6)=132         ! indicatorOfParameter
195    elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
196      isec1(6)=133         ! indicatorOfParameter
197! ESO Cloud water is in a) fields CLWC and CIWC, *or* b) field QC
198    elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
199      isec1(6)=246         ! indicatorOfParameter
200    elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
201      isec1(6)=247         ! indicatorOfParameter
202! ESO qc(=clwc+ciwc):
203    elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
204      isec1(6)=201031         ! indicatorOfParameter
205    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
206      isec1(6)=134         ! indicatorOfParameter
207    elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
208      isec1(6)=135         ! indicatorOfParameter
209    elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
210      isec1(6)=135         ! indicatorOfParameter
211    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
212      isec1(6)=151         ! indicatorOfParameter
213    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
214      isec1(6)=165         ! indicatorOfParameter
215    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
216      isec1(6)=166         ! indicatorOfParameter
217    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
218      isec1(6)=167         ! indicatorOfParameter
219    elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
220      isec1(6)=168         ! indicatorOfParameter
221    elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
222      isec1(6)=141         ! indicatorOfParameter
223      conversion_factor=1000.
224    elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
225      isec1(6)=164         ! indicatorOfParameter
226    elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
227      isec1(6)=142         ! indicatorOfParameter
228    elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
229      isec1(6)=143         ! indicatorOfParameter
230      conversion_factor=1000.
231    elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
232      isec1(6)=146         ! indicatorOfParameter
233    elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
234      isec1(6)=176         ! indicatorOfParameter
235!    elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS --wrong
236    elseif ((parCat.eq.2).and.(parNum.eq.38) .or. parId .eq. 180) then ! EWSS --correct
237      isec1(6)=180         ! indicatorOfParameter
238!    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS --wrong
239    elseif ((parCat.eq.2).and.(parNum.eq.37) .or. parId .eq. 181) then ! NSSS --correct
240      isec1(6)=181         ! indicatorOfParameter
241    elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
242      isec1(6)=129         ! indicatorOfParameter
243    elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
244      isec1(6)=160         ! indicatorOfParameter
245    elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
246         (typSurf.eq.1)) then ! LSM
247      isec1(6)=172         ! indicatorOfParameter
248    else
249      print*,'***WARNING: undefined GRiB2 message found!',discipl, &
250           parCat,parNum,typSurf
251    endif
252    if(parId .ne. isec1(6) .and. parId .ne. 77) then
253      write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
254!    stop
255    endif
256
257  endif
258
259!HSO  get the size and data of the values array
260  if (isec1(6).ne.-1) then
261    call grib_get_real4_array(igrib,'values',zsec4,iret)
262    call grib_check(iret,gribFunction,gribErrorMsg)
263  endif
264
265!HSO  get the required fields from section 2 in a gribex compatible manner
266  if (ifield.eq.1) then
267    call grib_get_int(igrib,'numberOfPointsAlongAParallel',isec2(2),iret)
268    call grib_check(iret,gribFunction,gribErrorMsg)
269    call grib_get_int(igrib,'numberOfPointsAlongAMeridian',isec2(3),iret)
270    call grib_check(iret,gribFunction,gribErrorMsg)
271    call grib_get_int(igrib,'numberOfVerticalCoordinateValues',isec2(12))
272    call grib_check(iret,gribFunction,gribErrorMsg)
273! CHECK GRID SPECIFICATIONS
274    if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
275    if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
276    if(isec2(12)/2-1.ne.nlev_ec) &
277         stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
278  endif ! ifield
279
280!HSO  get the second part of the grid dimensions only from GRiB1 messages
281  if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
282    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
283         xauxin,iret)
284    call grib_check(iret,gribFunction,gribErrorMsg)
285    call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
286         yauxin,iret)
287    call grib_check(iret,gribFunction,gribErrorMsg)
288    if (xauxin.gt.180.) xauxin=xauxin-360.0
289    if (xauxin.lt.-180.) xauxin=xauxin+360.0
290
291    xaux=xauxin+real(nxshift)*dx
292    yaux=yauxin
293    if (xaux.gt.180.) xaux=xaux-360.0
294    if(abs(xaux-xlon0).gt.eps) &
295         stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
296    if(abs(yaux-ylat0).gt.eps) &
297         stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
298    gotGrid=1
299  endif ! gotGrid
300
301  do j=0,nymin1
302    do i=0,nxfield-1
303      k=isec1(8)
304      if(isec1(6).eq.130) tth(i,j,nlev_ec-k+2,n)= &!! TEMPERATURE
305           zsec4(nxfield*(ny-j-1)+i+1)
306      if(isec1(6).eq.131) uuh(i,j,nlev_ec-k+2)= &!! U VELOCITY
307           zsec4(nxfield*(ny-j-1)+i+1)
308      if(isec1(6).eq.132) vvh(i,j,nlev_ec-k+2)= &!! V VELOCITY
309           zsec4(nxfield*(ny-j-1)+i+1)
310      if(isec1(6).eq.133) then                      !! SPEC. HUMIDITY
311        qvh(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
312        if (qvh(i,j,nlev_ec-k+2,n) .lt. 0.) &
313             qvh(i,j,nlev_ec-k+2,n) = 0.
314!        this is necessary because the gridded data may contain
315!        spurious negative values
316      endif
317      if(isec1(6).eq.134) ps(i,j,1,n)= &!! SURF. PRESS.
318           zsec4(nxfield*(ny-j-1)+i+1)
319
320      if(isec1(6).eq.135) wwh(i,j,nlev_ec-k+1)= &!! W VELOCITY
321           zsec4(nxfield*(ny-j-1)+i+1)
322      if(isec1(6).eq.141) sd(i,j,1,n)= &!! SNOW DEPTH
323           zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
324      if(isec1(6).eq.151) msl(i,j,1,n)= &!! SEA LEVEL PRESS.
325           zsec4(nxfield*(ny-j-1)+i+1)
326      if(isec1(6).eq.164) tcc(i,j,1,n)= &!! CLOUD COVER
327           zsec4(nxfield*(ny-j-1)+i+1)
328      if(isec1(6).eq.165) u10(i,j,1,n)= &!! 10 M U VELOCITY
329           zsec4(nxfield*(ny-j-1)+i+1)
330      if(isec1(6).eq.166) v10(i,j,1,n)= &!! 10 M V VELOCITY
331           zsec4(nxfield*(ny-j-1)+i+1)
332      if(isec1(6).eq.167) tt2(i,j,1,n)= &!! 2 M TEMPERATURE
333           zsec4(nxfield*(ny-j-1)+i+1)
334      if(isec1(6).eq.168) td2(i,j,1,n)= &!! 2 M DEW POINT
335           zsec4(nxfield*(ny-j-1)+i+1)
336      if(isec1(6).eq.142) then                      !! LARGE SCALE PREC.
337        lsprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
338        if (lsprec(i,j,1,n).lt.0.) lsprec(i,j,1,n)=0.
339      endif
340      if(isec1(6).eq.143) then                      !! CONVECTIVE PREC.
341        convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
342        if (convprec(i,j,1,n).lt.0.) convprec(i,j,1,n)=0.
343      endif
344      if(isec1(6).eq.146) sshf(i,j,1,n)= &!! SENS. HEAT FLUX
345           zsec4(nxfield*(ny-j-1)+i+1)
346      if((isec1(6).eq.146).and.(zsec4(nxfield*(ny-j-1)+i+1).ne.0.)) &
347           hflswitch=.true.    ! Heat flux available
348      if(isec1(6).eq.176) then                      !! SOLAR RADIATION
349        ssr(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
350        if (ssr(i,j,1,n).lt.0.) ssr(i,j,1,n)=0.
351      endif
352      if(isec1(6).eq.180) ewss(i,j)= &!! EW SURFACE STRESS
353           zsec4(nxfield*(ny-j-1)+i+1)
354      if(isec1(6).eq.181) nsss(i,j)= &!! NS SURFACE STRESS
355           zsec4(nxfield*(ny-j-1)+i+1)
356      if(((isec1(6).eq.180).or.(isec1(6).eq.181)).and. &
357           (zsec4(nxfield*(ny-j-1)+i+1).ne.0.)) strswitch=.true.    ! stress available
358!sec        strswitch=.true.
359      if(isec1(6).eq.129) oro(i,j)= &!! ECMWF OROGRAPHY
360           zsec4(nxfield*(ny-j-1)+i+1)/ga
361      if(isec1(6).eq.160) excessoro(i,j)= &!! STANDARD DEVIATION OF OROGRAPHY
362           zsec4(nxfield*(ny-j-1)+i+1)
363      if(isec1(6).eq.172) lsm(i,j)= &!! ECMWF LAND SEA MASK
364           zsec4(nxfield*(ny-j-1)+i+1)
365      if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
366      if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
367!ZHG READING CLOUD FIELDS ASWELL
368! ESO TODO: add check for if one of clwc/ciwc missing (error),
369! also if all 3 cw fields present, use qc and disregard the others
370      if(isec1(6).eq.246) then  !! CLWC  Cloud liquid water content [kg/kg]
371        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
372        readclouds=.true.
373        sumclouds=.false.
374      endif
375      if(isec1(6).eq.247) then  !! CIWC  Cloud ice water content
376        ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
377      endif
378!ZHG end
379!ESO read qc (=clwc+ciwc)
380      if(isec1(6).eq.201031) then  !! QC  Cloud liquid water content [kg/kg]
381        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
382        readclouds=.true.
383        sumclouds=.true.
384      endif
385
386    end do
387  end do
388
389  call grib_release(igrib)
390  goto 10                      !! READ NEXT LEVEL OR PARAMETER
391!
392! CLOSING OF INPUT DATA FILE
393!
394
39550 call grib_close_file(ifile)
396
397!error message if no fields found with correct first longitude in it
398  if (gotGrid.eq.0) then
399    print*,'***ERROR: input file needs to contain GRiB1 formatted'// &
400         'messages'
401    stop
402  endif
403
404  if(levdiff2.eq.0) then
405    iwmax=nlev_ec+1
406    do i=0,nxmin1
407      do j=0,nymin1
408        wwh(i,j,nlev_ec+1)=0.
409      end do
410    end do
411  endif
412
413! For global fields, assign the leftmost data column also to the rightmost
414! data column; if required, shift whole grid by nxshift grid points
415!*************************************************************************
416
417  if (xglobal) then
418    call shift_field_0(ewss,nxfield,ny)
419    call shift_field_0(nsss,nxfield,ny)
420    call shift_field_0(oro,nxfield,ny)
421    call shift_field_0(excessoro,nxfield,ny)
422    call shift_field_0(lsm,nxfield,ny)
423    call shift_field(ps,nxfield,ny,1,1,2,n)
424    call shift_field(sd,nxfield,ny,1,1,2,n)
425    call shift_field(msl,nxfield,ny,1,1,2,n)
426    call shift_field(tcc,nxfield,ny,1,1,2,n)
427    call shift_field(u10,nxfield,ny,1,1,2,n)
428    call shift_field(v10,nxfield,ny,1,1,2,n)
429    call shift_field(tt2,nxfield,ny,1,1,2,n)
430    call shift_field(td2,nxfield,ny,1,1,2,n)
431    call shift_field(lsprec,nxfield,ny,1,1,2,n)
432    call shift_field(convprec,nxfield,ny,1,1,2,n)
433    call shift_field(sshf,nxfield,ny,1,1,2,n)
434    call shift_field(ssr,nxfield,ny,1,1,2,n)
435    call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n)
436    call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n)
437    call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
438    call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
439    call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
440!ZHG
441    call shift_field(clwch,nxfield,ny,nuvzmax,nuvz,2,n)
442    if (.not.sumclouds) call shift_field(ciwch,nxfield,ny,nuvzmax,nuvz,2,n)
443!ZHG end
444
445  endif
446
447  do i=0,nxmin1
448    do j=0,nymin1
449      surfstr(i,j,1,n)=sqrt(ewss(i,j)**2+nsss(i,j)**2)
450    end do
451  end do
452
453  if ((.not.hflswitch).or.(.not.strswitch)) then
454    write(*,*) 'WARNING: No flux data contained in GRIB file ', &
455         wfname(indj)
456
457! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
458! As ECMWF has increased the model resolution, such that now the first model
459! level is at about 10 m (where 10-m wind is given), use the 2nd ECMWF level
460! (3rd model level in FLEXPART) for the profile method
461!***************************************************************************
462
463    do i=0,nxmin1
464      do j=0,nymin1
465        plev1=akz(3)+bkz(3)*ps(i,j,1,n)
466        pmean=0.5*(ps(i,j,1,n)+plev1)
467        tv=tth(i,j,3,n)*(1.+0.61*qvh(i,j,3,n))
468        fu=-r_air*tv/ga/pmean
469        hlev1=fu*(plev1-ps(i,j,1,n))   ! HEIGTH OF FIRST MODEL LAYER
470        ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2)
471        fflev1=sqrt(uuh(i,j,3)**2+vvh(i,j,3)**2)
472        call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
473             tt2(i,j,1,n),tth(i,j,3,n),ff10m,fflev1, &
474             surfstr(i,j,1,n),sshf(i,j,1,n))
475        if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
476        if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
477      end do
478    end do
479  endif
480
481
482! Assign 10 m wind to model level at eta=1.0 to have one additional model
483! level at the ground
484! Specific humidity is taken the same as at one level above
485! Temperature is taken as 2 m temperature
486!**************************************************************************
487
488  do i=0,nxmin1
489    do j=0,nymin1
490      uuh(i,j,1)=u10(i,j,1,n)
491      vvh(i,j,1)=v10(i,j,1,n)
492      qvh(i,j,1,n)=qvh(i,j,2,n)
493      tth(i,j,1,n)=tt2(i,j,1,n)
494    end do
495  end do
496
497  if(iumax.ne.nuvz-1) stop 'READWIND: NUVZ NOT CONSISTENT'
498  if(iwmax.ne.nwz)    stop 'READWIND: NWZ NOT CONSISTENT'
499
500  return
501
502888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
503  write(*,*) ' #### ',wfname(indj),'                    #### '
504  write(*,*) ' #### IS NOT GRIB FORMAT !!!                  #### '
505  stop 'Execution terminated'
506
507end subroutine readwind_ecmwf
508
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG