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.and.abs(cldch_cap-abs(field(i,j)))>small_cldchdiff !case when max in smartinit is cldch_cap
lcond3=amax<cldch_cap .or. & !case when the max assumed in smartinit does not matter
(amax>cldch_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.and.abs(cldch_cap-abs(field(i,j)))>small_cldchdiff !case when max in smartinit is cldch_cap
lcond3=amax<cldch_cap .or. & !case when the max assumed in smartinit does not matter
(amax>cldch_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