Skip to content

Commit

Permalink
make stochastics optional
Browse files Browse the repository at this point in the history
  • Loading branch information
@pjpegion
pjpegion committed Feb 2, 2021
1 parent 85023f8 commit 0b99c1f
Show file tree
Hide file tree
Showing 8 changed files with 145 additions and 103 deletions.
17 changes: 8 additions & 9 deletions config_src/coupled_driver/ocean_model_MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,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
Expand Down Expand Up @@ -187,7 +187,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.
Expand Down Expand Up @@ -566,7 +565,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.
Expand All @@ -581,12 +580,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)
Expand All @@ -608,16 +607,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)

Expand All @@ -634,7 +633,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
Expand Down
15 changes: 7 additions & 8 deletions config_src/mct_driver/mom_ocean_model_mct.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,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
Expand Down Expand Up @@ -185,7 +185,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.
Expand Down Expand Up @@ -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)
Expand All @@ -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)

Expand All @@ -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
Expand Down
78 changes: 51 additions & 27 deletions config_src/nuopc_driver/mom_ocean_model_nuopc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -176,6 +176,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, 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
Expand Down Expand Up @@ -426,20 +428,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

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))
! 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%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)
Expand Down Expand Up @@ -611,22 +631,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)
Expand All @@ -649,18 +669,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
Expand All @@ -675,9 +698,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

Expand Down
13 changes: 6 additions & 7 deletions config_src/solo_driver/MOM_driver.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ program MOM_main
use MOM_time_manager, only : 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_write_cputime, only : write_cputime, MOM_write_cputime_init
use MOM_write_cputime, only : write_cputime_start_clock, write_cputime_CS
Expand All @@ -84,7 +84,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
Expand Down Expand Up @@ -501,7 +500,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)
Expand All @@ -514,16 +513,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)

Expand All @@ -532,7 +531,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
Expand Down
Loading

0 comments on commit 0b99c1f

Please sign in to comment.