84 real(kind=4) :: xaux1,xaux2,yaux1,yaux2
85 real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
86 integer :: gribver,parcat,parnum,typsurf,valsurf,discipl
88 integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
89 real :: sizesouth,sizenorth,xauxa,pint
90 real :: akm_usort(nwzmax)
91 real,
parameter :: eps=0.0001
94 real :: pres(nwzmax), help
96 integer :: i179,i180,i181
100 integer :: isec1(8),isec2(3)
101 real(kind=4) :: zsec4(jpunp)
102 character(len=1) :: opt
105 character(len=24) :: griberrormsg =
'Error reading grib file'
106 character(len=20) :: gribfunction =
'gridcheckwind_gfs'
108 if (numbnests.ge.1)
then
109 write(*,*)
' ###########################################'
110 write(*,*)
' FLEXPART ERROR SUBROUTINE GRIDCHECK:'
111 write(*,*)
' NO NESTED WINDFIELDAS ALLOWED FOR GFS! '
112 write(*,*)
' ###########################################'
119 if(ideltas.gt.0)
then
127 5 call grib_open_file(ifile,path(3)(1:length(3)) &
128 //trim(wfname(ifn)),
'r',iret)
129 if (iret.ne.grib_success)
then
133 call grib_multi_support_on
140 call grib_new_from_file(ifile,igrib,iret)
141 if (iret.eq.grib_end_of_file )
then
143 elseif (iret.ne.grib_success)
then
148 call grib_get_int(igrib,
'editionNumber',gribver,iret)
149 call grib_check(iret,gribfunction,griberrormsg)
151 if (gribver.eq.1)
then
154 call grib_get_int(igrib,
'indicatorOfParameter',isec1(6),iret)
155 call grib_check(iret,gribfunction,griberrormsg)
156 call grib_get_int(igrib,
'indicatorOfTypeOfLevel',isec1(7),iret)
157 call grib_check(iret,gribfunction,griberrormsg)
158 call grib_get_int(igrib,
'level',isec1(8),iret)
159 call grib_check(iret,gribfunction,griberrormsg)
162 call grib_get_real4_array(igrib,
'values',zsec4,iret)
163 call grib_check(iret,gribfunction,griberrormsg)
168 call grib_get_int(igrib,
'discipline',discipl,iret)
169 call grib_check(iret,gribfunction,griberrormsg)
170 call grib_get_int(igrib,
'parameterCategory',parcat,iret)
171 call grib_check(iret,gribfunction,griberrormsg)
172 call grib_get_int(igrib,
'parameterNumber',parnum,iret)
173 call grib_check(iret,gribfunction,griberrormsg)
174 call grib_get_int(igrib,
'typeOfFirstFixedSurface',typsurf,iret)
175 call grib_check(iret,gribfunction,griberrormsg)
176 call grib_get_int(igrib,
'scaledValueOfFirstFixedSurface', &
178 call grib_check(iret,gribfunction,griberrormsg)
184 if ((parcat.eq.2).and.(parnum.eq.2).and.(typsurf.eq.100))
then
188 elseif ((parcat.eq.3).and.(parnum.eq.5).and.(typsurf.eq.1))
then
192 elseif ((parcat.eq.0).and.(parnum.eq.0).and.(typsurf.eq.1) &
193 .and.(discipl.eq.2))
then
199 if (isec1(6).ne.-1)
then
201 call grib_get_real4_array(igrib,
'values',zsec4,iret)
202 call grib_check(iret,gribfunction,griberrormsg)
211 call grib_get_int(igrib,
'numberOfPointsAlongAParallel', &
213 call grib_check(iret,gribfunction,griberrormsg)
214 call grib_get_int(igrib,
'numberOfPointsAlongAMeridian', &
216 call grib_check(iret,gribfunction,griberrormsg)
217 call grib_get_real8(igrib,
'longitudeOfFirstGridPointInDegrees', &
219 call grib_check(iret,gribfunction,griberrormsg)
220 call grib_get_real8(igrib,
'longitudeOfLastGridPointInDegrees', &
222 call grib_check(iret,gribfunction,griberrormsg)
223 call grib_get_real8(igrib,
'latitudeOfLastGridPointInDegrees', &
225 call grib_check(iret,gribfunction,griberrormsg)
226 call grib_get_real8(igrib,
'latitudeOfFirstGridPointInDegrees', &
228 call grib_check(iret,gribfunction,griberrormsg)
237 if((abs(xaux1).lt.eps).and.(xaux2.ge.359))
then
239 xaux2=-179.0+360.-360./
real(nxfield)
241 if (xaux1.gt.180) xaux1=xaux1-360.0
242 if (xaux2.gt.180) xaux2=xaux2-360.0
243 if (xaux1.lt.-180) xaux1=xaux1+360.0
244 if (xaux2.lt.-180) xaux2=xaux2+360.0
245 if (xaux2.lt.xaux1) xaux2=xaux2+360.
248 dx=(xaux2-xaux1)/
real(nxfield-1)
249 dy=(yaux2-yaux1)/
real(ny-1)
250 dxconst=180./(dx*r_earth*pi)
251 dyconst=180./(dy*r_earth*pi)
260 xauxa=abs(xaux2+dx-360.-xaux1)
261 if (xauxa.lt.0.001)
then
264 if (abs(nxshift).ge.nx) &
265 stop
'nxshift in file par_mod is too large'
266 xlon0=xlon0+
real(nxshift)*dx
271 stop
'nxshift (par_mod) must be zero for non-global domain'
275 if (xlon0.gt.180.) xlon0=xlon0-360.
277 if (xglobal.and.xauxa.lt.0.001)
then
281 sizesouth=6.*(switchsouth+90.)/dy
282 call
stlmbr(southpolemap,-90.,0.)
283 call
stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
284 sizesouth,switchsouth,180.)
285 switchsouthg=(switchsouth-ylat0)/dy
291 if (xglobal.and.xauxa.lt.0.001)
then
295 sizenorth=6.*(90.-switchnorth)/dy
296 call
stlmbr(northpolemap,90.,0.)
297 call
stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
298 sizenorth,switchnorth,180.)
299 switchnorthg=(switchnorth-ylat0)/dy
306 if (nxshift.lt.0) stop
'nxshift (par_mod) must not be negative'
307 if (nxshift.ge.nxfield) stop
'nxshift (par_mod) too large'
312 if((isec1(6).eq.33).and.(isec1(7).eq.100))
then
314 pres(iumax)=
real(isec1(8))*100.0
330 if((isec1(6).eq.007).and.(isec1(7).eq.001))
then
333 help=zsec4(nxfield*(ny-jy-1)+ix+1)
336 excessoro(i179+ix,jy)=0.0
339 excessoro(ix-i181,jy)=0.0
348 if((isec1(6).eq.081).and.(isec1(7).eq.001))
then
351 help=zsec4(nxfield*(ny-jy-1)+ix+1)
361 call grib_release(igrib)
370 call grib_close_file(ifile)
377 if (nx.gt.nxmax)
then
378 write(*,*)
'FLEXPART error: Too many grid points in x direction.'
379 write(*,*)
'Reduce resolution of wind fields.'
380 write(*,*)
'Or change parameter settings in file par_mod.'
385 if (ny.gt.nymax)
then
386 write(*,*)
'FLEXPART error: Too many grid points in y direction.'
387 write(*,*)
'Reduce resolution of wind fields.'
388 write(*,*)
'Or change parameter settings in file par_mod.'
393 if (nuvz.gt.nuvzmax)
then
394 write(*,*)
'FLEXPART error: Too many u,v grid points in z '// &
396 write(*,*)
'Reduce resolution of wind fields.'
397 write(*,*)
'Or change parameter settings in file par_mod.'
398 write(*,*) nuvz,nuvzmax
402 if (nwz.gt.nwzmax)
then
403 write(*,*)
'FLEXPART error: Too many w grid points in z '// &
405 write(*,*)
'Reduce resolution of wind fields.'
406 write(*,*)
'Or change parameter settings in file par_mod.'
407 write(*,*) nwz,nwzmax
425 write(*,
'(a,2i7)')
'# of vertical levels in NCEP data: ', &
428 write(*,
'(a)')
'Mother domain:'
429 write(*,
'(a,f10.2,a1,f10.2,a,f10.2)')
' Longitude range: ', &
430 xlon0,
' to ',xlon0+(nx-1)*dx,
' Grid distance: ',dx
431 write(*,
'(a,f10.2,a1,f10.2,a,f10.2)')
' Latitude range: ', &
432 ylat0,
' to ',ylat0+(ny-1)*dy,
' Grid distance: ',dy
443 k=nlev_ec+1+numskip+i
444 akm_usort(nwz-i+1)=pres(nwz-i+1)
452 if (akm_usort(1).gt.akm_usort(2))
then
455 akm(i)=akm_usort(nwz-i+1)
482 if (nz.gt.nzmax) stop
'nzmax too small'
506 pint=akz(i)+bkz(i)*101325.
507 if (pint.lt.5000.) goto 96
510 if (nconvlev.gt.nconvlevmax-1)
then
511 nconvlev=nconvlevmax-1
512 write(*,*)
'Attention, convection only calculated up to ', &
513 akz(nconvlev)+bkz(nconvlev)*1013.25,
' hPa'
519 write(*,*)
' ###########################################'// &
521 write(*,*)
' TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
522 write(*,*)
' CAN NOT OPEN INPUT DATA FILE '//wfname(ifn)
523 write(*,*)
' ###########################################'// &
526 write(*,
'(a)')
'!!! PLEASE INSERT A NEW CD-ROM AND !!!'
527 write(*,
'(a)')
'!!! PRESS ANY KEY TO CONTINUE... !!!'
528 write(*,
'(a)')
'!!! ...OR TERMINATE FLEXPART PRESSING!!!'
529 write(*,
'(a)')
'!!! THE "X" KEY... !!!'
subroutine, public stcm2p(strcmp, x1, y1, xlat1, xlong1, x2, y2, xlat2, xlong2)
subroutine, public stlmbr(strcmp, tnglat, xlong)
subroutine shift_field_0(field, nxf, nyf)