Skip to content

Commit

Permalink
Merge pull request #131 from KristenBathmann-NOAA/release/gfsda.v16.1.0
Browse files Browse the repository at this point in the history 
Commercial RO changes for Release/gfsda.v16.1.0
  • Loading branch information
@CatherineThomas-NOAA
CatherineThomas-NOAA authored Mar 16, 2021
2 parents 6974b87 + fcf1412 commit 1e49111
Show file tree
Hide file tree
Showing 4 changed files with 89 additions and 80 deletions.
2 changes: 1 addition & 1 deletion scripts/exglobal_atmos_analysis.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -749,7 +749,7 @@ cat > gsiparm.anl << EOF
iguess=-1,
tzr_qc=$TZR_QC,
oneobtest=.false.,retrieval=.false.,l_foto=.false.,
use_pbl=.false.,use_compress=.true.,nsig_ext=12,gpstop=50.,
use_pbl=.false.,use_compress=.true.,nsig_ext=12,gpstop=50.,commgpstop=50.,commgpserrinf=1.0,
use_gfs_nemsio=${use_gfs_nemsio},use_gfs_ncio=${use_gfs_ncio},sfcnst_comb=.true.,
use_readin_anl_sfcmask=${USE_READIN_ANL_SFCMASK},
lrun_subdirs=$lrun_subdirs,
Expand Down
6 changes: 4 additions & 2 deletions src/gsi/gsimod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -123,7 +123,7 @@ module gsimod
use_gfs_nemsio,sfcnst_comb,use_readin_anl_sfcmask,use_sp_eqspace,final_grid_vars,&
jcap_gfs,nlat_gfs,nlon_gfs,jcap_cut,wrf_mass_hybridcord,use_gfs_ncio,write_fv3_incr,&
use_fv3_aero
use guess_grids, only: ifact10,sfcmod_gfs,sfcmod_mm5,use_compress,nsig_ext,gpstop
use guess_grids, only: ifact10,sfcmod_gfs,sfcmod_mm5,use_compress,nsig_ext,gpstop,commgpstop,commgpserrinf
use gsi_io, only: init_io,lendian_in,verbose,print_obs_para
use regional_io_mod, only: regional_io_class
use wrf_params_mod, only: update_pint, preserve_restart_date
Expand Down Expand Up @@ -567,6 +567,8 @@ module gsimod
! use_compress - option to turn on the use of compressibility factors in geopotential heights
! nsig_ext - number of layers above the model top which are necessary to compute the bending angle for gpsro
! gpstop - maximum height for gpsro data assimilation. Reject anything above this height.
! commgpstop -Reject commercial ro above this height. Logic in setupbend assumes commgpstop <= gpstop.
! commgpserrinf - optional error inflation factor for commercial gpsro data
! use_gfs_nemsio - option to use nemsio to read global model NEMS/GFS first guess
! use_gfs_ncio - option to use netCDF to read global model FV3-GFS first guess
! use_fv3_aero - option to use FV3-Chem vs NGAC for global aerosol analysis
Expand Down Expand Up @@ -659,7 +661,7 @@ module gsimod
diag_rad,diag_pcp,diag_conv,diag_ozone,diag_aero,diag_co,diag_light,diag_radardbz,iguess, &
write_diag,reduce_diag, &
oneobtest,sfcmodel,dtbduv_on,ifact10,l_foto,offtime_data,&
use_pbl,use_compress,nsig_ext,gpstop,&
use_pbl,use_compress,nsig_ext,gpstop,commgpstop, commgpserrinf, &
perturb_obs,perturb_fact,oberror_tune,preserve_restart_date, &
crtm_coeffs_path,berror_stats, &
newpc4pred,adp_anglebc,angord,passive_bc,use_edges,emiss_bc,upd_pred,reset_bad_radbc,&
Expand Down
5 changes: 3 additions & 2 deletions src/gsi/guess_grids.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,7 +143,7 @@ module guess_grids
public :: geop_hgti,ges_lnprsi,ges_lnprsl,geop_hgtl,pbl_height,ges_geopi
public :: wgt_lcbas
public :: ges_qsat
public :: use_compress,nsig_ext,gpstop
public :: use_compress,nsig_ext,gpstop,commgpstop,commgpserrinf
public :: ges_tsen1,ges_q1
public :: ntguesaer,ifileaer,nfldaer,hrdifaer ! variables for external aerosol files

Expand Down Expand Up @@ -229,7 +229,8 @@ module guess_grids

real(r_kind):: gpstop=30.0_r_kind ! maximum gpsro height used in km
! geometric height for ref, impact height for bnd

real(r_kind):: commgpstop=30.0_r_kind
real(r_kind):: commgpserrinf=1.0_r_kind ! error inflation factor for commercial gnssro
real(r_kind):: ges_psfcavg ! average guess surface pressure
real(r_kind),allocatable,dimension(:):: ges_prslavg ! average guess pressure profile

Expand Down
156 changes: 81 additions & 75 deletions src/gsi/setupbend.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,7 @@ subroutine setupbend(obsLL,odiagLL, &

use gsi_4dvar, only: nobs_bins,hr_obsbin
use guess_grids, only: ges_lnprsi,hrdifsig,geop_hgti,nfldsig
use guess_grids, only: nsig_ext,gpstop
use guess_grids, only: nsig_ext,gpstop,commgpstop,commgpserrinf
use gridmod, only: nsig
use gridmod, only: get_ij,latlon11
use constants, only: fv,n_a,n_b,n_c,deg2rad,tiny_r_kind,r0_01,r18,r61,r63,r10000
Expand Down Expand Up @@ -257,6 +257,7 @@ subroutine setupbend(obsLL,odiagLL, &
real(r_kind),allocatable,dimension(:,:,:,:) :: ges_q

type(obsLList),pointer,dimension(:):: gpshead
logical:: commdat
gpshead => obsLL(:)

save_jacobian = conv_diagsave .and. jiter==jiterstart .and. lobsdiag_forenkf
Expand Down Expand Up @@ -679,6 +680,11 @@ subroutine setupbend(obsLL,odiagLL, &

repe_gps=exp(repe_gps) ! one/modified error in (rad-1*1E3)
repe_gps= r1em3*(one/abs(repe_gps)) ! modified error in rad
commdat=.false.
if (data(isatid,i)>=265 .and. data(isatid,i)<=269) then
commdat=.true.
repe_gps=commgpserrinf*repe_gps ! Inflate error for commercial data
endif
ratio_errors(i) = data(ier,i)/abs(repe_gps)

error(i)=one/data(ier,i) ! one/original error
Expand Down Expand Up @@ -768,85 +774,85 @@ subroutine setupbend(obsLL,odiagLL, &
rdiagbuf( 5,i) = (data(igps,i)-dbend)/data(igps,i) ! incremental bending angle (x100 %)

data(igps,i)=data(igps,i)-dbend !innovation vector

if (alt <= gpstop) then ! go into qc checks
cgrossuse=cgross(ikx)
cermaxuse=cermax(ikx)
cerminuse=cermin(ikx)
if (alt > five) then
cgrossuse=cgrossuse*r400
cermaxuse=cermaxuse*r400
cerminuse=cerminuse*r100
endif
! Gross error check
obserror = one/max(ratio_errors(i)*data(ier,i),tiny_r_kind)
obserrlm = max(cerminuse,min(cermaxuse,obserror))
residual = abs(data(igps,i))
ratio = residual/obserrlm

if (ratio > cgrossuse) then
if (luse(i)) then
awork(4) = awork(4)+one
endif
qcfail_gross(i)=one
data(ier,i) = zero
ratio_errors(i) = zero
muse(i)=.false.
else
! Statistics QC check if obs passed gross error check
cutoff=zero
if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then
cutoff1=(-4.725_r_kind+0.045_r_kind*alt+0.005_r_kind*alt**2)*one/two
else
cutoff1=(-4.725_r_kind+0.045_r_kind*alt+0.005_r_kind*alt**2)*two/three
end if
cutoff2=1.5_r_kind+one*cos(data(ilate,i)*deg2rad)
if(trefges<=r240) then
cutoff3=two
else
cutoff3=0.005_r_kind*trefges**2-2.3_r_kind*trefges+266_r_kind
endif
if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then
cutoff3=cutoff3*one/two
else
cutoff3=cutoff3*two/three
end if
if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then
cutoff4=(four+eight*cos(data(ilate,i)*deg2rad))*one/two
else
cutoff4=(four+eight*cos(data(ilate,i)*deg2rad))*two/three
end if
cutoff12=((36_r_kind-alt)/two)*cutoff2+&
((alt-34_r_kind)/two)*cutoff1
cutoff23=((eleven-alt)/two)*cutoff3+&
((alt-nine)/two)*cutoff2
cutoff34=((six-alt)/two)*cutoff4+&
((alt-four)/two)*cutoff3
if(alt>36_r_kind) cutoff=cutoff1
if((alt<=36_r_kind).and.(alt>34_r_kind)) cutoff=cutoff12
if((alt<=34_r_kind).and.(alt>eleven)) cutoff=cutoff2
if((alt<=eleven).and.(alt>nine)) cutoff=cutoff23
if((alt<=nine).and.(alt>six)) cutoff=cutoff3
if((alt<=six).and.(alt>four)) cutoff=cutoff34
if(alt<=four) cutoff=cutoff4

if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then
cutoff=two*cutoff*r0_01
else
cutoff=three*cutoff*r0_01
end if

if(abs(rdiagbuf(5,i)) > cutoff) then
qcfail(i)=.true.
if ((alt <= commgpstop) .or. (.not.commdat)) then
cgrossuse=cgross(ikx)
cermaxuse=cermax(ikx)
cerminuse=cermin(ikx)
if (alt > five) then
cgrossuse=cgrossuse*r400
cermaxuse=cermaxuse*r400
cerminuse=cerminuse*r100
endif
! Gross error check
obserror = one/max(ratio_errors(i)*data(ier,i),tiny_r_kind)
obserrlm = max(cerminuse,min(cermaxuse,obserror))
residual = abs(data(igps,i))
ratio = residual/obserrlm

if (ratio > cgrossuse) then
if (luse(i)) then
awork(4) = awork(4)+one
endif
qcfail_gross(i)=one
data(ier,i) = zero
ratio_errors(i) = zero
muse(i) = .false.
end if
end if ! gross qc check
muse(i)=.false.
else
! Statistics QC check if obs passed gross error check
cutoff=zero
if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then
cutoff1=(-4.725_r_kind+0.045_r_kind*alt+0.005_r_kind*alt**2)*one/two
else
cutoff1=(-4.725_r_kind+0.045_r_kind*alt+0.005_r_kind*alt**2)*two/three
end if
cutoff2=1.5_r_kind+one*cos(data(ilate,i)*deg2rad)
if(trefges<=r240) then
cutoff3=two
else
cutoff3=0.005_r_kind*trefges**2-2.3_r_kind*trefges+266_r_kind
endif
if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then
cutoff3=cutoff3*one/two
else
cutoff3=cutoff3*two/three
end if
if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then
cutoff4=(four+eight*cos(data(ilate,i)*deg2rad))*one/two
else
cutoff4=(four+eight*cos(data(ilate,i)*deg2rad))*two/three
end if
cutoff12=((36_r_kind-alt)/two)*cutoff2+&
((alt-34_r_kind)/two)*cutoff1
cutoff23=((eleven-alt)/two)*cutoff3+&
((alt-nine)/two)*cutoff2
cutoff34=((six-alt)/two)*cutoff4+&
((alt-four)/two)*cutoff3
if(alt>36_r_kind) cutoff=cutoff1
if((alt<=36_r_kind).and.(alt>34_r_kind)) cutoff=cutoff12
if((alt<=34_r_kind).and.(alt>eleven)) cutoff=cutoff2
if((alt<=eleven).and.(alt>nine)) cutoff=cutoff23
if((alt<=nine).and.(alt>six)) cutoff=cutoff3
if((alt<=six).and.(alt>four)) cutoff=cutoff34
if(alt<=four) cutoff=cutoff4

if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then
cutoff=two*cutoff*r0_01
else
cutoff=three*cutoff*r0_01
end if

if(abs(rdiagbuf(5,i)) > cutoff) then
qcfail(i)=.true.
data(ier,i) = zero
ratio_errors(i) = zero
muse(i) = .false.
end if
end if !gross qc check
end if ! commdat < commgpstop
end if ! qc checks (only below 50km)

! Remove obs above 50 km
if(alt > gpstop) then
if((alt > gpstop) .or. (commdat .and. (alt > commgpstop))) then
data(ier,i) = zero
ratio_errors(i) = zero
qcfail_high(i)=one
Expand Down

0 comments on commit 1e49111

Please sign in to comment.