source: trunk/src/readwind.f90 @ 20

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

move version 9.1.8 form branches to trunk. Contributions from HSO, saeck, pesei, NIK, RT, XKF, IP and others

File size: 19.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(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
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,xaux0,yaux0
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
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  !HSO conversion of ECMWF etadot to etadot*dp/deta
105  logical :: etacon=.false.
106  real,parameter :: p00=101325.
107  real :: dak,dbk
108
109  hflswitch=.false.
110  strswitch=.false.
111  levdiff2=nlev_ec-nwz+1
112  iumax=0
113  iwmax=0
114
115  !
116  ! OPENING OF DATA FILE (GRIB CODE)
117  !
1185   call grib_open_file(ifile,path(3)(1:length(3)) &
119         //trim(wfname(indj)),'r',iret)
120  if (iret.ne.GRIB_SUCCESS) then
121    goto 888   ! ERROR DETECTED
122  endif
123  !turn on support for multi fields messages */
124  !call grib_multi_support_on
125
126  gotGrid=0
127  ifield=0
12810   ifield=ifield+1
129  !
130  ! GET NEXT FIELDS
131  !
132  call grib_new_from_file(ifile,igrib,iret)
133  if (iret.eq.GRIB_END_OF_FILE)  then
134    goto 50    ! EOF DETECTED
135  elseif (iret.ne.GRIB_SUCCESS) then
136    goto 888   ! ERROR DETECTED
137  endif
138
139  !first see if we read GRIB1 or GRIB2
140  call grib_get_int(igrib,'editionNumber',gribVer,iret)
141  call grib_check(iret,gribFunction,gribErrorMsg)
142
143  if (gribVer.eq.1) then ! GRIB Edition 1
144
145  !print*,'GRiB Edition 1'
146  !read the grib2 identifiers
147  call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
148  call grib_check(iret,gribFunction,gribErrorMsg)
149  call grib_get_int(igrib,'level',isec1(8),iret)
150  call grib_check(iret,gribFunction,gribErrorMsg)
151
152  !change code for etadot to code for omega
153  if (isec1(6).eq.77) then
154    isec1(6)=135
155  endif
156
157  else
158
159  !print*,'GRiB Edition 2'
160  !read the grib2 identifiers
161  call grib_get_int(igrib,'discipline',discipl,iret)
162  call grib_check(iret,gribFunction,gribErrorMsg)
163  call grib_get_int(igrib,'parameterCategory',parCat,iret)
164  call grib_check(iret,gribFunction,gribErrorMsg)
165  call grib_get_int(igrib,'parameterNumber',parNum,iret)
166  call grib_check(iret,gribFunction,gribErrorMsg)
167  call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
168  call grib_check(iret,gribFunction,gribErrorMsg)
169  call grib_get_int(igrib,'level',valSurf,iret)
170  call grib_check(iret,gribFunction,gribErrorMsg)
171
172  !print*,discipl,parCat,parNum,typSurf,valSurf
173
174  !convert to grib1 identifiers
175  isec1(6)=-1
176  isec1(7)=-1
177  isec1(8)=-1
178  isec1(8)=valSurf     ! level
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.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
192    isec1(6)=151         ! indicatorOfParameter
193  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
194    isec1(6)=165         ! indicatorOfParameter
195  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
196    isec1(6)=166         ! indicatorOfParameter
197  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
198    isec1(6)=167         ! indicatorOfParameter
199  elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
200    isec1(6)=168         ! indicatorOfParameter
201  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
202    isec1(6)=141         ! indicatorOfParameter
203  elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC
204    isec1(6)=164         ! indicatorOfParameter
205  elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
206    isec1(6)=142         ! indicatorOfParameter
207  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
208    isec1(6)=143         ! indicatorOfParameter
209  elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
210    isec1(6)=146         ! indicatorOfParameter
211  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
212    isec1(6)=176         ! indicatorOfParameter
213  elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
214    isec1(6)=180         ! indicatorOfParameter
215  elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
216    isec1(6)=181         ! indicatorOfParameter
217  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
218    isec1(6)=129         ! indicatorOfParameter
219  elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
220    isec1(6)=160         ! indicatorOfParameter
221  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
222       (typSurf.eq.1)) then ! LSM
223    isec1(6)=172         ! indicatorOfParameter
224  else
225    print*,'***ERROR: undefined GRiB2 message found!',discipl, &
226         parCat,parNum,typSurf
227  endif
228
229  endif
230
231  !HSO  get the size and data of the values array
232  if (isec1(6).ne.-1) then
233    call grib_get_real4_array(igrib,'values',zsec4,iret)
234    call grib_check(iret,gribFunction,gribErrorMsg)
235  endif
236
237  !HSO  get the required fields from section 2 in a gribex compatible manner
238  if (ifield.eq.1) then
239  call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
240       isec2(2),iret)
241  call grib_check(iret,gribFunction,gribErrorMsg)
242  call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
243       isec2(3),iret)
244  call grib_check(iret,gribFunction,gribErrorMsg)
245  call grib_get_int(igrib,'numberOfVerticalCoordinateValues', &
246       isec2(12))
247  call grib_check(iret,gribFunction,gribErrorMsg)
248  ! CHECK GRID SPECIFICATIONS
249  if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
250  if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
251  if(isec2(12)/2-1.ne.nlev_ec) &
252       stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
253  endif ! ifield
254
255  !HSO  get the second part of the grid dimensions only from GRiB1 messages
256  if ((gribVer.eq.1).and.(gotGrid.eq.0)) then
257    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
258         xauxin,iret)
259    call grib_check(iret,gribFunction,gribErrorMsg)
260    call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
261         yauxin,iret)
262    call grib_check(iret,gribFunction,gribErrorMsg)
263    xaux=xauxin+real(nxshift)*dx
264    yaux=yauxin
265    xaux0=xlon0
266    yaux0=ylat0
267    if(xaux.lt.0.) xaux=xaux+360.
268    if(yaux.lt.0.) yaux=yaux+360.
269    if(xaux0.lt.0.) xaux0=xaux0+360.
270    if(yaux0.lt.0.) yaux0=yaux0+360.
271    if(abs(xaux-xaux0).gt.eps) &
272         stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
273    if(abs(yaux-yaux0).gt.eps) &
274         stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
275    gotGrid=1
276  endif ! gotGrid
277
278  do j=0,nymin1
279    do i=0,nxfield-1
280      k=isec1(8)
281      if(isec1(6).eq.130) tth(i,j,nlev_ec-k+2,n)= &!! TEMPERATURE
282           zsec4(nxfield*(ny-j-1)+i+1)
283      if(isec1(6).eq.131) uuh(i,j,nlev_ec-k+2)= &!! U VELOCITY
284           zsec4(nxfield*(ny-j-1)+i+1)
285      if(isec1(6).eq.132) vvh(i,j,nlev_ec-k+2)= &!! V VELOCITY
286           zsec4(nxfield*(ny-j-1)+i+1)
287      if(isec1(6).eq.133) then                      !! SPEC. HUMIDITY
288        qvh(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
289        if (qvh(i,j,nlev_ec-k+2,n) .lt. 0.) &
290             qvh(i,j,nlev_ec-k+2,n) = 0.
291  !        this is necessary because the gridded data may contain
292  !        spurious negative values
293      endif
294      if(isec1(6).eq.134) ps(i,j,1,n)= &!! SURF. PRESS.
295           zsec4(nxfield*(ny-j-1)+i+1)
296
297      if(isec1(6).eq.135) wwh(i,j,nlev_ec-k+1)= &!! W VELOCITY
298           zsec4(nxfield*(ny-j-1)+i+1)
299      if(isec1(6).eq.141) sd(i,j,1,n)= &!! SNOW DEPTH
300           zsec4(nxfield*(ny-j-1)+i+1)
301      if(isec1(6).eq.151) msl(i,j,1,n)= &!! SEA LEVEL PRESS.
302           zsec4(nxfield*(ny-j-1)+i+1)
303      if(isec1(6).eq.164) tcc(i,j,1,n)= &!! CLOUD COVER
304           zsec4(nxfield*(ny-j-1)+i+1)
305      if(isec1(6).eq.165) u10(i,j,1,n)= &!! 10 M U VELOCITY
306           zsec4(nxfield*(ny-j-1)+i+1)
307      if(isec1(6).eq.166) v10(i,j,1,n)= &!! 10 M V VELOCITY
308           zsec4(nxfield*(ny-j-1)+i+1)
309      if(isec1(6).eq.167) tt2(i,j,1,n)= &!! 2 M TEMPERATURE
310           zsec4(nxfield*(ny-j-1)+i+1)
311      if(isec1(6).eq.168) td2(i,j,1,n)= &!! 2 M DEW POINT
312           zsec4(nxfield*(ny-j-1)+i+1)
313      if(isec1(6).eq.142) then                      !! LARGE SCALE PREC.
314        lsprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
315        if (lsprec(i,j,1,n).lt.0.) lsprec(i,j,1,n)=0.
316      endif
317      if(isec1(6).eq.143) then                      !! CONVECTIVE PREC.
318        convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
319        if (convprec(i,j,1,n).lt.0.) convprec(i,j,1,n)=0.
320      endif
321      if(isec1(6).eq.146) sshf(i,j,1,n)= &!! SENS. HEAT FLUX
322           zsec4(nxfield*(ny-j-1)+i+1)
323      if((isec1(6).eq.146).and.(zsec4(nxfield*(ny-j-1)+i+1).ne.0.)) &
324           hflswitch=.true.    ! Heat flux available
325      if(isec1(6).eq.176) then                      !! SOLAR RADIATION
326        ssr(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
327        if (ssr(i,j,1,n).lt.0.) ssr(i,j,1,n)=0.
328      endif
329      if(isec1(6).eq.180) ewss(i,j)= &!! EW SURFACE STRESS
330           zsec4(nxfield*(ny-j-1)+i+1)
331      if(isec1(6).eq.181) nsss(i,j)= &!! NS SURFACE STRESS
332           zsec4(nxfield*(ny-j-1)+i+1)
333      if(((isec1(6).eq.180).or.(isec1(6).eq.181)).and. &
334           (zsec4(nxfield*(ny-j-1)+i+1).ne.0.)) strswitch=.true.    ! stress available
335  !sec        strswitch=.true.
336      if(isec1(6).eq.129) oro(i,j)= &!! ECMWF OROGRAPHY
337           zsec4(nxfield*(ny-j-1)+i+1)/ga
338      if(isec1(6).eq.160) excessoro(i,j)= &!! STANDARD DEVIATION OF OROGRAPHY
339           zsec4(nxfield*(ny-j-1)+i+1)
340      if(isec1(6).eq.172) lsm(i,j)= &!! ECMWF LAND SEA MASK
341           zsec4(nxfield*(ny-j-1)+i+1)
342      if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
343      if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
344
345    end do
346  end do
347
348  call grib_release(igrib)
349  goto 10                      !! READ NEXT LEVEL OR PARAMETER
350  !
351  ! CLOSING OF INPUT DATA FILE
352  !
353
35450   call grib_close_file(ifile)
355
356  !error message if no fields found with correct first longitude in it
357  if (gotGrid.eq.0) then
358    print*,'***ERROR: input file needs to contain GRiB1 formatted'// &
359         'messages'
360    stop
361  endif
362
363  if(levdiff2.eq.0) then
364    iwmax=nlev_ec+1
365    do i=0,nxmin1
366      do j=0,nymin1
367        wwh(i,j,nlev_ec+1)=0.
368      end do
369    end do
370  endif
371
372  ! convert from ECMWF etadot to etadot*dp/deta as needed by FLEXPART
373  if(etacon.eqv..true.) then
374    do k=1,nwzmax
375      dak=akm(k+1)-akm(k)
376      dbk=bkm(k+1)-bkm(k)
377      do i=0,nxmin1
378        do j=0,nymin1
379          wwh(i,j,k)=2*wwh(i,j,k)*ps(i,j,1,n)*(dak/ps(i,j,1,n)+dbk)/(dak/p00+dbk)
380          if (k.gt.1) then
381            wwh(i,j,k)=wwh(i,j,k)-wwh(i,j,k-1)
382          endif
383        end do
384      end do
385    end do
386  endif
387
388  ! For global fields, assign the leftmost data column also to the rightmost
389  ! data column; if required, shift whole grid by nxshift grid points
390  !*************************************************************************
391
392  if (xglobal) then
393    call shift_field_0(ewss,nxfield,ny)
394    call shift_field_0(nsss,nxfield,ny)
395    call shift_field_0(oro,nxfield,ny)
396    call shift_field_0(excessoro,nxfield,ny)
397    call shift_field_0(lsm,nxfield,ny)
398    call shift_field(ps,nxfield,ny,1,1,2,n)
399    call shift_field(sd,nxfield,ny,1,1,2,n)
400    call shift_field(msl,nxfield,ny,1,1,2,n)
401    call shift_field(tcc,nxfield,ny,1,1,2,n)
402    call shift_field(u10,nxfield,ny,1,1,2,n)
403    call shift_field(v10,nxfield,ny,1,1,2,n)
404    call shift_field(tt2,nxfield,ny,1,1,2,n)
405    call shift_field(td2,nxfield,ny,1,1,2,n)
406    call shift_field(lsprec,nxfield,ny,1,1,2,n)
407    call shift_field(convprec,nxfield,ny,1,1,2,n)
408    call shift_field(sshf,nxfield,ny,1,1,2,n)
409    call shift_field(ssr,nxfield,ny,1,1,2,n)
410    call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n)
411    call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n)
412    call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
413    call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
414    call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
415  endif
416
417  do i=0,nxmin1
418    do j=0,nymin1
419      surfstr(i,j,1,n)=sqrt(ewss(i,j)**2+nsss(i,j)**2)
420    end do
421  end do
422
423  if ((.not.hflswitch).or.(.not.strswitch)) then
424    write(*,*) 'WARNING: No flux data contained in GRIB file ', &
425         wfname(indj)
426
427  ! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
428  ! As ECMWF has increased the model resolution, such that now the first model
429  ! level is at about 10 m (where 10-m wind is given), use the 2nd ECMWF level
430  ! (3rd model level in FLEXPART) for the profile method
431  !***************************************************************************
432
433    do i=0,nxmin1
434      do j=0,nymin1
435        plev1=akz(3)+bkz(3)*ps(i,j,1,n)
436        pmean=0.5*(ps(i,j,1,n)+plev1)
437        tv=tth(i,j,3,n)*(1.+0.61*qvh(i,j,3,n))
438        fu=-r_air*tv/ga/pmean
439        hlev1=fu*(plev1-ps(i,j,1,n))   ! HEIGTH OF FIRST MODEL LAYER
440        ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2)
441        fflev1=sqrt(uuh(i,j,3)**2+vvh(i,j,3)**2)
442        call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
443             tt2(i,j,1,n),tth(i,j,3,n),ff10m,fflev1, &
444             surfstr(i,j,1,n),sshf(i,j,1,n))
445        if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
446        if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
447      end do
448    end do
449  endif
450
451
452  ! Assign 10 m wind to model level at eta=1.0 to have one additional model
453  ! level at the ground
454  ! Specific humidity is taken the same as at one level above
455  ! Temperature is taken as 2 m temperature
456  !**************************************************************************
457
458     do i=0,nxmin1
459        do j=0,nymin1
460           uuh(i,j,1)=u10(i,j,1,n)
461           vvh(i,j,1)=v10(i,j,1,n)
462           qvh(i,j,1,n)=qvh(i,j,2,n)
463           tth(i,j,1,n)=tt2(i,j,1,n)
464        end do
465     end do
466
467  if(iumax.ne.nuvz-1) stop 'READWIND: NUVZ NOT CONSISTENT'
468  if(iwmax.ne.nwz)    stop 'READWIND: NWZ NOT CONSISTENT'
469
470  return
471888   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
472  write(*,*) ' #### ',wfname(indj),'                    #### '
473  write(*,*) ' #### IS NOT GRIB FORMAT !!!                  #### '
474  stop 'Execution terminated'
475999   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
476  write(*,*) ' #### ',wfname(indj),'                    #### '
477  write(*,*) ' #### CANNOT BE OPENED !!!                    #### '
478  stop 'Execution terminated'
479
480end subroutine readwind
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG