67 real(kind=4) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
68 real(kind=4) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
69 real(kind=4) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
92 integer :: gribver,parcat,parnum,typsurf,valsurf,discipl
93 #if defined WITH_CTBTO_PATCHES
104 character(LEN=255),
parameter :: vtable_ecmwf_grib1_path = &
105 "Vtable_ecmwf_grib1", &
106 vtable_ecmwf_grib2_path = &
107 "Vtable_ecmwf_grib2", &
108 vtable_ecmwf_grib1_2_path = &
109 "Vtable_ecmwf_grib1_2"
111 integer :: gribfile_type
112 integer :: current_grib_level
113 character(len=255) :: gribfile_name
114 character(len=255) :: vtable_path
115 character(len=15) :: fpname
152 integer :: isec1(56),isec2(22+nxmax+nymax)
153 real(kind=4) :: zsec4(jpunp)
154 real(kind=4) :: xaux,yaux,xaux0,yaux0
155 #if defined WITH_CTBTO_PATCHES
156 real(kind=4) :: conversion_factor
158 real(kind=8) :: xauxin,yauxin
159 real(kind=4) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
160 real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
161 logical :: hflswitch,strswitch
170 integer :: i,j,k,levdiff2,ifield,iumax,iwmax
175 real,
parameter :: eps=1.e-4
177 character(len=24) :: griberrormsg =
'Error reading grib file'
178 character(len=20) :: gribfunction =
'readwind'
218 levdiff2=nlev_ec-nwz+1
225 gribfile_name = path(3)(1:length(3))//trim(wfname(indj))
226 print *,
'gribfile_name: ', gribfile_name
230 print *,
'gribfile_type: ', gribfile_type
232 if (gribfile_type .eq. vtable_gribfile_type_ecmwf_grib1)
then
233 vtable_path = vtable_ecmwf_grib1_path
234 else if (gribfile_type .eq. vtable_gribfile_type_ecmwf_grib2)
then
235 vtable_path = vtable_ecmwf_grib2_path
236 else if (gribfile_type .eq. vtable_gribfile_type_ecmwf_grib1_2)
then
237 vtable_path = vtable_ecmwf_grib1_2_path
239 print *,
'Unsupported gribfile_type: ', gribfile_type
245 print *,
'Loading Vtable: ', vtable_path
247 print *,
'Vtable Initialized: ', my_vtable%initialized
248 print *,
'Vtable num_entries: ', my_vtable%num_entries
259 5 call grib_open_file(ifile,path(3)(1:length(3)) &
260 //trim(wfname(indj)),
'r',iret)
261 if (iret.ne.grib_success)
then
273 call grib_new_from_file(ifile,igrib,iret)
274 if (iret.eq.grib_end_of_file)
then
276 elseif (iret.ne.grib_success)
then
290 #if defined WITH_CTBTO_PATCHES
291 conversion_factor=1.0
295 call grib_get_int(igrib,
'editionNumber',gribver,iret)
296 call grib_check(iret,gribfunction,griberrormsg)
299 if (gribver.eq.1)
then
301 call grib_get_int(igrib,
'level', current_grib_level,iret)
302 call grib_check(iret,gribfunction,griberrormsg)
307 #if defined WITH_CTBTO_PATCHES
309 call
grib2check(igrib, fpname, conversion_factor)
312 call grib_get_int(igrib,
'level', current_grib_level,iret)
313 call grib_check(iret,gribfunction,griberrormsg)
321 if (trim(fpname) .ne.
'None')
then
322 call grib_get_real4_array(igrib,
'values',zsec4,iret)
323 #if defined WITH_CTBTO_PATCHES
324 zsec4=zsec4/conversion_factor
326 call grib_check(iret,gribfunction,griberrormsg)
330 if (ifield.eq.1)
then
331 call grib_get_int(igrib,
'numberOfPointsAlongAParallel', &
333 call grib_check(iret,gribfunction,griberrormsg)
334 call grib_get_int(igrib,
'numberOfPointsAlongAMeridian', &
336 call grib_check(iret,gribfunction,griberrormsg)
337 call grib_get_int(igrib,
'numberOfVerticalCoordinateValues', &
339 call grib_check(iret,gribfunction,griberrormsg)
341 if(isec2(2).ne.nxfield) stop
'READWIND: NX NOT CONSISTENT'
342 if(isec2(3).ne.ny) stop
'READWIND: NY NOT CONSISTENT'
343 if(isec2(12)/2-1.ne.nlev_ec) &
344 stop
'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
348 #if defined WITH_CTBTO_PATCHES
353 call grib_get_real8(igrib,
'longitudeOfFirstGridPointInDegrees', &
355 call grib_check(iret,gribfunction,griberrormsg)
356 call grib_get_real8(igrib,
'latitudeOfLastGridPointInDegrees', &
358 call grib_check(iret,gribfunction,griberrormsg)
359 xaux=xauxin+
real(nxshift)*dx
363 if(xaux.lt.0.) xaux=xaux+360.
364 if(yaux.lt.0.) yaux=yaux+360.
365 if(xaux0.lt.0.) xaux0=xaux0+360.
366 if(yaux0.lt.0.) yaux0=yaux0+360.
367 if(abs(xaux-xaux0).gt.eps) &
368 stop
'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
369 if(abs(yaux-yaux0).gt.eps) &
370 stop
'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
375 k = current_grib_level
378 if(trim(fpname) .eq.
'TT')
then
381 tth(i,j,nlev_ec-k+2,n) = zsec4(nxfield*(ny-j-1)+i+1)
388 if(trim(fpname) .eq.
'UU')
then
389 iumax=max(iumax,nlev_ec-k+1)
392 uuh(i,j,nlev_ec-k+2) = zsec4(nxfield*(ny-j-1)+i+1)
399 if(trim(fpname) .eq.
'VV')
then
402 vvh(i,j,nlev_ec-k+2) = zsec4(nxfield*(ny-j-1)+i+1)
408 if(trim(fpname) .eq.
'ETADOT')
then
409 iwmax=max(iwmax,nlev_ec-k+1)
412 wwh(i,j,nlev_ec-k+1) = zsec4(nxfield*(ny-j-1)+i+1)
418 if(trim(fpname) .eq.
'QV')
then
421 qvh(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
422 if (qvh(i,j,nlev_ec-k+2,n) .lt. 0.) qvh(i,j,nlev_ec-k+2,n) = 0.
430 if(trim(fpname) .eq.
'PS')
then
433 ps(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
439 if(trim(fpname) .eq.
'SD')
then
442 sd(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
448 if(trim(fpname) .eq.
'MSL')
then
451 msl(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
457 if(trim(fpname) .eq.
'TCC')
then
460 tcc(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
466 if(trim(fpname) .eq.
'U10')
then
469 u10(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
475 if(trim(fpname) .eq.
'V10')
then
478 v10(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
484 if(trim(fpname) .eq.
'T2')
then
487 tt2(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
493 if(trim(fpname) .eq.
'TD2')
then
496 td2(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
502 if(trim(fpname) .eq.
'LSPREC')
then
505 lsprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
506 if (lsprec(i,j,1,n).lt.0.) lsprec(i,j,1,n)=0.
512 if(trim(fpname) .eq.
'CONVPREC')
then
515 convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
516 if (convprec(i,j,1,n).lt.0.) convprec(i,j,1,n)=0.
522 if(trim(fpname) .eq.
'SHF')
then
525 sshf(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
526 if (zsec4(nxfield*(ny-j-1)+i+1).ne.0.) hflswitch=.true.
532 if(trim(fpname) .eq.
'SR')
then
535 ssr(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
536 if (ssr(i,j,1,n).lt.0.) ssr(i,j,1,n)=0.
542 if(trim(fpname) .eq.
'EWSS')
then
545 ewss(i,j)=zsec4(nxfield*(ny-j-1)+i+1)
546 if (zsec4(nxfield*(ny-j-1)+i+1).ne.0.) strswitch=.true.
552 if(trim(fpname) .eq.
'NSSS')
then
555 nsss(i,j)=zsec4(nxfield*(ny-j-1)+i+1)
556 if (zsec4(nxfield*(ny-j-1)+i+1).ne.0.) strswitch=.true.
562 if(trim(fpname) .eq.
'ORO')
then
565 oro(i,j)=zsec4(nxfield*(ny-j-1)+i+1)/ga
571 if(trim(fpname) .eq.
'EXCESSORO')
then
574 excessoro(i,j)=zsec4(nxfield*(ny-j-1)+i+1)
580 if(trim(fpname) .eq.
'LSM')
then
583 lsm(i,j)=zsec4(nxfield*(ny-j-1)+i+1)
589 call grib_release(igrib)
595 50 call grib_close_file(ifile)
598 if (gotgrid.eq.0)
then
599 print*,
'***ERROR: no fields found with correct first longitude'// &
604 if(levdiff2.eq.0)
then
608 wwh(i,j,nlev_ec+1)=0.
644 surfstr(i,j,1,n)=sqrt(ewss(i,j)**2+nsss(i,j)**2)
648 if ((.not.hflswitch).or.(.not.strswitch))
then
649 write(*,*)
'WARNING: No flux data contained in GRIB file ', &
660 plev1=akz(3)+bkz(3)*ps(i,j,1,n)
661 pmean=0.5*(ps(i,j,1,n)+plev1)
662 tv=tth(i,j,3,n)*(1.+0.61*qvh(i,j,3,n))
663 fu=-r_air*tv/ga/pmean
664 hlev1=fu*(plev1-ps(i,j,1,n))
665 ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2)
666 fflev1=sqrt(uuh(i,j,3)**2+vvh(i,j,3)**2)
668 tt2(i,j,1,n),tth(i,j,3,n),ff10m,fflev1, &
669 surfstr(i,j,1,n),sshf(i,j,1,n))
670 if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
671 if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
685 uuh(i,j,1)=u10(i,j,1,n)
686 vvh(i,j,1)=v10(i,j,1,n)
687 qvh(i,j,1,n)=qvh(i,j,2,n)
688 tth(i,j,1,n)=tt2(i,j,1,n)
692 if(iumax.ne.nuvz-1) stop
'READWIND: NUVZ NOT CONSISTENT'
693 if(iwmax.ne.nwz) stop
'READWIND: NWZ NOT CONSISTENT'
696 888
write(*,*)
' #### FLEXPART MODEL ERROR! WINDFIELD #### '
697 write(*,*)
' #### ',wfname(indj),
' #### '
698 write(*,*)
' #### IS NOT GRIB FORMAT !!! #### '
699 stop
'Execution terminated'
700 999
write(*,*)
' #### FLEXPART MODEL ERROR! WINDFIELD #### '
701 write(*,*)
' #### ',wfname(indj),
' #### '
702 write(*,*)
' #### CANNOT BE OPENED !!! #### '
703 stop
'Execution terminated'
character(len=15) function, public vtable_get_fpname(igrib, vtable_object)
subroutine readwind_ecmwf(indj, n, uuh, vvh, wwh)
subroutine shift_field(field, nxf, nyf, nzfmax, nzf, nmax, n)
subroutine grib2check(igrib, fpname, conversion_factor)
subroutine pbl_profile(ps, td2m, zml1, t2m, tml1, u10m, uml1, stress, hf)
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)