diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index ee498f4184..652f9e5b47 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -97,6 +97,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 omp_lib , only : omp_set_num_threads implicit none; private @@ -1524,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 @@ -1738,14 +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) + 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 aab909e56e..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,6 +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 get_stochy_pattern_mod, only : write_stoch_restart_ocn use iso_fortran_env, only : int64 #include @@ -176,6 +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 :: 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 @@ -425,6 +430,17 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i endif 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 "//& + "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, & + "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 close_param_file(param_file) call diag_mediator_close_registration(OS%diag) @@ -686,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 @@ -694,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)) & @@ -733,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/core/MOM.F90 b/src/core/MOM.F90 index 1f1492f2e5..4072cf54a0 100644 --- a/src/core/MOM.F90 +++ b/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 @@ -396,6 +397,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 @@ -640,6 +642,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) & @@ -790,6 +794,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) @@ -1342,7 +1347,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)") @@ -2834,6 +2839,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call init_oda(Time, G, GV, CS%diag, 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) 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_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index a546bcdec0..c9df559583 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -68,6 +68,7 @@ module MOM_diabatic_driver use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units use MOM_wave_speed, only : wave_speeds use MOM_wave_interface, only : wave_parameters_CS +use MOM_stochastics, only : stochastic_CS implicit none ; private @@ -263,7 +264,7 @@ module MOM_diabatic_driver !> This subroutine imposes the diapycnal mass fluxes and the !! accompanying diapycnal advection of momentum and tracers. subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & - G, GV, US, CS, OBC, Waves) + G, GV, US, CS, stoch_CS, OBC, Waves) type(ocean_grid_type), intent(inout) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1] @@ -283,6 +284,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & type(time_type), intent(in) :: Time_end !< Time at the end of the interval type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(diabatic_CS), pointer :: CS !< module control structure + type(stochastic_CS), pointer :: stoch_CS !< stochastic control structure type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. type(Wave_parameters_CS), optional, pointer :: Waves !< Surface gravity waves @@ -295,6 +297,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 active + if (stoch_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 (stoch_CS%id_sppt_wts > 0) then + call post_data(stoch_CS%id_sppt_wts, stoch_CS%sppt_wts, CS%diag) + endif + endif + if (GV%ke == 1) return is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -378,10 +402,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (CS%useALEalgorithm .and. CS%use_legacy_diabatic) then call diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & - G, GV, US, CS, Waves) + G, GV, US, CS, stoch_CS, Waves) elseif (CS%useALEalgorithm) then call diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & - G, GV, US, CS, Waves) + G, GV, US, CS, stoch_CS, Waves) else call layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & G, GV, US, CS, Waves) @@ -447,13 +471,41 @@ 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 (stoch_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))*stoch_CS%sppt_wts(i,j) + t_tend = (tv%T(i,j,k)-t_in(i,j,k))*stoch_CS%sppt_wts(i,j) + s_tend = (tv%S(i,j,k)-s_in(i,j,k))*stoch_CS%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 !> Applies diabatic forcing and diapycnal mixing of temperature, salinity and other tracers for use !! with an ALE algorithm. This version uses an older set of algorithms compared with diabatic_ALE. subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & - G, GV, US, CS, Waves) + G, GV, US, CS, stoch_CS, Waves) type(ocean_grid_type), intent(inout) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -473,6 +525,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim real, intent(in) :: dt !< time increment [T ~> s] type(time_type), intent(in) :: Time_end !< Time at the end of the interval type(diabatic_CS), pointer :: CS !< module control structure + type(stochastic_CS), pointer :: stoch_CS !< stochastic control structure type(Wave_parameters_CS), optional, pointer :: Waves !< Surface gravity waves ! local variables @@ -774,7 +827,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, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, & + waves=waves) if (associated(Hml)) then call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, Hml(:,:), G, US) @@ -1037,7 +1091,7 @@ end subroutine diabatic_ALE_legacy !> This subroutine imposes the diapycnal mass fluxes and the !! accompanying diapycnal advection of momentum and tracers. subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & - G, GV, US, CS, Waves) + G, GV, US, CS, stoch_CS, Waves) type(ocean_grid_type), intent(inout) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -1057,6 +1111,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, real, intent(in) :: dt !< time increment [T ~> s] type(time_type), intent(in) :: Time_end !< Time at the end of the interval type(diabatic_CS), pointer :: CS !< module control structure + type(stochastic_CS), pointer :: stoch_CS !< stochastic control structure type(Wave_parameters_CS), optional, pointer :: Waves !< Surface gravity waves ! local variables @@ -1309,7 +1364,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, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, & + waves=waves) if (associated(Hml)) then call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, Hml(:,:), G, US) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 5a9e67bfd9..351f55984f 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -17,6 +17,7 @@ module MOM_energetic_PBL use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use MOM_wave_interface, only: wave_parameters_CS, Get_Langmuir_Number +use MOM_stochastics, only : stochastic_CS implicit none ; private @@ -168,7 +169,6 @@ module MOM_energetic_PBL real, allocatable, dimension(:,:) :: & ML_depth !< The mixed layer depth determined by active mixing in ePBL [Z ~> m]. - ! These are terms in the mixed layer TKE budget, all in [R Z3 T-3 ~> W m-2 = kg s-3]. real, allocatable, dimension(:,:) :: & diag_TKE_wind, & !< The wind source of TKE [R Z3 T-3 ~> W m-2]. @@ -245,9 +245,9 @@ module MOM_energetic_PBL !! mixed layer model. It assumes that heating, cooling and freshwater fluxes !! have already been applied. All calculations are done implicitly, and there !! is no stability limit on the time step. -subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS, & - dSV_dT, dSV_dS, TKE_forced, buoy_flux, dt_diag, last_call, & - dT_expected, dS_expected, Waves ) +subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS, & + stoch_CS, dSV_dT, dSV_dS, TKE_forced, buoy_flux, dt_diag, & + last_call, dT_expected, dS_expected, Waves) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -300,6 +300,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS !! diffusivities are applied [ppt]. type(wave_parameters_CS), & optional, pointer :: Waves !< Wave CS + type(stochastic_CS), pointer :: stoch_CS !< The control structure returned by a previous ! This subroutine determines the diffusivities from the integrated energetics ! mixed layer model. It assumes that heating, cooling and freshwater fluxes @@ -456,11 +457,17 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS ! Perhaps provide a first guess for MLD based on a stored previous value. MLD_io = -1.0 if (CS%MLD_iteration_guess .and. (CS%ML_Depth(i,j) > 0.0)) MLD_io = CS%ML_Depth(i,j) - - call ePBL_column(h, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, TKE_forcing, B_flux, absf, & - u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, & - US, CS, eCD, dt_diag=dt_diag, Waves=Waves, G=G, i=i, j=j) - + if (stoch_CS%pert_epbl) then ! stochastics are active + call ePBL_column(h, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, TKE_forcing, B_flux, absf, & + u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, & + US, CS, eCD, dt_diag=dt_diag, Waves=Waves, G=G, & + epbl1_wt=stoch_CS%epbl1_wts(i,j),epbl2_wt=stoch_CS%epbl2_wts(i,j), & + i=i, j=j) + else + call ePBL_column(h, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, TKE_forcing, B_flux, absf, & + u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, & + US, CS, eCD, dt_diag=dt_diag, Waves=Waves, G=G, i=i, j=j) + endif ! Copy the diffusivities to a 2-d array. do K=1,nz+1 @@ -531,6 +538,12 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS if (CS%id_LA > 0) call post_data(CS%id_LA, CS%LA, CS%diag) if (CS%id_LA_MOD > 0) call post_data(CS%id_LA_MOD, CS%LA_MOD, CS%diag) if (CS%id_MSTAR_LT > 0) call post_data(CS%id_MSTAR_LT, CS%MSTAR_LT, CS%diag) + ! only write random patterns if running with stochastic physics, otherwise the + ! array is unallocated and will give an error + if (stoch_CS%pert_epbl) then + if (stoch_CS%id_epbl1_wts > 0) call post_data(stoch_CS%id_epbl1_wts, stoch_CS%epbl1_wts, CS%diag) + if (stoch_CS%id_epbl2_wts > 0) call post_data(stoch_CS%id_epbl2_wts, stoch_CS%epbl2_wts, CS%diag) + endif endif if (debug) deallocate(eCD%dT_expect, eCD%dS_expect) @@ -543,7 +556,7 @@ end subroutine energetic_PBL !! mixed layer model for a single column of water. subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, absf, & u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, US, CS, eCD, & - dt_diag, Waves, G, i, j) + dt_diag, Waves, G, epbl1_wt, epbl2_wt, i, j) type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. @@ -587,6 +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 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) @@ -881,6 +896,8 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs else mech_TKE = MSTAR_total * (dt*GV%Rho0* u_star**3) endif + ! stochastically pertrub mech_TKE in the UFS + if (present(epbl1_wt)) mech_TKE=mech_TKE*epbl1_wt if (CS%TKE_diagnostics) then eCD%dTKE_conv = 0.0 ; eCD%dTKE_mixing = 0.0 @@ -963,7 +980,12 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs if (Idecay_len_TKE > 0.0) exp_kh = exp(-h(k-1)*Idecay_len_TKE) if (CS%TKE_diagnostics) & eCD%dTKE_mech_decay = eCD%dTKE_mech_decay + (exp_kh-1.0) * mech_TKE * I_dtdiag - mech_TKE = mech_TKE * exp_kh + if (present(epbl2_wt)) then ! perturb the TKE destruction + mech_TKE = mech_TKE * (1+(exp_kh-1) * epbl2_wt) + else + mech_TKE = mech_TKE * exp_kh + endif + !if ( i .eq. 10 .and. j .eq. 10 .and. k .eq. nz) print*,'mech TKE', mech_TKE ! Accumulate any convectively released potential energy to contribute ! to wstar and to drive penetrating convection. @@ -2373,7 +2395,6 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) Time, 'Velocity Scale that is used.', 'm s-1', conversion=US%Z_to_m*US%s_to_T) CS%id_MSTAR_mix = register_diag_field('ocean_model', 'MSTAR', diag%axesT1, & Time, 'Total mstar that is used.', 'nondim') - if (CS%use_LT) then CS%id_LA = register_diag_field('ocean_model', 'LA', diag%axesT1, & Time, 'Langmuir number.', 'nondim')