source: flexpart.git/flexpart_code/gridcheck.F90 @ 496c607

FPv9.3.1FPv9.3.1b_testingFPv9.3.2fp9.3.1-20161214-nc4grib2nc4_repair
Last change on this file since 496c607 was 496c607, checked in by Don Morton <Don.Morton@…>, 8 years ago

Initial commit of FPv9.3.1

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