src/layered/gotm.h

    General Ocean Turbulence Model interface for multilayer

    This module interfaces the multilayer solver with GOTM.

    Vertical diffusion of momentum, temperature and salinity is handled by GOTM.

    The temperature and salinity fields are undefined by default. If the calling solver defines them, the corresponding GOTM vertical diffusion routines will be called.

    scalar T = {-1}, S = {-1};

    The surface fluxes of momentum, heat and short-wave radiation are provided as surface fields by the calling solver. They are zero by default.

    // surface wind stress (eg. in Pa or N/m^2)
    (const) vector airsea_tau      = zerof;
    // surface heat flux (eg. in W/m^2)
    (const) scalar airsea_heat_flux = zeroc;
    // surface short-wave radiation flux (eg. in W/m^2)
    (const) scalar airsea_swr_flux  = zeroc;

    Implementation

    It uses the relevant headers from the C interface to GOTM.

    #include <gotm/common.h>
    #include <gotm/gotm/gotm.h>
    // #include <gotm/gotm/diagnostics.h>
    // #include <gotm/util/time.h>
    
    #include <gotm/meanflow/meanflow.h>
    // #include <gotm/meanflow/updategrid.h>
    // #include <gotm/meanflow/wequation.h>
    #include <gotm/meanflow/coriolis.h>
    #include <gotm/meanflow/uequation.h>
    #include <gotm/meanflow/vequation.h>
    // #include <gotm/meanflow/extpressure.h>
    // #include <gotm/meanflow/intpressure.h>
    #include <gotm/meanflow/friction.h>
    #include <gotm/meanflow/shear.h>
    #include <gotm/meanflow/stratification.h>
    #include <gotm/meanflow/salinity.h>
    #include <gotm/meanflow/temperature.h>
    #include <gotm/meanflow/convectiveadjustment.h>
    
    #include <gotm/util/convert_fluxes.h>
    #include <gotm/util/tridiagonal.h>
    #include <gotm/util/eqstate.h>
    
    #include <gotm/turbulence/turbulence.h>
    #include <gotm/turbulence/kpp.h>
    
    // #include <gotm/observations/observations.h>
    #include <gotm/airsea/airsea.h>

    The corresponding GOTM libraries need to be linked. Note that the GOTM environment variable need to point to the proper installation directory.

    #pragma autolink							\
      -L$GOTM/build/ -lgotm	-lairsea -lmeanflow -lobservations -lturbulence \
      -linput -lutil -loutput_manager  `nf-config --flibs`			\
      -lgfortran

    The timestepping routine is based on its Fortran counterpart.

    static void gotm_step (long n)
    {
      // prepare time and output
      //  time_update_time (&n);
    
      //     all observations/data
    
    #if 0  
      input_do_input (&time_julianday, &time_secondsofday,
    		  &nl, &meanflow_z);
      observations_get_all_obs (&time_julianday, &time_secondsofday,
    			    &nl, &meanflow_z);
    #endif
      
      //     external forcing
    #if 0  
      if (airsea_calc_fluxes) {
        // fixme: check indices
        airsea_set_sst (&meanflow_t.a[nl-1]); // check this
        airsea_set_ssuv (&meanflow_u.a[nl-1],
    		     &meanflow_v.a[nl-1]);
      }
    #endif
      //  airsea_do_airsea (&time_julianday, &time_secondsofday);
    
      //     reset some quantities
      airsea_tx /= meanflow_rho_0;
      airsea_ty /= meanflow_rho_0;
    
      // sum fluxes as a diagnostic
      //  airsea_integrated_fluxes (&dt);
    
      //     meanflow integration starts
    #if 0  
      realtype zeta = observations_get_zeta();
      meanflow_updategrid (&nl, &dt, &zeta);
      meanflow_wequation (&nl, &dt);
    #endif
    
    #if 0
      double omega = 2.*pi/86164.;
      double latitude = 50.; // 50.;
      meanflow_cori = 2.*omega*sin(2.*pi*latitude/360.);
      meanflow_coriolis (&nl, &dt);
    #endif
      
      //     update velocity
      integer ext_press_mode = 99; // no external pressure gradient
      meanflow_uequation (&nl, &dt, &gotm_cnpar, &airsea_tx,
    		      turbulence_num.a, turbulence_gamu.a,
    		      &ext_press_mode);
      meanflow_vequation (&nl, &dt, &gotm_cnpar, &airsea_ty,
    		      turbulence_num.a, turbulence_gamv.a,
    		      &ext_press_mode);
    #if 0
      meanflow_extpressure (&observations_ext_press_mode, &nl);
      meanflow_intpressure (&nl);
    #endif
      meanflow_friction (&turbulence_kappa, &meanflow_avmolu,
    		     &airsea_tx, &airsea_ty);
    
      //     update temperature and salinity
      if (S.i >= 0)
        meanflow_salinity (&nl, &dt, &gotm_cnpar,
    		       turbulence_nus.a, turbulence_gams.a);
      
      if (T.i >= 0)
        meanflow_temperature (&nl, &dt, &gotm_cnpar, &airsea_i_0, &airsea_heat,
    			  turbulence_nuh.a, turbulence_gamh.a, meanflow_rad.a);
      
      //     update shear and stratification
      meanflow_shear (&nl, &gotm_cnpar);
      meanflow_stratification (&nl, &gotm_buoy_method,
    			   &dt, &gotm_cnpar,
    			   turbulence_nuh.a, turbulence_gamh.a);
    
      //    compute turbulent mixing
      switch (turbulence_turb_method) {
    
      case 0:
        //        do convective adjustment
        meanflow_convectiveadjustment (&nl, turbulence_num.a, turbulence_nuh.a,
    				   &turbulence_const_num, &turbulence_const_nuh,
    				   &gotm_buoy_method, &meanflow_gravity,
    				   &meanflow_rho_0);
        break;
    
      case 99: {
        //        update KPP model
        realtype precip = airsea_precip + airsea_evap;
        realtype tFlux, btFlux, sFlux, bsFlux;
        realtype tRad[nl+1], bRad[nl+1]; // fixme: check this
        util_convert_fluxes (&nl, &meanflow_gravity, &meanflow_cp,
    			 &meanflow_rho_0,
    			 &airsea_heat, &precip,
    			 meanflow_rad.a, meanflow_t.a, meanflow_s.a,
    			 &tFlux, &sFlux, &btFlux, &bsFlux,
    			 tRad, bRad);
    
        kpp_do_kpp (&nl, &meanflow_depth,
    		meanflow_h.a, meanflow_rho.a, meanflow_u.a,
    		meanflow_v.a, meanflow_nn.a, meanflow_nnt.a, meanflow_nns.a,
    		meanflow_ss.a,
    		&meanflow_u_taus, &meanflow_u_taub,
    		&tFlux, &btFlux, &sFlux, &bsFlux,
    		tRad, bRad, &meanflow_cori);
        break;
      }
        
      default:
        //        update one-point models
        turbulence_do_turbulence (&nl, &dt, &meanflow_depth,
    			      &meanflow_u_taus,
    			      &meanflow_u_taub, &meanflow_z0s, &meanflow_z0b,
    			      meanflow_h.a, meanflow_nn.a, meanflow_ss.a,
    #ifdef SEAGRASS
    			      xP.a
    #else
    			      NULL
    #endif
    			      );
      }
    }

    The initialisation tries to call the relevant initialisation functions from GOTM, as they appear in the original gotm.F90 source code.

    event defaults (i = 0)
    {
      // cnpar       [float, minimum = 0, maximum = 1]
      //             Constant for the theta scheme used for time integration of
      //             diffusion-reaction components. \theta=0.5 for Cranck-Nicholson
      //             (second-order accurate), \theta=0 for Forward Euler (first-
      //             order accurate), \theta=1 for Backward Euler (first-order
      //             accurate). Note that only \theta=1 guarantees positive
      //             solutions for positive definite systems.
      gotm_cnpar = 1.;
    
      // buoy_method [integer]
      //               method to compute mean buoyancy
      //               1: from equation of state (i.e. from potential temperature
      //                                          and salinity)
      //               2: from prognostic equation
      gotm_buoy_method = 1;
    
      meanflow_gravity = G;
      
      mtridiagonal_init_tridiagonal (&nl);
      realtype latitude = 0.;
      meanflow_post_init_meanflow (&nl, &latitude);
      if (turbulence_turb_method == 99) {
        integer namlst = -1;
        kpp_init_kpp (&namlst, "", &nl, &meanflow_depth, meanflow_h.a,
    		  &G, &meanflow_rho_0);
      }
      else
        turbulence_post_init_turbulence (&nl);
    }

    Vertical diffusion is part of the “viscous terms” of the multilayer event loop. We just copy the relevant fields into the corresponding GOTM fields, call GOTM and retrieve the updated fields.

    event viscous_term (i++)
    {
      struct { realtype * x, * y; } gotm_u = { meanflow_u.a, meanflow_v.a };
      struct { realtype * x, * y; } gotm_t = { &airsea_tx, &airsea_ty };
      foreach() {
        memcpy (&meanflow_h.a[1], &h[], nl*sizeof(double));
        foreach_dimension() {
          *gotm_t.x = airsea_tau.x[];
          memcpy (&gotm_u.x[1], &u.x[], nl*sizeof(double));
        }
        if (T.i >= 0)
          memcpy (&meanflow_t.a[1], &T[], nl*sizeof(double));
        if (S.i >= 0)
          memcpy (&meanflow_s.a[1], &S[], nl*sizeof(double));
        
        double z = zb[];    
        foreach_layer() {
          meanflow_z.a[point.l + 1] = z + h[]/2.;
          z += h[];
        }
    
        airsea_heat = airsea_heat_flux[];
        airsea_i_0 = airsea_swr_flux[];
        gotm_step (i);
        
        foreach_dimension()
          memcpy (&u.x[], &gotm_u.x[1], nl*sizeof(double));
        if (T.i >= 0)
          memcpy (&T[], &meanflow_t.a[1], nl*sizeof(double));
        if (S.i >= 0)
          memcpy (&S[], &meanflow_s.a[1], nl*sizeof(double));
      }
    }

    The GOTM cleanup function deallocates Fortran fields etc.

    event cleanup (t = end)
    {
      gotm_clean_up();
    }

    Utility functions

    This function uses GOTM to initialize a vertical temperature profile corresponding to a squared buoyancy frequency NN for a given constant salinity and top temperature.

    void constant_NNT (double T_top, double S_const, double NN,
    		   scalar T)
    {
      foreach() {
        double z = zb[];
        foreach_layer()
          z += h[];
        T[0,0,nl-1] = T_top;
        for (int l = nl - 2; l >= 0; l--) {
          double pFace = G*(z - h[0,0,l]/2.);
          double alpha = eqstate_eos_alpha (&S_const, &T[0,0,l+1], &pFace, &G,
    					&meanflow_rho_0);
          T[0,0,l] = T[0,0,l+1] - 1./(G*alpha)*NN*h[0,0,l];
          z -= h[0,0,l];
        }
      }
      boundary ({T});
    }

    Usage

    Tests