program create_firstguess !*********************************************************************** ! prgmmr: pondeca org: np20 date: 2004-10-13 ! ! abstract: ! generate wrf format first guess for surface analysis. input is ! binary file contructed from first guess GRIB1 file using 'wgrib' ! ! program history log: ! 2004-10-13 pondeca ! 2007-06-11 pondeca - rewrite to handle not just the ndfd conus ! grid but also alaska, hawaii, guam and ! puerto rico ndfd grids ! 2010-04-08 pondeca - set XICE to zero. it was set to landmask ! 2011-02-03 zhu - add changes for pblh, gust, vis ! 2012-02-09 zhu - add changes for tcamt and lcbas ! 2016-02-09 sflamp - add inpaint subroutine for howv ! ! notes: (1) current code removes the bitmap for a single first guess. ! (2) if blending is in effect and some of the blending guesses ! contain a bitmap, they are "regularized" by the bitmap free ! first guess ! !*********************************************************************** implicit none character(60) cgrid include 'param.incl' integer(4),parameter:: nt0=17! # of horiz. slabs from wgrib integer(4),parameter:: nsig=1 real(4),parameter:: gravity=9.81 real(4),parameter:: qmin=1.e-06 real(4),parameter:: nhwrfmax=10 real(4),parameter:: ngesblendmax=30 !the maximum number of fcsts that !can be blended together to make the first guess. real(4),parameter:: ntmax=200 !maximum # of horiz. slabs from wgrib real(4),parameter::epsilon=1.e-3 real(4),parameter:: cldch_bitmapval=1.e+21 integer(4) nt !read in from namelist integer(4) nx,ny integer(4) i,j,n integer(4) iyear,imonth,iday,ihour,iminute,isecond integer(4) ihwrfyear,ihwrfmonth,ihwrfday,ihwrfhour, & ihwrfminute,ihwrfsecond integer(4) ihdrbuf(512) integer(4),dimension(5):: iadate, iadate2 integer(4) lun integer(4) ngesblend real(4) ds real(4) rmax,rmin real(4) gesfcsthour(ngesblendmax) real(4) dmax(ntmax) !the prescribed width of the blend zone in grid units for each slab character(60) gesmodelname(30) character(60) gesbinaryname(30) character(3) gesmodel_cldch_reflev(30) character(1),parameter::blank1=' ' integer(4),allocatable,dimension(:,:)::ifield,ivgtyp,isltyp real(4),allocatable,dimension(:,:)::glat,glon,field,sst,tslb, & tsk,vegfra,smois,qgrid,tdgrid, & psfcgrid,ugrid,vgrid,landmask, & dx,dy,mapfac,phbgrid,tgrid,gust, & cldch,vis,pblh,tcamt,lcbas,hgtgrid, & wspd,mxtm,mitm,pmsl,howv,uwnd,vwnd, & sfcrough, & !for similarity theory presgrid1,presgrid2,tmpgrid1,tmpgrid2,& qgrid1,qgrid2,ugrid1,vgrid1,hgtgrid1,& tggrid,presgsfc real(4),allocatable,dimension(:,:)::bias_qgrid,bias_psfcgrid, & bias_ugrid,bias_vgrid,bias_tgrid,bias_tdgrid, & bias_gust,bias_vis,bias_pblh,bias_tcamt,bias_lcbas, & bias_wspd,bias_mxtm,bias_mitm, & bias_pmsl,bias_howv,bias_cldch, & bias_uwnd,bias_vwnd real(4),allocatable,dimension(:,:)::undef_cldch real(4),allocatable,dimension(:,:)::terrain character(60) cbitfile, cnobitfile , cmodelwithbitmap character(3) cbitfile_cldch_reflev,cnobitfile_cldch_reflev character(10) var character(60) slabname(200) real(4) pbiascor,tbiascor,qbiascor,ubiascor,vbiascor,tdbiascor,& gustbiascor,visbiascor,pblhbiascor,tcamtbiascor,lcbasbiascor, & wspdbiascor,mxtmbiascor,mitmbiascor,pmslbiascor,howvbiascor,cldchbiascor, & uwndbiascor,vwndbiascor,sfcr0 character(3) cldch_reflev !used when there is no bitmap removal nor blending logical lbiascor logical mkrjlists logical fexist logical hwrfblend logical lbitmap logical lgesblend logical gesblend_ready logical use_sfcr0 logical neutral_stability_windfact_2dvar logical use_similarity_2dvar logical present_sfcrough_slab namelist/gridname/cgrid,lbiascor,pbiascor,tbiascor,qbiascor,ubiascor,vbiascor, & tdbiascor,gustbiascor,visbiascor,pblhbiascor,tcamtbiascor,lcbasbiascor, & wspdbiascor,mxtmbiascor,mitmbiascor,pmslbiascor,howvbiascor,cldchbiascor, & uwndbiascor,vwndbiascor, & cldch_reflev,mkrjlists,hwrfblend,use_sfcr0,sfcr0, & neutral_stability_windfact_2dvar,use_similarity_2dvar namelist/timeinfo/iyear,imonth,iday,ihour, & iminute,isecond namelist/hwrfinfo/ihwrfyear,ihwrfmonth,ihwrfday,ihwrfhour, & ihwrfminute,ihwrfsecond namelist/gesblendinfo/lgesblend,ngesblend,gesfcsthour, & gesmodelname,gesbinaryname,gesmodel_cldch_reflev namelist/bitmapinfo/lbitmap,cbitfile,cnobitfile,cmodelwithbitmap, & cbitfile_cldch_reflev,cnobitfile_cldch_reflev namelist/width_of_blendzone/dmax data slabname(1)/'PSFC'/ data slabname(2)/'HGT'/ data slabname(3)/'T2m'/ data slabname(4)/'TD2m'/ data slabname(5)/'U10m'/ data slabname(6)/'V10m'/ data slabname(7)/'Q2m'/ data slabname(8)/'GUST10m'/ data slabname(9)/'CLDCH'/ data slabname(10)/'VIS'/ data slabname(11)/'HPBL'/ data slabname(12)/'MXTM'/ data slabname(13)/'MITM'/ data slabname(14)/'PMSL'/ data slabname(15)/'HOWV'/ data slabname(16)/'TCDC'/ data slabname(17)/'LCBAS'/ data slabname(18)/'SFCR'/ data slabname(19)/'PRESGRID1'/ data slabname(20)/'PRESGRID2'/ data slabname(21)/'TMPGRID1'/ data slabname(22)/'TMPGRID2'/ data slabname(23)/'QGRID1'/ data slabname(24)/'QGRID2'/ data slabname(25)/'UGRID1'/ data slabname(26)/'VGRID1'/ data slabname(27)/'HGTGRID1'/ data slabname(28)/'TGGRID'/ data slabname(29)/'PRESGSFC'/ data cgrid/'conus'/ data lbiascor/.false./ data pbiascor/0.0/ data tbiascor/0.0/ data qbiascor/0.0/ data ubiascor/0.0/ data vbiascor/0.0/ data tdbiascor/0.0/ data gustbiascor/0.0/ data visbiascor/0.0/ data pblhbiascor/0.0/ data tcamtbiascor/0.0/ data lcbasbiascor/0.0/ data wspdbiascor/0.0/ data mxtmbiascor/0.0/ data mitmbiascor/0.0/ data pmslbiascor/0.0/ data howvbiascor/0.0/ data cldchbiascor/0.0/ data cldch_reflev/'ASL'/ data uwndbiascor/0.0/ data vwndbiascor/0.0/ data sfcr0/0.5/ data mkrjlists/.false./ data iyear/2005/ data imonth/6/ data iday/17/ data ihour/12/ data iminute/0/ data isecond/0/ data hwrfblend/.false./ data ihwrfyear/-9999/ data ihwrfmonth/-99/ data ihwrfday/-99/ data ihwrfhour/-99/ data ihwrfminute/0/ data ihwrfsecond/0/ data lgesblend/.false./ data ngesblend/0/ data gesmodel_cldch_reflev/30*'ASL'/ data lbitmap/.false./ data cbitfile_cldch_reflev /'ASL'/ data cnobitfile_cldch_reflev /'ASL'/ data use_sfcr0/.true./ data neutral_stability_windfact_2dvar/.true./ data use_similarity_2dvar/.false./ data dmax/200*50./ ! !==================================================================== !==================================================================== !==> find out what grid this is !==================================================================== open (55,file='gridname_input',form='formatted') read(55,gridname) close(55) print*,'in create_firstguess: cgrid is =',trim(cgrid) print*,'in create_firstguess: lbiascor=', lbiascor print*,'in create_firstguess: pbiascor=', pbiascor print*,'in create_firstguess: tbiascor=', tbiascor print*,'in create_firstguess: qbiascor=', qbiascor print*,'in create_firstguess: ubiascor=', ubiascor print*,'in create_firstguess: vbiascor=', vbiascor print*,'in create_firstguess: tdbiascor=', tdbiascor print*,'in create_firstguess: gustbiascor=', gustbiascor print*,'in create_firstguess: visbiascor=', visbiascor print*,'in create_firstguess: pblhbiascor=', pblhbiascor print*,'in create_firstguess: tcamtbiascor=',tcamtbiascor print*,'in create_firstguess: lcbasbiascor=',lcbasbiascor print*,'in create_firstguess: wspdbiascor=', wspdbiascor print*,'in create_firstguess: mxtmbiascor=', mxtmbiascor print*,'in create_firstguess: mitmbiascor=', mitmbiascor print*,'in create_firstguess: pmslbiascor=', pmslbiascor print*,'in create_firstguess: howvbiascor=', howvbiascor print*,'in create_firstguess: cldchbiascor=',cldchbiascor print*,'in create_firstguess: cldch_reflev=',cldch_reflev print*,'in create_firstguess: uwndbiascor=' ,uwndbiascor print*,'in create_firstguess: vwndbiascor=' ,vwndbiascor print*,'in create_firstguess: mkrjlists=', mkrjlists print*,'in create_firstguess: hwrfblend=', hwrfblend print*,'in create_firstguess: use_sfcr0=', use_sfcr0 print*,'in create_firstguess: sfcr0=',sfcr0 print*,'in create_firstguess: neutral_stability_windfact_2dvar=',neutral_stability_windfact_2dvar print*,'in create_firstguess: use_similarity_2dvar=',use_similarity_2dvar if (hwrfblend) then inquire(file='hwrfinfo_input',exist=fexist) if(fexist) then open (55,file='hwrfinfo_input',form='formatted') read (55,hwrfinfo) close(55) print*,'in create_firstguess: ihwrfyear=',ihwrfyear print*,'in create_firstguess: ihwrfmonth=',ihwrfmonth print*,'in create_firstguess: ihwrfday=',ihwrfday print*,'in create_firstguess: ihwrfhour=',ihwrfhour endif endif !==================================================================== !==================================================================== ! Determine total number of slabs nt=nt0 present_sfcrough_slab=.false. if(use_similarity_2dvar) then nt=nt0+12 present_sfcrough_slab=.true. else if (neutral_stability_windfact_2dvar .and. .not.use_sfcr0) then nt=nt0+1 present_sfcrough_slab=.true. endif print*,'in create_firstguess: present_sfcrough_slab=',present_sfcrough_slab print*,'in create_firstguess: nt=',nt !==================================================================== !==================================================================== !==> read in analysis time info from namelist !==================================================================== read(9,timeinfo) print*,'in create_firstguess: iyear,imonth,iday,ihour, & & iminute,isecond=',iyear,imonth,iday,ihour, & iminute,isecond print*,'in create_firstguess: recomputeq =',recomputeq call domain_dims(cgrid,nx,ny,ds) print*,'in create_firstguess: nx,ny,ds=',nx,ny,ds ! !==================================================================== !==> read in namelists for bitmap removal & guess blending !==================================================================== inquire(file='bitmapinfo_input',exist=fexist) if(fexist) then open (55,file='bitmapinfo_input',form='formatted') read(55,bitmapinfo) close(55) endif print*,'in create_firstguess: lbitmap=',lbitmap ; print* !------------------------------------------------------ inquire (file='gesblendinfo_input',exist=fexist) if (fexist) then open (55,file='gesblendinfo_input',form='formatted') read(55,gesblendinfo) close(55) endif print*,'in create_firstguess: lgesblend,ngesblend=',lgesblend,ngesblend !------------------------------------------------------ if (lbitmap .or. lgesblend) then inquire(file='width_of_blendzone_input',exist=fexist) if(fexist) then open (55,file='width_of_blendzone_input',form='formatted') read(55,width_of_blendzone) close(55) endif do n=1,nt print*,'in create_firstguess: n, dmax(n) =', n, dmax(n) enddo print* endif !==================================================================== !==> remove bitmap if necessary !==================================================================== if (lbitmap) then print*,'in create_firstguess: lbitmap,cbitfile,cnobitfile,cmodelwithbitmap=', & lbitmap,trim(cbitfile),blank1,trim(cnobitfile),blank1,trim(cmodelwithbitmap) print*,'in create_firstguess: cbitfile_cldch_reflev,cnobitfile_cldch_reflev=', & cbitfile_cldch_reflev,blank1,cnobitfile_cldch_reflev call remove_bitmap(cgrid,nx,ny,nt,cbitfile,cnobitfile,cmodelwithbitmap, & cbitfile_cldch_reflev,cnobitfile_cldch_reflev,dmax(1:nt)) endif !==================================================================== !==> use blended first guess if requested !==================================================================== gesblend_ready=.false. if (lgesblend .and. ngesblend > 1) then print*,'in create_firstguess: using multi-fcst blended first guess' do n=1,ngesblend print*,'n,gesfcsthour(n),gesbinaryname(n),gesmodelname(n)=',& n,gesfcsthour(n),trim(gesbinaryname(n)),blank1,trim(gesmodelname(n)) print*,'n,gesmodel_cldch_reflev(n)=',n,gesmodel_cldch_reflev(n) enddo call blend_modelfcsts(cgrid,ngesblend,gesfcsthour(1:ngesblend), & gesbinaryname(1:ngesblend),gesmodelname(1:ngesblend), & gesmodel_cldch_reflev(1:ngesblend), & nx,ny,nt,dmax(1:nt),gesblend_ready) else print*,'in create_firstguess: using single fcst first guess' endif print*,'in create_firstguess: gesblend_ready=',gesblend_ready !==================================================================== !==> allocate fields !==================================================================== allocate(ifield(nx,ny)) allocate(ivgtyp(nx,ny)) allocate(isltyp(nx,ny)) allocate(glat(nx,ny)) allocate(glon(nx,ny)) allocate(field(nx,ny)) allocate(sst(nx,ny)) allocate(tslb(nx,ny)) allocate(tsk(nx,ny)) allocate(vegfra(nx,ny)) allocate(smois(nx,ny)) allocate(sfcrough(nx,ny)) allocate(qgrid(nx,ny)) allocate(tdgrid(nx,ny)) allocate(psfcgrid(nx,ny)) allocate(ugrid(nx,ny)) allocate(vgrid(nx,ny)) allocate(tgrid(nx,ny)) allocate(landmask(nx,ny)) allocate(dx(nx,ny)) allocate(dy(nx,ny)) allocate(mapfac(nx,ny)) allocate(phbgrid(nx,ny)) allocate(gust(nx,ny)) allocate(cldch(nx,ny)) allocate(vis(nx,ny)) allocate(pblh(nx,ny)) allocate(tcamt(nx,ny)) allocate(lcbas(nx,ny)) allocate(hgtgrid(nx,ny)) allocate(wspd(nx,ny)) allocate(mxtm(nx,ny)) allocate(mitm(nx,ny)) allocate(pmsl(nx,ny)) allocate(howv(nx,ny)) allocate(uwnd(nx,ny)) allocate(vwnd(nx,ny)) allocate(undef_cldch(nx,ny)) allocate(terrain(nx,ny)) if (use_similarity_2dvar) then allocate(presgrid1(nx,ny)) allocate(presgrid2(nx,ny)) allocate(tmpgrid1(nx,ny)) allocate(tmpgrid2(nx,ny)) allocate(qgrid1(nx,ny)) allocate(qgrid2(nx,ny)) allocate(ugrid1(nx,ny)) allocate(vgrid1(nx,ny)) allocate(hgtgrid1(nx,ny)) allocate(tggrid(nx,ny)) allocate(presgsfc(nx,ny)) endif ! !==================================================================== !=> generate file for surface analysis in " wrf first-guess format" !==================================================================== open (48,file='rtma_terrain.dat',form='unformatted') read(48) terrain close(48) rewind(30) read(30) glat read(30) glon read(30) mapfac read(40) landmask dx=ds*mapfac dy=ds*mapfac ihdrbuf=0 write(88) ihdrbuf write(88) iyear,imonth,iday,ihour,iminute,isecond,nx,ny,nsig write(88) dx,dy write(88) glat write(88) glon if (gesblend_ready) then lun=57 open(lun,file='slabs.dat_blended',form='unformatted') print*,'n create_firstguess: using slabs.dat_blended' cldch_reflev='ASL' !consistent with output from blending code else if (lbitmap) then lun=57 open(lun,file=trim(cbitfile)//'_patched',form='unformatted') print*,'n create_firstguess: using ',trim(cbitfile)//'_patched' cldch_reflev='ASL' !consistent with output from bitmap removing code else lun=20 print*,'n create_firstguess: using slabs.dat' !cldch_reflev value is from namelist endif endif rewind(lun) do 300 n=1,nt print*,'in create_firstguess: read slab# ',n read (lun) field print*,'slab,rmin,rmax=',trim(slabname(n)),minval(field),maxval(field) if (n.eq.1) psfcgrid=field if (n.eq.2) then hgtgrid=terrain !field phbgrid=hgtgrid*gravity !field*gravity ! PHB (zsfc*g) endif if (n.eq.3) then tgrid=field ! T(k) ! TEMP (sensible) sst=field tslb=field tsk=field do j=1,ny do i=1,nx if (landmask(i,j).gt.0.5) sst(i,j)=0. enddo enddo endif if (n.eq.4) tdgrid=field if (n.eq.5) ugrid=field if (n.eq.6) vgrid=field if (n.eq.7) then qgrid=field if (recomputeq) call td_to_q(tdgrid,psfcgrid,qgrid,nx,ny) endif if (n.eq.8) gust=field if (n.eq.9) then if (cldch_reflev .ne. 'ASL') & call clch_reflev_conversion(field,'fmodel_terrain.dat_on_'//trim(cgrid),cldch_reflev,'ASL',nx,ny) print*,'cldch,min,max,Above Sea Level: ', minval(field),maxval(field) !================================ call clch_reflev_conversion(field,'rtma_terrain.dat','ASL','AGL',nx,ny) print*,'cldch,min,max,Above Ground Level: ', minval(field),maxval(field) !================================ cldch=field undef_cldch(:,:)=1. do j=1,ny do i=1,nx if (abs(field(i,j)) > 0.1*cldch_bitmapval) then !use "abs" because actual bitmap value could be negative undef_cldch(i,j)=-1. endif enddo enddo open (55,file='cldch_map_input.dat',form='unformatted') write(55) undef_cldch close(55) print*,'undef_cldch,min,max=',minval(undef_cldch),maxval(undef_cldch) endif if (n.eq.10) then print*,'vis,min,max,before capping: ', minval(field),maxval(field) do j=1,ny do i=1,nx vis(i,j)=min(20000.,field(i,j)) enddo enddo print*,'vis,min,max,after capping: ', minval(vis),maxval(vis) endif if (n.eq.11) then pblh=field call pblh_80km(pblh,nx,ny) end if if (n.eq.12) mxtm=field if (n.eq.13) mitm=field if (n.eq.14) pmsl=field if (n.eq.15) then print*,'howv,min,max,before inpaint=',minval(field),maxval(field) call inpaint(size (field,1), size (field,2), field, landmask) howv=field print*,'howv,min,max,after inpaint=',minval(field),maxval(field) end if if (n.eq.16) tcamt=field if (n.eq.17) then print*,'lcbas,min,max,before capping: ', minval(field),maxval(field) do j=1,ny do i=1,nx lcbas(i,j)=min(20000.,field(i,j)) enddo enddo print*,'lcbas,min,max,after capping: ', minval(lcbas),maxval(lcbas) end if if (n.eq.(nt0+1)) sfcrough=field if (n.eq.(nt0+2)) presgrid1=field if (n.eq.(nt0+3)) presgrid2=field if (n.eq.(nt0+4)) tmpgrid1=field if (n.eq.(nt0+5)) tmpgrid2=field if (n.eq.(nt0+6)) qgrid1=field if (n.eq.(nt0+7)) qgrid2=field if (n.eq.(nt0+8)) ugrid1=field if (n.eq.(nt0+9)) vgrid1=field if (n.eq.(nt0+10)) then hgtgrid1=max(0.,field) print*,'hgtgrid1,min,max,after setting min value to zero: ', minval(hgtgrid1),maxval(hgtgrid1) endif if (n.eq.(nt0+11)) tggrid=field if (n.eq.(nt0+12)) presgsfc=field if (n.eq.(nt0+11)) write(6,*)'tggrid sample =', tggrid(10,10) print*,'===================================================================' print*,'===================================================================' print* 300 continue close(lun) wspd=sqrt(ugrid*ugrid+vgrid*vgrid) uwnd=ugrid vwnd=vgrid if (use_sfcr0) then sfcrough(:,:)=sfcr0 print*,'sfcrough,min,max,adjusted: ', minval(sfcrough),maxval(sfcrough) endif if (present_sfcrough_slab) then if (trim(cgrid) == "conus" .or. trim(cgrid) == "cohres" .or. & trim(cgrid) == "conusext" .or. & trim(cgrid) == "cohresext" .or. trim(cgrid) == "cohreswexp" .or. & trim(cgrid) == "alaska" .or. trim(cgrid) == "akhres") then sfcrough(:,:)=max(0.0001,sfcrough(:,:)) endif endif if (use_similarity_2dvar) then sst(:,:)=tggrid(:,:) print*,'sst,min,max,adjusted: ', minval(sst),maxval(sst) tslb(:,:)=tggrid(:,:) print*,'tslb,min,max,adjusted: ', minval(tslb),maxval(tslb) tsk(:,:)=tggrid(:,:) print*,'tsk,min,max,adjusted: ', minval(tsk),maxval(tsk) endif open (70,file='slabs2_nobiasc.dat',form='unformatted') write(70) psfcgrid write(70) hgtgrid write(70) tgrid write(70) tdgrid write(70) ugrid write(70) vgrid write(70) qgrid write(70) gust write(70) vis write(70) pblh write(70) cldch write(70) wspd write(70) tdgrid !same as ctrtdgrid write(70) mxtm write(70) mitm write(70) pmsl write(70) howv write(70) tcamt write(70) lcbas write(70) uwnd write(70) vwnd if (present_sfcrough_slab) then write(70) sfcrough ! Surface Roughtness endif if (use_similarity_2dvar) then write(70) presgrid1 ! PRESSURE at lowest half sigma level write(70) presgrid2 ! PRESSURE at 2nd lowest half sigma level write(70) tmpgrid1 ! Temperature at lowest half sigma level write(70) tmpgrid2 ! Temperature at 2nd lowest half sigma level write(70) qgrid1 ! Mixing Ratio at lowest half sigma level write(70) qgrid2 ! Mixing Ratio at 2nd lowest half sigma level write(70) ugrid1 ! U at lowest half sigma level write(70) vgrid1 ! V at lowest half sigma level write(70) hgtgrid1 ! HGT at owest half sigma level write(70) tggrid ! Surface Temperature write(70) presgsfc ! model surface pressure endif close(70) if (lbiascor) then inquire(file='rtma_biascor_in.dat',exist=fexist) if(fexist) then allocate(bias_psfcgrid(nx,ny)) allocate(bias_tgrid(nx,ny)) allocate(bias_qgrid(nx,ny)) allocate(bias_ugrid(nx,ny)) allocate(bias_vgrid(nx,ny)) allocate(bias_gust(nx,ny)) allocate(bias_vis(nx,ny)) allocate(bias_pblh(nx,ny)) allocate(bias_cldch(nx,ny)) allocate(bias_wspd(nx,ny)) allocate(bias_tdgrid(nx,ny)) allocate(bias_mxtm(nx,ny)) allocate(bias_mitm(nx,ny)) allocate(bias_pmsl(nx,ny)) allocate(bias_howv(nx,ny)) allocate(bias_tcamt(nx,ny)) allocate(bias_lcbas(nx,ny)) allocate(bias_uwnd(nx,ny)) allocate(bias_vwnd(nx,ny)) bias_psfcgrid=0. bias_tgrid=0. bias_qgrid=0. bias_ugrid=0. bias_vgrid=0. bias_tdgrid=0. bias_gust=0. bias_vis=0. bias_pblh=0. bias_cldch=0. bias_wspd=0. bias_mxtm=0. bias_mitm=0. bias_pmsl=0. bias_howv=0. bias_tcamt=0. bias_lcbas=0. bias_uwnd=0. bias_vwnd=0. open (60,file='rtma_biascor_in.dat',form='unformatted') read(60) bias_psfcgrid read(60) bias_tgrid read(60) bias_qgrid read(60) bias_ugrid read(60) bias_vgrid read(60) bias_gust read(60) bias_vis read(60) bias_pblh read(60) bias_cldch read(60) bias_wspd read(60) bias_tdgrid read(60) bias_mxtm read(60) bias_mitm read(60) bias_pmsl read(60) bias_howv read(60) bias_tcamt read(60) bias_lcbas read(60) bias_uwnd read(60) bias_vwnd close(60) if (pbiascor > epsilon) psfcgrid = psfcgrid + bias_psfcgrid if (tbiascor > epsilon) tgrid = tgrid + bias_tgrid if (qbiascor > epsilon) qgrid = max(qmin, (qgrid + bias_qgrid)) if (ubiascor > epsilon) ugrid = ugrid + bias_ugrid if (vbiascor > epsilon) vgrid = vgrid + bias_vgrid if (gustbiascor > epsilon) gust = gust + bias_gust if (visbiascor > epsilon) vis = vis + bias_vis if (pblhbiascor > epsilon) pblh = pblh + bias_pblh if (cldchbiascor> epsilon) cldch = cldch + bias_cldch if (wspdbiascor > epsilon) wspd = wspd + bias_wspd if (tdbiascor > epsilon) tdgrid = tdgrid + bias_tdgrid if (mxtmbiascor > epsilon) mxtm = mxtm + bias_mxtm if (mitmbiascor > epsilon) mitm = mitm + bias_mitm if (pmslbiascor > epsilon) pmsl = pmsl + bias_pmsl if (howvbiascor > epsilon) howv = howv + bias_howv if (tcamtbiascor> epsilon) tcamt = tcamt + bias_tcamt if (lcbasbiascor> epsilon) lcbas = lcbas + bias_lcbas if (uwndbiascor> epsilon) uwnd = uwnd + bias_uwnd if (vwndbiascor> epsilon) vwnd = vwnd + bias_vwnd tdgrid = min(tgrid,tdgrid) ! if (recomputeq) call td_to_q(tdgrid,psfcgrid,qgrid,nx,ny) deallocate(bias_psfcgrid) deallocate(bias_tgrid) deallocate(bias_qgrid) deallocate(bias_ugrid) deallocate(bias_vgrid) deallocate(bias_gust) deallocate(bias_vis) deallocate(bias_pblh) deallocate(bias_cldch) deallocate(bias_wspd) deallocate(bias_tdgrid) deallocate(bias_mxtm) deallocate(bias_mitm) deallocate(bias_pmsl) deallocate(bias_howv) deallocate(bias_tcamt) deallocate(bias_lcbas) deallocate(bias_uwnd) deallocate(bias_vwnd) else print*,'in create_firstguess: lbias=.true., but found no input bias file.' print*,'assume zero input bias' endif endif if (hwrfblend) then open (70,file='slabs2_pre_hwrf_blending.dat',form='unformatted') write(70) psfcgrid write(70) hgtgrid write(70) tgrid write(70) tdgrid write(70) ugrid write(70) vgrid write(70) qgrid write(70) gust write(70) vis write(70) pblh write(70) cldch write(70) wspd write(70) tdgrid !same as ctrtdgrid write(70) mxtm write(70) mitm write(70) pmsl write(70) howv write(70) tcamt write(70) lcbas write(70) uwnd write(70) vwnd if (present_sfcrough_slab) then write(70) sfcrough endif if (use_similarity_2dvar) then write(70) presgrid1 ! PRESSURE at lowest half sigma level write(70) presgrid2 ! PRESSURE at 2nd lowest half sigma level write(70) tmpgrid1 ! Temperature at lowest half sigma level write(70) tmpgrid2 ! Temperature at 2nd lowest half sigma level write(70) qgrid1 ! Mixing Ratio at lowest half sigma level write(70) qgrid2 ! Mixing Ratio at 2nd lowest half sigma level write(70) ugrid1 ! U at lowest half sigma level write(70) vgrid1 ! V at lowest half sigma level write(70) hgtgrid1 ! HGT at owest half sigma level write(70) tggrid ! Surface Temperature write(70) presgsfc ! model surface temperature endif close(70) iadate(1)=iyear ; iadate2(1)=ihwrfyear iadate(2)=imonth ; iadate2(2)=ihwrfmonth iadate(3)=iday ; iadate2(3)=ihwrfday iadate(4)=ihour ; iadate2(4)=ihwrfhour iadate(5)=iminute ; iadate2(5)=ihwrfminute do n=1,nhwrfmax call blendwind(cgrid,glat,glon,iadate,iadate2,ugrid,vgrid,gust,nx,ny,n) ! call blendwind(cgrid,glat,glon,iadate,iadate2,ugrid,vgrid,wspd,gust,nx,ny,n) !work on including wspd enddo endif !gust cannot be less than wspd do j=1,ny do i=1,nx if (gust(i,j) < wspd(i,j)) gust(i,j)=wspd(i,j) enddo enddo print*,'gust,min,max,adjusted: ', minval(gust),maxval(gust) field=0. ifield=0 write(88) psfcgrid ! psfc0 write(88) phbgrid ! PHB (zsfc*g) write(88) tgrid ! T(k) ! TEMP (sensible) write(88) qgrid ! Q(k) write(88) ugrid ! U(K) write(88) vgrid ! V(K) write(88) landmask ! LANDMASK (0=water and >0.5 for land) write(88) field ! XICE write(88) sst ! SST write(88) ifield ! IVGTYP write(88) ifield ! ISLTYP write(88) field ! VEGFRA write(88) field ! SNOW write(88) field ! SMOIS write(88) tslb ! TSLB write(88) tsk ! TSK write(88) sfcrough ! SFCROUGH write(88) gust ! GUST write(88) vis ! VIS write(88) pblh ! PBLH write(88) cldch ! CLDCH write(88) wspd ! WSPD write(88) tdgrid ! TDGRID write(88) mxtm ! MXTM write(88) mitm ! MITM write(88) pmsl ! PMSL write(88) howv ! HOWV write(88) tcamt ! TCAMT write(88) lcbas ! LCBAS write(88) uwnd ! U(K) !for direct u-analysis write(88) vwnd ! V(K) !for direct v-analysis if (use_similarity_2dvar) then write(88) presgrid1 ! PRESSURE at lowest half sigma level write(88) presgrid2 ! PRESSURE at 2nd lowest half sigma level write(88) tmpgrid1 ! Temperature at lowest half sigma level write(88) tmpgrid2 ! Temperature at 2nd lowest half sigma level write(88) qgrid1 ! Mixing Ratio at lowest half sigma level write(88) qgrid2 ! Mixing Ratio at 2nd lowest half sigma level write(88) ugrid1 ! U at lowest half sigma level write(88) vgrid1 ! V at lowest half sigma level write(88) hgtgrid1 ! HGT at owest half sigma level write(88) presgsfc ! model surface pressure write(6,*)'xxxx firstguess presgsfc =',presgsfc(10,10) endif open (70,file='slabs2.dat',form='unformatted') write(70) psfcgrid write(70) hgtgrid write(70) tgrid write(70) tdgrid write(70) ugrid write(70) vgrid write(70) qgrid write(70) gust write(70) vis write(70) pblh write(70) cldch write(70) wspd write(70) tdgrid !change to ctrtdgrid later as it could be different from tdgrid if bias correction is on write(70) mxtm write(70) mitm write(70) pmsl write(70) howv write(70) tcamt write(70) lcbas write(70) uwnd write(70) vwnd if (present_sfcrough_slab) then write(70) sfcrough ; open (71,file=' sfcrough.dat',form='unformatted') ; write(71) sfcrough ; close(71) endif if (use_similarity_2dvar) then write(70) presgrid1 ; open (71,file='presgrid1.dat',form='unformatted') ; write(71) presgrid1 ; close(71) write(70) presgrid2 ; open (71,file='presgrid2.dat',form='unformatted') ; write(71) presgrid2 ; close(71) write(70) tmpgrid1 ; open (71,file='tmpgrid1.dat',form='unformatted') ; write(71) tmpgrid1 ; close(71) write(70) tmpgrid2 ; open (71,file='tmpgrid2.dat',form='unformatted') ; write(71) tmpgrid2 ; close(71) write(70) qgrid1 ; open (71,file='qgrid1.dat',form='unformatted') ; write(71) qgrid1 ; close(71) write(70) qgrid2 ; open (71,file='qgrid2.dat',form='unformatted') ; write(71) qgrid2 ; close(71) write(70) ugrid1 ; open (71,file='ugrid1.dat',form='unformatted') ; write(71) ugrid1 ; close(71) write(70) vgrid1 ; open (71,file='vgrid1.dat',form='unformatted') ; write(71) vgrid1 ; close(71) write(70) hgtgrid1 ; open (71,file='hgtgrid1.dat',form='unformatted') ; write(71) hgtgrid1 ; close(71) write(70) tggrid ; open (71,file='tggrid.dat',form='unformatted') ; write(71) tggrid ; close(71) write(70) presgsfc ; open (71,file='presgsfc.dat',form='unformatted') ; write(71) presgsfc ; close(71) endif close(70) deallocate(ifield) deallocate(ivgtyp) deallocate(isltyp) deallocate(glat) deallocate(glon) deallocate(field) deallocate(sst) deallocate(tslb) deallocate(tsk) deallocate(vegfra) deallocate(smois) deallocate(qgrid) deallocate(tdgrid) deallocate(psfcgrid) deallocate(ugrid) deallocate(vgrid) deallocate(landmask) deallocate(dx) deallocate(dy) deallocate(mapfac) deallocate(phbgrid) deallocate(tgrid) deallocate(gust) deallocate(cldch) deallocate(vis) deallocate(pblh) deallocate(hgtgrid) deallocate(wspd) deallocate(mxtm) deallocate(mitm) deallocate(pmsl) deallocate(howv) deallocate(tcamt) deallocate(lcbas) deallocate(uwnd) deallocate(vwnd) deallocate(sfcrough) deallocate(undef_cldch) deallocate(terrain) if (use_similarity_2dvar) then deallocate(presgrid1) deallocate(presgrid2) deallocate(tmpgrid1) deallocate(tmpgrid2) deallocate(qgrid1) deallocate(qgrid2) deallocate(ugrid1) deallocate(vgrid1) deallocate(hgtgrid1) deallocate(tggrid) deallocate(presgsfc) endif if (mkrjlists) then fexist=.false. inquire(file='bigrjlist.txt',exist=fexist) if (fexist) call separate_rjlists if (.not.fexist) call join_rjlists var='t' ; call process_allrjlists(var) var='q' ; call process_allrjlists(var) var='p' ; call process_allrjlists(var) var='w' ; call process_allrjlists(var) var='mass' ; call process_allrjlists(var) var='wspd10m'; call process_allrjlists(var) var='td2m' ; call process_allrjlists(var) var='mitm' ; call process_allrjlists(var) var='mxtm' ; call process_allrjlists(var) var='pmsl' ; call process_allrjlists(var) var='howv' ; call process_allrjlists(var) var='tcamt'; call process_allrjlists(var) var='lcbas'; call process_allrjlists(var) var='cldch'; call process_allrjlists(var) ! var='uwnd10m'; call process_allrjlists(var) ! var='vwnd10m'; call process_allrjlists(var) endif end !======================================================================= !======================================================================= !*********************************************************************** subroutine td_to_q(td,p,q,nx,ny) ! ! prgmmr: pondeca org: np20 date: 2004-10-13 ! ! abstract: given dewpt in K and pressure in Pa, compute ! specific humidity in kg/kg. this is done by inverting: ! ! qv=q(i,j)/(1.-q(i,j)) ! e=p(i,j)/100.*qv/(eps+qv) ! eln=alog(e) ! td(i,j) = (243.5*eln-440.8)/(19.48-eln)+273.15 ! ! program history log: ! 2004-10-13 pondeca !*********************************************************************** implicit none real(4),parameter::alpha=243.5 real(4),parameter::beta=440.8 real(4),parameter::gamma=19.48 real(4),parameter::eta=273.15 real(4),parameter::eps=0.62197 !=Rd/Rv integer(4) nx,ny,i,j real(4),dimension(nx,ny),intent(in)::td,p real(4),dimension(nx,ny),intent(inout)::q real(4) loge,e,qv do j=1,ny do i=1,nx loge=(gamma*(td(i,j)-eta)+beta)/(alpha+td(i,j)-eta) e=exp(loge) qv=100.*e*eps/(p(i,j)-100.*e) q(i,j)=qv/(1.+qv) enddo enddo return end !======================================================================= !======================================================================= subroutine pblh_80km(grd,nx,ny) implicit none integer(4) nx,ny,i,j,n,ii,jj real(4),dimension(nx,ny),intent(inout):: grd real(4),dimension(0:nx+1,0:ny+1):: grd2 real(4) corn,cent,side,temp,c2,c3,c4,rterm ! Set weights for nine point smoother corn=0.3 cent=1.0 side=0.5 c4=4.0 c3=3.0 c2=2.0 rterm = 1.0/(cent + c4*side + c4*corn) ! 30 passes of 9pt smoother do n=1,30 do i=0,nx+1 do j=0,ny+1 ii=min(max(i,1),nx) ; jj=min(max(j,1),ny) grd2(i,j)=grd(ii,jj) end do end do do i=1,nx do j=1,ny temp = cent*grd2(i,j) + side*(grd2(i+1,j) + & grd2(i-1,j) + grd2(i,j+1) + grd2(i,j-1)) + & corn*(grd2(i+1,j+1) + grd2(i+1,j-1) + grd2(i-1,j-1) + & grd2(i-1,j+1)) grd(i,j) = temp*rterm end do end do end do return end subroutine pblh_80km !=========================================================================================== !=========================================================================================== subroutine remove_bitmap(cgrid,nx,ny,nt,cbitfile,cnobitfile,cmodelwithbitmap, & cbitfile_cldch_reflev,cnobitfile_cldch_reflev,dmax) implicit none integer(4),intent(in)::nx,ny,nt character(60),intent(in)::cgrid,cmodelwithbitmap character(60),intent(in)::cbitfile,cnobitfile character(3),intent(in)::cbitfile_cldch_reflev,cnobitfile_cldch_reflev real(4),intent(in)::dmax(nt) real(4),allocatable,dimension(:,:)::gfield1,gfield2,rd2 real(4),allocatable,dimension(:,:)::terrain integer(4) n print*,'in remove_bitmap' print*,'cgrid=',trim(cgrid) print*,'cmodelwithbitmap=',trim(cmodelwithbitmap) print*,'cbitfile=',trim(cbitfile) print*,'cnobitfile=',trim(cnobitfile) print*,'cbitfile_cldch_reflev=',cbitfile_cldch_reflev print*,'cnobitfile_cldch_reflev=',cnobitfile_cldch_reflev if (trim(cgrid).ne.'akhres' .and. trim(cgrid).ne.'cohresext' .and. trim(cgrid).ne.'cohreswexp') then print*,'error!!!! ... grid ',trim(cgrid),' is not supported by remove_bitmap' call abort stop endif allocate(terrain(nx,ny)) open (48,file='rtma_terrain.dat',form='unformatted') read(48) terrain close(48) allocate(rd2(nx,ny)) if (trim(cgrid)=='akhres') then if (trim(cmodelwithbitmap)=='RAP_ON_AK3P0' .or. & trim(cmodelwithbitmap)=='rap_on_ak3p0') then open (48,file='rd2.out_rap_on_ak3p0',form='unformatted') elseif (trim(cmodelwithbitmap)=='NAMNEST_ON_AK3P0' .or. & trim(cmodelwithbitmap)=='namnest_on_ak3p0') then open (48,file='rd2.out_namnest_on_ak3p0',form='unformatted') elseif (trim(cmodelwithbitmap)=='HRRR_ON_AK3P0' .or. & trim(cmodelwithbitmap)=='hrrr_on_ak3p0') then open (48,file='rd2.out_hrrr_on_ak3p0',form='unformatted') else print*,'error!!!! ... cmodelwithbitmap,',& trim(cmodelwithbitmap),' is not supported by remove_bitmap' endif endif if (trim(cgrid)=='cohresext') then if (trim(cmodelwithbitmap)=='NAMNEST_ON_COHRESEXT' .or. & trim(cmodelwithbitmap)=='namnest_on_cohresext') then open (48,file='rd2.out_namnest_on_cohresext',form='unformatted') elseif (trim(cmodelwithbitmap)=='HRRRNATIVE_ON_COHRESEXT' .or. & !this is for hrrr interpolated to grid 187. No downscaling! trim(cmodelwithbitmap)=='hrrrnative_on_cohresext') then open (48,file='rd2.out_hrrrnative_on_cohresext',form='unformatted') elseif (trim(cmodelwithbitmap)=='HRRR_ON_COHRESEXT' .or. & !hrrr downscaled to grid 187 trim(cmodelwithbitmap)=='hrrr_on_cohresext') then open (48,file='rd2.out_hrrr_on_cohresext',form='unformatted') else print*,'error!!!! ... cmodelwithbitmap,',& trim(cmodelwithbitmap),' is not supported by remove_bitmap' endif endif if (trim(cgrid)=='cohreswexp') then if (trim(cmodelwithbitmap)=='NAMNEST_ON_COHRESWEXP' .or. & trim(cmodelwithbitmap)=='namnest_on_cohreswexp') then open (48,file='rd2.out_namnest_on_cohreswexp',form='unformatted') elseif (trim(cmodelwithbitmap)=='HRRRNATIVE_ON_COHRESWEXP' .or. & trim(cmodelwithbitmap)=='hrrrnative_on_cohreswexp') then open (48,file='rd2.out_hrrrnative_on_cohreswexp',form='unformatted') elseif (trim(cmodelwithbitmap)=='HRRR_ON_COHRESWEXP' .or. & trim(cmodelwithbitmap)=='hrrr_on_cohreswexp') then open (48,file='rd2.out_hrrr_on_cohreswexp',form='unformatted') else print*,'error!!!! ... cmodelwithbitmap,',& trim(cmodelwithbitmap),' is not supported by remove_bitmap' endif endif read(48) rd2 close(48) open(51,file=trim(cnobitfile),form='unformatted') open(52,file=trim(cbitfile),form='unformatted') open(53,file=trim(cbitfile)//'_patched',form='unformatted') allocate(gfield1(nx,ny)) allocate(gfield2(nx,ny)) do n=1,nt read(51) gfield1 read(52) gfield2 if (trim(cgrid)=='akhres') then ! call blend02(cgrid,cmodelwithbitmap,gfield1,gfield2,rd2,nx,ny,n) call blend03(cgrid,cmodelwithbitmap,gfield1,gfield2,rd2,dmax(n),nx,ny,n) else call blend03(cgrid,cmodelwithbitmap,gfield1,gfield2,rd2,dmax(n),nx,ny,n) endif if (n.eq.2) write(53) terrain if (n.ne.2) write(53) gfield1 enddo close(51) close(52) close(53) deallocate(rd2) deallocate(terrain) deallocate(gfield1) deallocate(gfield2) print*,'exit remove_bitmap ... cleanly' end subroutine remove_bitmap !======================================================================= !======================================================================= subroutine blend02(cgridhost,cmodelwithbitmap,field1,field2,rd2,nx,ny,kin) implicit none !Declare passed variables character(60) ,intent(in) :: cgridhost character(60) ,intent(in) :: cmodelwithbitmap integer(4),intent(in) ::nx,ny,kin real(4),intent(inout)::field1(nx,ny) real(4),intent(in) ::field2(nx,ny) real(4),intent(in) ::rd2(nx,ny) !Declare local parameters real(4),parameter::spval=-9999. real(4),parameter::spval2=1.e10 real(4),parameter::epsilon=1.e-3 !Declare local variables real(4),allocatable,dimension(:,:):: fldaux real(4) alpha,diff,diff2 real(4) s1,rl,amp integer(4) i,j,n real(4) rmin,rmax data s1/8./ data rl/7./ data amp/1./ if (trim(cgridhost)=='cohresext') then if (trim(cmodelwithbitmap)=='NAMNEST_ON_COHRESEXT' .or. trim(cmodelwithbitmap)=='namnest_on_cohresext') rl=8. if (trim(cmodelwithbitmap)=='HRRR_ON_COHRESEXT' .or. trim(cmodelwithbitmap)=='hrrr_on_cohresext') rl=8. if (trim(cmodelwithbitmap)=='HRRRNATIVE_ON_COHRESEXT' .or. trim(cmodelwithbitmap)=='hrrrnative_on_cohresext') rl=8. endif if (trim(cgridhost)=='cohreswexp') then if (trim(cmodelwithbitmap)=='NAMNEST_ON_COHRESWEXP' .or. trim(cmodelwithbitmap)=='namnest_on_cohreswexp') rl=8. if (trim(cmodelwithbitmap)=='HRRR_ON_COHRESWEXP' .or. trim(cmodelwithbitmap)=='hrrr_on_cohreswexp') rl=8. if (trim(cmodelwithbitmap)=='HRRRNATIVE_ON_COHRESWEXP' .or. trim(cmodelwithbitmap)=='hrrrnative_on_cohreswexp') rl=8. endif if (trim(cgridhost)=='akhres') then if (trim(cmodelwithbitmap)=='RAP_ON_AK3P0' .or. trim(cmodelwithbitmap)=='rap_on_ak3p0') rl=8. if (trim(cmodelwithbitmap)=='NAMNEST_ON_AK3P0' .or. trim(cmodelwithbitmap)=='namnest_on_ak3p0') rl=8. if (trim(cmodelwithbitmap)=='HRRR_ON_AK3P0' .or. trim(cmodelwithbitmap)=='hrrr_on_ak3p0') rl=8. endif print*,'in blend02: cgridhost,cmodelwithbitmap=',trim(cgridhost),trim(cmodelwithbitmap) print*,'in blend02: spval,spval2=',spval,spval2 print*,'in blend02: s1,rl,amp=',s1,rl,amp ; print* call rminmax_excludespval(field2,nx,ny,1.e19,rmin,rmax) print*,'in blend02: k,field2,min,max=',kin,rmin,rmax call rminmax_excludespval(field1,nx,ny,1.e19,rmin,rmax) print*,'in blend02: before blending k,field1,min,max=',kin,rmin,rmax if (kin==2) goto 2222 !do not change the terrain if (kin==10) goto 2222 !do not change the vis field if (kin==12) goto 2222 !do not change mxtm if (kin==13) goto 2222 !do not change mitm if (kin==15) goto 2222 !do not change howv if (kin==17) goto 2222 !do not change lcbas. Must fix this. lcbas from rap contains undefined values all over allocate(fldaux(nx,ny)) do j=1,ny do i=1,nx diff=abs(rd2(i,j)-spval) diff2=abs(rd2(i,j)-spval2) if (diff < epsilon) fldaux(i,j)=field1(i,j) if (diff2 < epsilon) fldaux(i,j)=field2(i,j) if (diff > epsilon .and. diff2 > epsilon ) then alpha=amp*exp(-rd2(i,j)*rd2(i,j)/(s1*rl*rl)) fldaux(i,j)=alpha*field1(i,j)+(1.-alpha)*field2(i,j) endif enddo enddo do j=1,ny do i=1,nx field1(i,j)=fldaux(i,j) enddo enddo deallocate(fldaux) 2222 continue call rminmax_excludespval(field1,nx,ny,1.e19,rmin,rmax) print*,'in blend02: after blending k,field1,min,max=',kin,rmin,rmax ; print* end subroutine blend02 !============================================================== !======================================================================= !======================================================================= subroutine blend03(cgridhost,cmodelwithbitmap,field1,field2,rd2,dmax0,nx,ny,kin) use pbend implicit none !Declare passed variables character(60) ,intent(in) :: cgridhost character(60) ,intent(in) :: cmodelwithbitmap integer(4),intent(in) ::nx,ny,kin real(4),intent(inout)::field1(nx,ny) real(4),intent(in) ::field2(nx,ny) real(4),intent(in) ::rd2(nx,ny) real(4),intent(in) ::dmax0 !Declare local parameters real(4),parameter::spval=-9999. real(4),parameter::spval2=1.e10 real(4),parameter::epsilon=1.e-3 !Declare local variables real(4),allocatable,dimension(:,:):: fldaux real(4) diff,diff2 real(8) x,alpha !double precision integer(4) i,j,n real(4) rmin,rmax logical visblend logical windblend logical gustblend logical lcbasblend logical cldchblend logical lcase1 logical lcase2 logical lcase3 lcase1=trim(cgridhost)=='cohresext' .and. & (trim(cmodelwithbitmap)=='HRRR_ON_COHRESEXT'.or. & trim(cmodelwithbitmap)=='hrrr_on_cohresext'.or. & trim(cmodelwithbitmap)=='HRRRNATIVE_ON_COHRESEXT'.or. & trim(cmodelwithbitmap)=='hrrrnative_on_cohresext') lcase2=trim(cgridhost)=='cohreswexp' .and. & (trim(cmodelwithbitmap)=='HRRR_ON_COHRESWEXP'.or. & trim(cmodelwithbitmap)=='hrrr_on_cohreswexp'.or. & trim(cmodelwithbitmap)=='HRRRNATIVE_ON_COHRESWEXP'.or. & trim(cmodelwithbitmap)=='hrrrnative_on_cohreswexp') lcase3=trim(cgridhost)=='akhres' .and. & (trim(cmodelwithbitmap)=='HRRR_ON_AK3P0'.or. & trim(cmodelwithbitmap)=='hrrr_on_ak3p0') visblend=.false. windblend=.false. gustblend=.false. lcbasblend=.false. cldchblend=.false. if ( lcase1 .or. lcase2 .or. lcase3 ) then visblend=.true. windblend=.true. gustblend=.true. lcbasblend=.true. cldchblend=.true. endif ! if (trim(cgridhost)=='akhres') then ! windblend=.true. ! gustblend=.true. ! endif print* print*,'===================================================================' print*,'in blend03: cgridhost=',trim(cgridhost) print*,'in blend03: cmodelwithbitmap=',trim(cmodelwithbitmap) print*,'in blend03: spval,spval2=',spval,spval2 print*,'in blend03: dmax0=',dmax0 ; print* print*,'in blend03: visblend=',visblend print*,'in blend03: windblend=',windblend print*,'in blend03: gustblend=',gustblend print*,'in blend03: lcbasblend=',lcbasblend print*,'in blend03: cldchblend=',cldchblend ; print* call rminmax_excludespval(field2,nx,ny,1.e19,rmin,rmax) print*,'in blend03: k,field2,min,max=',kin,rmin,rmax call rminmax_excludespval(field1,nx,ny,1.e19,rmin,rmax) print*,'in blend03: before blending k,field1,min,max=',kin,rmin,rmax if (kin==2) goto 2222 !do not change the terrain if (kin==5 .and. .not.windblend) goto 2222 !do not change u-wind if (kin==6 .and. .not.windblend) goto 2222 !do not change v-wind if (kin==8 .and. .not.gustblend) goto 2222 !do not change gust if (kin==9 .and. .not.cldchblend) goto 2222 !do not change cldch if (kin==10 .and. .not.visblend) goto 2222 !do not change the vis field if (kin==12) goto 2222 !do not change mxtm if (kin==13) goto 2222 !do not change mitm if (kin==15) goto 2222 !do not change howv if (kin==17 .and. .not.lcbasblend)goto 2222 !do not change lcbas !lcbas from rap contains undefined values all over allocate(fldaux(nx,ny)) do j=1,ny do i=1,nx diff=abs(rd2(i,j)-spval) diff2=abs(rd2(i,j)-spval2) if (diff < epsilon) fldaux(i,j)=field1(i,j) if (diff2 < epsilon) fldaux(i,j)=field2(i,j) if (diff >= epsilon .and. diff2 >= epsilon ) then x=min(1.,rd2(i,j)/dmax0) call wbend(x,alpha) fldaux(i,j)=alpha*field2(i,j)+(1.-alpha)*field1(i,j) endif enddo enddo do j=1,ny do i=1,nx field1(i,j)=fldaux(i,j) enddo enddo deallocate(fldaux) 2222 continue call rminmax_excludespval(field1,nx,ny,1.e19,rmin,rmax) print*,'in blend03: after blending k,field1,min,max=',kin,rmin,rmax ; print* end subroutine blend03 !======================================================================= !======================================================================= subroutine blend_modelfcsts(cgrid,ngesblend,gesfcsthour, & gesbinaryname,gesmodelname,gesmodel_cldch_reflev, & nx,ny,nt,dmax,gesblend_ready) implicit none integer(4),intent(in)::nx,ny,nt,ngesblend real(4),intent(in)::dmax(nt) real(4),intent(in)::gesfcsthour(ngesblend) character(60),intent(in)::cgrid character(60),intent(in)::gesbinaryname(ngesblend) character(60),intent(in)::gesmodelname(ngesblend) character(3),intent(in)::gesmodel_cldch_reflev(ngesblend) logical,intent(out)::gesblend_ready integer(4),parameter::lun=57 real(4),allocatable,dimension(:,:,:):: baseges real(4),allocatable,dimension(:,:):: field,fieldaux real(4),allocatable,dimension(:,:,:):: fieldsnew logical,allocatable,dimension(:)::gexist logical,allocatable,dimension(:)::bitmapyes real(4),allocatable,dimension(:,:):: terrain real(4),allocatable,dimension(:,:):: rd2,rd2aux1,rd2aux2,rd2aux3 integer(4) i,j integer(4) k,n,kflds real(4) weight,wsum character(60) cmodel logical fexist logical gexist0 logical lvalid1,lvalid2,lvalid3, lvalid4 print*,'===== entering blend_modelfcsts =====' if (trim(cgrid).ne.'akhres' .and. trim(cgrid).ne.'cohresext' .and. & trim(cgrid).ne.'cohreswexp' .and. trim(cgrid).ne.'prico') then print*,'error!!!! ... grid ',trim(cgrid),' is not supported by blend_modelfcsts' call abort stop endif print*,'cgrid=',trim(cgrid) ; print* do n=1,ngesblend cmodel=gesmodelname(n) lvalid1=trim(cmodel)=='RAP_ON_AK3P0'.or.trim(cmodel)=='rap_on_ak3p0'.or. & trim(cmodel)=='NAMNEST_ON_AK3P0'.or.trim(cmodel)=='namnest_on_ak3p0' .or. & trim(cmodel)=='HRRR_ON_AK3P0'.or.trim(cmodel)=='hrrr_on_ak3p0' lvalid2=trim(cmodel)=='RAP_ON_COHRESEXT'.or.trim(cmodel)=='rap_on_cohresext'.or. & trim(cmodel)=='NAMNEST_ON_COHRESEXT'.or.trim(cmodel)=='namnest_on_cohresext' lvalid3=trim(cmodel)=='RAP_ON_PR2P5'.or.trim(cmodel)=='rap_on_pr2p5'.or. & trim(cmodel)=='NAMNEST_ON_PR2P5'.or.trim(cmodel)=='namnest_on_pr2p5' lvalid4=trim(cmodel)=='RAP_ON_COHRESWEXP'.or.trim(cmodel)=='rap_on_cohreswexp'.or. & trim(cmodel)=='NAMNEST_ON_COHRESWEXP'.or.trim(cmodel)=='namnest_on_cohreswexp' if ( (trim(cgrid)=='akhres' .and. .not.lvalid1) .or. & (trim(cgrid)=='cohresext' .and. .not.lvalid2) .or. & (trim(cgrid)=='prico' .and. .not.lvalid3) .or. & (trim(cgrid)=='cohreswexp'.and. .not.lvalid4) ) then print*,'error!!!! ... cmodel ',trim(cmodel),' is not supported by blend_modelfcsts' call abort stop endif enddo allocate(terrain(nx,ny)) open (lun,file='rtma_terrain.dat',form='unformatted') read(lun) terrain close(lun) allocate(rd2(nx,ny)) allocate(rd2aux1(nx,ny)) allocate(rd2aux2(nx,ny)) allocate(rd2aux3(nx,ny)) if (trim(cgrid)=='akhres') then inquire(file='rd2.out_rap_on_ak3p0',exist=fexist) if (fexist) then open (48,file='rd2.out_rap_on_ak3p0',form='unformatted') read(48) rd2aux1 close(48) else print*,'WARNING: FILE rd2.out_rap_on_ak3p0 NOT PROVIDED. & &WILL ASSUME THAT RAP FIRST GUESS HAS NO BITMAP' rd2aux1(:,:)=1.e+10 endif open (48,file='rd2.out_namnest_on_ak3p0',form='unformatted') read(48) rd2aux2 close(48) open (48,file='rd2.out_hrrr_on_ak3p0',form='unformatted') read(48) rd2aux3 close(48) endif if (trim(cgrid)=='cohresext') then open (48,file='rd2.out_namnest_on_cohresext',form='unformatted') read(48) rd2 close(48) endif if (trim(cgrid)=='cohreswexp') then open (48,file='rd2.out_namnest_on_cohreswexp',form='unformatted') read(48) rd2 close(48) endif allocate(bitmapyes(ngesblend)) bitmapyes(:)=.false. allocate(gexist(ngesblend)) gexist(:)=.false. gesblend_ready=.false. allocate(field(nx,ny)) do n=1,ngesblend inquire(file=trim(gesbinaryname(n)),exist=gexist0) print*,'n,gesbinaryname(n),gexist0=',n,trim(gesbinaryname(n)),gexist0 gexist(n)=gexist0 if (gexist0) then open(lun,file=trim(gesbinaryname(n)),form='unformatted') kflds=0 do k=1,nt read(lun,end=1111) field kflds=kflds+1 if (k==1 .and. maxval(field) > 1.e18) bitmapyes(n)=.true. enddo 1111 continue close(lun) if (kflds < nt) then print*,'binary field does not contain enough records: kflds,nt=',kflds,nt print*,'reset gexist(n) to .false.' gexist(n)=.false. endif endif if (gexist(n)) gesblend_ready=.true. enddo if (.not.gesblend_ready) then print*,'missing all partial guesses. no blending& &will be performed' return endif print* do n=1,ngesblend print*,'n,gesmodelname(n),gesfcsthour(n)=', & n,trim(gesmodelname(n)),gesfcsthour(n) print*,'n,gesbinaryname(n),gexist(n),bitmapyes(n)=',& n,trim(gesbinaryname(n)),gexist(n),bitmapyes(n) ; print* if (trim(cgrid)=='prico' .and. gexist(n) .and. bitmapyes(n)) then print*,'error in blend_modelfcsts. there should be no bitmap for prico' call abort stop endif enddo ! important! the first binary file must be a nobitmap file. ! load its content to baseges and use to remove the bitmap from the other ! ges files allocate(baseges(nx,ny,nt)) open(lun,file=trim(gesbinaryname(1)),form='unformatted') do k=1,nt read(lun) field baseges(:,:,k)=field(:,:) enddo close(lun) ! do the blending allocate(fieldaux(nx,ny)) allocate(fieldsnew(nx,ny,nt)) fieldsnew(:,:,:)=0. wsum=0. do n=1,ngesblend if (.not.gexist(n)) cycle open(lun,file=trim(gesbinaryname(n)),form='unformatted') cmodel=gesmodelname(n) if (trim(cgrid)=='akhres') then if (trim(cmodel)=='RAP_ON_AK3P0'.or. & trim(cmodel)=='rap_on_ak3p0') rd2=rd2aux1 if (trim(cmodel)=='NAMNEST_ON_AK3P0' .or. & trim(cmodel)=='namnest_on_ak3p0') rd2=rd2aux2 if (trim(cmodel)=='HRRR_ON_AK3P0'.or. & trim(cmodel)=='hrrr_on_ak3p0') rd2=rd2aux3 endif IF (BITMAPYES(N)) OPEN (31,FILE=trim(gesbinaryname(n))//'_NOBITMAP',FORM='UNFORMATTED') weight=1./sqrt(max(1.,gesfcsthour(n))) wsum=wsum+weight print*,'n,gesbinaryname(n),gesfcsthour(n),weight=',n,trim(gesbinaryname(n)),gesfcsthour(n),weight print*,'n,gesmodel_cldch_reflev(n)=',n,gesmodel_cldch_reflev(n) ; print* do k=1,nt read(lun) fieldaux if (k==9) call clch_reflev_conversion(fieldaux,'terrain.dat_'//trim(cmodel),gesmodel_cldch_reflev(n),'ASL',nx,ny) if (bitmapyes(n)) then field(:,:)=baseges(:,:,k) if (trim(cgrid)=='akhres') then ! call blend02(cgrid,cmodel,field,fieldaux,rd2,nx,ny,k) call blend03(cgrid,cmodel,field,fieldaux,rd2,dmax(k),nx,ny,k) else call blend03(cgrid,cmodel,field,fieldaux,rd2,dmax(k),nx,ny,k) endif else field(:,:)=fieldaux(:,:) endif fieldsnew(:,:,k)=fieldsnew(:,:,k)+weight*field(:,:) IF (BITMAPYES(N)) WRITE (31) FIELD enddo close(lun) IF (BITMAPYES(N)) CLOSE (31) enddo print*,'wsum=',wsum open (lun,file='slabs.dat_blended',form='unformatted') do k=1,nt if (k.eq.2) then field(:,:)=terrain(:,:) ! elseif (k.eq.9) then !no blending for cldch for now ! field(:,:)=baseges(:,:,k) elseif (k.eq.12 .or. k.eq.13) then !no blending for maxT,minT, field(:,:)=baseges(:,:,k) !because all guess fields use the same maxT,minT elseif (k.eq.15 .or. k.eq.16) then !no blending for howv or tcamt field(:,:)=baseges(:,:,k) !because no valid aux fields available elseif (k.eq.17) then !no blending for lcbas for now field(:,:)=baseges(:,:,k) !similarity elseif (k.ge.18 .and. k .le.nt) then field(:,:)=baseges(:,:,k) else field(:,:)=fieldsnew(:,:,k)/wsum endif write(lun) field enddo close(lun) deallocate(bitmapyes) deallocate(gexist) deallocate(terrain) deallocate(rd2) deallocate(rd2aux1) deallocate(rd2aux2) deallocate(rd2aux3) deallocate(field) deallocate(fieldaux) deallocate(fieldsnew) deallocate(baseges) print*,'===== exiting blend_modelfcsts =====' return end subroutine blend_modelfcsts !======================================================================= !======================================================================= !======================================================================= subroutine clch_reflev_conversion(field,toponame,input_reflev,output_reflev,nx,ny) implicit none integer(4),intent(in):: nx,ny character(*),intent(in):: input_reflev character(*),intent(in):: output_reflev character(*),intent(in):: toponame real(4),intent(inout):: field(nx,ny) integer(4),parameter:: lun=27 real(4),parameter:: cldch_bitmapval=1.e+21 real(4),parameter:: cldch_cap=20000. real(4),parameter::small_cldchdiff=0.15 integer(4) i,j logical fexist logical lcond1,lcond2,lcond3 real(4),allocatable,dimension(:,:):: modelterrain real(4) amax print* ; print*,'===== entering clch_reflev_conversion =====' print*,'in clch_reflev_conversion: input_reflev=',input_reflev print*,'in clch_reflev_conversion: output_reflev=',output_reflev print*,'in clch_reflev_conversion: toponame=',toponame if ( (input_reflev.ne.'AGL' .and. input_reflev.ne.'ASL') .or. & (output_reflev.ne.'AGL' .and. output_reflev.ne.'ASL') ) then print*,'in clch_reflev_conversion: incorrect specification & &for input_reflev or output_reflev. they should be & &either AGL or ASL. ... aborting ...' call abort stop endif if (input_reflev .ne. output_reflev) then inquire(file=trim(toponame),exist=fexist) if(fexist) then allocate(modelterrain(nx,ny)) open (lun,file=trim(toponame),form='unformatted') read (lun) modelterrain close(lun) print*,'in clch_reflev_conversion:input cldch,min,max: ', minval(field),maxval(field) amax=maxval(abs(field)) print*,'in clch_reflev_conversion: amax=',amax if (output_reflev.eq.'ASL') then do j=1,ny do i=1,nx lcond1=amax>0.1*cldch_bitmapval .and. abs(field(i,j))<0.1*cldch_bitmapval lcond2=abs(amax-cldch_cap)small_cldchdiff !case when max in smartinit is cldch_cap lcond3=amaxcldch_cap.and.amax<0.1*cldch_bitmapval) !case when smartinit is coded to use cldch_bitmapval for clear-sky !but there is no clear-sky in the domain if (lcond1.or.lcond2.or.lcond3) then field(i,j)=field(i,j)+max(modelterrain(i,j),0.0) endif field(i,j)=max(min(cldch_cap,field(i,j)),0.1) enddo enddo else if (output_reflev.eq.'AGL') then do j=1,ny do i=1,nx lcond1=amax>0.1*cldch_bitmapval .and. abs(field(i,j))<0.1*cldch_bitmapval lcond2=abs(amax-cldch_cap)small_cldchdiff !case when max in smartinit is cldch_cap lcond3=amaxcldch_cap.and.amax<0.1*cldch_bitmapval) !case when smartinit is coded to use cldch_bitmapval for clear-sky !but there is no clear-sky in the domain if (lcond1.or.lcond2.or.lcond3) then field(i,j)=field(i,j)-max(modelterrain(i,j),0.0) endif field(i,j)=max(min(cldch_cap,field(i,j)),0.1) enddo enddo endif print*,'in clch_reflev_conversion:output cldch,min,max: ', minval(field),maxval(field) deallocate(modelterrain) else print*,'in clch_reflev_conversion: missing terrain file: ',toponame, '...aborting' call abort stop endif else print*,'in clch_reflev_conversion: no conversion necessary. input_reflev and output_reflev are the same' endif print*,'===== exiting clch_reflev_conversion =====' ; print* return end !======================================================================= !======================================================================= !======================================================================= !======================================================================= subroutine smoother_desmoother(g1,is,ie,js,je,smcf1,smcf2,one2one,ns) !$$$ subprogram documentation block ! . . . . ! subprogram: smoother_desmoother ! prgmmr: pondeca org: np22 date: 2006-08-01 ! ! abstract: apply smoother-desmoother to 2-d data slab ! ! program history log: ! 2011-09-15 pondeca ! ! input argument list: ! g1 - 2d array of field to be smoothed ! is,ie - first and last i values of data slab ! js,je - first and last j values of data slab ! ns - number of passes ! ! output argument list: ! g1 - smoothed 2d field ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block implicit none integer, parameter :: r_single = selected_real_kind(6) ! single precision integer, parameter :: i_kind = selected_int_kind(8) ! long integer integer(i_kind) , intent(in ) :: is, ie, js, je integer(i_kind) , intent(in ) :: ns real(r_single), dimension(is:ie, js:je), intent(inout) :: g1 ! on input: original data slab ! on ouput: filtered data slab real(r_single), intent(in) :: smcf1 ! smoothing coefficient real(r_single), intent(in) :: smcf2 ! desmoothing coefficient logical,intent(in) :: one2one ! if .true., perform 1-2-1 smoothing instead integer(i_kind) i,j,l,ip,im,jp,jm,istep real(r_single), allocatable:: g2(:,:) real(r_single) smcf0 real(r_single),parameter::half=0.5_r_single real(r_single),parameter::quarter=0.25_r_single allocate(g2(is:ie,js:je)) do l=1,ns if (.not.one2one) then do istep=1,2 if (istep==1) smcf0=+abs(smcf1) if (istep==2) smcf0=-abs(smcf2) do j=js,je do i=is,ie ip=min(i+1,ie) ; im=max(is,i-1) g2(i,j)=g1(i,j)+smcf0*(0.5*(g1(ip,j)+g1(im,j))-g1(i,j)) end do end do do i=is,ie do j=js,je jp=min(j+1,je) ; jm=max(js,j-1) g1(i,j)=g2(i,j)+smcf0*(0.5*(g2(i,jp)+g2(i,jm))-g2(i,j)) end do end do end do elseif (one2one) then do j=js,je do i=is,ie ip=min(i+1,ie) ; im=max(is,i-1) g2(i,j)=quarter*(g1(ip,j)+g1(im,j))+half*g1(i,j) end do end do do i=is,ie do j=js,je jp=min(j+1,je) ; jm=max(js,j-1) g1(i,j)=quarter*(g2(i,jp)+g2(i,jm))+half*g2(i,j) end do end do endif end do deallocate(g2) return end subroutine smoother_desmoother !======================================================================= !======================================================================= subroutine rminmax_excludespval(field,nx,ny,spval,rmin,rmax) implicit none integer(4),intent(in)::nx,ny real(4),intent(in)::spval real(4),intent(in)::field(nx,ny) real(4),intent(out)::rmin,rmax integer(4) i,j rmin=+huge(rmin) rmax=-huge(rmax) do j=1,ny do i=1,nx if (field(i,j) < spval ) then if (field(i,j) < rmin) rmin=field(i,j) if (field(i,j) > rmax) rmax=field(i,j) endif enddo enddo end subroutine rminmax_excludespval !> @brief Writes a 2D array to text file, column by column !! @param [in] fileName path to the output file !! @param [in] rda_A 2D array to write !! @details !< subroutine writeMatrix(fileName, rda_A) real(4), dimension(:, :), intent(in) :: rda_A character(len=*) , intent(in) :: fileName integer ib_i, ib_j, il_ios integer, parameter :: ip_fid = 41 ! logical :: iAmMaster ! ! iAmMaster = parallel_amIMaster() ! ! if(iAmMaster)then open( unit = ip_fid, file = fileName, status = 'replace' & , form = 'formatted', iostat = il_ios) if (il_ios /= 0) print*,'In writeMatrix : Error creating file'//fileName do ib_j = 1, size(rda_A,2) do ib_i = 1, size(rda_A,1) write(unit=ip_fid, fmt='(E18.8, $)') rda_A(ib_i,ib_j)!$ to avoid end of line end do write(unit=ip_fid, fmt=*)''! just for adding end of line end do close(ip_fid ) ! end if end subroutine writeMatrix !======================================================================= !======================================================================= ! Subroutine: inpaint ! Calling Functions : meanpos ! Author : stelios flampouris ! Comments, suggestions are welcomed at ! stylianos.flampouris@noaa.gov (gmail.com) ! License : Source code written by stelios is distributed under ! GNU General Public License v3.0 ! ! v0.5.prttp 11.04.2015 ! v0.7. 12.10.2015 ! v1.0 03.29.2017 ! ! ! ! Variables ! xlen :: length of x-dim (e.g longtitude) ! ylen :: length of y-dim (e.g latitude) ! arr :: 2D-array(lenght(x),length(y)) with data which has the gaps. The gaps have NaN as values ! subroutine inpaint(xlen, ylen, arr, arrmask) implicit none ! integer :: xlen,ylen real(4), dimension(xlen,ylen), intent(inout) :: arr real(4), dimension(xlen,ylen), intent(in) :: arrmask integer, parameter :: it_max = 50 ! !local logical, dimension(xlen,ylen) :: indNaN integer :: ix, iy, indx_min, indx_max, indy_min, indy_max, ext_it integer :: strt_ix, end_ix, strt_iy, end_iy,stp_ix, stp_iy real(4) :: max_val, Avg ! ! 1.Basic Input Manipulation max_val=30. ! indNaN = (arr>max_val).and.(arrmask.eq.0.) ! do ix = 1,xlen,1 do iy = 1,ylen,1 if (arr(ix,iy) > max_val ) then arr (ix,iy) = 0. end if enddo enddo ! ! 2. The hoop do ext_it = 1,it_max,1 do ix = 1,xlen,1 do iy = 1,ylen,1 if (.not.indNaN(ix,iy)) cycle indx_min = max(1,ix-1) indy_min = max(1,iy-1) indx_max = min(xlen,ix+1) indy_max = min(ylen,iy+1) if ((arr(ix,iy).eq.0.).and.(count(indNaN(indx_min:indx_max,indy_min:indy_max))>0) ) then arr(ix,iy) = meanpos(arr(indx_min:indx_max,indy_min:indy_max)) end if end do end do indNaN = (arr.eq.0.).and.(arrmask.eq.0.) if (count(indNaN)==0) exit end do ! ! 3. The remaining points Avg = 0.01*floor(100*meanpos(arr)) do ix = 1,xlen,1 do iy = 1,ylen,1 if (arr(ix,iy).eq.0) then arr (ix,iy) = Avg endif enddo enddo ! contains function meanpos(array) real :: meanpos real, intent(in), dimension(:,:) :: array meanpos = sum (array, array > 0) / max(1,count(array > 0 )) end function meanpos ! function mean(array) real :: mean real, intent(in), dimension(:,:) :: array mean = sum (array) / max(1,size(array)) end function mean ! end subroutine inpaint