make stochastics optional

NOAA-GFDL · Jan 28, 2022 · 0bf4ff0 · 0bf4ff0
1 parent 9e7029c
commit 0bf4ff0
Show file tree

Hide file tree

Showing 8 changed files with 181 additions and 137 deletions.
diff --git a/config_src/drivers/FMS_cap/ocean_model_MOM.F90 b/config_src/drivers/FMS_cap/ocean_model_MOM.F90
@@ -50,7 +50,7 @@ module ocean_model_mod
 use MOM_tracer_flow_control, only : call_tracer_register, tracer_flow_control_init
 use MOM_tracer_flow_control, only : call_tracer_flux_init
 use MOM_unit_scaling, only : unit_scale_type
-use MOM_variables,           only : surface, stochastic_pattern
+use MOM_variables, only : surface
 use MOM_verticalGrid, only : verticalGrid_type
 use MOM_ice_shelf, only : initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS
 use MOM_ice_shelf, only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart
@@ -186,7 +186,6 @@ module ocean_model_mod
                               !! timesteps are taken per thermodynamic step.
   type(surface)               :: sfc_state   !< A structure containing pointers to
                               !! the ocean surface state fields.
-  type(stochastic_pattern) :: stochastics !< A structure containing pointers to
   type(ocean_grid_type), pointer :: &
     grid => NULL()            !< A pointer to a grid structure containing metrics
                               !! and related information.
@@ -562,7 +561,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
     ! For now, the waves are only updated on the thermodynamics steps, because that is where
     ! the wave intensities are actually used to drive mixing.  At some point, the wave updates
     ! might also need to become a part of the ocean dynamics, according to B. Reichl.
-    call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves)
+    call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves, OS%forces)
   endif
 
   if ((OS%nstep==0) .and. (OS%nstep_thermo==0)) then ! This is the first call to update_ocean_model.
@@ -577,12 +576,12 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
     call step_offline(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp)
   elseif ((.not.do_thermo) .or. (.not.do_dyn)) then
     ! The call sequence is being orchestrated from outside of update_ocean_model.
-    call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dt_coupling, OS%MOM_CSp, &
+    call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp, &
                   Waves=OS%Waves, do_dynamics=do_dyn, do_thermodynamics=do_thermo, &
                   start_cycle=start_cycle, end_cycle=end_cycle, cycle_length=cycle_length, &
                   reset_therm=Ocn_fluxes_used)
   elseif (OS%single_step_call) then
-    call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves)
+    call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves)
   else  ! Step both the dynamics and thermodynamics with separate calls.
     n_max = 1 ; if (dt_coupling > OS%dt) n_max = ceiling(dt_coupling/OS%dt - 0.001)
     dt_dyn = dt_coupling / real(n_max)
@@ -604,16 +603,16 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
             "THERMO_SPANS_COUPLING and DIABATIC_FIRST.")
         if (modulo(n-1,nts)==0) then
           dtdia = dt_dyn*min(nts,n_max-(n-1))
-          call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dtdia, OS%MOM_CSp, &
+          call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dtdia, OS%MOM_CSp, &
                         Waves=OS%Waves, do_dynamics=.false., do_thermodynamics=.true., &
                         start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling)
         endif
 
-        call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dt_dyn, OS%MOM_CSp, &
+        call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_dyn, OS%MOM_CSp, &
                       Waves=OS%Waves, do_dynamics=.true., do_thermodynamics=.false., &
                       start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling)
       else
-        call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dt_dyn, OS%MOM_CSp, &
+        call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_dyn, OS%MOM_CSp, &
                       Waves=OS%Waves, do_dynamics=.true., do_thermodynamics=.false., &
                       start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling)
 
@@ -630,7 +629,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
         if (step_thermo) then
           ! Back up Time1 to the start of the thermodynamic segment.
           Time1 = Time1 - real_to_time(dtdia - dt_dyn)
-          call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dtdia, OS%MOM_CSp, &
+          call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dtdia, OS%MOM_CSp, &
                         Waves=OS%Waves, do_dynamics=.false., do_thermodynamics=.true., &
                         start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling)
         endif

diff --git a/config_src/drivers/mct_cap/mom_ocean_model_mct.F90 b/config_src/drivers/mct_cap/mom_ocean_model_mct.F90
@@ -47,7 +47,7 @@ module MOM_ocean_model_mct
 use MOM_tracer_flow_control,  only : call_tracer_register, tracer_flow_control_init
 use MOM_tracer_flow_control,  only : call_tracer_flux_init
 use MOM_unit_scaling,         only : unit_scale_type
-use MOM_variables,           only : surface, stochastic_pattern
+use MOM_variables,            only : surface
 use MOM_verticalGrid,         only : verticalGrid_type
 use MOM_ice_shelf,            only : initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS
 use MOM_ice_shelf,            only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart
@@ -187,7 +187,6 @@ module MOM_ocean_model_mct
                               !! timesteps are taken per thermodynamic step.
   type(surface)   :: sfc_state !< A structure containing pointers to
                               !! the ocean surface state fields.
-  type(stochastic_pattern) :: stochastics !< A structure containing pointers to
   type(ocean_grid_type), pointer :: &
     grid => NULL()            !< A pointer to a grid structure containing metrics
                               !! and related information.
@@ -587,12 +586,12 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
 
   elseif ((.not.do_thermo) .or. (.not.do_dyn)) then
     ! The call sequence is being orchestrated from outside of update_ocean_model.
-    call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dt_coupling, OS%MOM_CSp, &
+    call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp, &
                   Waves=OS%Waves, do_dynamics=do_thermo, do_thermodynamics=do_dyn, &
                   reset_therm=Ocn_fluxes_used)
 
   elseif (OS%single_step_call) then
-    call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves)
+    call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves)
 
   else
     n_max = 1 ; if (dt_coupling > OS%dt) n_max = ceiling(dt_coupling/OS%dt - 0.001)
@@ -616,16 +615,16 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
             "THERMO_SPANS_COUPLING and DIABATIC_FIRST.")
         if (modulo(n-1,nts)==0) then
           dtdia = dt_dyn*min(nts,n_max-(n-1))
-          call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time2, dtdia, OS%MOM_CSp, &
+          call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dtdia, OS%MOM_CSp, &
                         Waves=OS%Waves, do_dynamics=.false., do_thermodynamics=.true., &
                         start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling)
         endif
 
-        call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time2, dt_dyn, OS%MOM_CSp, &
+        call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dt_dyn, OS%MOM_CSp, &
                       Waves=OS%Waves, do_dynamics=.true., do_thermodynamics=.false., &
                       start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling)
       else
-        call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time2, dt_dyn, OS%MOM_CSp, &
+        call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dt_dyn, OS%MOM_CSp, &
                       Waves=OS%Waves, do_dynamics=.true., do_thermodynamics=.false., &
                       start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling)
 
@@ -642,7 +641,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
         if (step_thermo) then
           ! Back up Time2 to the start of the thermodynamic segment.
           Time2 = Time2 - set_time(int(floor((dtdia - dt_dyn) + 0.5)))
-          call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time2, dtdia, OS%MOM_CSp, &
+          call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dtdia, OS%MOM_CSp, &
                         Waves=OS%Waves, do_dynamics=.false., do_thermodynamics=.true., &
                         start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling)
         endif

diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90
@@ -177,10 +177,8 @@ module MOM_ocean_model_nuopc
                               !! steps can span multiple coupled time steps.
   logical :: diabatic_first   !< If true, apply diabatic and thermodynamic
                               !! processes before time stepping the dynamics.
-  logical :: do_sppt         !< If true, stochastically perturb the diabatic and
-                             !! write restarts
-  logical :: pert_epbl       !< If true, then randomly perturb the KE dissipation and
-                             !! genration termsand write restarts
+  logical :: do_sppt          !< If true, allocate array for SPPT
+  logical :: pert_epbl        !< If true, allocate arrays for energetic PBL perturbations
 
   real :: eps_omesh           !< Max allowable difference between ESMF mesh and MOM6
                               !! domain coordinates
@@ -435,20 +433,38 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
 
   endif
 
-  num_procs=num_PEs()
-  allocate(pelist(num_procs))
-  call Get_PElist(pelist,commID = mom_comm)
-  me=PE_here()
-  master=root_PE()
-
-  call init_stochastic_physics_ocn(OS%dt_therm,OS%grid%geoLonT,OS%grid%geoLatT,OS%grid%ied-OS%grid%isd+1,OS%grid%jed-OS%grid%jsd+1,OS%grid%ke,&
-                                   OS%stochastics%pert_epbl,OS%stochastics%do_sppt,master,mom_comm,iret)
-  print*,'after init_stochastic_physics_ocn',OS%stochastics%pert_epbl,OS%stochastics%do_sppt
+! get number of processors and PE list for stocasthci physics initialization
+  call get_param(param_file, mdl, "DO_SPPT", OS%do_sppt, &
+                 "If true, then stochastically perturb the thermodynamic "//&
+                 "tendemcies of T,S, amd h.  Amplitude and correlations are "//&
+                 "controlled by the nam_stoch namelist in the UFS model only.", &
+                 default=.false.)
+  call get_param(param_file, mdl, "PERT_EPBL", OS%pert_epbl, &
+                 "If true, then stochastically perturb the kinetic energy "//&
+                 "production and dissipation terms.  Amplitude and correlations are "//&
+                 "controlled by the nam_stoch namelist in the UFS model only.", &
+                 default=.false.)
+  if (OS%do_sppt .OR. OS%pert_epbl) then
+     num_procs=num_PEs()
+     allocate(pelist(num_procs))
+     call Get_PElist(pelist,commID = mom_comm)
+     me=PE_here()
+     master=root_PE()
+
+     call init_stochastic_physics_ocn(OS%dt_therm,OS%grid%geoLonT,OS%grid%geoLatT,OS%grid%ied-OS%grid%isd+1,OS%grid%jed-OS%grid%jsd+1,OS%grid%ke,&
+                                      OS%pert_epbl,OS%do_sppt,master,mom_comm,iret)
+     if (iret/=0)  then
+         write(6,*) 'call to init_stochastic_physics_ocn failed'
+         call MOM_error(FATAL, "stochastic physics in enambled in MOM6 but "// &
+                    "not activated in stochastic_physics namelist ")
+         return
+     endif
 
-  if (OS%stochastics%do_sppt) allocate(OS%stochastics%sppt_wts(OS%grid%isd:OS%grid%ied,OS%grid%jsd:OS%grid%jed))
-  if (OS%stochastics%pert_epbl) then
-    allocate(OS%stochastics%t_rp1(OS%grid%isd:OS%grid%ied,OS%grid%jsd:OS%grid%jed))
-    allocate(OS%stochastics%t_rp2(OS%grid%isd:OS%grid%ied,OS%grid%jsd:OS%grid%jed))
+     if (OS%do_sppt) allocate(OS%stochastics%sppt_wts(OS%grid%isd:OS%grid%ied,OS%grid%jsd:OS%grid%jed))
+     if (OS%pert_epbl) then
+       allocate(OS%stochastics%t_rp1(OS%grid%isd:OS%grid%ied,OS%grid%jsd:OS%grid%jed))
+       allocate(OS%stochastics%t_rp2(OS%grid%isd:OS%grid%ied,OS%grid%jsd:OS%grid%jed))
+     endif
   endif
   call close_param_file(param_file)
   call diag_mediator_close_registration(OS%diag)
@@ -622,22 +638,22 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
   Master_time = OS%Time ; Time1 = OS%Time
 
 ! update stochastic physics patterns before running next time-step
-  print*,'before call to stoch',OS%stochastics%do_sppt .OR.  OS%stochastics%pert_epbl
-  if (OS%stochastics%do_sppt .OR. OS%stochastics%pert_epbl ) then
+  if (OS%do_sppt .OR. OS%pert_epbl ) then
    call run_stochastic_physics_ocn(OS%stochastics%sppt_wts,OS%stochastics%t_rp1,OS%stochastics%t_rp2)
   endif
 
   if (OS%offline_tracer_mode) then
     call step_offline(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp)
   elseif ((.not.do_thermo) .or. (.not.do_dyn)) then
     ! The call sequence is being orchestrated from outside of update_ocean_model.
-    call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dt_coupling, OS%MOM_CSp, &
+    call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp, &
                   Waves=OS%Waves, do_dynamics=do_thermo, do_thermodynamics=do_dyn, &
-                  reset_therm=Ocn_fluxes_used)
+                  reset_therm=Ocn_fluxes_used, stochastics=OS%stochastics)
  !### What to do with these?   , start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling)
 
   elseif (OS%single_step_call) then
-    call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves)
+    call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves, &
+                  stochastics=OS%stochastics)
   else
     n_max = 1 ; if (dt_coupling > OS%dt) n_max = ceiling(dt_coupling/OS%dt - 0.001)
     dt_dyn = dt_coupling / real(n_max)
@@ -660,18 +676,21 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
             "THERMO_SPANS_COUPLING and DIABATIC_FIRST.")
         if (modulo(n-1,nts)==0) then
           dtdia = dt_dyn*min(nts,n_max-(n-1))
-          call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time2, dtdia, OS%MOM_CSp, &
+          call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dtdia, OS%MOM_CSp, &
                         Waves=OS%Waves, do_dynamics=.false., do_thermodynamics=.true., &
-                        start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling)
+                        start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling, &
+                        stochastics=OS%stochastics)
         endif
 
-        call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time2, dt_dyn, OS%MOM_CSp, &
+        call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dt_dyn, OS%MOM_CSp, &
                       Waves=OS%Waves, do_dynamics=.true., do_thermodynamics=.false., &
-                      start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling)
+                      start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling, &
+                      stochastics=OS%stochastics)
       else
-        call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time2, dt_dyn, OS%MOM_CSp, &
+        call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dt_dyn, OS%MOM_CSp, &
                       Waves=OS%Waves, do_dynamics=.true., do_thermodynamics=.false., &
-                      start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling)
+                      start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling, &
+                      stochastics=OS%stochastics)
 
         step_thermo = .false.
         if (thermo_does_span_coupling) then
@@ -686,9 +705,10 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
         if (step_thermo) then
           ! Back up Time2 to the start of the thermodynamic segment.
           Time2 = Time2 - set_time(int(floor((dtdia - dt_dyn) + 0.5)))
-          call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time2, dtdia, OS%MOM_CSp, &
+          call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dtdia, OS%MOM_CSp, &
                         Waves=OS%Waves, do_dynamics=.false., do_thermodynamics=.true., &
-                        start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling)
+                        start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling, &
+                        stochastics=OS%stochastics)
         endif
       endif
 

diff --git a/config_src/drivers/solo_driver/MOM_driver.F90 b/config_src/drivers/solo_driver/MOM_driver.F90
@@ -62,7 +62,7 @@ program MOM_main
   use MOM_time_manager,    only : JULIAN, GREGORIAN, NOLEAP, THIRTY_DAY_MONTHS, NO_CALENDAR
   use MOM_tracer_flow_control, only : tracer_flow_control_CS
   use MOM_unit_scaling,    only : unit_scale_type
-  use MOM_variables,           only : surface, stochastic_pattern
+  use MOM_variables,       only : surface
   use MOM_verticalGrid,    only : verticalGrid_type
   use MOM_wave_interface,  only : wave_parameters_CS, MOM_wave_interface_init
   use MOM_wave_interface,  only : Update_Surface_Waves
@@ -91,7 +91,6 @@ program MOM_main
   ! A structure containing pointers to the thermodynamic forcing fields
   ! at the ocean surface.
   type(forcing) :: fluxes
-  type(stochastic_pattern) :: stochastics !< A structure containing pointers to
 
   ! A structure containing pointers to the ocean surface state fields.
   type(surface) :: sfc_state
@@ -489,7 +488,7 @@ program MOM_main
     if (offline_tracer_mode) then
       call step_offline(forces, fluxes, sfc_state, Time1, dt_forcing, MOM_CSp)
     elseif (single_step_call) then
-      call step_MOM(forces, fluxes, sfc_state, stochastics, Time1, dt_forcing, MOM_CSp, Waves=Waves_CSP)
+      call step_MOM(forces, fluxes, sfc_state, Time1, dt_forcing, MOM_CSp, Waves=Waves_CSP)
     else
       n_max = 1 ; if (dt_forcing > dt) n_max = ceiling(dt_forcing/dt - 0.001)
       dt_dyn = dt_forcing / real(n_max)
@@ -502,16 +501,16 @@ program MOM_main
         if (diabatic_first) then
           if (modulo(n-1,nts)==0) then
             dtdia = dt_dyn*min(ntstep,n_max-(n-1))
-            call step_MOM(forces, fluxes, sfc_state, stochastics, Time2, dtdia, MOM_CSp, &
+            call step_MOM(forces, fluxes, sfc_state, Time2, dtdia, MOM_CSp, &
                           do_dynamics=.false., do_thermodynamics=.true., &
                           start_cycle=(n==1), end_cycle=.false., cycle_length=dt_forcing)
           endif
 
-          call step_MOM(forces, fluxes, sfc_state, stochastics, Time2, dt_dyn, MOM_CSp, &
+          call step_MOM(forces, fluxes, sfc_state, Time2, dt_dyn, MOM_CSp, &
                         do_dynamics=.true., do_thermodynamics=.false., &
                         start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_forcing)
         else
-          call step_MOM(forces, fluxes, sfc_state, stochastics, Time2, dt_dyn, MOM_CSp, &
+          call step_MOM(forces, fluxes, sfc_state, Time2, dt_dyn, MOM_CSp, &
                         do_dynamics=.true., do_thermodynamics=.false., &
                         start_cycle=(n==1), end_cycle=.false., cycle_length=dt_forcing)
 
@@ -520,7 +519,7 @@ program MOM_main
             ! Back up Time2 to the start of the thermodynamic segment.
             if (n > n_last_thermo+1) &
               Time2 = Time2 - real_to_time(dtdia - dt_dyn)
-            call step_MOM(forces, fluxes, sfc_state, stochastics, Time2, dtdia, MOM_CSp, &
+            call step_MOM(forces, fluxes, sfc_state, Time2, dtdia, MOM_CSp, &
                           do_dynamics=.false., do_thermodynamics=.true., &
                           start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_forcing)
             n_last_thermo = n