source: trunk/src/readwind_nests.f90 @ 28

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

version 9.2 beta. Changes from HH, AST, MC, NIK, IP. Changes in vert transform. New SPECIES input includes scavenging coefficients

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