subroutine read_guess(iyear,month,idd,mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    read_guess          read/compute various guess fields
!   prgmmr: parrish          org: np22                date: 1994-02-11
!
! abstract:  This routine performs various functions, all related in one
!            way or the other to the model guess.  Note that this routine
!            is for the global mode of the gsi.  Separate read_guess type
!            routines exist for the regional gsi.
!
!            Functions performed in this routine include the following
!              a) read atmospheric guess bias correction fields (optional)
!              b) read atmospheric guess fields (optionally update with bias correction)
!              c) read surface guess fields
!              d) compute average ozone at each level                           
!
!
! program history log:
!   1994-02-11  parrish
!   1998-04-03  weiyu yang
!   1999-08-24  derber, j., treadon, r., yang, w., first frozen mpp version
!   2004-06-16  treadon - update documentation
!   2004-08-02  treadon - add only to module use, add intent in/out
!   2004-12-15  treadon - remove variable mype from call load_geop_hgt
!   2005-01-27  treadon - replace inguesfc with rdgesfc
!   2005-02-23  wu - setup for qoption=2 and output stats of RH
!   2005-03-07  dee - support gmao model interface
!   2005-03-30  treadon - reformat code (cosmetic changes only)
!   2005-05-27  parrish - add call get_derivatives
!   2005-06-27  guo     - support of GMAO gridded fields
!   2005-07-28  guo     - added sfc component for GMAO grided fields
!   2005-09-29  kleist  - get derivatives of surface terrain for Jc term
!   2005-11-21  kleist  - expand calls to new genqsat and calctends
!   2005-11-21  derber  - modify qoption =1 to work consistently with =2
!   2005-11-29  derber - remove external iteration dependent calculations
!   2005-11-29  derber - add ozmz calculation                             
!   2005-12-09  guo     - remove GMAO derivative computation code.  Use
!                         unified NCEP compact_diff procedures.
!   2006-01-10  treadon - consolidate all read*guess calls into this routine
!   2006-02-02  treadon - load 3d pressure guess pressure and geopotential 
!                         height grids
!   2006-02-03  derber  - modify to increase reproducibility (ozmz)
!   2006-03-13  treadon - increase filename to 24 characters
!   2006-04-14  treadon - replace call read_gfsatm for bias with read_bias
!   2006-06-08  zhang,b - change "biascor>0" to "biascor>=0" for debug purpose
!   2006-07-28  derber  - include sensible temperature
!   2006-07-31  kleist  - use ges_ps instead of lnps
!   2006-09-28  treadon - add sfc_rough and load_fact10
!   2006-12-04  todling - merged NCEP & GMAO bias correction schemes
!   2007-03-15  todling - merged in da Silva/Cruz ESMF changes
!   2007-05-30  h.liu   - remove ozmz
!   2008-12-06  todling - some clean; generalized update biasa
!   2009-01-28  todling - remove original GMAO interface
!   2010-03-06  parrish - add option to read ozone from gfs
!   2010-03-15  parrish - add flag regional_ozone to turn on ozone in regional analysis
!   2010-03-31  treadon - replace read_gfsatm with read_gfs
!   2010-05-19  todling - pass year and month for read_gsf_chem; read_guess should never
!   2010-09-17  pagowski - add cmaq
!                         depend on obsmod - that's why idate not passed via common block
!   2010-10-18  hcHuang - add flag use_gfs_nemsio and link to read_nems and read_nems_chem
!   2010-10-21  r. yang - pass dd for read_gsf_chem
!   2012-02-21  wu      - remove regional_ozone--causes conflict with using gfs ozone
!   2012-12-21  s.liu   - add option to use reflectivity
!   2013-10-30  jung    - changed zero to qmin in sensible temp calc and re-compute sensible
!                         temperature after clipping supersaturation
!
!   input argument list:
!     mype     - mpi task id
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_kind,i_kind
  use jfunc, only: biascor,bcoption,clip_supersaturation
  use guess_grids, only:  nfldsig,ges_tv,ges_q,ges_tsen,ges_prsl,load_prsges,load_geop_hgt
  use m_gsiBiases,only : correct_bias,nbc
  use m_gsiBiases,only : bias_q,bias_tv,bias_cwmr,bias_oz,bias_ps,&
       bias_vor,bias_div,bias_tskin,bias_u,bias_v
  use gsi_io, only: read_bias
  use gridmod, only: lat2,lon2
  use gridmod, only: nsig
  use gridmod, only: wrf_mass_regional,wrf_nmm_regional,cmaq_regional,use_reflectivity,&
       twodvar_regional,netcdf,regional,nems_nmmb_regional,use_gfs_ozone
  use gridmod, only : raqms ! ajl
  use gridmod, only: use_gfs_nemsio
  use gfs_stratosphere, only: use_gfs_stratosphere

  use constants, only: zero,one,fv,qmin
  use ncepgfs_io, only: read_gfs,read_gfs_chem
  use ncepnems_io, only: read_nems,read_nems_chem

  implicit none

! Declare passed variables
  integer(i_kind),intent(in   ) :: iyear
  integer(i_kind),intent(in   ) :: month
  integer(i_kind),intent(in   ) :: idd
  integer(i_kind),intent(in   ) :: mype

! Declare local variables
  character(24) filename
  logical :: ice
  integer(i_kind) i,j,k,it,iret_bias
  integer(i_kind) :: iderivative

  real(r_kind) :: satval
  real(r_kind),dimension(lat2,lon2,nsig) :: satq
  real(r_kind),dimension(lat2,lon2):: work

!-----------------------------------------------------------------------------------
! Certain functions are only done once --> on the first outer iteration. 
! One-time functions include
!    a) read atmospheric guess fields (optionally add bias correction)
!    b) read surface guess fields
!

!    Handle regional interfaces
     if (regional)then
        if (wrf_nmm_regional) then
           if(netcdf) then
              call read_wrf_nmm_netcdf_guess(mype)
           else
              call read_wrf_nmm_binary_guess(mype)
           end if
        else if (wrf_mass_regional) then
           if(netcdf) then
              call read_wrf_mass_netcdf_guess(mype)
           else
              call read_wrf_mass_binary_guess(mype)
           end if
        else if(twodvar_regional) then
           call read_2d_guess(mype)
        else if (nems_nmmb_regional) then
           call read_nems_nmmb_guess(mype)
!          if(use_reflectivity) then
!          if(mype==0)write(6,*)'sliu in read_guess to read reflectivity',use_reflectivity
!          call read_nems_nmmb_guess_ref(mype)
!          end if
        else if (cmaq_regional) then
           call read_cmaq_guess(mype)
        end if
     

!    Otherwise, handle global interface (ie, NCEP GFS)
     else

!       If requested, read bias correction fields
        iret_bias=0
        if (biascor >= zero) then
           filename='biascor_in'
           call read_bias(filename,mype,nbc,work,bias_ps,bias_tskin,&
                bias_vor,bias_div,bias_u,bias_v,bias_tv,bias_q,&
                bias_cwmr,bias_oz,iret_bias)
        endif
        
!       Read atmospheric fields
#ifndef HAVE_ESMF
        if ( use_gfs_nemsio ) then
!!           WRITE(6,*)'WARNING :: you elect to read first guess field in NEMSIO format'
           call read_nems(mype)
           call read_nems_chem(iyear,month,idd)
        else
          if(raqms)then ! ajl
            call read_raqms(mype)
          else ! ajl
            call read_gfs(mype)
            call read_gfs_chem(iyear,month,idd)
          endif ! ajl
        end if
#endif

!    End of non-GMAO global interfaces
     endif
        
! If doing SBC, apply bias correction ...

  if(biascor>=zero .and. iret_bias==0 .and. bcoption==1 ) call correct_bias()

! Get sensible temperature (after bias correction's been applied)

  do it=1,nfldsig
      do k=1,nsig
         do j=1,lon2
            do i=1,lat2
!              ges_tsen(i,j,k,it)= ges_tv(i,j,k,it)/(one+fv*max(qmin,ges_q(i,j,k,it)))
               ges_tsen(i,j,k,it)= ges_tv(i,j,k,it)/(one+fv*max(zero,ges_q(i,j,k,it)))
            end do
         end do
      end do
  end do

! Load 3d subdomain pressure arrays from the guess fields
  call load_prsges

! recompute sensible temperature to remove supersaturation
  if ( clip_supersaturation ) then
    call tpause(mype,'pvoz')
    ice = .true.
    iderivative = 0
    do it=1,nfldsig
      call genqsat(satq,ges_tsen(1,1,1,it),ges_prsl(1,1,1,it),lat2,lon2, &
                 nsig,ice,iderivative)
      do k=1,nsig
         do j=1,lon2
            do i=1,lat2
               satval = min(ges_q(i,j,k,it),satq(i,j,k))
               satval = max(qmin,satval)
               ges_q(i,j,k,it) = satval
               ges_tsen(i,j,k,it)= ges_tv(i,j,k,it)/(one+fv*ges_q(i,j,k,it))
            end do
         end do
      end do
    end do
  endif   ! clip_supersaturation


! Compute 3d subdomain geopotential heights from the guess fields
  call load_geop_hgt

!  If this is a regional run and ozone is desired from the gfs model, bring it in here:
  if(regional.and.use_gfs_ozone.and..not.use_gfs_stratosphere) &
              call read_gfs_ozone_for_regional
  
  return
end subroutine read_guess
