source: flexpart.git/src/gridcheck_nests.f90 @ 02095e3

devrelease-10univie
Last change on this file since 02095e3 was db712a8, checked in by Espen Sollum ATMOS <eso@…>, 4 years ago

Completed handling of nested wind fields with cloud water (for wet deposition).

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