From 50190015445427a6671cf6060494ae50a8de1234 Mon Sep 17 00:00:00 2001 From: "john.derber" Date: Thu, 4 Nov 2021 14:44:13 +0000 Subject: [PATCH] GitHub Issue NOAA-EMC/GSI#219. Improve minimization and fix bug in vqc. --- src/gsi/berror.f90 | 147 ++++++---------------------- src/gsi/compute_derived.f90 | 7 +- src/gsi/compute_qvar3d.f90 | 4 +- src/gsi/constants.f90 | 1 - src/gsi/crtm_interface.f90 | 36 +++---- src/gsi/evalqlim.f90 | 8 +- src/gsi/genqsat.f90 | 3 +- src/gsi/gsimod.F90 | 8 +- src/gsi/intjcmod.f90 | 6 +- src/gsi/jfunc.f90 | 8 +- src/gsi/observer.F90 | 6 -- src/gsi/pcgsoi.f90 | 110 ++++++++++----------- src/gsi/prt_guess.f90 | 8 +- src/gsi/read_guess.F90 | 4 +- src/gsi/setupcldtot.F90 | 1 - src/gsi/setupq.f90 | 15 ++- src/gsi/setuprhsall.f90 | 5 - src/gsi/setupw.f90 | 29 +++--- src/gsi/stpcalc.f90 | 185 ++++++++++++++++++------------------ src/gsi/stpjcmod.f90 | 14 +-- src/gsi/stpw.f90 | 38 ++++---- src/gsi/update_guess.f90 | 4 +- 22 files changed, 268 insertions(+), 379 deletions(-) diff --git a/src/gsi/berror.f90 b/src/gsi/berror.f90 index 71bd4f1ea8..3a683c9b94 100644 --- a/src/gsi/berror.f90 +++ b/src/gsi/berror.f90 @@ -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 @@ -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 @@ -377,15 +375,14 @@ subroutine set_predictors_var ! !$$$ use constants, only: zero,one,two,one_tenth,r10 - use radinfo, only: ostats,varA,jpch_rad,npred,inew_rad,newpc4pred,biaspredvar - use aircraftinfo, only: aircraft_t_bc_pof,aircraft_t_bc,biaspredt,ntail,npredt,ostats_t,varA_t + use radinfo, only: varA,jpch_rad,npred,inew_rad,newpc4pred,biaspredvar + use aircraftinfo, only: aircraft_t_bc_pof,aircraft_t_bc,biaspredt,ntail,npredt,varA_t use gridmod, only: twodvar_regional use jfunc, only: nrclen, ntclen implicit none integer(i_kind) i,j,ii real(r_kind) stndev - real(r_kind) obs_count logical new_tail stndev = one/biaspredvar @@ -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 @@ -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 @@ -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 ! . . . . @@ -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 @@ -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 diff --git a/src/gsi/compute_derived.f90 b/src/gsi/compute_derived.f90 index a33dbe1f41..46413d9986 100644 --- a/src/gsi/compute_derived.f90 +++ b/src/gsi/compute_derived.f90 @@ -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 @@ -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') @@ -202,12 +201,12 @@ subroutine compute_derived(mype,init_pass) ! 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,ii)) + if(clip_supersaturation) ges_q(i,j,k) = min(ges_q(i,j,k),superfact*ges_qsat(i,j,k,ii)) end do end do end do end do - + ! Load guess cw for use in inner loop ! Get pointer to cloud water mixing ratio it=ntguessig diff --git a/src/gsi/compute_qvar3d.f90 b/src/gsi/compute_qvar3d.f90 index 54aa7721bf..851212d4f5 100644 --- a/src/gsi/compute_qvar3d.f90 +++ b/src/gsi/compute_qvar3d.f90 @@ -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 @@ -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 diff --git a/src/gsi/constants.f90 b/src/gsi/constants.f90 index da2368d70f..77aa027061 100644 --- a/src/gsi/constants.f90 +++ b/src/gsi/constants.f90 @@ -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 diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index d03715059f..3de679f7f4 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -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 @@ -326,7 +324,6 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo use radinfo, only: crtm_coeffs_path use radinfo, only: radjacindxs,radjacnames,jpch_rad,nusis,nuchan use aeroinfo, only: aerojacindxs - use guess_grids, only: ges_tsen,ges_prsl,nfldsig use gridmod, only: fv3_full_hydro use mpeu_util, only: getindex use constants, only: zero,max_varname_length @@ -347,9 +344,9 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo integer(i_kind), parameter :: length = 2621 ! lenth of GFL qsat table ! local variables - integer(i_kind) :: ier,ii,error_status,iderivative + integer(i_kind) :: ier,ii,error_status integer(i_kind) :: k, subset_start, subset_end - logical :: ice,Load_AerosolCoeff,Load_CloudCoeff + logical :: Load_AerosolCoeff,Load_CloudCoeff character(len=20),dimension(1) :: sensorlist integer(i_kind) :: indx,iii,icloud4crtm ! ...all "additional absorber" variables @@ -840,16 +837,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 @@ -905,7 +892,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) @@ -1048,7 +1034,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 @@ -1908,14 +1894,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 diff --git a/src/gsi/evalqlim.f90 b/src/gsi/evalqlim.f90 index 085fa87724..a5ad216cd8 100644 --- a/src/gsi/evalqlim.f90 +++ b/src/gsi/evalqlim.f90 @@ -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 @@ -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 diff --git a/src/gsi/genqsat.f90 b/src/gsi/genqsat.f90 index 2719ed28f8..ed0eb152e6 100644 --- a/src/gsi/genqsat.f90 +++ b/src/gsi/genqsat.f90 @@ -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)) @@ -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)) @@ -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. diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 48a6475eee..1e12efdd39 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -93,7 +93,7 @@ module gsimod pvis,pcldch,scale_cv,estvisoe,estcldchoe,vis_thres,cldch_thres,cao_check use qcmod, only: troflg,lat_c,nrand use pcpinfo, only: npredp,diag_pcp,dtphys,deltim,init_pcp - use jfunc, only: iout_iter,iguess,miter,factqmin,factqmax, & + use jfunc, only: iout_iter,iguess,miter,factqmin,factqmax,superfact,limitqobs, & factql,factqi,factqr,factqs,factqg, & factv,factl,factp,factg,factw10m,facthowv,factcldch,niter,niter_no_qc,biascor,& init_jfunc,qoption,cwoption,switch_on_derivatives,tendsflag,jiterstart,jiterend,R_option,& @@ -489,6 +489,8 @@ module gsimod ! gencode - source generation code ! factqmin - weighting factor for negative moisture constraint ! factqmax - weighting factor for supersaturated moisture constraint +! superfact- amount of supersaturation allowed 1.01 = 1% supersaturation +! limitqobs- limit q obs to be <= 100%RH based on model temperatures ! clip_supersaturation - flag to remove supersaturation during each outer loop default=.false. ! deltim - model timestep ! dtphys - physics timestep @@ -686,7 +688,7 @@ module gsimod ! NOTE: for now, if in regional mode, then iguess=-1 is forced internally. ! add use of guess file later for regional mode. - namelist/setup/gencode,factqmin,factqmax,clip_supersaturation, & + namelist/setup/gencode,factqmin,factqmax,superfact,limitqobs,clip_supersaturation, & factql,factqi,factqr,factqs,factqg, & factv,factl,factp,factg,factw10m,facthowv,factcldch,R_option,deltim,dtphys,& biascor,bcoption,diurnalbc,& @@ -1942,6 +1944,8 @@ subroutine gsimain_initialize dmesh=one factqmin=zero factqmax=zero + superfact=1._r_kind + limitqobs=.false. if (hilbert_curve) then write(6,*) 'Disabling hilbert_curve cross validation when oneobtest=.true.' hilbert_curve=.false. diff --git a/src/gsi/intjcmod.f90 b/src/gsi/intjcmod.f90 index d264cd4b5c..c0c23151ee 100644 --- a/src/gsi/intjcmod.f90 +++ b/src/gsi/intjcmod.f90 @@ -71,7 +71,7 @@ subroutine intlimq(rval,sval,itbin) ! !$$$ use gridmod, only: nsig,lat1,lon1,istart,wgtfactlats - use jfunc, only: factqmin,factqmax + use jfunc, only: factqmin,factqmax,superfact use gsi_metguess_mod, only: gsi_metguess_bundle use guess_grids, only: ges_qsat use mpimod, only: mype @@ -116,8 +116,8 @@ subroutine intlimq(rval,sval,itbin) /(ges_qsat(i,j,k,itbin)*ges_qsat(i,j,k,itbin)) ! Upper constraint limit - else if (q > ges_qsat(i,j,k,itbin)) then - rq(i,j,k) = rq(i,j,k) + (factqmax*wgtfactlats(ii))*(q-ges_qsat(i,j,k,itbin))/ & + else if (q > superfact*ges_qsat(i,j,k,itbin)) then + rq(i,j,k) = rq(i,j,k) + (factqmax*wgtfactlats(ii))*(q-superfact*ges_qsat(i,j,k,itbin))/ & (ges_qsat(i,j,k,itbin)*ges_qsat(i,j,k,itbin)) end if diff --git a/src/gsi/jfunc.f90 b/src/gsi/jfunc.f90 index 1f1adc5702..d65b4c8fd3 100644 --- a/src/gsi/jfunc.f90 +++ b/src/gsi/jfunc.f90 @@ -130,7 +130,7 @@ module jfunc public :: switch_on_derivatives,jiterend,jiterstart,jiter,iter,niter,miter public :: diurnalbc,bcoption,biascor,nval2d,xhatsave,first public :: factqmax,factqmin,clip_supersaturation,last,yhatsave,nvals_len,nval_levs,nval_levs_ens,iout_iter,nclen - public :: factql,factqi,factqr,factqs,factqg + public :: factql,factqi,factqr,factqs,factqg,superfact,limitqobs public :: niter_no_qc,print_diag_pcg,penorig,gnormorig,iguess public :: factg,factv,factp,factl,R_option,factw10m,facthowv,factcldch,diag_precon,step_start public :: pseudo_q2 @@ -139,7 +139,7 @@ module jfunc logical first,last,switch_on_derivatives,tendsflag,print_diag_pcg,tsensible,diag_precon logical clip_supersaturation,R_option - logical pseudo_q2 + logical pseudo_q2,limitqobs logical cnvw_option integer(i_kind) iout_iter,miter,iguess,nclen,qoption,cwoption integer(i_kind) jiter,jiterstart,jiterend,iter @@ -150,7 +150,7 @@ module jfunc integer(i_kind),dimension(0:50):: niter,niter_no_qc real(r_kind) factqmax,factqmin,gnormorig,penorig,biascor(2),diurnalbc,factg,factv,factp,factl,& - factw10m,facthowv,factcldch,step_start + factw10m,facthowv,factcldch,step_start,superfact real(r_kind) factql,factqi,factqr,factqs,factqg integer(i_kind) bcoption real(r_kind),allocatable,dimension(:,:):: varq @@ -202,6 +202,8 @@ subroutine init_jfunc factqmin=zero factqmax=zero + superfact=1.00_r_kind + limitqobs=.false. factql=zero factqi=zero factqr=zero diff --git a/src/gsi/observer.F90 b/src/gsi/observer.F90 index 7ac91e5eed..00f51448ac 100644 --- a/src/gsi/observer.F90 +++ b/src/gsi/observer.F90 @@ -173,12 +173,6 @@ subroutine guess_init_ #endif /* _LAG_MODEL_ */ endif -! Read output from previous min. - if (l4dvar.and.jiterstart>1) then - else - ! If requested and if available, read guess solution. - endif - ! Generate coefficients for compact differencing if(.not.regional)then if(.not.cdiff_created()) call create_cdiff_coefs() diff --git a/src/gsi/pcgsoi.f90 b/src/gsi/pcgsoi.f90 index bce3aa0835..0f7f3d32d8 100644 --- a/src/gsi/pcgsoi.f90 +++ b/src/gsi/pcgsoi.f90 @@ -170,15 +170,14 @@ subroutine pcgsoi() character(5) step(2) integer(i_kind) i,istep,iobs,ii,nprt real(r_kind) stp,b,converge - real(r_kind) gsave,small_step + real(r_kind) gsave,small_step,aindex real(r_kind) gnormx,penx,penalty,penaltynew - real(r_double) pennorm - real(r_quad) zjo - real(r_quad) :: zdla - real(r_quad),dimension(4):: dprod - real(r_kind),dimension(3):: gnorm real(r_kind) :: zgini,zfini,fjcost(4),fjcostnew(4),zgend,zfend real(r_kind) :: fjcost_e + real(r_kind),dimension(3):: gnorm + real(r_double) pennorm + real(r_quad) :: zdla,zjo + real(r_quad),dimension(4):: dprod type(control_vector) :: gradx,grady,dirx,diry,ydiff,xdiff type(gsi_bundle) :: sval(nobs_bins), rval(nobs_bins) type(gsi_bundle) :: eval(ntlevs_ens) @@ -322,7 +321,7 @@ subroutine pcgsoi() end if ! 2. Multiply by background error - call multb(lanlerr,gradx,grady) + call multb(gradx,grady) if(iorthomax>0) then ! save gradients @@ -338,34 +337,23 @@ subroutine pcgsoi() ! 3. Calculate new norm of gradients and factors going into b calculation dprod(1) = qdot_prod_sub(gradx,grady) - if(iter > 0)then - gsave=gnorm(3) - if (lanlerr) then -! xdiff used as a temporary array - do i=1,nclen - xdiff%values(i)=vprecond(i)*gradx%values(i) - end do - dprod(2) = qdot_prod_sub(xdiff,grady) - call mpl_allreduce(2,qpvals=dprod) - gnorm(2)=dprod(2) - gnorm(3)=dprod(2) - else - do i=1,nclen - xdiff%values(i)=vprecond(i)*(gradx%values(i)-xdiff%values(i)) - ydiff%values(i)=vprecond(i)*(grady%values(i)-ydiff%values(i)) - end do - dprod(2) = qdot_prod_sub(xdiff,grady) - dprod(3) = qdot_prod_sub(ydiff,gradx) + if(iter > 0 .and. .not. lanlerr)then + dprod(3) = qdot_prod_sub(xdiff,grady) + dprod(4) = qdot_prod_sub(ydiff,gradx) ! xdiff used as a temporary array - do i=1,nclen - xdiff%values(i)=vprecond(i)*gradx%values(i) - end do - dprod(4) = qdot_prod_sub(xdiff,grady) - call mpl_allreduce(4,qpvals=dprod) -! Two dot products in gnorm(2) should be same, but are slightly -! different due to round off, so use average. - gnorm(2)=0.5_r_quad*(dprod(2)+dprod(3)) - gnorm(3)=dprod(4) + do i=1,nclen + xdiff%values(i)=vprecond(i)*gradx%values(i) + end do + dprod(2) = qdot_prod_sub(xdiff,grady) + call mpl_allreduce(4,qpvals=dprod) +! Two dot products in dprod(3) and dprod(4) should be same, but are slightly +! different due to round off, so use average. + gnorm(2)=dprod(2)-0.5_r_quad*(dprod(3)+dprod(4)) + gnorm(3)=dprod(2) + if(mype == 0)then + aindex=abs(dprod(3)/dprod(2)) + write(iout_iter,*) 'NL Index ',aindex + if(aindex > 0.5_r_kind .or. print_verbose) write(iout_iter,*) 'NL Values ', dprod(3),dprod(2) end if else ! xdiff used as a temporary array @@ -375,38 +363,36 @@ subroutine pcgsoi() dprod(2) = qdot_prod_sub(xdiff,grady) call mpl_allreduce(2,qpvals=dprod) if(print_diag_pcg) call prt_control_norms(grady,'grady') + gnorm(2)=dprod(2) gnorm(3)=dprod(2) end if gnorm(1)=dprod(1) - if(mype == 0)write(iout_iter,*)'Minimization iteration',iter + if(mype == 0)write(iout_iter,*)'Minimization iteration',iter ! 4. Calculate b and new search direction b=zero if (.not. restart .or. iter > 0) then - if (gsave>1.e-16_r_kind .and. iter>0) b=gnorm(2)/gsave - if (b7.0_r_kind) then - if (mype==0) then - if (iout_6) write(6,105) gnorm(2),gsave,b - write(iout_iter,105) gnorm(2),gsave,b + if (iter > 1 .or. .not. read_success)then + if (gsave>1.e-16_r_kind) b=gnorm(2)/gsave + if (b10.0_r_kind) then + if (mype==0) then + if (iout_6) write(6,105) gnorm(2),gsave,b + write(iout_iter,105) gnorm(2),gsave,b + endif + b=zero endif - b=zero - endif - if (mype==0 .and. print_verbose) write(6,888)'pcgsoi: gnorm(1:3),b=',gnorm,b - if(read_success .and. iter == 1)b=zero - - if (.not. lanlerr) then - do i=1,nclen - xdiff%values(i)=gradx%values(i) - ydiff%values(i)=grady%values(i) - end do + if (mype==0 .and. print_verbose) write(6,888)'pcgsoi: gnorm(1:3),b=',gnorm,b end if -! Calculate new search direction + do i=1,nclen - dirx%values(i)=-vprecond(i)*grady%values(i)+b*dirx%values(i) - diry%values(i)=-vprecond(i)*gradx%values(i)+b*diry%values(i) +! Calculate new search direction + ydiff%values(i)=vprecond(i)*grady%values(i) + dirx%values(i)=-ydiff%values(i)+b*dirx%values(i) + xdiff%values(i)=vprecond(i)*gradx%values(i) + diry%values(i)=-xdiff%values(i)+b*diry%values(i) end do else ! If previous solution available, transfer into local arrays. @@ -416,9 +402,11 @@ subroutine pcgsoi() end do call read_guess_solution(diry,mype,read_success) ! Multiply by background error - call multb(lanlerr,diry,dirx) + call multb(diry,dirx) + restart=.false. endif - + gsave=gnorm(3) + ! 5. Calculate stepsize and update solution ! Convert search direction from control space to physical space do ii=1,nobs_bins @@ -568,7 +556,7 @@ subroutine pcgsoi() end do ! Multiply by background error - call multb(lanlerr,gradx,grady) + call multb(gradx,grady) ! Print final Jo table zgend=dot_product(gradx,grady,r_quad) @@ -630,8 +618,7 @@ subroutine pcgsoi() ! if (mype==0) write(6,*)'pcgsoi: Updating guess' if(iwrtinc<=0) call update_guess(sval,sbias) -! cloud analysis after iteration -! if(jiter == miter .and. i_gsdcldanal_type==1) then +! gsd cloud analysis after iteration if(jiter == miter) then if(i_gsdcldanal_type==2) then call gsdcloudanalysis4nmmb(mype) @@ -833,7 +820,7 @@ subroutine periodic_(gradx) end subroutine periodic_ -subroutine multb(lanlerr,vec1,vec2) +subroutine multb(vec1,vec2) !$$$ subprogram documentation block ! . . . . ! subprogram: multb multiply vec1 by background error to equal vec2 @@ -863,7 +850,6 @@ subroutine multb(lanlerr,vec1,vec2) type(control_vector),intent(inout) :: vec1 type(control_vector),intent(inout) :: vec2 - logical ,intent(in ) :: lanlerr if(periodic)call periodic_(vec1) ! start by setting vec2=vec1 and then operate on vec2 (unless gram_schmidt) @@ -892,7 +878,7 @@ end subroutine multb subroutine c2s(hat,val,bias,llprt,ltest) !$$$ subprogram documentation block ! . . . . -! subprogram: multb control2state for all options +! subprogram: c2s control2state for all options ! prgmmr: derber ! ! abstract: generalized control2state @@ -958,7 +944,7 @@ end subroutine c2s subroutine c2s_ad(hat,val,bias,llprt) !$$$ subprogram documentation block ! . . . . -! subprogram: multb control2state for all options +! subprogram: c2s_ad adjoint of control2state for all options ! prgmmr: derber ! ! abstract: generalized control2state diff --git a/src/gsi/prt_guess.f90 b/src/gsi/prt_guess.f90 index ad5cb80025..e1e251195a 100644 --- a/src/gsi/prt_guess.f90 +++ b/src/gsi/prt_guess.f90 @@ -165,7 +165,7 @@ subroutine prt_guess(sgrep) zloc(nvars+7) = minval(ges_cwmr_it(2:lat1+1,2:lon1+1,1:nsig)) zloc(nvars+8) = minval(ges_cf_it(2:lat1+1,2:lon1+1,1:nsig)) zloc(nvars+9) = minval(ges_div_it(2:lat1+1,2:lon1+1,1:nsig)) - zloc(nvars+10) = minval(ges_vor_it(2:lat1+1,2:lon1+1,1:nsig)) + zloc(nvars+10) = minval(ges_vor_it(2:lat1+1,2:lon1+1,1:nsig)) zloc(nvars+11) = minval(ges_prsl (2:lat1+1,2:lon1+1,1:nsig,ntsig)) zloc(nvars+12) = minval(ges_ps_it (2:lat1+1,2:lon1+1 )) zloc(nvars+13) = minval(sfct (2:lat1+1,2:lon1+1, ntsfc)) @@ -178,7 +178,7 @@ subroutine prt_guess(sgrep) zloc(2*nvars+7) = maxval(ges_cwmr_it(2:lat1+1,2:lon1+1,1:nsig)) zloc(2*nvars+8) = maxval(ges_cf_it(2:lat1+1,2:lon1+1,1:nsig)) zloc(2*nvars+9) = maxval(ges_div_it(2:lat1+1,2:lon1+1,1:nsig)) - zloc(2*nvars+10) = maxval(ges_vor_it(2:lat1+1,2:lon1+1,1:nsig)) + zloc(2*nvars+10) = maxval(ges_vor_it(2:lat1+1,2:lon1+1,1:nsig)) zloc(2*nvars+11) = maxval(ges_prsl (2:lat1+1,2:lon1+1,1:nsig,ntsig)) zloc(2*nvars+12) = maxval(ges_ps_it (2:lat1+1,2:lon1+1 )) zloc(2*nvars+13) = maxval(sfct (2:lat1+1,2:lon1+1, ntsfc)) @@ -188,8 +188,8 @@ subroutine prt_guess(sgrep) ! Gather contributions - call mpi_allgather(zloc,3*nvars+3,mpi_rtype, & - & zall,3*nvars+3,mpi_rtype, mpi_comm_world,ierror) + call mpi_gather(zloc,3*nvars+3,mpi_rtype, & + & zall,3*nvars+3,mpi_rtype,0, mpi_comm_world,ierror) if (mype==0) then zmin=zero diff --git a/src/gsi/read_guess.F90 b/src/gsi/read_guess.F90 index 364d4e305b..ecd85aa983 100644 --- a/src/gsi/read_guess.F90 +++ b/src/gsi/read_guess.F90 @@ -89,7 +89,7 @@ subroutine read_guess(iyear,month,idd,mype) !$$$ use kinds, only: r_kind,i_kind - use jfunc, only: bcoption,clip_supersaturation + use jfunc, only: bcoption,clip_supersaturation,superfact use guess_grids, only: nfldsig,ges_tsen,load_prsges,load_geop_hgt,ges_prsl,& ges_tsen1, geop_hgti, ges_geopi, ges_q1 use m_gsiBiases,only : bkg_bias_correction,nbc @@ -250,7 +250,7 @@ subroutine read_guess(iyear,month,idd,mype) do k=1,nsig do j=1,lon2 do i=1,lat2 - satval = min(ges_q(i,j,k),satq(i,j,k)) + satval = min(ges_q(i,j,k),superfact*satq(i,j,k)) satval = max(qmin,satval) ges_q(i,j,k) = satval ges_tsen(i,j,k,it)= ges_tv(i,j,k)/(one+fv*ges_q(i,j,k)) diff --git a/src/gsi/setupcldtot.F90 b/src/gsi/setupcldtot.F90 index fc3978e3d2..565075bec3 100755 --- a/src/gsi/setupcldtot.F90 +++ b/src/gsi/setupcldtot.F90 @@ -109,7 +109,6 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di external:: tintrp2a1,tintrp2a11 external:: tintrp31,tintrp3 external:: grdcrd1 - external:: genqsat external:: stop2 ! Declare local variables diff --git a/src/gsi/setupq.f90 b/src/gsi/setupq.f90 index 77c43ba121..8189b0bcb5 100644 --- a/src/gsi/setupq.f90 +++ b/src/gsi/setupq.f90 @@ -149,13 +149,13 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav use nc_diag_read_mod, only: nc_diag_read_init, nc_diag_read_get_dim, nc_diag_read_close use gsi_4dvar, only: nobs_bins,hr_obsbin,min_offset use oneobmod, only: oneobtest,maginnov,magoberr - use guess_grids, only: ges_lnprsl,hrdifsig,nfldsig,ges_tsen,ges_prsl,pbl_height + use guess_grids, only: ges_lnprsl,hrdifsig,nfldsig,ges_tsen,ges_prsl,pbl_height,ges_qsat use gridmod, only: lat2,lon2,nsig,get_ijk,twodvar_regional use constants, only: zero,one,r1000,r10,r100 use constants, only: huge_single,wgtlim,three use constants, only: tiny_r_kind,five,half,two,huge_r_kind,r0_01 use qcmod, only: npres_print,ptopq,pbotq,dfact,dfact1,njqc,vqc,nvqc - use jfunc, only: jiter,last,jiterstart,miter + use jfunc, only: jiter,last,jiterstart,miter,superfact,limitqobs use convinfo, only: nconvtype,cermin,cermax,cgross,cvar_b,cvar_pg,ictype use convinfo, only: ibeta,ikapa use convinfo, only: icsubtype @@ -417,7 +417,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav iderivative=0 do jj=1,nfldsig call genqsat(qg(1,1,1,jj),ges_tsen(1,1,1,jj),ges_prsl(1,1,1,jj),lat2,lon2,nsig,ice,iderivative) - call genqsat(qg2m(1,1,jj),ges_tsen(1,1,1,jj),ges_prsl(1,1,1,jj),lat2,lon2, 1,ice,iderivative) + qg2m(:,:,jj)=qg(:,:,1,jj) end do @@ -516,9 +516,15 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! Scale errors by guess saturation q - call tintrp31(qg,qsges,dlat,dlon,dpres,dtime,hrdifsig,& + qob = data(iqob,i) + if(limitqobs) then + call tintrp31(ges_qsat,qsges,dlat,dlon,dpres,dtime,hrdifsig,& mype,nfldsig) + qob=min(qob,superfact*qsges) + end if + call tintrp31(qg,qsges,dlat,dlon,dpres,dtime,hrdifsig,& + mype,nfldsig) ! Interpolate 2-m qs to obs locations/times if((i_use_2mq4b > 0) .and. ((itype > 179 .and. itype < 190) .or. itype == 199) & .and. .not.twodvar_regional)then @@ -527,7 +533,6 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! Load obs error and value into local variables obserror = max(cermin(ikx)*r0_01,min(cermax(ikx)*r0_01,data(ier,i))) - qob = data(iqob,i) rmaxerr=rmaxerr*qsges rmaxerr=max(small2,rmaxerr) diff --git a/src/gsi/setuprhsall.f90 b/src/gsi/setuprhsall.f90 index 8e34876201..3efcb69859 100644 --- a/src/gsi/setuprhsall.f90 +++ b/src/gsi/setuprhsall.f90 @@ -150,7 +150,6 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass) use lag_fields, only: lag_presetup,lag_state_write,lag_state_read,lag_destroy_uv use mpeu_util, only: getindex use mpl_allreducemod, only: mpl_allreduce - use berror, only: reset_predictors_var use rapidrefresh_cldsurf_mod, only: l_PBL_pseudo_SurfobsT,l_PBL_pseudo_SurfobsQ,& l_PBL_pseudo_SurfobsUV,i_gsdcldanal_type,& i_cloud_q_innovation @@ -579,10 +578,6 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass) call mpl_allreduce(npredt,ntail,rstats_t) end if - if (newpc4pred .or. aircraft_t_bc_pof .or. aircraft_t_bc) then - call reset_predictors_var - end if - ! Collect satellite and precip. statistics call mpi_reduce(aivals,aivals1,size(aivals1),mpi_rtype,mpi_sum,mype_rad, & mpi_comm_world,ierror) diff --git a/src/gsi/setupw.f90 b/src/gsi/setupw.f90 index aa6fe55094..33f0a62957 100644 --- a/src/gsi/setupw.f90 +++ b/src/gsi/setupw.f90 @@ -253,9 +253,9 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav real(r_double) rstation_id real(r_kind) qcu,qcv,trop5,tfact,fact real(r_kind) scale,ratio,obserror,obserrlm - real(r_kind) residual,ressw,ress,val,vals,val2,valqc2,dudiff,dvdiff - real(r_kind) valqc,valu,valv,dx10,rlow,rhgh,drpx,prsfc,var_jb - real(r_kind) cg_t,cvar,wgt,term,rat_err2,qcgross + real(r_kind) residual,ressw,ress,vals,val2,dudiff,dvdiff,rat_err2u + real(r_kind) valqc,valu,valv,dx10,rlow,rhgh,drpx,prsfc,var_jb,rat_err2v + real(r_kind) cg_t,cvar,wgt,term,qcgross,valqcu,valqcv real(r_kind) presw,factw,dpres,ugesin,vgesin,rwgt,dpressave real(r_kind) sfcchk,prsln2,error,dtime,dlon,dlat,r0_001,rsig,thirty,rsigp real(r_kind) ratio_errors,goverrd,spdges,spdob,ten,psges,zsges @@ -1222,8 +1222,6 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ratio_errors=ratio_errors*sqrt(data(ihil,i)) ! hilbert weight ! Compute penalty terms (linear & nonlinear qc). if(luse(i))then - val = valu*valu+valv*valv - vals=sqrt(val) if(vqc) then cg_t=cvar_b(ikx) cvar=cvar_pg(ikx) @@ -1238,17 +1236,23 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ibb=0 ikk=0 endif + vals=valu call vqc_setup(vals,ratio_errors,error,cvar,& - cg_t,ibb,ikk,var_jb,rat_err2,wgt,valqc) - rwgt = wgt/wgtlim + cg_t,ibb,ikk,var_jb,rat_err2u,wgt,valqcu) + rwgt = 0.5_r_kind*wgt/wgtlim + vals=valv + call vqc_setup(vals,ratio_errors,error,cvar,& + cg_t,ibb,ikk,var_jb,rat_err2v,wgt,valqcv) + rwgt = rwgt+0.5_r_kind*wgt/wgtlim + valqc=valqcu+valqcv ! Accumulate statistics for obs belonging to this task if (muse(i)) then if(rwgt < one) awork(21) = awork(21)+one jsig = dpres jsig=max(1,min(jsig,nsig)) - awork(4*nsig+jsig+100)=awork(4*nsig+jsig+100)+valu*valu*rat_err2 - awork(5*nsig+jsig+100)=awork(5*nsig+jsig+100)+valv*valv*rat_err2 + awork(4*nsig+jsig+100)=awork(4*nsig+jsig+100)+valu*valu*rat_err2u + awork(5*nsig+jsig+100)=awork(5*nsig+jsig+100)+valv*valv*rat_err2v awork(6*nsig+jsig+100)=awork(6*nsig+jsig+100)+one awork(3*nsig+jsig+100)=awork(3*nsig+jsig+100)+valqc end if @@ -1257,8 +1261,7 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! as a function of observation type. ress = scale*sqrt(dudiff**2+dvdiff**2) ressw = ress*ress - val2 = half*(valu*valu+valv*valv) - valqc2 = half*valqc + val2 = half*(valu*valu*rat_err2u+valv*valv*rat_err2v) nn=1 if (.not. muse(i)) then nn=2 @@ -1269,8 +1272,8 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav bwork(k,ikx,1,nn) = bwork(k,ikx,1,nn)+one ! count bwork(k,ikx,2,nn) = bwork(k,ikx,2,nn)+spdb ! speed bias bwork(k,ikx,3,nn) = bwork(k,ikx,3,nn)+ressw ! (o-g)**2 - bwork(k,ikx,4,nn) = bwork(k,ikx,4,nn)+val2*rat_err2 ! penalty - bwork(k,ikx,5,nn) = bwork(k,ikx,5,nn)+valqc2 ! nonlin qc penalty + bwork(k,ikx,4,nn) = bwork(k,ikx,4,nn)+val2 ! penalty + bwork(k,ikx,5,nn) = bwork(k,ikx,5,nn)+valqc ! nonlin qc penalty end if end do diff --git a/src/gsi/stpcalc.f90 b/src/gsi/stpcalc.f90 index 8b67e35742..80fac64d61 100644 --- a/src/gsi/stpcalc.f90 +++ b/src/gsi/stpcalc.f90 @@ -280,7 +280,7 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & real(r_kind),dimension(4)::sges real(r_kind),dimension(ioutpen):: outpen,outstp logical :: cxterm,change_dels,ifound - logical :: print_verbose + logical :: print_verbose,pjcalc !************************************************************************************ @@ -296,6 +296,8 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & outpen = zero nsteptot=0 istp_use=0 + kprt=3 + pjcalc=.false. pj=zero_quad ! Begin calculating contributions to penalty and stepsize for various terms @@ -386,12 +388,13 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & pstart=zero_quad pbc=zero_quad + if(iter == 0 .and. kprt >= 2)pjcalc=.true. ! penalty, b and c for background terms pstart(1,1) = qdot_prod_sub(xhatsave,yhatsave) - pj(1,1)=pstart(1,1) + if(pjcalc)pj(1,1)=pstart(1,1) ! two terms in next line should be the same, but roundoff makes average more accurate. @@ -406,7 +409,7 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & if (ljcdfi .and. nobs_bins>1) then call stpjcdfi(dval,sval,pstart(1,2),pstart(2,2),pstart(3,2)) - pj(2,1)=pstart(1,2) + if(pjcalc)pj(2,1)=pstart(1,2) end if ! Penalty, b, c for dry pressure @@ -416,13 +419,15 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & else call stpjcpdry(dval,sval,pstart(1,3),pstart(2,3),pstart(3,3),nobs_bins) end if - pj(3,1)=pstart(1,3) + if(pjcalc)pj(3,1)=pstart(1,3) end if ! iterate over number of stepsize iterations (istp_iter - currently set to a maximum of 5) dels = one_tenth_quad stepsize: do ii=1,istp_iter + pjcalc=.false. + if(iter == 0 .and. kprt >= 2 .and. ii == 1)pjcalc=.true. iis=ii ! Delta stepsize change_dels=.true. @@ -457,11 +462,76 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & ! penalties for moisture constraint if(.not. ltlint)then + if(.not.ljc4tlevs) then + call stplimq(dval(ibin_anl),sval(ibin_anl),sges,pbc(1,4),pbc(1,5),nstep,ntguessig) + if(pjcalc)then + pj(4,1)=pbc(1,4)+pbc(ipenloc,4) + pj(5,1)=pbc(1,5)+pbc(ipenloc,5) + end if + else + do ibin=1,nobs_bins + if (nobs_bins /= nfldsig) then + it=ntguessig + else + it=ibin + end if + call stplimq(dval(ibin),sval(ibin),sges,pbcqmin(1,ibin),pbcqmax(1,ibin),nstep,it) + end do + pbc(:,4)=zero_quad + pbc(:,5)=zero_quad + do ibin=1,nobs_bins + do j=1,nstep + pbc(j,4) = pbc(j,4)+pbcqmin(j,ibin) + pbc(j,5) = pbc(j,5)+pbcqmax(j,ibin) + end do + end do + if(pjcalc)then + do ibin=1,nobs_bins + pj(4,ibin)=pj(4,ibin)+pbcqmin(1,ibin)+pbcqmin(ipenloc,ibin) + pj(5,ibin)=pj(5,ibin)+pbcqmax(1,ibin)+pbcqmax(ipenloc,ibin) + end do + end if + end if +! penalties for gust constraint + if(getindex(cvars2d,'gust')>0) & + call stplimg(dval(1),sval(1),sges,pbc(1,6),nstep) + if(pjcalc)pj(6,1)=pbc(1,6)+pbc(ipenloc,6) + +! penalties for vis constraint + if(getindex(cvars2d,'vis')>0) & + call stplimv(dval(1),sval(1),sges,pbc(1,7),nstep) + if(pjcalc)pj(7,1)=pbc(1,7)+pbc(ipenloc,7) + +! penalties for pblh constraint + if(getindex(cvars2d,'pblh')>0) & + call stplimp(dval(1),sval(1),sges,pbc(1,8),nstep) + if(pjcalc)pj(8,1)=pbc(1,8)+pbc(ipenloc,8) + +! penalties for wspd10m constraint + if(getindex(cvars2d,'wspd10m')>0) & + call stplimw10m(dval(1),sval(1),sges,pbc(1,9),nstep) + if(pjcalc)pj(9,1)=pbc(1,9)+pbc(ipenloc,9) + +! penalties for howv constraint + if(getindex(cvars2d,'howv')>0) & + call stplimhowv(dval(1),sval(1),sges,pbc(1,10),nstep) + if(pjcalc)pj(10,1)=pbc(1,10)+pbc(ipenloc,10) + +! penalties for lcbas constraint + if(getindex(cvars2d,'lcbas')>0) & + call stpliml(dval(1),sval(1),sges,pbc(1,11),nstep) + if(pjcalc)pj(11,1)=pbc(1,11)+pbc(ipenloc,11) + +! penalties for cldch constraint + if(getindex(cvars2d,'cldch')>0) & + call stplimcldch(dval(1),sval(1),sges,pbc(1,12),nstep) + if(pjcalc)pj(12,1)=pbc(1,12)+pbc(ipenloc,12) + if (ljclimqc) then - if (getindex(cvars3d,'ql')>0) then + if (getindex(cvars3d,'ql')>0) then if(.not.ljc4tlevs) then call stplimqc(dval(ibin_anl),sval(ibin_anl),sges,pbc(1,13),nstep,ntguessig,'ql') - if(ii == 1) pj(13,1)=pbc(1,13)+pbc(ipenloc,13) + if(pjcalc) pj(13,1)=pbc(1,13)+pbc(ipenloc,13) else do ibin=1,nobs_bins if (nobs_bins /= nfldsig) then @@ -477,17 +547,17 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & pbc(j,13) = pbc(j,13)+pbcql(j,ibin) end do end do - if(ii == 1)then + if(pjcalc)then do ibin=1,nobs_bins pj(13,ibin)=pj(13,ibin)+pbcql(1,ibin)+pbcql(ipenloc,ibin) end do end if end if - end if - if (getindex(cvars3d,'qi')>0) then + end if + if (getindex(cvars3d,'qi')>0) then if(.not.ljc4tlevs) then call stplimqc(dval(ibin_anl),sval(ibin_anl),sges,pbc(1,14),nstep,ntguessig,'qi') - if(ii == 1) pj(14,1)=pbc(1,14)+pbc(ipenloc,14) + if(pjcalc) pj(14,1)=pbc(1,14)+pbc(ipenloc,14) else do ibin=1,nobs_bins if (nobs_bins /= nfldsig) then @@ -503,17 +573,17 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & pbc(j,14) = pbc(j,14)+pbcqi(j,ibin) end do end do - if(ii == 1)then + if(pjcalc)then do ibin=1,nobs_bins pj(14,ibin)=pj(14,ibin)+pbcqi(1,ibin)+pbcqi(ipenloc,ibin) end do end if end if - end if - if (getindex(cvars3d,'qr')>0) then + end if + if (getindex(cvars3d,'qr')>0) then if(.not.ljc4tlevs) then call stplimqc(dval(ibin_anl),sval(ibin_anl),sges,pbc(1,15),nstep,ntguessig,'qr') - if(ii == 1) pj(15,1)=pbc(1,15)+pbc(ipenloc,15) + if(pjcalc) pj(15,1)=pbc(1,15)+pbc(ipenloc,15) else do ibin=1,nobs_bins if (nobs_bins /= nfldsig) then @@ -529,17 +599,17 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & pbc(j,15) = pbc(j,15)+pbcqr(j,ibin) end do end do - if(ii == 1)then + if(pjcalc)then do ibin=1,nobs_bins pj(15,ibin)=pj(15,ibin)+pbcqr(1,ibin)+pbcqr(ipenloc,ibin) end do end if end if - end if - if (getindex(cvars3d,'qs')>0) then + end if + if (getindex(cvars3d,'qs')>0) then if(.not.ljc4tlevs) then call stplimqc(dval(ibin_anl),sval(ibin_anl),sges,pbc(1,16),nstep,ntguessig,'qs') - if(ii == 1) pj(16,1)=pbc(1,16)+pbc(ipenloc,16) + if(pjcalc) pj(16,1)=pbc(1,16)+pbc(ipenloc,16) else do ibin=1,nobs_bins if (nobs_bins /= nfldsig) then @@ -555,17 +625,17 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & pbc(j,16) = pbc(j,16)+pbcqs(j,ibin) end do end do - if(ii == 1)then + if(pjcalc)then do ibin=1,nobs_bins pj(16,ibin)=pj(16,ibin)+pbcqs(1,ibin)+pbcqs(ipenloc,ibin) end do end if end if - end if - if (getindex(cvars3d,'qg')>0) then + end if + if (getindex(cvars3d,'qg')>0) then if(.not.ljc4tlevs) then call stplimqc(dval(ibin_anl),sval(ibin_anl),sges,pbc(1,17),nstep,ntguessig,'qg') - if(ii == 1) pj(17,1)=pbc(1,17)+pbc(ipenloc,17) + if(pjcalc) pj(17,1)=pbc(1,17)+pbc(ipenloc,17) else do ibin=1,nobs_bins if (nobs_bins /= nfldsig) then @@ -581,78 +651,14 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & pbc(j,17) = pbc(j,17)+pbcqg(j,ibin) end do end do - if(ii == 1)then + if(pjcalc)then do ibin=1,nobs_bins pj(17,ibin)=pj(17,ibin)+pbcqg(1,ibin)+pbcqg(ipenloc,ibin) end do end if end if - end if + end if end if ! ljclimqc - if(.not.ljc4tlevs) then - call stplimq(dval(ibin_anl),sval(ibin_anl),sges,pbc(1,4),pbc(1,5),nstep,ntguessig) - if(ii == 1)then - pj(4,1)=pbc(1,4)+pbc(ipenloc,4) - pj(5,1)=pbc(1,5)+pbc(ipenloc,5) - end if - else - do ibin=1,nobs_bins - if (nobs_bins /= nfldsig) then - it=ntguessig - else - it=ibin - end if - call stplimq(dval(ibin),sval(ibin),sges,pbcqmin(1,ibin),pbcqmax(1,ibin),nstep,it) - end do - pbc(:,4)=zero_quad - pbc(:,5)=zero_quad - do ibin=1,nobs_bins - do j=1,nstep - pbc(j,4) = pbc(j,4)+pbcqmin(j,ibin) - pbc(j,5) = pbc(j,5)+pbcqmax(j,ibin) - end do - end do - if(ii == 1)then - do ibin=1,nobs_bins - pj(4,ibin)=pj(4,ibin)+pbcqmin(1,ibin)+pbcqmin(ipenloc,ibin) - pj(5,ibin)=pj(5,ibin)+pbcqmax(1,ibin)+pbcqmax(ipenloc,ibin) - end do - end if - end if -! penalties for gust constraint - if(getindex(cvars2d,'gust')>0) & - call stplimg(dval(1),sval(1),sges,pbc(1,6),nstep) - if(ii == 1)pj(6,1)=pbc(1,6)+pbc(ipenloc,6) - -! penalties for vis constraint - if(getindex(cvars2d,'vis')>0) & - call stplimv(dval(1),sval(1),sges,pbc(1,7),nstep) - if(ii == 1)pj(7,1)=pbc(1,7)+pbc(ipenloc,7) - -! penalties for pblh constraint - if(getindex(cvars2d,'pblh')>0) & - call stplimp(dval(1),sval(1),sges,pbc(1,8),nstep) - if(ii == 1)pj(8,1)=pbc(1,8)+pbc(ipenloc,8) - -! penalties for wspd10m constraint - if(getindex(cvars2d,'wspd10m')>0) & - call stplimw10m(dval(1),sval(1),sges,pbc(1,9),nstep) - if(ii == 1)pj(9,1)=pbc(1,9)+pbc(ipenloc,9) - -! penalties for howv constraint - if(getindex(cvars2d,'howv')>0) & - call stplimhowv(dval(1),sval(1),sges,pbc(1,10),nstep) - if(ii == 1)pj(10,1)=pbc(1,10)+pbc(ipenloc,10) - -! penalties for lcbas constraint - if(getindex(cvars2d,'lcbas')>0) & - call stpliml(dval(1),sval(1),sges,pbc(1,11),nstep) - if(ii == 1)pj(11,1)=pbc(1,11)+pbc(ipenloc,11) - -! penalties for cldch constraint - if(getindex(cvars2d,'cldch')>0) & - call stplimcldch(dval(1),sval(1),sges,pbc(1,12),nstep) - if(ii == 1)pj(12,1)=pbc(1,12)+pbc(ipenloc,12) end if ! penalties for Jo @@ -667,7 +673,7 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & end do end do enddo - if(ii == 1)then + if(pjcalc)then do ibin=1,size(pbcjoi,3) do j=1,size(pbcjoi,2) pj(n0+j,ibin)=pj(n0+j,ibin)+pbcjoi(ipenloc,j,ibin)+pbcjoi(1,j,ibin) @@ -866,7 +872,6 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, & exit stepsize end if end do stepsize - kprt=3 if(kprt >= 2 .and. iter == 0)then call mpl_allreduce(ipen,nobs_bins,pj) if(mype == 0)call prnt_j(pj,n0,ipen,kprt) diff --git a/src/gsi/stpjcmod.f90 b/src/gsi/stpjcmod.f90 index 289cc7672a..2c811912a0 100644 --- a/src/gsi/stpjcmod.f90 +++ b/src/gsi/stpjcmod.f90 @@ -77,7 +77,7 @@ subroutine stplimq(rval,sval,sges,outmin,outmax,nstep,itbin) ! !$$$ use gridmod, only: lat1,lon1,nsig,istart,wgtfactlats - use jfunc, only: factqmin,factqmax + use jfunc, only: factqmin,factqmax,superfact use guess_grids, only: ges_qsat use mpimod, only: mype implicit none @@ -125,8 +125,8 @@ subroutine stplimq(rval,sval,sges,outmin,outmax,nstep,itbin) /(ges_qsat(i,j,k,itbin)*ges_qsat(i,j,k,itbin)) else if(qx > ges_qsat(i,j,k,itbin))then - outmax(kk)=outmax(kk)+(factqmax*wgtfactlats(ii))*(qx-ges_qsat(i,j,k,itbin))* & - (qx-ges_qsat(i,j,k,itbin))/(ges_qsat(i,j,k,itbin)*ges_qsat(i,j,k,itbin)) + outmax(kk)=outmax(kk)+(factqmax*wgtfactlats(ii))*(qx-superfact*ges_qsat(i,j,k,itbin))**2 & + /ges_qsat(i,j,k,itbin)**2 end if end if end do @@ -143,9 +143,9 @@ subroutine stplimq(rval,sval,sges,outmin,outmax,nstep,itbin) if(q < zero)then outmin(1)=outmin(1)+(factqmin*wgtfactlats(ii))*q*q/(ges_qsat(i,j,k,itbin)*ges_qsat(i,j,k,itbin)) else - if(q > ges_qsat(i,j,k,itbin))then - outmax(1)=outmax(1)+(factqmax*wgtfactlats(ii))*(q-ges_qsat(i,j,k,itbin))*& - (q-ges_qsat(i,j,k,itbin))/(ges_qsat(i,j,k,itbin)*ges_qsat(i,j,k,itbin)) + if(q > superfact*ges_qsat(i,j,k,itbin))then + outmax(1)=outmax(1)+(factqmax*wgtfactlats(ii))*(q-superfact*ges_qsat(i,j,k,itbin))**2 & + /ges_qsat(i,j,k,itbin)**2 end if end if end do @@ -242,7 +242,7 @@ subroutine stplimqc(rval,sval,sges,out,nstep,itbin,cldtype) endif if (mype==0) write(6,*)'stplimqc: factqc = ', factqc if (mype==0) write(6,*)'stplimqc: ier ier1 = ', ier, ier1 - if ( factqc==0 ) return + if ( factqc <= zero) return if ( ier/=0 .or. ier1/=0 ) return ! Loop over interior of subdomain diff --git a/src/gsi/stpw.f90 b/src/gsi/stpw.f90 index ee1569d731..423019680d 100644 --- a/src/gsi/stpw.f90 +++ b/src/gsi/stpw.f90 @@ -105,7 +105,7 @@ subroutine stpw(whead,rval,sval,out,sges,nstep) real(r_kind) valu,facu,valv,facv,w1,w2,w3,w4,w5,w6,w7,w8 real(r_kind) cg_t,t_pg,var_jb real(r_kind) uu,vv - real(r_kind),dimension(max(1,nstep))::pen + real(r_kind),dimension(max(1,nstep))::pen,penu,penv real(r_kind),pointer,dimension(:):: ru,rv,su,sv type(wNode), pointer :: wptr @@ -157,41 +157,39 @@ subroutine stpw(whead,rval,sval,out,sges,nstep) do kk=1,nstep uu=facu+sges(kk)*valu vv=facv+sges(kk)*valv - pen(kk)= (uu*uu+vv*vv)*wptr%err2 + penu(kk)= (uu*uu)*wptr%err2 + penv(kk)= (vv*vv)*wptr%err2 end do else - pen(1)= (wptr%ures*wptr%ures+wptr%vres*wptr%vres)*wptr%err2 + penu(1)= (wptr%ures*wptr%ures)*wptr%err2 + penv(1)= (wptr%vres*wptr%vres)*wptr%err2 end if ! Modify penalty term if nonlinear QC - if (vqc .and. nlnqc_iter .and. wptr%pg > tiny_r_kind .and. & + t_pg=zero + cg_t=zero + var_jb=zero + ibb=0 + ikk=0 + if (vqc .and. nlnqc_iter .and. wptr%pg > tiny_r_kind .and. & wptr%b > tiny_r_kind) then t_pg=wptr%pg*varqc_iter cg_t=cg_term/wptr%b - else - t_pg=zero - cg_t=zero - endif - + else if(njqc .and. wptr%jb > tiny_r_kind .and. wptr%jb <10.0_r_kind) then ! for Dr. Jim purser' non liear quality control - if(njqc .and. wptr%jb > tiny_r_kind .and. wptr%jb <10.0_r_kind) then var_jb =wptr%jb - else - var_jb=zero - endif + else if(nvqc .and. wptr%ib >0) then ! mix model VQC - if(nvqc .and. wptr%ib >0) then ibb=wptr%ib ikk=wptr%ik - else - ibb=0 - ikk=0 endif - - call vqc_stp(pen,nstep,t_pg,cg_t,var_jb,ibb,ikk) - + call vqc_stp(penu,nstep,t_pg,cg_t,var_jb,ibb,ikk) + call vqc_stp(penv,nstep,t_pg,cg_t,var_jb,ibb,ikk) + do kk=1,nstep + pen(kk)=penu(kk)+penv(kk) + end do out(1) = out(1)+pen(1)*wptr%raterr2 do kk=2,nstep out(kk) = out(kk)+(pen(kk)-pen(1))*wptr%raterr2 diff --git a/src/gsi/update_guess.f90 b/src/gsi/update_guess.f90 index 041ad50e1d..a66495bc27 100644 --- a/src/gsi/update_guess.f90 +++ b/src/gsi/update_guess.f90 @@ -113,7 +113,7 @@ subroutine update_guess(sval,sbias) use mpimod, only: mype use constants, only: zero,one,fv,max_varname_length,qmin,qcmin,tgmin,& r100,one_tenth,tiny_r_kind - use jfunc, only: iout_iter,bcoption,tsensible,clip_supersaturation + use jfunc, only: iout_iter,bcoption,tsensible,clip_supersaturation,superfact use gridmod, only: lat2,lon2,nsig,& regional,twodvar_regional,regional_ozone,& l_reg_update_hydro_delz @@ -269,7 +269,7 @@ subroutine update_guess(sval,sbias) call gsi_bundlegetpointer (gsi_metguess_bundle(it),guess(ic),ptr3dges,istatus) if (trim(guess(ic))=='q') then call upd_positive_fldr3_(ptr3dges,ptr3dinc, qmin) - if(clip_supersaturation) ptr3dges(:,:,:) = min(ptr3dges(:,:,:),ges_qsat(:,:,:,it)) + if(clip_supersaturation) ptr3dges(:,:,:) = min(ptr3dges(:,:,:),superfact*ges_qsat(:,:,:,it)) cycle endif if (trim(guess(ic))=='oz') then