source: trunk/src/gridcheck_nests.f90 @ 28

Last change on this file since 28 was 27, checked in by hasod, 10 years ago
  • Implemented optional namelist input for COMMAND, RELEASES, SPECIES, AGECLASSES,OUTGRID,OUTGRID_NEST,RECEPTORS
  • Implemented com_mod switch nmlout to write input files as namelist to the output directory (.true. by default)
  • Proposed updated startup and runtime output (may change back to previous info if desired)
  • Property svn:executable set to *
File size: 17.1 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.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
173    isec1(6)=134         ! indicatorOfParameter
174  elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
175    isec1(6)=135         ! indicatorOfParameter
176  elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added bymc to make it consistent with new gridchek.f90
177    isec1(6)=135         ! indicatorOfParameter    !                     ! " " " " " " " " " " " " " " " " " " " " " " " "  " " "
178  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
179    isec1(6)=151         ! indicatorOfParameter
180  elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
181    isec1(6)=165         ! indicatorOfParameter
182  elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
183    isec1(6)=166         ! indicatorOfParameter
184  elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
185    isec1(6)=167         ! indicatorOfParameter
186  elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
187    isec1(6)=168         ! indicatorOfParameter
188  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
189    isec1(6)=141         ! indicatorOfParameter
190  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
191    isec1(6)=164         ! indicatorOfParameter
192 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
193    isec1(6)=142         ! indicatorOfParameter
194  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
195    isec1(6)=143         ! indicatorOfParameter
196  elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
197    isec1(6)=146         ! indicatorOfParameter
198  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
199    isec1(6)=176         ! indicatorOfParameter
200  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
201    isec1(6)=180         ! indicatorOfParameter
202  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
203    isec1(6)=181         ! indicatorOfParameter
204  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
205    isec1(6)=129         ! indicatorOfParameter
206  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
207    isec1(6)=160         ! indicatorOfParameter
208  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
209       (typSurf.eq.1)) then ! LSM
210    isec1(6)=172         ! indicatorOfParameter
211  else
212    print*,'***ERROR: undefined GRiB2 message found!',discipl, &
213         parCat,parNum,typSurf
214  endif
215  if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consistent with new gridchek.f90
216    write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
217!    stop
218  endif
219
220  endif
221
222  !get the size and data of the values array
223  if (isec1(6).ne.-1) then
224    call grib_get_real4_array(igrib,'values',zsec4,iret)
225    call grib_check(iret,gribFunction,gribErrorMsg)
226  endif
227
228  !HSO  get the required fields from section 2 in a gribex compatible manner
229  if (ifield.eq.1) then
230    call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
231         isec2(2),iret)
232    call grib_check(iret,gribFunction,gribErrorMsg)
233    call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
234         isec2(3),iret)
235    call grib_check(iret,gribFunction,gribErrorMsg)
236    call grib_get_int(igrib,'numberOfVerticalCoordinateValues', &
237         isec2(12),iret)
238    call grib_check(iret,gribFunction,gribErrorMsg)
239  !HSO    get the size and data of the vertical coordinate array
240    call grib_get_real4_array(igrib,'pv',zsec2,iret)
241    call grib_check(iret,gribFunction,gribErrorMsg)
242
243    nxn(l)=isec2(2)
244    nyn(l)=isec2(3)
245    nlev_ecn=isec2(12)/2-1
246  endif ! ifield
247
248  !HSO  get the second part of the grid dimensions only from GRiB1 messages
249  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!!
250    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
251         xaux1in,iret)
252    call grib_check(iret,gribFunction,gribErrorMsg)
253    call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
254         xaux2in,iret)
255    call grib_check(iret,gribFunction,gribErrorMsg)
256    call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
257         yaux1in,iret)
258    call grib_check(iret,gribFunction,gribErrorMsg)
259    call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees', &
260         yaux2in,iret)
261    call grib_check(iret,gribFunction,gribErrorMsg)
262    xaux1=xaux1in
263    xaux2=xaux2in
264    yaux1=yaux1in
265    yaux2=yaux2in
266    if(xaux1.gt.180.) xaux1=xaux1-360.0
267    if(xaux2.gt.180.) xaux2=xaux2-360.0
268    if(xaux1.lt.-180.) xaux1=xaux1+360.0
269    if(xaux2.lt.-180.) xaux2=xaux2+360.0
270    if (xaux2.lt.xaux1) xaux2=xaux2+360.0
271    xlon0n(l)=xaux1
272    ylat0n(l)=yaux1
273    dxn(l)=(xaux2-xaux1)/real(nxn(l)-1)
274    dyn(l)=(yaux2-yaux1)/real(nyn(l)-1)
275    gotGrib=1 !commetn by mc note tahthere gotGRIB is used instead of gotGrid!!!
276  endif ! ifield.eq.1
277
278  k=isec1(8)
279  if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
280  if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
281
282  if(isec1(6).eq.129) then
283    do j=0,nyn(l)-1
284      do i=0,nxn(l)-1
285        oron(i,j,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
286      end do
287    end do
288  endif
289  if(isec1(6).eq.172) then
290    do j=0,nyn(l)-1
291      do i=0,nxn(l)-1
292        lsmn(i,j,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
293      end do
294    end do
295  endif
296  if(isec1(6).eq.160) then
297    do j=0,nyn(l)-1
298      do i=0,nxn(l)-1
299        excessoron(i,j,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
300      end do
301    end do
302  endif
303
304  call grib_release(igrib)
305  goto 10                 !! READ NEXT LEVEL OR PARAMETER
306  !
307  ! CLOSING OF INPUT DATA FILE
308  !
309
31030   call grib_close_file(ifile)
311
312  !error message if no fields found with correct first longitude in it
313  if (gotGrib.eq.0) then
314    print*,'***ERROR: input file needs to contain GRiB1 formatted'// &
315         'messages'
316    stop
317  endif
318
319  nuvzn=iumax
320  nwzn=iwmax
321  if(nuvzn.eq.nlev_ec) nwzn=nlev_ecn+1
322
323  if (nxn(l).gt.nxmaxn) then
324  write(*,*) 'FLEXPART error: Too many grid points in x direction.'
325  write(*,*) 'Reduce resolution of wind fields (file GRIDSPEC)'
326  write(*,*) 'for nesting level ',l
327  write(*,*) 'Or change parameter settings in file par_mod.'
328  write(*,*) nxn(l),nxmaxn
329  stop
330  endif
331
332  if (nyn(l).gt.nymaxn) then
333  write(*,*) 'FLEXPART error: Too many grid points in y direction.'
334  write(*,*) 'Reduce resolution of wind fields (file GRIDSPEC)'
335  write(*,*) 'for nesting level ',l
336  write(*,*) 'Or change parameter settings in file par_mod.'
337  write(*,*) nyn(l),nymaxn
338  stop
339  endif
340
341  if ((nuvzn.gt.nuvzmax).or.(nwzn.gt.nwzmax)) then
342  write(*,*) 'FLEXPART error: Nested wind fields have too many'// &
343       'vertical levels.'
344  write(*,*) 'Problem was encountered for nesting level ',l
345  stop
346  endif
347
348
349  ! Output of grid info
350  !********************
351
352  write(*,'(a,i2,a)') ' Nested domain ',l,':'
353  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Longitude range: ', &
354       xlon0n(l),' to ',xlon0n(l)+(nxn(l)-1)*dxn(l), &
355       '   Grid distance: ',dxn(l)
356  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range : ', &
357       ylat0n(l),' to ',ylat0n(l)+(nyn(l)-1)*dyn(l), &
358       '   Grid distance: ',dyn(l)
359  write(*,*)
360
361  ! Determine, how much the resolutions in the nests are enhanced as
362  ! compared to the mother grid
363  !*****************************************************************
364
365  xresoln(l)=dx/dxn(l)
366  yresoln(l)=dy/dyn(l)
367
368  ! Determine the mother grid coordinates of the corner points of the
369  ! nested grids
370  ! Convert first to geographical coordinates, then to grid coordinates
371  !********************************************************************
372
373  xaux1=xlon0n(l)
374  xaux2=xlon0n(l)+real(nxn(l)-1)*dxn(l)
375  yaux1=ylat0n(l)
376  yaux2=ylat0n(l)+real(nyn(l)-1)*dyn(l)
377
378  xln(l)=(xaux1-xlon0)/dx
379  xrn(l)=(xaux2-xlon0)/dx
380  yln(l)=(yaux1-ylat0)/dy
381  yrn(l)=(yaux2-ylat0)/dy
382
383
384  if ((xln(l).lt.0.).or.(yln(l).lt.0.).or. &
385       (xrn(l).gt.real(nxmin1)).or.(yrn(l).gt.real(nymin1))) then
386    write(*,*) 'Nested domain does not fit into mother domain'
387    write(*,*) 'For global mother domain fields, you can shift'
388    write(*,*) 'shift the mother domain into x-direction'
389    write(*,*) 'by setting nxshift (file par_mod) to a'
390    write(*,*) 'positive value. Execution is terminated.'
391    stop
392  endif
393
394
395  ! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL
396  ! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
397
398  numskip=nlev_ecn-nuvzn ! number of ecmwf model layers not used by FLEXPART
399  do i=1,nwzn
400    j=numskip+i
401    k=nlev_ecn+1+numskip+i
402    akmn(nwzn-i+1)=zsec2(j)
403    bkmn(nwzn-i+1)=zsec2(k)
404  end do
405
406  !
407  ! CALCULATION OF AKZ, BKZ
408  ! AKZ,BKZ: model discretization parameters at the center of each model
409  !     layer
410  !
411  ! Assign the 10 m winds to an artificial model level with akz=0 and bkz=1.0,
412  ! i.e. ground level
413  !*****************************************************************************
414
415    akzn(1)=0.
416    bkzn(1)=1.0
417    do i=1,nuvzn
418      akzn(i+1)=0.5*(akmn(i+1)+akmn(i))
419      bkzn(i+1)=0.5*(bkmn(i+1)+bkmn(i))
420    end do
421    nuvzn=nuvzn+1
422
423  ! Check, whether the heights of the model levels of the nested
424  ! wind fields are consistent with those of the mother domain.
425  ! If not, terminate model run.
426  !*************************************************************
427
428    do i=1,nuvz
429      if ((akzn(i).ne.akz(i)).or.(bkzn(i).ne.bkz(i))) then
430  write(*,*) 'FLEXPART error: The wind fields of nesting level',l
431  write(*,*) 'are not consistent with the mother domain:'
432  write(*,*) 'Differences in vertical levels detected.'
433        stop
434      endif
435    end do
436
437    do i=1,nwz
438      if ((akmn(i).ne.akm(i)).or.(bkmn(i).ne.bkm(i))) then
439  write(*,*) 'FLEXPART error: The wind fields of nesting level',l
440  write(*,*) 'are not consistent with the mother domain:'
441  write(*,*) 'Differences in vertical levels detected.'
442        stop
443      endif
444    end do
445
446  end do
447
448  return
449
450999   write(*,*)
451  write(*,*) ' ###########################################'// &
452       '###### '
453  write(*,*) '                FLEXPART SUBROUTINE GRIDCHECK:'
454  write(*,*) ' CAN NOT OPEN INPUT DATA FILE '//wfnamen(l,ifn)
455  write(*,*) ' FOR NESTING LEVEL ',k
456  write(*,*) ' ###########################################'// &
457       '###### '
458  stop
459
460end subroutine gridcheck_nests
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG