From 182ef34031b24253d905bba72ba1feea0687fd6d Mon Sep 17 00:00:00 2001 From: Philip Pegion Date: Tue, 5 May 2020 09:41:38 -0500 Subject: [PATCH 01/21] additions for stochastic physics and ePBL perts --- .../mom_surface_forcing_nuopc.F90 | 2 ++ src/core/MOM.F90 | 25 ++++++++++++++++++- src/core/MOM_forcing_type.F90 | 13 ++++++++++ src/diagnostics/MOM_diagnostics.F90 | 13 +++++++++- src/framework/MOM_domains.F90 | 4 +-- .../vertical/MOM_energetic_PBL.F90 | 2 +- 6 files changed, 54 insertions(+), 5 deletions(-) diff --git a/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 b/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 index af59d7d6ea..868a281c62 100644 --- a/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 +++ b/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 @@ -285,6 +285,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, & call safe_alloc_ptr(fluxes%p_surf ,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%p_surf_full,isd,ied,jsd,jed) + print*,'allocate fluxes%t_rp' + call safe_alloc_ptr(fluxes%t_rp,isd,ied,jsd,jed) if (CS%use_limited_P_SSH) then fluxes%p_surf_SSH => fluxes%p_surf else diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 23c11cc05b..1709913b8a 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -27,6 +27,7 @@ module MOM use MOM_domains, only : To_All, Omit_corners, CGRID_NE, SCALAR_PAIR use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_domains, only : start_group_pass, complete_group_pass, Omit_Corners +use MOM_domains, only : root_PE,PE_here,Get_PElist,num_PEs use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe use MOM_error_handler, only : MOM_set_verbosity, callTree_showQuery use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint @@ -130,6 +131,7 @@ module MOM use MOM_offline_main, only : offline_fw_fluxes_into_ocean, offline_fw_fluxes_out_ocean use MOM_offline_main, only : offline_advection_layer, offline_transport_end use MOM_ALE, only : ale_offline_tracer_final, ALE_main_offline +use stochastic_physics, only : init_stochastic_physics_ocn,run_stochastic_physics_ocn implicit none ; private @@ -212,6 +214,8 @@ module MOM logical :: offline_tracer_mode = .false. !< If true, step_offline() is called instead of step_MOM(). !! This is intended for running MOM6 in offline tracer mode + logical :: do_stochy = .false. + !< If true, call stochastic physics pattern generator type(time_type), pointer :: Time !< pointer to the ocean clock real :: dt !< (baroclinic) dynamics time step [s] @@ -703,6 +707,9 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & enddo ; enddo endif + print*,'calling run_stochastic_physics_ocn',CS%do_stochy + if (CS%do_stochy) call run_stochastic_physics_ocn(forces%t_rp) + call step_MOM_dynamics(forces, CS%p_surf_begin, CS%p_surf_end, dt, & dt_therm_here, bbl_time_int, CS, & Time_local, Waves=Waves) @@ -843,7 +850,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & if (CS%time_in_thermo_cycle > 0.0) then call enable_averaging(CS%time_in_thermo_cycle, Time_local, CS%diag) call post_surface_thermo_diags(CS%sfc_IDs, G, GV, US, CS%diag, CS%time_in_thermo_cycle, & - sfc_state, CS%tv, ssh, CS%ave_ssh_ibc) + sfc_state, CS%tv, ssh, fluxes%t_rp, CS%ave_ssh_ibc) endif call disable_averaging(CS%diag) call cpu_clock_end(id_clock_diagnostics) @@ -1580,6 +1587,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & integer :: nkml, nkbl, verbosity, write_geom integer :: dynamics_stencil ! The computational stencil for the calculations ! in the dynamic core. + integer,allocatable :: pelist(:) ! list of pes for this instance of the ocean + integer :: num_procs +! model + integer :: me ! my pe + integer :: master ! root pe real :: conv2watt, conv2salt character(len=48) :: flux_units, S_flux_units @@ -2146,6 +2158,17 @@ 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) + num_procs=num_PEs() + allocate(pelist(num_procs)) + call Get_PElist(pelist) + me=PE_here() + master=root_PE() + + !call init_stochastic_physics_ocn(CS%dt_therm,G,me,master,pelist,CS%do_stochy) + print*,'callling init_stochastic_physics_ocn',maxval(G%geoLatT) + call init_stochastic_physics_ocn(CS%dt_therm,G%geoLonT,G%geoLatT,G%ied-G%isd+1,G%jed-G%jsd+1,nz,CS%do_stochy) + print*,'back from init_stochastic_physics_ocn',CS%do_stochy + ! Set a few remaining fields that are specific to the ocean grid type. call set_first_direction(G, first_direction) ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes. diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 7df4213a2f..07d567be43 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -123,6 +123,8 @@ module MOM_forcing_type !< Pressure at the top ocean interface [Pa] that is used in corrections to the sea surface !! height field that is passed back to the calling routines. !! p_surf_SSH may point to p_surf or to p_surf_full. + real, pointer, dimension(:,:) :: t_rp => NULL() + !< random pattern at t-points logical :: accumulate_p_surf = .false. !< If true, the surface pressure due to the atmosphere !! and various types of ice needs to be accumulated, and the !! surface pressure explicitly reset to zero at the driver level @@ -217,6 +219,8 @@ module MOM_forcing_type real, pointer, dimension(:,:) :: & rigidity_ice_u => NULL(), & !< Depth-integrated lateral viscosity of ice shelves or sea ice at u-points [m3 s-1] rigidity_ice_v => NULL() !< Depth-integrated lateral viscosity of ice shelves or sea ice at v-points [m3 s-1] + real, pointer, dimension(:,:) :: t_rp => NULL() + !< random pattern at t-points real :: dt_force_accum = -1.0 !< The amount of time over which the mechanical forcing fluxes !! have been averaged [s]. logical :: net_mass_src_set = .false. !< If true, an estimate of net_mass_src has been provided. @@ -2020,6 +2024,12 @@ subroutine copy_common_forcing_fields(forces, fluxes, G, skip_pres) do_pres = .true. ; if (present(skip_pres)) do_pres = .not.skip_pres + if (associated(forces%t_rp) .and. associated(fluxes%t_rp)) then + do j=js,je ; do i=is,ie + fluxes%t_rp(i,j) = forces%t_rp(i,j) + enddo ; enddo + endif + if (associated(forces%ustar) .and. associated(fluxes%ustar)) then do j=js,je ; do i=is,ie fluxes%ustar(i,j) = forces%ustar(i,j) @@ -2845,6 +2855,7 @@ subroutine allocate_mech_forcing(G, forces, stress, ustar, shelf, press, iceberg call myAlloc(forces%p_surf,isd,ied,jsd,jed, press) call myAlloc(forces%p_surf_full,isd,ied,jsd,jed, press) call myAlloc(forces%net_mass_src,isd,ied,jsd,jed, press) + call myAlloc(forces%t_rp,isd,ied,jsd,jed, press) call myAlloc(forces%rigidity_ice_u,IsdB,IedB,jsd,jed, shelf) call myAlloc(forces%rigidity_ice_v,isd,ied,JsdB,JedB, shelf) @@ -2908,6 +2919,7 @@ subroutine deallocate_forcing_type(fluxes) if (associated(fluxes%seaice_melt)) deallocate(fluxes%seaice_melt) if (associated(fluxes%salt_flux)) deallocate(fluxes%salt_flux) if (associated(fluxes%p_surf_full)) deallocate(fluxes%p_surf_full) + if (associated(fluxes%t_rp)) deallocate(fluxes%t_rp) if (associated(fluxes%p_surf)) deallocate(fluxes%p_surf) if (associated(fluxes%TKE_tidal)) deallocate(fluxes%TKE_tidal) if (associated(fluxes%ustar_tidal)) deallocate(fluxes%ustar_tidal) @@ -2932,6 +2944,7 @@ subroutine deallocate_mech_forcing(forces) if (associated(forces%ustar)) deallocate(forces%ustar) if (associated(forces%p_surf)) deallocate(forces%p_surf) if (associated(forces%p_surf_full)) deallocate(forces%p_surf_full) + if (associated(forces%t_rp)) deallocate(forces%t_rp) if (associated(forces%net_mass_src)) deallocate(forces%net_mass_src) if (associated(forces%rigidity_ice_u)) deallocate(forces%rigidity_ice_u) if (associated(forces%rigidity_ice_v)) deallocate(forces%rigidity_ice_v) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 54025a0ac0..4f951863f0 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -167,6 +167,8 @@ module MOM_diagnostics integer :: id_salt_deficit = -1 integer :: id_Heat_PmE = -1 integer :: id_intern_heat = -1 +! stochastic pattern + integer :: id_t_rp = -1 !!@} end type surface_diag_IDs @@ -1193,7 +1195,7 @@ end subroutine post_surface_dyn_diags !> This routine posts diagnostics of various ocean surface and integrated !! quantities at the time the ocean state is reported back to the caller subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv, & - ssh, ssh_ibc) + ssh, t_rp, ssh_ibc) type(surface_diag_IDs), intent(in) :: IDs !< A structure with the diagnostic IDs. type(ocean_grid_type), intent(in) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure @@ -1204,6 +1206,8 @@ subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: ssh !< Time mean surface height without corrections for ice displacement [m] + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: t_rp!< random pattern for stochastic proceeses real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh_ibc !< Time mean surface height with corrections !! for ice displacement and the inverse barometer [m] @@ -1328,6 +1332,11 @@ subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv call post_data(IDs%id_sss_sq, work_2d, diag, mask=G%mask2dT) endif + if (IDs%id_t_rp > 0) then + !call post_data(IDs%id_t_rp, t_rp, diag, mask=G%mask2dT) + call post_data(IDs%id_t_rp, t_rp, diag) + endif + call coupler_type_send_data(sfc_state%tr_fields, get_diag_time_end(diag)) end subroutine post_surface_thermo_diags @@ -1789,6 +1798,8 @@ subroutine register_surface_diags(Time, G, IDs, diag, tv) 'Heat flux into ocean from mass flux into ocean', 'W m-2') IDs%id_intern_heat = register_diag_field('ocean_model', 'internal_heat', diag%axesT1, Time,& 'Heat flux into ocean from geothermal or other internal sources', 'W m-2') + IDs%id_t_rp = register_diag_field('ocean_model', 'random_pattern', diag%axesT1, Time, & + 'random pattern for stochastics', 'None') end subroutine register_surface_diags diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index 64fddfe7fc..c28530fe65 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -3,7 +3,7 @@ module MOM_domains ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_coms, only : PE_here, root_PE, num_PEs, MOM_infra_init, MOM_infra_end +use MOM_coms, only : PE_here, root_PE, num_PEs, MOM_infra_init, MOM_infra_end, Get_PElist use MOM_coms, only : broadcast, sum_across_PEs, min_across_PEs, max_across_PEs use MOM_cpu_clock, only : cpu_clock_begin, cpu_clock_end use MOM_error_handler, only : MOM_error, MOM_mesg, NOTE, WARNING, FATAL, is_root_pe @@ -34,7 +34,7 @@ module MOM_domains public :: MOM_domains_init, MOM_infra_init, MOM_infra_end, get_domain_extent, get_domain_extent_dsamp2 public :: MOM_define_domain, MOM_define_io_domain, clone_MOM_domain -public :: pass_var, pass_vector, PE_here, root_PE, num_PEs +public :: pass_var, pass_vector, PE_here, root_PE, num_PEs, Get_PElist public :: pass_var_start, pass_var_complete, fill_symmetric_edges, broadcast public :: pass_vector_start, pass_vector_complete public :: global_field_sum, sum_across_PEs, min_across_PEs, max_across_PEs diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index b486e1e2ca..1f91d9549e 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -427,7 +427,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS do K=1,nz+1 ; Kd(K) = 0.0 ; enddo ! Make local copies of surface forcing and process them. - u_star = fluxes%ustar(i,j) + u_star = fluxes%ustar(i,j)*(fluxes%t_rp(i,j)) u_star_Mean = fluxes%ustar_gustless(i,j) B_flux = buoy_flux(i,j) if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then From 7de295c86cf9c2c88975a361418b945066bd11d0 Mon Sep 17 00:00:00 2001 From: pjpegion Date: Wed, 2 Dec 2020 15:35:05 +0000 Subject: [PATCH 02/21] cleanup of code and enhancement of ePBL perts --- .../mom_surface_forcing_nuopc.F90 | 2 - src/core/MOM.F90 | 28 ++++--- src/core/MOM_forcing_type.F90 | 37 ++++++--- src/diagnostics/MOM_diagnostics.F90 | 13 +--- .../vertical/MOM_diabatic_driver.F90 | 76 +++++++++++++++++-- .../vertical/MOM_energetic_PBL.F90 | 24 ++++-- 6 files changed, 129 insertions(+), 51 deletions(-) diff --git a/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 b/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 index 8279619fd1..3516ad3803 100644 --- a/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 +++ b/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 @@ -298,8 +298,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, call safe_alloc_ptr(fluxes%p_surf ,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%p_surf_full,isd,ied,jsd,jed) - print*,'allocate fluxes%t_rp' - call safe_alloc_ptr(fluxes%t_rp,isd,ied,jsd,jed) if (CS%use_limited_P_SSH) then fluxes%p_surf_SSH => fluxes%p_surf else diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index dd422ab445..6ed974708e 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -142,7 +142,7 @@ module MOM use MOM_offline_main, only : offline_fw_fluxes_into_ocean, offline_fw_fluxes_out_ocean use MOM_offline_main, only : offline_advection_layer, offline_transport_end use MOM_ALE, only : ale_offline_tracer_final, ALE_main_offline -use stochastic_physics, only : init_stochastic_physics_ocn,run_stochastic_physics_ocn +use stochastic_physics, only : init_stochastic_physics_ocn implicit none ; private @@ -761,9 +761,6 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS enddo ; enddo endif - print*,'calling run_stochastic_physics_ocn',CS%do_stochy - if (CS%do_stochy) call run_stochastic_physics_ocn(forces%t_rp) - call step_MOM_dynamics(forces, CS%p_surf_begin, CS%p_surf_end, dt, & dt_therm_here, bbl_time_int, CS, & Time_local, Waves=Waves) @@ -911,7 +908,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS if (CS%time_in_thermo_cycle > 0.0) then call enable_averages(CS%time_in_thermo_cycle, Time_local, CS%diag) call post_surface_thermo_diags(CS%sfc_IDs, G, GV, US, CS%diag, CS%time_in_thermo_cycle, & - sfc_state_diag, CS%tv, ssh,fluxes%t_rp, CS%ave_ssh_ibc) + sfc_state_diag, CS%tv, ssh, CS%ave_ssh_ibc) + !sfc_state_diag, CS%tv, ssh,fluxes%t_rp,fluxes%sppt_wts, CS%ave_ssh_ibc) endif call disable_averaging(CS%diag) call cpu_clock_end(id_clock_diagnostics) @@ -1672,6 +1670,7 @@ 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. This can be altered during the course @@ -1680,7 +1679,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & 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 :: num_procs + integer :: mom_comm ! list of pes for this instance of the ocean + integer :: num_procs,iret ! model integer :: me ! my pe integer :: master ! root pe @@ -2342,14 +2342,17 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & num_procs=num_PEs() allocate(pelist(num_procs)) - call Get_PElist(pelist) + call Get_PElist(pelist,commID = mom_comm) me=PE_here() master=root_PE() - !call init_stochastic_physics_ocn(CS%dt_therm,G,me,master,pelist,CS%do_stochy) - print*,'callling init_stochastic_physics_ocn',maxval(G%geoLatT) - call init_stochastic_physics_ocn(CS%dt_therm,G%geoLonT,G%geoLatT,G%ied-G%isd+1,G%jed-G%jsd+1,nz,CS%do_stochy) - print*,'back from init_stochastic_physics_ocn',CS%do_stochy + !print*,'callling init_stochastic_physics_ocn',maxval(G%geoLatT) + do_epbl=.false. + do_sppt=.false. + 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) + if (do_sppt .eq. .true.) CS%do_stochy=.true. + if (do_epbl .eq. .true.) CS%do_stochy=.true. + !print*,'back from init_stochastic_physics_ocn',CS%do_stochy ! Set a few remaining fields that are specific to the ocean grid type. call set_first_direction(G, first_direction) @@ -2763,6 +2766,9 @@ 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_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 5f48f4d7d6..930f54566d 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -141,8 +141,10 @@ module MOM_forcing_type !< Pressure at the top ocean interface [R L2 T-2 ~> Pa] that is used in corrections to the sea surface !! height field that is passed back to the calling routines. !! p_surf_SSH may point to p_surf or to p_surf_full. - real, pointer, dimension(:,:) :: t_rp => NULL() - !< random pattern at t-points +! real, pointer, dimension(:,:) :: t_rp => NULL() +! !< random pattern at t-points +! real, pointer, dimension(:,:) :: sppt_wts => NULL() +! !< random pattern at t-points logical :: accumulate_p_surf = .false. !< If true, the surface pressure due to the atmosphere !! and various types of ice needs to be accumulated, and the !! surface pressure explicitly reset to zero at the driver level @@ -240,8 +242,10 @@ module MOM_forcing_type !! u-points [L4 Z-1 T-1 ~> m3 s-1] rigidity_ice_v => NULL() !< Depth-integrated lateral viscosity of ice shelves or sea ice at !! v-points [L4 Z-1 T-1 ~> m3 s-1] - real, pointer, dimension(:,:) :: t_rp => NULL() - !< random pattern at t-points +! real, pointer, dimension(:,:) :: t_rp => NULL() +! !< random pattern at t-points +! real, pointer, dimension(:,:) :: sppt_wts => NULL() +! !< random pattern at t-points real :: dt_force_accum = -1.0 !< The amount of time over which the mechanical forcing fluxes !! have been averaged [s]. logical :: net_mass_src_set = .false. !< If true, an estimate of net_mass_src has been provided. @@ -2082,11 +2086,17 @@ subroutine copy_common_forcing_fields(forces, fluxes, G, skip_pres) do_pres = .true. ; if (present(skip_pres)) do_pres = .not.skip_pres - if (associated(forces%t_rp) .and. associated(fluxes%t_rp)) then - do j=js,je ; do i=is,ie - fluxes%t_rp(i,j) = forces%t_rp(i,j) - enddo ; enddo - endif +! if (associated(forces%t_rp) .and. associated(fluxes%t_rp)) then +! do j=js,je ; do i=is,ie +! fluxes%t_rp(i,j) = forces%t_rp(i,j) +! enddo ; enddo +! endif +! +! if (associated(forces%sppt_wts) .and. associated(fluxes%sppt_wts)) then +! do j=js,je ; do i=is,ie +! fluxes%sppt_wts(i,j) = forces%sppt_wts(i,j) +! enddo ; enddo +! endif if (associated(forces%ustar) .and. associated(fluxes%ustar)) then do j=js,je ; do i=is,ie @@ -3031,7 +3041,8 @@ subroutine allocate_mech_forcing_by_group(G, forces, stress, ustar, shelf, & call myAlloc(forces%p_surf,isd,ied,jsd,jed, press) call myAlloc(forces%p_surf_full,isd,ied,jsd,jed, press) call myAlloc(forces%net_mass_src,isd,ied,jsd,jed, press) - call myAlloc(forces%t_rp,isd,ied,jsd,jed, press) +! call myAlloc(forces%t_rp,isd,ied,jsd,jed, press) +! call myAlloc(forces%sppt_wts,isd,ied,jsd,jed, press) call myAlloc(forces%rigidity_ice_u,IsdB,IedB,jsd,jed, shelf) call myAlloc(forces%rigidity_ice_v,isd,ied,JsdB,JedB, shelf) @@ -3187,7 +3198,8 @@ subroutine deallocate_forcing_type(fluxes) if (associated(fluxes%seaice_melt)) deallocate(fluxes%seaice_melt) if (associated(fluxes%salt_flux)) deallocate(fluxes%salt_flux) if (associated(fluxes%p_surf_full)) deallocate(fluxes%p_surf_full) - if (associated(fluxes%t_rp)) deallocate(fluxes%t_rp) +! if (associated(fluxes%t_rp)) deallocate(fluxes%t_rp) +! if (associated(fluxes%sppt_wts)) deallocate(fluxes%sppt_wts) if (associated(fluxes%p_surf)) deallocate(fluxes%p_surf) if (associated(fluxes%TKE_tidal)) deallocate(fluxes%TKE_tidal) if (associated(fluxes%ustar_tidal)) deallocate(fluxes%ustar_tidal) @@ -3212,7 +3224,8 @@ subroutine deallocate_mech_forcing(forces) if (associated(forces%ustar)) deallocate(forces%ustar) if (associated(forces%p_surf)) deallocate(forces%p_surf) if (associated(forces%p_surf_full)) deallocate(forces%p_surf_full) - if (associated(forces%t_rp)) deallocate(forces%t_rp) +! if (associated(forces%t_rp)) deallocate(forces%t_rp) +! if (associated(forces%sppt_wts)) deallocate(forces%sppt_wts) if (associated(forces%net_mass_src)) deallocate(forces%net_mass_src) if (associated(forces%rigidity_ice_u)) deallocate(forces%rigidity_ice_u) if (associated(forces%rigidity_ice_v)) deallocate(forces%rigidity_ice_v) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index c8b5d6ac4e..50873863fd 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -169,8 +169,6 @@ module MOM_diagnostics integer :: id_salt_deficit = -1 integer :: id_Heat_PmE = -1 integer :: id_intern_heat = -1 -! stochastic pattern - integer :: id_t_rp = -1 !!@} end type surface_diag_IDs @@ -1199,7 +1197,7 @@ end subroutine post_surface_dyn_diags !> This routine posts diagnostics of various ocean surface and integrated !! quantities at the time the ocean state is reported back to the caller subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv, & - ssh, t_rp, ssh_ibc) + ssh, ssh_ibc) type(surface_diag_IDs), intent(in) :: IDs !< A structure with the diagnostic IDs. type(ocean_grid_type), intent(in) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure @@ -1210,8 +1208,6 @@ subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: ssh !< Time mean surface height without corrections for ice displacement [m] - real, dimension(SZI_(G),SZJ_(G)), & - intent(in) :: t_rp!< random pattern for stochastic proceeses real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh_ibc !< Time mean surface height with corrections !! for ice displacement and the inverse barometer [m] @@ -1336,11 +1332,6 @@ subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv call post_data(IDs%id_sss_sq, work_2d, diag, mask=G%mask2dT) endif - if (IDs%id_t_rp > 0) then - !call post_data(IDs%id_t_rp, t_rp, diag, mask=G%mask2dT) - call post_data(IDs%id_t_rp, t_rp, diag) - endif - call coupler_type_send_data(sfc_state%tr_fields, get_diag_time_end(diag)) end subroutine post_surface_thermo_diags @@ -1842,8 +1833,6 @@ subroutine register_surface_diags(Time, G, US, IDs, diag, tv) IDs%id_intern_heat = register_diag_field('ocean_model', 'internal_heat', diag%axesT1, Time,& 'Heat flux into ocean from geothermal or other internal sources', & 'W m-2', conversion=US%QRZ_T_to_W_m2) - IDs%id_t_rp = register_diag_field('ocean_model', 'random_pattern', diag%axesT1, Time, & - 'random pattern for stochastics', 'None') end subroutine register_surface_diags diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index d753afc97b..fc76e87a0c 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 stochastic_physics, only : run_stochastic_physics_ocn implicit none ; private @@ -175,7 +176,7 @@ module MOM_diabatic_driver 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_brine_lay = -1 + integer :: id_subMLN2 = -1, id_brine_lay = -1, id_sppt_wts = -1, id_t_rp1=-1,id_t_rp2=-1 ! diagnostic for fields prior to applying diapycnal physics integer :: id_u_predia = -1, id_v_predia = -1, id_h_predia = -1 @@ -207,6 +208,8 @@ 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 @@ -286,8 +289,35 @@ 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 + 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 + if (G%ke == 1) return + ! 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) + !print*,'in diabatic',CS%do_sppt,size(t_in,1),size(t_in,2),size(t_in,3),size(sppt_wts,1),size(sppt_wts,2) + !print*,'in diabatic',CS%do_sppt,minval(sppt_wts),maxval(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) + endif + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke if (.not. associated(CS)) call MOM_error(FATAL, "MOM_diabatic_driver: "// & @@ -369,10 +399,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, visc, ADp, CDp, dt, Time_end, & + call diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, t_rp, visc, ADp, CDp, dt, Time_end, & G, GV, US, CS, Waves) elseif (CS%useALEalgorithm) then - call diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & + call diabatic_ALE(u, v, h, tv, Hml, fluxes, t_rp, 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, & @@ -435,13 +465,37 @@ 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 + 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_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 + 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, & +subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, t_rp, 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 @@ -454,6 +508,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim 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(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 @@ -836,7 +891,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim 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, dt, Kd_ePBL, G, GV, US, & + call energetic_PBL(h, u_h, v_h, tv, fluxes, t_rp, dt, Kd_ePBL, G, GV, US, & CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) if (associated(Hml)) then @@ -1222,7 +1277,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, & +subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, t_rp, 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 @@ -1235,6 +1290,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, 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(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 @@ -1565,7 +1621,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, 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, dt, Kd_ePBL, G, GV, US, & + call energetic_PBL(h, u_h, v_h, tv, fluxes, t_rp, dt, Kd_ePBL, G, GV, US, & CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) if (associated(Hml)) then @@ -3385,6 +3441,12 @@ 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') 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 9fe95cac56..4f93c47e95 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -163,6 +163,7 @@ 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. @@ -245,7 +246,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, dt, Kd_int, G, GV, US, CS, & +subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, 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. @@ -282,6 +283,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS !! [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 real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: buoy_flux !< The surface buoyancy flux [Z2 T-3 ~> m2 s-3]. real, optional, intent(in) :: dt_diag !< The diagnostic time step, which may be less @@ -435,7 +438,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS do K=1,nz+1 ; Kd(K) = 0.0 ; enddo ! Make local copies of surface forcing and process them. - u_star = fluxes%ustar(i,j)*(fluxes%t_rp(i,j)) + !print*,'PJP EPBL',minval(t_rp),maxval(t_rp) + u_star = fluxes%ustar(i,j)!*t_rp(i,j) u_star_Mean = fluxes%ustar_gustless(i,j) B_flux = buoy_flux(i,j) if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then @@ -459,9 +463,10 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS 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) - + US, CS, eCD, t_rp(i,j,1),t_rp(i,j,2), dt_diag=dt_diag, Waves=Waves, G=G, i=i, j=j) + ! applly stochastic perturbation to TKE generation + ! Copy the diffusivities to a 2-d array. do K=1,nz+1 Kd_2d(i,K) = Kd(K) @@ -542,7 +547,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) + t_rp1,t_rp2, 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]. @@ -580,6 +585,8 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs type(energetic_PBL_CS), pointer :: CS !< The control structure returned by a previous !! call to mixedlayer_init. 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 production 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), & @@ -878,6 +885,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 + mech_TKE=mech_TKE*t_rp1 if (CS%TKE_diagnostics) then eCD%dTKE_conv = 0.0 ; eCD%dTKE_mixing = 0.0 @@ -959,8 +968,9 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs exp_kh = 1.0 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 + !eCD%dTKE_mech_decay = eCD%dTKE_mech_decay + (exp_kh-1.0) * mech_TKE * I_dtdiag + eCD%dTKE_mech_decay = exp_kh + mech_TKE = mech_TKE * (1+(exp_kh-1) * t_rp2) ! Accumulate any convectively released potential energy to contribute ! to wstar and to drive penetrating convection. From 7212400539c7e9678087c4087e8118192c883b5e Mon Sep 17 00:00:00 2001 From: Phil Pegion <38869668+pjpegion@users.noreply.github.com> Date: Wed, 2 Dec 2020 09:03:36 -0700 Subject: [PATCH 03/21] Update MOM_diabatic_driver.F90 remove conflict with dev/emc --- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index fc76e87a0c..fe88221f1a 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -176,7 +176,7 @@ module MOM_diabatic_driver 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_brine_lay = -1, id_sppt_wts = -1, id_t_rp1=-1,id_t_rp2=-1 + integer :: id_subMLN2 = -1, id_sppt_wts = -1, id_t_rp1=-1,id_t_rp2=-1 ! diagnostic for fields prior to applying diapycnal physics integer :: id_u_predia = -1, id_v_predia = -1, id_h_predia = -1 From bd477a951a2aacc9718eb167c8db7ab0d778d6fe Mon Sep 17 00:00:00 2001 From: Phil Pegion <38869668+pjpegion@users.noreply.github.com> Date: Wed, 2 Dec 2020 09:05:28 -0700 Subject: [PATCH 04/21] Update MOM_diabatic_driver.F90 further resolve conflict --- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index fe88221f1a..852cbe52fe 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -176,7 +176,7 @@ module MOM_diabatic_driver 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, id_t_rp1=-1,id_t_rp2=-1 + integer :: id_subMLN2 = -1 ! diagnostic for fields prior to applying diapycnal physics integer :: id_u_predia = -1, id_v_predia = -1, id_h_predia = -1 From 0c15f4c44c21b438da5afdacc867dc9815d76f6c Mon Sep 17 00:00:00 2001 From: Phil Pegion <38869668+pjpegion@users.noreply.github.com> Date: Wed, 2 Dec 2020 09:06:52 -0700 Subject: [PATCH 05/21] Update MOM_diabatic_driver.F90 put id_sppt_wts, etc back. --- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 9196b4b03a..15d41293c1 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -176,7 +176,7 @@ module MOM_diabatic_driver 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 + integer :: id_subMLN2 = -1, id_sppt_wts = -1, id_t_rp1=-1,id_t_rp2=-1 ! diagnostic for fields prior to applying diapycnal physics integer :: id_u_predia = -1, id_v_predia = -1, id_h_predia = -1 From a2a374bce160c71d63cf9e4e08147f988788ae89 Mon Sep 17 00:00:00 2001 From: pjpegion Date: Mon, 14 Dec 2020 15:47:16 +0000 Subject: [PATCH 06/21] add stochy_restart writing to mom_cap --- config_src/nuopc_driver/mom_cap.F90 | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/config_src/nuopc_driver/mom_cap.F90 b/config_src/nuopc_driver/mom_cap.F90 index 8d48607281..6880addfb2 100644 --- a/config_src/nuopc_driver/mom_cap.F90 +++ b/config_src/nuopc_driver/mom_cap.F90 @@ -91,6 +91,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 get_stochy_pattern_mod, only: write_stoch_restart_ocn implicit none; private @@ -1584,6 +1585,16 @@ subroutine ModelAdvance(gcomp, rc) ! write restart file(s) call ocean_model_restart(ocean_state, restartname=restartname) + + ! write stochastic physics restart file if active + if (ESMF_AlarmIsRinging(stop_alarm, rc=rc)) then + write(restartname,'(A)')"ocn_stoch.res.nc") + else + write(restartname,'(A,I4.4,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,A)') & + "oc_stoch.res.", year, month, day, hour, minute, seconds,".nc" + endif + call ESMF_LogWrite("MOM_cap: Writing restart : "//trim(restartname), ESMF_LOGMSG_INFO) + call write_stoch_restart_ocn('RESTART/'//trim(timestamp)//'.ocn_stoch.res.nc') endif if (is_root_pe()) then From 25ed5ef652fe30284192c9eb6ae4ef5882516863 Mon Sep 17 00:00:00 2001 From: pjpegion Date: Tue, 22 Dec 2020 18:40:13 +0000 Subject: [PATCH 07/21] additions for stochy restarts --- config_src/nuopc_driver/mom_cap.F90 | 11 ++++++++--- config_src/solo_driver/MOM_driver.F90 | 3 +++ src/core/MOM.F90 | 14 +++++++------- .../vertical/MOM_diabatic_driver.F90 | 16 +++++++++++++--- .../vertical/MOM_energetic_PBL.F90 | 19 ++++++++++--------- 5 files changed, 41 insertions(+), 22 deletions(-) diff --git a/config_src/nuopc_driver/mom_cap.F90 b/config_src/nuopc_driver/mom_cap.F90 index 6880addfb2..12d70fe4e7 100644 --- a/config_src/nuopc_driver/mom_cap.F90 +++ b/config_src/nuopc_driver/mom_cap.F90 @@ -91,7 +91,9 @@ 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 implicit none; private @@ -1587,14 +1589,17 @@ 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") + write(restartname,'(A)')"ocn_stoch.res.nc" else write(restartname,'(A,I4.4,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,"-",I2.2,A)') & - "oc_stoch.res.", year, month, day, hour, minute, seconds,".nc" + "ocn_stoch.res.", year, month, day, hour, minute, seconds,".nc" endif call ESMF_LogWrite("MOM_cap: Writing restart : "//trim(restartname), ESMF_LOGMSG_INFO) - call write_stoch_restart_ocn('RESTART/'//trim(timestamp)//'.ocn_stoch.res.nc') + if (is_root_pe()) print*,'calling write_stoch_restart_ocn ',trim(restartname) + call write_stoch_restart_ocn('RESTART/'//trim(restartname)) +#endif endif if (is_root_pe()) then diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index ba52d9c02a..81fb79207f 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -74,6 +74,9 @@ program MOM_main use MOM_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init use MOM_wave_interface, only: MOM_wave_interface_init_lite, Update_Surface_Waves +#ifdef UFS + use get_stochy_pattern_mod, only: write_stoch_restart_ocn +#endif implicit none diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index a4f2f81af2..5143ffbf77 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -142,7 +142,9 @@ 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 @@ -229,8 +231,6 @@ module MOM logical :: offline_tracer_mode = .false. !< If true, step_offline() is called instead of step_MOM(). !! This is intended for running MOM6 in offline tracer mode - logical :: do_stochy = .false. - !< If true, call stochastic physics pattern generator type(time_type), pointer :: Time !< pointer to the ocean clock real :: dt !< (baroclinic) dynamics time step [T ~> s] @@ -2361,6 +2361,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call copy_dyngrid_to_MOM_grid(dG_in, G_in, US) call destroy_dyn_horgrid(dG_in) + do_epbl=.false. + do_sppt=.false. +#ifdef UFS num_procs=num_PEs() allocate(pelist(num_procs)) call Get_PElist(pelist,commID = mom_comm) @@ -2368,12 +2371,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & master=root_PE() !print*,'callling init_stochastic_physics_ocn',maxval(G%geoLatT) - do_epbl=.false. - do_sppt=.false. + if (master) print*,'about to call init_stochastic_physics' 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) - if (do_sppt .eq. .true.) CS%do_stochy=.true. - if (do_epbl .eq. .true.) CS%do_stochy=.true. - !print*,'back from init_stochastic_physics_ocn',CS%do_stochy +#endif ! Set a few remaining fields that are specific to the ocean grid type. call set_first_direction(G, first_direction) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 15d41293c1..78068ba44d 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -69,7 +69,9 @@ 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 +#ifdef UFS use stochastic_physics, only : run_stochastic_physics_ocn +#endif implicit none ; private @@ -292,23 +294,28 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & 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 + print*,'calling run_stochastic_physics' call run_stochastic_physics_ocn(t_rp,sppt_wts) !print*,'in diabatic',CS%do_sppt,size(t_in,1),size(t_in,2),size(t_in,3),size(sppt_wts,1),size(sppt_wts,2) - !print*,'in diabatic',CS%do_sppt,minval(sppt_wts),maxval(sppt_wts) + print*,'in diabatic',CS%do_sppt,minval(sppt_wts),maxval(sppt_wts) + print*,'in diabatic',CS%do_sppt,minval(t_rp),maxval(t_rp) if (CS%id_t_rp1 > 0) then call post_data(CS%id_t_rp1, t_rp(:,:,1), CS%diag) endif @@ -318,6 +325,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (CS%id_sppt_wts > 0) then call post_data(CS%id_sppt_wts, sppt_wts, CS%diag) endif +#endif is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -465,6 +473,7 @@ 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 do k=1,nz do j=js,je @@ -488,6 +497,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo enddo endif +#endif end subroutine diabatic @@ -891,7 +901,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, dt, Kd_ePBL, G, GV, US, & + call energetic_PBL(h, u_h, v_h, tv, fluxes, t_rp, CS%do_epbl, dt, Kd_ePBL, G, GV, US, & CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) if (associated(Hml)) then @@ -1623,7 +1633,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, dt, Kd_ePBL, G, GV, US, & + call energetic_PBL(h, u_h, v_h, tv, fluxes, t_rp, CS%do_epbl, dt, Kd_ePBL, G, GV, US, & CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) if (associated(Hml)) then diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 277f714f2a..3c82049a8b 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -246,7 +246,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, dt, Kd_int, G, GV, US, CS, & +subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, stoch_epbl, 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. @@ -285,6 +285,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, dt, Kd_int, G, GV, !! 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]. real, optional, intent(in) :: dt_diag !< The diagnostic time step, which may be less @@ -438,8 +439,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, dt, Kd_int, G, GV, do K=1,nz+1 ; Kd(K) = 0.0 ; enddo ! Make local copies of surface forcing and process them. - !print*,'PJP EPBL',minval(t_rp),maxval(t_rp) - u_star = fluxes%ustar(i,j)!*t_rp(i,j) + u_star = fluxes%ustar(i,j) u_star_Mean = fluxes%ustar_gustless(i,j) B_flux = buoy_flux(i,j) if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then @@ -463,7 +463,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, dt, Kd_int, G, GV, 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), dt_diag=dt_diag, Waves=Waves, G=G, i=i, j=j) + 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) ! applly stochastic perturbation to TKE generation @@ -548,7 +548,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, & - t_rp1,t_rp2, dt_diag, Waves, G, i, j) + t_rp1,t_rp2, stoch_epbl, 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]. @@ -587,7 +587,8 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs !! call to mixedlayer_init. 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 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), & @@ -888,8 +889,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 - mech_TKE=mech_TKE*t_rp1 + ! stochastically pertrub mech_TKE in the UFS + if (stoch_epbl) mech_TKE=mech_TKE*t_rp1 if (CS%TKE_diagnostics) then eCD%dTKE_conv = 0.0 ; eCD%dTKE_mixing = 0.0 @@ -973,7 +974,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs if (CS%TKE_diagnostics) & !eCD%dTKE_mech_decay = eCD%dTKE_mech_decay + (exp_kh-1.0) * mech_TKE * I_dtdiag eCD%dTKE_mech_decay = exp_kh - mech_TKE = mech_TKE * (1+(exp_kh-1) * t_rp2) + if (stoch_epbl) mech_TKE = mech_TKE * (1+(exp_kh-1) * t_rp2) ! Accumulate any convectively released potential energy to contribute ! to wstar and to drive penetrating convection. From 4bd9b9ed4a5e45356578693487009975cb5668e6 Mon Sep 17 00:00:00 2001 From: pjpegion Date: Wed, 23 Dec 2020 21:46:27 +0000 Subject: [PATCH 08/21] clean up debug statements --- config_src/nuopc_driver/mom_cap.F90 | 1 - config_src/solo_driver/MOM_driver.F90 | 3 --- src/core/MOM.F90 | 2 -- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 4 ---- 4 files changed, 10 deletions(-) diff --git a/config_src/nuopc_driver/mom_cap.F90 b/config_src/nuopc_driver/mom_cap.F90 index 12d70fe4e7..cbc4e2d82c 100644 --- a/config_src/nuopc_driver/mom_cap.F90 +++ b/config_src/nuopc_driver/mom_cap.F90 @@ -1597,7 +1597,6 @@ subroutine ModelAdvance(gcomp, rc) "ocn_stoch.res.", year, month, day, hour, minute, seconds,".nc" endif call ESMF_LogWrite("MOM_cap: Writing restart : "//trim(restartname), ESMF_LOGMSG_INFO) - if (is_root_pe()) print*,'calling write_stoch_restart_ocn ',trim(restartname) call write_stoch_restart_ocn('RESTART/'//trim(restartname)) #endif endif diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index 81fb79207f..ba52d9c02a 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -74,9 +74,6 @@ program MOM_main use MOM_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init use MOM_wave_interface, only: MOM_wave_interface_init_lite, Update_Surface_Waves -#ifdef UFS - use get_stochy_pattern_mod, only: write_stoch_restart_ocn -#endif implicit none diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 5143ffbf77..2e749f5240 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2370,8 +2370,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & me=PE_here() master=root_PE() - !print*,'callling init_stochastic_physics_ocn',maxval(G%geoLatT) - if (master) print*,'about to call init_stochastic_physics' 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 diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 78068ba44d..f1dac73e7a 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -311,11 +311,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & t_in=tv%T s_in=tv%S endif - print*,'calling run_stochastic_physics' call run_stochastic_physics_ocn(t_rp,sppt_wts) - !print*,'in diabatic',CS%do_sppt,size(t_in,1),size(t_in,2),size(t_in,3),size(sppt_wts,1),size(sppt_wts,2) - print*,'in diabatic',CS%do_sppt,minval(sppt_wts),maxval(sppt_wts) - print*,'in diabatic',CS%do_sppt,minval(t_rp),maxval(t_rp) if (CS%id_t_rp1 > 0) then call post_data(CS%id_t_rp1, t_rp(:,:,1), CS%diag) endif From 1d7ffa370417126a4f6af2aa7cc82a70d521cae7 Mon Sep 17 00:00:00 2001 From: pjpegion Date: Wed, 6 Jan 2021 15:35:18 +0000 Subject: [PATCH 09/21] clean up code --- src/core/MOM.F90 | 7 +++-- src/core/MOM_forcing_type.F90 | 26 ------------------- src/diagnostics/MOM_diagnostics.F90 | 2 +- .../vertical/MOM_energetic_PBL.F90 | 6 ++--- 4 files changed, 6 insertions(+), 35 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 5c44b08c0b..4a9cb3dbb5 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -911,7 +911,6 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS call enable_averages(CS%time_in_thermo_cycle, Time_local, CS%diag) call post_surface_thermo_diags(CS%sfc_IDs, G, GV, US, CS%diag, CS%time_in_thermo_cycle, & sfc_state_diag, CS%tv, ssh, CS%ave_ssh_ibc) - !sfc_state_diag, CS%tv, ssh,fluxes%t_rp,fluxes%sppt_wts, CS%ave_ssh_ibc) endif call disable_averaging(CS%diag) call cpu_clock_end(id_clock_diagnostics) @@ -1699,9 +1698,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & 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,iret -! model + 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 diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 930f54566d..6720414b2b 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -141,10 +141,6 @@ module MOM_forcing_type !< Pressure at the top ocean interface [R L2 T-2 ~> Pa] that is used in corrections to the sea surface !! height field that is passed back to the calling routines. !! p_surf_SSH may point to p_surf or to p_surf_full. -! real, pointer, dimension(:,:) :: t_rp => NULL() -! !< random pattern at t-points -! real, pointer, dimension(:,:) :: sppt_wts => NULL() -! !< random pattern at t-points logical :: accumulate_p_surf = .false. !< If true, the surface pressure due to the atmosphere !! and various types of ice needs to be accumulated, and the !! surface pressure explicitly reset to zero at the driver level @@ -242,10 +238,6 @@ module MOM_forcing_type !! u-points [L4 Z-1 T-1 ~> m3 s-1] rigidity_ice_v => NULL() !< Depth-integrated lateral viscosity of ice shelves or sea ice at !! v-points [L4 Z-1 T-1 ~> m3 s-1] -! real, pointer, dimension(:,:) :: t_rp => NULL() -! !< random pattern at t-points -! real, pointer, dimension(:,:) :: sppt_wts => NULL() -! !< random pattern at t-points real :: dt_force_accum = -1.0 !< The amount of time over which the mechanical forcing fluxes !! have been averaged [s]. logical :: net_mass_src_set = .false. !< If true, an estimate of net_mass_src has been provided. @@ -2086,18 +2078,6 @@ subroutine copy_common_forcing_fields(forces, fluxes, G, skip_pres) do_pres = .true. ; if (present(skip_pres)) do_pres = .not.skip_pres -! if (associated(forces%t_rp) .and. associated(fluxes%t_rp)) then -! do j=js,je ; do i=is,ie -! fluxes%t_rp(i,j) = forces%t_rp(i,j) -! enddo ; enddo -! endif -! -! if (associated(forces%sppt_wts) .and. associated(fluxes%sppt_wts)) then -! do j=js,je ; do i=is,ie -! fluxes%sppt_wts(i,j) = forces%sppt_wts(i,j) -! enddo ; enddo -! endif - if (associated(forces%ustar) .and. associated(fluxes%ustar)) then do j=js,je ; do i=is,ie fluxes%ustar(i,j) = forces%ustar(i,j) @@ -3041,8 +3021,6 @@ subroutine allocate_mech_forcing_by_group(G, forces, stress, ustar, shelf, & call myAlloc(forces%p_surf,isd,ied,jsd,jed, press) call myAlloc(forces%p_surf_full,isd,ied,jsd,jed, press) call myAlloc(forces%net_mass_src,isd,ied,jsd,jed, press) -! call myAlloc(forces%t_rp,isd,ied,jsd,jed, press) -! call myAlloc(forces%sppt_wts,isd,ied,jsd,jed, press) call myAlloc(forces%rigidity_ice_u,IsdB,IedB,jsd,jed, shelf) call myAlloc(forces%rigidity_ice_v,isd,ied,JsdB,JedB, shelf) @@ -3198,8 +3176,6 @@ subroutine deallocate_forcing_type(fluxes) if (associated(fluxes%seaice_melt)) deallocate(fluxes%seaice_melt) if (associated(fluxes%salt_flux)) deallocate(fluxes%salt_flux) if (associated(fluxes%p_surf_full)) deallocate(fluxes%p_surf_full) -! if (associated(fluxes%t_rp)) deallocate(fluxes%t_rp) -! if (associated(fluxes%sppt_wts)) deallocate(fluxes%sppt_wts) if (associated(fluxes%p_surf)) deallocate(fluxes%p_surf) if (associated(fluxes%TKE_tidal)) deallocate(fluxes%TKE_tidal) if (associated(fluxes%ustar_tidal)) deallocate(fluxes%ustar_tidal) @@ -3224,8 +3200,6 @@ subroutine deallocate_mech_forcing(forces) if (associated(forces%ustar)) deallocate(forces%ustar) if (associated(forces%p_surf)) deallocate(forces%p_surf) if (associated(forces%p_surf_full)) deallocate(forces%p_surf_full) -! if (associated(forces%t_rp)) deallocate(forces%t_rp) -! if (associated(forces%sppt_wts)) deallocate(forces%sppt_wts) if (associated(forces%net_mass_src)) deallocate(forces%net_mass_src) if (associated(forces%rigidity_ice_u)) deallocate(forces%rigidity_ice_u) if (associated(forces%rigidity_ice_v)) deallocate(forces%rigidity_ice_v) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index c481243fcf..39c7d8b5f5 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -177,7 +177,7 @@ module MOM_diagnostics integer :: id_salt_deficit = -1 integer :: id_Heat_PmE = -1 integer :: id_intern_heat = -1 - !!@} + !>@} end type surface_diag_IDs diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 3c82049a8b..57ba388f45 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -465,8 +465,6 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, stoch_epbl, dt, Kd_ 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) - ! applly stochastic perturbation to TKE generation - ! Copy the diffusivities to a 2-d array. do K=1,nz+1 Kd_2d(i,K) = Kd(K) @@ -972,8 +970,8 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs exp_kh = 1.0 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 - eCD%dTKE_mech_decay = exp_kh + eCD%dTKE_mech_decay = eCD%dTKE_mech_decay + (exp_kh-1.0) * mech_TKE * I_dtdiag + mech_TKE = mech_TKE * exp_kh if (stoch_epbl) mech_TKE = mech_TKE * (1+(exp_kh-1) * t_rp2) ! Accumulate any convectively released potential energy to contribute From 6bb9d0bc84308dec6966dd1801a5b452c463a389 Mon Sep 17 00:00:00 2001 From: pjpegion Date: Thu, 7 Jan 2021 15:42:13 +0000 Subject: [PATCH 10/21] fix non stochastic ePBL calculation --- src/parameterizations/vertical/MOM_energetic_PBL.F90 | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 57ba388f45..94771e23d1 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -971,8 +971,11 @@ 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 (stoch_epbl) mech_TKE = mech_TKE * (1+(exp_kh-1) * t_rp2) + if (stoch_epbl) then ! perturb the TKE destruction + mech_TKE = mech_TKE * (1+(exp_kh-1) * t_rp2) + else + mech_TKE = mech_TKE * exp_kh + endif ! Accumulate any convectively released potential energy to contribute ! to wstar and to drive penetrating convection. From 1727d9ab59cc8708555de8d669b85d99817f2472 Mon Sep 17 00:00:00 2001 From: pjpegion Date: Fri, 29 Jan 2021 19:40:10 +0000 Subject: [PATCH 11/21] re-write of stochastic code to remove CPP directives --- config_src/coupled_driver/ocean_model_MOM.F90 | 15 ++-- config_src/mct_driver/mom_ocean_model_mct.F90 | 15 ++-- config_src/nuopc_driver/mom_cap.F90 | 4 -- .../nuopc_driver/mom_ocean_model_nuopc.F90 | 45 ++++++++++-- config_src/solo_driver/MOM_driver.F90 | 13 ++-- src/core/MOM.F90 | 46 +++--------- src/core/MOM_variables.F90 | 7 ++ .../vertical/MOM_diabatic_driver.F90 | 70 +++++++------------ .../vertical/MOM_energetic_PBL.F90 | 43 +++++++----- 9 files changed, 127 insertions(+), 131 deletions(-) diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 6cb358cdcb..223c179bfb 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -44,7 +44,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 @@ -187,6 +187,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. @@ -580,12 +581,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) @@ -607,16 +608,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) @@ -633,7 +634,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/mct_driver/mom_ocean_model_mct.F90 b/config_src/mct_driver/mom_ocean_model_mct.F90 index 5a04739971..d2d93d4fd2 100644 --- a/config_src/mct_driver/mom_ocean_model_mct.F90 +++ b/config_src/mct_driver/mom_ocean_model_mct.F90 @@ -46,7 +46,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 @@ -185,6 +185,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/nuopc_driver/mom_cap.F90 b/config_src/nuopc_driver/mom_cap.F90 index 966bdacf33..e078c0d74f 100644 --- a/config_src/nuopc_driver/mom_cap.F90 +++ b/config_src/nuopc_driver/mom_cap.F90 @@ -91,9 +91,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 implicit none; private @@ -1589,7 +1587,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 @@ -1598,7 +1595,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/nuopc_driver/mom_ocean_model_nuopc.F90 b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 index 493762f4bc..ef62430342 100644 --- a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 +++ b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 @@ -43,7 +43,7 @@ module MOM_ocean_model_nuopc 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 @@ -62,6 +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 MOM_domains, only : root_PE,PE_here,Get_PElist,num_PEs +use stochastic_physics, only : init_stochastic_physics_ocn, run_stochastic_physics_ocn #include @@ -187,6 +189,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. @@ -248,6 +251,13 @@ 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,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" @@ -416,6 +426,21 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i endif + 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) @@ -585,17 +610,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) @@ -618,16 +649,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) @@ -644,7 +675,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/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index ba52d9c02a..dd1cee96d6 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -57,7 +57,7 @@ program MOM_main use MOM_time_manager, only : 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_write_cputime, only : write_cputime, MOM_write_cputime_init use MOM_write_cputime, only : write_cputime_start_clock, write_cputime_CS @@ -84,6 +84,7 @@ 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 @@ -500,7 +501,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) @@ -513,16 +514,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) @@ -531,7 +532,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 4a9cb3dbb5..22abd99d25 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -29,7 +29,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 @@ -124,7 +123,7 @@ module MOM use MOM_variables, only : surface, allocate_surface_state, deallocate_surface_state use MOM_variables, only : thermo_var_ptrs, vertvisc_type 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 @@ -142,9 +141,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 @@ -420,13 +416,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), pointer :: CS !< control structure from initialize_MOM @@ -706,8 +703,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. @@ -804,8 +801,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS if (dtdia > dt) 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 @@ -1212,8 +1209,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 @@ -1226,6 +1223,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. @@ -1287,7 +1285,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, & + 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. @@ -1689,7 +1687,6 @@ 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. This can be altered during the course @@ -1697,12 +1694,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & integer :: nkml, nkbl, verbosity, write_geom integer :: dynamics_stencil ! The computational stencil for the calculations ! in the dynamic core. - 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 real :: RL2_T2_rescale, Z_rescale, QRZ_rescale ! Unit conversion factors character(len=48) :: flux_units, S_flux_units @@ -2360,18 +2351,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call copy_dyngrid_to_MOM_grid(dG_in, G_in, US) call destroy_dyn_horgrid(dG_in) - 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. call set_first_direction(G, first_direction) ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes. @@ -2784,9 +2763,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 0b225f0bf7..c4ee1eefb2 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -269,6 +269,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 f1dac73e7a..640b97ce95 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 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 @@ -210,8 +207,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 @@ -259,7 +254,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, & +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 @@ -271,6 +266,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & 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(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 @@ -292,36 +288,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 @@ -403,10 +387,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, & @@ -469,14 +453,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) @@ -493,7 +476,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo enddo endif -#endif end subroutine diabatic @@ -501,7 +483,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 @@ -513,8 +495,8 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, t_rp, visc, ADp, CDp, d !! unused have NULL ptrs 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 @@ -897,7 +879,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 @@ -1283,7 +1265,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 @@ -1296,7 +1278,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 @@ -1629,7 +1611,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 @@ -3444,12 +3426,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 94771e23d1..044dadec6a 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 @@ -163,13 +163,11 @@ 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. 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]. @@ -197,6 +195,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 @@ -246,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. @@ -277,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]. real, optional, intent(in) :: dt_diag !< The diagnostic time step, which may be less @@ -461,9 +459,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 @@ -534,6 +532,12 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, t_rp, stoch_epbl, dt, Kd_ 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 if (debug) deallocate(eCD%dT_expect, eCD%dS_expect) @@ -544,9 +548,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]. @@ -565,6 +569,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]. @@ -584,9 +589,6 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs type(energetic_PBL_CS), pointer :: CS !< The control structure returned by a previous !! call to mixedlayer_init. 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), & @@ -888,7 +890,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 @@ -971,8 +973,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 @@ -2385,7 +2387,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') From 5443f8effb07f2c06c812ed26588d2aaf90a974a Mon Sep 17 00:00:00 2001 From: pjpegion Date: Fri, 29 Jan 2021 20:34:07 +0000 Subject: [PATCH 12/21] remove blank link in MOM_diagnostics --- src/diagnostics/MOM_diagnostics.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 39c7d8b5f5..18c9fc96c4 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -1980,7 +1980,6 @@ subroutine register_surface_diags(Time, G, US, IDs, diag, tv) 'Heat flux into ocean from geothermal or other internal sources', & 'W m-2', conversion=US%QRZ_T_to_W_m2) - end subroutine register_surface_diags !> Register certain diagnostics related to transports From 80f9f44fa5abf829d2c0cc4bf58d055b5a95f2ba Mon Sep 17 00:00:00 2001 From: pjpegion Date: Fri, 29 Jan 2021 20:49:07 +0000 Subject: [PATCH 13/21] clean up MOM_domains --- src/framework/MOM_domains.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index 755507838d..06249daf6d 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -5,8 +5,6 @@ module MOM_domains use MOM_array_transform, only : rotate_array -use MOM_coms, only : PE_here, root_PE, num_PEs, MOM_infra_init, MOM_infra_end, Get_PElist - use MOM_coms, only : PE_here, root_PE, num_PEs, MOM_infra_init, MOM_infra_end, Get_PElist use MOM_coms, only : broadcast, sum_across_PEs, min_across_PEs, max_across_PEs use MOM_cpu_clock, only : cpu_clock_begin, cpu_clock_end From 0b99c1f82f562d4b0231d49353c6ee4b9204dc3d Mon Sep 17 00:00:00 2001 From: pjpegion Date: Tue, 2 Feb 2021 00:05:24 +0000 Subject: [PATCH 14/21] make stochastics optional --- config_src/coupled_driver/ocean_model_MOM.F90 | 17 ++-- config_src/mct_driver/mom_ocean_model_mct.F90 | 15 ++-- .../nuopc_driver/mom_ocean_model_nuopc.F90 | 78 ++++++++++++------- config_src/solo_driver/MOM_driver.F90 | 13 ++-- src/core/MOM.F90 | 25 +++--- src/core/MOM_variables.F90 | 2 - .../vertical/MOM_diabatic_driver.F90 | 60 ++++++++------ .../vertical/MOM_energetic_PBL.F90 | 38 +++++---- 8 files changed, 145 insertions(+), 103 deletions(-) diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index b091af9bb0..6cb358cdcb 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -44,7 +44,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, stochastic_pattern +use MOM_variables, only : surface 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,7 +187,6 @@ 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. @@ -566,7 +565,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda ! For now, the waves are only updated on the thermodynamics steps, because that is where ! the wave intensities are actually used to drive mixing. At some point, the wave updates ! might also need to become a part of the ocean dynamics, according to B. Reichl. - call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves) + call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves, OS%forces) endif if ((OS%nstep==0) .and. (OS%nstep_thermo==0)) then ! This is the first call to update_ocean_model. @@ -581,12 +580,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, OS%stochastics, Time1, dt_coupling, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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, OS%stochastics, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves) + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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) @@ -608,16 +607,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, OS%stochastics, Time1, dtdia, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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, OS%stochastics, Time1, dt_dyn, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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, OS%stochastics, Time1, dt_dyn, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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) @@ -634,7 +633,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, OS%stochastics, Time1, dtdia, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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/mct_driver/mom_ocean_model_mct.F90 b/config_src/mct_driver/mom_ocean_model_mct.F90 index d2d93d4fd2..5a04739971 100644 --- a/config_src/mct_driver/mom_ocean_model_mct.F90 +++ b/config_src/mct_driver/mom_ocean_model_mct.F90 @@ -46,7 +46,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, stochastic_pattern +use MOM_variables, only : surface 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 @@ -185,7 +185,6 @@ 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. @@ -587,12 +586,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, OS%stochastics, Time1, dt_coupling, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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, OS%stochastics, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves) + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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) @@ -616,16 +615,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, OS%stochastics, Time2, dtdia, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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, OS%stochastics, Time2, dt_dyn, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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, OS%stochastics, Time2, dt_dyn, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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) @@ -642,7 +641,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, OS%stochastics, Time2, dtdia, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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/nuopc_driver/mom_ocean_model_nuopc.F90 b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 index ef62430342..85624b95b8 100644 --- a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 +++ b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 @@ -176,6 +176,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 real :: eps_omesh !< Max allowable difference between ESMF mesh and MOM6 !! domain coordinates @@ -426,20 +428,38 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i endif - 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)) +! get number of processors and PE list for stocasthci physics initialization + call get_param(param_file, mdl, "DO_SPPT", OS%do_sppt, & + "If true, then stochastically perturb the thermodynamic "//& + "tendemcies of T,S, amd h. Amplitude and correlations are "//& + "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.) + if (OS%do_sppt .OR. OS%pert_epbl) then + 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%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%stochastics%sppt_wts(OS%grid%isd:OS%grid%ied,OS%grid%jsd:OS%grid%jed)) + if (OS%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 endif call close_param_file(param_file) call diag_mediator_close_registration(OS%diag) @@ -611,8 +631,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & 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 + if (OS%do_sppt .OR. OS%pert_epbl ) then call run_stochastic_physics_ocn(OS%stochastics%sppt_wts,OS%stochastics%t_rp1,OS%stochastics%t_rp2) endif @@ -620,13 +639,14 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & 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, OS%stochastics, Time1, dt_coupling, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=do_thermo, do_thermodynamics=do_dyn, & - reset_therm=Ocn_fluxes_used) + reset_therm=Ocn_fluxes_used, stochastics=OS%stochastics) !### 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, OS%stochastics, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves) + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves, & + stochastics=OS%stochastics) 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) @@ -649,18 +669,21 @@ 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, OS%stochastics, Time2, dtdia, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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) + start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling, & + stochastics=OS%stochastics) endif - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time2, dt_dyn, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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) + start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling, & + stochastics=OS%stochastics) else - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, OS%stochastics, Time2, dt_dyn, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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) + start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling, & + stochastics=OS%stochastics) step_thermo = .false. if (thermo_does_span_coupling) then @@ -675,9 +698,10 @@ 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, OS%stochastics, Time2, dtdia, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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) + start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling, & + stochastics=OS%stochastics) endif endif diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index dd1cee96d6..ba52d9c02a 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -57,7 +57,7 @@ program MOM_main use MOM_time_manager, only : 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, stochastic_pattern + use MOM_variables, only : surface use MOM_verticalGrid, only : verticalGrid_type use MOM_write_cputime, only : write_cputime, MOM_write_cputime_init use MOM_write_cputime, only : write_cputime_start_clock, write_cputime_CS @@ -84,7 +84,6 @@ 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 @@ -501,7 +500,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, stochastics, Time1, dt_forcing, MOM_CSp, Waves=Waves_CSP) + call step_MOM(forces, fluxes, sfc_state, 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) @@ -514,16 +513,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, stochastics, Time2, dtdia, MOM_CSp, & + call step_MOM(forces, fluxes, sfc_state, 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, stochastics, Time2, dt_dyn, MOM_CSp, & + call step_MOM(forces, fluxes, sfc_state, 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, stochastics, Time2, dt_dyn, MOM_CSp, & + call step_MOM(forces, fluxes, sfc_state, Time2, dt_dyn, MOM_CSp, & do_dynamics=.true., do_thermodynamics=.false., & start_cycle=(n==1), end_cycle=.false., cycle_length=dt_forcing) @@ -532,7 +531,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, stochastics, Time2, dtdia, MOM_CSp, & + call step_MOM(forces, fluxes, sfc_state, 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 22abd99d25..6798347ba5 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -416,14 +416,13 @@ 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, stochastics, Time_start, time_int_in, CS, & +subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS, & Waves, do_dynamics, do_thermodynamics, start_cycle, & - end_cycle, cycle_length, reset_therm) + end_cycle, cycle_length, reset_therm, stochastics) 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), pointer :: CS !< control structure from initialize_MOM @@ -444,6 +443,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, stochastics, Time_start, ti logical, optional, intent(in) :: reset_therm !< This indicates whether the running sums of !! thermodynamic quantities should be reset. !! If missing, this is like start_cycle. + type(stochastic_pattern), optional, intent(in) :: stochastics !< random patternss for stochastics ! local variables type(ocean_grid_type), pointer :: G => NULL() ! pointer to a structure containing @@ -703,8 +703,9 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, stochastics, Time_start, ti endif ! Apply diabatic forcing, do mixing, and regrid. - call step_MOM_thermo(CS, G, GV, US, u, v, h, CS%tv, fluxes, stochastics, & - dtdia, end_time_thermo, .true., Waves=Waves) + call step_MOM_thermo(CS, G, GV, US, u, v, h, CS%tv, fluxes, dtdia, & + end_time_thermo, .true., Waves=Waves, & + stochastics=stochastics) CS%time_in_thermo_cycle = CS%time_in_thermo_cycle + dtdia ! The diabatic processes are now ahead of the dynamics by dtdia. @@ -801,8 +802,9 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, stochastics, Time_start, ti if (dtdia > dt) 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, stochastics, & - dtdia, Time_local, .false., Waves=Waves) + call step_MOM_thermo(CS, G, GV, US, u, v, h, CS%tv, fluxes, dtdia, & + Time_local, .false., Waves=Waves, & + stochastics=stochastics) CS%time_in_thermo_cycle = CS%time_in_thermo_cycle + dtdia if ((CS%t_dyn_rel_thermo==0.0) .and. .not.do_dyn) then @@ -1209,8 +1211,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, stochastics, & - dtdia, Time_end_thermo, update_BBL, Waves) +subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & + Time_end_thermo, update_BBL, Waves, stochastics) 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 @@ -1285,8 +1287,9 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, stochastics, & call cpu_clock_begin(id_clock_diabatic) - 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) + 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, & + stochastics=stochastics) fluxes%fluxes_used = .true. if (showCallTree) call callTree_waypoint("finished diabatic (step_MOM_thermo)") diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index c4ee1eefb2..2881233767 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -270,8 +270,6 @@ 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 diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 640b97ce95..28d82e9143 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -141,6 +141,8 @@ module MOM_diabatic_driver !! 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 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,7 +177,7 @@ module MOM_diabatic_driver 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, id_t_rp1=-1,id_t_rp2=-1 + integer :: id_subMLN2 = -1, id_sppt_wts = -1, id_t_rp1 = -1, id_t_rp2 = -1 ! diagnostic for fields prior to applying diapycnal physics integer :: id_u_predia = -1, id_v_predia = -1, id_h_predia = -1 @@ -254,8 +256,8 @@ 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, stochastics, visc, ADp, CDp, dt, Time_end, & - G, GV, US, CS, OBC, WAVES) +subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & + G, GV, US, CS, OBC, WAVES, stochastics) 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] @@ -266,7 +268,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, stochastics, visc, ADp, CDp, dt, T 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(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 @@ -278,6 +279,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, stochastics, visc, ADp, CDp, dt, T 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(stochastic_pattern), optional, intent(in) :: stochastics !< random patterns ! local variables real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: & @@ -288,16 +290,16 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, stochastics, visc, ADp, CDp, dt, T integer :: i, j, k, m, is, ie, js, je, nz logical :: showCallTree ! If true, show the call tree - 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, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_in ! thickenss before thermodynamics + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: t_in ! temperature before thermodynamics + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: 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 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 (stochastics%do_sppt) then + if (CS%do_sppt) then h_in=h t_in=tv%T s_in=tv%S @@ -387,11 +389,11 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, stochastics, visc, ADp, CDp, dt, T 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, stochastics, visc, ADp, CDp, dt, Time_end, & - G, GV, US, CS, Waves) + call diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & + G, GV, US, CS, Waves,stochastics=stochastics) elseif (CS%useALEalgorithm) then - call diabatic_ALE(u, v, h, tv, Hml, fluxes, stochastics, visc, ADp, CDp, dt, Time_end, & - G, GV, US, CS, Waves) + call diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & + G, GV, US, CS, Waves, stochastics=stochastics) else call layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & G, GV, US, CS, Waves) @@ -453,7 +455,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, stochastics, visc, ADp, CDp, dt, T if (CS%debugConservation) call MOM_state_stats('leaving diabatic', u, v, h, tv%T, tv%S, G, GV, US) - if (stochastics%do_sppt) then + if (CS%do_sppt) then do k=1,nz do j=js,je do i=is,ie @@ -483,8 +485,8 @@ 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, stochastics, visc, ADp, CDp, dt, Time_end, & - G, GV, US, CS, WAVES) +subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & + G, GV, US, CS, WAVES, stochastics) 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 @@ -495,8 +497,6 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, stochastics, visc, ADp, !! unused have NULL ptrs 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 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 @@ -506,6 +506,8 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, stochastics, visc, ADp, 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(stochastic_pattern), optional, intent(in) :: stochastics !< random patterns for SPPT and + !! energetic PBL perturbations ! local variables real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & @@ -879,8 +881,9 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, stochastics, visc, ADp, 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, stochastics, dt, Kd_ePBL, G, GV, US, & - CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) + 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, stochastics=stochastics) if (associated(Hml)) then call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, Hml(:,:), G, US) @@ -1265,8 +1268,8 @@ 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, stochastics, visc, ADp, CDp, dt, Time_end, & - G, GV, US, CS, Waves) +subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & + G, GV, US, CS, Waves, stochastics) 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 @@ -1278,7 +1281,6 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, stochastics, visc, ADp, CDp, d 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(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 @@ -1288,6 +1290,8 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, stochastics, visc, ADp, CDp, d 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(stochastic_pattern), optional, intent(in) :: stochastics !< random patterns for SPPT and + !! energetic PBL perturbations ! local variables real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & @@ -1611,8 +1615,9 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, stochastics, 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, stochastics, dt, Kd_ePBL, G, GV, US, & - CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) + 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, stochastics=stochastics) if (associated(Hml)) then call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, Hml(:,:), G, US) @@ -3401,6 +3406,11 @@ 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.) ! 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 044dadec6a..856cf13598 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -56,6 +56,7 @@ 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. @@ -245,9 +246,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, stochastics, dt, Kd_int, G, GV, US, CS, & +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 ) + dT_expected, dS_expected, Waves, stochastics ) 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 @@ -276,8 +277,6 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, stochastics, dt, Kd_int, 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 @@ -302,6 +301,10 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, stochastics, dt, Kd_int, !! diffusivities are applied [ppt]. type(wave_parameters_CS), & optional, pointer :: Waves !< Wave CS + type(stochastic_pattern), optional, & + intent(in) :: stochastics !< A structure containing array to stochastic + !! patterns. Any unsued fields + !! are not allocated ! This subroutine determines the diffusivities from the integrated energetics ! mixed layer model. It assumes that heating, cooling and freshwater fluxes @@ -459,9 +462,10 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, stochastics, dt, Kd_int, 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, stochastics, & + 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) + GV, US, CS, eCD, dt_diag=dt_diag, Waves=Waves, G=G, & + stochastics=stochastics,i=i, j=j) ! Copy the diffusivities to a 2-d array. do K=1,nz+1 @@ -534,7 +538,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, stochastics, dt, Kd_int, 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%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 @@ -548,9 +552,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, stochastics, B_flux, absf, & +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, stochastics, 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]. @@ -569,7 +573,6 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, stochastics 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]. @@ -595,6 +598,8 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, stochastics optional, pointer :: Waves !< Wave CS for Langmuir turbulence type(ocean_grid_type), & optional, intent(inout) :: G !< The ocean's grid structure. + type(stochastic_pattern), & + optional, intent(in) :: stochastics !< stochastic patterns and logic controls 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) @@ -890,7 +895,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, stochastics mech_TKE = MSTAR_total * (dt*GV%Rho0* u_star**3) endif ! stochastically pertrub mech_TKE in the UFS - if (stochastics%pert_epbl) mech_TKE=mech_TKE*stochastics%t_rp1(i,j) + if (CS%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 @@ -973,7 +978,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, stochastics 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 (stochastics%pert_epbl) then ! perturb the TKE destruction + if (CS%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 @@ -2213,6 +2218,11 @@ 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. "//& @@ -2388,9 +2398,9 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) 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') + 'random pattern for KE generation', 'None') CS%id_t_rp2 = register_diag_field('ocean_model', 'random_pattern2', diag%axesT1, Time, & - 'random pattern2 for stochastics', 'None') + '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') From 6e3ea1b007fd798ca6c769054ffdad9f976cf781 Mon Sep 17 00:00:00 2001 From: pjpegion Date: Tue, 2 Feb 2021 00:13:24 +0000 Subject: [PATCH 15/21] correct coupled_driver/ocean_model_MOM.F90 and other cleanup --- config_src/coupled_driver/ocean_model_MOM.F90 | 2 +- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 4 ++-- src/parameterizations/vertical/MOM_energetic_PBL.F90 | 6 +++--- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 6cb358cdcb..082099158c 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -565,7 +565,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda ! For now, the waves are only updated on the thermodynamics steps, because that is where ! the wave intensities are actually used to drive mixing. At some point, the wave updates ! might also need to become a part of the ocean dynamics, according to B. Reichl. - call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves, OS%forces) + call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves) endif if ((OS%nstep==0) .and. (OS%nstep_thermo==0)) then ! This is the first call to update_ocean_model. diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 28d82e9143..aef8955ff8 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -290,7 +290,7 @@ 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),SZK_(G)) :: h_in ! thickenss before thermodynamics + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_in ! thickness before thermodynamics real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: t_in ! temperature before thermodynamics real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: s_in ! salinity before thermodynamics real :: t_tend,s_tend,h_tend ! holder for tendencey needed for SPPT @@ -390,7 +390,7 @@ 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,stochastics=stochastics) + G, GV, US, CS, Waves, stochastics=stochastics) elseif (CS%useALEalgorithm) then call diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & G, GV, US, CS, Waves, stochastics=stochastics) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 856cf13598..2dc9df57cc 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -462,9 +462,9 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS 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, & + 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, & stochastics=stochastics,i=i, j=j) ! Copy the diffusivities to a 2-d array. From eb88219af19e4624cdb379452a07821e276f0ed1 Mon Sep 17 00:00:00 2001 From: Philip Pegion Date: Tue, 2 Feb 2021 09:45:46 -0600 Subject: [PATCH 16/21] clean up of code for MOM6 coding standards --- .../nuopc_driver/mom_ocean_model_nuopc.F90 | 8 +++--- src/core/MOM.F90 | 2 +- .../vertical/MOM_diabatic_driver.F90 | 27 ++++++++++--------- .../vertical/MOM_energetic_PBL.F90 | 10 +++---- 4 files changed, 24 insertions(+), 23 deletions(-) diff --git a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 index 85624b95b8..2b19e7fc8f 100644 --- a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 +++ b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 @@ -455,10 +455,10 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i return endif - if (OS%do_sppt) allocate(OS%stochastics%sppt_wts(OS%grid%isd:OS%grid%ied,OS%grid%jsd:OS%grid%jed)) + if (OS%do_sppt) allocate(OS%stochastics%sppt_wts(OS%grid%isd:OS%grid%ied, OS%grid%jsd:OS%grid%jed)) if (OS%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)) + 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 endif call close_param_file(param_file) @@ -632,7 +632,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & ! update stochastic physics patterns before running next time-step if (OS%do_sppt .OR. OS%pert_epbl ) then - call run_stochastic_physics_ocn(OS%stochastics%sppt_wts,OS%stochastics%t_rp1,OS%stochastics%t_rp2) + call run_stochastic_physics_ocn(OS%stochastics%sppt_wts, OS%stochastics%t_rp1, OS%stochastics%t_rp2) endif if (OS%offline_tracer_mode) then diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 6798347ba5..c4facf7b46 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -123,7 +123,7 @@ module MOM use MOM_variables, only : surface, allocate_surface_state, deallocate_surface_state use MOM_variables, only : thermo_var_ptrs, vertvisc_type use MOM_variables, only : accel_diag_ptrs, cont_diag_ptrs, ocean_internal_state -use MOM_variables, only : rotate_surface_state,stochastic_pattern +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 diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index aef8955ff8..27f49a01c9 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -300,9 +300,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! save copy of the date for SPPT if (CS%do_sppt) then - 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, stochastics%sppt_wts, CS%diag) @@ -456,23 +456,24 @@ 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 + ! 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))*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) + 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) if (h_pert > GV%Angstrom_H) then - h(i,j,k)=h_pert + h(i,j,k) = h_pert else - h(i,j,k)=GV%Angstrom_H + h(i,j,k) = GV%Angstrom_H endif - tv%T(i,j,k)=t_pert + tv%T(i,j,k) = t_pert if (s_pert > 0.0) then - tv%S(i,j,k)=s_pert + tv%S(i,j,k) = s_pert endif enddo enddo diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 2dc9df57cc..f0cb692453 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -196,7 +196,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 + integer :: id_t_rp1=-1, id_t_rp2=-1 !>@} end type energetic_PBL_CS @@ -539,8 +539,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS ! 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_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) + 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 @@ -895,7 +895,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*stochastics%t_rp1(i,j) + if (CS%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 @@ -978,7 +978,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 (CS%pert_epbl) then ! perturb the TKE dissipation mech_TKE = mech_TKE * (1+(exp_kh-1) * stochastics%t_rp2(i,j)) else mech_TKE = mech_TKE * exp_kh From d984a7e1b1f3c2b4a8e146448791bf42ebff1f57 Mon Sep 17 00:00:00 2001 From: pjpegion Date: Thu, 4 Feb 2021 16:59:22 +0000 Subject: [PATCH 17/21] remove stochastics container --- .../nuopc_driver/mom_ocean_model_nuopc.F90 | 37 ++++++--------- src/core/MOM.F90 | 15 ++---- src/core/MOM_forcing_type.F90 | 4 ++ src/core/MOM_variables.F90 | 5 -- src/framework/MOM_domains.F90 | 2 +- .../vertical/MOM_diabatic_driver.F90 | 46 +++++++++---------- .../vertical/MOM_energetic_PBL.F90 | 43 ++++++++--------- 7 files changed, 69 insertions(+), 83 deletions(-) diff --git a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 index 85624b95b8..6b5a141a5e 100644 --- a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 +++ b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 @@ -43,7 +43,7 @@ module MOM_ocean_model_nuopc 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, stochastic_pattern +use MOM_variables, only : surface 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 @@ -51,7 +51,7 @@ module MOM_ocean_model_nuopc use coupler_types_mod, only : coupler_type_spawn, coupler_type_write_chksums use coupler_types_mod, only : coupler_type_initialized, coupler_type_copy_data use coupler_types_mod, only : coupler_type_set_diags, coupler_type_send_data -use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain +use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain,mpp_get_pelist use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain use fms_mod, only : stdout use mpp_mod, only : mpp_chksum @@ -62,7 +62,7 @@ module MOM_ocean_model_nuopc use MOM_surface_forcing_nuopc, only : convert_IOB_to_forces, ice_ocn_bnd_type_chksum use MOM_surface_forcing_nuopc, only : ice_ocean_boundary_type, surface_forcing_CS use MOM_surface_forcing_nuopc, only : forcing_save_restart -use MOM_domains, only : root_PE,PE_here,Get_PElist,num_PEs +use MOM_domains, only : root_PE,PE_here,num_PEs use stochastic_physics, only : init_stochastic_physics_ocn, run_stochastic_physics_ocn #include @@ -191,7 +191,6 @@ 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. @@ -254,7 +253,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i !! If HFrz <= 0 (default), melt potential will not be computed. 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 @@ -441,8 +439,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i 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) + call mpp_get_pelist(Ocean_sfc%domain, mom_comm) me=PE_here() master=root_PE() @@ -455,10 +452,10 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i return endif - if (OS%do_sppt) allocate(OS%stochastics%sppt_wts(OS%grid%isd:OS%grid%ied,OS%grid%jsd:OS%grid%jed)) + 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%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)) + 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) @@ -632,7 +629,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & ! update stochastic physics patterns before running next time-step if (OS%do_sppt .OR. OS%pert_epbl ) then - call run_stochastic_physics_ocn(OS%stochastics%sppt_wts,OS%stochastics%t_rp1,OS%stochastics%t_rp2) + call run_stochastic_physics_ocn(OS%fluxes%sppt_wts,OS%fluxes%epbl1_wts,OS%fluxes%epbl2_wts) endif if (OS%offline_tracer_mode) then @@ -641,12 +638,11 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & ! 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, & Waves=OS%Waves, do_dynamics=do_thermo, do_thermodynamics=do_dyn, & - reset_therm=Ocn_fluxes_used, stochastics=OS%stochastics) + 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, & - stochastics=OS%stochastics) + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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) @@ -671,19 +667,16 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & dtdia = dt_dyn*min(nts,n_max-(n-1)) call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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, & - stochastics=OS%stochastics) + 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, & Waves=OS%Waves, do_dynamics=.true., do_thermodynamics=.false., & - start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling, & - stochastics=OS%stochastics) + 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, & Waves=OS%Waves, do_dynamics=.true., do_thermodynamics=.false., & - start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling, & - stochastics=OS%stochastics) + start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling) step_thermo = .false. if (thermo_does_span_coupling) then @@ -700,8 +693,8 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & 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, & Waves=OS%Waves, do_dynamics=.false., do_thermodynamics=.true., & - start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling, & - stochastics=OS%stochastics) + start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling) + endif endif diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 6798347ba5..af56eb4c82 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -418,7 +418,7 @@ module MOM !! occur inside of diabatic. subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS, & Waves, do_dynamics, do_thermodynamics, start_cycle, & - end_cycle, cycle_length, reset_therm, stochastics) + 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 @@ -443,7 +443,6 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS logical, optional, intent(in) :: reset_therm !< This indicates whether the running sums of !! thermodynamic quantities should be reset. !! If missing, this is like start_cycle. - type(stochastic_pattern), optional, intent(in) :: stochastics !< random patternss for stochastics ! local variables type(ocean_grid_type), pointer :: G => NULL() ! pointer to a structure containing @@ -704,8 +703,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS ! 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, & - stochastics=stochastics) + 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. @@ -803,8 +801,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS ! 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, & - stochastics=stochastics) + 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 @@ -1212,7 +1209,7 @@ 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, stochastics) + 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 @@ -1225,7 +1222,6 @@ 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. @@ -1288,8 +1284,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, & - stochastics=stochastics) + 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)") diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index f0cc8f553c..2b0578ef49 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -166,6 +166,10 @@ 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] + ! stochastic patterns + real, allocatable :: sppt_wts(:,:) !< Random pattern for ocean SPPT + real, allocatable :: epbl1_wts(:,:) !< Random pattern for K.E. generation + real, allocatable :: epbl2_wts(:,:) !< Random pattern for K.E. dissipation ! 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/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 2881233767..0b225f0bf7 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -269,11 +269,6 @@ 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 - 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/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index 06249daf6d..33cb45814c 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -38,7 +38,7 @@ module MOM_domains public :: MOM_domains_init, MOM_infra_init, MOM_infra_end, get_domain_extent, get_domain_extent_dsamp2 public :: MOM_define_domain, MOM_define_io_domain, clone_MOM_domain -public :: pass_var, pass_vector, PE_here, root_PE, num_PEs, Get_PElist +public :: pass_var, pass_vector, PE_here, root_PE, num_PEs public :: pass_var_start, pass_var_complete, fill_symmetric_edges, broadcast public :: pass_vector_start, pass_vector_complete public :: global_field_sum, sum_across_PEs, min_across_PEs, max_across_PEs diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index aef8955ff8..88ee7a5dcb 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -65,7 +65,7 @@ module MOM_diabatic_driver use MOM_tracer_diabatic, only : tracer_vertdiff 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, stochastic_pattern +use MOM_variables, only : cont_diag_ptrs, MOM_thermovar_chksum, p3d use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units use MOM_wave_speed, only : wave_speeds use MOM_wave_interface, only : wave_parameters_CS @@ -177,7 +177,7 @@ module MOM_diabatic_driver 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, id_t_rp1 = -1, id_t_rp2 = -1 + integer :: id_subMLN2 = -1, id_sppt_wts = -1 ! diagnostic for fields prior to applying diapycnal physics integer :: id_u_predia = -1, id_v_predia = -1, id_h_predia = -1 @@ -257,7 +257,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, stochastics) + 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] @@ -279,7 +279,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & 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(stochastic_pattern), optional, intent(in) :: stochastics !< random patterns ! local variables real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: & @@ -290,9 +289,9 @@ 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),SZK_(G)) :: h_in ! thickness before thermodynamics - real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: t_in ! temperature before thermodynamics - real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: s_in ! salinity before thermodynamics + real, allocatable(:,:,:) :: h_in ! thickness before thermodynamics + real, allocatable(:,:,:) :: t_in ! temperature before thermodynamics + real, allocatable(:,:,:) :: 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 @@ -300,12 +299,15 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! save copy of the date for SPPT if (CS%do_sppt) then - h_in=h - t_in=tv%T - s_in=tv%S + 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 (CS%id_sppt_wts > 0) then - call post_data(CS%id_sppt_wts, stochastics%sppt_wts, CS%diag) + call post_data(CS%id_sppt_wts, fluxes%sppt_wts, CS%diag) endif endif @@ -390,10 +392,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, stochastics=stochastics) + G, GV, US, 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, stochastics=stochastics) + G, GV, US, CS, Waves) else call layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & G, GV, US, CS, Waves) @@ -459,9 +461,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & do k=1,nz do j=js,je do i=is,ie - 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_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_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) @@ -486,7 +488,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, visc, ADp, CDp, dt, Time_end, & - G, GV, US, CS, WAVES, stochastics) + 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 @@ -506,8 +508,6 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim 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(stochastic_pattern), optional, intent(in) :: stochastics !< random patterns for SPPT and - !! energetic PBL perturbations ! local variables real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & @@ -883,7 +883,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, & - waves=waves, stochastics=stochastics) + waves=waves) if (associated(Hml)) then call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, Hml(:,:), G, US) @@ -1269,7 +1269,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, stochastics) + 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 @@ -1290,8 +1290,6 @@ subroutine diabatic_ALE(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(diabatic_CS), pointer :: CS !< module control structure type(Wave_parameters_CS), optional, pointer :: Waves !< Surface gravity waves - type(stochastic_pattern), optional, intent(in) :: stochastics !< random patterns for SPPT and - !! energetic PBL perturbations ! local variables real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & @@ -1617,7 +1615,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, & - waves=waves, stochastics=stochastics) + 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 2dc9df57cc..b45f985a6b 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, stochastic_pattern +use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use MOM_wave_interface, only: wave_parameters_CS, Get_Langmuir_Number @@ -196,7 +196,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 + integer :: id_epbl1_wts=-1,id_epbl2_wts=-1 !>@} end type energetic_PBL_CS @@ -248,7 +248,7 @@ module MOM_energetic_PBL !! 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, stochastics ) + 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,10 +301,6 @@ 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_pattern), optional, & - intent(in) :: stochastics !< A structure containing array to stochastic - !! patterns. Any unsued fields - !! are not allocated ! This subroutine determines the diffusivities from the integrated energetics ! mixed layer model. It assumes that heating, cooling and freshwater fluxes @@ -461,11 +457,16 @@ 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, & - stochastics=stochastics,i=i, j=j) + if (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) + 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 @@ -539,8 +540,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS ! 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_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) + 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) endif endif @@ -554,7 +555,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, stochastics, 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]. @@ -598,8 +599,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. - type(stochastic_pattern), & - optional, intent(in) :: stochastics !< stochastic patterns and logic controls + 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) @@ -895,7 +896,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*stochastics%t_rp1(i,j) + if (CS%pert_epbl) mech_TKE=mech_TKE*epbl1_wt if (CS%TKE_diagnostics) then eCD%dTKE_conv = 0.0 ; eCD%dTKE_mixing = 0.0 @@ -979,7 +980,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs 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 - mech_TKE = mech_TKE * (1+(exp_kh-1) * stochastics%t_rp2(i,j)) + mech_TKE = mech_TKE * (1+(exp_kh-1) * epbl2_wt) else mech_TKE = mech_TKE * exp_kh endif @@ -2397,9 +2398,9 @@ 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, & + CS%id_epbl1_wts = register_diag_field('ocean_model', 'epbl1_wts', diag%axesT1, Time, & 'random pattern for KE generation', 'None') - CS%id_t_rp2 = register_diag_field('ocean_model', 'random_pattern2', diag%axesT1, Time, & + 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, & From 25ed4fc31f0c9c6dafacfe146f0159ee712edf11 Mon Sep 17 00:00:00 2001 From: pjpegion Date: Fri, 5 Feb 2021 20:18:31 +0000 Subject: [PATCH 18/21] revert MOM_domains.F90 --- src/framework/MOM_domains.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index 33cb45814c..f578df61c9 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -3,9 +3,8 @@ module MOM_domains ! This file is part of MOM6. See LICENSE.md for the license. - use MOM_array_transform, only : rotate_array -use MOM_coms, only : PE_here, root_PE, num_PEs, MOM_infra_init, MOM_infra_end, Get_PElist +use MOM_coms, only : PE_here, root_PE, num_PEs, MOM_infra_init, MOM_infra_end use MOM_coms, only : broadcast, sum_across_PEs, min_across_PEs, max_across_PEs use MOM_cpu_clock, only : cpu_clock_begin, cpu_clock_end use MOM_error_handler, only : MOM_error, MOM_mesg, NOTE, WARNING, FATAL, is_root_pe From 8afe969206a402b596700b09ebc0a342329e5214 Mon Sep 17 00:00:00 2001 From: pjpegion Date: Fri, 5 Feb 2021 20:50:56 +0000 Subject: [PATCH 19/21] clean up of mom_ocean_model_nuopc.F90 --- config_src/nuopc_driver/mom_ocean_model_nuopc.F90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 index aa40f1a1d1..6e57f66a95 100644 --- a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 +++ b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 @@ -51,7 +51,7 @@ module MOM_ocean_model_nuopc use coupler_types_mod, only : coupler_type_spawn, coupler_type_write_chksums use coupler_types_mod, only : coupler_type_initialized, coupler_type_copy_data use coupler_types_mod, only : coupler_type_set_diags, coupler_type_send_data -use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain,mpp_get_pelist +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 fms_mod, only : stdout use mpp_mod, only : mpp_chksum @@ -442,7 +442,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i if (OS%do_sppt .OR. OS%pert_epbl) then num_procs=num_PEs() allocate(pelist(num_procs)) - !call mpp_get_pelist(pelist, commID=mom_comm) call Get_PElist(pelist,commID = mom_comm) me=PE_here() master=root_PE() @@ -698,7 +697,6 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, 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 endif From 689a73ffdc87582a595e496aedfd1d898f862f9b Mon Sep 17 00:00:00 2001 From: pjpegion Date: Fri, 5 Feb 2021 20:54:05 +0000 Subject: [PATCH 20/21] remove PE_here from mom_ocean_model_nuopc.F90 --- config_src/nuopc_driver/mom_ocean_model_nuopc.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 index 6e57f66a95..b24e5d493c 100644 --- a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 +++ b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 @@ -62,7 +62,7 @@ module MOM_ocean_model_nuopc use MOM_surface_forcing_nuopc, only : convert_IOB_to_forces, ice_ocn_bnd_type_chksum use MOM_surface_forcing_nuopc, only : ice_ocean_boundary_type, surface_forcing_CS use MOM_surface_forcing_nuopc, only : forcing_save_restart -use MOM_domains, only : root_PE,PE_here,num_PEs +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 @@ -443,7 +443,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i 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,& From 565e0bb43842f69e6afdf1a88d30da90d643f86b Mon Sep 17 00:00:00 2001 From: pjpegion Date: Fri, 26 Feb 2021 17:43:50 +0000 Subject: [PATCH 21/21] remove debug statements --- src/parameterizations/vertical/MOM_energetic_PBL.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 9f7c465fc9..10ee9bad07 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -985,6 +985,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs 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.