85 real(kind=4) :: xaux1,xaux2,yaux1,yaux2
86 real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
87 integer :: gribver,parcat,parnum,typsurf,valsurf,discipl
89 integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
90 real :: sizesouth,sizenorth,xauxa,pint
91 real :: akm_usort(nwzmax)
92 real,
parameter :: eps=0.0001
95 real :: pres(nwzmax), help
97 integer :: i179,i180,i181
102 real(kind=4) :: zsec4(jpunp)
103 character(len=1) :: opt
106 character(len=24) :: griberrormsg =
'Error reading grib file'
107 character(len=20) :: gribfunction =
'gridcheckwind_gfs'
117 character(LEN=255),
parameter :: vtable_ncep_grib1_path = &
118 "Vtable_ncep_grib1", &
119 vtable_ncep_grib2_path = &
124 integer :: gribfile_type
125 integer :: current_grib_level
126 character(len=255) :: gribfile_name
127 character(len=255) :: vtable_path
128 character(len=15) :: fpname
140 if (numbnests.ge.1)
then
141 write(*,*)
' ###########################################'
142 write(*,*)
' FLEXPART ERROR SUBROUTINE GRIDCHECK:'
143 write(*,*)
' NO NESTED WINDFIELDAS ALLOWED FOR GFS! '
144 write(*,*)
' ###########################################'
151 if(ideltas.gt.0)
then
162 gribfile_name = path(3)(1:length(3))//trim(wfname(ifn))
163 print *,
'gribfile_name: ', gribfile_name
167 print *,
'gribfile_type: ', gribfile_type
169 if (gribfile_type .eq. vtable_gribfile_type_ncep_grib1)
then
170 vtable_path = vtable_ncep_grib1_path
171 else if (gribfile_type .eq. vtable_gribfile_type_ncep_grib2)
then
172 vtable_path = vtable_ncep_grib2_path
174 print *,
'Unsupported gribfile_type: ', gribfile_type
180 print *,
'Loading Vtable: ', vtable_path
182 print *,
'Vtable Initialized: ', my_vtable%initialized
183 print *,
'Vtable num_entries: ', my_vtable%num_entries
196 5 call grib_open_file(ifile,path(3)(1:length(3)) &
197 //trim(wfname(ifn)),
'r',iret)
198 if (iret.ne.grib_success)
then
202 call grib_multi_support_on
209 call grib_new_from_file(ifile,igrib,iret)
210 if (iret.eq.grib_end_of_file )
then
212 elseif (iret.ne.grib_success)
then
230 call grib_get_int(igrib,
'editionNumber',gribver,iret)
231 call grib_check(iret,gribfunction,griberrormsg)
233 if (gribver.eq.1)
then
235 call grib_get_int(igrib,
'level',current_grib_level,iret)
236 call grib_check(iret,gribfunction,griberrormsg)
240 current_grib_level = current_grib_level*100.0
244 call grib_get_int(igrib,
'scaledValueOfFirstFixedSurface', &
245 current_grib_level,iret)
246 call grib_check(iret,gribfunction,griberrormsg)
253 if ( trim(fpname) .ne.
'None' )
then
255 call grib_get_real4_array(igrib,
'values',zsec4,iret)
256 call grib_check(iret,gribfunction,griberrormsg)
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_real8(igrib,
'longitudeOfLastGridPointInDegrees', &
274 call grib_check(iret,gribfunction,griberrormsg)
275 call grib_get_real8(igrib,
'latitudeOfLastGridPointInDegrees', &
277 call grib_check(iret,gribfunction,griberrormsg)
278 call grib_get_real8(igrib,
'latitudeOfFirstGridPointInDegrees', &
280 call grib_check(iret,gribfunction,griberrormsg)
286 if (xaux2in .lt. 0) xaux2in = 359.0
298 if((abs(xaux1).lt.eps).and.(xaux2.ge.359))
then
300 xaux2=-179.0+360.-360./
real(nxfield)
302 if (xaux1.gt.180) xaux1=xaux1-360.0
303 if (xaux2.gt.180) xaux2=xaux2-360.0
304 if (xaux1.lt.-180) xaux1=xaux1+360.0
305 if (xaux2.lt.-180) xaux2=xaux2+360.0
306 if (xaux2.lt.xaux1) xaux2=xaux2+360.
309 dx=(xaux2-xaux1)/
real(nxfield-1)
310 dy=(yaux2-yaux1)/
real(ny-1)
311 dxconst=180./(dx*r_earth*pi)
312 dyconst=180./(dy*r_earth*pi)
321 xauxa=abs(xaux2+dx-360.-xaux1)
322 if (xauxa.lt.0.001)
then
325 if (abs(nxshift).ge.nx) &
326 stop
'nxshift in file par_mod is too large'
327 xlon0=xlon0+
real(nxshift)*dx
332 stop
'nxshift (par_mod) must be zero for non-global domain'
336 if (xlon0.gt.180.) xlon0=xlon0-360.
338 if (xglobal.and.xauxa.lt.0.001)
then
342 sizesouth=6.*(switchsouth+90.)/dy
343 call
stlmbr(southpolemap,-90.,0.)
344 call
stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
345 sizesouth,switchsouth,180.)
346 switchsouthg=(switchsouth-ylat0)/dy
352 if (xglobal.and.xauxa.lt.0.001)
then
356 sizenorth=6.*(90.-switchnorth)/dy
357 call
stlmbr(northpolemap,90.,0.)
358 call
stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
359 sizenorth,switchnorth,180.)
360 switchnorthg=(switchnorth-ylat0)/dy
367 if (nxshift.lt.0) stop
'nxshift (par_mod) must not be negative'
368 if (nxshift.ge.nxfield) stop
'nxshift (par_mod) too large'
373 if( trim(fpname) .eq.
'UU' )
then
375 pres(iumax) =
real(current_grib_level)
391 if( trim(fpname) .eq.
'ORO' )
then
394 help=zsec4(nxfield*(ny-jy-1)+ix+1)
397 excessoro(i179+ix,jy)=0.0
400 excessoro(ix-i181,jy)=0.0
409 if( trim(fpname) .eq.
'LSM' )
then
412 help=zsec4(nxfield*(ny-jy-1)+ix+1)
422 call grib_release(igrib)
431 call grib_close_file(ifile)
438 if (nx.gt.nxmax)
then
439 write(*,*)
'FLEXPART error: Too many grid points in x direction.'
440 write(*,*)
'Reduce resolution of wind fields.'
441 write(*,*)
'Or change parameter settings in file par_mod.'
446 if (ny.gt.nymax)
then
447 write(*,*)
'FLEXPART error: Too many grid points in y direction.'
448 write(*,*)
'Reduce resolution of wind fields.'
449 write(*,*)
'Or change parameter settings in file par_mod.'
454 if (nuvz.gt.nuvzmax)
then
455 write(*,*)
'FLEXPART error: Too many u,v grid points in z '// &
457 write(*,*)
'Reduce resolution of wind fields.'
458 write(*,*)
'Or change parameter settings in file par_mod.'
459 write(*,*) nuvz,nuvzmax
463 if (nwz.gt.nwzmax)
then
464 write(*,*)
'FLEXPART error: Too many w grid points in z '// &
466 write(*,*)
'Reduce resolution of wind fields.'
467 write(*,*)
'Or change parameter settings in file par_mod.'
468 write(*,*) nwz,nwzmax
486 write(*,
'(a,2i7)')
'# of vertical levels in NCEP data: ', &
489 write(*,
'(a)')
'Mother domain:'
490 write(*,
'(a,f10.2,a1,f10.2,a,f10.2)')
' Longitude range: ', &
491 xlon0,
' to ',xlon0+(nx-1)*dx,
' Grid distance: ',dx
492 write(*,
'(a,f10.2,a1,f10.2,a,f10.2)')
' Latitude range: ', &
493 ylat0,
' to ',ylat0+(ny-1)*dy,
' Grid distance: ',dy
504 k=nlev_ec+1+numskip+i
505 akm_usort(nwz-i+1)=pres(nwz-i+1)
513 if (akm_usort(1).gt.akm_usort(2))
then
516 akm(i)=akm_usort(nwz-i+1)
543 if (nz.gt.nzmax) stop
'nzmax too small'
567 pint=akz(i)+bkz(i)*101325.
568 if (pint.lt.5000.) goto 96
571 if (nconvlev.gt.nconvlevmax-1)
then
572 nconvlev=nconvlevmax-1
573 write(*,*)
'Attention, convection only calculated up to ', &
574 akz(nconvlev)+bkz(nconvlev)*1013.25,
' hPa'
580 write(*,*)
' ###########################################'// &
582 write(*,*)
' TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
583 write(*,*)
' CAN NOT OPEN INPUT DATA FILE '//wfname(ifn)
584 write(*,*)
' ###########################################'// &
587 write(*,
'(a)')
'!!! PLEASE INSERT A NEW CD-ROM AND !!!'
588 write(*,
'(a)')
'!!! PRESS ANY KEY TO CONTINUE... !!!'
589 write(*,
'(a)')
'!!! ...OR TERMINATE FLEXPART PRESSING!!!'
590 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, 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)