source: flexpart.git/src/gridcheck_ecmwf.f90 @ 3481cc1

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file since 3481cc1 was 3481cc1, checked in by Ignacio Pisso <ip@…>, 4 years ago

move license from headers to a different file
  • Property mode set to 100644
File size: 20.5 KB
Line 
1subroutine gridcheck_ecmwf
2
3  !**********************************************************************
4  !                                                                     *
5  !             FLEXPART MODEL SUBROUTINE GRIDCHECK                     *
6  !                                                                     *
7  !**********************************************************************
8  !                                                                     *
9  !             AUTHOR:      G. WOTAWA                                  *
10  !             DATE:        1997-08-06                                 *
11  !             LAST UPDATE: 1997-10-10                                 *
12  !                                                                     *
13  !             Update:      1999-02-08, global fields allowed, A. Stohl*
14  !             CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
15  !                                 ECMWF grib_api                      *
16  !             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
17  !                                 ECMWF grib_api                      *
18  !                                                                     *
19  !   Unified ECMWF and GFS builds                                      *
20  !   Marian Harustak, 12.5.2017                                        *
21  !     - Renamed from gridcheck to gridcheck_ecmwf                     *
22  !                                                                     *
23  !**********************************************************************
24  !                                                                     *
25  ! DESCRIPTION:                                                        *
26  !                                                                     *
27  ! THIS SUBROUTINE DETERMINES THE GRID SPECIFICATIONS (LOWER LEFT      *
28  ! LONGITUDE, LOWER LEFT LATITUDE, NUMBER OF GRID POINTS, GRID DIST-   *
29  ! ANCE AND VERTICAL DISCRETIZATION OF THE ECMWF MODEL) FROM THE       *
30  ! GRIB HEADER OF THE FIRST INPUT FILE. THE CONSISTANCY (NO CHANGES    *
31  ! WITHIN ONE FLEXPART RUN) IS CHECKED IN THE ROUTINE "READWIND" AT    *
32  ! ANY CALL.                                                           *
33  !                                                                     *
34  ! XLON0                geographical longitude of lower left gridpoint *
35  ! YLAT0                geographical latitude of lower left gridpoint  *
36  ! NX                   number of grid points x-direction              *
37  ! NY                   number of grid points y-direction              *
38  ! DX                   grid distance x-direction                      *
39  ! DY                   grid distance y-direction                      *
40  ! NUVZ                 number of grid points for horizontal wind      *
41  !                      components in z direction                      *
42  ! NWZ                  number of grid points for vertical wind        *
43  !                      component in z direction                       *
44  ! sizesouth, sizenorth give the map scale (i.e. number of virtual grid*
45  !                      points of the polar stereographic grid):       *
46  !                      used to check the CFL criterion                *
47  ! UVHEIGHT(1)-         heights of gridpoints where u and v are        *
48  ! UVHEIGHT(NUVZ)       given                                          *
49  ! WHEIGHT(1)-          heights of gridpoints where w is given         *
50  ! WHEIGHT(NWZ)                                                        *
51  !                                                                     *
52  !**********************************************************************
53
54  use grib_api
55  use par_mod
56  use com_mod
57  use conv_mod
58  use cmapf_mod, only: stlmbr,stcm2p
59
60  implicit none
61
62  !HSO  parameters for grib_api
63  integer :: ifile
64  integer :: iret
65  integer :: igrib
66  integer :: gotGrid
67  real(kind=4) :: xaux1,xaux2,yaux1,yaux2
68  real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
69  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parId
70  !HSO  end
71  integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
72  real :: sizesouth,sizenorth,xauxa,pint,conversion_factor
73
74  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
75
76  ! dimension of isec2 at least (22+n), where n is the number of parallels or
77  ! meridians in a quasi-regular (reduced) Gaussian or lat/long grid
78
79  ! dimension of zsec2 at least (10+nn), where nn is the number of vertical
80  ! coordinate parameters
81
82  integer :: isec1(56),isec2(22+nxmax+nymax)
83  real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
84  character(len=1) :: opt
85
86  !HSO  grib api error messages
87  character(len=24) :: gribErrorMsg = 'Error reading grib file'
88  character(len=20) :: gribFunction = 'gridcheck'
89
90
91  iumax=0
92  iwmax=0
93
94  if(ideltas.gt.0) then
95    ifn=1
96  else
97    ifn=numbwf
98  endif
99  !
100  ! OPENING OF DATA FILE (GRIB CODE)
101  !
1025 call grib_open_file(ifile,path(3)(1:length(3)) &
103       //trim(wfname(ifn)),'r',iret)
104  if (iret.ne.GRIB_SUCCESS) then
105    goto 999   ! ERROR DETECTED
106  endif
107  !turn on support for multi fields messages
108  !call grib_multi_support_on
109
110  gotGrid=0
111  ifield=0
11210 ifield=ifield+1
113
114  !
115  ! GET NEXT FIELDS
116  !
117  call grib_new_from_file(ifile,igrib,iret)
118  if (iret.eq.GRIB_END_OF_FILE )  then
119    goto 30    ! EOF DETECTED
120  elseif (iret.ne.GRIB_SUCCESS) then
121    goto 999   ! ERROR DETECTED
122  endif
123
124  !first see if we read GRIB1 or GRIB2
125  call grib_get_int(igrib,'editionNumber',gribVer,iret)
126  call grib_check(iret,gribFunction,gribErrorMsg)
127
128  if (gribVer.eq.1) then ! GRIB Edition 1
129
130    !print*,'GRiB Edition 1'
131    !read the grib2 identifiers
132    call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
133    call grib_check(iret,gribFunction,gribErrorMsg)
134    call grib_get_int(igrib,'level',isec1(8),iret)
135    call grib_check(iret,gribFunction,gribErrorMsg)
136
137    !change code for etadot to code for omega
138    if (isec1(6).eq.77) then
139      isec1(6)=135
140    endif
141
142    !print*,isec1(6),isec1(8)
143
144  else
145
146    !print*,'GRiB Edition 2'
147    !read the grib2 identifiers
148    call grib_get_int(igrib,'discipline',discipl,iret)
149    call grib_check(iret,gribFunction,gribErrorMsg)
150    call grib_get_int(igrib,'parameterCategory',parCat,iret)
151    call grib_check(iret,gribFunction,gribErrorMsg)
152    call grib_get_int(igrib,'parameterNumber',parNum,iret)
153    call grib_check(iret,gribFunction,gribErrorMsg)
154    call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
155    call grib_check(iret,gribFunction,gribErrorMsg)
156    call grib_get_int(igrib,'level',valSurf,iret)
157    call grib_check(iret,gribFunction,gribErrorMsg)
158    call grib_get_int(igrib,'paramId',parId,iret)
159    call grib_check(iret,gribFunction,gribErrorMsg)
160
161    !print*,discipl,parCat,parNum,typSurf,valSurf
162
163    !convert to grib1 identifiers
164    isec1(6)=-1
165    isec1(7)=-1
166    isec1(8)=-1
167    isec1(8)=valSurf     ! level
168    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
169      isec1(6)=130         ! indicatorOfParameter
170    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
171      isec1(6)=131         ! indicatorOfParameter
172    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
173      isec1(6)=132         ! indicatorOfParameter
174    elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
175      isec1(6)=133         ! indicatorOfParameter
176      !ZHG FOR CLOUDS FROM GRIB
177    elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
178      isec1(6)=246         ! indicatorOfParameter
179    elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
180      isec1(6)=247         ! indicatorOfParameter
181      !ZHG end
182      ! ESO qc(=clwc+ciwc)
183    elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
184      isec1(6)=201031      ! indicatorOfParameter
185    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
186      isec1(6)=134         ! indicatorOfParameter
187    elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
188      isec1(6)=135         ! indicatorOfParameter
189    elseif ((parCat.eq.128).and.(parNum.eq.77)) 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) .or. parId .eq. 164) then ! CC
204      isec1(6)=164         ! indicatorOfParameter
205    elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) 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) .or. parId .eq. 180) then ! EWSS
214      isec1(6)=180         ! indicatorOfParameter
215    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) 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) .or. parId .eq. 160) 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    if(parId .ne. isec1(6) .and. parId .ne. 77) then
229      write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
230      !    stop
231    endif
232
233  endif
234
235  CALL grib_get_int(igrib,'numberOfPointsAlongAParallel', &
236       isec2(2),iret)
237  ! nx=isec2(2)
238  ! WRITE(*,*) nx,nxmax
239  IF (isec2(2).GT.nxmax) THEN
240    WRITE(*,*) 'FLEXPART error: Too many grid points in x direction.'
241    WRITE(*,*) 'Reduce resolution of wind fields.'
242    WRITE(*,*) 'Or change parameter settings in file ecmwf_mod.'
243    WRITE(*,*) nx,nxmax
244!    STOP
245  ENDIF
246
247  !get the size and data of the values array
248  if (isec1(6).ne.-1) then
249    call grib_get_real4_array(igrib,'values',zsec4,iret)
250    call grib_check(iret,gribFunction,gribErrorMsg)
251  endif
252
253  if (ifield.eq.1) then
254
255    !HSO  get the required fields from section 2 in a gribex compatible manner
256    call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
257         isec2(2),iret)
258    call grib_check(iret,gribFunction,gribErrorMsg)
259    call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
260         isec2(3),iret)
261    call grib_check(iret,gribFunction,gribErrorMsg)
262    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
263         xaux1in,iret)
264    call grib_check(iret,gribFunction,gribErrorMsg)
265    call grib_get_int(igrib,'numberOfVerticalCoordinateValues', &
266         isec2(12),iret)
267    call grib_check(iret,gribFunction,gribErrorMsg)
268
269    !  get the size and data of the vertical coordinate array
270    call grib_get_real4_array(igrib,'pv',zsec2,iret)
271    call grib_check(iret,gribFunction,gribErrorMsg)
272
273    nxfield=isec2(2)
274    ny=isec2(3)
275    nlev_ec=isec2(12)/2-1
276  endif
277
278  !HSO  get the second part of the grid dimensions only from GRiB1 messages
279  if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
280    call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
281         xaux2in,iret)
282    call grib_check(iret,gribFunction,gribErrorMsg)
283    call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
284         yaux1in,iret)
285    call grib_check(iret,gribFunction,gribErrorMsg)
286    call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees', &
287         yaux2in,iret)
288    call grib_check(iret,gribFunction,gribErrorMsg)
289    xaux1=xaux1in
290    xaux2=xaux2in
291    yaux1=yaux1in
292    yaux2=yaux2in
293    if (xaux1.gt.180.) xaux1=xaux1-360.0
294    if (xaux2.gt.180.) xaux2=xaux2-360.0
295    if (xaux1.lt.-180.) xaux1=xaux1+360.0
296    if (xaux2.lt.-180.) xaux2=xaux2+360.0
297    if (xaux2.lt.xaux1) xaux2=xaux2+360.0
298    xlon0=xaux1
299    ylat0=yaux1
300    dx=(xaux2-xaux1)/real(nxfield-1)
301    dy=(yaux2-yaux1)/real(ny-1)
302    dxconst=180./(dx*r_earth*pi)
303    dyconst=180./(dy*r_earth*pi)
304    gotGrid=1
305    ! Check whether fields are global
306    ! If they contain the poles, specify polar stereographic map
307    ! projections using the stlmbr- and stcm2p-calls
308    !***********************************************************
309
310    xauxa=abs(xaux2+dx-360.-xaux1)
311    if (xauxa.lt.0.001) then
312      nx=nxfield+1                 ! field is cyclic
313      xglobal=.true.
314      if (abs(nxshift).ge.nx) &
315           stop 'nxshift in file par_mod is too large'
316      xlon0=xlon0+real(nxshift)*dx
317    else
318      nx=nxfield
319      xglobal=.false.
320      if (nxshift.ne.0) &
321           stop 'nxshift (par_mod) must be zero for non-global domain'
322    endif
323    nxmin1=nx-1
324    nymin1=ny-1
325    if (xlon0.gt.180.) xlon0=xlon0-360.
326    xauxa=abs(yaux1+90.)
327    if (xglobal.and.xauxa.lt.0.001) then
328      sglobal=.true.               ! field contains south pole
329      ! Enhance the map scale by factor 3 (*2=6) compared to north-south
330      ! map scale
331      sizesouth=6.*(switchsouth+90.)/dy
332      call stlmbr(southpolemap,-90.,0.)
333      call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
334           sizesouth,switchsouth,180.)
335      switchsouthg=(switchsouth-ylat0)/dy
336    else
337      sglobal=.false.
338      switchsouthg=999999.
339    endif
340    xauxa=abs(yaux2-90.)
341    if (xglobal.and.xauxa.lt.0.001) then
342      nglobal=.true.               ! field contains north pole
343      ! Enhance the map scale by factor 3 (*2=6) compared to north-south
344      ! map scale
345      sizenorth=6.*(90.-switchnorth)/dy
346      call stlmbr(northpolemap,90.,0.)
347      call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
348           sizenorth,switchnorth,180.)
349      switchnorthg=(switchnorth-ylat0)/dy
350    else
351      nglobal=.false.
352      switchnorthg=999999.
353    endif
354    if (nxshift.lt.0) &
355         stop 'nxshift (par_mod) must not be negative'
356    if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large'
357  endif ! gotGrid
358
359  if (nx.gt.nxmax) then
360    write(*,*) 'FLEXPART error: Too many grid points in x direction.'
361    write(*,*) 'Reduce resolution of wind fields.'
362    write(*,*) 'Or change parameter settings in file par_mod.'
363    write(*,*) nx,nxmax
364    stop
365  endif
366
367  if (ny.gt.nymax) then
368    write(*,*) 'FLEXPART error: Too many grid points in y direction.'
369    write(*,*) 'Reduce resolution of wind fields.'
370    write(*,*) 'Or change parameter settings in file par_mod.'
371    write(*,*) ny,nymax
372    stop
373  endif
374
375  k=isec1(8)
376  if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
377  if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
378
379  if(isec1(6).eq.129) then
380    do jy=0,ny-1
381      do ix=0,nxfield-1
382        oro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)/ga
383      end do
384    end do
385  endif
386  if(isec1(6).eq.172) then
387    do jy=0,ny-1
388      do ix=0,nxfield-1
389        lsm(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
390      end do
391    end do
392  endif
393  if(isec1(6).eq.160) then
394    do jy=0,ny-1
395      do ix=0,nxfield-1
396        excessoro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
397      end do
398    end do
399  endif
400
401  call grib_release(igrib)
402  goto 10                      !! READ NEXT LEVEL OR PARAMETER
403  !
404  ! CLOSING OF INPUT DATA FILE
405  !
406
40730 call grib_close_file(ifile)
408
409  !error message if no fields found with correct first longitude in it
410  if (gotGrid.eq.0) then
411    print*,'***ERROR: input file needs to contain GRiB1 formatted'// &
412         'messages'
413    stop
414  endif
415
416  nuvz=iumax
417  nwz =iwmax
418  if(nuvz.eq.nlev_ec) nwz=nlev_ec+1
419
420  if (nuvz+1.gt.nuvzmax) then
421    write(*,*) 'FLEXPART error: Too many u,v grid points in z '// &
422         'direction.'
423    write(*,*) 'Reduce resolution of wind fields.'
424    write(*,*) 'Or change parameter settings in file par_mod.'
425    write(*,*) nuvz+1,nuvzmax
426    stop
427  endif
428
429  if (nwz.gt.nwzmax) then
430    write(*,*) 'FLEXPART error: Too many w grid points in z '// &
431         'direction.'
432    write(*,*) 'Reduce resolution of wind fields.'
433    write(*,*) 'Or change parameter settings in file par_mod.'
434    write(*,*) nwz,nwzmax
435    stop
436  endif
437
438  ! If desired, shift all grids by nxshift grid cells
439  !**************************************************
440
441  if (xglobal) then
442    call shift_field_0(oro,nxfield,ny)
443    call shift_field_0(lsm,nxfield,ny)
444    call shift_field_0(excessoro,nxfield,ny)
445  endif
446
447  ! Output of grid info
448  !********************
449
450  if (lroot) then
451    write(*,'(a,2i7)') ' Vertical levels in ECMWF data: ', &
452         nuvz+1,nwz
453    write(*,*)
454    write(*,'(a)') ' Mother domain:'
455    write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Longitude range: ', &
456         xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
457    write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range : ', &
458         ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
459    write(*,*)
460  end if
461
462  ! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL
463  ! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
464
465  numskip=nlev_ec-nuvz  ! number of ecmwf model layers not used
466  ! by trajectory model
467  !do 8940 i=1,244
468  !   write (*,*) 'zsec2:',i,ifield,zsec2(i),numskip
469  !940  continue
470  !   stop
471  ! SEC SEC SEC
472  ! for unknown reason zsec 1 to 10 is filled in this version
473  ! compared to the old one
474  ! SEC SEC SE
475  do i=1,nwz
476    j=numskip+i
477    k=nlev_ec+1+numskip+i
478    akm(nwz-i+1)=zsec2(j)
479    !   write (*,*) 'ifield:',ifield,k,j,zsec2(10+j)
480    bkm(nwz-i+1)=zsec2(k)
481  end do
482
483  !
484  ! CALCULATION OF AKZ, BKZ
485  ! AKZ,BKZ: model discretization parameters at the center of each model
486  !     layer
487  !
488  ! Assign the 10 m winds to an artificial model level with akz=0 and bkz=1.0,
489  ! i.e. ground level
490  !*****************************************************************************
491
492  akz(1)=0.
493  bkz(1)=1.0
494  do i=1,nuvz
495    akz(i+1)=0.5*(akm(i+1)+akm(i))
496    bkz(i+1)=0.5*(bkm(i+1)+bkm(i))
497  end do
498  nuvz=nuvz+1
499
500  ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
501  ! upon the transformation to z levels. In order to save computer memory, this is
502  ! not done anymore in the standard version. However, this option can still be
503  ! switched on by replacing the following lines with those below, that are
504  ! currently commented out. For this, similar changes are necessary in
505  ! verttransform.f and verttranform_nests.f
506  !*****************************************************************************
507
508  nz=nuvz
509  if (nz.gt.nzmax) stop 'nzmax too small'
510  do i=1,nuvz
511    aknew(i)=akz(i)
512    bknew(i)=bkz(i)
513  end do
514
515  ! Switch on following lines to use doubled vertical resolution
516  !*************************************************************
517  !nz=nuvz+nwz-1
518  !if (nz.gt.nzmax) stop 'nzmax too small'
519  !do 100 i=1,nwz
520  !  aknew(2*(i-1)+1)=akm(i)
521  !00     bknew(2*(i-1)+1)=bkm(i)
522  !do 110 i=2,nuvz
523  !  aknew(2*(i-1))=akz(i)
524  !10     bknew(2*(i-1))=bkz(i)
525  ! End doubled vertical resolution
526
527
528  ! Determine the uppermost level for which the convection scheme shall be applied
529  ! by assuming that there is no convection above 50 hPa (for standard SLP)
530  !*****************************************************************************
531
532  do i=1,nuvz-2
533    pint=akz(i)+bkz(i)*101325.
534    if (pint.lt.5000.) goto 96
535  end do
53696 nconvlev=i
537  if (nconvlev.gt.nconvlevmax-1) then
538    nconvlev=nconvlevmax-1
539    write(*,*) 'Attention, convection only calculated up to ', &
540         akz(nconvlev)+bkz(nconvlev)*1013.25,' hPa'
541  endif
542
543  return
544
545999 write(*,*)
546  write(*,*) ' ###########################################'// &
547       '###### '
548  write(*,*) '       TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
549  write(*,*) ' CAN NOT OPEN INPUT DATA FILE '//wfname(ifn)
550  write(*,*) ' ###########################################'// &
551       '###### '
552  write(*,*)
553  write(*,'(a)') '!!! PLEASE INSERT A NEW CD-ROM AND   !!!'
554  write(*,'(a)') '!!! PRESS ANY KEY TO CONTINUE...     !!!'
555  write(*,'(a)') '!!! ...OR TERMINATE FLEXPART PRESSING!!!'
556  write(*,'(a)') '!!! THE "X" KEY...                   !!!'
557  write(*,*)
558  read(*,'(a)') opt
559  if(opt.eq.'X') then
560    stop
561  else
562    goto 5
563  endif
564
565end subroutine gridcheck_ecmwf
566
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG