source: flexpart.git/src/gridcheck_gfs.f90

bugfixes+enhancements
Last change on this file was dba4221, checked in by Pirmin Kaufmann <pirmin.kaufmann@…>, 18 months ago

Bugfixes:

  • options/SPECIES/SPECIES_009: corrected wrong number format, replaced comma by decimal point
  • options/SPECIES/SPECIES_028: corrected wrong number format, moved sign of exponent to after the E
  • options/SPECIES/specoverview.f90: added namelist parameters that appear in SPECIES files but were missing here
  • src/FLEXPART.f90: replaced compiler-specific command line argument routines by standard Fortran intrinsic routines
  • src/FLEXPART_MPI.f90: ditto
  • src/gridcheck_ecmwf.f90: corrected handling of vertical levels when input files do not contain uppermost layers
  • src/gridcheck_nests.f90: ditto
  • src/readwind_ecmwf.f90: corrected handling of vertical levels when input files do not contain uppermost layers
  • readwind_ecmwf_mpi.f90: ditto

Code enhancements:

  • options/OUTGRID: added comments describing contents
  • options/SPECIES/SPECIES_*: aligned comments
  • options/SPECIES/specoverview.f90: removed commented lines, rectified lines indenting
  • src/FLEXPART.f90: rectified lines indenting, updated date in version string
  • src/FLEXPART_MPI.f90: ditto, and realigned code with src/FLEXPART.f90
  • src/gridcheck_*.f90: added code to write out name of file before it is opened (helps a lot when an input file causes troubles)
  • src/par_mod.f90: added comment explaining relevance of nuvzmax for GRIB input
  • src/readreleases.f90: write out warning if too few particles are used to randomize release
  • src/readspecies.f90: write out name of SPECIES file before it is read
  • src/readwind_*.f90: write out name of input file before opening it
  • src/writeheader_txt.f90: removed wrong comment
  • Property mode set to 100644
File size: 17.7 KB
Line 
1! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
2! SPDX-License-Identifier: GPL-3.0-or-later
3
4subroutine gridcheck_gfs
5
6  !**********************************************************************
7  !                                                                     *
8  !             FLEXPART MODEL SUBROUTINE GRIDCHECK                     *
9  !                                                                     *
10  !**********************************************************************
11  !                                                                     *
12  !             AUTHOR:      G. WOTAWA                                  *
13  !             DATE:        1997-08-06                                 *
14  !             LAST UPDATE: 1997-10-10                                 *
15  !                                                                     *
16  !             Update:      1999-02-08, global fields allowed, A. Stohl*
17  !             CHANGE: 17/11/2005, Caroline Forster, GFS data          *
18  !             CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
19  !                                 ECMWF grib_api                      *
20  !             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
21  !                                 ECMWF grib_api                      *
22  !                                                                     *
23  !   Unified ECMWF and GFS builds                                      *
24  !   Marian Harustak, 12.5.2017                                        *
25  !     - Renamed routine from gridcheck to gridcheck_gfs               *
26  !                                                                     *
27  !**********************************************************************
28  !                                                                     *
29  ! DESCRIPTION:                                                        *
30  !                                                                     *
31  ! THIS SUBROUTINE DETERMINES THE GRID SPECIFICATIONS (LOWER LEFT      *
32  ! LONGITUDE, LOWER LEFT LATITUDE, NUMBER OF GRID POINTS, GRID DIST-   *
33  ! ANCE AND VERTICAL DISCRETIZATION OF THE ECMWF MODEL) FROM THE       *
34  ! GRIB HEADER OF THE FIRST INPUT FILE. THE CONSISTANCY (NO CHANGES    *
35  ! WITHIN ONE FLEXPART RUN) IS CHECKED IN THE ROUTINE "READWIND" AT    *
36  ! ANY CALL.                                                           *
37  !                                                                     *
38  ! XLON0                geographical longitude of lower left gridpoint *
39  ! YLAT0                geographical latitude of lower left gridpoint  *
40  ! NX                   number of grid points x-direction              *
41  ! NY                   number of grid points y-direction              *
42  ! DX                   grid distance x-direction                      *
43  ! DY                   grid distance y-direction                      *
44  ! NUVZ                 number of grid points for horizontal wind      *
45  !                      components in z direction                      *
46  ! NWZ                  number of grid points for vertical wind        *
47  !                      component in z direction                       *
48  ! sizesouth, sizenorth give the map scale (i.e. number of virtual grid*
49  !                      points of the polar stereographic grid):       *
50  !                      used to check the CFL criterion                *
51  ! UVHEIGHT(1)-         heights of gridpoints where u and v are        *
52  ! UVHEIGHT(NUVZ)       given                                          *
53  ! WHEIGHT(1)-          heights of gridpoints where w is given         *
54  ! WHEIGHT(NWZ)                                                        *
55  !                                                                     *
56  !**********************************************************************
57
58  use grib_api
59  use par_mod
60  use com_mod
61  use conv_mod
62  use cmapf_mod, only: stlmbr,stcm2p
63
64  implicit none
65
66  !HSO  parameters for grib_api
67  integer :: ifile
68  integer :: iret
69  integer :: igrib
70  real(kind=4) :: xaux1,xaux2,yaux1,yaux2
71  real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
72  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
73  !HSO  end
74  integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
75  real :: sizesouth,sizenorth,xauxa,pint
76  real :: akm_usort(nwzmax)
77  real,parameter :: eps=0.0001
78
79  ! NCEP GFS
80  real :: pres(nwzmax), help
81
82  integer :: i179,i180,i181
83
84  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
85
86  integer :: isec1(8),isec2(3)
87  real(kind=4) :: zsec4(jpunp)
88  character(len=1) :: opt
89
90  !HSO  grib api error messages
91  character(len=24) :: gribErrorMsg = 'Error reading grib file'
92  character(len=20) :: gribFunction = 'gridcheckwind_gfs'
93  !
94  if (numbnests.ge.1) then
95  write(*,*) ' ###########################################'
96  write(*,*) ' FLEXPART ERROR SUBROUTINE GRIDCHECK:'
97  write(*,*) ' NO NESTED WINDFIELDAS ALLOWED FOR GFS!      '
98  write(*,*) ' ###########################################'
99  stop
100  endif
101
102  iumax=0
103  iwmax=0
104
105  if(ideltas.gt.0) then
106    ifn=1
107  else
108    ifn=numbwf
109  endif
110  !
111  ! OPENING OF DATA FILE (GRIB CODE)
112  !
113  write(*,*)     'Reading: '//path(3)(1:length(3)) &
114         //trim(wfname(ifn))
1155   call grib_open_file(ifile,path(3)(1:length(3)) &
116         //trim(wfname(ifn)),'r',iret)
117  if (iret.ne.GRIB_SUCCESS) then
118    goto 999   ! ERROR DETECTED
119  endif
120  !turn on support for multi fields messages
121  call grib_multi_support_on
122
123  ifield=0
12410   ifield=ifield+1
125  !
126  ! GET NEXT FIELDS
127  !
128  call grib_new_from_file(ifile,igrib,iret)
129  if (iret.eq.GRIB_END_OF_FILE )  then
130    goto 30    ! EOF DETECTED
131  elseif (iret.ne.GRIB_SUCCESS) then
132    goto 999   ! ERROR DETECTED
133  endif
134
135  !first see if we read GRIB1 or GRIB2
136  call grib_get_int(igrib,'editionNumber',gribVer,iret)
137  call grib_check(iret,gribFunction,gribErrorMsg)
138
139  if (gribVer.eq.1) then ! GRIB Edition 1
140
141  !read the grib1 identifiers
142  call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
143  call grib_check(iret,gribFunction,gribErrorMsg)
144  call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
145  call grib_check(iret,gribFunction,gribErrorMsg)
146  call grib_get_int(igrib,'level',isec1(8),iret)
147  call grib_check(iret,gribFunction,gribErrorMsg)
148
149  !get the size and data of the values array
150  call grib_get_real4_array(igrib,'values',zsec4,iret)
151  call grib_check(iret,gribFunction,gribErrorMsg)
152
153  else ! GRIB Edition 2
154
155  !read the grib2 identifiers
156  call grib_get_int(igrib,'discipline',discipl,iret)
157  call grib_check(iret,gribFunction,gribErrorMsg)
158  call grib_get_int(igrib,'parameterCategory',parCat,iret)
159  call grib_check(iret,gribFunction,gribErrorMsg)
160  call grib_get_int(igrib,'parameterNumber',parNum,iret)
161  call grib_check(iret,gribFunction,gribErrorMsg)
162  call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
163  call grib_check(iret,gribFunction,gribErrorMsg)
164  call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
165       valSurf,iret)
166  call grib_check(iret,gribFunction,gribErrorMsg)
167
168  !convert to grib1 identifiers
169  isec1(6)=-1
170  isec1(7)=-1
171  isec1(8)=-1
172  if ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
173    isec1(6)=33          ! indicatorOfParameter
174    isec1(7)=100         ! indicatorOfTypeOfLevel
175    isec1(8)=valSurf/100 ! level, convert to hPa
176  elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
177    isec1(6)=7           ! indicatorOfParameter
178    isec1(7)=1           ! indicatorOfTypeOfLevel
179    isec1(8)=0
180  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
181       .and.(discipl.eq.2)) then ! LSM
182    isec1(6)=81          ! indicatorOfParameter
183    isec1(7)=1           ! indicatorOfTypeOfLevel
184    isec1(8)=0
185  endif
186
187  if (isec1(6).ne.-1) then
188  !  get the size and data of the values array
189    call grib_get_real4_array(igrib,'values',zsec4,iret)
190    call grib_check(iret,gribFunction,gribErrorMsg)
191  endif
192
193  endif ! gribVer
194
195  if(ifield.eq.1) then
196
197  !get the required fields from section 2
198  !store compatible to gribex input
199  call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
200       isec2(2),iret)
201  call grib_check(iret,gribFunction,gribErrorMsg)
202  call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
203       isec2(3),iret)
204  call grib_check(iret,gribFunction,gribErrorMsg)
205  call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
206       xaux1in,iret)
207  call grib_check(iret,gribFunction,gribErrorMsg)
208  call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
209       xaux2in,iret)
210  call grib_check(iret,gribFunction,gribErrorMsg)
211  call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
212       yaux1in,iret)
213  call grib_check(iret,gribFunction,gribErrorMsg)
214  call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees', &
215       yaux2in,iret)
216  call grib_check(iret,gribFunction,gribErrorMsg)
217
218  ! Fix for flexpart.eu ticket #48
219  if (xaux2in.lt.0) xaux2in = 359.0
220
221  xaux1=xaux1in
222  xaux2=xaux2in
223  yaux1=yaux1in
224  yaux2=yaux2in
225
226  nxfield=isec2(2)
227  ny=isec2(3)
228  if((abs(xaux1).lt.eps).and.(xaux2.ge.359)) then ! NCEP DATA FROM 0 TO
229    xaux1=-179.0                             ! 359 DEG EAST ->
230    xaux2=-179.0+360.-360./real(nxfield)    ! TRANSFORMED TO -179
231  endif                                      ! TO 180 DEG EAST
232  if (xaux1.gt.180) xaux1=xaux1-360.0
233  if (xaux2.gt.180) xaux2=xaux2-360.0
234  if (xaux1.lt.-180) xaux1=xaux1+360.0
235  if (xaux2.lt.-180) xaux2=xaux2+360.0
236  if (xaux2.lt.xaux1) xaux2=xaux2+360.
237  xlon0=xaux1
238  ylat0=yaux1
239  dx=(xaux2-xaux1)/real(nxfield-1)
240  dy=(yaux2-yaux1)/real(ny-1)
241  dxconst=180./(dx*r_earth*pi)
242  dyconst=180./(dy*r_earth*pi)
243  !HSO end edits
244
245
246  ! Check whether fields are global
247  ! If they contain the poles, specify polar stereographic map
248  ! projections using the stlmbr- and stcm2p-calls
249  !***********************************************************
250
251    xauxa=abs(xaux2+dx-360.-xaux1)
252    if (xauxa.lt.0.001) then
253      nx=nxfield+1                 ! field is cyclic
254      xglobal=.true.
255      if (abs(nxshift).ge.nx) &
256           stop 'nxshift in file par_mod is too large'
257      xlon0=xlon0+real(nxshift)*dx
258    else
259      nx=nxfield
260      xglobal=.false.
261      if (nxshift.ne.0) &
262           stop 'nxshift (par_mod) must be zero for non-global domain'
263    endif
264    nxmin1=nx-1
265    nymin1=ny-1
266    if (xlon0.gt.180.) xlon0=xlon0-360.
267    xauxa=abs(yaux1+90.)
268    if (xglobal.and.xauxa.lt.0.001) then
269      sglobal=.true.               ! field contains south pole
270  ! Enhance the map scale by factor 3 (*2=6) compared to north-south
271  ! map scale
272      sizesouth=6.*(switchsouth+90.)/dy
273      call stlmbr(southpolemap,-90.,0.)
274      call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
275           sizesouth,switchsouth,180.)
276      switchsouthg=(switchsouth-ylat0)/dy
277    else
278      sglobal=.false.
279      switchsouthg=999999.
280    endif
281    xauxa=abs(yaux2-90.)
282    if (xglobal.and.xauxa.lt.0.001) then
283      nglobal=.true.               ! field contains north pole
284  ! Enhance the map scale by factor 3 (*2=6) compared to north-south
285  ! map scale
286      sizenorth=6.*(90.-switchnorth)/dy
287      call stlmbr(northpolemap,90.,0.)
288      call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
289           sizenorth,switchnorth,180.)
290      switchnorthg=(switchnorth-ylat0)/dy
291    else
292      nglobal=.false.
293      switchnorthg=999999.
294    endif
295  endif ! ifield.eq.1
296
297  if (nxshift.lt.0) stop 'nxshift (par_mod) must not be negative'
298  if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large'
299
300  ! NCEP ISOBARIC LEVELS
301  !*********************
302
303  if((isec1(6).eq.33).and.(isec1(7).eq.100)) then ! check for U wind
304    iumax=iumax+1
305    pres(iumax)=real(isec1(8))*100.0
306  endif
307
308
309  i179=nint(179./dx)
310  if (dx.lt.0.7) then
311    i180=nint(180./dx)+1    ! 0.5 deg data
312  else
313    i180=nint(179./dx)+1    ! 1 deg data
314  endif
315  i181=i180+1
316
317
318  ! NCEP TERRAIN
319  !*************
320
321  if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
322    do jy=0,ny-1
323      do ix=0,nxfield-1
324        help=zsec4(nxfield*(ny-jy-1)+ix+1)
325        if(ix.le.i180) then
326          oro(i179+ix,jy)=help
327          excessoro(i179+ix,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
328        else
329          oro(ix-i181,jy)=help
330          excessoro(ix-i181,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
331        endif
332      end do
333    end do
334  endif
335
336  ! NCEP LAND SEA MASK
337  !*******************
338
339  if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
340    do jy=0,ny-1
341      do ix=0,nxfield-1
342        help=zsec4(nxfield*(ny-jy-1)+ix+1)
343        if(ix.le.i180) then
344          lsm(i179+ix,jy)=help
345        else
346          lsm(ix-i181,jy)=help
347        endif
348      end do
349    end do
350  endif
351
352  call grib_release(igrib)
353
354  goto 10                      !! READ NEXT LEVEL OR PARAMETER
355  !
356  ! CLOSING OF INPUT DATA FILE
357  !
358
359  ! HSO
36030   continue
361  call grib_close_file(ifile)
362  ! HSO end edits
363
364  nuvz=iumax
365  nwz =iumax
366  nlev_ec=iumax
367
368  if (nx.gt.nxmax) then
369   write(*,*) 'FLEXPART error: Too many grid points in x direction.'
370    write(*,*) 'Reduce resolution of wind fields.'
371    write(*,*) 'Or change parameter settings in file par_mod.'
372    write(*,*) nx,nxmax
373    stop
374  endif
375
376  if (ny.gt.nymax) then
377   write(*,*) 'FLEXPART error: Too many grid points in y direction.'
378    write(*,*) 'Reduce resolution of wind fields.'
379    write(*,*) 'Or change parameter settings in file par_mod.'
380    write(*,*) ny,nymax
381    stop
382  endif
383
384  if (nuvz.gt.nuvzmax) then
385    write(*,*) 'FLEXPART error: Too many u,v grid points in z '// &
386         'direction.'
387    write(*,*) 'Reduce resolution of wind fields.'
388    write(*,*) 'Or change parameter settings in file par_mod.'
389    write(*,*) nuvz,nuvzmax
390    stop
391  endif
392
393  if (nwz.gt.nwzmax) then
394    write(*,*) 'FLEXPART error: Too many w grid points in z '// &
395         'direction.'
396    write(*,*) 'Reduce resolution of wind fields.'
397    write(*,*) 'Or change parameter settings in file par_mod.'
398    write(*,*) nwz,nwzmax
399    stop
400  endif
401
402  ! If desired, shift all grids by nxshift grid cells
403  !**************************************************
404
405  if (xglobal) then
406    call shift_field_0(oro,nxfield,ny)
407    call shift_field_0(lsm,nxfield,ny)
408    call shift_field_0(excessoro,nxfield,ny)
409  endif
410
411  ! Output of grid info
412  !********************
413
414  if (lroot) then
415    write(*,*)
416    write(*,*)
417    write(*,'(a,2i7)') 'Vertical levels in NCEP data: ', &
418         nuvz,nwz
419    write(*,*)
420    write(*,'(a)') 'Mother domain:'
421    write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
422         xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
423    write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range : ', &
424         ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
425    write(*,*)
426  end if
427
428  ! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL
429  ! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
430
431  numskip=nlev_ec-nuvz  ! number of ecmwf model layers not used
432                        ! by trajectory model
433  do i=1,nwz
434    j=numskip+i
435    k=nlev_ec+1+numskip+i
436    akm_usort(nwz-i+1)=pres(nwz-i+1)
437    bkm(nwz-i+1)=0.0
438  end do
439
440  !******************************
441  ! change Sabine Eckhardt: akm should always be in descending order ... readwind adapted!
442  !******************************
443      do i=1,nwz
444         if (akm_usort(1).gt.akm_usort(2)) then
445            akm(i)=akm_usort(i)
446         else
447            akm(i)=akm_usort(nwz-i+1)
448         endif
449      end do
450
451  !
452  ! CALCULATION OF AKZ, BKZ
453  ! AKZ,BKZ: model discretization parameters at the center of each model
454  !     layer
455  !
456  ! Assign the 10 m winds to an artificial model level with akz=0 and bkz=1.0,
457  ! i.e. ground level
458  !*****************************************************************************
459
460  do i=1,nuvz
461     akz(i)=akm(i)
462     bkz(i)=bkm(i)
463  end do
464
465  ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
466  ! upon the transformation to z levels. In order to save computer memory, this is
467  ! not done anymore in the standard version. However, this option can still be
468  ! switched on by replacing the following lines with those below, that are
469  ! currently commented out. For this, similar changes are necessary in
470  ! verttransform.f and verttranform_nests.f
471  !*****************************************************************************
472
473  nz=nuvz
474  if (nz.gt.nzmax) stop 'nzmax too small'
475  do i=1,nuvz
476    aknew(i)=akz(i)
477    bknew(i)=bkz(i)
478  end do
479
480  ! Switch on following lines to use doubled vertical resolution
481  !*************************************************************
482  !nz=nuvz+nwz-1
483  !if (nz.gt.nzmax) stop 'nzmax too small'
484  !do 100 i=1,nwz
485  !  aknew(2*(i-1)+1)=akm(i)
486  !00     bknew(2*(i-1)+1)=bkm(i)
487  !do 110 i=2,nuvz
488  !  aknew(2*(i-1))=akz(i)
489  !10     bknew(2*(i-1))=bkz(i)
490  ! End doubled vertical resolution
491
492
493  ! Determine the uppermost level for which the convection scheme shall be applied
494  ! by assuming that there is no convection above 50 hPa (for standard SLP)
495  !*****************************************************************************
496
497  do i=1,nuvz-2
498    pint=akz(i)+bkz(i)*101325.
499    if (pint.lt.5000.) goto 96
500  end do
50196   nconvlev=i
502  if (nconvlev.gt.nconvlevmax-1) then
503    nconvlev=nconvlevmax-1
504    write(*,*) 'Attention, convection only calculated up to ', &
505         akz(nconvlev)+bkz(nconvlev)*1013.25,' hPa'
506  endif
507
508  return
509
510999   write(*,*)
511  write(*,*) ' ###########################################'// &
512       '###### '
513  write(*,*) '       TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
514  write(*,*) ' CAN NOT OPEN INPUT DATA FILE '//wfname(ifn)
515  write(*,*) ' ###########################################'// &
516       '###### '
517  write(*,*)
518  write(*,'(a)') '!!! PLEASE INSERT A NEW CD-ROM AND   !!!'
519  write(*,'(a)') '!!! PRESS ANY KEY TO CONTINUE...     !!!'
520  write(*,'(a)') '!!! ...OR TERMINATE FLEXPART PRESSING!!!'
521  write(*,'(a)') '!!! THE "X" KEY...                   !!!'
522  write(*,*)
523  read(*,'(a)') opt
524  if(opt.eq.'X') then
525    stop
526  else
527    goto 5
528  endif
529
530end subroutine gridcheck_gfs
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG