50 integer :: gribver,parcat,parnum,typsurf,valsurf,discipl
51 #if defined WITH_CTBTO_PATCHES
56 integer :: i,j,k,l,ifn,ifield,iumax,iwmax,numskip,nlev_ecn
58 real :: akmn(nwzmax),bkmn(nwzmax),akzn(nuvzmax),bkzn(nuvzmax)
59 real(kind=4) :: xaux1,xaux2,yaux1,yaux2
60 real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
70 integer :: isec1(56),isec2(22+nxmaxn+nymaxn)
71 real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
72 #if defined WITH_CTBTO_PATCHES
73 real(kind=4) ::conversion_factor
77 character(len=24) :: griberrormsg =
'Error reading grib file'
78 character(len=20) :: gribfunction =
'gridcheck_nests'
86 character(LEN=255),
parameter :: vtable_ecmwf_grib1_path = &
87 "Vtable_ecmwf_grib1", &
88 vtable_ecmwf_grib2_path = &
89 "Vtable_ecmwf_grib2", &
90 vtable_ecmwf_grib1_2_path = &
91 "Vtable_ecmwf_grib1_2"
93 integer :: gribfile_type
94 integer :: current_grib_level
95 character(len=255) :: gribfile_name
96 character(len=255) :: vtable_path
97 character(len=15) :: fpname
116 if(ideltas.gt.0)
then
126 gribfile_name = path(3)(1:length(3))//trim(wfname(ifn))
127 print *,
'gribfile_name: ', gribfile_name
131 print *,
'gribfile_type: ', gribfile_type
133 if (gribfile_type .eq. vtable_gribfile_type_ecmwf_grib1)
then
134 vtable_path = vtable_ecmwf_grib1_path
135 else if (gribfile_type .eq. vtable_gribfile_type_ecmwf_grib2)
then
136 vtable_path = vtable_ecmwf_grib2_path
137 else if (gribfile_type .eq. vtable_gribfile_type_ecmwf_grib1_2)
then
138 vtable_path = vtable_ecmwf_grib1_2_path
140 print *,
'Unsupported gribfile_type: ', gribfile_type
146 print *,
'Loading Vtable: ', vtable_path
148 print *,
'Vtable Initialized: ', my_vtable%initialized
149 print *,
'Vtable num_entries: ', my_vtable%num_entries
164 5 call grib_open_file(ifile,path(numpath+2*(l-1)+1) &
165 (1:length(numpath+2*(l-1)+1))//trim(wfnamen(l,ifn)),
'r',iret)
166 if (iret.ne.grib_success)
then
179 call grib_new_from_file(ifile,igrib,iret)
180 if (iret.eq.grib_end_of_file)
then
182 elseif (iret.ne.grib_success)
then
202 #if defined WITH_CTBTO_PATCHES
203 conversion_factor=1.0
207 call grib_get_int(igrib,
'editionNumber',gribver,iret)
208 call grib_check(iret,gribfunction,griberrormsg)
210 if (gribver.eq.1)
then
216 call grib_get_int(igrib,
'level',current_grib_level,iret)
217 call grib_check(iret,gribfunction,griberrormsg)
221 #if defined WITH_CTBTO_PATCHES
222 call
grib2check(igrib, fpname, conversion_factor)
225 call grib_get_int(igrib,
'level', current_grib_level, iret)
226 call grib_check(iret, gribfunction, griberrormsg)
233 if (trim(fpname) .ne.
'None')
then
234 call grib_get_real4_array(igrib,
'values',zsec4,iret)
235 #if defined WITH_CTBTO_PATCHES
236 zsec4=zsec4/conversion_factor
238 call grib_check(iret,gribfunction,griberrormsg)
242 if (ifield.eq.1)
then
243 call grib_get_int(igrib,
'numberOfPointsAlongAParallel', &
245 call grib_check(iret,gribfunction,griberrormsg)
246 call grib_get_int(igrib,
'numberOfPointsAlongAMeridian', &
248 call grib_check(iret,gribfunction,griberrormsg)
249 call grib_get_int(igrib,
'numberOfVerticalCoordinateValues', &
251 call grib_check(iret,gribfunction,griberrormsg)
253 call grib_get_real4_array(igrib,
'pv',zsec2,iret)
254 call grib_check(iret,gribfunction,griberrormsg)
258 nlev_ecn=isec2(12)/2-1
262 #if defined WITH_CTBTO_PATCHES
267 call grib_get_real8(igrib,
'longitudeOfFirstGridPointInDegrees', &
269 call grib_check(iret,gribfunction,griberrormsg)
270 call grib_get_real8(igrib,
'longitudeOfLastGridPointInDegrees', &
272 call grib_check(iret,gribfunction,griberrormsg)
273 call grib_get_real8(igrib,
'latitudeOfLastGridPointInDegrees', &
275 call grib_check(iret,gribfunction,griberrormsg)
276 call grib_get_real8(igrib,
'latitudeOfFirstGridPointInDegrees', &
278 call grib_check(iret,gribfunction,griberrormsg)
283 if(xaux1.gt.180) xaux1=xaux1-360.0
284 if(xaux2.gt.180) xaux2=xaux2-360.0
285 if(xaux1.lt.-180) xaux1=xaux1+360.0
286 if(xaux2.lt.-180) xaux2=xaux2+360.0
287 if (xaux2.lt.xaux1) xaux2=xaux2+360.
290 dxn(l)=(xaux2-xaux1)/
real(nxn(l)-1)
291 dyn(l)=(yaux2-yaux1)/
real(nyn(l)-1)
296 if(trim(fpname) .eq.
'UU') iumax=max(iumax,nlev_ec-k+1)
297 if(trim(fpname) .eq.
'ETADOT') iwmax=max(iwmax,nlev_ec-k+1)
299 if(trim(fpname) .eq.
'ORO')
then
302 oron(i,j,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
306 if(trim(fpname) .eq.
'LSM')
then
309 lsmn(i,j,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
313 if(trim(fpname) .eq.
'EXCESSORO')
then
316 excessoron(i,j,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
321 call grib_release(igrib)
327 30 call grib_close_file(ifile)
330 if (gotgrib.eq.0)
then
331 print*,
'***ERROR: input file needs to contain GRiB1 formatted'// &
338 if(nuvzn.eq.nlev_ec) nwzn=nlev_ecn+1
340 if (nxn(l).gt.nxmaxn)
then
341 write(*,*)
'FLEXPART error: Too many grid points in x direction.'
342 write(*,*)
'Reduce resolution of wind fields (file GRIDSPEC)'
343 write(*,*)
'for nesting level ',l
344 write(*,*)
'Or change parameter settings in file par_mod.'
345 write(*,*) nxn(l),nxmaxn
349 if (nyn(l).gt.nymaxn)
then
350 write(*,*)
'FLEXPART error: Too many grid points in y direction.'
351 write(*,*)
'Reduce resolution of wind fields (file GRIDSPEC)'
352 write(*,*)
'for nesting level ',l
353 write(*,*)
'Or change parameter settings in file par_mod.'
354 write(*,*) nyn(l),nymaxn
358 if ((nuvzn.gt.nuvzmax).or.(nwzn.gt.nwzmax))
then
359 write(*,*)
'FLEXPART error: Nested wind fields have too many'// &
361 write(*,*)
'Problem was encountered for nesting level ',l
369 write(*,
'(a,i2)')
'Nested domain #: ',l
370 write(*,
'(a,f10.2,a1,f10.2,a,f10.2)')
' Longitude range: ', &
371 xlon0n(l),
' to ',xlon0n(l)+(nxn(l)-1)*dxn(l), &
372 ' Grid distance: ',dxn(l)
373 write(*,
'(a,f10.2,a1,f10.2,a,f10.2)')
' Latitude range: ', &
374 ylat0n(l),
' to ',ylat0n(l)+(nyn(l)-1)*dyn(l), &
375 ' Grid distance: ',dyn(l)
391 xaux2=xlon0n(l)+
real(nxn(l)-1)*dxn(l)
393 yaux2=ylat0n(l)+
real(nyn(l)-1)*dyn(l)
395 xln(l)=(xaux1-xlon0)/dx
396 xrn(l)=(xaux2-xlon0)/dx
397 yln(l)=(yaux1-ylat0)/dy
398 yrn(l)=(yaux2-ylat0)/dy
401 if ((xln(l).lt.0.).or.(yln(l).lt.0.).or. &
402 (xrn(l).gt.
real(nxmin1)).or.(yrn(l).gt.
real(nymin1))) then
403 write(*,*)
'Nested domain does not fit into mother domain'
404 write(*,*)
'For global mother domain fields, you can shift'
405 write(*,*)
'shift the mother domain into x-direction'
406 write(*,*)
'by setting nxshift (file par_mod) to a'
407 write(*,*)
'positive value. Execution is terminated.'
415 numskip=nlev_ecn-nuvzn
418 k=nlev_ecn+1+numskip+i
419 akmn(nwzn-i+1)=zsec2(j)
420 bkmn(nwzn-i+1)=zsec2(k)
435 akzn(i+1)=0.5*(akmn(i+1)+akmn(i))
436 bkzn(i+1)=0.5*(bkmn(i+1)+bkmn(i))
446 if ((akzn(i).ne.akz(i)).or.(bkzn(i).ne.bkz(i)))
then
447 write(*,*)
'FLEXPART error: The wind fields of nesting level',l
448 write(*,*)
'are not consistent with the mother domain:'
449 write(*,*)
'Differences in vertical levels detected.'
455 if ((akmn(i).ne.akm(i)).or.(bkmn(i).ne.bkm(i)))
then
456 write(*,*)
'FLEXPART error: The wind fields of nesting level',l
457 write(*,*)
'are not consistent with the mother domain:'
458 write(*,*)
'Differences in vertical levels detected.'
468 write(*,*)
' ###########################################'// &
470 write(*,*)
' FLEXPART SUBROUTINE GRIDCHECK:'
471 write(*,*)
' CAN NOT OPEN INPUT DATA FILE '//wfnamen(l,ifn)
472 write(*,*)
' FOR NESTING LEVEL ',k
473 write(*,*)
' ###########################################'// &
character(len=15) function, public vtable_get_fpname(igrib, vtable_object)
subroutine grib2check(igrib, fpname, conversion_factor)
subroutine gridcheck_nests
subroutine, public vtable_load_by_name(vtable_name, the_vtable_data)
integer function, public vtable_detect_gribfile_type(gribfilename)