diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 0bad9638f4..652f9e5b47 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -98,7 +98,6 @@ module MOM_cap_mod 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 @@ -1526,7 +1525,7 @@ subroutine ModelAdvance(gcomp, rc) integer :: nc type(ESMF_Time) :: MyTime integer :: seconds, day, year, month, hour, minute - character(ESMF_MAXSTR) :: restartname, cvalue + character(ESMF_MAXSTR) :: restartname, cvalue, stoch_restartname character(240) :: msgString character(ESMF_MAXSTR) :: casename integer :: iostat @@ -1740,26 +1739,19 @@ subroutine ModelAdvance(gcomp, rc) ! write the final restart without a timestamp if (ESMF_AlarmIsRinging(stop_alarm, rc=rc)) then write(restartname,'(A)')"MOM.res" + write(stoch_restartname,'(A)')"ocn_stoch.res.nc" else write(restartname,'(A,I4.4,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2)') & "MOM.res.", year, month, day, hour, minute, seconds + write(stoch_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 restart : "//trim(restartname), ESMF_LOGMSG_INFO) ! write restart file(s) - call ocean_model_restart(ocean_state, restartname=restartname) + call ocean_model_restart(ocean_state, restartname=restartname, & + stoch_restartname=stoch_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 diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index 1064c13d85..448f23140e 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -62,8 +62,7 @@ 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 get_stochy_pattern_mod, only : write_stoch_restart_ocn use iso_fortran_env, only : int64 #include @@ -178,8 +177,10 @@ 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 + logical :: do_sppt !< If true, stochastically perturb the diabatic and + !! write restarts + logical :: pert_epbl !< If true, then randomly perturb the KE dissipation and + !! genration termsand write restarts real :: eps_omesh !< Max allowable difference between ESMF mesh and MOM6 !! domain coordinates @@ -428,10 +429,11 @@ 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 extract_surface_state(OS%MOM_CSp, OS%sfc_state) +! 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 "//& + "tendencies of T,S, and 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, & @@ -439,7 +441,6 @@ 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.) - call extract_surface_state(OS%MOM_CSp, OS%sfc_state) call close_param_file(param_file) call diag_mediator_close_registration(OS%diag) @@ -701,7 +702,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & end subroutine update_ocean_model !> This subroutine writes out the ocean model restart file. -subroutine ocean_model_restart(OS, timestamp, restartname, num_rest_files) +subroutine ocean_model_restart(OS, timestamp, restartname, stoch_restartname, num_rest_files) type(ocean_state_type), pointer :: OS !< A pointer to the structure containing the !! internal ocean state being saved to a restart file character(len=*), optional, intent(in) :: timestamp !< An optional timestamp string that should be @@ -709,6 +710,9 @@ subroutine ocean_model_restart(OS, timestamp, restartname, num_rest_files) character(len=*), optional, intent(in) :: restartname !< Name of restart file to use !! This option distinguishes the cesm interface from the !! non-cesm interface + character(len=*), optional, intent(in) :: stoch_restartname !< Name of restart file to use + !! This option distinguishes the cesm interface from the + !! non-cesm interface integer, optional, intent(out) :: num_rest_files !< number of restart files written if (.not.MOM_state_is_synchronized(OS%MOM_CSp)) & @@ -748,6 +752,11 @@ subroutine ocean_model_restart(OS, timestamp, restartname, num_rest_files) endif endif endif + if (present(stoch_restartname)) then + if (OS%do_sppt .OR. OS%pert_epbl) then + call write_stoch_restart_ocn('RESTART/'//trim(stoch_restartname)) + endif + endif end subroutine ocean_model_restart ! NAME="ocean_model_restart" diff --git a/config_src/external/stochastic_physics/stochastic_physics.F90 b/config_src/external/stochastic_physics/stochastic_physics.F90 new file mode 100644 index 0000000000..df62aa1591 --- /dev/null +++ b/config_src/external/stochastic_physics/stochastic_physics.F90 @@ -0,0 +1,68 @@ +! The are stubs for ocean stochastic physics +! the fully functional code is available at +! http://github.com/noaa-psd/stochastic_physics +module stochastic_physics + +implicit none + +private + +public :: init_stochastic_physics_ocn +public :: run_stochastic_physics_ocn + +contains + +!!!!!!!!!!!!!!!!!!!! +subroutine init_stochastic_physics_ocn(delt,geoLonT,geoLatT,nx,ny,nz,pert_epbl_in,do_sppt_in, & + mpiroot, mpicomm, iret) +implicit none +real,intent(in) :: delt !< timestep in seconds between calls to run_stochastic_physics_ocn +integer,intent(in) :: nx !< number of gridpoints in the x-direction of the compute grid +integer,intent(in) :: ny !< number of gridpoints in the y-direction of the compute grid +integer,intent(in) :: nz !< number of gridpoints in the z-direction of the compute grid +real,intent(in) :: geoLonT(nx,ny) !< Longitude in degrees +real,intent(in) :: geoLatT(nx,ny) !< Latitude in degrees +logical,intent(in) :: pert_epbl_in !< logical flag, if true generate random pattern for ePBL perturbations +logical,intent(in) :: do_sppt_in !< logical flag, if true generate random pattern for SPPT perturbations +integer,intent(in) :: mpiroot !< root processor +integer,intent(in) :: mpicomm !< mpi communicator +integer, intent(out) :: iret !< return code + +iret=0 +if (pert_epbl_in .EQV. .true. ) then + print*,'pert_epbl needs to be false if using the stub' + iret=-1 +endif +if (do_sppt_in.EQV. .true. ) then + print*,'do_sppt needs to be false if using the stub' + iret=-1 +endif +return +end subroutine init_stochastic_physics_ocn + +subroutine run_stochastic_physics_ocn(sppt_wts,t_rp1,t_rp2) +implicit none +real, intent(inout) :: sppt_wts(:,:) !< array containing random weights for SPPT range [0,2] +real, intent(inout) :: t_rp1(:,:) !< array containing random weights for ePBL + !! perturbations (KE generation) range [0,2] +real, intent(inout) :: t_rp2(:,:) !< array containing random weights for ePBL + !! perturbations (KE dissipation) range [0,2] +return +end subroutine run_stochastic_physics_ocn + +end module stochastic_physics + +module get_stochy_pattern_mod + +private + +public :: write_stoch_restart_ocn + +contains +subroutine write_stoch_restart_ocn(sfile) + +character(len=*) :: sfile !< name of restart file +return +end subroutine write_stoch_restart_ocn + +end module get_stochy_pattern_mod diff --git a/src/parameterizations/stochastic/MOM_stochastics.F90 b/src/parameterizations/stochastic/MOM_stochastics.F90 new file mode 100644 index 0000000000..21a22a222e --- /dev/null +++ b/src/parameterizations/stochastic/MOM_stochastics.F90 @@ -0,0 +1,144 @@ +!> 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 + +#include + +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 + integer :: id_sppt_wts = -1 !< Diagnostic id for SPPT + integer :: id_epbl1_wts=-1 !< Diagnostic id for epbl generation perturbation + integer :: id_epbl2_wts=-1 !< Diagnostic id for epbl dissipation perturbation + ! 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 !< stochastic control structure + 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 :: pe_zero ! 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) + pe_zero=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,pe_zero,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 + +end module MOM_stochastics + diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 3828e9483b..351f55984f 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -600,8 +600,8 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs optional, pointer :: Waves !< Wave CS for Langmuir turbulence type(ocean_grid_type), & optional, intent(inout) :: G !< The ocean's grid structure. - real, optional, intent(in) :: epbl1_wt ! random number to perturb KE generation - real, optional, intent(in) :: epbl2_wt ! random number to perturb KE dissipation + real, optional, intent(in) :: epbl1_wt !< random number to perturb KE generation + real, optional, intent(in) :: epbl2_wt !< random number to perturb KE dissipation integer, optional, intent(in) :: i !< The i-index to work on (used for Waves) integer, optional, intent(in) :: j !< The i-index to work on (used for Waves)