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