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
!   2020-03-09  pondeca - deallocate arrays used in similarity theory-based 10-wind
!                         adjustment for twodvar_regional
!
!   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_similarity_2dvar
  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 aux2dvarflds, only: destroy_aux2dvarflds

  use m_prad, only: prad_updatePredx    ! was -- prad_bias()
  use m_obsdiags, only: obsdiags_write
  use gsi_io,only: verbose
  use m_berror_stats,only: inquire_berror

  implicit none

! Declare passed variables

! Declare local variables
  logical laltmin

  integer(i_kind) jiterlast,lunix,lunit
  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

  lunit=0
  lunix=0
  call isetprm('MXMSGL',400000)
  call isetprm('MAXSS',250000)
!  Initialize bufr read on all processors (so that exitbufr works, lunit and
!  lunix are dummys)
  call openbf(lunit,'FIRST',lunix)
! Read observations and scatter
  call observer_set
!  release bufr memory
  call exitbufr

! 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
     lunit=22
     call inquire_berror(lunit,mype)
     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

  if (twodvar_regional) then
     if (use_similarity_2dvar .and. jiter == miter) call destroy_aux2dvarflds
  endif

! Write updated bias correction coefficients
  if (.not.twodvar_regional) then
     if (l4dvar) then
        call radinfo_write(0)
        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
           call radinfo_write(0)
           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