additions for stochastic physics and ePBL perts

config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90

@@ -315,6 +315,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
 
     call safe_alloc_ptr(fluxes%p_surf     ,isd,ied,jsd,jed)
     call safe_alloc_ptr(fluxes%p_surf_full,isd,ied,jsd,jed)
+    print*,'allocate fluxes%t_rp'
+    call safe_alloc_ptr(fluxes%t_rp,isd,ied,jsd,jed)
     if (CS%use_limited_P_SSH) then
        fluxes%p_surf_SSH => fluxes%p_surf
     else

src/core/MOM.F90

@@ -30,6 +30,7 @@ module MOM
 use MOM_domains,              only : To_All, Omit_corners, CGRID_NE, SCALAR_PAIR
 use MOM_domains,              only : create_group_pass, do_group_pass, group_pass_type
 use MOM_domains,              only : start_group_pass, complete_group_pass, Omit_Corners
+use MOM_domains,              only : root_PE,PE_here,Get_PElist,num_PEs
 use MOM_error_handler,        only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe
 use MOM_error_handler,        only : MOM_set_verbosity, callTree_showQuery
 use MOM_error_handler,        only : callTree_enter, callTree_leave, callTree_waypoint
@@ -155,8 +156,7 @@ module MOM
 use MOM_offline_main,          only : offline_fw_fluxes_into_ocean, offline_fw_fluxes_out_ocean
 use MOM_offline_main,          only : offline_advection_layer, offline_transport_end
 use MOM_ALE,                   only : ale_offline_tracer_final, ALE_main_offline
-use MOM_ice_shelf,             only : ice_shelf_CS, ice_shelf_query, initialize_ice_shelf
-use MOM_particles_mod,         only : particles, particles_init, particles_run, particles_save_restart, particles_end
+use stochastic_physics,        only : init_stochastic_physics_ocn,run_stochastic_physics_ocn
 
 implicit none ; private
 
@@ -249,6 +249,8 @@ module MOM
   logical :: offline_tracer_mode = .false.
                     !< If true, step_offline() is called instead of step_MOM().
                     !! This is intended for running MOM6 in offline tracer mode
+  logical :: do_stochy = .false.
+                    !< If true, call stochastic physics pattern generator
 
   type(time_type), pointer :: Time   !< pointer to the ocean clock
   real    :: dt                      !< (baroclinic) dynamics time step [T ~> s]
@@ -827,6 +829,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
         enddo ; enddo
       endif
 
+      print*,'calling run_stochastic_physics_ocn',CS%do_stochy
+      if (CS%do_stochy) call run_stochastic_physics_ocn(forces%t_rp)
 
       call step_MOM_dynamics(forces, CS%p_surf_begin, CS%p_surf_end, dt, &
                              dt_therm_here, bbl_time_int, CS, &
@@ -979,7 +983,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
     if (CS%time_in_thermo_cycle > 0.0) then
       call enable_averages(CS%time_in_thermo_cycle, Time_local, CS%diag)
       call post_surface_thermo_diags(CS%sfc_IDs, G, GV, US, CS%diag, CS%time_in_thermo_cycle, &
-                                     sfc_state_diag, CS%tv, ssh, CS%ave_ssh_ibc)
+                                    sfc_state, CS%tv, ssh, fluxes%t_rp, CS%ave_ssh_ibc)
     endif
     call disable_averaging(CS%diag)
     call cpu_clock_end(id_clock_diagnostics)
@@ -1812,10 +1816,12 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   integer :: nkml, nkbl, verbosity, write_geom
   integer :: dynamics_stencil  ! The computational stencil for the calculations
                                ! in the dynamic core.
-  real :: conv2watt            ! A conversion factor from temperature fluxes to heat
-                               ! fluxes [J m-2 H-1 degC-1 ~> J m-3 degC-1 or J kg-1 degC-1]
-  real :: conv2salt            ! A conversion factor for salt fluxes [m H-1 ~> 1] or [kg m-2 H-1 ~> 1]
-  real :: RL2_T2_rescale, Z_rescale, QRZ_rescale ! Unit conversion factors
+  integer,allocatable :: pelist(:) ! list of pes for this instance of the ocean
+  integer             :: num_procs
+!  model
+  integer :: me                !  my pe
+  integer :: master            !  root pe
+  real :: conv2watt, conv2salt
   character(len=48) :: flux_units, S_flux_units
 
   type(vardesc) :: vd_T, vd_S  ! Structures describing temperature and salinity variables.
@@ -2492,6 +2498,25 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
     call callTree_waypoint("returned from ALE_init() (initialize_MOM)")
   endif
 
+  !   Shift from using the temporary dynamic grid type to using the final
+  ! (potentially static) ocean-specific grid type.
+  !   The next line would be needed if G%Domain had not already been init'd above:
+  !     call clone_MOM_domain(dG%Domain, G%Domain)
+  call MOM_grid_init(G, param_file, US, HI, bathymetry_at_vel=bathy_at_vel)
+  call copy_dyngrid_to_MOM_grid(dG, G, US)
+  call destroy_dyn_horgrid(dG)
+
+  num_procs=num_PEs()
+  allocate(pelist(num_procs))
+  call Get_PElist(pelist)
+  me=PE_here()
+  master=root_PE()
+
+  !call init_stochastic_physics_ocn(CS%dt_therm,G,me,master,pelist,CS%do_stochy)
+  print*,'callling init_stochastic_physics_ocn',maxval(G%geoLatT)
+  call init_stochastic_physics_ocn(CS%dt_therm,G%geoLonT,G%geoLatT,G%ied-G%isd+1,G%jed-G%jsd+1,nz,CS%do_stochy)
+  print*,'back from init_stochastic_physics_ocn',CS%do_stochy
+
   ! Set a few remaining fields that are specific to the ocean grid type.
   if (CS%rotate_index) then
     call set_first_direction(G, modulo(first_direction + turns, 2))

src/core/MOM_forcing_type.F90

@@ -145,6 +145,8 @@ module MOM_forcing_type
                 !< Pressure at the top ocean interface [R L2 T-2 ~> Pa] that is used in corrections to the sea surface
                 !! height field that is passed back to the calling routines.
                 !! p_surf_SSH may point to p_surf or to p_surf_full.
+  real, pointer, dimension(:,:) :: t_rp => NULL()
+                !< random pattern at t-points
   logical :: accumulate_p_surf = .false. !< If true, the surface pressure due to the atmosphere
                                  !! and various types of ice needs to be accumulated, and the
                                  !! surface pressure explicitly reset to zero at the driver level
@@ -250,10 +252,10 @@ module MOM_forcing_type
                 !! nondimensional from 0 to 1 [nondim]. This is only associated if ice shelves are enabled,
                 !! and is exactly 0 away from shelves or on land.
   real, pointer, dimension(:,:) :: &
-    rigidity_ice_u => NULL(), & !< Depth-integrated lateral viscosity of ice shelves or sea ice at
-                                !! u-points [L4 Z-1 T-1 ~> m3 s-1]
-    rigidity_ice_v => NULL()    !< Depth-integrated lateral viscosity of ice shelves or sea ice at
-                                !! v-points [L4 Z-1 T-1 ~> m3 s-1]
+    rigidity_ice_u => NULL(), & !< Depth-integrated lateral viscosity of ice shelves or sea ice at u-points [m3 s-1]
+    rigidity_ice_v => NULL()    !< Depth-integrated lateral viscosity of ice shelves or sea ice at v-points [m3 s-1]
+  real, pointer, dimension(:,:) :: t_rp => NULL()
+                                 !< random pattern at t-points
   real :: dt_force_accum = -1.0 !< The amount of time over which the mechanical forcing fluxes
                                 !! have been averaged [s].
   logical :: net_mass_src_set = .false. !< If true, an estimate of net_mass_src has been provided.
@@ -2135,6 +2137,12 @@ subroutine copy_common_forcing_fields(forces, fluxes, G, skip_pres)
 
   do_pres = .true. ; if (present(skip_pres)) do_pres = .not.skip_pres
 
+  if (associated(forces%t_rp) .and. associated(fluxes%t_rp)) then
+    do j=js,je ; do i=is,ie
+      fluxes%t_rp(i,j) = forces%t_rp(i,j)
+    enddo ; enddo
+  endif
+
   if (associated(forces%ustar) .and. associated(fluxes%ustar)) then
     do j=js,je ; do i=is,ie
       fluxes%ustar(i,j) = forces%ustar(i,j)
@@ -3109,6 +3117,7 @@ subroutine allocate_mech_forcing_by_group(G, forces, stress, ustar, shelf, &
   call myAlloc(forces%p_surf,isd,ied,jsd,jed, press)
   call myAlloc(forces%p_surf_full,isd,ied,jsd,jed, press)
   call myAlloc(forces%net_mass_src,isd,ied,jsd,jed, press)
+  call myAlloc(forces%t_rp,isd,ied,jsd,jed, press)
 
   call myAlloc(forces%rigidity_ice_u,IsdB,IedB,jsd,jed, shelf)
   call myAlloc(forces%rigidity_ice_v,isd,ied,JsdB,JedB, shelf)
@@ -3273,6 +3282,7 @@ subroutine deallocate_forcing_type(fluxes)
   if (associated(fluxes%netMassIn))            deallocate(fluxes%netMassIn)
   if (associated(fluxes%salt_flux))            deallocate(fluxes%salt_flux)
   if (associated(fluxes%p_surf_full))          deallocate(fluxes%p_surf_full)
+  if (associated(fluxes%t_rp))                 deallocate(fluxes%t_rp)
   if (associated(fluxes%p_surf))               deallocate(fluxes%p_surf)
   if (associated(fluxes%TKE_tidal))            deallocate(fluxes%TKE_tidal)
   if (associated(fluxes%ustar_tidal))          deallocate(fluxes%ustar_tidal)
@@ -3303,6 +3313,7 @@ subroutine deallocate_mech_forcing(forces)
   if (associated(forces%ustar)) deallocate(forces%ustar)
   if (associated(forces%p_surf))         deallocate(forces%p_surf)
   if (associated(forces%p_surf_full))    deallocate(forces%p_surf_full)
+  if (associated(forces%t_rp))           deallocate(forces%t_rp)
   if (associated(forces%net_mass_src))   deallocate(forces%net_mass_src)
   if (associated(forces%rigidity_ice_u)) deallocate(forces%rigidity_ice_u)
   if (associated(forces%rigidity_ice_v)) deallocate(forces%rigidity_ice_v)

src/diagnostics/MOM_diagnostics.F90

@@ -139,7 +139,9 @@ module MOM_diagnostics
   integer :: id_salt_deficit = -1
   integer :: id_Heat_PmE     = -1
   integer :: id_intern_heat  = -1
-  !>@}
+!  stochastic pattern
+  integer :: id_t_rp         = -1
+  !!@}
 end type surface_diag_IDs
 
 
@@ -1277,7 +1279,7 @@ end subroutine post_surface_dyn_diags
 !> This routine posts diagnostics of various ocean surface and integrated
 !! quantities at the time the ocean state is reported back to the caller
 subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv, &
-                                    ssh, ssh_ibc)
+                                    ssh, t_rp, ssh_ibc)
   type(surface_diag_IDs),   intent(in) :: IDs !< A structure with the diagnostic IDs.
   type(ocean_grid_type),    intent(in) :: G   !< ocean grid structure
   type(verticalGrid_type),  intent(in) :: GV  !< ocean vertical grid structure
@@ -1286,8 +1288,10 @@ subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv
   real,                     intent(in) :: dt_int !< total time step associated with these diagnostics [T ~> s].
   type(surface),            intent(in) :: sfc_state !< structure describing the ocean surface state
   type(thermo_var_ptrs),    intent(in) :: tv  !< A structure pointing to various thermodynamic variables
-  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh !< Time mean surface height without corrections
-                                              !! for ice displacement [Z ~> m]
+  real, dimension(SZI_(G),SZJ_(G)), &
+                            intent(in) :: ssh !< Time mean surface height without corrections for ice displacement [m]
+  real, dimension(SZI_(G),SZJ_(G)), &
+                            intent(in) :: t_rp!< random pattern for stochastic proceeses
   real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh_ibc !< Time mean surface height with corrections
                                               !! for ice displacement and the inverse barometer [Z ~> m]
 
@@ -1409,6 +1413,11 @@ subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv
     call post_data(IDs%id_sss_sq, work_2d, diag, mask=G%mask2dT)
   endif
 
+  if (IDs%id_t_rp > 0) then
+     !call post_data(IDs%id_t_rp, t_rp, diag, mask=G%mask2dT)
+     call post_data(IDs%id_t_rp, t_rp, diag)
+  endif
+
   call coupler_type_send_data(sfc_state%tr_fields, get_diag_time_end(diag))
 
 end subroutine post_surface_thermo_diags
@@ -1896,8 +1905,9 @@ subroutine register_surface_diags(Time, G, US, IDs, diag, tv)
          'Heat flux into ocean from mass flux into ocean', &
          'W m-2', conversion=US%QRZ_T_to_W_m2)
   IDs%id_intern_heat = register_diag_field('ocean_model', 'internal_heat', diag%axesT1, Time,&
-         'Heat flux into ocean from geothermal or other internal sources', &
-         'W m-2', conversion=US%QRZ_T_to_W_m2)
+         'Heat flux into ocean from geothermal or other internal sources', 'W m-2')
+  IDs%id_t_rp = register_diag_field('ocean_model', 'random_pattern', diag%axesT1, Time, &
+      'random pattern for stochastics', 'None')
 
 end subroutine register_surface_diags
 

src/framework/MOM_domains.F90

@@ -3,56 +3,23 @@ module MOM_domains
 
 ! This file is part of MOM6. See LICENSE.md for the license.
 
-use MOM_coms_infra,       only : MOM_infra_init, MOM_infra_end
-use MOM_coms_infra,       only : PE_here, root_PE, num_PEs, broadcast
-use MOM_coms_infra,       only : sum_across_PEs, min_across_PEs, max_across_PEs
-use MOM_domain_infra,     only : MOM_domain_type, domain2D, domain1D, group_pass_type
-use MOM_domain_infra,     only : create_MOM_domain, clone_MOM_domain, deallocate_MOM_domain
-use MOM_domain_infra,     only : get_domain_extent, get_domain_components, same_domain
-use MOM_domain_infra,     only : compute_block_extent, get_global_shape
-use MOM_domain_infra,     only : pass_var, pass_vector, fill_symmetric_edges, global_field_sum
-use MOM_domain_infra,     only : pass_var_start, pass_var_complete
-use MOM_domain_infra,     only : pass_vector_start, pass_vector_complete
-use MOM_domain_infra,     only : create_group_pass, do_group_pass
-use MOM_domain_infra,     only : start_group_pass, complete_group_pass
-use MOM_domain_infra,     only : rescale_comp_data, global_field, redistribute_array, broadcast_domain
-use MOM_domain_infra,     only : MOM_thread_affinity_set, set_MOM_thread_affinity
-use MOM_domain_infra,     only : AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR, BITWISE_EXACT_SUM
-use MOM_domain_infra,     only : CORNER, CENTER, NORTH_FACE, EAST_FACE
-use MOM_domain_infra,     only : To_East, To_West, To_North, To_South, To_All, Omit_Corners
-use MOM_error_handler,    only : MOM_error, MOM_mesg, NOTE, WARNING, FATAL
-use MOM_file_parser,      only : get_param, log_param, log_version, param_file_type
-use MOM_io_infra,         only : file_exists
+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
+use MOM_file_parser, only : get_param, log_param, log_version
+use MOM_file_parser, only : param_file_type
 use MOM_string_functions, only : slasher
 
 implicit none ; private
 
-public :: MOM_infra_init, MOM_infra_end
-!  Domain types and creation and destruction routines
-public :: MOM_domain_type, domain2D, domain1D
-public :: MOM_domains_init, create_MOM_domain, clone_MOM_domain, deallocate_MOM_domain
-public :: MOM_thread_affinity_set, set_MOM_thread_affinity
-!  Domain query routines
-public :: get_domain_extent, get_domain_components, get_global_shape, same_domain
-public :: PE_here, root_PE, num_PEs
-!  Blocks are not actively used in MOM6, so this routine could be deprecated.
-public :: compute_block_extent
-!  Single call communication routines
-public :: pass_var, pass_vector, fill_symmetric_edges, broadcast
-!  Non-blocking communication routines
-public :: pass_var_start, pass_var_complete, pass_vector_start, pass_vector_complete
-!  Multi-variable group communication routines and type
-public :: create_group_pass, do_group_pass, group_pass_type, start_group_pass, complete_group_pass
-!  Global reduction routines
-public :: sum_across_PEs, min_across_PEs, max_across_PEs
-public :: global_field, redistribute_array, broadcast_domain
-!  Simple index-convention-invariant array manipulation routine
-public :: rescale_comp_data
-!> These encoding constants are used to indicate the staggering of scalars and vectors
-public :: AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR
-!> These encoding constants are used to indicate the discretization position of a variable
-public :: CORNER, CENTER, NORTH_FACE, EAST_FACE
-!> These encoding constants indicate communication patterns.  In practice they can be added.
+public :: MOM_domains_init, MOM_infra_init, MOM_infra_end, get_domain_extent, get_domain_extent_dsamp2
+public :: MOM_define_domain, MOM_define_io_domain, clone_MOM_domain
+public :: pass_var, pass_vector, PE_here, root_PE, num_PEs, Get_PElist
+public :: pass_var_start, pass_var_complete, fill_symmetric_edges, broadcast
+public :: pass_vector_start, pass_vector_complete
+public :: global_field_sum, sum_across_PEs, min_across_PEs, max_across_PEs
+public :: AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR, BITWISE_EXACT_SUM, CORNER, CENTER
 public :: To_East, To_West, To_North, To_South, To_All, Omit_Corners
 ! These are no longer used by MOM6 because the reproducing sum works so well, but they are
 ! still referenced by some of the non-GFDL couplers.

src/parameterizations/vertical/MOM_energetic_PBL.F90

@@ -401,7 +401,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS
       do K=1,nz+1 ; Kd(K) = 0.0 ; enddo
 
       ! Make local copies of surface forcing and process them.
-      u_star = fluxes%ustar(i,j)
+      u_star = fluxes%ustar(i,j)*(fluxes%t_rp(i,j))
       u_star_Mean = fluxes%ustar_gustless(i,j)
       B_flux = buoy_flux(i,j)
       if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then