module intallmod
!$$$ mdoule documentation block
! . . . .
! module: intallmod module for intall and its tangent linear intall_tl
! prgmmr:
!
! abstract: module for intall and its tangent linear intall_tl
!
! program history log:
! 2005-05-18 Yanqiu zhu - wrap intall and its tangent linear intall_tl into one module
! 2005-11-21 Derber - remove interface and clean up code
! 2008-11-26 Todling - remove intall_tl
! 2009-08-13 lueken - update documentation
! 2012-02-08 kleist - changes related to 4d-ensemble-var additions and consolidation of
! int... individual modules to one intjcmod
! 2015-09-03 guo - obsmod::yobs has been replaced with m_obsHeadBundle,
! where yobs is created and destroyed when and where it
! is needed.
!
! subroutines included:
! sub intall
!
! variable definitions:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
PRIVATE
PUBLIC intall
contains
subroutine intall(sval,sbias,rval,rbias)
!$$$ subprogram documentation block
! . . . .
! subprogram: intall calculate RHS for analysis equation
! prgmmr: derber org: np23 date: 2003-12-18
!
! abstract: calculate RHS for all variables (nonlinear qc version)
!
! A description of nonlinear qc follows:
!
! The observation penalty Jo is defined as
!
! Jo = - (sum over obs) 2*log(Po)
!
! where,
!
! Po = Wnotgross*exp(-.5*(Hn(x+xb) - yo)**2 ) + Wgross
! with
! Hn = the forward model (possibly non-linear) normalized by
! observation error
! x = the current estimate of the analysis increment
! xb = the background state
! yo = the observation normalized by observation error
!
! Note: The factor 2 in definition of Jo is present because the
! penalty Jo as used in this code is 2*(usual definition
! of penalty)
!
! Wgross = Pgross*cg
!
! Wnotgross = 1 - Wgross
!
! Pgross = probability of gross error for observation (assumed
! here to have uniform distribution over the possible
! range of values)
!
! cg = sqrt(2*pi)/2b
!
! b = possible range of variable for gross errors, normalized by
! observation error
!
! The values for the above parameters that Bill Collins used in the
! eta 3dvar are:
!
! cg = cg_term/b, where cg_term = sqrt(2*pi)/2
!
! b = 10. ! range for gross errors, normalized by obs error
!
! pg_q=.002 ! probability of gross error for specific humidity
! pg_pw=.002 ! probability of gross error for precipitable water
! pg_p=.002 ! probability of gross error for pressure
! pg_w=.005 ! probability of gross error for wind
! pg_t=.007 ! probability of gross error for temperature
! pg_rad=.002 ! probability of gross error for radiances
!
!
! Given the above Jo, the gradient of Jo is as follows:
!
! T
! gradx(Jo) = - (sum over observations) 2*H (Hn(x+xb)-yo)*(Po - Wgross)/Po
!
! where,
!
! H = tangent linear model of Hn about x+xb
!
!
! Note that if Pgross = 0.0, then Wnotgross=1.0 and Wgross=0.0. That is,
! the code runs as though nonlinear quality control were not present
! (which is indeed the case since the gross error probability is 0).
!
! As a result the same int* routines may be used for use with or without
! nonlinear quality control.
!
!
! program history log:
! 2003-12-18 derber
! 2004-07-23 derber - modify to include conventional sst
! 2004-07-28 treadon - add only to module use, add intent in/out
! 2004-10-06 parrish - add nonlinear qc option
! 2004-10-06 kleist - separate control vector for u,v, & convert int
! for wind components into int for st,vp
! 2004-11-30 treadon - add brightness temperatures to nonlinear
! quality control
! 2004-12-03 treadon - replace mpe_iallreduce (IBM extension) with
! standard mpi_allreduce
! 2005-01-20 okamoto - add u,v to intrad
! 2005-02-23 wu - changes related to normalized rh option
! 2005-04-11 treadon - rename intall_qc as intall
! 2005-05-18 yanqiu zhu - add 'use int*mod',and modify call interfaces for using these modules
! 2005-05-24 pondeca - take into consideration that npred=npredp=0
! for 2dvar only surface analysis option
! 2005-06-03 parrish - add horizontal derivatives
! 2005-07-10 kleist - add dynamic constraint term
! 2005-09-29 kleist - expand Jc term, include time derivatives vector
! 2005-11-21 kleist - separate tendencies from Jc term, add call to calctends adjoint
! 2005-12-01 cucurull - add code for GPS local bending angle, add use obsmod for ref_obs
! 2005-12-20 parrish - add arguments to call to intt to allow for option of using boundary
! layer forward tlm.
! 2006-02-03 derber - modify to increase reproducibility
! 2006-03-17 park - correct error in call to intt--rval,sval --> rvaluv,svaluv
! in order to correctly pass wind variables.
! 2006-04-06 kleist - include both Jc formulations
! 2006-07-26 parrish - correct inconsistency in computation of space and time derivatives of q
! currently, if derivatives computed, for q it is normalized q, but
! should be mixing ratio.
! 2006-07-26 parrish - add strong constraint initialization option
! 2007-03-19 tremolet - binning of observations
! 2007-04-13 tremolet - split Jo and 3dvar components into intjo and int3dvar
! 2007-10-01 todling - add timers
! 2011-10-20 todling - observation operators refer to state- not control-vec (cvars->svars)
! 2014-03-19 pondeca - Add RHS calculation for wspd10m constraint
! 2014-05-07 pondeca - Add RHS calculation for howv constraint
! 2014-06-17 carley/zhu - Add RHS calculation for lcbas constraint
! 2015-07-10 pondeca - Add RHS calculation for cldch constraint
!
! input argument list:
! sval - solution on grid
! sbias
! rval
! rbias
!
! output argument list:
! rval - RHS on grid
! rbias
!
! remarks:
! if strong initialization, then svalt, svalp, svaluv, sval_x, sval_y, sval_t
! are all grid fields after strong initialization.
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: i_kind,r_quad
use gsi_4dvar, only: nobs_bins,ltlint,ibin_anl
use constants, only: zero,zero_quad
use jcmod, only: ljcpdry,ljc4tlevs,ljcdfi
use jfunc, only: nrclen,nsclen,npclen,ntclen
use intradmod, only: setrad
use intjomod, only: intjo
use bias_predictors, only : predictors,assignment(=)
use state_vectors, only: allocate_state,deallocate_state
use intjcmod, only: intlimq,intlimg,intlimv,intlimp,intlimw10m,intlimhowv,intlimcldch,&
intliml,intjcpdry1,intjcpdry2,intjcdfi
use timermod, only: timer_ini,timer_fnl
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: assignment(=)
use state_vectors, only: svars2d
use mpeu_util, only: getindex
use guess_grids, only: ntguessig,nfldsig
use mpl_allreducemod, only: mpl_allreduce
use m_obsHeadBundle, only: obsHeadBundle
use m_obsHeadBundle, only: obsHeadBundle_create
use m_obsHeadBundle, only: obsHeadBundle_destroy
implicit none
! Declare passed variables
type(gsi_bundle), intent(in ) :: sval(nobs_bins)
type(predictors), intent(in ) :: sbias
type(gsi_bundle), intent(inout) :: rval(nobs_bins)
type(predictors), intent(inout) :: rbias
real(r_quad),dimension(max(1,nrclen),nobs_bins) :: qpred_bin
real(r_quad),dimension(max(1,nrclen)) :: qpred
real(r_quad),dimension(2*nobs_bins) :: mass
! Declare local variables
integer(i_kind) :: ibin,ii,it,i
type(obsHeadBundle),pointer,dimension(:):: yobs
!******************************************************************************
! Initialize timer
call timer_ini('intall')
! Zero gradient arrays
do ii=1,nobs_bins
rval(ii)=zero
enddo
! Compute RHS in physical space
call setrad(sval(1))
qpred_bin=zero_quad
call obsHeadBundle_create(yobs,nobs_bins)
! RHS for Jo
!$omp parallel do schedule(dynamic,1) private(ibin)
do ibin=1,size(yobs) ! == nobs_bins
call intjo(yobs(ibin),rval(ibin),qpred_bin(:,ibin),sval(ibin),sbias,ibin)
end do
qpred=zero_quad
do ibin=1,size(yobs) ! == nobs_bins
do i=1,nrclen
qpred(i)=qpred(i)+qpred_bin(i,ibin)
end do
end do
call obsHeadBundle_destroy(yobs)
if(.not.ltlint)then
! RHS for moisture constraint
if (.not.ljc4tlevs) then
call intlimq(rval(ibin_anl),sval(ibin_anl),ntguessig)
else
do ibin=1,nobs_bins
if (nobs_bins /= nfldsig) then
it=ntguessig
else
it=ibin
end if
call intlimq(rval(ibin),sval(ibin),it)
end do
end if
! RHS for gust constraint
if (getindex(svars2d,'gust')>0)call intlimg(rval(1),sval(1))
! RHS for vis constraint
if (getindex(svars2d,'vis')>0) call intlimv(rval(1),sval(1))
! RHS for pblh constraint
if (getindex(svars2d,'pblh')>0) call intlimp(rval(1),sval(1))
! RHS for wspd10m constraint
if (getindex(svars2d,'wspd10m')>0) call intlimw10m(rval(1),sval(1))
! RHS for howv constraint
if (getindex(svars2d,'howv')>0) call intlimhowv(rval(1),sval(1))
! RHS for lcbas constraint
if (getindex(svars2d,'lcbas')>0) call intliml(rval(1),sval(1))
! RHS for cldch constraint
if (getindex(svars2d,'cldch')>0) call intlimcldch(rval(1),sval(1))
end if
! RHS for dry ps constraint: part 1
if(ljcpdry)then
if (.not.ljc4tlevs) then
call intjcpdry1(sval(ibin_anl),1,mass)
else
call intjcpdry1(sval,nobs_bins,mass)
end if
! Put reduces together to minimize wait time
! First, use MPI to get global mean increment
call mpl_allreduce(2*nobs_bins,qpvals=mass)
end if
! Take care of background error for bias correction terms
call mpl_allreduce(nrclen,qpvals=qpred)
! RHS for dry ps constraint: part 2
if(ljcpdry)then
if (.not.ljc4tlevs) then
call intjcpdry2(rval(ibin_anl),1,mass)
else
call intjcpdry2(rval,nobs_bins,mass)
end if
end if
! RHS for Jc DFI
if (ljcdfi .and. nobs_bins>1) call intjcdfi(rval,sval)
if(nsclen > 0)then
do i=1,nsclen
rbias%predr(i)=qpred(i)
end do
end if
if(npclen > 0)then
do i=1,npclen
rbias%predp(i)=qpred(nsclen+i)
end do
end if
if (ntclen>0) then
do i=1,ntclen
rbias%predt(i)=qpred(nsclen+npclen+i)
end do
end if
! Finalize timer
call timer_fnl('intall')
return
end subroutine intall
end module intallmod