source: flexpart.git/src/gridcheck_gfs.f90 @ e2132cf

dev
Last change on this file since e2132cf was e2132cf, checked in by Ignacio Pisso <ip@…>, 3 years ago

fix permissions
  • Property mode set to 100644
File size: 17.6 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 eccodes
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,xsec18
76  real :: akm_usort(nwzmax)
77  real,parameter :: eps=spacing(2.0_4*360.0_4)
78
79  ! NCEP GFS
80  real :: pres(nwzmax), help
81
82  integer :: i180
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  !
1135   call grib_open_file(ifile,path(3)(1:length(3)) &
114         //trim(wfname(ifn)),'r',iret)
115  if (iret.ne.GRIB_SUCCESS) then
116    goto 999   ! ERROR DETECTED
117  endif
118  !turn on support for multi fields messages
119  call grib_multi_support_on
120
121  ifield=0
12210   ifield=ifield+1
123  !
124  ! GET NEXT FIELDS
125  !
126  call grib_new_from_file(ifile,igrib,iret)
127  if (iret.eq.GRIB_END_OF_FILE )  then
128    goto 30    ! EOF DETECTED
129  elseif (iret.ne.GRIB_SUCCESS) then
130    goto 999   ! ERROR DETECTED
131  endif
132
133  !first see if we read GRIB1 or GRIB2
134  call grib_get_int(igrib,'editionNumber',gribVer,iret)
135  call grib_check(iret,gribFunction,gribErrorMsg)
136
137  if (gribVer.eq.1) then ! GRIB Edition 1
138
139  !read the grib1 identifiers
140  call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
141  call grib_check(iret,gribFunction,gribErrorMsg)
142  call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),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  !get the size and data of the values array
148  call grib_get_real4_array(igrib,'values',zsec4,iret)
149  call grib_check(iret,gribFunction,gribErrorMsg)
150
151  xsec18 = real(isec1(8))
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  xsec18 = -1.0
173  if ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
174    isec1(6)=33          ! indicatorOfParameter
175    isec1(7)=100         ! indicatorOfTypeOfLevel
176    xsec18=valSurf/100.0 ! level, convert to hPa
177  elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
178    isec1(6)=7           ! indicatorOfParameter
179    isec1(7)=1           ! indicatorOfTypeOfLevel
180    xsec18=real(0)
181  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
182       .and.(discipl.eq.2)) then ! LSM
183    isec1(6)=81          ! indicatorOfParameter
184    isec1(7)=1           ! indicatorOfTypeOfLevel
185    xsec18=real(0)
186  endif
187
188  if (isec1(6).ne.-1) then
189  !  get the size and data of the values array
190    call grib_get_real4_array(igrib,'values',zsec4,iret)
191    call grib_check(iret,gribFunction,gribErrorMsg)
192  endif
193
194  endif ! gribVer
195
196  if(ifield.eq.1) then
197
198  !get the required fields from section 2
199  !store compatible to gribex input
200  call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
201       isec2(2),iret)
202  call grib_check(iret,gribFunction,gribErrorMsg)
203  call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
204       isec2(3),iret)
205  call grib_check(iret,gribFunction,gribErrorMsg)
206  call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
207       xaux1in,iret)
208  call grib_check(iret,gribFunction,gribErrorMsg)
209  call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
210       xaux2in,iret)
211  call grib_check(iret,gribFunction,gribErrorMsg)
212  call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
213       yaux1in,iret)
214  call grib_check(iret,gribFunction,gribErrorMsg)
215  call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees', &
216       yaux2in,iret)
217  call grib_check(iret,gribFunction,gribErrorMsg)
218
219  ! Fix for flexpart.eu ticket #48
220  if (xaux2in.lt.0) xaux2in = 359.0
221
222  xaux1=xaux1in
223  xaux2=xaux2in
224  yaux1=yaux1in
225  yaux2=yaux2in
226
227  nxfield=isec2(2)
228  ny=isec2(3)
229  if((abs(xaux1).lt.eps).and.(xaux2.ge.359.)) then ! NCEP DATA FROM 0 TO
230    xaux1=-180.0                             ! 359 DEG EAST ->
231    xaux2=-180.0+360.-360./real(nxfield)    ! TRANSFORMED TO -179
232  endif                                      ! TO 180 DEG EAST
233  if (xaux1.gt.180) xaux1=xaux1-360.0
234  if (xaux2.gt.180) xaux2=xaux2-360.0
235  if (xaux1.lt.-180) xaux1=xaux1+360.0
236  if (xaux2.lt.-180) xaux2=xaux2+360.0
237  if (xaux2.lt.xaux1) xaux2=xaux2+360.
238  xlon0=xaux1
239  ylat0=yaux1
240  dx=(xaux2-xaux1)/real(nxfield-1)
241  dy=(yaux2-yaux1)/real(ny-1)
242  dxconst=180./(dx*r_earth*pi)
243  dyconst=180./(dy*r_earth*pi)
244  !HSO end edits
245
246
247  ! Check whether fields are global
248  ! If they contain the poles, specify polar stereographic map
249  ! projections using the stlmbr- and stcm2p-calls
250  !***********************************************************
251
252    xauxa=abs(xaux2+dx-360.-xaux1)
253    if (xauxa.lt.0.001) then
254      nx=nxfield+1                 ! field is cyclic
255      xglobal=.true.
256      if (abs(nxshift).ge.nx) &
257           stop 'nxshift in file par_mod is too large'
258      xlon0=xlon0+real(nxshift)*dx
259    else
260      nx=nxfield
261      xglobal=.false.
262      if (nxshift.ne.0) &
263           stop 'nxshift (par_mod) must be zero for non-global domain'
264    endif
265    nxmin1=nx-1
266    nymin1=ny-1
267    if (xlon0.gt.180.) xlon0=xlon0-360.
268    xauxa=abs(yaux1+90.)
269    if (xglobal.and.xauxa.lt.0.001) then
270      sglobal=.true.               ! field contains south pole
271  ! Enhance the map scale by factor 3 (*2=6) compared to north-south
272  ! map scale
273      sizesouth=6.*(switchsouth+90.)/dy
274      call stlmbr(southpolemap,-90.,0.)
275      call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
276           sizesouth,switchsouth,180.)
277      switchsouthg=(switchsouth-ylat0)/dy
278    else
279      sglobal=.false.
280      switchsouthg=999999.
281    endif
282    xauxa=abs(yaux2-90.)
283    if (xglobal.and.xauxa.lt.0.001) then
284      nglobal=.true.               ! field contains north pole
285  ! Enhance the map scale by factor 3 (*2=6) compared to north-south
286  ! map scale
287      sizenorth=6.*(90.-switchnorth)/dy
288      call stlmbr(northpolemap,90.,0.)
289      call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
290           sizenorth,switchnorth,180.)
291      switchnorthg=(switchnorth-ylat0)/dy
292    else
293      nglobal=.false.
294      switchnorthg=999999.
295    endif
296  endif ! ifield.eq.1
297
298  if (nxshift.lt.0) stop 'nxshift (par_mod) must not be negative'
299  if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large'
300
301  ! NCEP ISOBARIC LEVELS
302  !*********************
303
304  if((isec1(6).eq.33).and.(isec1(7).eq.100)) then ! check for U wind
305    iumax=iumax+1
306    pres(iumax)=xsec18*100.0
307  endif
308
309
310  i180=nint(180./dx)    ! 0.5 deg data
311
312
313  ! NCEP TERRAIN
314  !*************
315
316  if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
317    do jy=0,ny-1
318      do ix=0,nxfield-1
319        help=zsec4(nxfield*(ny-jy-1)+ix+1)
320        if(ix.lt.i180) then
321          oro(i180+ix,jy)=help
322          excessoro(i180+ix,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
323        else
324          oro(ix-i180,jy)=help
325          excessoro(ix-i180,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
326        endif
327      end do
328    end do
329  endif
330
331  ! NCEP LAND SEA MASK
332  !*******************
333
334  if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
335    do jy=0,ny-1
336      do ix=0,nxfield-1
337        help=zsec4(nxfield*(ny-jy-1)+ix+1)
338        if(ix.lt.i180) then
339          lsm(i180+ix,jy)=help
340        else
341          lsm(ix-i180,jy)=help
342        endif
343      end do
344    end do
345  endif
346
347  call grib_release(igrib)
348
349  goto 10                      !! READ NEXT LEVEL OR PARAMETER
350  !
351  ! CLOSING OF INPUT DATA FILE
352  !
353
354  ! HSO
35530   continue
356  call grib_close_file(ifile)
357  ! HSO end edits
358
359  nuvz=iumax
360  nwz =iumax
361  nlev_ec=iumax
362
363  if (nx.gt.nxmax) then
364   write(*,*) 'FLEXPART error: Too many grid points in x direction.'
365    write(*,*) 'Reduce resolution of wind fields.'
366    write(*,*) 'Or change parameter settings in file par_mod.'
367    write(*,*) nx,nxmax
368    stop
369  endif
370
371  if (ny.gt.nymax) then
372   write(*,*) 'FLEXPART error: Too many grid points in y direction.'
373    write(*,*) 'Reduce resolution of wind fields.'
374    write(*,*) 'Or change parameter settings in file par_mod.'
375    write(*,*) ny,nymax
376    stop
377  endif
378
379  if (nuvz.gt.nuvzmax) then
380    write(*,*) 'FLEXPART error: Too many u,v grid points in z '// &
381         'direction.'
382    write(*,*) 'Reduce resolution of wind fields.'
383    write(*,*) 'Or change parameter settings in file par_mod.'
384    write(*,*) nuvz,nuvzmax
385    stop
386  endif
387
388  if (nwz.gt.nwzmax) then
389    write(*,*) 'FLEXPART error: Too many w grid points in z '// &
390         'direction.'
391    write(*,*) 'Reduce resolution of wind fields.'
392    write(*,*) 'Or change parameter settings in file par_mod.'
393    write(*,*) nwz,nwzmax
394    stop
395  endif
396
397  ! If desired, shift all grids by nxshift grid cells
398  !**************************************************
399
400  if (xglobal) then
401    call shift_field_0(oro,nxfield,ny)
402    call shift_field_0(lsm,nxfield,ny)
403    call shift_field_0(excessoro,nxfield,ny)
404  endif
405
406  ! Output of grid info
407  !********************
408
409  if (lroot) then
410    write(*,*)
411    write(*,*)
412  write(*,'(a,2i7)') '# of vertical levels in NCEP data: ', &
413         nuvz,nwz
414    write(*,*)
415    write(*,'(a)') 'Mother domain:'
416    write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
417         xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
418    write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range : ', &
419         ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
420    write(*,*)
421  end if
422
423  ! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL
424  ! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
425
426  numskip=nlev_ec-nuvz  ! number of ecmwf model layers not used
427                        ! by trajectory model
428  do i=1,nwz
429    j=numskip+i
430    k=nlev_ec+1+numskip+i
431    akm_usort(nwz-i+1)=pres(nwz-i+1)
432    bkm(nwz-i+1)=0.0
433  end do
434
435  !******************************
436  ! change Sabine Eckhardt: akm should always be in descending order ... readwind adapted!
437  !******************************
438      do i=1,nwz
439         if (akm_usort(1).gt.akm_usort(2)) then
440            akm(i)=akm_usort(i)
441         else
442            akm(i)=akm_usort(nwz-i+1)
443         endif
444      end do
445
446  !
447  ! CALCULATION OF AKZ, BKZ
448  ! AKZ,BKZ: model discretization parameters at the center of each model
449  !     layer
450  !
451  ! Assign the 10 m winds to an artificial model level with akz=0 and bkz=1.0,
452  ! i.e. ground level
453  !*****************************************************************************
454
455  do i=1,nuvz
456     akz(i)=akm(i)
457     bkz(i)=bkm(i)
458  end do
459
460  ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
461  ! upon the transformation to z levels. In order to save computer memory, this is
462  ! not done anymore in the standard version. However, this option can still be
463  ! switched on by replacing the following lines with those below, that are
464  ! currently commented out. For this, similar changes are necessary in
465  ! verttransform.f and verttranform_nests.f
466  !*****************************************************************************
467
468  nz=nuvz
469  if (nz.gt.nzmax) stop 'nzmax too small'
470  do i=1,nuvz
471    aknew(i)=akz(i)
472    bknew(i)=bkz(i)
473  end do
474
475  ! Switch on following lines to use doubled vertical resolution
476  !*************************************************************
477  !nz=nuvz+nwz-1
478  !if (nz.gt.nzmax) stop 'nzmax too small'
479  !do 100 i=1,nwz
480  !  aknew(2*(i-1)+1)=akm(i)
481  !00     bknew(2*(i-1)+1)=bkm(i)
482  !do 110 i=2,nuvz
483  !  aknew(2*(i-1))=akz(i)
484  !10     bknew(2*(i-1))=bkz(i)
485  ! End doubled vertical resolution
486
487
488  ! Determine the uppermost level for which the convection scheme shall be applied
489  ! by assuming that there is no convection above 50 hPa (for standard SLP)
490  !*****************************************************************************
491
492  do i=1,nuvz-2
493    pint=akz(i)+bkz(i)*101325.
494    if (pint.lt.5000.) goto 96
495  end do
49696   nconvlev=i
497  if (nconvlev.gt.nconvlevmax-1) then
498    nconvlev=nconvlevmax-1
499    write(*,*) 'Attention, convection only calculated up to ', &
500         akz(nconvlev)+bkz(nconvlev)*1013.25,' hPa'
501  endif
502
503  return
504
505999   write(*,*)
506  write(*,*) ' ###########################################'// &
507       '###### '
508  write(*,*) '       TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
509  write(*,*) ' CAN NOT OPEN INPUT DATA FILE '//wfname(ifn)
510  write(*,*) ' ###########################################'// &
511       '###### '
512  write(*,*)
513  write(*,'(a)') '!!! PLEASE INSERT A NEW CD-ROM AND   !!!'
514  write(*,'(a)') '!!! PRESS ANY KEY TO CONTINUE...     !!!'
515  write(*,'(a)') '!!! ...OR TERMINATE FLEXPART PRESSING!!!'
516  write(*,'(a)') '!!! THE "X" KEY...                   !!!'
517  write(*,*)
518  read(*,'(a)') opt
519  if(opt.eq.'X') then
520    stop
521  else
522    goto 5
523  endif
524
525end subroutine gridcheck_gfs
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG