source: trunk/src/readwind.f90 @ 24

Last change on this file since 24 was 24, checked in by igpis, 8 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