Skip to content

Commit

Permalink
Add FED EnVar DA Capability (NOAA-EMC#632)
Browse files Browse the repository at this point in the history
- This PR supports RRFS_B GSI FED assimilation. 
- This PR adds a new GSI EnVar FED assimilation capability.

The summary of the changes:
- Read FED background and ensemble from restart phy files
- Add new control/state variable of fed ( in anavinfo, section: metguess, state and control variable)
- Create intfed.f90 and sfpfed.f90 for minimization.
- Other related codes. For example, update hydrometers when either dbz or fed is assimilated, or both are assimilated. Previously the update of hydrometers is done only when dbz is assimilated.

Please see  Fixes NOAA-EMC#622  

This PR was tested with:
- One FED obs DA test
- Real FED DA with pseudo ensemble for code development and debug
- Real FED DA with real  ensemble 
---------
Co-authored-by: David Dowell <David.Dowell@noaa.gov>
hongli-wang authored Oct 31, 2023
1 parent 2cb0f5b commit acfe56d
Showing 13 changed files with 997 additions and 487 deletions.
235 changes: 179 additions & 56 deletions src/gsi/cplr_get_fv3_regional_ensperts.f90

Large diffs are not rendered by default.

10 changes: 10 additions & 0 deletions src/gsi/gsi_fedOper.F90
Original file line number Diff line number Diff line change
@@ -9,6 +9,7 @@ module gsi_fedOper
! 2023-07-10 D. Dowell - created new module for FED (flash extent
! density); gsi_dbzOper.F90 code used as a
! starting point for developing this new module
! 2023-08-24 H. Wang - Turned on intfed and stpfed
!
! input argument list: see Fortran 90 style document below
!
@@ -128,6 +129,7 @@ subroutine setup_(self, lunin, mype, is, nobs, init_pass, last_pass)
end subroutine setup_

subroutine intjo1_(self, ibin, rval,sval, qpred,sbias)
use intfedmod, only: intjo => intfed
use gsi_bundlemod , only: gsi_bundle
use bias_predictors, only: predictors
use m_obsNode , only: obsNode
@@ -145,9 +147,14 @@ subroutine intjo1_(self, ibin, rval,sval, qpred,sbias)
character(len=*),parameter:: myname_=myname//"::intjo1_"
class(obsNode),pointer:: headNode

headNode => obsLList_headNode(self%obsLL(ibin))
call intjo(headNode, rval,sval)
headNode => null()

end subroutine intjo1_

subroutine stpjo1_(self, ibin, dval,xval,pbcjo,sges,nstep,dbias,xbias)
use stpfedmod, only: stpjo => stpfed
use gsi_bundlemod, only: gsi_bundle
use bias_predictors, only: predictors
use m_obsNode , only: obsNode
@@ -169,6 +176,9 @@ subroutine stpjo1_(self, ibin, dval,xval,pbcjo,sges,nstep,dbias,xbias)
character(len=*),parameter:: myname_=myname//"::stpjo1_"
class(obsNode),pointer:: headNode

headNode => obsLList_headNode(self%obsLL(ibin))
call stpjo(headNode,dval,xval,pbcjo(:),sges,nstep)
headNode => null()
end subroutine stpjo1_

end module gsi_fedOper
2 changes: 2 additions & 0 deletions src/gsi/gsi_files.cmake
Original file line number Diff line number Diff line change
@@ -274,6 +274,7 @@ intaod.f90
intcldch.f90
intco.f90
intdbz.f90
intfed.f90
intdw.f90
intgps.f90
intgust.f90
@@ -594,6 +595,7 @@ stpcalc.f90
stpcldch.f90
stpco.f90
stpdbz.f90
stpfed.f90
stpdw.f90
stpgps.f90
stpgust.f90
89 changes: 55 additions & 34 deletions src/gsi/gsi_rfv3io_mod.f90

Large diffs are not rendered by default.

39 changes: 32 additions & 7 deletions src/gsi/gsimod.F90
Original file line number Diff line number Diff line change
@@ -23,12 +23,13 @@ module gsimod
use gsi_dbzOper, only: diag_radardbz
use gsi_fedOper, only: diag_fed

use obsmod, only: doradaroneob,oneoblat,oneoblon,oneobheight,oneobvalue,oneobddiff,oneobradid,&
use obsmod, only: doradaroneob,dofedoneob,oneoblat,oneoblon,oneobheight,oneobvalue,oneobddiff,oneobradid,&
radar_no_thinning,ens_hx_dbz_cut,static_gsi_nopcp_dbz,rmesh_dbz,&
rmesh_vr,zmesh_dbz,zmesh_vr,if_vterminal, if_model_dbz,if_vrobs_raw,if_use_w_vr,&
rmesh_vr,zmesh_dbz,zmesh_vr,if_vterminal, if_model_dbz,if_model_fed,innov_use_model_fed,if_vrobs_raw,if_use_w_vr,&
minobrangedbz,maxobrangedbz,maxobrangevr,maxtiltvr,missing_to_nopcp,&
ntilt_radarfiles,whichradar,&
minobrangevr,maxtiltdbz,mintiltvr,mintiltdbz,l2rwthin,hurricane_radar
minobrangevr,maxtiltdbz,mintiltvr,mintiltdbz,l2rwthin,hurricane_radar,&
r_hgt_fed

use obsmod, only: lwrite_predterms, &
lwrite_peakwt,use_limit,lrun_subdirs,l_foreaft_thin,lobsdiag_forenkf,&
@@ -202,7 +203,7 @@ module gsimod
use gsi_nstcouplermod, only: gsi_nstcoupler_init_nml
use gsi_nstcouplermod, only: nst_gsi,nstinfo,zsea1,zsea2,fac_dtl,fac_tsl
use ncepnems_io, only: init_nems,imp_physics,lupp
use wrf_vars_mod, only: init_wrf_vars
use wrf_vars_mod, only: init_wrf_vars,fed_exist,dbz_exist
use gsi_rfv3io_mod,only : fv3sar_bg_opt
use radarz_cst, only: mphyopt, MFflg
use radarz_iface, only: init_mphyopt
@@ -510,6 +511,15 @@ module gsimod
! 2023-07-30 Zhao - added namelist options for analysis of significant wave height
! (aka howv in GSI code): corp_howv, hwllp_howv
! (in namelist session rapidrefresh_cldsurf)
!
! 2023-09-14 H. Wang - add namelist option for FED EnVar DA.
! - if_model_fed=.true. : FED in background and ens. If
! perform FED DA, this has to be true along with fed in
! control/analysis and metguess vectors. If only run GSI observer,
! it can be false.
! - innov_use_model_fed=.true. : Use FED from BG to calculate innovation.
! this requires if_model_fed=.true.
! it works either an EnVar DA run or a GSI observer run.
!
!EOP
!-------------------------------------------------------------------------
@@ -770,16 +780,17 @@ module gsimod
use_sp_eqspace,lnested_loops,lsingleradob,thin4d,use_readin_anl_sfcmask,&
luse_obsdiag,id_drifter,id_ship,verbose,print_obs_para,lsingleradar,singleradar,lnobalance, &
missing_to_nopcp,minobrangedbz,minobrangedbz,maxobrangedbz,&
maxobrangevr,maxtiltvr,whichradar,doradaroneob,oneoblat,&
maxobrangevr,maxtiltvr,whichradar,doradaroneob,dofedoneob,oneoblat,&
oneoblon,oneobheight,oneobvalue,oneobddiff,oneobradid,&
rmesh_vr,zmesh_dbz,zmesh_vr, ntilt_radarfiles, whichradar,&
radar_no_thinning,ens_hx_dbz_cut,static_gsi_nopcp_dbz,rmesh_dbz,&
minobrangevr, maxtiltdbz, mintiltvr,mintiltdbz,if_vterminal,if_vrobs_raw,if_use_w_vr,&
if_model_dbz,imp_physics,lupp,netcdf_diag,binary_diag,l_wcp_cwm,aircraft_recon,diag_version,&
if_model_dbz,if_model_fed,innov_use_model_fed,imp_physics,lupp,netcdf_diag,binary_diag,l_wcp_cwm,aircraft_recon,diag_version,&
write_fv3_incr,incvars_to_zero,incvars_zero_strat,incvars_efold,diag_version,&
cao_check,lcalc_gfdl_cfrac,tau_fcst,efsoi_order,lupdqc,lqcoef,cnvw_option,l2rwthin,hurricane_radar,&
l_reg_update_hydro_delz, l_obsprvdiag,&
l_use_dbz_directDA, l_use_rw_columntilt, ta2tb, optconv
l_use_dbz_directDA, l_use_rw_columntilt, ta2tb, optconv, &
r_hgt_fed

! GRIDOPTS (grid setup variables,including regional specific variables):
! jcap - spectral resolution
@@ -1978,6 +1989,20 @@ subroutine gsimain_initialize
endif
endif

if (innov_use_model_fed .and. .not.if_model_fed) then
if(mype==0) write(6,*)' GSIMOD: invalid innov_use_model_fed=.true. but if_model_fed=.false.'
call die(myname_,'invalid innov_use_model_fed,if_model_fed, check namelist settings',330)
end if

if (.not. (miter == 0 .or. lobserver) .and. if_model_fed .and. .not. fed_exist) then
if(mype==0) write(6,*)' GSIMOD: .not. (miter == 0 .or. lobserver) and if_model_fed=.true. but fed is not in anavinfo file'
call die(myname_,'Please check namelist parameters and/or add fed in anavinfo (contro/state_vector and met_guess) when miter > 0 and if_model_fed=.true.',332)
end if

if (.not. (miter == 0 .or. lobserver) .and. if_model_dbz .and. .not. dbz_exist) then
if(mype==0) write(6,*)' GSIMOD: .not. (miter == 0 .or. lobserver) and if_model_dbz=.true. but dbz is not in anavinfo file'
call die(myname_,'Please check namelist parameters and/or add dbz in anavinfo (contro/state_vector and met_guess) when miter > 0 and if_model_fed=.true.',334)
end if

! Ensure valid number of horizontal scales
if (nhscrf<0 .or. nhscrf>3) then
187 changes: 187 additions & 0 deletions src/gsi/intfed.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,187 @@
module intfedmod
!$$$ module documentation block
! . . . .
! module: intfedmod module for intfed and its tangent linear intfed_tl
! prgmmr:
!
! abstract: module for intfed and its tangent linear intfed_tl
!
! program history log:
! 2023-08-24 H. Wang - add tangent linear of fed operator to directly assimilate FED
!
! subroutines included:
! sub intfed_
!
! variable definitions:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block

use m_obsNode, only: obsNode
use m_fedNode, only: fedNode
use m_fedNode, only: fedNode_typecast
use m_fedNode, only: fedNode_nextcast
use m_obsdiagNode, only: obsdiagNode_set
implicit none

PRIVATE
PUBLIC intfed

interface intfed; module procedure &
intfed_
end interface

contains

subroutine intfed_(fedhead,rval,sval)
!$$$ subprogram documentation block
! . . . .
! subprogram: intfed apply nonlin qc operator for GLM FED
!
! abstract: apply observation operator for radar winds
! with nonlinear qc operator
!
! program history log:
! 2023-08-24 H.Wang - modified based on intdbz.f90
! - using tangent linear fed operator

!
! input argument list:
! fedhead - obs type pointer to obs structure
! sfed - current fed solution increment
!
! output argument list:
! rfed - fed results from fed observation operator
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use constants, only: half,one,tiny_r_kind,cg_term,r3600
use obsmod, only: lsaveobsens,l_do_adjoint,luse_obsdiag
use qcmod, only: nlnqc_iter,varqc_iter
use jfunc, only: jiter
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
use gsi_4dvar, only: ladtest_obs
use wrf_vars_mod, only : fed_exist
implicit none

! Declare passed variables
class(obsNode), pointer, intent(in ) :: fedhead
type(gsi_bundle), intent(in ) :: sval
type(gsi_bundle), intent(inout) :: rval

! Declare local variables
integer(i_kind) j1,j2,j3,j4,j5,j6,j7,j8,ier,istatus
! real(r_kind) penalty
real(r_kind) val,w1,w2,w3,w4,w5,w6,w7,w8,valfed
real(r_kind) cg_fed,p0,grad,wnotgross,wgross,pg_fed
real(r_kind),pointer,dimension(:) :: sfed
real(r_kind),pointer,dimension(:) :: rfed
type(fedNode), pointer :: fedptr

! If no fed obs type data return
if(.not. associated(fedhead))return

! Retrieve pointers
! Simply return if any pointer not found
ier=0
if(fed_exist)then
call gsi_bundlegetpointer(sval,'fed',sfed,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'fed',rfed,istatus);ier=istatus+ier
else
return
end if

if(ier/=0)return


fedptr => fedNode_typecast(fedhead)
do while (associated(fedptr))
j1=fedptr%ij(1)
j2=fedptr%ij(2)
j3=fedptr%ij(3)
j4=fedptr%ij(4)
j5=fedptr%ij(5)
j6=fedptr%ij(6)
j7=fedptr%ij(7)
j8=fedptr%ij(8)
w1=fedptr%wij(1)
w2=fedptr%wij(2)
w3=fedptr%wij(3)
w4=fedptr%wij(4)
w5=fedptr%wij(5)
w6=fedptr%wij(6)
w7=fedptr%wij(7)
w8=fedptr%wij(8)


! Forward model
if( fed_exist )then
val = w1* sfed(j1)+w2* sfed(j2)+w3* sfed(j3)+w4* sfed(j4)+ &
w5* sfed(j5)+w6* sfed(j6)+w7* sfed(j7)+w8* sfed(j8)
end if

if(luse_obsdiag)then
if (lsaveobsens) then
grad = val*fedptr%raterr2*fedptr%err2
!-- fedptr%diags%obssen(jiter) = grad
call obsdiagNode_set(fedptr%diags,jiter=jiter,obssen=grad)

else
!-- if (fedptr%luse) fedptr%diags%tldepart(jiter)=val
if (fedptr%luse) call obsdiagNode_set(fedptr%diags,jiter=jiter,tldepart=val)
endif
endif

if (l_do_adjoint) then
if (.not. lsaveobsens) then
if( .not. ladtest_obs ) val=val-fedptr%res

! gradient of nonlinear operator
if (nlnqc_iter .and. fedptr%pg > tiny_r_kind .and. &
fedptr%b > tiny_r_kind) then
pg_fed=fedptr%pg*varqc_iter
cg_fed=cg_term/fedptr%b
wnotgross= one-pg_fed
wgross = pg_fed*cg_fed/wnotgross
p0 = wgross/(wgross+exp(-half*fedptr%err2*val**2))
val = val*(one-p0)
endif

if( ladtest_obs) then
grad = val
else
grad = val*fedptr%raterr2*fedptr%err2
end if

endif

! Adjoint
if(fed_exist)then
valfed = grad
rfed(j1)=rfed(j1)+w1*valfed
rfed(j2)=rfed(j2)+w2*valfed
rfed(j3)=rfed(j3)+w3*valfed
rfed(j4)=rfed(j4)+w4*valfed
rfed(j5)=rfed(j5)+w5*valfed
rfed(j6)=rfed(j6)+w6*valfed
rfed(j7)=rfed(j7)+w7*valfed
rfed(j8)=rfed(j8)+w8*valfed
end if

endif

!fedptr => fedptr%llpoint
fedptr => fedNode_nextcast(fedptr)
end do
return
end subroutine intfed_

end module intfedmod
13 changes: 12 additions & 1 deletion src/gsi/m_berror_stats_reg.f90
Original file line number Diff line number Diff line change
@@ -400,7 +400,7 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt
integer(i_kind) :: nrf2_td2m,nrf2_mxtm,nrf2_mitm,nrf2_pmsl,nrf2_howv,nrf2_tcamt,nrf2_lcbas,nrf2_cldch
integer(i_kind) :: nrf2_uwnd10m,nrf2_vwnd10m
integer(i_kind) :: nrf3_sfwter,nrf3_vpwter
integer(i_kind) :: nrf3_dbz
integer(i_kind) :: nrf3_dbz,nrf3_fed
integer(i_kind) :: nrf3_ql,nrf3_qi,nrf3_qr,nrf3_qs,nrf3_qg,nrf3_qnr,nrf3_w
integer(i_kind) :: inerr,istat
integer(i_kind) :: nsigstat,nlatstat,isig
@@ -624,6 +624,7 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt
nrf3_sf =getindex(cvars3d,'sf')
nrf3_vp =getindex(cvars3d,'vp')
nrf3_dbz=getindex(cvars3d,'dbz')
nrf3_fed=getindex(cvars3d,'fed')
nrf2_sst=getindex(cvars2d,'sst')
nrf2_gust=getindex(cvars2d,'gust')
nrf2_vis=getindex(cvars2d,'vis')
@@ -671,6 +672,16 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt
vz(:,:,nrf3_dbz)=vz(:,:,nrf3_t)
endif

if( nrf3_fed>0 )then
if(.not. nrf3_t>0) then
write(6,*)'not as expect,stop'
stop
endif
corz(:,:,nrf3_fed)=10.0_r_kind
hwll(:,:,nrf3_fed)=hwll(:,:,nrf3_t)
vz(:,:,nrf3_fed)=vz(:,:,nrf3_t)
endif

if (nrf3_oz>0) then
factoz = 0.0002_r_kind*r25
corz(:,:,nrf3_oz)=factoz
22 changes: 19 additions & 3 deletions src/gsi/obsmod.F90
Original file line number Diff line number Diff line change
@@ -161,6 +161,7 @@ module obsmod
! observation provider and sub-provider information into
! obsdiags files (used for AutoObsQC)
! 2023-07-10 Y. Wang, D. Dowell - add variables for flash extent density
! 2023-10-10 H. Wang (GSL) - add variables for flash extent density EnVar DA
!
! Subroutines Included:
! sub init_obsmod_dflts - initialize obs related variables to default values
@@ -188,6 +189,12 @@ module obsmod
! def diag_radardbz- namelist logical to compute/write (=true) radar
! reflectiivty diag files
! def diag_fed - namelist logical to compute/write (=true) flash extent density diag files
! def innov_use_model_fed - namelist logical. True: use (the FEB in background to calculate innovation
! False: calculate innvation use
! the obs operator in GSI
! def if_model_fed - namelist logical. True: Read in FED from background
! including from ensemble.
! def r_hgt_fed - height of fed observations
! def reduce_diag - namelist logical to produce reduced radiance diagnostic files
! def use_limit - parameter set equal to -1 if diag files produced or 0 if not diag files or reduce_diag
! def obs_setup - prefix for files passing pe relative obs data to setup routines
@@ -487,7 +494,12 @@ module obsmod
public :: iout_dbz, mype_dbz
! --- DBZ DA ---

! ==== FED DA ===
public :: if_model_fed, innov_use_model_fed
public :: r_hgt_fed
public :: iout_fed, mype_fed
public :: dofedoneob
! --- FED DA ---

public :: obsmod_init_instr_table
public :: obsmod_final_instr_table
@@ -577,7 +589,7 @@ module obsmod

real(r_kind) perturb_fact,time_window_max,time_offset,time_window_rad
real(r_kind),dimension(50):: dmesh

real(r_kind) r_hgt_fed
integer(i_kind) nchan_total,ianldate
integer(i_kind) ndat,ndat_types,ndat_times,nprof_gps
integer(i_kind) lunobs_obs,nloz_v6,nloz_v8,nobskeep,nloz_omi
@@ -621,8 +633,8 @@ module obsmod
integer(i_kind) ntilt_radarfiles,tcp_posmatch,tcp_box

logical :: ta2tb
logical :: doradaroneob
logical :: vr_dealisingopt, if_vterminal, if_model_dbz, inflate_obserr, if_vrobs_raw, if_use_w_vr, l2rwthin
logical :: doradaroneob,dofedoneob
logical :: vr_dealisingopt, if_vterminal, if_model_dbz,if_model_fed, innov_use_model_fed,inflate_obserr, if_vrobs_raw, if_use_w_vr, l2rwthin
character(4) :: whichradar,oneobradid
real(r_kind) :: oneoblat,oneoblon,oneobddiff,oneobvalue,oneobheight
logical :: radar_no_thinning
@@ -755,11 +767,15 @@ subroutine init_obsmod_dflts
if_vrobs_raw=.false.
if_use_w_vr=.true.
if_model_dbz=.false.
if_model_fed=.false.
innov_use_model_fed=.false.
inflate_obserr=.false.
whichradar="KKKK"

oneobradid="KKKK"
doradaroneob=.false.
r_hgt_fed=6500_r_kind
dofedoneob=.false.
oneoblat=-999_r_kind
oneoblon=-999_r_kind
oneobddiff=-999_r_kind
3 changes: 2 additions & 1 deletion src/gsi/read_dbz_netcdf.f90
Original file line number Diff line number Diff line change
@@ -526,7 +526,8 @@ subroutine read_dbz_mrms_netcdf(nread,ndata,nodata,infile,obstype,lunout,sis,nob

if(thislon>=r360) thislon=thislon-r360
if(thislon<zero ) thislon=thislon+r360

if(thislon>=r360 .or. thislat >90.0_r_kind) cycle

!-Convert back to radians

thislat = thislat*deg2rad
445 changes: 162 additions & 283 deletions src/gsi/read_fed.f90

Large diffs are not rendered by default.

229 changes: 129 additions & 100 deletions src/gsi/setupfed.f90

Large diffs are not rendered by default.

171 changes: 171 additions & 0 deletions src/gsi/stpfed.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,171 @@
module stpfedmod

!$$$ module documentation block
! . . . .
! module: stpfedmod module for stpfed and its tangent linear stpfed_tl
! prgmmr:
!
! abstract: module for stpfed and its tangent linear stpfed_tl
!
! program history log:
! 2023-08-23 H. Wang - Modified based on sftdbzmod
! - add adjoint of fed operator
!
! subroutines included:
! sub stpfed
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block

implicit none

PRIVATE
PUBLIC stpfed

contains

subroutine stpfed(fedhead,rval,sval,out,sges,nstep)
!$$$ subprogram documentation block
! . . . .
! subprogram: stpfed calculate penalty and contribution to
! stepsize with nonlinear qc added.
! prgmmr: derber org: np23 date: 1991-02-26
!
!
! program history log:
! 2019-08-23 H.Wang - added for FED DA
! input argument list:
! fedhead
! sges - step size estimates (nstep)
! nstep - number of step sizes (== 0 means use outer iteration value)
!
! output argument list - output for step size calculation
! out(1:nstep) - penalty from fed sges(1:nstep)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind,r_quad
use qcmod, only: nlnqc_iter,varqc_iter
use constants, only: half,one,two,tiny_r_kind,cg_term,zero_quad,r3600
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
use gridmod, only: wrf_mass_regional, fv3_regional
use wrf_vars_mod, only : fed_exist
use m_obsNode, only: obsNode
use m_fedNode , only: fedNode
use m_fedNode , only: fedNode_typecast
use m_fedNode , only: fedNode_nextcast
! use directDA_radaruse_mod, only: l_use_fed_directDA
use radarz_cst, only: mphyopt

implicit none

! Declare passed variables
class(obsNode), pointer ,intent(in ) :: fedhead
integer(i_kind) ,intent(in ) :: nstep
real(r_quad),dimension(max(1,nstep)),intent(inout) :: out
type(gsi_bundle) ,intent(in ) :: rval,sval
real(r_kind),dimension(max(1,nstep)),intent(in ) :: sges

! Declare local variables
integer(i_kind) ier,istatus
integer(i_kind) j1,j2,j3,j4,j5,j6,j7,j8,kk
real(r_kind) w1,w2,w3,w4,w5,w6,w7,w8
real(r_kind) valfed
real(r_kind) fedcur
real(r_kind) cg_fed,fed,wgross,wnotgross
real(r_kind),dimension(max(1,nstep))::pen
real(r_kind) pg_fed
real(r_kind),pointer,dimension(:) :: sfed
real(r_kind),pointer,dimension(:) :: rfed
type(fedNode), pointer :: fedptr

out=zero_quad

! If no fed data return
if(.not. associated(fedhead))return

! Retrieve pointers
! Simply return if any pointer not found
ier=0
if(fed_exist)then
call gsi_bundlegetpointer(sval,'fed',sfed,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'fed',rfed,istatus);ier=istatus+ier
else
return
end if

if(ier/=0)return

fedptr => fedNode_typecast(fedhead)
do while (associated(fedptr))
if(fedptr%luse)then
if(nstep > 0)then
j1=fedptr%ij(1)
j2=fedptr%ij(2)
j3=fedptr%ij(3)
j4=fedptr%ij(4)
j5=fedptr%ij(5)
j6=fedptr%ij(6)
j7=fedptr%ij(7)
j8=fedptr%ij(8)
w1=fedptr%wij(1)
w2=fedptr%wij(2)
w3=fedptr%wij(3)
w4=fedptr%wij(4)
w5=fedptr%wij(5)
w6=fedptr%wij(6)
w7=fedptr%wij(7)
w8=fedptr%wij(8)

if( fed_exist )then
valfed= w1* rfed(j1)+w2*rfed(j2)+w3*rfed(j3)+w4*rfed(j4)+ &
w5* rfed(j5)+w6*rfed(j6)+w7*rfed(j7)+w8*rfed(j8)

fedcur= w1* sfed(j1)+w2* sfed(j2)+w3* sfed(j3)+w4*sfed(j4)+ &
w5* sfed(j5)+w6* sfed(j6)+w7* sfed(j7)+w8* sfed(j8)- &
fedptr%res
end if


do kk=1,nstep
fed=fedcur+sges(kk)*valfed
pen(kk)=fed*fed*fedptr%err2
end do
else
pen(1)=fedptr%res*fedptr%res*fedptr%err2
end if

! Modify penalty term if nonlinear QC
if (nlnqc_iter .and. fedptr%pg > tiny_r_kind .and. &
fedptr%b > tiny_r_kind) then

pg_fed=fedptr%pg*varqc_iter
cg_fed=cg_term/fedptr%b
wnotgross= one-pg_fed
wgross = pg_fed*cg_fed/wnotgross
do kk=1,max(1,nstep)
pen(kk)= -two*log((exp(-half*pen(kk)) + wgross)/(one+wgross))
end do
end if

out(1) = out(1)+pen(1)*fedptr%raterr2
kk=1
do kk=2,nstep
out(kk) = out(kk)+(pen(kk)-pen(1))*fedptr%raterr2
end do
end if

fedptr => fedNode_nextcast(fedptr)

end do
return
end subroutine stpfed

end module stpfedmod
39 changes: 37 additions & 2 deletions src/gsi/wrf_vars_mod.f90
Original file line number Diff line number Diff line change
@@ -39,28 +39,63 @@ module wrf_vars_mod
use mpimod, only: mype
use control_vectors, only: nc3d,cvars3d
use kinds, only: i_kind
use gsi_metguess_mod, only: gsi_metguess_get
use constants, only: max_varname_length
implicit none
private
! public methods
public :: init_wrf_vars
! common block variables
public :: w_exist
public :: dbz_exist
public :: fed_exist

logical,save :: w_exist, dbz_exist
logical,save :: w_exist, dbz_exist, fed_exist
contains

subroutine init_wrf_vars
integer(i_kind) ii
integer(i_kind) ii,istatus
character(max_varname_length),allocatable,dimension(:) :: cloud
integer(i_kind) ncloud
logical :: dbz_cloud_exist,fed_cloud_exist

w_exist=.false.
dbz_exist=.false.
fed_exist=.false.
dbz_cloud_exist=.false.
fed_cloud_exist=.false.

do ii=1,nc3d
if(mype == 0 ) write(6,*)"anacv cvars3d is ",cvars3d(ii)
if(trim(cvars3d(ii)) == 'w'.or.trim(cvars3d(ii))=='W') w_exist=.true.
if(trim(cvars3d(ii))=='dbz'.or.trim(cvars3d(ii))=='DBZ') dbz_exist=.true.
if(trim(cvars3d(ii))=='fed'.or.trim(cvars3d(ii))=='FED') fed_exist=.true.
enddo

! Inquire about clouds

call gsi_metguess_get('clouds::3d',ncloud,istatus)
if (ncloud>0) then
allocate(cloud(ncloud))
call gsi_metguess_get('clouds::3d',cloud,istatus)
endif

do ii=1,ncloud
if(mype == 0 ) write(6,*)"metguess cloud3d is ",cloud(ii)
if(trim(cloud(ii))=='fed'.or.trim(cloud(ii))=='FED')fed_cloud_exist=.true.
if(trim(cloud(ii))=='dbz'.or.trim(cloud(ii))=='DBZ')dbz_cloud_exist=.true.
end do

if(.not.fed_exist .or. .not.fed_cloud_exist )then
fed_exist=.false.
endif

if(.not.dbz_exist .or. .not.dbz_cloud_exist )then
dbz_exist=.false.
endif

if(ncloud>0) deallocate(cloud)

end subroutine init_wrf_vars

end module wrf_vars_mod

0 comments on commit acfe56d

Please sign in to comment.