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 1ef5b07c06..9ffe4cd794 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -64,7 +64,6 @@ module MOM_ocean_model_nuopc use MOM_surface_forcing_nuopc, only : forcing_save_restart 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 @@ -178,8 +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 :: do_sppt !< If true, allocate array for SPPT - logical :: pert_epbl !< If true, allocate arrays for energetic PBL perturbations + 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 @@ -255,12 +254,6 @@ 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 :: 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" @@ -432,7 +425,7 @@ 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 + ! 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 "//& @@ -443,27 +436,7 @@ 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.) - if (OS%do_sppt .OR. OS%pert_epbl) then - num_procs=num_PEs() - allocate(pelist(num_procs)) - call Get_PElist(pelist,commID = mom_comm) - 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) @@ -635,11 +608,6 @@ 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 diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 477d185e96..a0594290ee 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -403,16 +403,6 @@ 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(porous_barrier_ptrs) :: pbv !< porous barrier fractional cell metrics - real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) & - :: por_face_areaU !< fractional open area of U-faces [nondim] - real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) & - :: por_face_areaV !< fractional open area of V-faces [nondim] - real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NK_INTERFACE_) & - :: por_layer_widthU !< fractional open width of U-faces [nondim] - real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NK_INTERFACE_) & - :: por_layer_widthV !< fractional open width of V-faces [nondim] - type(particles), pointer :: particles => NULL() ! NULL() !< a pointer to the stochastics control structure end type MOM_control_struct @@ -673,7 +663,7 @@ 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 @@ -825,6 +815,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) @@ -1391,7 +1382,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)") @@ -2500,6 +2491,16 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & else call set_first_direction(G, modulo(first_direction, 2)) endif + call MOM_grid_init(G_in, param_file, US, HI_in, bathymetry_at_vel=bathy_at_vel) + call copy_dyngrid_to_MOM_grid(dG_in, G_in, US) + call destroy_dyn_horgrid(dG_in) + + if (.not. CS%rotate_index) & + G => G_in + ! initialize stochastic physics + !call stochastics_init(CS%dt_therm, CS%G, CS%GV, CS%stoch_CS, param_file, diag, Time) + ! Set a few remaining fields that are specific to the ocean grid type. + call set_first_direction(G, first_direction) ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes. if (CS%debug .or. G%symmetric) then call clone_MOM_domain(G%Domain, G%Domain_aux, symmetric=.false.) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index bdf4823f81..4a702d4c0e 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -170,8 +170,6 @@ 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] - real, pointer, dimension(:,:) :: shelf_sfc_mass_flux => NULL() !< Ice shelf surface mass flux - !! deposition from the atmosphere. [R Z T-1 ~> kg m-2 s-1] ! Scalars set by surface forcing modules real :: vPrecGlobalAdj = 0. !< adjustment to restoring vprec to zero out global net [kg m-2 s-1] diff --git a/src/parameterizations/stochastic/MOM_stochastics.F90 b/src/parameterizations/stochastic/MOM_stochastics.F90 index 21a22a222e..5bcf158f7e 100644 --- a/src/parameterizations/stochastic/MOM_stochastics.F90 +++ b/src/parameterizations/stochastic/MOM_stochastics.F90 @@ -4,10 +4,13 @@ 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 +! for initialization, termination and update of ocean model state. This ! particular version wraps all of the calls for MOM6 in the calls that had ! been used for MOM4. ! +! This code is a stop-gap wrapper of the MOM6 code to enable it to be called +! in the same way as 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 @@ -26,13 +29,11 @@ module MOM_stochastics 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 + 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 @@ -42,14 +43,22 @@ module MOM_stochastics type(time_type), pointer :: Time !< Pointer to model time (needed for sponges) end type stochastic_CS +!> This type is used for communication with other components via the FMS coupler. +!! The element names and types can be changed only with great deliberation, hence +!! the persistnce of things like the cutsy element name "avg_kount". contains -!! This subroutine initializes the stochastics physics control structure. +!> ocean_model_init initializes the ocean model, including registering fields +!! for restarts and reading restart files if appropriate. +!! +!! This subroutine initializes both the ocean state and the ocean surface type. +!! Because of the way that indicies and domains are handled, Ocean_sfc must have +!! been used in a previous call to cean_type. 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(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 @@ -59,7 +68,7 @@ subroutine stochastics_init(dt, grid, GV, CS, param_file, diag, Time) 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 :: master ! root pe integer :: nx ! number of x-points including halo integer :: ny ! number of x-points including halo @@ -95,13 +104,15 @@ subroutine stochastics_init(dt, grid, GV, CS, param_file, diag, Time) 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 + 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,pe_zero,mom_comm,iret) + 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") + 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 @@ -122,7 +133,6 @@ subroutine stochastics_init(dt, grid, GV, CS, param_file, diag, Time) 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 @@ -135,9 +145,9 @@ subroutine update_stochastics(CS) 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) + call run_stochastic_physics_ocn(CS%sppt_wts,CS%epbl1_wts,CS%epbl2_wts) + print*,'in update_stoch',minval(CS%sppt_wts),maxval(CS%sppt_wts),minval(CS%epbl1_wts),maxval(CS%epbl1_wts) - return end subroutine update_stochastics end module MOM_stochastics diff --git a/src/parameterizations/stochastic/MOM_stochastics_stub.F90 b/src/parameterizations/stochastic/MOM_stochastics_stub.F90 new file mode 100644 index 0000000000..f03f5283d3 --- /dev/null +++ b/src/parameterizations/stochastic/MOM_stochastics_stub.F90 @@ -0,0 +1,64 @@ +!> 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, termination and update of ocean model state. This +! particular version wraps all of the calls for MOM6 in the calls that had +! been used for MOM4. +! +! This code is a stop-gap wrapper of the MOM6 code to enable it to be called +! in the same way as 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 + +#include + +implicit none ; private + +public stochastics_init, update_stochastics + +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 + +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 + +end module MOM_stochastics + diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 9ae391815b..319a6f2be6 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -69,6 +69,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 @@ -142,12 +143,12 @@ module MOM_diabatic_driver integer :: NKBL !< The number of buffer layers (if bulk_mixed_layer) logical :: massless_match_targets !< If true (the default), keep the T & S !! consistent with the target values. - logical :: mix_boundary_tracers !< If true, mix the passive tracers in massless - !! 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 + logical :: mix_boundary_tracers !< If true, mix the passive tracers in massless layers at the + !! bottom into the interior as though a diffusivity of + !! Kd_min_tr (see below) were operating. + logical :: mix_boundary_tracer_ALE !< If true, in ALE mode mix the passive tracers in massless + !! layers at the bottom into the interior as though a + !! diffusivity of Kd_min_tr (see below) were operating. 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,15 +176,20 @@ module MOM_diabatic_driver real :: MLD_EN_VALS(3) !< Energy values for energy mixed layer diagnostics !>@{ Diagnostic IDs - integer :: id_cg1 = -1 ! diag handle for mode-1 speed (BDM) - integer, allocatable, dimension(:) :: id_cn ! diag handle for all mode speeds (BDM) - integer :: id_wd = -1, id_ea = -1, id_eb = -1 ! used by layer diabatic - integer :: id_dudt_dia = -1, id_dvdt_dia = -1, id_ea_s = -1, id_eb_s = -1 - integer :: id_ea_t = -1, id_eb_t = -1 - 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, id_sppt_wts = -1 + integer :: id_cg1 = -1 ! diag handle for mode-1 speed + integer, allocatable, dimension(:) :: id_cn ! diag handle for all mode speeds + integer :: id_ea = -1, id_eb = -1 ! used by layer diabatic + integer :: id_ea_t = -1, id_eb_t = -1, id_ea_s = -1, id_eb_s = -1 + integer :: id_Kd_heat = -1, id_Kd_salt = -1, id_Kd_int = -1, id_Kd_ePBL = -1 + integer :: id_Tdif = -1, id_Sdif = -1, id_Tadv = -1, id_Sadv = -1 + ! These are handles to diagnostics related to the mixed layer properties. + integer :: id_MLD_003 = -1, id_MLD_0125 = -1, id_MLD_user = -1, id_mlotstsq = -1 + integer :: id_MLD_EN1 = -1, id_MLD_EN2 = -1, id_MLD_EN3 = -1, id_subMLN2 = -1 + + ! These are handles to diatgnostics that are only available in non-ALE layered mode. + integer :: id_wd = -1 + integer :: id_dudt_dia = -1, id_dvdt_dia = -1 + integer :: id_hf_dudt_dia_2d = -1, id_hf_dvdt_dia_2d = -1 ! diagnostic for fields prior to applying diapycnal physics integer :: id_u_predia = -1, id_v_predia = -1, id_h_predia = -1 @@ -266,13 +272,13 @@ 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) - 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_(G)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1] - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1] - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< thickness [H ~> m or kg m-2] - type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields + 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] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< thickness [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields !! unused have NULL ptrs real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m] type(forcing), intent(inout) :: fluxes !< points to forcing fields @@ -281,13 +287,14 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & type(accel_diag_ptrs), intent(inout) :: ADp !< related points to accelerations in momentum !! equations, to enable the later derived !! diagnostics, like energy budgets - type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations - real, intent(in) :: dt !< time increment [T ~> s] - 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(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. - type(Wave_parameters_CS), optional, pointer :: Waves !< Surface gravity waves + type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations + real, intent(in) :: dt !< time increment [T ~> s] + 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 ! local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & @@ -307,16 +314,16 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (G%ke == 1) return ! save copy of the date for SPPT - if (CS%do_sppt) then + 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(:,:) + 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) + if (stoch_CS%id_sppt_wts > 0) then + call post_data(stoch_CS%id_sppt_wts, stoch_CS%sppt_wts, CS%diag) endif endif @@ -405,10 +412,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) @@ -474,14 +481,14 @@ 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 + 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))*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_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) @@ -507,14 +514,14 @@ 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) - 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 - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1] - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1] - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< thickness [H ~> m or kg m-2] - type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields + 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 + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< thickness [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields !! unused have NULL ptrs real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m] type(forcing), intent(inout) :: fluxes !< points to forcing fields @@ -522,11 +529,12 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim type(accel_diag_ptrs), intent(inout) :: ADp !< related points to accelerations in momentum !! equations, to enable the later derived !! diagnostics, like energy budgets - type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations - 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(Wave_parameters_CS), optional, pointer :: Waves !< Surface gravity waves + type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations + 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 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & @@ -827,7 +835,7 @@ 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, & + CS%energetic_PBL_CSp, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, & waves=waves) if (associated(Hml)) then @@ -1091,14 +1099,14 @@ 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) - 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 - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1] - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1] - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< thickness [H ~> m or kg m-2] - type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields + 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 + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< thickness [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields !! unused have NULL ptrs real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m] type(forcing), intent(inout) :: fluxes !< points to forcing fields @@ -1107,11 +1115,12 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, type(accel_diag_ptrs), intent(inout) :: ADp !< related points to accelerations in momentum !! equations, to enable the later derived !! diagnostics, like energy budgets - type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations - 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(Wave_parameters_CS), optional, pointer :: Waves !< Surface gravity waves + type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations + 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 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & @@ -1363,7 +1372,7 @@ 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, & + CS%energetic_PBL_CSp, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, & waves=waves) if (associated(Hml)) then @@ -3086,11 +3095,9 @@ 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.) + if (CS%use_energetic_PBL .and. .not.CS%useALEalgorithm) & + call MOM_error(FATAL, "diabatic_driver_init: "//& + "ENERGETICS_SFC_PBL = True is only coded to work when USE_REGRIDDING = True.") ! Register all available diagnostics for this module. thickness_units = get_thickness_units(GV) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 7cf8d48ae0..36e92c2850 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -58,7 +58,6 @@ module MOM_energetic_PBL !! self-consistent mixed layer depth. Otherwise use the false position !! after a maximum and minimum bound have been evaluated and the !! returned value from the previous guess or bisection before this. - logical :: pert_epbl !! If true, then randomly perturb the KE dissipation and genration terms integer :: max_MLD_its !< The maximum number of iterations that can be used to find a !! self-consistent mixed layer depth with Use_MLD_iteration. real :: MixLenExponent !< Exponent in the mixing length shape-function. @@ -198,7 +197,6 @@ module MOM_energetic_PBL integer :: id_TKE_mech_decay = -1, id_TKE_conv_decay = -1 integer :: id_Mixing_Length = -1, id_Velocity_Scale = -1 integer :: id_MSTAR_mix = -1, id_LA_mod = -1, id_LA = -1, id_MSTAR_LT = -1 - integer :: id_epbl1_wts=-1,id_epbl2_wts=-1 !>@} end type energetic_PBL_CS @@ -246,9 +244,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 @@ -301,6 +299,8 @@ 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 @@ -440,11 +440,12 @@ 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) - if (CS%pert_epbl) then ! stochastics are active + 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=epbl1_wts(i,j),epbl2_wt=epbl2_wts(i,j),i=i, j=j) + 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, & @@ -508,9 +509,9 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS 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 (CS%pert_epbl) then - if (CS%id_epbl1_wts > 0) call post_data(CS%id_epbl1_wts, stochastics%epbl1_wts, CS%diag) - if (CS%id_epbl2_wts > 0) call post_data(CS%id_epbl2_wts, stochastics%epbl2_wts, CS%diag) + 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 end subroutine energetic_PBL @@ -857,7 +858,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs mech_TKE = MSTAR_total * (dt*GV%Rho0* u_star**3) endif ! stochastically pertrub mech_TKE in the UFS - if (CS%pert_epbl) mech_TKE=mech_TKE*epbl1_wt + 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 @@ -940,7 +941,7 @@ 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 - if (CS%pert_epbl) then ! perturb the TKE destruction + 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 @@ -2174,11 +2175,6 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "This is only used if USE_MLD_ITERATION is True.", & units="nondim", default=2.0) - 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.) !/ Turbulent velocity scale in mixing coefficient call get_param(param_file, mdl, "EPBL_VEL_SCALE_SCHEME", tmpstr, & "Selects the method for translating TKE into turbulent velocities. "//& @@ -2353,10 +2349,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') - CS%id_epbl1_wts = register_diag_field('ocean_model', 'epbl1_wts', diag%axesT1, Time, & - 'random pattern for KE generation', 'None') - CS%id_epbl2_wts = register_diag_field('ocean_model', 'epbl2_wts', diag%axesT1, Time, & - 'random pattern for KE dissipation', 'None') if (CS%use_LT) then CS%id_LA = register_diag_field('ocean_model', 'LA', diag%axesT1, & Time, 'Langmuir number.', 'nondim')