diff --git a/config_src/drivers/FMS_cap/ocean_model_MOM.F90 b/config_src/drivers/FMS_cap/ocean_model_MOM.F90 index 97fb869ad4..6292c32469 100644 --- a/config_src/drivers/FMS_cap/ocean_model_MOM.F90 +++ b/config_src/drivers/FMS_cap/ocean_model_MOM.F90 @@ -50,7 +50,7 @@ module ocean_model_mod use MOM_tracer_flow_control, only : call_tracer_register, tracer_flow_control_init use MOM_tracer_flow_control, only : call_tracer_flux_init use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : surface +use MOM_variables, only : surface, stochastic_pattern use MOM_verticalGrid, only : verticalGrid_type use MOM_ice_shelf, only : initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS use MOM_ice_shelf, only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart @@ -186,6 +186,7 @@ module ocean_model_mod !! timesteps are taken per thermodynamic step. type(surface) :: sfc_state !< A structure containing pointers to !! the ocean surface state fields. + type(stochastic_pattern) :: stochastics !< A structure containing pointers to type(ocean_grid_type), pointer :: & grid => NULL() !< A pointer to a grid structure containing metrics !! and related information. @@ -576,12 +577,12 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda 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 ! The call sequence is being orchestrated from outside of update_ocean_model. - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dt_coupling, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=do_dyn, do_thermodynamics=do_thermo, & start_cycle=start_cycle, end_cycle=end_cycle, cycle_length=cycle_length, & reset_therm=Ocn_fluxes_used) elseif (OS%single_step_call) then - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves) + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves) else ! Step both the dynamics and thermodynamics with separate calls. n_max = 1 ; if (dt_coupling > OS%dt) n_max = ceiling(dt_coupling/OS%dt - 0.001) dt_dyn = dt_coupling / real(n_max) @@ -603,16 +604,16 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda "THERMO_SPANS_COUPLING and DIABATIC_FIRST.") if (modulo(n-1,nts)==0) then dtdia = dt_dyn*min(nts,n_max-(n-1)) - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dtdia, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dtdia, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=.false., do_thermodynamics=.true., & start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling) endif - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_dyn, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dt_dyn, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=.true., do_thermodynamics=.false., & start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling) else - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_dyn, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dt_dyn, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=.true., do_thermodynamics=.false., & start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling) @@ -629,7 +630,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda if (step_thermo) then ! Back up Time1 to the start of the thermodynamic segment. Time1 = Time1 - real_to_time(dtdia - dt_dyn) - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dtdia, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dtdia, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=.false., do_thermodynamics=.true., & start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling) endif diff --git a/config_src/drivers/mct_cap/mom_ocean_model_mct.F90 b/config_src/drivers/mct_cap/mom_ocean_model_mct.F90 index 3bd0e1e28d..4a4f6eee05 100644 --- a/config_src/drivers/mct_cap/mom_ocean_model_mct.F90 +++ b/config_src/drivers/mct_cap/mom_ocean_model_mct.F90 @@ -47,7 +47,7 @@ module MOM_ocean_model_mct use MOM_tracer_flow_control, only : call_tracer_register, tracer_flow_control_init use MOM_tracer_flow_control, only : call_tracer_flux_init use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : surface +use MOM_variables, only : surface, stochastic_pattern use MOM_verticalGrid, only : verticalGrid_type use MOM_ice_shelf, only : initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS use MOM_ice_shelf, only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart @@ -187,6 +187,7 @@ module MOM_ocean_model_mct !! timesteps are taken per thermodynamic step. type(surface) :: sfc_state !< A structure containing pointers to !! the ocean surface state fields. + type(stochastic_pattern) :: stochastics !< A structure containing pointers to type(ocean_grid_type), pointer :: & grid => NULL() !< A pointer to a grid structure containing metrics !! and related information. @@ -586,12 +587,12 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & elseif ((.not.do_thermo) .or. (.not.do_dyn)) then ! The call sequence is being orchestrated from outside of update_ocean_model. - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dt_coupling, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=do_thermo, do_thermodynamics=do_dyn, & reset_therm=Ocn_fluxes_used) elseif (OS%single_step_call) then - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves) + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves) else n_max = 1 ; if (dt_coupling > OS%dt) n_max = ceiling(dt_coupling/OS%dt - 0.001) @@ -615,16 +616,16 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & "THERMO_SPANS_COUPLING and DIABATIC_FIRST.") if (modulo(n-1,nts)==0) then dtdia = dt_dyn*min(nts,n_max-(n-1)) - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dtdia, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time2, dtdia, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=.false., do_thermodynamics=.true., & start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling) endif - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dt_dyn, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time2, dt_dyn, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=.true., do_thermodynamics=.false., & start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling) else - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dt_dyn, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time2, dt_dyn, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=.true., do_thermodynamics=.false., & start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling) @@ -641,7 +642,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & if (step_thermo) then ! Back up Time2 to the start of the thermodynamic segment. Time2 = Time2 - set_time(int(floor((dtdia - dt_dyn) + 0.5))) - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dtdia, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time2, dtdia, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=.false., do_thermodynamics=.true., & start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling) endif diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 8778df067a..3de56c0511 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -97,9 +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 -#ifdef UFS use get_stochy_pattern_mod, only: write_stoch_restart_ocn -#endif !$use omp_lib , only : omp_set_num_threads @@ -1755,7 +1753,6 @@ subroutine ModelAdvance(gcomp, rc) call ocean_model_restart(ocean_state, restartname=restartname) ! write stochastic physics restart file if active -#ifdef UFS if (ESMF_AlarmIsRinging(stop_alarm, rc=rc)) then write(restartname,'(A)')"ocn_stoch.res.nc" else @@ -1764,7 +1761,6 @@ subroutine ModelAdvance(gcomp, rc) endif call ESMF_LogWrite("MOM_cap: Writing restart : "//trim(restartname), ESMF_LOGMSG_INFO) call write_stoch_restart_ocn('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 448f23140e..290e6b30df 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -43,7 +43,7 @@ module MOM_ocean_model_nuopc use time_interp_external_mod,only : time_interp_external_init use MOM_tracer_flow_control, only : call_tracer_flux_init use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : surface +use MOM_variables, only : surface, stochastic_pattern use MOM_verticalGrid, only : verticalGrid_type use MOM_ice_shelf, only : initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS use MOM_ice_shelf, only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart @@ -62,8 +62,8 @@ 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 +use MOM_domains, only : root_PE,PE_here,Get_PElist,num_PEs +use stochastic_physics, only : init_stochastic_physics_ocn, run_stochastic_physics_ocn #include @@ -194,6 +194,7 @@ module MOM_ocean_model_nuopc !! timesteps are taken per thermodynamic step. type(surface) :: sfc_state !< A structure containing pointers to !! the ocean surface state fields. + type(stochastic_pattern) :: stochastics !< A structure containing pointers to type(ocean_grid_type), pointer :: & grid => NULL() !< A pointer to a grid structure containing metrics !! and related information. @@ -255,9 +256,14 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i !! The actual depth over which melt potential is computed will !! 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 - logical :: use_CFC !< If true, allocated arrays for surface CFCs. - + logical :: use_melt_pot!< If true, allocate melt_potential array +! stochastic physics + integer,allocatable :: pelist(:) ! list of pes for this instance of the ocean + integer :: mom_comm ! list of pes for this instance of the ocean + integer :: num_procs ! number of processors to pass to stochastic physics + integer :: iret ! return code from stochastic physics + integer :: me ! my pe + integer :: master ! root pe ! This include declares and sets the variable "version". #include "version_variable.h" @@ -429,19 +435,21 @@ 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.) + num_procs=num_PEs() + allocate(pelist(num_procs)) + call Get_PElist(pelist,commID = mom_comm) + me=PE_here() + master=root_PE() + + call init_stochastic_physics_ocn(OS%dt_therm,OS%grid%geoLonT,OS%grid%geoLatT,OS%grid%ied-OS%grid%isd+1,OS%grid%jed-OS%grid%jsd+1,OS%grid%ke,& + OS%stochastics%pert_epbl,OS%stochastics%do_sppt,master,mom_comm,iret) + print*,'after init_stochastic_physics_ocn',OS%stochastics%pert_epbl,OS%stochastics%do_sppt + if (OS%stochastics%do_sppt) allocate(OS%stochastics%sppt_wts(OS%grid%isd:OS%grid%ied,OS%grid%jsd:OS%grid%jed)) + if (OS%stochastics%pert_epbl) then + allocate(OS%stochastics%t_rp1(OS%grid%isd:OS%grid%ied,OS%grid%jsd:OS%grid%jed)) + allocate(OS%stochastics%t_rp2(OS%grid%isd:OS%grid%ied,OS%grid%jsd:OS%grid%jed)) + endif call close_param_file(param_file) call diag_mediator_close_registration(OS%diag) @@ -613,17 +621,23 @@ 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 + print*,'before call to stoch',OS%stochastics%do_sppt .OR. OS%stochastics%pert_epbl + if (OS%stochastics%do_sppt .OR. OS%stochastics%pert_epbl ) then + call run_stochastic_physics_ocn(OS%stochastics%sppt_wts,OS%stochastics%t_rp1,OS%stochastics%t_rp2) + 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 ! The call sequence is being orchestrated from outside of update_ocean_model. - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dt_coupling, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=do_thermo, do_thermodynamics=do_dyn, & reset_therm=Ocn_fluxes_used) !### What to do with these? , start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling) elseif (OS%single_step_call) then - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves) + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves) else n_max = 1 ; if (dt_coupling > OS%dt) n_max = ceiling(dt_coupling/OS%dt - 0.001) dt_dyn = dt_coupling / real(n_max) @@ -646,16 +660,16 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & "THERMO_SPANS_COUPLING and DIABATIC_FIRST.") if (modulo(n-1,nts)==0) then dtdia = dt_dyn*min(nts,n_max-(n-1)) - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dtdia, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time2, dtdia, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=.false., do_thermodynamics=.true., & start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling) endif - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dt_dyn, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time2, dt_dyn, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=.true., do_thermodynamics=.false., & start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling) else - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dt_dyn, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time2, dt_dyn, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=.true., do_thermodynamics=.false., & start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling) @@ -672,7 +686,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & if (step_thermo) then ! Back up Time2 to the start of the thermodynamic segment. Time2 = Time2 - set_time(int(floor((dtdia - dt_dyn) + 0.5))) - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dtdia, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time2, dtdia, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=.false., do_thermodynamics=.true., & start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling) endif diff --git a/config_src/drivers/solo_driver/MOM_driver.F90 b/config_src/drivers/solo_driver/MOM_driver.F90 index 5bfe0bb8a2..ee5e92e0f4 100644 --- a/config_src/drivers/solo_driver/MOM_driver.F90 +++ b/config_src/drivers/solo_driver/MOM_driver.F90 @@ -62,7 +62,7 @@ program MOM_main use MOM_time_manager, only : JULIAN, GREGORIAN, NOLEAP, THIRTY_DAY_MONTHS, NO_CALENDAR use MOM_tracer_flow_control, only : tracer_flow_control_CS use MOM_unit_scaling, only : unit_scale_type - use MOM_variables, only : surface + use MOM_variables, only : surface, stochastic_pattern use MOM_verticalGrid, only : verticalGrid_type use MOM_wave_interface, only : wave_parameters_CS, MOM_wave_interface_init use MOM_wave_interface, only : Update_Surface_Waves @@ -91,6 +91,8 @@ program MOM_main ! A structure containing pointers to the thermodynamic forcing fields ! at the ocean surface. type(forcing) :: fluxes + type(stochastic_pattern) :: stochastics !< A structure containing pointers to + ! A structure containing pointers to the ocean surface state fields. type(surface) :: sfc_state @@ -487,7 +489,7 @@ program MOM_main if (offline_tracer_mode) then call step_offline(forces, fluxes, sfc_state, Time1, dt_forcing, MOM_CSp) elseif (single_step_call) then - call step_MOM(forces, fluxes, sfc_state, Time1, dt_forcing, MOM_CSp, Waves=Waves_CSP) + call step_MOM(forces, fluxes, sfc_state, stochastics, Time1, dt_forcing, MOM_CSp, Waves=Waves_CSP) else n_max = 1 ; if (dt_forcing > dt) n_max = ceiling(dt_forcing/dt - 0.001) dt_dyn = dt_forcing / real(n_max) @@ -500,16 +502,16 @@ program MOM_main if (diabatic_first) then if (modulo(n-1,nts)==0) then dtdia = dt_dyn*min(ntstep,n_max-(n-1)) - call step_MOM(forces, fluxes, sfc_state, Time2, dtdia, MOM_CSp, & + call step_MOM(forces, fluxes, sfc_state, stochastics, Time2, dtdia, MOM_CSp, & do_dynamics=.false., do_thermodynamics=.true., & start_cycle=(n==1), end_cycle=.false., cycle_length=dt_forcing) endif - call step_MOM(forces, fluxes, sfc_state, Time2, dt_dyn, MOM_CSp, & + call step_MOM(forces, fluxes, sfc_state, stochastics, Time2, dt_dyn, MOM_CSp, & do_dynamics=.true., do_thermodynamics=.false., & start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_forcing) else - call step_MOM(forces, fluxes, sfc_state, Time2, dt_dyn, MOM_CSp, & + call step_MOM(forces, fluxes, sfc_state, stochastics, Time2, dt_dyn, MOM_CSp, & do_dynamics=.true., do_thermodynamics=.false., & start_cycle=(n==1), end_cycle=.false., cycle_length=dt_forcing) @@ -518,7 +520,7 @@ program MOM_main ! Back up Time2 to the start of the thermodynamic segment. if (n > n_last_thermo+1) & Time2 = Time2 - real_to_time(dtdia - dt_dyn) - call step_MOM(forces, fluxes, sfc_state, Time2, dtdia, MOM_CSp, & + call step_MOM(forces, fluxes, sfc_state, stochastics, Time2, dtdia, MOM_CSp, & do_dynamics=.false., do_thermodynamics=.true., & start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_forcing) n_last_thermo = n diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 0c5a1a9788..50e58b8a63 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -30,7 +30,6 @@ 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 @@ -133,7 +132,7 @@ module MOM use MOM_variables, only : surface, allocate_surface_state, deallocate_surface_state use MOM_variables, only : thermo_var_ptrs, vertvisc_type, porous_barrier_ptrs use MOM_variables, only : accel_diag_ptrs, cont_diag_ptrs, ocean_internal_state -use MOM_variables, only : rotate_surface_state +use MOM_variables, only : rotate_surface_state,stochastic_pattern use MOM_verticalGrid, only : verticalGrid_type, verticalGridInit, verticalGridEnd use MOM_verticalGrid, only : fix_restart_scaling use MOM_verticalGrid, only : get_thickness_units, get_flux_units, get_tr_flux_units @@ -155,9 +154,6 @@ 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 -#ifdef UFS -use stochastic_physics, only : init_stochastic_physics_ocn -#endif implicit none ; private @@ -455,13 +451,14 @@ module MOM !! The action of lateral processes on tracers occur in calls to !! advect_tracer and tracer_hordiff. Vertical mixing and possibly remapping !! occur inside of diabatic. -subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS, & +subroutine step_MOM(forces_in, fluxes_in, sfc_state, stochastics, Time_start, time_int_in, CS, & Waves, do_dynamics, do_thermodynamics, start_cycle, & end_cycle, cycle_length, reset_therm) type(mech_forcing), target, intent(inout) :: forces_in !< A structure with the driving mechanical forces type(forcing), target, intent(inout) :: fluxes_in !< A structure with pointers to themodynamic, !! tracer and mass exchange forcing fields type(surface), target, intent(inout) :: sfc_state !< surface ocean state + type(stochastic_pattern), intent(in) :: stochastics !< surface ocean state type(time_type), intent(in) :: Time_start !< starting time of a segment, as a time type real, intent(in) :: time_int_in !< time interval covered by this run segment [s]. type(MOM_control_struct), intent(inout), target :: CS !< control structure from initialize_MOM @@ -770,8 +767,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS endif ! Apply diabatic forcing, do mixing, and regrid. - call step_MOM_thermo(CS, G, GV, US, u, v, h, CS%tv, fluxes, dtdia, & - end_time_thermo, .true., Waves=Waves) + call step_MOM_thermo(CS, G, GV, US, u, v, h, CS%tv, fluxes, stochastics, & + dtdia, end_time_thermo, .true., Waves=Waves) CS%time_in_thermo_cycle = CS%time_in_thermo_cycle + dtdia ! The diabatic processes are now ahead of the dynamics by dtdia. @@ -871,8 +868,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS CS%Time = CS%Time - real_to_time(0.5*US%T_to_s*(dtdia-dt)) ! Apply diabatic forcing, do mixing, and regrid. - call step_MOM_thermo(CS, G, GV, US, u, v, h, CS%tv, fluxes, dtdia, & - Time_local, .false., Waves=Waves) + call step_MOM_thermo(CS, G, GV, US, u, v, h, CS%tv, fluxes, stochastics, & + dtdia, Time_local, .false., Waves=Waves) CS%time_in_thermo_cycle = CS%time_in_thermo_cycle + dtdia if ((CS%t_dyn_rel_thermo==0.0) .and. .not.do_dyn) then @@ -1304,8 +1301,8 @@ end subroutine step_MOM_tracer_dyn !> MOM_step_thermo orchestrates the thermodynamic time stepping and vertical !! remapping, via calls to diabatic (or adiabatic) and ALE_main. -subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & - Time_end_thermo, update_BBL, Waves) +subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, stochastics, & + dtdia, Time_end_thermo, update_BBL, Waves) type(MOM_control_struct), intent(inout) :: CS !< Master MOM control structure type(ocean_grid_type), intent(inout) :: G !< ocean grid structure type(verticalGrid_type), intent(inout) :: GV !< ocean vertical grid structure @@ -1318,6 +1315,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & intent(inout) :: h !< layer thickness [H ~> m or kg m-2] type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic variables type(forcing), intent(inout) :: fluxes !< pointers to forcing fields + type(stochastic_pattern), intent(in) :: stochastics !< surface ocean state real, intent(in) :: dtdia !< The time interval over which to advance [T ~> s] type(time_type), intent(in) :: Time_end_thermo !< End of averaging interval for thermo diags logical, intent(in) :: update_BBL !< If true, calculate the bottom boundary layer properties. @@ -1393,8 +1391,8 @@ 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, CS%stoch_CS, CS%OBC, Waves) + call diabatic(u, v, h, tv, CS%Hml, fluxes, stochastics, CS%visc, CS%ADp, CS%CDp, dtdia, & + Time_end_thermo, G, GV, US, CS%diabatic_CSp, OBC=CS%OBC, Waves=Waves) fluxes%fluxes_used = .true. if (showCallTree) call callTree_waypoint("finished diabatic (step_MOM_thermo)") @@ -1806,19 +1804,12 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & logical :: calc_dtbt ! Indicates whether the dynamically adjusted barotropic ! time step needs to be updated before it is used. logical :: debug_truncations ! If true, turn on diagnostics useful for debugging truncations. - logical :: do_epbl,do_sppt integer :: first_direction ! An integer that indicates which direction is to be ! updated first in directionally split parts of the ! calculation. integer :: nkml, nkbl, verbosity, write_geom integer :: dynamics_stencil ! The computational stencil for the calculations ! in the dynamic core. - integer,allocatable :: pelist(:) ! list of pes for this instance of the ocean - integer :: mom_comm ! list of pes for this instance of the ocean - integer :: num_procs ! number of processors to pass to stochastic physics - integer :: iret ! return code from stochastic physics - integer :: me ! my pe - integer :: master ! root pe real :: conv2watt, conv2salt character(len=48) :: flux_units, S_flux_units @@ -2503,18 +2494,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call copy_dyngrid_to_MOM_grid(dG, G, US) call destroy_dyn_horgrid(dG) - do_epbl=.false. - do_sppt=.false. -#ifdef UFS - num_procs=num_PEs() - allocate(pelist(num_procs)) - call Get_PElist(pelist,commID = mom_comm) - me=PE_here() - master=root_PE() - - call init_stochastic_physics_ocn(CS%dt_therm,G%geoLonT,G%geoLatT,G%ied-G%isd+1,G%jed-G%jsd+1,nz,do_epbl,do_sppt,master,mom_comm,iret) -#endif - ! 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)) @@ -2969,9 +2948,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ! call fix_restart_scaling(GV) ! call fix_restart_unit_scaling(US) - CS%diabatic_CSp%do_epbl=do_epbl - CS%diabatic_CSp%do_sppt=do_sppt - call callTree_leave("initialize_MOM()") call cpu_clock_end(id_clock_init) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index a9bf6c3dcf..4020721075 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -276,6 +276,13 @@ module MOM_variables !> Container for information about the summed layer transports !! and how they will vary as the barotropic velocity is changed. +type, public :: stochastic_pattern + logical :: do_sppt = .false. + logical :: pert_epbl = .false. + real, allocatable :: sppt_wts(:,:) !< Random pattern for ocean SPPT + real, allocatable :: t_rp1(:,:) !< Random pattern for K.E. generation + real, allocatable :: t_rp2(:,:) !< Random pattern for K.E. dissipation +end type stochastic_pattern type, public :: BT_cont_type real, allocatable :: FA_u_EE(:,:) !< The effective open face area for zonal barotropic transport !! drawing from locations far to the east [H L ~> m2 or kg m-1]. diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 01ae5811af..df4dd4c453 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -65,13 +65,10 @@ module MOM_diabatic_driver use MOM_tracer_diabatic, only : tracer_vertdiff, tracer_vertdiff_Eulerian use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs, vertvisc_type, accel_diag_ptrs -use MOM_variables, only : cont_diag_ptrs, MOM_thermovar_chksum, p3d +use MOM_variables, only : cont_diag_ptrs, MOM_thermovar_chksum, p3d, stochastic_pattern use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units use MOM_wave_speed, only : wave_speeds use MOM_wave_interface, only : wave_parameters_CS -#ifdef UFS -use stochastic_physics, only : run_stochastic_physics_ocn -#endif implicit none ; private @@ -218,8 +215,6 @@ module MOM_diabatic_driver logical :: diabatic_diff_tendency_diag = .false. !< If true calculate diffusive tendency diagnostics logical :: boundary_forcing_tendency_diag = .false. !< If true calculate frazil diagnostics logical :: frazil_tendency_diag = .false. !< If true calculate frazil tendency diagnostics - logical,public :: do_epbl = .false. !< If true pertrub u_start in ePBL calculation - logical,public :: do_sppt = .false. !< If true perturb all physics tendenceies in MOM_diabatic_driver real, allocatable, dimension(:,:,:) :: frazil_heat_diag !< diagnose 3d heat tendency from frazil real, allocatable, dimension(:,:,:) :: frazil_temp_diag !< diagnose 3d temp tendency from frazil @@ -270,21 +265,21 @@ 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, 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 +subroutine diabatic(u, v, h, tv, Hml, fluxes, stochastics, 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 !! unused have NULL ptrs real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m] type(forcing), intent(inout) :: fluxes !< points to forcing fields !! unused fields have NULL ptrs - type(vertvisc_type), intent(inout) :: visc !< Structure with vertical viscosities, - !! BBL properties and related fields - type(accel_diag_ptrs), intent(inout) :: ADp !< Points to accelerations in momentum + type(stochastic_pattern), intent(in) :: stochastics !< random patterns + type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, BBL properies, and + type(accel_diag_ptrs), intent(inout) :: ADp !< related points to accelerations in momentum !! equations, to enable the later derived !! diagnostics, like energy budgets type(cont_diag_ptrs), intent(inout) :: CDp !< points to terms in continuity equations @@ -305,36 +300,24 @@ 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, dimension(SZI_(G),SZJ_(G)) :: sppt_wts - real, dimension(SZI_(G),SZJ_(G),2) :: t_rp -#ifdef UFS real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_in real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: t_in !< thickness [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: s_in !< thickness [H ~> m or kg m-2] real :: t_tend,s_tend,h_tend ! holder for tendencey needed for SPPT real :: t_pert,s_pert,h_pert ! holder for tendencey needed for SPPT -#endif if (G%ke == 1) return -#ifdef UFS ! save copy of the date for SPPT - if (CS%do_sppt) then - h_in=h - t_in=tv%T - s_in=tv%S - endif - call run_stochastic_physics_ocn(t_rp,sppt_wts) - if (CS%id_t_rp1 > 0) then - call post_data(CS%id_t_rp1, t_rp(:,:,1), CS%diag) - endif - if (CS%id_t_rp2 > 0) then - call post_data(CS%id_t_rp2, t_rp(:,:,2), CS%diag) - endif - if (CS%id_sppt_wts > 0) then - call post_data(CS%id_sppt_wts, sppt_wts, CS%diag) + if (stochastics%do_sppt) then + 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, stochastics%sppt_wts, CS%diag) + endif endif -#endif is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -420,10 +403,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & endif ! end CS%use_int_tides if (CS%useALEalgorithm .and. CS%use_legacy_diabatic) then - call diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, t_rp, visc, ADp, CDp, dt, Time_end, & + call diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, stochastics, visc, ADp, CDp, dt, Time_end, & G, GV, US, CS, Waves) elseif (CS%useALEalgorithm) then - call diabatic_ALE(u, v, h, tv, Hml, fluxes, t_rp, visc, ADp, CDp, dt, Time_end, & + call diabatic_ALE(u, v, h, tv, Hml, fluxes, stochastics, visc, ADp, CDp, dt, Time_end, & G, GV, US, CS, Waves) else call layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & @@ -490,14 +473,13 @@ 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) -#ifdef UFS - if (CS%do_sppt) then + if (stochastics%do_sppt) then do k=1,nz do j=js,je do i=is,ie - h_tend = (h(i,j,k)-h_in(i,j,k))*sppt_wts(i,j) - t_tend = (tv%T(i,j,k)-t_in(i,j,k))*sppt_wts(i,j) - s_tend = (tv%S(i,j,k)-s_in(i,j,k))*sppt_wts(i,j) + h_tend = (h(i,j,k)-h_in(i,j,k))*stochastics%sppt_wts(i,j) + t_tend = (tv%T(i,j,k)-t_in(i,j,k))*stochastics%sppt_wts(i,j) + s_tend = (tv%S(i,j,k)-s_in(i,j,k))*stochastics%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) @@ -514,7 +496,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo enddo endif -#endif end subroutine diabatic @@ -523,7 +504,7 @@ 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, t_rp, visc, ADp, CDp, dt, Time_end, & +subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, stochastics, 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 @@ -533,10 +514,10 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, t_rp, visc, ADp, CDp, d 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 !! unused have NULL ptrs - real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m] - type(forcing), intent(inout) :: fluxes !< points to forcing fields + real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m] + type(forcing), intent(inout) :: fluxes !< points to forcing fields + type(stochastic_pattern), intent(in) :: stochastics !< points to forcing fields !! unused fields have NULL ptrs - real, dimension(SZI_(G),SZJ_(G),2), intent(in) :: t_rp !< random pattern type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, BBL properies, and type(accel_diag_ptrs), intent(inout) :: ADp !< related points to accelerations in momentum !! equations, to enable the later derived @@ -846,7 +827,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, t_rp, visc, ADp, CDp, d endif 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, t_rp, CS%do_epbl, dt, Kd_ePBL, G, GV, US, & + call energetic_PBL(h, u_h, v_h, tv, fluxes, stochastics, dt, Kd_ePBL, G, GV, US, & CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) if (associated(Hml)) then @@ -1109,7 +1090,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, t_rp, visc, ADp, CDp, dt, Time_end, & +subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, stochastics, 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 @@ -1122,7 +1103,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, t_rp, visc, ADp, CDp, dt, Time real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m] type(forcing), intent(inout) :: fluxes !< points to forcing fields !! unused fields have NULL ptrs - real, dimension(SZI_(G),SZJ_(G),2), intent(in) :: t_rp !< random pattern + type(stochastic_pattern), intent(in) :: stochastics !< random patterns type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, BBL properies, and type(accel_diag_ptrs), intent(inout) :: ADp !< related points to accelerations in momentum !! equations, to enable the later derived @@ -1383,7 +1364,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, t_rp, visc, ADp, CDp, dt, Time endif 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, t_rp, CS%do_epbl, dt, Kd_ePBL, G, GV, US, & + call energetic_PBL(h, u_h, v_h, tv, fluxes, stochastics, dt, Kd_ePBL, G, GV, US, & CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) if (associated(Hml)) then @@ -3151,12 +3132,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di 'Zonal Acceleration from Diapycnal Mixing', 'm s-2', conversion=US%L_T2_to_m_s2) CS%id_dvdt_dia = register_diag_field('ocean_model', 'dvdt_dia', diag%axesCvL, Time, & 'Meridional Acceleration from Diapycnal Mixing', 'm s-2', conversion=US%L_T2_to_m_s2) - CS%id_t_rp1 = register_diag_field('ocean_model', 'random_pattern1', diag%axesT1, Time, & - 'random pattern1 for stochastics', 'None') - CS%id_t_rp2 = register_diag_field('ocean_model', 'random_pattern2', diag%axesT1, Time, & - 'random pattern2 for stochastics', 'None') CS%id_sppt_wts = register_diag_field('ocean_model', 'sppt_pattern', diag%axesT1, Time, & - 'random pattern for sppt', 'None') + 'random pattern for sppt', 'None') if (CS%use_int_tides) then CS%id_cg1 = register_diag_field('ocean_model', 'cn1', diag%axesT1, & diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 1cf089d5ed..a0b9ee0b51 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -14,7 +14,7 @@ module MOM_energetic_PBL use MOM_grid, only : ocean_grid_type use MOM_string_functions, only : uppercase use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : thermo_var_ptrs +use MOM_variables, only : thermo_var_ptrs, stochastic_pattern use MOM_verticalGrid, only : verticalGrid_type use MOM_wave_interface, only: wave_parameters_CS, Get_Langmuir_Number use MOM_stochastics, only : stochastic_CS @@ -165,7 +165,6 @@ module MOM_energetic_PBL !! potential energy change code. Otherwise, it uses a newer version !! that can work with successive increments to the diffusivity in !! upward or downward passes. - logical :: do_epbl type(diag_ctrl), pointer :: diag=>NULL() !< A structure that is used to regulate the !! timing of diagnostic output. @@ -198,6 +197,7 @@ 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_t_rp1=-1,id_t_rp2=-1 !>@} end type energetic_PBL_CS @@ -245,7 +245,7 @@ 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, t_rp, stoch_epbl, dt, Kd_int, G, GV, US, CS, & +subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, stochastics, dt, Kd_int, G, GV, US, 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. @@ -276,15 +276,14 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, stoch_epbl, dt, Kd_ type(forcing), intent(inout) :: fluxes !< A structure containing pointers to any !! possible forcing fields. Unused fields have !! NULL ptrs. + type(stochastic_pattern), intent(in) :: stochastics !< A structure containing array to any unsued fields + !! are not allocated real, intent(in) :: dt !< Time increment [T ~> s]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & intent(out) :: Kd_int !< The diagnosed diffusivities at interfaces !! [Z2 s-1 ~> m2 s-1]. type(energetic_PBL_CS), pointer :: CS !< The control structure returned by a previous !! call to mixedlayer_init. - real, dimension(SZI_(G),SZJ_(G),2), & - intent(in) :: t_rp !< random pattern to perturb wind - logical, intent(in) :: stoch_epbl !< flag to pertrub production and dissipation of TKE real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: buoy_flux !< The surface buoyancy flux [Z2 T-3 ~> m2 s-3]. type(wave_parameters_CS), pointer :: Waves !< Waves control structure for Langmuir turbulence @@ -429,9 +428,9 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, stoch_epbl, dt, Kd_ 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, t_rp(i,j,1),t_rp(i,j,2), stoch_epbl, dt_diag=dt_diag, Waves=Waves, G=G, i=i, j=j) + call ePBL_column(h, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, TKE_forcing, stochastics, & + 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) ! Copy the diffusivities to a 2-d array. do K=1,nz+1 @@ -470,26 +469,30 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, stoch_epbl, dt, Kd_ enddo ! j-loop - if (CS%id_ML_depth > 0) call post_data(CS%id_ML_depth, CS%ML_depth, CS%diag) - if (CS%id_hML_depth > 0) call post_data(CS%id_hML_depth, CS%ML_depth, CS%diag) - if (CS%id_TKE_wind > 0) call post_data(CS%id_TKE_wind, CS%diag_TKE_wind, CS%diag) - if (CS%id_TKE_MKE > 0) call post_data(CS%id_TKE_MKE, CS%diag_TKE_MKE, CS%diag) - if (CS%id_TKE_conv > 0) call post_data(CS%id_TKE_conv, CS%diag_TKE_conv, CS%diag) - if (CS%id_TKE_forcing > 0) call post_data(CS%id_TKE_forcing, CS%diag_TKE_forcing, CS%diag) - if (CS%id_TKE_mixing > 0) call post_data(CS%id_TKE_mixing, CS%diag_TKE_mixing, CS%diag) - if (CS%id_TKE_mech_decay > 0) & - call post_data(CS%id_TKE_mech_decay, CS%diag_TKE_mech_decay, CS%diag) - if (CS%id_TKE_conv_decay > 0) & - call post_data(CS%id_TKE_conv_decay, CS%diag_TKE_conv_decay, CS%diag) - if (CS%id_Mixing_Length > 0) call post_data(CS%id_Mixing_Length, CS%Mixing_Length, CS%diag) - if (CS%id_Velocity_Scale >0) call post_data(CS%id_Velocity_Scale, CS%Velocity_Scale, CS%diag) - if (CS%id_MSTAR_MIX > 0) call post_data(CS%id_MSTAR_MIX, CS%MSTAR_MIX, CS%diag) - 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) - 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) + if (write_diags) then + if (CS%id_ML_depth > 0) call post_data(CS%id_ML_depth, CS%ML_depth, CS%diag) + if (CS%id_hML_depth > 0) call post_data(CS%id_hML_depth, CS%ML_depth, CS%diag) + if (CS%id_TKE_wind > 0) call post_data(CS%id_TKE_wind, CS%diag_TKE_wind, CS%diag) + if (CS%id_TKE_MKE > 0) call post_data(CS%id_TKE_MKE, CS%diag_TKE_MKE, CS%diag) + if (CS%id_TKE_conv > 0) call post_data(CS%id_TKE_conv, CS%diag_TKE_conv, CS%diag) + if (CS%id_TKE_forcing > 0) call post_data(CS%id_TKE_forcing, CS%diag_TKE_forcing, CS%diag) + if (CS%id_TKE_mixing > 0) call post_data(CS%id_TKE_mixing, CS%diag_TKE_mixing, CS%diag) + if (CS%id_TKE_mech_decay > 0) & + call post_data(CS%id_TKE_mech_decay, CS%diag_TKE_mech_decay, CS%diag) + if (CS%id_TKE_conv_decay > 0) & + call post_data(CS%id_TKE_conv_decay, CS%diag_TKE_conv_decay, CS%diag) + if (CS%id_Mixing_Length > 0) call post_data(CS%id_Mixing_Length, CS%Mixing_Length, CS%diag) + if (CS%id_Velocity_Scale >0) call post_data(CS%id_Velocity_Scale, CS%Velocity_Scale, CS%diag) + if (CS%id_MSTAR_MIX > 0) call post_data(CS%id_MSTAR_MIX, CS%MSTAR_MIX, CS%diag) + 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 (stochastics%pert_epbl) then + if (CS%id_t_rp1 > 0) call post_data(CS%id_t_rp1, stochastics%t_rp1, CS%diag) + if (CS%id_t_rp2 > 0) call post_data(CS%id_t_rp2, stochastics%t_rp2, CS%diag) + endif endif end subroutine energetic_PBL @@ -497,9 +500,9 @@ end subroutine energetic_PBL !> This subroutine determines the diffusivities from the integrated energetics !! 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, & +subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, stochastics, B_flux, absf, & u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, US, CS, eCD, & - t_rp1,t_rp2, stoch_epbl, dt_diag, Waves, G, i, j) + dt_diag, Waves, G, 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]. @@ -518,6 +521,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs real, dimension(SZK_(GV)), intent(in) :: TKE_forcing !< The forcing requirements to homogenize the !! forcing that has been applied to each layer !! [R Z3 T-2 ~> J m-2]. + type(stochastic_pattern), intent(in) :: stochastics !< stochastic patterns and logic controls real, intent(in) :: B_flux !< The surface buoyancy flux [Z2 T-3 ~> m2 s-3] real, intent(in) :: absf !< The absolute value of the Coriolis parameter [T-1 ~> s-1]. real, intent(in) :: u_star !< The surface friction velocity [Z T-1 ~> m s-1]. @@ -536,9 +540,6 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs intent(out) :: mixlen !< The mixing length scale used in Kd [Z ~> m]. type(energetic_PBL_CS), intent(inout) :: CS !< Energetic PBL control struct type(ePBL_column_diags), intent(inout) :: eCD !< A container for passing around diagnostics. - real, intent(in) :: t_rp1 !< random value to perturb TKE production - real, intent(in) :: t_rp2 !< random value to perturb TKE dissipation - logical, intent(in) :: stoch_epbl !< flag to pertrub production and dissipation of TKE real, optional, intent(in) :: dt_diag !< The diagnostic time step, which may be less !! than dt if there are two calls to mixedlayer [T ~> s]. type(wave_parameters_CS), & @@ -836,7 +837,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 (stoch_epbl) mech_TKE=mech_TKE*t_rp1 + if (stochastics%pert_epbl) mech_TKE=mech_TKE*stochastics%t_rp1(i,j) if (CS%TKE_diagnostics) then eCD%dTKE_conv = 0.0 ; eCD%dTKE_mixing = 0.0 @@ -919,8 +920,8 @@ 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 (stoch_epbl) then ! perturb the TKE destruction - mech_TKE = mech_TKE * (1+(exp_kh-1) * t_rp2) + if (stochastics%pert_epbl) then ! perturb the TKE destruction + mech_TKE = mech_TKE * (1+(exp_kh-1) * stochastics%t_rp2(i,j)) else mech_TKE = mech_TKE * exp_kh endif @@ -2326,6 +2327,10 @@ 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_t_rp1 = register_diag_field('ocean_model', 'random_pattern1', diag%axesT1, Time, & + 'random pattern1 for stochastics', 'None') + CS%id_t_rp2 = register_diag_field('ocean_model', 'random_pattern2', diag%axesT1, Time, & + 'random pattern2 for stochastics', 'None') if (CS%use_LT) then CS%id_LA = register_diag_field('ocean_model', 'LA', diag%axesT1, & Time, 'Langmuir number.', 'nondim')