Initial stochastic physics implementation

config_src/nuopc_driver/mom_cap.F90

@@ -92,6 +92,7 @@ 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 get_stochy_pattern_mod, only: write_stoch_restart_ocn
 
 implicit none; private
 
@@ -1378,6 +1379,8 @@ 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: ")
@@ -1585,6 +1588,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_stoch_restart_ocn('RESTART/'//trim(restartname))
+        endif
      endif
 
      if (is_root_pe()) then

config_src/nuopc_driver/mom_ocean_model_nuopc.F90

@@ -51,7 +51,7 @@ module MOM_ocean_model_nuopc
 use coupler_types_mod,       only : coupler_type_spawn, coupler_type_write_chksums
 use coupler_types_mod,       only : coupler_type_initialized, coupler_type_copy_data
 use coupler_types_mod,       only : coupler_type_set_diags, coupler_type_send_data
-use mpp_domains_mod,         only : domain2d, mpp_get_layout, mpp_get_global_domain
+use mpp_domains_mod,         only : domain2d, mpp_get_layout, mpp_get_global_domain,mpp_get_pelist
 use mpp_domains_mod,         only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain
 use fms_mod,                 only : stdout
 use mpp_mod,                 only : mpp_chksum
@@ -62,6 +62,9 @@ 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,PE_here,num_PEs
+use MOM_coms,                  only : Get_PElist
+use stochastic_physics,        only : init_stochastic_physics_ocn, run_stochastic_physics_ocn
 
 #include <MOM_memory.h>
 
@@ -174,6 +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,public :: do_sppt   !< If true, allocate array for SPPT
+  logical,public :: pert_epbl !< If true, allocate arrays for energetic PBL perturbations
 
   real :: eps_omesh           !< Max allowable difference between ESMF mesh and MOM6
                               !! domain coordinates
@@ -248,6 +253,13 @@ 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"
@@ -416,6 +428,40 @@ 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
+  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 mpp_get_pelist(pelist, commID=mom_comm)
+     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%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)
 
@@ -585,6 +631,11 @@ 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
@@ -647,6 +698,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
           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
       endif
 

src/core/MOM_forcing_type.F90

@@ -166,6 +166,10 @@ 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]

src/framework/MOM_domains.F90

@@ -3,8 +3,9 @@ module MOM_domains
 
 ! This file is part of MOM6. See LICENSE.md for the license.
 
+
 use MOM_array_transform, only : rotate_array
-use MOM_coms, only : PE_here, root_PE, num_PEs, MOM_infra_init, MOM_infra_end
+use MOM_coms, only : PE_here, root_PE, num_PEs, MOM_infra_init, MOM_infra_end, Get_PElist
 use MOM_coms, only : broadcast, sum_across_PEs, min_across_PEs, max_across_PEs
 use MOM_cpu_clock, only : cpu_clock_begin, cpu_clock_end
 use MOM_error_handler, only : MOM_error, MOM_mesg, NOTE, WARNING, FATAL, is_root_pe

src/parameterizations/vertical/MOM_diabatic_driver.F90

@@ -141,6 +141,8 @@ module MOM_diabatic_driver
                                      !! layers at the bottom into the interior as though
                                      !! a diffusivity of Kd_min_tr (see below) were
                                      !! operating.
+  logical :: do_sppt                 !< If true, stochastically perturb the diabatic
+                                     !! tendencies with a number between 0 and 2
   real    :: Kd_BBL_tr               !< A bottom boundary layer tracer diffusivity that
                                      !! will allow for explicitly specified bottom fluxes
                                      !! [Z2 T-1 ~> m2 s-1].  The entrainment at the bottom is at
@@ -175,7 +177,7 @@ module MOM_diabatic_driver
   integer :: id_Kd_heat  = -1, id_Kd_salt  = -1, id_Kd_interface = -1, id_Kd_ePBL  = -1
   integer :: id_Tdif     = -1, id_Tadv     = -1, id_Sdif         = -1, id_Sadv     = -1
   integer :: id_MLD_003  = -1, id_MLD_0125  = -1, id_MLD_user     = -1, id_mlotstsq = -1
-  integer :: id_subMLN2  = -1
+  integer :: id_subMLN2  = -1, id_sppt_wts  = -1
 
   ! diagnostic for fields prior to applying diapycnal physics
   integer :: id_u_predia = -1, id_v_predia = -1, id_h_predia = -1
@@ -287,8 +289,28 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
   integer :: i, j, k, m, is, ie, js, je, nz
   logical :: showCallTree ! If true, show the call tree
 
+  real, allocatable, dimension(:,:,:)    :: h_in  ! thickness before thermodynamics
+  real, allocatable, dimension(:,:,:)    :: t_in  ! temperature before thermodynamics
+  real, allocatable, dimension(:,:,:)    :: s_in  ! salinity before thermodynamics
+  real :: t_tend,s_tend,h_tend                  ! holder for tendencey needed for SPPT
+  real :: t_pert,s_pert,h_pert                  ! holder for perturbations needed for SPPT
+
   if (G%ke == 1) return
 
+   ! save  copy of the date for SPPT
+  if (CS%do_sppt) then
+    allocate(h_in(G%isd:G%ied, G%jsd:G%jed,G%ke))
+    allocate(t_in(G%isd:G%ied, G%jsd:G%jed,G%ke))
+    allocate(s_in(G%isd:G%ied, G%jsd:G%jed,G%ke))
+    h_in(:,:,:)=h(:,:,:)
+    t_in(:,:,:)=tv%T(:,:,:)
+    s_in(:,:,:)=tv%S(:,:,:)
+
+    if (CS%id_sppt_wts > 0) then
+      call post_data(CS%id_sppt_wts, fluxes%sppt_wts, CS%diag)
+    endif
+  endif
+
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
 
   if (.not. associated(CS)) call MOM_error(FATAL, "MOM_diabatic_driver: "// &
@@ -435,6 +457,34 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
 
   if (CS%debugConservation) call MOM_state_stats('leaving diabatic', u, v, h, tv%T, tv%S, G, GV, US)
 
+  if (CS%do_sppt) then
+  ! perturb diabatic tendecies
+    do k=1,nz
+      do j=js,je
+        do i=is,ie
+          h_tend = (h(i,j,k)-h_in(i,j,k))*fluxes%sppt_wts(i,j)
+          t_tend = (tv%T(i,j,k)-t_in(i,j,k))*fluxes%sppt_wts(i,j)
+          s_tend = (tv%S(i,j,k)-s_in(i,j,k))*fluxes%sppt_wts(i,j)
+          h_pert=h_tend+h_in(i,j,k)
+          t_pert=t_tend+t_in(i,j,k)
+          s_pert=s_tend+s_in(i,j,k)
+          if (h_pert > GV%Angstrom_H) then
+            h(i,j,k) = h_pert
+          else
+            h(i,j,k) = GV%Angstrom_H
+          endif
+          tv%T(i,j,k) = t_pert
+          if (s_pert > 0.0) then
+            tv%S(i,j,k) = s_pert
+          endif
+        enddo
+      enddo
+    enddo
+    deallocate(h_in)
+    deallocate(t_in)
+    deallocate(s_in)
+  endif
+
 end subroutine diabatic
 
 
@@ -453,7 +503,6 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim
                                                                         !! unused have NULL ptrs
   real, dimension(:,:),                      pointer       :: Hml       !< Active mixed layer depth [Z ~> m]
   type(forcing),                             intent(inout) :: fluxes    !< points to forcing fields
-                                                                        !! unused fields have NULL ptrs
   type(vertvisc_type),                       intent(inout) :: visc      !< vertical viscosities, BBL properies, and
   type(accel_diag_ptrs),                     intent(inout) :: ADp       !< related points to accelerations in momentum
                                                                         !! equations, to enable the later derived
@@ -837,7 +886,8 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim
 
     call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
     call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, US, &
-                       CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves)
+                      CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, &
+                      waves=waves)
 
     if (associated(Hml)) then
       call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, Hml(:,:), G, US)
@@ -1568,7 +1618,8 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end,
 
     call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US)
     call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, US, &
-                       CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves)
+                       CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, &
+                       waves=waves)
 
     if (associated(Hml)) then
       call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, Hml(:,:), G, US)
@@ -3357,6 +3408,11 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
                  "mass loss is passed down through the column.", &
                  units="nondim", default=0.8)
 
+  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.)
 
   ! Register all available diagnostics for this module.
   thickness_units = get_thickness_units(GV)
@@ -3382,6 +3438,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
       'Zonal Acceleration from Diapycnal Mixing', 'm s-2', conversion=US%L_T2_to_m_s2)
   CS%id_dvdt_dia = register_diag_field('ocean_model', 'dvdt_dia', diag%axesCvL, Time, &
       'Meridional Acceleration from Diapycnal Mixing', 'm s-2', conversion=US%L_T2_to_m_s2)
+  CS%id_sppt_wts = register_diag_field('ocean_model', 'sppt_pattern', diag%axesT1, Time, &
+     'random pattern for sppt', 'None')
 
   if (CS%use_int_tides) then
     CS%id_cg1 = register_diag_field('ocean_model', 'cn1', diag%axesT1, &