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

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 02095e3 was 61e07ba, checked in by Espen Sollum ATMOS <eso@…>, 7 years ago

Adding files from CTBTO project wo_17

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