module atm_comp_mct 2,31
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_mod
use seq_flds_mod
use seq_flds_indices
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_comp
use cam_control_mod
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_types , only: surface_state, srfflx_state
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
!
! !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(srfflx_state) , pointer :: cam_in(:)
type(surface_state), 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
! (%c = caseid, $y = year, $m = month, $d = day, $s = seconds in day, %t = tape number)
!
character(len=*), parameter :: rsfilename_spec_cam = '%c.cam2.rs.%y-%m-%d-%s.nc' ! cam srf restarts
integer :: nrg = -1 ! Logical unit number for cam srf restart dataset
!
! 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
integer :: id_isop, id_c10h16
!
!================================================================================
CONTAINS
!================================================================================
subroutine atm_init_mct( EClock, cdata_a, x2a_a, a2x_a, NLFilename ) 2,50
use constituents, only: cnst_get_ind
!-----------------------------------------------------------------------
!
! 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)
logical :: single_column
real(r8):: scmlat,scmlon
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.
!-----------------------------------------------------------------------
!
! 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
! Redirect share output to cam log
call spmdinit(mpicom_atm)
if (masterproc) then
inquire(file='atm_modelio.nml',exist=exists)
if (exists) then
iulog = shr_file_getUnit()
call shr_file_setIO('atm_modelio.nml',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_Faxx_fco2_ocn /= 0) then
write(iulog,*)'error co2_readFlux_ocn and index_x2a_Faxx_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)
!
! call shr_orb_print( iyear_AD, eccen, obliq, mvelp )
!
! 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
!
call read_namelist(single_column_in=single_column, scmlat_in=scmlat, scmlon_in=scmlon)
!
! 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 cnst_get_ind ( 'ISOP', id_isop, abort=.false.)
call cnst_get_ind ( 'C10H16', id_c10h16, abort=.false.)
end subroutine atm_init_mct
!================================================================================
subroutine atm_run_mct( EClock, cdata_a, x2a_a, a2x_a) 1,39
!-----------------------------------------------------------------------
!
! Uses
!
use time_manager, only: advance_timestep, get_curr_date, get_curr_calday, &
get_nstep, get_step_size
use scamMod, only: single_column
! use iop, only: scam_use_iop_srf
use pmgrid, only: plev, plevp
use constituents, only: pcnst
use shr_sys_mod, only: shr_sys_flush
!
! 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)
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)
! 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
dosend = (seq_timemgr_EClockDateInSync( EClock, ymd, tod))
! 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
call t_startf ('CAM_adv_timestep')
call advance_timestep()
call t_stopf ('CAM_adv_timestep')
! 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 ) )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( ) 1,1
call cam_final( cam_out, cam_in )
end subroutine atm_final_mct
!================================================================================
subroutine atm_SetgsMap_mct( mpicom_atm, ATMID, GSMap_atm ) 1,5
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 ) 3,11
!-----------------------------------------------------------------------
!
! 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(srfflx_state), 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.
! factors used to convert MEGAN source units [ug C/m2/hr] to CAM srf emis units [kg/m2/sec]
real(r8), parameter :: isop_factor = 1._r8/3600._r8/1.e9_r8* 68._r8/ 60._r8
real(r8), parameter :: c10h16_factor = 1._r8/3600._r8/1.e9_r8*136._r8/120._r8
#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)
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)%icefrac(i) = x2a_a%rAttr(index_x2a_Sx_ifrac, ig)
cam_in(c)%ocnfrac(i) = x2a_a%rAttr(index_x2a_Sx_ofrac, ig)
cam_in(c)%landfrac(i) = x2a_a%rAttr(index_x2a_Sx_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)
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
! heald: add MEGAN VOC fluxes (convert from units C to units X)
if ( id_isop > 0 .and. index_x2a_Faxx_flxvoc1 > 0 ) then
cam_in(c)%cflx(i,id_isop) = -x2a_a%rAttr(index_x2a_Faxx_flxvoc1, ig)*isop_factor
endif
if ( id_c10h16 > 0 .and. index_x2a_Faxx_flxvoc2 > 0 ) then
cam_in(c)%cflx(i,id_c10h16) = -x2a_a%rAttr(index_x2a_Faxx_flxvoc2, ig)*c10h16_factor
endif
if ( index_x2a_Sx_ddvel/=0 .and. n_drydep>0 ) then
cam_in(c)%depvel(i,:n_drydep) = &
x2a_a%rAttr(index_x2a_Sx_ddvel:index_x2a_Sx_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_Faxx_fco2_lnd /= 0) then
cam_in(c)%fco2_lnd(i) = -x2a_a%rAttr(index_x2a_Faxx_fco2_lnd,ig)
end if
if (index_x2a_Faxx_fco2_ocn /= 0) then
cam_in(c)%fco2_ocn(i) = -x2a_a%rAttr(index_x2a_Faxx_fco2_ocn,ig)
end if
if (index_x2a_Faxx_fdms /= 0) then
cam_in(c)%fdms(i) = -x2a_a%rAttr(index_x2a_Faxx_fdms,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_Faxx_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_Faxx_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 ) 3,1
!-------------------------------------------------------------------
!
! Arguments
!
type(surface_state), intent(in) :: cam_out(begchunk:endchunk)
type(mct_aVect) , intent(out) :: 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 ) 1,7
!-------------------------------------------------------------------
!
! 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) 1,7
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, & 1,3
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