85 real(kind=4) :: xaux1,xaux2,yaux1,yaux2
86 real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
87 integer :: gribver,parcat,parnum,typsurf,valsurf,discipl
88 #if defined WITH_CTBTO_PATCHES
92 integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
93 real :: sizesouth,sizenorth,xauxa,pint
103 integer :: isec1(56),isec2(22+nxmax+nymax)
104 real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
105 #if defined WITH_CTBTO_PATCHES
106 real(kind=4) :: conversion_factor
108 character(len=1) :: opt
111 character(len=24) :: griberrormsg =
'Error reading grib file'
112 character(len=20) :: gribfunction =
'gridcheck'
121 character(LEN=255),
parameter :: vtable_ecmwf_grib1_path = &
122 "Vtable_ecmwf_grib1", &
123 vtable_ecmwf_grib2_path = &
124 "Vtable_ecmwf_grib2", &
125 vtable_ecmwf_grib1_2_path = &
126 "Vtable_ecmwf_grib1_2"
128 integer :: gribfile_type
129 integer :: current_grib_level
130 character(len=255) :: gribfile_name
131 character(len=255) :: vtable_path
132 character(len=15) :: fpname
146 if(ideltas.gt.0)
then
156 gribfile_name = path(3)(1:length(3))//trim(wfname(ifn))
157 print *,
'gribfile_name: ', gribfile_name
161 print *,
'gribfile_type: ', gribfile_type
163 if (gribfile_type .eq. vtable_gribfile_type_ecmwf_grib1)
then
164 vtable_path = vtable_ecmwf_grib1_path
165 else if (gribfile_type .eq. vtable_gribfile_type_ecmwf_grib2)
then
166 vtable_path = vtable_ecmwf_grib2_path
167 else if (gribfile_type .eq. vtable_gribfile_type_ecmwf_grib1_2)
then
168 vtable_path = vtable_ecmwf_grib1_2_path
170 print *,
'Unsupported gribfile_type: ', gribfile_type
176 print *,
'Loading Vtable: ', vtable_path
178 print *,
'Vtable Initialized: ', my_vtable%initialized
179 print *,
'Vtable num_entries: ', my_vtable%num_entries
189 5 call grib_open_file(ifile,path(3)(1:length(3)) &
190 //trim(wfname(ifn)),
'r',iret)
191 if (iret.ne.grib_success)
then
204 call grib_new_from_file(ifile,igrib,iret)
205 if (iret.eq.grib_end_of_file )
then
207 elseif (iret.ne.grib_success)
then
224 #if defined WITH_CTBTO_PATCHES
225 conversion_factor=1.0
229 call grib_get_int(igrib,
'editionNumber',gribver,iret)
230 call grib_check(iret,gribfunction,griberrormsg)
233 if (gribver.eq.1)
then
237 call grib_get_int(igrib,
'level', current_grib_level, iret)
238 call grib_check(iret,gribfunction,griberrormsg)
242 #if defined WITH_CTBTO_PATCHES
243 call
grib2check(igrib, fpname, conversion_factor)
246 call grib_get_int(igrib,
'level', current_grib_level, iret)
247 call grib_check(iret,gribfunction,griberrormsg)
252 if (trim(fpname) .ne.
'None')
then
253 call grib_get_real4_array(igrib,
'values',zsec4,iret)
254 #if defined WITH_CTBTO_PATCHES
255 zsec4=zsec4/conversion_factor
257 call grib_check(iret,gribfunction,griberrormsg)
260 if (ifield.eq.1)
then
263 call grib_get_int(igrib,
'numberOfPointsAlongAParallel', &
265 call grib_check(iret,gribfunction,griberrormsg)
266 call grib_get_int(igrib,
'numberOfPointsAlongAMeridian', &
268 call grib_check(iret,gribfunction,griberrormsg)
269 call grib_get_real8(igrib,
'longitudeOfFirstGridPointInDegrees', &
271 call grib_check(iret,gribfunction,griberrormsg)
272 call grib_get_int(igrib,
'numberOfVerticalCoordinateValues', &
274 call grib_check(iret,gribfunction,griberrormsg)
277 call grib_get_real4_array(igrib,
'pv',zsec2,iret)
278 call grib_check(iret,gribfunction,griberrormsg)
282 nlev_ec=isec2(12)/2-1
286 #if defined WITH_CTBTO_PATCHES
291 call grib_get_real8(igrib,
'longitudeOfLastGridPointInDegrees', &
293 call grib_check(iret,gribfunction,griberrormsg)
294 call grib_get_real8(igrib,
'latitudeOfLastGridPointInDegrees', &
296 call grib_check(iret,gribfunction,griberrormsg)
297 call grib_get_real8(igrib,
'latitudeOfFirstGridPointInDegrees', &
299 call grib_check(iret,gribfunction,griberrormsg)
304 #if defined WITH_CTBTO_PATCHES
305 if (xaux1.ge.180.) xaux1=xaux1-360.0
307 if (xaux1.gt.180.) xaux1=xaux1-360.0
309 if (xaux2.gt.180.) xaux2=xaux2-360.0
310 if (xaux1.lt.-180.) xaux1=xaux1+360.0
311 if (xaux2.lt.-180.) xaux2=xaux2+360.0
312 if (xaux2.lt.xaux1) xaux2=xaux2+360.0
315 dx=(xaux2-xaux1)/
real(nxfield-1)
316 dy=(yaux2-yaux1)/
real(ny-1)
317 dxconst=180./(dx*r_earth*pi)
318 dyconst=180./(dy*r_earth*pi)
326 xauxa=abs(xaux2+dx-360.-xaux1)
327 if (xauxa.lt.0.001)
then
330 if (abs(nxshift).ge.nx) &
331 stop
'nxshift in file par_mod is too large'
332 xlon0=xlon0+
real(nxshift)*dx
337 stop
'nxshift (par_mod) must be zero for non-global domain'
341 if (xlon0.gt.180.) xlon0=xlon0-360.
343 if (xglobal.and.xauxa.lt.0.001)
then
347 sizesouth=6.*(switchsouth+90.)/dy
348 call
stlmbr(southpolemap,-90.,0.)
349 call
stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
350 sizesouth,switchsouth,180.)
351 switchsouthg=(switchsouth-ylat0)/dy
357 if (xglobal.and.xauxa.lt.0.001)
then
361 sizenorth=6.*(90.-switchnorth)/dy
362 call
stlmbr(northpolemap,90.,0.)
363 call
stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
364 sizenorth,switchnorth,180.)
365 switchnorthg=(switchnorth-ylat0)/dy
371 stop
'nxshift (par_mod) must not be negative'
372 if (nxshift.ge.nxfield) stop
'nxshift (par_mod) too large'
376 if(trim(fpname) .eq.
'UU') iumax=max(iumax,nlev_ec-k+1)
377 if(trim(fpname) .eq.
'ETADOT') iwmax=max(iwmax,nlev_ec-k+1)
379 if(trim(fpname) .eq.
'ORO')
then
382 oro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)/ga
386 if(trim(fpname) .eq.
'LSM')
then
389 lsm(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
393 if(trim(fpname) .eq.
'EXCESSORO')
then
396 excessoro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
401 call grib_release(igrib)
407 30 call grib_close_file(ifile)
410 if (gotgrid.eq.0)
then
411 print*,
'***gridcheck() ERROR: no fields found with correct first longitude'// &
418 if(nuvz.eq.nlev_ec) nwz=nlev_ec+1
420 if (nx.gt.nxmax)
then
421 write(*,*)
'FLEXPART error: Too many grid points in x direction.'
422 write(*,*)
'Reduce resolution of wind fields.'
423 write(*,*)
'Or change parameter settings in file par_mod.'
428 if (ny.gt.nymax)
then
429 write(*,*)
'FLEXPART error: Too many grid points in y direction.'
430 write(*,*)
'Reduce resolution of wind fields.'
431 write(*,*)
'Or change parameter settings in file par_mod.'
436 if (nuvz+1.gt.nuvzmax)
then
437 write(*,*)
'FLEXPART error: Too many u,v grid points in z '// &
439 write(*,*)
'Reduce resolution of wind fields.'
440 write(*,*)
'Or change parameter settings in file par_mod.'
441 write(*,*) nuvz+1,nuvzmax
445 if (nwz.gt.nwzmax)
then
446 write(*,*)
'FLEXPART error: Too many w grid points in z '// &
448 write(*,*)
'Reduce resolution of wind fields.'
449 write(*,*)
'Or change parameter settings in file par_mod.'
450 write(*,*) nwz,nwzmax
468 write(*,
'(a,2i7)')
'# of vertical levels in ECMWF data: ', &
471 write(*,
'(a)')
'Mother domain:'
472 write(*,
'(a,f10.2,a1,f10.2,a,f10.2)')
' Longitude range: ', &
473 xlon0,
' to ',xlon0+(nx-1)*dx,
' Grid distance: ',dx
474 write(*,
'(a,f10.2,a1,f10.2,a,f10.2)')
' Latitude range: ', &
475 ylat0,
' to ',ylat0+(ny-1)*dy,
' Grid distance: ',dy
494 k=nlev_ec+1+numskip+i
495 akm(nwz-i+1)=zsec2(j)
497 bkm(nwz-i+1)=zsec2(k)
512 akz(i+1)=0.5*(akm(i+1)+akm(i))
513 bkz(i+1)=0.5*(bkm(i+1)+bkm(i))
526 if (nz.gt.nzmax) stop
'nzmax too small'
550 pint=akz(i)+bkz(i)*101325.
551 if (pint.lt.5000.) goto 96
554 if (nconvlev.gt.nconvlevmax-1)
then
555 nconvlev=nconvlevmax-1
556 write(*,*)
'Attention, convection only calculated up to ', &
557 akz(nconvlev)+bkz(nconvlev)*1013.25,
' hPa'
563 write(*,*)
' ###########################################'// &
565 write(*,*)
' TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
566 write(*,*)
' CAN NOT OPEN INPUT DATA FILE '//wfname(ifn)
567 write(*,*)
' ###########################################'// &
570 write(*,
'(a)')
'!!! PLEASE INSERT A NEW CD-ROM AND !!!'
571 write(*,
'(a)')
'!!! PRESS ANY KEY TO CONTINUE... !!!'
572 write(*,
'(a)')
'!!! ...OR TERMINATE FLEXPART PRESSING!!!'
573 write(*,
'(a)')
'!!! THE "X" KEY... !!!'
subroutine, public stcm2p(strcmp, x1, y1, xlat1, xlong1, x2, y2, xlat2, xlong2)
character(len=15) function, public vtable_get_fpname(igrib, vtable_object)
subroutine gridcheck_ecmwf
subroutine grib2check(igrib, fpname, conversion_factor)
subroutine, public stlmbr(strcmp, tnglat, xlong)
subroutine, public vtable_load_by_name(vtable_name, the_vtable_data)
subroutine shift_field_0(field, nxf, nyf)
integer function, public vtable_detect_gribfile_type(gribfilename)