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 #if defined WITH_CTBTO_PATCHES
91 integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
92 real :: sizesouth,sizenorth,xauxa,pint
102 integer :: isec1(56),isec2(22+nxmax+nymax)
103 real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
104 #if defined WITH_CTBTO_PATCHES
105 real(kind=4) :: conversion_factor
107 character(len=1) :: opt
110 character(len=24) :: griberrormsg =
'Error reading grib file'
111 character(len=20) :: gribfunction =
'gridcheck'
116 if(ideltas.gt.0)
then
124 5 call grib_open_file(ifile,path(3)(1:length(3)) &
125 //trim(wfname(ifn)),
'r',iret)
126 if (iret.ne.grib_success)
then
139 call grib_new_from_file(ifile,igrib,iret)
140 if (iret.eq.grib_end_of_file )
then
142 elseif (iret.ne.grib_success)
then
146 #if defined WITH_CTBTO_PATCHES
147 conversion_factor=1.0
151 call grib_get_int(igrib,
'editionNumber',gribver,iret)
152 call grib_check(iret,gribfunction,griberrormsg)
154 if (gribver.eq.1)
then
158 call grib_get_int(igrib,
'indicatorOfParameter',isec1(6),iret)
159 call grib_check(iret,gribfunction,griberrormsg)
160 call grib_get_int(igrib,
'level',isec1(8),iret)
161 call grib_check(iret,gribfunction,griberrormsg)
164 #if defined WITH_CTBTO_PATCHES
165 if (isec1(6).eq.77.or.isec1(6).eq.87)
then
167 if (isec1(6).eq.77)
then
176 #if defined WITH_CTBTO_PATCHES
177 call
grib2check(igrib,isec1,conversion_factor)
181 call grib_get_int(igrib,
'discipline',discipl,iret)
182 call grib_check(iret,gribfunction,griberrormsg)
183 call grib_get_int(igrib,
'parameterCategory',parcat,iret)
184 call grib_check(iret,gribfunction,griberrormsg)
185 call grib_get_int(igrib,
'parameterNumber',parnum,iret)
186 call grib_check(iret,gribfunction,griberrormsg)
187 call grib_get_int(igrib,
'typeOfFirstFixedSurface',typsurf,iret)
188 call grib_check(iret,gribfunction,griberrormsg)
189 call grib_get_int(igrib,
'level',valsurf,iret)
190 call grib_check(iret,gribfunction,griberrormsg)
199 if ((parcat.eq.0).and.(parnum.eq.0).and.(typsurf.eq.105))
then
201 elseif ((parcat.eq.2).and.(parnum.eq.2).and.(typsurf.eq.105))
then
203 elseif ((parcat.eq.2).and.(parnum.eq.3).and.(typsurf.eq.105))
then
205 elseif ((parcat.eq.1).and.(parnum.eq.0).and.(typsurf.eq.105))
then
207 elseif ((parcat.eq.3).and.(parnum.eq.0).and.(typsurf.eq.1))
then
209 elseif ((parcat.eq.2).and.(parnum.eq.32))
then
211 elseif ((parcat.eq.3).and.(parnum.eq.0).and.(typsurf.eq.101))
then
213 elseif ((parcat.eq.2).and.(parnum.eq.2).and.(typsurf.eq.103))
then
215 elseif ((parcat.eq.2).and.(parnum.eq.3).and.(typsurf.eq.103))
then
217 elseif ((parcat.eq.0).and.(parnum.eq.0).and.(typsurf.eq.103))
then
219 elseif ((parcat.eq.0).and.(parnum.eq.6).and.(typsurf.eq.103))
then
221 elseif ((parcat.eq.1).and.(parnum.eq.11).and.(typsurf.eq.1))
then
223 elseif ((parcat.eq.6).and.(parnum.eq.1))
then
225 elseif ((parcat.eq.1).and.(parnum.eq.9))
then
227 elseif ((parcat.eq.1).and.(parnum.eq.10))
then
229 elseif ((parcat.eq.0).and.(parnum.eq.11).and.(typsurf.eq.1))
then
231 elseif ((parcat.eq.4).and.(parnum.eq.9).and.(typsurf.eq.1))
then
233 elseif ((parcat.eq.2).and.(parnum.eq.17))
then
235 elseif ((parcat.eq.2).and.(parnum.eq.18))
then
237 elseif ((parcat.eq.3).and.(parnum.eq.4))
then
239 elseif ((parcat.eq.3).and.(parnum.eq.7))
then
241 elseif ((discipl.eq.2).and.(parcat.eq.0).and.(parnum.eq.0).and. &
245 print*,
'***ERROR: undefined GRiB2 message found!',discipl, &
246 parcat,parnum,typsurf
253 if (isec1(6).ne.-1)
then
254 call grib_get_real4_array(igrib,
'values',zsec4,iret)
255 #if defined WITH_CTBTO_PATCHES
256 zsec4=zsec4/conversion_factor
258 call grib_check(iret,gribfunction,griberrormsg)
261 if (ifield.eq.1)
then
264 call grib_get_int(igrib,
'numberOfPointsAlongAParallel', &
266 call grib_check(iret,gribfunction,griberrormsg)
267 call grib_get_int(igrib,
'numberOfPointsAlongAMeridian', &
269 call grib_check(iret,gribfunction,griberrormsg)
270 call grib_get_real8(igrib,
'longitudeOfFirstGridPointInDegrees', &
272 call grib_check(iret,gribfunction,griberrormsg)
273 call grib_get_int(igrib,
'numberOfVerticalCoordinateValues', &
275 call grib_check(iret,gribfunction,griberrormsg)
278 call grib_get_real4_array(igrib,
'pv',zsec2,iret)
279 call grib_check(iret,gribfunction,griberrormsg)
283 nlev_ec=isec2(12)/2-1
287 #if defined WITH_CTBTO_PATCHES
288 if (gotgrid.eq.0)
then
290 if ((gribver.eq.1).and.(gotgrid.eq.0))
then
292 call grib_get_real8(igrib,
'longitudeOfLastGridPointInDegrees', &
294 call grib_check(iret,gribfunction,griberrormsg)
295 call grib_get_real8(igrib,
'latitudeOfLastGridPointInDegrees', &
297 call grib_check(iret,gribfunction,griberrormsg)
298 call grib_get_real8(igrib,
'latitudeOfFirstGridPointInDegrees', &
300 call grib_check(iret,gribfunction,griberrormsg)
305 #if defined WITH_CTBTO_PATCHES
306 if (xaux1.ge.180.) xaux1=xaux1-360.0
308 if (xaux1.gt.180.) xaux1=xaux1-360.0
310 if (xaux2.gt.180.) xaux2=xaux2-360.0
311 if (xaux1.lt.-180.) xaux1=xaux1+360.0
312 if (xaux2.lt.-180.) xaux2=xaux2+360.0
313 if (xaux2.lt.xaux1) xaux2=xaux2+360.0
316 dx=(xaux2-xaux1)/
real(nxfield-1)
317 dy=(yaux2-yaux1)/
real(ny-1)
318 dxconst=180./(dx*r_earth*pi)
319 dyconst=180./(dy*r_earth*pi)
327 xauxa=abs(xaux2+dx-360.-xaux1)
328 if (xauxa.lt.0.001)
then
331 if (abs(nxshift).ge.nx) &
332 stop
'nxshift in file par_mod is too large'
333 xlon0=xlon0+
real(nxshift)*dx
338 stop
'nxshift (par_mod) must be zero for non-global domain'
342 if (xlon0.gt.180.) xlon0=xlon0-360.
344 if (xglobal.and.xauxa.lt.0.001)
then
348 sizesouth=6.*(switchsouth+90.)/dy
349 call
stlmbr(southpolemap,-90.,0.)
350 call
stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
351 sizesouth,switchsouth,180.)
352 switchsouthg=(switchsouth-ylat0)/dy
358 if (xglobal.and.xauxa.lt.0.001)
then
362 sizenorth=6.*(90.-switchnorth)/dy
363 call
stlmbr(northpolemap,90.,0.)
364 call
stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
365 sizenorth,switchnorth,180.)
366 switchnorthg=(switchnorth-ylat0)/dy
372 stop
'nxshift (par_mod) must not be negative'
373 if (nxshift.ge.nxfield) stop
'nxshift (par_mod) too large'
377 if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
378 if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
380 if(isec1(6).eq.129)
then
383 oro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)/ga
387 if(isec1(6).eq.172)
then
390 lsm(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
394 if(isec1(6).eq.160)
then
397 excessoro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
402 call grib_release(igrib)
408 30 call grib_close_file(ifile)
411 if (gotgrid.eq.0)
then
412 print*,
'***ERROR: input file needs to contain GRiB1 formatted'// &
419 if(nuvz.eq.nlev_ec) nwz=nlev_ec+1
421 if (nx.gt.nxmax)
then
422 write(*,*)
'FLEXPART error: Too many grid points in x direction.'
423 write(*,*)
'Reduce resolution of wind fields.'
424 write(*,*)
'Or change parameter settings in file par_mod.'
429 if (ny.gt.nymax)
then
430 write(*,*)
'FLEXPART error: Too many grid points in y direction.'
431 write(*,*)
'Reduce resolution of wind fields.'
432 write(*,*)
'Or change parameter settings in file par_mod.'
437 if (nuvz+1.gt.nuvzmax)
then
438 write(*,*)
'FLEXPART error: Too many u,v grid points in z '// &
440 write(*,*)
'Reduce resolution of wind fields.'
441 write(*,*)
'Or change parameter settings in file par_mod.'
442 write(*,*) nuvz+1,nuvzmax
446 if (nwz.gt.nwzmax)
then
447 write(*,*)
'FLEXPART error: Too many w grid points in z '// &
449 write(*,*)
'Reduce resolution of wind fields.'
450 write(*,*)
'Or change parameter settings in file par_mod.'
451 write(*,*) nwz,nwzmax
469 write(*,
'(a,2i7)')
'# of vertical levels in ECMWF data: ', &
472 write(*,
'(a)')
'Mother domain:'
473 write(*,
'(a,f10.2,a1,f10.2,a,f10.2)')
' Longitude range: ', &
474 xlon0,
' to ',xlon0+(nx-1)*dx,
' Grid distance: ',dx
475 write(*,
'(a,f10.2,a1,f10.2,a,f10.2)')
' Latitude range: ', &
476 ylat0,
' to ',ylat0+(ny-1)*dy,
' Grid distance: ',dy
495 k=nlev_ec+1+numskip+i
496 akm(nwz-i+1)=zsec2(j)
498 bkm(nwz-i+1)=zsec2(k)
513 akz(i+1)=0.5*(akm(i+1)+akm(i))
514 bkz(i+1)=0.5*(bkm(i+1)+bkm(i))
527 if (nz.gt.nzmax) stop
'nzmax too small'
551 pint=akz(i)+bkz(i)*101325.
552 if (pint.lt.5000.) goto 96
555 if (nconvlev.gt.nconvlevmax-1)
then
556 nconvlev=nconvlevmax-1
557 write(*,*)
'Attention, convection only calculated up to ', &
558 akz(nconvlev)+bkz(nconvlev)*1013.25,
' hPa'
564 write(*,*)
' ###########################################'// &
566 write(*,*)
' TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
567 write(*,*)
' CAN NOT OPEN INPUT DATA FILE '//wfname(ifn)
568 write(*,*)
' ###########################################'// &
571 write(*,
'(a)')
'!!! PLEASE INSERT A NEW CD-ROM AND !!!'
572 write(*,
'(a)')
'!!! PRESS ANY KEY TO CONTINUE... !!!'
573 write(*,
'(a)')
'!!! ...OR TERMINATE FLEXPART PRESSING!!!'
574 write(*,
'(a)')
'!!! THE "X" KEY... !!!'
subroutine, public stcm2p(strcmp, x1, y1, xlat1, xlong1, x2, y2, xlat2, xlong2)
subroutine gridcheck_ecmwf
subroutine, public stlmbr(strcmp, tnglat, xlong)
subroutine shift_field_0(field, nxf, nyf)
subroutine grib2check(igrib, isec1, conversion_factor)