!<CONTACT EMAIL="petri@pik-potsdam.de"> Stefan Petri
!</CONTACT>
!<REVIEWER EMAIL="">
!
!</REVIEWER>
!-----------------------------------------------------------------------
!<OVERVIEW>
!  Driver for the atmospheric model Aeolus2, contains routines to advance the
!  atmospheric model state by one time step.
!</OVERVIEW>

!<DESCRIPTION> This version of atmos_model_mod is going to implement
!     a coupling interface to the atmosphere model Aeolus2, which is
!     basedon a pseudo-spectral moist-convective Thermal Rotating Shallow
!     Water model (mcTRSW). 
!     It requires two routines to advance
!     the atmospheric model one time step into the future, up and and down.

!     The boundary variables needed by other component models for coupling
!     are contained in a derived data type. A variable of this derived type
!     is returned when initializing the atmospheric model. It is used by other
!     routines in this module and by coupling routines. The contents of
!     this derived type should only be modified by the atmospheric model.

!</DESCRIPTION>

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                   !!
!!                   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                      !!
!!                                                                   !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#include <fms_platform.h>

module atmosphere_mod

  use time_manager_mod, only: time_type, get_time, set_time, get_date, set_date, operator(+), &
       operator(/), &
       days_in_year, get_calendar_type, THIRTY_DAY_MONTHS, NOLEAP

  use fms_mod, only:                      &
       file_exist, open_namelist_file,    &
       error_mesg, FATAL,                 &
       check_nml_error, stdlog, stdout,   &
       write_version_number,              &
       mpp_pe, mpp_root_pe,               &
       close_file, set_domain,            &
       mpp_clock_id, mpp_clock_begin,     &
       mpp_clock_end, CLOCK_SUBCOMPONENT, &
       clock_flag_default,                &
       mpp_npes

  !use fms_io_mod, only:                   &
  !     return_domain
  
  use mpp_domains_mod, only: domain1D, domain2D,        &
       mpp_domains_init, mpp_define_domains,            &
       mpp_get_compute_domain, mpp_get_compute_domains, &
       mpp_get_domain_components, mpp_get_pelist
  use  field_manager_mod, only: MODEL_ATMOS
  use      constants_mod, only: omega, cp_air, rdgas, grav, rvgas, kappa, radius, pstd_mks, PI, stefan
  use stock_constants_mod, only : ISTOCK_WATER, ISTOCK_HEAT, ISTOCK_SALT

  use diag_manager_mod, only: diag_axis_init, register_diag_field
  use diag_manager_mod, only: register_static_field, send_data
  use diag_manager_mod, only: diag_manager_end

  use mpp_mod, only: mpp_error, mpp_get_current_pelist, mpp_sum
  use tracer_manager_mod, only: get_number_tracers, get_tracer_index, NO_TRACER
  use xgrid_mod, only: grid_box_type
  use grid_mod, only: get_grid_cell_area, get_grid_cell_centers

  use topography_mod, only: get_topog_mean, get_topog_stdev

  ! we now have two alternate ways to get the cosine of the zenith angle
  use astronomy_mod, only: astronomy_init, daily_mean_solar !, annual_mean_solar, get_period
  

  !! To enable strict floating point exceptions inside this module.
  !! Taken from the example in the Intel Fortran Compiler Reference Manual.
  !! Most probably will not work with any other compiler
  use ifcore

  use aeolus2fms_mod, only: aeolus2fms_interface_check
  use aeolus2fms_mod, only: aeolus2fms_init_grid, aeolus2fms_get_axes_coords, aeolus2fms_init
  use aeolus2fms_mod, only: aeolus2fms_atmosphere_down, aeolus2fms_atmosphere_up
  use aeolus2fms_mod, only: aeolus2fms_finish, aeolus2fms_restart
  !use aeolus2fms_mod, only: aeolus2fms_get_boundary
  use aeolus2fms_mod, only: aeolus2fms_get_bottom_mass, aeolus2fms_get_bottom_wind, aeolus2fms_get_stock_pe

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

  implicit none
  private

  public atmosphere_down,     atmosphere_up,       &
       atmosphere_init,       atmosphere_end,      &
       atmosphere_resolution, atmosphere_boundary, &
       atmosphere_cell_area,  atmosphere_restart,  &
       get_atmosphere_axes,   atmosphere_domain

  public  get_bottom_mass,  get_bottom_wind,  get_stock_pe

  public  surf_diff_type

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

  !<PUBLICTYPE>
  type surf_diff_type
     real, pointer, dimension(:,:) :: dtmass  => NULL() ! dt/mass, where dt = atmospheric time step (sec)
                                                        ! mass = mass per unit area of lowest atmospheric layer (Kg/m2)
     real, pointer, dimension(:,:) :: dflux_t => NULL() ! derivative of implicit part of downward temperature flux at the top of the lowest
                                                        ! atmospheric layer with respect to the temperature
                                                        ! of that layer  (J/(m2 K))
     real, pointer, dimension(:,:,:) :: dflux_tr => NULL() ! derivative of implicit part of downward flux of specific humidity
                                                        ! at the top of the lowest atmospheric layer with respect to
                                                        ! the specific humidity of that layer  (kg/(m2 s))
     real, pointer, dimension(:,:) :: delta_t => NULL() ! the increment in temperature in the lowest atmospheric
                                                        ! layer (((i+1)-(i-1) if atmos model is leapfrog) (K)
                                                        ! (defined in  gcm_vert_diff_down as the increment computed up
                                                        ! to this point in model, including effect of vertical
                                                        ! diffusive flux at top of lowest model level, presumed
                                                        ! to be modified to include effects of surface fluxes
                                                        ! outside of this module, then used to start upward
                                                        ! tridiagonal sweep,
     real, pointer, dimension(:,:,:):: delta_tr => NULL() ! similarly for the increment in specific humidity
                                                        ! (non-dimensional  = Kg/Kg)
     real, pointer, dimension(:,:) :: delta_u => NULL()
     real, pointer, dimension(:,:) :: delta_v => NULL()
     real, pointer, dimension(:,:) :: sst_miz => NULL() ! only used if /flux_exchange_nml/do_forecast is .true.
  end type surf_diff_type
  !</PUBLICTYPE>

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

  character(len=128) :: version = '$Id: atmosphere.F90 3862 2021-12-06 15:51:04Z petri $'
  character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $'

  !---- namelist (saved in file input.nml) ----
  !<NAMELIST NAME="atmos_aeolus2_nml">
  ! <DATA NAME="NoNx" TYPE="integer">
  !  Longitudinal resolution of atmosphere grid (at the equator).
  !  Must match the data in the grid specification file.
  ! </DATA>
  ! <DATA NAME="NoNy" TYPE="integer">
  !  Latitudinal resolution of atmosphere grid.
  !  Must match the data in the grid specification file.
  ! </DATA>
  !! <DATA NAME="NoNz" TYPE="integer">
  !!  Number of vertical layers of atmosphere grid in the Troposphere
  !! </DATA>
  !! <DATA NAME="NoNz_Strato" TYPE="integer">
  !!  Number of vertical layers of atmosphere grid in the Stratosphere
  !! </DATA>
  ! <DATA NAME="without_topo" TYPE="logical">
  !  If true, use flat topography in Aeolus2, else use real topography heights
  ! </DATA>
  ! <DATA NAME="update_land_frac_always" TYPE="logical">
  !  If true, update Aeolus2 topo.land_frac in every _down and _up step.
  !  Also update the surface type fractions.
  !  If false, do it only in the very first update step.
  !  To update in every step might be useful if the land or ocean might decide to change the
  !  land sea mask. With current implementations, that seems never to be the case.
  ! </DATA>
  ! <DATA NAME="fix_astro" TYPE="logical" DEFAULT=".false.">
  !  Freeze time for obtaining astronomical data.
  !  This can be useful to study spin-up behaviour / stability / equilibrium without interference of seasonal (or diurnal) cycle.
  ! </DATA>
  ! <DATA NAME="astro_date" TYPE="integer, dimension(6)"  DEFAULT="0">
  !  If fix_astro is true, this specifies the time to use
  ! </DATA>
  ! <DATA NAME="solar_constant" TYPE="real" UNIT="W/m2" DEFAULT="1365.0">
  !  Solar constant
  ! </DATA>
  !</NAMELIST>
  namelist /atmos_aeolus2_nml/ &
       NoNx, NoNy,          & ! Number of grid cells in lon, lat
       !NoNz, NoNz_Strato,  & ! Number of grid cells in height (Tropo+Strato)
       lon_0_360,           & ! whether longitudes run 0 to 360 in the grid spec file
       without_topo,        & ! choose between flat and real topography
       update_land_frac_always, & ! wether to update land_frac (and with it the surface type fractions)
                                  ! in every _down and _up step, or only in the first one
       fix_astro,          & ! freeze astronomical time, i.e. keep time for obtaining astronomical data constant
       astro_date,         & ! if fix_astro, this gives the time to use
       solar_constant

       

  !-----------------------------------------------------------------------
  !---- private data ----

  integer :: NoNx, NoNy !, NoNz, NoNz_Strato   ! values from namelist
  logical :: lon_0_360 = .false. ! .true. was the default for old Aeolus
  logical :: without_topo = .false.
  logical :: update_land_frac_always = .false.

  type    (time_type) :: Time_step_atmos
  real                :: dt_atmos        ! Time_step_atmos converted to seconds
  integer, dimension(4)              :: atmos_axes ! dimension(4) is required by atmos_coupled_mod

  logical             :: fix_astro = .false. ! freeze astronomical time, i.e. keep time for solar data constant
  integer, dimension(6) :: astro_date = (/ 0, 0, 0, 0, 0, 0 /)          ! year, month, day, hour, minute, second
  type    (time_type) :: Time_fix_astro             ! if fix_astro, this gives the time to use
  integer             :: fixdayoftheyear = -1       ! if fix_astro, this the day of the year to use for sinsol

  real :: solar_constant = 1365.0

  ! clock ids for timing
  integer :: clock_id_aeolus2
  integer :: clockid_init, clockid_up, clockid_down, clockid_diag

  logical :: module_is_initialized =.false.

  integer :: nr_tracers = 0      ! total number of tracers in surf_diff_type fields
  integer :: q_ind = NO_TRACER   ! index of humidity in tracers array
  integer :: co2_ind = NO_TRACER ! index of CO2 in tracers array

  type(domain2D) :: grid_domain
  integer :: is, ie, js, je  ! local domain boundary indices

  real, dimension(:), allocatable :: lonb, latb, lon, lat ! , level, levelb ! for the _global_ axes

  ! diag var ids
  integer :: id_t = -1, id_astro_cosz = -1, id_astro_solar = -1
  integer :: atmos_stock_water_id = -1, atmos_stock_heat_id = -1
  integer :: id_flux_lw_sf_up = -1 ! to diagnose problems with interpolation to reduced grid

  ! for the FMS astronomy_mod
  integer :: DaysPerYear
  real, dimension(:), allocatable :: rad_lat ! grid-cell midpoint coordinates in meridional direction, in radians

  !------------------------- axis names ----------------------------------
  character(len=8) :: axiset   = 'aeolus2'
  character(len=8) :: mod_name = 'aeolus2'

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

  !! To enable strict floating point exceptions inside this module.
  !! Taken from the example in the Intel Fortran Compiler Reference Manual.
  !! Most probably will not work with any other compiler
  integer(4) :: original_fpe_flags
  integer(4) :: strict_fpe_flags = FPE_M_TRAP_UND + FPE_M_TRAP_OVF + FPE_M_TRAP_DIV0 &
       + FPE_M_TRAP_INV + FPE_M_ABRUPT_UND + FPE_M_ABRUPT_DMZ
  !! Some Fortran-operatiosn apparently abuse floating point registers for fast (?) transfer
  !! of non-floating-point data. This can trigger unwanted exceptions.
  !! Lets try suppress some of that.
  integer(4) :: loose_fpe_flags = FPE_M_TRAP_DIV0 + FPE_M_ABRUPT_UND + FPE_M_ABRUPT_DMZ

contains

  ! ==================================================================================
  
  ! <SUBROUTINE NAME="atmosphere_init">
  !
  ! <OVERVIEW>
  !  public routine required for atmospheric components of coupled models
  !  Read in restart files and initialize arrays
  ! </OVERVIEW>
  !
  ! <DESCRIPTION>
  !  public routine required for atmospheric components of coupled models
  !  Read in restart files and initialize arrays
  ! </DESCRIPTION>
  !
  ! <TEMPLATE>
  !   call atmosphere_init(Time_init, Time_in, Time_step_in, Surf_diff, Grid_box)
  ! </TEMPLATE>
  !
  ! <IN NAME="Time_init" TYPE="type(time_type)" >
  !   The base (or initial) time of the experiment.
  ! </IN>
  
  ! <IN NAME="Time" TYPE="type(time_type)" >
  !   The current time.
  ! </IN>

  ! <IN NAME="Time_step" TYPE="type(time_type)" >
  !   The atmospheric model/physics time step.
  ! </IN>

  ! <INOUT NAME="Surf_diff" TYPE="type(surf_diff_type)">
  !   The surface terms for vertical diffusion that are exchanged
  !   with other component models. On input fields have not been
  !   allocated, on output all arrays are allocated.
  ! </INOUT>

  ! <INOUT NAME="Grid_box" TYPE="type(grid_box_type)">
  !   Contains information about the grid cells 
  !   In the 18dec2009 release, all other atmos models set the fields
  !   of this parameter simply to zero, or even just ignore it entirely.
  !   So do we, for now.
  ! </INOUT>

  subroutine atmosphere_init (Time_init, Time, Time_step, Surf_diff, Grid_box)

    type (time_type),     intent(in)    :: Time_init, Time, Time_step
    type(surf_diff_type), intent(inout) :: Surf_diff
    type(grid_box_type),  intent(inout) :: Grid_box

    integer :: Time_init_year, Time_init_month, Time_init_day, Time_init_hours, Time_init_minutes, Time_init_seconds
    integer :: Time_year, Time_month, Time_days, Time_hours, Time_seconds
    integer :: Time_step_days, Time_step_seconds
    integer :: unit, ierr, io
    integer :: log_unit

    integer :: pe, npes, layout(2)
    integer :: atmos_comm ! the MPI communicator for the atmosphere
    integer, allocatable, dimension(:) :: atmos_pelist ! needed to obtain the peset

    integer :: id_lonb, id_latb !, id_levelb
    integer, dimension(1:2) :: axes_2d
    logical :: used

    ! global axis info retrieved from exchange grid, for comparison with
    ! the axis info retrieved from the python side.
    real, dimension(:), allocatable :: exglon, exglat, londiff, latdiff

    ! 2-D grid cell boundary arrays, for calls to topography_mod 
    real, dimension(:,:), allocatable :: lon_boundaries, lat_boundaries
    logical :: got_topo
    real, _ALLOCATABLE :: SIGORO(:,:) ! standard deviation of topography on local domain
    real, _ALLOCATABLE :: HORO(:,:)   ! Orography height, averaged over grid cell
    real, _ALLOCATABLE :: land_frac(:,:) ! provisional land fraction mask. The real one is delivered by the coupler in every down and up step.
    real, _ALLOCATABLE :: area(:,:)   ! grid cell area (obtained from the exchange grid)

    integer id_horo, id_sigoro, id_area

    !character(len=64) :: filename

    integer(4) :: previous_fpe_flags

    if (module_is_initialized) return

    !! To enable strict floating point exceptions inside this module.
    !! Taken from the example in the Intel Fortran Compiler Reference Manual.
    !! Most probably will not work with any other compiler
    call clearstatusfpqq() ! do not inherit signalled FP exceptions from FMS
    original_fpe_flags = FOR_SET_FPE(strict_fpe_flags)

    clock_id_aeolus2 = mpp_clock_id('Aeolus2', &
         flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT )
    clockid_init = mpp_clock_id('Aeolus2_init',   &
         flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT )
    clockid_up   = mpp_clock_id('Aeolus2_update_up',   &
         flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT )
    clockid_down = mpp_clock_id('Aeolus2_update_down',   &
         flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT )
    clockid_diag = mpp_clock_id('Aeolus2_diagnostic',   &
         flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT )

    


    call mpp_clock_begin(clockid_init)
    
    !----- read namelist -----
    
    if (file_exist('input.nml')) then
       unit = open_namelist_file ( )
       ierr=1;
       do while (ierr /= 0)
          read (unit, nml=atmos_aeolus2_nml, iostat=io, end=10)
          ierr = check_nml_error (io, 'atmos_aeolus2_nml')
       enddo
10     call close_file (unit)
    endif

    !----- write version and namelist to log file -----
    call write_version_number ( version, tagname )
    if ( mpp_pe() == mpp_root_pe() ) then
       log_unit = stdlog()
       write (log_unit, nml=atmos_aeolus2_nml)
       write (*, nml=atmos_aeolus2_nml)
    endif

    ! some sanity checks
    if (0 == NoNx .or. 0 == NoNy) &
         write(*,*) 'Warning: strange grid size in atmos_aeolus2_nml:',NoNx,NoNy

    !---- compute physics/atmos time step in seconds ----
    Time_step_atmos = Time_step ! remember Time_step
    call get_time (Time_step, Time_step_seconds)
    dt_atmos = real(Time_step_seconds)

    ! initialize domain decomposition
    ! heavily inspired by atmos_spectral/tools/spec_mpp.F90

    call mpp_domains_init()
    pe = mpp_pe()
    npes = mpp_npes()

    !grid domain: atmos_fv and EBM do, by default, 1D decomposition along Y
    layout = (/1,npes/)
    !if( PRESENT(grid_layout) ) layout = grid_layout
    call mpp_define_domains( (/1,NoNx,1,NoNy/), layout, grid_domain )
    if(pe == mpp_root_pe()) call print_decomp (npes, layout, grid_domain )

    !!requirement of equal domains: can be generalized to retain mirror symmetry between N/S if unequal.
    !!the equal-domains requirement permits us to eliminate one buffer/unbuffer in the transpose_fourier routines.
    !if( mod(NoNy,layout(2)).NE.0 ) then
    !   !       call mpp_error( FATAL, 'atmos_aeolus2fms: currently requires equal grid domains on all PEs.' )
    !   write(chtmp1,'(i4)') layout(2)
    !   write(chtmp2,'(i4)') NoNy
    !   call mpp_error( FATAL, 'atmos_aeolus2fms: Requires num_lat_rows/num_pes=int;num_pes='&
    !        &//chtmp1//';num_lat_rows='//chtmp2 )
    !endif
    
    call mpp_get_compute_domain(grid_domain, is, ie, js, je)

    ! --- now we can allocate variables ---
    allocate(atmos_pelist(npes))
    call mpp_get_current_pelist( atmos_pelist, commID=atmos_comm )
    !write(*,*) 'current pelist ', atmos_pelist
    !write(*,'(a,z)') 'current atmos_comm ',atmos_comm
    deallocate(atmos_pelist)

    ! get the number of prognostic tracers
    call get_number_tracers( MODEL_ATMOS, num_prog=nr_tracers)
    q_ind = get_tracer_index( MODEL_ATMOS, 'sphum' )
    if (q_ind==NO_TRACER) &
         q_ind = get_tracer_index( MODEL_ATMOS, 'mix_rat' )
    if (q_ind==NO_TRACER) &
         call error_mesg('atmosphere_init', 'Cannot find water vapor in ATM tracer table', FATAL)

    co2_ind = get_tracer_index( MODEL_ATMOS, 'co2' )

    allocate (Surf_diff%dtmass  (is:ie, js:je) )
    allocate (Surf_diff%dflux_t (is:ie, js:je) )
    allocate (Surf_diff%dflux_tr (is:ie, js:je, nr_tracers) )
    allocate (Surf_diff%delta_t (is:ie, js:je) )
    allocate (Surf_diff%delta_tr (is:ie, js:je, nr_tracers) )
    allocate (Surf_diff%delta_u (is:ie, js:je) )
    allocate (Surf_diff%delta_v (is:ie, js:je) )
    
    Surf_diff%dtmass  = 0.0
    Surf_diff%dflux_t = 0.0
    Surf_diff%dflux_tr = 0.0
    Surf_diff%delta_t = 0.0
    Surf_diff%delta_tr = 0.0
    Surf_diff%delta_u = 0.0
    Surf_diff%delta_v = 0.0

    ! --- read restart ---
    !filename = 'INPUT/atmosphere.res.nc'
    !if (file_exist(filename)) then
    !  !call read_data(filename, 'tg', tg, grid_domain)
    !  !call read_data(filename, 'qg', qg, grid_domain)
    !else
    !  ! tg-initial atmospheric temperature default value 273 deg Kelvin
    !  ! qg-initial atmospheric specific humidity
    !  !tg = tg_init
    !  !qg = rh*sat_mixing_ratio(tg)
    !endif

    ! possibly set up some more tracers
    ! ...


    call get_date(Time_init, Time_init_year, Time_init_month, Time_init_day, &
         Time_init_hours, Time_init_minutes, Time_init_seconds)
    call get_time(Time, Time_seconds, Time_days)
    call get_time(Time_step, Time_step_seconds, Time_step_days)

    call aeolus2fms_interface_check(atmos_comm)

    call aeolus2fms_init_grid(    &
         NoNx, NoNy,       &
         !NoNz, NoNz_Strato, &
         lon_0_360,        &
         atmos_comm, pe, mpp_root_pe(), npes, &
         is, ie, js, je   &
         )

    !----- initialize atmos_axes and diagnostics

    ! If we want to use the FMS diagnostics manager infrastructure,
    ! we can set it up here. Even if we wanted to handle that from
    ! the python side only, the flux exchange module of the coupler wants to
    ! register some diagnostic fields using the atmosphere axes.
    ! Thus we really need to set those up here.
    ! The diagnostic axes require _global_ information

    allocate(lonb(1:NoNx+1), latb(1:NoNy+1), lon(1:NoNx), lat(1:NoNy))
    !allocate(level(NoNz+NoNz_Strato), levelb(NoNz+NoNz_Strato+1))
    ! Fill arrays with well-known values to ease debugging
    !do i=1,NoNx
    !   lonb(i) = i
    !   lon(i) = i
    !end do
    !lonb(NoNx+1) = NoNx+1
    !do i=1,NoNy
    !   latb(i) = i
    !   lat(i) = i
    !end do
    !latb(NoNy+1) = NoNy+1
    !do i=1,NoNz+NoNz_Strato
    !   level(i) = i
    !   levelb(i) = i
    !end do
    !levelb(NoNz+NoNz_Strato+1) = NoNz+NoNz_Strato+1

    ! receive global(!) axis labels from the python side
    call aeolus2fms_get_axes_coords(lon, lonb, lat, latb)
    write (*,*) 'lon', lon
    write (*,*) 'lonb', lonb
    write (*,*) 'lat', lat
    write (*,*) 'latb', latb
    !write (*,*) 'level', level
    !write (*,*) 'levelb', levelb

    ! ... and also retrieve global(!) axis information from the exchange grid definition,
    ! to countercheck that python and FMS agree in their grid definitions.
    allocate(exglon(1:NoNx), exglat(1:NoNy), londiff(1:NoNx), latdiff(1:NoNy))
    call get_grid_cell_centers('ATM', tile=1, glon=exglon, glat=exglat)
    londiff = exglon-lon
    if (any(abs(londiff) > 1e-13)) then
       write(*,*) 'Error: longitude axis info non-trivial mismatch between FMS and python'
       write(*,*) 'FMS', exglon
       write(*,*) 'python', lon
       write(*,*) 'FMS-python', exglon-lon
    endif
    latdiff = exglat-lat
    if (any(abs(latdiff) > 1e-13)) then
       write(*,*) 'Error: latitude axis info non-trivial mismatch between FMS and python'
       write(*,*) 'FMS', exglat
       write(*,*) 'python', lat
       write(*,*) 'FMS-python', exglat-lat
    endif
    deallocate(exglon, exglat, londiff, latdiff)

    id_lonb       = diag_axis_init('lonb', lonb, 'degrees_E', 'x', 'longitude edges', set_name=axiset, Domain2=grid_domain)
    id_latb       = diag_axis_init('latb', latb, 'degrees_N', 'y', 'latitude edges',  set_name=axiset, Domain2=grid_domain)
    atmos_axes(1) = diag_axis_init('lon', lon, 'degrees_E', 'x', 'longitude', set_name=axiset, Domain2=grid_domain, edges=id_lonb)
    atmos_axes(2) = diag_axis_init('lat', lat, 'degrees_N', 'y', 'latitude',  set_name=axiset, Domain2=grid_domain, edges=id_latb)
    !id_levelb     = diag_axis_init('levelb', levelb, 'm', 'z', 'level edges', set_name=axiset)
    !atmos_axes(3) = diag_axis_init('level', level, 'm', 'z', 'level height', set_name=axiset, edges=id_levelb)

    ! init topography
    allocate(lon_boundaries(ie-is+2, je-js+2))
    allocate(lat_boundaries(ie-is+2, je-js+2))
    write(*,*) 'allocate lon_boundaries to sizes ', ie-is+2, je-js+2
    call atmosphere_boundary(lon_boundaries, lat_boundaries, global=.false.)
    !write(*,*) 'lon_boundaries', lon_boundaries*180.0/PI ! output in degrees for better readability
    !write(*,*) 'lat_boundaries', lat_boundaries*180.0/PI

    previous_fpe_flags = FOR_SET_FPE(loose_fpe_flags)

    allocate(HORO(ie-is+1, je-js+1))
    got_topo = get_topog_mean(lon_boundaries, lat_boundaries, HORO)
    if (.not. got_topo) &
         call error_mesg ('atmosphere_init', 'could not get topography mean', FATAL)

    allocate(SIGORO(ie-is+1, je-js+1))
    got_topo = get_topog_stdev(lon_boundaries, lat_boundaries, SIGORO)
    if (.not. got_topo) &
         call error_mesg ('atmosphere_init', 'could not get topography stdev', FATAL)

    allocate(land_frac(ie-is+1, je-js+1))
    
    ! This is not needed during init
    !where (HORO > 0)
    !   land_frac = 1.0
    !elsewhere
    !   land_frac = 0.0
    !end where

    allocate(area(ie-is+1, je-js+1))
    call atmosphere_cell_area(area)

    if (THIRTY_DAY_MONTHS /= get_calendar_type() .and. NOLEAP /= get_calendar_type()) then
       call mpp_error(FATAL,'atmosphere_init: need integral number of days per year, i.e. calendar_type THIRTY_DAY_MONTHS or NOLEAP')
    endif
    DaysPerYear = days_in_year(Time_init)

    ! prepare for the use of astronomy_mod
    call astronomy_init()
    allocate(rad_lat(js:je))
    !write(*,*)'size rad_lat ',size(rad_lat),' size lat ',size(lat)
    !write(*,*) 'bounds rad_lat ', lbound(rad_lat), ubound(rad_lat), ' rad ',lbound(lat),ubound(lat)
    rad_lat(:) = lat(js:je)*PI/180.
    !write(*,*)'size rad_lat ',size(rad_lat),' size lat ',size(lat)
    !write(*,*) 'bounds rad_lat ', lbound(rad_lat), ubound(rad_lat), ' rad ',lbound(lat),ubound(lat)

    if (fix_astro) then
       Time_fix_astro = set_date(astro_date(1), astro_date(2), astro_date(3), astro_date(4), astro_date(5), astro_date(6))
       call get_time(Time_fix_astro, astro_date(6), fixdayoftheyear)
       fixdayoftheyear = MODULO(fixdayoftheyear, DaysPerYear)+1
    end if

    call aeolus2fms_init(     &
         Time_init_year, Time_init_month, Time_init_day, Time_init_seconds, &
         Time_days, Time_seconds,           &
         Time_step_days, Time_step_seconds, &
         without_topo,     &
         HORO,             &
         SIGORO,           &
         land_frac,        &
         area,             &
         update_land_frac_always, &
         nr_tracers,       &
         q_ind,            &
         co2_ind           &
         )


    ! here we could register some diagnostic fields
    axes_2d = atmos_axes(1:2)
    
    !id_astro_cosz = register_diag_field(mod_name, 'fms_astro_cosz'  , (/atmos_axes(2)/), Time, long_name='cosz is cosine of an effective zenith angle that would produce the correct daily solar flux if the sun were fixed at that position for the period of daylight; from FMS astronomy_mod', units='1') 
    !id_astro_solar = register_diag_field(mod_name, 'fms_astro_solar'  , (/atmos_axes(2)/), Time, long_name='shortwave flux factor: cosine of zenith angle * daylight fraction / (earth-sun distance squared); from FMS astronomy_mod', units='1') 


    !liq_wat = get_tracer_index(MODEL_ATMOS, 'liq_wat')
    !ice_wat = get_tracer_index(MODEL_ATMOS, 'ice_wat')

    id_horo = register_static_field(mod_name, "horo", axes_2d, "Orographic Height from FMS", "m")
    id_sigoro = register_static_field(mod_name, "sigoro", axes_2d, "Standard Deviation of Orographic Height from FMS", "m")
    if (id_horo > 0) used = send_data(id_horo, HORO, Time)
    if (id_sigoro > 0) used = send_data(id_sigoro, SIGORO, Time)

    id_area = register_static_field(mod_name, "area", axes_2d, "Cell Area", "m**2")
    if (id_area > 0) used = send_data(id_area, area, Time)

    !deallocate(lonb, latb, lon, lat)
    !deallocate(level, levelb)
    deallocate(HORO)
    deallocate(SIGORO)
    deallocate(land_frac)
    deallocate(area)

    id_flux_lw_sf_up = register_diag_field(mod_name, 'flux_lwr_sf_up', axes_2d, Time, long_name='lwr upward flux calculated in atmos wrapper', units='W/m2')

    call diag_manager_end(Time) ! to get static fields onto disk before something crashes

    !call register_Aeolus2_diag_vars(Time)

    atmos_stock_water_id = register_diag_field(mod_name, 'waterstock', Time, long_name = 'Total water content of atmosphere', units='')
    atmos_stock_heat_id = register_diag_field(mod_name, 'heatstock', Time, long_name = 'Total heat content of atmosphere', units='J')

    !call send_Aeolus2_static_diag_vars(Time)

    call mpp_clock_end(clockid_init)
    
    module_is_initialized = .true.

    !! To enable strict floating point exceptions inside this module.
    !! Taken from the example in the Intel Fortran Compiler Reference Manual.
    !! Most probably will not work with any other compiler
    original_fpe_flags = FOR_SET_FPE(original_fpe_flags)

    write (*,*) 'atmosphere_init finished'
  end subroutine atmosphere_init
  ! </SUBROUTINE>


  ! ==================================================================================
  ! <SUBROUTINE NAME="atmosphere_down">
  ! <OVERVIEW>
  !  public routine required for atmospheric components of coupled models
  ! </OVERVIEW>

  ! <DESCRIPTION>
  !   atmosphere_down
  !  This routine calls the dynamical core and the
  !  "downward pass" of the atmospheric physics.
  !  It should only be called once per time step and before
  !  calling atmosphere_up.
  ! </DESCRIPTION>

  ! <TEMPLATE>
  !     call atmosphere_down(Time, frac_land, t_surf, albedo, rough_mom, &
  !                          u_star, b_star, q_star, dtau_du, dtau_dv, tau_x, tau_y,  &
  !                          gust, coszen, coszen_radwt, net_surf_sw_down, surf_lw_down, &
  !                          Surf_diff)
  ! </TEMPLATE>

  ! <IN NAME="Time" TYPE="type(time_type)">
  !  time at the current time level (tau)
  ! </IN>

  ! <IN NAME="frac_land" TYPE="real, dimension(is:ie,js:je)">
  !  fraction (0. to 1.) of underlying surface which covered by land
  ! </IN>

  ! <IN NAME="t_surf" TYPE="real, dimension(is:ie,js:je)">
  !  surface (skin) temperature (in deg k)
  ! </IN>

  ! <IN NAME="albedo" TYPE="real, dimension(is:ie,js:je)">
  !  surface albedo
  ! </IN>

  ! <IN NAME="rough_mom" TYPE="real, dimension(is:ie,js:je)">
  !  momentum roughness length (units=m)
  ! </IN>

  ! <IN NAME="u_star" TYPE=" real, dimension(is:ie,js:je)">
  !  friction velocity
  ! </IN>

  ! <IN NAME="b_star" TYPE="real, dimension(is:ie,js:je) ">
  !  buoyancy scale
  ! </IN>

  ! <IN NAME="q_star" TYPE="real, dimension(is:ie,js:je) ">
  !  moisture scale
  ! </IN>

  ! <IN NAME="dtau_du" TYPE=" real, dimension(is:ie,js:je)">
  !  derivative of zonal wind stress w.r.t. the lowest level wind speed
  ! </IN>

  ! <IN NAME="dtau_dv" TYPE=" real, dimension(is:ie,js:je)">
  !  derivative of meridional wind stress w.r.t. the lowest level wind speed
  ! </IN>

  ! <IN NAME="frac_open_sea" TYPE=" real, dimension(is:ie,js:je)">
  !  fraction of oopen sea without ice
  ! </IN>

  ! <OUT NAME="gust" TYPE="real, dimension(is:ie,js:je)" >
  !  wind gustiness
  ! </OUT>

  ! <OUT NAME="coszen" TYPE="real, dimension(is:ie,js:je)" >
  !  cosine of the zenith angle
  ! </OUT>

  ! <OUT NAME="coszen_radwt" TYPE="real, dimension(is:ie,js:je)" >
  !  cosine of the zenith angle, radiation-weighted, Climber-2 style.
  !  Ignored with Aeolus2.
  ! </OUT>

  ! <OUT NAME="flux_sw" TYPE="real, dimension(is:ie,js:je)" >
  !  net shortwave surface flux (down minus up) (in watts/m**2)
  ! </OUT>

  ! <OUT NAME="flux_lw" TYPE="real, dimension(is:ie,js:je)" >
  !  downward longwave surface flux (in watts/m**2)
  ! </OUT>

  ! <INOUT NAME="Surf_diff" TYPE="type(surf_diff_type)" >
  !  Surface diffusion terms computed by the vertical diffusion scheme
  !  In the GFDL-provided example atmosphere implementations, all of the
  !  fields of this variable are asigned new values here.
  ! </INOUT>

  subroutine atmosphere_down (           &
       ! input parameters                &
       Time,                             &
       frac_land, t_surf,  albedo,       &
       albedo_vis_dir, albedo_nir_dir,   &
       albedo_vis_dif, albedo_nir_dif,   &
       rough_mom,                        &
       u_star,  b_star, q_star,          &
       dtau_du, dtau_dv,                 & ! derivatives of wind stress w.r.t. the lowest level wind speed
       ! inout parameters                &
       u_flux, v_flux,                   & ! wind stress
       frac_open_sea,                    & ! non-seaice fraction (%)
       ! output parameters               &
       gust, coszen, coszen_radwt,       &
       flux_sw,                          &
       flux_sw_dir, flux_sw_dif,         &
       flux_sw_down_vis_dir,             &
       flux_sw_down_vis_dif,             &
       flux_sw_down_total_dir,           &
       flux_sw_down_total_dif,           &
       flux_sw_vis,                      &
       flux_sw_vis_dir,                  &
       flux_sw_vis_dif,                  &
       flux_lw,                          &
       Surf_diff                         )

    type(time_type),  intent(in)                 :: Time
    real, intent(in),     dimension(is:ie,js:je) :: frac_land, t_surf, albedo
    real, intent(in),     dimension(is:ie,js:je) :: albedo_vis_dir, albedo_nir_dir
    real, intent(in),     dimension(is:ie,js:je) :: albedo_vis_dif, albedo_nir_dif
    real, intent(in),     dimension(is:ie,js:je) :: rough_mom
    real, intent(in),     dimension(is:ie,js:je) :: u_star, b_star, q_star, dtau_du, dtau_dv
    real, intent(inout),  dimension(is:ie,js:je) :: u_flux, v_flux
    real, intent(in),     dimension(is:ie,js:je) :: frac_open_sea
    real, intent(out),    dimension(is:ie,js:je) :: gust, coszen
    real, intent(inout),  dimension(is:ie,js:je) :: coszen_radwt ! declare as inout to suppress compiler warning
    real, intent(out),    dimension(is:ie,js:je) :: flux_sw
    real, intent(out),    dimension(is:ie,js:je) :: flux_sw_dir, flux_sw_dif, flux_sw_down_vis_dir
    real, intent(out),    dimension(is:ie,js:je) :: flux_sw_down_vis_dif, flux_sw_down_total_dir
    real, intent(out),    dimension(is:ie,js:je) :: flux_sw_down_total_dif, flux_sw_vis
    real, intent(out),    dimension(is:ie,js:je) :: flux_sw_vis_dir, flux_sw_vis_dif
    real, intent(out),    dimension(is:ie,js:je) :: flux_lw
    type(surf_diff_type),   intent(inout)        :: Surf_diff

    integer :: Time_seconds, Time_days
    integer :: dayoftheyear
    real, dimension(js:je) :: cosz, solar, cosz_radwt
    real, dimension(is:ie,js:je) :: flux_lw_sf_up


    integer :: j, used

    call mpp_clock_begin (clockid_down)

    !! To enable strict floating point exceptions inside this module.
    !! Taken from the example in the Intel Fortran Compiler Reference Manual.
    !! Most probably will not work with any other compiler
    call clearstatusfpqq() ! do not inherit signalled FP exceptions from FMS
    original_fpe_flags = FOR_SET_FPE(strict_fpe_flags)

    call get_time(Time, Time_seconds, Time_days)
    write (*,*) 'atmosphere_down starting at', Time_days, Time_seconds

    !type(time_type) :: Time_prev, Time_next
    !
    !Time_prev = Time                       ! two time-level scheme
    !Time_next = Time + Time_step_atmos
   
    !---- dynamics -----

    !call timing_on('fv_dynamics')
    !call mpp_clock_begin (id_dynam)
    !call mpp_clock_end (id_dynam)
    !call timing_off('fv_dynamics')

    !---- call physics -----

    !call timing_on('fv_physics_down')
    !call mpp_clock_begin (id_phys_down)
    !call mpp_clock_end (id_phys_down)
    !call timing_off('fv_physics_down')

    !call fv_print_chksums( 'Exiting  atmosphere_down' )

    ! TODO: possible optimization: only re-calculate solar data once per day
    if (fix_astro) then
       dayoftheyear = fixdayoftheyear
    else
       dayoftheyear = MODULO(Time_days, DaysPerYear)+1
    end if
    !write(*,*)'size rad_lat ',size(rad_lat),' size cosz ',size(cosz),' size solar ', size(solar)
    if (id_astro_cosz > 0 .or. id_astro_solar > 0) then
       if (fix_astro) then
          call daily_mean_solar  (rad_lat, Time_fix_astro, cosz, solar)
       else
          !write(*,*)'size rad_lat ',size(rad_lat)
          !write(*,*) 'bounds rad_lat ', lbound(rad_lat), ubound(rad_lat)
          !write(*,*)'size cosz ',size(cosz)
          !write(*,*) 'bounds cosz ', lbound(cosz), ubound(cosz)
          call daily_mean_solar  (rad_lat, Time,           cosz, solar)
       end if
       if (id_astro_cosz > 0) used = send_data(id_astro_cosz, cosz, Time)
       if (id_astro_solar > 0) used = send_data(id_astro_solar, solar, Time)
    end if
    ! the return value solar from FMS-astronomy must be multiplied by solar constant to
    ! obtain the same meaning as solarm from sinsol
    solar = solar * solar_constant
    do j = js,je
       coszen(:,j) = cosz(j) ! pass back to coupler the non-radiation-weighted coszen
    enddo
    !cosz_radwt = 0.0 was already set during initialisation, and is not used with Aeolus2

    !write(*,*) 'aeolus2fms_atmosphere_down()', 'delta_t ', size(Surf_diff%delta_t), &
    !     lbound(Surf_diff%delta_t), ubound(Surf_diff%delta_t),                    &
    !                                         'delta_tr', size(Surf_diff%delta_tr), &
    !     lbound(Surf_diff%delta_tr), ubound(Surf_diff%delta_tr),                   &
    !                                         'dflux_tr', size(Surf_diff%dflux_tr), &
    !     lbound(Surf_diff%dflux_tr), ubound(Surf_diff%dflux_tr)


    ! same as in surface_flux(), before interpolation to reduced grid
    flux_lw_sf_up = stefan*t_surf**4
    if (id_flux_lw_sf_up > 0) used = send_data(id_flux_lw_sf_up, flux_lw_sf_up, Time)

    call aeolus2fms_atmosphere_down (        &
         ! input parameters                &
         Time_seconds, Time_days,          &
         dayoftheyear,                     &
         frac_land, t_surf,  albedo,       &
         albedo_vis_dir, albedo_nir_dir,   &
         albedo_vis_dif, albedo_nir_dif,   &
         rough_mom,                        &
         u_star,  b_star, q_star,          &
         !frac_open_sea,                   & ! non-seaice fraction (%)
         dtau_du, dtau_dv,                 & ! derivatives of wind stress w.r.t. the lowest level wind speed
         ! inout parameters                &
         u_flux, v_flux,                   & ! wind stress
         ! output parameters               &
         gust,                             &
         ! coszen, cosz, solar were already calculated above, and thus are now input-parameters for Aeolus
         coszen,                           &
         !cosz_radwt,                      &
         solar,                            & ! here pass the 1-D variables.
         flux_lw_sf_up,                    & 
         ! output parameters
         flux_sw,                          &
         flux_sw_dir, flux_sw_dif,         &
         flux_sw_down_vis_dir,             &
         flux_sw_down_vis_dif,             &
         flux_sw_down_total_dir,           &
         flux_sw_down_total_dif,           &
         flux_sw_vis,                      &
         flux_sw_vis_dir,                  &
         flux_sw_vis_dif,                  &
         flux_lw,                          &
         Surf_diff%delta_tr,               &
         Surf_diff%dflux_tr,               &
         Surf_diff%delta_t,                &
         Surf_diff%dflux_t,                &
         Surf_diff%dtmass,                 &
         Surf_diff%delta_u,                &
         Surf_diff%delta_v                 &
         )

    !! To enable strict floating point exceptions inside this module.
    !! Taken from the example in the Intel Fortran Compiler Reference Manual.
    !! Most probably will not work with any other compiler
    original_fpe_flags = FOR_SET_FPE(original_fpe_flags)

    call mpp_clock_end(clockid_down)
  end subroutine atmosphere_down

  
  ! ==================================================================================
  ! <SUBROUTINE NAME="atmosphere_up">
  ! <OVERVIEW>
  !  public routine required for atmospheric components of coupled models
  !</OVERVIEW>
  !
  !<DESCRIPTION>
  !   atmosphere_up
  !  This routine calls the "upward pass" of the atmospheric physics,
  !  spherical diagnostics, and time differencing.  The prognostic
  !  variables are advanced to the next time step.  It should only be
  !  called once per time step and after calling atmosphere_down.
  !</DESCRIPTION>

  !   <TEMPLATE>
  !     call  atmosphere_up(Time, frac_land, Surf_diff, lprec, fprec, gust)
  !   </TEMPLATE>

  ! <IN NAME="Time" TYPE="type(time_type)">
  !  time at the current time level (tau)
  ! </IN>

  ! <IN NAME="frac_land" TYPE="real, dimension(is:ie,js:je)">
  !  fraction (0. to 1.) of underlying surface which covered by land
  ! </IN>

  ! <OUT NAME="lprec" TYPE="real, dimension(is:ie,js:je)" >
  !  liquid precipitation rate (rain) in kg/m2/s
  ! </OUT>

  ! <OUT NAME="fprec" TYPE="real, dimension(is:ie,js:je)" >
  !  frozen precipitation rate (snow) in kg/m2/s
  ! </OUT>

  ! <OUT NAME="gust" TYPE="real, dimension(is:ie,js:je)" >
  !  wind gustiness
  ! </OUT>

  ! <INOUT NAME=" Surf_diff" TYPE="type(surf_diff_type)" >
  !  surface diffusion terms computed by the vertical diffusion scheme
  ! </INOUT>

  ! <IN NAME="u_star" TYPE=" real, dimension(is:ie,js:je)">
  !  friction velocity
  ! </IN>

  ! <IN NAME="b_star" TYPE="real, dimension(is:ie,js:je) ">
  !  buoyancy scale
  ! </IN>

  ! <IN NAME="q_star" TYPE="real, dimension(is:ie,js:je) ">
  !  moisture scale
  ! </IN>

  subroutine atmosphere_up (Time,  frac_land, Surf_diff, lprec, fprec, gust, u_star, b_star, q_star )

    type(time_type),      intent(in)                          :: Time
    real,                 intent(in),  dimension(is:ie,js:je) :: frac_land
    type(surf_diff_type), intent(in)                          :: Surf_diff
    real,                 intent(out), dimension(is:ie,js:je) :: lprec, fprec, gust
    ! u_star, b_star, q_star are not used in EBM atmosphere_up
    real,                 intent(in),  dimension(is:ie,js:je) :: u_star, b_star, q_star

    ! Surf_diff%delta_t and Surf_diff%delta_tr have just been updated
    ! in the caller update_atmos_model_up(), and can be used now.
    ! None of the GFDL-provided atmosphere example models
    ! uses any of the other Surf_diff fields in the upward routines.

    integer :: Time_year, Time_month   ! 
    integer :: Time_seconds, Time_days ! seconds-of-day and total days since start of time axis, not day-of-year!
    integer :: used
    type(time_type) :: Time_next ! to get time stamp of diag output right
    real    :: waterstock, heatstock

    call mpp_clock_begin(clockid_up)

    !! To enable strict floating point exceptions inside this module.
    !! Taken from the example in the Intel Fortran Compiler Reference Manual.
    !! Most probably will not work with any other compiler
    call clearstatusfpqq() ! do not inherit signalled FP exceptions from FMS
    original_fpe_flags = FOR_SET_FPE(strict_fpe_flags)

    ! here we need only year and month
    call get_date(Time, Time_year, Time_month, used, used, used, used)
    call get_time(Time, Time_seconds, Time_days)

    !type(time_type) :: !Time_prev, Time_next
    !
    !Time_prev = Time                       ! two time-level scheme
    !Time_next = Time + Time_step_atmos

    !call timing_on('fv_physics_up')
    !call mpp_clock_begin (id_phys_up)
    !call mpp_clock_end (id_phys_up)
    !call timing_off('fv_physics_up')

    !---- diagnostics for FV dynamics -----
    !call timing_on('FV_DIAG')
    !call mpp_clock_begin(id_fv_diag)
    !fv_time = Time + Time_step_atmos

    !call fv_diag(fv_time, nlon, mlat, nlev, beglat, endlat, &
    !     ncnst, zvir, dt_atmos, .false.)
    !call timing_off('FV_DIAG')
    !call mpp_clock_end(id_fv_diag)

    call aeolus2fms_atmosphere_up ( &
         ! input parameters
         Time_year, Time_month,   &
         Time_seconds, Time_days, & ! seconds-of-day and days-since-time-axis-start
         frac_land,               &
         Surf_diff%delta_t,       &
         Surf_diff%delta_tr,      &
         ! output parameters
         lprec, fprec, gust,      &
         ! input parameters again
         u_star, b_star, q_star   &
         )

    Time_next = Time + Time_step_atmos
    !call send_Aeolus2fms_diag_vars(Time_next)

    if (atmos_stock_water_id > 0) then
       call get_stock_atm_global(ISTOCK_WATER, waterstock)
       used = send_data(atmos_stock_water_id, waterstock, Time)
    end if
    if (atmos_stock_heat_id > 0) then
       call get_stock_atm_global(ISTOCK_HEAT, heatstock)
       used = send_data(atmos_stock_heat_id, heatstock, Time)
    end if
    
    !! To enable strict floating point exceptions inside this module.
    !! Taken from the example in the Intel Fortran Compiler Reference Manual.
    !! Most probably will not work with any other compiler
    original_fpe_flags = FOR_SET_FPE(original_fpe_flags)

    call mpp_clock_end(clockid_up)
  end subroutine atmosphere_up
  ! </SUBROUTINE>


  ! =====================================================================
  ! <SUBROUTINE NAME="atmosphere_end">
  !
  ! <OVERVIEW>
  !  write out restart file
  !  Termination routine for atmosphere_mod.
  ! </OVERVIEW>
  ! <DESCRIPTION>
  !  public routine required for atmospheric components of coupled models
  !  write out restart file
  !  Termination routine for atmosphere_mod.
  ! </DESCRIPTION>
  !
  ! <TEMPLATE>
  !     call atmosphere_end(Time, Grid_box)
  ! </TEMPLATE>
  !
  ! <IN NAME="Time" TYPE="type(time_type)">
  !  time at the current time level (tau)
  ! </IN>
  !
  ! <INOUT NAME="Grid_box" TYPE="type(grid_box_type)">
  !   Contains information about the grid cells 
  ! </INOUT>
  subroutine atmosphere_end (Time, Grid_box)
    type (time_type), intent(in) :: Time
    type(grid_box_type),  intent(inout) :: Grid_box ! unused

    integer :: Time_days, Time_seconds

    if (.NOT. module_is_initialized) return

    !----- initialize domains for writing global physics data -----

    !call set_domain ( fv_domain )

    ! Note that in atmos_coupled/atmos_model.F90:atmos_model_local_restart()
    ! already the fields lprec,fprec,gust,t_bot,q_bot are saved, and restored in
    ! atmos_model_init()
    call get_time(Time, Time_seconds, Time_days)
    call aeolus2fms_finish(Time_seconds, Time_days)

    ! de-allocate memory
    
    ! sigh. the surf_diff_type cannot de-allocated in any atmosphere implementation,
    ! because it is never passed to any end-subroutine
    !deallocate (Surf_diff%dtmass)
    !deallocate (Surf_diff%dflux_t)
    !deallocate (Surf_diff%dflux_tr)
    !deallocate (Surf_diff%delta_t)
    !deallocate (Surf_diff%delta_tr)
    !deallocate (Surf_diff%delta_u)
    !deallocate (Surf_diff%delta_v)

    deallocate(lonb, latb, lon, lat)
    deallocate(rad_lat)

    !deallocate(Aeolus2_diag_ids)
    !deallocate(Aeolus2_diag_types)

    module_is_initialized = .false.
  end subroutine atmosphere_end
  ! </SUBROUTINE>


  !#######################################################################
  ! <SUBROUTINE NAME="atmosphere_restart">
  ! <DESCRIPTION>
  !  Write out restart data
  ! </DESCRIPTION>
  subroutine atmosphere_restart(timestamp)
    ! format yyyymmdd.hhmmss, from date_to_string()
    character(len=*),  intent(in) :: timestamp 

    ! sigh deeply. For the time stamp in the aeolus restart file, we want
    ! days since 0000-01-01. But all we get is this timestamp string,
    ! which even cannot be converted back into a time_type, because
    ! the FMS-provided set_date() expects a different format than
    ! what comes out of date_to_string().
    ! On the other hand, set_date() already has all the calendar-handling built in,
    ! so we really want to use that and not re-invent the wheel.
    ! First apporach: lets try some string manipulation in Fortran...
    ! Convert timestamp to the format 1980-01-01 00:00:00

    type(time_type) :: Time
    integer :: Time_days, Time_seconds
    character(len=50) :: timebuf
    timebuf = timestamp(1:4)//'-'//timestamp(5:6)//'-'//timestamp(7:8)&
         //' '//timestamp(10:11)//':'//timestamp(12:13)//':'//timestamp(14:15)
    write(*,*) 'atmosphere_restart: timebuf>', timebuf, '<'


    ! TODO: if necessary, save things into atmosphere.res.nc

    Time = set_date(timebuf)
    call get_time(Time, Time_seconds, Time_days)

    ! Note that in atmos_coupled/atmos_model.F90:atmos_model_local_restart()
    ! already the fields lprec,fprec,gust,t_bot,q_bot are saved, and restored in
    ! atmos_model_init()
    call aeolus2fms_restart(Time_seconds, Time_days)

  end subroutine atmosphere_restart
  ! </SUBROUTINE>


  ! ==================================================================================
  ! <SUBROUTINE NAME="atmosphere_resolution">
  ! <OVERVIEW>
  !    returns the number of longitude and latitude grid points
  !    for either the local PEs grid (default) or the global grid
  ! </OVERVIEW>

  ! <DESCRIPTION>
  !  public routine required for atmospheric components of coupled models
  !    returns the number of longitude and latitude grid points
  !    for either the local PEs grid (default) or the global grid
  ! </DESCRIPTION>

  ! <TEMPLATE>
  !   call  atmosphere_resolution(num_lon, num_lat, global)
  ! </TEMPLATE>
  
  ! <OUT NAME="num_lon" TYPE="integer">
  !  The number of longitude points in the compute domain.
  ! </OUT>

  ! <OUT NAME="num_lat" TYPE="integer">
  !  The number of latitude points in the compute domain.
  ! </OUT>

  ! <IN NAME="global" TYPE="logical,optional">
  !  Flag that specifies whether the returned compute domain size is
  !  for the global grid (TRUE) or for the current processor (FALSE).
  ! </IN>

  subroutine atmosphere_resolution (num_lon, num_lat, global)
    integer, intent(out)          :: num_lon, num_lat
    logical, intent(in), optional :: global
    logical :: local

    local = .TRUE.
    if (PRESENT(global)) local = .NOT.global

    if (local) then
       num_lon = ie - is + 1
       num_lat = je - js + 1
    else
       num_lon = NoNx
       num_lat = NoNy
    end if
  end subroutine atmosphere_resolution
  ! </SUBROUTINE>


  ! ==================================================================================
  ! <SUBROUTINE NAME="atmosphere_cell_area">
  !
  ! <OVERVIEW>
  !   returns grid cell areas for the domain
  ! </OVERVIEW>
  !
  ! <DESCRIPTION>
  !   returns grid cell areas for the domain
  ! </DESCRIPTION>

  subroutine atmosphere_cell_area  (area_out)
    real, dimension(:,:),  intent(out) :: area_out

    integer :: xsize, ysize
    xsize = ie-is+1
    ysize = je-js+1
    if(any(shape(area_out) /= (/xsize,ysize/))) then
       call mpp_error(FATAL,'atmosphere_cell_area: argument has wrong shape. Its shape is ', &
            shape(area_out),'  It should be ',(/xsize,ysize/))
    endif

    call get_grid_cell_area('ATM', tile=1, cellarea=area_out, domain=grid_domain)

    !do xsize=is,ie
    !   do ysize=js,je
    !      write(*,*) 'i ', xsize, ' j ', ysize, ' area ', area_out(xsize, ysize)
    !   end do
    !end do

  end subroutine atmosphere_cell_area
  ! </SUBROUTINE>


  ! ==================================================================================
  ! <SUBROUTINE NAME="atmosphere_boundary">
  ! <OVERVIEW>
  !    returns the longitude and latitude grid box edges
  !    for either the local PEs grid (default) or the global grid
  ! </OVERVIEW>
  ! <DESCRIPTION>
  !  public routine required for atmospheric components of coupled models
  !    returns the longitude and latitude grid box edges
  !    for either the local PEs grid (default) or the global grid
  ! </DESCRIPTION>

  ! <TEMPLATE>
  !   call atmosphere_boundary(lon_boundaries, lat_boundaries, global)
  ! </TEMPLATE>

  ! <OUT NAME = "lon_boundaries" TYPE = "real,dimension(:,:)"> 
  !  The west-to-east longitude edges of grid boxes (in radians).
  ! </OUT>

  ! <OUT NAME = "lat_boundaries" TYPE = "real,dimension(:,:)"> 
  !  The south-to-north latitude edges of grid boxes (in radians).
  ! </OUT>

  ! <IN NAME="global" TYPE="logical, optional">
  !  Flag that specifies whether the returned grid box edges are
  !  for the global grid (TRUE) or for the current processor (FALSE).
  ! </IN>

  subroutine atmosphere_boundary (lon_boundaries, lat_boundaries, global)
    real,    intent(out)          :: lon_boundaries(:,:), lat_boundaries(:,:)   ! radians !!
    logical, intent(in), optional :: global

    logical :: local
    integer :: i, i1, i2, j1, j2

    i1=is; i2=ie+1; j1=js; j2=je+1
    local = .TRUE.
    if( PRESENT(global) )local = .NOT.global
    if (local) then
       write(*,*) 'atmosphere_boundary local'
    else
       write(*,*) 'atmosphere_boundary global'
       i1=1; i2=size(lonb); j1=1; j2=size(latb)
    end if

    write(*,*) 'shape of lon_boundaries', shape(lon_boundaries)
    write(*,*) 'shape of lonb', shape(lonb)
    call flush(6)
    do i = 1, size(lon_boundaries,2)
       lon_boundaries(:,i) = lonb(i1:i2)
    end do

    do i = 1, size(lon_boundaries,1)
       lat_boundaries(i,:) = latb(j1:j2)
    end do

    lon_boundaries = lon_boundaries/180.0*PI
    lat_boundaries = lat_boundaries/180.0*PI
  end subroutine atmosphere_boundary
  ! </SUBROUTINE>



  ! ==================================================================================
  ! <SUBROUTINE NAME="atmosphere_domain">
  !
  ! <OVERVIEW>
  !    returns the domain2D variable associated with the coupling grid
  ! </OVERVIEW>
  !
  ! <DESCRIPTION>
  !  public routine required for atmospheric components of coupled models
  !    returns the domain2D variable associated with the coupling grid
  !  note: coupling is done using the mass/temperature grid with no halos
  !
  !OUTPUT
  !  Domain   The domain2D variable describing the grid used for coupling.
  !           For the B-grid, this corresponds to the temperature grid
  !           without halos.
  ! </DESCRIPTION>
 
  ! <TEMPLATE>
  !     call atmosphere_domain(domain)
  ! </TEMPLATE>
 
  ! <OUT NAME="domain" TYPE="type (domain2D)">
  !  The domain2D variable describing the grid used for coupling.
  !  For the B-grid, this corresponds to the temperature grid
  !  without halos.
  ! </OUT>
 
  subroutine atmosphere_domain (domain)
    type(domain2D), intent(out) :: domain
    domain = grid_domain
  end subroutine atmosphere_domain
  !</SUBROUTINE>

  ! ==================================================================================
  ! <SUBROUTINE NAME="get_atmosphere_axes">
  !
  ! <OVERVIEW>
  !    returns the axis indices associated with the coupling grid
  ! </OVERVIEW>
 
  ! <DESCRIPTION>
  !  public routine required for atmospheric components of coupled models
  !    returns the axis indices associated with the coupling grid
  ! </DESCRIPTION>
  
  ! <TEMPLATE>
  !     call  get_atmosphere_axes(axes)
  ! </TEMPLATE>

  ! <OUT NAME="axes" TYPE="integer,dimension(:)" >
  !  The axis identifiers for the atmospheric grids.
  !  The size of axes must be least 1 but not greater than 4.
  !  The axes are returned in the order (/ x, y, p_full, p_half /)
  !</OUT>
  subroutine get_atmosphere_axes ( axes )
    integer, intent(out), dimension(:) :: axes

    !----- returns the axis indices for the atmospheric (mass) grid -----

    if ( size(axes(:)) < 1 .or. size(axes(:)) > 4 ) call error_mesg (    &
         'get_atmosphere_axes in atmosphere_mod', &
         'size of argument is incorrect', FATAL   )
    
    axes (1:size(axes(:))) = atmos_axes (1:size(axes(:)))
    
  end subroutine get_atmosphere_axes
  !</SUBROUTINE>


  ! ==================================================================================
  ! <SUBROUTINE NAME="get_bottom_mass">
  !
  ! <OVERVIEW>
  ! returns temp, sphum, pres, height at the lowest model level
  !         and surface pressure
  ! </OVERVIEW>
  
  ! <DESCRIPTION>
  !  public routine required for atmospheric components of coupled models
  !  returns temp, sphum, pres, height at the lowest model level
  !         and surface pressure
  ! </DESCRIPTION>

  !   <TEMPLATE>
  !     call get_bottom_mass (t_bot, tr_bot, p_bot, z_bot, p_surf, slp)
  !   </TEMPLATE>

  !   <OUT NAME="t_bot" TYPE="real">
  !                near surface temperature in degrees Kelvin
  !   </OUT>

  !   <OUT NAME="tr_bot" TYPE="real">
  !                array of tracers, 1st tracer is
  !                near surface mixing  ratio
  !   </OUT>

  !   <OUT NAME="p_bot" TYPE="real">
  !                 pressure at which atmos near surface values are assumed to be defined
  !   </OUT>

  !   <OUT NAME="z_bot" TYPE="real">
  !                 height at which atmos near surface values are assumed to be defined
  !   </OUT>

  !   <OUT NAME="p_surf" TYPE="real">
  !                 surface pressure
  !   </OUT>

  !   <OUT NAME="slp" TYPE="real">
  !                 sea level pressure
  !   </OUT>

  subroutine get_bottom_mass (t_bot, tr_bot, p_bot, z_bot, p_surf, slp)
    real, intent(out), dimension(is:ie,js:je) :: t_bot, p_bot, z_bot, p_surf, slp
    real, intent(out), dimension(:,:,:) :: tr_bot

    !! To enable strict floating point exceptions inside this module.
    !! Taken from the example in the Intel Fortran Compiler Reference Manual.
    !! Most probably will not work with any other compiler
    call clearstatusfpqq() ! do not inherit signalled FP exceptions from FMS
    original_fpe_flags = FOR_SET_FPE(strict_fpe_flags)

    ! for now, just set some dummy values which prevent the coupler from
    ! division by zero
    !t_bot = 274.0
    !tr_bot = 
    !p_bot = 1
    !z_bot = 1
    !p_surf = 1
    !slp = 1

    call aeolus2fms_get_bottom_mass(t_bot, tr_bot, p_bot, z_bot, p_surf, slp)

    !! To enable strict floating point exceptions inside this module.
    !! Taken from the example in the Intel Fortran Compiler Reference Manual.
    !! Most probably will not work with any other compiler
    original_fpe_flags = FOR_SET_FPE(original_fpe_flags)
    
  end subroutine get_bottom_mass



  ! ==================================================================================
  ! <SUBROUTINE NAME="get_bottom_wind">
  !
  ! <OVERVIEW>
  ! returns u and v on the mass grid at the lowest model level
  ! </OVERVIEW>
  
  ! <DESCRIPTION>
  !  public routine required for atmospheric components of coupled models
  !  returns u and v on the mass grid at the lowest model level
  ! </DESCRIPTION>

  !   <TEMPLATE>
  !     call get_bottom_wind (u_bot_out, v_bot_out)
  !   </TEMPLATE>

  !   <OUT NAME="u_bot"  TYPE="real" >
  !    near surface zonal      wind
  !   </OUT>

  !   <OUT NAME="v_bot"  TYPE="real" >
  !    near surface meridional wind
  !   </OUT>

  subroutine get_bottom_wind (u_bot, v_bot)
    real, intent(out), dimension(is:ie,js:je) :: u_bot, v_bot
    call aeolus2fms_get_bottom_wind(u_bot, v_bot)
  end subroutine get_bottom_wind

  ! =====================================================================
  ! <SUBROUTINE NAME="get_stock_pe">
  !
  ! <OVERVIEW>
  !  Return for this pe the total content of a kind of stock (water, salt, heat) in the atmosphere.
  ! </OVERVIEW>
  ! <DESCRIPTION>
  !  Return for this pe the total content of a kind of stock (water, salt, heat) in the atmosphere.
  ! </DESCRIPTION>

  ! <TEMPLATE>
  !     call get_stock_pe(index, value)
  ! </TEMPLATE>

  ! <IN NAME="idx" TYPE="integer">
  !   index to indicate wich kind of stock to get
  ! </IN>
  ! <OUT NAME="val" TYPE="real">
  !   value of requested stock
  ! </OUT>
  subroutine get_stock_pe(idx, val)
    integer, intent(in) :: idx
    real, intent(out)   :: val

    call aeolus2fms_get_stock_pe(idx, val)
  end subroutine get_stock_pe
  !</SUBROUTINE>

  ! =====================================================================
  ! <SUBROUTINE NAME="get_stock_atm_global">
  !
  ! <OVERVIEW>
  !  Return global total content of a kind of stock (water, salt, heat) in the atmosphere.
  ! </OVERVIEW>
  ! <DESCRIPTION>
  !  Return global total content of a kind of stock (water, salt, heat) in the atmosphere.
  ! </DESCRIPTION>

  ! <TEMPLATE>
  !     call get_stock_atm_global(index, value)
  ! </TEMPLATE>

  ! <IN NAME="idx" TYPE="integer">
  !   index to indicate wich kind of stock to get
  ! </IN>
  ! <OUT NAME="val" TYPE="real">
  !   value of requested stock
  ! </OUT>
  subroutine get_stock_atm_global(idx, val)
    integer, intent(in) :: idx
    real, intent(out)   :: val

    call aeolus2fms_get_stock_pe(idx, val)
    call mpp_sum(val)
  end subroutine get_stock_atm_global
  !</SUBROUTINE>


  ! This private subroutine is simply copied from atmos_spectral/tools/spec_mpp.F90
  ! because it is not accesible in there from here
  subroutine print_decomp (npes, layout, Domain)
    integer, intent(in) :: npes, layout(2)
    type(domain2D), intent(in) :: Domain
    integer, dimension(0:npes-1) :: xsize, ysize
    integer :: i, j, xlist(layout(1)), ylist(layout(2))
    type (domain1D) :: Xdom, Ydom
    
    call mpp_get_compute_domains   ( Domain, xsize=xsize, ysize=ysize )
    call mpp_get_domain_components ( Domain, Xdom, Ydom )
    call mpp_get_pelist ( Xdom, xlist )
    call mpp_get_pelist ( Ydom, ylist )

    write (*,100)
    write (*,110) (xsize(xlist(i)),i=1,layout(1))
    write (*,120) (ysize(ylist(j)),j=1,layout(2))
    
100 format ('ATMOS MODEL DOMAIN DECOMPOSITION')
110 format ('  X-AXIS = ',24i4,/,(11x,24i4))
120 format ('  Y-AXIS = ',24i4,/,(11x,24i4))

  end subroutine print_decomp

end module atmosphere_mod

! And now one routine that is called back from the python side.
! It is placed outside the module, to make cross-language-calling a little more portable.
! When the python code exits or aborts, the FMS diag_manager buffers are lost.
! Here, we try to flush them before submitting to the crash.

subroutine abort_with_flush

  use diag_manager_mod, only: diag_manager_end
  use time_manager_mod, only: time_type, set_time

  implicit none

  type(time_type) :: t

  t = set_time(0, 999999)
  call diag_manager_end(t)
  return

end subroutine abort_with_flush
