Skip to content

Commit

Permalink
stochastic physics re-write
Browse files Browse the repository at this point in the history
  • Loading branch information
pjpegion committed Jul 26, 2021
1 parent 202cbd4 commit e4bc007
Show file tree
Hide file tree
Showing 8 changed files with 264 additions and 90 deletions.
2 changes: 0 additions & 2 deletions config_src/drivers/nuopc_cap/mom_cap.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1527,8 +1527,6 @@ subroutine ModelAdvance(gcomp, rc)
character(len=*),parameter :: subname='(MOM_cap:ModelAdvance)'
character(len=8) :: suffix
integer :: num_rest_files
logical :: do_sppt = .false.
logical :: pert_epbl = .false.

rc = ESMF_SUCCESS
if(profile_memory) call ESMF_VMLogMemInfo("Entering MOM Model_ADVANCE: ")
Expand Down
41 changes: 4 additions & 37 deletions config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,6 @@ module MOM_ocean_model_nuopc
use MOM_surface_forcing_nuopc, only : forcing_save_restart
use MOM_domains, only : root_PE,num_PEs
use MOM_coms, only : Get_PElist
use stochastic_physics, only : init_stochastic_physics_ocn, run_stochastic_physics_ocn

#include <MOM_memory.h>

Expand Down Expand Up @@ -177,8 +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,public :: do_sppt !< If true, allocate array for SPPT
logical,public :: pert_epbl !< If true, allocate arrays for energetic PBL perturbations
logical,public :: do_sppt !< If true, write stochastic physics restarts
logical,public :: pert_epbl !< If true, write stochastic physics restarts

real :: eps_omesh !< Max allowable difference between ESMF mesh and MOM6
!! domain coordinates
Expand Down Expand Up @@ -253,13 +252,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
!! min(HFrz, OBLD), where OBLD is the boundary layer depth.
!! If HFrz <= 0 (default), melt potential will not be computed.
logical :: use_melt_pot!< If true, allocate melt_potential array
! stochastic physics
integer,allocatable :: pelist(:) ! list of pes for this instance of the ocean
integer :: mom_comm ! list of pes for this instance of the ocean
integer :: num_procs ! number of processors to pass to stochastic physics
integer :: iret ! return code from stochastic physics
integer :: me ! my pe
integer :: master ! root pe

! This include declares and sets the variable "version".
#include "version_variable.h"
Expand Down Expand Up @@ -428,7 +420,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i

endif

! get number of processors and PE list for stocasthci physics initialization
! check to see if stochastic physics is active
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 "//&
Expand All @@ -439,27 +431,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
"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)
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%fluxes%sppt_wts(OS%grid%isd:OS%grid%ied,OS%grid%jsd:OS%grid%jed))
if (OS%pert_epbl) then
allocate(OS%fluxes%epbl1_wts(OS%grid%isd:OS%grid%ied,OS%grid%jsd:OS%grid%jed))
allocate(OS%fluxes%epbl2_wts(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 @@ -629,11 +601,6 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
call disable_averaging(OS%diag)
Master_time = OS%Time ; Time1 = OS%Time

! update stochastic physics patterns before running next time-step
if (OS%do_sppt .OR. OS%pert_epbl ) then
call run_stochastic_physics_ocn(OS%fluxes%sppt_wts,OS%fluxes%epbl1_wts,OS%fluxes%epbl2_wts)
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
Expand Down
13 changes: 11 additions & 2 deletions src/core/MOM.F90
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,7 @@ module MOM
use MOM_coord_initialization, only : MOM_initialize_coord
use MOM_diabatic_driver, only : diabatic, diabatic_driver_init, diabatic_CS, extract_diabatic_member
use MOM_diabatic_driver, only : adiabatic, adiabatic_driver_init, diabatic_driver_end
use MOM_stochastics, only : stochastics_init, update_stochastics, stochastic_CS
use MOM_diagnostics, only : calculate_diagnostic_fields, MOM_diagnostics_init
use MOM_diagnostics, only : register_transport_diags, post_transport_diagnostics
use MOM_diagnostics, only : register_surface_diags, write_static_fields
Expand Down Expand Up @@ -380,6 +381,7 @@ module MOM
type(ODA_CS), pointer :: odaCS => NULL() !< a pointer to the control structure for handling
!! ensemble model state vectors and data assimilation
!! increments and priors
type(stochastic_CS), pointer :: stoch_CS => NULL() !< a pointer to the stochastics control structure
end type MOM_control_struct

public initialize_MOM, finish_MOM_initialization, MOM_end
Expand Down Expand Up @@ -624,6 +626,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
call disable_averaging(CS%diag)
endif
endif

if (CS%stoch_CS%do_sppt .OR. CS%stoch_CS%pert_epbl) call update_stochastics(CS%stoch_CS)

if (do_dyn) then
if (G%nonblocking_updates) &
Expand Down Expand Up @@ -774,6 +778,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
enddo ; enddo
endif


call step_MOM_dynamics(forces, CS%p_surf_begin, CS%p_surf_end, dt, &
dt_therm_here, bbl_time_int, CS, &
Time_local, Waves=Waves)
Expand Down Expand Up @@ -1305,7 +1310,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
call cpu_clock_begin(id_clock_diabatic)

call diabatic(u, v, h, tv, CS%Hml, fluxes, CS%visc, CS%ADp, CS%CDp, dtdia, &
Time_end_thermo, G, GV, US, CS%diabatic_CSp, OBC=CS%OBC, Waves=Waves)
Time_end_thermo, G, GV, US, CS%diabatic_CSp, CS%stoch_CS,OBC=CS%OBC, Waves=Waves)
fluxes%fluxes_used = .true.

if (showCallTree) call callTree_waypoint("finished diabatic (step_MOM_thermo)")
Expand Down Expand Up @@ -2370,7 +2375,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &

if (.not. CS%rotate_index) &
G => G_in

! initialize stochastic physics
!call stochastics_init(CS%dt_therm, CS%G, CS%GV, CS%stoch_CS, param_file, diag, Time)
! Set a few remaining fields that are specific to the ocean grid type.
call set_first_direction(G, first_direction)
! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes.
Expand Down Expand Up @@ -2799,6 +2805,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
call init_oda(Time, G, GV, CS%odaCS)
endif

! initialize stochastic physics
call stochastics_init(CS%dt_therm, CS%G, CS%GV, CS%stoch_CS, param_file, diag, Time)

!### This could perhaps go here instead of in finish_MOM_initialization?
! call fix_restart_scaling(GV)
! call fix_restart_unit_scaling(US)
Expand Down
4 changes: 0 additions & 4 deletions src/core/MOM_forcing_type.F90
Original file line number Diff line number Diff line change
Expand Up @@ -164,10 +164,6 @@ module MOM_forcing_type
!! exactly 0 away from shelves or on land.
real, pointer, dimension(:,:) :: iceshelf_melt => NULL() !< Ice shelf melt rate (positive)
!! or freezing (negative) [R Z T-1 ~> kg m-2 s-1]
! stochastic patterns
real, allocatable :: sppt_wts(:,:) !< Random pattern for ocean SPPT
real, allocatable :: epbl1_wts(:,:) !< Random pattern for K.E. generation
real, allocatable :: epbl2_wts(:,:) !< Random pattern for K.E. dissipation

! Scalars set by surface forcing modules
real :: vPrecGlobalAdj = 0. !< adjustment to restoring vprec to zero out global net [kg m-2 s-1]
Expand Down
154 changes: 154 additions & 0 deletions src/parameterizations/stochastic/MOM_stochastics.F90
Original file line number Diff line number Diff line change
@@ -0,0 +1,154 @@
!> Top-level module for the MOM6 ocean model in coupled mode.
module MOM_stochastics

! This file is part of MOM6. See LICENSE.md for the license.

! This is the top level module for the MOM6 ocean model. It contains routines
! for initialization, termination and update of ocean model state. This
! particular version wraps all of the calls for MOM6 in the calls that had
! been used for MOM4.
!
! This code is a stop-gap wrapper of the MOM6 code to enable it to be called
! in the same way as MOM4.

use MOM_diag_mediator, only : register_diag_field, diag_ctrl, time_type
use MOM_grid, only : ocean_grid_type
use MOM_verticalGrid, only : verticalGrid_type
use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe
use MOM_error_handler, only : callTree_enter, callTree_leave
use MOM_file_parser, only : get_param, log_version, close_param_file, param_file_type
use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain
use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain
use MOM_domains, only : root_PE,num_PEs
use MOM_coms, only : Get_PElist
use stochastic_physics, only : init_stochastic_physics_ocn, run_stochastic_physics_ocn

#include <MOM_memory.h>

implicit none ; private

public stochastics_init, update_stochastics

type, public:: stochastic_CS
logical :: do_sppt !< If true, stochastically perturb the diabatic
logical :: pert_epbl !! If true, then randomly perturb the KE dissipation and genration terms
integer :: id_sppt_wts = -1
integer :: id_epbl1_wts=-1,id_epbl2_wts=-1
! stochastic patterns
real, allocatable :: sppt_wts(:,:) !< Random pattern for ocean SPPT
!! tendencies with a number between 0 and 2
real, allocatable :: epbl1_wts(:,:) !< Random pattern for K.E. generation
real, allocatable :: epbl2_wts(:,:) !< Random pattern for K.E. dissipation
type(diag_ctrl), pointer :: diag !< structure used to regulate timing of diagnostic output
type(time_type), pointer :: Time !< Pointer to model time (needed for sponges)
end type stochastic_CS

!> This type is used for communication with other components via the FMS coupler.
!! The element names and types can be changed only with great deliberation, hence
!! the persistnce of things like the cutsy element name "avg_kount".
contains

!> ocean_model_init initializes the ocean model, including registering fields
!! for restarts and reading restart files if appropriate.
!!
!! This subroutine initializes both the ocean state and the ocean surface type.
!! Because of the way that indicies and domains are handled, Ocean_sfc must have
!! been used in a previous call to cean_type.
subroutine stochastics_init(dt, grid, GV, CS, param_file, diag, Time)
real, intent(in) :: dt !< time step [T ~> s]
type(ocean_grid_type), intent(in) :: grid ! horizontal grid information
type(verticalGrid_type), intent(in) :: GV ! vertical grid structure
type(stochastic_CS), pointer, intent(inout):: CS
type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
type(diag_ctrl), target, intent(inout) :: diag !< structure to regulate diagnostic output
type(time_type), target :: Time !< model time
! Local variables
integer,allocatable :: pelist(:) ! list of pes for this instance of the ocean
integer :: mom_comm ! list of pes for this instance of the ocean
integer :: num_procs ! number of processors to pass to stochastic physics
integer :: iret ! return code from stochastic physics
integer :: me ! my pe
integer :: master ! root pe
integer :: nx ! number of x-points including halo
integer :: ny ! number of x-points including halo

! This include declares and sets the variable "version".
#include "version_variable.h"
character(len=40) :: mdl = "ocean_stochastics_init" ! This module's name.

call callTree_enter("ocean_model_stochastic_init(), MOM_stochastics.F90")
if (associated(CS)) then
call MOM_error(WARNING, "MOM_stochastics_init called with an "// &
"associated control structure.")
return
else ; allocate(CS) ; endif

CS%diag => diag
CS%Time => Time

! Read all relevant parameters and write them to the model log.
call log_version(param_file, mdl, version, "")

! get number of processors and PE list for stocasthci physics initialization
call get_param(param_file, mdl, "DO_SPPT", CS%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", CS%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 (CS%do_sppt .OR. CS%pert_epbl) then
num_procs=num_PEs()
allocate(pelist(num_procs))
call Get_PElist(pelist,commID = mom_comm)
master=root_PE()
nx=grid%ied-grid%isd+1
ny=grid%jed-grid%jsd+1
call init_stochastic_physics_ocn(dt,grid%geoLonT,grid%geoLatT,nx,ny,GV%ke, &
CS%pert_epbl,CS%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 (CS%do_sppt) allocate(CS%sppt_wts(grid%isd:grid%ied,grid%jsd:grid%jed))
if (CS%pert_epbl) then
allocate(CS%epbl1_wts(grid%isd:grid%ied,grid%jsd:grid%jed))
allocate(CS%epbl2_wts(grid%isd:grid%ied,grid%jsd:grid%jed))
endif
endif
CS%id_sppt_wts = register_diag_field('ocean_model', 'sppt_pattern', CS%diag%axesT1, Time, &
'random pattern for sppt', 'None')
CS%id_epbl1_wts = register_diag_field('ocean_model', 'epbl1_wts', CS%diag%axesT1, Time, &
'random pattern for KE generation', 'None')
CS%id_epbl2_wts = register_diag_field('ocean_model', 'epbl2_wts', CS%diag%axesT1, Time, &
'random pattern for KE dissipation', 'None')

if (is_root_pe()) &
write(*,'(/12x,a/)') '=== COMPLETED MOM STOCHASTIC INITIALIZATION ====='

call callTree_leave("ocean_model_init(")
end subroutine stochastics_init

!> update_ocean_model uses the forcing in Ice_ocean_boundary to advance the
!! ocean model's state from the input value of Ocean_state (which must be for
!! time time_start_update) for a time interval of Ocean_coupling_time_step,
!! returning the publicly visible ocean surface properties in Ocean_sfc and
!! storing the new ocean properties in Ocean_state.
subroutine update_stochastics(CS)
type(stochastic_CS), intent(inout) :: CS !< diabatic control structure
call callTree_enter("update_stochastics(), MOM_stochastics.F90")

! update stochastic physics patterns before running next time-step
call run_stochastic_physics_ocn(CS%sppt_wts,CS%epbl1_wts,CS%epbl2_wts)
print*,'in update_stoch',minval(CS%sppt_wts),maxval(CS%sppt_wts),minval(CS%epbl1_wts),maxval(CS%epbl1_wts)

end subroutine update_stochastics

end module MOM_stochastics

64 changes: 64 additions & 0 deletions src/parameterizations/stochastic/MOM_stochastics_stub.F90
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
!> Top-level module for the MOM6 ocean model in coupled mode.
module MOM_stochastics

! This file is part of MOM6. See LICENSE.md for the license.

! This is the top level module for the MOM6 ocean model. It contains routines
! for initialization, termination and update of ocean model state. This
! particular version wraps all of the calls for MOM6 in the calls that had
! been used for MOM4.
!
! This code is a stop-gap wrapper of the MOM6 code to enable it to be called
! in the same way as MOM4.

use MOM_diag_mediator, only : register_diag_field, diag_ctrl, time_type
use MOM_grid, only : ocean_grid_type
use MOM_verticalGrid, only : verticalGrid_type
use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe
use MOM_error_handler, only : callTree_enter, callTree_leave
use MOM_file_parser, only : get_param, log_version, close_param_file, param_file_type
use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain
use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain
use MOM_domains, only : root_PE,num_PEs
use MOM_coms, only : Get_PElist

#include <MOM_memory.h>

implicit none ; private

public stochastics_init, update_stochastics

type, public:: stochastic_CS
logical :: do_sppt !< If true, stochastically perturb the diabatic
logical :: pert_epbl !! If true, then randomly perturb the KE dissipation and genration terms
integer :: id_sppt_wts = -1
integer :: id_epbl1_wts=-1,id_epbl2_wts=-1
! stochastic patterns
real, allocatable :: sppt_wts(:,:) !< Random pattern for ocean SPPT
!! tendencies with a number between 0 and 2
real, allocatable :: epbl1_wts(:,:) !< Random pattern for K.E. generation
real, allocatable :: epbl2_wts(:,:) !< Random pattern for K.E. dissipation
type(diag_ctrl), pointer :: diag !< structure used to regulate timing of diagnostic output
type(time_type), pointer :: Time !< Pointer to model time (needed for sponges)
end type stochastic_CS

contains

subroutine stochastics_init(dt, grid, GV, CS, param_file, diag, Time)
real, intent(in) :: dt !< time step [T ~> s]
type(ocean_grid_type), intent(in) :: grid ! horizontal grid information
type(verticalGrid_type), intent(in) :: GV ! vertical grid structure
type(stochastic_CS), pointer, intent(inout):: CS
type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
type(diag_ctrl), target, intent(inout) :: diag !< structure to regulate diagnostic output
type(time_type), target :: Time !< model time
return
end subroutine stochastics_init

subroutine update_stochastics(CS)
type(stochastic_CS), intent(inout) :: CS !< diabatic control structure
return
end subroutine update_stochastics

end module MOM_stochastics

Loading

0 comments on commit e4bc007

Please sign in to comment.