#ifdef USE_LPJ
module lpj_driver_mod

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

! <CONTACT EMAIL="petri@pik-potsdam.de">
!   Stefan Petri
! </CONTACT> 

! <REVIEWER EMAIL="">
! 
! </REVIEWER>

! <OVERVIEW>
!   Wrapper for the vegetation model LPJ
! </OVERVIEW>

! <DESCRIPTION>
!   This module is the Fortran90 side of the wrapper for the
!   vegetation model LPJ. It wraps the replacement for the main program of LPJ.
! </DESCRIPTION>

!#include <fms_platform.h>

#define __ERROR__(message) \
call my_error(module_name,message,FATAL,thisfile,__LINE__)
#define __NOTE__(message) \
call my_error(module_name,message,NOTE,thisfile,__LINE__)
! CHECK is like ASSERT, but independent of NDEBUG
#define __CHECK__(x, message) \
if(.NOT.(x))call my_error(module_name,message,FATAL,thisfile,__LINE__)
#ifndef NDEBUG
#define __ASSERT__(x, message) \
if(.NOT.(x))call my_error(module_name,message,FATAL,thisfile,__LINE__)
#else
#define __ASSERT(x, message) ! empty
#endif

  use time_manager_mod, only: time_type, operator(+), operator(==), assignment(=), get_time, get_date, increment_time

  use mpp_domains_mod,  only: domain1D, domain2d, mpp_get_layout, mpp_define_layout, &
                              mpp_define_domains, mpp_get_compute_domain,            &
                              CYCLIC_GLOBAL_DOMAIN, mpp_get_compute_domains,         &
                              mpp_get_domain_components, mpp_get_pelist,             &
                              mpp_define_mosaic, mpp_global_field

  use fms_mod,          only: stdlog, open_namelist_file, check_nml_error,   &
                              write_version_number, error_mesg, FATAL, NOTE, &
                              mpp_pe, mpp_npes, mpp_root_pe, close_file,     &
                              mpp_clock_id, mpp_clock_begin, mpp_clock_end,  &
                              clock_flag_default, CLOCK_COMPONENT, CLOCK_SUBCOMPONENT, &
                              field_size, read_data, field_exist, get_mosaic_tile_grid

  use mpp_mod,          only: mpp_error, FATAL, mpp_get_current_pelist, stdout, mpp_chksum, mpp_sum

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

  use field_manager_mod,  only: MODEL_LAND
  use tracer_manager_mod, only: register_tracers, get_tracer_index, NO_TRACER

  use mosaic_mod,       only: get_mosaic_ntiles, calc_mosaic_grid_area, &
                              get_mosaic_xgrid_size, get_mosaic_xgrid
  
  use constants_mod,    only: PI, SECONDS_PER_DAY, KELVIN, RAD_TO_DEG, DEG_TO_RAD, &
       WTMAIR, & ! 28.96440 [AMU] molecular weight of air
       WTMC, &   ! 12.00000 [AMU] molecular weight of carbon
       WTMCO2    ! 44.00995 [AMU] molecular weight of carbon dioxide

  use land_types_mod,   only: land_data_type

  use numerics_mod,     only: is_latlon

  use topography_mod, only: get_topog_mean

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

implicit none

private

! ==== public interfaces =====================================================

public lpj_driver_init
public lpj_driver_update_fast ! , lpj_driver_update_slow
public update_lpj_bnd_fast, update_lpj_bnd_slow
public lpj_driver_end
!public lpj_driver_restart


! ==== some names, for information only ======================================
character(len=*),   parameter :: &
     module_name = 'lpj_driver_mod', &
     version     = '$Id$', &
     tagname     = '$Name$', &
     thisfile    = __FILE__
logical                       :: module_is_initialized = .FALSE.


! ==== public data type =====================================================

! <TYPE NAME="lpj_wrapper_type">
!   <DESCRIPTION>
!     The type describing the state of the LPJ model wrapper. There is only one
!     variable of this type, called "lpj_state", which is used in the model to
!     hold the data.
!     Strictly speaking, this type is not necessary, but there is no harm in
!     having it and it is possible to image a situation where a single entity
!     describing the lpj model wrapper state would be useful.
!   </DESCRIPTION>
type lpj_wrapper_type

   private  ! it's an opaque type: all data members are private

!   <DATA NAME="domain" TYPE="domain2d" DIM="2">
!     Domain of calculations
!   </DATA>
   type(domain2d)         :: domain     ! domain of calculations

   ! The update subroutines of the land model do not get the current time as arguments.
   ! Instead they must save the start time and time steps and advance the current time
   ! on each invocation.
   type(time_type) :: current_time, dt_fast, dt_slow
   integer         :: dt_fast_seconds, dt_slow_seconds  ! unpacked into seconds

   ! LPJ requires daily values. land_lad typically uses sub-daily time steps.
   ! Thus we aggregate here over one day's time steps
   integer :: yesterday = -1

   ! The following data arrays are on the local compute domain, unless noted otherwise.
   ! In these, the input data from land_lad is accumulated over the time steps of one day
   real, dimension(:,:), pointer :: &
        !lprec => NULL(),            & ! [kg/m2] == [mm] daily sum liquid precipitation
        !fprec => NULL(),            & ! [kg/m2] == [mm] daily sum frozen precipitation
        ! Currently, LPJ requires total precipitation and computes liquid and snow fractions
        ! on its own, depending on temperature.
        ! It is planned that LPJ can also make use of separate Xprec input directly.
        ! But that still needs to be implemented by th LPJ team.
        prec => NULL(),             & ! [kg/m2] == [mm] daily sum of lprec and fprec
        temp_mean => NULL(),        & ! [degC] daily mean
        temp_max => NULL(),         & ! [degC] daily max
        temp_min => NULL(),         & ! [degC] daily min
        swdown => NULL(),           & ! [W/m2] daily mean
        lwnet => NULL(),            & ! [W/m2] daily mean
        t_flux => NULL(),           & ! [W/m2] Sensible heat flux
        qflux => NULL(),            & ! [kg/m2/s] water vapor flux near surface
        dedq => NULL(),             & ! [kg/m2/s] Moisture flux humidity sensitivity
        dhdt => NULL(),             & ! derivative of sens over T
        drdt => NULL(),             & ! derivative of LW radiation over T
        drag_q => NULL(),           & ! [m/s] Product of drag coefficient for water vapor by wind
        p_surf => NULL(),           & ! [N/m2] == [Pa] == [kg/m/s2] Surface pressure
        wind => NULL(),             & ! [m/s] daily mean
        co2 => NULL()                 ! CO2 concentration [kg CO2/kg wet air]; Note: LPJ wants [ppm] daily mean dry volume mixing ratio
   ! TODO: lightning

   real, dimension(:,:), pointer :: lad_land_mask => NULL() ! land fraction, 0 = ocean
   real, dimension(:,:), pointer :: lad_warea => NULL() ! Area covered by water per grid cell

   real, dimension(:,:), pointer :: lpj_land_mask => NULL()

   ! LPJ output data fields
   real, dimension(:,:), pointer ::         &
        co2_flux_out => NULL(),         & ! daily co2 flux from LPJ to atmosphere [kg CO2/m2/sec]
        q_ca_out => NULL(),             & !
        runoff_out => NULL(),           & !
        roughness_length_out => NULL(), & !
        surface_temp_out => NULL(),     & !
        albedo_out => NULL(),           & !
        evap1_out => NULL(),            & !
        dedt_out => NULL(),             & !
        gc_out => NULL(),               & !
        exlpj_out => NULL()
   ! TODO: folia protective cover [%] 
   ! TODO: LAI 

   ! the local runoff_out maps are assembled into a global map
   ! TODO: possible optimisations:
   ! 1. each task stores only the coast cell information for those cells
   !    which pour out into ocean cells within its own compute domain.
   ! 2. instead of global runoff map, it might be sufficient
   !    to store in each task a data domain with a halo width corresponding
   !    to the maximum redirection distance used in subroutine
   !    prepare_coastal_runoff_routing(), currently 4.
   !    The FMS domains code has some flexibility in the halowidth,
   !    but it is yet unclear wether this really increases performance.
   real, dimension(:,:), pointer :: discharge => NULL() ! local map of runoff _into_ the ocean

   ! list of coast cells, global grid indices for runoff source and destination
   integer :: ncoastcells = 0
   integer, dimension(:), pointer :: &
        coast_src_i => NULL(), &
        coast_src_j => NULL(), &
        coast_dst_i => NULL(), &
        coast_dst_j => NULL()

   real, dimension(:,:), pointer :: inland_sinks => NULL()

end type lpj_wrapper_type
! </TYPE>


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

! ==== local module variables ================================================
type(lpj_wrapper_type),save    :: lpj_state  ! internal state of the model ! do we really need save here?

integer :: lpj_init_clock_id, lpj_update_fast_clock_id, lpj_update_slow_clock_id, lpj_update_cpl_clock_id
! wrapper function in fmsclocks.F90, to pull in the object file, so that it is available
! in the LPJ C interface
INTERFACE
   FUNCTION FMSCLOCK_ID(NAME)
     CHARACTER(*), INTENT(IN) :: NAME
     INTEGER(KIND=4) :: FMSCLOCK_ID
   END FUNCTION FMSCLOCK_ID
END INTERFACE

integer :: isphum ! number of specific humidity in the tracer array, or NO_TRACER
integer :: ico2   ! number of CO2 in the tracer array, or NO_TRACER

! cold start CO2 concentration [kg CO2 / kg wet air]
! 434.5e-06 [kg/kg] corresponds to 285.9 [ppmv] in dry air, i.e. q_ca = 0
real :: co2_init = 434.5e-06 ! [kg/kg wet air] also used as constant CO2 concentration if no CO2 tracer is defined in the field table

! Taken from land_lad2/canopy/cana_tile.F90
! The value 10.0 is set in the namelist of the ESM2M_pi-control_C2 example setup
real :: canopy_air_mass_for_tracers = 10.0 ! mass of wet air in the canopy air
                                           ! space for tracers other than water vapor, [kg wet air/m^2]

! Number of local and global grid cells in each dimensions.
! We require that FMS/LAD and LPJ run on the same grid resolution.
integer :: nlons, nlats, nlons_global, nlats_global
integer :: is, ie, js, je  ! local domain boundary indices

logical :: &
     do_lpj = .true.,                          & ! invoke LPJ at all?
     do_check_land_masks = .true.,             & ! set to false for production runs
     use_co2_from_lpj = .false.,               &
     use_q_ca_from_lpj = .false.,              &
     use_runoff_from_lpj = .false.,            &
     ! use_runoff_antarctica_from_lpj is obsolete
     use_runoff_antarctica_from_lpj = .false., & 
     use_roughness_length_from_lpj = .false.,  &
     use_surface_temp_from_lpj = .false.,      &
     use_albedo_from_lpj = .false.,            &
     enable_fpe = .false.


real    :: min_water_frac = 1e-4 ! min water area fraction for the point to be valid discharge destination
real    :: min_land_frac  = 1e-4 ! min land fraction for the cell to be considered            
!                                  "land" for discharge destination calculations

integer :: & ! diag manager field ids for input data on land_lad resolution
     id_lad_mask, &
     id_fprec, id_lprec, id_prec, &
     id_temp, id_temp_mean, id_temp_min, id_temp_max, &
     id_swdown, id_lwnet, id_t_flux, id_wind, id_qflux, id_dedq, id_dhdt, id_drdt,&
     id_co2, id_fco2, id_Dfco2Dq, id_co2_ppmv

integer :: & ! diag manager field ids for output data, on LPJ units and daily time step
     id_lpj_domainmask, id_lpj_mask, &
     id_carbon_flux_out, id_co2_flux_out, id_q_ca_out, id_roughness_length_out, id_runoff_out, id_surface_temp_out, id_albedo_out, id_evap1_out, id_dedt_out, id_gc_out, id_exlpj_out
! ... and on land_lad units and time step
integer :: id_co2_out_lad, id_q_ca_out_lad, id_roughness_length_out_lad, id_runoff_out_lad, id_surface_temp_out_lad, id_albedo_out_lad, id_discharge_out_lad, id_discharge_out_lad_aw


! <NAMELIST NAME="lpj_driver_nml">
!   <DATA NAME="use_co2_from_lpj" TYPE="logical" DEFAULT=".false.">
!     Override land_lad-produced data with LPJ-produced data
!   </DATA>
!   <DATA NAME="use_q_ca_from_lpj" TYPE="logical" DEFAULT=".false.">
!     Override land_lad-produced data with LPJ-produced data
!   </DATA>
!   <DATA NAME="use_runoff_from_lpj" TYPE="logical" DEFAULT=".false.">
!     Override land_lad-produced data with LPJ-produced data
!   </DATA>
!   <DATA NAME="use_runoff_antarctica_from_lpj" TYPE="logical" DEFAULT=".false.">
!    obsolete, just to keep old namelists valid.
!   </DATA>
!   <DATA NAME="use_roughness_length_from_lpj" TYPE="logical" DEFAULT=".false.">
!     Override land_lad-produced data with LPJ-produced data
!   </DATA>
!   <DATA NAME="use_surface_temp_from_lpj" TYPE="logical" DEFAULT=".false.">
!     Override land_lad-produced data with LPJ-produced data
!   </DATA>
!   <DATA NAME="use_albedo_from_lpj" TYPE="logical" DEFAULT=".false.">
!     Override land_lad-produced data with LPJ-produced data
!   </DATA>
!   <DATA NAME="canopy_air_mass_for_tracers" TYPE="real" DEFAULT="10.0" UNIT="kg/m2">
!    Mass of wet air in the canopy air space for tracers other than water vapor.
!    The default value 10.0 is taken from the ESM2M_pi-control_C2 example setup input.nml
!   </DATA>
!   <DATA NAME="co2_init" TYPE="real" DEFAULT="434.5e-06" UNIT="kg/kg wet air">
!    Initial CO2 concentration.
!    Unfortunately, this is not passed from the atmosphere.
!    434.5e-06 [kg/kg] corresponds to 285.9 [ppmv] in dry air.
!   </DATA>
!   <DATA NAME="min_water_frac" TYPE="real" DEFAULT="1e-4">
!     Minimum water area fraction for the point to be a
!     valid discharge destination
!   </DATA>
!   <DATA NAME="min_land_frac" TYPE="real" DEFAULT="1e-4">
!     "Land" for discharge destination calculations
!   </DATA>

namelist /lpj_driver_nml/ &
     do_lpj, &
     use_co2_from_lpj, &
     use_q_ca_from_lpj, &
     use_runoff_from_lpj, &
     use_runoff_antarctica_from_lpj, &
     use_roughness_length_from_lpj, &
     use_surface_temp_from_lpj, &
     use_albedo_from_lpj, &
     canopy_air_mass_for_tracers, &
     co2_init, &
     enable_fpe, &
     min_water_frac, &
     min_land_frac
! </NAMELIST>



contains

  ! ============================================================================
  ! <SUBROUTINE NAME="lpj_driver_init">
  !   <OVERVIEW>
  !     Initializes the LPJ driver
  !   </OVERVIEW>
  !   <DESCRIPTION>
  !     ...
  ! Sigh. land_lad uses two very similar named variables,
  ! n_tiles and ntiles,
  ! but to my understanding they are indeed very independant from
  ! each other. One seems to reflect the grid type (1 = spherical lat/lon,
  ! 6 = so-called cubed sphere, else unknown).
  ! The other has yet unknown meaning.
  ! Here we require that both must be 1.
  !   </DESCRIPTION>
  !   <PUBLICROUTINE>
  subroutine lpj_driver_init &
       ( time_init, time, dt_fast, dt_slow, &
       lad_domain, lad_layout, &
       glon, glat, &
       n_tiles, ntiles, &
       lad_id_lon, lad_id_lat, &
       lad_blon, lad_blat, lad_garea, lad_gfrac, isphum_in, ico2_in, &
       bnd)
    !USE, INTRINSIC :: IEEE_EXCEPTIONS  
    !type(atmos_land_boundary_type), intent(inout) :: atmos2land ! land boundary data
    type(time_type), intent(in)    :: time_init    ! initial time of simulation (?)
    type(time_type), intent(in)    :: time         ! current time
    type(time_type), intent(in)    :: dt_fast      ! fast time step
    type(time_type), intent(in)    :: dt_slow      ! slow time step
    type(domain2d),  intent(in)    :: lad_domain ! domain of computations
    integer, dimension(1:2), intent(in) :: lad_layout ! PE layout of lad domain
    real, intent(in)               :: glon(:,:), glat(:,:) ! global cell centers, used for mapping between LPJ and FMS domains
    integer, intent(in)            :: n_tiles, ntiles ! sigh. am I confused, or was the land_lad programmer? S.P.
    integer, intent(in)            :: lad_id_lon, lad_id_lat ! diag_manager axis ids for land_lad coordinates
    real, intent(in)               :: lad_blon(:,:)   ! lon corners of the local grid cells in radians
    real, intent(in)               :: lad_blat(:,:)   ! lat corners of the local grid cells in radians
    real, intent(in)               :: lad_garea(:,:)  ! global cell area
    real, intent(in)               :: lad_gfrac(:,:)  ! global map of fraction of land in the grid cell, 0 for ocean - or should we use land area instead?
    integer, intent(in)            :: isphum_in, ico2_in ! index in tracer array for specific humidity and CO2, or NO_TRACER
    type(land_data_type), intent(in) :: bnd    ! land boundary data, to share the initial conditions from land_lad
    !   </PUBLICROUTINE>

    integer :: i,j,k, unit, ierr, io, used

    ! some magical values for testing the interface between Fortran90 and C++ code.
    ! we pass them from Fortran90 to C++ functions and test in there if they were
    ! received correctly.
    integer, parameter :: deadbeef = X'eadbeef' ! X'deadbeef' yields overflow error on 32bit machines
    integer, parameter :: fse = 4711
    real,    parameter :: onetwothree = 1234567809.0987654321D0, &
         pitest = PI ! from constants_mod ! 3.1415926535897932384626433832795029D0 ! from linux <math.h>
    logical :: true = .true., false = .false.
    integer :: kind_i = KIND(deadbeef), kind_r = KIND(pitest), kind_l = KIND(true)
    real, dimension(2,3,4) :: testarray3d

    integer :: atmos_comm ! the MPI communicator for the atmosphere
    integer,allocatable,dimension(:) :: lad_pelist ! needed to obtain the peset

    integer :: Time_init_year, Time_init_month, Time_init_day, Time_init_hour, Time_init_minute, Time_init_second
    integer :: Time_year, Time_month, Time_day, Time_hour, Time_minute, Time_second
    !integer :: Time_step_days, Time_step_seconds
    
    integer :: log_unit
    integer :: pe, npes

    integer :: lad_nlons, lad_nlats ! number of cells in lad resolution
    real, dimension(:), allocatable :: lad_lons, lad_lats ! lad cell center coordinates

    ! info for assembling the global runoff routing map
    integer :: global_number_lpj_cells = 0 ! number of LPJ cells (land cells)
    real, allocatable, dimension(:) :: lpj_cell_lons, lpj_cell_lats
    integer, allocatable, dimension(:) :: lpj_cell_next

    integer :: id_runoff_map = -1, id_runoff_coastcells = -1
    integer :: id_runoff_dest_lon = -1, id_runoff_dest_lat = -1
    integer :: id_runoff_dcount = -1, id_runoff_basin = -1
    integer :: id_runin_coastcells = -1
    integer :: id_inland_sinks = -1  ! map of inland sink cells
    !! Water content of all rivers is written to lpj_runoff.
    !! Water content of inland sinks can be computed by
    !! applying inland_sinks as mask to lpj_runoff.
    
    integer :: id_warea = -1

    ! just for diagnostic output, as helper to prepare runoff routing map
    integer :: id_horo = -1
    real, allocatable, dimension(:,:) :: HORO  ! Orography height, averaged over grid cell

    integer :: id_lad_domainmask
    real, allocatable, dimension(:,:) :: lad_domainmask

    integer(4) :: previous_fpe_flags
    !CALL IEEE_SET_HALTING_MODE(IEEE_DIVIDE_BY_ZERO, .FALSE.)
    if (module_is_initialized) return

    !----- write version and namelist to log file -----
    call write_version_number ( version, tagname )
    write(*,*) 'lpj_driver_init', version
    unit = open_namelist_file()
    ierr=1;
    do while (ierr /= 0)
       read  (unit, nml=lpj_driver_nml, iostat=io, end=10)
       ierr = check_nml_error(io,'lpj_driver_nml') ! also initializes nml error codes
    enddo
10  call close_file (unit)

    if ( mpp_pe() == mpp_root_pe() ) then
       log_unit = stdlog()
       write (log_unit, nml=lpj_driver_nml)
       write (*, nml=lpj_driver_nml)
    endif

    ! for preparation of LPJ configuration we need the topography map on
    ! LPJ resolution and land-sea-mask even before we can actually start and use LPJ.
    ! Thus write it out here, before checking for do_lpj .

    call mpp_get_compute_domain(lad_domain, is, ie, js, je)

    id_horo = register_static_field('lpj', "horo", (/lad_id_lon,lad_id_lat/), "Orographic Height from FMS", "m")
    if (id_horo > 0) then
       allocate(HORO(ie-is+1, je-js+1))
       write(*,*) 'for get_topog_mean(): shape(HORO)', shape(HORO)
       write(*,*) '  vs. (/size(lad_blon,1)-1,size(lad_blon,2)-1/): ', size(lad_blon,1)-1, size(lad_blon,2)-1
       used = get_topog_mean(lad_blon, lad_blat, HORO)
       if (.not. used) &
            call error_mesg (module_name, 'could not get topography mean HORO', FATAL)
       used = send_data(id_horo, HORO, Time)
       deallocate(HORO)
    end if

    if (.not. do_lpj) then
       write (*,*) 'do_lpj = .false. -- Not invoking LPJ!'
       return
    endif

    if (enable_fpe) then
       !! 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)
    endif

    lpj_init_clock_id = mpp_clock_id( 'LPJ_init', flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT )
    lpj_update_slow_clock_id = mpp_clock_id( 'LPJ_update_slow', flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT )
    lpj_update_fast_clock_id = mpp_clock_id( 'LPJ_update_fast', flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT )
    lpj_update_cpl_clock_id = fmsclock_id( 'LPJ_update_cpl' ) ! to force the linker to load the fmsclock wrapper object file

    call mpp_clock_begin(lpj_init_clock_id)

    if (n_tiles /= 1) then
       call error_mesg('lpj_driver_init', 'n_tiles /= 1', FATAL)
    endif
    if (ntiles /= 1) then
       call error_mesg('lpj_driver_init', 'ntiles /= 1', FATAL)
    endif

    if (.not. is_latlon(lad_blon, lad_blat)) then
       call error_mesg('lpj_driver_init', 'FMS grid is not lat-lon type', FATAL)
    endif

    ! TODO: collect LPJ parameters from namelist and pass them to C-side of LPJ driver

    pe = mpp_pe()
    npes = mpp_npes()
    allocate(lad_pelist(npes))
    call mpp_get_pelist(lad_domain, lad_pelist)
    call mpp_get_current_pelist(lad_pelist, commID=atmos_comm)
    deallocate(lad_pelist)

    lpj_state%domain = lad_domain

    lad_nlons = size(lad_blon, 1)-1 ! local domain
    lad_nlats = size(lad_blat, 2)-1

    !write(*,*) 'lad_blon lbound', lbound(lad_blon), 'ubound', ubound(lad_blon)
    !write(*,*) 'lad_blat lbound', lbound(lad_blat), 'ubound', ubound(lad_blat)
    write(*,'(a,i4,a,i4,a,i4,a,i4,a,i5,a,i5)') 'is ', is, ' ie ', ie, ' js ', js, ' je ', je, ' lad_nlons ', lad_nlons, ' lad_nlats ', lad_nlats
    if (ie-is+1 /= lad_nlons) &
         call error_mesg('lpj_driver_init', 'dimension of land_lad lon boundaries does not match boundary indices', FATAL)
    if (je-js+1 /= lad_nlats) &
         call error_mesg('lpj_driver_init', 'dimension of land_lad lat boundaries does not match boundary indices', FATAL)
    write(*,*) 'lad_gfrac bound ',lbound(lad_gfrac),ubound(lad_gfrac)

    allocate(lad_lons(is:ie), lad_lats(js:je))
    ! sigh. The array boundary values are lost during the subroutine call.
    ! Thus, lad_blon and lad_blat always start at 1 here, even though with domain decomposition
    ! they were allocated with different boundaries by our caller.
    lad_lons = ((lad_blon(1:lad_nlons,1)+lad_blon(2:lad_nlons+1,1))/2)*RAD_TO_DEG ! lad cell center longitudes in degrees
    lad_lats = ((lad_blat(1,1:lad_nlats)+lad_blat(1,2:lad_nlats+1))/2)*RAD_TO_DEG ! lad cell center latitudes in degrees

    lpj_state%current_time = Time
    lpj_state%dt_fast = dt_fast
    lpj_state%dt_slow = dt_slow
    
    call get_date(Time_init, Time_init_year, Time_init_month, Time_init_day, &
         Time_init_hour, Time_init_minute, Time_init_second)
    call get_date(Time, Time_year, Time_month, Time_day, &
         Time_hour, Time_minute, Time_second)
    call get_time(dt_fast, lpj_state%dt_fast_seconds)
    !lpj_state%dt_fast_seconds = Time_step_seconds
    call get_time(dt_slow, lpj_state%dt_slow_seconds)
    !lpj_state%dt_slow_seconds = Time_step_seconds
    
    allocate(                                      & ! data fields 
         !lpj_state%lprec(is:ie,js:je),            &
         !lpj_state%fprec(is:ie,js:je),            &
         lpj_state%prec(is:ie,js:je),              &
         lpj_state%temp_mean(is:ie,js:je),         &
         lpj_state%temp_max(is:ie,js:je),          &
         lpj_state%temp_min(is:ie,js:je),          &
         lpj_state%swdown(is:ie,js:je),            &
         lpj_state%lwnet(is:ie,js:je),             &
         lpj_state%t_flux(is:ie,js:je),            &
         lpj_state%qflux(is:ie,js:je),             &
         lpj_state%dedq(is:ie,js:je),              &
         lpj_state%dhdt(is:ie,js:je),             &
         lpj_state%drdt(is:ie,js:je),              &
         lpj_state%drag_q(is:ie,js:je),            &
         lpj_state%p_surf(is:ie,js:je),            &
         lpj_state%wind(is:ie,js:je),              &
         lpj_state%co2(is:ie,js:je),               &
         !lpj_state%fco2(is:ie,js:je),             &
         !lpj_state%Dfco2Dq(is:ie,js:je)           &
         STAT=ierr)

    if (0 /= ierr) &
         call error_mesg('lpj_driver_init', 'allocate input fields failed', FATAL)

    allocate(                                      &
         lpj_state%co2_flux_out(is:ie,js:je),      &
         lpj_state%q_ca_out(is:ie,js:je),          &
         lpj_state%runoff_out(is:ie,js:je),        &
         lpj_state%roughness_length_out(is:ie,js:je), &
         lpj_state%surface_temp_out(is:ie,js:je),  &
         lpj_state%albedo_out(is:ie,js:je),        &
         lpj_state%evap1_out(is:ie,js:je),         &
         lpj_state%dedt_out(is:ie,js:je),          &
         lpj_state%gc_out(is:ie,js:je),            &
         lpj_state%exlpj_out(is:ie,js:je),         &
         STAT=ierr)

    if (0 /= ierr) &
         call error_mesg('lpj_driver_init', 'allocate output fields failed', FATAL)

    ! TODO: check which of the lpj_state variables need to be saved in RESTART files

    ! Make sure that we have nice values also over ocean and antarctica!
    ! else netcdf output will complain bitterly.
    lpj_state%co2_flux_out = 0
    lpj_state%q_ca_out = 0
    lpj_state%runoff_out = 0
    lpj_state%roughness_length_out = bnd%rough_mom(:,:,1)
    !lpj_state%surface_temp_out = 288.0 - KELVIN ! default initial condition from soil.F90 ???
    lpj_state%surface_temp_out = bnd%t_surf(:,:,1) ! copy the initial conditions that were set in soil.F90
    lpj_state%albedo_out = 0
    lpj_state%evap1_out = 0
    lpj_state%dedt_out = 0
    lpj_state%gc_out = 0
    lpj_state%exlpj_out = 0


    allocate(lpj_state%lad_land_mask(is:ie, js:je))
    lpj_state%lad_land_mask = lad_gfrac(is:ie,js:je)
    allocate(lpj_state%lad_warea(is:ie, js:je))
    lpj_state%lad_warea = lad_garea(is:ie,js:je) * (1-lad_gfrac(is:ie,js:je))

    ! now prepare to call out into the C - implemented heart of LPJ
    ! first assemble some parameters with well-known values to verify
    ! that the F90<->C++ interface is OK.
    do i=1,2
       do j=1,3
          do k=1,4
             testarray3d(i,j,k) = i*100.0+j+k/10.0
          end do
       end do
    end do
    ! The arrays in lpj_state%XXX are allocate()d, they might behave different than
    ! testarray3d. We put remarkable values at the corners of such array to test
    ! that we can find those values on the C side in the expected locations.
    lpj_state%co2(is,js) = 42.42
    lpj_state%co2(ie,js) = 43.43
    lpj_state%co2(is,je) = 44.44
    lpj_state%co2(ie,je) = 45.45
    !write(*,*) 'co2 test array', lpj_state%co2

    write(*,*) 'MAXEXPONENT ', MAXEXPONENT(pitest)

    ! Get back from LPJ information about its grid resolution.
    ! TODO: Check if it matches the land_lad grid resolution.
    ! We pass the boundary box (cell borders) of the FMS compute domain to LPJ.
    ! On the LPJ side we then can compute a grid in LPJ-resolution inside those boxes,
    ! and can compute the mappings between LPJ grid cell indices and FMS lat-lon-array indices.
    ! This slightly complicated approach is necessary because in LPJ only the land cells
    ! are represented, while the boundary box of the FMS domain might also contain ocean-only
    ! cells along its borders. 
    call lpj_init(         &
         ! first some arguments for checking the F90->C++ interface
         deadbeef, pitest, &
         true, false,      &
         testarray3d,      &
         lpj_state%co2,    &
         kind_i, kind_r, kind_l,               &
         ! now the ``real'' arguments
         ! the FMS time stamps are broken out into days and seconds
         Time_init_year, Time_init_month, Time_init_day,     &
         Time_init_hour, Time_init_minute, Time_init_second, &
         Time_year, Time_month, Time_day,                    &
         Time_hour, Time_minute, Time_second,                &
         lpj_state%dt_fast_seconds, lpj_state%dt_slow_seconds, &
         ! information about domain partitioning
         atmos_comm, pe, mpp_root_pe(), npes, &
         lad_layout, & ! of domain decomposition, number of PEs in lon and lat directions
         ! cell centers of corners of global grid
         glon(1,1)*RAD_TO_DEG, glon(size(glon,1),size(glon,2))*RAD_TO_DEG, &
         glat(1,1)*RAD_TO_DEG, glat(size(glat,1),size(glat,2))*RAD_TO_DEG, &
         ! the boundary indices of the local compute domain
         is, ie, js, je,   &
         ! grid axis info of local domain
         lad_nlons, lad_nlats,     &
         !!blon, blat,       & pass only 1-D arrays here because we require rectilinear lat-lon grid
         lad_blon(1,1)*RAD_TO_DEG, lad_blon(lad_nlons+1,lad_nlats+1)*RAD_TO_DEG, &
         lad_blat(1,1)*RAD_TO_DEG, lad_blat(lad_nlons+1,lad_nlats+1)*RAD_TO_DEG, & ! boundary box of FMS compute domain
         !lons, lats,       &
         lpj_update_cpl_clock_id, &
         ! LPJ-specific options TODO
         1, & !argc,
         !argv, &
         ! output parameters
         global_number_lpj_cells, &
         ! finally interface safety parameters.
         onetwothree, fse &
         )

    
    ! now, after the co2 array has been (ab)used for the F90->C interface test,
    ! we can set the initial value for the CO2 concentration
    lpj_state%co2 = co2_init ! [kg/kg wet air]

    isphum = isphum_in
    ico2 = ico2_in
    if (ico2 /= NO_TRACER) bnd%tr(:,:,1,ico2) = co2_init ! remains [kg/kg wet air]
    
    ! get information to set up the horizontal interpolation
    ! (not really needed until we allow different resolution between land_lad and LPJ)
    nlons = lad_nlons
    nlats = lad_nlats
    nlons_global = size(glon,1)
    nlats_global = size(glon,2)

    write(*,'(a,i4,a,i4,a,i4,a,i4)') 'nlons ',nlons,' nlats ',nlats, &
         ' nlons_global ',nlons_global,' nlats_global ',nlats_global

    allocate(lpj_state%lpj_land_mask(is:ie,js:je))

    allocate(lpj_state%discharge(is:ie,js:je), STAT=ierr)
    if (0 /= ierr) &
         call error_mesg('lpj_driver_init', 'allocate discharge failed', FATAL)
    lpj_state%discharge = 0.0

    ! The arrays in lpj_state%XXX are allocate()d, they might behave different than
    ! testarray3d. We put remarkable values at the corners of such array to test
    ! that we can find those values on the C side in the expected locations.
    lpj_state%lpj_land_mask(is,js) = 4242.42
    lpj_state%lpj_land_mask(ie,js) = 4343.43
    lpj_state%lpj_land_mask(is,je) = 4444.44
    lpj_state%lpj_land_mask(ie,je) = 4545.45

    allocate(lpj_cell_lons(global_number_lpj_cells))
    allocate(lpj_cell_lats(global_number_lpj_cells))
    allocate(lpj_cell_next(global_number_lpj_cells))
    lpj_cell_lons = 0;
    lpj_cell_lats = 0;
    lpj_cell_next = 0;

    !allocate(lpj_state%global_gfrac(nlons_global, nlats_global))
    !lpj_state%global_gfrac = lad_gfrac

    call lpj_init_grid_interpolation(  &
         ! output parameters
         lpj_state%lpj_land_mask, &
         lpj_cell_lons, &
         lpj_cell_lats, &
         lpj_cell_next, &
         ! interface safety parameters
         fse &
         )

    !if (do_check_land_masks) & ! here do it unconditionally
    call check_land_masks(lpj_state%lad_land_mask, lpj_state%lpj_land_mask, 'LPJ', 'land_lad', 'LPJ')

    id_lad_mask = register_static_field('lpj', 'lad_land_mask', &
         (/lad_id_lon,lad_id_lat/), &
         long_name='land fraction as obtained from land_lad (gfrac)', &
         missing_value = 0.0 &
         )
    id_lad_domainmask = register_static_field('lpj', 'lad_domain_mask', &
         (/lad_id_lon,lad_id_lat/), long_name='MPI task numbers of lad cells', &
         units="lad task number", missing_value = -1.0 &
         )
    id_lpj_domainmask = register_static_field('lpj', 'lpj_domain_mask', &
         (/lad_id_lon,lad_id_lat/), long_name='MPI task numbers of LPJ cells', &
         units="lpj task number", missing_value = -1.0 &
         )
    id_lpj_mask = register_static_field('lpj', 'lpj_land_mask', &
         (/lad_id_lon,lad_id_lat/), long_name='land mask as obtained from LPJ', &
         missing_value = -1.0 &
         )
    id_warea = register_static_field('lpj', 'lad_warea', &
         (/lad_id_lon,lad_id_lat/), long_name='ocean area per grid cell', &
         units = 'm^2', missing_value = -1.0 &
         )
    ! Note that the static fields are actually written out into the file only in
    ! diag_manager_end()
    if (id_lad_mask > 0) &
         used = send_data(id_lad_mask, lpj_state%lad_land_mask)
    if (id_lad_domainmask > 0) then
       allocate(lad_domainmask(is:ie,js:je))
       lad_domainmask(:,:) = pe
       used = send_data(id_lad_domainmask, lad_domainmask)
       deallocate(lad_domainmask)
    end if
    if (id_lpj_domainmask > 0) &
         used = send_data(id_lpj_domainmask, lpj_state%lpj_land_mask)
    if (id_lpj_mask > 0) then
       where (lpj_state%lpj_land_mask > 0)
          lpj_state%lpj_land_mask = 1
       elsewhere
          lpj_state%lpj_land_mask = 0
       end where
       used = send_data(id_lpj_mask, lpj_state%lpj_land_mask)
    end if
    if (id_warea > 0) &
         used = send_data(id_warea, lpj_state%lad_warea)

    id_runoff_map = register_static_field('lpj', 'lpj_runoffmap', &
         (/lad_id_lon,lad_id_lat/), &
         long_name='runoff routing map (destination cell index)', &
         missing_value = -10.0 & ! FMS ocean
         )

    id_runoff_dest_lon = register_static_field('lpj', 'lpj_runoff_dest_lon', &
         (/lad_id_lon,lad_id_lat/), &
         long_name='runoff routing map (destination cell longitude)', &
         units = 'degrees_E', &
         missing_value = -1000.0 & ! FMS ocean
         )
    id_runoff_dest_lat = register_static_field('lpj', 'lpj_runoff_dest_lat', &
         (/lad_id_lon,lad_id_lat/), &
         long_name='runoff routing map (destination cell latitude)', &
         units = 'degrees_N', &
         missing_value = -1000.0 & ! FMS ocean
         )

    id_runoff_coastcells = register_static_field('lpj', 'lpj_coastcells', &
         (/lad_id_lon,lad_id_lat/), &
         long_name='land coastcells which route into ocean cells', &
         missing_value = 0.0 &
         )

    id_runoff_dcount = register_static_field('lpj', 'lpj_dcount', &
         (/lad_id_lon,lad_id_lat/), &
         long_name='count of cells which route into this destination point', &
         missing_value = 0.0 &
         )
    id_runoff_basin = register_static_field('lpj', 'lpj_basin', &
         (/lad_id_lon,lad_id_lat/), &
         long_name='river basins', &
         missing_value = 0.0 &
         )

    id_runin_coastcells = register_static_field('lpj', 'lpj_beachcells', &
         (/lad_id_lon,lad_id_lat/), &
         long_name='ocean coastcells showing count of cells which route into this destination point', &
         missing_value = -1.0 &
         )

    id_inland_sinks = register_static_field('lpj', 'lpj_inland_sinks', &
         (/lad_id_lon,lad_id_lat/), &
         long_name='inland river sinks', &
         missing_value = 0.0 &
         )

    call prepare_coastal_runoff_routing(global_number_lpj_cells, &
         lpj_cell_lons, lpj_cell_lats, lpj_cell_next, &
         lad_gfrac, &
         id_runoff_map, id_runoff_coastcells, &
         id_runoff_dest_lon, id_runoff_dest_lat, &
         id_runoff_dcount, id_runoff_basin, &
         id_runin_coastcells, id_inland_sinks)

    deallocate(lpj_cell_lons)
    deallocate(lpj_cell_lats)
    deallocate(lpj_cell_next)

    !!! to get the above-written static fields out into the file now,
    !!! and not only after the entire program ended successfully.
    !call diag_manager_end(Time) !!!
    !flush(stdout())
    

    !! init diagnostic output

    ! input data on land_lad grid
    id_prec = register_diag_field('lpj', 'prec', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='daily sum of liquid and frozen precipitation', &
         !standard_name='', &
         units='mm', missing_value=-999.0 &
         )
    id_temp_mean = register_diag_field('lpj', 'temp_mean', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='daily mean temperature', &
         !standard_name='', &
         units='degC', missing_value=-999.0 &
         )
    id_temp_max = register_diag_field('lpj', 'temp_max', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='daily maximum temperature', &
         !standard_name='', &
         units='degC', missing_value=-999.0 &
         )
    id_temp_min = register_diag_field('lpj', 'temp_min', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='daily minimum temperature', &
         units='degC', missing_value=-999.0 &
         )
    id_swdown = register_diag_field('lpj', 'swdown', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='daily mean shortwave downwelling radiation', &
         !standard_name='', &
         units='W/m2', missing_value=-999.0 &
         )
    id_lwnet = register_diag_field('lpj', 'lwnet', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='daily mean longwave net radiation', &
         !standard_name='', &
         units='W/m2', missing_value=-999.0 &
         )
    id_t_flux = register_diag_field('lpj', 't_flux', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='daily mean sensible heat flux', &
         !standard_name='', &
         units='W/m2', missing_value=-999.0 &
         )
    id_qflux = register_diag_field('lpj', 'qflux', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='water vapor flux near surface', &
         !standard_name='', &
         units='kg/m2/s', missing_value=-999.0 &
         )
    id_dedq = register_diag_field('lpj', 'dedq', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='Moisture flux humidity sensitivity', &
         !standard_name='', &
         units='kg/m2/s', missing_value=-999.0 &
         )
    id_dhdt = register_diag_field('lpj', 'dhdt', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='derivative of sens over T', &
         !standard_name='', &
         units='kg/m2/s', missing_value=-999.0 &
         )
    id_drdt = register_diag_field('lpj', 'drdt', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='derivative of LW radiation over T', &
         !standard_name='', &
         units='kg/m2/s', missing_value=-999.0 &
         )
    id_wind = register_diag_field('lpj', 'wind', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='daily mean absolute wind speed', &
         !standard_name='', &
         units='m/s', missing_value=-999.0 &
         )
    id_fco2 = register_diag_field('lpj', 'fco2', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='daily mean atmospheric CO2 flux', &
         !standard_name='', &
         units='kg CO2/(m2 s)', missing_value=-999.0 &
         )
    id_Dfco2Dq = register_diag_field('lpj', 'Dfco2Dq', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='daily mean atmospheric CO2 d(tracer flux)/d(surf tracer)', &
         !standard_name='', &
         units='kg wet air/(m2 s)', missing_value=-999.0 &
         )
    id_co2 = register_diag_field('lpj', 'co2', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='daily integrated atmospheric CO2 concentration', &
         !standard_name='', &
         units='kg CO2/kg wet air', missing_value=-999.0 &
         )
    id_co2_ppmv = register_diag_field('lpj', 'co2_ppmv', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='daily inegrated atmospheric CO2 concentration', &
         !standard_name='', &
         units='ppmv', missing_value=-999.0 &
         )

    ! daily output data in LPJ units
    id_carbon_flux_out = register_diag_field('lpj', 'carbon_flux_out', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='LPJ daily carbon flux to atmosphere output', &
         !standard_name='', &
         units='g C/(m^2 day)', missing_value=-999.0 &
         )
    id_co2_flux_out = register_diag_field('lpj', 'co2_flux_out', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='LPJ daily CO2 flux to atmosphere output', &
         !standard_name='', &
         units='kg CO2/(m^2 sec)', missing_value=-999.0 &
         )
    id_q_ca_out = register_diag_field('lpj', 'q_ca_out', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='surface/canopy humidity', &
         !standard_name='', &
         units='kg/kg', missing_value=-999.0 &
         )
    id_roughness_length_out = register_diag_field('lpj', 'roughness_length_out', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='LPJ roughness length output', &
         !standard_name='', &
         units='??', missing_value=-999.0 &
         )
    id_runoff_out = register_diag_field('lpj', 'runoff_out', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='LPJ runoff', &
         !standard_name='', &
         units='kg/s', missing_value=-999.0 &
         )
    id_surface_temp_out = register_diag_field('lpj', 'surface_temp_out', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='LPJ surface temperature output', &
         !standard_name='', &
         units='degC', missing_value=-999.0 &
         )
    id_albedo_out = register_diag_field('lpj', 'albedo_out', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='LPJ albedo output', &
         !standard_name='', &
         units='1', missing_value=-999.0 &
         )
    id_evap1_out = register_diag_field('lpj', 'evap1_out', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='LPJ evapotranspiration output', &
         !standard_name='', &
         units='kg/m2/s', missing_value=-999.0 &
         )   
    id_dedt_out = register_diag_field('lpj', 'dedt_out', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='LPJ evapotranspiration output', &
         !standard_name='', &
         units='kg/m2/s', missing_value=-999.0 &
         )
    id_gc_out = register_diag_field('lpj', 'gc_out', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='LPJ evapotranspiration output', &
         !standard_name='', &
         units='kg/m2/s', missing_value=-999.0 &
         )
    id_exlpj_out = register_diag_field('lpj', 'exlpj_out', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='LPJ evapotranspiration output', &
         !standard_name='', &
         units='kg/m2/s', missing_value=-999.0 &
         )


    ! output data interpolated onto sub-daily time steps and converted to land_lad units
    ! should be output into a different file than above
    id_co2_out_lad = register_diag_field('lpj_lad', 'co2_out', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='LPJ CO2 output on land_lad resolution', &
         !standard_name='', &
         units='kg/kg wet air', missing_value=-999.0 &
         )
    id_q_ca_out_lad = register_diag_field('lpj_lad', 'q_ca_out', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='surface/canopy humidity', &
         !standard_name='', &
         units='kg/kg', missing_value=-999.0 &
         )
    id_roughness_length_out_lad = register_diag_field('lpj_lad', 'roughness_length_out', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='LPJ roughness length output on land_lad resolution', &
         !standard_name='', &
         units='??', missing_value=-999.0 &
         )
    id_runoff_out_lad = register_diag_field('lpj_lad', 'runoff_out', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='LPJ runoff on land_lad resolution', &
         !standard_name='', &
         units='kg/s', missing_value=-999.0 &
         )
    id_surface_temp_out_lad = register_diag_field('lpj_lad', 'surface_temp_out', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='LPJ surface temperature output on land_lad resolution', &
         !standard_name='', &
         units='K', missing_value=-999.0 &
         )
    id_albedo_out_lad = register_diag_field('lpj_lad', 'albedo_out', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='LPJ albedo output on land_lad resolution', &
         !standard_name='', &
         units='1', missing_value=-999.0 &
         )

    id_discharge_out_lad = register_diag_field('lpj_lad', 'discharge_out', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='LPJ runoff discharged into the ocean on land_lad resolution', &
         !standard_name='', &
         units='kg/s', missing_value=-999.0 &
         )
    id_discharge_out_lad_aw = register_diag_field('lpj_lad', 'discharge_out_aw', (/lad_id_lon,lad_id_lat/), Time, &
         long_name='LPJ area-weighted runoff discharged into the ocean on land_lad resolution', &
         !standard_name='', &
         units='kg/m2/s', missing_value=-999.0 &
         )

    deallocate(lad_lons, lad_lats)

    !{ begin init_output
    !
    ! Init daily output from values that are saved in the LPJ restart file.
    ! We must pass back meaningful values to land_lad/FMS already in the first time step,
    ! even before the invocation of lpj_update.
    ! This can be done only after all else has been initialised above.
    ! Note that the same code for writing out diagnostics and units conversion
    ! appears in lpj_update below.
    call lpj_init_output( &
         ! get back output from LPJ
         lpj_state%co2_flux_out, & ! returned as carbon flux in [g C/ (m^2/day)], must be converted to [kg CO2/(m^2 sec)]
         lpj_state%q_ca_out, &
         lpj_state%runoff_out, &       ! returned in [m3/s]
         lpj_state%roughness_length_out, &
         lpj_state%surface_temp_out, & ! returned in [degC]
         lpj_state%albedo_out, &
         lpj_state%evap1_out, &
         lpj_state%dedt_out, &
         lpj_state%gc_out, &
         lpj_state%exlpj_out, &

         ! input from FMS to LPJ
         lpj_state%inland_sinks, &

         ! size of the arrays
         size(lpj_state%q_ca_out, 1), &
         size(lpj_state%q_ca_out, 2), &
         ! interface safety parameters
         fse &
         )

    deallocate(lpj_state%inland_sinks)
    lpj_state%inland_sinks => NULL()

    ! experimental, to prevent instabilities
    where (lpj_state%q_ca_out>0.4)
       lpj_state%q_ca_out=0.4
    endwhere
    where (lpj_state%q_ca_out<0)
       lpj_state%q_ca_out=0
    endwhere

    write(*,'(a,i4,i4,a,i3,i3,i4,i4)') 'lpj_init_output returned size q_ca_out ', size(lpj_state%q_ca_out,1), size(lpj_state%q_ca_out,2),' bounds q_ca_out ', lbound(lpj_state%q_ca_out,1), lbound(lpj_state%q_ca_out,2), ubound(lpj_state%q_ca_out,1), ubound(lpj_state%q_ca_out,2)
    write(*,'(a,f8.3,a,f8.3,a,f8.3,a,Z20)') '   q_ca_out min ',minval(lpj_state%q_ca_out(:,:)),' max ',maxval(lpj_state%q_ca_out(:,:)),' mean ',sum(lpj_state%q_ca_out(:,:))/size(lpj_state%q_ca_out(:,:)), ' checksum', mpp_chksum(lpj_state%q_ca_out(:,:))

    call pour_into_ocean ! for runoff -> discharge

    !
    ! LPJ output
    !
    !if (id_carbon_flux_out > 0) & ! [g Carbon/(m^2*day)]
    !     used = send_data(id_carbon_flux_out, lpj_state%co2_flux_out, lpj_state%current_time, &
    !     rmask=lpj_state%lpj_land_mask)
    write(*,*) 'sending to id_q_ca_out ', id_q_ca_out, ' at ', Time_year, Time_month, Time_day, Time_hour, Time_minute, Time_second
    if (id_q_ca_out > 0) &
         used = send_data(id_q_ca_out, lpj_state%q_ca_out, lpj_state%current_time &
         !, rmask=lpj_state%lpj_land_mask &
         )
    !if (maxval(lpj_state%runoff_out) == 0.0) __NOTE__('lpj_init_output: runoff_out is 0')
    if (id_runoff_out > 0) &
         used = send_data(id_runoff_out, lpj_state%runoff_out, lpj_state%current_time)
    !if (id_roughness_length_out > 0) &
    !     used = send_data(id_roughness_length_out, lpj_state%roughness_length_out, lpj_state%current_time, &
    !     rmask=lpj_state%lpj_land_mask)
    !if (id_surface_temp_out > 0) &
    !     used = send_data(id_surface_temp_out, lpj_state%surface_temp_out, lpj_state%current_time, &
    !     rmask=lpj_state%lpj_land_mask)
    !if (id_albedo_out > 0) &
    !     used = send_data(id_albedo_out, lpj_state%albedo_out, lpj_state%current_time, &
    !     rmask=lpj_state%lpj_land_mask)
    
    ! now convert to correct units for land_lad
    !lpj_state%co2_flux_out = lpj_state%co2_flux_out / SECONDS_PER_DAY * WTMCO2/WTMC / 1000 ! from [g Carbon/(m^2*day)] to [kg CO2/(m^2*sec)]
    !if (id_co2_flux_out > 0) &
    !     used = send_data(id_co2_flux_out, lpj_state%co2_flux_out, lpj_state%current_time, &
    !     rmask=lpj_state%lpj_land_mask)
    lpj_state%q_ca_out = lpj_state%q_ca_out  ! from daily sum [kg/m^2/day]==[liter/day] to rate [kg/m^2/sec]
    write(*,'(a)') 'lpj_init_output converted to FMS units '
    write(*,'(a,g11.3,a,g11.3,a,g11.3,a,Z20)') '   q_ca_out min ',minval(lpj_state%q_ca_out(:,:)),' max ',maxval(lpj_state%q_ca_out(:,:)),' mean ',sum(lpj_state%q_ca_out(:,:))/size(lpj_state%q_ca_out(:,:)), ' checksum', mpp_chksum(lpj_state%q_ca_out(:,:))
    !lpj_state%runoff_out = lpj_state%runoff_out / SECONDS_PER_DAY ! from daily sum [kg/m^2/day] to rate [kg/m^2/sec]
    ! roughness_length_out
    !lpj_state%surface_temp_out = lpj_state%surface_temp_out + KELVIN ! from [C] to [K]
    ! albedo_out

    !} end init_output
    
    !!! to get the above-written static fields out into the file now,
    !!! and not only after the entire program ended successfully.
    !call diag_manager_end(Time) !!!
    module_is_initialized = .true.

    if (enable_fpe) then
       !! 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)
    endif

    call mpp_clock_end(lpj_init_clock_id)

  end subroutine lpj_driver_init
  ! </SUBROUTINE>


  ! ============================================================================
  ! <SUBROUTINE NAME="lpj_driver_end">
  !   <OVERVIEW>
  !     clean-up for the LPJ driver
  !   </OVERVIEW>
  !   <DESCRIPTION>
  !     ...
  !   </DESCRIPTION>
  !   <PUBLICROUTINE>
  subroutine lpj_driver_end( &
       )
    if (.not. do_lpj) return

    call lpj_end()

    !deallocate(lpj_state%lprec)
    !deallocate(lpj_state%fprec)
    deallocate(lpj_state%prec)
    deallocate(lpj_state%temp_mean)
    deallocate(lpj_state%temp_max)
    deallocate(lpj_state%temp_min)
    deallocate(lpj_state%swdown)
    deallocate(lpj_state%lwnet)
    deallocate(lpj_state%t_flux)
    deallocate(lpj_state%qflux)
    deallocate(lpj_state%dedq)
    deallocate(lpj_state%dhdt)
    deallocate(lpj_state%drdt)
    deallocate(lpj_state%drag_q)
    deallocate(lpj_state%p_surf)
    deallocate(lpj_state%wind)
    deallocate(lpj_state%co2)
    !deallocate(lpj_state%fco2)
    !deallocate(lpj_state%Dfco2Dq)

    deallocate(lpj_state%co2_flux_out)
    deallocate(lpj_state%q_ca_out)
    deallocate(lpj_state%runoff_out)
    deallocate(lpj_state%roughness_length_out)
    deallocate(lpj_state%surface_temp_out)
    deallocate(lpj_state%albedo_out)
    deallocate(lpj_state%evap1_out)    
    deallocate(lpj_state%dedt_out)
    deallocate(lpj_state%gc_out)
    deallocate(lpj_state%exlpj_out)

    deallocate(lpj_state%discharge)

    deallocate(lpj_state%lad_land_mask)
    deallocate(lpj_state%lad_warea)
    deallocate(lpj_state%lpj_land_mask)

    deallocate( &
         lpj_state%coast_src_i, &
         lpj_state%coast_src_j, &
         lpj_state%coast_dst_i, &
         lpj_state%coast_dst_j)

    module_is_initialized = .false.

  end subroutine lpj_driver_end
  ! </SUBROUTINE>  


  subroutine check_land_masks(regridded_mask, native_mask, resname, srcname, dstname)
    real, dimension(:,:), intent(in) :: regridded_mask, native_mask
    character(*), intent(in) :: resname, srcname, dstname

    real, dimension(size(regridded_mask, 1),size(regridded_mask,2)) :: mask_check
    integer :: n

    ! Check matching of regridded interpolation output mask with ``real'' native destination land mask
    mask_check = regridded_mask
    where(native_mask > 0) mask_check = 0.0 ! remaining are cells where regridded src cells fall into dst native ocean cells
    n = count(mask_check > 0)
    if (n > 0) &
         write(stdout(),*) 'WARNING: for ', n, ' cells on ', resname, &
         ' resolution the input from ', srcname, ' vanishes into ocean/antarctica'
       
    ! check if there are native dst cells which do not get input because
    ! the interpolated-from-src grid assumes they blong to the ocean
    mask_check = native_mask
    where (regridded_mask > 0) mask_check = 0.0 ! remaining are dst cells that do not get input
    n = count(mask_check > 0)
    if (n > 0) &
         write(stdout(),'(a,i5,a,a,a,a,a)') 'WARNING: ', n, ' ', dstname, ' cells do not get input from ', srcname, ' because of horizontal interpolation errors'
    
    !! TODO: make this an one-time static field
    !if (id_lpj_mask_out > 0) &
    !     used = send_data(id_lpj_mask_out, regridded_mask, lpj_state%current_time)

  end subroutine check_land_masks




  ! ============================================================================
  ! <SUBROUTINE NAME="lpj_driver_update_fast">
  !   <OVERVIEW>
  !     Accumulates the FMS data of one fast time step as input data for LPJ
  !   </OVERVIEW>
  !   <DESCRIPTION>
  !     CO2 input data is delivered by the coupler on every fast time step.
  !     Thus we must accumulated it here, and save it for lpj_driver_update_slow()
  !
  !     FMS land model does not pass the current time to the update functions,
  !     thus we must use start time and time step information to keep track
  !     of the current time here.
  !   </DESCRIPTION>
  !   <PUBLICROUTINE>
  subroutine lpj_driver_update_fast( &
       ! input variables
       time, & ! from soil
       ! on lad resolution
       q_ca, &                  ! mandatory
       qflux, dedq, &           ! mandatory
       evap1, &
       dedt, &
       gc, &
       exlpj, &
       dhdt, drdt, &            ! mandatory
       drag_q, &                ! mandatory
       p_surf, &                ! mandatory
       lprec, fprec, temp, swdown, lwnet, t_flux, wind, &
       tr_flux_co2, dfdtr_co2 & ! optional
       ! no output variables here
       )
    type(time_type), intent(in) :: time
    ! instant values in FMS units. must be aggregated over each day for LPJ,
    ! and converted into the units expected by LPJ
    real, dimension(:,:), intent(out)           :: & ! on lad resolution
         evap1, &
         dedt, &
         gc, &
         exlpj
    real, dimension(:,:), intent(in)           :: & ! on lad resolution
         q_ca, &            ! [kg/kg] Canopy air specific humidity, needed to convert between CO2 dry ppm and kg CO2/kg wet air
         qflux, &           ! [kg/m2/s] water vapor flux near surface
         dedq, &            ! [kg/m2/s] Moisture flux humidity sensitivity
         dhdt, &            ! derivative of sens over T
         drdt, &            ! derivative of LW radiation over T
         drag_q, &          ! [m/s] Product of drag coefficient for water vapor by wind
         p_surf             ! [N/m2] == [Pa] == [kg/m/s2] Surface pressure
    real, dimension(:,:), intent(in) :: &
         lprec,     & ! [kg/m2/s] liquid precipitation; LPJ needs daily sum [mm] == [kg/m2]
         fprec,     & ! [kg/m2/s] frozen precipitation; LPJ needs daily sum [mm] == [kg/m2]
         temp,      & ! [K] temperature; LPJ needs [C]
         swdown,    & ! [W/m2] shortwave downwelling; LPJ needs same
         lwnet,     & ! [W/m2] longwave net; LPJ needs same
         t_flux,    & ! [W/m2] Sensible heat flux; LPJ needs the same
         wind         ! [m/s] wind speed; LPJ needs same
    real, dimension(:,:), intent(in), optional :: & ! on lad resolution
         tr_flux_co2,     & ! [kg CO2/m2/s]
         dfdtr_co2          ! [kg wet air/m2/s] d(tracer_flux)/d(surface tracer)


    integer :: Time_year, Time_month, Time_day, Time_hour, Time_minute, Time_second, fortytwo = 42
    logical         :: goodmorning, goodnight
    integer :: used
    ! dont refer to optional argument. strange crashes occur if tr_flux_co2 isnot present
    !real, dimension(size(tr_flux_co2,1),size(tr_flux_co2,2)) :: tfco2, dco2
    real, dimension(size(wind,1),size(wind,2)) :: co2_ppmv ! CO2 concentration in ppmv

    if (.not. do_lpj) return
    call mpp_clock_begin(lpj_update_fast_clock_id)

    if (enable_fpe) then
       !! 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)
    endif

    !lpj_state%current_time = lpj_state%current_time + lpj_state%dt_fast
    lpj_state%current_time = time
    call get_date(lpj_state%current_time, Time_year, Time_month, Time_day, Time_hour, Time_minute, Time_second)

    goodmorning = lpj_state%yesterday /= Time_day;
    goodnight = Time_hour*3600+Time_minute*60+Time_second + lpj_state%dt_fast_seconds >= SECONDS_PER_DAY
   ! write(*,*)
   ! write(*,*)
   ! write(*,*)
    !write(*,'(a,i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,l2,a,l2)') &
    !     'lpj_driver_update_fast ',Time_year,'.',Time_month,'.',Time_day, &
    !     ' ',Time_hour,':',Time_minute,':',Time_Second, &
    !     '  goodmorning',goodmorning,' goodnight', goodnight

    if (goodmorning) then
     !   write(*,'(a,i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)') &
     !    'Model time:',Time_year,'.',Time_month,'.',Time_day, &
     !    ' ',Time_hour,':',Time_minute,':',Time_Second
       !write(*,*) '   - goodmorning: zero aggregations'
       !lpj_state%nfaststepstoday = 0
       !lpj_state%lprec = 0.0
       !lpj_state%fprec = 0.0
       lpj_state%prec = 0.0
       !lpj_state%lpj_temp_mean = 288.0-KELVIN ! default initial condition from soil.F90
       !lpj_state%lpj_temp_max = 288.0-KELVIN ! default initial condition from soil.F90
       !lpj_state%lpj_temp_min = 288.0-KELVIN ! default initial condition from soil.F90
       lpj_state%temp_mean = 0.0
       lpj_state%temp_max = -1000.0
       lpj_state%temp_min = 1000.0
       lpj_state%swdown = 0.0
       lpj_state%lwnet = 0.0
       lpj_state%t_flux = 0.0
       lpj_state%qflux = 0.0
       lpj_state%dedq = 0.0
       lpj_state%dhdt = 0.0
       lpj_state%drdt = 0.0
       lpj_state%drag_q = 0.0
       lpj_state%p_surf = 0.0
       lpj_state%wind = 0.0
       ! just continue the integration of flux values over the coming day
       !lpj_state%co2 = 0.0
    end if ! goodmorning

    ! aggregate in FMS-provided units
    !lpj_state%nfaststepstoday = lpj_state%nfaststepstoday + 1

    if (present(tr_flux_co2)) then
       !write(*,*) '   - collecting aggregation for CO2'
       ! This yields kg CO2/kg wet air. LPJ needs [ppm]
       !tfco2 = tr_flux_co2(:,:,1)
       !dco2 = dfdtr_co2(:,:,1)
       !lpj_state%co2 = lpj_state%co2 + tfco2 / (canopy_air_mass_for_tracers/lpj_state%dt_fast_seconds + dco2)
       lpj_state%co2 = lpj_state%co2 &
            + ( lpj_state%co2_flux_out - tr_flux_co2) &
            / (canopy_air_mass_for_tracers/lpj_state%dt_fast_seconds + dfdtr_co2)
       if (id_fco2 > 0) &
            used = send_data(id_fco2, tr_flux_co2, lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
       if (id_Dfco2Dq > 0) &
            used = send_data(id_Dfco2Dq, dfdtr_co2, lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
    end if ! present(tr_flux_co2)

       ! init output variables with default / missing values
       ! TODO: missing values and masks
       lpj_state%co2_flux_out = 0
       lpj_state%q_ca_out = 0
       lpj_state%albedo_out = 0
       if (goodnight) lpj_state%runoff_out = 0 ! recomputed only on goodnight steps
       lpj_state%roughness_length_out = 0
       lpj_state%surface_temp_out = 288.0-KELVIN ! default initial condition from soil.F90 ???
       
    !write(*,*) 'before lpj update'
    ! tracer fluxes (humidity and CO2) are aggregated on fast time step
     call lpj_update(Time_year, Time_month, Time_day, Time_hour, Time_minute, Time_second, &
            ! pass climate status information
            !lpj_state%lpj_fprec, lpj_state%fpj_fprec, &
            (fprec + lprec)* lpj_state%dt_fast_seconds, &
            temp - KELVIN, &
            swdown, lwnet, &
            t_flux, &
            qflux, dedq, &
            dhdt, drdt, &
            drag_q, p_surf, &
            wind, &
            co2_ppmv, &
            ! get back output from LPJ
            lpj_state%co2_flux_out, & ! returned as carbon flux in [g C/ (m^2/day)], must be converted to [kg CO2/(m^2 sec)]
            lpj_state%q_ca_out, &
            lpj_state%runoff_out, & ! recomputed only on goodnight steps, returned in [l/s]
            lpj_state%roughness_length_out, &
            lpj_state%surface_temp_out, & ! returned in [degC]
            lpj_state%albedo_out, &
            lpj_state%evap1_out, &
            lpj_state%dedt_out, &
            lpj_state%gc_out, &
            lpj_state%exlpj_out, &
            !lpj_state%dt_fast_seconds, &
            ! last arg, minimal safety against interface changes
            fortytwo &
            )
       !write(*,*) 'update lpj done'
       !if (goodnight) &
       !write(*,*) 'SURFTEMP_LPJOUT=', maxval(lpj_state%surface_temp_out), 'min=', minval(lpj_state%surface_temp_out), 'mean=', sum(lpj_state%surface_temp_out)/size(lpj_state%surface_temp_out)

       evap1 = lpj_state%evap1_out
       dedt = lpj_state%dedt_out
       gc = lpj_state%gc_out
       exlpj = lpj_state%exlpj_out

       ! diagnostic output
       !
       ! input to LPJ
       !
        if (id_prec > 0) &
            used = send_data(id_prec, lpj_state%prec, lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
       if (id_temp_mean > 0) &
            used = send_data(id_temp_mean, lpj_state%temp_mean, lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
       if (id_temp_min > 0) &
            used = send_data(id_temp_min, lpj_state%temp_min, lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
       if (id_temp_max > 0) &
            used = send_data(id_temp_max, lpj_state%temp_max, lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
       if (id_swdown > 0) &
            used = send_data(id_swdown, lpj_state%swdown, lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
       if (id_lwnet > 0) &
            used = send_data(id_lwnet, lpj_state%lwnet, lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
       if (id_t_flux > 0) &
            used = send_data(id_t_flux, lpj_state%t_flux, lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
       if (id_qflux > 0) &
            used = send_data(id_qflux, lpj_state%qflux, lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
       if (id_dedq > 0) &
            used = send_data(id_dedq, lpj_state%dedq, lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
       if (id_dhdt > 0) &
            used = send_data(id_dhdt, lpj_state%dhdt, lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
       if (id_drdt > 0) &
            used = send_data(id_drdt, lpj_state%drdt, lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
       if (id_wind > 0) &
            used = send_data(id_wind, lpj_state%wind, lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
       ! CO2 is updated on every time step.
       !if (id_co2 > 0) & ! write out CO2 input from coupler in [kg CO2/kg wet air]
       !     used = send_data(id_co2, lpj_state%co2, lpj_state%current_time, &
       !     rmask=lpj_state%lad_land_mask)
       if (id_co2_ppmv > 0) & ! write out CO2 input in ppmv dry air
            used = send_data(id_co2_ppmv, co2_ppmv, lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
       !
       ! LPJ output
       !
       if (id_carbon_flux_out > 0) & ! [g Carbon/(m^2*day)]
            used = send_data(id_carbon_flux_out, lpj_state%co2_flux_out, lpj_state%current_time, &
            rmask=lpj_state%lpj_land_mask)
       !write(*,*) 'sending to id_q_ca_out ', id_q_ca_out, ' at ', Time_year, Time_month, Time_day, Time_hour, Time_minute, Time_second
       if (id_q_ca_out > 0) &
            used = send_data(id_q_ca_out, lpj_state%q_ca_out, lpj_state%current_time &
            !, rmask=lpj_state%lpj_land_mask &
            )
       !if (maxval(lpj_state%runoff_out) == 0.0) __NOTE__('lpj_update (fast): runoff_out is 0')
       if (id_runoff_out > 0) &
            used = send_data(id_runoff_out, lpj_state%runoff_out, lpj_state%current_time)
       if (id_roughness_length_out > 0) &
            used = send_data(id_roughness_length_out, lpj_state%roughness_length_out, lpj_state%current_time, &
            rmask=lpj_state%lpj_land_mask)
       if (id_surface_temp_out > 0) &
            used = send_data(id_surface_temp_out, lpj_state%surface_temp_out, lpj_state%current_time, &
            rmask=lpj_state%lpj_land_mask)
       if (id_albedo_out > 0) &
            used = send_data(id_albedo_out, lpj_state%albedo_out, lpj_state%current_time, &
            rmask=lpj_state%lpj_land_mask)
       if (id_evap1_out > 0) &
            used = send_data(id_evap1_out, lpj_state%evap1_out, lpj_state%current_time, &
            rmask=lpj_state%lpj_land_mask)
       if (id_dedt_out > 0) &
            used = send_data(id_dedt_out, lpj_state%dedt_out, lpj_state%current_time, &
            rmask=lpj_state%lpj_land_mask)
       if (id_gc_out > 0) &
            used = send_data(id_gc_out, lpj_state%gc_out, lpj_state%current_time, &
            rmask=lpj_state%lpj_land_mask)
       if (id_exlpj_out > 0) &
            used = send_data(id_exlpj_out, lpj_state%exlpj_out, lpj_state%current_time, &
            rmask=lpj_state%lpj_land_mask)

       
       ! Convert to land_lad / FMS conforming units

       !lpj_state%co2_flux_out = lpj_state%co2_flux_out / SECONDS_PER_DAY * WTMCO2/WTMC / 1000 ! from [g Carbon/(m^2*day)] to [kg CO2/(m^2*sec)]
       !if (id_co2_flux_out > 0) &
       !     used = send_data(id_co2_flux_out, lpj_state%co2_flux_out, lpj_state%current_time, &
       !     rmask=lpj_state%lpj_land_mask)

       ! runoff is recomputed only at goodnight timesteps.
       if (goodnight) then
          call pour_into_ocean
       end if ! goodnight

       ! roughness_length_out
       lpj_state%surface_temp_out = lpj_state%surface_temp_out + KELVIN ! from [C] to [K]
       ! albedo_out


!    if (.false.) then  !from old version, can be deleted
!
!    ! output CO2 concentration on every fast timestep,
!    ! because it is (potentially) updated on every fast timestep.
!    if (id_co2 > 0) & ! write out CO2 input to LPJ from coupler in [kg CO2/kg wet air]
!         used = send_data(id_co2, lpj_state%co2, lpj_state%current_time, &
!         rmask=lpj_state%lad_land_mask)
!       write(*,*) 'dtslowsteps', lpj_state%dt_slow_seconds
!       write(*,*) '   - goodnight: finish aggregations'
!       write(*,*) 'nslowtimestep=',lpj_state%nslowstepstoday
!
!       write(*,*) 'dtslowsteps', lpj_state%dt_slow_seconds
!       write(*,*) '   - goodnight: finish aggregations'
!       write(*,*) 'nslowtimestep=',lpj_state%nslowstepstoday
!
!! LPJ wants ppm based on volume, not on mass.
!       ! modelled after cold start in land_lad2/canopy_air/canopy_air.F90
!       co2_ppmv = lpj_state%co2 * WTMAIR / (WTMCO2 * (1-q_ca)) * 1e6 ! [kg CO2/kg wet air] -> ppmv dry air
!
!       if (do_check_land_masks) &
!            call check_land_masks(lpj_state%lad_land_mask, lpj_state%lpj_land_mask, 'LPJ', 'land_lad', 'LPJ')
!
!       call lpj_update(...)
!
!       if (do_check_land_masks) &
!            call check_land_masks(lpj_state%lpj_land_mask, lpj_state%lad_land_mask, 'land_lad', 'LPJ', 'land_lad')
!
!
!
!!!! XXX for debugging
!    !call diag_manager_end(time)
!    !stop
!
!
!    else ! not goodnight
!       ! Invoke C++ side to advance time and things, but avoid the overhead for aggregation and interpolation.
!       ! The actual climate info is not accessed for those maintenance works,
!       ! thus just pass NULL.
!       ! Sigh deeply. Passing NULL() is rejected by ifort version 15.
!       ! Workaround: pass the actual arrays; in this invocation they will not be used inside lpj_update() anyway.
!
!       if (.false.) &
!            call lpj_update(...)
!       
!    endif ! goodnight

    lpj_state%yesterday = Time_day

    if (enable_fpe) then
       !! 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)
    endif
    where (lpj_state%q_ca_out>0.2)
       lpj_state%q_ca_out=0.2
    endwhere
    where (lpj_state%q_ca_out<0)
       lpj_state%q_ca_out=0
    endwhere

    call mpp_clock_end(lpj_update_fast_clock_id)

  end subroutine lpj_driver_update_fast
  ! </SUBROUTINE>


  ! ============================================================================
  ! <SUBROUTINE NAME="lpj_driver_update_slow">
  !   <OVERVIEW>
  !     Advances the LPJ model by one (slow!) time step
  !   </OVERVIEW>
  !   <DESCRIPTION>
  !   </DESCRIPTION>
  !   <PUBLICROUTINE>
!  subroutine lpj_driver_update_slow( &
!       ! input variables
!       !time, &
!       ! climate status information on lad resolution
!       !!lprec, fprec, temp, swdown, lwnet, t_flux, wind &
!       ! qflux and dedq are accumulated on fast time step
!       ! , co2 & ! CO2 is accumulated on fast time step
!       ! output variables on lad resolution
!       ! no output variables here, output is written back in the _update_bnd_ subroutines.
!       !co2_flux_out, q_ca_out, runoff_out, roughness_length_out, surface_temp_out &
!       ! TODO: folia protective cover [%] &
!       ! TODO: LAI &
!       )
!
!    !type(time_type), intent(in) :: time
!    ! instant values in FMS units. must be aggregated over each day for LPJ,
!    ! and converted into the units expected by LPJ
!    !!real, dimension(:,:), intent(in) :: &
!    !!     lprec,     & ! [kg/m2/s] liquid precipitation; LPJ needs daily sum [mm] == [kg/m2]
!    !!     fprec,     & ! [kg/m2/s] frozen precipitation; LPJ needs daily sum [mm] == [kg/m2]
!    !!     temp,      & ! [K] temperature; LPJ needs [C]
!    !!     swdown,    & ! [W/m2] shortwave downwelling; LPJ needs same
!    !!     lwnet,     & ! [W/m2] longwave net; LPJ needs same
!    !!     t_flux,    & ! [W/m2] Sensible heat flux; LPJ needs the same
!    !!     wind         ! [m/s] wind speed; LPJ needs same
!
!    if (.not. do_lpj) return
!    call mpp_clock_begin(lpj_update_slow_clock_id)
!
!    if (enable_fpe) then
!       !! 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)
!    endif
!    !write(*,*) 'swdown=', sum(swdown)/size(swdown)
!    write(*,*) 'lpj_driver_update_slow - collecting aggregation'
!    !write(*,'(a,f7.2,a,f7.2,a,f7.2)') '   temp min ',minval(temp),' max ',maxval(temp),' mean ',sum(temp)/size(temp)
!
!    ! Accumulate the climate status information for LPJ over all slow time steps in this day.
!    ! Should be done here and not on the C-side to minimize interpolation overhead.
!    ! aggregate in FMS-provided units
!    ! re-set to 0 is done at goodmorning time steps in lpj_driver_update_fast()
!    lpj_state%nslowstepstoday = lpj_state%nslowstepstoday + 1
!    !lpj_state%lprec = lpj_state%lprec + lprec * lpj_state%dt_slow_seconds
!    !lpj_state%fprec = lpj_state%fprec + fprec * lpj_state%dt_slow_seconds
!
!
!!    lpj_state%prec = lpj_state%prec + (fprec + lprec) ! * lpj_state%dt_slow_seconds
!!    lpj_state%temp_mean = lpj_state%temp_mean + temp
!!    lpj_state%temp_max = max(lpj_state%temp_max, temp)
!!    lpj_state%temp_min = min(lpj_state%temp_min, temp)
!!    lpj_state%swdown = lpj_state%swdown + swdown
!!    lpj_state%lwnet = lpj_state%lwnet + lwnet
!!    lpj_state%t_flux = lpj_state%t_flux + t_flux
!!    lpj_state%wind = lpj_state%wind + wind
!    
!
!!lpj_state%co2 = lpj_state%co2 + co2 ! on fast time step
!
!    if (enable_fpe) then
!       !! 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)
!    endif
!
!    call mpp_clock_end(lpj_update_slow_clock_id)
!
!  end subroutine lpj_driver_update_slow
  ! </SUBROUTINE>  


  ! ============================================================================
  ! <SUBROUTINE NAME="update_lpj_bnd_fast">
  !   <OVERVIEW>
  !     Override land_lad land boundary data data with LPJ data
  !   </OVERVIEW>
  !   <DESCRIPTION>
  !     Updates land boundary data (the ones that atmosphere sees) on the
  !     fast time scale.
  !     The results of the most recent lpj_update() invocation are time-interpolated
  !     onto the land_lad fast time steps and written into the land_data_type structure.
  !   </DESCRIPTION>
  !   <PUBLICROUTINE>
  !   </PUBLICROUTINE>
  subroutine update_lpj_bnd_fast(bnd)
    type(land_data_type), intent(inout) :: bnd        ! state to update
    integer :: used

    !integer :: Time_year, Time_month, Time_day, Time_hour, Time_minute, Time_second
    !call get_date(lpj_state%current_time, Time_year, Time_month, Time_day, Time_hour, Time_minute, Time_second)
    !write(*,'(a,i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)') &
    !     'update_lpj_bnd_fast at', Time_year,'.',Time_month,'.',Time_day, &
    !     ' ',Time_hour,':',Time_minute,':',Time_second

    if (.not. do_lpj) return

    if (enable_fpe) then
       !! 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)
    endif

    ! now time-interpolate the most recent daily output onto the individual time steps
    ! of the following (!!) day
    if (use_co2_from_lpj .and. ico2 /= NO_TRACER) then
       ! put CO2 into tracer
       !!!!! TODO: according to the comments for land_lad2/canopy/cana_tile.F90::cana_prog_type%co2,
       !!!!!       we must return co2 concentration in [kg CO2/kg of wet air] .
       !!!!!       That conforms to the document unit of evaporation as humidity tracer in [kg/kg] .
       bnd%tr(:,:,1,ico2) = lpj_state%co2(:,:)
    end if
    if (id_co2_out_lad > 0 .and. ico2 /= NO_TRACER) then
       used = send_data(id_co2_out_lad, lpj_state%co2(:,:), lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
    end if
    !!
    !! override values which are set in update_vegetation_bnd_fast()
    !!
    !write(*,*) 'computed from land_lad'
    !write(*,'(a,i4,i4,a,i3,i3,i4,i4)') ' size tr(:,:,1,isphum) ', size(bnd%tr(:,:,1,isphum),1), size(bnd%tr(:,:,1,isphum),2), &
    !     ' bounds tr(:,:,1,isphum) ', lbound(bnd%tr(:,:,1,isphum),1), lbound(bnd%tr(:,:,1,isphum),2), &
    !     ubound(bnd%tr(:,:,1,isphum),1), ubound(bnd%tr(:,:,1,isphum),2)
    !write(*,'(a,g11.3,a,g11.3,a,g11.3,a,Z20)') '   tr(:,:,1,isphum) min ',minval(bnd%tr(:,:,1,isphum)),' max ',maxval(bnd%tr(:,:,1,isphum)),' mean ',sum(bnd%tr(:,:,1,isphum))/size(bnd%tr(:,:,1,isphum)), ' checksum', mpp_chksum(bnd%tr(:,:,1,isphum))
    if (use_q_ca_from_lpj) then
       ! put evaporation into humidity tracer
       ! 3rd dimension n_tiles must be fixed to 1
       !write(*,'(a,i4,i4,a,i3,i3,i4,i4)') 'overriding from lpj with size q_ca_out ', size(lpj_state%q_ca_out(:,:),1), size(lpj_state%q_ca_out(:,:),2),' bounds q_ca_out ', lbound(lpj_state%q_ca_out(:,:),1), lbound(lpj_state%q_ca_out(:,:),2), ubound(lpj_state%q_ca_out(:,:),1), ubound(lpj_state%q_ca_out(:,:),2)
       !write(*,'(a,g11.3,a,g11.3,a,g11.3,a,Z20)') '   q_ca_out min ',minval(lpj_state%q_ca_out(:,:)),' max ',maxval(lpj_state%q_ca_out(:,:)),' mean ',sum(lpj_state%q_ca_out(:,:))/size(lpj_state%q_ca_out(:,:)), ' checksum', mpp_chksum(lpj_state%q_ca_out(:,:))

!       where(lpj_state%q_ca_out(:,:)<0)
!         lpj_state%q_ca_out(:,:)=0;
!       endwhere
!       where(lpj_state%q_ca_out(:,:)>0.1)
!         lpj_state%q_ca_out(:,:)=0.1
!       endwhere
       bnd%tr(:,:,1,isphum) = lpj_state%q_ca_out(:,:)

       !write(*,'(a,g11.3,a,g11.3,a,g11.3,a,Z20)') '   tr(:,:,1,isphum) min ',minval(bnd%tr(:,:,1,isphum)),' max ',maxval(bnd%tr(:,:,1,isphum)),' mean ',sum(bnd%tr(:,:,1,isphum))/size(bnd%tr(:,:,1,isphum)), ' checksum', mpp_chksum(bnd%tr(:,:,1,isphum))
    else
       !write(*,*) 'NOT overriding from lpj with'
       !write(*,'(a,g11.3,a,g11.3,a,g11.3,a,Z20)') '   q_ca_out min ',minval(lpj_state%q_ca_out(:,:)),' max ',maxval(lpj_state%q_ca_out(:,:)),' mean ',sum(lpj_state%q_ca_out(:,:))/size(lpj_state%q_ca_out(:,:)), ' checksum', mpp_chksum(lpj_state%q_ca_out(:,:))
    end if
    if (id_q_ca_out_lad > 0) then
       used = send_data(id_q_ca_out_lad, lpj_state%q_ca_out, lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
    end if
    !!
    !! end of override values which are set in update_vegetation_bnd_fast()
    !!


    !!
    !! override values which are set in update_soil_bnd_fast()
    !!
    !write(*,*) 'computed from land_lad'
    !write(*,'(a,i3,i3,a,i3,i3,i3,i3)') ' size t_surf ', size(bnd%t_surf(:,:,1),1), size(bnd%t_surf(:,:,1),2), &
    !     ' bounds t_surf ', lbound(bnd%t_surf(:,:,1),1), lbound(bnd%t_surf(:,:,1),2), &
    !     ubound(bnd%t_surf(:,:,1),1), ubound(bnd%t_surf(:,:,1),2)
    !write(*,'(a,i3,i3,a,i3,i3,i3,i3)')' size t_ca   ', size(bnd%t_ca(:,:,1),1), size(bnd%t_ca(:,:,1),2), &
    !     ' bounds t_ca   ', lbound(bnd%t_ca(:,:,1),1), lbound(bnd%t_ca(:,:,1),2), &
    !     ubound(bnd%t_ca(:,:,1),1), ubound(bnd%t_ca(:,:,1),2)
    !write(*,'(a,f7.2,a,f7.2,a,f7.2,a,Z20)') '   t_surf min ',minval(bnd%t_surf(:,:,1)),' max ',maxval(bnd%t_surf(:,:,1)),' mean ',sum(bnd%t_surf(:,:,1))/size(bnd%t_surf(:,:,1)), ' checksum', mpp_chksum(bnd%t_surf(:,:,1))
    !write(*,'(a,f7.2,a,f7.2,a,f7.2,a,Z20)') '   t_ca   min ',minval(bnd%t_ca(:,:,1)),' max ',maxval(bnd%t_ca(:,:,1)),' mean ',sum(bnd%t_ca(:,:,1))/size(bnd%t_ca(:,:,1)), ' checksum', mpp_chksum(bnd%t_ca(:,:,1))

    if (use_surface_temp_from_lpj) then
       !lpj_state%surface_temp_out(:,:) = bnd%t_surf(:,:,1) ! for testing

       where (lpj_state%surface_temp_out > 400)
          lpj_state%surface_temp_out = 400
       endwhere
       where (lpj_state%surface_temp_out < 190)
          lpj_state%surface_temp_out = 190
       endwhere

       bnd%t_surf(:,:,1) = lpj_state%surface_temp_out(:,:)
       bnd%t_ca(:,:,1) = lpj_state%surface_temp_out(:,:)
      ! write(*,'(a,i3,i3,a,i3,i3,i3,i3)') 'overriding from lpj with size s_t_o ', size(lpj_state%surface_temp_out(:,:)), ' bounds s_t_o ', lbound(lpj_state%surface_temp_out(:,:)), ubound(lpj_state%surface_temp_out(:,:))
      ! write(*,'(a,f7.2,a,f7.2,a,f7.2,a,Z20)') '   temp_out min ',minval(lpj_state%surface_temp_out(:,:)),' max ',maxval(lpj_state%surface_temp_out(:,:)),' mean ',sum(lpj_state%surface_temp_out(:,:))/size(lpj_state%surface_temp_out(:,:)), ' checksum', mpp_chksum(lpj_state%surface_temp_out(:,:))
      ! write(*,'(a,f7.2,a,f7.2,a,f7.2,a,Z20)') '   t_surf min ',minval(bnd%t_surf(:,:,1)),' max ',maxval(bnd%t_surf(:,:,1)),' mean ',sum(bnd%t_surf(:,:,1))/size(bnd%t_surf(:,:,1)), ' checksum', mpp_chksum(bnd%t_surf(:,:,1))
      ! write(*,'(a,f7.2,a,f7.2,a,f7.2,a,Z20)') '   t_ca   min ',minval(bnd%t_ca(:,:,1)),' max ',maxval(bnd%t_ca(:,:,1)),' mean ',sum(bnd%t_ca(:,:,1))/size(bnd%t_ca(:,:,1)), ' checksum', mpp_chksum(bnd%t_ca(:,:,1))
    else
      ! write(*,*) 'NOT overriding from lpj with'
      ! write(*,'(a,f7.2,a,f7.2,a,f7.2,a,Z20)') '   temp_out min ',minval(lpj_state%surface_temp_out(:,:)),' max ',maxval(lpj_state%surface_temp_out(:,:)),' mean ',sum(lpj_state%surface_temp_out(:,:))/size(lpj_state%surface_temp_out(:,:)), ' checksum', mpp_chksum(lpj_state%surface_temp_out(:,:))
    end if

    if (id_surface_temp_out_lad > 0) then
       used = send_data(id_surface_temp_out_lad, lpj_state%surface_temp_out, lpj_state%current_time &
!            , rmask=lpj_state%lad_land_mask &
            )
    end if
    !bnd % albedo     = soil % albedo
    !bnd % albedo_vis_dir = soil%albedo_vis_dir
    !bnd % albedo_nir_dir = soil%albedo_nir_dir
    !bnd % albedo_vis_dif = soil%albedo_vis_dif
    !bnd % albedo_nir_dif = soil%albedo_nir_dif
    !TODO: albedo values from LPJ must be merged/blended with glacier / ice / snow albedos from land_lad
    if (use_albedo_from_lpj) then
       bnd%albedo(:,:,1) = lpj_state%albedo_out(:,:)
       !for now set all to the same value, TODO: caculate seperate values within lpjml
       !it looks like land lad also sets one value for all components
       bnd % albedo_vis_dir(:,:,1) = lpj_state%albedo_out(:,:)
       bnd % albedo_nir_dir(:,:,1) = lpj_state%albedo_out(:,:)
       bnd % albedo_vis_dif(:,:,1) = lpj_state%albedo_out(:,:)
       bnd % albedo_nir_dif(:,:,1) = lpj_state%albedo_out(:,:)

    end if
    if (id_albedo_out_lad > 0) then
       used = send_data(id_albedo_out_lad, lpj_state%albedo_out, lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
       !!<TODO: albedo_X_Y !!
    end if


    !!
    !! end of override values which are set in update_soil_bnd_fast()
    !!

    !! do we really need to pass these to land_lad? what for?
    ! TODO: folia protective cover [%]
    ! TODO: LAI

!!! XXX for debugging
    !call diag_manager_end(lpj_state%current_time)

    if (enable_fpe) then
       !! 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)
    endif

  end subroutine update_lpj_bnd_fast
  ! </SUBROUTINE>


  ! ============================================================================
  ! <SUBROUTINE NAME="pour_into_ocean">
  !   <OVERVIEW>
  !     Pour the water from LPJ runoff into the ocean cells.
  !   </OVERVIEW>
  !   <DESCRIPTION>
  !     Parts of this code was taken from / inspired by land_lad/soil/river.F90
  !
  !     Actually pour the water from LPJ cells into ocean cells,
  !     using the prepared lpj_runoff_coast_map_i/j .
  !     lpj_state%runoff_out was updated at goodnight timestep in lpj_update_fast()
  !     Here, it is assembled into a global map, and then forwarded into discharge,
  !     which then is copied into the corresponding FMS data structure in
  !     subroutine update_lpj_bnd_slow(). There also the division by
  !     per-cell ocean area is done, as in land_lad update_rivers_bnd_slow()
  !     
  !     This subroutine is invoked during init to pour out the runoff that
  !     is read from the LPJ restart file,
  !     and at goodnight timesteps.
  !   </DESCRIPTION>
  subroutine pour_into_ocean
    integer :: i, cellspoured
    real, dimension(1:nlons_global,1:nlats_global) :: global_runoff_out

    !lpj_state%runoff_out(:,:) = lpj_state%runoff_out(:,:) / SECONDS_PER_DAY ! from daily sum [kg/m^2/day] to rate [kg/m^2/sec]
    global_runoff_out(:,:) = 0.0
    !write(*,'(a,4(i4),2(2(i4)))') 'is,ie,js,je,bounds(runoff_out)', &
    !     is, ie, js, je, &
    !     lbound(lpj_state%runoff_out), ubound(lpj_state%runoff_out)
    !write(*,'(a,4(i4),2(2(i4)))') 'is,ie,js,je,bounds(discharge)', &
    !     is, ie, js, je, &
    !     lbound(lpj_state%discharge), ubound(lpj_state%discharge)
    !write(*,'(a,2(2(i4)))') 'bounds(global_runoff_out)', lbound(global_runoff_out), ubound(global_runoff_out)
    !if (maxval(lpj_state%runoff_out) == 0.0) __NOTE__('pour_into_ocean: runoff_out is 0')
    global_runoff_out(is:ie,js:je) = lpj_state%runoff_out(is:ie,js:je)
    !call mpp_sum(global_runoff_out, size(global_runoff_out))
    call mpp_global_field(lpj_state%domain, lpj_state%runoff_out, global_runoff_out)

    ! now pour the water into the oceans
    lpj_state%discharge(:,:) = 0
    cellspoured = 0
    do i=1,lpj_state%ncoastcells
       if (lpj_state%coast_dst_i(i) >= is &
            .and. lpj_state%coast_dst_i(i) <= ie &
            .and. lpj_state%coast_dst_j(i) >= js &
            .and. lpj_state%coast_dst_j(i) <= je) then
          cellspoured = cellspoured + 1
          lpj_state%discharge(lpj_state%coast_dst_i(i),lpj_state%coast_dst_j(i)) &
               = lpj_state%discharge(lpj_state%coast_dst_i(i),lpj_state%coast_dst_j(i)) &
               + global_runoff_out(lpj_state%coast_src_i(i),lpj_state%coast_src_j(i))
          !write(*,'(a,i5,2i4,a,2i4,2g12.5)') &
          !     'pour ',i,lpj_state%coast_src_i(i),lpj_state%coast_src_j(i), &
          !     ' -> ',lpj_state%coast_dst_i(i),lpj_state%coast_dst_j(i), &
          !     global_runoff_out(lpj_state%coast_src_i(i),lpj_state%coast_src_j(i)), &
          !     lpj_state%discharge(lpj_state%coast_dst_i(i),lpj_state%coast_dst_j(i))
       end if
    end do
    !write(*,*) 'poured from ',cellspoured, ' cells'
    !if (maxval(lpj_state%discharge) == 0.0) __NOTE__('pour_into_ocean: discharge is 0')
  end subroutine pour_into_ocean
  ! </SUBROUTINE>


  ! ============================================================================
  ! <SUBROUTINE NAME="update_lpj_bnd_slow">
  !   <OVERVIEW>
  !     Override land_lad land boundary data data with LPJ data
  !   </OVERVIEW>
  !   <DESCRIPTION>
  !     Updates land boundary data (the ones that atmosphere sees) on the
  !     slow time scale.
  !     The results of the most recent lpj_update() invocation are time-interpolated
  !     onto the land_lad slow time steps and written into the land_data_type structure.
  !   </DESCRIPTION>
  !   <PUBLICROUTINE>
  !   </PUBLICROUTINE>
  subroutine update_lpj_bnd_slow(bnd)
    type(land_data_type), intent(inout) :: bnd        ! state to update
    integer :: used

    if (.not. do_lpj) return

    ! now time-interpolate the most recent daily output onto the individual time steps
    ! of the following (!!) day

    if (enable_fpe) then
       !! 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)
    endif

    !!
    !! override values which are set in update_vegetation_bnd_slow()
    !!
    ! no operation
    !!
    !! end of override values which are set in update_vegetation_bnd_slow()
    !!

    !!
    !! override values which are set in update_soil_bnd_slow()
    !!
    if (use_roughness_length_from_lpj) then
       !! TODO: is this the correct variable pairing ?
       !write(*,*) 'roughness_lad=', maxval(bnd%rough_mom(:,:,1)), 'min=', minval(bnd%rough_mom(:,:,1)), 'mean=', sum(bnd%rough_mom(:,:,1))/size(bnd%rough_mom(:,:,1))
       !write(*,*) 'roughness_lad_heat=', maxval(bnd%rough_heat(:,:,1)), 'min=', minval(bnd%rough_heat(:,:,1)), 'mean=', sum(bnd%rough_heat(:,:,1))/size(bnd%rough_heat(:,:,1))
       !write(*,*) 'roughness_lad_scale=', maxval(bnd%rough_scale(:,:,1)), 'min=', minval(bnd%rough_scale(:,:,1)), 'mean=', sum(bnd%rough_scale(:,:,1))/size(bnd%rough_scale(:,:,1))

       bnd%rough_mom(:,:,1) = lpj_state%roughness_length_out(:,:)
       !set the same value for now. TODO. calculate all values within LPJ
       bnd % rough_heat(:,:,1) =  lpj_state%roughness_length_out(:,:) 
       bnd % rough_scale(:,:,1) =  lpj_state%roughness_length_out(:,:)
        !write(*,*) 'roughness=', maxval(bnd%rough_mom(:,:,1)), 'min=', minval(bnd%rough_mom(:,:,1)), 'mean=', sum(bnd%rough_mom(:,:,1))/size(bnd%rough_mom(:,:,1))

    end if
    if (id_roughness_length_out_lad > 0) then
       used = send_data(id_roughness_length_out_lad, lpj_state%roughness_length_out, lpj_state%current_time, &
            rmask=lpj_state%lad_land_mask)
    end if
    !!<TODO: bnd % rough_heat
    !!<TODO: bnd % rough_scale


    !!
    !! override values which are set in update_rivers_bnd_slow()
    !!
    if (use_runoff_from_lpj) then
       ! cut out local compute domain from global map
       !if (maxval(lpj_state%discharge) == 0.0) __NOTE__('update_lpj_bnd_slow: lpj_state%discharge is 0')
       !if (maxval(      bnd%discharge) == 0.0) __NOTE__('update_lpj_bnd_slow: bnd%discharge is 0')
       bnd%discharge(:,:) = lpj_state%discharge(is:ie,js:je) ! [l/s] == [kg/s]
       where ( lpj_state%lad_warea > 0.0 )
          ! adjust partial-ocean-cell discharge values for area fraction
          ! and recalculate to kg/(m2 s)
          bnd % discharge      = bnd % discharge      / lpj_state%lad_warea ! [l/s] -> [l/m^2/s] = [kg/(m2 s)]
          !bnd % discharge_snow = bnd % discharge_snow / rivers%
       elsewhere
          bnd % discharge      = 0
          !bnd % discharge_snow = 0
       endwhere
       if (id_discharge_out_lad_aw > 0) &
            used = send_data(id_discharge_out_lad_aw, bnd%discharge, lpj_state%current_time)
    end if

#if 0
    ! for debugging, pass on in parallel to land_lad runoff
    bnd%lpj_discharge(:,:) = lpj_state%discharge(is:ie,js:je)
    where ( lpj_state%lad_warea > 0.0 )
       ! adjust partial-ocean-cell discharge values for area fraction
       ! and recalculate to kg/(m2 s)
       bnd % lpj_discharge      = bnd % lpj_discharge      / lpj_state%lad_warea
    elsewhere
       bnd % lpj_discharge      = 0
    endwhere
#endif

    ! dont use land mask here, because runoff into ocean is just outside land boundaries
    !if (maxval(lpj_state%runoff_out) == 0.0) __NOTE__('update_lpj_bnd_slow: runoff_out is 0')
    if (id_runoff_out_lad > 0) &
       used = send_data(id_runoff_out_lad, lpj_state%runoff_out, lpj_state%current_time)
    if (id_discharge_out_lad > 0) &
         used = send_data(id_discharge_out_lad, lpj_state%discharge, lpj_state%current_time)
    
    !!<TODO: bnd%discharge_heat
    !!<TODO: bnd%discharge_snow
    !!<TODO: bnd%discharge_snow_heat
    !!<bnd%discharge_heat and bnd%discharge_snow_heat are initialized as 0 and never changed in land_lad

    !!
    !! end of override values which are set in update_rivers_bnd_slow()
    !!

    if (enable_fpe) then
       !! 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)
    endif

  end subroutine update_lpj_bnd_slow
  ! </SUBROUTINE>


  ! ============================================================================
  ! <SUBROUTINE NAME="prepare_coastal_runoff_routing">
  !   <OVERVIEW>
  !     Prepare coastal runoff routing map from LPJmL river routing map
  !   </OVERVIEW>
  !   <DESCRIPTION>
  !     Parts of this code was taken from / inspired by land_lad/soil/river.F90
  !
  !     The input arrays contain coordinates, and the LPJ cell index
  !     of the next cell in the river routing chain.
  !     -1 marks dead ends, i.e. inland sinks, or coastal river mouths.
  !     -9 means no data, i.e. ocean cells. Should not occur in this data.
  !     TODO: check that!!
  !     Follow the method used in the land_lad code to prepare
  !     the coastal runoff routing map:
  !     Iterate over the map.
  !     First look at dead ends of runoff routes, i.e. cells with
  !     value -1. For each such sink cell, search its vicinity for an
  !     ocean cell. If there are ocean cells nearby, store the i,j indices
  !     of the nearest ocean cell as runoff destination in a
  !     global coast cell map.
  !     Further, check for cells where the runoff destination already is
  !     an ocean cell. Copy the values of those cells also into the
  !     coast cell map.
  !     The computed routing information is stored in module-global arrays.
  !   </DESCRIPTION>
  subroutine prepare_coastal_runoff_routing(global_number_lpj_cells, &
         lpj_cell_lons, lpj_cell_lats, lpj_cell_next, &
         lad_gfrac, &
         id_runoff_map, id_runoff_coastcells, &
         id_runoff_dest_lon, id_runoff_dest_lat, &
         id_runoff_dcount, id_runoff_basin, &
         id_runin_coastcells, id_inland_sinks)
    integer, intent(in) :: global_number_lpj_cells
    real, intent(inout), dimension(global_number_lpj_cells) :: lpj_cell_lons, lpj_cell_lats
    integer, intent(inout), dimension(global_number_lpj_cells) :: lpj_cell_next
    real, dimension(:,:), intent(in) :: lad_gfrac ! global fraction of the cell covered by land
    integer, intent(in) :: id_runoff_map, id_runoff_coastcells
    integer, intent(in) :: id_runoff_dest_lon, id_runoff_dest_lat
    integer, intent(in) :: id_runoff_dcount, id_runoff_basin
    integer, intent(in) :: id_runin_coastcells
    integer, intent(in) :: id_inland_sinks

    character(len=*), parameter :: module_name = 'prepare_coastal_runoff_routing'
    integer :: i, j, cell, ii, this_i, this_j, next_i, next_j, used
    real :: dlon, dlat, maxbasin
    ! count statistics of obtained mappings
    integer :: n_10=0, n_9=0, n_1=0, n_other=0, n0=0, nvalid=0
    integer :: ncoastcells = 0, nsinks_already = 0, nsinks_redirected = 0, nsinks_remain = 0

    ! Global runoff routing maps.
    ! For each coast cell, the i and j indices of the target ocean cell are stored,
    ! in Fortran numbering 1..nlons, 1..nlats
    ! Ocean cells are -9, inland sinks are -1.
    integer, dimension(1:nlons_global,1:nlats_global) :: lpj_runoff_coast_map_i
    integer, dimension(1:nlons_global,1:nlats_global) :: lpj_runoff_coast_map_j

    ! for static diagnostic output (dcount, basin)
    real, dimension(nlons_global, nlats_global) :: map, diag
    !real, dimension(is:ie,js:je) :: inland_sinks

    ! path to find redirection points in the vicinity
    integer, parameter :: n_didj = 81
    integer, parameter :: di(n_didj) = &
         (/ -4, 4,-4, 4,-4, 4,-3, 3, 4,-4, 3,-3,-2, 2,-4, 4,-4, 4,-2, 2, &
         3,-3, 3,-3,-1, 1, 4,-4,-4, 4,-1, 1, 0, 0,-4, 4,                 &
         2,-3, 3,-3, 3,-3, 2,-2,                                         &
         1,-1, 3,-3,-3, 3, 1,-1, 3,-3, 0, 0,                             &
         -2, 2,-2, 2,-1, 1, 2,-2,-2, 2,-1, 1, 0, 0,-2, 2,                &
         1,-1, 1,-1, 1,-1, 0, 0, 0/)

    integer, parameter :: dj(n_didj) = &
         (/  4,-4,-4, 4, 3,-3,-4, 4, 3,-3,-4, 4, 4,-4,-2, 2,-2,-2,-4, 4, &
         3,-3,-3, 3,-4, 4,-1, 1,-1, 1, 4,-4,-4, 4, 0, 0,                 &
         3,-2,-2, 2, 2,-2,-3, 3,                                         &
         -3, 3, 1,-1, 1,-1, 3,-3, 0, 0, 3,-3,                            &
         2,-2,-2, 2,-2, 2,-1, 1,-1, 1, 2,-2,-2, 2, 0, 0,                 &
         1,-1,-1, 1, 0, 0, 1,-1, 0/)

    character(len=nlons_global) :: buffer
    logical :: changed
    integer :: riverlength
    integer :: riveradd  ! how much dcount to add to new cells in the path
    integer :: riverincr ! increment from cell cell, 1 or 0

    __NOTE__("starting")

    ! first print map before the calls to mpp_sum. 
    ! LPJ counts cells from 0, thus subtract 1 in the printout
    !write(*,*) 'discharge map as obtained from LPJ, only local domain filled in:'
    !do cell = 1, global_number_lpj_cells
    !   write(*,'(a,i5,f8.2,f8.2,a,i5)') 'discharge map A ',cell-1,lpj_cell_lons(cell),lpj_cell_lats(cell),' -> ',lpj_cell_next(cell)
    !end do

    dlon = 360.0/nlons_global
    dlat = 180.0/nlats_global
    write(*,'(a,f5.2,a,f5.2)') 'dlon ',dlon,' dlat ',dlat

    ! we (ab)use global sum. This is safe, because each task
    ! contributes non-zero values only for its local cells.
    ! This was verified using the above and below write statements.
    ! It works, despite strange warning from compiler and from
    ! valgrind about using uninitialized values. Sigh.
    call mpp_sum(lpj_cell_lons, global_number_lpj_cells)
    call mpp_sum(lpj_cell_lats, global_number_lpj_cells)
    call mpp_sum(lpj_cell_next, global_number_lpj_cells)

    ! print map resulting from mpp_sum.
    !write(*,*) 'discharge map after mpp_sum, now global domain should be filled in:'
    !do cell = 1, global_number_lpj_cells
    !   write(*,'(a,i5,f8.2,f8.2,a,i5)') 'discharge map B ',cell-1,lpj_cell_lons(cell),lpj_cell_lats(cell),' -> ',lpj_cell_next(cell)
    !end do

    lpj_runoff_coast_map_i = -10
    lpj_runoff_coast_map_j = -10

    write(*,*) 'convert the LPJ next cell indices to FMS global grid indices i,j'
    do cell = 1, global_number_lpj_cells
       ! LPJ counts cells from 0, thus subtract 1 in the printout
       this_i = lon2i(lpj_cell_lons(cell), dlon)
       this_j = lat2j(lpj_cell_lats(cell), dlat)
       if (is_coastal(this_i,this_j,lad_gfrac)) ncoastcells = ncoastcells+1
       !write(*,'(a,i5,f8.2,f8.2,i5,i5,a,i5$)') ' discharge map ',cell-1,lpj_cell_lons(cell),lpj_cell_lats(cell),this_i,this_j,' -> ',lpj_cell_next(cell)
       __CHECK__(this_i >= 1 .and. this_i <= nlons_global, "this_i out of range")
       __CHECK__(this_j >= 1 .and. this_j <= nlats_global, "this_j out of range")
       if (lpj_cell_next(cell) >= 0) then
          __CHECK__(lpj_cell_next(cell) < global_number_lpj_cells, "next out of range")
          nvalid = nvalid+1
          next_i = lon2i(lpj_cell_lons(lpj_cell_next(cell)+1), dlon)
          next_j = lat2j(lpj_cell_lats(lpj_cell_next(cell)+1), dlat)
          !write(*,'(a,f8.2,f8.2,i5,i5$)') ' ',lpj_cell_lons(lpj_cell_next(cell)+1),lpj_cell_lats(lpj_cell_next(cell)+1),next_i,next_j
          __CHECK__(next_i >= 1 .and. next_i <= nlons_global, "next_i out of range")
          __CHECK__(next_j >= 1 .and. next_j <= nlats_global, "next_j out of range")
          !if (this_i == next_i .and. this_j == next_j) write(*,'(a$)') ' loop to self'
       else if (lpj_cell_next(cell) == 0) then
          n0 = n0+1
          next_i = 0 ! ?
          next_j = 0 ! ?
       else if (lpj_cell_next(cell) == -1) then ! dead end
          next_i = -1
          next_j = -1
          n_1 = n_1+1
       else if (lpj_cell_next(cell) == -9) then ! ocean cell in land?
          n_9 = n_9+1
          next_i = lpj_cell_next(cell)
          next_j = lpj_cell_next(cell)
       else ! some other negative value
          n_other = n_other+1
          next_i = lpj_cell_next(cell)
          next_j = lpj_cell_next(cell)
          !write(*,'(a,i7,a,i7,a,f7.2,a,f7.2)') &
          !     'Warning: '//module_name//' found river routing next == ', &
          !     lpj_cell_next(cell), &
          !     ' in cell ', cell, ' at ', &
          !     lpj_cell_lons(cell),',',lpj_cell_lats(cell)
       end if
       !write(*,*) ' ' ! just end the line
       lpj_runoff_coast_map_i(this_i, this_j) = next_i
       lpj_runoff_coast_map_j(this_i, this_j) = next_j
    end do

    write(*,*) 'got from LPJ -10/-9/-1/-X/0/valid/other', &
         n_10, n_9, n_1, n_other, n0, nvalid, &
         global_number_lpj_cells-(n_10+n_9+n_1+n_other+n0+nvalid)
    write(*,*) 'river route map has ', ncoastcells, ' coast cells'
         
    ! Now the remaining -9 entries in lpj_runoff_coast_map_i (and _j)
    ! should match the ocean part of the lad_land_mask. Check that.
    ! (this could be done in above loop already,
    ! but the code would look uglier)
    do j=1,nlats_global
       do i=1,nlons_global
          if (lad_gfrac(i,j) > 0) then
             if (-9 == lpj_runoff_coast_map_i(i,j)) then
                !write(*,'(a,i5,a,i5)') &
                !     'Warning: '//module_name &
                !     //' found land cell in lad_gfrac' &
                !     //' but LPJ ocean cell in river routing next' &
                !     //' at index ', i,',',j
             end if
             if (-10 == lpj_runoff_coast_map_i(i,j)) then
                !write(*,'(a,i5,a,i5)') &
                !     'Warning: '//module_name &
                !     //' found land cell in lad_gfrac' &
                !     //' but lad ocean cell in river routing next' &
                !     //' at index ', i,',',j
             end if
          else
             if (-9 /= lpj_runoff_coast_map_i(i,j)) then
                !write(*,'(a,i5,a,i5)') &
                !     'Warning: '//module_name &
                !     //' found LPJ ocean cell in lad_gfrac,' &
                !     //' but land cell in river routing next' &
                !     //' at index ', i,',',j
             end if
             if (-10 /= lpj_runoff_coast_map_i(i,j)) then
                !write(*,'(a,i5,a,i5)') &
                !     'Warning: '//module_name &
                !     //' found lad ocean cell in lad_gfrac,' &
                !     //' but land cell in river routing next' &
                !     //' at index ', i,',',j
             end if
          end if
       end do
    end do

    !write(*,*) 'print a nice map'
    !do j=nlats_global,1,-1 ! print N to S
    !   do i=1,nlons_global
    !      select case (lpj_runoff_coast_map_i(i,j))
    !      case (-10)
    !         buffer(i:i) = '.'
    !      case (-9)
    !         buffer(i:i) = '+'
    !      case (-1)
    !         buffer(i:i) = 'X'
    !      case (0)
    !         buffer(i:i) = '0'
    !      case default
    !         if (lpj_runoff_coast_map_i(i,j) == i .and. &
    !              lpj_runoff_coast_map_j(i,j) == j) then
    !            buffer(i:i) = '@' ! route into self
    !         else
    !            buffer(i:i) = '*'
    !         end if
    !      end select
    !   end do
    !   write(*,'(a)') buffer
    !end do

    __NOTE__("Now try to re-route coastal sink cells")

    ! In the first pass, just count the number of successes,
    ! so that we know how much space to allocate for the final map.
    ! TODO: Optimization: in each task, store and process only those parts
    ! of the coast runoff map, that pours into the local compute domain.
    nsinks_redirected = 0
    nsinks_already = 0
    do j=1,nlats_global
       do i=1,nlons_global
          next_i = lpj_runoff_coast_map_i(i,j)
          next_j = lpj_runoff_coast_map_j(i,j)
          if (next_i > 0) then
             if (is_coastal(next_i, next_j, lad_gfrac)) then
                ! cell i,j already routes into coast cell.
                ! Just count it.
                nsinks_already = nsinks_already + 1 
                cycle
             end if
          end if
          if (-1 == lpj_runoff_coast_map_i(i,j)) then
             ! try to find a nearby ocean cell to redirect runoff to.
             ! to avoid branching out of loop, currently search all
             ! possibilities, over-writing from most distant to least distant.
             ! backward search direction, from closest to farthest,
             ! and break out of loop.
             do ii = size(di(:)), 1, -1
                next_i = modulo(i+di(ii)-1, nlons_global) + 1
                next_j = max(min(j+dj(ii), nlats_global), 1)
                if (is_coastal(next_i, next_j, lad_gfrac)) then
                   ! found a new destination that is OK
                   nsinks_redirected = nsinks_redirected + 1
                   exit ! from loop over ii
                end if
             end do
          end if
       end do
    end do
    nsinks_remain = n_1 - nsinks_redirected
    write(*,*) nsinks_already, ' cells already routing into ocean cell'
    write(*,*) 'redirected routing for ',nsinks_redirected,' coast cells.'
    write(*,*) 'remaining inland sink cells: ',nsinks_remain

    ! Prepare a compressed coastal runoff map in the shape of 4 1-D arrays,
    ! on the assumption that this is cheaper to handle than full global maps.
    !lpj_state%ncoastcells = nsinks_already+nsinks_redirected
    lpj_state%ncoastcells = nsinks_redirected+nsinks_already
    allocate( &
         lpj_state%coast_src_i(nsinks_redirected+nsinks_already), &
         lpj_state%coast_src_j(nsinks_redirected+nsinks_already), &
         lpj_state%coast_dst_i(nsinks_redirected+nsinks_already), &
         lpj_state%coast_dst_j(nsinks_redirected+nsinks_already))

    __NOTE__("2nd pass")

    ! 2nd pass
    map = 0
    diag = -1
    allocate(lpj_state%inland_sinks(is:ie,js:je))
    lpj_state%inland_sinks = 0
    cell = 0
    do j=1,nlats_global
       do i=1,nlons_global
          next_i = lpj_runoff_coast_map_i(i,j)
          next_j = lpj_runoff_coast_map_j(i,j)
          if (next_i > 0) then
             map(i, j) = (next_j-1)*nlons+next_i+1
             if (is_coastal(next_i, next_j, lad_gfrac)) then
                !write(*,'(a,i5,i5,a,i5,i5,a,f6.4)') &
                !     'cell ',i,j,' already routes into coast cell ',next_i, next_j,' gfrac ',lad_gfrac(i,j)
                cell = cell + 1
                lpj_state%coast_src_i(cell) = i
                lpj_state%coast_src_j(cell) = j
                lpj_state%coast_dst_i(cell) = next_i
                lpj_state%coast_dst_j(cell) = next_j
                diag(i, j) = 0
             end if
          else
             map(i, j) = next_i ! copy missing value
          end if
          if (-1 == next_i) then
             ! try to find a nearby ocean cell to redirect runoff to.
             ! to avoid branching out of loop, currently search all
             ! possibilities, over-writing from most distant to least distant.
             ! backward search direction, from closest to farthest,
             ! and break out of loop.
             do ii = size(di(:)), 1, -1
                next_i = modulo(i+di(ii)-1, nlons_global) + 1
                next_j = max(min(j+dj(ii), nlats_global), 1)
                if (is_coastal(next_i, next_j, lad_gfrac)) then
                   ! found a new destination that is OK
                   lpj_runoff_coast_map_i(i,j) = next_i
                   lpj_runoff_coast_map_j(i,j) = next_j
                   !write(*,'(a,i5,a,i5,a,i5,a,i5,a,f6.4$)') 'redirect coastal runoff from ', &
                   !     i,',',j,' to ', next_i,',',next_j,' gfrac ',lad_gfrac(i,j)
                   !if (i == next_i .and. j == next_j) write(*,'(a$)') ' redirected loop to self'
                   !write(*,*) ' '
                   cell = cell + 1
                   lpj_state%coast_src_i(cell) = i
                   lpj_state%coast_src_j(cell) = j
                   lpj_state%coast_dst_i(cell) = next_i
                   lpj_state%coast_dst_j(cell) = next_j
                   diag(i, j) = ii
                   map(i, j) = (next_j-1)*nlons+next_i+1
                   exit ! from loop over ii
                end if
             end do ! loop over ii
             if (0 == ii) then
                write(*,'(a,i5,a,i5)') 'no coastal redirection found from ',i,',',j
                if (i>=is.and.i<=ie.and.j>=js.and.j<=je) lpj_state%inland_sinks(i,j) = 1 ! prepare diagnostic output of inland sink locations
             end if
          end if
       end do
    end do
    if (id_runoff_map > 0) &
         used = send_data(id_runoff_map, map(is:ie,js:je))
    if (id_runoff_coastcells > 0) &
         used = send_data(id_runoff_coastcells, diag(is:ie,js:je))
    if (id_inland_sinks > 0) &
         used = send_data(id_inland_sinks, lpj_state%inland_sinks(is:ie,js:je))
    __CHECK__(cell == nsinks_redirected+nsinks_already, "2nd redirection pass got inconsistent result")

    if (id_runoff_dest_lon > 0) then
       diag = -1000.0
       do j=js,je
          do i=is,ie
             if (lpj_runoff_coast_map_i(i,j) > 0) &
                  diag(i,j) = i2lon(lpj_runoff_coast_map_i(i,j), dlon)
          end do
       end do
       used = send_data(id_runoff_dest_lon, diag(is:ie,js:je))
    end if

    if (id_runoff_dest_lat > 0) then
       diag = -1000.0
       do j=js,je
          do i=is,ie
             if (lpj_runoff_coast_map_j(i,j) > 0) &
                  diag(i,j) = j2lat(lpj_runoff_coast_map_j(i,j), dlat)
          end do
       end do
       used = send_data(id_runoff_dest_lat, diag(is:ie,js:je))
    end if
    

    ! static diagnostics:
    if (id_runoff_dcount > 0 .or. id_runin_coastcells > 0) then
       __NOTE__("dcount diagnostics")
       !TODO: use array map to improve cycle detection by coloring it with
       ! the number of the currently followed river flow. If we encounter a map cell
       ! which already is colored with the current number, we have a cycle.
       ! That saves counting up to the limit.
       diag = -10 ! ocean cells
       map = 0
       ii = 0 ! use here for river number
       do j=1,nlats_global
          do i=1,nlons_global
             riverlength = 0
             if (diag(i,j) /= -10) &
                  cycle ! We have been here before. No sense to add 0 all along downstream path
             if (lpj_runoff_coast_map_i(i,j) <= 0) cycle
             ! follow down-river
             !write(*,'(a,i5,i5,a,i5,i5)') 'follow from ',i,j,' to ',lpj_runoff_coast_map_i(i,j), lpj_runoff_coast_map_j(i,j)
             this_i = i
             this_j = j
             ii = ii+1
             riveradd = 1 ! how much to add to new cells in the path
             riverincr = 1 ! increment from cell to cell, might be set to 0 below
             riverlength = 1
             !write(*,'(a,i5,i5$)') ' start ',i,j
             do while (.true.) ! sigh. F90 does not have && shortcut evaluation like C
                next_i = lpj_runoff_coast_map_i(this_i,this_j)
                next_j = lpj_runoff_coast_map_j(this_i,this_j)
                !write(*,'(a,i5,i5$)') ' -> ',next_i, next_j
                if (next_i <= 0 .or. next_j <= 0) then ! sink
                   !write(*,*) ' sink nodata len ',riverlength
                   exit
                end if
                if (next_i == this_i .and. next_j == this_j) then ! sink to self
                   !write(*,*) ' sink self len ',riverlength
                   exit
                end if
                if (map(this_i,this_j) == ii & ! routing loop, or bug
                     .or. riverlength == 1000 &
                     .or. diag(next_i, next_j) >= global_number_lpj_cells &
                     ) then
                   !write(*,*) ' loop len ',riverlength
                   !write(*,'(a,i5,a,i5,a,f7.2,f7.2,a,i7)') 'Warning: river routing loop at ',next_i,',',next_j,' = ',i2lon(next_i,dlon),j2lat(next_j,dlat),' len ',riverlength
                   exit
                end if
                map(this_i,this_j) = ii ! color it
                !if (diag(next_i, next_j) >= global_number_lpj_cells) exit ! loop or bug
                if (diag(next_i, next_j) == -10) then
                   diag(next_i, next_j) = riveradd
                else
                   diag(next_i, next_j) = diag(next_i, next_j) + riveradd
                   ! We have been here before. Stop increasing riveradd,
                   ! in order to not count downstream cells multiple times.
                   riverincr = 0
                end if
                this_i = next_i ! and move on
                this_j = next_j
                riveradd = riveradd + riverincr
                riverlength = riverlength + 1
             end do
          end do
       end do
       if (id_runoff_dcount > 0) &
            used = send_data(id_runoff_dcount, diag(is:ie,js:je))
       !write(*,*) 'print a nice dcount map'
       !do j=nlats_global,1,-1 ! print N to S
       !   write(*,'(*(i6))') int(diag(1:nlons_global,j))
       !end do

       if (id_runin_coastcells > 0) then
          __NOTE__("beach cell diagnostics")
          ! beach cells := ocean cells with number of river cells routing into them
          map = -1.0
          do cell=1,nsinks_redirected
             this_i = lpj_state%coast_src_i(cell)
             this_j = lpj_state%coast_src_j(cell)
             next_i = lpj_state%coast_dst_i(cell)
             next_j = lpj_state%coast_dst_j(cell)
             if (map(next_i, next_j) < 0) map(next_i, next_j) = 0
             map(next_i, next_j) = map(next_i, next_j) + diag(this_i, this_j)
          end do
          used = send_data(id_runin_coastcells, map(is:ie,js:je))
       end if ! id_runin_coastcells > 0
    end if ! id_runoff_dcount > 0 .or. id_runin_coastcells > 0


    if (id_runoff_basin > 0) then
       __NOTE__("basin diagnostics")
    
       ! basins of size 0, i.e. cells without outflow and without inflow get 0.
       ! All non-trivial basins get a unique number.
       ! Ocean cells get -10
       diag = -10
       maxbasin = 1
       changed = .true. ! Termination flag. If anything was changed, we need another loop over the cells.
       do while (changed)
          changed = .false.
          do j=1,nlats_global
             do i=1,nlons_global
                !write(*,*) 'looking at ',i,j, lpj_runoff_coast_map_i(i,j), lpj_runoff_coast_map_j(i,j)
                if (lpj_runoff_coast_map_i(i,j) > 0) then
                   this_i = i
                   this_j = j
                   next_i = lpj_runoff_coast_map_i(this_i,this_j)
                   next_j = lpj_runoff_coast_map_j(this_i,this_j)
                   if (diag(this_i, this_j) == -10 & ! new cell
                        .and. (next_i <= 0 & ! no outflow
                        .or. (next_i == this_i .and. next_j == this_j))) then ! outflow to self
                      !write(*,*) '  new cell without outflow -> basin size 0'
                      diag(this_i, this_j) = 0 ! might be a basin of size 0
                      changed = .true.
                      cycle
                   end if
                   if (diag(this_i,this_j) == -10) then ! found a new cell with outflow
                      !write(*,*) '  new cell with outflow -> basin',maxbasin
                      diag(this_i,this_j) = maxbasin
                      maxbasin = maxbasin+1
                      changed = .true.
                   end if
                   ! follow down-river
                   riverlength = 1
                   !write(*,'(a,i5,i5,a,i5,i5)') 'follow from ',i,j,' to ',lpj_runoff_coast_map_i(i,j), lpj_runoff_coast_map_j(i,j)
                   do while (.true.) ! sigh. F90 does not have && shortcut evaluation like C
                      if (next_i <= 0 .or. next_j <= 0) exit
                      if (next_i == this_i .and. next_j == this_j) exit
                      if (riverlength >= global_number_lpj_cells) exit
                      if (diag(next_i,next_j) == -10 & ! new cell
                           .or. diag(next_i,next_j) == 0) then ! cell without outflow
                         diag(next_i,next_j) = diag(this_i,this_j) ! now we found an inflow
                         !write(*,*) '  found inflow to basin 0 from', diag(this_i,this_j)
                         changed = .true.
                      else if (diag(next_i,next_j) == diag(this_i,this_j)) then
                         !write(*,*) ' stop follow downriver'
                         exit ! apparently we have been along this path before, no more changes to do downstream
                      else
                         ! Found inflow into a previously found different basin.
                         ! Propagate the lower basin number.
                         diag(next_i,next_j) = min(diag(next_i,next_j), diag(this_i,this_j))
                         diag(this_i,this_j) = diag(next_i,next_j)
                         !write(*,*) ' connecting basins',diag(next_i,next_j), diag(this_i,this_j)
                         changed = .true.
                      end if
                      riverlength = riverlength + 1
                      !if (riverlength == global_number_lpj_cells) &
                      !     write(*,*) 'Warning: river routing loop at ',next_i,',',next_j
                      this_i = next_i ! and move on
                      this_j = next_j
                      next_i = lpj_runoff_coast_map_i(this_i,this_j)
                      next_j = lpj_runoff_coast_map_j(this_i,this_j)
                   end do
                end if
             end do
          end do
       end do

       write(*,*) 'maxbasin ', maxbasin
       used = send_data(id_runoff_basin, diag(is:ie,js:je))
    end if ! id_runoff_basin > 0

  end subroutine prepare_coastal_runoff_routing
  ! </SUBROUTINE>


  ! <FUNCTION NAME="lon2i">
  ! <OVERVIEW>Convert from longitude to grid index</OVERVIEW>
  integer function lon2i(lon_in, dlon)
    real, intent(in) :: lon_in, dlon
    real :: lon
    lon = lon_in
    if (lon_in < 0.0) lon = lon + 360.0 ! transform -180 .. 180 to 0 .. 360
    lon2i = lon / dlon + 1
  end function lon2i
  ! </FUNCTION>


  ! <FUNCTION NAME="lat2j">
  ! <OVERVIEW>Convert from latitude to grid index</OVERVIEW>
  integer function lat2j(lat, dlat)
    real, intent(in) :: lat, dlat
    lat2j = (lat+90.0) / dlat + 1
  end function lat2j
  ! </FUNCTION>

  ! <FUNCTION NAME="i2lon">
  ! <OVERVIEW>Convert from grid index to longitude</OVERVIEW>
  real function i2lon(i_in, dlon)
    integer, intent(in) :: i_in
    real, intent(in) :: dlon
    i2lon = dlon/2.0 + (i_in-1) * dlon
  end function i2lon
  ! </FUNCTION>

  ! <FUNCTION NAME="j2lat">
  ! <OVERVIEW>Convert from grid index to longitude</OVERVIEW>
  real function j2lat(j_in, dlat)
    integer, intent(in) :: j_in
    real, intent(in) :: dlat
    j2lat = -90.0 + dlat/2.0 + (j_in-1) * dlat
  end function j2lat
  ! </FUNCTION>



  ! <FUNCTION NAME="is_coastal">
  !   <OVERVIEW>
  !     Tests whether a point is a coastal point.
  !   </OVERVIEW>
  !
  !   <DESCRIPTION>
  !     Returns true if and only if point (i,j) is a coastal point on the map gfrac.
  !   </DESCRIPTION>
  !
  !   <TEMPLATE>
  !     value=is_coastal(i, j, gfrac)
  !   </TEMPLATE>
  logical function is_coastal(i, j, gfrac)
    integer, intent(in) :: i, j       ! coordinates of the point in question
    real,    intent(in) :: gfrac(:,:) ! fractional area of the land

    integer ip,in, jp,jn ! coordinates of surrounding points

    ! find coordinates of surrounding points
    ip = i-1; if (ip<1) ip = size(gfrac, 1)
    in = i+1; if (in>size(gfrac,1)) in = 1
    jp = max(1, j-1); 
    jn = min(size(gfrac,2), j+1)

    !write(*,*) 'is_coastal bounds(gfrac) ',lbound(gfrac), ubound(gfrac)
    !write(*,'(a,i4,i4,i4,i4,i4,i4)') '  i,j,ip,jp,in,jn ',i,j,ip,jp,in,jn

    __ASSERT__(i >= lbound(gfrac,1) .and. i <= ubound(gfrac,1), "i out of range")
    __ASSERT__(j >= lbound(gfrac,2) .and. j <= ubound(gfrac,2), "j out of range")
    __ASSERT__(ip >= lbound(gfrac,1) .and. ip <= ubound(gfrac,1), "ip out of range")
    __ASSERT__(jp >= lbound(gfrac,2) .and. jp <= ubound(gfrac,2), "jp out of range")
    __ASSERT__(in >= lbound(gfrac,1) .and. in <= ubound(gfrac,1), "in out of range")
    __ASSERT__(jn >= lbound(gfrac,2) .and. jn <= ubound(gfrac,2), "jn out of range")
    
    is_coastal = gfrac(i,j) < (1.0-min_water_frac) .and.( &
         gfrac(i ,j ) > min_land_frac .or. &
         gfrac(in,j ) > min_land_frac .or. &
         gfrac(ip,j ) > min_land_frac .or. &
         gfrac(i ,jn) > min_land_frac .or. &
         gfrac(i ,jp) > min_land_frac      )
  end function is_coastal
  ! </FUNCTION>



! ==============================================================================
! Reports error, including file name and line.
subroutine my_error(mod_name, message, mode, file, line)

  character(len=*), intent(in) :: mod_name
  character(len=*), intent(in) :: message
  integer,          intent(in) :: mode
  character(len=*), intent(in), optional :: file
  integer,          intent(in), optional :: line

  ! ---- local vars ----------------------------------------------------------
  character(len=512) :: mesg
  if(present(file)) then ! assume that file and line are either both present or not
     write(mesg,'("File ",a," Line ",i4.4," :: ",a)')&
       file, line, trim(message)
  else
    mesg = trim(message)
  endif
  call error_mesg(mod_name, mesg, mode)
end subroutine

end module lpj_driver_mod
#endif
