66 real(kind=4) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
67 real(kind=4) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
68 real(kind=4) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
91 integer :: gribver,parcat,parnum,typsurf,valsurf,discipl
92 #if defined WITH_CTBTO_PATCHES
126 integer :: isec1(56),isec2(22+nxmax+nymax)
127 real(kind=4) :: zsec4(jpunp)
128 real(kind=4) :: xaux,yaux,xaux0,yaux0
129 #if defined WITH_CTBTO_PATCHES
130 real(kind=4) :: conversion_factor
132 real(kind=8) :: xauxin,yauxin
133 real(kind=4) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
134 real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
135 logical :: hflswitch,strswitch
144 integer :: i,j,k,levdiff2,ifield,iumax,iwmax
149 real,
parameter :: eps=1.e-4
151 character(len=24) :: griberrormsg =
'Error reading grib file'
152 character(len=20) :: gribfunction =
'readwind'
192 levdiff2=nlev_ec-nwz+1
199 5 call grib_open_file(ifile,path(3)(1:length(3)) &
200 //trim(wfname(indj)),
'r',iret)
201 if (iret.ne.grib_success)
then
213 call grib_new_from_file(ifile,igrib,iret)
214 if (iret.eq.grib_end_of_file)
then
216 elseif (iret.ne.grib_success)
then
220 #if defined WITH_CTBTO_PATCHES
221 conversion_factor=1.0
225 call grib_get_int(igrib,
'editionNumber',gribver,iret)
226 call grib_check(iret,gribfunction,griberrormsg)
228 if (gribver.eq.1)
then
232 call grib_get_int(igrib,
'indicatorOfParameter',isec1(6),iret)
233 call grib_check(iret,gribfunction,griberrormsg)
234 call grib_get_int(igrib,
'level',isec1(8),iret)
235 call grib_check(iret,gribfunction,griberrormsg)
238 #if defined WITH_CTBTO_PATCHES
239 if (isec1(6).eq.77.or.isec1(6).eq.87)
then
241 if (isec1(6).eq.77)
then
248 #if defined WITH_CTBTO_PATCHES
249 call
grib2check(igrib,isec1,conversion_factor)
253 call grib_get_int(igrib,
'discipline',discipl,iret)
254 call grib_check(iret,gribfunction,griberrormsg)
255 call grib_get_int(igrib,
'parameterCategory',parcat,iret)
256 call grib_check(iret,gribfunction,griberrormsg)
257 call grib_get_int(igrib,
'parameterNumber',parnum,iret)
258 call grib_check(iret,gribfunction,griberrormsg)
259 call grib_get_int(igrib,
'typeOfFirstFixedSurface',typsurf,iret)
260 call grib_check(iret,gribfunction,griberrormsg)
261 call grib_get_int(igrib,
'level',valsurf,iret)
262 call grib_check(iret,gribfunction,griberrormsg)
271 if ((parcat.eq.0).and.(parnum.eq.0).and.(typsurf.eq.105))
then
273 elseif ((parcat.eq.2).and.(parnum.eq.2).and.(typsurf.eq.105))
then
275 elseif ((parcat.eq.2).and.(parnum.eq.3).and.(typsurf.eq.105))
then
277 elseif ((parcat.eq.1).and.(parnum.eq.0).and.(typsurf.eq.105))
then
279 elseif ((parcat.eq.3).and.(parnum.eq.0).and.(typsurf.eq.1))
then
281 elseif ((parcat.eq.2).and.(parnum.eq.32))
then
283 elseif ((parcat.eq.3).and.(parnum.eq.0).and.(typsurf.eq.101))
then
285 elseif ((parcat.eq.2).and.(parnum.eq.2).and.(typsurf.eq.103))
then
287 elseif ((parcat.eq.2).and.(parnum.eq.3).and.(typsurf.eq.103))
then
289 elseif ((parcat.eq.0).and.(parnum.eq.0).and.(typsurf.eq.103))
then
291 elseif ((parcat.eq.0).and.(parnum.eq.6).and.(typsurf.eq.103))
then
293 elseif ((parcat.eq.1).and.(parnum.eq.11).and.(typsurf.eq.1))
then
295 elseif ((parcat.eq.6).and.(parnum.eq.1))
then
297 elseif ((parcat.eq.1).and.(parnum.eq.9))
then
299 elseif ((parcat.eq.1).and.(parnum.eq.10))
then
301 elseif ((parcat.eq.0).and.(parnum.eq.11).and.(typsurf.eq.1))
then
303 elseif ((parcat.eq.4).and.(parnum.eq.9).and.(typsurf.eq.1))
then
305 elseif ((parcat.eq.2).and.(parnum.eq.17))
then
307 elseif ((parcat.eq.2).and.(parnum.eq.18))
then
309 elseif ((parcat.eq.3).and.(parnum.eq.4))
then
311 elseif ((parcat.eq.3).and.(parnum.eq.7))
then
313 elseif ((discipl.eq.2).and.(parcat.eq.0).and.(parnum.eq.0).and. &
317 print*,
'***ERROR: undefined GRiB2 message found!',discipl, &
318 parcat,parnum,typsurf
326 if (isec1(6).ne.-1)
then
327 call grib_get_real4_array(igrib,
'values',zsec4,iret)
328 #if defined WITH_CTBTO_PATCHES
329 zsec4=zsec4/conversion_factor
331 call grib_check(iret,gribfunction,griberrormsg)
335 if (ifield.eq.1)
then
336 call grib_get_int(igrib,
'numberOfPointsAlongAParallel', &
338 call grib_check(iret,gribfunction,griberrormsg)
339 call grib_get_int(igrib,
'numberOfPointsAlongAMeridian', &
341 call grib_check(iret,gribfunction,griberrormsg)
342 call grib_get_int(igrib,
'numberOfVerticalCoordinateValues', &
344 call grib_check(iret,gribfunction,griberrormsg)
346 if(isec2(2).ne.nxfield) stop
'READWIND: NX NOT CONSISTENT'
347 if(isec2(3).ne.ny) stop
'READWIND: NY NOT CONSISTENT'
348 if(isec2(12)/2-1.ne.nlev_ec) &
349 stop
'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
353 #if defined WITH_CTBTO_PATCHES
354 if (gotgrid.eq.0)
then
356 if ((gribver.eq.1).and.(gotgrid.eq.0))
then
358 call grib_get_real8(igrib,
'longitudeOfFirstGridPointInDegrees', &
360 call grib_check(iret,gribfunction,griberrormsg)
361 call grib_get_real8(igrib,
'latitudeOfLastGridPointInDegrees', &
363 call grib_check(iret,gribfunction,griberrormsg)
364 xaux=xauxin+
real(nxshift)*dx
368 if(xaux.lt.0.) xaux=xaux+360.
369 if(yaux.lt.0.) yaux=yaux+360.
370 if(xaux0.lt.0.) xaux0=xaux0+360.
371 if(yaux0.lt.0.) yaux0=yaux0+360.
372 if(abs(xaux-xaux0).gt.eps) &
373 stop
'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
374 if(abs(yaux-yaux0).gt.eps) &
375 stop
'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
382 if(isec1(6).eq.130) tth(i,j,nlev_ec-k+2,n)= &
383 zsec4(nxfield*(ny-j-1)+i+1)
384 if(isec1(6).eq.131) uuh(i,j,nlev_ec-k+2)= &
385 zsec4(nxfield*(ny-j-1)+i+1)
386 if(isec1(6).eq.132) vvh(i,j,nlev_ec-k+2)= &
387 zsec4(nxfield*(ny-j-1)+i+1)
388 if(isec1(6).eq.133)
then
389 qvh(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
390 if (qvh(i,j,nlev_ec-k+2,n) .lt. 0.) &
391 qvh(i,j,nlev_ec-k+2,n) = 0.
395 if(isec1(6).eq.134) ps(i,j,1,n)= &
396 zsec4(nxfield*(ny-j-1)+i+1)
398 if(isec1(6).eq.135) wwh(i,j,nlev_ec-k+1)= &
399 zsec4(nxfield*(ny-j-1)+i+1)
400 if(isec1(6).eq.141) sd(i,j,1,n)= &
401 zsec4(nxfield*(ny-j-1)+i+1)
402 if(isec1(6).eq.151) msl(i,j,1,n)= &
403 zsec4(nxfield*(ny-j-1)+i+1)
404 if(isec1(6).eq.164) tcc(i,j,1,n)= &
405 zsec4(nxfield*(ny-j-1)+i+1)
406 if(isec1(6).eq.165) u10(i,j,1,n)= &
407 zsec4(nxfield*(ny-j-1)+i+1)
408 if(isec1(6).eq.166) v10(i,j,1,n)= &
409 zsec4(nxfield*(ny-j-1)+i+1)
410 if(isec1(6).eq.167) tt2(i,j,1,n)= &
411 zsec4(nxfield*(ny-j-1)+i+1)
412 if(isec1(6).eq.168) td2(i,j,1,n)= &
413 zsec4(nxfield*(ny-j-1)+i+1)
414 if(isec1(6).eq.142)
then
415 lsprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
416 if (lsprec(i,j,1,n).lt.0.) lsprec(i,j,1,n)=0.
418 if(isec1(6).eq.143)
then
419 convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
420 if (convprec(i,j,1,n).lt.0.) convprec(i,j,1,n)=0.
422 if(isec1(6).eq.146) sshf(i,j,1,n)= &
423 zsec4(nxfield*(ny-j-1)+i+1)
424 if((isec1(6).eq.146).and.(zsec4(nxfield*(ny-j-1)+i+1).ne.0.)) &
426 if(isec1(6).eq.176)
then
427 ssr(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
428 if (ssr(i,j,1,n).lt.0.) ssr(i,j,1,n)=0.
430 if(isec1(6).eq.180) ewss(i,j)= &
431 zsec4(nxfield*(ny-j-1)+i+1)
432 if(isec1(6).eq.181) nsss(i,j)= &
433 zsec4(nxfield*(ny-j-1)+i+1)
434 if(((isec1(6).eq.180).or.(isec1(6).eq.181)).and. &
435 (zsec4(nxfield*(ny-j-1)+i+1).ne.0.)) strswitch=.true.
437 if(isec1(6).eq.129) oro(i,j)= &
438 zsec4(nxfield*(ny-j-1)+i+1)/ga
439 if(isec1(6).eq.160) excessoro(i,j)= &
440 zsec4(nxfield*(ny-j-1)+i+1)
441 if(isec1(6).eq.172) lsm(i,j)= &
442 zsec4(nxfield*(ny-j-1)+i+1)
443 if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
444 if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
449 call grib_release(igrib)
455 50 call grib_close_file(ifile)
458 if (gotgrid.eq.0)
then
459 print*,
'***ERROR: input file needs to contain GRiB1 formatted'// &
464 if(levdiff2.eq.0)
then
468 wwh(i,j,nlev_ec+1)=0.
504 surfstr(i,j,1,n)=sqrt(ewss(i,j)**2+nsss(i,j)**2)
508 if ((.not.hflswitch).or.(.not.strswitch))
then
509 write(*,*)
'WARNING: No flux data contained in GRIB file ', &
520 plev1=akz(3)+bkz(3)*ps(i,j,1,n)
521 pmean=0.5*(ps(i,j,1,n)+plev1)
522 tv=tth(i,j,3,n)*(1.+0.61*qvh(i,j,3,n))
523 fu=-r_air*tv/ga/pmean
524 hlev1=fu*(plev1-ps(i,j,1,n))
525 ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2)
526 fflev1=sqrt(uuh(i,j,3)**2+vvh(i,j,3)**2)
528 tt2(i,j,1,n),tth(i,j,3,n),ff10m,fflev1, &
529 surfstr(i,j,1,n),sshf(i,j,1,n))
530 if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
531 if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
545 uuh(i,j,1)=u10(i,j,1,n)
546 vvh(i,j,1)=v10(i,j,1,n)
547 qvh(i,j,1,n)=qvh(i,j,2,n)
548 tth(i,j,1,n)=tt2(i,j,1,n)
552 if(iumax.ne.nuvz-1) stop
'READWIND: NUVZ NOT CONSISTENT'
553 if(iwmax.ne.nwz) stop
'READWIND: NWZ NOT CONSISTENT'
556 888
write(*,*)
' #### FLEXPART MODEL ERROR! WINDFIELD #### '
557 write(*,*)
' #### ',wfname(indj),
' #### '
558 write(*,*)
' #### IS NOT GRIB FORMAT !!! #### '
559 stop
'Execution terminated'
560 999
write(*,*)
' #### FLEXPART MODEL ERROR! WINDFIELD #### '
561 write(*,*)
' #### ',wfname(indj),
' #### '
562 write(*,*)
' #### CANNOT BE OPENED !!! #### '
563 stop
'Execution terminated'
subroutine readwind_ecmwf(indj, n, uuh, vvh, wwh)
subroutine shift_field(field, nxf, nyf, nzfmax, nzf, nmax, n)
subroutine pbl_profile(ps, td2m, zml1, t2m, tml1, u10m, uml1, stress, hf)
subroutine shift_field_0(field, nxf, nyf)
subroutine grib2check(igrib, isec1, conversion_factor)