subroutine control2model(xhat,sval,bval)
!$$$ subprogram documentation block
! . . . .
! subprogram: control2model
! prgmmr: tremolet
!
! abstract: Converts control variable to physical space
!
! program history log:
! 2007-04-13 tremolet - initial code
! 2007-04-27 tremolet - multiply by sqrt(B) (from ckerror D. Parrish)
! 2008-12-04 todling - update interface to ckgcov; add tsen/p3d
! 2008-12-29 todling - add call to strong balance contraint
! 2009-04-21 derber - modify call to getuv to getuv(*,0)
! 2009-08-14 lueken - update documentation
! 2010-03-15 zhu - make changes to ckbcov, add assign_cs2array
! 2010-04-28 todling - update to use gsi_bundle
! 2010-05-31 todling - better consistency checks; add co/co2
! - ready to bypass analysis of (any) meteorological fields
! 2010-06-15 todling - generalized handling of chemistry
! 2011-05-15 auligne/todling - generalized cloud handling
! 2011-07-12 zhu - add cw_to_hydro and cwhydromod
! 2011-12-14 mkim - changed clouds4crtm to clouds in metguess
! 2012-06-25 parrish - modify wbundle by adding motley variables to control vector
! so will be in form expected by new version of ckgcov which uses general_grid2sub.
! 2012-10-09 wgu - update interface to tbalance (fut2ps)
! 2013-10-25 todling - nullify work pointers
! 2013-10-28 todling - rename p3d to prse
! 2013-05-23 zhu - add ntclen and predt for aircraft temperature bias correction
!
! input argument list:
! xhat - Control variable
! sval - State variable
! bval - Bias predictors
!
! output argument list:
! sval - State variable
! bval - Bias predictors
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind
use control_vectors, only: control_vector
use control_vectors, only: cvars3d,cvars2d,cvarsmd
use control_vectors, only: nc2d,mvars
use bias_predictors, only: predictors
use gsi_4dvar, only: nsubwin, l4dvar, lsqrtb
use gridmod, only: lat2,lon2,nsig
use jfunc, only: nsclen,npclen,ntclen
use berror, only: varprd,fpsproj,fut2ps
use balmod, only: balance
use cwhydromod, only: cw2hydro_tl
use gsi_bundlemod, only: gsi_bundlecreate
use gsi_bundlemod, only: gsi_gridcreate
use gsi_bundlemod, only: gsi_grid
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
use gsi_bundlemod, only: gsi_bundlegetvar
use gsi_bundlemod, only: gsi_bundledestroy
use gsi_chemguess_mod, only: gsi_chemguess_get
use gsi_metguess_mod, only: gsi_metguess_get
use mpeu_util, only: getindex
use constants, only: max_varname_length
implicit none
! Declare passed variables
type(control_vector), intent(in ) :: xhat
type(gsi_bundle) , intent(inout) :: sval(nsubwin)
type(predictors) , intent(inout) :: bval
! Declare local variables
character(len=*),parameter:: myname ='control2model'
real(r_kind),dimension(lat2,lon2,nsig) :: workst,workvp,workrh
type(gsi_bundle) :: wbundle
type(gsi_grid) :: grid
integer(i_kind) :: ii,jj,i,ic,id,ngases,nclouds,istatus
character(len=10),allocatable,dimension(:) :: gases
character(len=max_varname_length),allocatable,dimension(:) :: clouds
! Note: The following does not aim to get all variables in
! the state and control vectors, but rather the ones
! explicitly needed by this routine.
! Declare required local state variables
logical :: ls_u,ls_v,ls_prse,ls_q,ls_tsen,ls_tv,ls_ps,ls_ql,ls_qi
integer(i_kind), parameter :: nsvars = 9
integer(i_kind) :: isps(nsvars)
character(len=4), parameter :: mysvars(nsvars) = (/ & ! vars from SV needed here
'u ', 'v ', 'prse', 'q ', 'tsen', 'tv ', 'ps ','ql ', 'qi ' /)
character(len=max_varname_length),allocatable,dimension(:) :: cvars2dpm ! names of 2d fields including
! motley vars (if any)
real(r_kind),pointer,dimension(:,:) :: sv_ps=>NULL()
real(r_kind),pointer,dimension(:,:) :: sv_sst=>NULL()
real(r_kind),pointer,dimension(:,:,:) :: sv_u=>NULL()
real(r_kind),pointer,dimension(:,:,:) :: sv_v=>NULL()
real(r_kind),pointer,dimension(:,:,:) :: sv_prse=>NULL()
real(r_kind),pointer,dimension(:,:,:) :: sv_q=>NULL()
real(r_kind),pointer,dimension(:,:,:) :: sv_tsen=>NULL()
real(r_kind),pointer,dimension(:,:,:) :: sv_tv=>NULL()
real(r_kind),pointer,dimension(:,:,:) :: sv_oz=>NULL()
real(r_kind),pointer,dimension(:,:,:) :: sv_rank3=>NULL()
real(r_kind),pointer,dimension(:,:) :: sv_rank2=>NULL()
logical :: do_balance,do_getprs_tl,do_normal_rh_to_q,do_tv_to_tsen,do_getuv,cw_to_hydro
!******************************************************************************
if (.not.lsqrtb) then
write(6,*)trim(myname),': assumes lsqrtb'
call stop2(104)
end if
if (nsubwin/=1 .and. .not.l4dvar) then
write(6,*)trim(myname),': error 3dvar',nsubwin
call stop2(105)
end if
! create extended list of 2d variable names to include motley variables.
! NOTE: if mvars=0, there are no motley variables so cvars2dpm=cvars2d.
allocate(cvars2dpm(nc2d+mvars))
do i=1,nc2d
cvars2dpm(i)=cvars2d(i)
end do
do i=1,mvars
cvars2dpm(nc2d+i)=cvarsmd(i)
end do
call gsi_gridcreate(grid,lat2,lon2,nsig)
! Inquire about chemistry
call gsi_chemguess_get('dim',ngases,istatus)
if (ngases>0) then
allocate(gases(ngases))
call gsi_chemguess_get('gsinames',gases,istatus)
endif
! Since each internal vector of xhat has the same structure, pointers are
! the same independent of the subwindow jj
call gsi_bundlegetpointer (sval(1),mysvars,isps,istatus)
ls_u =isps(1)>0; ls_v =isps(2)>0; ls_prse=isps(3)>0
ls_q =isps(4)>0; ls_tsen=isps(5)>0; ls_tv =isps(6)>0
ls_ps =isps(7)>0; ls_ql =isps(8)>0; ls_qi =isps(9)>0
! Define what to do depending on what's in SV and
! what's explictly needed in this routine
do_balance =ls_tv.and.ls_ps
do_getprs_tl =ls_ps.and.ls_tv .and.ls_prse
do_normal_rh_to_q=ls_tv.and.ls_prse.and.ls_q
do_tv_to_tsen =ls_tv.and.ls_q .and.ls_tsen
do_getuv =ls_u .and.ls_v
! Inquire about cloud-variables
cw_to_hydro=.false.
call gsi_metguess_get ('clouds::3d',nclouds,istatus)
if (nclouds>0) then
allocate(clouds(nclouds))
call gsi_metguess_get ('clouds::3d',clouds,istatus)
if (getindex(cvars3d,'cw')>0 .and. ls_ql .and. ls_qi) cw_to_hydro=.true.
end if
! Loop over control steps
do jj=1,nsubwin
! create an internal structure w/ the same vars as those in the control vector, including motley vars
call gsi_bundlecreate (wbundle,grid,'control2model work',istatus,names2d=cvars2dpm,names3d=cvars3d)
if(istatus/=0) then
write(6,*) trim(myname), ': trouble creating work bundle'
call stop2(999)
endif
! Multiply by sqrt of background error (ckerror)
! -----------------------------------------------------------------------------
! Apply sqrt of variance, as well as vertical & horizontal parts of background
! error
call ckgcov(xhat%step(jj)%values(:),wbundle,size(xhat%step(jj)%values(:)))
! Get pointers to required state variables
call gsi_bundlegetpointer (sval(jj),'u' ,sv_u, istatus)
call gsi_bundlegetpointer (sval(jj),'v' ,sv_v, istatus)
call gsi_bundlegetpointer (sval(jj),'ps' ,sv_ps, istatus)
call gsi_bundlegetpointer (sval(jj),'prse',sv_prse,istatus)
call gsi_bundlegetpointer (sval(jj),'tv' ,sv_tv, istatus)
call gsi_bundlegetpointer (sval(jj),'tsen',sv_tsen,istatus)
call gsi_bundlegetpointer (sval(jj),'q' ,sv_q , istatus)
call gsi_bundlegetpointer (sval(jj),'oz' ,sv_oz , istatus)
call gsi_bundlegetpointer (sval(jj),'sst' ,sv_sst, istatus)
! Copy variables from CV to SV
call gsi_bundlegetvar ( wbundle, 'sf' , workst, istatus )
call gsi_bundlegetvar ( wbundle, 'vp' , workvp, istatus )
call gsi_bundlegetvar ( wbundle, 'q' , workrh, istatus )
call gsi_bundlegetvar ( wbundle, 't' , sv_tv, istatus )
call gsi_bundlegetvar ( wbundle, 'oz' , sv_oz, istatus )
call gsi_bundlegetvar ( wbundle, 'ps' , sv_ps, istatus )
call gsi_bundlegetvar ( wbundle, 'sst', sv_sst, istatus )
! Same one-to-one map for chemistry-vars; take care of them together
do ic=1,ngases
id=getindex(cvars3d,gases(ic))
if (id>0) then
call gsi_bundlegetpointer (sval(jj),gases(ic),sv_rank3,istatus)
call gsi_bundlegetvar (wbundle, gases(ic),sv_rank3,istatus)
endif
id=getindex(cvars2d,gases(ic))
if (id>0) then
call gsi_bundlegetpointer (sval(jj),gases(ic),sv_rank2,istatus)
call gsi_bundlegetvar (wbundle, gases(ic),sv_rank2,istatus)
endif
enddo
! Balance equation
if(do_balance) call balance(sv_tv,sv_ps,workst,workvp,fpsproj,fut2ps)
! -----------------------------------------------------------------------------
! Get 3d pressure
if(do_getprs_tl) call getprs_tl(sv_ps,sv_tv,sv_prse)
! Convert input normalized RH to q
if(do_normal_rh_to_q) call normal_rh_to_q(workrh,sv_tv,sv_prse,sv_q)
! Calculate sensible temperature
if(do_tv_to_tsen) call tv_to_tsen(sv_tv,sv_q,sv_tsen)
! Convert streamfunction and velocity potential to u,v
if(do_getuv) call getuv(sv_u,sv_v,workst,workvp,0)
if (nclouds>0) then
if (cw_to_hydro) then
! Case when cw is generated from hydrometeors
call cw2hydro_tl(sval(jj),wbundle,clouds,nclouds)
else
! Case when cloud-vars map one-to-one, take care of them together
do ic=1,nclouds
id=getindex(cvars3d,clouds(ic))
if (id>0) then
call gsi_bundlegetpointer (sval(jj),clouds(ic),sv_rank3,istatus)
call gsi_bundlegetvar (wbundle, clouds(ic),sv_rank3,istatus)
endif
enddo
end if
end if
! destroy temporary bundle
call gsi_bundledestroy(wbundle,istatus)
if(istatus/=0) then
write(6,*) trim(myname), ': trouble destroying work bundle'
call stop2(999)
endif
end do
! Bias correction terms
do ii=1,nsclen
bval%predr(ii)=xhat%predr(ii)*sqrt(varprd(ii))
enddo
do ii=1,npclen
bval%predp(ii)=xhat%predp(ii)*sqrt(varprd(nsclen+ii))
enddo
if (ntclen>0) then
do ii=1,ntclen
bval%predt(ii)=xhat%predt(ii)*sqrt(varprd(nsclen+npclen+ii))
enddo
end if
! Clean up
if (ngases>0) then
deallocate(gases)
endif
if (nclouds>0) deallocate(clouds)
return
end subroutine control2model