#ifdef AIX
  @PROCESS ALIAS_SIZE(805306368)
#endif
module datm_comp_mod

  ! !USES:

  use mct_mod
  use esmf
  use perf_mod
  use shr_const_mod
  use shr_sys_mod
  use shr_kind_mod   , only: IN=>SHR_KIND_IN, R8=>SHR_KIND_R8, CS=>SHR_KIND_CS, CL=>SHR_KIND_CL
  use shr_file_mod   , only: shr_file_getunit, shr_file_freeunit
  use shr_cal_mod    , only: shr_cal_date2julian, shr_cal_ymdtod2string
  use shr_mpi_mod    , only: shr_mpi_bcast
  use shr_precip_mod , only: shr_precip_partition_rain_snow_ramp
  use shr_strdata_mod, only: shr_strdata_type, shr_strdata_pioinit, shr_strdata_init
  use shr_strdata_mod, only: shr_strdata_setOrbs, shr_strdata_print, shr_strdata_restRead
  use shr_strdata_mod, only: shr_strdata_advance, shr_strdata_restWrite
  use shr_dmodel_mod , only: shr_dmodel_gsmapcreate, shr_dmodel_rearrGGrid
  use shr_dmodel_mod , only: shr_dmodel_translate_list, shr_dmodel_translateAV_list
  use seq_timemgr_mod, only: seq_timemgr_EClockGetData, seq_timemgr_RestartAlarmIsOn

  use datm_shr_mod   , only: datm_shr_getNextRadCDay, datm_shr_esat, datm_shr_CORE2getFactors
  use datm_shr_mod   , only: datamode       ! namelist input
  use datm_shr_mod   , only: decomp         ! namelist input
  use datm_shr_mod   , only: wiso_datm      ! namelist input
  use datm_shr_mod   , only: rest_file      ! namelist input
  use datm_shr_mod   , only: rest_file_strm ! namelist input
  use datm_shr_mod   , only: factorfn       ! namelist input
  use datm_shr_mod   , only: iradsw         ! namelist input
  use datm_shr_mod   , only: nullstr

  ! !PUBLIC TYPES:

  implicit none
  private ! except

  !--------------------------------------------------------------------------
  ! Public interfaces
  !--------------------------------------------------------------------------

  public :: datm_comp_init
  public :: datm_comp_run
  public :: datm_comp_final

  !--------------------------------------------------------------------------
  ! Private data
  !--------------------------------------------------------------------------

  character(CS) :: myModelName = 'atm'   ! user defined model name
  logical       :: firstcall = .true.    ! first call logical
  real(R8)      :: tbotmax               ! units detector
  real(R8)      :: tdewmax               ! units detector
  real(R8)      :: anidrmax              ! existance detector

  character(len=*),parameter :: rpfile = 'rpointer.atm'

  real(R8),parameter :: tKFrz  = SHR_CONST_TKFRZ
  real(R8),parameter :: degtorad = SHR_CONST_PI/180.0_R8
  real(R8),parameter :: pstd   = SHR_CONST_PSTD     ! standard pressure ~ Pa
  real(R8),parameter :: stebol = SHR_CONST_STEBOL   ! Stefan-Boltzmann constant ~ W/m^2/K^4
  real(R8),parameter :: rdair  = SHR_CONST_RDAIR    ! dry air gas constant   ~ J/K/kg
  real(R8),parameter :: avg_c0 =  61.846_R8
  real(R8),parameter :: avg_c1 =   1.107_R8
  real(R8),parameter :: amp_c0 = -21.841_R8
  real(R8),parameter :: amp_c1 =  -0.447_R8
  real(R8),parameter :: phs_c0 =   0.298_R8
  real(R8),parameter :: dLWarc =  -5.000_R8

  real(R8)           :: dTarc(12)
  data   dTarc      / 0.49_R8, 0.06_R8,-0.73_R8,  -0.89_R8,-0.77_R8,-1.02_R8, &
       -1.99_R8,-0.91_R8, 1.72_R8,   2.30_R8, 1.81_R8, 1.06_R8/

  integer(IN) :: kz,ktopo,ku,kv,ktbot,kptem,kshum,kdens,kpbot,kpslv,klwdn
  integer(IN) :: krc,krl,ksc,ksl,kswndr,kswndf,kswvdr,kswvdf,kswnet
  integer(IN) :: kanidr,kanidf,kavsdr,kavsdf
  integer(IN) :: stbot,swind,sz,spbot,sshum,stdew,srh,slwdn,sswdn,sswdndf,sswdndr
  integer(IN) :: sprecc,sprecl,sprecn,sco2p,sco2d,sswup,sprec,starcf

  ! water isotopes / tracer input
  integer(IN) :: kshum_16O, kshum_18O, kshum_HDO
  integer(IN) :: krc_18O, krc_HDO
  integer(IN) :: krl_18O, krl_HDO
  integer(IN) :: ksc_18O, ksc_HDO
  integer(IN) :: ksl_18O, ksl_HDO
  integer(IN) :: srh_16O, srh_18O, srh_HDO
  integer(IN) :: sprecn_16O, sprecn_18O, sprecn_HDO

  ! anomaly forcing
  integer(IN) :: sprecsf
  integer(IN) :: sprec_af,su_af,sv_af,stbot_af,sshum_af,spbot_af,slwdn_af,sswdn_af

  type(mct_rearr)      :: rearr
  type(mct_avect)      :: avstrm   ! av of data from stream
  integer(IN), pointer :: imask(:)
  real(R8), pointer    :: yc(:)
  real(R8), pointer    :: windFactor(:)
  real(R8), pointer    :: winddFactor(:)
  real(R8), pointer    :: qsatFactor(:)

  integer(IN),parameter :: ktrans  = 77

  character(16),parameter  :: avofld(1:ktrans) = &
       (/"Sa_z            ","Sa_topo         ", &
       "Sa_u            ","Sa_v            ","Sa_tbot         ", &
       "Sa_ptem         ","Sa_shum         ","Sa_dens         ","Sa_pbot         ", &
       "Sa_pslv         ","Faxa_lwdn       ","Faxa_rainc      ","Faxa_rainl      ", &
       "Faxa_snowc      ","Faxa_snowl      ","Faxa_swndr      ","Faxa_swvdr      ", &
       "Faxa_swndf      ","Faxa_swvdf      ","Faxa_swnet      ","Sa_co2prog      ", &
       "Sa_co2diag      ","Faxa_bcphidry   ","Faxa_bcphodry   ","Faxa_bcphiwet   ", &
       "Faxa_ocphidry   ","Faxa_ocphodry   ","Faxa_ocphiwet   ","Faxa_dstwet1    ", &
       "Faxa_dstwet2    ","Faxa_dstwet3    ","Faxa_dstwet4    ","Faxa_dstdry1    ", &
       "Faxa_dstdry2    ","Faxa_dstdry3    ","Faxa_dstdry4    ",                    &
       "Sx_tref         ","Sx_qref         ","Sx_avsdr        ","Sx_anidr        ", &
       "Sx_avsdf        ","Sx_anidf        ","Sx_t            ","So_t            ", &
       "Sl_snowh        ","Sf_lfrac        ","Sf_ifrac        ","Sf_ofrac        ", &
       "Faxx_taux       ","Faxx_tauy       ","Faxx_lat        ","Faxx_sen        ", &
       "Faxx_lwup       ","Faxx_evap       ","Fall_fco2_lnd   ","Faoo_fco2_ocn   ", &
       "Faoo_fdms_ocn   ",  &
                                ! add values for bias correction / anomaly forcing
       "Sa_precsf       ", &
       "Sa_prec_af      ","Sa_u_af         ","Sa_v_af         ","Sa_tbot_af      ",&
       "Sa_pbot_af      ","Sa_shum_af      ","Sa_swdn_af      ","Sa_lwdn_af      ",&
                                ! isotopic forcing
       "Faxa_rainc_18O  ","Faxa_rainc_HDO  ","Faxa_rainl_18O  ","Faxa_rainl_HDO  ",&
       "Faxa_snowc_18O  ","Faxa_snowc_HDO  ","Faxa_snowl_18O  ","Faxa_snowl_HDO  ",&
       "Sa_shum_16O     ","Sa_shum_18O     ","Sa_shum_HDO     " &
       /)

  character(16),parameter  :: avifld(1:ktrans) = &
       (/"z               ","topo            ", &
       "u               ","v               ","tbot            ", &
       "ptem            ","shum            ","dens            ","pbot            ", &
       "pslv            ","lwdn            ","rainc           ","rainl           ", &
       "snowc           ","snowl           ","swndr           ","swvdr           ", &
       "swndf           ","swvdf           ","swnet           ","co2prog         ", &
       "co2diag         ","bcphidry        ","bcphodry        ","bcphiwet        ", &
       "ocphidry        ","ocphodry        ","ocphiwet        ","dstwet1         ", &
       "dstwet2         ","dstwet3         ","dstwet4         ","dstdry1         ", &
       "dstdry2         ","dstdry3         ","dstdry4         ",                    &
       "tref            ","qref            ","avsdr           ","anidr           ", &
       "avsdf           ","anidf           ","ts              ","to              ", &
       "snowhl          ","lfrac           ","ifrac           ","ofrac           ", &
       "taux            ","tauy            ","lat             ","sen             ", &
       "lwup            ","evap            ","co2lnd          ","co2ocn          ", &
       "dms             ", &
                                ! add values for bias correction / anomaly forcing (add Sa_precsf for precip scale factor)
       "precsf          ", &
       "prec_af         ","u_af            ","v_af            ","tbot_af         ", &
       "pbot_af         ","shum_af         ","swdn_af         ","lwdn_af         ", &
                                ! isotopic forcing
       "rainc_18O       ","rainc_HDO       ","rainl_18O       ","rainl_HDO       ", &
       "snowc_18O       ","snowc_HDO       ","snowl_18O       ","snowl_HDO       ", &
       "shum_16O        ","shum_18O        ","shum_HDO        " &
       /)

  ! The stofld and stifld lists are used for fields that are read but not passed to the
  ! coupler (e.g., they are used to compute fields that are passed to the coupler), and
  ! other fields used in calculations. Fields that are simply read and passed directly to
  ! the coupler do not need to be in these lists.

  integer(IN),parameter :: ktranss = 33

  character(16),parameter  :: stofld(1:ktranss) = &
       (/"strm_tbot       ","strm_wind       ","strm_z          ","strm_pbot       ", &
       "strm_shum       ","strm_tdew       ","strm_rh         ","strm_lwdn       ", &
       "strm_swdn       ","strm_swdndf     ","strm_swdndr     ","strm_precc      ", &
       "strm_precl      ","strm_precn      ","strm_co2prog    ","strm_co2diag    ", &
       "strm_swup       ","strm_prec       ","strm_tarcf      ", &
                                ! add bias correction / anomaly forcing streams
       "strm_precsf     ", &
       "strm_prec_af    ","strm_u_af       ","strm_v_af       ","strm_tbot_af    ", &
       "strm_pbot_af    ","strm_shum_af    ","strm_swdn_af    ","strm_lwdn_af    ", &
       "strm_rh_18O     ","strm_rh_HDO     ", &
       "strm_precn_16O  ","strm_precn_18O  ","strm_precn_HDO  "  &
       /)

  character(16),parameter  :: stifld(1:ktranss) = &
       (/"tbot            ","wind            ","z               ","pbot            ", &
       "shum            ","tdew            ","rh              ","lwdn            ", &
       "swdn            ","swdndf          ","swdndr          ","precc           ", &
       "precl           ","precn           ","co2prog         ","co2diag         ", &
                                ! add precsf
       "swup            ","prec            ","tarcf           ","precsf          ", &
                                ! add anomaly forcing streams
       "prec_af         ","u_af            ","v_af            ","tbot_af         ", &
       "pbot_af         ","shum_af         ","swdn_af         ","lwdn_af         ", &
                                ! isotopic forcing
       "rh_18O          ","rh_HDO          ", &
       "precn_16O       ","precn_18O       ","precn_HDO       "  &
       /)

  character(CL), pointer :: ilist_av(:)     ! input list for translation
  character(CL), pointer :: olist_av(:)     ! output list for translation
  character(CL), pointer :: ilist_st(:)     ! input list for translation
  character(CL), pointer :: olist_st(:)     ! output list for translation
  integer(IN)  , pointer :: count_av(:)
  integer(IN)  , pointer :: count_st(:)

  save

  !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
CONTAINS
  !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  !===============================================================================
  subroutine datm_comp_init(Eclock, x2a, a2x, &
       seq_flds_x2a_fields, seq_flds_a2x_fields, &
       SDATM, gsmap, ggrid, mpicom, compid, my_task, master_task, &
       inst_suffix, inst_name, logunit, read_restart, &
       scmMode, scmlat, scmlon, &
       orbEccen, orbMvelpp, orbLambm0, orbObliqr, phase, nextsw_cday)

    ! !DESCRIPTION: initialize data atm model
    implicit none

    ! !INPUT/OUTPUT PARAMETERS:
    type(ESMF_Clock)       , intent(in)    :: EClock
    type(mct_aVect)        , intent(inout) :: x2a, a2x
    character(len=*)       , intent(in)    :: seq_flds_x2a_fields ! fields from mediator
    character(len=*)       , intent(in)    :: seq_flds_a2x_fields ! fields to mediator
    type(shr_strdata_type) , intent(inout) :: SDATM               ! model shr_strdata instance (output)
    type(mct_gsMap)        , pointer       :: gsMap               ! model global sep map (output)
    type(mct_gGrid)        , pointer       :: ggrid               ! model ggrid (output)
    integer(IN)            , intent(in)    :: mpicom              ! mpi communicator
    integer(IN)            , intent(in)    :: compid              ! mct comp id
    integer(IN)            , intent(in)    :: my_task             ! my task in mpi communicator mpicom
    integer(IN)            , intent(in)    :: master_task         ! task number of master task
    character(len=*)       , intent(in)    :: inst_suffix         ! char string associated with instance
    character(len=*)       , intent(in)    :: inst_name           ! fullname of current instance (ie. "lnd_0001")
    integer(IN)            , intent(in)    :: logunit             ! logging unit number
    logical                , intent(in)    :: read_restart        ! start from restart
    logical                , intent(in)    :: scmMode             ! single column mode
    real(R8)               , intent(in)    :: scmLat              ! single column lat
    real(R8)               , intent(in)    :: scmLon              ! single column lon
    real(R8)               , intent(in)    :: orbEccen            ! orb eccentricity (unit-less)
    real(R8)               , intent(in)    :: orbMvelpp           ! orb moving vernal eq (radians)
    real(R8)               , intent(in)    :: orbLambm0           ! orb mean long of perhelion (radians)
    real(R8)               , intent(in)    :: orbObliqr           ! orb obliquity (radians)
    integer                , intent(in)    :: phase               ! initialization phase index
    real(R8)               , intent(out)   :: nextsw_cday         ! calendar of next atm sw

    !--- local variables ---
    integer(IN)   :: n,k         ! generic counters
    integer(IN)   :: lsize     ! local size
    integer(IN)   :: kmask       ! field reference
    integer(IN)   :: klat        ! field reference
    integer(IN)   :: kfld        ! fld index
    integer(IN)   :: cnt         ! counter
    integer(IN)   :: idt         ! integer timestep

    logical       :: exists      ! filename existance
    integer(IN)   :: nu          ! unit number
    integer(IN)   :: CurrentYMD  ! model date
    integer(IN)   :: CurrentTOD  ! model sec into model date
    integer(IN)   :: stepno      ! step number
    character(CL) :: calendar    ! calendar type
    character(CL) :: flds_strm

    !--- formats ---
    character(*), parameter :: F00   = "('(datm_comp_init) ',8a)"
    character(*), parameter :: F0L   = "('(datm_comp_init) ',a, l2)"
    character(*), parameter :: F01   = "('(datm_comp_init) ',a,5i8)"
    character(*), parameter :: F02   = "('(datm_comp_init) ',a,4es13.6)"
    character(*), parameter :: F03   = "('(datm_comp_init) ',a,i8,a)"
    character(*), parameter :: F04   = "('(datm_comp_init) ',2a,2i8,'s')"
    character(*), parameter :: F05   = "('(datm_comp_init) ',a,2f10.4)"
    character(*), parameter :: F90   = "('(datm_comp_init) ',73('='))"
    character(*), parameter :: F91   = "('(datm_comp_init) ',73('-'))"
    character(*), parameter :: subName = "(datm_comp_init) "
    !-------------------------------------------------------------------------------

    call t_startf('DATM_INIT')

    if (phase == 1) then
       call t_startf('datm_strdata_init')

       !----------------------------------------------------------------------------
       ! Initialize PIO
       !----------------------------------------------------------------------------

       call shr_strdata_pioinit(SDATM, COMPID)

       !----------------------------------------------------------------------------
       ! Initialize SDATM
       !----------------------------------------------------------------------------

       call seq_timemgr_EClockGetData( EClock, dtime=idt, calendar=calendar )

       ! NOTE: shr_strdata_init calls shr_dmodel_readgrid which reads the data model
       ! grid and from that computes SDATM%gsmap and SDATM%ggrid. DATM%gsmap is created
       ! using the decomp '2d1d' (1d decomp of 2d grid)
       if (scmmode) then
          if (my_task == master_task) then
             write(logunit,F05) ' scm lon lat = ',scmlon,scmlat
          end if
          call shr_strdata_init(SDATM,&
               mpicom, compid, name='atm', &
               scmmode=scmmode,scmlon=scmlon,scmlat=scmlat, &
               calendar=calendar)
       else
          call shr_strdata_init(SDATM,&
               mpicom, compid, name='atm', &
               calendar=calendar)
       endif

       !--- overwrite mask and frac ---
       k = mct_aVect_indexRA(SDATM%grid%data,'mask')
       SDATM%grid%data%rAttr(k,:) = 1.0_R8

       k = mct_aVect_indexRA(SDATM%grid%data,'frac')
       SDATM%grid%data%rAttr(k,:) = 1.0_R8

       if (my_task == master_task) then
          call shr_strdata_print(SDATM,'ATM data')
       endif

       call t_stopf('datm_strdata_init')

       !----------------------------------------------------------------------------
       ! Initialize MCT global seg map, 1d decomp, gsmap
       !----------------------------------------------------------------------------

       call t_startf('datm_initgsmaps')
       if (my_task == master_task) write(logunit,F00) ' initialize gsmaps'
       call shr_sys_flush(logunit)

       ! create a data model global seqmap (gsmap) given the data model global grid sizes
       ! NOTE: gsmap is initialized using the decomp read in from the datm_in namelist
       ! (which by default is "1d")
       call shr_dmodel_gsmapcreate(gsmap, SDATM%nxg*SDATM%nyg, compid, mpicom, decomp)
       lsize = mct_gsmap_lsize(gsmap,mpicom)

       ! create a rearranger from the data model SDATM%gsmap to gsmap
       call mct_rearr_init(SDATM%gsmap, gsmap, mpicom, rearr)
       call t_stopf('datm_initgsmaps')

       !----------------------------------------------------------------------------
       ! Initialize MCT domain
       !----------------------------------------------------------------------------

       call t_startf('datm_initmctdom')
       if (my_task == master_task) write(logunit,F00) 'copy domains'
       call shr_sys_flush(logunit)

       call shr_dmodel_rearrGGrid(SDATM%grid, ggrid, gsmap, rearr, mpicom)
       call t_stopf('datm_initmctdom')

       !----------------------------------------------------------------------------
       ! Initialize MCT attribute vectors
       !----------------------------------------------------------------------------

       call t_startf('datm_initmctavs')
       if (my_task == master_task) write(logunit,F00) 'allocate AVs'
       call shr_sys_flush(logunit)

       call mct_aVect_init(a2x, rList=seq_flds_a2x_fields, lsize=lsize)
       call mct_aVect_zero(a2x)

       kz    = mct_aVect_indexRA(a2x,'Sa_z')
       ktopo = mct_aVect_indexRA(a2x,'Sa_topo')
       ku    = mct_aVect_indexRA(a2x,'Sa_u')
       kv    = mct_aVect_indexRA(a2x,'Sa_v')
       ktbot = mct_aVect_indexRA(a2x,'Sa_tbot')
       kptem = mct_aVect_indexRA(a2x,'Sa_ptem')
       kshum = mct_aVect_indexRA(a2x,'Sa_shum')
       kdens = mct_aVect_indexRA(a2x,'Sa_dens')
       kpbot = mct_aVect_indexRA(a2x,'Sa_pbot')
       kpslv = mct_aVect_indexRA(a2x,'Sa_pslv')
       klwdn = mct_aVect_indexRA(a2x,'Faxa_lwdn')
       krc   = mct_aVect_indexRA(a2x,'Faxa_rainc')
       krl   = mct_aVect_indexRA(a2x,'Faxa_rainl')
       ksc   = mct_aVect_indexRA(a2x,'Faxa_snowc')
       ksl   = mct_aVect_indexRA(a2x,'Faxa_snowl')
       kswndr= mct_aVect_indexRA(a2x,'Faxa_swndr')
       kswndf= mct_aVect_indexRA(a2x,'Faxa_swndf')
       kswvdr= mct_aVect_indexRA(a2x,'Faxa_swvdr')
       kswvdf= mct_aVect_indexRA(a2x,'Faxa_swvdf')
       kswnet= mct_aVect_indexRA(a2x,'Faxa_swnet')

       if (wiso_datm) then  ! water isotopic forcing
          kshum_16O = mct_aVect_indexRA(a2x,'Sa_shum_16O')
          kshum_18O = mct_aVect_indexRA(a2x,'Sa_shum_18O')
          kshum_HDO = mct_aVect_indexRA(a2x,'Sa_shum_HDO')
          krc_18O   = mct_aVect_indexRA(a2x,'Faxa_rainc_18O')
          krc_HDO   = mct_aVect_indexRA(a2x,'Faxa_rainc_HDO')
          krl_18O   = mct_aVect_indexRA(a2x,'Faxa_rainl_18O')
          krl_HDO   = mct_aVect_indexRA(a2x,'Faxa_rainl_HDO')
          ksc_18O   = mct_aVect_indexRA(a2x,'Faxa_snowc_18O')
          ksc_HDO   = mct_aVect_indexRA(a2x,'Faxa_snowc_HDO')
          ksl_18O   = mct_aVect_indexRA(a2x,'Faxa_snowl_18O')
          ksl_HDO   = mct_aVect_indexRA(a2x,'Faxa_snowl_HDO')
       end if

       call mct_aVect_init(x2a, rList=seq_flds_x2a_fields, lsize=lsize)
       call mct_aVect_zero(x2a)

       kanidr = mct_aVect_indexRA(x2a,'Sx_anidr')
       kanidf = mct_aVect_indexRA(x2a,'Sx_anidf')
       kavsdr = mct_aVect_indexRA(x2a,'Sx_avsdr')
       kavsdf = mct_aVect_indexRA(x2a,'Sx_avsdf')

       !--- figure out what's on the standard streams ---
       cnt = 0
       flds_strm = ''
       do n = 1,SDATM%nstreams
          do k = 1,ktranss
             kfld = mct_aVect_indexRA(SDATM%avs(n),trim(stifld(k)),perrWith='quiet')
             if (kfld > 0) then
                cnt = cnt + 1
                if (cnt == 1) then
                   flds_strm = trim(stofld(k))
                else
                   flds_strm = trim(flds_strm)//':'//trim(stofld(k))
                endif
             endif
          enddo
       enddo

       if (my_task == master_task) write(logunit,F00) ' flds_strm = ',trim(flds_strm)
       call shr_sys_flush(logunit)

       call mct_aVect_init(avstrm, rList=flds_strm, lsize=lsize)
       call mct_aVect_zero(avstrm)

       stbot  = mct_aVect_indexRA(avstrm,'strm_tbot'   ,perrWith='quiet')
       swind  = mct_aVect_indexRA(avstrm,'strm_wind'   ,perrWith='quiet')
       sz     = mct_aVect_indexRA(avstrm,'strm_z'      ,perrWith='quiet')
       spbot  = mct_aVect_indexRA(avstrm,'strm_pbot'   ,perrWith='quiet')
       sshum  = mct_aVect_indexRA(avstrm,'strm_shum'   ,perrWith='quiet')
       stdew  = mct_aVect_indexRA(avstrm,'strm_tdew'   ,perrWith='quiet')
       srh    = mct_aVect_indexRA(avstrm,'strm_rh'     ,perrWith='quiet')
       slwdn  = mct_aVect_indexRA(avstrm,'strm_lwdn'   ,perrWith='quiet')
       sswdn  = mct_aVect_indexRA(avstrm,'strm_swdn'   ,perrWith='quiet')
       sswdndf= mct_aVect_indexRA(avstrm,'strm_swdndf' ,perrWith='quiet')
       sswdndr= mct_aVect_indexRA(avstrm,'strm_swdndr' ,perrWith='quiet')
       sprecc = mct_aVect_indexRA(avstrm,'strm_precc'  ,perrWith='quiet')
       sprecl = mct_aVect_indexRA(avstrm,'strm_precl'  ,perrWith='quiet')
       sprecn = mct_aVect_indexRA(avstrm,'strm_precn'  ,perrWith='quiet')
       sco2p  = mct_aVect_indexRA(avstrm,'strm_co2p'   ,perrWith='quiet')
       sco2d  = mct_aVect_indexRA(avstrm,'strm_co2d'   ,perrWith='quiet')
       sswup  = mct_aVect_indexRA(avstrm,'strm_swup'   ,perrWith='quiet')
       sprec  = mct_aVect_indexRA(avstrm,'strm_prec'   ,perrWith='quiet')
       starcf = mct_aVect_indexRA(avstrm,'strm_tarcf'  ,perrWith='quiet')

       ! anomaly forcing
       sprecsf  = mct_aVect_indexRA(avstrm,'strm_precsf'  ,perrWith='quiet')
       sprec_af = mct_aVect_indexRA(avstrm,'strm_prec_af' ,perrWith='quiet')
       su_af    = mct_aVect_indexRA(avstrm,'strm_u_af'    ,perrWith='quiet')
       sv_af    = mct_aVect_indexRA(avstrm,'strm_v_af'    ,perrWith='quiet')
       stbot_af = mct_aVect_indexRA(avstrm,'strm_tbot_af' ,perrWith='quiet')
       spbot_af = mct_aVect_indexRA(avstrm,'strm_pbot_af' ,perrWith='quiet')
       sshum_af = mct_aVect_indexRA(avstrm,'strm_shum_af' ,perrWith='quiet')
       sswdn_af = mct_aVect_indexRA(avstrm,'strm_swdn_af' ,perrWith='quiet')
       slwdn_af = mct_aVect_indexRA(avstrm,'strm_lwdn_af' ,perrWith='quiet')

       if(wiso_datm) then
          ! isotopic forcing
          sprecn_16O = mct_aVect_indexRA(avstrm,'strm_precn_16O',perrWith='quiet')
          sprecn_18O = mct_aVect_indexRA(avstrm,'strm_precn_18O',perrWith='quiet')
          sprecn_HDO = mct_aVect_indexRA(avstrm,'strm_precn_HDO',perrWith='quiet')
          ! Okay here to just use srh_18O and srh_HDO, because the forcing is (should)
          ! just be deltas, applied in lnd_comp_mct to the base tracer
          srh_16O    = mct_aVect_indexRA(avstrm,'strm_rh_16O',perrWith='quiet')
          srh_18O    = mct_aVect_indexRA(avstrm,'strm_rh_18O',perrWith='quiet')
          srh_HDO    = mct_aVect_indexRA(avstrm,'strm_rh_HDO',perrWith='quiet')
       end if

       allocate(imask(lsize))
       allocate(yc(lsize))
       allocate(windFactor(lsize))
       allocate(winddFactor(lsize))
       allocate(qsatFactor(lsize))

       kmask = mct_aVect_indexRA(ggrid%data,'mask')
       imask(:) = nint(ggrid%data%rAttr(kmask,:))
       klat = mct_aVect_indexRA(ggrid%data,'lat')
       yc(:) = ggrid%data%rAttr(klat,:)

       call t_stopf('datm_initmctavs')

       !----------------------------------------------------------------------------
       ! Read restart
       !----------------------------------------------------------------------------

       if (read_restart) then
          if (trim(rest_file)      == trim(nullstr) .and. &
               trim(rest_file_strm) == trim(nullstr)) then
             if (my_task == master_task) then
                write(logunit,F00) ' restart filenames from rpointer'
                call shr_sys_flush(logunit)
                inquire(file=trim(rpfile)//trim(inst_suffix),exist=exists)
                if (.not.exists) then
                   write(logunit,F00) ' ERROR: rpointer file does not exist'
                   call shr_sys_abort(trim(subname)//' ERROR: rpointer file missing')
                endif
                nu = shr_file_getUnit()
                open(nu,file=trim(rpfile)//trim(inst_suffix),form='formatted')
                read(nu,'(a)') rest_file
                read(nu,'(a)') rest_file_strm
                close(nu)
                call shr_file_freeUnit(nu)
                inquire(file=trim(rest_file_strm),exist=exists)
             endif
             call shr_mpi_bcast(rest_file,mpicom,'rest_file')
             call shr_mpi_bcast(rest_file_strm,mpicom,'rest_file_strm')
          else
             ! use namelist already read
             if (my_task == master_task) then
                write(logunit,F00) ' restart filenames from namelist '
                call shr_sys_flush(logunit)
                inquire(file=trim(rest_file_strm),exist=exists)
             endif
          endif

          call shr_mpi_bcast(exists,mpicom,'exists')

          if (exists) then
             if (my_task == master_task) write(logunit,F00) ' reading ',trim(rest_file_strm)
             call shr_strdata_restRead(trim(rest_file_strm),SDATM,mpicom)
          else
             if (my_task == master_task) write(logunit,F00) ' file not found, skipping ',trim(rest_file_strm)
          endif
          call shr_sys_flush(logunit)
       endif

       if (read_restart) then
          call seq_timemgr_EClockGetData( EClock, curr_ymd=CurrentYMD, curr_tod=CurrentTOD)
          call seq_timemgr_EClockGetData( EClock, stepno=stepno, dtime=idt )
          call seq_timemgr_EClockGetData( EClock, calendar=calendar )
          nextsw_cday = datm_shr_getNextRadCDay( CurrentYMD, CurrentTOD, stepno, idt, iradsw, calendar )
       else
          call seq_timemgr_EClockGetData( EClock, curr_cday=nextsw_cday, stepno=stepno )
       endif

    else  ! phase = 2

       call seq_timemgr_EClockGetData( EClock, curr_ymd=CurrentYMD, curr_tod=CurrentTOD)
       call seq_timemgr_EClockGetData( EClock, stepno=stepno, dtime=idt)
       call seq_timemgr_EClockGetData( EClock, calendar=calendar )
       nextsw_cday = datm_shr_getNextRadCDay( CurrentYMD, CurrentTOD, stepno, idt, iradsw, calendar )

    endif

    !----------------------------------------------------------------------------
    ! Set initial atm state
    !----------------------------------------------------------------------------

    call t_adj_detailf(+2)
    call datm_comp_run( &
         EClock = EClock, &
         x2a = x2a, &
         a2x = a2x, &
         SDATM = SDATM, &
         gsmap = gsmap, &
         ggrid = ggrid, &
         mpicom = mpicom, &
         compid = compid, &
         my_task = my_task, &
         master_task = master_task, &
         inst_suffix = inst_suffix, &
         logunit = logunit, &
         orbEccen = orbEccen, &
         orbMvelpp = orbMvelpp, &
         orbLambm0 = orbLambm0, &
         orbObliqr = orbObliqr, &
         nextsw_cday = nextsw_cday)
    call t_adj_detailf(-2)

    call t_stopf('DATM_INIT')

  end subroutine datm_comp_init

  !===============================================================================
  subroutine datm_comp_run(EClock, x2a, a2x, &
       SDATM, gsmap, ggrid, mpicom, compid, my_task, master_task, &
       inst_suffix, logunit, &
       orbEccen, orbMvelpp, orbLambm0, orbObliqr, &
       nextsw_cday, case_name)

    ! !DESCRIPTION: run method for datm model

    implicit none

    ! !INPUT/OUTPUT PARAMETERS:
    type(ESMF_Clock)       , intent(in)    :: EClock
    type(mct_aVect)        , intent(inout) :: x2a
    type(mct_aVect)        , intent(inout) :: a2x
    type(shr_strdata_type) , intent(inout) :: SDATM
    type(mct_gsMap)        , pointer       :: gsMap
    type(mct_gGrid)        , pointer       :: ggrid
    integer(IN)            , intent(in)    :: mpicom           ! mpi communicator
    integer(IN)            , intent(in)    :: compid           ! mct comp id
    integer(IN)            , intent(in)    :: my_task          ! my task in mpi communicator mpicom
    integer(IN)            , intent(in)    :: master_task      ! task number of master task
    character(len=*)       , intent(in)    :: inst_suffix      ! char string associated with instance
    integer(IN)            , intent(in)    :: logunit          ! logging unit number
    real(R8)               , intent(in)    :: orbEccen         ! orb eccentricity (unit-less)
    real(R8)               , intent(in)    :: orbMvelpp        ! orb moving vernal eq (radians)
    real(R8)               , intent(in)    :: orbLambm0        ! orb mean long of perhelion (radians)
    real(R8)               , intent(in)    :: orbObliqr        ! orb obliquity (radians)
    real(R8)               , intent(out)   :: nextsw_cday      ! calendar of next atm sw
    character(CL)          , intent(in), optional :: case_name ! case name

    !--- local ---
    integer(IN)   :: CurrentYMD        ! model date
    integer(IN)   :: CurrentTOD        ! model sec into model date
    integer(IN)   :: yy,mm,dd          ! year month day
    integer(IN)   :: n                 ! indices
    integer(IN)   :: lsize             ! size of attr vect
    integer(IN)   :: idt               ! integer timestep
    real(R8)      :: dt                ! timestep
    logical       :: write_restart     ! restart now
    character(CL) :: rest_file         ! restart_file
    character(CL) :: rest_file_strm    ! restart_file
    integer(IN)   :: nu                ! unit number
    integer(IN)   :: stepno            ! step number
    real(R8)      :: rday              ! elapsed day
    real(R8)      :: cosFactor         ! cosine factor
    real(R8)      :: factor            ! generic/temporary correction factor
    real(R8)      :: avg_alb           ! average albedo
    real(R8)      :: tMin              ! minimum temperature
    character(CL) :: calendar          ! calendar type

    character(len=18) :: date_str
    !--- temporaries
    real(R8)      :: uprime,vprime,swndr,swndf,swvdr,swvdf,ratio_rvrf
    real(R8)      :: tbot,pbot,rtmp,vp,ea,e,qsat,frac

    character(*), parameter :: F00   = "('(datm_comp_run) ',8a)"
    character(*), parameter :: F04   = "('(datm_comp_run) ',2a,2i8,'s')"
    character(*), parameter :: subName = "(datm_comp_run) "
    !-------------------------------------------------------------------------------

    call t_startf('DATM_RUN')

    call t_startf('datm_run1')

    call seq_timemgr_EClockGetData( EClock, curr_ymd=CurrentYMD, curr_tod=CurrentTOD)
    call seq_timemgr_EClockGetData( EClock, curr_yr=yy, curr_mon=mm, curr_day=dd)
    call seq_timemgr_EClockGetData( EClock, stepno=stepno, dtime=idt)
    call seq_timemgr_EClockGetData( EClock, calendar=calendar)
    dt = idt * 1.0_r8
    write_restart = seq_timemgr_RestartAlarmIsOn(EClock)

    call t_stopf('datm_run1')

    !--------------------
    ! ADVANCE ATM
    !--------------------

    call t_barrierf('datm_BARRIER',mpicom)
    call t_startf('datm')

    nextsw_cday = datm_shr_getNextRadCDay( CurrentYMD, CurrentTOD, stepno, idt, iradsw, calendar )

    !--- set data needed for cosz t-interp method ---
    call shr_strdata_setOrbs(SDATM,orbEccen,orbMvelpp,orbLambm0,orbObliqr,idt)

    !--- copy all fields from streams to a2x as default ---

    call t_startf('datm_strdata_advance')
    call shr_strdata_advance(SDATM,currentYMD,currentTOD,mpicom,'datm')
    call t_stopf('datm_strdata_advance')

    call t_barrierf('datm_scatter_BARRIER',mpicom)

    call t_startf('datm_scatter')
    if (firstcall) then
       allocate(ilist_av(SDATM%nstreams))
       allocate(olist_av(SDATM%nstreams))
       allocate(ilist_st(SDATM%nstreams))
       allocate(olist_st(SDATM%nstreams))
       allocate(count_av(SDATM%nstreams))
       allocate(count_st(SDATM%nstreams))
    end if
    do n = 1,SDATM%nstreams
       if (firstcall) then
          call shr_dmodel_translate_list(SDATM%avs(n),a2x,&
               avifld(1:ktrans),avofld,ilist_av(n),olist_av(n),count_av(n))
       end if
       if (count_av(n) > 0) then
          call shr_dmodel_translateAV_list(SDATM%avs(n),a2x,&
               ilist_av(n),olist_av(n),rearr)
       end if
    enddo
    do n = 1,SDATM%nstreams
       if (firstcall) then
          call shr_dmodel_translate_list(SDATM%avs(n),avstrm,&
               stifld(1:ktranss),stofld,ilist_st(n),olist_st(n),count_st(n))
       end if
       if (count_st(n) > 0) then
          call shr_dmodel_translateAV_list(SDATM%avs(n),avstrm,&
               ilist_st(n),olist_st(n),rearr)
       end if
    enddo
    call t_stopf('datm_scatter')

    !-------------------------------------------------
    ! Determine data model behavior based on the mode
    !-------------------------------------------------

    call t_startf('datm_datamode')
    select case (trim(datamode))

    case('COPYALL')
       ! do nothing extra

    case('CORE2_NYF','CORE2_IAF')
       if (firstcall) then
          if (sprec < 1 .or. sswdn < 1) then
             write(logunit,F00) 'ERROR: prec and swdn must be in streams for CORE2'
             call shr_sys_abort(trim(subname)//'ERROR: prec and swdn must be in streams for CORE2')
          endif
          if (trim(datamode) == 'CORE2_IAF' ) then
             if (starcf < 1 ) then
                write(logunit,F00) 'ERROR: tarcf must be in an input stream for CORE2_IAF'
                call shr_sys_abort(trim(subname)//'tarcf must be in an input stream for CORE2_IAF')
             endif
          endif
          call datm_shr_CORE2getFactors(factorFn,windFactor,winddFactor,qsatFactor, &
               mpicom,compid,gsmap,ggrid,SDATM%nxg,SDATM%nyg)
       endif
       call shr_cal_date2julian(currentYMD,currentTOD,rday,calendar)
       rday = mod((rday - 1.0_R8),365.0_R8)
       cosfactor = cos((2.0_R8*SHR_CONST_PI*rday)/365 - phs_c0)

       lsize = mct_avect_lsize(a2x)
       do n = 1,lsize
          a2x%rAttr(kz,n) = 10.0_R8

          !--- correction to NCEP winds based on QSCAT ---
          uprime    = a2x%rAttr(ku,n)*windFactor(n)
          vprime    = a2x%rAttr(kv,n)*windFactor(n)
          a2x%rAttr(ku,n) = uprime*cos(winddFactor(n)*degtorad)- &
               vprime*sin(winddFactor(n)*degtorad)
          a2x%rAttr(kv,n) = uprime*sin(winddFactor(n)*degtorad)+ &
               vprime*cos(winddFactor(n)*degtorad)

          !--- density, tbot, & pslv taken directly from input stream, set pbot ---
          a2x%rAttr(kpbot,n) = a2x%rAttr(kpslv,n)

          !--- correction to NCEP Arctic & Antarctic air T & potential T ---
          if      ( yc(n) < -60.0_R8 ) then
             tMin = (avg_c0 + avg_c1*yc(n)) + (amp_c0 + amp_c1*yc(n))*cosFactor + tKFrz
             a2x%rAttr(ktbot,n) = max(a2x%rAttr(ktbot,n), tMin)
          else if ( yc(n) > 60.0_R8 ) then
             factor = MIN(1.0_R8, 0.1_R8*(yc(n)-60.0_R8) )
             a2x%rAttr(ktbot,n) = a2x%rAttr(ktbot,n) + factor * dTarc(mm)
          endif
          a2x%rAttr(kptem,n) = a2x%rAttr(ktbot,n)

          !---  correction to NCEP relative humidity for heat budget balance ---
          a2x%rAttr(kshum,n) = a2x%rAttr(kshum,n) + qsatFactor(n)

          !--- Dupont correction to NCEP Arctic air T  ---
          !--- don't correct during summer months (July-September)
          !--- ONLY correct when forcing year is 1997->2004
          if (trim(datamode) == 'CORE2_IAF' ) then
             a2x%rAttr(ktbot,n) = a2x%rAttr(ktbot,n) +  avstrm%rAttr(starcf,n)
             a2x%rAttr(kptem,n) = a2x%rAttr(ktbot,n)
          end if

          !-------------------------------------------------------------------------
          ! PRECIPITATION DATA
          !-------------------------------------------------------------------------

          avstrm%rAttr(sprec,n) = avstrm%rAttr(sprec,n)/86400.0_R8        ! convert mm/day to kg/m^2/s

          !  only correct satellite products, do not correct Serreze Arctic data
          if ( yc(n) < 58. ) then
             avstrm%rAttr(sprec,n) = avstrm%rAttr(sprec,n)*1.14168_R8
          endif
          if ( yc(n) >= 58. .and. yc(n) < 68. ) then
             factor = MAX(0.0_R8, 1.0_R8 - 0.1_R8*(yc(n)-58.0_R8) )
             avstrm%rAttr(sprec,n) = avstrm%rAttr(sprec,n)*(factor*(1.14168_R8 - 1.0_R8) + 1.0_R8)
          endif

          a2x%rAttr(krc,n) = 0.0_R8                    ! default zero
          a2x%rAttr(ksc,n) = 0.0_R8
          if (a2x%rAttr(ktbot,n) < tKFrz ) then        ! assign precip to rain/snow components
             a2x%rAttr(krl,n) = 0.0_R8
             a2x%rAttr(ksl,n) = avstrm%rAttr(sprec,n)
          else
             a2x%rAttr(krl,n) = avstrm%rAttr(sprec,n)
             a2x%rAttr(ksl,n) = 0.0_R8
          endif

          !-------------------------------------------------------------------------
          ! RADIATION DATA
          !-------------------------------------------------------------------------

          !--- fabricate required swdn components from net swdn ---
          a2x%rAttr(kswvdr,n) = avstrm%rAttr(sswdn,n)*(0.28_R8)
          a2x%rAttr(kswndr,n) = avstrm%rAttr(sswdn,n)*(0.31_R8)
          a2x%rAttr(kswvdf,n) = avstrm%rAttr(sswdn,n)*(0.24_R8)
          a2x%rAttr(kswndf,n) = avstrm%rAttr(sswdn,n)*(0.17_R8)

          !--- compute net short-wave based on LY08 latitudinally-varying albedo ---
          avg_alb = ( 0.069 - 0.011*cos(2.0_R8*yc(n)*degtorad ) )
          a2x%rAttr(kswnet,n) = avstrm%rAttr(sswdn,n)*(1.0_R8 - avg_alb)

          !--- corrections to GISS sswdn for heat budget balancing ---
          factor = 1.0_R8
          if      ( -60.0_R8 < yc(n) .and. yc(n) < -50.0_R8 ) then
             factor = 1.0_R8 - (yc(n) + 60.0_R8)*(0.05_R8/10.0_R8)
          else if ( -50.0_R8 < yc(n) .and. yc(n) <  30.0_R8 ) then
             factor = 0.95_R8
          else if (  30.0_R8 < yc(n) .and. yc(n) <  40._R8 ) then
             factor = 1.0_R8 - (40.0_R8 - yc(n))*(0.05_R8/10.0_R8)
          endif
          a2x%rAttr(kswnet,n) = a2x%rAttr(kswnet,n)*factor
          a2x%rAttr(kswvdr,n) = a2x%rAttr(kswvdr,n)*factor
          a2x%rAttr(kswndr,n) = a2x%rAttr(kswndr,n)*factor
          a2x%rAttr(kswvdf,n) = a2x%rAttr(kswvdf,n)*factor
          a2x%rAttr(kswndf,n) = a2x%rAttr(kswndf,n)*factor

          !--- correction to GISS lwdn in Arctic ---
          if ( yc(n) > 60._R8 ) then
             factor = MIN(1.0_R8, 0.1_R8*(yc(n)-60.0_R8) )
             a2x%rAttr(klwdn,n) = a2x%rAttr(klwdn,n) + factor * dLWarc
          endif

       enddo   ! lsize

    case('CORE_IAF_JRA', 'CORE_RYF_JRA', 'IAF_JRA_1p5')
       if (firstcall) then
          if (sprec < 1 .or. sswdn < 1) then
             write(logunit,F00) 'ERROR: prec and swdn must be in streams for IAF/RYF JRA'
             call shr_sys_abort(trim(subname)//'ERROR: prec and swdn must be in streams for IAF/RYF JRA')
          endif
          if (trim(factorFn) == 'null') then
             windFactor = 1.0_R8
             winddFactor = 1.0_R8
             qsatFactor = 1.0_R8
          else
             call datm_shr_CORE2getFactors(factorFn,windFactor,winddFactor,qsatFactor, &
                  mpicom,compid,gsmap,ggrid,SDATM%nxg,SDATM%nyg)
          endif
       endif
       call shr_cal_date2julian(currentYMD,currentTOD,rday,calendar)
       rday = mod((rday - 1.0_R8),365.0_R8)
       cosfactor = cos((2.0_R8*SHR_CONST_PI*rday)/365 - phs_c0)

       lsize = mct_avect_lsize(a2x)
       do n = 1,lsize
          a2x%rAttr(kz,n) = 10.0_R8

          !--- density, tbot, & pslv taken directly from input stream, set pbot ---
          a2x%rAttr(kpbot,n) = a2x%rAttr(kpslv,n)

          a2x%rAttr(kptem,n) = a2x%rAttr(ktbot,n)

          !--- density computation for JRA55 forcing ---
          a2x%rAttr(kdens,n) = a2x%rAttr(kpbot,n)/(rdair*a2x%rAttr(ktbot,n) &
               *(1+0.608* a2x%rAttr(kshum,n)))

          !-------------------------------------------------------------------------
          ! PRECIPITATION DATA
          !-------------------------------------------------------------------------

          a2x%rAttr(krc,n) = 0.0_R8                    ! default zero
          a2x%rAttr(ksc,n) = 0.0_R8
          if (a2x%rAttr(ktbot,n) < tKFrz ) then        ! assign precip to rain/snow components
             a2x%rAttr(krl,n) = 0.0_R8
             a2x%rAttr(ksl,n) = avstrm%rAttr(sprec,n)
          else
             a2x%rAttr(krl,n) = avstrm%rAttr(sprec,n)
             a2x%rAttr(ksl,n) = 0.0_R8
          endif

          !-------------------------------------------------------------------------
          ! RADIATION DATA
          !-------------------------------------------------------------------------

          !--- fabricate required swdn components from net swdn ---
          a2x%rAttr(kswvdr,n) = avstrm%rAttr(sswdn,n)*(0.28_R8)
          a2x%rAttr(kswndr,n) = avstrm%rAttr(sswdn,n)*(0.31_R8)
          a2x%rAttr(kswvdf,n) = avstrm%rAttr(sswdn,n)*(0.24_R8)
          a2x%rAttr(kswndf,n) = avstrm%rAttr(sswdn,n)*(0.17_R8)

          !--- compute net short-wave based on LY08 latitudinally-varying albedo ---
          avg_alb = ( 0.069 - 0.011*cos(2.0_R8*yc(n)*degtorad ) )
          a2x%rAttr(kswnet,n) = avstrm%rAttr(sswdn,n)*(1.0_R8 - avg_alb)

       enddo   ! lsize

    case('CLMNCEP')
       if (firstcall) then
          if (swind < 1 .or. stbot < 1) then
             write(logunit,F00) ' ERROR: wind and tbot must be in streams for CLMNCEP'
             call shr_sys_abort(trim(subname)//' ERROR: wind and tbot must be in streams for CLMNCEP')
          endif
          rtmp = maxval(a2x%rAttr(ktbot,:))
          call shr_mpi_max(rtmp,tbotmax,mpicom,'datm_tbot',all=.true.)
          rtmp = maxval(x2a%rAttr(kanidr,:))
          call shr_mpi_max(rtmp,anidrmax,mpicom,'datm_ani',all=.true.)
          if (stdew > 0) then
             rtmp = maxval(avstrm%rAttr(stdew,:))
             call shr_mpi_max(rtmp,tdewmax,mpicom,'datm_tdew',all=.true.)
          endif
          if (my_task == master_task) &
               write(logunit,*) trim(subname),' max values = ',tbotmax,tdewmax,anidrmax
       endif
       lsize = mct_avect_lsize(a2x)
       do n = 1,lsize
          !--- bottom layer height ---
          if (sz < 1) a2x%rAttr(kz,n) = 30.0_R8

          !--- temperature ---
          if (tbotmax < 50.0_R8) a2x%rAttr(ktbot,n) = a2x%rAttr(ktbot,n) + tkFrz
          ! Limit very cold forcing to 180K
          a2x%rAttr(ktbot,n) = max(180._r8, a2x%rAttr(ktbot,n))
          a2x%rAttr(kptem,n) = a2x%rAttr(ktbot,n)

          !--- pressure ---
          if (spbot < 1) a2x%rAttr(kpbot,n) = pstd
          a2x%rAttr(kpslv,n) = a2x%rAttr(kpbot,n)

          !--- u, v wind velocity ---
          a2x%rAttr(ku,n) = avstrm%rAttr(swind,n)/sqrt(2.0_R8)
          a2x%rAttr(kv,n) = a2x%rAttr(ku,n)

          !--- specific humidity ---
          tbot = a2x%rAttr(ktbot,n)
          pbot = a2x%rAttr(kpbot,n)
          if (sshum > 0) then
             e = datm_shr_esat(tbot,tbot)
             qsat = (0.622_R8 * e)/(pbot - 0.378_R8 * e)
             if (qsat < a2x%rAttr(kshum,n)) then
                a2x%rAttr(kshum,n) = qsat
             endif
          else if (srh > 0) then
             e = avstrm%rAttr(srh,n) * 0.01_R8 * datm_shr_esat(tbot,tbot)
             qsat = (0.622_R8 * e)/(pbot - 0.378_R8 * e)
             a2x%rAttr(kshum,n) = qsat
             if(wiso_datm) then
                ! isotopic forcing
                ! For tracer specific humidity, lnd_import_mct expects a delta, so
                ! just keep the delta from the input file - TW
                a2x%rAttr(kshum_16O,n) = avstrm%rAttr(srh_16O,n)
                a2x%rAttr(kshum_18O,n) = avstrm%rAttr(srh_18O,n)
                a2x%rAttr(kshum_HDO,n) = avstrm%rAttr(srh_HDO,n)
             end if
          else if (stdew > 0) then
             if (tdewmax < 50.0_R8) avstrm%rAttr(stdew,n) = avstrm%rAttr(stdew,n) + tkFrz
             e = datm_shr_esat(avstrm%rAttr(stdew,n),tbot)
             qsat = (0.622_R8 * e)/(pbot - 0.378_R8 * e)
             a2x%rAttr(kshum,n) = qsat
          else
             call shr_sys_abort(subname//'ERROR: cannot compute shum')
          endif

          !--- density ---
          vp = (a2x%rAttr(kshum,n)*pbot) / (0.622_R8 + 0.378_R8 * a2x%rAttr(kshum,n))
          a2x%rAttr(kdens,n) = (pbot - 0.378_R8 * vp) / (tbot*rdair)

          !--- downward longwave ---
          if (slwdn < 1) then
             e  = a2x%rAttr(kpslv,n) * a2x%rAttr(kshum,n) / (0.622_R8 + 0.378_R8 * a2x%rAttr(kshum,n))
             ea = 0.70_R8 + 5.95e-05_R8 * 0.01_R8 * e * exp(1500.0_R8/tbot)
             a2x%rAttr(klwdn,n) = ea * stebol * tbot**4
          endif

          !--- shortwave radiation ---
          if (sswdndf > 0 .and. sswdndr > 0) then
             a2x%rAttr(kswndr,n) = avstrm%rAttr(sswdndr,n) * 0.50_R8
             a2x%rAttr(kswvdr,n) = avstrm%rAttr(sswdndr,n) * 0.50_R8
             a2x%rAttr(kswndf,n) = avstrm%rAttr(sswdndf,n) * 0.50_R8
             a2x%rAttr(kswvdf,n) = avstrm%rAttr(sswdndf,n) * 0.50_R8
          elseif (sswdn > 0) then
             ! relationship between incoming NIR or VIS radiation and ratio of
             ! direct to diffuse radiation calculated based on one year's worth of
             ! hourly CAM output from CAM version cam3_5_55
             swndr = avstrm%rAttr(sswdn,n) * 0.50_R8
             ratio_rvrf =   min(0.99_R8,max(0.29548_R8 + 0.00504_R8*swndr  &
                  -1.4957e-05_R8*swndr**2 + 1.4881e-08_R8*swndr**3,0.01_R8))
             a2x%rAttr(kswndr,n) = ratio_rvrf*swndr
             swndf = avstrm%rAttr(sswdn,n) * 0.50_R8
             a2x%rAttr(kswndf,n) = (1._R8 - ratio_rvrf)*swndf

             swvdr = avstrm%rAttr(sswdn,n) * 0.50_R8
             ratio_rvrf =   min(0.99_R8,max(0.17639_R8 + 0.00380_R8*swvdr  &
                  -9.0039e-06_R8*swvdr**2 + 8.1351e-09_R8*swvdr**3,0.01_R8))
             a2x%rAttr(kswvdr,n) = ratio_rvrf*swvdr
             swvdf = avstrm%rAttr(sswdn,n) * 0.50_R8
             a2x%rAttr(kswvdf,n) = (1._R8 - ratio_rvrf)*swvdf
          else
             call shr_sys_abort(subName//'ERROR: cannot compute short-wave down')
          endif

          !--- swnet: a diagnostic quantity ---
          if (anidrmax < 1.0e-8 .or. anidrmax > SHR_CONST_SPVAL * 0.9_R8) then
             a2x%rAttr(kswnet,n) = 0.0_R8
          else
             a2x%rAttr(kswnet,n) = (1.0_R8-x2a%rAttr(kanidr,n))*a2x%rAttr(kswndr,n) + &
                  (1.0_R8-x2a%rAttr(kavsdr,n))*a2x%rAttr(kswvdr,n) + &
                  (1.0_R8-x2a%rAttr(kanidf,n))*a2x%rAttr(kswndf,n) + &
                  (1.0_R8-x2a%rAttr(kavsdf,n))*a2x%rAttr(kswvdf,n)
          endif

          !--- rain and snow ---
          if (sprecc > 0 .and. sprecl > 0) then
             a2x%rAttr(krc,n) = avstrm%rAttr(sprecc,n)
             a2x%rAttr(krl,n) = avstrm%rAttr(sprecl,n)
          elseif (sprecn > 0) then
             a2x%rAttr(krc,n) = avstrm%rAttr(sprecn,n)*0.1_R8
             a2x%rAttr(krl,n) = avstrm%rAttr(sprecn,n)*0.9_R8
          else
             call shr_sys_abort(subName//'ERROR: cannot compute rain and snow')
          endif

          !--- split precip between rain & snow ---
          call shr_precip_partition_rain_snow_ramp(tbot, frac)
          a2x%rAttr(ksc,n) = max(0.0_R8, a2x%rAttr(krc,n)*(1.0_R8 - frac) )
          a2x%rAttr(ksl,n) = max(0.0_R8, a2x%rAttr(krl,n)*(1.0_R8 - frac) )
          a2x%rAttr(krc,n) = max(0.0_R8, a2x%rAttr(krc,n)*(         frac) )
          a2x%rAttr(krl,n) = max(0.0_R8, a2x%rAttr(krl,n)*(         frac) )

       enddo

    end select

    call t_stopf('datm_datamode')

    !----------------------------------------------------------
    ! bias correction / anomaly forcing ( start block )
    !----------------------------------------------------------

    ! modify atmospheric input fields if streams exist

    lsize = mct_avect_lsize(avstrm)

    ! bias correct precipitation relative to observed
    ! (via bias_correct nameslist option)
    if (sprecsf > 0) then
       do n = 1,lsize
          a2x%rAttr(ksc,n) = a2x%rAttr(ksc,n)*min(1.e2_r8,avstrm%rAttr(sprecsf,n))
          a2x%rAttr(ksl,n) = a2x%rAttr(ksl,n)*min(1.e2_r8,avstrm%rAttr(sprecsf,n))
          a2x%rAttr(krc,n) = a2x%rAttr(krc,n)*min(1.e2_r8,avstrm%rAttr(sprecsf,n))
          a2x%rAttr(krl,n) = a2x%rAttr(krl,n)*min(1.e2_r8,avstrm%rAttr(sprecsf,n))

       end do
    endif

    ! adjust atmospheric input fields if anomaly forcing streams exist
    ! (via anomaly_forcing namelist option)

    ! wind
    if (su_af > 0 .and. sv_af > 0) then
       do n = 1,lsize
          a2x%rAttr(ku,n) = a2x%rAttr(ku,n) + avstrm%rAttr(su_af,n)
          a2x%rAttr(kv,n) = a2x%rAttr(kv,n) + avstrm%rAttr(sv_af,n)
       end do
    endif

    ! specific humidity
    if (sshum_af > 0) then
       do n = 1,lsize
          a2x%rAttr(kshum,n) = a2x%rAttr(kshum,n) + avstrm%rAttr(sshum_af,n)

          ! avoid possible negative q values
          if(a2x%rAttr(kshum,n) < 0._r8) then
             a2x%rAttr(kshum,n) = 1.e-6_r8
          endif

       end do
    endif

    ! pressure
    if (spbot_af > 0) then
       do n = 1,lsize
          a2x%rAttr(kpbot,n) = a2x%rAttr(kpbot,n) + avstrm%rAttr(spbot_af,n)
       end do
    endif

    ! temperature
    if (stbot_af > 0) then
       do n = 1,lsize
          a2x%rAttr(ktbot,n) = a2x%rAttr(ktbot,n) + avstrm%rAttr(stbot_af,n)
       end do
    endif

    ! longwave
    if (slwdn_af > 0) then
       do n = 1,lsize
          a2x%rAttr(klwdn,n) = a2x%rAttr(klwdn,n)*avstrm%rAttr(slwdn_af,n)
       end do
    endif

    ! precipitation
    if (sprec_af > 0) then
       do n = 1,lsize
          a2x%rAttr(ksc,n) = a2x%rAttr(ksc,n)*avstrm%rAttr(sprec_af,n)
          a2x%rAttr(ksl,n) = a2x%rAttr(ksl,n)*avstrm%rAttr(sprec_af,n)
          a2x%rAttr(krc,n) = a2x%rAttr(krc,n)*avstrm%rAttr(sprec_af,n)
          a2x%rAttr(krl,n) = a2x%rAttr(krl,n)*avstrm%rAttr(sprec_af,n)
       enddo
    endif
    ! shortwave
    if (sswdn_af > 0) then
       do n = 1,lsize
          a2x%rAttr(kswndr,n) = a2x%rAttr(kswndr,n)*avstrm%rAttr(sswdn_af,n)
          a2x%rAttr(kswvdr,n) = a2x%rAttr(kswvdr,n)*avstrm%rAttr(sswdn_af,n)
          a2x%rAttr(kswndf,n) = a2x%rAttr(kswndf,n)*avstrm%rAttr(sswdn_af,n)
          a2x%rAttr(kswvdf,n) = a2x%rAttr(kswvdf,n)*avstrm%rAttr(sswdn_af,n)
       enddo
    endif
    !----------------------------------------------------------
    ! bias correction / anomaly forcing ( end block )
    !----------------------------------------------------------

    !--------------------
    ! Write restart
    !--------------------

    if (write_restart) then
       call t_startf('datm_restart')
       call shr_cal_ymdtod2string(date_str, yy,mm,dd,currentTOD)

       write(rest_file,"(6a)") &
            trim(case_name), '.datm',trim(inst_suffix),'.r.', trim(date_str), '.nc'
       write(rest_file_strm,"(6a)") &
            trim(case_name), '.datm',trim(inst_suffix),'.rs1.', trim(date_str), '.bin'
       if (my_task == master_task) then
          nu = shr_file_getUnit()
          open(nu,file=trim(rpfile)//trim(inst_suffix),form='formatted')
          write(nu,'(a)') rest_file
          write(nu,'(a)') rest_file_strm
          close(nu)
          call shr_file_freeUnit(nu)
       endif

       if (my_task == master_task) write(logunit,F04) ' writing ',trim(rest_file_strm),currentYMD,currentTOD
       call shr_strdata_restWrite(trim(rest_file_strm),SDATM,mpicom,trim(case_name),'SDATM strdata')
       call shr_sys_flush(logunit)
       call t_stopf('datm_restart')
    endif

    call t_stopf('datm')

    !----------------------------------------------------------------------------
    ! Log output for model date
    ! Reset shr logging to original values
    !----------------------------------------------------------------------------

    call t_startf('datm_run2')
    if (my_task == master_task) then
       write(logunit,F04) trim(myModelName),': model date ', CurrentYMD,CurrentTOD
       call shr_sys_flush(logunit)
    end if

    firstcall = .false.
    call t_stopf('datm_run2')

    call t_stopf('DATM_RUN')

  end subroutine datm_comp_run

  !===============================================================================
  subroutine datm_comp_final(my_task, master_task, logunit)

    ! !DESCRIPTION: finalize method for datm model
    implicit none

    ! !INPUT/OUTPUT PARAMETERS:
    integer(IN) , intent(in) :: my_task     ! my task in mpi communicator mpicom
    integer(IN) , intent(in) :: master_task ! task number of master task
    integer(IN) , intent(in) :: logunit     ! logging unit number

    !--- formats ---
    character(*), parameter :: F00   = "('(datm_comp_final) ',8a)"
    character(*), parameter :: F91   = "('(datm_comp_final) ',73('-'))"
    character(*), parameter :: subName = "(datm_comp_final) "
    !-------------------------------------------------------------------------------

    call t_startf('DATM_FINAL')

    if (my_task == master_task) then
       write(logunit,F91)
       write(logunit,F00) trim(myModelName),': end of main integration loop'
       write(logunit,F91)
    end if

    call t_stopf('DATM_FINAL')

  end subroutine datm_comp_final
  !===============================================================================

end module datm_comp_mod
