Skip to content

Commit

Permalink
GitHub Issue NOAA-EMC#219 Improve Minimization and fix bug in vqc
Browse files Browse the repository at this point in the history
  • Loading branch information
@jderber-NOAA
jderber-NOAA committed Oct 22, 2021
1 parent 059e402 commit f938842
Show file tree
Hide file tree
Showing 23 changed files with 268 additions and 372 deletions.
2 changes: 1 addition & 1 deletion regression/regression_var.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -131,7 +131,7 @@ case $machine in

# On Hera, there are no scrubbers to remove old contents from stmp* directories.
# After completion of regression tests, will remove the regression test subdirecories
export clean=".true."
export clean=".false."
;;
Jet)

Expand Down
143 changes: 29 additions & 114 deletions src/gsi/berror.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,6 @@ module berror
! 2011-04-07 todling - move newpc4pred to radinfo
! 2012-10-09 Gu - add fut2ps to project unbalanced temp to surface pressure in static B modeling
! 2013-05-27 zhu - add background error variances for aircraft temperature bias correction coefficients
! 2013-10-02 zhu - add reset_predictors_var
!
! subroutines included:
! sub init_berror - initialize background error related variables
Expand Down Expand Up @@ -120,7 +119,6 @@ module berror
public :: create_berror_vars
public :: destroy_berror_vars
public :: set_predictors_var
public :: reset_predictors_var
public :: init_rftable
public :: initable
public :: create_berror_vars_reg
Expand Down Expand Up @@ -385,7 +383,6 @@ subroutine set_predictors_var

integer(i_kind) i,j,ii
real(r_kind) stndev
real(r_kind) obs_count
logical new_tail

stndev = one/biaspredvar
Expand All @@ -407,17 +404,12 @@ subroutine set_predictors_var
do j=1,npred
ii=ii+1
if (inew_rad(i)) then
varprd(ii)=10000.0_r_kind
varA(j,i)=r10
else
if (ostats(i)<=20.0_r_kind) then
varA(j,i)=two*varA(j,i)+1.0e-6_r_kind
varprd(ii)=varA(j,i)
else
varprd(ii)=1.1_r_kind*varA(j,i)+1.0e-6_r_kind
end if
if (varprd(ii)>r10) varprd(ii)=r10
if (varA(j,i)>10000.0_r_kind) varA(j,i)=10000.0_r_kind
varA(j,i)=1.1_r_kind*varA(j,i)+1.0e-6_r_kind
varA(j,i)= min(r10,varA(j,i))
end if
varprd(ii)=varA(j,i)
end do
end do

Expand All @@ -427,50 +419,33 @@ subroutine set_predictors_var
do j=1,npredt
ii=ii+1

if (aircraft_t_bc_pof) then
obs_count = ostats_t(j,i)
new_tail = varA_t(j,i)==zero
end if
new_tail = varA_t(j,i)==zero
if (aircraft_t_bc) then
obs_count = ostats_t(1,i)
new_tail = .true.
if (any(varA_t(:,i)/=zero)) new_tail = .false.
end if

if (new_tail) then
varprd(ii)=one_tenth*one_tenth
if (aircraft_t_bc .and. j==2) varprd(ii)=1.0e-4_r_kind
if (aircraft_t_bc .and. j==3) varprd(ii)=1.0e-5_r_kind
else
if (obs_count<=10.0_r_kind) then
if (aircraft_t_bc .and. j==2) then
varA_t(j,i)=1.01_r_kind*varA_t(j,i)+1.0e-6_r_kind
else if (aircraft_t_bc .and. j==3) then
varA_t(j,i)=1.01_r_kind*varA_t(j,i)+1.0e-7_r_kind
else
varA_t(j,i)=1.01_r_kind*varA_t(j,i)+1.0e-5_r_kind
end if
varprd(ii)=varA_t(j,i)
if (aircraft_t_bc .and. j==2) then
varA_t(j,i)=1.0e-4_r_kind
else if (aircraft_t_bc .and. j==3) then
varA_t(j,i)=1.0e-5_r_kind
else
if (aircraft_t_bc .and. j==2) then
varprd(ii)=1.005_r_kind*varA_t(j,i)+1.0e-6_r_kind
else if (aircraft_t_bc .and. j==3) then
varprd(ii)=1.005_r_kind*varA_t(j,i)+1.0e-7_r_kind
else
varprd(ii)=1.005_r_kind*varA_t(j,i)+1.0e-5_r_kind
end if
varA_t(j,i)=one_tenth*one_tenth
end if
if (varprd(ii)>one_tenth) varprd(ii)=one_tenth
if (varA_t(j,i)>one_tenth) varA_t(j,i)=one_tenth
else
if (aircraft_t_bc .and. j==2) then
if (varprd(ii)>1.0e-3_r_kind) varprd(ii)=1.0e-3_r_kind
if (varA_t(j,i)>1.0e-3_r_kind) varA_t(j,i)=1.0e-3_r_kind
end if
if (aircraft_t_bc .and. j==3) then
if (varprd(ii)>1.0e-4_r_kind) varprd(ii)=1.0e-4_r_kind
if (varA_t(j,i)>1.0e-4_r_kind) varA_t(j,i)=1.0e-4_r_kind
varA_t(j,i)=1.005_r_kind*varA_t(j,i)+1.0e-6_r_kind
varA_t(j,i)=min(varA_t(j,i),1.0e-3_r_kind)
else if (aircraft_t_bc .and. j==3) then
varA_t(j,i)=1.005_r_kind*varA_t(j,i)+1.0e-7_r_kind
varA_t(j,i)=min(varA_t(j,i),1.0e-4_r_kind)
else
varA_t(j,i)=1.005_r_kind*varA_t(j,i)+1.0e-5_r_kind
varA_t(j,i)=min(varA_t(j,i),one_tenth)
end if
end if
varprd(ii)=varA_t(j,i)
end do
end do
end if
Expand All @@ -479,70 +454,6 @@ subroutine set_predictors_var
return
end subroutine set_predictors_var


subroutine reset_predictors_var
!$$$ subprogram documentation block
! . . . .
! subprogram: reset_predictors_var sets variances for bias correction predictors
! prgmmr: yanqiu org: np20 date: 2013-10-01
!
! abstract: resets variances for bias correction predictors
!
! program history log:
! output argument list:
! 2013-10-01 zhu
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use constants, only: one,one_tenth
use radinfo, only: newpc4pred,jpch_rad,npred,ostats,inew_rad,iuse_rad
use aircraftinfo, only: aircraft_t_bc_pof,aircraft_t_bc,biaspredt,ntail,npredt,ostats_t
use gridmod, only: twodvar_regional
use jfunc, only: nrclen, ntclen
implicit none

integer(i_kind) i,j,ii,obs_count
real(r_kind) stndev

stndev = one/biaspredt

! reset variances for bias predictor coeff. based on current data count
if (.not. twodvar_regional .and. newpc4pred) then
ii=0
do i=1,jpch_rad
do j=1,npred
ii=ii+1
if (.not.inew_rad(i) .and. iuse_rad(i)>0 .and. ostats(i)<=20.0_r_kind) then
varprd(ii)=1.0e-6_r_kind
end if
end do
end do

if ((aircraft_t_bc_pof .or. aircraft_t_bc) .and. ntclen>0) then
ii=nrclen-ntclen
do i=1,ntail
do j=1,npredt
ii=ii+1

if (aircraft_t_bc_pof) obs_count = ostats_t(j,i)
if (aircraft_t_bc) obs_count = ostats_t(1,i)

if (obs_count<=10.0_r_kind .and. varprd(ii)>stndev) then
varprd(ii)=stndev
if (aircraft_t_bc .and. j==2) varprd(ii)=one_tenth*stndev
if (aircraft_t_bc .and. j==3) varprd(ii)=one_tenth*one_tenth*stndev
end if
end do
end do
end if
end if

return
end subroutine reset_predictors_var

subroutine pcinfo
!$$$ subprogram documentation block
! . . . .
Expand Down Expand Up @@ -590,13 +501,12 @@ subroutine pcinfo
do i=1,jpch_rad
do j=1,npred
ii=ii+1
! if (ostats(i)>zero) vprecond(nclen1+ii)=vprecond(nclen1+ii)/(one+rstats(j,i)*varprd(ii))
if (ostats(i)>zero) vprecond(nclen1+ii)=one/(one+rstats(j,i)*varprd(ii))
if (ostats(i)>20.0_r_kind) then
if (rstats(j,i)>zero) then
varA(j,i)=one/(one/varprd(ii)+rstats(j,i))
else
varA(j,i)=10000.0_r_kind
if(varA(j,i) <= zero)varA(j,i)=10000.0_r_kind
end if
end if
end do
Expand All @@ -612,13 +522,18 @@ subroutine pcinfo
ii=ii+1
jj=jj+1

if (aircraft_t_bc_pof) obs_count = ostats_t(j,i)
if (aircraft_t_bc) obs_count = ostats_t(1,i)
obs_count=0
if (aircraft_t_bc_pof) then
obs_count = ostats_t(j,i)
else if (aircraft_t_bc) then
obs_count = ostats_t(1,i)
end if

! if (obs_count>zero) vprecond(nclen1+ii)=vprecond(nclen1+ii)/(one+rstats_t(j,i)*varprd(jj))
if (obs_count>zero) vprecond(nclen1+ii)=one/(one+rstats_t(j,i)*varprd(jj))
if (obs_count>3.0_r_kind) then
varA_t(j,i)=one/(one/varprd(jj)+rstats_t(j,i))
else
if(varA_t(j,i) <= zero)varA_t(j,i)=10000.0_r_kind
end if
end do
end do
Expand Down
5 changes: 2 additions & 3 deletions src/gsi/compute_derived.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,7 @@ subroutine compute_derived(mype,init_pass)
use kinds, only: r_kind,i_kind
use jfunc, only: jiter,jiterstart,&
qoption,switch_on_derivatives,&
tendsflag,clip_supersaturation
tendsflag,superfact,clip_supersaturation
use control_vectors, only: cvars3d
use control_vectors, only: nrf_var
use control_vectors, only: an_amp0
Expand Down Expand Up @@ -178,7 +178,6 @@ subroutine compute_derived(mype,init_pass)
if(init_pass .and. (ntguessig<1 .or. ntguessig>nfldsig)) &
call die(myname,'invalid init_pass, ntguessig =',ntguessig)


! Get required indexes from control vector names
nrf3_q=getindex(cvars3d,'q')
iq_loc=getindex(nrf_var,'q')
Expand Down Expand Up @@ -207,7 +206,7 @@ subroutine compute_derived(mype,init_pass)
end do
end do
end do

! Load guess cw for use in inner loop
! Get pointer to cloud water mixing ratio
it=ntguessig
Expand Down
4 changes: 2 additions & 2 deletions src/gsi/compute_qvar3d.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ subroutine compute_qvar3d
!$$$
use kinds, only: r_kind,i_kind,r_single
use berror, only: dssv
use jfunc, only: varq,qoption,varcw,cwoption,clip_supersaturation
use jfunc, only: varq,qoption,varcw,cwoption,clip_supersaturation,superfact
use derivsmod, only: qsatg,qgues
use control_vectors, only: cvars3d
use gridmod, only: lat2,lon2,nsig
Expand Down Expand Up @@ -94,7 +94,7 @@ subroutine compute_qvar3d
! Limit q to be >= qmin
ges_q(i,j,k)=max(ges_q(i,j,k),qmin)
! Limit q to be <= ges_qsat
if(clip_supersaturation) ges_q(i,j,k)=min(ges_q(i,j,k),ges_qsat(i,j,k,it))
if(clip_supersaturation) ges_q(i,j,k)=min(ges_q(i,j,k),superfact*ges_qsat(i,j,k,it))
end do
end do
end do
Expand Down
1 change: 0 additions & 1 deletion src/gsi/constants.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -238,7 +238,6 @@ module constants
real(r_kind),parameter:: ke2 = 0.00002_r_kind
real(r_kind),parameter:: row = r1000
real(r_kind),parameter:: rrow = one/row
! real(r_kind),parameter:: qmin = 1.e-7_r_kind !lower bound on ges_q

! Constant used to process ozone
real(r_kind),parameter:: constoz = 604229.0_r_kind
Expand Down
31 changes: 9 additions & 22 deletions src/gsi/crtm_interface.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -171,8 +171,6 @@ module crtm_interface
real(r_kind) , save ,allocatable,dimension(:,:) :: cloud_efr ! effective radius of cloud type in CRTM
real(r_kind) , save ,allocatable,dimension(:) :: cf ! effective radius of cloud type in CRTM
real(r_kind) , save ,allocatable,dimension(:) :: hwp_guess ! column total for each hydrometeor

real(r_kind) , save ,allocatable,dimension(:,:,:,:) :: gesqsat ! qsat to calc rh for aero particle size estimate
real(r_kind) , save ,allocatable,dimension(:) :: table,table2,tablew ! GFDL saturation water vapor pressure tables
real(r_kind) , save ,allocatable,dimension(:) :: des2,desw ! GFDL saturation water vapor presure
real(r_kind) , save ,allocatable,dimension(:) :: lcloud4crtm_wk ! cloud info usage index for each channel
Expand Down Expand Up @@ -840,16 +838,6 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo
endif ! nvege_type
endif ! regional or IGBP

! Calculate RH when aerosols are present and/or cloud-fraction used
if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0) then
allocate(gesqsat(lat2,lon2,nsig,nfldsig))
ice=.true.
iderivative=0
do ii=1,nfldsig
call genqsat(gesqsat(1,1,1,ii),ges_tsen(1,1,1,ii),ges_prsl(1,1,1,ii),lat2,lon2,nsig,ice,iderivative)
end do
endif

! Initial GFDL saturation water vapor pressure tables
if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0 .and. imp_physics==11) then

Expand Down Expand Up @@ -905,7 +893,6 @@ subroutine destroy_crtm
if (error_status /= success) &
write(6,*)myname_,': ***ERROR*** error_status=',error_status
if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0) then
deallocate(gesqsat)
if (imp_physics==11) then
deallocate(table)
deallocate(table2)
Expand Down Expand Up @@ -1048,7 +1035,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
use radinfo, only: nsigradjac
use gsi_nstcouplermod, only: nst_gsi
use guess_grids, only: ges_tsen,&
ges_prsl,ges_prsi,tropprs,dsfct,add_rtm_layers, &
ges_prsl,ges_prsi,ges_qsat,tropprs,dsfct,add_rtm_layers, &
hrdifsig,nfldsig,hrdifsfc,nfldsfc,ntguessfc,isli2,sno2, &
hrdifaer,nfldaer ! for separate aerosol input file
use cloud_efr_mod, only: efr_ql,efr_qi,efr_qr,efr_qs,efr_qg,efr_qh
Expand Down Expand Up @@ -1908,14 +1895,14 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
end if ! lread_ext_aerosol
end if ! n_actual_aerosols_wk > 0
do k=1,nsig
qs(k) = (gesqsat(ix ,iy ,k,itsig )*w00+ &
gesqsat(ixp,iy ,k,itsig )*w10+ &
gesqsat(ix ,iyp,k,itsig )*w01+ &
gesqsat(ixp,iyp,k,itsig )*w11)*dtsig + &
(gesqsat(ix ,iy ,k,itsigp)*w00+ &
gesqsat(ixp,iy ,k,itsigp)*w10+ &
gesqsat(ix ,iyp,k,itsigp)*w01+ &
gesqsat(ixp,iyp,k,itsigp)*w11)*dtsigp
qs(k) = (ges_qsat(ix ,iy ,k,itsig )*w00+ &
ges_qsat(ixp,iy ,k,itsig )*w10+ &
ges_qsat(ix ,iyp,k,itsig )*w01+ &
ges_qsat(ixp,iyp,k,itsig )*w11)*dtsig + &
(ges_qsat(ix ,iy ,k,itsigp)*w00+ &
ges_qsat(ixp,iy ,k,itsigp)*w10+ &
ges_qsat(ix ,iyp,k,itsigp)*w01+ &
ges_qsat(ixp,iyp,k,itsigp)*w11)*dtsigp
rh(k) = q(k)/qs(k)
end do
endif
Expand Down
8 changes: 4 additions & 4 deletions src/gsi/evalqlim.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ subroutine evalqlim(sval,pbc,rval)
use kinds, only: r_kind,i_kind,r_quad
use constants, only: zero,one,zero_quad
use gridmod, only: lat1,lon1,nsig,istart,wgtfactlats
use jfunc, only: factqmin,factqmax
use jfunc, only: factqmin,factqmax,superfact
use derivsmod, only: qgues,qsatg
use mpl_allreducemod, only: mpl_allreduce
use gsi_bundlemod, only: gsi_bundle
Expand Down Expand Up @@ -83,9 +83,9 @@ subroutine evalqlim(sval,pbc,rval)
endif
! Compute penalty for excess q
if (q>qsatg(i,j,k)) then
term=(factqmax*wgtfactlats(ii))*(q-qsatg(i,j,k))&
/(qsatg(i,j,k)*qsatg(i,j,k))
zbc(2) = zbc(2) + term*(q-qsatg(i,j,k))
term=(factqmax*wgtfactlats(ii))*((q-superfact*qsatg(i,j,k))&
/qsatg(i,j,k))**2
zbc(2) = zbc(2) + term
! Adjoint
rq(i,j,k) = rq(i,j,k) + term
endif
Expand Down
3 changes: 1 addition & 2 deletions src/gsi/genqsat.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,7 +121,6 @@ subroutine genqsat(qsat,tsen,prsl,lat2,lon2,nsig,ice,iderivative)
end do
do i=1,lat2
tdry = mint(i)
if( abs(tdry) < 1.0e-8_r_kind ) tdry = 1.0e-8_r_kind
tr = ttp/tdry
if (tdry >= ttp .or. .not. ice) then
estmax(i) = psat * (tr**xa) * exp(xb*(one-tr))
Expand All @@ -137,7 +136,6 @@ subroutine genqsat(qsat,tsen,prsl,lat2,lon2,nsig,ice,iderivative)
do k = 1,nsig
do i = 1,lat2
tdry = tsen(i,j,k)
if( abs(tdry) < 1.0e-8_r_kind ) tdry = 1.0e-8_r_kind
tr = ttp/tdry
if (tdry >= ttp .or. .not. ice) then
es = psat * (tr**xa) * exp(xb*(one-tr))
Expand All @@ -164,6 +162,7 @@ subroutine genqsat(qsat,tsen,prsl,lat2,lon2,nsig,ice,iderivative)
qsat(i,j,k) = max(qmin,qsat(i,j,k))

if(iderivative > 0)then
! if(es <= esmax .and. iderivative == 2 .and. qsat(i,j,k) > qmin )then
if(es <= esmax .and. iderivative == 2)then
idpupdate=.true.
idtupdate=.true.
Expand Down
Loading

0 comments on commit f938842

Please sign in to comment.