diff --git a/config_src/external/stochastic_physics/stochastic_physics.F90 b/config_src/external/stochastic_physics/stochastic_physics.F90 index 14fa1bf289..fdfd701892 100644 --- a/config_src/external/stochastic_physics/stochastic_physics.F90 +++ b/config_src/external/stochastic_physics/stochastic_physics.F90 @@ -16,7 +16,7 @@ module stochastic_physics !> Initializes the stochastic physics perturbations. subroutine init_stochastic_physics_ocn(delt, geoLonT, geoLatT, nx, ny, nz, pert_epbl_in, do_sppt_in, & - mpiroot, mpicomm, iret) + do_skeb_in,mpiroot, mpicomm, iret) real, intent(in) :: delt !< timestep in seconds between calls to run_stochastic_physics_ocn [s] integer, intent(in) :: nx !< number of gridpoints in the x-direction of the compute grid integer, intent(in) :: ny !< number of gridpoints in the y-direction of the compute grid @@ -25,6 +25,7 @@ subroutine init_stochastic_physics_ocn(delt, geoLonT, geoLatT, nx, ny, nz, pert_ real, intent(in) :: geoLatT(nx,ny) !< Latitude in degrees logical, intent(in) :: pert_epbl_in !< logical flag, if true generate random pattern for ePBL perturbations logical, intent(in) :: do_sppt_in !< logical flag, if true generate random pattern for SPPT perturbations + logical, intent(in) :: do_skeb_in !< logical flag, if true generate random pattern for SKEB perturbations integer, intent(in) :: mpiroot !< root processor integer, intent(in) :: mpicomm !< mpi communicator integer, intent(out) :: iret !< return code @@ -38,14 +39,20 @@ subroutine init_stochastic_physics_ocn(delt, geoLonT, geoLatT, nx, ny, nz, pert_ call MOM_error(WARNING, 'init_stochastic_physics_ocn: do_sppt needs to be false if using the stub') iret=-1 endif + if (do_skeb_in) then + call MOM_error(WARNING, 'init_stochastic_physics_ocn: do_skeb needs to be false if using the stub') + iret=-1 + endif ! This stub function does not actually do anything. return end subroutine init_stochastic_physics_ocn + !> Determines the stochastic physics perturbations. -subroutine run_stochastic_physics_ocn(sppt_wts, t_rp1, t_rp2) +subroutine run_stochastic_physics_ocn(sppt_wts, skeb_wts, t_rp1, t_rp2) real, intent(inout) :: sppt_wts(:,:) !< array containing random weights for SPPT range [0,2] + real, intent(inout) :: skeb_wts(:,:) !< array containing random weights for SKEB real, intent(inout) :: t_rp1(:,:) !< array containing random weights for ePBL !! perturbations (KE generation) range [0,2] real, intent(inout) :: t_rp2(:,:) !< array containing random weights for ePBL diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 77ee1192a2..600439d5b2 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -97,6 +97,9 @@ module MOM_ALE !! values result in the use of more robust and accurate forms of !! mathematically equivalent expressions. + logical :: conserve_ke !< Apply a correction to the baroclinic velocity after remapping to + !! conserve KE. + logical :: debug !< If true, write verbose checksums for debugging purposes. logical :: show_call_tree !< For debugging @@ -117,6 +120,8 @@ module MOM_ALE integer :: id_e_preale = -1 !< diagnostic id for interface heights before ALE. integer :: id_vert_remap_h = -1 !< diagnostic id for layer thicknesses used for remapping integer :: id_vert_remap_h_tendency = -1 !< diagnostic id for layer thickness tendency due to ALE + integer :: id_remap_delta_integ_u2 = -1 !< Change in depth-integrated rho0*u**2/2 + integer :: id_remap_delta_integ_v2 = -1 !< Change in depth-integrated rho0*v**2/2 end type @@ -298,6 +303,11 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) if (CS%use_hybgen_unmix) & call init_hybgen_unmix(CS%hybgen_unmixCS, GV, US, param_file, hybgen_regridCS) + call get_param(param_file, mdl, "REMAP_VEL_CONSERVE_KE", CS%conserve_ke, & + "If true, a correction is applied to the baroclinic component of velocity "//& + "after remapping so that total KE is conserved. KE may not be conserved "//& + "when (CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)", & + default=.false.) call get_param(param_file, "MOM", "DEBUG", CS%debug, & "If true, write out verbose debugging data.", & default=.false., debuggingParam=.true.) @@ -341,13 +351,23 @@ subroutine ALE_register_diags(Time, G, GV, US, diag, CS) CS%id_dzRegrid = register_diag_field('ocean_model', 'dzRegrid', diag%axesTi, Time, & 'Change in interface height due to ALE regridding', 'm', conversion=GV%H_to_m) - cs%id_vert_remap_h = register_diag_field('ocean_model', 'vert_remap_h', diag%axestl, Time, & + CS%id_vert_remap_h = register_diag_field('ocean_model', 'vert_remap_h', diag%axestl, Time, & 'layer thicknesses after ALE regridding and remapping', & thickness_units, conversion=GV%H_to_MKS, v_extensive=.true.) - cs%id_vert_remap_h_tendency = register_diag_field('ocean_model', & + CS%id_vert_remap_h_tendency = register_diag_field('ocean_model', & 'vert_remap_h_tendency', diag%axestl, Time, & 'Layer thicknesses tendency due to ALE regridding and remapping', & trim(thickness_units)//" s-1", conversion=GV%H_to_MKS*US%s_to_T, v_extensive=.true.) + CS%id_remap_delta_integ_u2 = register_diag_field('ocean_model', 'ale_u2', diag%axesCu1, Time, & + 'Rate of change in half rho0 times depth integral of squared zonal'//& + ' velocity by remapping. If REMAP_VEL_CONSERVE_KE is .true. then '//& + ' this measures the change before the KE-conserving correction is applied.', & + 'W m-2', conversion=US%RZ3_T3_to_W_m2 * US%L_to_Z**2) + CS%id_remap_delta_integ_v2 = register_diag_field('ocean_model', 'ale_v2', diag%axesCv1, Time, & + 'Rate of change in half rho0 times depth integral of squared meridional'//& + ' velocity by remapping. If REMAP_VEL_CONSERVE_KE is .true. then '//& + ' this measures the change before the KE-conserving correction is applied.', & + 'W m-2', conversion=US%RZ3_T3_to_W_m2 * US%L_to_Z**2) end subroutine ALE_register_diags @@ -1020,7 +1040,8 @@ end subroutine ALE_remap_set_h_vel_OBC !! This routine may be called during initialization of the model at time=0, to !! remap initial conditions to the model grid. It is also called during a !! time step to update the state. -subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u, v, debug) +subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u, v, debug, & + dt, allow_preserve_variance) type(ALE_CS), intent(in) :: CS !< ALE control structure type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure @@ -1041,6 +1062,9 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1] logical, optional, intent(in) :: debug !< If true, show the call tree + real, optional, intent(in) :: dt !< time step for diagnostics [T ~> s] + logical, optional, intent(in) :: allow_preserve_variance !< If true, enables ke-conserving + !! correction ! Local variables real :: h_mask_vel ! A depth below which the thicknesses at a velocity point are masked out [H ~> m or kg m-2] @@ -1051,6 +1075,16 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u real :: h1(GV%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2] real :: h2(GV%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2] real :: h_neglect, h_neglect_edge ! Tiny thicknesses used in remapping [H ~> m or kg m-2] + real :: rescale_coef ! Factor that scales the baroclinic velocity to conserve ke [nondim] + real :: u_bt, v_bt ! Depth-averaged velocity components [L T-1 ~> m s-1] + real :: ke_c_src, ke_c_tgt ! \int [u_c or v_c]^2 dz on src and tgt grids [H L2 T-2 ~> m3 s-2] + real, dimension(SZIB_(G),SZJ_(G)) :: du2h_tot ! The rate of change of vertically integrated + ! 0.5 * rho0 * u**2 [R Z L2 T-3 ~> W m-2] + real, dimension(SZI_(G),SZJB_(G)) :: dv2h_tot ! The rate of change of vertically integrated + ! 0.5 * rho0 * v**2 [R Z L2 T-3 ~> W m-2] + real :: u2h_tot, v2h_tot ! The vertically integrated u**2 and v**2 [H L2 T-2 ~> m3 s-2 or kg s-2] + real :: I_dt ! 1 / dt [T-1 ~> s-1] + logical :: variance_option ! Contains the value of allow_preserve_variance when present, else false logical :: show_call_tree integer :: i, j, k, nz @@ -1058,6 +1092,17 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u if (present(debug)) show_call_tree = debug if (show_call_tree) call callTree_enter("ALE_remap_velocities()") + ! Setup related to KE conservation + variance_option = .false. + if (present(allow_preserve_variance)) variance_option=allow_preserve_variance + if (present(dt)) I_dt = 1.0 / dt + + if (CS%id_remap_delta_integ_u2>0) du2h_tot(:,:) = 0. + if (CS%id_remap_delta_integ_v2>0) dv2h_tot(:,:) = 0. + + if (((CS%id_remap_delta_integ_u2>0) .or. (CS%id_remap_delta_integ_v2>0)) .and. .not.present(dt))& + call MOM_error(FATAL, "ALE KE diagnostics requires passing dt into ALE_remap_velocities") + if (CS%answer_date >= 20190101) then h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff elseif (GV%Boussinesq) then @@ -1070,7 +1115,9 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u ! --- Remap u profiles from the source vertical grid onto the new target grid. - !$OMP parallel do default(shared) private(h1,h2,u_src,h_mask_vel,u_tgt) + !$OMP parallel do default(shared) private(h1,h2,u_src,h_mask_vel,u_tgt, & + !$OMP u_bt,ke_c_src,ke_c_tgt,rescale_coef, & + !$OMP u2h_tot,v2h_tot) do j=G%jsc,G%jec ; do I=G%IscB,G%IecB ; if (G%mask2dCu(I,j)>0.) then ! Make a 1-d copy of the start and final grids and the source velocity do k=1,nz @@ -1079,9 +1126,47 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u u_src(k) = u(I,j,k) enddo + if (CS%id_remap_delta_integ_u2>0) then + u2h_tot = 0. + do k=1,nz + u2h_tot = u2h_tot - h1(k) * (u_src(k)**2) + enddo + endif + call remapping_core_h(CS%vel_remapCS, nz, h1, u_src, nz, h2, u_tgt, & h_neglect, h_neglect_edge) + if (variance_option .and. CS%conserve_ke) then + ! Conserve ke_u by correcting baroclinic component. + ! Assumes total depth doesn't change during remap, and + ! that \int u(z) dz doesn't change during remap. + ! First get barotropic component + u_bt = 0.0 + do k=1,nz + u_bt = u_bt + h2(k) * u_tgt(k) ! Dimensions [H L T-1] + enddo + u_bt = u_bt / (sum(h2(1:nz)) + h_neglect) ! Dimensions return to [L T-1] + ! Next get baroclinic ke = \int (u-u_bt)^2 from source and target + ke_c_src = 0.0 + ke_c_tgt = 0.0 + do k=1,nz + ke_c_src = ke_c_src + h1(k) * (u_src(k) - u_bt)**2 + ke_c_tgt = ke_c_tgt + h2(k) * (u_tgt(k) - u_bt)**2 + enddo + ! Next rescale baroclinic component on target grid to conserve ke + rescale_coef = min(1.25, sqrt(ke_c_src / (ke_c_tgt + 1.E-19))) + do k=1,nz + u_tgt(k) = u_bt + rescale_coef * (u_tgt(k) - u_bt) + enddo + endif + + if (CS%id_remap_delta_integ_u2>0) then + do k=1,nz + u2h_tot = u2h_tot + h2(k) * (u_tgt(k)**2) + enddo + du2h_tot(I,j) = GV%H_to_RZ * u2h_tot * I_dt + endif + if ((CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)) & call mask_near_bottom_vel(u_tgt, h2, CS%BBL_h_vel_mask, CS%h_vel_mask, nz) @@ -1091,12 +1176,16 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u enddo !k endif ; enddo ; enddo + if (CS%id_remap_delta_integ_u2>0) call post_data(CS%id_remap_delta_integ_u2, du2h_tot, CS%diag) + if (show_call_tree) call callTree_waypoint("u remapped (ALE_remap_velocities)") ! --- Remap v profiles from the source vertical grid onto the new target grid. - !$OMP parallel do default(shared) private(h1,h2,v_src,h_mask_vel,v_tgt) + !$OMP parallel do default(shared) private(h1,h2,v_src,h_mask_vel,v_tgt, & + !$OMP v_bt,ke_c_src,ke_c_tgt,rescale_coef, & + !$OMP u2h_tot,v2h_tot) do J=G%JscB,G%JecB ; do i=G%isc,G%iec ; if (G%mask2dCv(i,J)>0.) then do k=1,nz @@ -1105,9 +1194,47 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u v_src(k) = v(i,J,k) enddo + if (CS%id_remap_delta_integ_v2>0) then + v2h_tot = 0. + do k=1,nz + v2h_tot = v2h_tot - h1(k) * (v_src(k)**2) + enddo + endif + call remapping_core_h(CS%vel_remapCS, nz, h1, v_src, nz, h2, v_tgt, & h_neglect, h_neglect_edge) + if (variance_option .and. CS%conserve_ke) then + ! Conserve ke_v by correcting baroclinic component. + ! Assumes total depth doesn't change during remap, and + ! that \int v(z) dz doesn't change during remap. + ! First get barotropic component + v_bt = 0.0 + do k=1,nz + v_bt = v_bt + h2(k) * v_tgt(k) ! Dimensions [H L T-1] + enddo + v_bt = v_bt / (sum(h2(1:nz)) + h_neglect) ! Dimensions return to [L T-1] + ! Next get baroclinic ke = \int (u-u_bt)^2 from source and target + ke_c_src = 0.0 + ke_c_tgt = 0.0 + do k=1,nz + ke_c_src = ke_c_src + h1(k) * (v_src(k) - v_bt)**2 + ke_c_tgt = ke_c_tgt + h2(k) * (v_tgt(k) - v_bt)**2 + enddo + ! Next rescale baroclinic component on target grid to conserve ke + rescale_coef = min(1.25, sqrt(ke_c_src / (ke_c_tgt + 1.E-19))) + do k=1,nz + v_tgt(k) = v_bt + rescale_coef * (v_tgt(k) - v_bt) + enddo + endif + + if (CS%id_remap_delta_integ_v2>0) then + do k=1,nz + v2h_tot = v2h_tot + h2(k) * (v_tgt(k)**2) + enddo + dv2h_tot(I,j) = GV%H_to_RZ * v2h_tot * I_dt + endif + if ((CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)) then call mask_near_bottom_vel(v_tgt, h2, CS%BBL_h_vel_mask, CS%h_vel_mask, nz) endif @@ -1118,6 +1245,8 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u enddo !k endif ; enddo ; enddo + if (CS%id_remap_delta_integ_v2>0) call post_data(CS%id_remap_delta_integ_v2, dv2h_tot, CS%diag) + if (show_call_tree) call callTree_waypoint("v remapped (ALE_remap_velocities)") if (show_call_tree) call callTree_leave("ALE_remap_velocities()") diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 52941944a9..7b080a5537 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -64,7 +64,7 @@ module MOM use MOM_diabatic_driver, only : diabatic, diabatic_driver_init, diabatic_CS, extract_diabatic_member use MOM_diabatic_driver, only : adiabatic, adiabatic_driver_init, diabatic_driver_end use MOM_diabatic_driver, only : register_diabatic_restarts -use MOM_stochastics, only : stochastics_init, update_stochastics, stochastic_CS +use MOM_stochastics, only : stochastics_init, update_stochastics, stochastic_CS, apply_skeb use MOM_diagnostics, only : calculate_diagnostic_fields, MOM_diagnostics_init use MOM_diagnostics, only : register_transport_diags, post_transport_diagnostics use MOM_diagnostics, only : register_surface_diags, write_static_fields @@ -741,7 +741,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS endif endif ! advance the random pattern if stochastic physics is active - if (CS%stoch_CS%do_sppt .OR. CS%stoch_CS%pert_epbl) call update_stochastics(CS%stoch_CS) + if (CS%stoch_CS%do_sppt .OR. CS%stoch_CS%pert_epbl .OR. CS%stoch_CS%do_skeb) & + call update_stochastics(CS%stoch_CS) if (do_dyn) then if (nonblocking_p_surf_update) & @@ -1151,7 +1152,8 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & if (CS%VarMix%use_variable_mixing) & call calc_slope_functions(h, CS%tv, dt, G, GV, US, CS%VarMix, OBC=CS%OBC) call thickness_diffuse(h, CS%uhtr, CS%vhtr, CS%tv, dt_thermo, G, GV, US, & - CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp) + CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp, & + CS%stoch_CS) call cpu_clock_end(id_clock_thick_diff) call pass_var(h, G%Domain, clock=id_clock_pass, halo=max(2,CS%cont_stencil)) if (showCallTree) call callTree_waypoint("finished thickness_diffuse_first (step_MOM)") @@ -1223,7 +1225,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & call step_MOM_dyn_split_RK2(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, & CS%eta_av_bc, G, GV, US, CS%dyn_split_RK2_CSp, calc_dtbt, CS%VarMix, & - CS%MEKE, CS%thickness_diffuse_CSp, CS%pbv, waves=waves) + CS%MEKE, CS%thickness_diffuse_CSp, CS%pbv, CS%stoch_CS, waves=waves) if (showCallTree) call callTree_waypoint("finished step_MOM_dyn_split (step_MOM)") elseif (CS%do_dynamics) then ! ------------------------------------ not SPLIT @@ -1237,11 +1239,13 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & if (CS%use_RK2) then call step_MOM_dyn_unsplit_RK2(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, & - CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_RK2_CSp, CS%VarMix, CS%MEKE, CS%pbv) + CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_RK2_CSp, CS%VarMix, CS%MEKE, CS%pbv, & + CS%stoch_CS) else call step_MOM_dyn_unsplit(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, & - CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_CSp, CS%VarMix, CS%MEKE, CS%pbv, Waves=Waves) + CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_CSp, CS%VarMix, CS%MEKE, CS%pbv, & + CS%stoch_CS, Waves=Waves) endif if (showCallTree) call callTree_waypoint("finished step_MOM_dyn_unsplit (step_MOM)") @@ -1282,7 +1286,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & if (CS%VarMix%use_variable_mixing) & call calc_slope_functions(h, CS%tv, dt, G, GV, US, CS%VarMix, OBC=CS%OBC) call thickness_diffuse(h, CS%uhtr, CS%vhtr, CS%tv, dt, G, GV, US, & - CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp) + CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp, CS%stoch_CS) if (CS%debug) call hchksum(h,"Post-thickness_diffuse h", G%HI, haloshift=1, scale=GV%H_to_MKS) call cpu_clock_end(id_clock_thick_diff) @@ -1584,6 +1588,10 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & Time_end_thermo, G, GV, US, CS%diabatic_CSp, CS%stoch_CS, CS%OBC, Waves) fluxes%fluxes_used = .true. + if (CS%stoch_CS%do_skeb) then + call apply_skeb(CS%G,CS%GV,CS%stoch_CS,CS%u,CS%v,CS%h,CS%tv,dtdia,Time_end_thermo) + endif + if (showCallTree) call callTree_waypoint("finished diabatic (step_MOM_thermo)") ! Regridding/remapping is done here, at end of thermodynamics time step @@ -1641,7 +1649,8 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & endif ! Remap the velocity components. - call ALE_remap_velocities(CS%ALE_CSp, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u, v, showCallTree) + call ALE_remap_velocities(CS%ALE_CSp, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u, v, showCallTree, & + dtdia, allow_preserve_variance=.true.) if (allocated(tv%SpV_avg)) tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid. diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index debc63cb46..14b0942009 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -48,7 +48,7 @@ module MOM_dynamics_split_RK2 use MOM_debugging, only : check_redundant use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type -use MOM_hor_visc, only : horizontal_viscosity, hor_visc_CS +use MOM_hor_visc, only : horizontal_viscosity, hor_visc_CS, hor_visc_vel_stencil use MOM_hor_visc, only : hor_visc_init, hor_visc_end use MOM_interface_heights, only : thickness_to_dz, find_col_avg_SpV use MOM_lateral_mixing_coeffs, only : VarMix_CS @@ -59,6 +59,7 @@ module MOM_dynamics_split_RK2 use MOM_PressureForce, only : PressureForce, PressureForce_CS use MOM_PressureForce, only : PressureForce_init use MOM_set_visc, only : set_viscous_ML, set_visc_CS +use MOM_stochastics, only : stochastic_CS use MOM_thickness_diffuse, only : thickness_diffuse_CS use MOM_self_attr_load, only : SAL_CS use MOM_self_attr_load, only : SAL_init, SAL_end @@ -286,7 +287,7 @@ module MOM_dynamics_split_RK2 !> RK2 splitting for time stepping MOM adiabatic dynamics subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_surf_begin, p_surf_end, & uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, calc_dtbt, VarMix, & - MEKE, thickness_diffuse_CSp, pbv, Waves) + MEKE, thickness_diffuse_CSp, pbv, STOCH, 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 @@ -326,6 +327,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s type(thickness_diffuse_CS), intent(inout) :: thickness_diffuse_CSp !< Pointer to a structure containing !! interface height diffusivities type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics + type(stochastic_CS), intent(inout) :: STOCH !< Stochastic control structure type(wave_parameters_CS), optional, pointer :: Waves !< A pointer to a structure containing !! fields related to the surface wave conditions @@ -399,7 +401,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s logical :: showCallTree, sym integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz - integer :: cont_stencil, obc_stencil + integer :: cont_stencil, obc_stencil, vel_stencil is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -466,19 +468,20 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s if (associated(CS%OBC)) then if (CS%OBC%oblique_BCs_exist_globally) obc_stencil = 3 endif + vel_stencil = max(2, obc_stencil, hor_visc_vel_stencil(CS%hor_visc)) call cpu_clock_begin(id_clock_pass) call create_group_pass(CS%pass_eta, eta, G%Domain, halo=1) call create_group_pass(CS%pass_visc_rem, CS%visc_rem_u, CS%visc_rem_v, G%Domain, & To_All+SCALAR_PAIR, CGRID_NE, halo=max(1,cont_stencil)) call create_group_pass(CS%pass_uvp, up, vp, G%Domain, halo=max(1,cont_stencil)) call create_group_pass(CS%pass_hp_uv, hp, G%Domain, halo=2) - call create_group_pass(CS%pass_hp_uv, u_av, v_av, G%Domain, halo=max(2,obc_stencil)) - call create_group_pass(CS%pass_hp_uv, uh(:,:,:), vh(:,:,:), G%Domain, halo=max(2,obc_stencil)) + call create_group_pass(CS%pass_hp_uv, u_av, v_av, G%Domain, halo=vel_stencil) + call create_group_pass(CS%pass_hp_uv, uh(:,:,:), vh(:,:,:), G%Domain, halo=vel_stencil) call create_group_pass(CS%pass_uv, u, v, G%Domain, halo=max(2,cont_stencil)) call create_group_pass(CS%pass_h, h, G%Domain, halo=max(2,cont_stencil)) - call create_group_pass(CS%pass_av_uvh, u_av, v_av, G%Domain, halo=max(2,obc_stencil)) - call create_group_pass(CS%pass_av_uvh, uh(:,:,:), vh(:,:,:), G%Domain, halo=max(2,obc_stencil)) + call create_group_pass(CS%pass_av_uvh, u_av, v_av, G%Domain, halo=vel_stencil) + call create_group_pass(CS%pass_av_uvh, uh(:,:,:), vh(:,:,:), G%Domain, halo=vel_stencil) call cpu_clock_end(id_clock_pass) !--- end set up for group halo pass @@ -839,7 +842,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s if (CS%debug) then call MOM_state_chksum("Predictor ", up, vp, hp, uh, vh, G, GV, US, symmetric=sym) - call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s) + call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=vel_stencil, symmetric=sym, scale=US%L_T_to_m_s) call hchksum(h_av, "Predictor avg h", G%HI, haloshift=0, scale=GV%H_to_MKS) ! call MOM_state_chksum("Predictor avg ", u_av, v_av, h_av, uh, vh, G, GV, US) call check_redundant("Predictor up ", up, vp, G, unscale=US%L_T_to_m_s) @@ -851,7 +854,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, & MEKE, Varmix, G, GV, US, CS%hor_visc, & OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, & - ADp=CS%ADp) + ADp=CS%ADp, STOCH=STOCH) call cpu_clock_end(id_clock_horvisc) if (showCallTree) call callTree_wayPoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2)") diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index c87e6e9958..f1d3311a89 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -87,8 +87,9 @@ module MOM_dynamics_unsplit use MOM_open_boundary, only : open_boundary_zero_normal_flow use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS use MOM_set_visc, only : set_viscous_ML, set_visc_CS -use MOM_self_attr_load, only : SAL_init, SAL_end, SAL_CS +use MOM_stochastics, only : stochastic_CS use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_end, tidal_forcing_CS +use MOM_self_attr_load, only : SAL_init, SAL_end, SAL_CS use MOM_unit_scaling, only : unit_scale_type use MOM_vert_friction, only : vertvisc, vertvisc_coef, vertvisc_init, vertvisc_CS use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units @@ -189,7 +190,7 @@ module MOM_dynamics_unsplit !! 3rd order (for the inviscid momentum equations) order scheme subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, & - VarMix, MEKE, pbv, Waves) + VarMix, MEKE, pbv, STOCH, 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 @@ -223,6 +224,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & type(VarMix_CS), intent(inout) :: VarMix !< Variable mixing control structure type(MEKE_type), intent(inout) :: MEKE !< MEKE fields type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics + type(stochastic_CS), intent(inout) :: STOCH !< Stochastic control structure type(wave_parameters_CS), optional, pointer :: Waves !< A pointer to a structure containing !! fields related to the surface wave conditions @@ -263,7 +265,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! diffu = horizontal viscosity terms (u,h) call enable_averages(dt, Time_local, CS%diag) call cpu_clock_begin(id_clock_horvisc) - call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc) + call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc, STOCH=STOCH) call cpu_clock_end(id_clock_horvisc) call disable_averaging(CS%diag) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index b515229566..3c71fac0e4 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -87,6 +87,7 @@ module MOM_dynamics_unsplit_RK2 use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS use MOM_set_visc, only : set_viscous_ML, set_visc_CS use MOM_self_attr_load, only : SAL_init, SAL_end, SAL_CS +use MOM_stochastics, only : stochastic_CS use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_end, tidal_forcing_CS use MOM_unit_scaling, only : unit_scale_type use MOM_vert_friction, only : vertvisc, vertvisc_coef, vertvisc_init, vertvisc_CS @@ -192,7 +193,7 @@ module MOM_dynamics_unsplit_RK2 !> Step the MOM6 dynamics using an unsplit quasi-2nd order Runge-Kutta scheme subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, & - VarMix, MEKE, pbv) + VarMix, MEKE, pbv, STOCH) 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 @@ -237,6 +238,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, !! fields related to the Mesoscale !! Eddy Kinetic Energy. type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics + type(stochastic_CS), intent(inout) :: STOCH !< Stochastic control structure ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_av ! Averaged layer thicknesses [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: hp ! Predicted layer thicknesses [H ~> m or kg m-2] @@ -276,7 +278,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, call enable_averages(dt,Time_local, CS%diag) call cpu_clock_begin(id_clock_horvisc) call horizontal_viscosity(u_in, v_in, h_in, CS%diffu, CS%diffv, MEKE, VarMix, & - G, GV, US, CS%hor_visc) + G, GV, US, CS%hor_visc, STOCH=STOCH) call cpu_clock_end(id_clock_horvisc) call disable_averaging(CS%diag) call pass_vector(CS%diffu, CS%diffv, G%Domain, clock=id_clock_pass) diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index 46d0448699..40e7661e7f 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -11,7 +11,7 @@ module MOM_shared_initialization use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, log_param, param_file_type, log_version -use MOM_io, only : create_MOM_file, file_exists, field_size +use MOM_io, only : create_MOM_file, file_exists, field_size, get_filename_appendix use MOM_io, only : MOM_infra_file, MOM_field use MOM_io, only : MOM_read_data, MOM_read_vector, read_variable, stdout use MOM_io, only : open_file_to_read, close_file_to_read, SINGLE_FILE, MULTIPLE @@ -1348,6 +1348,7 @@ subroutine write_ocean_geometry_file(G, param_file, directory, US, geom_file) ! Local variables. character(len=240) :: filepath ! The full path to the file to write character(len=40) :: mdl = "write_ocean_geometry_file" + character(len=32) :: filename_appendix = '' ! Appendix to geom filename for ensemble runs type(vardesc), dimension(:), allocatable :: & vars ! Types with metadata about the variables and their staggering type(MOM_field), dimension(:), allocatable :: & @@ -1355,6 +1356,7 @@ subroutine write_ocean_geometry_file(G, param_file, directory, US, geom_file) type(MOM_infra_file) :: IO_handle ! The I/O handle of the fileset integer :: nFlds ! The number of variables in this file integer :: file_threading + integer :: geom_file_len ! geometry file name length logical :: multiple_files call callTree_enter('write_ocean_geometry_file()') @@ -1408,6 +1410,17 @@ subroutine write_ocean_geometry_file(G, param_file, directory, US, geom_file) filepath = trim(directory) // "ocean_geometry" endif + ! Append ensemble run number to filename if it is an ensemble run + call get_filename_appendix(filename_appendix) + if (len_trim(filename_appendix) > 0) then + geom_file_len = len_trim(filepath) + if (filepath(geom_file_len-2:geom_file_len) == ".nc") then + filepath = filepath(1:geom_file_len-3) // '.' // trim(filename_appendix) // ".nc" + else + filepath = filepath // '.' // trim(filename_appendix) + endif + endif + call get_param(param_file, mdl, "PARALLEL_RESTARTFILES", multiple_files, & "If true, the IO layout is used to group processors that write to the same "//& "restart file or each processor writes its own (numbered) restart file. "//& diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 9b1d81348e..6d188990a1 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -20,6 +20,7 @@ module MOM_hor_visc use MOM_MEKE_types, only : MEKE_type use MOM_open_boundary, only : ocean_OBC_type, OBC_DIRECTION_E, OBC_DIRECTION_W use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S, OBC_NONE +use MOM_stochastics, only : stochastic_CS use MOM_unit_scaling, only : unit_scale_type use MOM_verticalGrid, only : verticalGrid_type use MOM_variables, only : accel_diag_ptrs @@ -30,7 +31,7 @@ module MOM_hor_visc #include -public horizontal_viscosity, hor_visc_init, hor_visc_end +public horizontal_viscosity, hor_visc_init, hor_visc_end, hor_visc_vel_stencil !> Control structure for horizontal viscosity type, public :: hor_visc_CS ; private @@ -241,7 +242,7 @@ module MOM_hor_visc !! v[is-2:ie+2,js-2:je+2] !! h[is-1:ie+1,js-1:je+1] subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, & - CS, OBC, BT, TD, ADp) + CS, OBC, BT, TD, ADp, STOCH) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & @@ -265,6 +266,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, type(barotropic_CS), intent(in), optional :: BT !< Barotropic control structure type(thickness_diffuse_CS), intent(in), optional :: TD !< Thickness diffusion control structure type(accel_diag_ptrs), intent(in), optional :: ADp !< Acceleration diagnostics + type(stochastic_CS), intent(inout), optional :: STOCH !< Stochastic control structure ! Local variables real, dimension(SZIB_(G),SZJ_(G)) :: & @@ -395,6 +397,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, logical :: apply_OBC = .false. logical :: use_MEKE_Ku logical :: use_MEKE_Au + logical :: skeb_use_frict integer :: is_vort, ie_vort, js_vort, je_vort ! Loop ranges for vorticity terms integer :: is_Kh, ie_Kh, js_Kh, je_Kh ! Loop ranges for thickness point viscosities integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz @@ -426,6 +429,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, inv_PI2 = 1.0/((4.0*atan(1.0))**2) inv_PI6 = inv_PI3 * inv_PI3 + skeb_use_frict = .false. + if (present(STOCH)) skeb_use_frict = STOCH%skeb_use_frict + m_leithy(:,:) = 0.0 ! Initialize if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then @@ -588,12 +594,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP CS, G, GV, US, OBC, VarMix, MEKE, u, v, h, & !$OMP is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, is_vort, ie_vort, js_vort, je_vort, & !$OMP is_Kh, ie_Kh, js_Kh, je_Kh, apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, & - !$OMP use_MEKE_Ku, use_MEKE_Au, u_smooth, v_smooth, & + !$OMP use_MEKE_Ku, use_MEKE_Au, u_smooth, v_smooth, skeb_use_frict, & !$OMP backscat_subround, GME_effic_h, GME_effic_q, & !$OMP h_neglect, h_neglect3, inv_PI3, inv_PI6, & !$OMP diffu, diffv, Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_GME, & !$OMP div_xx_h, sh_xx_h, vort_xy_q, sh_xy_q, GME_coeff_h, GME_coeff_q, & - !$OMP KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah, NoSt, ShSt & + !$OMP KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah, NoSt, ShSt, STOCH & !$OMP ) & !$OMP private( & !$OMP i, j, k, n, & @@ -1192,10 +1198,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if ((CS%Smagorinsky_Ah) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) then if (CS%Smagorinsky_Ah) then if (CS%bound_Coriolis) then - do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh + do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh AhSm = Shear_mag(i,j) * (CS%Biharm_const_xx(i,j) & - + CS%Biharm_const2_xx(i,j) * Shear_mag(i,j) & - ) + + CS%Biharm_const2_xx(i,j) * Shear_mag(i,j)) Ah(i,j) = max(Ah(i,j), AhSm) enddo ; enddo else @@ -1426,10 +1431,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! Pass the velocity gradients and thickness to ZB2020 if (CS%use_ZB2020) then - call ZB2020_copy_gradient_and_thickness( & - sh_xx, sh_xy, vort_xy, & - hq, & - G, GV, CS%ZB2020, k) + call ZB2020_copy_gradient_and_thickness(sh_xx, sh_xy, vort_xy, hq, G, GV, CS%ZB2020, k) endif if (CS%Laplacian) then @@ -1569,8 +1571,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%bound_Coriolis) then do J=js-1,Jeq ; do I=is-1,Ieq AhSm = Shear_mag(I,J) * (CS%Biharm_const_xy(I,J) & - + CS%Biharm_const2_xy(I,J) * Shear_mag(I,J) & - ) + + CS%Biharm_const2_xy(I,J) * Shear_mag(I,J)) Ah(I,J) = max(Ah(I,J), AhSm) enddo ; enddo else @@ -1599,8 +1600,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! *Add* the MEKE contribution do J=js-1,Jeq ; do I=is-1,Ieq Ah(I,J) = Ah(I,J) + 0.25 * ( & - (MEKE%Au(i,j) + MEKE%Au(i+1,j+1)) + (MEKE%Au(i+1,j) + MEKE%Au(i,j+1)) & - ) + (MEKE%Au(i,j) + MEKE%Au(i+1,j+1)) + (MEKE%Au(i+1,j) + MEKE%Au(i,j+1)) ) enddo ; enddo endif @@ -1786,6 +1786,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, + (v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)) ) ) ) enddo ; enddo ; endif + if (skeb_use_frict) then ; do j=js,je ; do i=is,ie + ! Note that the sign convention is FrictWork < 0 means energy dissipation. + STOCH%skeb_diss(i,j,k) = STOCH%skeb_diss(i,j,k) - STOCH%skeb_frict_coef * & + FrictWork(i,j,k) / (GV%H_to_RZ * (h(i,j,k) + h_neglect)) + enddo ; enddo ; endif + ! Make a similar calculation as for FrictWork above but accumulating into ! the vertically integrated MEKE source term, and adjusting for any ! energy loss seen as a reduction in the (biharmonic) frictional source term. @@ -1885,11 +1891,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%debug) then if (CS%Laplacian) then + ! In symmetric memory mode, Kh_h should also be valid with a haloshift of 1. call hchksum(Kh_h, "Kh_h", G%HI, haloshift=0, scale=US%L_to_m**2*US%s_to_T) - call Bchksum(Kh_q, "Kh_q", G%HI, haloshift=0, scale=US%L_to_m**2*US%s_to_T) + call Bchksum(Kh_q, "Kh_q", G%HI, haloshift=0, symmetric=.true., scale=US%L_to_m**2*US%s_to_T) + endif + if (CS%biharmonic) then + ! In symmetric memory mode, Ah_h should also be valid with a haloshift of 1. + call hchksum(Ah_h, "Ah_h", G%HI, haloshift=0, scale=US%L_to_m**4*US%s_to_T) + call Bchksum(Ah_q, "Ah_q", G%HI, haloshift=0, symmetric=.true., scale=US%L_to_m**4*US%s_to_T) endif - if (CS%biharmonic) call hchksum(Ah_h, "Ah_h", G%HI, haloshift=0, scale=US%L_to_m**4*US%s_to_T) - if (CS%biharmonic) call Bchksum(Ah_q, "Ah_q", G%HI, haloshift=0, scale=US%L_to_m**4*US%s_to_T) endif if (CS%id_FrictWorkIntz > 0) then @@ -2391,14 +2401,31 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) ALLOC_(CS%m_leithy_max(isd:ied,jsd:jed)) ; CS%m_leithy_max(:,:) = 0.0 endif if (CS%Re_Ah > 0.0) then - ALLOC_(CS%Re_Ah_const_xx(isd:ied,jsd:jed)); CS%Re_Ah_const_xx(:,:) = 0.0 - ALLOC_(CS%Re_Ah_const_xy(IsdB:IedB,JsdB:JedB)); CS%Re_Ah_const_xy(:,:) = 0.0 + ALLOC_(CS%Re_Ah_const_xx(isd:ied,jsd:jed)) ; CS%Re_Ah_const_xx(:,:) = 0.0 + ALLOC_(CS%Re_Ah_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Re_Ah_const_xy(:,:) = 0.0 endif endif do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 CS%dx2q(I,J) = G%dxBu(I,J)*G%dxBu(I,J) ; CS%dy2q(I,J) = G%dyBu(I,J)*G%dyBu(I,J) - CS%DX_dyBu(I,J) = G%dxBu(I,J)*G%IdyBu(I,J) ; CS%DY_dxBu(I,J) = G%dyBu(I,J)*G%IdxBu(I,J) enddo ; enddo + + if (((CS%Leith_Kh) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) .and. & + ((G%isc-G%isd < 3) .or. (G%isc-G%isd < 3))) call MOM_error(FATAL, & + "The minimum halo size is 3 when a Leith viscosity is being used.") + if (CS%use_Leithy) then + do J=js-3,Jeq+2 ; do I=is-3,Ieq+2 + CS%DX_dyBu(I,J) = G%dxBu(I,J)*G%IdyBu(I,J) ; CS%DY_dxBu(I,J) = G%dyBu(I,J)*G%IdxBu(I,J) + enddo ; enddo + elseif ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then + do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2 + CS%DX_dyBu(I,J) = G%dxBu(I,J)*G%IdyBu(I,J) ; CS%DY_dxBu(I,J) = G%dyBu(I,J)*G%IdxBu(I,J) + enddo ; enddo + else + do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 + CS%DX_dyBu(I,J) = G%dxBu(I,J)*G%IdyBu(I,J) ; CS%DY_dxBu(I,J) = G%dyBu(I,J)*G%IdxBu(I,J) + enddo ; enddo + endif + do j=js-2,Jeq+2 ; do i=is-2,Ieq+2 CS%dx2h(i,j) = G%dxT(i,j)*G%dxT(i,j) ; CS%dy2h(i,j) = G%dyT(i,j)*G%dyT(i,j) CS%DX_dyT(i,j) = G%dxT(i,j)*G%IdyT(i,j) ; CS%DY_dxT(i,j) = G%dyT(i,j)*G%IdxT(i,j) @@ -2529,12 +2556,12 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) endif endif if (CS%Leith_Ah) then - CS%biharm6_const_xx(i,j) = Leith_bi_const * (grid_sp_h3 * grid_sp_h3) + CS%biharm6_const_xx(i,j) = Leith_bi_const * (grid_sp_h3 * grid_sp_h3) endif if (CS%use_Leithy) then - CS%biharm6_const_xx(i,j) = Leith_bi_const * max(G%dxT(i,j),G%dyT(i,j))**6 - CS%m_const_leithy(i,j) = 0.5 * sqrt(CS%c_K) * max(G%dxT(i,j),G%dyT(i,j)) - CS%m_leithy_max(i,j) = 4. / max(G%dxT(i,j),G%dyT(i,j))**2 + CS%biharm6_const_xx(i,j) = Leith_bi_const * max(G%dxT(i,j),G%dyT(i,j))**6 + CS%m_const_leithy(i,j) = 0.5 * sqrt(CS%c_K) * max(G%dxT(i,j),G%dyT(i,j)) + CS%m_leithy_max(i,j) = 4. / max(G%dxT(i,j),G%dyT(i,j))**2 endif CS%Ah_bg_xx(i,j) = MAX(Ah, Ah_vel_scale * grid_sp_h2 * sqrt(grid_sp_h2)) if (CS%Re_Ah > 0.0) CS%Re_Ah_const_xx(i,j) = grid_sp_h3 / CS%Re_Ah @@ -2559,12 +2586,12 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) endif endif if ((CS%Leith_Ah) .or. (CS%use_Leithy))then - CS%biharm6_const_xy(I,J) = Leith_bi_const * (grid_sp_q3 * grid_sp_q3) + CS%biharm6_const_xy(I,J) = Leith_bi_const * (grid_sp_q3 * grid_sp_q3) endif CS%Ah_bg_xy(I,J) = MAX(Ah, Ah_vel_scale * grid_sp_q2 * sqrt(grid_sp_q2)) if (CS%Re_Ah > 0.0) CS%Re_Ah_const_xy(i,j) = grid_sp_q3 / CS%Re_Ah if (Ah_time_scale > 0.) CS%Ah_bg_xy(i,j) = & - MAX(CS%Ah_bg_xy(i,j), (grid_sp_q2 * grid_sp_q2) / Ah_time_scale) + MAX(CS%Ah_bg_xy(i,j), (grid_sp_q2 * grid_sp_q2) / Ah_time_scale) if (CS%bound_Ah .and. .not.CS%better_bound_Ah) then CS%Ah_Max_xy(I,J) = Ah_Limit * (grid_sp_q2 * grid_sp_q2) CS%Ah_bg_xy(I,J) = MIN(CS%Ah_bg_xy(I,J), CS%Ah_Max_xy(I,J)) @@ -2810,6 +2837,18 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) end subroutine hor_visc_init +!> hor_visc_vel_stencil returns the horizontal viscosity input velocity stencil size +function hor_visc_vel_stencil(CS) result(stencil) + type(hor_visc_CS), intent(in) :: CS !< Control structure for horizontal viscosity + integer :: stencil !< The horizontal viscosity velocity stencil size with the current settings. + + stencil = 2 + + if ((CS%Leith_Kh) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) then + stencil = 3 + endif +end function hor_visc_vel_stencil + !> Calculates factors in the anisotropic orientation tensor to be align with the grid. !! With n1=1 and n2=0, this recovers the approach of Large et al, 2001. subroutine align_aniso_tensor_to_grid(CS, n1, n2) diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 2638ca71e1..52b55ad252 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -18,6 +18,7 @@ module MOM_thickness_diffuse use MOM_isopycnal_slopes, only : vert_fill_TS use MOM_lateral_mixing_coeffs, only : VarMix_CS use MOM_MEKE_types, only : MEKE_type +use MOM_stochastics, only : stochastic_CS use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs, cont_diag_ptrs use MOM_verticalGrid, only : verticalGrid_type @@ -119,7 +120,7 @@ module MOM_thickness_diffuse !> Calculates isopycnal height diffusion coefficients and applies isopycnal height diffusion !! by modifying to the layer thicknesses, h. Diffusivities are limited to ensure stability. !! Also returns along-layer mass fluxes used in the continuity equation. -subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp, CS) +subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp, CS, STOCH) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -134,6 +135,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp type(VarMix_CS), target, intent(in) :: VarMix !< Variable mixing coefficients type(cont_diag_ptrs), intent(inout) :: CDp !< Diagnostics for the continuity equation type(thickness_diffuse_CS), intent(inout) :: CS !< Control structure for thickness_diffuse + type(stochastic_CS), intent(inout) :: STOCH !< Stochastic control structure ! Local variables real :: e(SZI_(G),SZJ_(G),SZK_(GV)+1) ! heights of interfaces, relative to mean ! sea level [Z ~> m], positive up. @@ -477,12 +479,23 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp endif ! Calculate uhD, vhD from h, e, KH_u, KH_v, tv%T/S - if (use_stored_slopes) then - call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, CS, & - int_slope_u, int_slope_v, VarMix%slope_x, VarMix%slope_y) + if (STOCH%skeb_use_gm) then + if (use_stored_slopes) then + call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, CS, & + int_slope_u, int_slope_v, VarMix%slope_x, VarMix%slope_y, & + STOCH=STOCH, VarMix=VarMix) + else + call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, CS, & + int_slope_u, int_slope_v, STOCH=STOCH, VarMix=VarMix) + endif else - call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, CS, & - int_slope_u, int_slope_v) + if (use_stored_slopes) then + call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, CS, & + int_slope_u, int_slope_v, VarMix%slope_x, VarMix%slope_y) + else + call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, CS, & + int_slope_u, int_slope_v) + endif endif if (VarMix%use_variable_mixing) then @@ -593,7 +606,7 @@ end subroutine thickness_diffuse !! Fluxes are limited to give positive definite thicknesses. !! Called by thickness_diffuse(). subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, & - CS, int_slope_u, int_slope_v, slope_x, slope_y) + CS, int_slope_u, int_slope_v, slope_x, slope_y, STOCH, VarMix) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -622,6 +635,8 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV !! density gradients [nondim]. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), optional, intent(in) :: slope_x !< Isopyc. slope at u [Z L-1 ~> nondim] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), optional, intent(in) :: slope_y !< Isopyc. slope at v [Z L-1 ~> nondim] + type(stochastic_CS), optional, intent(inout) :: STOCH !< Stochastic control structure + type(VarMix_CS), target, optional, intent(in) :: VarMix !< Variable mixing coefficents ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & @@ -759,6 +774,11 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! [H L2 T-1 ~> m3 s-1 or kg s-1] real :: diag_sfn_unlim_y(SZI_(G),SZJB_(G),SZK_(GV)+1) ! Diagnostic of the y-face streamfunction before ! applying limiters [Z L2 T-1 ~> m3 s-1] + ! applying limiters [H L2 T-1 ~> m3 s-1 or kg s-1] + real, allocatable :: skeb_gm_work(:,:) ! Temp array to hold GM work for SKEB + real, allocatable :: skeb_ebt_norm2(:,:) ! Used to normalize EBT for SKEB + real :: h_tot ! total depth [H ~> m] + logical :: present_slope_x, present_slope_y, calc_derivatives integer, dimension(2) :: EOSdom_u ! The shifted I-computational domain to use for equation of ! state calculations at u-points. @@ -766,7 +786,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! state calculations at v-points. integer, dimension(2) :: EOSdom_h1 ! The shifted i-computational domain to use for equation of ! state calculations at h points with 1 extra halo point - logical :: use_stanley + logical :: use_stanley, skeb_use_gm integer :: is, ie, js, je, nz, IsdB, halo integer :: i, j, k is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke ; IsdB = G%IsdB @@ -786,6 +806,13 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV use_stanley = CS%use_stanley_gm + skeb_use_gm = .false. + if (present(STOCH)) skeb_use_gm = STOCH%skeb_use_gm + if (skeb_use_gm) then + allocate(skeb_gm_work(is:ie,js:je), source=0.) + allocate(skeb_ebt_norm2(is:ie,js:je), source=0.) + endif + nk_linear = max(GV%nkml, 1) Slope_x_PE(:,:,:) = 0.0 @@ -795,6 +822,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV find_work = allocated(MEKE%GM_src) find_work = (allocated(CS%GMwork) .or. find_work) + find_work = (skeb_use_gm .or. find_work) if (use_EOS) then halo = 1 ! Default halo to fill is 1 @@ -1548,8 +1576,23 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (.not. CS%GM_src_alt) then ; if (allocated(MEKE%GM_src)) then MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + Work_h endif ; endif + if (skeb_use_gm) then + h_tot = sum(h(i,j,1:nz)) + skeb_gm_work(i,j) = STOCH%skeb_gm_coef * Work_h + skeb_ebt_norm2(i,j) = GV%H_to_RZ * & + (sum(h(i,j,1:nz) * VarMix%ebt_struct(i,j,1:nz)**2) + h_neglect) + endif enddo ; enddo ; endif + if (skeb_use_gm) then + ! This block spreads the GM work down through the column using the ebt vertical structure, squared. + ! Note the sign convention. + do k=1,nz ; do j=js,je ; do i=is,ie + STOCH%skeb_diss(i,j,k) = STOCH%skeb_diss(i,j,k) - skeb_gm_work(i,j) * & + VarMix%ebt_struct(i,j,k)**2 / skeb_ebt_norm2(i,j) + enddo ; enddo ; enddo + endif + if (find_work .and. CS%GM_src_alt) then ; if (allocated(MEKE%GM_src)) then do j=js,je ; do i=is,ie ; do k=nz,1,-1 PE_release_h = -0.25 * (GV%H_to_RZ*US%L_to_Z**2) * & diff --git a/src/parameterizations/stochastic/MOM_stochastics.F90 b/src/parameterizations/stochastic/MOM_stochastics.F90 index 04a29019fa..1bc42660d3 100644 --- a/src/parameterizations/stochastic/MOM_stochastics.F90 +++ b/src/parameterizations/stochastic/MOM_stochastics.F90 @@ -8,8 +8,12 @@ module MOM_stochastics ! particular version wraps all of the calls for MOM6 in the calls that had ! been used for MOM4. ! -use MOM_diag_mediator, only : register_diag_field, diag_ctrl, time_type +use MOM_debugging, only : hchksum, uvchksum, qchksum +use MOM_diag_mediator, only : register_diag_field, diag_ctrl, time_type, post_data +use MOM_diag_mediator, only : register_static_field, enable_averages, disable_averaging use MOM_grid, only : ocean_grid_type +use MOM_variables, only : thermo_var_ptrs +use MOM_domains, only : pass_var, pass_vector, CORNER, SCALAR_PAIR use MOM_verticalGrid, only : verticalGrid_type use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave @@ -18,28 +22,56 @@ module MOM_stochastics use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain use MOM_domains, only : root_PE, num_PEs use MOM_coms, only : Get_PElist +use MOM_EOS, only : calculate_density, EOS_domain use stochastic_physics, only : init_stochastic_physics_ocn, run_stochastic_physics_ocn #include implicit none ; private -public stochastics_init, update_stochastics +public stochastics_init, update_stochastics, apply_skeb !> This control structure holds parameters for the MOM_stochastics module type, public:: stochastic_CS logical :: do_sppt !< If true, stochastically perturb the diabatic + logical :: do_skeb !< If true, stochastically perturb the diabatic + logical :: skeb_use_gm !< If true, adds GM work to the amplitude of SKEBS + logical :: skeb_use_frict !< If true, adds viscous dissipation rate to the amplitude of SKEBS logical :: pert_epbl !< If true, then randomly perturb the KE dissipation and genration terms - integer :: id_sppt_wts = -1 !< Diagnostic id for SPPT - integer :: id_epbl1_wts = -1 !< Diagnostic id for epbl generation perturbation - integer :: id_epbl2_wts = -1 !< Diagnostic id for epbl dissipation perturbation + integer :: id_sppt_wts = -1 !< Diagnostic id for SPPT + integer :: id_skeb_wts = -1 !< Diagnostic id for SKEB + integer :: id_skebu = -1 !< Diagnostic id for SKEB + integer :: id_skebv = -1 !< Diagnostic id for SKEB + integer :: id_diss = -1 !< Diagnostic id for SKEB + integer :: skeb_npass = -1 !< number of passes of the 9-point smoother for the dissipation estimate + integer :: id_psi = -1 !< Diagnostic id for SPPT + integer :: id_epbl1_wts = -1 !< Diagnostic id for epbl generation perturbation + integer :: id_epbl2_wts = -1 !< Diagnostic id for epbl dissipation perturbation + integer :: id_skeb_taperu = -1 !< Diagnostic id for u taper of SKEB velocity increment + integer :: id_skeb_taperv = -1 !< Diagnostic id for v taper of SKEB velocity increment + real :: skeb_gm_coef !< If skeb_use_gm is true, then skeb_gm_coef * GM_work is added to the + !! dissipation rate used to set the amplitude of SKEBS [nondim] + real :: skeb_frict_coef !< If skeb_use_frict is true, then skeb_gm_coef * GM_work is added to the + !! dissipation rate used to set the amplitude of SKEBS [nondim] + real, allocatable :: skeb_diss(:,:,:) !< Dissipation rate used to set amplitude of SKEBS [L2 T-3 ~> m2 s-2] + !! Index into this at h points. ! stochastic patterns real, allocatable :: sppt_wts(:,:) !< Random pattern for ocean SPPT !! tendencies with a number between 0 and 2 + real, allocatable :: skeb_wts(:,:) !< Random pattern for ocean SKEB real, allocatable :: epbl1_wts(:,:) !< Random pattern for K.E. generation real, allocatable :: epbl2_wts(:,:) !< Random pattern for K.E. dissipation - type(diag_ctrl), pointer :: diag !< structure used to regulate timing of diagnostic output type(time_type), pointer :: Time !< Pointer to model time (needed for sponges) + type(diag_ctrl), pointer :: diag=>NULL() !< A structure that is used to regulate the + + ! Taper array to smoothly zero out the SKEBS velocity increment near land + real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: taperCu !< Taper applied to u component of + !! stochastic velocity increment + !! range [0,1], [nondim] + real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: taperCv !< Taper applied to v component of + !! stochastic velocity increment + !! range [0,1], [nondim] + end type stochastic_CS contains @@ -62,20 +94,24 @@ subroutine stochastics_init(dt, grid, GV, CS, param_file, diag, Time) integer :: pe_zero ! root pe integer :: nx ! number of x-points including halo integer :: ny ! number of x-points including halo + integer :: i, j, k ! loop indices + real :: tmp(grid%isdB:grid%iedB,grid%jsdB:grid%jedB) ! Used to construct tapers + integer :: taper_width ! Width (in cells) of the taper that brings the stochastic velocity + ! increments to 0 at the boundary. ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "ocean_stochastics_init" ! This module's name. - call callTree_enter("ocean_model_stochastic_init(), MOM_stochastics.F90") + call callTree_enter("stochastic_init(), MOM_stochastics.F90") if (associated(CS)) then call MOM_error(WARNING, "MOM_stochastics_init called with an "// & "associated control structure.") return else ; allocate(CS) ; endif - CS%diag => diag CS%Time => Time + CS%diag => diag ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") @@ -83,48 +119,130 @@ subroutine stochastics_init(dt, grid, GV, CS, param_file, diag, Time) ! get number of processors and PE list for stochastic physics initialization call get_param(param_file, mdl, "DO_SPPT", CS%do_sppt, & "If true, then stochastically perturb the thermodynamic "//& - "tendemcies of T,S, amd h. Amplitude and correlations are "//& + "tendencies 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, "DO_SKEB", CS%do_skeb, & + "If true, then stochastically perturb the currents "//& + "using the stochastic kinetic energy backscatter scheme.",& + default=.false.) + call get_param(param_file, mdl, "SKEB_NPASS", CS%skeb_npass, & + "number of passes of a 9-point smoother of the "//& + "dissipation estimate.", default=3, do_not_log=.not.CS%do_skeb) + call get_param(param_file, mdl, "SKEB_TAPER_WIDTH", taper_width, & + "number of cells over which the stochastic velocity increment "//& + "is tapered to zero.", default=4, do_not_log=.not.CS%do_skeb) + call get_param(param_file, mdl, "SKEB_USE_GM", CS%skeb_use_gm, & + "If true, adds GM work rate to the SKEBS amplitude.", & + default=.false., do_not_log=.not.CS%do_skeb) + if ((.not. CS%do_skeb) .and. (CS%skeb_use_gm)) call MOM_error(FATAL, "If SKEB_USE_GM is True "//& + "then DO_SKEB must also be True.") + call get_param(param_file, mdl, "SKEB_GM_COEF", CS%skeb_gm_coef, & + "Fraction of GM work that is added to backscatter rate.", & + units="nondim", default=0.0, do_not_log=.not.CS%skeb_use_gm) + call get_param(param_file, mdl, "SKEB_USE_FRICT", CS%skeb_use_frict, & + "If true, adds horizontal friction dissipation rate "//& + "to the SKEBS amplitude.", default=.false., do_not_log=.not.CS%do_skeb) + if ((.not. CS%do_skeb) .and. (CS%skeb_use_frict)) call MOM_error(FATAL, "If SKEB_USE_FRICT is "//& + "True then DO_SKEB must also be True.") + call get_param(param_file, mdl, "SKEB_FRICT_COEF", CS%skeb_frict_coef, & + "Fraction of horizontal friction work that is added to backscatter rate.", & + units="nondim", default=0.0, do_not_log=.not.CS%skeb_use_frict) call get_param(param_file, mdl, "PERT_EPBL", CS%pert_epbl, & "If true, then stochastically perturb the kinetic energy "//& "production and dissipation terms. Amplitude and correlations are "//& "controlled by the nam_stoch namelist in the UFS model only.", & default=.false.) - if (CS%do_sppt .OR. CS%pert_epbl) then - num_procs = num_PEs() - allocate(pelist(num_procs)) - call Get_PElist(pelist,commID = mom_comm) - pe_zero = root_PE() - nx = grid%ied - grid%isd + 1 - ny = grid%jed - grid%jsd + 1 - call init_stochastic_physics_ocn(dt,grid%geoLonT,grid%geoLatT,nx,ny,GV%ke, & - CS%pert_epbl,CS%do_sppt,pe_zero,mom_comm,iret) - if (iret/=0) then - call MOM_error(FATAL, "call to init_stochastic_physics_ocn failed") - endif - - if (CS%do_sppt) allocate(CS%sppt_wts(grid%isd:grid%ied,grid%jsd:grid%jed), source=0.0) - if (CS%pert_epbl) then - allocate(CS%epbl1_wts(grid%isd:grid%ied,grid%jsd:grid%jed), source=0.0) - allocate(CS%epbl2_wts(grid%isd:grid%ied,grid%jsd:grid%jed), source=0.0) - endif - endif - if (CS%do_sppt) then - CS%id_sppt_wts = register_diag_field('ocean_model', 'sppt_pattern', CS%diag%axesT1, Time, & - 'random pattern for sppt', 'None') + + if (CS%do_sppt .OR. CS%pert_epbl .OR. CS%do_skeb) then + num_procs = num_PEs() + allocate(pelist(num_procs)) + call Get_PElist(pelist,commID = mom_comm) + pe_zero = root_PE() + nx = grid%iedB - grid%isdB + 1 + ny = grid%jedB - grid%jsdB + 1 + call init_stochastic_physics_ocn(dt,grid%geoLonBu,grid%geoLatBu,nx,ny,GV%ke, & + CS%pert_epbl,CS%do_sppt,CS%do_skeb,pe_zero,mom_comm,iret) + if (iret/=0) then + call MOM_error(FATAL, "call to init_stochastic_physics_ocn failed") + return + endif + + if (CS%do_sppt) allocate(CS%sppt_wts(grid%isdB:grid%iedB,grid%jsdB:grid%jedB)) + if (CS%do_skeb) allocate(CS%skeb_wts(grid%isdB:grid%iedB,grid%jsdB:grid%jedB)) + if (CS%do_skeb) allocate(CS%skeb_diss(grid%isd:grid%ied,grid%jsd:grid%jed,GV%ke), source=0.) + if (CS%pert_epbl) then + allocate(CS%epbl1_wts(grid%isdB:grid%iedB,grid%jsdB:grid%jedB)) + allocate(CS%epbl2_wts(grid%isdB:grid%iedB,grid%jsdB:grid%jedB)) + endif endif - if (CS%pert_epbl) then - CS%id_epbl1_wts = register_diag_field('ocean_model', 'epbl1_wts', CS%diag%axesT1, Time, & - 'random pattern for KE generation', 'None') - CS%id_epbl2_wts = register_diag_field('ocean_model', 'epbl2_wts', CS%diag%axesT1, Time, & - 'random pattern for KE dissipation', 'None') + + CS%id_sppt_wts = register_diag_field('ocean_model', 'sppt_pattern', CS%diag%axesB1, Time, & + 'random pattern for sppt', 'None') + CS%id_skeb_wts = register_diag_field('ocean_model', 'skeb_pattern', CS%diag%axesB1, Time, & + 'random pattern for skeb', 'None') + CS%id_epbl1_wts = register_diag_field('ocean_model', 'epbl1_wts', CS%diag%axesB1, Time, & + 'random pattern for KE generation', 'None') + CS%id_epbl2_wts = register_diag_field('ocean_model', 'epbl2_wts', CS%diag%axesB1, Time, & + 'random pattern for KE dissipation', 'None') + CS%id_skebu = register_diag_field('ocean_model', 'skebu', CS%diag%axesCuL, Time, & + 'zonal current perts', 'None') + CS%id_skebv = register_diag_field('ocean_model', 'skebv', CS%diag%axesCvL, Time, & + 'zonal current perts', 'None') + CS%id_diss = register_diag_field('ocean_model', 'skeb_amp', CS%diag%axesTL, Time, & + 'SKEB amplitude', 'm s-1') + CS%id_psi = register_diag_field('ocean_model', 'psi', CS%diag%axesBL, Time, & + 'stream function', 'None') + CS%id_skeb_taperu = register_static_field('ocean_model', 'skeb_taper_u', CS%diag%axesCu1, & + 'SKEB taper u', 'None', interp_method='none') + CS%id_skeb_taperv = register_static_field('ocean_model', 'skeb_taper_v', CS%diag%axesCv1, & + 'SKEB taper v', 'None', interp_method='none') + + ! Initialize the "taper" fields. These fields multiply the components of the stochastic + ! velocity increment in such a way as to smoothly taper them to zero at land boundaries. + if ((CS%do_skeb) .or. (CS%id_skeb_taperu > 0) .or. (CS%id_skeb_taperv > 0)) then + ALLOC_(CS%taperCu(grid%IsdB:grid%IedB,grid%jsd:grid%jed)) + ALLOC_(CS%taperCv(grid%isd:grid%ied,grid%JsdB:grid%JedB)) + ! Initialize taper from land mask + do j=grid%jsd,grid%jed ; do I=grid%isdB,grid%iedB + CS%taperCu(I,j) = grid%mask2dCu(I,j) + enddo ; enddo + do J=grid%jsdB,grid%jedB ; do i=grid%isd,grid%ied + CS%taperCv(i,J) = grid%mask2dCv(i,J) + enddo ; enddo + ! Extend taper land + do k=1,(taper_width / 2) + do j=grid%jsc-1,grid%jec+1 ; do I=grid%iscB-1,grid%iecB+1 + tmp(I,j) = minval(CS%taperCu(I-1:I+1,j-1:j+1)) + enddo ; enddo + do j=grid%jsc,grid%jec ; do I=grid%iscB,grid%iecB + CS%taperCu(I,j) = minval(tmp(I-1:I+1,j-1:j+1)) + enddo ; enddo + do J=grid%jscB-1,grid%jecB+1 ; do i=grid%isc-1,grid%iec+1 + tmp(i,J) = minval(CS%taperCv(i-1:i+1,J-1:J+1)) + enddo ; enddo + do J=grid%jscB,grid%jecB ; do i=grid%isc,grid%iec + CS%taperCv(i,J) = minval(tmp(i-1:i+1,J-1:J+1)) + enddo ; enddo + ! Update halo + call pass_vector(CS%taperCu, CS%taperCv, grid%Domain, SCALAR_PAIR) + enddo + ! Smooth tapers. Each call smooths twice. + do k=1,(taper_width - (taper_width/2)) + call smooth_x9_uv(grid, CS%taperCu, CS%taperCv, zero_land=.true.) + call pass_vector(CS%taperCu, CS%taperCv, grid%Domain, SCALAR_PAIR) + enddo endif - if (CS%do_sppt .OR. CS%pert_epbl) & + !call uvchksum("SKEB taper [uv]", CS%taperCu, CS%taperCv, grid%HI) + + if (CS%id_skeb_taperu > 0) call post_data(CS%id_skeb_taperu, CS%taperCu, CS%diag, .true.) + if (CS%id_skeb_taperv > 0) call post_data(CS%id_skeb_taperv, CS%taperCv, CS%diag, .true.) + + if (CS%do_sppt .OR. CS%pert_epbl .OR. CS%do_skeb) & call MOM_mesg(' === COMPLETED MOM STOCHASTIC INITIALIZATION =====') - call callTree_leave("ocean_model_init(") + call callTree_leave("stochastic_init(), MOM_stochastics.F90") end subroutine stochastics_init @@ -138,10 +256,202 @@ subroutine update_stochastics(CS) call callTree_enter("update_stochastics(), MOM_stochastics.F90") ! update stochastic physics patterns before running next time-step - call run_stochastic_physics_ocn(CS%sppt_wts,CS%epbl1_wts,CS%epbl2_wts) + call run_stochastic_physics_ocn(CS%sppt_wts,CS%skeb_wts,CS%epbl1_wts,CS%epbl2_wts) + + call callTree_leave("update_stochastics(), MOM_stochastics.F90") - return end subroutine update_stochastics +subroutine apply_skeb(grid,GV,CS,uc,vc,thickness,tv,dt,Time_end) + + type(ocean_grid_type), intent(in) :: grid !< ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid + type(stochastic_CS), intent(inout) :: CS !< stochastic control structure + + real, dimension(SZIB_(grid),SZJ_(grid),SZK_(GV)), intent(inout) :: uc !< zonal velocity [L T-1 ~> m s-1] + real, dimension(SZI_(grid),SZJB_(grid),SZK_(GV)), intent(inout) :: vc !< meridional velocity [L T-1 ~> m s-1] + real, dimension(SZI_(grid),SZJ_(grid),SZK_(GV)), intent(in) :: thickness !< thickness [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(in) :: tv !< points to thermodynamic fields + real, intent(in) :: dt !< time increment [T ~> s] + type(time_type), intent(in) :: Time_end !< Time at the end of the interval +! locals + + real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_,NKMEM_) :: psi + real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: ustar + real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: vstar + real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: diss_tmp + + real, dimension(3,3) :: local_weights + + real :: shr,ten,tot,kh + integer :: i,j,k,iter + integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state + + call callTree_enter("apply_skeb(), MOM_stochastics.F90") + ALLOC_(diss_tmp(grid%isd:grid%ied,grid%jsd:grid%jed)) + ALLOC_(psi(grid%isdB:grid%iedB,grid%jsdB:grid%jedB,GV%ke)) + ALLOC_(ustar(grid%isdB:grid%iedB,grid%jsd:grid%jed,GV%ke)) + ALLOC_(vstar(grid%isd:grid%ied,grid%jsdB:grid%jedB,GV%ke)) + + if ((.not. CS%skeb_use_gm) .and. (.not. CS%skeb_use_frict)) then + ! fill in halos with zeros + do k=1,GV%ke + do j=grid%jsd,grid%jed ; do i=grid%isd,grid%ied + CS%skeb_diss(i,j,k) = 0.0 + enddo ; enddo + enddo + + !kh needs to be scaled + + kh=1!(120*111)**2 + do k=1,GV%ke + do j=grid%jsc,grid%jec ; do i=grid%isc,grid%iec + ! Shear + shr = (vc(i,J,k)-vc(i-1,J,k))*grid%mask2dCv(i,J)*grid%mask2dCv(i-1,J)*grid%IdxCv(i,J)+& + (uc(I,j,k)-uc(I,j-1,k))*grid%mask2dCu(I,j)*grid%mask2dCu(I,j-1)*grid%IdyCu(I,j) + ! Tension + ten = (vc(i,J,k)-vc(i-1,J,k))*grid%mask2dCv(i,J)*grid%mask2dCv(i-1,J)*grid%IdyCv(i,J)+& + (uc(I,j,k)-uc(I,j-1,k))*grid%mask2dCu(I,j)*grid%mask2dCu(I,j-1)*grid%IdxCu(I,j) + + tot = sqrt( shr**2 + ten**2 ) * grid%mask2dT(i,j) + CS%skeb_diss(i,j,k) = tot**3 * kh * grid%areaT(i,j)!!**2 + enddo ; enddo + enddo + endif ! Sets CS%skeb_diss without GM or FrictWork + + ! smooth dissipation skeb_npass times + do iter=1,CS%skeb_npass + if (mod(iter,2) == 1) call pass_var(CS%skeb_diss, grid%domain) + do k=1,GV%ke + do j=grid%jsc-1,grid%jec+1 ; do i=grid%isc-1,grid%iec+1 + ! This does not preserve rotational symmetry + local_weights = grid%mask2dT(i-1:i+1,j-1:j+1)*grid%areaT(i-1:i+1,j-1:j+1) + diss_tmp(i,j) = sum(local_weights*CS%skeb_diss(i-1:i+1,j-1:j+1,k)) / & + (sum(local_weights) + 1.E-16) + enddo ; enddo + do j=grid%jsc-1,grid%jec+1 ; do i=grid%isc-1,grid%iec+1 + if (grid%mask2dT(i,j)==0.) cycle + CS%skeb_diss(i,j,k) = diss_tmp(i,j) + enddo ; enddo + enddo + enddo + call pass_var(CS%skeb_diss, grid%domain) + + ! call hchksum(CS%skeb_diss, "SKEB DISS", grid%HI, haloshift=2) + ! call qchksum(CS%skeb_wts, "SKEB WTS", grid%HI, haloshift=1) + + do k=1,GV%ke + do J=grid%jscB-1,grid%jecB ; do I=grid%iscB-1,grid%iecB + psi(I,J,k) = sqrt(0.25 * dt * max((CS%skeb_diss(i ,j ,k) + CS%skeb_diss(i+1,j+1,k)) + & + (CS%skeb_diss(i ,j+1,k) + CS%skeb_diss(i+1,j ,k)), 0.) ) & + * CS%skeb_wts(I,J) + enddo ; enddo + enddo + !call qchksum(psi,"SKEB PSI", grid%HI, haloshift=1) + !call pass_var(psi, grid%domain, position=CORNER) + do k=1,GV%ke + do j=grid%jsc,grid%jec ; do I=grid%iscB,grid%iecB + ustar(I,j,k) = - (psi(I,J,k) - psi(I,J-1,k)) * CS%taperCu(I,j) * grid%IdyCu(I,j) + uc(I,j,k) = uc(I,j,k) + ustar(I,j,k) + enddo ; enddo + do J=grid%jscB,grid%jecB ; do i=grid%isc,grid%iec + vstar(i,J,k) = (psi(I,J,k) - psi(I-1,J,k)) * CS%taperCv(i,J) * grid%IdxCv(i,J) + vc(i,J,k) = vc(i,J,k) + vstar(i,J,k) + enddo ; enddo + enddo + + !call uvchksum("SKEB increment [uv]", ustar, vstar, grid%HI) + + call enable_averages(dt, Time_end, CS%diag) + if (CS%id_diss > 0) then + call post_data(CS%id_diss, sqrt(dt * max(CS%skeb_diss(:,:,:), 0.)), CS%diag) + endif + if (CS%id_skeb_wts > 0) then + call post_data(CS%id_skeb_wts, CS%skeb_wts, CS%diag) + endif + if (CS%id_skebu > 0) then + call post_data(CS%id_skebu, ustar(:,:,:), CS%diag) + endif + if (CS%id_skebv > 0) then + call post_data(CS%id_skebv, vstar(:,:,:), CS%diag) + endif + if (CS%id_psi > 0) then + call post_data(CS%id_psi, psi(:,:,:), CS%diag) + endif + call disable_averaging(CS%diag) + DEALLOC_(diss_tmp) + DEALLOC_(ustar) + DEALLOC_(vstar) + DEALLOC_(psi) + CS%skeb_diss(:,:,:) = 0.0 ! Must zero before next time step. + + call callTree_leave("apply_skeb(), MOM_stochastics.F90") + +end subroutine apply_skeb + +!> Apply a 9-point smoothing filter twice to a pair of velocity components to reduce +!! horizontal two-grid-point noise. +!! Note that this subroutine does not conserve angular momentum, so don't use it +!! in situations where you need conservation. Also note that it assumes that the +!! input fields have valid values in the first two halo points upon entry. +subroutine smooth_x9_uv(G, field_u, field_v, zero_land) + type(ocean_grid_type), intent(in) :: G !< Ocean grid + real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: field_u !< u-point field to be smoothed[arbitrary] + real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: field_v !< v-point field to be smoothed [arbitrary] + logical, optional, intent(in) :: zero_land !< If present and false, return the average + !! of the surrounding ocean points when + !! smoothing, otherwise use a value of 0 for + !! land points and include them in the averages. + + ! Local variables. + real :: fu_prev(SZIB_(G),SZJ_(G)) ! The value of the u-point field at the previous iteration [arbitrary] + real :: fv_prev(SZI_(G),SZJB_(G)) ! The value of the v-point field at the previous iteration [arbitrary] + real :: Iwts ! The inverse of the sum of the weights [nondim] + logical :: zero_land_val ! The value of the zero_land optional argument or .true. if it is absent. + integer :: i, j, s, is, ie, js, je, Isq, Ieq, Jsq, Jeq + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + + zero_land_val = .true. ; if (present(zero_land)) zero_land_val = zero_land + + do s=1,0,-1 + fu_prev(:,:) = field_u(:,:) + ! apply smoothing on field_u using rotationally symmetric expressions. + do j=js-s,je+s ; do I=Isq-s,Ieq+s ; if (G%mask2dCu(I,j) > 0.0) then + Iwts = 0.0625 + if (.not. zero_land_val) & + Iwts = 1.0 / ( (4.0*G%mask2dCu(I,j) + & + ( 2.0*((G%mask2dCu(I-1,j) + G%mask2dCu(I+1,j)) + & + (G%mask2dCu(I,j-1) + G%mask2dCu(I,j+1))) + & + ((G%mask2dCu(I-1,j-1) + G%mask2dCu(I+1,j+1)) + & + (G%mask2dCu(I-1,j+1) + G%mask2dCu(I+1,j-1))) ) ) + 1.0e-16 ) + field_u(I,j) = Iwts * ( 4.0*G%mask2dCu(I,j) * fu_prev(I,j) & + + (2.0*((G%mask2dCu(I-1,j) * fu_prev(I-1,j) + G%mask2dCu(I+1,j) * fu_prev(I+1,j)) + & + (G%mask2dCu(I,j-1) * fu_prev(I,j-1) + G%mask2dCu(I,j+1) * fu_prev(I,j+1))) & + + ((G%mask2dCu(I-1,j-1) * fu_prev(I-1,j-1) + G%mask2dCu(I+1,j+1) * fu_prev(I+1,j+1)) + & + (G%mask2dCu(I-1,j+1) * fu_prev(I-1,j+1) + G%mask2dCu(I+1,j-1) * fu_prev(I-1,j-1))) )) + endif ; enddo ; enddo + + fv_prev(:,:) = field_v(:,:) + ! apply smoothing on field_v using rotationally symmetric expressions. + do J=Jsq-s,Jeq+s ; do i=is-s,ie+s ; if (G%mask2dCv(i,J) > 0.0) then + Iwts = 0.0625 + if (.not. zero_land_val) & + Iwts = 1.0 / ( (4.0*G%mask2dCv(i,J) + & + ( 2.0*((G%mask2dCv(i-1,J) + G%mask2dCv(i+1,J)) + & + (G%mask2dCv(i,J-1) + G%mask2dCv(i,J+1))) + & + ((G%mask2dCv(i-1,J-1) + G%mask2dCv(i+1,J+1)) + & + (G%mask2dCv(i-1,J+1) + G%mask2dCv(i+1,J-1))) ) ) + 1.0e-16 ) + field_v(i,J) = Iwts * ( 4.0*G%mask2dCv(i,J) * fv_prev(i,J) & + + (2.0*((G%mask2dCv(i-1,J) * fv_prev(i-1,J) + G%mask2dCv(i+1,J) * fv_prev(i+1,J)) + & + (G%mask2dCv(i,J-1) * fv_prev(i,J-1) + G%mask2dCv(i,J+1) * fv_prev(i,J+1))) & + + ((G%mask2dCv(i-1,J-1) * fv_prev(i-1,J-1) + G%mask2dCv(i+1,J+1) * fv_prev(i+1,J+1)) + & + (G%mask2dCv(i-1,J+1) * fv_prev(i-1,J+1) + G%mask2dCv(i+1,J-1) * fv_prev(i-1,J-1))) )) + endif ; enddo ; enddo + enddo + +end subroutine smooth_x9_uv + end module MOM_stochastics diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 8e95edd563..832d8bf4b1 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -1342,47 +1342,60 @@ subroutine KPP_smooth_BLD(CS, G, GV, US, dz) type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: dz !< Layer thicknesses [Z ~> m] - ! local + ! local variables real, dimension(SZI_(G),SZJ_(G)) :: OBLdepth_prev ! OBLdepth before s.th smoothing iteration [Z ~> m] + real, dimension(SZI_(G),SZJ_(G)) :: total_depth ! The total depth of the water column, adjusted + ! for the minimum layer thickness [Z ~> m] real, dimension( GV%ke ) :: cellHeight ! Cell center heights referenced to surface [Z ~> m] ! (negative in the ocean) real, dimension( GV%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface [Z ~> m] ! (negative in the ocean) real :: wc, ww, we, wn, ws ! averaging weights for smoothing [nondim] real :: dh ! The local thickness used for calculating interface positions [Z ~> m] + real :: h_cor(SZI_(G)) ! A cumulative correction arising from inflation of vanished layers [Z ~> m] real :: hcorr ! A cumulative correction arising from inflation of vanished layers [Z ~> m] - integer :: i, j, k, s + integer :: i, j, k, s, halo call cpu_clock_begin(id_clock_KPP_smoothing) - ! Update halos + ! Find the total water column thickness first, as it is reused for each smoothing pass. + total_depth(:,:) = 0.0 + + !$OMP parallel do default(shared) private(dh, h_cor) + do j = G%jsc, G%jec + h_cor(:) = 0. + do k=1,GV%ke + do i=G%isc,G%iec ; if (G%mask2dT(i,j) > 0.0) then + ! This code replicates the interface height calculations below. It could be simpler, as shown below. + dh = dz(i,j,k) ! Nominal thickness to use for increment + dh = dh + h_cor(i) ! Take away the accumulated error (could temporarily make dh<0) + h_cor(i) = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0 + dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness + total_depth(i,j) = total_depth(i,j) + dh + endif ; enddo + enddo + enddo + ! A much simpler (but answer changing) version of the total_depth calculation would be + ! do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + ! total_depth(i,j) = total_depth(i,j) + dz(i,j,k) + ! enddo ; enddo ; enddo + + ! Update halos once, then march inward for each iteration + if (CS%n_smooth > 1) call pass_var(total_depth, G%Domain, halo=CS%n_smooth, complete=.false.) call pass_var(CS%OBLdepth, G%Domain, halo=CS%n_smooth) - if (CS%id_OBLdepth_original > 0) CS%OBLdepth_original = CS%OBLdepth + if (CS%id_OBLdepth_original > 0) CS%OBLdepth_original(:,:) = CS%OBLdepth(:,:) do s=1,CS%n_smooth - OBLdepth_prev = CS%OBLdepth + OBLdepth_prev(:,:) = CS%OBLdepth(:,:) + halo = CS%n_smooth - s ! apply smoothing on OBL depth - !$OMP parallel do default(none) shared(G, GV, US, CS, dz, OBLdepth_prev) & - !$OMP private(wc, ww, we, wn, ws, dh, hcorr, cellHeight, iFaceHeight) - do j = G%jsc, G%jec - do i = G%isc, G%iec ; if (G%mask2dT(i,j) > 0.0) then - - iFaceHeight(1) = 0.0 ! BBL is all relative to the surface - hcorr = 0. - do k=1,GV%ke - - ! cell center and cell bottom in meters (negative values in the ocean) - dh = dz(i,j,k) ! Nominal thickness to use for increment - dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0) - hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0 - dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness - cellHeight(k) = iFaceHeight(k) - 0.5 * dh - iFaceHeight(k+1) = iFaceHeight(k) - dh - enddo - + !$OMP parallel do default(none) shared(G, GV, CS, OBLdepth_prev, total_depth, halo) & + !$OMP private(wc, ww, we, wn, ws) + do j = G%jsc-halo, G%jec+halo + do i = G%isc-halo, G%iec+halo ; if (G%mask2dT(i,j) > 0.0) then ! compute weights ww = 0.125 * G%mask2dT(i-1,j) we = 0.125 * G%mask2dT(i+1,j) @@ -1400,19 +1413,37 @@ subroutine KPP_smooth_BLD(CS, G, GV, US, dz) if (CS%deepen_only) CS%OBLdepth(i,j) = max(CS%OBLdepth(i,j), OBLdepth_prev(i,j)) ! prevent OBL depths deeper than the bathymetric depth - CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -iFaceHeight(GV%ke+1) ) ! no deeper than bottom - CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) ) + CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), total_depth(i,j) ) ! no deeper than bottom endif ; enddo enddo enddo ! s-loop + ! Determine the fractional index of the bottom of the boundary layer. + !$OMP parallel do default(none) shared(G, GV, CS, dz) & + !$OMP private(dh, hcorr, cellHeight, iFaceHeight) + do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (G%mask2dT(i,j) > 0.0) then + + iFaceHeight(1) = 0.0 ! BBL is all relative to the surface + hcorr = 0. + do k=1,GV%ke + ! cell center and cell bottom in meters (negative values in the ocean) + dh = dz(i,j,k) ! Nominal thickness to use for increment + dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0) + hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0 + dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness + cellHeight(k) = iFaceHeight(k) - 0.5 * dh + iFaceHeight(k+1) = iFaceHeight(k) - dh + enddo + + CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) ) + endif ; enddo ; enddo + call cpu_clock_end(id_clock_KPP_smoothing) end subroutine KPP_smooth_BLD - !> Copies KPP surface boundary layer depth into BLD, in units of [Z ~> m] unless other units are specified. subroutine KPP_get_BLD(CS, BLD, G, US, m_to_BLD_units) type(KPP_CS), pointer :: CS !< Control structure for diff --git a/src/parameterizations/vertical/MOM_opacity.F90 b/src/parameterizations/vertical/MOM_opacity.F90 index 61a7a0c7d0..831607d2db 100644 --- a/src/parameterizations/vertical/MOM_opacity.F90 +++ b/src/parameterizations/vertical/MOM_opacity.F90 @@ -40,7 +40,18 @@ module MOM_opacity real :: PenSW_flux_absorb !< A heat flux that is small enough to be completely absorbed in the next !! sufficiently thick layer [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]. real :: PenSW_absorb_Invlen !< The inverse of the thickness that is used to absorb the remaining - !! shortwave heat flux when it drops below PEN_SW_FLUX_ABSORB [H ~> m or kg m-2]. + !! shortwave heat flux when it drops below PEN_SW_FLUX_ABSORB [H ~> m or kg m-2]. + + !! Lookup tables for Ohlmann solar penetration scheme + !! These would naturally exist as private module variables but that is prohibited in MOM6 + real :: dlog10chl !< Chl increment within lookup table + real :: log10chl_min !< Lower bound of Chl in lookup table + real :: log10chl_max !< Upper bound of Chl in lookup table + real, allocatable, dimension(:) :: a1_lut,& !< Coefficient for band 1 + & a2_lut,& !< Coefficient for band 2 + & b1_lut,& !< Exponential decay scale for band 1 + & b2_lut !< Exponential decay scale for band 2 + integer :: answer_date !< The vintage of the order of arithmetic and expressions in the optics !! calculations. Values below 20190101 recover the answers from the !! end of 2018, while higher values use updated and more robust @@ -77,11 +88,13 @@ module MOM_opacity end type opacity_CS !>@{ Coded integers to specify the opacity scheme -integer, parameter :: NO_SCHEME = 0, MANIZZA_05 = 1, MOREL_88 = 2, SINGLE_EXP = 3, DOUBLE_EXP = 4 +integer, parameter :: NO_SCHEME = 0, MANIZZA_05 = 1, MOREL_88 = 2, SINGLE_EXP = 3, DOUBLE_EXP = 4,& + & OHLMANN_03 = 5 !>@} character*(10), parameter :: MANIZZA_05_STRING = "MANIZZA_05" !< String to specify the opacity scheme character*(10), parameter :: MOREL_88_STRING = "MOREL_88" !< String to specify the opacity scheme +character*(10), parameter :: OHLMANN_03_STRING = "OHLMANN_03" !< String to specify the opacity scheme character*(10), parameter :: SINGLE_EXP_STRING = "SINGLE_EXP" !< String to specify the opacity scheme character*(10), parameter :: DOUBLE_EXP_STRING = "DOUBLE_EXP" !< String to specify the opacity scheme @@ -254,6 +267,16 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir ! use the "blue" band in the parameterizations to determine the e-folding ! depth of the incoming shortwave attenuation. The red portion is lumped ! into the net heating at the surface. +! Adding Ohlmann scheme. Needs sw_total and chl as inputs. Produces 2 penetrating bands. +! This implementation follows that in CESM-POP using a lookup table in log10(chl) space. +! The table is initialized in subroutine init_ohlmann and the coefficients are recovered +! with routines lookup_ohlmann_swpen and lookup_ohlmann_opacity. +! Note that this form treats the IR solar input implicitly: the sum of partioning +! coefficients < 1.0. The remainder is non-penetrating and is deposited in first layer +! irrespective of thickness. The Ohlmann (2003) paper states that the scheme is not valid +! for vertcal grids with first layer thickness < 2.0 meters. +! +! Ohlmann, J.C. Ocean radiant heating in climate models. J. Climate, 16, 1337-1351, 2003. ! ! Morel, A., Optical modeling of the upper ocean in relation to its biogenous ! matter content (case-i waters)., J. Geo. Res., {93}, 10,749--10,768, 1988. @@ -353,13 +376,44 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir do n=1,nbands optics%sw_pen_band(n,i,j) = Inv_nbands*sw_pen_tot enddo - enddo ; enddo + enddo; enddo + case (OHLMANN_03) + ! want exactly two penetrating bands. If not, throw an error. + if ( nbands /= 2 ) then + call MOM_error(FATAL, "opacity_from_chl: CS%opacity_scheme requires nbands==2.") + endif + !$OMP parallel do default(shared) private(SW_vis_tot) + do j=js,je ; do i=is,ie + SW_vis_tot = 0.0 ! Ohlmann does not classify as vis/nir. Using vis to add up total + if (G%mask2dT(i,j) < 0.5) then + optics%sw_pen_band(1:2,i,j) = 0. ! Make sure there is a valid value for land points + else + if (multiband_vis_input ) then ! If multiband_vis_input is true then so is multiband_nir_input + SW_vis_tot = sw_vis_dir(i,j) + sw_vis_dif(i,j) + & + & sw_nir_dir(i,j) + sw_nir_dif(i,j) + elseif (total_sw_input) then + SW_vis_tot = sw_total(i,j) + else + call MOM_error(FATAL, "No shortwave input was provided.") + endif + + ! Bands 1-2 (Ohlmann factors A with coefficients for Table 1a) + optics%sw_pen_band(1:2,i,j) = lookup_ohlmann_swpen(chl_data(i,j),optics)*SW_vis_tot + endif + enddo; enddo case default call MOM_error(FATAL, "opacity_from_chl: CS%opacity_scheme is not valid.") end select !$OMP parallel do default(shared) firstprivate(chl_data) do k=1,nz + !! FOB + !!! I don't think this is what we want to do with Ohlmann. + !!! The surface CHL is used in developing the parameterization. + !!! Only the surface CHL is used above in setting optics%sw_pen_band for all schemes. + !!! Seems inconsistent to use depth dependent CHL in opacity calculation. + !!! Nevertheless, leaving as is for now. + !! FOB if (present(chl_3d)) then do j=js,je ; do i=is,ie ; chl_data(i,j) = chl_3d(i,j,k) ; enddo ; enddo endif @@ -389,14 +443,22 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir do n=2,optics%nbands optics%opacity_band(n,i,j,k) = optics%opacity_band(1,i,j,k) enddo - enddo ; enddo - + enddo; enddo + case (OHLMANN_03) + !! not testing for 2 bands since we did it above + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) <= 0.5) then + optics%opacity_band(1:2,i,j,k) = CS%opacity_land_value + else + ! Bands 1-2 (Ohlmann factors B with coefficients for Table 1a + optics%opacity_band(1:2,i,j,k) = lookup_ohlmann_opacity(chl_data(i,j),optics) * US%Z_to_m + endif + enddo; enddo case default call MOM_error(FATAL, "opacity_from_chl: CS%opacity_scheme is not valid.") end select enddo - end subroutine opacity_from_chl !> This sets the blue-wavelength opacity according to the scheme proposed by @@ -998,7 +1060,8 @@ subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics) "concentrations are translated into opacities. Currently "//& "valid options include:\n"//& " \t\t MANIZZA_05 - Use Manizza et al., GRL, 2005. \n"//& - " \t\t MOREL_88 - Use Morel, JGR, 1988.", & + " \t\t MOREL_88 - Use Morel, JGR, 1988. \n"//& + " \t\t OHLMANN_03 - Use Ohlmann, J Clim, 2003.", & default=MANIZZA_05_STRING) if (len_trim(tmpstr) > 0) then tmpstr = uppercase(tmpstr) @@ -1007,6 +1070,8 @@ subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics) CS%opacity_scheme = MANIZZA_05 ; scheme_string = MANIZZA_05_STRING case (MOREL_88_STRING) CS%opacity_scheme = MOREL_88 ; scheme_string = MOREL_88_STRING + case (OHLMANN_03_STRING) + CS%opacity_scheme = OHLMANN_03 ; scheme_string = OHLMANN_03_STRING case default call MOM_error(FATAL, "opacity_init: #DEFINE OPACITY_SCHEME "//& trim(tmpstr) // "in input file is invalid.") @@ -1072,6 +1137,9 @@ subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics) elseif (CS%Opacity_scheme == SINGLE_EXP ) then if (optics%nbands /= 1) call MOM_error(FATAL, & "set_opacity: \Cannot use a single_exp opacity scheme with nbands!=1.") + elseif (CS%Opacity_scheme == OHLMANN_03 ) then + if (optics%nbands /= 2) call MOM_error(FATAL, & + "set_opacity: \OHLMANN_03 scheme requires nbands==2") endif call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & @@ -1143,8 +1211,175 @@ subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics) longname, 'm-1', conversion=US%m_to_Z) enddo + !! FOB + if (CS%opacity_scheme == OHLMANN_03) then + ! Set up the lookup table + call init_ohlmann_table(optics) + endif + !! FOB + end subroutine opacity_init +!> Initialize the lookup table for Ohlmann solar penetration scheme. +!! Step size in Chl is a constant in log-space to make lookups easy. +!! Step size is fine enough that nearest neighbor lookup is sufficiently +!! accurate. +subroutine init_ohlmann_table(optics) + + implicit none + + type(optics_type), intent(inout) :: optics + + ! Local variables + + !! These are the data from Ohlmann (2003) Table 1a with additional + !! values provided by C. Ohlmann and implemented in CESM-POP by B. Briegleb + integer, parameter :: nval_tab1a = 31 + real, parameter, dimension(nval_tab1a) :: & + chl_tab1a = (/ & + .001, .005, .01, .02, & + .03, .05, .10, .15, & + .20, .25, .30, .35, & + .40, .45, .50, .60, & + .70, .80, .90, 1.00, & + 1.50, 2.00, 2.50, 3.00, & + 4.00, 5.00, 6.00, 7.00, & + 8.00, 9.00, 10.00 /) + + real, parameter, dimension(nval_tab1a) :: & + a1_tab1a = (/ & + 0.4421, 0.4451, 0.4488, 0.4563, & + 0.4622, 0.4715, 0.4877, 0.4993, & + 0.5084, 0.5159, 0.5223, 0.5278, & + 0.5326, 0.5369, 0.5408, 0.5474, & + 0.5529, 0.5576, 0.5615, 0.5649, & + 0.5757, 0.5802, 0.5808, 0.5788, & + 0.56965, 0.55638, 0.54091, 0.52442, & + 0.50766, 0.49110, 0.47505 /) + + real, parameter, dimension(nval_tab1a) :: & + a2_tab1a = (/ & + 0.2981, 0.2963, 0.2940, 0.2894, & + 0.2858, 0.2800, 0.2703, 0.2628, & + 0.2571, 0.2523, 0.2481, 0.2444, & + 0.2411, 0.2382, 0.2356, 0.2309, & + 0.2269, 0.2235, 0.2206, 0.2181, & + 0.2106, 0.2089, 0.2113, 0.2167, & + 0.23357, 0.25504, 0.27829, 0.30274, & + 0.32698, 0.35056, 0.37303 /) + + real, parameter, dimension(nval_tab1a) :: & + b1_tab1a = (/ & + 0.0287, 0.0301, 0.0319, 0.0355, & + 0.0384, 0.0434, 0.0532, 0.0612, & + 0.0681, 0.0743, 0.0800, 0.0853, & + 0.0902, 0.0949, 0.0993, 0.1077, & + 0.1154, 0.1227, 0.1294, 0.1359, & + 0.1640, 0.1876, 0.2082, 0.2264, & + 0.25808, 0.28498, 0.30844, 0.32932, & + 0.34817, 0.36540, 0.38132 /) + + real, parameter, dimension(nval_tab1a) :: & + b2_tab1a = (/ & + 0.3192, 0.3243, 0.3306, 0.3433, & + 0.3537, 0.3705, 0.4031, 0.4262, & + 0.4456, 0.4621, 0.4763, 0.4889, & + 0.4999, 0.5100, 0.5191, 0.5347, & + 0.5477, 0.5588, 0.5682, 0.5764, & + 0.6042, 0.6206, 0.6324, 0.6425, & + 0.66172, 0.68144, 0.70086, 0.72144, & + 0.74178, 0.76190, 0.78155 /) + + !! Make the table big enough so step size is smaller + !! in log-space that any increment in Table 1a + integer, parameter :: nval_lut=401 + real :: chl, log10chl_lut, w1, w2 + integer :: n,m,mm1,err + + allocate(optics%a1_lut(nval_lut),optics%b1_lut(nval_lut),& + & optics%a2_lut(nval_lut),optics%b2_lut(nval_lut),& + & stat=err) + if ( err /= 0 ) then + call MOM_error(FATAL,"init_ohlmann: Cannot allocate lookup table") + endif + + optics%log10chl_min = log10(chl_tab1a(1)) + optics%log10chl_max = log10(chl_tab1a(nval_tab1a)) + optics%dlog10chl = (optics%log10chl_max - optics%log10chl_min)/(nval_lut-1) + + ! step through the lookup table + m = 2 + do n=1,nval_lut + log10chl_lut = optics%log10chl_min + (n-1)*optics%dlog10chl + chl = 10.0**log10chl_lut + chl = max(chl_tab1a(1),min(chl,chl_tab1a(nval_tab1a))) + + ! find interval in Table 1a (m-1,m] + do while (chl > chl_tab1a(m)) + m = m + 1 + enddo + mm1 = m-1 + + ! interpolation weights + w2 = (chl - chl_tab1a(mm1))/(chl_tab1a(m) - chl_tab1a(mm1)) + w1 = 1. - w2 + + ! fill in the tables + optics%a1_lut(n) = w1*a1_tab1a(mm1) + w2*a1_tab1a(m) + optics%a2_lut(n) = w1*a2_tab1a(mm1) + w2*a2_tab1a(m) + optics%b1_lut(n) = w1*b1_tab1a(mm1) + w2*b1_tab1a(m) + optics%b2_lut(n) = w1*b2_tab1a(mm1) + w2*b2_tab1a(m) + enddo + + return +end subroutine init_ohlmann_table + +!> Get the partion of total solar into bands from Ohlmann lookup table +function lookup_ohlmann_swpen(chl,optics) result(A) + + implicit none + + real, intent(in) :: chl + type(optics_type), intent(in) :: optics + real, dimension(2) :: A + + ! Local variables + + real :: log10chl + integer :: n + + ! Make sure we are in the table + log10chl = max(optics%log10chl_min,min(log10(chl),optics%log10chl_max)) + ! Do a nearest neighbor lookup + n = nint( (log10chl - optics%log10chl_min)/optics%dlog10chl ) + 1 + + A(1) = optics%a1_lut(n) + A(2) = optics%a2_lut(n) + +end function lookup_ohlmann_swpen + +!> Get the opacity (decay scale) from Ohlmann lookup table +function lookup_ohlmann_opacity(chl,optics) result(B) + + implicit none + real, intent(in) :: chl + type(optics_type), intent(in) :: optics + real, dimension(2) :: B + + ! Local variables + real :: log10chl + integer :: n + + ! Make sure we are in the table + log10chl = max(optics%log10chl_min,min(log10(chl),optics%log10chl_max)) + ! Do a nearest neighbor lookup + n = nint( (log10chl - optics%log10chl_min)/optics%dlog10chl ) + 1 + + B(1) = optics%b1_lut(n) + B(2) = optics%b2_lut(n) + + return +end function lookup_ohlmann_opacity subroutine opacity_end(CS, optics) type(opacity_CS) :: CS !< Opacity control structure @@ -1159,7 +1394,11 @@ subroutine opacity_end(CS, optics) if (allocated(optics%max_wavelength_band)) & deallocate(optics%max_wavelength_band) if (allocated(optics%min_wavelength_band)) & - deallocate(optics%min_wavelength_band) + deallocate(optics%min_wavelength_band) + if (allocated(optics%a1_lut)) deallocate(optics%a1_lut) + if (allocated(optics%a2_lut)) deallocate(optics%a2_lut) + if (allocated(optics%b1_lut)) deallocate(optics%b1_lut) + if (allocated(optics%b2_lut)) deallocate(optics%b2_lut) end subroutine opacity_end !> \namespace mom_opacity @@ -1179,4 +1418,7 @@ end subroutine opacity_end !! and sea-ice in a global model, Geophys. Res. Let., 32, L05603, !! doi:10.1029/2004GL020778. +!! Ohlmann, J.C., 2003: Ocean radiant heating in climate models. +!! J. Climate, 16, 1337-1351, 2003. + end module MOM_opacity diff --git a/src/tracer/MOM_hor_bnd_diffusion.F90 b/src/tracer/MOM_hor_bnd_diffusion.F90 index bb5e48fad9..a81a42b428 100644 --- a/src/tracer/MOM_hor_bnd_diffusion.F90 +++ b/src/tracer/MOM_hor_bnd_diffusion.F90 @@ -32,7 +32,7 @@ module MOM_hor_bnd_diffusion public boundary_k_range, hor_bnd_diffusion_end ! Private parameters to avoid doing string comparisons for bottom or top boundary layer -integer, public, parameter :: SURFACE = -1 !< Set a value that corresponds to the surface bopundary +integer, public, parameter :: SURFACE = -1 !< Set a value that corresponds to the surface boundary integer, public, parameter :: BOTTOM = 1 !< Set a value that corresponds to the bottom boundary #include @@ -146,10 +146,11 @@ logical function hor_bnd_diffusion_init(Time, G, GV, US, param_file, diag, diaba call initialize_remapping( CS%remap_CS, string, boundary_extrapolation = boundary_extrap ,& check_reconstruction=.false., check_remapping=.false.) call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg) - call get_param(param_file, mdl, "DEBUG", debug, default=.false., do_not_log=.true.) + call get_param(param_file, mdl, "DEBUG", debug, & + default=.false., debuggingParam=.true., do_not_log=.true.) call get_param(param_file, mdl, "HBD_DEBUG", CS%debug, & "If true, write out verbose debugging data in the HBD module.", & - default=debug) + default=debug, debuggingParam=.true.) id_clock_hbd = cpu_clock_id('(Ocean HBD)', grain=CLOCK_MODULE) @@ -208,7 +209,7 @@ subroutine hor_bnd_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) tracer => Reg%tr(m) if (CS%debug) then - call hchksum(tracer%t, "before HBD "//tracer%name,G%HI) + call hchksum(tracer%t, "before HBD "//tracer%name, G%HI, scale=tracer%conc_scale) endif ! for diagnostics @@ -264,10 +265,10 @@ subroutine hor_bnd_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) endif if (CS%debug) then - call hchksum(tracer%t, "after HBD "//tracer%name,G%HI) + call hchksum(tracer%t, "after HBD "//tracer%name, G%HI, scale=tracer%conc_scale) ! tracer (native grid) integrated tracer amounts before and after HBD - tracer_int_prev = global_mass_integral(h, G, GV, tracer_old) - tracer_int_end = global_mass_integral(h, G, GV, tracer%t) + tracer_int_prev = global_mass_integral(h, G, GV, tracer_old, scale=tracer%conc_scale) + tracer_int_end = global_mass_integral(h, G, GV, tracer%t, scale=tracer%conc_scale) write(mesg,*) 'Total '//tracer%name//' before/after HBD:', tracer_int_prev, tracer_int_end call MOM_mesg(mesg) endif @@ -1213,7 +1214,7 @@ end subroutine hor_bnd_diffusion_end !! !! \subsection section_harmonic_mean Harmonic Mean !! -!! The harmonic mean (HM) betwen h1 and h2 is defined as: +!! The harmonic mean (HM) between h1 and h2 is defined as: !! !! \f[ HM = \frac{2 \times h1 \times h2}{h1 + h2} \f] !!