source: flexpart.git/src/readwind_nests.f90

bugfixes+enhancements
Last change on this file was dba4221, checked in by Pirmin Kaufmann <pirmin.kaufmann@…>, 17 months ago

Bugfixes:

  • options/SPECIES/SPECIES_009: corrected wrong number format, replaced comma by decimal point
  • options/SPECIES/SPECIES_028: corrected wrong number format, moved sign of exponent to after the E
  • options/SPECIES/specoverview.f90: added namelist parameters that appear in SPECIES files but were missing here
  • src/FLEXPART.f90: replaced compiler-specific command line argument routines by standard Fortran intrinsic routines
  • src/FLEXPART_MPI.f90: ditto
  • src/gridcheck_ecmwf.f90: corrected handling of vertical levels when input files do not contain uppermost layers
  • src/gridcheck_nests.f90: ditto
  • src/readwind_ecmwf.f90: corrected handling of vertical levels when input files do not contain uppermost layers
  • readwind_ecmwf_mpi.f90: ditto

Code enhancements:

  • options/OUTGRID: added comments describing contents
  • options/SPECIES/SPECIES_*: aligned comments
  • options/SPECIES/specoverview.f90: removed commented lines, rectified lines indenting
  • src/FLEXPART.f90: rectified lines indenting, updated date in version string
  • src/FLEXPART_MPI.f90: ditto, and realigned code with src/FLEXPART.f90
  • src/gridcheck_*.f90: added code to write out name of file before it is opened (helps a lot when an input file causes troubles)
  • src/par_mod.f90: added comment explaining relevance of nuvzmax for GRIB input
  • src/readreleases.f90: write out warning if too few particles are used to randomize release
  • src/readspecies.f90: write out name of SPECIES file before it is read
  • src/readwind_*.f90: write out name of input file before opening it
  • src/writeheader_txt.f90: removed wrong comment
  • Property mode set to 100644
File size: 18.5 KB
Line 
1! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
2! SPDX-License-Identifier: GPL-3.0-or-later
3
4subroutine readwind_nests(indj,n,uuhn,vvhn,wwhn)
5  !                           i   i  o    o    o
6  !*****************************************************************************
7  !                                                                            *
8  !     This routine reads the wind fields for the nested model domains.       *
9  !     It is similar to subroutine readwind, which reads the mother domain.   *
10  !                                                                            *
11  !     Authors: A. Stohl, G. Wotawa                                           *
12  !                                                                            *
13  !     8 February 1999                                                        *
14  !                                                                            *
15  !     Last update: 17 October 2000, A. Stohl                                 *
16  !                                                                            *
17  !*****************************************************************************
18  !  Changes, Bernd C. Krueger, Feb. 2001:                                     *
19  !        Variables tthn and qvhn (on eta coordinates) in common block        *
20  !  CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with ECMWF grib_api    *
21  !  CHANGE: 03/12/2008, Harald Sodemann, update to f90 with ECMWF grib_api    *
22  !*****************************************************************************
23
24  use grib_api
25  use par_mod
26  use com_mod
27
28  implicit none
29
30  !HSO  parameters for grib_api
31  integer :: ifile
32  integer :: iret
33  integer :: igrib
34  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
35  integer :: parId !!added by mc for making it consistent with new readwind.f90
36  integer :: gotGrid
37  !HSO  end
38
39  real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
40  real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
41  real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
42  integer :: indj,i,j,k,n,levdiff2,ifield,iumax,iwmax,l
43
44  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
45
46  ! dimension of isec2 at least (22+n), where n is the number of parallels or
47  ! meridians in a quasi-regular (reduced) Gaussian or lat/long grid
48
49  ! dimension of zsec2 at least (10+nn), where nn is the number of vertical
50  ! coordinate parameters
51
52  integer :: isec1(56),isec2(22+nxmaxn+nymaxn)
53  real(kind=4) :: zsec4(jpunp)
54  real(kind=4) :: xaux,yaux
55  real(kind=8) :: xauxin,yauxin
56  real,parameter :: eps=1.e-4
57  real :: ewss(0:nxmaxn-1,0:nymaxn-1),nsss(0:nxmaxn-1,0:nymaxn-1)
58  real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
59  real :: conversion_factor !added by mc to make it consistent with new gridchek.f90
60
61  logical :: hflswitch,strswitch
62
63  !HSO  grib api error messages
64  character(len=24) :: gribErrorMsg = 'Error reading grib file'
65  character(len=20) :: gribFunction = 'readwind_nests'
66
67  do l=1,numbnests
68    hflswitch=.false.
69    strswitch=.false.
70    levdiff2=nlev_ec-nwz+1
71    iumax=0
72    iwmax=0
73
74    ifile=0
75    igrib=0
76    iret=0
77
78  !
79  ! OPENING OF DATA FILE (GRIB CODE)
80  !
81    write(*,*)   'Reading: '//path(numpath+2*(l-1)+1) &
82         (1:length(numpath+2*(l-1)+1))//trim(wfnamen(l,indj))
835   call grib_open_file(ifile,path(numpath+2*(l-1)+1) &
84         (1:length(numpath+2*(l-1)+1))//trim(wfnamen(l,indj)),'r')
85  if (iret.ne.GRIB_SUCCESS) then
86    goto 888   ! ERROR DETECTED
87  endif
88  !turn on support for multi fields messages */
89  !call grib_multi_support_on
90
91    gotGrid=0
92    ifield=0
9310   ifield=ifield+1
94  !
95  ! GET NEXT FIELDS
96  !
97  call grib_new_from_file(ifile,igrib,iret)
98  if (iret.eq.GRIB_END_OF_FILE)  then
99    goto 50    ! EOF DETECTED
100  elseif (iret.ne.GRIB_SUCCESS) then
101    goto 888   ! ERROR DETECTED
102  endif
103
104  !first see if we read GRIB1 or GRIB2
105  call grib_get_int(igrib,'editionNumber',gribVer,iret)
106  call grib_check(iret,gribFunction,gribErrorMsg)
107
108  if (gribVer.eq.1) then ! GRIB Edition 1
109
110  !print*,'GRiB Edition 1'
111  !read the grib2 identifiers
112  call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
113  call grib_check(iret,gribFunction,gribErrorMsg)
114  call grib_get_int(igrib,'level',isec1(8),iret)
115  call grib_check(iret,gribFunction,gribErrorMsg)
116
117  !change code for etadot to code for omega
118  if (isec1(6).eq.77) then
119    isec1(6)=135
120  endif
121
122  conversion_factor=1.
123
124
125  else
126
127  !print*,'GRiB Edition 2'
128  !read the grib2 identifiers
129  call grib_get_int(igrib,'discipline',discipl,iret)
130  call grib_check(iret,gribFunction,gribErrorMsg)
131  call grib_get_int(igrib,'parameterCategory',parCat,iret)
132  call grib_check(iret,gribFunction,gribErrorMsg)
133  call grib_get_int(igrib,'parameterNumber',parNum,iret)
134  call grib_check(iret,gribFunction,gribErrorMsg)
135  call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
136  call grib_check(iret,gribFunction,gribErrorMsg)
137  call grib_get_int(igrib,'level',valSurf,iret)
138  call grib_check(iret,gribFunction,gribErrorMsg)
139  call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consisitent with new readwind.f90
140  call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consisitent with new readwind.f90
141
142  !print*,discipl,parCat,parNum,typSurf,valSurf
143
144  !convert to grib1 identifiers
145  isec1(6)=-1
146  isec1(7)=-1
147  isec1(8)=-1
148  isec1(8)=valSurf     ! level
149   conversion_factor=1.
150  if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
151    isec1(6)=130         ! indicatorOfParameter
152  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
153    isec1(6)=131         ! indicatorOfParameter
154  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
155    isec1(6)=132         ! indicatorOfParameter
156  elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
157    isec1(6)=133         ! indicatorOfParameter
158! ESO Cloud water is in a) fields CLWC and CIWC, *or* b) field QC
159    elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
160      isec1(6)=246         ! indicatorOfParameter
161    elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
162      isec1(6)=247         ! indicatorOfParameter
163! ESO qc(=clwc+ciwc):
164    elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
165      isec1(6)=201031         ! indicatorOfParameter
166  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
167    isec1(6)=134         ! indicatorOfParameter
168  elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot !
169    isec1(6)=135         ! indicatorOfParameter
170  elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added by mc to make it consisitent with new readwind.f90
171    isec1(6)=135         ! indicatorOfParameter    !added by mc to make it consisitent with new readwind.f90
172  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
173    isec1(6)=151         ! indicatorOfParameter
174  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
175    isec1(6)=165         ! indicatorOfParameter
176  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
177    isec1(6)=166         ! indicatorOfParameter
178  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
179    isec1(6)=167         ! indicatorOfParameter
180  elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
181    isec1(6)=168         ! indicatorOfParameter
182  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
183    isec1(6)=141         ! indicatorOfParameter
184    conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90
185  elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC !added by mc to make it consisitent with new readwind.f90
186    isec1(6)=164         ! indicatorOfParameter
187  elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP !added by mc to make it consisitent with new readwind.f90
188    isec1(6)=142         ! indicatorOfParameter
189  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
190    isec1(6)=143         ! indicatorOfParameter
191    conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90
192  elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
193    isec1(6)=146         ! indicatorOfParameter
194  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
195    isec1(6)=176         ! indicatorOfParameter
196  elseif ((parCat.eq.2).and.(parNum.eq.38) .or. parId .eq. 180) then ! EWSS !added by mc to make it consisitent with new readwind.f90
197    isec1(6)=180         ! indicatorOfParameter
198  elseif ((parCat.eq.2).and.(parNum.eq.37) .or. parId .eq. 181) then ! NSSS !added by mc to make it consisitent with new readwind.f90
199    isec1(6)=181         ! indicatorOfParameter
200  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
201    isec1(6)=129         ! indicatorOfParameter
202   elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO !added by mc to make it consisitent with new readwind.f90
203    isec1(6)=160         ! indicatorOfParameter
204  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
205       (typSurf.eq.1)) then ! LSM
206    isec1(6)=172         ! indicatorOfParameter
207  elseif (parNum.eq.152) then
208      isec1(6)=152         ! avoid warning for lnsp       
209  else
210    print*,'***WARNING: undefined GRiB2 message found!',discipl, &
211         parCat,parNum,typSurf
212  endif
213  if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consisitent with new readwind.f90
214    write(*,*) 'parId',parId, 'isec1(6)',isec1(6)  !
215!    stop
216  endif
217
218  endif
219
220  !HSO  get the size and data of the values array
221  if (isec1(6).ne.-1) then
222    call grib_get_real4_array(igrib,'values',zsec4,iret)
223    call grib_check(iret,gribFunction,gribErrorMsg)
224  endif
225
226  !HSO  get the required fields from section 2 in a gribex compatible manner
227  if(ifield.eq.1) then
228  call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
229       isec2(2),iret)
230  call grib_check(iret,gribFunction,gribErrorMsg)
231  call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
232       isec2(3),iret)
233  call grib_check(iret,gribFunction,gribErrorMsg)
234  call grib_get_int(igrib,'numberOfVerticalCoordinateValues', &
235       isec2(12))
236  call grib_check(iret,gribFunction,gribErrorMsg)
237  ! CHECK GRID SPECIFICATIONS
238  if(isec2(2).ne.nxn(l)) stop &
239  'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'
240  if(isec2(3).ne.nyn(l)) stop &
241  'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'
242  if(isec2(12)/2-1.ne.nlev_ec) stop 'READWIND: VERTICAL DISCRET&
243       &IZATION NOT CONSISTENT FOR A NESTING LEVEL'
244  endif ! ifield
245
246  !HSO  get the second part of the grid dimensions only from GRiB1 messages
247 if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then ! !added by mc to make it consisitent with new readwind.f90
248    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
249         xauxin,iret)
250    call grib_check(iret,gribFunction,gribErrorMsg)
251    call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
252         yauxin,iret)
253    call grib_check(iret,gribFunction,gribErrorMsg)
254    if (xauxin.gt.180.) xauxin=xauxin-360.0
255    if (xauxin.lt.-180.) xauxin=xauxin+360.0
256
257    xaux=xauxin
258    yaux=yauxin
259    if (abs(xaux-xlon0n(l)).gt.eps) &
260    stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT FOR A NESTING LEVEL'
261    if (abs(yaux-ylat0n(l)).gt.eps) &
262    stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT FOR A NESTING LEVEL'
263    gotGrid=1
264  endif
265
266    do j=0,nyn(l)-1
267      do i=0,nxn(l)-1
268        k=isec1(8)
269        if(isec1(6).eq.130) tthn(i,j,nlev_ec-k+2,n,l)= &!! TEMPERATURE
270             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
271        if(isec1(6).eq.131) uuhn(i,j,nlev_ec-k+2,l)= &!! U VELOCITY
272             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
273        if(isec1(6).eq.132) vvhn(i,j,nlev_ec-k+2,l)= &!! V VELOCITY
274             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
275        if(isec1(6).eq.133) then                         !! SPEC. HUMIDITY
276          qvhn(i,j,nlev_ec-k+2,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
277          if (qvhn(i,j,nlev_ec-k+2,n,l) .lt. 0.) &
278               qvhn(i,j,nlev_ec-k+2,n,l) = 0.
279  !          this is necessary because the gridded data may contain
280  !          spurious negative values
281        endif
282        if(isec1(6).eq.134) psn(i,j,1,n,l)= &!! SURF. PRESS.
283             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
284        if(isec1(6).eq.135) wwhn(i,j,nlev_ec-k+1,l)= &!! W VELOCITY
285             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
286        if(isec1(6).eq.141) sdn(i,j,1,n,l)= &!! SNOW DEPTH
287             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90!
288        if(isec1(6).eq.151) msln(i,j,1,n,l)= &!! SEA LEVEL PRESS.
289             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
290        if(isec1(6).eq.164) tccn(i,j,1,n,l)= &!! CLOUD COVER
291             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
292        if(isec1(6).eq.165) u10n(i,j,1,n,l)= &!! 10 M U VELOCITY
293             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
294        if(isec1(6).eq.166) v10n(i,j,1,n,l)= &!! 10 M V VELOCITY
295             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
296        if(isec1(6).eq.167) tt2n(i,j,1,n,l)= &!! 2 M TEMPERATURE
297             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
298        if(isec1(6).eq.168) td2n(i,j,1,n,l)= &!! 2 M DEW POINT
299             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
300        if(isec1(6).eq.142) then                         !! LARGE SCALE PREC.
301          lsprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
302          if (lsprecn(i,j,1,n,l).lt.0.) lsprecn(i,j,1,n,l)=0.
303        endif
304        if(isec1(6).eq.143) then                         !! CONVECTIVE PREC.
305          convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90
306          if (convprecn(i,j,1,n,l).lt.0.) convprecn(i,j,1,n,l)=0.
307        endif
308        if(isec1(6).eq.146) sshfn(i,j,1,n,l)= &!! SENS. HEAT FLUX
309             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
310        if((isec1(6).eq.146).and. &
311             (zsec4(nxn(l)*(nyn(l)-j-1)+i+1).ne.0.)) hflswitch=.true.    ! Heat flux available
312        if(isec1(6).eq.176) then                         !! SOLAR RADIATION
313          ssrn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
314          if (ssrn(i,j,1,n,l).lt.0.) ssrn(i,j,1,n,l)=0.
315        endif
316        if(isec1(6).eq.180) ewss(i,j)= &!! EW SURFACE STRESS
317             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
318        if(isec1(6).eq.181) nsss(i,j)= &!! NS SURFACE STRESS
319             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
320        if(((isec1(6).eq.180).or.(isec1(6).eq.181)).and. &
321             (zsec4(nxn(l)*(nyn(l)-j-1)+i+1).ne.0.)) strswitch=.true.    ! stress available
322        if(isec1(6).eq.129) oron(i,j,l)= &!! ECMWF OROGRAPHY
323             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
324        if(isec1(6).eq.160) excessoron(i,j,l)= &!! STANDARD DEVIATION OF OROGRAPHY
325             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
326        if(isec1(6).eq.172) lsmn(i,j,l)= &!! ECMWF LAND SEA MASK
327             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
328        if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
329        if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
330
331! ESO TODO:
332! -add check for if one of clwc/ciwc missing (error),
333!    also if all 3 cw fields present, use qc and disregard the others
334        if(isec1(6).eq.246) then  !! CLWC  Cloud liquid water content [kg/kg]
335          clwchn(i,j,nlev_ec-k+2,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
336          readclouds_nest(l)=.true.
337          sumclouds_nest(l)=.false.
338        endif
339        if(isec1(6).eq.247) then  !! CIWC  Cloud ice water content
340          ciwchn(i,j,nlev_ec-k+2,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
341        endif
342!ZHG end
343!ESO read qc (=clwc+ciwc)
344        if(isec1(6).eq.201031) then  !! QC  Cloud liquid water content [kg/kg]
345          clwchn(i,j,nlev_ec-k+2,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
346          readclouds_nest(l)=.true.
347          sumclouds_nest(l)=.true.
348        endif
349
350
351      end do
352    end do
353
354  call grib_release(igrib)
355  goto 10                      !! READ NEXT LEVEL OR PARAMETER
356  !
357  ! CLOSING OF INPUT DATA FILE
358  !
35950   call grib_close_file(ifile)
360
361  !error message if no fields found with correct first longitude in it
362  if (gotGrid.eq.0) then
363    print*,'***ERROR: input file needs to contain GRiB1 formatted'// &
364         'messages'
365    stop
366  endif
367
368!  if(levdiff2.eq.0) then
369    iwmax=iwmax+1
370    do i=0,nxn(l)-1
371      do j=0,nyn(l)-1
372        wwhn(i,j,iwmax,l)=0.
373      end do
374    end do
375!  endif
376
377  do i=0,nxn(l)-1
378    do j=0,nyn(l)-1
379      surfstrn(i,j,1,n,l)=sqrt(ewss(i,j)**2+nsss(i,j)**2)
380    end do
381  end do
382
383  if ((.not.hflswitch).or.(.not.strswitch)) then
384    write(*,*) 'WARNING: No flux data contained in GRIB file ', &
385         wfnamen(l,indj)
386
387  ! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
388  ! As ECMWF has increased the model resolution, such that now the first model
389  ! level is at about 10 m (where 10-m wind is given), use the 2nd ECMWF level
390  ! (3rd model level in FLEXPART) for the profile method
391  !***************************************************************************
392
393    do i=0,nxn(l)-1
394      do j=0,nyn(l)-1
395        plev1=akz(3)+bkz(3)*psn(i,j,1,n,l)
396        pmean=0.5*(psn(i,j,1,n,l)+plev1)
397        tv=tthn(i,j,3,n,l)*(1.+0.61*qvhn(i,j,3,n,l))
398        fu=-r_air*tv/ga/pmean
399        hlev1=fu*(plev1-psn(i,j,1,n,l))   ! HEIGTH OF FIRST MODEL LAYER
400        ff10m= sqrt(u10n(i,j,1,n,l)**2+v10n(i,j,1,n,l)**2)
401        fflev1=sqrt(uuhn(i,j,3,l)**2+vvhn(i,j,3,l)**2)
402        call pbl_profile(psn(i,j,1,n,l),td2n(i,j,1,n,l),hlev1, &
403             tt2n(i,j,1,n,l),tthn(i,j,3,n,l),ff10m,fflev1, &
404             surfstrn(i,j,1,n,l),sshfn(i,j,1,n,l))
405        if(sshfn(i,j,1,n,l).gt.200.) sshfn(i,j,1,n,l)=200.
406        if(sshfn(i,j,1,n,l).lt.-400.) sshfn(i,j,1,n,l)=-400.
407      end do
408    end do
409  endif
410
411
412  ! Assign 10 m wind to model level at eta=1.0 to have one additional model
413  ! level at the ground
414  ! Specific humidity is taken the same as at one level above
415  ! Temperature is taken as 2 m temperature
416  !**************************************************************************
417
418    do i=0,nxn(l)-1
419      do j=0,nyn(l)-1
420        uuhn(i,j,1,l)=u10n(i,j,1,n,l)
421        vvhn(i,j,1,l)=v10n(i,j,1,n,l)
422        qvhn(i,j,1,n,l)=qvhn(i,j,2,n,l)
423        tthn(i,j,1,n,l)=tt2n(i,j,1,n,l)
424      end do
425    end do
426
427    if(iumax.ne.nuvz-1) stop &
428         'READWIND: NUVZ NOT CONSISTENT FOR A NESTING LEVEL'
429    if(iwmax.ne.nwz) stop &
430         'READWIND: NWZ NOT CONSISTENT FOR A NESTING LEVEL'
431
432  end do
433
434  return
435888   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
436  write(*,*) ' #### ',wfnamen(l,indj),' FOR NESTING LEVEL  #### '
437  write(*,*) ' #### ',l,' IS NOT GRIB FORMAT !!!           #### '
438  stop 'Execution terminated'
439
440
441999   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
442  write(*,*) ' #### ',wfnamen(l,indj),'                    #### '
443  write(*,*) ' #### CANNOT BE OPENED FOR NESTING LEVEL ',l,'####'
444
445end subroutine readwind_nests
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG