subroutine glbsoi
!$$$ subprogram documentation block
! . . . .
! subprogram: glbsoi driver for gridpoint statistical
! interpolation
! prgmmr: parrish org: np22 date: 1990-10-06
!
! abstract: executes gridpoint spectral interpolation analysis.
!
! program history log:
! 1990-10-06 parrish
! 1994-02-11 parrish
! 1998-05-15 weiyu yang mpp version
! 1999-08-24 derber, j., treadon, r., yang, w., first frozen mpp version
! 1993-12-22 kleist,d., treadon, r.,derber, j. modules, updates and comments
! 2004-06-21 treadon - update documentation
! 2004-07-08 treadon - fix vertical indexing bug in call set_ozone_var
! 2004-07-24 treadon - add only to module use, add intent in/out
! 2004-11-22 parrish - add code to handle regional netcdf i/o
! 2004-12-13 treadon - limit runtime output from CRTM & IRSSE initilizations
! 2004-12-22 treadon - add optional code to compute/write out innovation
! information following all outer loops
! 2005-01-22 parrish - add balmod, compact_diffs
! 2005-02-23 wu - setup norm rh
! 2005-03-01 parrish - add parts of regional anisotropic background option
! (not complete yet)
! 2005-04-14 yanqiu zhu - support for observation sensitivity study
! 2005-05-24 pondeca - add 2dvar only surface analysis option
! 2005-05-27 kleist - add call to destroy grids from patch interpolation
! 2005-06-08 treadon - add switch_on_derivatives for create/destroy_ges_bias_grids
! 2005-06-13 li/treadon - move destroy_sfc_grids after write_all
! 2005-07-12 kleist - add Jc term
! 2005-08-16 guo - add gmao surface interface
! 2005-09-28 parrish - add logic to only allow NCEP global or WRF_NMM to use jc routines
! 2005-09-29 derber - simplify observation file handling
! 2005-11-21 kleist - use tendency module, add call to get moisture diagnostics
! 2005-11-28 derber - move read_obs to glbsoi
! 2005-11-29 derber - move read guess and background error calculations outside
! external iteration
! 2006-01-09 derber - absorb set_nrh_var into compute_derived
! 2006-01-10 treadon - move read*files into gesinfo, consolidate read*guess calls
! 2006-01-12 treadon - replace pCRTM with CRTM
! 2006-02-03 derber - modify for new obs control and obs count
! 2006-02-24 derber - modify to take advantage of convinfo module
! 2006-04-12 treadon - remove mpimod (not used)
! 2006-04-21 kleist - remove call to setupjc, no longer exists
! 2006-07-28 derber - remove creation of obs inner loop data file name
! 2006-08-15 parrish - add call to create_vtrans (get vertical modes if nvmodes_keep > 0)
! 2006-12-04 todling - split bias and guess init/final (rename ges_bias routines)
! 2006-12-15 todling - no need to protect destroy call to tendvars
! 2007-03-15 todling - merged in da Silva/Cruz ESMF changes
! 2007-03-15 su - to delocate the converr arrays
! 2007-05-08 kleist - add preliminary verions of mp_compact_diffs module
! 2007-06-08 kleist/treadon - add prefix (task id or path) to obs_setup
! 2007-06-21 rancic - add pbl code
! 2007-06-27 tremolet- observation sensitivity
! 2007-06-29 jung - update CRTM interface
! 2007-07-05 todling - skip calculating o-a when running in 4d-var mode
! 2007-08-27 tremolet- changed outer loop control for 4dvar
! 2007-09-30 todling - add timer
! 2007-10-25 todling - obsdiag files now written by observer
! 2007-11-12 todling - write sat bias moved from write_all here (ESMF-interface support)
! 2008-01-04 tremolet- outer loop for sensitivity computations
! 2008-11-03 sato - enable use of global anisotropic mode
! 2008-12-02 todling - remove references to pcgsoi_tl and old obs_sen
! 2009-01-28 todling - move write_all to pcgsoi (for consistency w/ 4dvar branch of code)
! - use observer to avoid redundant code
! 2009-08-19 guo - changed setuprhsall() interface for multi-pass observer.
! 2009-08-31 parrish - add call to fmg_initialize_e when tlnmc_type=4. Initializes
! alternative regional tangent linear normal mode constraint which
! allows for variation of coriolis parameter and map factor.
! 2009-09-12 parrish - add call to hybrid_ensemble_setup. if l_hyb_ens=.true., then
! subroutine hybrid_ensemble_setup is called, which creates
! everything needed for a hybrid ensemble 3dvar analysis.
! 2009-09-14 guo - move compact_diff related statements to observer_init() in observer.F90
! 2010-02-20 parrish - move hybrid_ensemble_setup to beginning of code and read
! in ensemble perturbations where hybrid_ensemble_setup was previously located.
! 2010-04-27 zhu - add call to pcinfo. when newpc4pred=.true., pcinfo is called
! to calculate additional preconditioner; add newpc4pred in
! radinfo_write's interface
! 2010-05-12 zhu - add option passive_bc for radiance bias correction for monitored channels
! 2010-10-01 el akkraoui/todling - add Bi-CG as optional minimization scheme
! 2011-04-07 todling - newpc4pred now in radinfo
! 2011-08-01 lueken - replaced F90 with f90 (no machine logic)
! 2011-09-05 todling - repositioned initialization of ensemble to enable sqrt-ens feature
! 2012-07-12 todling - read yhatsave as well as xhatsave in 4dvar mode; knob for nested resolution
! 2012-09-14 Syed RH Rizvi, NCAR/NESL/MMM/DAS - implemented obs adjoint test
! 2013-05-19 zhu - add aircraft temperature bias correction
! 2013-07-02 parrish - remove references to init_strongvars_1, init_strongvars_2
! 2014-02-03 todling - move cost function create/destroy from observer to this routine;
! reposition load of ens due to init of sqrt(ens) dims dependences
! 2014-02-05 todling - update interface to prebal
! 2014-06-19 carley/zhu - Modify for R_option: optional variable correlation length twodvar_regional
! lcbas analysis variable
! 2015-10-01 guo - omb at full res; support via obsdiags
! 2015-12-08 el akkraoui - Y. Zhu sat-bias-corr now works with BiCG option
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use mpimod, only: mype,npe
use adjtest_obs, only: adtest_obs
use jfunc, only: miter,jiter,jiterstart,jiterend,iguess,&
write_guess_solution,R_option,&
xhatsave,yhatsave,create_jfunc,destroy_jfunc
use anberror, only: anisotropic, &
create_anberror_vars_reg,destroy_anberror_vars_reg,&
create_anberror_vars,destroy_anberror_vars
use anisofilter, only: anprewgt_reg
use anisofilter_glb, only: anprewgt
use berror, only: create_berror_vars_reg,create_berror_vars,&
set_predictors_var,destroy_berror_vars_reg,destroy_berror_vars,&
bkgv_flowdep,pcinfo,fut2ps,cwcoveqqcov
use balmod, only: create_balance_vars_reg,create_balance_vars, &
destroy_balance_vars_reg,destroy_balance_vars,prebal,prebal_reg
use compact_diffs, only: create_cdiff_coefs,inisph
use gridmod, only: regional,twodvar_regional
use guess_grids, only: nfldsig
use obsmod, only: write_diag,perturb_obs,ditype,iadate
use qcmod,only: njqc
use turblmod, only: create_turblvars,destroy_turblvars
use obs_sensitivity, only: lobsensfc, iobsconv, lsensrecompute, &
init_fc_sens, save_fc_sens, lobsensincr, lobsensjb
use smooth_polcarf, only: norsp,destroy_smooth_polcas
use jcmod, only: ljcdfi
use gsi_4dvar, only: l4dvar, lsqrtb, lbicg, lanczosave, lnested_loops, ladtest_obs
use pcgsoimod, only: pcgsoi
use control_vectors, only: dot_product
use radinfo, only: radinfo_write,passive_bc,newpc4pred
use pcpinfo, only: pcpinfo_write
use converr, only: converr_destroy
use converr_ps, only: converr_ps_destroy
use converr_q, only: converr_q_destroy
use converr_t, only: converr_t_destroy
use converr_uv, only: converr_uv_destroy
use converr_pw, only: converr_pw_destroy
use convb_ps, only: convb_ps_destroy
use convb_q, only: convb_q_destroy
use convb_t, only: convb_t_destroy
use convb_uv, only: convb_uv_destroy
use zrnmi_mod, only: zrnmi_initialize
use observermod, only: observer_init,observer_set,observer_finalize,ndata
use timermod, only: timer_ini, timer_fnl
use hybrid_ensemble_parameters, only: l_hyb_ens,destroy_hybens_localization_parameters
use hybrid_ensemble_isotropic, only: create_ensemble,load_ensemble,destroy_ensemble, &
hybens_localization_setup,hybens_grid_setup
use gfs_stratosphere, only: destroy_nmmb_vcoords,use_gfs_stratosphere
use aircraftinfo, only: aircraftinfo_write,aircraft_t_bc_pof,aircraft_t_bc,mype_airobst
use rapidrefresh_cldsurf_mod, only: i_gsdcldanal_type
use m_prad, only: prad_updatePredx ! was -- prad_bias()
use m_obsdiags, only: obsdiags_write
use gsi_io,only: verbose
implicit none
! Declare passed variables
! Declare local variables
logical laltmin
integer(i_kind) jiterlast
real(r_kind) :: zgg,zxy
character(len=12) :: clfile
logical print_verbose
print_verbose=.false.
if(verbose)print_verbose=.true.
!*******************************************************************************************
!
! Initialize timer for this procedure
call timer_ini('glbsoi')
if(mype==0) write(6,*) 'glbsoi: starting ...'
! If l_hyb_ens is true, then initialize machinery for hybrid ensemble 3dvar
if(l_hyb_ens) then
call hybens_grid_setup
call create_ensemble
end if
! Check for alternative minimizations
laltmin = lsqrtb.or.lbicg
! Initialize observer
call observer_init
! Check GSI options against available number of guess time levels
if (nfldsig == 1) then
if (bkgv_flowdep) then
bkgv_flowdep = .false.
if (mype==0) &
write(6,*)'GLBSOI: ***WARNING*** reset bkgv_flowdep=',bkgv_flowdep,&
', because only ',nfldsig,' guess time level available'
endif
endif
! Set cost function
call create_jfunc
! Read observations and scatter
call observer_set
! cloud analysis
if(i_gsdcldanal_type==6 .or. i_gsdcldanal_type==3) then
call gsdcloudanalysis(mype)
! Write output analysis files
call write_all(-1,mype)
call prt_guess('analysis')
! Finalize observer
call observer_finalize
! Finalize timer for this procedure
call timer_fnl('glbsoi')
return
endif
! Create/setup background error and background error balance
if (regional)then
call create_balance_vars_reg(mype)
if(anisotropic) then
call create_anberror_vars_reg(mype)
else
call create_berror_vars_reg
end if
call prebal_reg(cwcoveqqcov)
if (.not. R_option) then
if(anisotropic) then
call anprewgt_reg(mype)
else
call prewgt_reg(mype)
end if
end if
else
call create_balance_vars
if(anisotropic) then
call create_anberror_vars(mype)
else
call create_berror_vars
end if
call prebal(fut2ps,cwcoveqqcov)
! Load background error arrays used by recursive filters
if(anisotropic) then
call anprewgt(mype)
else
call prewgt(mype)
end if
end if
! If l_hyb_ens is true, then read in ensemble perturbations
if(l_hyb_ens) then
call load_ensemble
call hybens_localization_setup
end if
! Set error (variance) for predictors (only use guess)
call set_predictors_var
! Set errors and create variables for dynamical constraint
if (ljcdfi) call init_jcdfi
! Read output from previous min.
if (l4dvar.and.jiterstart>1) then
clfile='xhatsave.ZZZ'
write(clfile(10:12),'(I3.3)') jiterstart-1
call view_cv_ad(xhatsave,iadate,clfile,.not.lnested_loops)
zgg=dot_product(xhatsave,xhatsave)
if (mype==0) write(6,*)'Norm xhatsave=',sqrt(zgg)
if (.not.lsqrtb) then
clfile='yhatsave.ZZZ'
write(clfile(10:12),'(I3.3)') jiterstart-1
call view_cv_ad(yhatsave,iadate,clfile,.not.lnested_loops)
zgg=dot_product(yhatsave,yhatsave)
zxy=dot_product(xhatsave,yhatsave)
if (mype==0) then
write(6,*)'Norm yhatsave=',sqrt(zgg)
write(6,*)'Norm x,yhatsave=',zxy
endif
endif
endif
jiterlast=miter
if (lsensrecompute) jiterlast=jiterend
if (l4dvar) jiterlast=jiterstart
if (ladtest_obs) jiterlast=jiterstart
! Main outer analysis loop
do jiter=jiterstart,jiterlast
if (mype==0) write(6,'(a44,4i5)')'GLBSOI: jiter,jiterstart,jiterlast,jiterend=', &
jiter,jiterstart,jiterlast,jiterend
! Set up right hand side of analysis equation
call setuprhsall(ndata,mype,.true.,.true.)
! Estimate correlation length for lcbas if R_option==.true.
! For this to work we need to have run setuplcbas first to get the weights.
! Note that R_option is only applicable for lcbas and not for any
! another twodvar_regional variables. All other variables are handled
! as they would have been otherwise.
if (R_option .and. twodvar_regional .and. jiter==jiterstart) then
if(anisotropic) then
call anprewgt_reg(mype)
else
call prewgt_reg(mype)
end if
end if
! implement obs adjoint test and return
if( ladtest_obs) then
call adtest_obs
return
end if
if (jiter<=miter) then
! Set up right hand side of adjoint of analysis equation
if (lsensrecompute) lobsensfc=(jiter==jiterend)
if (lobsensfc.or.iobsconv>0) call init_fc_sens
! Call inner minimization loop
if (laltmin) then
if (lsqrtb) then
if (newpc4pred) then
if (mype==0) write(6,*)'GLBSOI: newpc4pred is not available for lsqrtb'
call stop2(334)
end if
if (mype==0) write(6,*)'GLBSOI: Using sqrt(B), jiter=',jiter
call sqrtmin
endif
if (lbicg) then
if (mype==0) write(6,*)'GLBSOI: Using bicg, jiter=',jiter
call pcinfo ! Set up additional preconditioning information
call bicg
endif
else
! Set up additional preconditioning information
call pcinfo
! Standard run
if (mype==0 .and. print_verbose) write(6,*)'GLBSOI: START pcgsoi jiter=',jiter
call pcgsoi
end if
! Save information for next minimization
if (lobsensfc) then
clfile='obsdiags.ZZZ'
write(clfile(10:12),'(I3.3)') 100+jiter
call obsdiags_write(clfile) ! replacing write_obsdiags()
if (lobsensincr .or. lobsensjb) then
clfile='xhatsave.ZZZ'
write(clfile(10:12),'(I3.3)') jiter
call view_cv(xhatsave,iadate,clfile,.not.lnested_loops)
endif
elseif (l4dvar.or.lanczosave) then
if (l4dvar) then
clfile='obsdiags.ZZZ'
write(clfile(10:12),'(I3.3)') jiter
call obsdiags_write(clfile) ! replacing write_obsdiags()
endif
clfile='xhatsave.ZZZ'
write(clfile(10:12),'(I3.3)') jiter
call view_cv(xhatsave,iadate,clfile,.not.lnested_loops)
zgg=dot_product(xhatsave,xhatsave)
if (mype==0) write(6,*)'Norm xhatsave=',sqrt(zgg)
if(.not.lsqrtb) then
clfile='yhatsave.ZZZ'
write(clfile(10:12),'(I3.3)') jiter
call view_cv(yhatsave,iadate,clfile,.not.lnested_loops)
zgg=dot_product(yhatsave,yhatsave)
zxy=dot_product(xhatsave,yhatsave)
if (mype==0) write(6,*)'Norm yhatsave=',sqrt(zgg)
if (mype==0) write(6,*)'Norm x,yhatsave=',zxy
endif
endif
! Save output of adjoint of analysis equation
if (lobsensfc.or.iobsconv>0) call save_fc_sens
endif
! End of outer iteration loop
end do
if (.not.l4dvar) then
jiter=jiterlast+1
! If requested, write obs-anl information to output files
if (write_diag(jiter)) then
call setuprhsall(ndata,mype,.true.,.true.)
if (.not. lsqrtb) call pcinfo
if (any(ditype=='rad') .and. passive_bc) call prad_updatePredx() ! was -- call prad_bias
end if
! Write xhat- and yhat-save for use as a guess for the solution
if (iguess==0 .or. iguess==1) call write_guess_solution(mype)
endif
! Deallocate arrays
if(perturb_obs) then
if(njqc) then
call converr_ps_destroy
call converr_q_destroy
call converr_t_destroy
call converr_uv_destroy
call converr_pw_destroy
call convb_ps_destroy
call convb_q_destroy
call convb_t_destroy
call convb_uv_destroy
else
call converr_destroy
endif
endif
if (regional) then
if(anisotropic) then
call destroy_anberror_vars_reg
else
call destroy_berror_vars_reg
end if
call destroy_balance_vars_reg
if(use_gfs_stratosphere) call destroy_nmmb_vcoords
else
if(anisotropic) then
call destroy_anberror_vars
else
call destroy_berror_vars
end if
call destroy_balance_vars
if (norsp > 0) call destroy_smooth_polcas
endif
if (l_hyb_ens) call destroy_hybens_localization_parameters
! Write updated bias correction coefficients
if (.not.twodvar_regional) then
if (l4dvar) then
if(mype == 0) call radinfo_write
if(mype == npe-1) call pcpinfo_write
if(mype==mype_airobst .and. (aircraft_t_bc_pof .or. aircraft_t_bc)) call aircraftinfo_write
else
if (jiter==miter+1 ) then
if(mype == 0) call radinfo_write
if(mype == npe-1) call pcpinfo_write
if(mype==mype_airobst .and. (aircraft_t_bc_pof .or. aircraft_t_bc)) call aircraftinfo_write
endif
endif
endif
! Finalize cost function
call destroy_jfunc
! Finalize observer
call observer_finalize
! When applicable, finalize ensemble
if(l_hyb_ens) then
call destroy_ensemble
endif
if(mype==0) write(6,*) 'glbsoi: complete'
! Finalize timer for this procedure
call timer_fnl('glbsoi')
! End of routine
end subroutine glbsoi