Merge pull request mom-ocean#72 from pjpegion/ocn_stoch_july2021

Ocean stochastic physics

config_src/drivers/nuopc_cap/mom_cap.F90

@@ -97,6 +97,8 @@ module MOM_cap_mod
 use NUOPC_Model, only: model_label_SetRunClock    => label_SetRunClock
 use NUOPC_Model, only: model_label_Finalize       => label_Finalize
 use NUOPC_Model, only: SetVM
+
+use MOM_stochastics,         only : write_mom_restart_stoch
 !$use omp_lib             , only : omp_set_num_threads
 
 implicit none; private
@@ -1746,6 +1748,18 @@ subroutine ModelAdvance(gcomp, rc)
 
         ! write restart file(s)
         call ocean_model_restart(ocean_state, restartname=restartname)
+
+        if (ocean_state%do_sppt .OR. ocean_state%pert_epbl) then
+          if (ESMF_AlarmIsRinging(stop_alarm, rc=rc)) then
+             write(restartname,'(A)')"ocn_stoch.res.nc"
+          else
+             write(restartname,'(A,I4.4,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,A)') &
+                  "ocn_stoch.res.", year, month, day, hour, minute, seconds,".nc"
+          endif
+          call ESMF_LogWrite("MOM_cap: Writing stoch restart :  "//trim(restartname), &
+                             ESMF_LOGMSG_INFO)
+          call write_mom_restart_stoch('RESTART/'//trim(restartname))
+        endif
      endif
 
      if (is_root_pe()) then

config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90

@@ -62,6 +62,8 @@ module MOM_ocean_model_nuopc
 use MOM_surface_forcing_nuopc, only : convert_IOB_to_forces, ice_ocn_bnd_type_chksum
 use MOM_surface_forcing_nuopc, only : ice_ocean_boundary_type, surface_forcing_CS
 use MOM_surface_forcing_nuopc, only : forcing_save_restart
+use MOM_domains,               only : root_PE,num_PEs
+use MOM_coms,                  only : Get_PElist
 use iso_fortran_env,           only : int64
 
 #include <MOM_memory.h>
@@ -176,6 +178,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, 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
@@ -424,6 +428,17 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
 
   endif
 
+  ! 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 "//&
+                 "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.)
   call extract_surface_state(OS%MOM_CSp, OS%sfc_state)
 
   call close_param_file(param_file)

config_src/external/OCEAN_stochastic_phyiscs/MOM_stochastics.F90

@@ -0,0 +1,155 @@
+!> 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, update, and writing restart of stochastic physics. This
+! particular version wraps all of the calls for MOM6 in the calls that had
+! been used for 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
+use get_stochy_pattern_mod, only: write_stoch_restart_ocn
+
+#include <MOM_memory.h>
+
+implicit none ; private
+
+public stochastics_init, update_stochastics, write_mom_restart_stoch
+
+!> This control structure holds parameters for the MOM_stochastics module
+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
+
+!!   This subroutine initializes the stochastics physics control structure.
+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
+         call MOM_error(FATAL, "call to init_stochastic_physics_ocn failed")
+         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(")
+  return
+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)
+
+  return
+end subroutine update_stochastics
+
+!< wrapper to write ocean stochastic restarts
+subroutine write_mom_restart_stoch(filename)
+  character(len=*) :: filename
+
+  call callTree_enter("write_mom_restart_stoch(), MOM_stochastics.F90")
+
+  call write_stoch_restart_ocn(filename)
+
+  return
+end subroutine write_mom_restart_stoch
+
+end module MOM_stochastics
+

config_src/external/OCEAN_stochastic_phyiscs/MOM_stochastics_stub.F90

@@ -0,0 +1,68 @@
+!> 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
+! placeholder for initialization, update, and writing restarts of ocean stochastic physics.
+! The actualy stochastic physics is available at
+! https://github.com/ufs-community/ufs-weather-model
+! 
+
+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
+
+!> This control structure holds parameters for the MOM_stochastics module
+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
+  !>@{ Diagnostic IDs
+  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
+subroutine write_mom_restart_stoch(filename)
+  character(len=*) :: filename
+  return
+end subroutine write_mom_restart_stoch
+end module MOM_stochastics
+

src/core/MOM.F90

@@ -59,6 +59,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
@@ -391,6 +392,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
@@ -635,6 +637,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
       call disable_averaging(CS%diag)
     endif
   endif
+  ! advance the random pattern if stochastic physics is active
+  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) &
@@ -785,6 +789,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)
@@ -1316,7 +1321,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)")
@@ -2814,6 +2819,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)