source: flexpart.git/src/gridcheck_nests.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@…>, 5 years ago

move license from headers to a different file
  • Property mode set to 100644
File size: 16.0 KB
RevLine 
[e200b7a]1subroutine gridcheck_nests
2
3  !*****************************************************************************
4  !                                                                            *
5  !     This routine checks the grid specification for the nested model        *
6  !     domains. It is similar to subroutine gridcheck, which checks the       *
7  !     mother domain.                                                         *
8  !                                                                            *
9  !     Authors: A. Stohl, G. Wotawa                                           *
10  !                                                                            *
11  !     8 February 1999                                                        *
12  !                                                                            *
13  !*****************************************************************************
14  !  CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with ECMWF grib_api    *
15  !  CHANGE: 03/12/2008, Harald Sodemann, change to f90 grib_api               *
16  !*****************************************************************************
17
18  use grib_api
19  use par_mod
20  use com_mod
21
22  implicit none
23
24  !HSO  parameters for grib_api
25  integer :: ifile
26  integer :: iret
27  integer :: igrib
28  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
[4fbe7a5]29  integer :: parID !added by mc for making it consistent with new gridcheck.f90
[e200b7a]30  integer :: gotGrib
31  !HSO  end
32  integer :: i,j,k,l,ifn,ifield,iumax,iwmax,numskip,nlev_ecn
33  integer :: nuvzn,nwzn
34  real :: akmn(nwzmax),bkmn(nwzmax),akzn(nuvzmax),bkzn(nuvzmax)
35  real(kind=4) :: xaux1,xaux2,yaux1,yaux2
36  real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
[4fbe7a5]37  real :: conversion_factor !added by mc to make it consistent with new gridchek.f90
[e200b7a]38
39  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
40
41  ! dimension of isec2 at least (22+n), where n is the number of parallels or
42  ! meridians in a quasi-regular (reduced) Gaussian or lat/long grid
43
44  ! dimension of zsec2 at least (10+nn), where nn is the number of vertical
45  ! coordinate parameters
46
47  integer :: isec1(56),isec2(22+nxmaxn+nymaxn)
48  real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
49
50  !HSO  grib api error messages
51  character(len=24) :: gribErrorMsg = 'Error reading grib file'
52  character(len=20) :: gribFunction = 'gridcheck_nests'
53
54  xresoln(0)=1.       ! resolution enhancement for mother grid
55  yresoln(0)=1.       ! resolution enhancement for mother grid
56
57  ! Loop about all nesting levels
58  !******************************
59
60  do l=1,numbnests
61
62    iumax=0
63    iwmax=0
64
65    if(ideltas.gt.0) then
66      ifn=1
67    else
68      ifn=numbwf
69    endif
70  !
71  ! OPENING OF DATA FILE (GRIB CODE)
72  !
73  ifile=0
74  igrib=0
75  iret=0
76
775   call grib_open_file(ifile,path(numpath+2*(l-1)+1) &
78         (1:length(numpath+2*(l-1)+1))//trim(wfnamen(l,ifn)),'r',iret)
79  if (iret.ne.GRIB_SUCCESS) then
80    goto 999   ! ERROR DETECTED
81  endif
82  !turn on support for multi fields messages
83  !call grib_multi_support_on
84
85  gotGrib=0
86  ifield=0
8710   ifield=ifield+1
88
89  !
90  ! GET NEXT FIELDS
91  !
92  call grib_new_from_file(ifile,igrib,iret)
93  if (iret.eq.GRIB_END_OF_FILE)  then
94    goto 30    ! EOF DETECTED
95  elseif (iret.ne.GRIB_SUCCESS) then
96    goto 999   ! ERROR DETECTED
97  endif
98
99  !first see if we read GRIB1 or GRIB2
100  call grib_get_int(igrib,'editionNumber',gribVer,iret)
101  call grib_check(iret,gribFunction,gribErrorMsg)
102
103  if (gribVer.eq.1) then ! GRIB Edition 1
104
105  !print*,'GRiB Edition 1'
106  !read the grib2 identifiers
107  call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
108  call grib_check(iret,gribFunction,gribErrorMsg)
109  call grib_get_int(igrib,'level',isec1(8),iret)
110  call grib_check(iret,gribFunction,gribErrorMsg)
111
112  !change code for etadot to code for omega
113  if (isec1(6).eq.77) then
114    isec1(6)=135
115  endif
116
117  !print*,isec1(6),isec1(8)
118
119  else
120
121  !print*,'GRiB Edition 2'
122  !read the grib2 identifiers
123  call grib_get_int(igrib,'discipline',discipl,iret)
124  call grib_check(iret,gribFunction,gribErrorMsg)
125  call grib_get_int(igrib,'parameterCategory',parCat,iret)
126  call grib_check(iret,gribFunction,gribErrorMsg)
127  call grib_get_int(igrib,'parameterNumber',parNum,iret)
128  call grib_check(iret,gribFunction,gribErrorMsg)
129  call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
130  call grib_check(iret,gribFunction,gribErrorMsg)
131  call grib_get_int(igrib,'level',valSurf,iret)
132  call grib_check(iret,gribFunction,gribErrorMsg)
[4fbe7a5]133  call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consisitent with new grid_check.f90
134  call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consisitent with new  grid_check.f90
[e200b7a]135
136  !print*,discipl,parCat,parNum,typSurf,valSurf
137
138  !convert to grib1 identifiers
139  isec1(6)=-1
140  isec1(7)=-1
141  isec1(8)=-1
142  isec1(8)=valSurf     ! level
143  if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
144    isec1(6)=130         ! indicatorOfParameter
145  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
146    isec1(6)=131         ! indicatorOfParameter
147  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
148    isec1(6)=132         ! indicatorOfParameter
149  elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
150    isec1(6)=133         ! indicatorOfParameter
[b0434e1]151  elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
152    isec1(6)=246         ! indicatorOfParameter
153  elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
154    isec1(6)=247         ! indicatorOfParameter
155!ZHG end
156! ESO qc(=clwc+ciwc)
157  elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
158    isec1(6)=201031      ! indicatorOfParameter
[e200b7a]159  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
160    isec1(6)=134         ! indicatorOfParameter
161  elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
162    isec1(6)=135         ! indicatorOfParameter
[db712a8]163  elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added bymc to make it consistent with new gridcheck.f90
164    isec1(6)=135         ! indicatorOfParameter    !
[e200b7a]165  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
166    isec1(6)=151         ! indicatorOfParameter
167  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
168    isec1(6)=165         ! indicatorOfParameter
169  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
170    isec1(6)=166         ! indicatorOfParameter
171  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
172    isec1(6)=167         ! indicatorOfParameter
173  elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
174    isec1(6)=168         ! indicatorOfParameter
175  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
176    isec1(6)=141         ! indicatorOfParameter
[4fbe7a5]177  elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC !added by mc to make it consistent with new gridchek.f90
[e200b7a]178    isec1(6)=164         ! indicatorOfParameter
[4fbe7a5]179 elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP !added by mc to make it consistent with new gridchek.f90
[e200b7a]180    isec1(6)=142         ! indicatorOfParameter
181  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
182    isec1(6)=143         ! indicatorOfParameter
183  elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
184    isec1(6)=146         ! indicatorOfParameter
185  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
186    isec1(6)=176         ! indicatorOfParameter
[4fbe7a5]187  elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS !added by mc to make it consistent with new gridchek.f90
[e200b7a]188    isec1(6)=180         ! indicatorOfParameter
[4fbe7a5]189  elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS !added by mc to make it consistent with new gridchek.f90
[e200b7a]190    isec1(6)=181         ! indicatorOfParameter
191  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
192    isec1(6)=129         ! indicatorOfParameter
[4fbe7a5]193  elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO !added by mc to make it consistent with new gridchek.f90
[e200b7a]194    isec1(6)=160         ! indicatorOfParameter
195  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
196       (typSurf.eq.1)) then ! LSM
197    isec1(6)=172         ! indicatorOfParameter
198  else
199    print*,'***ERROR: undefined GRiB2 message found!',discipl, &
200         parCat,parNum,typSurf
201  endif
[4fbe7a5]202  if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consistent with new gridchek.f90
203    write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
204!    stop
205  endif
[e200b7a]206
207  endif
208
209  !get the size and data of the values array
210  if (isec1(6).ne.-1) then
211    call grib_get_real4_array(igrib,'values',zsec4,iret)
212    call grib_check(iret,gribFunction,gribErrorMsg)
213  endif
214
215  !HSO  get the required fields from section 2 in a gribex compatible manner
216  if (ifield.eq.1) then
217    call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
218         isec2(2),iret)
219    call grib_check(iret,gribFunction,gribErrorMsg)
220    call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
221         isec2(3),iret)
222    call grib_check(iret,gribFunction,gribErrorMsg)
223    call grib_get_int(igrib,'numberOfVerticalCoordinateValues', &
224         isec2(12),iret)
225    call grib_check(iret,gribFunction,gribErrorMsg)
226  !HSO    get the size and data of the vertical coordinate array
227    call grib_get_real4_array(igrib,'pv',zsec2,iret)
228    call grib_check(iret,gribFunction,gribErrorMsg)
229
230    nxn(l)=isec2(2)
231    nyn(l)=isec2(3)
232    nlev_ecn=isec2(12)/2-1
233  endif ! ifield
234
[db712a8]235  if (nxn(l).gt.nxmaxn) then
236  write(*,*) 'FLEXPART error: Too many grid points in x direction.'
237  write(*,*) 'Reduce resolution of wind fields (file GRIDSPEC)'
238  write(*,*) 'for nesting level ',l
239  write(*,*) 'Or change parameter settings in file par_mod.'
240  write(*,*) nxn(l),nxmaxn
241  stop
242  endif
243
244  if (nyn(l).gt.nymaxn) then
245  write(*,*) 'FLEXPART error: Too many grid points in y direction.'
246  write(*,*) 'Reduce resolution of wind fields (file GRIDSPEC)'
247  write(*,*) 'for nesting level ',l
248  write(*,*) 'Or change parameter settings in file par_mod.'
249  write(*,*) nyn(l),nymaxn
250  stop
251  endif
252
[e200b7a]253  !HSO  get the second part of the grid dimensions only from GRiB1 messages
[4fbe7a5]254  if (isec1(6) .eq. 167 .and. (gotGrib.eq.0)) then !added by mc to make it consistent with new gridchek.f90 note that gotGrid must be changed in gotGrib!!
255    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', & !comment by mc: note that this was in the (if (ifield.eq.1) ..end above in gridchek.f90 see line 257
[e200b7a]256         xaux1in,iret)
257    call grib_check(iret,gribFunction,gribErrorMsg)
258    call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
259         xaux2in,iret)
260    call grib_check(iret,gribFunction,gribErrorMsg)
261    call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
262         yaux1in,iret)
263    call grib_check(iret,gribFunction,gribErrorMsg)
264    call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees', &
265         yaux2in,iret)
266    call grib_check(iret,gribFunction,gribErrorMsg)
267    xaux1=xaux1in
268    xaux2=xaux2in
269    yaux1=yaux1in
270    yaux2=yaux2in
[4fbe7a5]271    if(xaux1.gt.180.) xaux1=xaux1-360.0
272    if(xaux2.gt.180.) xaux2=xaux2-360.0
273    if(xaux1.lt.-180.) xaux1=xaux1+360.0
274    if(xaux2.lt.-180.) xaux2=xaux2+360.0
275    if (xaux2.lt.xaux1) xaux2=xaux2+360.0
[e200b7a]276    xlon0n(l)=xaux1
277    ylat0n(l)=yaux1
278    dxn(l)=(xaux2-xaux1)/real(nxn(l)-1)
279    dyn(l)=(yaux2-yaux1)/real(nyn(l)-1)
[4fbe7a5]280    gotGrib=1 !commetn by mc note tahthere gotGRIB is used instead of gotGrid!!!
[e200b7a]281  endif ! ifield.eq.1
282
283  k=isec1(8)
284  if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
285  if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
286
287  if(isec1(6).eq.129) then
288    do j=0,nyn(l)-1
289      do i=0,nxn(l)-1
290        oron(i,j,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
291      end do
292    end do
293  endif
294  if(isec1(6).eq.172) then
295    do j=0,nyn(l)-1
296      do i=0,nxn(l)-1
297        lsmn(i,j,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
298      end do
299    end do
300  endif
301  if(isec1(6).eq.160) then
302    do j=0,nyn(l)-1
303      do i=0,nxn(l)-1
304        excessoron(i,j,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
305      end do
306    end do
307  endif
308
309  call grib_release(igrib)
310  goto 10                 !! READ NEXT LEVEL OR PARAMETER
311  !
312  ! CLOSING OF INPUT DATA FILE
313  !
314
31530   call grib_close_file(ifile)
316
317  !error message if no fields found with correct first longitude in it
318  if (gotGrib.eq.0) then
319    print*,'***ERROR: input file needs to contain GRiB1 formatted'// &
320         'messages'
321    stop
322  endif
323
324  nuvzn=iumax
325  nwzn=iwmax
326  if(nuvzn.eq.nlev_ec) nwzn=nlev_ecn+1
327
328  if ((nuvzn.gt.nuvzmax).or.(nwzn.gt.nwzmax)) then
329  write(*,*) 'FLEXPART error: Nested wind fields have too many'// &
330       'vertical levels.'
331  write(*,*) 'Problem was encountered for nesting level ',l
332  stop
333  endif
334
335
336  ! Output of grid info
337  !********************
338
[b4d29ce]339  write(*,'(a,i2,a)') ' Nested domain ',l,':'
[4fbe7a5]340  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Longitude range: ', &
[e200b7a]341       xlon0n(l),' to ',xlon0n(l)+(nxn(l)-1)*dxn(l), &
342       '   Grid distance: ',dxn(l)
[b4d29ce]343  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range : ', &
[e200b7a]344       ylat0n(l),' to ',ylat0n(l)+(nyn(l)-1)*dyn(l), &
345       '   Grid distance: ',dyn(l)
346  write(*,*)
347
348  ! Determine, how much the resolutions in the nests are enhanced as
349  ! compared to the mother grid
350  !*****************************************************************
351
352  xresoln(l)=dx/dxn(l)
353  yresoln(l)=dy/dyn(l)
354
355  ! Determine the mother grid coordinates of the corner points of the
356  ! nested grids
357  ! Convert first to geographical coordinates, then to grid coordinates
358  !********************************************************************
359
360  xaux1=xlon0n(l)
361  xaux2=xlon0n(l)+real(nxn(l)-1)*dxn(l)
362  yaux1=ylat0n(l)
363  yaux2=ylat0n(l)+real(nyn(l)-1)*dyn(l)
364
365  xln(l)=(xaux1-xlon0)/dx
366  xrn(l)=(xaux2-xlon0)/dx
367  yln(l)=(yaux1-ylat0)/dy
368  yrn(l)=(yaux2-ylat0)/dy
369
370
371  if ((xln(l).lt.0.).or.(yln(l).lt.0.).or. &
372       (xrn(l).gt.real(nxmin1)).or.(yrn(l).gt.real(nymin1))) then
373    write(*,*) 'Nested domain does not fit into mother domain'
374    write(*,*) 'For global mother domain fields, you can shift'
375    write(*,*) 'shift the mother domain into x-direction'
376    write(*,*) 'by setting nxshift (file par_mod) to a'
377    write(*,*) 'positive value. Execution is terminated.'
378    stop
379  endif
380
381
382  ! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL
383  ! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
384
385  numskip=nlev_ecn-nuvzn ! number of ecmwf model layers not used by FLEXPART
386  do i=1,nwzn
387    j=numskip+i
388    k=nlev_ecn+1+numskip+i
389    akmn(nwzn-i+1)=zsec2(j)
390    bkmn(nwzn-i+1)=zsec2(k)
391  end do
392
393  !
394  ! CALCULATION OF AKZ, BKZ
395  ! AKZ,BKZ: model discretization parameters at the center of each model
396  !     layer
397  !
398  ! Assign the 10 m winds to an artificial model level with akz=0 and bkz=1.0,
399  ! i.e. ground level
400  !*****************************************************************************
401
402    akzn(1)=0.
403    bkzn(1)=1.0
404    do i=1,nuvzn
405      akzn(i+1)=0.5*(akmn(i+1)+akmn(i))
406      bkzn(i+1)=0.5*(bkmn(i+1)+bkmn(i))
407    end do
408    nuvzn=nuvzn+1
409
410  ! Check, whether the heights of the model levels of the nested
411  ! wind fields are consistent with those of the mother domain.
412  ! If not, terminate model run.
413  !*************************************************************
414
415    do i=1,nuvz
416      if ((akzn(i).ne.akz(i)).or.(bkzn(i).ne.bkz(i))) then
417  write(*,*) 'FLEXPART error: The wind fields of nesting level',l
418  write(*,*) 'are not consistent with the mother domain:'
419  write(*,*) 'Differences in vertical levels detected.'
420        stop
421      endif
422    end do
423
424    do i=1,nwz
425      if ((akmn(i).ne.akm(i)).or.(bkmn(i).ne.bkm(i))) then
426  write(*,*) 'FLEXPART error: The wind fields of nesting level',l
427  write(*,*) 'are not consistent with the mother domain:'
428  write(*,*) 'Differences in vertical levels detected.'
429        stop
430      endif
431    end do
432
433  end do
434
435  return
436
437999   write(*,*)
438  write(*,*) ' ###########################################'// &
439       '###### '
440  write(*,*) '                FLEXPART SUBROUTINE GRIDCHECK:'
441  write(*,*) ' CAN NOT OPEN INPUT DATA FILE '//wfnamen(l,ifn)
442  write(*,*) ' FOR NESTING LEVEL ',k
443  write(*,*) ' ###########################################'// &
444       '###### '
445  stop
446
447end subroutine gridcheck_nests
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG