/**
 * \file Aeolus.cpp
 *
 * $Id$
 *
 * This is the C++ side of the the interface between
 * MOM4/FMS and the Potsdam-3 atmosphere model
 *
 * The official name of the stand-alone atmosphere model is now Aeolus.
 * It started its life under the nick-name potsdam3,
 * thus some functions still carry that nick-name.
 */
#if defined(USE_MPI) || defined(USE_MPI_DOMAIN_DECOMPOSE)
// Intel MPI requires mpi.h before stdio.h
#include <mpi.h>
#else
typedef int MPI_Comm;
#endif

#include <cstdlib>
#include <cmath>
#include <cstdio>
#include <assert.h>
#include <cstring>
#include <unistd.h>
#include <mcheck.h>

#include "SphericalGrid.h"
//#include "Variable.h"
//#include "CellVariable.h"
#include "HorizontalVariable.h"
#include "ScalarVariable.h"
#include "NetCDFOutput.h"
#include "Restart.h"
#include "NetCDFInput.h"
#include "AuxiliaryFunctionality.h"
#include "PlanWave.h"
#include "SynopticScaleDynamics.h"
#include "HeatTransfer.h"
#include "LargeScaleDynamics.h"
#include "PlanetaryBoundaryLayerPhysics.h"
#include "Cloudiness_2.h"
#include "Cloudiness_3Layers.h"
#include "RadiativeTransfer.h"
#include "HeatTransfer.h"
#include "WaterVaporTransfer.h"
#include "Atmosphere.h"
#include "FMSAstronomy.h"

#ifdef WITH_SIMENV
#include "simenv_mod_c.inc"
#include "SimenvOutput.h"
int simenv_run_int;
static SimenvOutput *simenv_output = NULL;
#endif

#ifdef _OPENMP
#include "omp.h"
#endif

// Define the following symbol to enable full tracing of
// all variables that are exchanged between FMS and Aeolus.
// The same facility is already in atmos_coupled/atmos_model.F90 ,
// using the FMS diag_manager . However, that diag_manager buffers
// all output before writing, thus the most interesting data is lost
// if the program crashes.
// In contrast, Aeolus NetCDFOutput writes out everything immediately.
// To avoid overhead in time, memory, and disk space, we make this
// facility a compile-time conditional.
//#define TRACE_AEOLUS_UPDATE

// Define the following symbol to invoke stock taking inbetween the
// single compute steps inside the update functions.
// This is expensive and only useful for debugging (if at all),
// thus we make it a compile-time conditional.
//#define DEBUG_AEOLUS_STOCKS


//#ifndef __FUNCTION__
// // __FUNCTION__ seems to be GNU specific.
// // __func__ is supposed to be ISO C99 conformant
//#define __FUNCTION__ __func__
//#endif


static int is, ie, js, je;    // boundary indices of local compute domain in Fortran notation 1..n
static MPI_Comm comm;         // MPI communicator handle for the atmosphere domain
static int pe, root_pe, npes; // our own MPI rank, root MPI rank, total number of atmosphere ranks

static int nr_tracers;  // total number of tracers in surf_diff_type fields
static int q_ind;       // index of humidity in surf_diff_type tracer fields
static int co2_ind;     // index of CO2 in surf_diff_type tracer fiels, or -1 for none

static double dt_atm_sec;   // timestep of atmosphere in seconds
#ifndef NDEBUG
static int first_day;       // day of first timestep of this run, for human-readable timestamps
#endif

//static int debug_out_counter;        // counter for DEBUG output
//static int count_out;                // counter for normal Aeolus output
static const double epsilon = 1e-10;




// Atmosphere variables
// ...


static SphericalGrid *grid = NULL;
#ifdef ONLY_TRANSPORT_ON_REDUCED_GRID
static SphericalGrid *reducedgrid = NULL;
#endif
//static NetCDFOutput *diag_ncout = NULL;

// data storage (naming convention: lowercase)
static Topography                      *topo           = NULL;
static Temperature                     *temperature    = NULL;
static Humidity                        *humidity       = NULL;
static RadiativeFluxes                 *radiation      = NULL;
static LargeScaleWind                  *wind           = NULL;
static SynopticScaleField              *synoptic       = NULL;
static PlanetaryWave                   *standing_wave  = NULL;
static Clouds_2                        *clouds_2       = NULL;
static Clouds_3Layers                  *clouds_3       = NULL;
static SurfaceLayer                    *sf_layer       = NULL;
static PlanetaryBoundaryLayer          *pbl            = NULL;
//static EmpiricalData                 *era40          = NULL;
static Atm                             *atm            = NULL;

// computation schemes (naming convention: CAPITALS)
// dynamics:
static LargeScaleDynamics              *WIND           = NULL;
static PlanetaryBoundaryLayerPhysics   *PBL_PHYSICS    = NULL;
static PlanWave                        *PLANWAVE       = NULL;
static SynopticScaleDynamics           *SYNOPTIC       = NULL;
// conservation:
static HeatTransfer                    *ENERGY_CYCLE   = NULL;
static WaterVaporTransfer              *VAPOR_CYCLE    = NULL;
// physics:
static RadiativeTransfer               *RADIATION      = NULL;
static Cloudiness_2                    *CLOUDS_2       = NULL;
static Cloudiness_3Layers              *CLOUDS_3       = NULL;

// for handling heat and moisture flux with land surface, required for strict mass/energy conservation
static HorizontalVariable              *dT             = NULL;  ///< change in temperature of the lowest atmospheric layer 
static HorizontalVariable              *dFlux_dT       = NULL;  ///< derivative of temperature flux at the top of the lowest atmospheric layer
                                                                /// with respect to the temperature in that layer.
static HorizontalVariable              *dq             = NULL;  ///< change in spec. humidity of the lowest atmospheric layer 
static HorizontalVariable              *dFlux_dq       = NULL;  ///< derivative of spec. humidity flux at the top of the lowest atmospheric layer
                                                                /// with respect to the spec. humidity in that layer.
static HorizontalVariable              *dt_mass        = NULL;  ///< dt/mass, where dt = atmospheric time step (sec)
                                                                /// mass = mass of lowest atmospheric layer (Kg/m2)
static HorizontalVariable              *du             = NULL;  ///< change in u in PBL, just needed to calculate the implicit fluxes
static HorizontalVariable              *dv             = NULL;  ///< change in v in PBL, just needed to calculate the implicit fluxes
static HorizontalVariable              *qS_old         = NULL;  ///< stores qS (surface specific humidity) of previous timestep

// for handling function aeolus_get_bottom_mass_
//static HorizontalVariable              *T_bot          = NULL;  ///< near surface temperature (at lowest model level) [deg K]
//static HorizontalVariable              *q_bot          = NULL;  ///< near surface spec. humidity (at lowest model level) [deg K]
//static HorizontalVariable              *P_bot          = NULL;  ///< pressure at which atmos near surface values are assumed to be defined (at lowest model level)
static HorizontalVariable              *Z_bot          = NULL;  ///< height above the surface for the lowest model level
// for handling function aeolus_get_bottom_wind_
//static HorizontalVariable              *U_bot          = NULL;  ///< zonal wind component at lowest model level
//static HorizontalVariable              *V_bot          = NULL;  ///< meridional wind component at lowest model level


// Flag to indicate that the longitude range should be 0 to 360 deg instead of -180 to 180.
// in most GFDL-provided example data sets, the longitudes run from 
// 0 to 360 deg in the grid spec file for the atmos grid.
// Atmosphere_III currently uses -180 to 180.
// It might be possible to create a transfer grid spec also with -180 to 180 range,
// but for now we prefer to just re-use GFDL-provded example data sets unmodified.
// Thus we use this flag to indicate that all accesses to longitude via grid index i
// must be rotated by one half.
static bool lon_0_360    = true; ///< whether longitudes run from 0 to 360 in the grid spec file  DC: this is an Aeolus input parameter, should be loaded from file
//static bool without_topo = true; // run atmosphere with/without topography. DC: commented-out, this is an Aeolus input parameter

static bool use_sinsol = true;
static bool update_land_frac_always = false;
static double dt_mass_atm; // [s*m^2/kg]
static bool ocean_albedo_climber2 = false; ///< where ocean, override albedo from FMS with values as calculated in Climber-2 svat.f


// For tuning experiments, we need to partially decouple Aeolus from FMS/MOM,
// by using data from forcing files instead from the FMS coupler.
// The FMS data override mechanism is intended for such types
// of experiments.
// However, specifically of interest is surface humidity,
// which the coupler does not provide directly, but only via fluxes.
// And it is unclear how to set up data files for such flux values.
// Instead, we implement some kind of data override mechanism on our own here.
//static const char *override_qS_filename = NULL;
static NetCDFInput *override_qS_input = NULL;
static long override_qS_recordcount = -1;
static long override_qS_rec = 0;
Parameter<const char *> *override_qS = NULL;
void override_qS_setup();

// For the clouds module, we need synoptic scale vertical wind ww.
// This is currently not computed, but set to constant 0.0 .
// As workaround, we implement an override mechanism here.
static NetCDFInput *override_ww_input = NULL;
static long override_ww_recordcount = -1;
static long override_ww_rec = 0;
Parameter<const char *> *override_ww = NULL;
void override_ww_setup();

// Sigh.
// To send data from Aeolus reduced grid to FMS coupler regular grid, we must do
// some kind of interpolation.
// Since alle quantities are temperatures, speeds, or fractions, or fluxes (defined per m^2),
// for all of them the simplest way to handle that is Variable::INTERP_OUT_REPLICATE .
// However, we have observed strange zonal jumpy patterns in the polar regions.
// So, the next sophisiticate alternate is Variable::INTERP_OUT_INTERPOL .
// But with that, we experienced even more unexpected instabilities.
// Thus, here we can define for each of the passed-to-FMS variables a flag
// to switch it from INTERP_OUT_INTERPOL to INTERP_OUT_REPLICATE
// true : force INTERP_OUT_REPLICATE ; false : use INTERP_OUT_INTERPOL

Parameter<bool> *INTERP_gust = NULL;
Parameter<bool> *INTERP_flux_sw = NULL;
Parameter<bool> *INTERP_flux_sw_dir = NULL;
Parameter<bool> *INTERP_flux_sw_dif = NULL;
Parameter<bool> *INTERP_flux_sw_vis = NULL;
Parameter<bool> *INTERP_flux_sw_vis_dir = NULL;
Parameter<bool> *INTERP_flux_sw_vis_dif = NULL;
Parameter<bool> *INTERP_flux_sw_down_vis_dir = NULL;
Parameter<bool> *INTERP_flux_sw_down_vis_dif = NULL;
Parameter<bool> *INTERP_flux_sw_down_total_dir = NULL;
Parameter<bool> *INTERP_flux_sw_down_total_dif = NULL;
Parameter<bool> *INTERP_flux_lw = NULL;
Parameter<bool> *INTERP_surf_diff_delta_t = NULL;
Parameter<bool> *INTERP_surf_diff_dflux_t = NULL;
Parameter<bool> *INTERP_surf_diff_delta_tr = NULL;
Parameter<bool> *INTERP_surf_diff_dflux_tr  = NULL;
Parameter<bool> *INTERP_surf_diff_dtmass = NULL;
Parameter<bool> *INTERP_surf_diff_delta_u = NULL;
Parameter<bool> *INTERP_surf_diff_delta_v = NULL;
//Parameter<bool> *INTERP_surf_diff_delta_co2 = NULL;
//Parameter<bool> *INTERP_surf_diff_dflux_co2 = NULL;
Parameter<bool> *INTERP_lprec = NULL;
Parameter<bool> *INTERP_fprec = NULL;
// in aeolus_get_bottom_mass_()
Parameter<bool> *INTERP_z_bot = NULL;
Parameter<bool> *INTERP_t_bot = NULL;
Parameter<bool> *INTERP_tr_bot = NULL;
Parameter<bool> *INTERP_p_bot = NULL;
Parameter<bool> *INTERP_p_surf = NULL;
Parameter<bool> *INTERP_slp = NULL;
// in aeolus_get_bottom_wind_()
Parameter<bool> *INTERP_u_bot = NULL;
Parameter<bool> *INTERP_v_bot = NULL;


// We want to use the facilities of the FMS diagnostics manager for Aeolus-generated Variables.
// During initialization, Aeolus Variables are first registered for FMS diag output with FMSAddVariable().
// The number of registered Variables is passed back to the invoking atmosphere module.
// After the registrations, the atmosphere module iteratively calls aeolus_FMS_Variable_metadata_()
// to query the metadata of those Variables (and invokes the actual registration with the FMS diaga_manager module).
// During the simulation run, aeolus_FMS_Variable_data_() is invoked iteratively to
// obtain the current data of those Variables which are actually requested through the FMS diag_table.
// Advantages: Only calls F90 -> C++, not vice versa, thus assumed to be fairly robust.
//             Simple
//             About as run-time efficient as it could possibly get
// Disadvantage: works only with Variables which are constructed during init and live throughout the simulation run,
//               does not work with temporary / local Variables.
static vector<const Variable *> FMS_registered_vars;
static bool FMS_open_for_registration = true;

static void FMSAddVariable(const Variable& v);
static void FMSAddVariable(const LargeScaleWind *data );
static void FMSAddVariable(const SynopticScaleField *data );
static void FMSAddVariable(const Temperature *data );
static void FMSAddVariable(const Humidity *data );
static void FMSAddVariable(const Clouds_2 *data );
static void FMSAddVariable(const Clouds_3Layers *data );
static void FMSAddVariable(const PlanetaryWave *data );
static void FMSAddVariable(const PlanetaryBoundaryLayer *data );
static void FMSAddVariable(const Topography *data );
//static void FMSAddVariable(const EmpiricalData *data );
static void FMSAddVariable(const RadiativeFluxes *data );
static void FMSAddVariable(const Atm *data );
static void FMSAddVariable(const SurfaceLayer *data );


#ifdef TRACE_AEOLUS_UPDATE
static NetCDFOutput *Trace_Aeolus_Update_Down_In = NULL;
// Input to aeolus_atmosphere_down_()
static HorizontalVariable *u_star_down = NULL;
static HorizontalVariable *b_star_down = NULL;
static HorizontalVariable *q_star_down = NULL;
static HorizontalVariable *t_down = NULL;
static HorizontalVariable *dtaudu_down = NULL;
static HorizontalVariable *dtaudv_down = NULL;
static HorizontalVariable *u_flux_in_down = NULL;
static HorizontalVariable *v_flux_in_down = NULL;

// diagnostics from aeolus_atmosphere_down_()
static double stockflux_down_heat_in = 0.0;
/* no water stock is moved in downward sweep */

static NetCDFOutput *Trace_Aeolus_Update_Down_Out = NULL;
// Output from aeolus_atmosphere_down_()
static HorizontalVariable *u_flux_out_down = NULL;
static HorizontalVariable *v_flux_out_down = NULL;
static HorizontalVariable *gust_down = NULL;
static MeridionalVariable1D *coszen_down = NULL;
static const HorizontalVariable *flux_sw_sf_down = NULL;
static const HorizontalVariable *flux_sw_toa_down = NULL;
static HorizontalVariable *flux_sw_dir_down = NULL;
static HorizontalVariable *flux_sw_dif_down = NULL;
static HorizontalVariable *flux_sw_down_vis_dir_down = NULL;
static HorizontalVariable *flux_sw_down_vis_dif_down = NULL;
static HorizontalVariable *flux_sw_down_total_dir_down = NULL;
static HorizontalVariable *flux_sw_down_total_dif_down = NULL;
static HorizontalVariable *flux_sw_vis_down = NULL;
static HorizontalVariable *flux_sw_vis_dir_down = NULL;
static HorizontalVariable *flux_sw_vis_dif_down = NULL;
static HorizontalVariable *flux_lw_sf_down = NULL;
static const HorizontalVariable *flux_lw_toa_down = NULL;
static HorizontalVariable *dtmass_down = NULL;
static HorizontalVariable *dflux_t_down = NULL;
static HorizontalVariable *dflux_sphum_down = NULL;
static HorizontalVariable *delta_t_down = NULL;
static HorizontalVariable *delta_sphum_down = NULL;
static HorizontalVariable *delta_u_down = NULL;
static HorizontalVariable *delta_v_down = NULL;

// diagnostics from aeolus_atmosphere_down_()
static double stockflux_down_heat_out = 0.0;
/* no water stock is moved in downward sweep */

static NetCDFOutput *Trace_Aeolus_Update_Up_In = NULL;
// Input to aeolus_atmosphere_up_()
static HorizontalVariable *v_flux_up = NULL;
static HorizontalVariable *u_star_up = NULL;
static HorizontalVariable *b_star_up = NULL;
static HorizontalVariable *q_star_up = NULL;
static HorizontalVariable *dt_t_up = NULL;
static HorizontalVariable *dt_sphum_up = NULL;

// diagnostics from aeolus_atmosphere_up_()
static double stockflux_up_heat_in = 0.0;
static double stockflux_up_water_in = 0.0;

static NetCDFOutput *Trace_Aeolus_Update_Up_Out = NULL;
// Output from aeolus_atmosphere_up_()
//static HorizontalVariable *lprec_up = NULL;
//static HorizontalVariable *fprec_up = NULL;
static HorizontalVariable *gust_up = NULL;
static const HorizontalVariable *t_bot_up = NULL;
static const HorizontalVariable *sphum_bot_up = NULL;
static const HorizontalVariable *p_bot_up = NULL;
static const HorizontalVariable *z_bot_up = NULL;
static const HorizontalVariable *p_surf_up = NULL;
static const HorizontalVariable *slp_up = NULL;
static const HorizontalVariable *u_bot_up = NULL;
static const HorizontalVariable *v_bot_up = NULL;

// diagnostics from aeolus_atmosphere_up_()
static double stockflux_up_heat_out = 0.0;
static double stockflux_up_water_out = 0.0;

#endif


#ifdef DEBUG_AEOLUS_STOCKS
static ScalarVariable *stockmove_down_heat_toa = NULL;
static ScalarVariable *stockmove_down_heat_out  = NULL;
/* no water stock is moved in downward sweep */

static ScalarVariable *stockmove_up_heat_in   = NULL;
static ScalarVariable *stockmove_up_water_in  = NULL;
static ScalarVariable *stockmove_up_heat_out  = NULL;
static ScalarVariable *stockmove_up_water_out = NULL;
#endif

#ifdef TRACE_AEOLUS_INIT
NetCDFOutput *Trace_Aeolus_Init = NULL;
double trace_aeolus_init_counter = 0;
#endif

#ifdef USE_CLOCK
#include "Clock.h"
static Clock clock_Aeolus("Aeolus");
static Clock clock_Aeolus_down("Aeolus_down");
static Clock clock_Aeolus_up("Aeolus_up");
#endif

#include "fmsclock.h"
static int clock_id_aeolus = -1;
static int clock_id_aeolus_down = -1;
static int clock_id_aeolus_up = -1;
//static int clock_id_aeolus_getstocks = -1;

// declare functions:
void LoadTopologyFromF90( const double *const landfrac, Topography& topo);
void LoadAlbedosFromF90  (  const double *const albedo_vis_dir, 
                            const double *const albedo_nir_dir,
                            const double *const albedo_vis_dif,
                            const double *const albedo_nir_dif, 
                            const double *const roughness,
                            SurfaceLayer& sf,
                            const double *const cosz = NULL);
void InitSurfaceTypesFromLandOceanMask(    const HorizontalVariable& fraction_ocean, const HorizontalVariable& fraction_land, SurfaceLayer& sf );
void LoadTemperatureFromF90( const double *const ts, Temperature& temperature);
//void LoadSurfaceFluxesFromF90( const double *const delta_t, const double *const delta_q, SurfaceLayer& sf);
void LoadSurfaceFluxesFromF90( const HorizontalVariable& delta_t, const HorizontalVariable& delta_q, SurfaceLayer& sf);

// Compare double values for equality.
// They are considered equal if the relative error is smaller than epsilon.
inline static bool double_eq(const double v1, const double v2, const double e) {
  return fabs((v1-v2)/((v1+v2)/2)) < e;
}


/* we try to override abort() to get FMS diag_manager output flushed to disk */
#ifndef __THROW
#define __THROW /* empty */
#endif
extern "C" void abort_with_flush_(void);
void abort(void) __THROW {
#if defined(WITH_SIMENV) && defined(WITH_COSTFUNCTION)
  /* avoid Missing-Values in the cost function output */
  /* For SimEnv, we need to output the worst possible score, which is 1 */
  /* But on stdout we print out the worst possible taylor skill value, which is 0 */
  double one = 1.0;
  simenv_put_c("skillscore",(char *)&one);
  cout << "SkillScore " << 0 << " aborted " << simenv_run_int << 1 << endl;
#endif

  // lets try to call back into Fortran90
  fprintf(stderr, "abort overriden\n");
  abort_with_flush_();
  exit(42);
}

#ifdef WITH_COSTFUNCTION
/**
 * #ifdef WITH_COSTFUNCTION a pice of code named Aeolus_SkillScore.cpp is included.
 * It must provide the functions
 * void init_Aeolus_Score()
 * void Aeolus_Score_Add_Tstep(const int year, const int month)
 * void Aeolus_Score_Calc()
 * void Aeolus_Score_End(void)
 * There are several implementations in the SVN repository.
 * To choose one for that is suitable for your intended purpose,
 * create a symbolic link named Aeolus_SkillScore.cpp which to points to your choice.
 */
#include "Aeolus_SkillScore.cpp"
#endif

extern "C" {

  // This block contains those functions that are called by the F90-implemented FMS.
  // It is tested with our two main target platforms GNU gfortran with g++ and Intel Fortran
  // with Intel C++, both on Linux.
  // For other Fortran/C++ combos, different naming and argument passing conventions might
  // be needed.
  // GNU and Intel Fortran convert names of external functions to lowercase and append one
  // underscore. (This their default behaviour, can be changed with command line options and/or
  // special directive comments in the code.


  /// aeolus_init_grid_()  first of all is used to check for
  /// compatibility of the procedure call interface between
  /// the Fortran90 and C++ modules.
  /// Possible pitfalls include:
  /// - size of Fortran integer vs. sizeof C++ int
  /// - size of Fortran real vs. sizeof C++ float or double
  /// - assumed size arrays and pointers might be passed as
  ///   descriptors (alias dope vectors)
  /// - alignment of elements in derived types
  ///   => to avoid this problem, we try to avoid passing
  ///      derived types, and use the single elements as call
  ///      arguments (we shall see if this will always be possible)
  ///
  /// Our primary target platforms are
  /// Linux with Intel Compilers,
  /// and Linux with GNU compilers.
  void aeolus_init_grid_
  ( /// first some arguments for checking the F90->C++ interface
   const int *const deadbeef, const double *const pitest,
   const int *const true_in, const int *const false_in,
   const double *const testarray3d,
   const int *const kind_i, const int *const kind_r, const int *const kind_l,

   /// now the ``real'' arguments

   // inout: the fields of the surf_diff_type
   // TODO: do we really need to pass them here? - does not seem so.
 
   /// the model parameter from the namelist
   const int *const NoNx, const int *const NoNy, const int *const NoNz, const int *const NoNz_Strato,
   const int *const polar_merging,
   const int *const lon_0_360_in,

   /// information about domain partitioning
   const int *const peset_in,   ///< MPI communicator for the set of atmoshpere pes
   const int *const pe_in,      ///< our own MPI rank == processing element number
   const int *const root_pe_in, ///< MPI rank of atmosphere master process
   const int *const npes_in,    ///< number of total atmosphere processes
   /// the boundary indices of the local compute domain
   const int *const is_in, const int *const ie_in, const int *const js_in, const int *const je_in,

   /// clock id
   const int *const clock_id_aeolus_in,

   /// and finally some interface-checking args again
   const double *const onetwothree, const int *const fse
    ) { // begin aeolus_init_grid_()

    const double myonetwothree = 1234567809.0987654321L;

    //When mtrace() is called, it checks the value of the environment variable MALLOC_TRACE, which  should  contain
    //the  pathname  of a file in which the tracing information is to be recorded.  If the pathname is successfully
    //opened, it is truncated to zero length.
    //
    //If MALLOC_TRACE is not set, or the pathname it specifies is invalid or not writable, then no  hook  functions
    //are  installed,  and  mtrace()  has  no  effect.
    mtrace();
#ifdef USE_CLOCK
    clock_Aeolus.start();
#endif
    clock_id_aeolus = *clock_id_aeolus_in;
    fmsclock_begin_(&clock_id_aeolus);

    printf("%s\n", "$Id$");
    // TODO: perform this test and all the printout only on the atmosphere ``root'' process
    // to avoid cluttering of the output logs
    printf("%s %s F90->C++ interface test\n", __FILE__, __FUNCTION__); fflush(stdout);
    printf("Fortran KIND(integer) %d KIND(real) %d KIND(logical) %d C++ sizeof int %d long %d float %d double %d bool %d\n",
	   *kind_i, *kind_r, *kind_l, (int)sizeof(int), (int)sizeof(long), (int)sizeof(float), (int)sizeof(double),
	   (int)sizeof(bool));
    fflush(stdout);
    if (*kind_i != sizeof(int)) { printf("Error: Fortran KIND(integer) != C++ sizeof(int)\n"); fflush(stdout); abort(); }
    if (*kind_r != sizeof(double)) { printf("Error: Fortran KIND(real) != C++ sizeof(double)\n"); fflush(stdout); abort(); }
    if (*kind_l != sizeof(int)) { printf("Error: Fortran KIND(logical) != C++ sizeof(int)\n"); fflush(stdout); abort(); }
    if (*deadbeef != 0xeadbeef /* 0xdeadbeef yields overflow error on 32bit machines */) {
      printf("Error: test1 failed: 0xdeadbeef != %x\n", *deadbeef);
      fflush(stdout);
      abort();
    }
    if (! double_eq(*pitest, M_PIl, epsilon)) {
      printf("Error: test2 failed: %30.28f != %30.28f\n", M_PI, *pitest); fflush(stdout);
      abort();
    }
    if (! (*true_in)) { printf("Error: test3 failed: true_in %x is not true\n", *true_in);
      fflush(stdout);
      abort();
    }
    if (*false_in) { printf("Error: test4 failed: false_in %x is not false\n", *false_in);
      fflush(stdout);
      abort();
    }

    for (int k = 1; k <= 4; k++) {
      for (int j = 1; j <= 3; j++) {
	for (int i = 1; i <= 2; i++) {
	  const int idx = (k-1)*3*2 + (j-1)*2 + (i-1);
	  const double testval = testarray3d[idx];
	  const double myval = i*100.0+j+k/10.0;
	  //printf("testarray3d(%2d,%2d,%2d) = testarray3d[%2d] = %10.9g vs. %10.9g\n",
	  //	 i, j, k, idx, testval, myval);
	  if (! double_eq(testval, myval, epsilon)) {
	    printf("Error: test5 failed: testarray3d(%d,%d,%d) is %10.9g != %10.9g\n", i, j, k, testval, myval);
	    fflush(stdout);
	    abort();
	  }
	}
      }
    }

    if (! double_eq(*onetwothree, myonetwothree, epsilon)) {
      printf("Error: test6 failed: %24.12f != %24.12f\n", myonetwothree, *onetwothree);
      fflush(stdout);
      abort();
    }
    if (*fse != 4711) { printf("Error: test7 failed: 4711 != %d\n", *fse); fflush(stdout); abort(); }
    //printf("double HUGE_VAL is %g\n", HUGE_VAL);
    printf("... interface test finished\n"); fflush(stdout);

#ifdef _OPENMP
    int tid, nthreads;
#pragma omp parallel shared(nthreads) private(tid)
    {
      tid = omp_get_thread_num();
      if (tid == 0) { /* Only master thread does this */
        nthreads = omp_get_num_threads();
        printf("Running with %d OpenMP threads\n", nthreads);
      }
    }
#else
    printf("Compiled without OpenMP\n");
#endif

    clock_id_aeolus_down = fmsclock_id("Aeolus_down");
    clock_id_aeolus_up   = fmsclock_id("Aeolus_up");
    //clock_id_aeolus_getstocks = fmsclock_id("Aeolus_getstocks");

#ifdef WITH_SIMENV
    { int simenv_sts;
      char  simenv_run_char[6];

      simenv_sts = simenv_ini_safe_c();
      if (simenv_sts != 0) { fprintf(stderr, "simenv_ini_c() returned %d\n", simenv_sts); exit(simenv_sts); }
      /* get current run number */ 
      simenv_get_run_c(&simenv_run_int,simenv_run_char);
      printf("SimEnv model run %d\n", simenv_run_int);

      /* to avoid (slow) creation of symlinks for every SimEnv run */
      //NetCDFInput::filenameprefix = "../";
    }
#endif

    //
    // now look at the actual model parameters
    //

    is = *is_in; ie = *ie_in; js = *js_in; je = *je_in;
    comm = *peset_in; pe = *pe_in; root_pe = *root_pe_in; npes = *npes_in;
    lon_0_360 = *lon_0_360_in;
    printf("Grid resolution: lon %d lat %d vertical %d+%d %s\n", *NoNx, *NoNy, *NoNz, *NoNz_Strato,
	   (*polar_merging?" Merged Cells":" regular grid"));
    printf("grid_spec longitudes run from 0 to 360? -> %s\n", lon_0_360 ? "Yes" : "No");
    printf("MPI: comm %d 0x%x pe %d root_pe %d npes %d\n", comm, comm, pe, root_pe, npes);
#ifdef USE_MPI_DOMAIN_DECOMPOSE
    printf("MPI_COMM_WORLD 0x%x MPI_COMM_SELF 0x%x MPI_COMM_NULL 0x%x\n", MPI_COMM_WORLD, MPI_COMM_SELF, MPI_COMM_NULL);
    printf("MPI_GROUP_EMPTY 0x%x MPI_GROUP_NULL 0x%x\n", MPI_GROUP_EMPTY, MPI_GROUP_NULL);
#else
    if (npes > 1) {
      printf("Error: %s was not compiled with USE_MPI_DOMAIN_DECOMPOSE #defined, but is running on %d MPI tasks\n",
              __FILE__, npes);
      fflush(stdout);
      abort();
    }
#endif

    printf("Domain: local boundary indices: is %d ie %d js %d je %d\n", is, ie, js, je);

#ifdef USE_MPI_DOMAIN_DECOMPOSE
    Aeolus_Communicator = comm; // implicit invocation of copy constructor which converts handle to object
#endif
    grid = new SphericalGrid(*NoNx, *NoNy, *NoNz, *NoNz_Strato, *polar_merging, lon_0_360 /* use 0to360 deg */
#ifdef USE_MPI_DOMAIN_DECOMPOSE
                             , Aeolus_Communicator
#endif
                             );
    //grid->OutputToScreen(); // very verbose

#ifdef ONLY_TRANSPORT_ON_REDUCED_GRID
    if (*polar_merging) {
      fprintf(stderr, "WARNING: ONLY_TRANSPORT_ON_REDUCED_GRID is #defined, but polar_merging is true - are you sure?\n");
      fprintf(stderr, "   No separate reduced grid is generated, and Advection is done on the ``main'' grid of Aeolus\n");
    } else {
      reducedgrid = new SphericalGrid(*NoNx, *NoNy, *NoNz, *NoNz_Strato, true, lon_0_360 /* use 0to360 deg */, comm);
      //grid->OutputToScreen(); // very verbose
    }
#endif

#ifdef USE_CLOCK
    clock_Aeolus.stop();
#endif
    fmsclock_end_(&clock_id_aeolus);
  } // end aeolus_init_grid_()


  void aeolus_init_
  (
   /// the FMS time stamps are broken out into days and seconds
   const int *const Time_init_year, const int *const Time_init_month,
   const int *const Time_init_day, const int *const Time_init_seconds,
   const int *const Time_days, const int *const Time_seconds,
   const int *const Time_step_days, const int *const Time_step_seconds,

   /// the model parameter from the namelist
   const int *const without_topo_in,

   const double *const HORO_in,   ///< Orography height, averaged over grid cell
   const double *const SIGORO_in, ///< standard deviation of topography on local domain
   //const double *const land_frac_in, ///< preliminary land fraction per grid cell
   // the actual land_frac is only passed in the downward sweep, and not needed during init

   /// the model parameter from the namelist
   const int *use_sinsol_in,
   const int *update_land_frac_always_in,
   const int *ocean_albedo_climber2_in,

   /// information about tracer indices in surf_diff_type
   const int *const nr_tracers_in, ///< total number of tracers in surf_diff_type
   // TODO: if more than one tracer, we probably need to pass an array of tracer indices
   const int *const q_ind_in, ///< index of humidity in tracer fields
   const int *const co_ind_in, ///< index of CO2 in tracer fields

   /// output variable: number of Aeolus variables that will be registered with the FMS
   /// FMS diagnostics manager lateron
   int *const nr_Aeolus_diag_vars,

   /// interface checking
   const int *const fse
    ) { // begin Aeolus_init_()
#ifdef USE_CLOCK
    clock_Aeolus.start();
#endif
    fmsclock_begin_(&clock_id_aeolus);
    if (*fse != 4711) { printf("Error: test8 failed: 4711 != %d\n", *fse); fflush(stdout); abort(); }
    //printf("double HUGE_VAL is %g\n", HUGE_VAL);
    printf("... interface test finished\n"); fflush(stdout);

    printf("Time_init %05d-%02d-%02d:%05dsec\n", *Time_init_year, *Time_init_month, *Time_init_day, *Time_init_seconds);
    printf("Time      %05dY%02dD:%05dsec\n", (*Time_days)/356, (*Time_days)%356, *Time_seconds);
    printf("Time_step %d:%05d days:sec\n", *Time_step_days, *Time_step_seconds);
#ifndef NDEBUG
    first_day = *Time_days;
#endif

    nr_tracers = *nr_tracers_in;
    q_ind = *q_ind_in;
    co2_ind = *co_ind_in;
    printf("total number of tracers in surf_diff_type fields: %d; index of humidity: %d; of CO2: %d\n",
	   nr_tracers, q_ind, co2_ind);
    fflush(stdout);

    use_sinsol = (*use_sinsol_in != 0);
    printf("Using solar data from %s module\n", use_sinsol ? "Climber-2 sinsol" : "FMS astronomy");
    if (!use_sinsol) FMSAstronomy::Init(*grid);

    update_land_frac_always = (*update_land_frac_always_in != 0);
    ocean_albedo_climber2 = (*ocean_albedo_climber2_in != 0);

    NetCDFInput::filenameprefix = "INPUT/";

    NML *Aeolus_nml = DefineNML("Aeolus_nml");
    
    override_qS = new Parameter<const char *>(NULL, "override_qS", Aeolus_nml, "name of forcing file, default NULL");
    override_ww = new Parameter<const char *>(NULL, "override_ww", Aeolus_nml, "name of forcing file, default NULL");
    INTERP_gust = new Parameter<bool>(false, "INTERP_gust", Aeolus_nml, "const 1");
    INTERP_flux_sw = new Parameter<bool>(false, "INTERP_flux_sw", Aeolus_nml);
    INTERP_flux_sw_dir = new Parameter<bool>(false, "INTERP_flux_sw_dir", Aeolus_nml);
    INTERP_flux_sw_dif = new Parameter<bool>(false, "INTERP_flux_sw_dif", Aeolus_nml);
    INTERP_flux_sw_vis = new Parameter<bool>(false, "INTERP_flux_sw_vis", Aeolus_nml);
    INTERP_flux_sw_vis_dir = new Parameter<bool>(false, "INTERP_flux_sw_vis_dir", Aeolus_nml);
    INTERP_flux_sw_vis_dif = new Parameter<bool>(false, "INTERP_flux_sw_vis_dif", Aeolus_nml);
    INTERP_flux_sw_down_vis_dir = new Parameter<bool>(false, "INTERP_flux_sw_down_vis_dir", Aeolus_nml, "const 0");
    INTERP_flux_sw_down_vis_dif = new Parameter<bool>(false, "INTERP_flux_sw_down_vis_dif", Aeolus_nml, "const 0");
    INTERP_flux_sw_down_total_dir = new Parameter<bool>(false, "INTERP_flux_sw_down_total_dir", Aeolus_nml, "const 0");
    INTERP_flux_sw_down_total_dif = new Parameter<bool>(false, "INTERP_flux_sw_down_total_dif", Aeolus_nml, "const 0");
    INTERP_flux_lw = new Parameter<bool>(false, "INTERP_flux_lw", Aeolus_nml);
    INTERP_surf_diff_delta_t = new Parameter<bool>(false, "INTERP_surf_diff_delta_t", Aeolus_nml);
    INTERP_surf_diff_dflux_t = new Parameter<bool>(false, "INTERP_surf_diff_dflux_t", Aeolus_nml, "const 0");
    INTERP_surf_diff_delta_tr = new Parameter<bool>(false, "INTERP_surf_diff_delta_tr", Aeolus_nml);
    INTERP_surf_diff_dflux_tr = new Parameter<bool>(false, "INTERP_surf_diff_dflux_tr", Aeolus_nml, "const 0");
    INTERP_surf_diff_dtmass = new Parameter<bool>(false, "INTERP_surf_diff_dtmass", Aeolus_nml, "const uniform");
    INTERP_surf_diff_delta_u = new Parameter<bool>(false, "INTERP_surf_diff_delta_u", Aeolus_nml, "const 0");
    INTERP_surf_diff_delta_v = new Parameter<bool>(false, "INTERP_surf_diff_delta_v", Aeolus_nml, "const 0");
    //INTERP_surf_diff_delta_co2 = new Parameter<bool>(false, "INTERP_surf_diff_delta_co2", Aeolus_nml);
    //INTERP_surf_diff_dflux_co2 = new Parameter<bool>(false, "INTERP_surf_diff_dflux_co2", Aeolus_nml);
    INTERP_lprec = new Parameter<bool>(false, "INTERP_lprec", Aeolus_nml);
    INTERP_fprec = new Parameter<bool>(false, "INTERP_fprec", Aeolus_nml);
    // in aeolus_get_bottom_mass_()
    INTERP_z_bot = new Parameter<bool>(false, "INTERP_z_bot", Aeolus_nml, "const uniform");
    INTERP_t_bot = new Parameter<bool>(false, "INTERP_t_bot", Aeolus_nml);
    INTERP_tr_bot = new Parameter<bool>(false, "INTERP_tr_bot", Aeolus_nml);
    INTERP_p_bot = new Parameter<bool>(false, "INTERP_p_bot", Aeolus_nml);
    INTERP_p_surf = new Parameter<bool>(false, "INTERP_p_surf", Aeolus_nml);
    INTERP_slp = new Parameter<bool>(false, "INTERP_slp", Aeolus_nml);
    // in aeolus_get_bottom_wind_()
    INTERP_u_bot = new Parameter<bool>(false, "INTERP_u_bot", Aeolus_nml);
    INTERP_v_bot = new Parameter<bool>(false, "INTERP_v_bot", Aeolus_nml);
    ReadNML (Aeolus_nml, "input.nml"  );

    Variable::SetInvalidIsNotFatal();

    // Now decide whether to read initial data from files,
    // or restore previously saved state from restart file.
    // Check if there is a valid restart file to read saved state from
    unsigned int restNoNx, restNoNy, restNoNz, restNoNz_Strato;
    bool restpolar_merging, restlon_0_360;
    double restTime;
    const char *readresfilename = "INPUT/aeolus.res.nc";
    if (access("RESTART/aeolus.res.nc", R_OK) < 0
        // maybe we still have a restart file name with old name?
        && 0 == access("INPUT/aeolus.res.nc", R_OK)) {
      // read from old file name
      readresfilename = "INPUT/aeolus.res.nc";
    }
    bool restart_from_saved_state =
      restart_init("RESTART/aeolus.res.nc", // new state to be saved
                   "Restart file for Aeolus Atmosphere coupled to FMS $Id$",
                   readresfilename,  // previous state to read
                   restNoNx, restNoNy, restNoNz, restNoNz_Strato,
                   restpolar_merging, restlon_0_360,
                   restTime);

    if (restart_from_saved_state) {
      // check that the grid parameters from restart file match those from input namelist
      if (grid->ZonalCells() != restNoNx || grid->MeridionalCellsGlobal() != restNoNy
          || grid->TroposphereLevels() != restNoNz || grid->StratosphereLevels() != restNoNz_Strato
          /* Fortran logical might use -1 for .true. Both C and Fortran use 0 for .false. */
	  || !(grid->IsRegular()) != restpolar_merging
          || grid->Is0to360() != restlon_0_360) {
	printf("Error: some Grid parameters from restart file %s do not match those from FMS input namelist\n",
	       readresfilename);
	printf("  %u / %u ; %u / %u ; %u / %u ; %u / %u ; %d / %d ; %d / %d\n",
	       grid->ZonalCells(), restNoNx, grid->MeridionalCellsGlobal(), restNoNy,
               grid->TroposphereLevels(), restNoNz, grid->StratosphereLevels(), restNoNz_Strato,
	       !(grid->IsRegular()), restpolar_merging, grid->Is0to360(), restlon_0_360);
	fflush(stdout);
	abort();
      }

      // the timestamp read from restart file must match the current time passed as parameter by FMS
      if (restTime != (double)(*Time_days) + (double)(*Time_seconds)/(double)(24*60*60)) {
	printf("Error: Timestamp read from restart file %s does not match that passed from FMS: %15.8g != %15.8g\n",
	       readresfilename, restTime, (double)(*Time_days) + (double)(*Time_seconds)/(double)(24*60*60));
	fflush(stdout);
	abort();
      }
      // Everything looks fine so far.
    }

    // allocate atmosphere variables

    // constant timestepping
    dt_atm_sec = static_cast<double>(*Time_step_seconds);
    //debug_out_counter = 0;

    // print out the maximum permissible wind speeds for CFL stability for this grid and dt_atm_sec
    // sigh. with the new advection_dt options, we would use not the atmos time step dt_atm_sec,
    // but the advection_dt.
    // However, that is currently local to the HeatTransfer and WaterVaporTransfer modules.
    double u_max, v_max, w_max;
    grid->Max_Winds(dt_atm_sec, &u_max, &v_max, &w_max);

    // Initialise memory for all variables, NST = number of surface types
    // all created classes contain several 2D or 3D variables capturing certain physics
    // convention: use regular fonts for datastorage classes and CAPITALS for computation classes
    //topo           = new Topography             ( *grid );
    topo           = new Topography             ( *grid );
    temperature    = new Temperature            ( *grid );
    humidity       = new Humidity               ( *grid );
    radiation      = new RadiativeFluxes        ( *grid/*, NST*/ );
    // set validity ranges according to time step and grid resolution
    // Default in DataStorage: 300, 300, 10
    // Regular grid, 96x60 cells, dt=60 sec
    //maximum permissible wind speeds [m/s] for CFL stability at timestep length 60 [s]:
    //u  181.822 at  1.88 -88.50 1000.00
    //v 5556.692 at  1.88 -88.50 1000.00
    //w   33.333 at  1.88 -88.50 1000.00
    wind           = new LargeScaleWind         ( *grid, 180 /*, v_max, w_max*/ );
    synoptic       = new SynopticScaleField     ( *grid );
    standing_wave  = new PlanetaryWave          ( *grid );
    if (false) // TODO: depend on namelist parameter
      clouds_2     = new Clouds_2               ( *grid );
    else
      clouds_3     = new Clouds_3Layers         ( *grid );
    sf_layer       = new SurfaceLayer           ( *grid/*, NST*/ );
    pbl            = new PlanetaryBoundaryLayer ( *grid );
    //era40        = new EmpiricalData          ( *grid ); // needed only for initial testing
    atm            = new Atm                    ( *grid );
    
    // atmosphere variables are registered for writing restart information
    // in the c-tors of the DataStorage classes directly

#ifdef TRACE_AEOLUS_INIT
    Trace_Aeolus_Init = new NetCDFOutput(*grid, "trace_aeolus_init.nc", false,
                                         "Tracing computation of Variables during initialisation of Aeolus coupled to FMS $Id$",
                                         "seconds since 0000-01-01 00:00");
#endif

    // The FMS topography is then used for Aeolus.
    // Unfortunately, land_fraction is not passed by FMS to the atmosphere during init,
    // but in every up and down update step anew.
    // Another purpose here is to check interface consistency on both sides,
    // that the topography matches and is not rotated / mirrored / etc.
    LoadTopographyFromF90(HORO_in, SIGORO_in, *without_topo_in != 0, *topo);
    //!TODO: make smoothing of topography selectable by namelist parameter
    //LoadLandfracFromF90(land_frac_in, *topo);
    { // open a new scope for toponc
      NetCDFOutput toponc(*grid, "aeolus_diag_topo", false,
                          "Topography of Aeolus Atmosphere coupled to FMS $Id$");
      // for comparison, write out the Atmosphere-Internally generated topography
      toponc.AddVariable(topo->real_height_());
      toponc.AddVariable(topo->real_sigma_());
      toponc.AddVariable(topo->height_());
      toponc.AddVariable(topo->sigma_());
      //toponc.AddVariable(topo->fraction_land_());
      //toponc.AddVariable(topo->fraction_ocean_());
    
      toponc.OutputToFile(0);
      // end of life for toponc - should be closed and destroyed automagically here
    }

    // adapt grid for topography. If not using topography, adapt for height 0
    //!TODO: make this adaptation selectable via namelist parameter
    /*if (*without_topo_in == 0)*/ //grid->AdaptForTopography(topo->height_());
    /*else*/ printf("NOT adapting Cells and Interfaces for topography\n");

    // prepare diagnostic output
    char days_since[100];
    snprintf(days_since, sizeof(days_since)-1, "days since %04d-%02d-%02d 00:00",
             *Time_init_year, *Time_init_month, *Time_init_day);
    //diag_ncout = new NetCDFOutput(*grid, "aeolus_diag.nc", false /* not a restart file */,
    //				  "Diagnostic output for Aeolus Atmosphere coupled to FMS $Id$",
    //                              days_since, "NOLEAP");

    // Initialise required computations
    // dynamics
    PLANWAVE       = new PlanWave		               ( *grid, *temperature, *wind, *synoptic, *pbl, *atm, *topo, *standing_wave );
    PBL_PHYSICS    = new PlanetaryBoundaryLayerPhysics ( *grid, false, *wind, *synoptic, *sf_layer, *temperature, *topo, *atm, *humidity, *pbl );
    SYNOPTIC       = new SynopticScaleDynamics         ( *grid, *temperature, *humidity, *wind, *atm, *pbl, *topo, *clouds_3, *radiation, *synoptic );
    // physics
    if (clouds_2) {
      printf("Error: class LargeScaleDynamics now requires 3-layer clouds, but Climber-2 clouds are selected\n");
      fflush(stdout);
      abort();
      //CLOUDS_2     = new Cloudiness_2                  ( *grid, *wind, *synoptic, *temperature, *humidity, *clouds_2 );
    } else {
      WIND         = new LargeScaleDynamics            ( *grid, *standing_wave, *temperature, *synoptic, *pbl, *topo, *atm, *clouds_3, *PBL_PHYSICS, *wind );
      CLOUDS_3     = new Cloudiness_3Layers            (*grid, *topo, *wind, *synoptic, *standing_wave, *temperature, *humidity, *pbl, *radiation, *clouds_3); // radiation is not used in constructor
    }
    RADIATION      = new RadiativeTransfer             ( *grid, *temperature, *humidity, *wind, clouds_2, clouds_3, *topo, *pbl, *sf_layer, *radiation, use_sinsol );
    VAPOR_CYCLE    = new WaterVaporTransfer            ( *grid, *wind, *synoptic, *radiation, *temperature, clouds_2, clouds_3, *atm, *topo, *sf_layer, *humidity
#ifdef ONLY_TRANSPORT_ON_REDUCED_GRID
                                                         , reducedgrid
#endif
                                                         );
    ENERGY_CYCLE   = new HeatTransfer                  ( *grid, *wind, *synoptic, *radiation, *humidity, clouds_2, clouds_3, *atm, *topo, *pbl, *sf_layer, *temperature
#ifdef ONLY_TRANSPORT_ON_REDUCED_GRID
                                                         , reducedgrid
#endif
                                                         );

    if (restart_from_saved_state)
    {
      // read the values for the registered variables
      restart_do_restore();
      // set fraction of surface types: "open water" (0) set to "fraction ocean" and "trees" (2) set to fraction land. All others set to 0
      //InitSurfaceTypesFromLandOceanMask( topo->fraction_ocean_(), topo->fraction_land_(), *sf_layer );
      VAPOR_CYCLE->CalculateHe();
      WIND->InterpolateWindFieldToInterfaces(); // needed since interfaces cannot be stored in restart-file
      ENERGY_CYCLE->DerivedVariables();
    } 
    else
    {
      // Cold start.
      // reading TS and qs from file
      // Careful: from the CM2M restart file, we get data at 993mbar height, i.e. sea level.
      // But currently, Aeolus seems to expect surface height data here.
      //! TODO: rename the input variables TS and QS in the NetCDF file tos T_sl and Q_sl, to clarify their semantic
      VAPOR_CYCLE->LoadQSFromFileQSL("AEOLUS_input_1976-2001_yseasmean.cdf", "QS", 0, topo->height_() ); // reads sealevel humidity from and converts it to surface humidity in humidity->qS
      //! TODO: rename the input variable TS in the NetCDF file to T_sl, to clarify its semantic
      // This function reads sealevel temperature in K from a file,
      // converts it to C, computes Lapse Rate, and then surface temperature TSa in K.
      ENERGY_CYCLE->LoadT_slFromFile("AEOLUS_input_1976-2001_yseasmean.cdf", "TS", 0, false ); // true = smooth, else do not smooth

      // now initialize all compute modules using TS and qs loaded from file
      InitPresentDayClimate ( *grid, temperature->TSa_(), humidity->qS_(), *PBL_PHYSICS, *PLANWAVE, *SYNOPTIC, *WIND, *ENERGY_CYCLE, *VAPOR_CYCLE, *RADIATION, CLOUDS_2, CLOUDS_3, false, true );
      // set fraction of surface types: "open water" (0) set to "fraction ocean" and "trees" (2) set to fraction land. All others set to 0
      InitSurfaceTypesFromLandOceanMask( topo->fraction_ocean_(), topo->fraction_land_(), *sf_layer );
    }

    double z_level (10.); // level at which near-surface variables are defined
    dt_mass_atm = dt_atm_sec/Atmosphere.H0()/Atmosphere.rho_0(); // timestep divided by mass of atm column [s*m^2/kg]

    // initialise FMS flux exchange terms
    dT       = new HorizontalVariable ( *grid, SF, "dT", "K", "change in temperature of the lowest atmospheric layer", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_INTERPOL);
    dFlux_dT = new HorizontalVariable ( *grid, SF, "dFlux_dT", "?", "derivative of temperature flux", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_INTERPOL);
    dq       = new HorizontalVariable ( *grid, SF, "dq", "", "change in spec. humidity of the lowest atmospheric layer", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_INTERPOL /*?*/);
    dFlux_dq = new HorizontalVariable ( *grid, SF, "dFlux_dq", "?", "derivative of spec. humidity flux", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_INTERPOL /*?*/);
    dt_mass  = new HorizontalVariable ( *grid, SF, "dt_mass", "s m2/Kg", "dt / mass of lowest atmospheric layer", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_INTERPOL /*?*/); 
    du       = new HorizontalVariable ( *grid, PBL, "du", "m/s", "change of u in PBL", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_INTERPOL); 
    dv       = new HorizontalVariable ( *grid, PBL, "dv", "m/s", "change of v in PBL", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_INTERPOL); 
    Z_bot    = new HorizontalVariable ( *grid, SF, "Z", "m", "height of lowest model level", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_INTERPOL);
    qS_old   = new HorizontalVariable ( *grid, SF, "qS_old", "", "spec. humidity of the lowest atmospheric layer of previous timestep", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE /*?*/);

    // fill FMS flux exchange terms with initial values
    dT->CopyDataFrom(temperature->realTSa_()); // set equal to 2-meter temperature
    dq->CopyDataFrom(humidity->qS_() );    // set equal to q in bottom layer [kg/kg]
    // following remain constant throughout simulation
    dFlux_dT->UniformValue(0.);            // set to zero: bottom flux is added to complete atm. column. Fluxes at the top are thus zero.
    dFlux_dq->UniformValue(0.);            // assumed zero, like dFlux_dT
    // dt/mass, where dt = atmospheric time step (sec) and mass = mass of atmospheric layer (Kg/m2) to which surface fluxes are added
    dt_mass->UniformValue(dt_mass_atm);    // assume fluxes are added to full atmospheric column [s m2/Kg]

    // du / dv = u / v stress on atmosphere 
    // only needed to calculate the implicit fluxes
    // in EBL: set to 0.
    du->UniformValue(0.);
    dv->UniformValue(0.);
    Z_bot->UniformValue( z_level );


    // set up data override facility
    if ((*override_qS)()) {
      override_qS_setup();
      override_qS_rec = (*Time_days)%365 /*+ 1*/; // records are numbered from 0
    }
    if ((*override_ww)()) {
      override_ww_setup();
      // our forcing file for ww has seasonal values
      override_ww_rec = -1; // force values to be read on first time step
    }

    // add all standard fields to output (this is A LOT (!), should be reduced after initial testing)
    // -> we prefer to do diag output via the FMS diagnostics manager
    //AddOutputVariables(*diag_ncout, *topo); static field, in separate file
    //AddOutputVariables(*diag_ncout, *temperature); 
    //AddOutputVariables(*diag_ncout, *humidity); 
    //AddOutputVariables(*diag_ncout, *radiation); 
    //AddOutputVariables(*diag_ncout, *wind);
    //AddOutputVariables(*diag_ncout, *synoptic); 
    //AddOutputVariables(*diag_ncout, *standing_wave); 
    //if (clouds_2) AddOutputVariables(*diag_ncout, *clouds_2);
    //if (clouds_3) AddOutputVariables(*diag_ncout, *clouds_3);
    //AddOutputVariables(*diag_ncout, *sf_layer);
    //AddOutputVariables(*diag_ncout, *pbl);
    //AddOutputVariables(*diag_ncout, *atm);

    FMSAddVariable(topo);
    FMSAddVariable(temperature); 
    FMSAddVariable(humidity); 
    FMSAddVariable(radiation); 
    FMSAddVariable(wind);
    FMSAddVariable(synoptic); 
    FMSAddVariable(standing_wave); 
    if (clouds_2) FMSAddVariable(clouds_2);
    if (clouds_3) FMSAddVariable(clouds_3);
    FMSAddVariable(sf_layer);
    FMSAddVariable(pbl);
    FMSAddVariable(atm);

    // set counter to zero
    //count_out = 0;
    // FMS exchange variables
    //diag_ncout->AddVariable( *dT );
    //diag_ncout->AddVariable( *dq );
    //diag_ncout->AddVariable( *dFlux_dT );
    //diag_ncout->AddVariable( *dFlux_dq );
    //diag_ncout->AddVariable( *dt_mass );
    //diag_ncout->AddVariable( *du );
    //diag_ncout->AddVariable( *dv );
    //diag_ncout->AddVariable( *Z_bot );

    FMSAddVariable(*dT);
    FMSAddVariable(*dq);
    FMSAddVariable(*dFlux_dT);
    FMSAddVariable(*dFlux_dq);
    FMSAddVariable(*dt_mass);
    FMSAddVariable(*du);
    FMSAddVariable(*dv);
    FMSAddVariable(*Z_bot);

#ifdef DEBUG_AEOLUS_STOCKS
    stockmove_down_heat_toa = new ScalarVariable ( *grid, "stockmove_down_heat_toa", "J", "net heat stock moved into Aeolus at TOA (irradiation - upward SWR - upward LWR)");
    stockmove_down_heat_out = new ScalarVariable ( *grid, "stockmove_down_heat_out",  "J", "Heat stock moved out of Aeolus in Downward sweep");
    //stockmove_down_water_out= new ScalarVariable ( *grid, "stockmove_down_water", "kg", "Water stock moved down from Aeolus");
    stockmove_up_heat_in = new ScalarVariable ( *grid, "stockmove_up_heat_in",  "J", "Heat stock moved into Aeolus in Upward sweep");
    stockmove_up_water_in= new ScalarVariable ( *grid, "stockmove_up_water_in", "kg", "Water stock moved into Aeolus in Upward sweep");
    stockmove_up_heat_out = new ScalarVariable ( *grid, "stockmove_up_heat_out",  "J", "Heat stock moved out of Aeolus in Upward sweep");
    stockmove_up_water_out= new ScalarVariable ( *grid, "stockmove_up_water_out", "kg", "Water stock moved out of Aeolus in Upward sweep");
    FMSAddVariable(*stockmove_down_heat_toa);
    FMSAddVariable(*stockmove_down_heat_out);
    //FMSAddVariable(*stockmove_down_water);
    FMSAddVariable(*stockmove_up_heat_in);
    FMSAddVariable(*stockmove_up_water_in);
    FMSAddVariable(*stockmove_up_heat_out);
    FMSAddVariable(*stockmove_up_water_out);
#endif


#ifdef TRACE_AEOLUS_UPDATE
    Trace_Aeolus_Update_Down_In = new NetCDFOutput(*grid, "aeolus_trace_down_in.nc", false, 
                                                   "instant trace of input parameters of aeolus_atmosphere_down_()",
                                                   days_since, "NOLEAP");
    // Input to aeolus_atmosphere_down_()
    //u_star_down = new HorizontalVariable(*grid, SF, "u_star", "unit", "friction velocity (in)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE); 
    //Trace_Aeolus_Update_Down_In->AddVariable(*u_star_down);
    //b_star_down = new HorizontalVariable(*grid, SF, "b_star", "unit", "bouyancy scale (in)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    //Trace_Aeolus_Update_Down_In->AddVariable(*b_star_down);
    //q_star_down = new HorizontalVariable(*grid, SF, "q_star", "unit", "moisture scale (in)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    //Trace_Aeolus_Update_Down_In->AddVariable(*q_star_down);
    t_down = new HorizontalVariable(*grid, SF, "t", "K", "surface temperature (in)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    Trace_Aeolus_Update_Down_In->AddVariable(*t_down);
    //dtaudu_down = new HorizontalVariable(*grid, SF, "dtaudu", "unit", "derivative of zonal wind stress w.r.t. the lowest zonal level wind speed (in)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    //Trace_Aeolus_Update_Down_In->AddVariable(*dtaudu_down);
    //dtaudv_down = new HorizontalVariable(*grid, SF, "dtaudv", "unit", "derivative of meridional wind stress w.r.t. the lowest zonal level wind speed (in)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    //Trace_Aeolus_Update_Down_In->AddVariable(*dtaudv_down);
    //u_flux_in_down = new HorizontalVariable(*grid, SF, "u_flux_in", "pa", "zonal wind stress (in)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    //Trace_Aeolus_Update_Down_In->AddVariable(*u_flux_in_down);
    //v_flux_in_down = new HorizontalVariable(*grid, SF, "v_flux_in", "pa", "meridional wind stress (in)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    //Trace_Aeolus_Update_Down_In->AddVariable(*v_flux_in_down);
    // will be modified in LWRM(). write out the input values here.
    Trace_Aeolus_Update_Down_In->AddVariable(radiation->lwr_up_sf_());

    Trace_Aeolus_Update_Down_In->AddVariable(stockflux_down_heat_in, "stockflux_down_heat_in", "J", "Heat stock moved into Aeolus in Downward sweep");

    // add internal state of Aeolus to get more insights
    Trace_Aeolus_Update_Down_In->AddVariable(temperature->T_()); // 3-D temperature
    Trace_Aeolus_Update_Down_In->AddVariable(temperature->TSa_()); // "virtual" surface temperature
    Trace_Aeolus_Update_Down_In->AddVariable(temperature->realTSa_()); // "virtual" surface temperature with real topo
    Trace_Aeolus_Update_Down_In->AddVariable(temperature->QT_()); // heat content atmospheric column
    Trace_Aeolus_Update_Down_In->AddVariable(temperature->dTdx_()); // 
    Trace_Aeolus_Update_Down_In->AddVariable(temperature->dTdy_()); // 
    Trace_Aeolus_Update_Down_In->AddVariable(temperature->LR_());
    Trace_Aeolus_Update_Down_In->AddVariable(pbl->TS_());
    Trace_Aeolus_Update_Down_In->AddVariable(humidity->q_()); // 3-D humidity
    Trace_Aeolus_Update_Down_In->AddVariable(humidity->qS_()); // Surface humidity
    Trace_Aeolus_Update_Down_In->AddVariable(humidity->Qq_()); // Water content per column
    Trace_Aeolus_Update_Down_In->AddVariable(humidity->dqdx_()); //
    Trace_Aeolus_Update_Down_In->AddVariable(humidity->dqdy_()); // 
    
    Trace_Aeolus_Update_Down_Out = new NetCDFOutput(*grid, "aeolus_trace_down_out.nc", false, 
                                                    "instant trace of output parameters of aeolus_atmosphere_down_()",
                                                    days_since, "NOLEAP");
    // Output from aeolus_atmosphere_down_()
    //u_flux_out_down = new HorizontalVariable(*grid, SF, "u_flux_out", "pa", "zonal wind stress (out)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    //Trace_Aeolus_Update_Down_Out->AddVariable(*u_flux_out_down);
    //v_flux_out_down = new HorizontalVariable(*grid, SF, "v_flux_out", "pa", "meridional wind stress (out)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    //Trace_Aeolus_Update_Down_Out->AddVariable(*v_flux_out_down);
    // gust is kept constant at 1 in Aeolus
    //gust_down = new HorizontalVariable(*grid, SF, "gust", "1", "gustiness factor (downward) (out)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    // coszen is input to Aeolus, but output from F90 side of the atmosphere to FMS, thus listed here as output
    coszen_down = new MeridionalVariable1D(*grid, SF, "cosz", "", "cosine of the zenith angle (out)");
    Trace_Aeolus_Update_Down_Out->AddVariable(*coszen_down);
    dtmass_down = dt_mass; // remains const however throughout simulation...
    Trace_Aeolus_Update_Down_Out->AddVariable(*dtmass_down);
    dflux_t_down = dFlux_dT; // remains const however throughout simulation...
    Trace_Aeolus_Update_Down_Out->AddVariable(*dflux_t_down);
    dflux_sphum_down = dFlux_dq; // remains const however throughout simulation...
    Trace_Aeolus_Update_Down_Out->AddVariable(*dflux_sphum_down);
    delta_t_down = dT;
    Trace_Aeolus_Update_Down_Out->AddVariable(*delta_t_down);
    delta_sphum_down = dq;
    Trace_Aeolus_Update_Down_Out->AddVariable(*delta_sphum_down);
    delta_u_down = du; // remains const however throughout simulation...
    Trace_Aeolus_Update_Down_Out->AddVariable(*delta_u_down);
    delta_v_down = dv;  // remains const however throughout simulation...
    Trace_Aeolus_Update_Down_Out->AddVariable(*delta_v_down);
    // create radiative flux HorizontalVariables
    // fluxes: NOTE: "_down" in names indicate that these are calculated in downward sweep (they are not downward fluxes!)
//    flux_sw_down = new HorizontalVariable(*grid, SF, "flux_sw", "W/m2", "Shortwave radiation flux at surface (out)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    flux_sw_sf_down = &(radiation->swr_flux_sf_());
    flux_sw_toa_down = &(radiation->swr_flux_toa_());
    flux_sw_dir_down = new HorizontalVariable(*grid, SF, "flux_sw_dir", "W/m2", "Direct shortwave radiation flux at surface (out)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    flux_sw_dif_down = new HorizontalVariable(*grid, SF, "flux_sw_dif", "W/m2", "Diffusive shortwave radiation flux at surface (out)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    // downward component
    flux_sw_down_vis_dir_down = new HorizontalVariable(*grid, SF, "flux_sw_down_vis_dir", "W/m2", "Downward direct shortwave radiation at surface - visible (out)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    flux_sw_down_vis_dif_down = new HorizontalVariable(*grid, SF, "flux_sw_down_vis_dif", "W/m2", "Downward diffusive shortwave radiation at surface - visible (out)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    flux_sw_down_total_dir_down = new HorizontalVariable(*grid, SF, "flux_sw_down_total_dir", "W/m2", "Downward direct shortwave radiation at surface - total (out)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    flux_sw_down_total_dif_down = new HorizontalVariable(*grid, SF, "flux_sw_down_total_dif", "W/m2", "Downward diffusive shortwave radiation at surface - total (out)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    // visible fluxes
    flux_sw_vis_down = new HorizontalVariable(*grid, SF, "flux_sw_vis", "W/m2", "Shortwave radiation flux at surface - visible (out)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    flux_sw_vis_dir_down = new HorizontalVariable(*grid, SF, "flux_sw_vis_dir", "W/m2", "Shortwave direct radiation flux at surface - visible (out)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    flux_sw_vis_dif_down = new HorizontalVariable(*grid, SF, "flux_sw_vis_dif", "W/m2", "Shortwave diffusive radiation flux at surface - visible (out)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    // longwave fluxes at SF
    flux_lw_sf_down = new HorizontalVariable(*grid, SF, "flux_lw", "W/m2", "Longwave radiation flux at surface (out)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    flux_lw_toa_down = &(radiation->lwr_up_toa_());

    // add to the netcdf file handler
    Trace_Aeolus_Update_Down_Out->AddVariable(*flux_sw_sf_down);
    Trace_Aeolus_Update_Down_Out->AddVariable(*flux_sw_toa_down);
    Trace_Aeolus_Update_Down_Out->AddVariable(*flux_sw_dir_down);
    Trace_Aeolus_Update_Down_Out->AddVariable(*flux_sw_dif_down);
    Trace_Aeolus_Update_Down_Out->AddVariable(*flux_sw_down_vis_dir_down);
    Trace_Aeolus_Update_Down_Out->AddVariable(*flux_sw_down_vis_dif_down);
    Trace_Aeolus_Update_Down_Out->AddVariable(*flux_sw_down_total_dir_down);
    Trace_Aeolus_Update_Down_Out->AddVariable(*flux_sw_down_total_dif_down);
    Trace_Aeolus_Update_Down_Out->AddVariable(*flux_sw_vis_down);
    Trace_Aeolus_Update_Down_Out->AddVariable(*flux_sw_vis_dir_down);
    Trace_Aeolus_Update_Down_Out->AddVariable(*flux_sw_vis_dif_down);
    Trace_Aeolus_Update_Down_Out->AddVariable(*flux_lw_sf_down);
    Trace_Aeolus_Update_Down_Out->AddVariable(*flux_lw_toa_down);
    Trace_Aeolus_Update_Down_Out->AddVariable(radiation->lwr_up_sf_());

    Trace_Aeolus_Update_Down_Out->AddVariable(stockflux_down_heat_out, "stockflux_down_heat_out", "J", "Heat stock moved out of Aeolus in Downward sweep");

    // add internal state of Aeolus to get more insights
    Trace_Aeolus_Update_Down_Out->AddVariable(temperature->T_()); // 3-D temperature
    Trace_Aeolus_Update_Down_Out->AddVariable(temperature->TSa_()); // "virtual" surface temperature
    Trace_Aeolus_Update_Down_Out->AddVariable(temperature->realTSa_()); // "virtual" surface temperature with real topo
    Trace_Aeolus_Update_Down_Out->AddVariable(temperature->QT_()); // heat content atmospheric column
    Trace_Aeolus_Update_Down_Out->AddVariable(temperature->dTdx_()); // 
    Trace_Aeolus_Update_Down_Out->AddVariable(temperature->dTdy_()); // 
    Trace_Aeolus_Update_Down_Out->AddVariable(temperature->LR_());
    Trace_Aeolus_Update_Down_Out->AddVariable(pbl->TS_());
    Trace_Aeolus_Update_Down_Out->AddVariable(humidity->q_()); // 3-D humidity
    Trace_Aeolus_Update_Down_Out->AddVariable(humidity->qS_()); // Surface humidity
    Trace_Aeolus_Update_Down_Out->AddVariable(humidity->Qq_()); // Water content per column
    Trace_Aeolus_Update_Down_Out->AddVariable(humidity->dqdx_()); // Surface humidity
    Trace_Aeolus_Update_Down_Out->AddVariable(humidity->dqdy_()); // Surface humidity

    Trace_Aeolus_Update_Down_Out->AddVariable(temperature->large_transport_x_());
    Trace_Aeolus_Update_Down_Out->AddVariable(temperature->large_transport_y_());
    Trace_Aeolus_Update_Down_Out->AddVariable(temperature->synop_transport_x_());
    Trace_Aeolus_Update_Down_Out->AddVariable(temperature->synop_transport_y_());
    Trace_Aeolus_Update_Down_Out->AddVariable(temperature->meso_transport_x_());
    Trace_Aeolus_Update_Down_Out->AddVariable(temperature->meso_transport_y_());
    Trace_Aeolus_Update_Down_Out->AddVariable(humidity->large_transport_x_());
    Trace_Aeolus_Update_Down_Out->AddVariable(humidity->large_transport_y_());
    Trace_Aeolus_Update_Down_Out->AddVariable(humidity->synop_transport_x_());
    Trace_Aeolus_Update_Down_Out->AddVariable(humidity->synop_transport_y_());
    Trace_Aeolus_Update_Down_Out->AddVariable(humidity->meso_transport_x_());
    Trace_Aeolus_Update_Down_Out->AddVariable(humidity->meso_transport_y_());


    Trace_Aeolus_Update_Up_In = new NetCDFOutput(*grid, "aeolus_trace_up_in.nc", false,
                                                 "instant trace of input parameters of aeolus_atmosphere_up_()",
                                                 days_since, "NOLEAP");
    // Input to aeolus_atmosphere_up_()
    //v_flux_up = NULL;
    //u_star_up = new HorizontalVariable(*grid, SF, "u_star", "unit", "friction velocity (in)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    //b_star_up = new HorizontalVariable(*grid, SF, "b_star", "unit", "bouyancy scale (in)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    //q_star_up = new HorizontalVariable(*grid, SF, "q_star", "unit", "moisture scale (in)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    dt_t_up = new HorizontalVariable(*grid, SF, "dt_t", "K", "increment in temperature in the lowest atmospheric layer (in)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    Trace_Aeolus_Update_Up_In->AddVariable(*dt_t_up);
    dt_sphum_up = new HorizontalVariable(*grid, SF, "dt_sphum", "kg/kg", "increment in specific humidity (tracer) in the lowest atmospheric layer (in)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    Trace_Aeolus_Update_Up_In->AddVariable(*dt_sphum_up);

    Trace_Aeolus_Update_Up_In->AddVariable(stockflux_up_heat_in, "stockflux_up_heat_in", "J", "Heat stock moved into Aeolus in Upward sweep");
    Trace_Aeolus_Update_Up_In->AddVariable(stockflux_up_water_in, "stockflux_up_water_in", "kg", "Water stock moved into Aeolus in Upward sweep");

    // add internal state of Aeolus to get more insights
    Trace_Aeolus_Update_Up_In->AddVariable(humidity->q_()); // 3-D humidity
    Trace_Aeolus_Update_Up_In->AddVariable(humidity->qS_()); // Surface humidity

    Trace_Aeolus_Update_Up_Out = new NetCDFOutput(*grid, "aeolus_trace_up_out.nc", false,
                                                  "instant trace of output parameters of aeolus_atmosphere_up_()",
                                                  days_since, "NOLEAP");
    // Output from aeolus_atmosphere_up_()
    //lprec_up = new HorizontalVariable(*grid, SF, "lprec", "kg/m2/s", "mass of liquid precipitation since last time step (out)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    Trace_Aeolus_Update_Up_Out->AddVariable(humidity->lprec_());
    //fprec_up = new HorizontalVariable(*grid, SF, "fprec", "kg/m2/s", "mass of frozen precipitation since last time step (out)", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_REPLICATE);
    Trace_Aeolus_Update_Up_Out->AddVariable(humidity->fprec_());
    Trace_Aeolus_Update_Up_Out->AddVariable(humidity->precip_());
    //Trace_Aeolus_Update_Up_Out->AddVariable(humidity->preciptest_()); // to check writing via FMS diag_manager
    //gust_up = NULL; gust is kept constant at 1 in Aeolus
    // actually returned in aeolus_get_bottom_mass_()
    t_bot_up = &(temperature->realTSa_());
    Trace_Aeolus_Update_Up_Out->AddVariable(*t_bot_up);
    sphum_bot_up = &(humidity->qS_());
    Trace_Aeolus_Update_Up_Out->AddVariable(*sphum_bot_up);
    p_bot_up = &(wind->sfp_());
    Trace_Aeolus_Update_Up_Out->AddVariable(*p_bot_up);
    z_bot_up = Z_bot; // remains constant ...
    Trace_Aeolus_Update_Up_Out->AddVariable(*z_bot_up);
    p_surf_up = &(wind->sfp_());
    Trace_Aeolus_Update_Up_Out->AddVariable(*p_surf_up);
    slp_up = &(wind->slp_());
    Trace_Aeolus_Update_Up_Out->AddVariable(*slp_up);
    // actually returned in get_bottom_wind()
    u_bot_up = &(pbl->us_());
    Trace_Aeolus_Update_Up_Out->AddVariable(*u_bot_up);
    v_bot_up = &(pbl->vs_());
    Trace_Aeolus_Update_Up_Out->AddVariable(*v_bot_up);

    Trace_Aeolus_Update_Up_Out->AddVariable(stockflux_up_heat_out, "stockflux_up_heat_out", "J", "Heat stock moved out of Aeolus in Upward sweep");
    Trace_Aeolus_Update_Up_Out->AddVariable(stockflux_up_water_out, "stockflux_up_water_out", "kg", "Water stock moved out of Aeolus in Upward sweep");

    // add internal state of Aeolus to get more insights
    Trace_Aeolus_Update_Up_Out->AddVariable(humidity->q_()); // 3-D humidity
    Trace_Aeolus_Update_Up_Out->AddVariable(humidity->qS_()); // Surface humidity

#ifdef TRACE_AEOLUS_INIT
    /* write state at end of init with timestamp 0 */
    Trace_Aeolus_Update_Down_In->OutputToFile(0);
    Trace_Aeolus_Update_Down_Out->OutputToFile(0);
    Trace_Aeolus_Update_Up_In->OutputToFile(0);
    Trace_Aeolus_Update_Up_Out->OutputToFile(0);
#endif

#endif /* TRACE_AEOLUS_UPDATE */


#ifdef WITH_SIMENV
    simenv_output = new SimenvOutput(*grid);
    simenv_output->AddVariable(pbl->TS_());
    /// Add more Variables
#endif
#ifdef WITH_COSTFUNCTION
    init_Aeolus_Score();
#endif

    *nr_Aeolus_diag_vars = FMS_registered_vars.size();

#ifdef TRACE_AEOLUS_INIT
    //delete Trace_Aeolus_Init; Trace_Aeolus_Init = NULL;
#endif
#ifdef USE_CLOCK
    clock_Aeolus.stop();
#endif
    fmsclock_end_(&clock_id_aeolus);
  } // end aeolus_init_()


  void aeolus_finish_(const int *const Time_seconds, const int *const Time_days) {

    double days = (double)(*Time_days) + (double)(*Time_seconds) / (double)(24*60*60);

#ifdef TRACE_AEOLUS_INIT
    if (Trace_Aeolus_Init) { delete Trace_Aeolus_Init; Trace_Aeolus_Init = NULL; }
#endif

#ifdef WITH_COSTFUNCTION
    Aeolus_Score_Calc();
    Aeolus_Score_End(days);
#endif
#ifdef WITH_SIMENV
    simenv_end_c();
#endif

    // save last diagnostic output
    // diag_ncout.OutputToFile(dt_atm_sec);

    // clean up ...
    // save state into restart file
    restart_do_save(days);

    // deallocate atmosphere variables
#define AWAY(_ptr) delete(_ptr); (_ptr) = NULL

    // preferably do it in reverse order of allocation
    FMS_registered_vars.resize(0);
    FMS_open_for_registration = true;

#ifdef TRACE_AEOLUS_UPDATE
    AWAY(Trace_Aeolus_Update_Down_In);
    AWAY(Trace_Aeolus_Update_Down_Out);
    AWAY(Trace_Aeolus_Update_Up_In);
    AWAY(Trace_Aeolus_Update_Up_Out);
#endif
    //AWAY(diag_ncout);

    AWAY(ENERGY_CYCLE);
    AWAY(VAPOR_CYCLE);
    printf("deleting RADIATION\n");
    AWAY(RADIATION);
    printf("deleted RADIATION\n");
    if (CLOUDS_3) AWAY(CLOUDS_3);
    if (CLOUDS_2) AWAY(CLOUDS_2);
    AWAY(WIND);
    AWAY(SYNOPTIC);
    AWAY(PBL_PHYSICS);
    AWAY(PLANWAVE);

    AWAY(atm);
    //AWAY(era40);
    AWAY(pbl);
    AWAY(sf_layer);
    if (clouds_3) AWAY(clouds_3);
    //if (clouds_2) AWAY(clouds_2);
    AWAY(standing_wave);
    AWAY(synoptic);
    AWAY(wind);
    AWAY(radiation);
    AWAY(humidity);
    AWAY(temperature);
    AWAY(topo);

    AWAY(dT);
    AWAY(dFlux_dT);
    AWAY(dq);
    AWAY(dFlux_dq);
    AWAY(dt_mass);
    AWAY(du);
    AWAY(dv);
    AWAY(Z_bot);

    AWAY(INTERP_gust);
    AWAY(INTERP_flux_sw);
    AWAY(INTERP_flux_sw_dir);
    AWAY(INTERP_flux_sw_dif);
    AWAY(INTERP_flux_sw_vis);
    AWAY(INTERP_flux_sw_vis_dir);
    AWAY(INTERP_flux_sw_vis_dif);
    AWAY(INTERP_flux_sw_down_vis_dir);
    AWAY(INTERP_flux_sw_down_vis_dif);
    AWAY(INTERP_flux_sw_down_total_dir);
    AWAY(INTERP_flux_sw_down_total_dif);
    AWAY(INTERP_flux_lw);
    AWAY(INTERP_surf_diff_delta_t);
    AWAY(INTERP_surf_diff_dflux_t);
    AWAY(INTERP_surf_diff_delta_tr);
    AWAY(INTERP_surf_diff_dflux_tr );
    AWAY(INTERP_surf_diff_dtmass);
    AWAY(INTERP_surf_diff_delta_u);
    AWAY(INTERP_surf_diff_delta_v);
    //AWAY(INTERP_surf_diff_delta_co2);
    //AWAY(INTERP_surf_diff_dflux_co2);
    AWAY(INTERP_lprec);
    AWAY(INTERP_fprec);
    AWAY(INTERP_z_bot);
    AWAY(INTERP_t_bot);
    AWAY(INTERP_tr_bot);
    AWAY(INTERP_p_bot);
    AWAY(INTERP_p_surf);
    AWAY(INTERP_slp);
    AWAY(INTERP_u_bot);
    AWAY(INTERP_v_bot);

#ifdef ONLY_TRANSPORT_ON_REDUCED_GRID
    AWAY(reducedgrid);
#endif
    AWAY(grid);

#if 0 && defined(USE_CLOCK)
    // they are printed in implicit destructor invocation after exit()
    clock_Aeolus.print();
    clock_Aeolus_down.print();
    clock_Aeolus_up.print();
#endif
  }


  void aeolus_restart_(const int *const Time_seconds, const int *const Time_days) {
    double days = (double)(*Time_days) + (double)(*Time_seconds) / (double)(24*60);
    restart_do_save(days);
  }




  /** void aeolus_get_axes_coords_
   *
   * Return longitudes and latitudes of cell boundaries and cell midpoints
   * in degrees, and of level heights in m to be used as axes descriptions in the FMS diagnostics code.
   * The parameters are 1-D arrays, defined on the global domain.
   * The returned values always describe a regular grid, regardless if the
   * atmosphere actually uses a reduced grid or not.
   * With reduced grids, it is hard to find the correct southward interface
   * latitude of a cell. Thus we hardcode here the assumption that the grid is
   * regularly spaced in latitudes.
   */
  void aeolus_get_axes_coords_(double *const lonb, double *const latb,
                                 double *const lon, double *const lat,
                                 double *const level, double *const levelb,
                                 int *const nlons, int *const nlats, int *const nlevels) {
    //printf("%s\n", __FUNCTION__);
    //for (int i = 0; i < *nlons; i++) printf(" lon[%d] %g lonb[%d] %g\n", i, lon[i], i, lonb[i]);
    //printf("                lonb[%d] %g\n", *nlons, lonb[*nlons]);
    //for (int i = 0; i < *nlats; i++) printf(" lat[%d] %g latb[%d] %g\n", i, lat[i], i, latb[i]);
    //printf("                latb[%d] %g\n", *nlats, latb[*nlats]);
    //for (int i = 0; i < *nlevels; i++) printf(" level[%d] %g levelb[%d] %g\n", i, level[i], i, levelb[i]);
    //printf("                levelb[%d] %g\n", *nlevels, levelb[*nlevels]);
    //fflush(stdout);
    grid->get_axes_coords(lonb, latb, lon, lat, level, levelb, *nlons, *nlats, *nlevels, /*global*/ true);
#if 0
    const int meridionalcells = grid->MeridionalCells();
    const int zonalcells = grid->ZonalCells();
    const int j_equator = meridionalcells/2;
    const double dlat2 = grid->MeridionalResolution()/2;

    //printf("%s %s\n", __FILE__, __FUNCTION__);
    // for j and i we use the Fortran index range, for code-reading-compatibility
    // with the Fortran side of FMS. When accessing C/C++ arrays,
    // we write index-1.

    assert(nlons == zonalcells);
    for (int i = 1; i <= zonalcells; i++) {
        const Cell *const cell = grid->GetCell(i-1, j_equator, 0);
        lonb[i-1] = cell->FirstConnectedXInterface()->lon();
        lon [i-1] = cell->lon();
    }
    double lonbmax = grid->GetCell(zonalcells-1, j_equator, 0)->SecondConnectedXInterface()->lon();
    lonb[zonalcells] = (epsilon > lonbmax) ? 360 : lonbmax;

    assert(nlats == meridionalcells);
    for (int j = 1; j <= meridionalcells; j++) {
      const Cell *const cell = grid->GetCell(0, j-1, 0);
      const double l = cell->lat();
      lat [j-1] = l;
      latb[j-1] = l - dlat2;
    }
    latb[meridionalcells] = grid->GetCell(0, meridionalcells-1, 0)->lat() + dlat2;

    unsigned int k;
    Cell *c;
    assert(nlevels == grid->VerticalCells());
    for (k=0; k != grid->VerticalCells(); k++) {
      c = grid->GetCell(0, 0, k);
      level[k] = c->z();
      levelb[k] = c->FirstConnectedZInterface()->z();
    }
    levelb[k] = c->SecondConnectedZInterface()->z();
#endif
  }


  /** void aeolus_get_boundary_(double *lon_boundaries, double *lat_boundaries, const int *local)
   *
   * Return the grid box boundaries.
   * lon_boundaries: The west-to-east longitude edges of grid boxes (in radians).
   *   The example data sets provided by GFDL all use grid specifications
   *   with longitudes 0..360.
   * lat_boundaries: The south-to-north latitude edges of grid boxes (in radians).
   * local: a flag that specifies wether boundaries for local or global domain
   *   should be returned.
   * {lon,lat}_{u,l}_{1,2} are the lower and upper bounds of the first and second
   * dimensions of the output array parameters.
   * 
   * The output parameters lon_boundaries and lat_boundaries are
   * 2D-arrays.
   * The dimensions of the output arrays are one more than the number of cells,
   * in particluar (NoNx+1)*(NoNy+1) for global domain and (ie-is+2)*(je-js+2) for local domain.
   * Since we use rectilinear grid, the actual information
   * is 1D only, and thus filled in redundantly. However, if this rectilinearity
   * should change one day, we are already prepared here.
   *
   * We simply use the same method as aeolus_get_axes_coords_() to compute
   * 1-D boundary info, and replicate that into all rows/columns.
   *
   * The returned values always describe a regular grid, regardless if the
   * atmosphere actually uses a reduced grid or not.
   */
  void aeolus_get_boundary_(double *const lon_boundaries, double *const lat_boundaries, const int *const local,
			      const int *const lon_l_1, const int *const lon_u_1,
			      const int *const lon_l_2, const int *const lon_u_2,
			      const int *const lat_l_1, const int *const lat_u_1,
			      const int *const lat_l_2, const int *const lat_u_2) {

    assert(NULL != grid);

    assert(*lon_l_1 == *lat_l_1);
    assert(*lon_l_2 == *lat_l_2);
    assert(*lon_u_1 == *lat_u_1);
    assert(*lon_u_2 == *lat_u_2);

    const int meridionalcellsglobal = grid->MeridionalCellsGlobal();
    const int zonalcells = grid->ZonalCells();
    const double dphi = grid->dPHI();
    //const double dphi2 = dphi/2.0;
    const double dlambda = grid->dLABDA();
    //const double dlambda2 = dlambda/2.0;

    const int ilow  = *local ? is : 1;
    const int ihigh = (*local ? ie : zonalcells)+1; // NoNx+1
    const int jlow  = *local ? js : 1;
    const int jhigh = (*local ? je : meridionalcellsglobal)+1;      // NoNy+1
    const int nlons = ihigh-ilow+1;
    const int nlats = jhigh-jlow+1;

    printf("%s %s output array bounds: (%d:%d, %d:%d) (%d:%d, %d:%d) local %d ilow %d ihigh %d jlow %d jhigh %d, nlons %d, nlats %d\n",
	   __FILE__, __FUNCTION__,
	   *lon_l_1, *lon_u_1, *lon_l_2, *lon_u_2, *lat_l_1, *lat_u_1, *lat_l_2, *lat_u_2,
           *local,
	   ilow, ihigh, jlow, jhigh, nlons, nlats);
    //assert(*lon_l_1 == ilow);
    //assert(*lon_l_2 == jlow);
    //assert(*lon_u_1 == ihigh);
    //assert(*lon_u_2 == jhigh);

    // for j and i we use the Fortran index range, for code-reading-compatibility
    // with the Fortran side of FMS. When accessing C/C++ arrays,
    // we write index-1.
    for (int i = ilow; i <= ihigh; i++) {
      const double phi = (i-1)*dphi;
      for (int j = jlow; j <= jhigh; j++) {
        //printf("i %d j %d index %d phi %g %g\n", i, j, (j-jlow)*nlons+(i-ilow), phi, phi*180.0/M_PI);
        lon_boundaries[(j-jlow)*nlons+(i-ilow)] = phi;
      }
    }
    for (int j  = jlow; j <= jhigh; j++) {
      const double lambda = (j-1)*dlambda-M_PI_2;
      for (int i = ilow; i <= ihigh; i++) {
        //printf("i %d j %d index %d lambda %g %g\n", i, j, (j-jlow)*nlons+(i-ilow), lambda, lambda*180.0/M_PI);
        lat_boundaries[(j-jlow)*nlons+(i-ilow)] = lambda;
      }
    }
    fflush(stdout);
#if 0
    for (int i = ilow; i <= ihigh; i++) {
        const Cell *const cell = grid->GetCell(i-1, j_equator, 0);
        lon_boundaries[i-ilow] = cell->FirstConnectedXInterface()->labda();
	//printf("lon_boundaries[%d] = %f = %f\n",
	//       i-ilow,
	//       cell->FirstConnectedXInterface()->lon(),
	//       cell->FirstConnectedXInterface()->labda());
    }
    double lonbmax = grid->GetCell(ihigh-1, j_equator, 0)->SecondConnectedXInterface()->labda();
    lon_boundaries[ihigh-ilow+1] = (epsilon > lonbmax) ? 2*M_PI : lonbmax;
    //printf("lon_boundaries[%d] = %f = %f\n",
    //	   ihigh-ilow+1,
    //	   grid->GetCell(ihigh-1, j_equator, 0)->SecondConnectedXInterface()->lon(),
    //	   lon_boundaries[ihigh-ilow+1]);

    // now replicate this 1-D result into all other rows
    for (int j = jlow+1; j <= jhigh+1; j++) {
      for (int i = ilow; i <= ihigh+1; i++) {
	lon_boundaries[(j-jlow) * (ihigh+1-ilow+1) + (i-ilow)] = lon_boundaries[i-ilow];
	//printf("lon_boundaries[%d,%d = %d] = l_b[%d] = %f\n",
	//       j, i, (j-jlow) * (ihigh+1-ilow+1) + (i-ilow),
	//       i-ilow,
	//       lon_boundaries[i-ilow]);
      }
    }

    for (int j = jlow; j <= jhigh; j++) {
      const Cell *const cell = grid->GetCell(0, j-1, 0);
      const double l = cell->phi();
      lat_boundaries[(j-jlow) * (ihigh+1-ilow+1)] = l - dphi2;
      //printf("lat_boundaries[%d] = %f = %f\n",
      //	     (j-jlow) * (ihigh+1-ilow+1),
      //	     cell->lat(),
      //	     l - dphi2);
    }
    lat_boundaries[(jhigh-jlow+1) * (ihigh+1-ilow+1)] = grid->GetCell(0, jhigh-1, 0)->phi() + dphi2;
    //printf("lat_boundaries[%d] = %f = %f\n",
    //	   (jhigh-jlow+1) * (ihigh+1-ilow+1),
    //	   grid->GetCell(0, jhigh-1, 0)->lat(),
    //	   grid->GetCell(0, jhigh-1, 0)->phi() + dphi2);

    // now replicate this 1-D result into all other columns
    for (int i = ilow+1; i <= ihigh+1; i++) {
      for (int j = jlow; j <= jhigh+1; j++) {
	lat_boundaries[(j-jlow) * (ihigh+1-ilow+1) + (i-ilow)] = lat_boundaries[(j-jlow) * (ihigh+1-ilow+1)];
      }
    }

    //fflush(stdout);
#endif
  }




  /** aeolus_cell_area_(double *area_out)
   * area_out: 2D array defined on local domain
   */
  void aeolus_cell_area_(double *const area_out) {

    // When using a grid with merged cells towards the pole,
    // we must split up our grid cells into corresponding cells of the
    // regular lat/lon grid that FMS expects. Thus we use for each longitudinal
    // belt the merge_factor, which gives the number of regular cells that were
    // merged into our cells. For regular, non-merged, grids, the merge factor
    // will always be 1.
    // Note that for the merged grid case different i,j indices will
    // yield pointers to the same cell.

    const int zonalcells = grid->ZonalCells(); // at equator

    // for j and i we use the Fortran index range, for code-reading-compatibility
    // with the Fortran side of FMS. When accessing C/C++ arrays,
    // we write index-1.
    //printf("%s is %d ie %d js %d je %d\n", __FUNCTION__, is, ie, js, je);
    //fprintf(stderr, "%s is %d ie %d js %d je %d\n", __FUNCTION__, is, ie, js, je);
    for (int j = js; j <= je; j++) {
      const int merge_factor = zonalcells / grid->ZonalCells(j-js);
      //printf("  j %d merge_factor %d\n", j, merge_factor);
      for (int i = is; i <= ie; i++) {
        //fprintf(stderr, "GetCell(%d, %d, 0)\n", i-is, j-js); fflush(stderr);
	const Cell *cell = grid->GetCell(i-is, j-js, 0);
	const Interface *iface = cell->FirstConnectedZInterface();
	assert(iface->AtBoundary() && (iface->BoundaryFlag() & BOTTOM));
	const double area = iface->Area() / merge_factor;
	//printf("i %3d j %3d index %d area %f\n", i, j, (j-js) * (ie-is+1) + (i-is), area);
	area_out[(j-js) * (ie-is+1) + (i-is)] = area;
      }
    }
  }



  /** enum stock_idx mimics definitions in shared/exchange/stock_constants.F90 */
  enum stock_idx { ISTOCK_WATER = 1, ISTOCK_HEAT = 2, ISTOCK_SALT = 3 };

  // Compute current stock of water/heat/salt present in the local data
  // domain of this processing element.
  // If checking of mass/energy conservation is switched on in the model
  // configuration, this is called once after initialisation, and
  // then after each time step.
  // idx: index specifying which stock to compute
  // val: scalar output parameter
  void aeolus_get_stock_pe_(const int *const idx, double *const val) {
    // Sigh. This function is sometimes called from atmosphere tasks,
    // and sometimes from the coupler for all tasks.
    // I.e. from different pe_set contexts. However, an FMS clock
    // can must be invoked in the pe_set (i.e. Communicator) context
    // in which it was defined.
    //fmsclock_begin_(&clock_id_aeolus_getstocks);
    switch (*idx) {
    case ISTOCK_WATER:
        // calculate total water vapor in atmosphere
        *val = 0.; 
        for(std::vector<Cell>::const_iterator it=grid->SurfaceCellsBegin(); it!=grid->SurfaceCellsEnd(); it++)
            *val += humidity->Qq_().Read(it->ID()) * it->FirstConnectedZInterface()->Area(); // kg/m2 * m2 = kg
        break;
    case ISTOCK_HEAT:
        // calculate total heat content in atmosphere
        *val = 0.; 
        for(std::vector<Cell>::const_iterator it=grid->SurfaceCellsBegin(); it!=grid->SurfaceCellsEnd(); it++)
            *val += temperature->QT_().Read(it->ID()) * it->FirstConnectedZInterface()->Area(); // J/m2 * m2 = J
        break;
    case ISTOCK_SALT:
      *val = 0;   // atmosphere does not have salt stock
      break;
    default:
      *val = 0;   // don't know how to compute any other stock
    }
    //fmsclock_end_(&clock_id_aeolus_getstocks);
  }


  // Compute current stock of water/heat/salt present globally.
  // If checking of mass/energy conservation is switched on in the model
  // idx: index specifying which stock to compute
  // val: scalar output parameter
  void aeolus_get_stock_atm_global_(const int *const idx, double *const val) {
    aeolus_get_stock_pe_(idx, val);
#ifdef USE_MPI_DOMAIN_DECOMPOSE
    double global;
    grid->MPI_Comm().Allreduce(val, &global, 1, MPI::DOUBLE, MPI::SUM);
    *val = global;
#endif
  }

} // extern "C"


#if 0
// Compute amount of stock moved upward into the atmosphere.
// On 16-May-2014, Dim Coumou wrote:
//What Aeolus *receives***is stored in Aeolus DataStorage class SurfaceLayer variables:
//HorizontalVariable   F_sens;           //< sensible heat flux from surface: J/s/m2
//HorizontalVariable     evapo;            //< evaporation / evapotranspiration: kg/m2/s
//To convert this to absolute "stocks" rather than fluxes, multiply with timestep and cell-surface area (using "cell->FirstConnectedZInterface()->Area()").
//  heat_flux = F_sens * cell->FirstConnectedZInterface()->Area() * timestep    // JOULES
//  water_flux = evapo * cell->FirstConnectedZInterface()->Area() * timestep    // kg
//  (So this can be implemented similar to aeolus_get_stock_pe_)
// Note that hose are passed as surf_diff_delta_t and surf_diff_delta_tr[:,:,q_ind] from the coupler
double Aeolus_get_flux_stock_up_pe(const int idx) {
  double val = 0.0;
  switch (idx) {
  case ISTOCK_WATER:
    for(std::vector<Cell>::const_iterator it=grid->SurfaceCellsBegin(); it!=grid->SurfaceCellsEnd(); it++)
      val += sf_layer->evapo_().Read(it->ID()) * it->FirstConnectedZInterface()->Area();
    break;
  case ISTOCK_HEAT:
    for(std::vector<Cell>::const_iterator it=grid->SurfaceCellsBegin(); it!=grid->SurfaceCellsEnd(); it++)
      val += sf_layer->F_sens_().Read(it->ID()) * it->FirstConnectedZInterface()->Area();
    break;
  case ISTOCK_SALT:
    val = 0;   // atmosphere does not receive salt stock
    break;
  default:
    val = 0;   // don't know how to compute any other stock flux
  }
  val *= dt_atm_sec;
  return val;
}
#endif


#if 0
// Compute amount of stock moved downward out of the atmosphere.
// On 16-May-2014, Dim Coumou wrote:
//What Aeolus *gives***is stored in FMS variables:
//- radiative fluxes:
//  double *const flux_sw,               // net shortwave flux (W/m2) at the surface
//  double *const flux_lw,                // DC: downward longwave flux (W/m2) at the surface
//  Again, to convert:
//  heat_flux = (flux_sw + flux_lw) * cell->FirstConnectedZInterface()->Area() * timestep    // JOULES
// flux_sw is taken in aeolus_atmosphere_down_() from radiation->swr_flux_sf_(),
// flux_lw is taken in aeolus_atmosphere_down_() from radiation->lwr_down_sf_(),
//    positive upward in Aeolus, downward in FMS
//- precipitation
//    double *const lprec,  // mass of liquid precipitation since last time step (Kg/m2/s)
//    double *const fprec,  // mass of frozen precipitation since last time step (Kg/m2/s): set zero for now
//    Again, to convert:
//    water_flux = (lprec + fprec) * cell->FirstConnectedZInterface()->Area() * timestep    // kg
// lprec is copied in aeolus_atmosphere_up_() from humidity->precip_()
// fprec is currently always 0.0

double Aeolus_get_flux_stock_down_pe(const int idx) {
  double val = 0.0;
  switch (idx) {
  case ISTOCK_WATER:
    for(std::vector<Cell>::const_iterator it=grid->SurfaceCellsBegin(); it!=grid->SurfaceCellsEnd(); it++)
      val += (humidity->precip_().Read(it->ID()) /* + fprec */) * it->FirstConnectedZInterface()->Area();
    break;
  case ISTOCK_HEAT:
    for(std::vector<Cell>::const_iterator it=grid->SurfaceCellsBegin(); it!=grid->SurfaceCellsEnd(); it++)
      val += (radiation->swr_flux_sf_().Read(it->ID()) - radiation->lwr_down_sf_().Read(it->ID())) * it->FirstConnectedZInterface()->Area();
    break;
  case ISTOCK_SALT:
    val = 0;   // atmosphere does not receive salt stock
    break;
  default:
    val = 0;   // don't know how to compute any other stock flux
  }
  val *= dt_atm_sec;
              return val;
}
#endif

#ifdef DEBUG_AEOLUS_STOCKS
#define DEBUG_STOCK_CHECK(_fun,_place) \
  aeolus_get_stock_pe_(&iwater, &water2); \
  if (water1 != water2) { printf("%s:%s: water stock changed from %22.15g to %22.15g  diff %22.15g\n", (_fun), (_place), water1, water2, water2-water1); water1 = water2; } \
  aeolus_get_stock_pe_(&iheat,  &heat2); \
  if (heat1 != heat2) { printf("%s:%s: heat stock changed from %22.15g to %22.15g  diff %22.15g\n", (_fun), (_place), heat1, heat2, heat2-heat1); heat1 = heat2; }
#else
#define DEBUG_STOCK_CHECK(_fun,_place) {}
#endif


// Some helper functions wich are declared as friend functions to the storage classes.
// For safety and simplicity, we define them as C++ linkage
// They are defined before use to ease inlining and optimisation for the compiler.
void LoadTopographyFromF90(const double *const height, const double *const stddev, bool without_topography, Topography& t)
{
  t.real_height.from_F90(grid, height);
  t.real_height.UpdateHalos(*grid);
  t.real_sigma.from_F90(grid, stddev);
  if (without_topography) {
    t.sigma.UniformValue(0.);
    t.height.UniformValue(0.);
  } else {
    t.height.from_F90(grid, height);
    t.sigma.from_F90(grid, stddev);
  }
  t.height.UpdateHalos(*grid);
  t.real_sigma.PRINTSUM;
  t.real_height.PRINTSUM;
  t.sigma.PRINTSUM;
  t.height.PRINTSUM;
}

void LoadLandfracFromF90( const double *const landfrac, Topography& t)
{
  t.fraction_land.from_F90(grid, landfrac);
  t.fraction_land.PRINTSUM;
  t.fraction_ocean  = t.fraction_land;
  t.fraction_ocean *= -1.;
  t.fraction_ocean += 1.0; // f_ocean = 1.0 - f_land
}

void LoadAlbedosFromF90  (  const double *const albedo_vis_dir, 
                            const double *const albedo_nir_dir,
                            const double *const albedo_vis_dif,
                            const double *const albedo_nir_dif, 
                            const double *const roughness,
                            SurfaceLayer& sf,
                            const double *const cosz_radwt)
{
  // TO DO:   - adapt SWRM interface to handle averaged albedos instead of for each surface type
  //          - adapt SurfaceLayer accordingly: now only storing averaged albedos
  //          - uncomment lines below...

  // this assigns the same albedos to each of the NST surface types
  // at some stage this should be removed: atmosphere & radiation module do not need to know about surface types
  for(unsigned int sf_type=0; sf_type!=sf.SurfaceTypes(); sf_type++) {
    sf.alb_vis_clear[sf_type]->from_F90(grid, albedo_vis_dir);
    sf.alb_vis_cloudy[sf_type]->from_F90(grid, albedo_vis_dif);
    sf.alb_ir_clear[sf_type]->from_F90(grid, albedo_nir_dir);
    sf.alb_ir_cloudy[sf_type]->from_F90(grid, albedo_nir_dif);
  }
  sf.roughness_length.from_F90(grid, roughness);

  sf.alb_vis_clear_(0).PRINTSUM;
  sf.alb_vis_cloudy_(0).PRINTSUM;
  sf.alb_ir_clear_(0).PRINTSUM;
  sf.alb_ir_cloudy_(0).PRINTSUM;
  sf.roughness_length_().PRINTSUM;

  if (ocean_albedo_climber2 && !(use_sinsol) && cosz_radwt) { // here, cosz is defined only when !use_sinsol
    for(std::vector<Cell>::const_iterator it=grid->SurfaceCellsBegin(); it!=grid->SurfaceCellsEnd(); it++) {
      if (sf.fraction_sftype[0]->Read(it->ID()) >= 0.99999) { // ocean-only. Maybe use threshold of 0.5 or .9
        // However, here we cannot distinguish between ice-free ocean and sea-ice-covered ocean.
        const int j = grid->GetIndexFromPHI(it->phi());
        const double COSZi = cosz_radwt[j];
        const double ALB1 = COSZi > 0.0001 ? 0.05/COSZi : 0.05;
        const double AMIN1 = ALB1 < 0.20 ? ALB1 : 0.20;
        sf.alb_vis_clear[0]->Write(it->ID(), AMIN1);
        sf.alb_vis_cloudy[0]->Write(it->ID(), 0.07);
        sf.alb_ir_clear[0]->Write(it->ID(), AMIN1);    // use same values for Near-Infrared albedos
        sf.alb_ir_cloudy[0]->Write(it->ID(), 0.07);
      }
    }
  }
}

void InitSurfaceTypesFromLandOceanMask( const HorizontalVariable& fraction_ocean, const HorizontalVariable& fraction_land, SurfaceLayer& sf )
{
  // initialise to zero
  for (unsigned int sf_type = 0; sf_type != sf.SurfaceTypes(); sf_type++)
    sf.fraction_sftype[sf_type]->UniformValue(0.);
  // set "open water" to fraction ocean and "trees" to fraction land. All others remain 0
  *sf.fraction_sftype[0] = fraction_ocean;
  *sf.fraction_sftype[2] = fraction_land;
}

void LoadSkinTemperatureFromF90( const double *const t_skin, SurfaceLayer& sf) {
    sf.T_skin.from_F90(grid, t_skin);
    sf.T_skin.PRINTSUM;
}
/*
void LoadSurfaceFluxesFromF90( const double *const delta_t, const double *const delta_q, SurfaceLayer& sf)
{
  // heat flux
  sf.F_sens.from_F90(grid, delta_t); // F_sens = dT [K]
  sf.F_sens /= dt_mass_atm;       // [K/s kg/m2]
  sf.F_sens *= Atmosphere.cp(); // J/s/m2
  // vapor flux
  sf.evapo.from_F90(grid, delta_q); // evapo = delta_q [-]
  sf.evapo /= dt_mass_atm;       // [kg/m2/s]
}
*/
void LoadSurfaceFluxesFromF90( const HorizontalVariable& delta_t, const HorizontalVariable& delta_q, SurfaceLayer& sf)
{
  // heat flux
  sf.F_sens.CopyDataFrom(delta_t); // F_sens = dT [K]
  sf.F_sens.PRINTSUM;
  sf.F_sens /= dt_mass_atm;           // [K/s kg/m2]
  sf.F_sens *= Atmosphere.cp();    // J/s/m2
  // vapor flux
  sf.evapo.CopyDataFrom(delta_q);  // evapo = delta_q [-]
  sf.evapo.PRINTSUM;
  sf.evapo /= dt_mass_atm;            // [kg/m2/s]
}

void  LoadFluxLWRSFUp(const double *const flux_lwr_sf_up, class RadiativeFluxes *r) {
  r->lwr_up_sf.from_F90(grid, flux_lwr_sf_up);
}

void override_qS_setup() {
  override_qS_input = new NetCDFInput(*grid, (*override_qS)());
  override_qS_recordcount = override_qS_input->RecordCount();
  if (override_qS_recordcount < 1) {
    fflush(stdout);
    fprintf(stderr, "ERROR: override_qS: file %s must have at least one record\n",
            (*override_qS)());
    fflush(stderr);
    abort();
  }
  // for now, we require daily climatology
  if (override_qS_recordcount != 365) {
    fflush(stdout);
    fprintf(stderr, "ERROR: override_qS: file %s must have 365 timesteps for daily climatology\n",
            (*override_qS)());
    fflush(stderr);
    abort();
  }
  // sigh. qS is private.
  override_qS_input->GetVariable(humidity->qS, "q");
}

void override_ww_setup() {
  override_ww_input = new NetCDFInput(*grid, (*override_ww)());
  override_ww_recordcount = override_ww_input->RecordCount();
  if (override_ww_recordcount < 1) {
    fflush(stdout);
    fprintf(stderr, "ERROR: override_ww: file %s must have at least one record\n",
            (*override_ww)());
    fflush(stderr);
    abort();
  }
  // for now, we require seasonal climatology
  if (override_ww_recordcount != 4) {
    fflush(stdout);
    fprintf(stderr, "ERROR: override_ww: file %s must have 4 timesteps for seasonal climatology\n",
            (*override_ww)());
    fflush(stderr);
    abort();
  }
  // sigh. ww is private.
  override_ww_input->GetVariable(synoptic->ww, "WW");
}

extern "C" {

  void aeolus_atmosphere_down_
  (
   const int *const Time_seconds, const int *const Time_days,
   const int *const dayoftheyear,

   /// most of the following parameters are 2-D arrays,
   /// - except the tracer arrays denoted as 3D arrays -
   /// defined on the local data domain with index range (is:ie,js:je)

   /// input parameters
   /// quantities going from land+ice to atmos
   const double *const land_frac,       ///< fraction amount of land in a grid box
   const double *const t_surf,          ///< surface temperature for radiation calculations
   const double *const albedo,          ///< surface albedo for radiation calculations
   const double *const albedo_vis_dir,
   const double *const albedo_nir_dir,
   const double *const albedo_vis_dif,
   const double *const albedo_nir_dif,
   const double *const rough_mom,       ///< surface roughness (used for momentum)

   /// DC: following block of input parameters are not required: surface winds are directly calculated from roughness_length only.
   /// these are not used
   const double *const u_star,          ///< friction velocity
   const double *const b_star,          ///< bouyancy scale
   const double *const q_star,          ///< moisture scale
   const double *const dtau_du,         ///< derivative of zonal wind stress w.r.t. the lowest zonal level wind speed
   const double *const dtau_dv,         ///< derivative of meridional wind stress w.r.t. the lowest meridional level wind speed
   //const double *const frac_open_sea,   ///< non-seaice fraction
   // inout parameters
   double *const u_flux,                ///< zonal wind stress       DC: output only!
   double *const v_flux,                ///< meridional wind stress  DC: output only!

   // output parameters
   double *const gust,                  ///< gustiness factor
   // cosz and solar were already calculated in the F90-side, and thus are now input-parameters for Aeolus
   const double *const cosz_radwt,      ///< cosine of the zenith angle, per latitude, 1-D array (js:je)
                                        ///< Note that the SWR code in Aelus needs radiation-weighted cosz!
   const double *const solar,           ///< solar irradiation at top of atmosphere, per latitude, 1-D array (js:je)
   // to avoid non-linear interference with regular-to-reduced grid interpolation, we calculate
   // the upgoing long wave radiation at the surface, according to Stefan-Boltzmann-Law,
   // from t_surf already in the Fortran interface, and pass it to here.
   double *const flux_lwr_sf_up,        ///< upgoing lwr at surface, calculated as in surface_flux()

   // output parameters again
   double *const flux_sw,               ///< net shortwave flux (W/m2) at the surface
   double *const flux_sw_dir,
   double *const flux_sw_dif,
   double *const flux_sw_down_vis_dir,
   double *const flux_sw_down_vis_dif,
   double *const flux_sw_down_total_dir,
   double *const flux_sw_down_total_dif,
   double *const flux_sw_vis,
   double *const flux_sw_vis_dir,
   double *const flux_sw_vis_dif,
   double *const flux_lw,               ///< DC: downward longwave flux (W/m2) at the surface. The FMS coupler computes the upgoing lwr flux on its own, and then the net LWR flux at surface
   // output parameters from surf_diff_type fields
   double *const surf_diff_delta_tr, ///< 3D, 3rd dim for tracer id.
                                     /// - similarly for the increment in specific humidity
                                     ///   (non-dimensional  = Kg/Kg)
   double *const surf_diff_dflux_tr, ///< 3D, 3rd dim for tracer id.
                                     /// - derivative of the flux of specific humidity
                                     ///   at the top of the lowest atmospheric layer with respect to
                                     ///   the specific humidity of that layer  (--/(m2 K))
   double *const surf_diff_delta_t,  ///< the increment in temperature in the lowest atmospheric
                                     /// layer (((i+1)-(i-1) if atmos model is leapfrog) (K)
                                     /// (defined in  gcm_vert_diff_down as the increment computed up
                                     /// to this point in model, including effect of vertical
                                     /// diffusive flux at top of lowest model level, presumed
                                     /// to be modified to include effects of surface fluxes
	                             /// outside of this module, then used to start upward
	                             /// tridiagonal sweep,
   double *const surf_diff_dflux_t,  ///< derivative of the temperature flux at the top of the lowest
                                     /// atmospheric layer with respect to the temperature
                                     /// of that layer  (J/(m2 K))
   double *const surf_diff_dtmass,   ///< dt/mass, where dt = atmospheric time step (sec)
                                     /// mass = mass of lowest atmospheric layer (Kg/m2)
   double *const surf_diff_delta_u,  //
   double *const surf_diff_delta_v,  //
   const int *const fortytwo ///< minimal interface check
   ) { // begin aeolus_atmosphere_down_()
    // ----------------------------
    // update Atmospheric variables
    // ----------------------------

#ifdef USE_CLOCK
    clock_Aeolus.start();
    clock_Aeolus_down.start();
#endif
    fmsclock_begin_(&clock_id_aeolus);
    fmsclock_begin_(&clock_id_aeolus_down);
    static bool do_update_land_frac = true;
    assert(42 == *fortytwo);
#ifndef NDEBUG
    static unsigned int n_atmos_down_steps = 0;
    n_atmos_down_steps++;
    if (pe == root_pe) {
      cout << __FUNCTION__ << " timestep nr " << n_atmos_down_steps << " at " << *Time_days << ":" << *Time_seconds << " day:sec, "
           << (*Time_days-first_day)/365 << "Y" << (*Time_days-first_day)%356 << "d" << *Time_seconds << "s after start of run" << endl;
    }
#endif

#ifdef DEBUG_AEOLUS_STOCKS
    double water1, water2, heat1, heat2;
    const int iwater = ISTOCK_WATER, iheat = ISTOCK_HEAT;
    aeolus_get_stock_pe_(&iwater, &water1);
    aeolus_get_stock_pe_(&iheat,  &heat1);
#endif

    // FOR NOW: assume regular atmospheric grid, and hence use F90_to_SurfaceVariable_simple
    // TODO:   make more flexible to handle reduce grids, hint: use F90_to_SurfaceVariable_avg or appropriate other exchange function

    // map FMS input variables on Aeolus variables
    if (do_update_land_frac) {
      LoadLandfracFromF90      ( land_frac, *topo );
      InitSurfaceTypesFromLandOceanMask(topo->fraction_ocean_(), topo->fraction_land_(), *sf_layer);
      if (!update_land_frac_always) do_update_land_frac = false;
    }
    // if ocean_albedo_climber2 is true, the albedo for ice-free ocean is overridden with the equation
    // from Climber-2 SVAT inside LoadAlbedosFromF90()
    LoadAlbedosFromF90         ( albedo_vis_dir, albedo_nir_dir, albedo_vis_dif, albedo_nir_dif, rough_mom, *sf_layer, cosz_radwt );

    LoadSkinTemperatureFromF90 ( t_surf, *sf_layer );      // storing t_surf in sf_layer.T_skin

    LoadFluxLWRSFUp(flux_lwr_sf_up, radiation);

    if ((*override_qS)()) {
      override_qS_input->InputFromFile(override_qS_rec, false);
      override_qS_rec++;
      printf("overriding qS from record %ld\n", override_qS_rec);
      if (override_qS_rec >= override_qS_recordcount) override_qS_rec = 0;
    }

    if ((*override_ww)()) {
      int season = ((*Time_days)%365)/4 /*+ 1*/; // records are numbered from 0
      if (override_ww_rec != season) { // read data on change of season, not on every time step
        override_ww_rec = season;
        override_ww_input->InputFromFile(override_ww_rec, true /* convert pressure levels to meters */);
        printf("overriding ww from record %ld\n", override_ww_rec);
      }
    }

    // Before any calculation is done: 
    // set FMS flux variables to values of previous timestep: T(n), qS(n)
    dT->CopyDataFrom( temperature->realTSa_()); // set equal to 2-meter temperature
    dq->CopyDataFrom( humidity->qS_() );    // set equal to q in bottom layer [kg/kg]

    if (!use_sinsol) {
      FMSAstronomy::setCOSZM(cosz_radwt);
      FMSAstronomy::setSOLARM(solar);
    }

    // output before calculations
    //debug_out_counter = 0;

#ifdef TRACE_AEOLUS_UPDATE
    t_down->from_F90(grid, t_surf);
    if (!use_sinsol) coszen_down->from_F90(cosz_radwt);
    Trace_Aeolus_Update_Down_In->OutputToFile(*Time_days + (*Time_seconds / (24.0 * 60.0 * 60.0)));    
#endif

    // ----------------------------
    // take Atmospheric timestep
    // ----------------------------

    VAPOR_CYCLE->CalculateHe();
    VAPOR_CYCLE->Calculate_rh_sf(); // needed for PBL::CalculatePBLThickness()
    // Dynamics: PBL -> PlanWave -> Synoptics -> LargeScaleField
    PBL_PHYSICS -> Calculate(dt_atm_sec);
    DEBUG_STOCK_CHECK("aeolus_atmosphere_down_", "A")
    PLANWAVE    -> Calculate(dt_atm_sec);
    DEBUG_STOCK_CHECK("aeolus_atmosphere_down_", "B")
    SYNOPTIC    -> Calculate(dt_atm_sec);
    DEBUG_STOCK_CHECK("aeolus_atmosphere_down_", "C")
    WIND        -> Calculate(dt_atm_sec);
    DEBUG_STOCK_CHECK("aeolus_atmosphere_down_", "D")

    // Physics: Clouds -> Radiation ->  WaterVaporTransfer -> HeatTransfer
    if (CLOUDS_2)
      CLOUDS_2  -> Calculate(dt_atm_sec);
    else {
      VAPOR_CYCLE->Calculate_rhum();
      CLOUDS_3  -> Calculate(dt_atm_sec);
    }
    DEBUG_STOCK_CHECK("aeolus_atmosphere_down_", "E")

    RADIATION   -> Calculate(*dayoftheyear);
    DEBUG_STOCK_CHECK("aeolus_atmosphere_down_", "F")
    VAPOR_CYCLE -> CalculateDownwardSweep(dt_atm_sec);    // calculating vapor transport. No precipitation, no surface fluxes
    DEBUG_STOCK_CHECK("aeolus_atmosphere_down_", "G")
    ENERGY_CYCLE-> CalculateDownwardSweep(dt_atm_sec);    // calculating energy transport & latent heating & radiative heating. Excluding surface fluxes
    DEBUG_STOCK_CHECK("aeolus_atmosphere_down_", "H")

    // ----------------------------
    // map return variables back on FMS pointers
    // ----------------------------
    // output: gust, flux_sw, flux_sw_dir, flux_sw_dif, flux_sw_down_vis_dir, flux_sw_down_vis_dif, flux_sw_down_total_dir,
    //         flux_sw_down_total_dif, flux_sw_vis, flux_sw_vis_dir, flux_sw_vis_dif, flux_lw
    //         surf_diff_delta_tr,  surf_diff_dflux_tr, surf_diff_delta_t, surf_diff_dflux_t, surf_diff_dtmass,
    //         surf_diff_delta_u, surf_diff_delta_v

      static HorizontalVariable dummy (*grid, SF, "DUMMY_DOWN", "", "dummy variable used in aeolus_atmosphere_down", Variable::INTERP_OUT_INTERPOL|Variable::INTERP_IN_AVERAGE);


    // wind stress
    (pbl->tau_x_()).to_F90(grid, u_flux);
    (pbl->tau_y_()).to_F90(grid, v_flux);
    
    // gustiness:
    dummy.UniformValue(1.); // (EBM treats gustiness this way)
    dummy.to_F90(grid, gust,  (*INTERP_gust)());

    // SW radiation. Current implementation takes same approach as EBM 
    (radiation->swr_flux_sf_()).to_F90(grid, flux_sw, (*INTERP_flux_sw)());     // net swr flux at surface

#ifdef TRACE_AEOLUS_UPDATE
    //flux_sw_down == radiation->swr_flux_sf_()
#endif

    // copied from atmos_ebm/atmosphere.F90
    //
    //           VIS                NIR          VIS+NIR
    //-----------------------------------------------------
    //       |                 |              |
    // DIR   |flux_sw_vis_dir  |              | flux_sw_dir
    //       |                 |              |
    //-----------------------------------------------------
    //       |                 |              |
    // DIF   |flux_sw_vis_dif  |              | flux_sw_dif
    //-----------------------------------------------------
    //       |                 |              |
    //DIR+DIF|flux_sw_vis      |              | flux_sw        
    //       |                 |              |(diagnostic only)
    //-----------------------------------------------------

    // bls-total Short wave is equally divided into VIS and NIR
    // flux_sw_dif - total Diffusive Short wave : VIS + NIR

    // bls- Visible Direct and Diffusive short wave will be assumed to be
    // half of the total VIS+DIR components

    // NIR components are determined in the ice model

    dummy.CopyDataFrom(radiation->swr_flux_sf_());
    dummy *= 0.5; // assume visible, direct and diffusive components are half of total (EBM handles fluxes this way)
    dummy.to_F90(grid, flux_sw_dir, (*INTERP_flux_sw_dir)()); // net swr direct flux at surface
    dummy.to_F90(grid, flux_sw_dif, (*INTERP_flux_sw_dif)()); // net swr diffusive flux at surface
    dummy.to_F90(grid, flux_sw_vis, (*INTERP_flux_sw_vis)()); // net visible swr flux at surface

#ifdef TRACE_AEOLUS_UPDATE
    flux_sw_dir_down->CopyDataFrom(dummy); ///< TODO: improve naming of trace variables to avoid confusion with radiation directions
    flux_sw_dif_down->CopyDataFrom(dummy);
    flux_sw_vis_down->CopyDataFrom(dummy);
#endif

    dummy *= 0.5; // assume direct and diffusive parts of visible component are half of total visible (EBM handles fluxes this way)
    dummy.to_F90(grid, flux_sw_vis_dir, (*INTERP_flux_sw_vis_dir)()); // net direct visible swr flux at surface
    dummy.to_F90(grid, flux_sw_vis_dif, (*INTERP_flux_sw_vis_dif)()); // net diffusive visible swr flux at surface

#ifdef TRACE_AEOLUS_UPDATE
    flux_sw_vis_dir_down->CopyDataFrom(dummy);
    flux_sw_vis_dif_down->CopyDataFrom(dummy);
#endif

    dummy.UniformValue(0.); // assume all other partial components to be zero (EBM handles fluxes this way)
    dummy.to_F90(grid, flux_sw_down_vis_dir,   (*INTERP_flux_sw_down_vis_dir)());   // net direct visible swr flux at surface
    dummy.to_F90(grid, flux_sw_down_vis_dif,   (*INTERP_flux_sw_down_vis_dif)());   // net diffusive visible swr flux at surface
    dummy.to_F90(grid, flux_sw_down_total_dir, (*INTERP_flux_sw_down_total_dir)()); // net direct total swr flux at surface
    dummy.to_F90(grid, flux_sw_down_total_dif, (*INTERP_flux_sw_down_total_dif)()); // net diffusive total swr flux at surface

#ifdef TRACE_AEOLUS_UPDATE
    flux_sw_down_vis_dir_down->CopyDataFrom(dummy);
    flux_sw_down_vis_dif_down->CopyDataFrom(dummy);
    flux_sw_down_total_dir_down->CopyDataFrom(dummy);
    flux_sw_down_total_dif_down->CopyDataFrom(dummy);
#endif
    
    // LW radiation:
    dummy.CopyDataFrom(radiation->lwr_down_sf_());
    dummy *= -1.;                 // LWR is defined positive upwards in Aeolus
    dummy.to_F90(grid, flux_lw, (*INTERP_flux_lw)());  // downward lwr flux surface
#ifdef TRACE_AEOLUS_UPDATE
    flux_lw_sf_down->CopyDataFrom(dummy);
#endif
    
    // ----------------------------
    // compute FMS flux conserving exchange variables
    // ----------------------------
 
    // calculate the change in near-surface temperature and humidity 
    (*dT) *= -1.;                  // dT = - TSa(n)
    (*dT) += temperature->realTSa_();// dT = TSa(n+1) - TSa(n)
    (*dq) *= -1.;                  // dq = - qS(n)
    (*dq) += humidity->qS_();      // dq = qS(n+1) - qS(n)

    // dt_mass remains constant with time
    // flux-derivatives are zero
    // du / dv = u / v stress on atmosphere, are set to zero

    // update the FMS pointers
    dT      ->to_F90(grid, surf_diff_delta_t, (*INTERP_surf_diff_delta_t)());
    dFlux_dT->to_F90(grid, surf_diff_dflux_t, (*INTERP_surf_diff_dflux_t)());    // remains const 0 however throughout simulation...
    dq      ->to_F90(grid, surf_diff_delta_tr+(ie-is+1)*(je-js+1)*(q_ind-1), (*INTERP_surf_diff_delta_tr)()); // INDEX ??
    dFlux_dq->to_F90(grid, surf_diff_dflux_tr+(ie-is+1)*(je-js+1)*(q_ind-1), (*INTERP_surf_diff_dflux_tr)());   // remains const 0 however throughout simulation... // INDEX ??
    dt_mass ->to_F90(grid, surf_diff_dtmass, (*INTERP_surf_diff_dtmass)());     // remains const however throughout simulation...
    du      ->to_F90(grid, surf_diff_delta_u, (*INTERP_surf_diff_delta_u)());    // remains const 0 however throughout simulation...
    dv      ->to_F90(grid, surf_diff_delta_v, (*INTERP_surf_diff_delta_v)());    // remains const 0 however throughout simulation...

    //! TODO: send CO2 tracer changes to coupler, if CO2 tracer was defined in field_table
    //if (co2_ind > 0) {
    //  XX      ->to_F90(grid, surf_diff_delta_tr+(ie-is+1)*(je-js+1)*(co2_ind-1), (*INTERP_surf_diff_delta_co2); // INDEX ??
    //  dFlux_XX->to_F90(grid, surf_diff_dflux_tr+(ie-is+1)*(je-js+1)*(co2_ind-1), (*INTERP_surf_diff_dflux_co2); // remains const however throughout simulation... // INDEX ??
    //}

    // check if qS undershoots (dummy stores dq since previous timestep)
    assert(humidity->qS_().Min() >= 0.);

    // output after calculations
#if defined(DEBUG_AEOLUS_STOCKS)
    stockmove_down_heat_toa->Write(0.0);
    stockmove_down_heat_out->Write(0.0);
    /* no water stock is moved in downward sweep */
    for(std::vector<Cell>::const_iterator it=grid->SurfaceCellsBegin(); it!=grid->SurfaceCellsEnd(); it++) {
      (*stockmove_down_heat_toa) += (radiation->swr_flux_toa_().Read(it->ID())
                                     - radiation->lwr_up_toa_().Read(it->ID())
                                     ) * it->FirstConnectedZInterface()->Area();
      (*stockmove_down_heat_out)  += (radiation->swr_flux_sf_().Read(it->ID())
                                      - radiation->lwr_flux_sf_().Read(it->ID())
                                      ) * it->FirstConnectedZInterface()->Area();
    }
#ifdef USE_MPI_DOMAIN_DECOMPOSE
    { double globalstock;
      grid->MPI_Comm().Allreduce(&(*stockmove_down_heat_toa)[0], &globalstock, 1, MPI::DOUBLE, MPI::SUM);
      stockmove_down_heat_toa->Write(globalstock);
      grid->MPI_Comm().Allreduce(&(*stockmove_down_heat_out)[0], &globalstock, 1, MPI::DOUBLE, MPI::SUM);
      stockmove_down_heat_out->Write(globalstock);
    }
#endif
    (*stockmove_down_heat_toa)  *= dt_atm_sec;
    (*stockmove_down_heat_out)  *= dt_atm_sec;
#if defined(TRACE_AEOLUS_UPDATE)
    stockflux_down_heat_out = stockmove_down_heat_out->Read();
    printf("%s: heat  stock moved toa  into Aeolus %22g\n", __FUNCTION__, stockmove_down_heat_toa->Read());
    printf("%s: heat  stock moved down from Aeolus %22g\n", __FUNCTION__, stockflux_down_heat_out->Read());
#endif
#endif
    //diag_ncout->OutputToFile( *Time_days + (*Time_seconds + debug_out_counter++) / (24.0 * 60.0 * 60.0));
#ifdef TRACE_AEOLUS_UPDATE
    //dtmass_down == dt_mass
    //dflux_t_down == dFlux_dT
    //dflux_sphum_down == dFlux_dq
    //delta_t_down == dT
    //delta_sphum_down == dq
    //delta_u_down == du
    //delta_v_down == dv
    Trace_Aeolus_Update_Down_Out->OutputToFile(*Time_days + (*Time_seconds / (24.0 * 60.0 * 60.0)));
#endif
    DEBUG_STOCK_CHECK("aeolus_atmosphere_down_", "I")

#ifdef USE_CLOCK
    clock_Aeolus_down.stop();
    clock_Aeolus.stop();
#endif
    fmsclock_end_(&clock_id_aeolus_down);
    fmsclock_end_(&clock_id_aeolus);
  } // end aeolus_atmosphere_down_()
} // extern "C"



extern "C" {

  void aeolus_atmosphere_up_
  (
   const int *const Time_year,
   const int *const Time_month,
   const int *const Time_seconds, ///< seconds of day
   const int *const Time_days,    ///< days since start of time axis, not day of year!

   // all following parameters are 2-D arrays,
   // defined on the local data domain with index range (is:ie,js:je)
   // input parameter
   const double *const land_frac,  // fraction amount of land in a grid box
   // two input parameter fields from Surf_diff_type
   const double *const surf_diff_delta_t,  // the increment in temperature in the lowest atmospheric
                                           // layer (((i+1)-(i-1) if atmos model is leapfrog) (K)
                                           // (defined in  gcm_vert_diff_down as the increment computed up
                                           // to this point in model, including effect of vertical
                                           // diffusive flux at top of lowest model level, presumed
                                           // to be modified to include effects of surface fluxes
	                                   // outside of this module, then used to start upward
	                                   // tridiagonal sweep,
   const double *const surf_diff_delta_tr, // 3D, 3rd dim for tracer id.
                                           // - similarly for the increment in specific humidity
                                           //   (non-dimensional  = Kg/Kg)
   // output parameters
   double *const lprec,  // mass of liquid precipitation since last time step (Kg/m2/s)
   double *const fprec,  // mass of frozen precipitation since last time step (Kg/m2/s)
   double *const gust,   // gustiness factor
   // input parameters again
   // u_star, b_star, q_star are not used in EBM atmosphere_up
   const double *const u_star, // friction velocity
   const double *const b_star, // bouyancy scale
   const double *const q_star, // moisture scale
   const int *const fortytwo   // minimal interface check
   ) { // begin aeolus_atmosphere_up
#ifdef USE_CLOCK
    clock_Aeolus.start();
    clock_Aeolus_up.start();
#endif
    fmsclock_begin_(&clock_id_aeolus);
    fmsclock_begin_(&clock_id_aeolus_up);
    assert(42 == *fortytwo);
#ifndef NDEBUG
    static unsigned int n_atmos_up_steps = 0;
    n_atmos_up_steps++;
    if (pe == root_pe) {
      cout << __FUNCTION__ << " timestep nr " << n_atmos_up_steps << " at " << *Time_days << ":" << *Time_seconds << " day:sec, "
           << (*Time_days-first_day)/365 << "Y" << (*Time_days-first_day)%356 << "d" << *Time_seconds << "s after start of run" << endl;
    }
#endif

#ifdef DEBUG_AEOLUS_STOCKS
    double water1, water2, heat1, heat2;
    const int iwater = ISTOCK_WATER, iheat = ISTOCK_HEAT;
    aeolus_get_stock_pe_(&iwater, &water1);
    aeolus_get_stock_pe_(&iheat,  &heat1);
#endif

    // from delta_t (K) the sensible heat flux needs to be calculated:
    // this should be done using dt_mass = dt/mass lowest atm. layer = [s m2 / kg]
    // from delta_tr[:,:,q_ind] (kg/kg) the amount of evaporation needs to be calculated in [kg/m2/s]
    // again using dt_mass

    // dT and dq still store change in temperature and humidity due to downward sweep
    // surf_diff_delta_t & surf_diff_delta_tr[:,:,q_ind] store change in temperature and humidity since previous timestep
    // needed: change in temperature and humidity due to surface flux only: 
    // = surf_diff_delta_t - dT   and    surf_diff_delta_tr[:,:,q_ind] - dq
    static HorizontalVariable dummy(*grid, SF, "dummy","","just for I/O", Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_INTERPOL);
    dummy.from_F90(grid, surf_diff_delta_t);
    dummy.PRINTSUM1("surf_diff_delta_t");
#ifdef TRACE_AEOLUS_UPDATE
    dt_t_up->CopyDataFrom(dummy);
#endif
    (*dT) -= dummy;  // dT - surf_diff_delta_t
    (*dT) *= -1.;    // surf_diff_delta_t - dT
    //(*dT) = dummy;  // dT - surf_diff_delta_t // Levke 5-Aug-14
    dT->DebugIsValid(grid, __FUNCTION__);

    dummy.from_F90(grid, surf_diff_delta_tr+(ie-is+1)*(je-js+1)*(q_ind-1));    // [kg/kg] // INDEX ??
    dummy.PRINTSUM1("surf_diff_delta_tr[q]");
#ifdef TRACE_AEOLUS_UPDATE
    dt_sphum_up->CopyDataFrom(dummy);
#endif
    (*dq) -= dummy;  // dq - surf_diff_delta_tr[:,:,q_ind]
    (*dq) *= -1.;    // surf_diff_delta_tr - dq
    //(*dq) = dummy;  // dq - surf_diff_delta_tr[:,:,q_ind] // Levke 5-Aug-14
    dq->DebugIsValid(grid, __FUNCTION__);
    LoadSurfaceFluxesFromF90( *dT, *dq, *sf_layer ); // convert dT and dq to F_sens and evapo
#if defined(DEBUG_AEOLUS_STOCKS)
    stockmove_up_heat_in->Write(0.0); stockmove_up_water_in->Write(0.0);
    for(std::vector<Cell>::const_iterator it=grid->SurfaceCellsBegin(); it!=grid->SurfaceCellsEnd(); it++) {
      // latent heat stock move can be calculated only after ENERGY_CYCLE-> CalculateUpwardSweep(),
      // because precipitation is computed in there
      (*stockmove_up_heat_in)  += sf_layer->F_sens_().Read(it->ID()) * it->FirstConnectedZInterface()->Area(); // sensible heat
      (*stockmove_up_water_in) += sf_layer->evapo_().Read(it->ID()) * it->FirstConnectedZInterface()->Area(); // evaporation
    }
#ifdef USE_MPI_DOMAIN_DECOMPOSE
    { double globalstock;
      grid->MPI_Comm().Allreduce(&(*stockmove_up_water_in)[0], &globalstock, 1, MPI::DOUBLE, MPI::SUM);
      stockmove_up_water_in->Write(globalstock);
    }
#endif
    // (*stockmove_up_heat_in)  *= dt_atm_sec; do it only below, after adding latent heat
    (*stockmove_up_water_in) *= dt_atm_sec;
#if defined(TRACE_AEOLUS_UPDATE)
    stockflux_up_heat_in = stockmove_up_heat_in->Read();
    stockflux_up_water_in = stockmove_up_water_in->Read();
    printf("%s: water stock moved up into Aeolus %22g\n", __FUNCTION__, stockflux_up_water_in);
    printf("%s: heat  stock moved up into Aeolus %22g\n", __FUNCTION__, stockflux_up_heat_in);
#endif
#endif


    //! TODO: get CO2 flux from coupler, if CO2 tracer was defined in field_table
    //if (co2_ind > 0) {
    //  dummy.from_F90(grid, surf_diff_delta_tr+(ie-is+1)*(je-js+1)*(co2_ind-1)); // INDEX ??
    //  dummy.PRINTSUM1("surf_diff_delta_tr[co2]");
    //}

#ifdef TRACE_AEOLUS_UPDATE
    Trace_Aeolus_Update_Up_In->OutputToFile(*Time_days + (*Time_seconds / (24.0 * 60.0 * 60.0)));
#endif

    // update Atmospheric heat and water content
    VAPOR_CYCLE -> CalculateUpwardSweep(dt_atm_sec);    // add evaporation to moist balance, update all affected variables, subtract precipitation
    DEBUG_STOCK_CHECK("aeolus_atmosphere_up_", "A")
    // here we can implement the split-up between liquid and frozen precipitation
    // However, this function is not allowed to change humidity->precip .
    // Thus, we do it inside VAPOR_CYCLE -> CalculateUpwardSweep()

    ENERGY_CYCLE-> CalculateUpwardSweep(dt_atm_sec);    // add sensible heat flux to energy balance, update all affected variables
    DEBUG_STOCK_CHECK("aeolus_atmosphere_up_", "B")

    // output before calculations
    //diag_ncout->OutputToFile( *Time_days + ( *Time_seconds + debug_out_counter++) / (24.0 * 60.0 * 60.0) );

    // return precip variables
    //dummy.CopyDataFrom(humidity->precip_()); // precip is in kg/m2/s
    //dummy *= dt_atm_sec;                     // dummy (for FMS) is in kg/m2: NO, should be per-seconds!
    //dummy.to_F90(grid, lprec, INTERP_lprec);               // store in lprec
    //since we do not need to multiply with timestep anymore, we can avoid the copy in/out of dummy
    humidity->lprec_().to_F90(grid, lprec, (*INTERP_lprec)()); // liquid precip is in kg/m2/s
#ifdef TRACE_AEOLUS_UPDATE
    //lprec_up->CopyDataFrom(dummy);
    //lprec_up->CopyDataFrom(humidity->lprec_());
#endif

    //dummy.UniformValue(0.);
    //dummy.to_F90(grid, fprec, INTERP_fprec); // not implemented yet: frozen precipitation is zero
    humidity->fprec_().to_F90(grid, fprec, (*INTERP_fprec)()); // frozen precip is in kg/m2/s
#ifdef TRACE_AEOLUS_UPDATE
    //fprec_up->CopyDataFrom(dummy);
    //fprec_up->CopyDataFrom(dummy);
#endif

#ifdef DEBUG_AEOLUS_STOCKS
    stockmove_up_heat_out->Write(0.0);
    stockmove_up_water_out->Write(0.0);
    for(std::vector<Cell>::const_iterator it=grid->SurfaceCellsBegin(); it!=grid->SurfaceCellsEnd(); it++) {
      //(*stockmove_up_heat_in) += humidity->precip_().Read(it->ID()) * Atmosphere.L() * it->FirstConnectedZInterface()->Area(); // latent heat
      (*stockmove_up_heat_in) += humidity->lprec_().Read(it->ID()) * Atmosphere.L() * it->FirstConnectedZInterface()->Area()
        + humidity->fprec_().Read(it->ID()) * Atmosphere.HLF() * it->FirstConnectedZInterface()->Area(); // latent heat
      //(*stockmove_up_heat_out)  += dT->Read(it->ID()) * it->FirstConnectedZInterface()->Area(); leave it at zero
      (*stockmove_up_water_out) += humidity->precip_().Read(it->ID()) * it->FirstConnectedZInterface()->Area();
    }
#ifdef USE_MPI_DOMAIN_DECOMPOSE
    { double globalstock;
      grid->MPI_Comm().Allreduce(&(*stockmove_up_heat_in)[0], &globalstock, 1, MPI::DOUBLE, MPI::SUM);
      stockmove_up_heat_in->Write(globalstock);
      grid->MPI_Comm().Allreduce(&(*stockmove_up_heat_out)[0], &globalstock, 1, MPI::DOUBLE, MPI::SUM);
      stockmove_up_heat_out->Write(globalstock);
      grid->MPI_Comm().Allreduce(&(*stockmove_up_water_out)[0], &globalstock, 1, MPI::DOUBLE, MPI::SUM);
      stockmove_up_water_out->Write(globalstock);
    }
#endif
    (*stockmove_up_heat_in)  *= dt_atm_sec;
    (*stockmove_up_heat_out)  *= dt_atm_sec;
    (*stockmove_up_water_out) *= dt_atm_sec;
#if defined(TRACE_AEOLUS_UPDATE)
    stockflux_up_heat_out = stockmove_up_heat_out->Read();
    stockflux_up_water_out = stockmove_up_water_out->Read();
    printf("%s: water stock moved up outof Aeolus %22g\n", __FUNCTION__, stockflux_up_water_out);
    printf("%s: heat  stock moved up outof Aeolus %22g\n", __FUNCTION__, stockflux_up_heat_out);
#endif
#endif

    // output to file (all atmosphere calculations done for this timestep)
    //count_out++;
    //if (count_out == 24) // hardcoded end of day output
    //{
    //  diag_ncout->OutputToFile( *Time_days + (*Time_seconds + debug_out_counter++) / (24.0 * 60.0 * 60.0));
    //  count_out = 0;
    //}

#ifdef TRACE_AEOLUS_UPDATE
    Trace_Aeolus_Update_Up_Out->OutputToFile(*Time_days + (*Time_seconds / (24.0 * 60.0 * 60.0)));
#endif
#ifdef USE_CLOCK
    clock_Aeolus.stop();
    clock_Aeolus_up.stop();
#endif
    fmsclock_end_(&clock_id_aeolus_up);
    fmsclock_end_(&clock_id_aeolus);

#ifdef WITH_SIMENV
    simenv_output->OutputToSimenv();
#endif

#ifdef WITH_COSTFUNCTION
    Aeolus_Score_Add_Tstep(*Time_year, *Time_month,
                           *Time_days, // days since start of time axis, not day of year!
                           *Time_seconds); // seconds of day
#endif

  } // end of aeolus_atmosphere_up_()



  /*
   * aeolus_get_bottom_mass_() and aeolus_get_bottom_wind_() are called
   * - during init just after aeolus_init_()
   * - during upward update, just after aeolus_atmosphere_up_()
   */

  // t_bot, p_bot, z_bot, p_surf, slp: output parameters, 2-D arrays
  // tr_bot: output parameter, 3-D array, 3rd dimension is for tracer index
  // defined on the local data domain with index range (is:ie,js:je)
  void aeolus_get_bottom_mass_
  (
   double *const t_bot,  // near surface temperature (at lowest model level) [deg K]
   double *const tr_bot, // tracers at lowest model level, 1st tracer is near surface mixing  ratio
   double *const p_bot,  // pressure at which atmos near surface values are assumed to be defined (at lowest model level)
   double *const z_bot,  // height above the surface for the lowest model level
   double *const p_surf, // surface pressure
   double *const slp     // sea level pressure
   )
  {
    fmsclock_begin_(&clock_id_aeolus);
    // load Horizontal Variables from 3D aeolus variables. Disadvantage: lots of memory copying!
    /*T_bot->ImportFrom( temperature->T_() );            // loads lowest level of T on T_bot
     *q_bot->ImportFrom( humidity->q_() );               // loads lowest level of q on q_bot
     *P_bot->ImportFrom( wind->p_() );                   // loads lowest level of p on p_bot
     *
     * // convert to FMS double*
     *T_bot->to_F90(grid, t_bot);
     *q_bot->to_F90(grid, tr_bot); // tracer index !!
     *P_bot->to_F90(grid, p_bot);
     */
    Z_bot->PRINTSUM;
    Z_bot->to_F90(grid, z_bot, (*INTERP_z_bot)());  // remains constant ...

    (temperature->realTSa_()).PRINTSUM;
    (temperature->realTSa_()).to_F90(grid, t_bot, (*INTERP_t_bot)());
    (humidity->qS_()).PRINTSUM;
    (humidity->qS_()).to_F90(grid, tr_bot+(ie-is+1)*(je-js+1)*(q_ind-1), (*INTERP_tr_bot)());    // INDEX ??
    //! TODO: the same for CO2 tracer tr_bot+(ie-is+1)*(je-js+1)*(co2_ind-1)
    (wind->sfp_()).PRINTSUM;
    (wind->sfp_()).to_F90(grid, p_bot, (*INTERP_p_bot)()); // old: wind->p2m_()
    (wind->sfp_()).to_F90(grid, p_surf, (*INTERP_p_surf)());
    (wind->slp_()).PRINTSUM;
    (wind->slp_()).to_F90(grid, slp, (*INTERP_slp)());
    fmsclock_end_(&clock_id_aeolus);
  }


  /// \param u_bot and v_bot: output parameters, 2-D arrays
  /// defined on the local data domain with index range (is:ie,js:je)
  void aeolus_get_bottom_wind_
  (
   double *const u_bot, ///< zonal wind component at lowest model level
   double *const v_bot  ///< meridional wind component at lowest model level
   )
  { 
    fmsclock_begin_(&clock_id_aeolus);
    (pbl->us_eff_()).PRINTSUM;
    (pbl->us_eff_()).to_F90(grid, u_bot, (*INTERP_u_bot)());
    (pbl->vs_eff_()).PRINTSUM;
    (pbl->vs_eff_()).to_F90(grid, v_bot, (*INTERP_v_bot)());
    fmsclock_end_(&clock_id_aeolus);
  }


#ifdef DEBUG

  /// a function just for testing the F90 / C++ interface
  // frac_land: 2-D array, input parameter
  // t_bot, p_bot: output parameters, 2-D arrays
  // tr_bot: output parameter, 3-D array, 3rd dimension is for tracer index
  // defined on the local data domain with index range (is:ie,js:je)
  void aeolus_test_var_access_
  (
   const double *const frac_land,
   double *const t_bot, double *const tr_bot, double *const p_bot				
   )
  {

    // allocate atmosphere variables,
    HorizontalVariable *hor = new HorizontalVariable(*grid, SF, "testvar", "1", "dummy test variable",
						     Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_DIVIDE);
    HorizontalVariable *hor2 = new HorizontalVariable(*grid, SF, "testvar", "1", "dummy test variable",
						     Variable::INTERP_IN_SUM|Variable::INTERP_OUT_DIVIDE);
    HorizontalVariable *hor3 = new HorizontalVariable(*grid, SF, "testvar", "1", "dummy test variable",
						     Variable::INTERP_IN_AVERAGE|Variable::INTERP_OUT_INTERPOL);

    // initial data
    for (unsigned int i = 0; i < grid->ZonalCells(); i++) {
      for (unsigned int j = 0; j < grid->MeridionalCells(); j++) {
	hor->Write(grid->GetCell(i, j, 0)->ID(), i*1000+j);
	printf("%u %u %g\n", i, j, hor->Read(grid->GetCell(i, j, 0)->ID()));
      }
    }

    // assign some variables
    hor->to_F90(grid, t_bot, true);
    hor->to_F90(grid, p_bot);
    hor->to_F90(grid, tr_bot+((q_ind-1) * (je-js+1) * (ie-is+1)));

    // assign some variables in different ways
    hor->from_F90(grid, frac_land, true);
    hor->from_F90(grid, frac_land);
    hor2->from_F90(grid, frac_land);

    delete hor;
    delete hor2;

    fflush(stdout);
  }
#endif /* DEBUG */

} /* extern "C" */


static void FMSAddVariable(const Variable& v) {
  if (!FMS_open_for_registration) {
    fprintf(stderr, "%s: Error: registration of Aeolus Variables for FMS diagnostic manager is closed already.\n"
            "Attempt to add Variable %s after first sending of data\n",
            __FUNCTION__, v.Name().c_str());
    fflush(stdout);
    fflush(stderr);
    abort();
  }
  if (FMS_registered_vars.size() > 0) {
    std::vector<const Variable*>::iterator it;
    it = find(FMS_registered_vars.begin(), FMS_registered_vars.end(), &v);
    if (it != FMS_registered_vars.end()) { // variable exists in vector already
      fprintf(stderr, "%s: Warning: variable %s already registered for FMS diagnostic manager, nothing is done ...\n",
              __FUNCTION__, v.Name().c_str());
      return;
    }
    for (it = FMS_registered_vars.begin(); it != FMS_registered_vars.end(); it++) {
      if ((*it)->Name() == v.Name()) {
        fprintf(stderr, "%s: Warning: a different Variable but with same name %s was already registered for FMS diagnostic manager, nothing is done ...\n",
                __FUNCTION__, v.Name().c_str());
        return;
      }
    }
  } 
  FMS_registered_vars.push_back(&v);
  fprintf(stderr, "%s: Info: %lu variable %s registered for FMS diagnostic manager\n",
          __FUNCTION__, FMS_registered_vars.size(), v.Name().c_str());
}


static void FMSAddVariable(const LargeScaleWind *data ) {
  FMSAddVariable( data->u_() );	FMSAddVariable( data->v_() );	FMSAddVariable( data->w_() ); 
  FMSAddVariable( data->dudx_() ); FMSAddVariable( data->dudy_() ); FMSAddVariable( data->dudz_() );
  FMSAddVariable( data->dvdx_() ); FMSAddVariable( data->dvdy_() ); FMSAddVariable( data->dvdz_() );
  FMSAddVariable( data->u_za_() ); FMSAddVariable( data->v_za_() ); FMSAddVariable( data->w_za_() ); 
  FMSAddVariable( data->v_za_unsmoothed_() );
  FMSAddVariable( data->p_() ); FMSAddVariable( data->slp_() ); FMSAddVariable( data->sfp_() );
  FMSAddVariable( data->p_za_() ); 
  FMSAddVariable( data->dslp_dx_() );FMSAddVariable( data->dslp_dy_() );FMSAddVariable( data->dslp_dy_za_() );
  FMSAddVariable( data->PSI_za_() );
  FMSAddVariable( data->u_nt_() ); FMSAddVariable( data->v_nt_() );
  FMSAddVariable( data->v_za_conv_() );
  FMSAddVariable( data->denominator_1_() ); FMSAddVariable( data->denominator_2_() );
  FMSAddVariable( data->denominator_3_() ); FMSAddVariable( data->denominator_4_() );
  FMSAddVariable( data->nominator_1_() ); FMSAddVariable( data->nominator_2_() );
  FMSAddVariable( data->nominator_3_() );
}
static void FMSAddVariable(const SynopticScaleField *data ) {
  FMSAddVariable( data->Tu_() );	FMSAddVariable( data->Tv_() );	/*FMSAddVariable( data->Tw_() );*/
  FMSAddVariable( data->qu_() );	FMSAddVariable( data->qv_() );	/*FMSAddVariable( data->qw_() );*/
  FMSAddVariable( data->uu_() );	FMSAddVariable( data->uv_() );	FMSAddVariable( data->uw_() );
  FMSAddVariable( data->vv_() );	FMSAddVariable( data->ww_() );    FMSAddVariable( data->vw_() );
  FMSAddVariable( data->Ksyn_() );	FMSAddVariable( data->fF_ag_() );FMSAddVariable( data->L_Ro_() );
  FMSAddVariable( data->alfa_() );  FMSAddVariable( data->KS_() );    FMSAddVariable( data->Ah_() );
}
static void FMSAddVariable(const Temperature *data ) {
  FMSAddVariable( data->T_() );	FMSAddVariable( data->Theta_() ); 
  FMSAddVariable( data->dThetadx_() ); FMSAddVariable( data->dThetady_() );
  FMSAddVariable( data->LR_() );	FMSAddVariable( data->LR_wa_() ); FMSAddVariable( data->LR_za_() );
  FMSAddVariable( data->dTdx_() );	FMSAddVariable( data->dTdy_() ); 
  FMSAddVariable( data->Ta_() );	FMSAddVariable( data->T_za_() ); FMSAddVariable( data->TSa_() );
  FMSAddVariable( data->realTSa_() );
  FMSAddVariable( data->T_sl_() );  FMSAddVariable( data->LR_pot_() );
  FMSAddVariable( data->QT_() ); FMSAddVariable( data->Divergence_() ); 
  FMSAddVariable( data->large_transport_x_() );  FMSAddVariable( data->large_transport_y_() );
  FMSAddVariable( data->synop_transport_x_() );  FMSAddVariable( data->synop_transport_y_() );
  FMSAddVariable( data->meso_transport_x_() );  FMSAddVariable( data->meso_transport_y_() );
  FMSAddVariable( data->dTdphi_() );
  //FMSAddVariable( data->thermal_equator_(), "thermal_equator", "deg" );
}
static void FMSAddVariable(const Humidity *data ) {
  FMSAddVariable( data->q_() ); FMSAddVariable( data->qS_() ); FMSAddVariable( data->Divergence_() ); 
  FMSAddVariable( data->He_() );
  FMSAddVariable( data->dqdx_() );	FMSAddVariable( data->dqdy_() ); // FMSAddVariable( data->dqdz_() ); 
  FMSAddVariable( data->Qq_() ); FMSAddVariable( data->precip_() );
  //FMSAddVariable( data->preciptest_() );
  FMSAddVariable( data->lprec_() ); FMSAddVariable( data->fprec_() );
  FMSAddVariable( data->q_sat_() ); FMSAddVariable( data->P_over_sat_() ); 
  FMSAddVariable( data->large_transport_x_() );  FMSAddVariable( data->large_transport_y_() );
  FMSAddVariable( data->synop_transport_x_() );  FMSAddVariable( data->synop_transport_y_() );
  FMSAddVariable( data->meso_transport_x_() );  FMSAddVariable( data->meso_transport_y_() );
  FMSAddVariable( data->rhum_() ); FMSAddVariable( data->rh_sf_() ); FMSAddVariable( data->qS_sat_() );
}
static void FMSAddVariable(const Clouds_2 *data ) {
  //FMSAddVariable( data->cumulus_() ); FMSAddVariable( data->stratus_() ); FMSAddVariable( data->total_() ); 
  //FMSAddVariable( data->Fc_() );  FMSAddVariable( data->w_eff_() ); 
}
static void FMSAddVariable(const Clouds_3Layers *data ) {
  FMSAddVariable( data->iconvect_() );
  FMSAddVariable( data->prectot_()  );
  FMSAddVariable( data->precls_()   );
  FMSAddVariable( data->precconv_() );
  for(unsigned int i=0;i!=NCLOUD;i++) {
    FMSAddVariable( data->clam_(i) );
    FMSAddVariable( data->clwattot_(i) );
    FMSAddVariable( data->hbase_(i) );
    FMSAddVariable( data->htop_(i) );
    FMSAddVariable( data->w_eff_(i) );
  }
  FMSAddVariable( data->clam_(NCLOUD) );
  FMSAddVariable( data->clwattot_(NCLOUD) );
  FMSAddVariable( data->stratus_() );
}
static void FMSAddVariable(const PlanetaryWave *data ) {
  FMSAddVariable( data->PSI_() ); FMSAddVariable( data->PSI_therm_() ); FMSAddVariable( data->PSI_oro_() ); 
  FMSAddVariable( data->w_oro_() );
  FMSAddVariable( data->ua_() ); FMSAddVariable( data->va_() ); FMSAddVariable( data->wa_() );
  FMSAddVariable( data->uava_() );
  FMSAddVariable( data->ua_ageos_() ); FMSAddVariable( data->va_ageos_() );
  FMSAddVariable( data->ua_geos_() ); FMSAddVariable( data->va_geos_() );
  FMSAddVariable( data->pa_() ); FMSAddVariable( data->am_() );
  FMSAddVariable( data->INTv_() );
}
static void FMSAddVariable(const PlanetaryBoundaryLayer *data ) {
  FMSAddVariable( data->alfa_() ); FMSAddVariable( data->alfa_za_() ); /*FMSAddVariable( data->Ts_pbl_() );*/
  /*FMSAddVariable( data->T_skin_() );*/ FMSAddVariable( data->us_() ); FMSAddVariable( data->vs_() );
  FMSAddVariable( data->us_eff_() ); FMSAddVariable( data->vs_eff_() );
  FMSAddVariable( data->ua_ageos_() ); FMSAddVariable( data->va_ageos_() ); FMSAddVariable( data->H_PBL_() );
  FMSAddVariable( data->H_PBL_diag_() );
  FMSAddVariable( data->C_drag_() ); FMSAddVariable( data->tau_x_() ); FMSAddVariable( data->tau_y_() );
  FMSAddVariable( data->uu_() ); FMSAddVariable( data->vv_() );
  FMSAddVariable( data->TS_() ); FMSAddVariable( data->ug_top_pbl_() ); FMSAddVariable( data->vg_top_pbl_() );
}
static void FMSAddVariable(const Topography *data ) {
  FMSAddVariable( data->real_height_() );
  FMSAddVariable( data->real_sigma_() );
  FMSAddVariable( data->height_() );
  FMSAddVariable( data->sigma_() );
  FMSAddVariable( data->fraction_land_() );
}
//static void FMSAddVariable(const EmpiricalData *data ) {
//  FMSAddVariable( data->Tu_() );FMSAddVariable( data->Tv_() );FMSAddVariable( data->Tw_() );
//  FMSAddVariable( data->qu_() );FMSAddVariable( data->qv_() );FMSAddVariable( data->qw_() );
//  FMSAddVariable( data->uu_() );FMSAddVariable( data->vv_() );
//  FMSAddVariable( data->uv_() );FMSAddVariable( data->uw_() );FMSAddVariable( data->vw_() );
//}
static void FMSAddVariable(const RadiativeFluxes *data ) {
  FMSAddVariable( data->swr_flux_sf_() );FMSAddVariable( data->swr_flux_toa_() );FMSAddVariable( data->lwr_up_sf_() );
  FMSAddVariable( data->lwr_flux_sf_() );FMSAddVariable( data->lwr_down_sf_() );FMSAddVariable( data->lwr_up_toa_() );
  FMSAddVariable( data->plan_albedo_() );FMSAddVariable( data->RBTP_() );FMSAddVariable( data->RBSR_() );
  FMSAddVariable( data->T_c_() ); FMSAddVariable( data->hstrat_() );
  FMSAddVariable( data->ITF_() ); FMSAddVariable( data->Htrop_() ); FMSAddVariable( data->rad_flux_toa_() );
  FMSAddVariable( data->swr_down_toa_() );
  for (int i = 0; i < NST; i++) FMSAddVariable( *(data->totrad_down_sf_()[i]) );
#ifdef DEBUG_ITF
  FMSAddVariable(data->D_LWR_()); FMSAddVariable(data->D_HTR_());
#endif
}
static void FMSAddVariable(const Atm *data ) {
  FMSAddVariable( data->rho_() ); FMSAddVariable( data->F0_() ); FMSAddVariable( data->f_() ); FMSAddVariable( data->f_beta_() );
}
static void FMSAddVariable(const SurfaceLayer *data ) {
  FMSAddVariable( data->roughness_length_() );
  FMSAddVariable( *(data->alb_vis_clear_()[0]) );  // surface type 0 == ocean
  FMSAddVariable( *(data->alb_vis_cloudy_()[0]) );
  FMSAddVariable( *(data->alb_ir_clear_()[0]) );
  FMSAddVariable( *(data->alb_ir_cloudy_()[0]) );
  FMSAddVariable( *(data->alb_vis_clear_()[2]) ); // surface type 2 == land (originally wood in Climber-2) 
  FMSAddVariable( *(data->alb_vis_cloudy_()[2]) );
  FMSAddVariable( *(data->alb_ir_clear_()[2]) );
  FMSAddVariable( *(data->alb_ir_cloudy_()[2]) );
  /*FMSAddVariable( data->Vs_() ); */
  FMSAddVariable( data->F_sens_() ); 
  FMSAddVariable( data->evapo_() );
  FMSAddVariable( data->T_skin_() );
  FMSAddVariable( data->q_sat_() );
  FMSAddVariable( data->q_soil_() );
}


/* a simple implementation of strncpy(),
 * BUT it fills up the destination with blanks, while
 * the real strncpy() fills it up with null bytes.
 * The blank-filling is needed for inter-operation with Fortran character type
 */
static void strncpyblank(char *dest, const char *src, size_t n) {
  size_t i;

  for (i = 0 ; i < n && src[i] != '\0' ; i++)
    dest[i] = src[i];
  for ( ; i < n ; i++)
    dest[i] = ' ';
}

extern "C"
void aeolus_fms_variable_metadata_(const int *n_in,
                                     char *name, char *long_name, char *units,
                                     double *missing_value, int *Variable_type,
                                     int *timedependent,
                                     int namelen, int long_namelen, int unitslen) {
  const int n = (*n_in) - 1; // Fortran to C index
  if ( n < 0 || n >= FMS_registered_vars.size()) {
    fprintf(stderr, "%s: Error: n == %d is outside 0 ... %lu\n",
            __FUNCTION__, n, FMS_registered_vars.size()-1);
    fflush(stderr);
    fflush(stdout);
    abort();
  }
  const Variable *v = FMS_registered_vars[n];
  //printf("%s %d %s %d %d %d\n", __FUNCTION__, *n_in, v->Name(), namelen, long_namelen, unitslen);
  strncpyblank(name, v->Name().c_str(), namelen);
  strncpyblank(long_name, v->Longname().c_str(), long_namelen);
  strncpyblank(units, v->Unit().c_str(), unitslen);
  *missing_value = v->MissingValue();
  *Variable_type = v->Type();
  *timedependent = v->IsTimedependent();
}

/// Called from the Fortran side to access that current content of the to-be-diagnosed Variables.
extern "C"
void aeolus_fms_variable_data_(const int *n_in, double *f) {
  const int n = (*n_in) - 1; // Fortran to C index
  FMS_open_for_registration = false;
  if ( n < 0 || n >= FMS_registered_vars.size()) {
    fprintf(stderr, "%s: Error: n == %d is outside 0 ... %lu\n",
            __FUNCTION__, n, FMS_registered_vars.size()-1);
    fflush(stderr);
    fflush(stdout);
    abort();
  }
  const Variable *v = FMS_registered_vars[n];
  //printf("%s: %d getting data for %s type %d\n", __FUNCTION__, *n_in, v->Name().c_str(), v->Type()); fflush(stdout);
  switch(v->Type()) {
  case CELL: // CellVariable
    {
      const CellVariable *cv = static_cast<const CellVariable *>(v);
      cv->to_F90(grid, f, true, is, ie, js, je);
    }
    break;
  case LVL:   // VerticalVariable
  case MV:    // Meridional Variable
  case MV_1D: // 1D Meridional Variable
    for (unsigned int i=0, dit = v->DataBegin(); dit!=v->DataEnd(); i++, dit++)
      f[i] = (*v)[dit];
    break;
    // there are 4 types of HorizontalVariable
  case HOR:
  case SF:
  case EBL:
  case PBL:
    {
      const HorizontalVariable *hv = static_cast<const HorizontalVariable *>(v);
      hv->to_F90(grid, f, true, is, ie, js, je);
      //if (!strcmp(v->Name().c_str(), "PRECIPITATION")) {
      //  double min = hv->Min();
      //  double max = hv->Max();
      //  printf("  Min %g Max %g\n", min, max);
      //  for (unsigned int i = 0; i < 6; i++) printf(" %g", hv->Read(i));
      //  printf("\n");
      //  for (unsigned int i = hv->Size()-6; i < hv->Size(); i++) printf(" %g", hv->Read(i));
      //  printf("\n");
      //}
    } // end of scope for HorizontalVariable *hv
    break;
    // values of InterfaceVariable are interpolated to CellVariable before writing them out.
  case X:
  case Y:
  case Z:
    {
      const InterfaceVariable *iv = static_cast<const InterfaceVariable *>(v);
      static CellVariable cv(*grid, "InterfaceVariableOutput_dummy_cv", "", "", Variable::INTERP_OUT_REPLICATE|Variable::INTERP_IN_SIMPLE, 0, -1);
      InterpolateFromTo(*grid, *iv, cv);
      cv.to_F90(grid, f, true, is, ie, js, je);
    }
    break;
  case SCALAR:
    f[0] = (*v)[0];
    break;
  default:
    cerr << "Error: " << __FUNCTION__ << ": dont know how to handle Variable " <<v->Name()
         <<" with Type " << v->Type() << endl << flush;
  }
}