subroutine anbkerror(gradx,grady)
!$$$ subprogram documentation block
! . . . .
! subprogram: anbkerror apply anisotropic background error covariance
! prgmmr: parrish org: np22 date: 2005-02-03
!
! abstract: apply regional anisotropic background error.
!
! program history log:
! 2005-02-08 parrish
! 2005-04-29 parrish - replace coarse2fine with fgrid2agrid;
! remove ansmoothrf_reg_d
! 2006-11-30 todling - add fpsproj as arg to (t)balance routine(s)
! 2008-10-10 derber - add strong constraint to background error
! 2008-12-29 todling - update interface to strong_bk/bk_ad
! 2009-04-13 derber - move strong_bk into balance
! 2009-07-01 sato - update for global mode
! 2010-05-05 todling - update to use gsi_bundle
! 2010-06-22 todling - update to better handle bundle pointers
! 2010-06-29 lueken - replaced tv with t in call to gsi_bundlegetpointer
! 2010-08-19 lueken - add only to module use
! 2012-10-09 Gu - add fut2ps as arg to (t)balance routine(s)
! 2013-05-23 zhu - add ntclen for aircraft temperature bias correction
! 2014-02-07 pondeca - update to handle motley variables. rename p_st to p_sf
! 2014-02-14 pondeca - update to handle optional separation of sf and vp control variables
! into land-only and water-only parts
! 2015-05-02 parrish - add subroutine ansmoothrf_reg_sub2slab_option, and
! parameter rtma_bkerr_sub2slab.
! rtma_bkerr_sub2slab = F, then use ansmoothrf_reg_subdomain_option
! rtma_bkerr_sub2slab = T, then use ansmoothrf_reg_sub2slab_option
! This allows dual resolution for the anisotropic recursive filter, which
! currently can only be used with full horizontal domain (slab) storage.
! 2015-07-02 pondeca - update slab mode option to work with any number of control variables
!
! input argument list:
! gradx - input field
!
! output
! grady - background structure * gradx
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind
use gridmod, only: lat2,lon2
use jfunc, only: nsclen,npclen,ntclen
use balmod, only: balance,tbalance,ke_vp,bvk
use berror, only: varprd,fpsproj,fut2ps
use constants, only: zero
use control_vectors, only: control_vector,assignment(=)
use control_vectors, only: mvars,nrf,nrf_var,nrf_3d
use gsi_4dvar, only: nsubwin
use timermod, only: timer_ini,timer_fnl
use gsi_bundlemod, only: gsi_bundlegetpointer,gsi_bundlemerge,gsi_bundle,gsi_bundledup,gsi_bundledestroy
implicit none
! Declare passed variables
type(control_vector),intent(inout) :: gradx
type(control_vector),intent(inout) :: grady
! Declare local variables
integer(i_kind) i,j,k,ii,istatus
real(r_kind),dimension(:,:,:),pointer::p_t =>NULL()
real(r_kind),dimension(:,:,:),pointer::p_sf =>NULL()
real(r_kind),dimension(:,:,:),pointer::p_vp =>NULL()
real(r_kind),dimension(:,: ),pointer::p_ps =>NULL()
real(r_kind),dimension(:,:,:),pointer::p_sfwter =>NULL()
real(r_kind),dimension(:,:,:),pointer::p_vpwter =>NULL()
real(r_kind),pointer::rank2a(:,:) =>NULL()
real(r_kind),pointer::rank2b(:,:) =>NULL()
real(r_kind),pointer::rank3a(:,:,:)=>NULL()
real(r_kind),pointer::rank3b(:,:,:)=>NULL()
logical do_balance, do_balancewter
integer(i_kind), parameter :: myvars = 6
integer(i_kind) :: ipnts(myvars)
type(gsi_bundle) :: mbundle
character(len=6), parameter :: myvnames(myvars) = (/ &
'sf ', 'vp ', 'ps ', 't ', 'sfwter', 'vpwter'/)
! Initialize timer
call timer_ini('anbkerror')
! Put things in grady first since operations change input variables
grady=gradx
! Since each internal vector [step(jj)] of grad has the same structure, pointers
! are the same independent of the subwindow jj
call gsi_bundlegetpointer (grady%step(1),myvnames,ipnts,istatus)
! Define what to do depending on what's in CV and SV
do_balance=ipnts(1)>0.and.ipnts(2)>0.and.ipnts(3)>0.and.ipnts(4)>0
do_balancewter=ipnts(5)>0.and.ipnts(6)>0
! Loop on control steps
do ii=1,nsubwin
! Create temporary bundle which merges grady%step(ii) with grady%motley(ii)
if(mvars>0) then
call gsi_bundlemerge(mbundle,grady%step(ii),grady%motley(ii),' add motley to step',istatus)
else
call gsi_bundledup(grady%step(ii),mbundle,' copy of step ',istatus)
end if
! Transpose of balance equation
if(do_balancewter) then
call gsi_bundlegetpointer (mbundle,'sfwter',p_sfwter, istatus)
call gsi_bundlegetpointer (mbundle,'vpwter',p_vpwter, istatus)
!Adjoint of contribution to velocity potential from streamfunction.
do k=1,ke_vp
do j=1,lon2
do i=1,lat2
p_sfwter(i,j,k)=p_sfwter(i,j,k)+bvk(i,j,k)*p_vpwter(i,j,k)
end do
end do
end do
endif
if(do_balance) then
call gsi_bundlegetpointer (mbundle,'sf',p_sf, istatus)
call gsi_bundlegetpointer (mbundle,'vp',p_vp, istatus)
call gsi_bundlegetpointer (mbundle,'ps',p_ps, istatus)
call gsi_bundlegetpointer (mbundle,'t ',p_t, istatus)
call gsi_bundlegetpointer (mbundle,'sfwter',p_sfwter, istatus)
call gsi_bundlegetpointer (mbundle,'vpwter',p_vpwter, istatus)
call tbalance(p_t,p_ps,p_sf,p_vp,fpsproj,fut2ps)
endif
! Apply variances, as well as vertical & horizontal parts of background error
call anbkgcov(mbundle)
! Balance equation
if(do_balance) call balance(p_t,p_ps,p_sf,p_vp,fpsproj,fut2ps)
if(do_balancewter) then
do k=1,ke_vp
do j=1,lon2
do i=1,lat2
p_vpwter(i,j,k)=p_vpwter(i,j,k)+bvk(i,j,k)*p_sfwter(i,j,k)
end do
end do
end do
endif
! Transfer step part of mbundle back to grady%step(ii)
do i=1,nrf
if(nrf_3d(i)) then
call gsi_bundlegetpointer(mbundle,trim(nrf_var(i)),rank3a,istatus)
call gsi_bundlegetpointer(grady%step(ii),trim(nrf_var(i)),rank3b,istatus)
rank3b=rank3a
else
call gsi_bundlegetpointer(mbundle,trim(nrf_var(i)),rank2a,istatus)
call gsi_bundlegetpointer(grady%step(ii),trim(nrf_var(i)),rank2b,istatus)
rank2b=rank2a
end if
end do
! Transfer motley part of mbundle back to grady%motley(ii)
do i=1,mvars
call gsi_bundlegetpointer(mbundle,trim(nrf_var(nrf+i)),rank2a,istatus)
call gsi_bundlegetpointer(grady%motley(ii),trim(nrf_var(nrf+i)),rank2b,istatus)
rank2b=rank2a
end do
end do
! clean work space
call gsi_bundledestroy(mbundle,istatus)
if(istatus/=0) then
write(6,*) ' in anbkerror: trouble destroying work mbundle'
call stop2(999)
endif
! Take care of background error for bias correction terms
do i=1,nsclen
grady%predr(i)=grady%predr(i)*varprd(i)
end do
do i=1,npclen
grady%predp(i)=grady%predp(i)*varprd(nsclen+i)
end do
if(ntclen>0)then
do i=1,ntclen
grady%predt(i)=grady%predt(i)*varprd(nsclen+npclen+i)
end do
end if
! Finalize timer
call timer_fnl('anbkerror')
end subroutine anbkerror
subroutine anbkgcov(bundle)
!$$$ subprogram documentation block
! . . . .
! subprogram: anbkgcov apply anisotropic background error covar
! prgmmr: parrish org: np22 date: 2005-02-14
!
! abstract: apply regional anisotropic background error covariance
!
! program history log:
! 2005-02-14 parrish
! 2009-07-01 sato - update for global mode
! 2010-05-20 todling - update fit interface to sug2grid/grid2sub (bundle)
! 2010-06-22 todling - update interface (remove cwmr since it's in bunlde)
! 2010-06-29 lueken - added if(ipnts(2)>0) to second call of anbkgvar
! 2011-02-22 zhu - replace the argument list of ansmoothrf_reg_subdomain_option by a bundle
! 2012-06-25 parrish - replace sub2grid and grid2sub calls with general_sub2grid, general_grid2sub.
! NOTE: This will not work with sst and the motley variables stl,sti. However
! this is not currently used in this version of RTMA.
! 2014-03-19 pondeca - add wspd10m
! 2014-04-10 pondeca - add td2m,mxtm,mitm,pmsl
! 2014-05-07 pondeca - add howv
! 2016-05-03 pondeca - add uwnd10m, vwnd10m
!
! input argument list:
! t - t on subdomain
! p - p surface pressure on subdomain
! q - q on subdomain
! oz - ozone on subdomain
! skint - skin temperature on subdomain
! sst - sea surface temperature on subdomain
! stl - land surface temperature on subdomain
! sti - ice surface temperature on subdomain
! gust - 10-m gust on subdomain
! vis - surface visibility on subdomain
! wspd10m - 10m-wind speed on subdomain
! td2m - td on subdomain
! mxtm - daily maximum temperature
! mitm - daily minimum temperature
! pmsl - pressure at mean sea level
! howv - significant wave height
! tcamt - total cloud amount
! lcbas - lowest cloud base height
! cldch - cloud ceiling height
! uwnd10m - 10m-uwind on subdomain
! vwnd10m - 10m-vwind on subdomain
! pswter - water surface pressure on subdomain
! twter - water 2m-temperature on subdomain
! qwter - water 2m-specific humidity on subdomain
! gustwter - water 10m-gust on subdomain
! wspd10mwter - water 10m-wind speed on subdomain
! td2mwter - water td on subdomain
! mxtmwter - water daily maximum temperature
! mitmwter - water daily minimum temperature
! uwnd10mwter - water 10m-uwind on subdomain
! vwnd10mwter - water 10m-vwind on subdomain
!
! output argument list:
! all after smoothing, combining scales
! t - t on subdomain
! p - p surface pressure on subdomain
! q - q on subdomain
! oz - ozone on subdomain
! skint - skin temperature on subdomain
! sst - sea surface temperature on subdomain
! stl - land surface temperature on subdomain
! sti - ice surface temperature on subdomain
! gust - 10-m gust on subdomain
! vis - surface visibility on subdomain
! wspd10m - 10m-wind speed on subdomain
! td2m - td on subdomain
! mxtm - daily maximum temperature
! mitm - daily minimum temperature
! pmsl - pressure at mean sea level
! howv - significant wave height
! tcamt - total cloud amount
! lcbas - lowest cloud base height
! cldch - cloud ceiling height
! uwnd10m - 10m-uwind on subdomain
! vwnd10m - 10m-vwind on subdomain
! pswter - water surface pressure on subdomain
! twter - water 2m-temperature on subdomain
! qwter - water 2m-specific humidity on subdomain
! gustwter - water 10m-gust on subdomain
! wspd10mwter - water 10m-wind speed on subdomain
! td2mwter - water td on subdomain
! mxtmwter - water daily maximum temperature
! mitmwter - water daily minimum temperature
! uwnd10mwter - water 10m-uwind on subdomain
! vwnd10mwter - water 10m-vwind on subdomain
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!$$$
use kinds, only: r_kind,i_kind
use gridmod, only: lat2,lon2,twodvar_regional
use anberror, only: rtma_subdomain_option,nsmooth, nsmooth_shapiro,rtma_bkerr_sub2slab
use constants, only: zero
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
use general_sub2grid_mod, only: general_sub2grid,general_grid2sub
implicit none
! Passed Variables
type(gsi_bundle), intent(inout) :: bundle
! Local Variables
integer(i_kind) n,istatus
integer(i_kind) i_sst,i_stl,i_sti,i_ps,i_t,i_q,i_gust,i_wspd10m, &
i_td2m,i_mxtm,i_mitm,i_uwnd10m,i_vwnd10m, &
i_pswter,i_twter,i_qwter,i_gustwter,i_wspd10mwter, &
i_td2mwter,i_mxtmwter,i_mitmwter,i_uwnd10mwter,i_vwnd10mwter
real(r_kind),dimension(lat2,lon2) :: skint,sst,stl,sti
real(r_kind),dimension(lat2,lon2) :: field,fld,fldwter
real(r_kind),dimension(:,:), pointer :: ptrsst=>NULL()
real(r_kind),dimension(:,:), pointer :: ptrstl=>NULL()
real(r_kind),dimension(:,:), pointer :: ptrsti=>NULL()
real(r_kind),dimension(:,:), pointer :: ptrps=>NULL()
real(r_kind),dimension(:,:,:), pointer :: ptrt=>NULL()
real(r_kind),dimension(:,:,:), pointer :: ptrq=>NULL()
real(r_kind),dimension(:,:), pointer :: ptrgust=>NULL()
real(r_kind),dimension(:,:), pointer :: ptrwspd10m=>NULL()
real(r_kind),dimension(:,:), pointer :: ptrtd2m=>NULL()
real(r_kind),dimension(:,:), pointer :: ptrmxtm=>NULL()
real(r_kind),dimension(:,:), pointer :: ptrmitm=>NULL()
real(r_kind),dimension(:,:), pointer :: ptruwnd10m=>NULL()
real(r_kind),dimension(:,:), pointer :: ptrvwnd10m=>NULL()
real(r_kind),dimension(:,:,:), pointer :: ptr3d=>NULL()
real(r_kind),dimension(:,:), pointer :: ptrpswter=>NULL()
real(r_kind),dimension(:,:), pointer :: ptrtwter=>NULL()
real(r_kind),dimension(:,:), pointer :: ptrqwter=>NULL()
real(r_kind),dimension(:,:), pointer :: ptrgustwter=>NULL()
real(r_kind),dimension(:,:), pointer :: ptrwspd10mwter=>NULL()
real(r_kind),dimension(:,:), pointer :: ptrtd2mwter=>NULL()
real(r_kind),dimension(:,:), pointer :: ptrmxtmwter=>NULL()
real(r_kind),dimension(:,:), pointer :: ptrmitmwter=>NULL()
real(r_kind),dimension(:,:), pointer :: ptruwnd10mwter=>NULL()
real(r_kind),dimension(:,:), pointer :: ptrvwnd10mwter=>NULL()
! Perform simple vertical smoothing while fields are in sudomain mode.
! The accompanying smoothing in the horizontal is performed inside the
! recursive filter. Motivation: Reduce possible high frequency noise in
! the analysis that would arise from the use of a "non-blending" RF algorithm.
do n=1,bundle%n3d
call gsi_bundlegetpointer ( bundle,bundle%r3(n)%shortname,ptr3d,istatus )
call vert_smther(ptr3d,nsmooth,nsmooth_shapiro)
end do
! Break up skin temp into components
! First get pointers
!SST
call gsi_bundlegetpointer (bundle, 'sst', i_sst, istatus)
call gsi_bundlegetpointer (bundle, 'stl', i_stl, istatus)
call gsi_bundlegetpointer (bundle, 'sti', i_sti, istatus)
if (i_sst>0) call gsi_bundlegetpointer (bundle, 'sst', ptrsst, istatus)
if (i_stl>0) call gsi_bundlegetpointer (bundle, 'stl', ptrstl, istatus)
if (i_sti>0) call gsi_bundlegetpointer (bundle, 'sti', ptrsti, istatus)
if(twodvar_regional) then
!SURFACE PRESSURE
call gsi_bundlegetpointer (bundle, 'ps', i_ps, istatus)
call gsi_bundlegetpointer (bundle, 'pswter', i_pswter, istatus)
if (i_ps>0) call gsi_bundlegetpointer (bundle, 'ps', ptrps, istatus)
if (i_pswter>0) call gsi_bundlegetpointer (bundle, 'pswter', ptrpswter, istatus)
!TEMPERATURE
call gsi_bundlegetpointer (bundle, 't', i_t, istatus)
call gsi_bundlegetpointer (bundle, 'twter', i_twter, istatus)
if (i_t>0) call gsi_bundlegetpointer (bundle, 't', ptrt, istatus)
if (i_twter>0) call gsi_bundlegetpointer (bundle, 'twter', ptrtwter, istatus)
!SPECIFIC HUMIDITY
call gsi_bundlegetpointer (bundle, 'q', i_q, istatus)
call gsi_bundlegetpointer (bundle, 'qwter', i_qwter, istatus)
if (i_q>0) call gsi_bundlegetpointer (bundle, 'q', ptrq, istatus)
if (i_qwter>0) call gsi_bundlegetpointer (bundle, 'qwter', ptrqwter, istatus)
!WIND GUST
call gsi_bundlegetpointer (bundle, 'gust', i_gust, istatus)
call gsi_bundlegetpointer (bundle, 'gustwter', i_gustwter, istatus)
if (i_gust>0) call gsi_bundlegetpointer (bundle, 'gust', ptrgust, istatus)
if (i_gustwter>0) call gsi_bundlegetpointer (bundle, 'gustwter', ptrgustwter, istatus)
!10-m WIND SPEED
call gsi_bundlegetpointer (bundle, 'wspd10m', i_wspd10m, istatus)
call gsi_bundlegetpointer (bundle, 'wspd10mwter', i_wspd10mwter, istatus)
if (i_wspd10m>0) call gsi_bundlegetpointer (bundle, 'wspd10m', ptrwspd10m, istatus)
if (i_wspd10mwter>0) call gsi_bundlegetpointer (bundle, 'wspd10mwter', ptrwspd10mwter, istatus)
!2-m DEW POINT
call gsi_bundlegetpointer (bundle, 'td2m', i_td2m, istatus)
call gsi_bundlegetpointer (bundle, 'td2mwter', i_td2mwter, istatus)
if (i_td2m>0) call gsi_bundlegetpointer (bundle, 'td2m', ptrtd2m, istatus)
if (i_td2mwter>0) call gsi_bundlegetpointer (bundle, 'td2mwter', ptrtd2mwter, istatus)
!MAXIMUM TEMPERATURE
call gsi_bundlegetpointer (bundle, 'mxtm', i_mxtm, istatus)
call gsi_bundlegetpointer (bundle, 'mxtmwter', i_mxtmwter, istatus)
if (i_mxtm>0) call gsi_bundlegetpointer (bundle, 'mxtm', ptrmxtm, istatus)
if (i_mxtmwter>0) call gsi_bundlegetpointer (bundle, 'mxtmwter', ptrmxtmwter, istatus)
!MINIMUM TEMPERATURE
call gsi_bundlegetpointer (bundle, 'mitm', i_mitm, istatus)
call gsi_bundlegetpointer (bundle, 'mitmwter', i_mitmwter, istatus)
if (i_mitm>0) call gsi_bundlegetpointer (bundle, 'mitm', ptrmitm, istatus)
if (i_mitmwter>0) call gsi_bundlegetpointer (bundle, 'mitmwter', ptrmitmwter, istatus)
!10-m UWIND
call gsi_bundlegetpointer (bundle, 'uwnd10m', i_uwnd10m, istatus)
call gsi_bundlegetpointer (bundle, 'uwnd10mwter', i_uwnd10mwter, istatus)
if (i_uwnd10m>0) call gsi_bundlegetpointer (bundle, 'uwnd10m', ptruwnd10m, istatus)
if (i_uwnd10mwter>0) call gsi_bundlegetpointer (bundle, 'uwnd10mwter', ptruwnd10mwter, istatus)
!10-m VWIND
call gsi_bundlegetpointer (bundle, 'vwnd10m', i_vwnd10m, istatus)
call gsi_bundlegetpointer (bundle, 'vwnd10mwter', i_vwnd10mwter, istatus)
if (i_vwnd10m>0) call gsi_bundlegetpointer (bundle, 'vwnd10m', ptrvwnd10m, istatus)
if (i_vwnd10mwter>0) call gsi_bundlegetpointer (bundle, 'vwnd10mwter', ptrvwnd10mwter, istatus)
endif
! Break up skin temp
if(i_sst>0 .and. i_stl>0 .and. i_sti>0) then
stl=zero
sst=zero
sti=zero
skint=ptrsst
call anbkgvar(skint,sst,stl,sti,0)
ptrsst=sst
ptrstl=stl
ptrsti=sti
endif
if(twodvar_regional) then
if(i_ps>0.and.i_pswter>0) then
fld=zero
fldwter=zero
field=ptrps
call anbkgvar_lw(field,fld,fldwter,0)
ptrps=fld
ptrpswter=fldwter
endif
if(i_t>0.and.i_twter>0) then
fld=zero
fldwter=zero
field(:,:)=ptrt(:,:,1)
call anbkgvar_lw(field,fld,fldwter,0)
ptrt(:,:,1)=fld(:,:)
ptrtwter=fldwter
endif
if(i_q>0.and.i_qwter>0) then
fld=zero
fldwter=zero
field(:,:)=ptrq(:,:,1)
call anbkgvar_lw(field,fld,fldwter,0)
ptrq(:,:,1)=fld(:,:)
ptrqwter=fldwter
endif
if(i_gust>0.and.i_gustwter>0) then
fld=zero
fldwter=zero
field=ptrgust
call anbkgvar_lw(field,fld,fldwter,0)
ptrgust=fld
ptrgustwter=fldwter
endif
if(i_wspd10m>0.and.i_wspd10mwter>0) then
fld=zero
fldwter=zero
field=ptrwspd10m
call anbkgvar_lw(field,fld,fldwter,0)
ptrwspd10m=fld
ptrwspd10mwter=fldwter
endif
if(i_td2m>0.and.i_td2mwter>0) then
fld=zero
fldwter=zero
field=ptrtd2m
call anbkgvar_lw(field,fld,fldwter,0)
ptrtd2m=fld
ptrtd2mwter=fldwter
endif
if(i_mxtm>0.and.i_mxtmwter>0) then
fld=zero
fldwter=zero
field=ptrmxtm
call anbkgvar_lw(field,fld,fldwter,0)
ptrmxtm=fld
ptrmxtmwter=fldwter
endif
if(i_mitm>0.and.i_mitmwter>0) then
fld=zero
fldwter=zero
field=ptrmitm
call anbkgvar_lw(field,fld,fldwter,0)
ptrmitm=fld
ptrmitmwter=fldwter
endif
if(i_uwnd10m>0.and.i_uwnd10mwter>0) then
fld=zero
fldwter=zero
field=ptruwnd10m
call anbkgvar_lw(field,fld,fldwter,0)
ptruwnd10m=fld
ptruwnd10mwter=fldwter
endif
if(i_vwnd10m>0.and.i_vwnd10mwter>0) then
fld=zero
fldwter=zero
field=ptrvwnd10m
call anbkgvar_lw(field,fld,fldwter,0)
ptrvwnd10m=fld
ptrvwnd10mwter=fldwter
endif
endif
! Apply auto-covariance
if(rtma_subdomain_option) then
if(rtma_bkerr_sub2slab) then
call ansmoothrf_reg_sub2slab_option(bundle)
else
call ansmoothrf_reg_subdomain_option(bundle)
end if
else
call ansmoothrf(bundle)
end if
!Adjoint of simple vertical smoothing
do n=bundle%n3d,1,-1
call gsi_bundlegetpointer ( bundle,bundle%r3(n)%shortname,ptr3d,istatus )
call tvert_smther(ptr3d,nsmooth,nsmooth_shapiro)
end do
if(twodvar_regional) then
if(i_vwnd10m>0.and.i_vwnd10mwter>0) then
fld=ptrvwnd10m
fldwter=ptrvwnd10mwter
field=zero
call anbkgvar_lw(field,fld,fldwter,1)
ptrvwnd10m=field
! ignore content of remaining arrays
endif
if(i_uwnd10m>0.and.i_uwnd10mwter>0) then
fld=ptruwnd10m
fldwter=ptruwnd10mwter
field=zero
call anbkgvar_lw(field,fld,fldwter,1)
ptruwnd10m=field
! ignore content of remaining arrays
endif
if(i_mitm>0.and.i_mitmwter>0) then
fld=ptrmitm
fldwter=ptrmitmwter
field=zero
call anbkgvar_lw(field,fld,fldwter,1)
ptrmitm=field
! ignore content of remaining arrays
endif
if(i_mxtm>0.and.i_mxtmwter>0) then
fld=ptrmxtm
fldwter=ptrmxtmwter
field=zero
call anbkgvar_lw(field,fld,fldwter,1)
ptrmxtm=field
! ignore content of remaining arrays
endif
if(i_td2m>0.and.i_td2mwter>0) then
fld=ptrtd2m
fldwter=ptrtd2mwter
field=zero
call anbkgvar_lw(field,fld,fldwter,1)
ptrtd2m=field
! ignore content of remaining arrays
endif
if(i_wspd10m>0.and.i_wspd10mwter>0) then
fld=ptrwspd10m
fldwter=ptrwspd10mwter
field=zero
call anbkgvar_lw(field,fld,fldwter,1)
ptrwspd10m=field
! ignore content of remaining arrays
endif
if(i_gust>0.and.i_gustwter>0) then
fld=ptrgust
fldwter=ptrgustwter
field=zero
call anbkgvar_lw(field,fld,fldwter,1)
ptrgust=field
! ignore content of remaining arrays
endif
if(i_q>0.and.i_qwter>0) then
fld(:,:)=ptrq(:,:,1)
fldwter=ptrqwter
field=zero
call anbkgvar_lw(field,fld,fldwter,1)
ptrq(:,:,1)=field(:,:)
! ignore content of remaining arrays
endif
if(i_t>0.and.i_twter>0) then
fld(:,:)=ptrt(:,:,1)
fldwter=ptrtwter
field(:,:)=zero
call anbkgvar_lw(field,fld,fldwter,1)
ptrt(:,:,1)=field(:,:)
! ignore content of remaining arrays
endif
if(i_ps>0.and.i_pswter>0) then
fld=ptrps
fldwter=ptrpswter
field=zero
call anbkgvar_lw(field,fld,fldwter,1)
ptrps=field
! ignore content of remaining arrays
endif
endif
!==> combine sst,stl, and sti into skin temperature field
if(i_sst>0 .and. i_stl>0 .and. i_sti>0) then
sst=ptrsst
stl=ptrstl
sti=ptrsti
skint=zero
call anbkgvar(skint,sst,stl,sti,1)
ptrsst=skint
! ignore contents of remaining arrays
endif
end subroutine anbkgcov
subroutine anbkgvar(skint,sst,stl,sti,iflg)
!$$$ subprogram documentation block
! . . . .
! subprogram: anbkgvar_reg manipulate skin temp
! prgmmr: parrish org: np22 date: 1990-10-06
!
! abstract: manipulate skin temp <--> sst,sfc temp, and ice temp fields
!
! program history log:
! 2005-01-22 parrish
! 2008-06-05 safford - rm unused uses
! 2012-06-25 parrish - remove _i_kind from integer constants
!
! input argument list:
! skint - skin temperature grid values
! sst - sst grid values
! stl - land surface temperature grid values
! sti - snow/ice covered surface temperature grid values
! iflg - flag for skin temperature manipulation
! 0: skint --> sst,stl,sti
! 1: sst,stl,sti --> skint
!
! output argument list:
! skint - skin temperature grid values
! sst - sst grid values
! stl - land surface temperature grid values
! sti - snow/ice covered surface temperature grid values
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use gridmod, only: lat2,lon2
use guess_grids, only: isli2
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: iflg
real(r_kind),dimension(lat2,lon2),intent(inout) :: skint,sst,stl,sti
! Declare local variables
integer(i_kind) i,j
do j=1,lon2
do i=1,lat2
if(iflg == 0) then
! Break skin temperature into components
! If land point
if(isli2(i,j) == 1) then
stl(i,j)=skint(i,j)
! If ice
else if(isli2(i,j) == 2) then
sti(i,j)=skint(i,j)
! Else treat as a water point
else
sst(i,j)=skint(i,j)
end if
else if (iflg==1) then
! Combine sst,stl, and sti into skin temperature field
! Land point, load land sfc t into skint
if(isli2(i,j) == 1) then
skint(i,j)=stl(i,j)
! Ice, load ice temp into skint
else if(isli2(i,j) == 2) then
skint(i,j)=sti(i,j)
! Treat as a water point, load sst into skint
else
skint(i,j)=sst(i,j)
end if
end if
end do
end do
return
end subroutine anbkgvar
subroutine anbkgvar_lw_original(field,fld,fldwter,iflg)
!$$$ subprogram documentation block
! . . .
! subprogram: anbkgvar_lw land/water field manipulation
! prgmmr: pondeca org: np22 date: 2013-04-02
!
! abstract: manipulate field <--> fld,fldwter
! (based on anbkgvar)
!
! program history log:
! 2013-04-02 pondeca
!
! input argument list:
! field - field grid values
! fld - field grid values over land
! fldwter - field grid values over water
! iflg - flag for field manipulation
! 0: field --> fld,fldwter
! 1: fld,fldwter --> field
!
! output argument list:
! field - field grid values
! fld - field grid values over land
! fldwter - field grid values over water
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use gridmod, only: lat2,lon2
use guess_grids, only: isli2
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: iflg
real(r_kind),dimension(lat2,lon2),intent(inout) :: field, &
fld,fldwter
! Declare local variables
integer(i_kind) i,j
do j=1,lon2
do i=1,lat2
if(iflg == 0) then
! Break field into components
! If land point
if(isli2(i,j) == 1) then
fld(i,j)=field(i,j)
! Else treat as a water point
else
fldwter(i,j)=field(i,j)
end if
else if (iflg==1) then
! Combine fld,fldwter into field
! Land point
if(isli2(i,j) == 1) then
field(i,j)=fld(i,j)
! Treat as a water point,
else
field(i,j)=fldwter(i,j)
end if
end if
end do
end do
return
end subroutine anbkgvar_lw_original
subroutine anbkgvar_lw(field,fld,fldwter,iflg)
!$$$ subprogram documentation block
! . . .
! subprogram: anbkgvar_lw land/water field manipulation
! prgmmr: pondeca org: np22 date: 2013-04-02
!
! abstract: manipulate field <--> fld,fldwter
! (based on anbkgvar)
!
! program history log:
! 2013-04-02 pondeca
!
! input argument list:
! field - field grid values
! fld - field grid values over land
! fldwter - field grid values over water
! iflg - flag for field manipulation
! 0: field --> fld,fldwter
! 1: fld,fldwter --> field
!
! output argument list:
! field - field grid values
! fld - field grid values over land
! fldwter - field grid values over water
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use gridmod, only: lat2,lon2,region_lat,region_lon, &
nlon_regional,nlat_regional,istart,jstart
use guess_grids, only: isli2
use mpimod, only: mype
use constants, only: rad2deg
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: iflg
real(r_kind),dimension(lat2,lon2),intent(inout) :: field, &
fld,fldwter
! Declare local variables
integer(i_kind) i,j
integer(i_kind) iglob,jglob
integer(i_kind) mm1
real(r_kind) flon,flat
logical glerlarea
! Declare local parameters
! Great Lakes
real(r_kind),parameter::flon1=-93._r_kind
real(r_kind),parameter::flon2=-75._r_kind
real(r_kind),parameter::flat1=40.5_r_kind
real(r_kind),parameter::flat2=49.5_r_kind
! Great Salt Lake
real(r_kind),parameter::slon1=-113._r_kind
real(r_kind),parameter::slon2=-112._r_kind
real(r_kind),parameter::slat1=40.6_r_kind
real(r_kind),parameter::slat2=41.7_r_kind
mm1=mype+1
do j=1,lon2
jglob=max(1,min(j+jstart(mm1)-2,nlon_regional))
do i=1,lat2
iglob=max(1,min(i+istart(mm1)-2,nlat_regional))
flat=region_lat(iglob,jglob)*rad2deg
flon=region_lon(iglob,jglob)*rad2deg
glerlarea=(flat>=flat1.and.flat<=flat2).and.(flon>=flon1.and.flon<=flon2)
glerlarea=glerlarea.or.((flat>=slat1.and.flat<=slat2).and.(flon>=slon1.and.flon<=slon2))
if(iflg == 0) then
! Break field into components
! If land point
if(isli2(i,j) == 1 .or. .not.glerlarea) then
fld(i,j)=field(i,j)
! Else treat as a water point
else
fldwter(i,j)=field(i,j)
end if
else if (iflg==1) then
! Combine fld,fldwter into field
! Land point
if(isli2(i,j) == 1 .or. .not.glerlarea) then
field(i,j)=fld(i,j)
! Treat as a water point,
else
field(i,j)=fldwter(i,j)
end if
end if
end do
end do
return
end subroutine anbkgvar_lw
subroutine ansmoothrf(cstate)
!$$$ subprogram documentation block
! . . . .
! subprogram: ansmoothrf anisotropic rf for regional mode
! prgmmr: parrish org: np22 date: 2005-02-14
!
! abstract: apply anisotropic rf for regional mode
!
! program history log:
! 2005-02-14 parrish
! 2008-12-04 sato - update for global mode
! 2009-01-02 todling - get mype from mpimod directly
! 2015-07-01 pondeca - rewrite to handle motley variables
!
! input argument list:
! cstate - bundle containing horizontal fields to be smoothed
!
! output argument list:
! cstate - bundle containing smoothed horizontal fields
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!$$$
use kinds, only: r_kind,i_kind,r_single
use anberror, only: indices,indices_p,ngauss,pf2aP1, &
filter_all,filter_p2,filter_p3
use patch2grid_mod, only: patch2grid, tpatch2grid
use mpimod, only: npe
use constants, only: zero
use gridmod, only: nlat,nlon,lat2,lon2,nsig,regional
use control_vectors, only: nvars,nrf,nrf_var,nrf_3d
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
use general_commvars_mod, only: s2g_raf
use general_sub2grid_mod, only: general_sub2grid,general_grid2sub
use fgrid2agrid_mod, only: fgrid2agrid,tfgrid2agrid
use raflib, only: raf4_ad,raf4
use raflib, only: raf4_ad_wrap,raf4_wrap
implicit none
! Declare passed variables
type(gsi_bundle),intent(inout) :: cstate
! Declare local variables
integer(i_kind):: ips,ipe,jps,jpe,kps,kpe,kds,kde
integer(i_kind):: p_ips,p_ipe,p_jps,p_jpe,p_kps,p_kpe
integer(i_kind) i,igauss,j,k,kk,n,istatus
real(r_kind),allocatable:: fields(:,:,:,:)
real(r_kind),allocatable:: work(:,:,:,:)
real(r_kind),allocatable:: worka(:,:,:)
real(r_single),allocatable:: workb(:,:,:,:)
real(r_kind),allocatable,dimension(:,:,:) :: workanp,workasp
real(r_single),allocatable,dimension(:,:,:,:):: workbnp,workbsp
real(r_kind),pointer::rank2(:,:)
real(r_kind),pointer::rank3(:,:,:)
! Convert from subdomain to full horizontal field distributed among processors
kds=1
kde=s2g_raf%num_fields
kps=s2g_raf%kbegin_loc
kpe=s2g_raf%kend_loc
ips=indices%ips
ipe=indices%ipe
jps=indices%jps
jpe=indices%jpe
p_kps=s2g_raf%kbegin_loc !indices_p%kps
p_kpe=s2g_raf%kend_loc !indices_p%kpe
p_ips=indices_p%ips
p_ipe=indices_p%ipe
p_jps=indices_p%jps
p_jpe=indices_p%jpe
!WHY THE "1" IN THE ALLOCATE STATEMENTS? / MPondeca
allocate(fields(1,lat2,lon2,kds:kde)) !MUST DEALLOCATE / MPondeca
allocate(work(1,nlat,nlon,kps:max(kps,kpe))) !MUST DEALLOCATE / MPondeca
allocate(worka(ips:ipe,jps:jpe,kps:max(kps,kpe)))
allocate(workb(ngauss,ips:ipe,jps:jpe,kps:max(kps,kpe)))
if(.not.regional) then
allocate(workanp(p_ips:p_ipe, p_jps:p_jpe, p_kps:p_kpe))
allocate(workasp(p_ips:p_ipe ,p_jps:p_jpe ,p_kps:p_kpe))
allocate(workbnp(ngauss, p_ips:p_ipe, p_jps:p_jpe, p_kps:p_kpe))
allocate(workbsp(ngauss, p_ips:p_ipe, p_jps:p_jpe, p_kps:p_kpe ))
end if
fields=zero
kk=0
do n=1,nvars
if (n<=nrf .and. nrf_3d(min(n,nrf))) then
call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank3,istatus)
if(istatus==0) then
do k=1,nsig
kk=kk+1
do j=1,lon2
do i=1,lat2
fields(1,i,j,kk)=rank3(i,j,k)
end do
end do
end do
endif
else !2d flds including motley flds
call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank2,istatus)
if(istatus==0) then
kk=kk+1
do j=1,lon2
do i=1,lat2
fields(1,i,j,kk)=rank2(i,j)
end do
end do
endif
endif
end do
! Convert from subdomain to full horizontal fields distributed among processors
call general_sub2grid(s2g_raf,fields,work)
! adjoint of coarse to fine grid
do k=kps,kpe
if(regional) then
call tfgrid2agrid(pf2aP1,work(1,1,1,k),worka(ips,jps,k))
else
call tpatch2grid(work(1,1,1,k), &
worka (ips,jps,k), &
workanp(p_ips,p_jps,k), &
workasp(p_ips,p_jps,k))
end if
end do
! transfer coarse grid fields to ngauss copies
do k=kps,kpe
do j=jps,jpe
do i=ips,ipe
do igauss=1,ngauss
workb(igauss,i,j,k)=worka(i,j,k)
end do
end do
end do
end do
if(.not.regional) then
do k=p_kps,p_kpe
do j=p_jps,p_jpe
do i=p_ips,p_ipe
do igauss=1,ngauss
workbnp(igauss,i,j,k)=workanp(i,j,k)
workbsp(igauss,i,j,k)=workasp(i,j,k)
end do
end do
end do
end do
end if
! apply recursive filter
! call raf4_wrap( workb,filter_all,ngauss,indices,npe) /MPondeca
! call raf4_ad_wrap(workb,filter_all,ngauss,indices,npe) /MPondeca
call raf4 (workb,filter_all,ngauss,ips,ipe,jps,jpe,ips,ipe,jps,jpe,kps,kpe,npe)
call raf4_ad(workb,filter_all,ngauss,ips,ipe,jps,jpe,ips,ipe,jps,jpe,kps,kpe,npe)
if(.not.regional) then
call raf4_wrap( workbnp,filter_p2,ngauss,indices_p,npe) !WORK ON THIS / MPondeca
call raf4_ad_wrap(workbnp,filter_p2,ngauss,indices_p,npe) !WORK ON THIS / MPondeca
call raf4_wrap( workbsp,filter_p3,ngauss,indices_p,npe) !WORK ON THIS / MPondeca
call raf4_ad_wrap(workbsp,filter_p3,ngauss,indices_p,npe) !WORK ON THIS / MPondeca
end if
! add together ngauss copies
worka=zero
do k=kps,kpe
do j=jps,jpe
do i=ips,ipe
do igauss=1,ngauss
worka(i,j,k)=worka(i,j,k)+workb(igauss,i,j,k)
end do
end do
end do
end do
deallocate(workb)
if(.not.regional) then
workanp=zero
workasp=zero
do k=p_kps,p_kpe
do j=p_jps,p_jpe
do i=p_ips,p_ipe
do igauss=1,ngauss
workanp(i,j,k)=workanp(i,j,k)+workbnp(igauss,i,j,k)
workasp(i,j,k)=workasp(i,j,k)+workbsp(igauss,i,j,k)
end do
end do
end do
end do
end if
! coarse to fine grid
do k=kps,kpe
if(regional) then
call fgrid2agrid(pf2aP1,worka(ips,jps,k),work(1,1,1,k))
else
call patch2grid(work(1,1,1,k), &
worka (ips,jps,k),&
workanp(p_ips,p_jps,k),&
workasp(p_ips,p_jps,k))
end if
end do
deallocate(worka)
call general_grid2sub(s2g_raf,work,fields)
deallocate(work)
if(.not.regional) then
deallocate(workanp)
deallocate(workbnp)
deallocate(workasp)
deallocate(workbsp)
end if
! transfer from work back to bundle
kk=0
do n=1,nvars
if (n<=nrf .and. nrf_3d(min(n,nrf))) then
call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank3,istatus)
if(istatus==0) then
do k=1,nsig
kk=kk+1
do j=1,lon2
do i=1,lat2
rank3(i,j,k)=fields(1,i,j,kk)
end do
end do
end do
endif
else !2d flds including motley flds
call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank2,istatus)
if(istatus==0) then
kk=kk+1
do j=1,lon2
do i=1,lat2
rank2(i,j)=fields(1,i,j,kk)
end do
end do
endif
endif
end do
deallocate(fields)
end subroutine ansmoothrf
subroutine vert_smther(g,nsmooth,nsmooth_shapiro)
!$$$ subprogram documentation block
! . . . .
! subprogram: vert_smther
! prgmmr:
!
! abstract:
!
! program history log:
! 2009-09-15 lueken - added subprogram doc block
!
! input argument list:
! nsmooth
! nsmooth_shapiro
! g
!
! output argument list:
! g
!
! Notes: nsmooth > 0 ==> apply 1-2-1 smoother
! nsmooth_shapiro 0 ==> apply second moment preserving
! "shapiro smoother"
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
use constants, only: quarter,half
use kinds, only: r_kind,i_kind
use gridmod, only: lat2,lon2,nsig
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: nsmooth,nsmooth_shapiro
real(r_kind),dimension(1:lat2,1:lon2,1:nsig),intent(inout) :: g
! Declare local variables
integer(i_kind) i,j,l,k,kp,km,kp3,km3
real(r_kind), allocatable:: gaux(:)
if (nsig==1)return
allocate(gaux(1:nsig))
if (nsmooth > 0 ) then
do i=1,lat2
do j=1,lon2
do l=1,nsmooth
gaux(1:nsig)=g(i,j,1:nsig)
do k=1,nsig
kp=min(k+1,nsig) ; km=max(1,k-1)
g(i,j,k)=quarter*(gaux(kp)+gaux(km))+half*gaux(k)
enddo
enddo
enddo
enddo
endif
if (nsmooth_shapiro > 0 .and. nsmooth <= 0) then
do i=1,lat2
do j=1,lon2
do l=1,nsmooth_shapiro
gaux(1:nsig)=g(i,j,1:nsig)
do k=1,nsig
kp=min(k+1,nsig) ; km=max(1,k-1)
kp3=min(k+3,nsig) ; km3=max(1,k-3)
g(i,j,k)=.28125_r_kind*(gaux(kp)+gaux(km))+half*gaux(k)-.03125_r_kind*(gaux(kp3)+gaux(km3))
enddo
enddo
enddo
enddo
endif
deallocate(gaux)
return
end subroutine vert_smther
subroutine tvert_smther(g,nsmooth,nsmooth_shapiro)
!$$$ subprogram documentation block
! . . . .
! subprogram: tvert_smther
! prgmmr:
!
! abstract:
!
! program history log:
! 2009-09-15 lueken - added subprogram doc block
!
! input argument list:
! nsmooth
! nsmooth_shapiro
! g
!
! output argument list:
! g
!
! Notes: nsmooth > 0 ==> 1-2-1 smoother
! nsmooth_shapiro 0 ==> second moment preserving
! "shapiro smoother"
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
use constants, only: zero,quarter,half
use kinds, only: r_kind,i_kind
use gridmod, only: lat2,lon2,nsig
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: nsmooth,nsmooth_shapiro
real(r_kind),dimension(1:lat2,1:lon2,1:nsig),intent(inout) :: g
! Declare local variables
integer(i_kind) i,j,l,k,kp,km,kp3,km3
real(r_kind), allocatable:: gaux(:)
if (nsig==1)return
allocate(gaux(1:nsig))
if (nsmooth > 0 ) then
do i=1,lat2
do j=1,lon2
do l=1,nsmooth
gaux(1:nsig)=zero
do k=1,nsig
kp=min(k+1,nsig) ; km=max(1,k-1)
gaux(k)=gaux(k)+half*g(i,j,k)
gaux(km)=gaux(km)+quarter*g(i,j,k)
gaux(kp)=gaux(kp)+quarter*g(i,j,k)
enddo
g(i,j,1:nsig)=gaux(1:nsig)
enddo
enddo
enddo
endif
if (nsmooth_shapiro > 0 .and. nsmooth <= 0) then
do i=1,lat2
do j=1,lon2
do l=1,nsmooth_shapiro
gaux(1:nsig)=zero
do k=1,nsig
kp=min(k+1,nsig) ; km=max(1,k-1)
kp3=min(k+3,nsig) ; km3=max(1,k-3)
gaux(km3)=gaux(km3)-.03125_r_kind*g(i,j,k)
gaux(kp3)=gaux(kp3)-.03125_r_kind*g(i,j,k)
gaux(k)=gaux(k)+half*g(i,j,k)
gaux(km)=gaux(km)+.28125_r_kind*g(i,j,k)
gaux(kp)=gaux(kp)+.28125_r_kind*g(i,j,k)
enddo
g(i,j,1:nsig)=gaux(1:nsig)
enddo
enddo
enddo
endif
deallocate(gaux)
return
end subroutine tvert_smther
subroutine ansmoothrf_reg_sub2slab_option(cstate)
!$$$ subprogram documentation block
! . . . .
! subprogram: ansmoothrf_reg_subdomain_option anisotropic rf for regional mode
! prgmmr: parrish org: np22 date: 2015-05-01
!
! abstract: apply anisotropic rf for regional mode (input subdomains but slab for filter)
!
! program history log:
! 2005-02-14 parrish
! 2011-02-22 zhu - use cstate to replace argument list such as p,t,q,vp,st
! 2014-02-07 pondeca - update to include motley variables as well in 2d set of filtered variables
! 2015-05-02 parrish - make copy of ansmoothrf_reg_subdomain_option and
! modify to allow conversion from subdomains to slabs.
!
! input argument list:
! t,p,q,oz,st,stl,sti,cwmr,st,vp - fields to be smoothed
!
! output argument list:
! t,p,q,oz,st,stl,sti,cwmr,st,vp - smoothed fields
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!$$$
use kinds, only: r_kind,i_kind,r_single
use anberror, only: filter_all,ngauss
use anberror, only: pf2aP1
use mpimod, only: npe
use constants, only: zero,zero_single
use gridmod, only: lat2,lon2,nsig
use raflib, only: raf4_ad,raf4
use control_vectors, only: nvars,nrf,nrf_var,nrf_3d
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
use general_commvars_mod, only: s2g_raf
use general_sub2grid_mod, only: general_sub2grid,general_grid2sub
use fgrid2agrid_mod, only: fgrid2agrid,tfgrid2agrid
implicit none
! Declare passed variables
type(gsi_bundle),intent(inout) :: cstate
! Declare local variables
integer(i_kind) i,igauss,j,k,kk,n,istatus
real(r_kind),allocatable:: worka(:,:,:,:)
real(r_single),allocatable:: workb(:,:,:,:)
real(r_kind),pointer::rank2(:,:)
real(r_kind),pointer::rank3(:,:,:)
real(r_kind),allocatable::slaba(:,:,:,:)
real(r_kind),allocatable::slabb(:,:,:)
integer(i_kind):: kds,kde,kps,kpe,nlatf,nlonf,nlat,nlon
kds=1
kde=s2g_raf%num_fields
kps=s2g_raf%kbegin_loc
kpe=s2g_raf%kend_loc
nlat=pf2ap1%nlata
nlon=pf2ap1%nlona
nlatf=pf2ap1%nlatf
nlonf=pf2ap1%nlonf
allocate(worka(1,lat2,lon2,kds:kde))
allocate(slaba(1,nlat,nlon,kps:max(kps,kpe)))
allocate(slabb(nlatf,nlonf,kps:max(kps,kpe)))
allocate(workb(ngauss,nlatf,nlonf,kps:max(kps,kpe)))
! transfer from bundle to worka
worka=zero
kk=0
do n=1,nvars
if (n<=nrf .and. nrf_3d(min(n,nrf))) then
call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank3,istatus)
if(istatus==0) then
do k=1,nsig
kk=kk+1
do j=1,lon2
do i=1,lat2
worka(1,i,j,kk)=rank3(i,j,k)
end do
end do
end do
endif
else !2d flds including motley flds
call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank2,istatus)
if(istatus==0) then
kk=kk+1
do j=1,lon2
do i=1,lat2
worka(1,i,j,kk)=rank2(i,j)
end do
end do
endif
endif
end do
! transfer from subdomains to slabs
call general_sub2grid(s2g_raf,worka,slaba)
! adjoint of coarse to fine grid interpolation:
do k=kps,kpe
call tfgrid2agrid(pf2aP1,slaba(1,1,1,k),slabb(1,1,k))
end do
! transfer to ngauss copies
do k=kps,kpe
do j=1,nlonf
do i=1,nlatf
do igauss=1,ngauss
workb(igauss,i,j,k)=slabb(i,j,k)
end do
end do
end do
end do
! apply recursive filter
call raf4 (workb,filter_all,ngauss,1,nlatf,1,nlonf,1,nlatf,1,nlonf,kps,kpe,npe)
call raf4_ad(workb,filter_all,ngauss,1,nlatf,1,nlonf,1,nlatf,1,nlonf,kps,kpe,npe)
! add together ngauss copies
slabb=zero
do k=kps,kpe
do j=1,nlonf
do i=1,nlatf
do igauss=1,ngauss
slabb(i,j,k)=workb(igauss,i,j,k)+slabb(i,j,k)
end do
end do
end do
end do
deallocate(workb)
! coarse to fine grid interpolation:
do k=kps,kpe
call fgrid2agrid(pf2aP1,slabb(1,1,k),slaba(1,1,1,k))
end do
deallocate(slabb)
call general_grid2sub(s2g_raf,slaba,worka)
deallocate(slaba)
! transfer from worka back to bundle
kk=0
do n=1,nvars
if (n<=nrf .and. nrf_3d(min(n,nrf))) then
call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank3,istatus)
if(istatus==0) then
do k=1,nsig
kk=kk+1
do j=1,lon2
do i=1,lat2
rank3(i,j,k)=worka(1,i,j,kk)
end do
end do
end do
endif
else !2d flds including motley flds
call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank2,istatus)
if(istatus==0) then
kk=kk+1
do j=1,lon2
do i=1,lat2
rank2(i,j)=worka(1,i,j,kk)
end do
end do
endif
endif
end do
deallocate(worka)
end subroutine ansmoothrf_reg_sub2slab_option
subroutine ansmoothrf_reg_subdomain_option(cstate)
!$$$ subprogram documentation block
! . . . .
! subprogram: ansmoothrf_reg_subdomain_option anisotropic rf for regional mode
! prgmmr: parrish org: np22 date: 2005-02-14
!
! abstract: apply anisotropic rf for regional mode (using subdomains instead of slabs)
! NOTE: only works if filter grid is same as analysis grid
!
! program history log:
! 2005-02-14 parrish
! 2011-02-22 zhu - use cstate to replace argument list such as p,t,q,vp,st
! 2014-02-07 pondeca - update to include motley variables as well in 2d set of filtered variables
!
! input argument list:
! t,p,q,oz,st,stl,sti,cwmr,st,vp - fields to be smoothed
!
! output argument list:
! t,p,q,oz,st,stl,sti,cwmr,st,vp - smoothed fields
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!$$$
use kinds, only: r_kind,i_kind,r_single
use anberror, only: indices, filter_all,ngauss,halo_update_reg
use mpimod, only: mype,npe
use constants, only: zero,zero_single
use gridmod, only: istart,jstart,nsig
use raflib, only: raf4_ad_wrap,raf4_wrap
use control_vectors, only: nvars,nrf,nrf_var,nrf_3d
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
implicit none
! Declare passed variables
type(gsi_bundle),intent(inout) :: cstate
! Declare local variables
integer(i_kind) i,igauss,iloc,j,jloc,k,kk,mm1,n,istatus
real(r_single),dimension(ngauss, &
indices%ips:indices%ipe,&
indices%jps:indices%jpe,&
indices%kps:indices%kpe):: workb
real(r_kind),pointer::rank2(:,:)
real(r_kind),pointer::rank3(:,:,:)
integer(i_kind):: ids,ide,jds,jde,kds,kde,ips,ipe,jps,jpe,kps,kpe
ids=indices%ids; ide=indices%ide
jds=indices%jds; jde=indices%jde
kds=indices%kds; kde=indices%kde
ips=indices%ips; ipe=indices%ipe
jps=indices%jps; jpe=indices%jpe
kps=indices%kps; kpe=indices%kpe
mm1=mype+1
workb=zero_single
! transfer variables to ngauss copies
kk=0
do n=1,nvars
if (n<=nrf .and. nrf_3d(min(n,nrf))) then
call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank3,istatus)
if(istatus==0) then
do k=1,nsig
kk=kk+1
do j=jps,jpe
jloc=j-jstart(mm1)+2
do i=ips,ipe
iloc=i-istart(mm1)+2
do igauss=1,ngauss
workb(igauss,i,j,kk)=rank3(iloc,jloc,k)
end do
end do
end do
end do
endif
else !2d flds including motley flds
call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank2,istatus)
if(istatus==0) then
kk=kk+1
do j=jps,jpe
jloc=j-jstart(mm1)+2
do i=ips,ipe
iloc=i-istart(mm1)+2
do igauss=1,ngauss
workb(igauss,i,j,kk)=rank2(iloc,jloc)
end do
end do
end do
endif
endif
end do
! apply recursive filter
call raf4_wrap(workb,filter_all,ngauss,indices,npe)
call raf4_ad_wrap(workb,filter_all,ngauss,indices,npe)
! add together ngauss copies
kk=0
do n=1,nvars
if (n<=nrf .and. nrf_3d(min(n,nrf))) then
call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank3,istatus)
if(istatus==0) then
do k=1,nsig
kk=kk+1
do j=jps,jpe
jloc=j-jstart(mm1)+2
do i=ips,ipe
iloc=i-istart(mm1)+2
rank3(iloc,jloc,k)=zero
do igauss=1,ngauss
rank3(iloc,jloc,k)=rank3(iloc,jloc,k)+workb(igauss,i,j,kk)
end do
end do
end do
end do
call halo_update_reg(rank3,nsig)
endif
else
call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank2,istatus)
if(istatus==0) then
kk=kk+1
do j=jps,jpe
jloc=j-jstart(mm1)+2
do i=ips,ipe
iloc=i-istart(mm1)+2
rank2(iloc,jloc)=zero
do igauss=1,ngauss
rank2(iloc,jloc)=rank2(iloc,jloc)+workb(igauss,i,j,kk)
end do
end do
end do
call halo_update_reg(rank2,1)
endif
endif
end do
end subroutine ansmoothrf_reg_subdomain_option