!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                   !!
!!                   GNU General Public License                      !!
!!                                                                   !!
!! This file is part of the Flexible Modeling System (FMS).          !!
!!                                                                   !!
!! FMS is free software; you can redistribute it and/or modify       !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation.             !!
!!                                                                   !!
!! FMS is distributed in the hope that it will be useful,            !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of    !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the     !!
!! GNU General Public License for more details.                      !!
!!                                                                   !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to:                                 !!
!!          Free Software Foundation, Inc.                           !!
!!          59 Temple Place, Suite 330                               !!
!!          Boston, MA  02111-1307  USA                              !!
!! or see:                                                           !!
!!          http://www.gnu.org/licenses/gpl.txt                      !!
!!                                                                   !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  coupler_main couples component models and controls the time integration
!
program coupler_main
!-----------------------------------------------------------------------
!                   GNU General Public License                        
!                                                                      
! This program is free software; you can redistribute it and/or modify it and  
! are expected to follow the terms of the GNU General Public License  
! as published by the Free Software Foundation; either version 2 of   
! the License, or (at your option) any later version.                 
!                                                                      
! MOM is distributed in the hope that it will be useful, but WITHOUT    
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  
! or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    
! License for more details.                                           
!                                                                      
! For the full text of the GNU General Public License,                
! write to: Free Software Foundation, Inc.,                           
!           675 Mass Ave, Cambridge, MA 02139, USA.                   
! or see:   http://www.gnu.org/licenses/gpl.html                      
!-----------------------------------------------------------------------
! <CONTACT EMAIL="Bruce.Wyman@noaa.gov"> Bruce Wyman </CONTACT>
! <CONTACT EMAIL="V.Balaji@noaa.gov"> V. Balaji </CONTACT>

! <HISTORY SRC="http://www.gfdl.noaa.gov/fms-cgi-bin/cvsweb.cgi/FMS/"/>

! <OVERVIEW>
!  A main program that couples component models for atmosphere, ocean, land, 
!  and sea ice on independent grids. 
! </OVERVIEW>

! <DESCRIPTION>
!  This version couples model components representing atmosphere, ocean, land 
!  and sea ice on independent grids. Each model component is represented by a 
!  data type giving the instantaneous model state.
!
!  The component models are coupled to allow implicit vertical diffusion of 
!  heat and moisture at the interfaces of the atmosphere, land, and ice models. 
!  As a result, the atmosphere, land, and ice models all use the same time step. 
!  The atmospheric model has been separated into down and up calls that 
!  correspond to the down and up sweeps of the standard tridiagonal elimination.
!
!  The ocean interface uses explicit mixing. Fluxes to and from the ocean must
!  be passed through the ice model. This includes atmospheric fluxes as well as 
!  fluxes from the land to the ocean (runoff).
!
!  This program contains the model's main time loop. Each iteration of the 
!  main time loop is one coupled (slow) time step. Within this slow time step 
!  loop is a fast time step loop, using the atmospheric time step, where the
!  tridiagonal vertical diffusion equations are solved. Exchange between sea 
!  ice and ocean occurs once every slow timestep.
!
! <PRE>
!      MAIN PROGRAM EXAMPLE
!      --------------------
!
!         DO slow time steps (ocean)
!
!              call flux_ocean_to_ice
!
!              call ICE_SLOW_UP
!
!              DO fast time steps (atmos)
!
!                   call flux_calculation
!
!                   call ATMOS_DOWN
!
!                   call flux_down_from_atmos
!
!                   call LAND_FAST
!
!                   call ICE_FAST
!
!                   call flux_up_to_atmos
!
!                   call ATMOS_UP
!
!              END DO
!
!              call ICE_SLOW_DN
!
!              call flux_ice_to_ocean
!
!              call OCEAN
!
!         END DO

!  </PRE>

! </DESCRIPTION>
! <INFO>
!   <NOTE>
!     <PRE>
!   1.If no value is set for current_date, start_date, or calendar (or default value 
!     specified) then the value from restart file "INPUT/coupler.res" will be used. 
!     If neither a namelist value or restart file value exist the program will fail. 
!   2.The actual run length will be the sum of months, days, hours, minutes, and 
!     seconds. A run length of zero is not a valid option. 
!   3.The run length must be an intergal multiple of the coupling timestep dt_cpld. 
!     </PRE>
!   </NOTE>

!   <ERROR MSG="no namelist value for current_date " STATUS="FATAL">
!     A namelist value for current_date must be given if no restart file for
!     coupler_main (INPUT/coupler.res) is found. 
!   </ERROR>
!   <ERROR MSG="invalid namelist value for calendar" STATUS="FATAL">
!     The value of calendar must be 'julian', 'noleap', or 'thirty_day'. 
!     See the namelist documentation. 
!   </ERROR>
!   <ERROR MSG="no namelist value for calendar" STATUS="FATAL">
!     If no restart file is present, then a namelist value for calendar 
!     must be specified. 
!   </ERROR>
!   <ERROR MSG="initial time is greater than current time" STATUS="FATAL">
!     If a restart file is present, then the namelist value for either 
!     current_date or start_date was incorrectly set. 
!   </ERROR>
!   <ERROR MSG="run length must be multiple of ocean time step " STATUS="FATAL">
!     There must be an even number of ocean time steps for the requested run length. 
!   </ERROR>
!   <ERROR MSG="final time does not match expected ending time " STATUS="WARNING">
!     This error should probably not occur because of checks done at initialization time. 
!   </ERROR>

! </INFO>

  use constants_mod,           only: constants_init

  use time_manager_mod,        only: time_type, set_calendar_type, set_time
  use time_manager_mod,        only: set_date, get_date, days_in_month, month_name
  use time_manager_mod,        only: operator(+), operator(-), operator (<)
  use time_manager_mod,        only: operator (>), operator ( /= ), operator ( / )
  use time_manager_mod,        only: operator (*), THIRTY_DAY_MONTHS, JULIAN
  use time_manager_mod,        only: NOLEAP, NO_CALENDAR, INVALID_CALENDAR
  use time_manager_mod,        only: date_to_string, increment_date
  use time_manager_mod,        only: operator(>=), operator(<=), operator(==)

  use fms_mod,                 only: open_namelist_file, field_exist, file_exist, check_nml_error
  use fms_mod,                 only: uppercase, error_mesg, write_version_number
  use fms_mod,                 only: fms_init, fms_end, stdout
  use fms_mod,                 only: read_data, write_data

  use fms_io_mod,              only: fms_io_exit
  use fms_io_mod,              only: restart_file_type, register_restart_field, save_restart

  use diag_manager_mod,        only: diag_manager_init, diag_manager_end, diag_grid_end
  use diag_manager_mod,        only: DIAG_OCEAN, DIAG_OTHER, DIAG_ALL, get_base_date

  use field_manager_mod,       only: MODEL_ATMOS, MODEL_LAND, MODEL_ICE

  use tracer_manager_mod,      only: tracer_manager_init, get_tracer_index
  use tracer_manager_mod,      only: get_number_tracers, get_tracer_names, NO_TRACER

  use coupler_types_mod,       only: coupler_types_init

  use data_override_mod,       only: data_override_init

!
! model interfaces used to couple the component models:
!               atmosphere, land, ice, and ocean
!

  use atmos_model_mod,         only: atmos_model_init, atmos_model_end
  use atmos_model_mod,         only: update_atmos_model_down
  use atmos_model_mod,         only: update_atmos_model_up
  use atmos_model_mod,         only: atmos_data_type
  use atmos_model_mod,         only: land_ice_atmos_boundary_type
  use atmos_model_mod,         only: atmos_model_restart

  use land_model_mod,          only: land_model_init, land_model_end
  use land_model_mod,          only: land_data_type, atmos_land_boundary_type
  use land_model_mod,          only: update_land_model_fast, update_land_model_slow
  use land_model_mod,          only: land_model_restart

  use ice_model_mod,           only: ice_model_init, ice_model_end
  use ice_model_mod,           only: update_ice_model_slow_up
  use ice_model_mod,           only: update_ice_model_fast
  use ice_model_mod,           only: update_ice_model_slow_dn
  use ice_model_mod,           only: ice_data_type, land_ice_boundary_type
  use ice_model_mod,           only: ocean_ice_boundary_type, atmos_ice_boundary_type
  use ice_model_mod,           only: ice_model_restart

  use ocean_model_mod,         only: update_ocean_model, ocean_model_init
  use ocean_model_mod,         only: ocean_model_end, ocean_public_type, ocean_state_type, ice_ocean_boundary_type
  use ocean_model_mod,         only: ocean_model_restart
!
! flux_ calls translate information between model grids - see flux_exchange.f90
!

  use flux_exchange_mod,       only: flux_exchange_init
  use flux_exchange_mod,       only: sfc_boundary_layer
  use flux_exchange_mod,       only: generate_sfc_xgrid
  use flux_exchange_mod,       only: flux_down_from_atmos
  use flux_exchange_mod,       only: flux_up_to_atmos
  use flux_exchange_mod,       only: flux_land_to_ice
  use flux_exchange_mod,       only: flux_ice_to_ocean
  use flux_exchange_mod,       only: flux_ocean_to_ice
  use flux_exchange_mod,       only: flux_check_stocks, flux_init_stocks, flux_ice_to_ocean_stocks, flux_ocean_from_ice_stocks

  use atmos_tracer_driver_mod, only: atmos_tracer_driver_gather_data

  use mpp_mod,                 only: mpp_clock_id, mpp_clock_begin, mpp_clock_end, mpp_chksum
  use mpp_mod,                 only: mpp_init, mpp_pe, mpp_npes, mpp_root_pe, MAXPES
  use mpp_mod,                 only: stderr, stdlog, mpp_error, NOTE, FATAL, WARNING
  use mpp_mod,                 only: mpp_set_current_pelist, mpp_declare_pelist

  use mpp_io_mod,              only: mpp_open, mpp_close, mpp_io_clock_on
  use mpp_io_mod,              only: MPP_NATIVE, MPP_RDONLY, MPP_DELETE

  use mpp_domains_mod,         only: mpp_broadcast_domain

  use memutils_mod,            only: print_memuse_stats

  implicit none

!-----------------------------------------------------------------------

  character(len=128) :: version = '$Id: coupler_main.F90,v 17.0.4.2.2.1 2009/11/30 20:17:46 z1l Exp $'
  character(len=128) :: tag = '$Name: mom4p1_pubrel_dec2009_nnz $'

!-----------------------------------------------------------------------
!---- model defined-types ----

  type (atmos_data_type) :: Atm
  type  (land_data_type) :: Land
  type   (ice_data_type) :: Ice
  ! allow members of ocean type to be aliased (ap)
  type (ocean_public_type), target :: Ocean
  type (ocean_state_type),  pointer :: Ocean_state => NULL()

  type(atmos_land_boundary_type)     :: Atmos_land_boundary
  type(atmos_ice_boundary_type)      :: Atmos_ice_boundary
  type(land_ice_atmos_boundary_type) :: Land_ice_atmos_boundary
  type(land_ice_boundary_type)       :: Land_ice_boundary
  type(ice_ocean_boundary_type)      :: Ice_ocean_boundary
  type(ocean_ice_boundary_type)      :: Ocean_ice_boundary

!-----------------------------------------------------------------------
! ----- coupled model time -----

  type (time_type) :: Time, Time_init, Time_end, &
                      Time_step_atmos, Time_step_cpld
  type(time_type) :: Time_atmos, Time_ocean
  integer :: num_atmos_calls, na
  integer :: num_cpld_calls, nc

!------ for intermediate restart
  type(restart_file_type), allocatable :: Ice_bc_restart(:), Ocn_bc_restart(:)
  character(len=64),       allocatable :: ice_bc_restart_file(:), ocn_bc_restart_file(:) 
  integer                              :: num_ice_bc_restart=0, num_ocn_bc_restart=0
  type(time_type)                      :: Time_restart, Time_restart_current, Time_start
  character(len=32)                    :: timestamp

! ----- coupled model initial date -----

  integer :: date_init(6)
  integer :: calendar_type = INVALID_CALENDAR

!-----------------------------------------------------------------------
!------ namelist interface -------

! <NAMELIST NAME="coupler_nml">
!   <DATA NAME="current_date"  TYPE="integer, dimension(6)"  DEFAULT="0">
!     The date that the current integration starts with. 
!   </DATA>
!   <DATA NAME="force_date_from_namelist"  TYPE="logical"  DEFAULT=".false.">
!     Flag that determines whether the namelist variable current_date should 
!     override the date in the restart file INPUT/coupler.res. If the restart 
!     file does not exist then force_date_from_namelist has not effect, the value of current_date 
!     will be used.
!   </DATA>
!   <DATA NAME="calendar"  TYPE="character(maxlen=17)"  DEFAULT="''">
!     The calendar type used by the current integration. Valid values are consistent 
!     with the time_manager module: 'julian', 'noleap', or 'thirty_day'. The value 
!     'no_calendar' can not be used because the time_manager's date  function are used. 
!     All values must be lowercase.
!   </DATA>
!   <DATA NAME="months "  TYPE="integer"  DEFAULT="0">
!     The number of months that the current integration will be run for. 
!   </DATA>
!   <DATA NAME="days "  TYPE="integer"  DEFAULT="0">
!     The number of days that the current integration will be run for. 
!   </DATA>
!   <DATA NAME="hours"  TYPE="integer"  DEFAULT="0">
!     The number of hours that the current integration will be run for. 
!   </DATA>
!   <DATA NAME="minutes "  TYPE="integer"  DEFAULT="0">
!     The number of minutes that the current integration will be run for. 
!   </DATA>
!   <DATA NAME="seconds"  TYPE="integer"  DEFAULT="0">
!     The number of seconds that the current integration will be run for. 
!   </DATA>
!   <DATA NAME="dt_atmos"  TYPE="integer"  DEFAULT="0">
!     Atmospheric model time step in seconds, including the fast coupling with 
!     land and sea ice. 
!   </DATA>
!   <DATA NAME="dt_cpld"  TYPE="integer"  DEFAULT="0">
!     Time step in seconds for coupling between ocean and atmospheric models: 
!     must be an integral multiple of dt_atmos and dt_ocean. This is the "slow" timestep.
!   </DATA>
!  <DATA NAME="do_atmos, do_ocean, do_ice, do_land, do_flux" TYPE="logical">
!  If true (default), that particular model component (atmos, etc.) is run.
!  If false, the execution of that component is skipped. This is used when
!  ALL the output fields sent by that component to the coupler have been
!  overridden using the data_override feature. For advanced users only:
!  if you're not sure, you should leave these values at TRUE.
!  </DATA> 
!  <DATA NAME="concurrent" TYPE="logical">
!  If true, the ocean executes concurrently with the atmosphere-land-ocean
!   on a separate set of PEs.
!  If false (default), the execution is serial: call atmos... followed by
!  call ocean...
!  If using concurrent execution, you must set one of
!   atmos_npes and ocean_npes, see below.
!  </DATA> 
!  <DATA NAME="atmos_npes, ocean_npes" TYPE="integer">
!  If concurrent is set to true, we use these to set the list of PEs on which
!   each component runs.
!  At least one of them must be set to a number between 0 and NPES.
!  If exactly one of these two is set non-zero, the other is set to the
!   remainder from NPES.
!  If both are set non-zero they must add up to NPES.
!  </DATA> 
!  <DATA NAME="use_lag_fluxes" TYPE="logical">
!  If true, then mom4 is forced with SBCs from one coupling timestep ago
!  If false, then mom4 is forced with most recent SBCs.
!  For a leapfrog MOM coupling with dt_cpld=dt_ocean, lag fluxes
!  can be shown to be stable and current fluxes to be unconditionally unstable.
!  For dt_cpld>dt_ocean there is probably sufficient damping.
!  use_lag_fluxes is set to TRUE by default.
!  </DATA>
!  <DATA NAME="n_mask" TYPE="integer">
!    number of region to be masked out. Its value should be less than MAX_PES.
!  </DATA>
!  <DATA NAME="mask_list(2,MAXPES)" TYPE="integer, dimension(2,MAX_MASK_REGION)">
!    The position of the region to be masked out. mask_list(1,:) is the x-layout position
!    and mask_list(2,:) is y-layout position.  
!  </DATA>
!  <DATA NAME="layout_mask" TYPE="integer, dimension(2)">
!   Processor domain layout for all the component model. layout_mask need to be set when and only 
!   when n_mask is greater than 0 ( some domain region is masked out ). When this namelist is set,
!   it will overload the layout in each component model. The default value is (0,0).
!   Currently we require all the component model has the same layout and same grid size.
!  </DATA>
!  <DATA NAME="restart_interval" TYPE="integer, dimension(6)"  DEFAULT="0">
!     The time interval that write out intermediate restart file. The format is (yr,mo,day,hr,min,sec).
!     When restart_interval is all zero, no intermediate restart file will be written out.
!   </DATA>
!   <NOTE>
!     <PRE>
!     1.If no value is set for current_date, start_date, or calendar (or default value specified) then the value from restart
!       file "INPUT/coupler.res" will be used. If neither a namelist value or restart file value exist the program will fail. 
!     2.The actual run length will be the sum of months, days, hours, minutes, and seconds. A run length of zero is not a
!       valid option. 
!     3.The run length must be an intergal multiple of the coupling timestep dt_cpld. 
!     </PRE>
!   </NOTE>
! </NAMELIST>

  integer, dimension(6) :: restart_interval = (/ 0, 0, 0, 0, 0, 0/)
  integer, dimension(6) :: current_date     = (/ 0, 0, 0, 0, 0, 0 /)
  character(len=17) :: calendar = '                 '
  logical :: force_date_from_namelist = .false.  ! override restart values for date
  integer :: months=0, days=0, hours=0, minutes=0, seconds=0
  integer :: dt_atmos = 0  ! fluxes passed between atmosphere & ice/land
  integer :: dt_cpld  = 0  ! fluxes passed between ice & ocean


  integer ::atmos_npes=0, ocean_npes=0, ice_npes=0, land_npes=0
  logical :: do_atmos =.true., do_land =.true., do_ice =.true., do_ocean=.true.
  logical :: do_flux =.true.
  logical :: concurrent=.FALSE.
  logical :: use_lag_fluxes=.TRUE.
  logical :: do_chksum=.FALSE.
  integer :: layout_mask(2) = (/0 , 0/)
  integer :: n_mask = 0
  integer :: mask_list(2, MAXPES), n, m 
  integer :: check_stocks = 0 ! -1: never 0: at end of run only n>0: every n coupled steps
  integer, parameter :: mp = 2*MAXPES
  data ((mask_list(n,m),n=1, 2),m=1,MAXPES) /mp*0/

  namelist /coupler_nml/ current_date, calendar, force_date_from_namelist, months, days, hours,      &
                         minutes, seconds, dt_cpld, dt_atmos, do_atmos,              &
                         do_land, do_ice, do_ocean, do_flux, atmos_npes, ocean_npes, &
                         ice_npes, land_npes, concurrent, use_lag_fluxes, do_chksum, &
                         n_mask, layout_mask, mask_list, check_stocks, restart_interval

  integer :: initClock, mainClock, termClock

  character(len=80) :: text
  character(len=48), parameter                    :: mod_name = 'coupler_main_mod'
 
  integer :: ensemble_id = 1 , outunit
  integer, allocatable :: ensemble_pelist(:, :) 

!
!-----------------------------------------------------------------------
!     local parameters
!-----------------------------------------------------------------------
!

character(len=64), parameter    :: sub_name = 'coupler_main'
character(len=256), parameter   :: error_header =                               &
     '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):'
character(len=256), parameter   :: warn_header =                                &
     '==>Warning from ' // trim(mod_name) // '(' // trim(sub_name) // '):'
character(len=256), parameter   :: note_header =                                &
     '==>Note from ' // trim(mod_name) // '(' // trim(sub_name) // '):'

!#######################################################################

  call mpp_init()
!these clocks are on the global pelist
  initClock = mpp_clock_id( 'Initialization' )
  mainClock = mpp_clock_id( 'Main loop' )
  termClock = mpp_clock_id( 'Termination' )
  call mpp_clock_begin(initClock)
  
  call fms_init
  call constants_init

  call coupler_init
  if(do_chksum) call coupler_chksum('coupler_init+', 0)

  call mpp_set_current_pelist()

  call mpp_clock_end (initClock) !end initialization

  call mpp_clock_begin(mainClock) !begin main loop

!-----------------------------------------------------------------------
!------ ocean/slow-ice integration loop ------

     if(check_stocks >= 0) then
        call mpp_set_current_pelist()
        call flux_init_stocks(Time, Atm, Land, Ice, Ocean_state)
     endif

  do nc = 1, num_cpld_calls
     if( Atm%pe )then
        call mpp_set_current_pelist(Atm%pelist)
        call generate_sfc_xgrid( Land, Ice )
     end if
     call mpp_set_current_pelist()

     ! Calls to flux_ocean_to_ice and flux_ice_to_ocean are all PE communication
     ! points when running concurrently. The calls are placed next to each other in
     ! concurrent mode to avoid multiple synchronizations within the main loop.
     ! This is only possible in the serial case when use_lag_fluxes.
     call flux_ocean_to_ice( Time, Ocean, Ice, Ocean_ice_boundary )

     ! Update Ice_ocean_boundary; first iteration is supplied by restart     
     if( use_lag_fluxes )then
        call flux_ice_to_ocean( Time, Ice, Ocean, Ice_ocean_boundary )
     end if

     ! To print the value of frazil heat flux at the right time the following block
     ! needs to sit here rather than at the end of the coupler loop.
     if(check_stocks > 0) then
        if(check_stocks*((nc-1)/check_stocks) == nc-1 .AND. nc > 1) then
           call mpp_set_current_pelist()
           call flux_check_stocks(Time=Time, Atm=Atm, Lnd=Land, Ice=Ice, Ocn_state=Ocean_state)
        endif
     endif

     if( Atm%pe )then
        call mpp_set_current_pelist(Atm%pelist)
        if (do_ice) call update_ice_model_slow_up( Ocean_ice_boundary, Ice )

        !-----------------------------------------------------------------------
        !   ------ atmos/fast-land/fast-ice integration loop -------


        do na = 1, num_atmos_calls

           Time_atmos = Time_atmos + Time_step_atmos

           if (do_atmos) then
              call atmos_tracer_driver_gather_data(Atm%fields, Atm%tr_bot)
           endif

           if (do_flux) then
              !if(do_chksum) call coupler_chksum('sfc-', (nc-1)*num_atmos_calls+na)
              call sfc_boundary_layer( REAL(dt_atmos), Time_atmos, &
                   Atm, Land, Ice, Land_ice_atmos_boundary )
              !if(do_chksum) call coupler_chksum('sfc+', (nc-1)*num_atmos_calls+na)
           end if

           !      ---- atmosphere down ----

           if (do_atmos) &
                call update_atmos_model_down( Land_ice_atmos_boundary, Atm )

           call flux_down_from_atmos( Time_atmos, Atm, Land, Ice, &
                Land_ice_atmos_boundary, &
                Atmos_land_boundary, &
                Atmos_ice_boundary )


           !      --------------------------------------------------------------

           !      ---- land model ----

           if (do_land) &
                call update_land_model_fast( Atmos_land_boundary, Land )

           !      ---- ice model ----
           if (do_ice) &
                call update_ice_model_fast( Atmos_ice_boundary, Ice )

           !      --------------------------------------------------------------
           !      ---- atmosphere up ----

           call flux_up_to_atmos( Time_atmos, Land, Ice, Land_ice_atmos_boundary, &
                & Atmos_land_boundary, Atmos_ice_boundary )

           if (do_atmos) &
                call update_atmos_model_up( Land_ice_atmos_boundary, Atm )

           !--------------

        enddo

        !   ------ end of atmospheric time step loop -----
        if (do_land) call update_land_model_slow(Atmos_land_boundary,Land)
        !-----------------------------------------------------------------------

        !
        !     need flux call to put runoff and p_surf on ice grid
        !
        call flux_land_to_ice( Time, Land, Ice, Land_ice_boundary )

        Atmos_ice_boundary%p = 0.0 ! call flux_atmos_to_ice_slow ?

        !   ------ slow-ice model ------

        if (do_ice) then 
           call update_ice_model_slow_dn( Atmos_ice_boundary, &
                & Land_ice_boundary, Ice )
           call flux_ice_to_ocean_stocks(Ice)
        endif
        Time = Time_atmos
     end if                     !Atm%pe block

     if( .NOT.use_lag_fluxes )then !this will serialize
        call mpp_set_current_pelist()
        call flux_ice_to_ocean( Time, Ice, Ocean, Ice_ocean_boundary )
     end if

     if( Ocean%is_ocean_pe )then
        call mpp_set_current_pelist(Ocean%pelist)

        ! update_ocean_model since fluxes don't change here

        if (do_ocean) &
          call update_ocean_model( Ice_ocean_boundary, Ocean_state,  Ocean, &
                                   Time_ocean, Time_step_cpld )

        ! Get stocks from "Ice_ocean_boundary" and add them to Ocean stocks.
        ! This call is just for record keeping of stocks transfer and
        ! does not modify either Ocean or Ice_ocean_boundary
        call flux_ocean_from_ice_stocks(Ocean_state, Ocean, Ice_ocean_boundary)

        Time_ocean = Time_ocean +  Time_step_cpld

        !-----------------------------------------------------------------------
        Time = Time_ocean

     end if

     !--- write out intermediate restart file when needed.
     if( Time >= Time_restart ) then
        Time_restart_current = Time
        Time_restart = increment_date(Time, restart_interval(1), restart_interval(2), &
             restart_interval(3), restart_interval(4), restart_interval(5), restart_interval(6) )
        timestamp = date_to_string(time_restart_current)
        outunit= stdout()
        write(outunit,*) '=> NOTE from program coupler: intermediate restart file is written and ', &
             trim(timestamp),' is appended as prefix to each restart file name'
        if( Atm%pe )then        
           call atmos_model_restart(Atm, timestamp)
           call land_model_restart(timestamp)
           call ice_model_restart(timestamp)
        endif
        if( Ocean%is_ocean_pe) then
           call ocean_model_restart(Ocean_state, timestamp)
        endif
        call coupler_restart(Time, Time_restart_current, timestamp)
     end if

     !--------------
     if(do_chksum) call coupler_chksum('MAIN_LOOP+', nc)
     write( text,'(a,i4)' )'Main loop at coupling timestep=', nc
     call print_memuse_stats(text)


  enddo

  if(check_stocks >= 0) then
     call mpp_set_current_pelist()
     call flux_check_stocks(Time=Time, Atm=Atm, Lnd=Land, Ice=Ice, Ocn_state=Ocean_state)
  endif

! Need final update of Ice_ocean_boundary for concurrent restart
!  if( concurrent )then
!      call mpp_set_current_pelist()
!      call flux_ice_to_ocean( Time, Ice, Ocean, Ice_ocean_boundary )
!  endif

  call mpp_set_current_pelist()
!-----------------------------------------------------------------------
  call mpp_clock_end(mainClock)
  call mpp_clock_begin(termClock)

  if(do_chksum) call coupler_chksum('coupler_end-', nc)
  call coupler_end

  call mpp_clock_end(termClock)

  call print_memuse_stats( 'Memory HiWaterMark', always=.TRUE. )
  call fms_end

!-----------------------------------------------------------------------

contains

!#######################################################################

  subroutine coupler_init

    use ensemble_manager_mod, only : ensemble_manager_init, get_ensemble_id,ensemble_pelist_setup
    use ensemble_manager_mod, only : get_ensemble_size, get_ensemble_pelist

!-----------------------------------------------------------------------
!   initialize all defined exchange grids and all boundary maps
!-----------------------------------------------------------------------
 
!
!-----------------------------------------------------------------------
!     local parameters
!-----------------------------------------------------------------------
!

    character(len=64), parameter    :: sub_name = 'coupler_init'
    character(len=256), parameter   :: error_header =                               &
         '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):'
    character(len=256), parameter   :: warn_header =                                &
         '==>Warning from ' // trim(mod_name) // '(' // trim(sub_name) // '):'
    character(len=256), parameter   :: note_header =                                &
         '==>Note from ' // trim(mod_name) // '(' // trim(sub_name) // '):'

    integer :: unit,  ierr, io,    m, i, outunit, logunit
    integer :: date(6)
    type (time_type) :: Run_length
    character(len=9) :: month
    integer :: pe, npes

    integer :: ens_siz(4), ensemble_size

    integer :: atmos_pe_start=0, atmos_pe_end=0, &
               ocean_pe_start=0, ocean_pe_end=0
    integer :: n
    integer :: diag_model_subset=DIAG_ALL
    logical :: other_fields_exist
    logical, allocatable :: maskmap(:,:)
    character(len=256) :: err_msg
    integer :: date_restart(6)
    character(len=64)  :: filename, fieldname
    integer :: id_restart, l
!-----------------------------------------------------------------------

!----- write version to logfile -------
    call write_version_number(version, tag)

!----- read namelist -------

    unit = open_namelist_file()
    ierr=1; do while (ierr /= 0)
       read  (unit, nml=coupler_nml, iostat=io, end=10)
       ierr = check_nml_error (io, 'coupler_nml')
    enddo
10  call mpp_close(unit)
    outunit = stdout()
    logunit = stdlog()
    
!---- when concurrent is set true and mpp_io_nml io_clock_on is set true, the model
!---- will crash with error message "MPP_CLOCK_BEGIN: cannot change pelist context of a clock",
!---- so need to make sure it will not happen
    if(concurrent) then
       if(mpp_io_clock_on()) then
          call error_mesg ('program coupler', 'when coupler_nml variable concurrent is set to true, '// &
              'mpp_io_nml variable io_clock_non can not be set to true.', FATAL )
       endif
    endif
!----- read date and calendar type from restart file -----

    if( file_exist('INPUT/coupler.res') )then
!Balaji: currently written in binary, needs form=MPP_NATIVE
        call mpp_open( unit, 'INPUT/coupler.res', action=MPP_RDONLY )
        read( unit,*,err=999 )calendar_type
        read( unit,* )date_init
        read( unit,* )date
        goto 998 !back to fortran-4
!read old-style coupler.res
999     call mpp_close(unit)
        call mpp_open( unit, 'INPUT/coupler.res', action=MPP_RDONLY, form=MPP_NATIVE )
        read(unit)calendar_type
        read(unit)date
998     call mpp_close(unit)
    else
        force_date_from_namelist = .true.
    endif

!----- use namelist value (either no restart or override flag on) ---

    if ( force_date_from_namelist ) then

        if ( sum(current_date) <= 0 ) then
            call error_mesg ('program coupler',  &
                 'no namelist value for base_date or current_date', FATAL)
        else
            date      = current_date
        endif

!----- override calendar type with namelist value -----

        select case( uppercase(trim(calendar)) )
        case( 'JULIAN' )
            calendar_type = JULIAN
        case( 'NOLEAP' )
            calendar_type = NOLEAP
        case( 'THIRTY_DAY' )
            calendar_type = THIRTY_DAY_MONTHS
        case( 'NO_CALENDAR' )
            calendar_type = NO_CALENDAR
        end select

    endif

    call set_calendar_type (calendar_type, err_msg)
    if(err_msg /= '') then
      call mpp_error(FATAL, 'ERROR in coupler_init: '//trim(err_msg))
    endif

    if( concurrent .AND. .NOT.use_lag_fluxes )call mpp_error( WARNING, &
            'coupler_init: you have set concurrent=TRUE and use_lag_fluxes=FALSE &
            & in coupler_nml. When not using lag fluxes, components &
            & will synchronize at two points, and thus run serially.' )


    !Check with the ensemble_manager module for the size of ensemble
    !and PE counts for each member of the ensemble.
    !
    !NOTE: ensemble_manager_init renames all the output files (restart and diagnostics)
    !      to show which ensemble member they are coming from.
    !      There also need to be restart files for each member of the ensemble in INPUT.
    !
    !NOTE: if the ensemble_size=1 the input/output files will not be renamed.
    !
    
    call ensemble_manager_init() ! init pelists for ensembles
    ens_siz = get_ensemble_size()   
    ensemble_size = ens_siz(1)      
    npes = ens_siz(2)              

    !Check for the consistency of PE counts
    if( concurrent )then
!atmos_npes + ocean_npes must equal npes
        if( atmos_npes.EQ.0 )atmos_npes = npes - ocean_npes
        if( ocean_npes.EQ.0 )ocean_npes = npes - atmos_npes
!both must now be non-zero
        if( atmos_npes.EQ.0 .OR. ocean_npes.EQ.0 ) &
             call mpp_error( FATAL, 'coupler_init: atmos_npes or ocean_npes must be specified for concurrent coupling.' )
        if( atmos_npes+ocean_npes.NE.npes ) &
             call mpp_error( FATAL, 'coupler_init: atmos_npes+ocean_npes must equal npes for concurrent coupling.' )
    else                        !serial timestepping
        if( atmos_npes.EQ.0 )atmos_npes = npes
        if( ocean_npes.EQ.0 )ocean_npes = npes
        if( max(atmos_npes,ocean_npes).EQ.npes )then !overlapping pelists
            ! do nothing
        else                    !disjoint pelists
            if( atmos_npes+ocean_npes.NE.npes ) call mpp_error( FATAL,  &
                 'coupler_init: atmos_npes+ocean_npes must equal npes for serial coupling on disjoint pelists.' )
        end if
    end if    

    allocate( Atm%pelist  (atmos_npes) )
    allocate( Ocean%pelist(ocean_npes) )

    !Set up and declare all the needed pelists
    call ensemble_pelist_setup(concurrent, atmos_npes, ocean_npes, Atm%pelist, Ocean%pelist)
    ensemble_id = get_ensemble_id() 
 
    allocate(ensemble_pelist(1:ensemble_size,1:npes))   
    call get_ensemble_pelist(ensemble_pelist) 

    Atm%pe            = ANY(Atm%pelist   .EQ. mpp_pe()) 
    Ocean%is_ocean_pe = ANY(Ocean%pelist .EQ. mpp_pe())  
    Ice%pe  = Atm%pe
    Land%pe = Atm%pe
 
    !Why is the following needed?
    if( Atm%pe )            call mpp_set_current_pelist( Atm%pelist   )
    if( Ocean%is_ocean_pe ) call mpp_set_current_pelist( Ocean%pelist )
    
    !Write out messages on root PEs
    if(mpp_pe().EQ.mpp_root_pe() )then
       write( text,'(a,2i6,a,i2.2)' )'Atmos PE range: ', Atm%pelist(1)  , Atm%pelist(atmos_npes)  ,&
            ' ens_', ensemble_id
       call mpp_error( NOTE, 'coupler_init: '//trim(text) )
       write( text,'(a,2i6,a,i2.2)' )'Ocean PE range: ', Ocean%pelist(1), Ocean%pelist(ocean_npes), &
            ' ens_', ensemble_id
       call mpp_error( NOTE, 'coupler_init: '//trim(text) )
       if( concurrent )then
          call mpp_error( NOTE, 'coupler_init: Running with CONCURRENT coupling.' )

          write( logunit,'(a)' )'Using concurrent coupling...'
          write( logunit,'(a,4i4)' ) &
               'atmos_pe_start, atmos_pe_end, ocean_pe_start, ocean_pe_end=', &
               Atm%pelist(1)  , Atm%pelist(atmos_npes), Ocean%pelist(1), Ocean%pelist(ocean_npes) 
       else
          call mpp_error( NOTE, 'coupler_init: Running with SERIAL coupling.' )
       end if
       if( use_lag_fluxes )then
          call mpp_error( NOTE, 'coupler_init: Sending LAG fluxes to ocean.' )
       else
          call mpp_error( NOTE, 'coupler_init: Sending most recent fluxes to ocean.' )
       end if
    endif

    if( ice_npes.NE.0 ) &
         call mpp_error( WARNING, 'coupler_init: pelists not yet implemented for ice.' )
    if( land_npes.NE.0 ) &
         call mpp_error( WARNING, 'coupler_init: pelists not yet implemented for land.' )

!----- write namelist to logfile -----
    if( mpp_pe() == mpp_root_pe() )write( logunit, nml=coupler_nml )

!----- write current/initial date actually used to logfile file -----

    if ( mpp_pe().EQ.mpp_root_pe() ) &
         write( logunit, 16 )date(1),trim(month_name(date(2))),date(3:6)
16  format ('  current date used = ',i4,1x,a,2i3,2(':',i2.2),' gmt') 

!----- check the value of layout and setup the maskmap for domain layout.
    if( n_mask > 0 ) then
       if(do_atmos .OR. do_land) call mpp_error(FATAL, &
            'program coupler: do_atmos and do_land should be false when n_mask > 0')

       if(concurrent) call mpp_error(FATAL, &
            'program coupler: can not run concurrent run when some regions are masked out')
       if( layout_mask(1)*layout_mask(2) - n_mask .NE. npes ) call mpp_error(FATAL, &
            'program coupler: layout(1)*layout(2) - n_mask should equal to npes when n_mask>0')
       call mpp_error(NOTE, 'program coupler: layout_mask and mask_list is set in coupler_nml, ' // &
                            'the value of layout_mask will override the layout specified in each component model')

       allocate(maskmap(layout_mask(1), layout_mask(2)) )
       maskmap = .TRUE.
       do n=1, n_mask
          if (mask_list(1,n) .gt. layout_mask(1) ) &
             call mpp_error( FATAL, 'program coupler: mask_list elements outside layout defines.' )
          if (mask_list(2,n) .gt. layout_mask(2) ) &
             call mpp_error( FATAL, 'program coupler: mask_list elements outside layout defines.' )
          maskmap(mask_list(1,n),mask_list(2,n)) = .false.
       enddo
       !--- copy maskmap value to each model data type
       allocate(Atm%maskmap(layout_mask(1), layout_mask(2)), Land%maskmap(layout_mask(1), layout_mask(2)) )
       allocate(Ice%maskmap(layout_mask(1), layout_mask(2)), Ocean%maskmap(layout_mask(1), layout_mask(2)))
       Atm%maskmap = maskmap;  Land%maskmap = maskmap
       Ice%maskmap = maskmap;  Ocean%maskmap = maskmap
       deallocate(maskmap)
    else
       if( layout_mask(1)*layout_mask(2) .NE. 0 ) call mpp_error(NOTE, &
            'program coupler: when no region is masked out, layout_mask need not be set' )
    end if

!-----------------------------------------------------------------------
!------ initialize diagnostics manager ------

!jwd Fork here is somewhat dangerous. It relies on "no side effects" from
!    diag_manager_init. diag_manager_init or this section should be 
!    re-architected to guarantee this or remove this assumption.
!    For instance, what follows assumes that get_base_date has the same
!    time for both Atm and Ocean pes. While this should be the case, the
!    possible error condition needs to be checked

    if( Atm%pe )then
        call mpp_set_current_pelist(Atm%pelist)
        if(atmos_npes /= npes)diag_model_subset = DIAG_OTHER  ! change diag_model_subset from DIAG_ALL
    elseif( Ocean%is_ocean_pe )then  ! Error check above for disjoint pelists should catch any problem
        call mpp_set_current_pelist(Ocean%pelist)
        if(ocean_npes /= npes)diag_model_subset = DIAG_OCEAN  ! change diag_model_subset from DIAG_ALL
    end if
    call diag_manager_init(DIAG_MODEL_SUBSET=diag_model_subset)   ! initialize diag_manager for processor subset output
    call print_memuse_stats( 'diag_manager_init' )
!-----------------------------------------------------------------------
!------ reset pelist to "full group" ------

    call mpp_set_current_pelist()
!----- always override initial/base date with diag_manager value -----

    call get_base_date ( date_init(1), date_init(2), date_init(3), &
         date_init(4), date_init(5), date_init(6)  )

!----- use current date if no base date ------

    if ( date_init(1) == 0 ) date_init = date

!----- set initial and current time types ------

    Time_init = set_date (date_init(1), date_init(2), date_init(3), &
         date_init(4), date_init(5), date_init(6))

    Time      = set_date (date(1), date(2), date(3),  &
         date(4), date(5), date(6))

    Time_start = Time

!----- compute the ending time -----

    Time_end = Time
    do m=1,months
       Time_end = Time_end + set_time(0,days_in_month(Time_end))
    end do
    Time_end   = Time_end + set_time(hours*3600+minutes*60+seconds, days)
    Run_length = Time_end - Time

!--- get the time that last intermediate restart file was written out.
    if (file_exist('INPUT/coupler.intermediate.res')) then
       call mpp_open(unit,'INPUT/coupler.intermediate.res',action=MPP_RDONLY)
       read(unit,*) date_restart
       call mpp_close(unit)
    else
       date_restart = date
    endif

    Time_restart_current = Time
    if(ALL(restart_interval ==0)) then
       Time_restart = increment_date(Time_end, 1, 0, 0, 0, 0, 0)   ! no intermediate restart
    else
       Time_restart = set_date(date_restart(1), date_restart(2), date_restart(3),  &
                               date_restart(4), date_restart(5), date_restart(6) )
       Time_restart = increment_date(Time_restart, restart_interval(1), restart_interval(2), &
            restart_interval(3), restart_interval(4), restart_interval(5), restart_interval(6) )
       if(Time_restart <= Time) call mpp_error(FATAL, &
            '==>Error from program coupler: The first intermediate restart time is no larger than the start time')
    end if

!-----------------------------------------------------------------------
!----- write time stamps (for start time and end time) ------

    call mpp_open( unit, 'time_stamp.out', nohdrs=.TRUE. )

    month = month_name(date(2))
    if ( mpp_pe().EQ.mpp_root_pe() ) write (unit,20) date, month(1:3)

    call get_date (Time_end, date(1), date(2), date(3),  &
         date(4), date(5), date(6))
    month = month_name(date(2))
    if ( mpp_pe().EQ.mpp_root_pe() ) write (unit,20) date, month(1:3)

    call mpp_close(unit)

20  format (6i4,2x,a3)

!-----------------------------------------------------------------------
!----- compute the time steps ------

    Time_step_cpld  = set_time (dt_cpld ,0)
    Time_step_atmos = set_time (dt_atmos,0)

!----- determine maximum number of iterations per loop ------

    num_cpld_calls  = Run_length      / Time_step_cpld
    num_atmos_calls = Time_step_cpld  / Time_step_atmos

!-----------------------------------------------------------------------
!------------------- some error checks ---------------------------------

!----- initial time cannot be greater than current time -------

    if ( Time_init > Time ) call error_mesg ('program coupler',  &
         'initial time is greater than current time', FATAL)

!----- make sure run length is a multiple of ocean time step ------

    if ( num_cpld_calls * Time_step_cpld  /= Run_length )  &
         call error_mesg ('program coupler',  &
         'run length must be multiple of coupled time step', FATAL)

! ---- make sure cpld time step is a multiple of atmos time step ----

    if ( num_atmos_calls * Time_step_atmos /= Time_step_cpld )  &
         call error_mesg ('program coupler',   &
         'cpld time step is not a multiple of the atmos time step', FATAL)

!
!       Initialize the tracer manager. This needs to be done on all PEs,
!       before the individual models are initialized.
!

    call tracer_manager_init
!
!       Initialize the coupler types
!

    call coupler_types_init

!-----------------------------------------------------------------------
!------ initialize component models ------
!------ grid info now comes from grid_spec file

    if( Atm%pe )then
        call mpp_set_current_pelist(Atm%pelist)
!---- atmosphere ----
        call atmos_model_init( Atm, Time_init, Time, Time_step_atmos )
        call print_memuse_stats( 'atmos_model_init' )

!---- land ----------
        call land_model_init( Atmos_land_boundary, Land, Time_init, Time, &
             Time_step_atmos, Time_step_cpld )
        call print_memuse_stats( 'land_model_init' )

!---- ice -----------
        call ice_model_init( Ice, Time_init, Time, Time_step_atmos, Time_step_cpld )
        call print_memuse_stats( 'ice_model_init' )
        call data_override_init(Atm_domain_in = Atm%domain, Ice_domain_in = Ice%domain, Land_domain_in=Land%domain)
    end if
    if( Ocean%is_ocean_pe )then
        call mpp_set_current_pelist(Ocean%pelist)
!---- ocean ---------
        call ocean_model_init( Ocean, Ocean_state, Time_init, Time )
        call print_memuse_stats( 'ocean_model_init' )
        call data_override_init(Ocean_domain_in = Ocean%domain )
    end if
    call mpp_set_current_pelist(ensemble_pelist(ensemble_id,:))

    call mpp_broadcast_domain(Ice%domain)
    call mpp_broadcast_domain(Ocean%domain)
!-----------------------------------------------------------------------
!---- initialize flux exchange module ----
    call flux_exchange_init ( Time, Atm, Land, Ice, Ocean, Ocean_state,&
         atmos_ice_boundary, land_ice_atmos_boundary, &
         land_ice_boundary, ice_ocean_boundary, ocean_ice_boundary, &
         dt_atmos=dt_atmos, dt_cpld=dt_cpld)

    Time_atmos = Time
    Time_ocean = Time

!
!       read in extra fields for the air-sea gas fluxes
!

    if ( Atm%pe ) then
      call mpp_set_current_pelist(Atm%pelist)
      allocate(Ice_bc_restart(Ice%ocean_fluxes%num_bcs))
      allocate(ice_bc_restart_file(Ice%ocean_fluxes%num_bcs))
      do n = 1, Ice%ocean_fluxes%num_bcs  !{
        if(Ice%ocean_fluxes%bc(n)%num_fields .LE. 0) cycle
        filename = trim(Ice%ocean_fluxes%bc(n)%ice_restart_file)
        do l = 1, num_ice_bc_restart
           if(trim(filename) == ice_bc_restart_file(l)) exit
        end do
        if(l>num_ice_bc_restart) then
           num_ice_bc_restart = num_ice_bc_restart + 1
           ice_bc_restart_file(l) = trim(filename)
        end if
        filename = 'INPUT/'//trim(filename)
        other_fields_exist = .false.
        do m = 1, Ice%ocean_fluxes%bc(n)%num_fields  !{
          fieldname = trim(Ice%ocean_fluxes%bc(n)%field(m)%name)
          id_restart = register_restart_field(Ice_bc_restart(l), ice_bc_restart_file(l), &
                       fieldname, Ice%ocean_fluxes%bc(n)%field(m)%values, Ice%domain    )
          if (field_exist(filename, fieldname, Ice%domain) ) then
            other_fields_exist = .true.
            write (outunit,*) trim(note_header), ' Reading restart info for ',         &
                 trim(fieldname), ' from ',  trim(filename)
            call read_data(filename, fieldname, Ice%ocean_fluxes%bc(n)%field(m)%values, Ice%domain)
          elseif (other_fields_exist) then
            call mpp_error(FATAL, trim(error_header) // ' Couldn''t find field ' //     &
                 trim(fieldname) // ' in file ' //trim(filename))
          endif
        enddo  !} m
      enddo  !} n
    endif
    if ( Ocean%is_ocean_pe ) then
      call mpp_set_current_pelist(Ocean%pelist)
      allocate(Ocn_bc_restart(Ocean%fields%num_bcs))
      allocate(ocn_bc_restart_file(Ocean%fields%num_bcs))
      do n = 1, Ocean%fields%num_bcs  !{
        if(Ocean%fields%bc(n)%num_fields .LE. 0) cycle
        filename = trim(Ocean%fields%bc(n)%ocean_restart_file)
        do l = 1, num_ocn_bc_restart
           if(trim(filename) == ocn_bc_restart_file(l)) exit
        end do
        if(l>num_ocn_bc_restart) then
           num_ocn_bc_restart = num_ocn_bc_restart + 1
           ocn_bc_restart_file(l) = trim(filename)
        end if
        filename = 'INPUT/'//trim(filename)
        other_fields_exist = .false.
        do m = 1, Ocean%fields%bc(n)%num_fields  !{
          fieldname = trim(Ocean%fields%bc(n)%field(m)%name)
          id_restart = register_restart_field(Ocn_bc_restart(l), Ocn_bc_restart_file(l), &
                       fieldname, Ocean%fields%bc(n)%field(m)%values, Ocean%domain    )
          if (field_exist(filename, fieldname, Ocean%domain) ) then
            other_fields_exist = .true.
            write (outunit,*) trim(note_header), ' Reading restart info for ',         &
                 trim(fieldname), ' from ', trim(filename)
            call read_data(filename, fieldname, Ocean%fields%bc(n)%field(m)%values, Ocean%domain)
          elseif (other_fields_exist) then
            call mpp_error(FATAL, trim(error_header) // ' Couldn''t find field ' //     &
                 trim(fieldname) // ' in file ' //trim(filename))
          endif
        enddo  !} m
      enddo  !} n
    endif

    call mpp_set_current_pelist()

!-----------------------------------------------------------------------
!---- open and close dummy file in restart dir to check if dir exists --

    call mpp_open( unit, 'RESTART/file' )
    call mpp_close(unit, MPP_DELETE)

    ! Call to daig_grid_end to free up memory used during regional
    ! output setup
    CALL diag_grid_end()

!-----------------------------------------------------------------------
    call print_memuse_stats('coupler_init')
  end subroutine coupler_init

!#######################################################################

  subroutine coupler_end

!-----------------------------------------------------------------------

    call mpp_set_current_pelist()

!----- check time versus expected ending time ----

    if (Time /= Time_end) call error_mesg ('program coupler',  &
         'final time does not match expected ending time', WARNING)

!-----------------------------------------------------------------------
!the call to fms_io_exit has been moved here
!this will work for serial code or concurrent (disjoint pelists)
!but will fail on overlapping but unequal pelists
    if( Ocean%is_ocean_pe )then
        call mpp_set_current_pelist(Ocean%pelist)
        call ocean_model_end (Ocean, Ocean_state, Time)
    end if
    if( Atm%pe )then
        call mpp_set_current_pelist(Atm%pelist)
        call atmos_model_end (Atm)
        call  land_model_end (Atmos_land_boundary, Land)
        call   ice_model_end (Ice)
    end if

    !----- write restart file ------
    call coupler_restart(Time, Time_restart_current)

    call fms_io_exit
    call diag_manager_end (Time)
    call mpp_set_current_pelist()

!-----------------------------------------------------------------------

  end subroutine coupler_end

  !--- writing restart file that contains running time and restart file writing time.
  subroutine coupler_restart(Time_run, Time_res, time_stamp)
    type(time_type),   intent(in)           :: Time_run, Time_res
    character(len=*), intent(in),  optional :: time_stamp
    character(len=128)                      :: file_run, file_res
    integer :: yr, mon, day, hr, min, sec, date(6), unit

    call mpp_set_current_pelist()

    ! write restart file
    if(present(time_stamp)) then
       file_run = 'RESTART/'//trim(time_stamp)//'.coupler.res'
       file_res = 'RESTART/'//trim(time_stamp)//'.coupler.intermediate.res'
    else
       file_run = 'RESTART/coupler.res'
       file_res = 'RESTART/coupler.intermediate.res'
    endif

    !----- compute current date ------
    call get_date (Time_run, date(1), date(2), date(3),  &
                   date(4), date(5), date(6))
    call mpp_open( unit, file_run, nohdrs=.TRUE. )
    if ( mpp_pe().EQ.mpp_root_pe() )then
        write( unit, '(i6,8x,a)' )calendar_type, &
             '(Calendar: no_calendar=0, thirty_day_months=1, julian=2, gregorian=3, noleap=4)'

        write( unit, '(6i6,8x,a)' )date_init, &
             'Model start time:   year, month, day, hour, minute, second'
        write( unit, '(6i6,8x,a)' )date, &
             'Current model time: year, month, day, hour, minute, second'
    end if
    call mpp_close(unit)

    if(Time_res > Time_start) then
       call mpp_open( unit, file_res, nohdrs=.TRUE. )
       if ( mpp_pe().EQ.mpp_root_pe() )then
          call get_date(Time_res ,yr,mon,day,hr,min,sec)
          write( unit, '(6i6,8x,a)' )yr,mon,day,hr,min,sec, &
               'Current intermediate restart time: year, month, day, hour, minute, second'
       end if
       call mpp_close(unit)     
    end if

    if( Ocean%is_ocean_pe )then
        call mpp_set_current_pelist(Ocean%pelist)
        do n = 1, num_ocn_bc_restart
           call save_restart(Ocn_bc_restart(n), time_stamp)
        enddo
    endif
    if( Atm%pe )then
        call mpp_set_current_pelist(Atm%pelist)
        do n = 1, num_ice_bc_restart
           call save_restart(Ice_bc_restart(n), time_stamp)
        enddo
    endif


  end subroutine coupler_restart

!--------------------------------------------------------------------------

  subroutine coupler_chksum(id, timestep)

    character(len=*), intent(in) :: id
    integer         , intent(in) :: timestep

    type :: tracer_ind_type
       integer :: atm, ice, lnd ! indices of the tracer in the respective models
    end type tracer_ind_type
    integer                            :: n_atm_tr, n_lnd_tr, n_exch_tr
    integer                            :: n_atm_tr_tot, n_lnd_tr_tot
    integer                            :: i, tr, n, m, outunit
    type(tracer_ind_type), allocatable :: tr_table(:)
    character(32) :: tr_name

    call get_number_tracers (MODEL_ATMOS, num_tracers=n_atm_tr_tot, &
                             num_prog=n_atm_tr)
    call get_number_tracers (MODEL_LAND, num_tracers=n_lnd_tr_tot, &
                             num_prog=n_lnd_tr)

    ! assemble the table of tracer number translation by matching names of
    ! prognostic tracers in the atmosphere and surface models; skip all atmos.
    ! tracers that have no corresponding surface tracers.
    allocate(tr_table(n_atm_tr))
    n = 1
    do i = 1,n_atm_tr
       call get_tracer_names( MODEL_ATMOS, i, tr_name )
       tr_table(n)%atm = i
       tr_table(n)%ice = get_tracer_index ( MODEL_ICE,  tr_name )
       tr_table(n)%lnd = get_tracer_index ( MODEL_LAND, tr_name )
       if(tr_table(n)%ice/=NO_TRACER.or.tr_table(n)%lnd/=NO_TRACER) &
            n = n+1
    enddo
    n_exch_tr = n-1

100 FORMAT("CHECKSUM::",A32," = ",Z20)
101 FORMAT("CHECKSUM::",A16,a,'%',a," = ",Z20)

    if( Atm%pe )then
       call mpp_set_current_pelist(Atm%pelist)

       outunit = stdout()
       write(outunit,*) 'BEGIN CHECKSUM(Atm):: ', id, timestep
       write(outunit,100) 'atm%t_bot', mpp_chksum(atm%t_bot)
       write(outunit,100) 'atm%z_bot', mpp_chksum(atm%z_bot)
       write(outunit,100) 'atm%p_bot', mpp_chksum(atm%p_bot)
       write(outunit,100) 'atm%u_bot', mpp_chksum(atm%u_bot)
       write(outunit,100) 'atm%v_bot', mpp_chksum(atm%v_bot)
       write(outunit,100) 'atm%p_surf', mpp_chksum(atm%p_surf)
       write(outunit,100) 'atm%gust', mpp_chksum(atm%gust)
       do tr = 1,n_exch_tr
          n = tr_table(tr)%atm
          if(n /= NO_TRACER ) then
             call get_tracer_names( MODEL_ATMOS, tr_table(tr)%atm, tr_name )
             write(outunit,100) 'atm%'//trim(tr_name), mpp_chksum(Atm%tr_bot(:,:,n))
          endif
       enddo
   
       write(outunit,100) 'land%t_surf', mpp_chksum(land%t_surf)
       write(outunit,100) 'land%t_ca', mpp_chksum(land%t_ca)
       write(outunit,100) 'land%rough_mom', mpp_chksum(land%rough_mom)
       write(outunit,100) 'land%rough_heat', mpp_chksum(land%rough_heat)
       write(outunit,100) 'land%rough_scale', mpp_chksum(land%rough_scale)
       do tr = 1,n_exch_tr
          n = tr_table(tr)%lnd
          if(n /= NO_TRACER ) then
             call get_tracer_names( MODEL_ATMOS, tr_table(tr)%atm, tr_name )
             write(outunit,100) 'land%'//trim(tr_name), mpp_chksum(Land%tr(:,:,:,n))
          endif
       enddo
   
       write(outunit,100) 'ice%t_surf', mpp_chksum(ice%t_surf)
       write(outunit,100) 'ice%rough_mom', mpp_chksum(ice%rough_mom)
       write(outunit,100) 'ice%rough_heat', mpp_chksum(ice%rough_heat)
       write(outunit,100) 'ice%rough_moist', mpp_chksum(ice%rough_moist)
       write(outunit,*) 'STOP CHECKSUM(Atm):: ', id, timestep
   
    !endif

    !if( Ocean%is_ocean_pe )then
        !call mpp_set_current_pelist(Ocean%pelist)

       write(outunit,*) 'BEGIN CHECKSUM(Ice):: ', id, timestep
       do n = 1, ice%ocean_fields%num_bcs  !{
          do m = 1, ice%ocean_fields%bc(n)%num_fields  !{
             !write(outunit,101) 'ice%', m, n, mpp_chksum(Ice%ocean_fields%bc(n)%field(m)%values)
             write(outunit,101) 'ice%',trim(ice%ocean_fields%bc(n)%name), &
                  trim(ice%ocean_fields%bc(n)%field(m)%name), mpp_chksum(Ice%ocean_fields%bc(n)%field(m)%values)
          enddo  !} m
       enddo  !} n
       write(outunit,*) 'STOP CHECKSUM(Ice):: ', id, timestep

    endif

    deallocate(tr_table)

    call mpp_set_current_pelist()

  end subroutine coupler_chksum

  !#######################################################################

  end program coupler_main