source: trunk/src/gridcheck_fnl.f90 @ 28

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