atm_comp_mct.F90

Go to the documentation of this file.

module atm_comp_mct

  use pio              , only: file_desc_t, io_desc_t, var_desc_t, pio_double, pio_def_dim, &
                               pio_put_att, pio_enddef, pio_initdecomp, pio_read_darray, pio_freedecomp, &
                               pio_closefile, pio_write_darray, pio_def_var, pio_inq_varid, &
                          pio_noerr, pio_bcast_error, pio_internal_error, pio_seterrorhandling 
  use mct_mod
  use esmf
  use seq_flds_mod
  use seq_cdata_mod
  use seq_infodata_mod
  use seq_timemgr_mod

  use shr_kind_mod     , only: r8 => shr_kind_r8, cl=>shr_kind_cl
  use shr_file_mod     , only: shr_file_getunit, shr_file_freeunit, &
                               shr_file_setlogunit, shr_file_setloglevel, &
                               shr_file_getlogunit, shr_file_getloglevel, &
                     shr_file_setio
  use shr_sys_mod      , only: shr_sys_flush, shr_sys_abort

  use cam_cpl_indices
  use cam_comp
  use cam_instance     , only: cam_instance_init, inst_suffix
  use cam_control_mod  , only: nsrest, adiabatic, ideal_phys, aqua_planet, eccen, obliqr, lambm0, mvelpp
  use radiation        , only: radiation_get, radiation_do, radiation_nextsw_cday
  use phys_grid        , only: get_ncols_p, get_gcol_all_p, & 
                               ngcols, get_gcol_p, get_rlat_all_p, &
                          get_rlon_all_p, get_area_all_p
  use ppgrid           , only: pcols, begchunk, endchunk       
  use dyn_grid         , only: get_horiz_grid_dim_d
  use camsrfexch       , only: cam_out_t, cam_in_t     
  use cam_restart      , only: get_restcase, get_restartdir
  use cam_history      , only: outfld, ctitle
  use abortutils       , only: endrun
  use filenames        , only: interpret_filename_spec, caseid, brnch_retain_casename
#ifdef SPMD
  use spmd_utils       , only: spmdinit, masterproc, iam
  use mpishorthand     , only: mpicom
#else
  use spmd_utils       , only: spmdinit, masterproc, mpicom, iam
#endif
  use time_manager     , only: get_curr_calday, advance_timestep, get_curr_date, get_nstep, &
                               is_first_step, get_step_size, timemgr_init, timemgr_check_restart
  use iofilemod             
  use perf_mod
  use cam_logfile      , only: iulog
  use co2_cycle        , only: c_i, co2_readflux_ocn, co2_readflux_fuel, co2_transport, &
                               co2_time_interp_ocn, co2_time_interp_fuel, data_flux_ocn, data_flux_fuel
  use physconst       ,  only: mwco2
  use runtime_opts     , only: read_namelist
  use phys_control     , only: cam_chempkg_is
  use scammod          , only: single_column,scmlat,scmlon

!
! !PUBLIC TYPES:
  implicit none
  save
  private ! except

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

  public :: atm_init_mct
  public :: atm_run_mct
  public :: atm_final_mct

!--------------------------------------------------------------------------
! Private interfaces
!--------------------------------------------------------------------------

  private :: atm_setgsmap_mct
  private :: atm_import_mct
  private :: atm_export_mct
  private :: atm_domain_mct
  private :: atm_read_srfrest_mct
  private :: atm_write_srfrest_mct

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

  type(cam_in_t) , pointer :: cam_in(:)
  type(cam_out_t), pointer :: cam_out(:)

  type(mct_avect)   :: a2x_a_snap
  type(mct_avect)   :: a2x_a_sum

  integer, parameter  :: nlen = 256     ! Length of character strings
  character(len=nlen) :: fname_srf_cam  ! surface restart filename
  character(len=nlen) :: pname_srf_cam  ! surface restart full pathname

  ! Filename specifier for restart surface file
  character(len=cl) :: rsfilename_spec_cam
!
! Time averaged counter for flux fields
!
  integer :: avg_count
!
! Time averaged flux fields
!  
  character(*), parameter :: a2x_avg_flds = "Faxa_rainc:Faxa_rainl:Faxa_snowc:Faxa_snowl"  
!
! Are all surface types present   
!
  logical :: lnd_present ! if true => land is present
  logical :: ocn_present ! if true => ocean is present

!
!================================================================================
CONTAINS
!================================================================================

  subroutine atm_init_mct( EClock, cdata_a, x2a_a, a2x_a, NLFilename )

    !-----------------------------------------------------------------------
    !
    ! Arguments
    !
    type(esmf_clock),intent(in)                 :: eclock
    type(seq_cdata), intent(inout)              :: cdata_a
    type(mct_avect), intent(inout)              :: x2a_a
    type(mct_avect), intent(inout)              :: a2x_a   
    character(len=*), optional,   intent(IN)    :: nlfilename ! Namelist filename
    !
    ! Locals
    !
    type(mct_gsmap), pointer   :: gsmap_atm
    type(mct_ggrid), pointer   :: dom_a
    type(seq_infodata_type),pointer :: infodata
    integer :: atmid
    integer :: mpicom_atm
    integer :: lsize
    integer :: iradsw
    logical :: exists           ! true if file exists
    real(r8):: nextsw_cday      ! calendar of next atm shortwave
    integer :: stepno           ! time step         
    integer :: dtime_sync       ! integer timestep size
    integer :: currentymd       ! current year-month-day
    integer :: dtime            ! time step increment (sec)
    integer :: atm_cpl_dt       ! driver atm coupling time step 
    integer :: nstep            ! CAM nstep
    real(r8):: caldayp1         ! CAM calendar day for for next cam time step
    integer :: dtime_cam        ! Time-step increment (sec)
    integer :: ymd              ! CAM current date (YYYYMMDD)
    integer :: yr               ! CAM current year
    integer :: mon              ! CAM current month
    integer :: day              ! CAM current day
    integer :: tod              ! CAM current time of day (sec)
    integer :: start_ymd        ! Start date (YYYYMMDD)
    integer :: start_tod        ! Start time of day (sec)
    integer :: ref_ymd          ! Reference date (YYYYMMDD)
    integer :: ref_tod          ! Reference time of day (sec)
    integer :: stop_ymd         ! Stop date (YYYYMMDD)
    integer :: stop_tod         ! Stop time of day (sec)
    logical :: perpetual_run    ! If in perpetual mode or not
    integer :: perpetual_ymd    ! Perpetual date (YYYYMMDD)
    integer :: shrlogunit,shrloglev ! old values
    logical :: first_time = .true.
    character(len=SHR_KIND_CS) :: calendar  ! Calendar type
    character(len=SHR_KIND_CS) :: starttype ! infodata start type
    integer :: lbnum
    integer :: hdim1_d, hdim2_d ! dimensions of rectangular horizontal grid
                                ! data structure, If 1D data structure, then
                                ! hdim2_d == 1.
    character(len=64) :: filein ! Input namelist filename
    !-----------------------------------------------------------------------
    !
    ! Determine cdata points
    !
#if (defined _MEMTRACE)
    if(masterproc) then
      lbnum=1
      call memmon_dump_fort('memmon.out','atm_init_mct:start::',lbnum)
    endif                      
#endif                         
    call seq_cdata_setptrs(cdata_a, id=atmid, mpicom=mpicom_atm, &
         gsmap=gsmap_atm, dom=dom_a, infodata=infodata)

    if (first_time) then
       
       call cam_instance_init(atmid)

       ! Set filename specifier for restart surface file
       ! (%c=caseid, $y=year, $m=month, $d=day, $s=seconds in day)
       rsfilename_spec_cam = '%c.cam' // trim(inst_suffix) // '.rs.%y-%m-%d-%s.nc' 

       ! Determine attribute vector indices

       call cam_cpl_indices_set()

       ! Redirect share output to cam log
       
       call spmdinit(mpicom_atm)
       
       if (masterproc) then
          inquire(file='atm_modelio.nml'//trim(inst_suffix), exist=exists)
          if (exists) then
             iulog = shr_file_getunit()
             call shr_file_setio('atm_modelio.nml'//trim(inst_suffix), iulog)
          endif
          write(iulog,*) "CAM atmosphere model initialization"
       endif
       
       call shr_file_getlogunit(shrlogunit)
       call shr_file_getloglevel(shrloglev)
       call shr_file_setlogunit(iulog)
       ! 
       ! Consistency check                              
       !
       if (co2_readflux_ocn .and. index_x2a_faoo_fco2_ocn /= 0) then
          write(iulog,*)'error co2_readFlux_ocn and index_x2a_Faoo_fco2_ocn cannot both be active'
          call shr_sys_abort()
       end if
       ! 
       ! Get data from infodata object
       !
       call seq_infodata_getdata( infodata,                                           &
            case_name=caseid, case_desc=ctitle,                                       &
            start_type=starttype,                                                     &
            atm_adiabatic=adiabatic,                                                  &
            atm_ideal_phys=ideal_phys,                                                &
            aqua_planet=aqua_planet,                                                  &
            brnch_retain_casename=brnch_retain_casename,                              &
            single_column=single_column, scmlat=scmlat, scmlon=scmlon,                &
            orb_eccen=eccen, orb_mvelpp=mvelpp, orb_lambm0=lambm0, orb_obliqr=obliqr, &
            lnd_present=lnd_present, ocn_present=ocn_present,                         & 
            perpetual=perpetual_run, perpetual_ymd=perpetual_ymd)
       !
       ! Get nsrest from startup type methods
       !
       if (     trim(starttype) == trim(seq_infodata_start_type_start)) then
          nsrest = 0
       else if (trim(starttype) == trim(seq_infodata_start_type_cont) ) then
          nsrest = 1
       else if (trim(starttype) == trim(seq_infodata_start_type_brnch)) then
          nsrest = 3
       else
          write(iulog,*) 'atm_comp_mct: ERROR: unknown starttype'
          call shr_sys_abort()
       end if
       !
       ! Initialize time manager.
       !
       call seq_timemgr_eclockgetdata(eclock, &
                                      start_ymd=start_ymd, start_tod=start_tod, &
                                      ref_ymd=ref_ymd, ref_tod=ref_tod,         &
                                      stop_ymd=stop_ymd, stop_tod=stop_tod,     &
                                      calendar=calendar )
       !
       ! Read namelist
       !
       filein = "atm_in" // trim(inst_suffix)
       call read_namelist(single_column_in=single_column, scmlat_in=scmlat, &
            scmlon_in=scmlon, nlfilename_in=filein)
       !
       ! Initialize cam time manager
       !
       if ( nsrest == 0 )then
          call timemgr_init( calendar_in=calendar, start_ymd=start_ymd, &
                             start_tod=start_tod, ref_ymd=ref_ymd,      &
                             ref_tod=ref_tod, stop_ymd=stop_ymd,        &
                             stop_tod=stop_tod,                         &
                             perpetual_run=perpetual_run,               &
                             perpetual_ymd=perpetual_ymd )
       end if
       !
       ! First phase of cam initialization 
       ! Initialize mpicom_atm, allocate cam_in and cam_out and determine 
       ! atm decomposition (needed to initialize gsmap) 
       ! for an initial run, cam_in and cam_out are allocated in cam_initial
       ! for a restart/branch run, cam_in and cam_out are allocated in restart 
       ! Set defaults then override with user-specified input and initialize time manager
       ! Note that the following arguments are needed to cam_init for timemgr_restart only
       !
       call cam_init( cam_out, cam_in, mpicom_atm, &
                      start_ymd, start_tod, ref_ymd, ref_tod, stop_ymd, stop_tod, &
                      perpetual_run, perpetual_ymd, calendar)
       !
       ! Check consistency of restart time information with input clock
       !
       if (nsrest /= 0) then
          dtime_cam = get_step_size()
          call timemgr_check_restart( calendar, start_ymd, start_tod, ref_ymd, &
                                      ref_tod, dtime_cam, perpetual_run, perpetual_ymd)
       end if
       !
       ! Initialize MCT gsMap, domain and attribute vectors
       !
       call atm_setgsmap_mct( mpicom_atm, atmid, gsmap_atm )
       lsize = mct_gsmap_lsize(gsmap_atm, mpicom_atm)
       !
       ! Initialize MCT domain 
       !
       call atm_domain_mct( lsize, gsmap_atm, dom_a )
       !
       ! Initialize MCT attribute vectors
       !
       call mct_avect_init(a2x_a, rlist=seq_flds_a2x_fields, lsize=lsize)
       call mct_avect_zero(a2x_a)
       
       call mct_avect_init(x2a_a, rlist=seq_flds_x2a_fields, lsize=lsize) 
       call mct_avect_zero(x2a_a)
       
       call mct_avect_init(a2x_a_snap, rlist=a2x_avg_flds, lsize=lsize)
       call mct_avect_zero(a2x_a_snap)
       
       call mct_avect_init(a2x_a_sum , rlist=a2x_avg_flds, lsize=lsize)
       call mct_avect_zero(a2x_a_sum )
       !
       ! Initialize averaging counter
       !
       avg_count = 0
       !
       ! Create initial atm export state
       !
       call atm_export_mct( cam_out, a2x_a )
       !
       ! Set flag to specify that an extra albedo calculation is to be done (i.e. specify active)
       !
       call seq_infodata_putdata(infodata, atm_prognostic=.true.)
       call get_horiz_grid_dim_d(hdim1_d, hdim2_d)
       call seq_infodata_putdata(infodata, atm_nx=hdim1_d, atm_ny=hdim2_d)

       ! Set flag to indicate that CAM will provide carbon and dust deposition fluxes.
       ! This is now hardcoded to .true. since the ability of CICE to read these
       ! fluxes from a file has been removed.
       call seq_infodata_putdata(infodata, atm_aero=.true.)

       !
       ! Set time step of radiation computation as the current calday
       ! This will only be used on the first timestep of an initial run
       !
       if (nsrest == 0) then
          nextsw_cday = get_curr_calday()
          call seq_infodata_putdata( infodata, nextsw_cday=nextsw_cday )
       end if
       
       ! End redirection of share output to cam log
       
       call shr_file_setlogunit(shrlogunit)
       call shr_file_setloglevel(shrloglev)

       first_time = .false.

    else
       
       ! For initial run, run cam radiation/clouds and return
       ! For restart run, read restart x2a_a
       ! Note - a2x_a is computed upon the completion of the previous run - cam_run1 is called
       ! only for the purposes of finishing the flux averaged calculation to compute a2x_a
       ! Note - cam_run1 is called on restart only to have cam internal state consistent with the 
       ! a2x_a state sent to the coupler

       ! Redirect share output to cam log

       call shr_file_getlogunit(shrlogunit)
       call shr_file_getloglevel(shrloglev)
       call shr_file_setlogunit(iulog)

       call seq_timemgr_eclockgetdata(eclock,curr_ymd=currentymd, stepno=stepno, dtime=dtime_sync )
       if (stepno == 0) then
          call atm_import_mct( x2a_a, cam_in )
          call cam_run1( cam_in, cam_out ) 
          call atm_export_mct( cam_out, a2x_a )
       else
          call atm_read_srfrest_mct( eclock, cdata_a, x2a_a, a2x_a )
          call atm_import_mct( x2a_a, cam_in )
          call cam_run1( cam_in, cam_out ) 
       end if

       ! Compute time of next radiation computation, like in run method for exact restart

! tcx was
!       nextsw_cday = radiation_nextsw_cday() 

       call seq_timemgr_eclockgetdata(eclock,dtime=atm_cpl_dt)
       dtime = get_step_size()          
       nstep = get_nstep()
       if (nstep < 1 .or. dtime < atm_cpl_dt) then
          nextsw_cday = radiation_nextsw_cday() 
       else if (dtime == atm_cpl_dt) then
          caldayp1 = get_curr_calday(offset=int(dtime))
          nextsw_cday = radiation_nextsw_cday() 
          if (caldayp1 /= nextsw_cday) nextsw_cday = -1._r8
       else
          call shr_sys_abort('dtime must be less than or equal to atm_cpl_dt')
       end if
       call seq_infodata_putdata( infodata, nextsw_cday=nextsw_cday ) 

       ! End redirection of share output to cam log
       
       call shr_file_setlogunit(shrlogunit)
       call shr_file_setloglevel(shrloglev)
       
    end if

#if (defined _MEMTRACE )
    if(masterproc) then
      lbnum=1
      call memmon_dump_fort('memmon.out','atm_init_mct:end::',lbnum)
      call memmon_reset_addr()
    endif
#endif
    
    call shr_sys_flush(iulog)

 end subroutine atm_init_mct

!================================================================================

  subroutine atm_run_mct( EClock, cdata_a, x2a_a, a2x_a)

    !-----------------------------------------------------------------------
    !
    ! Uses
    !
    use time_manager,    only: advance_timestep, get_curr_date, get_curr_calday, &
                          get_nstep, get_step_size
!   use iop,             only: scam_use_iop_srf
    use pmgrid,          only: plev, plevp
    use constituents,    only: pcnst
    use shr_sys_mod,     only: shr_sys_flush
    use chemistry,       only: chem_reset_fluxes
    use offline_driver,  only: offline_driver_dorun, offline_driver_end_of_data

    ! 
    ! Arguments
    !
    type(esmf_clock)            ,intent(in)    :: eclock
    type(seq_cdata)             ,intent(inout) :: cdata_a
    type(mct_avect)             ,intent(inout) :: x2a_a
    type(mct_avect)             ,intent(inout) :: a2x_a
    !
    ! Local variables
    !
    type(seq_infodata_type),pointer :: infodata
    integer :: lsize           ! size of attribute vector
    integer :: stepno          ! time step          
    integer :: dtime_sync      ! integer timestep size
    integer :: currentymd      ! current year-month-day
    integer :: iradsw          ! shortwave radation frequency (time steps) 
    logical :: dosend          ! true => send data back to driver
    integer :: dtime           ! time step increment (sec)
    integer :: atm_cpl_dt      ! driver atm coupling time step 
    integer :: ymd_sync        ! Sync date (YYYYMMDD)
    integer :: yr_sync         ! Sync current year
    integer :: mon_sync        ! Sync current month
    integer :: day_sync        ! Sync current day
    integer :: tod_sync        ! Sync current time of day (sec)
    integer :: ymd             ! CAM current date (YYYYMMDD)
    integer :: yr              ! CAM current year
    integer :: mon             ! CAM current month
    integer :: day             ! CAM current day
    integer :: tod             ! CAM current time of day (sec)
    integer :: nstep           ! CAM nstep
    integer :: shrlogunit,shrloglev ! old values
    real(r8):: caldayp1        ! CAM calendar day for for next cam time step
    real(r8):: nextsw_cday     ! calendar of next atm shortwave
    logical :: rstwr           ! .true. ==> write restart file before returning
    logical :: nlend           ! Flag signaling last time-step
    logical :: rstwr_sync      ! .true. ==> write restart file before returning
    logical :: nlend_sync      ! Flag signaling last time-step
    logical :: first_time = .true.
    character(len=*), parameter :: subname="atm_run_mct"
    !-----------------------------------------------------------------------
    integer :: lbnum

#if (defined _MEMTRACE)
    if(masterproc) then
      lbnum=1
      call memmon_dump_fort('memmon.out',subname //':start::',lbnum)
    endif
#endif

    ! Redirect share output to cam log
    
    call shr_file_getlogunit(shrlogunit)
    call shr_file_getloglevel(shrloglev)
    call shr_file_setlogunit(iulog)
    
    ! Note that sync clock time should match cam time at end of time step/loop not beginning
    
    call seq_cdata_setptrs(cdata_a, infodata=infodata)
    call seq_timemgr_eclockgetdata(eclock,curr_ymd=ymd_sync,curr_tod=tod_sync, &
       curr_yr=yr_sync,curr_mon=mon_sync,curr_day=day_sync)

    !load orbital parameters
    call seq_infodata_getdata( infodata,                                           &
       orb_eccen=eccen, orb_mvelpp=mvelpp, orb_lambm0=lambm0, orb_obliqr=obliqr)

    nlend_sync = seq_timemgr_stopalarmison(eclock)
    rstwr_sync = seq_timemgr_restartalarmison(eclock)

    ! Map input from mct to cam data structure

    call t_startf('CAM_import')
    call atm_import_mct( x2a_a, cam_in )
    call t_stopf('CAM_import')
    
    ! Cycle over all time steps in the atm coupling interval
    
    dosend = .false.
    do while (.not. dosend)

       ! (re)set surface fluxes of chem tracers here to MEGAN fluxes (from CLM)
       ! or to zero so that fluxes read from file can be added to MEGAN
       call chem_reset_fluxes( x2a_a%rAttr, cam_in )

       ! Determine if dosend
       ! When time is not updated at the beginning of the loop - then return only if
       ! are in sync with clock before time is updated
       
       call get_curr_date( yr, mon, day, tod )
       ymd = yr*10000 + mon*100 + day
       tod = tod
       if( offline_driver_dorun ) then
          dosend = offline_driver_end_of_data()
       else
          dosend = (seq_timemgr_eclockdateinsync( eclock, ymd, tod))
       endif
       
       ! Determine if time to write cam restart and stop
       
       rstwr = .false.
       if (rstwr_sync .and. dosend) rstwr = .true.
       nlend = .false.
       if (nlend_sync .and. dosend) nlend = .true.
       
       ! Single column specific input 
       
       if (single_column) then
          call scam_use_iop_srf( cam_in )
       endif

       ! Run CAM (run2, run3, run4)
       
       call t_startf('CAM_run2')
       call cam_run2( cam_out, cam_in )
       call t_stopf('CAM_run2')

       call t_startf('CAM_run3')
       call cam_run3( cam_out )
       call t_stopf('CAM_run3')
       
       call t_startf('CAM_run4')
       call cam_run4( cam_out, cam_in, rstwr, nlend, &
            yr_spec=yr_sync, mon_spec=mon_sync, day_spec=day_sync, sec_spec=tod_sync)
       call t_stopf('CAM_run4')
       
       ! Advance cam time step 
       if( .not.offline_driver_dorun ) then
          call t_startf('CAM_adv_timestep')
          call advance_timestep()
          call t_stopf('CAM_adv_timestep')
       endif
       
       ! Run cam radiation/clouds (run1)
          
       call t_startf('CAM_run1')
       call cam_run1( cam_in, cam_out ) 
       call t_stopf('CAM_run1')
       
       ! Map output from cam to mct data structures
       
       call t_startf('CAM_export')
       call atm_export_mct( cam_out, a2x_a )
       call t_stopf('CAM_export')
       
       ! Compute snapshot attribute vector for accumulation
       
! don't accumulate on first coupling freq ts1 and ts2
! for consistency with ccsm3 when flxave is off
       nstep = get_nstep()
       if (nstep <= 2) then
          call mct_avect_copy( a2x_a, a2x_a_sum )
          avg_count = 1
       else
          call mct_avect_copy( a2x_a, a2x_a_snap )
          call mct_avect_accum( avin=a2x_a_snap, avout=a2x_a_sum )
          avg_count = avg_count + 1
       endif
       
    end do

    ! Finish accumulation of attribute vector and average and copy accumulation 
    ! field into output attribute vector
    
    call mct_avect_avg( a2x_a_sum, avg_count)
    call mct_avect_copy( a2x_a_sum, a2x_a )
    call mct_avect_zero( a2x_a_sum) 
    avg_count = 0                   
    
    ! Get time of next radiation calculation - albedos will need to be 
    ! calculated by each surface model at this time
    
    call seq_timemgr_eclockgetdata(eclock,dtime=atm_cpl_dt)
    dtime = get_step_size()          
    if (dtime < atm_cpl_dt) then
       nextsw_cday = radiation_nextsw_cday() 
    else if (dtime == atm_cpl_dt) then
       caldayp1 = get_curr_calday(offset=int(dtime))
       nextsw_cday = radiation_nextsw_cday() 
       if (caldayp1 /= nextsw_cday) nextsw_cday = -1._r8
    else
       call shr_sys_abort('dtime must be less than or equal to atm_cpl_dt')
    end if
    call seq_infodata_putdata( infodata, nextsw_cday=nextsw_cday ) 
    
    ! Write merged surface data restart file if appropriate
    
    if (rstwr_sync) then
       call atm_write_srfrest_mct( cdata_a, x2a_a, a2x_a, &
            yr_spec=yr_sync, mon_spec=mon_sync, day_spec=day_sync, sec_spec=tod_sync)
    end if
    
    ! Check for consistency of internal cam clock with master sync clock 
    
    dtime = get_step_size()
    call get_curr_date( yr, mon, day, tod, offset=-dtime )
    ymd = yr*10000 + mon*100 + day
    tod = tod
    if ((.not.seq_timemgr_eclockdateinsync( eclock, ymd, tod )) .and. (.not.offline_driver_dorun))then
       call seq_timemgr_eclockgetdata(eclock, curr_ymd=ymd_sync, curr_tod=tod_sync )
       write(iulog,*)' cam ymd=',ymd     ,'  cam tod= ',tod
       write(iulog,*)'sync ymd=',ymd_sync,' sync tod= ',tod_sync
       call shr_sys_abort( subname//': CAM clock is not in sync with master Sync Clock' )
    end if
    
    ! End redirection of share output to cam log

    call shr_file_setlogunit(shrlogunit)
    call shr_file_setloglevel(shrloglev)

#if (defined _MEMTRACE)
    if(masterproc) then
      lbnum=1
      call memmon_dump_fort('memmon.out',subname //':end::',lbnum)
      call memmon_reset_addr()
    endif
#endif

  end subroutine atm_run_mct

!================================================================================

  subroutine atm_final_mct( EClock, cdata_a, x2a_a, a2x_a)
    type(esmf_clock)            ,intent(in)    :: eclock
    type(seq_cdata)             ,intent(inout) :: cdata_a
    type(mct_avect)             ,intent(inout) :: x2a_a
    type(mct_avect)             ,intent(inout) :: a2x_a

    call cam_final( cam_out, cam_in )

  end subroutine atm_final_mct

!================================================================================

  subroutine atm_setgsmap_mct( mpicom_atm, ATMID, GSMap_atm )
    use phys_grid, only : get_nlcols_p
    !-------------------------------------------------------------------
    !
    ! Arguments
    !
    integer        , intent(in)  :: mpicom_atm
    integer        , intent(in)  :: atmid
    type(mct_gsmap), intent(out) :: gsmap_atm
    !
    ! Local variables
    !
    integer, allocatable :: gindex(:)
    integer :: i, n, c, ncols, sizebuf, nlcols
    integer :: ier            ! error status
    !-------------------------------------------------------------------

    ! Build the atmosphere grid numbering for MCT
    ! NOTE:  Numbering scheme is: West to East and South to North
    ! starting at south pole.  Should be the same as what's used in SCRIP
    
    ! Determine global seg map

    sizebuf=0
    do c = begchunk, endchunk
       ncols = get_ncols_p(c)
       do i = 1,ncols
          sizebuf = sizebuf+1
       end do
    end do

    allocate(gindex(sizebuf))

    n=0
    do c = begchunk, endchunk
       ncols = get_ncols_p(c)
       do i = 1,ncols
          n=n+1
          gindex(n) = get_gcol_p(c,i)
       end do
    end do

    nlcols = get_nlcols_p()
    call mct_gsmap_init( gsmap_atm, gindex, mpicom_atm, atmid, nlcols, ngcols)

    deallocate(gindex)

  end subroutine atm_setgsmap_mct

!===============================================================================

  subroutine atm_import_mct( x2a_a, cam_in )

    !-----------------------------------------------------------------------
    !
    ! Uses  
    !
    use dust_intr,     only: dust_idx1
#if (defined MODAL_AERO)
    use mo_chem_utls,  only: get_spc_ndx
#endif
    use shr_const_mod, only: shr_const_stebol
    use seq_drydep_mod,only: n_drydep
    !
    ! Arguments
    !
    type(mct_avect),    intent(inout) :: x2a_a
    type(cam_in_t),     intent(inout) :: cam_in(begchunk:endchunk)
    !
    ! Local variables
    !    
    integer  :: i,lat,n,c,ig  ! indices
    integer  :: ncols         ! number of columns
    integer  :: dust_ndx
    logical, save :: first_time = .true.

#if (defined MODAL_AERO)
    integer, parameter:: ndst =2
    integer, target   :: spc_ndx(ndst)
#if (defined MODAL_AERO_7MODE)
    integer, pointer  :: dst_a5_ndx, dst_a7_ndx
#elif (defined MODAL_AERO_3MODE)
    integer, pointer  :: dst_a1_ndx, dst_a3_ndx
#endif
#endif
    !-----------------------------------------------------------------------
    !
#if (defined MODAL_AERO)
#if (defined MODAL_AERO_7MODE)
    dst_a5_ndx => spc_ndx(1)
    dst_a7_ndx => spc_ndx(2)
    dst_a5_ndx = get_spc_ndx( 'dst_a5' )
    dst_a7_ndx = get_spc_ndx( 'dst_a7' )
#elif (defined MODAL_AERO_3MODE)
    dst_a1_ndx => spc_ndx(1)
    dst_a3_ndx => spc_ndx(2)
    dst_a1_ndx = get_spc_ndx( 'dst_a1' )
    dst_a3_ndx = get_spc_ndx( 'dst_a3' )
#endif
#endif

    ! ccsm sign convention is that fluxes are positive downward

    ig=1
    do c=begchunk,endchunk
       ncols = get_ncols_p(c) 

       ! initialize constituent surface fluxes to zero
       cam_in(c)%cflx(:,:) = 0._r8 
                                               
       do i =1,ncols                                                               
          cam_in(c)%wsx(i)       = -x2a_a%rAttr(index_x2a_faxx_taux,ig)     
          cam_in(c)%wsy(i)       = -x2a_a%rAttr(index_x2a_faxx_tauy,ig)     
          cam_in(c)%lhf(i)       = -x2a_a%rAttr(index_x2a_faxx_lat, ig)     
          cam_in(c)%shf(i)       = -x2a_a%rAttr(index_x2a_faxx_sen, ig)     
          cam_in(c)%lwup(i)      = -x2a_a%rAttr(index_x2a_faxx_lwup,ig)    
          cam_in(c)%cflx(i,1)    = -x2a_a%rAttr(index_x2a_faxx_evap,ig)                
          cam_in(c)%asdir(i)     =  x2a_a%rAttr(index_x2a_sx_avsdr, ig)  
          cam_in(c)%aldir(i)     =  x2a_a%rAttr(index_x2a_sx_anidr, ig)  
          cam_in(c)%asdif(i)     =  x2a_a%rAttr(index_x2a_sx_avsdf, ig)  
          cam_in(c)%aldif(i)     =  x2a_a%rAttr(index_x2a_sx_anidf, ig)
          cam_in(c)%ts(i)        =  x2a_a%rAttr(index_x2a_sx_t,     ig)  
          cam_in(c)%sst(i)       =  x2a_a%rAttr(index_x2a_so_t,     ig)             
          cam_in(c)%snowhland(i) =  x2a_a%rAttr(index_x2a_sl_snowh, ig)  
          cam_in(c)%snowhice(i)  =  x2a_a%rAttr(index_x2a_si_snowh, ig)  
          cam_in(c)%tref(i)      =  x2a_a%rAttr(index_x2a_sx_tref,  ig)  
          cam_in(c)%qref(i)      =  x2a_a%rAttr(index_x2a_sx_qref,  ig)
          cam_in(c)%u10(i)       =  x2a_a%rAttr(index_x2a_sx_u10,   ig)
          cam_in(c)%icefrac(i)   =  x2a_a%rAttr(index_x2a_sf_ifrac, ig)  
          cam_in(c)%ocnfrac(i)   =  x2a_a%rAttr(index_x2a_sf_ofrac, ig)
     cam_in(c)%landfrac(i)  =  x2a_a%rAttr(index_x2a_sf_lfrac, ig)
          if ( associated(cam_in(c)%ram1) ) &
               cam_in(c)%ram1(i) =  x2a_a%rAttr(index_x2a_sl_ram1 , ig)
          if ( associated(cam_in(c)%fv) ) &
               cam_in(c)%fv(i)   =  x2a_a%rAttr(index_x2a_sl_fv   , ig)
          if ( associated(cam_in(c)%soilw) ) &
               cam_in(c)%soilw(i) =  x2a_a%rAttr(index_x2a_sl_soilw, ig)
          dust_ndx = dust_idx1()
          ! check that dust constituents are actually in the simulation
          if (dust_ndx>0) then
#if (defined MODAL_AERO)
#if (defined MODAL_AERO_7MODE)
            cam_in(c)%cflx(i,dust_ndx   )  = 0.13_r8  &  ! 1st mode, based on Zender et al (2003) Table 1
#elif (defined MODAL_AERO_3MODE)
            cam_in(c)%cflx(i,dust_ndx   )  = 0.032_r8  &  ! 1st mode, based on Zender et al (2003) Table 1
#endif
                                           * (-x2a_a%rAttr(index_x2a_fall_flxdst1, ig) &
                                              -x2a_a%rAttr(index_x2a_fall_flxdst2, ig) &
                                              -x2a_a%rAttr(index_x2a_fall_flxdst3, ig) &
                                              -x2a_a%rAttr(index_x2a_fall_flxdst4, ig))
#if (defined MODAL_AERO_7MODE)
            cam_in(c)%cflx(i,dust_ndx-spc_ndx(1)+spc_ndx(2))  = 0.87_r8 &  ! 2nd mode
#elif (defined MODAL_AERO_3MODE)
            cam_in(c)%cflx(i,dust_ndx-spc_ndx(1)+spc_ndx(2))  = 0.968_r8 &  ! 2nd mode
#endif
                                           * (-x2a_a%rAttr(index_x2a_fall_flxdst1, ig) &
                                              -x2a_a%rAttr(index_x2a_fall_flxdst2, ig) &
                                              -x2a_a%rAttr(index_x2a_fall_flxdst3, ig) &
                                              -x2a_a%rAttr(index_x2a_fall_flxdst4, ig))
#else
       cam_in(c)%cflx(i,dust_ndx   )  = -x2a_a%rAttr(index_x2a_fall_flxdst1, ig)
       cam_in(c)%cflx(i,dust_ndx +1)  = -x2a_a%rAttr(index_x2a_fall_flxdst2, ig)
       cam_in(c)%cflx(i,dust_ndx +2)  = -x2a_a%rAttr(index_x2a_fall_flxdst3, ig)
       cam_in(c)%cflx(i,dust_ndx +3)  = -x2a_a%rAttr(index_x2a_fall_flxdst4, ig)
#endif
          endif

          ! dry dep velocities
          if ( index_x2a_sl_ddvel/=0 .and. n_drydep>0 ) then
             cam_in(c)%depvel(i,:n_drydep) = &
                  x2a_a%rAttr(index_x2a_sl_ddvel:index_x2a_sl_ddvel+n_drydep-1, ig)
          endif
          !
          ! fields needed to calculate water isotopes to ocean evaporation processes
          !
          cam_in(c)%ustar(i) = x2a_a%rAttr(index_x2a_so_ustar,ig)
          cam_in(c)%re(i)    = x2a_a%rAttr(index_x2a_so_re   ,ig)
          cam_in(c)%ssq(i)   = x2a_a%rAttr(index_x2a_so_ssq  ,ig)
          !
          ! bgc scenarios
          !
          if (index_x2a_fall_fco2_lnd /= 0) then
             cam_in(c)%fco2_lnd(i) = -x2a_a%rAttr(index_x2a_fall_fco2_lnd,ig)
          end if
          if (index_x2a_faoo_fco2_ocn /= 0) then
             cam_in(c)%fco2_ocn(i) = -x2a_a%rAttr(index_x2a_faoo_fco2_ocn,ig)
          end if
          if (index_x2a_faoo_fdms_ocn /= 0) then
             cam_in(c)%fdms(i)     = -x2a_a%rAttr(index_x2a_faoo_fdms_ocn,ig)
          end if

          ig=ig+1

       end do
    end do

    ! Get total co2 flux from components,
    ! Note - co2_transport determines if cam_in(c)%cflx(i,c_i(1:4)) is allocated

    if (co2_transport()) then

       ! Interpolate in time for flux data read in
       if (co2_readflux_ocn) then
          call co2_time_interp_ocn
       end if
       if (co2_readflux_fuel) then
          call co2_time_interp_fuel
       end if
       
       ! from ocn : data read in or from coupler or zero
       ! from fuel: data read in or zero
       ! from lnd : through coupler or zero
       do c=begchunk,endchunk
          ncols = get_ncols_p(c)                                                 
          do i=1,ncols                                                               
             
             ! all co2 fluxes in unit kgCO2/m2/s ! co2 flux from ocn 
             if (index_x2a_faoo_fco2_ocn /= 0) then
                cam_in(c)%cflx(i,c_i(1)) = cam_in(c)%fco2_ocn(i)
             else if (co2_readflux_ocn) then 
                ! convert from molesCO2/m2/s to kgCO2/m2/s
                cam_in(c)%cflx(i,c_i(1)) = &
                     -data_flux_ocn%co2flx(i,c)*(1._r8- cam_in(c)%landfrac(i)) &
                     *mwco2*1.0e-3_r8
             else
                cam_in(c)%cflx(i,c_i(1)) = 0._r8
             end if
             
             ! co2 flux from fossil fuel
             if (co2_readflux_fuel) then
                cam_in(c)%cflx(i,c_i(2)) = data_flux_fuel%co2flx(i,c)
             else
                cam_in(c)%cflx(i,c_i(2)) = 0._r8
             end if
             
             ! co2 flux from land (cpl already multiplies flux by land fraction)
             if (index_x2a_fall_fco2_lnd /= 0) then
                cam_in(c)%cflx(i,c_i(3)) = cam_in(c)%fco2_lnd(i)
             else
                cam_in(c)%cflx(i,c_i(3)) = 0._r8
             end if
             
             ! merged co2 flux
             cam_in(c)%cflx(i,c_i(4)) = cam_in(c)%cflx(i,c_i(1)) + &
                                        cam_in(c)%cflx(i,c_i(2)) + &
                                        cam_in(c)%cflx(i,c_i(3))
          end do
       end do
    end if
    !
    ! if first step, determine longwave up flux from the surface temperature 
    !
    if (first_time) then
       if (is_first_step()) then
          do c=begchunk, endchunk
             ncols = get_ncols_p(c)
             do i=1,ncols
                cam_in(c)%lwup(i) = shr_const_stebol*(cam_in(c)%ts(i)**4)
             end do
          end do
       end if
       first_time = .false.
    end if

  end subroutine atm_import_mct

!===============================================================================

  subroutine atm_export_mct( cam_out, a2x_a )

    !-------------------------------------------------------------------
    !
    ! Arguments
    !
    type(cam_out_t),     intent(in)  :: cam_out(begchunk:endchunk) 
    type(mct_avect),     intent(inout) :: a2x_a
    !
    ! Local variables
    !
    integer :: avsize, avnat
    integer :: i,m,c,n,ig       ! indices
    integer :: ncols            ! Number of columns
    !-----------------------------------------------------------------------

    ! Copy from component arrays into chunk array data structure
    ! Rearrange data from chunk structure into lat-lon buffer and subsequently
    ! create attribute vector

    ig=1
    do c=begchunk, endchunk
       ncols = get_ncols_p(c)
       do i=1,ncols
          a2x_a%rAttr(index_a2x_sa_pslv   ,ig) = cam_out(c)%psl(i)
          a2x_a%rAttr(index_a2x_sa_z      ,ig) = cam_out(c)%zbot(i)   
          a2x_a%rAttr(index_a2x_sa_u      ,ig) = cam_out(c)%ubot(i)   
          a2x_a%rAttr(index_a2x_sa_v      ,ig) = cam_out(c)%vbot(i)   
          a2x_a%rAttr(index_a2x_sa_tbot   ,ig) = cam_out(c)%tbot(i)   
          a2x_a%rAttr(index_a2x_sa_ptem   ,ig) = cam_out(c)%thbot(i)  
          a2x_a%rAttr(index_a2x_sa_pbot   ,ig) = cam_out(c)%pbot(i)   
          a2x_a%rAttr(index_a2x_sa_shum   ,ig) = cam_out(c)%qbot(i,1) 
     a2x_a%rAttr(index_a2x_sa_dens   ,ig) = cam_out(c)%rho(i)
          a2x_a%rAttr(index_a2x_faxa_swnet,ig) = cam_out(c)%netsw(i)      
          a2x_a%rAttr(index_a2x_faxa_lwdn ,ig) = cam_out(c)%flwds(i)  
          a2x_a%rAttr(index_a2x_faxa_rainc,ig) = (cam_out(c)%precc(i)-cam_out(c)%precsc(i))*1000._r8
          a2x_a%rAttr(index_a2x_faxa_rainl,ig) = (cam_out(c)%precl(i)-cam_out(c)%precsl(i))*1000._r8
          a2x_a%rAttr(index_a2x_faxa_snowc,ig) = cam_out(c)%precsc(i)*1000._r8
          a2x_a%rAttr(index_a2x_faxa_snowl,ig) = cam_out(c)%precsl(i)*1000._r8
          a2x_a%rAttr(index_a2x_faxa_swndr,ig) = cam_out(c)%soll(i)   
          a2x_a%rAttr(index_a2x_faxa_swvdr,ig) = cam_out(c)%sols(i)   
          a2x_a%rAttr(index_a2x_faxa_swndf,ig) = cam_out(c)%solld(i)  
          a2x_a%rAttr(index_a2x_faxa_swvdf,ig) = cam_out(c)%solsd(i)  

          ! aerosol deposition fluxes
          a2x_a%rAttr(index_a2x_faxa_bcphidry,ig) = cam_out(c)%bcphidry(i)
          a2x_a%rAttr(index_a2x_faxa_bcphodry,ig) = cam_out(c)%bcphodry(i)
          a2x_a%rAttr(index_a2x_faxa_bcphiwet,ig) = cam_out(c)%bcphiwet(i)
          a2x_a%rAttr(index_a2x_faxa_ocphidry,ig) = cam_out(c)%ocphidry(i)
          a2x_a%rAttr(index_a2x_faxa_ocphodry,ig) = cam_out(c)%ocphodry(i)
          a2x_a%rAttr(index_a2x_faxa_ocphiwet,ig) = cam_out(c)%ocphiwet(i)
          a2x_a%rAttr(index_a2x_faxa_dstwet1,ig)  = cam_out(c)%dstwet1(i)
          a2x_a%rAttr(index_a2x_faxa_dstdry1,ig)  = cam_out(c)%dstdry1(i)
          a2x_a%rAttr(index_a2x_faxa_dstwet2,ig)  = cam_out(c)%dstwet2(i)
          a2x_a%rAttr(index_a2x_faxa_dstdry2,ig)  = cam_out(c)%dstdry2(i)
          a2x_a%rAttr(index_a2x_faxa_dstwet3,ig)  = cam_out(c)%dstwet3(i)
          a2x_a%rAttr(index_a2x_faxa_dstdry3,ig)  = cam_out(c)%dstdry3(i)
          a2x_a%rAttr(index_a2x_faxa_dstwet4,ig)  = cam_out(c)%dstwet4(i)
          a2x_a%rAttr(index_a2x_faxa_dstdry4,ig)  = cam_out(c)%dstdry4(i)

          if (index_a2x_sa_co2prog /= 0) then
             a2x_a%rAttr(index_a2x_sa_co2prog,ig) = cam_out(c)%co2prog(i) ! atm prognostic co2
          end if
          if (index_a2x_sa_co2diag /= 0) then
             a2x_a%rAttr(index_a2x_sa_co2diag,ig) = cam_out(c)%co2diag(i) ! atm diagnostic co2
          end if

          ig=ig+1
       end do
    end do
    
  end subroutine atm_export_mct

!===============================================================================

  subroutine atm_domain_mct( lsize, gsMap_a, dom_a )

    !-------------------------------------------------------------------
    !
    ! Arguments
    !
    integer        , intent(in)   :: lsize
    type(mct_gsmap), intent(in)   :: gsmap_a
    type(mct_ggrid), intent(inout):: dom_a  
    !
    ! Local Variables
    !
    integer  :: n,i,c,ncols           ! indices 
    real(r8) :: lats(pcols)           ! array of chunk latitudes
    real(r8) :: lons(pcols)           ! array of chunk longitude
    real(r8) :: area(pcols)           ! area in radians squared for each grid point
    real(r8), pointer  :: data(:)     ! temporary
    integer , pointer  :: idata(:)    ! temporary
    real(r8), parameter:: radtodeg = 180.0_r8/shr_const_pi
    !-------------------------------------------------------------------
    !
    ! Initialize mct atm domain
    !
    call mct_ggrid_init( ggrid=dom_a, coordchars=trim(seq_flds_dom_coord), otherchars=trim(seq_flds_dom_other), lsize=lsize )
    !
    ! Allocate memory
    !
    allocate(data(lsize))
    !
    ! Initialize attribute vector with special value
    !
    call mct_gsmap_orderedpoints(gsmap_a, iam, idata)
    call mct_ggrid_importiattr(dom_a,'GlobGridNum',idata,lsize)
    !
    ! Determine domain (numbering scheme is: West to East and South to North to South pole)
    ! Initialize attribute vector with special value
    !
    data(:) = -9999.0_r8 
    call mct_ggrid_importrattr(dom_a,"lat"  ,data,lsize) 
    call mct_ggrid_importrattr(dom_a,"lon"  ,data,lsize) 
    call mct_ggrid_importrattr(dom_a,"area" ,data,lsize) 
    call mct_ggrid_importrattr(dom_a,"aream",data,lsize) 
    data(:) = 0.0_r8     
    call mct_ggrid_importrattr(dom_a,"mask" ,data,lsize) 
    data(:) = 1.0_r8
    call mct_ggrid_importrattr(dom_a,"frac" ,data,lsize)
    !
    ! Fill in correct values for domain components
    !
    n=0
    do c = begchunk, endchunk
       ncols = get_ncols_p(c)
       call get_rlat_all_p(c, ncols, lats)
       do i=1,ncols
          n = n+1
          data(n) = lats(i)*radtodeg
       end do
    end do
    call mct_ggrid_importrattr(dom_a,"lat",data,lsize) 

    n=0
    do c = begchunk, endchunk
       ncols = get_ncols_p(c)
       call get_rlon_all_p(c, ncols, lons)
       do i=1,ncols
          n = n+1
          data(n) = lons(i)*radtodeg
       end do
    end do
    call mct_ggrid_importrattr(dom_a,"lon",data,lsize) 

    n=0
    do c = begchunk, endchunk
       ncols = get_ncols_p(c)
       call get_area_all_p(c, ncols, area)
       do i=1,ncols
          n = n+1
          data(n) = area(i) 
       end do
    end do
    call mct_ggrid_importrattr(dom_a,"area",data,lsize) 

    n=0
    do c = begchunk, endchunk
       ncols = get_ncols_p(c)
       do i=1,ncols
          n = n+1
          data(n) = 1._r8 ! mask
       end do
    end do
    call mct_ggrid_importrattr(dom_a,"mask"   ,data,lsize) 
    deallocate(data)

  end subroutine atm_domain_mct

!===========================================================================================
!
  subroutine atm_read_srfrest_mct( EClock, cdata_a, x2a_a, a2x_a)
    use cam_pio_utils
    !-----------------------------------------------------------------------
    !
    ! Arguments
    !
    type(esmf_clock),intent(in)    :: eclock
    type(seq_cdata), intent(inout) :: cdata_a
    type(mct_avect), intent(inout) :: x2a_a
    type(mct_avect), intent(inout) :: a2x_a
    ! 
    ! Local variables
    !
    integer         :: npts         ! array size
    integer         :: rcode        ! return error code
    type(mct_avect) :: gdata        ! global/gathered bundle data
    integer         :: yr_spec      ! Current year
    integer         :: mon_spec     ! Current month
    integer         :: day_spec     ! Current day
    integer         :: sec_spec     ! Current time of day (sec)
    !-----------------------------------------------------------------------
    !
    ! Determine and open surface restart dataset
    !
    integer, pointer :: dof(:)
    integer :: lnx, nf_x2a, nf_a2x, k
    real(r8), allocatable :: tmp(:)
    type(file_desc_t) :: file
    type(io_desc_t) :: iodesc
    type(var_desc_t) :: varid
    character(CL)    :: itemc       ! string converted to char
    type(mct_string) :: mstring     ! mct char type



    call seq_timemgr_eclockgetdata( eclock, curr_yr=yr_spec,curr_mon=mon_spec, &
         curr_day=day_spec, curr_tod=sec_spec ) 
    fname_srf_cam = interpret_filename_spec( rsfilename_spec_cam, case=get_restcase(), &
         yr_spec=yr_spec, mon_spec=mon_spec, day_spec=day_spec, sec_spec= sec_spec )
    pname_srf_cam = trim(get_restartdir() )//fname_srf_cam
    call getfil(pname_srf_cam, fname_srf_cam)
    
    call cam_pio_openfile(file, fname_srf_cam, 0)
    call mct_gsmap_orderedpoints(cdata_a%gsmap, iam, dof)
    lnx = mct_gsmap_gsize(cdata_a%gsmap)
    call pio_initdecomp(pio_subsystem, pio_double, (/lnx/), dof, iodesc)
    allocate(tmp(size(dof)))
    deallocate(dof)
    
    nf_x2a = mct_avect_nrattr(x2a_a)

    do k=1,nf_x2a
       call mct_avect_getrlist(mstring,k,x2a_a)
       itemc = mct_string_tochar(mstring)
       call mct_string_clean(mstring)

       call pio_seterrorhandling(file, pio_bcast_error)
       rcode = pio_inq_varid(file,'x2a_'//trim(itemc) ,varid)
       if (rcode == pio_noerr) then
          call pio_read_darray(file, varid, iodesc, tmp, rcode)
          x2a_a%rattr(k,:) = tmp(:)
       else
     if (masterproc) then
             write(iulog,*)'srfrest warning: field ',trim(itemc),' is not on restart file'
             write(iulog,*)'for backwards compatibility will set it to 0'
          end if
          x2a_a%rattr(k,:) = 0._r8
       end if
       call pio_seterrorhandling(file, pio_internal_error)
    end do

    nf_a2x = mct_avect_nrattr(a2x_a)

    do k=1,nf_a2x
       call mct_avect_getrlist(mstring,k,a2x_a)
       itemc = mct_string_tochar(mstring)
       call mct_string_clean(mstring)

       rcode = pio_inq_varid(file,'a2x_'//trim(itemc) ,varid)
       call pio_read_darray(file, varid, iodesc, tmp, rcode)
       a2x_a%rattr(k,:) = tmp(:)
    end do

    call pio_freedecomp(file,iodesc)
    call pio_closefile(file)
    deallocate(tmp)

  end subroutine atm_read_srfrest_mct
!
!===========================================================================================
!
  subroutine atm_write_srfrest_mct( cdata_a, x2a_a, a2x_a, & 
       yr_spec, mon_spec, day_spec, sec_spec)
    use cam_pio_utils
    !-----------------------------------------------------------------------
    !
    ! Arguments
    !
    type(seq_cdata), intent(in) :: cdata_a
    type(mct_avect), intent(in) :: x2a_a
    type(mct_avect), intent(in) :: a2x_a
    integer        , intent(in) :: yr_spec         ! Simulation year
    integer        , intent(in) :: mon_spec        ! Simulation month
    integer        , intent(in) :: day_spec        ! Simulation day
    integer        , intent(in) :: sec_spec        ! Seconds into current simulation day
    !
    ! Local variables
    !
    integer         :: rcode        ! return error code
    type(mct_avect) :: gdata        ! global/gathered bundle data
    !-----------------------------------------------------------------------
    !
    ! Determine and open surface restart dataset
    !

    integer, pointer :: dof(:)
    integer :: nf_x2a, nf_a2x, lnx, dimid(1), k
    type(file_desc_t) :: file
    type(var_desc_t), pointer :: varid_x2a(:), varid_a2x(:)
    type(io_desc_t)  :: iodesc
    character(CL)    :: itemc       ! string converted to char
    type(mct_string) :: mstring     ! mct char type


    fname_srf_cam = interpret_filename_spec( rsfilename_spec_cam, &
         yr_spec=yr_spec, mon_spec=mon_spec, day_spec=day_spec, sec_spec= sec_spec )
    call cam_pio_createfile(file, fname_srf_cam, 0)

    call mct_gsmap_orderedpoints(cdata_a%gsmap, iam, dof)
    lnx = mct_gsmap_gsize(cdata_a%gsmap)
    call pio_initdecomp(pio_subsystem, pio_double, (/lnx/), dof, iodesc)

    deallocate(dof)
    
    nf_x2a = mct_avect_nrattr(x2a_a)
    allocate(varid_x2a(nf_x2a))
    
    rcode = pio_def_dim(file,'x2a_nx',lnx,dimid(1))
    do k = 1,nf_x2a
       call mct_avect_getrlist(mstring,k,x2a_a)
       itemc = mct_string_tochar(mstring)
       call mct_string_clean(mstring)
       rcode = pio_def_var(file,'x2a_'//trim(itemc),pio_double,dimid,varid_x2a(k))
       rcode = pio_put_att(file,varid_x2a(k),"_fillvalue",fillvalue)
    enddo

    nf_a2x = mct_avect_nrattr(a2x_a)
    allocate(varid_a2x(nf_a2x))
    
    rcode = pio_def_dim(file,'a2x_nx',lnx,dimid(1))
    do k = 1,nf_a2x
       call mct_avect_getrlist(mstring,k,a2x_a)
       itemc = mct_string_tochar(mstring)
       call mct_string_clean(mstring)
       rcode = pio_def_var(file,'a2x_'//trim(itemc),pio_double,dimid,varid_a2x(k))
       rcode = pio_put_att(file,varid_a2x(k),"_fillvalue",fillvalue)
    enddo

    rcode = pio_enddef(file)  ! don't check return code, might be enddef already


    do k=1,nf_x2a
       call pio_write_darray(file, varid_x2a(k), iodesc, x2a_a%rattr(k,:), rcode)
    end do

    do k=1,nf_a2x
       call pio_write_darray(file, varid_a2x(k), iodesc, a2x_a%rattr(k,:), rcode)       
    end do

    deallocate(varid_x2a, varid_a2x)

    call pio_freedecomp(file,iodesc)
    call pio_closefile(file)


  end subroutine atm_write_srfrest_mct

!================================================================================

end module atm_comp_mct