diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 3c1fe0d151..989d12e378 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -113,9 +113,11 @@ module MOM use MOM_obsolete_diagnostics, only : register_obsolete_diagnostics use MOM_open_boundary, only : ocean_OBC_type, OBC_registry_type use MOM_open_boundary, only : register_temp_salt_segments, update_segment_tracer_reservoirs +use MOM_open_boundary, only : setup_OBC_tracer_reservoirs, fill_temp_salt_segments use MOM_open_boundary, only : open_boundary_register_restarts, remap_OBC_fields use MOM_open_boundary, only : open_boundary_setup_vert, update_OBC_segment_data -use MOM_open_boundary, only : rotate_OBC_config, rotate_OBC_init, write_OBC_info, chksum_OBC_segments +use MOM_open_boundary, only : rotate_OBC_config, rotate_OBC_init, open_boundary_halo_update +use MOM_open_boundary, only : write_OBC_info, chksum_OBC_segments use MOM_porous_barriers, only : porous_widths_layer, porous_widths_interface, porous_barriers_init use MOM_porous_barriers, only : porous_barrier_CS use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_CS @@ -777,7 +779,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS if (CS%VarMix%use_variable_mixing) then call enable_averages(cycle_time, Time_start + real_to_time(US%T_to_s*cycle_time), CS%diag) - call calc_resoln_function(h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, dt) + call calc_resoln_function(h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, CS%OBC, dt) call calc_depth_function(G, CS%VarMix) call disable_averaging(CS%diag) endif @@ -1284,7 +1286,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_tr_adv, & endif endif endif - if (CS%debug_OBCs .and. associated(CS%OBC)) call chksum_OBC_segments(CS%OBC, G, GV, US, 3) + ! if (CS%debug_OBCs .and. associated(CS%OBC)) call chksum_OBC_segments(CS%OBC, G, GV, US, 3) if (CS%do_dynamics .and. CS%split) then !--------------------------- start SPLIT ! This section uses a split time stepping scheme for the dynamic equations, @@ -2067,7 +2069,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS if (.not. skip_diffusion) then if (CS%VarMix%use_variable_mixing) then call pass_var(CS%h, G%Domain) - call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, dt_offline) + call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, CS%OBC, dt_offline) call calc_depth_function(G, CS%VarMix) call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC) endif @@ -2094,7 +2096,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS if (.not. skip_diffusion) then if (CS%VarMix%use_variable_mixing) then call pass_var(CS%h, G%Domain) - call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, dt_offline) + call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, CS%OBC, dt_offline) call calc_depth_function(G, CS%VarMix) call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC) endif @@ -3260,7 +3262,13 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & endif if (associated(OBC_in)) then - call rotate_OBC_init(OBC_in, G, GV, US, param_file, CS%tv, restart_CSp, CS%OBC) + if (use_temperature) then + call pass_var(CS%tv%T, G%Domain) + call pass_var(CS%tv%S, G%Domain) + call fill_temp_salt_segments(G, GV, US, CS%OBC, CS%tv) + endif + + call rotate_OBC_init(OBC_in, G, CS%OBC) call calc_derived_thermo(CS%tv, CS%h, G, GV, US) call update_OBC_segment_data(G, GV, US, CS%OBC, CS%tv, CS%h, Time) endif @@ -3274,7 +3282,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & endif if (use_ice_shelf) & deallocate(frac_shelf_in,mass_shelf_in) - else + else ! The model is being run without grid rotation. This is true of all production runs. if (use_ice_shelf) then call initialize_ice_shelf(param_file, G, Time, ice_shelf_CSp, diag_ptr, Time_init, & dirs%output_directory, calve_ice_shelf_bergs=point_calving) @@ -3343,8 +3351,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & if (ALE_remap_init_conds(CS%ALE_CSp) .and. .not. query_initialized(CS%h,"h",restart_CSp)) then ! This block is controlled by the ALE parameter REMAP_AFTER_INITIALIZATION. - ! \todo This block exists for legacy reasons and we should phase it out of - ! all examples. !### + ! \todo This block exists for legacy reasons and we should phase it out of all examples. !### if (CS%debug) then call uvchksum("Pre ALE adjust init cond [uv]", CS%u, CS%v, G%HI, haloshift=1, unscale=US%L_T_to_m_s) call hchksum(CS%h,"Pre ALE adjust init cond h", G%HI, haloshift=1, unscale=GV%H_to_MKS) @@ -3657,6 +3664,14 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & call register_diags_offline_transport(Time, CS%diag, CS%offline_CSp, GV, US) endif + if (associated(CS%OBC)) then + ! At this point any information related to the tracer reservoirs has either been read from + ! the restart file or has been specified in the segments. Initialize the tracer reservoir + ! values from the segments if they have not been set via the restart file. + call setup_OBC_tracer_reservoirs(G, GV, CS%OBC, restart_CSp) + call open_boundary_halo_update(G, CS%OBC) + endif + call register_obsolete_diagnostics(param_file, CS%diag) if (use_frazil) then diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index cdba3e0ba9..b6bbd0a0ba 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -28,8 +28,8 @@ module MOM_isopycnal_slopes !> Calculate isopycnal slopes, and optionally return other stratification dependent functions such as N^2 !! and dz*S^2*g-prime used, or calculable from factors used, during the calculation. -subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stanley, & - slope_x, slope_y, N2_u, N2_v, dzu, dzv, dzSxN, dzSyN, halo, OBC) +subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stanley, slope_x, slope_y, & + N2_u, N2_v, dzu, dzv, dzSxN, dzSyN, halo, OBC, OBC_N2) type(ocean_grid_type), intent(in) :: 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 @@ -61,6 +61,9 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan !! Eady growth rate at v-points. [Z T-1 ~> m s-1] integer, optional, intent(in) :: halo !< Halo width over which to compute type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. + logical, optional, intent(in) :: OBC_N2 !< If present and true, use interior data + !! to calculate stratification at open boundary + !! condition faces. ! Local variables real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: & @@ -127,6 +130,8 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan logical :: present_N2_u, present_N2_v logical :: local_open_u_BC, local_open_v_BC ! True if u- or v-face OBCs exist anywhere in the global domain. + logical :: OBC_friendly ! If true, open boundary conditions are in use and only interior data should + ! be used to calculate N2 at OBC faces. integer, dimension(2) :: EOSdom_u ! The shifted I-computational domain to use for equation of ! state calculations at u-points. integer, dimension(2) :: EOSdom_v ! The shifted i-computational domain to use for equation of @@ -155,9 +160,11 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan local_open_u_BC = .false. local_open_v_BC = .false. + OBC_friendly = .false. if (present(OBC)) then ; if (associated(OBC)) then local_open_u_BC = OBC%open_u_BCs_exist_globally local_open_v_BC = OBC%open_v_BCs_exist_globally + if (present(OBC_N2)) OBC_friendly = OBC_N2 endif ; endif use_EOS = associated(tv%eqn_of_state) @@ -241,12 +248,12 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan enddo do I=is-1,ie - GxSpV_u(I) = G_Rho0 !This will be changed if both use_EOS and allocated(tv%SpV_avg) are true + GxSpV_u(I) = G_Rho0 ! This will be changed if both use_EOS and allocated(tv%SpV_avg) are true enddo !$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv,h,e, & !$OMP h_neglect,dz_neglect,h_neglect2, & !$OMP present_N2_u,G_Rho0,N2_u,slope_x,dzSxN,EOSdom_u,EOSdom_h1, & - !$OMP local_open_u_BC,dzu,OBC,use_stanley) & + !$OMP local_open_u_BC,dzu,OBC,use_stanley,OBC_friendly) & !$OMP private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u, & !$OMP drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, & !$OMP drho_dT_dT_h,scrap,pres_h,T_h,S_h, & @@ -266,6 +273,22 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan T_u(I) = 0.25*((T(i,j,k) + T(i+1,j,k)) + (T(i,j,k-1) + T(i+1,j,k-1))) S_u(I) = 0.25*((S(i,j,k) + S(i+1,j,k)) + (S(i,j,k-1) + S(i+1,j,k-1))) enddo + if (OBC_friendly) then + do I=is-1,ie + l_seg = OBC%segnum_u(I,j) + if (l_seg /= OBC_NONE) then + if (OBC%segment(l_seg)%direction == OBC_DIRECTION_E) then + pres_u(I) = pres(i,j,K) + T_u(I) = 0.5*(T(i,j,k) + T(i,j,k-1)) + S_u(I) = 0.5*(S(i,j,k) + S(i,j,k-1)) + elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_W) then + pres_u(I) = pres(i+1,j,K) + T_u(I) = 0.5*(T(i+1,j,k) + T(i+1,j,k-1)) + S_u(I) = 0.5*(S(i+1,j,k) + S(i+1,j,k-1)) + endif + endif + enddo + endif call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u, drho_dS_u, & tv%eqn_of_state, EOSdom_u) if (present_N2_u .or. (present(dzSxN))) then @@ -338,8 +361,21 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan ! The expression for drdz above is mathematically equivalent to: ! drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / & ! ((hg2L/haL) + (hg2R/haR)) - ! This is the gradient of density along geopotentials. + ! which is an estimate of the gradient of density across geopotentials. if (present_N2_u) then + if (OBC_friendly) then ; if (OBC%segnum_u(I,j) /= OBC_NONE) then + l_seg = OBC%segnum_u(I,j) + if (OBC%segment(l_seg)%direction == OBC_DIRECTION_E) then + drdz = drdkL / dzaL ! Note that drdz is not used for slopes at OBC faces. + if (use_EOS .and. allocated(tv%SpV_avg)) & + GxSpV_u(I) = GV%g_Earth * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j,k-1)) + elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_W) then + drdz = drdkR / dzaR + if (use_EOS .and. allocated(tv%SpV_avg)) & + GxSpV_u(I) = GV%g_Earth * 0.5 * (tv%SpV_avg(i+1,j,k) + tv%SpV_avg(i+1,j,k-1)) + endif + endif ; endif + N2_u(I,j,K) = GxSpV_u(I) * drdz * G%mask2dCu(I,j) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2] endif @@ -375,6 +411,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan endif slope = slope * max(G%mask2dT(i,j), G%mask2dT(i+1,j)) endif + slope_x(I,j,K) = slope if (present(dzSxN)) & dzSxN(I,j,K) = sqrt( GxSpV_u(I) * max(0., (wtL * ( dzaL * drdkL )) & @@ -391,7 +428,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan !$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv, & !$OMP h,h_neglect,e,dz_neglect, & !$OMP h_neglect2,present_N2_v,G_Rho0,N2_v,slope_y,dzSyN,EOSdom_v, & - !$OMP dzv,local_open_v_BC,OBC,use_stanley) & + !$OMP dzv,local_open_v_BC,OBC,use_stanley,OBC_friendly) & !$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v, & !$OMP drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, & !$OMP drho_dT_dT_h,scrap,pres_h,T_h,S_h, & @@ -411,6 +448,22 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan T_v(i) = 0.25*((T(i,j,k) + T(i,j+1,k)) + (T(i,j,k-1) + T(i,j+1,k-1))) S_v(i) = 0.25*((S(i,j,k) + S(i,j+1,k)) + (S(i,j,k-1) + S(i,j+1,k-1))) enddo + if (OBC_friendly) then + do i=is,ie + l_seg = OBC%segnum_v(i,J) + if (l_seg /= OBC_NONE) then + if (OBC%segment(l_seg)%direction == OBC_DIRECTION_N) then + pres_v(i) = pres(i,j,K) + T_v(i) = 0.5*(T(i,j,k) + T(i,j,k-1)) + S_v(i) = 0.5*(S(i,j,k) + S(i,j,k-1)) + elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_S) then + pres_v(i) = pres(i,j+1,K) + T_v(i) = 0.5*(T(i,j+1,k) + T(i,j+1,k-1)) + S_v(i) = 0.5*(S(i,j+1,k) + S(i,j+1,k-1)) + endif + endif + enddo + endif call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, drho_dS_v, & tv%eqn_of_state, EOSdom_v) @@ -490,8 +543,23 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan ! The expression for drdz above is mathematically equivalent to: ! drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / & ! ((hg2L/haL) + (hg2R/haR)) - ! This is the gradient of density along geopotentials. - if (present_N2_v) N2_v(i,J,K) = GxSpV_v(i) * drdz * G%mask2dCv(i,J) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2] + ! which is an estimate of the gradient of density across geopotentials. + if (present_N2_v) then + if (OBC_friendly) then ; if (OBC%segnum_v(i,J) /= OBC_NONE) then + l_seg = OBC%segnum_v(i,J) + if (OBC%segment(l_seg)%direction == OBC_DIRECTION_N) then + drdz = drdkL / dzaL ! Note that drdz is not used for slopes at OBC faces. + if (use_EOS .and. allocated(tv%SpV_avg)) & + GxSpV_v(i) = GV%g_Earth * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j,k-1)) + elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_S) then + drdz = drdkL / dzaL + if (use_EOS .and. allocated(tv%SpV_avg)) & + GxSpV_v(i) = GV%g_Earth * 0.5 * (tv%SpV_avg(i,j+1,k) + tv%SpV_avg(i,j+1,k-1)) + endif + endif ; endif + + N2_v(i,J,K) = GxSpV_v(i) * drdz * G%mask2dCv(i,J) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2] + endif if (use_EOS) then drdy = ((wtA * drdjA + wtB * drdjB) / (wtA + wtB) - & diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 5c949112cb..75fc9cf193 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -25,7 +25,7 @@ module MOM_open_boundary use MOM_remapping, only : remappingSchemesDoc, remappingDefaultScheme, remapping_CS use MOM_remapping, only : initialize_remapping, remapping_core_h, end_remapping use MOM_restart, only : register_restart_field, register_restart_pair -use MOM_restart, only : query_initialized, MOM_restart_CS +use MOM_restart, only : query_initialized, set_initialized, MOM_restart_CS use MOM_string_functions, only : extract_word, remove_spaces, uppercase, lowercase use MOM_tidal_forcing, only : astro_longitudes, astro_longitudes_init, eq_phase, nodal_fu, tidal_frequency use MOM_time_manager, only : set_date, time_type, time_type_to_real, operator(-) @@ -41,7 +41,7 @@ module MOM_open_boundary public open_boundary_apply_normal_flow public open_boundary_config public open_boundary_setup_vert -public open_boundary_init +public open_boundary_halo_update public open_boundary_query public open_boundary_end public open_boundary_impose_normal_slope @@ -68,6 +68,7 @@ module MOM_open_boundary public setup_OBC_tracer_reservoirs public open_boundary_register_restarts public update_segment_tracer_reservoirs +public set_initialized_OBC_tracer_reservoirs public update_OBC_ramp public remap_OBC_fields public rotate_OBC_config @@ -375,6 +376,8 @@ module MOM_open_boundary real, pointer :: tres_y(:,:,:,:) => Null() !< Array storage of tracer reservoirs for restarts, !! in unscaled units [conc] logical :: debug !< If true, write verbose checksums for debugging purposes. + integer :: nk_OBC_debug = 0 !< The number of layers of OBC segment data to write out + !! in full when DEBUG_OBCS is true. real :: silly_h !< A silly value of thickness outside of the domain that can be used to test !! the independence of the OBCs to this external data [Z ~> m]. real :: silly_u !< A silly value of velocity outside of the domain that can be used to test @@ -398,6 +401,10 @@ module MOM_open_boundary type(group_pass_type) :: pass_oblique !< Structure for group halo pass logical :: exterior_OBC_bug !< If true, use incorrect form of tracers exterior to OBCs. logical :: hor_index_bug !< If true, recover set of a horizontal indexing bugs in the OBC code. + logical :: reservoir_init_bug !< If true, set the OBC tracer reservoirs at the startup of a new + !! run from the interior tracer concentrations regardless of + !! properties that may be explicitly specified for the reservoir + !! concentrations. end type ocean_OBC_type !> Control structure for open boundaries that read from files. @@ -563,6 +570,10 @@ subroutine open_boundary_config(G, US, param_file, OBC) "If true, do additional calls resetting certain values to help verify the correctness "//& "of the open boundary condition code.", & default=.false., old_name="DEBUG_OBC", debuggingParam=.true.) + call get_param(param_file, mdl, "NK_OBC_DEBUG", OBC%nk_OBC_debug, & + "The number of layers of OBC segment data to write out in full "//& + "when DEBUG_OBCS is true.", & + default=0, debuggingParam=.true., do_not_log=.not.OBC%debug) call get_param(param_file, mdl, "OBC_SILLY_THICK", OBC%silly_h, & "A silly value of thicknesses used outside of open boundary "//& @@ -579,6 +590,11 @@ subroutine open_boundary_config(G, US, param_file, OBC) call get_param(param_file, mdl, "OBC_HOR_INDEXING_BUG", OBC%hor_index_bug, & "If true, recover set of a horizontal indexing bugs in the OBC code.", & default=.true.) + call get_param(param_file, mdl, "OBC_RESERVOIR_INIT_BUG", OBC%reservoir_init_bug, & + "If true, set the OBC tracer reservoirs at the startup of a new run from the "//& + "interior tracer concentrations regardless of properties that may be explicitly "//& + "specified for the reservoir concentrations.", default=.true., do_not_log=.true.) + !### Change the default of OBC_RESERVOIR_INIT_BUG to false. reentrant_x = .false. call get_param(param_file, mdl, "REENTRANT_X", reentrant_x, default=.true.) reentrant_y = .false. @@ -1945,21 +1961,13 @@ subroutine parse_for_tracer_reservoirs(OBC, PF, use_temperature) end subroutine parse_for_tracer_reservoirs -!> Initialize open boundary control structure and do any necessary rescaling of OBC -!! fields that have been read from a restart file. -subroutine open_boundary_init(G, GV, US, param_file, OBC, restart_CS) +!> Do any necessary halo updates on OBC-related fields. +subroutine open_boundary_halo_update(G, OBC) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - type(verticalGrid_type), intent(in) :: GV !< Container for vertical grid information - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(param_file_type), intent(in) :: param_file !< Parameter file handle type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure - type(MOM_restart_CS), intent(in) :: restart_CS !< Restart structure, data intent(inout) ! Local variables - integer :: i, j, k, isd, ied, jsd, jed, nz, m - integer :: IsdB, IedB, JsdB, JedB - isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke - IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + integer :: m if (.not.associated(OBC)) return @@ -1989,7 +1997,7 @@ subroutine open_boundary_init(G, GV, US, param_file, OBC, restart_CS) enddo endif -end subroutine open_boundary_init +end subroutine open_boundary_halo_update logical function open_boundary_query(OBC, apply_open_OBC, apply_specified_OBC, apply_Flather_OBC, & apply_nudged_OBC, needs_ext_seg_data) @@ -2232,50 +2240,99 @@ subroutine open_boundary_impose_land_mask(OBC, G, areaCu, areaCv, US) end subroutine open_boundary_impose_land_mask -!> Make sure the OBC tracer reservoirs are initialized. -subroutine setup_OBC_tracer_reservoirs(G, GV, OBC) - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure - type(ocean_OBC_type), target, intent(inout) :: OBC !< Open boundary control structure +!> Initialize the tracer reservoirs values, perhaps only if they have not been set via a restart file. +subroutine setup_OBC_tracer_reservoirs(G, GV, OBC, restart_CS) + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(ocean_OBC_type), target, intent(inout) :: OBC !< Open boundary control structure + type(MOM_restart_CS), optional, intent(in) :: restart_CS !< MOM restart control structure + ! Local variables type(OBC_segment_type), pointer :: segment => NULL() real :: I_scale ! The inverse of the scaling factor for the tracers. ! For salinity the units would be [ppt S-1 ~> 1] + logical :: set_tres_x, set_tres_y + character(len=12) :: x_var_name, y_var_name integer :: i, j, k, m, n - do n=1,OBC%number_of_segments - segment=>OBC%segment(n) - if (associated(segment%tr_Reg)) then - if (segment%is_E_or_W) then - I = segment%HI%IsdB - do m=1,OBC%ntr - I_scale = 1.0 ; if (segment%tr_Reg%Tr(m)%scale /= 0.0) I_scale = 1.0 / segment%tr_Reg%Tr(m)%scale - if (allocated(segment%tr_Reg%Tr(m)%tres)) then - do k=1,GV%ke - do j=segment%HI%jsd,segment%HI%jed - OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%t(i,j,k) - enddo - enddo - endif - enddo + do m=1,OBC%ntr + + set_tres_x = associated(OBC%tres_x) .and. OBC%tracer_x_reservoirs_used(m) + set_tres_y = associated(OBC%tres_y) .and. OBC%tracer_y_reservoirs_used(m) + + if (present(restart_CS)) then + ! Set the names of the reservoirs for this tracer in the restart file, and inquire whether + ! they have been initialized + if (modulo(G%HI%turns, 2) == 0) then + write(x_var_name,'("tres_x_",I3.3)') m + write(y_var_name,'("tres_y_",I3.3)') m else - J = segment%HI%JsdB - do m=1,OBC%ntr - I_scale = 1.0 ; if (segment%tr_Reg%Tr(m)%scale /= 0.0) I_scale = 1.0 / segment%tr_Reg%Tr(m)%scale - if (allocated(segment%tr_Reg%Tr(m)%tres)) then - do k=1,GV%ke - do i=segment%HI%isd,segment%HI%ied - OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%t(i,J,k) - enddo - enddo - endif - enddo + write(x_var_name,'("tres_y_",I3.3)') m + write(y_var_name,'("tres_x_",I3.3)') m endif + if (set_tres_x) set_tres_x = .not.query_initialized(OBC%tres_x, x_var_name, restart_CS) + if (set_tres_y) set_tres_y = .not.query_initialized(OBC%tres_y, y_var_name, restart_CS) endif + + do n=1,OBC%number_of_segments + segment => OBC%segment(n) + if (associated(segment%tr_Reg)) then ; if (allocated(segment%tr_Reg%Tr(m)%tres)) then + I_scale = 1.0 ; if (segment%tr_Reg%Tr(m)%scale /= 0.0) I_scale = 1.0 / segment%tr_Reg%Tr(m)%scale + + if (segment%is_E_or_W .and. set_tres_x) then + I = segment%HI%IsdB + if (segment%tr_Reg%Tr(m)%is_initialized) then + do k=1,GV%ke ; do j=segment%HI%jsd,segment%HI%jed + OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,j,k) + enddo ; enddo + else + do k=1,GV%ke ; do j=segment%HI%jsd,segment%HI%jed + OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%t(i,j,k) + enddo ; enddo + endif + elseif (segment%is_N_or_S .and. set_tres_y) then + J = segment%HI%JsdB + if (segment%tr_Reg%Tr(m)%is_initialized) then + do k=1,GV%ke ; do i=segment%HI%isd,segment%HI%ied + OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,J,k) + enddo ; enddo + else + do k=1,GV%ke ; do i=segment%HI%isd,segment%HI%ied + OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%t(i,J,k) + enddo ; enddo + endif + endif + endif ; endif + enddo enddo end subroutine setup_OBC_tracer_reservoirs +!> Record that the tracer reservoirs have been initialized so that their values are not reset later. +subroutine set_initialized_OBC_tracer_reservoirs(G, OBC, restart_CS) + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(ocean_OBC_type), intent(in) :: OBC !< Open boundary control structure + type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure + character(len=12) :: x_var_name, y_var_name + integer :: i, j, k, m, n + + do m=1,OBC%ntr + ! Set the names of the reservoirs for this tracer in the restart file + if (modulo(G%HI%turns, 2) == 0) then + write(x_var_name,'("tres_x_",I3.3)') m + write(y_var_name,'("tres_y_",I3.3)') m + else + write(x_var_name,'("tres_y_",I3.3)') m + write(y_var_name,'("tres_x_",I3.3)') m + endif + + if (OBC%tracer_x_reservoirs_used(m)) call set_initialized(OBC%tres_x, x_var_name, restart_CS) + if (OBC%tracer_y_reservoirs_used(m)) call set_initialized(OBC%tres_y, y_var_name, restart_CS) + enddo + +end subroutine set_initialized_OBC_tracer_reservoirs + + !> Apply radiation conditions to 3D u,v at open boundaries subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, dt) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure @@ -2326,6 +2383,9 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, if (.not.(OBC%open_u_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally)) & return + if (OBC%debug) call chksum_OBC_segments(OBC, G, GV, US, OBC%nk_OBC_debug) + + eps = 1.0e-20*US%m_s_to_L_T**2 !! Copy previously calculated phase velocity from global arrays into segments @@ -4894,8 +4954,7 @@ subroutine segment_tracer_registry_end(Reg) endif end subroutine segment_tracer_registry_end -!> Register temperature and salinity tracers that are active on an OBC segment. -!! tracer inflow values are specified. +!> Registers the temperature and salinity in the segment tracer registry. subroutine register_temp_salt_segments(GV, US, OBC, tr_Reg, param_file) type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< Unit scaling type @@ -4903,7 +4962,7 @@ subroutine register_temp_salt_segments(GV, US, OBC, tr_Reg, param_file) type(tracer_registry_type), pointer :: tr_Reg !< Tracer registry type(param_file_type), intent(in) :: param_file !< file to parse for model parameter values -! Local variables + ! Local variables integer :: n, ntr_id character(len=32) :: name type(OBC_segment_type), pointer :: segment => NULL() ! pointer to segment type list @@ -4970,12 +5029,13 @@ subroutine get_obgc_segments_props(node, tr_name,obc_src_file_name,obc_src_field node => node%next end subroutine get_obgc_segments_props +!> Registers a named tracer in the segment tracer registries for the OBC segments on which it is active. subroutine register_obgc_segments(GV, OBC, tr_Reg, param_file, tr_name) type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(ocean_OBC_type), pointer :: OBC !< Open boundary structure type(tracer_registry_type), pointer :: tr_Reg !< Tracer registry type(param_file_type), intent(in) :: param_file !< file to parse for model parameter values - character(len=*), intent(in) :: tr_name!< Tracer name + character(len=*), intent(in) :: tr_name !< Tracer name ! Local variables integer :: isd, ied, IsdB, IedB, jsd, jed, JsdB, JedB, nz, nf, ntr_id, fd_id integer :: i, j, k, n, m @@ -4998,6 +5058,7 @@ subroutine register_obgc_segments(GV, OBC, tr_Reg, param_file, tr_name) end subroutine register_obgc_segments +!> Stores the interior tracer values on the segment, and in some cases also sets the tracer reservoir values. subroutine fill_obgc_segments(G, GV, OBC, tr_ptr, tr_name) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure @@ -5017,7 +5078,7 @@ subroutine fill_obgc_segments(G, GV, OBC, tr_ptr, tr_name) do n=1, OBC%number_of_segments segment => OBC%segment(n) if (.not. segment%on_pe) cycle - nt=get_tracer_index(segment,tr_name) + nt = get_tracer_index(segment, tr_name) if (nt < 0) then call MOM_error(FATAL,"fill_obgc_segments: Did not find tracer "// tr_name) endif @@ -5025,40 +5086,66 @@ subroutine fill_obgc_segments(G, GV, OBC, tr_ptr, tr_name) jsd = segment%HI%jsd ; jed = segment%HI%jed IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB - I_scale = 1.0 - if (segment%tr_Reg%Tr(nt)%scale /= 0.0) I_scale = 1.0 / segment%tr_Reg%Tr(nt)%scale - ! Fill with Tracer values - if (segment%is_E_or_W) then - I=segment%HI%IsdB + + ! Fill segments with Tracer values + if (segment%direction == OBC_DIRECTION_W) then + I = segment%HI%IsdB do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed - if (segment%direction == OBC_DIRECTION_W) then - segment%tr_Reg%Tr(nt)%t(I,j,k) = tr_ptr(i+1,j,k) - else - segment%tr_Reg%Tr(nt)%t(I,j,k) = tr_ptr(i,j,k) - endif - OBC%tres_x(I,j,k,nt) = I_scale * segment%tr_Reg%Tr(nt)%t(I,j,k) + segment%tr_Reg%Tr(nt)%t(I,j,k) = tr_ptr(i+1,j,k) enddo ; enddo - else - J=segment%HI%JsdB + elseif (segment%direction == OBC_DIRECTION_E) then + I = segment%HI%IsdB + do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed + segment%tr_Reg%Tr(nt)%t(I,j,k) = tr_ptr(i,j,k) + enddo ; enddo + elseif (segment%direction == OBC_DIRECTION_S) then + J = segment%HI%JsdB do k=1,nz ; do i=segment%HI%isd,segment%HI%ied - if (segment%direction == OBC_DIRECTION_S) then - segment%tr_Reg%Tr(nt)%t(i,J,k) = tr_ptr(i,j+1,k) - else - segment%tr_Reg%Tr(nt)%t(i,J,k) = tr_ptr(i,j,k) - endif - OBC%tres_y(i,J,k,nt) = I_scale * segment%tr_Reg%Tr(nt)%t(i,J,k) + segment%tr_Reg%Tr(nt)%t(i,J,k) = tr_ptr(i,j+1,k) + enddo ; enddo + elseif (segment%direction == OBC_DIRECTION_N) then + J = segment%HI%JsdB + do k=1,nz ; do i=segment%HI%isd,segment%HI%ied + segment%tr_Reg%Tr(nt)%t(i,J,k) = tr_ptr(i,j,k) enddo ; enddo endif - segment%tr_Reg%Tr(nt)%tres(:,:,:) = segment%tr_Reg%Tr(nt)%t(:,:,:) - enddo + + if (.not.segment%tr_Reg%Tr(nt)%is_initialized) & + segment%tr_Reg%Tr(nt)%tres(:,:,:) = segment%tr_Reg%Tr(nt)%t(:,:,:) + + if (OBC%reservoir_init_bug) then + ! OBC%tres_x and OBC%tres_y should not be set here, but in a subsequent call to setup_OBC_tracer_reservoirs. + ! Note that fill_obgc_segments is not called for runs that start from a restart file. + I_scale = 1.0 + if (segment%tr_Reg%Tr(nt)%scale /= 0.0) I_scale = 1.0 / segment%tr_Reg%Tr(nt)%scale + if (segment%is_E_or_W) then + if (associated(OBC%tres_x)) then + I = segment%HI%IsdB + do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed + OBC%tres_x(I,j,k,nt) = I_scale * segment%tr_Reg%Tr(nt)%tres(I,j,k) + enddo ; enddo + endif + else ! segment%is_N_or_S + if (associated(OBC%tres_y)) then + J = segment%HI%JsdB + do k=1,nz ; do i=segment%HI%isd,segment%HI%ied + OBC%tres_y(i,J,k,nt) = I_scale * segment%tr_Reg%Tr(nt)%tres(i,J,k) + enddo ; enddo + endif + endif + endif + + enddo ! End of loop over segments. + end subroutine fill_obgc_segments +!> Set the value of temperatures and salinities on OBC segments subroutine fill_temp_salt_segments(G, GV, US, OBC, tv) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< Unit scaling type(ocean_OBC_type), pointer :: OBC !< Open boundary structure - type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure integer :: isd, ied, IsdB, IedB, jsd, jed, JsdB, JedB, n, nz integer :: i, j, k @@ -5068,9 +5155,6 @@ subroutine fill_temp_salt_segments(G, GV, US, OBC, tv) if (.not. associated(tv%T) .and. associated(tv%S)) return ! Both temperature and salinity fields - call pass_var(tv%T, G%Domain) - call pass_var(tv%S, G%Domain) - nz = GV%ke do n=1, OBC%number_of_segments @@ -5106,11 +5190,12 @@ subroutine fill_temp_salt_segments(G, GV, US, OBC, tv) endif enddo ; enddo endif - segment%tr_Reg%Tr(1)%tres(:,:,:) = segment%tr_Reg%Tr(1)%t(:,:,:) - segment%tr_Reg%Tr(2)%tres(:,:,:) = segment%tr_Reg%Tr(2)%t(:,:,:) + if (.not.segment%tr_Reg%Tr(1)%is_initialized) & + segment%tr_Reg%Tr(1)%tres(:,:,:) = segment%tr_Reg%Tr(1)%t(:,:,:) + if (.not.segment%tr_Reg%Tr(2)%is_initialized) & + segment%tr_Reg%Tr(2)%tres(:,:,:) = segment%tr_Reg%Tr(2)%t(:,:,:) enddo - call setup_OBC_tracer_reservoirs(G, GV, OBC) end subroutine fill_temp_salt_segments !> Find the region outside of all open boundary segments and @@ -6076,9 +6161,12 @@ subroutine rotate_OBC_config(OBC_in, G_in, OBC, G, turns) OBC%remap_h_CS = OBC_in%remap_h_CS endif - ! TODO: The OBC registry seems to be a list of "registered" OBC types. + ! TODO: The OBC registry is a list of "registered" OBC types that is set up + ! on the rotated grid after rotate_OBC_config is called. ! It does not appear to be used, so for now we skip this record. - !OBC%OBC_Reg => OBC_in%OBC_Reg + ! OBC%OBC_Reg => OBC_in%OBC_Reg + ! OBC%num_obgc_tracers = OBC_in%num_obgc_tracers + ! if (associated(OBC_in%OBC_Reg) .and. (.not.associated(OBC%OBC_Reg))) allocate(OBC%OBC_Reg) end subroutine rotate_OBC_config !> Rotate the OBC segment configuration data from the input to model index map. @@ -6279,36 +6367,20 @@ end function rotate_OBC_segment_direction !> Initialize the segments and field-related data of a rotated OBC. -subroutine rotate_OBC_init(OBC_in, G, GV, US, param_file, tv, restart_CS, OBC) +subroutine rotate_OBC_init(OBC_in, G, OBC) type(ocean_OBC_type), intent(in) :: OBC_in !< OBC on input map type(ocean_grid_type), intent(in) :: G !< Rotated grid metric - type(verticalGrid_type), intent(in) :: GV !< Vertical grid - type(unit_scale_type), intent(in) :: US !< Unit scaling - type(param_file_type), intent(in) :: param_file !< Input parameters - type(thermo_var_ptrs), intent(inout) :: tv !< Tracer fields - type(MOM_restart_CS), intent(in) :: restart_CS !< Restart CS type(ocean_OBC_type), pointer, intent(inout) :: OBC !< Rotated OBC - logical :: use_temperature - integer :: l + integer :: l_seg ! update_OBC may have been updated during initialization. OBC%update_OBC = OBC_in%update_OBC - call get_param(param_file, "MOM", "ENABLE_THERMODYNAMICS", use_temperature, & - "If true, Temperature and salinity are used as state "//& - "variables.", default=.true., do_not_log=.true.) - - if (use_temperature) & - call fill_temp_salt_segments(G, GV, US, OBC, tv) - - do l = 1, OBC%number_of_segments - call rotate_OBC_segment_data(OBC_in%segment(l), OBC%segment(l), G%HI%turns) + do l_seg = 1, OBC%number_of_segments + call rotate_OBC_segment_data(OBC_in%segment(l_seg), OBC%segment(l_seg), G%HI%turns) enddo - ! There is already a call to setup_OBC_tracer_reservoirs in fill_temp_salt_segments - call setup_OBC_tracer_reservoirs(G, GV, OBC) - call open_boundary_init(G, GV, US, param_file, OBC, restart_CS) end subroutine rotate_OBC_init @@ -6456,11 +6528,46 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) if (allocated(segment_in%nudged_tangential_grad)) & call rotate_array(segment_in%nudged_tangential_grad, turns, segment%nudged_tangential_grad) if (associated(segment_in%tr_Reg)) then + isd = segment%HI%isd ; ied = segment%HI%ied ; IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB + jsd = segment%HI%jsd ; jed = segment%HI%jed ; JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB + + if (.not.associated(segment%tr_Reg)) allocate(segment%tr_Reg) + segment%tr_Reg%ntseg = segment_in%tr_Reg%ntseg + do n = 1, segment_in%tr_Reg%ntseg - call rotate_array(segment_in%tr_Reg%tr(n)%t, turns, segment%tr_Reg%tr(n)%t) - call rotate_array(segment_in%tr_Reg%tr(n)%tres, turns, segment%tr_Reg%tr(n)%tres) - ! Testing this to see if it works for contant tres values. Probably wrong otherwise. - segment%tr_Reg%Tr(n)%is_initialized=.true. + ! segment_in already points to the rotated tracer fields in the registry. + if (associated(segment_in%tr_Reg%Tr(n)%Tr)) & + segment%tr_Reg%Tr(n)%Tr => segment_in%tr_Reg%Tr(n)%Tr + segment%tr_Reg%Tr(n)%name = segment_in%tr_Reg%Tr(n)%name + segment%tr_Reg%Tr(n)%ntr_index = segment_in%tr_Reg%Tr(n)%ntr_index + segment%tr_Reg%Tr(n)%fd_index = segment_in%tr_Reg%Tr(n)%fd_index + segment%tr_Reg%Tr(n)%scale = segment_in%tr_Reg%Tr(n)%scale + segment%tr_Reg%Tr(n)%OBC_inflow_conc = segment_in%tr_Reg%Tr(n)%OBC_inflow_conc + + if (allocated(segment_in%tr_Reg%tr(n)%t)) then + if (.not.allocated(segment%tr_Reg%tr(n)%t)) then + ke = size(segment_in%tr_Reg%tr(n)%t, 3) + if (segment%is_E_or_W) then + allocate(segment%tr_Reg%Tr(n)%t(IsdB:IedB,jsd:jed,1:ke), source=0.0) + elseif (segment%is_N_or_S) then + allocate(segment%tr_Reg%Tr(n)%t(isd:ied,JsdB:JedB,1:ke), source=0.0) + endif + endif + call rotate_array(segment_in%tr_Reg%tr(n)%t, turns, segment%tr_Reg%tr(n)%t) + endif + + if (allocated(segment_in%tr_Reg%tr(n)%tres)) then + if (.not.allocated(segment%tr_Reg%tr(n)%tres)) then + ke = size(segment_in%tr_Reg%tr(n)%tres, 3) + if (segment%is_E_or_W) then + allocate(segment%tr_Reg%Tr(n)%tres(IsdB:IedB,jsd:jed,1:ke), source=0.0) + elseif (segment%is_N_or_S) then + allocate(segment%tr_Reg%Tr(n)%tres(isd:ied,JsdB:JedB,1:ke), source=0.0) + endif + endif + call rotate_array(segment_in%tr_Reg%tr(n)%tres, turns, segment%tr_Reg%tr(n)%tres) + endif + segment%tr_Reg%Tr(n)%is_initialized = segment_in%tr_Reg%Tr(n)%is_initialized enddo endif diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 2f1b2209af..d8f6dda94d 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -20,13 +20,13 @@ module MOM_state_initialization use MOM_interface_heights, only : find_eta, dz_to_thickness, dz_to_thickness_simple use MOM_interface_heights, only : calc_derived_thermo use MOM_io, only : file_exists, field_size, MOM_read_data, MOM_read_vector, slasher -use MOM_open_boundary, only : ocean_OBC_type, open_boundary_init, set_tracer_data +use MOM_open_boundary, only : ocean_OBC_type, open_boundary_halo_update, set_tracer_data use MOM_open_boundary, only : OBC_NONE use MOM_open_boundary, only : open_boundary_query use MOM_open_boundary, only : set_tracer_data, initialize_segment_data use MOM_open_boundary, only : open_boundary_test_extern_h -use MOM_open_boundary, only : fill_temp_salt_segments -use MOM_open_boundary, only : update_OBC_segment_data +use MOM_open_boundary, only : fill_temp_salt_segments, setup_OBC_tracer_reservoirs +use MOM_open_boundary, only : update_OBC_segment_data, set_initialized_OBC_tracer_reservoirs !use MOM_open_boundary, only : set_3D_OBC_data use MOM_grid_initialize, only : initialize_masks, set_grid_metrics use MOM_restart, only : restore_state, is_new_run, copy_restart_var, copy_restart_vector @@ -164,6 +164,9 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & logical :: new_sim, rotate_index logical :: use_temperature, use_sponge, use_OBC, use_oda_incupd logical :: verify_restart_time + logical :: OBC_reservoir_init_bug ! If true, set the OBC tracer reservoirs at the startup of a new + ! run from the interior tracer concentrations regardless of properties that + ! may be explicitly specified for the reservoir concentrations. logical :: use_EOS ! If true, density is calculated from T & S using an equation of state. logical :: depress_sfc ! If true, remove the mass that would be displaced ! by a large surface pressure by squeezing the column. @@ -437,8 +440,22 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & end select endif endif ! not from_Z_file. - if (use_temperature .and. use_OBC) & - call fill_temp_salt_segments(G, GV, US, OBC, tv) + + if (use_temperature .and. use_OBC) then + ! Log this parameter later with the other OBC parameters. + call get_param(PF, mdl, "OBC_RESERVOIR_INIT_BUG", OBC_reservoir_init_bug, & + "If true, set the OBC tracer reservoirs at the startup of a new run from the "//& + "interior tracer concentrations regardless of properties that may be explicitly "//& + "specified for the reservoir concentrations.", default=.true., do_not_log=.true.) + !### Change the default of OBC_RESERVOIR_INIT_BUG to false. + if (OBC_reservoir_init_bug) then + ! These calls should be moved down to join the OBC code, but doing so changes answers because + ! the temperatures and salinities can change due to the remapping and reading from the restarts. + call pass_var(tv%T, G%Domain, complete=.false.) + call pass_var(tv%S, G%Domain, complete=.true.) + call fill_temp_salt_segments(G, GV, US, OBC, tv) + endif + endif ! Convert thicknesses from geometric distances in depth units to thickness units or mass-per-unit-area. if (new_sim .and. convert) call dz_to_thickness(dz, tv, h, G, GV, US) @@ -487,6 +504,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & if (new_sim .and. debug) & call hchksum(h, "Pre-ALE_regrid: h ", G%HI, haloshift=1, unscale=GV%H_to_MKS) + ! In this call, OBC is only used for the directions of OBCs when setting thicknesses at + ! velocity points. call ALE_regrid_accelerated(ALE_CSp, G, GV, US, h, tv, regrid_iterations, u, v, OBC, tracer_Reg, & dt=dt, initial=.true.) endif @@ -621,11 +640,32 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & end select endif - ! Reads OBC parameters not pertaining to the location of the boundaries - call open_boundary_init(G, GV, US, PF, OBC, restart_CS) - - ! This controls user code for setting open boundary data if (associated(OBC)) then + call get_param(PF, mdl, "OBC_RESERVOIR_INIT_BUG", OBC_reservoir_init_bug, & + "If true, set the OBC tracer reservoirs at the startup of a new run from the "//& + "interior tracer concentrations regardless of properties that may be explicitly "//& + "specified for the reservoir concentrations.", default=.true.) + !### Change the default of OBC_RESERVOIR_INIT_BUG to false. + if (use_temperature) then + if (OBC_reservoir_init_bug) then + if (new_sim) then + ! Set up OBC%trex_x and OBC%tres_y as they have not been read from a restart file. + call setup_OBC_tracer_reservoirs(G, GV, OBC) + ! Ensure that the values of the tracer reservoirs that have just been set will not be revised. + call set_initialized_OBC_tracer_reservoirs(G, OBC, restart_CS) + endif + else + ! Store the updated temperatures and salinities at the open boundaries, noting that they may + ! still be updated by the calls in the next 50 lines, so the code setting the tracer + ! reservoir values will come later in the calling routine. + call fill_temp_salt_segments(G, GV, US, OBC, tv) + endif + endif + + ! This does halo updates for OBC-related fields, but it is not needed? + ! call open_boundary_halo_update(G, OBC) + + ! This allocates the arrays on the segments for open boundary data call initialize_segment_data(G, GV, US, OBC, PF) ! call open_boundary_config(G, US, PF, OBC) @@ -633,6 +673,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & ! Call this during initialization to fill boundary arrays from fixed values call update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) + ! This controls user code for setting open boundary data call get_param(PF, mdl, "OBC_USER_CONFIG", config, & "A string that sets how the user code is invoked to set open boundary data: \n"//& " DOME - specified inflow on northern boundary\n"//& @@ -664,22 +705,24 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & call MOM_error(FATAL, "The open boundary conditions specified by "//& "OBC_USER_CONFIG = "//trim(config)//" have not been fully implemented.") endif + if (open_boundary_query(OBC, apply_open_OBC=.true.)) then call set_tracer_data(OBC, tv, h, G, GV, PF) endif - endif -! if (open_boundary_query(OBC, apply_nudged_OBC=.true.)) then -! call set_3D_OBC_data(OBC, tv, h, G, PF, tracer_Reg) -! endif - ! Still need a way to specify the boundary values - if (debug.and.associated(OBC)) then - call hchksum(G%mask2dT, 'MOM_initialize_state: mask2dT ', G%HI) - call uvchksum('MOM_initialize_state: mask2dC[uv]', G%mask2dCu, & - G%mask2dCv, G%HI) - call qchksum(G%mask2dBu, 'MOM_initialize_state: mask2dBu ', G%HI) + + ! if (open_boundary_query(OBC, apply_nudged_OBC=.true.)) then + ! call set_3D_OBC_data(OBC, tv, h, G, PF, tracer_Reg) + ! endif + ! Still need a way to specify the boundary values + + if (debug) then + call hchksum(G%mask2dT, 'MOM_initialize_state: mask2dT ', G%HI) + call uvchksum('MOM_initialize_state: mask2dC[uv]', G%mask2dCu, G%mask2dCv, G%HI) + call qchksum(G%mask2dBu, 'MOM_initialize_state: mask2dBu ', G%HI) + endif + if (debug_OBC) call open_boundary_test_extern_h(G, GV, OBC, h) endif - if (debug_OBC) call open_boundary_test_extern_h(G, GV, OBC, h) call callTree_leave('MOM_initialize_state()') ! Set-up of data Assimilation with incremental update diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index f2f476b0c8..e9dae96e92 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -17,10 +17,10 @@ module MOM_lateral_mixing_coeffs use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use MOM_wave_speed, only : wave_speed, wave_speed_CS, wave_speed_init -use MOM_open_boundary, only : ocean_OBC_type +use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE +use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S use MOM_MEKE_types, only : MEKE_type - implicit none ; private #include @@ -80,10 +80,13 @@ module MOM_lateral_mixing_coeffs !! in its denominator, rather than just the nominal depth of !! the bathymetry. This only applies when using the model !! interface heights as a proxy for isopycnal slopes. + logical :: OBC_friendly !< If true, use only interior data for thickness weighting and + !! to calculate stratification and other fields at open boundary + !! condition faces. real :: cropping_distance !< Distance from surface or bottom to filter out outcropped or !! incropped interfaces for the Eady growth rate calc [Z ~> m] real :: h_min_N2 !< The minimum vertical distance to use in the denominator of the - !! bouyancy frequency used in the slope calculation [H ~> m or kg m-2] + !! buoyancy frequency used in the slope calculation [H ~> m or kg m-2] real, allocatable :: SN_u(:,:) !< S*N at u-points [T-1 ~> s-1] real, allocatable :: SN_v(:,:) !< S*N at v-points [T-1 ~> s-1] @@ -232,7 +235,7 @@ subroutine calc_depth_function(G, CS) end subroutine calc_depth_function !> Calculates and stores the non-dimensional resolution functions -subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE, dt) +subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE, OBC, dt) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] @@ -240,6 +243,7 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE, dt) type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure type(MEKE_type), intent(in) :: MEKE !< MEKE struct + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure real, intent(in) :: dt !< Time increment [T ~> s] ! Local variables @@ -284,7 +288,7 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE, dt) endif if (CS%BS_use_sqg_struct .or. CS%khth_use_sqg_struct .or. CS%khtr_use_sqg_struct & .or. CS%kdgl90_use_sqg_struct .or. CS%id_sqg_struct>0) then - call calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE) + call calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE, OBC) call pass_var(CS%sqg_struct, G%Domain) endif @@ -536,21 +540,21 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE, dt) end subroutine calc_resoln_function !> Calculates and stores functions of SQG mode -subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE) +subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE, OBC) type(ocean_grid_type), intent(inout) :: 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 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] - type(thermo_var_ptrs), intent(in) :: tv ! s] + type(thermo_var_ptrs), intent(in) :: tv ! s] type(VarMix_CS), intent(inout) :: CS !< Variable mixing control struct - type(MEKE_type), intent(in) :: MEKE !< MEKE struct + type(MEKE_type), intent(in) :: MEKE !< MEKE struct + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure ! Local variables - real, dimension(SZI_(G), SZJ_(G),SZK_(GV)+1) :: & - e ! The interface heights relative to mean sea level [Z ~> m]. - real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: N2_u ! Square of Brunt-Vaisala freq at u-points [L2 Z-2 T-2 ~> s-2] - real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: N2_v ! Square of Brunt-Vaisala freq at v-points [L2 Z-2 T-2 ~> s-2] + real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: e ! The interface heights relative to mean sea level [Z ~> m] + real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: N2_u ! Square of buoyancy frequency at u-points [L2 Z-2 T-2 ~> s-2] + real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: N2_v ! Square of buoyancy frequency at v-points [L2 Z-2 T-2 ~> s-2] real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: dzu ! Z-thickness at u-points [Z ~> m] real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: dzv ! Z-thickness at v-points [Z ~> m] real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: dzSxN ! |Sx| N times dz at u-points [Z T-1 ~> m s-1] @@ -570,8 +574,8 @@ subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE) call find_eta(h, tv, G, GV, US, e, halo_size=2) call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, & - CS%slope_x, CS%slope_y, N2_u=N2_u, N2_v=N2_v,dzu=dzu, dzv=dzv, & - dzSxN=dzSxN, dzSyN=dzSyN, halo=1) + CS%slope_x, CS%slope_y, N2_u=N2_u, N2_v=N2_v, dzu=dzu, dzv=dzv, & + dzSxN=dzSxN, dzSyN=dzSyN, halo=1, OBC=OBC, OBC_N2=CS%OBC_friendly) if (CS%sqg_expo<=0.) then CS%sqg_struct(:,:,:) = 1. @@ -620,10 +624,11 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC) real, intent(in) :: dt !< Time increment [T ~> s] type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure + ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! The interface heights relative to mean sea level [Z ~> m] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: N2_u ! Square of Brunt-Vaisala freq at u-points [L2 Z-2 T-2 ~> s-2] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: N2_v ! Square of Brunt-Vaisala freq at v-points [L2 Z-2 T-2 ~> s-2] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: N2_u ! Square of buoyancy frequency at u-points [L2 Z-2 T-2 ~> s-2] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: N2_v ! Square of buoyancy frequency at v-points [L2 Z-2 T-2 ~> s-2] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: dzu ! Z-thickness at u-points [Z ~> m] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: dzv ! Z-thickness at v-points [Z ~> m] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: dzSxN ! |Sx| N times dz at u-points [Z T-1 ~> m s-1] @@ -637,12 +642,13 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC) if (CS%use_simpler_Eady_growth_rate) then call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, & CS%slope_x, CS%slope_y, N2_u=N2_u, N2_v=N2_v, dzu=dzu, dzv=dzv, & - dzSxN=dzSxN, dzSyN=dzSyN, halo=1, OBC=OBC) + dzSxN=dzSxN, dzSyN=dzSyN, halo=1, OBC=OBC, OBC_N2=CS%OBC_friendly) call calc_Eady_growth_rate_2D(CS, G, GV, US, h, e, dzu, dzv, dzSxN, dzSyN, CS%SN_u, CS%SN_v) elseif (CS%use_stored_slopes) then call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, & - CS%slope_x, CS%slope_y, N2_u=N2_u, N2_v=N2_v, halo=1, OBC=OBC) - call calc_Visbeck_coeffs_old(h, CS%slope_x, CS%slope_y, N2_u, N2_v, G, GV, US, CS) + CS%slope_x, CS%slope_y, N2_u=N2_u, N2_v=N2_v, halo=1, OBC=OBC, & + OBC_N2=CS%OBC_friendly) + call calc_Visbeck_coeffs_old(h, CS%slope_x, CS%slope_y, N2_u, N2_v, G, GV, US, CS, OBC) else call calc_slope_functions_using_just_e(h, G, GV, US, CS, e) endif @@ -666,7 +672,7 @@ end subroutine calc_slope_functions !> Calculates factors used when setting diffusivity coefficients similar to Visbeck et al., 1997. !! This is on older implementation that is susceptible to large values of Eady growth rate !! for incropping layers. -subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, CS) +subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, CS, OBC) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] @@ -679,6 +685,7 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C !! at v-points [L2 Z-2 T-2 ~> s-2] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. ! Local variables real :: S2 ! Interface slope squared [Z2 L-2 ~> nondim] @@ -697,7 +704,14 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C real :: S2_v(SZI_(G),SZJB_(G)) ! At first the thickness-weighted depth integral of the squared ! slope [H Z2 L-2 ~> m or kg m-2] and then the average of the ! squared slope [Z2 L-2 ~> nondim] at v points. - + integer :: OBC_dir_u(SZIB_(G),SZJ_(G)) ! An integer indicating where there are u OBCs: +1 for + ! eastern OBCs, -1 for western OBCs and 0 at points with no OBCs. + integer :: OBC_dir_v(SZI_(G),SZJB_(G)) ! An integer indicating where there are v OBCs: +1 for + ! northern OBCs, -1 for southern OBCs and 0 at points with no OBCs. + real :: h4_u(SZIB_(G),SZJ_(G),SZK_(GV)+1) ! The product of the 4 thicknesses surrounding a u-point + ! interface or the inward equivalent with OBCs [H4 ~> m4 or kg2 m-4] + real :: h4_v(SZI_(G),SZJB_(G),SZK_(GV)+1) ! The product of the 4 thicknesses surrounding a v-point + ! interface or the inward equivalent with OBCs [H4 ~> m4 or kg2 m-4] integer :: i, j, k, is, ie, js, je, nz, l_seg if (.not. CS%initialized) call MOM_error(FATAL, "calc_Visbeck_coeffs_old: "// & @@ -713,11 +727,61 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C S2max = CS%Visbeck_S_max**2 - !$OMP parallel do default(shared) - do j=js-1,je+1 ; do i=is-1,ie+1 - CS%SN_u(i,j) = 0.0 - CS%SN_v(i,j) = 0.0 - enddo ; enddo + CS%SN_u(:,:) = 0.0 + CS%SN_v(:,:) = 0.0 + + ! These settings apply where there are not open boundary conditions. + OBC_dir_u(:,:) = 0 ; OBC_dir_v(:,:) = 0 + + if (associated(OBC).and. CS%OBC_friendly) then + ! Store the direction of any OBC faces. + !$OMP parallel do default(shared) private(l_seg) + do j=js-1,je+1 ; do I=is-1,ie ; if (OBC%segnum_u(I,j) /= OBC_NONE) then + l_seg = OBC%segnum_u(I,j) + if (OBC%segment(l_seg)%direction == OBC_DIRECTION_E) OBC_dir_u(I,j) = 1 + if (OBC%segment(l_seg)%direction == OBC_DIRECTION_W) OBC_dir_u(I,j) = -1 + endif ; enddo ; enddo + !$OMP parallel do default(shared) + do J=js-1,je ; do i=is-1,ie+1 ; if (OBC%segnum_v(i,J) /= OBC_NONE) then + l_seg = OBC%segnum_v(i,J) + if (OBC%segment(l_seg)%direction == OBC_DIRECTION_N) OBC_dir_v(i,J) = 1 + if (OBC%segment(l_seg)%direction == OBC_DIRECTION_S) OBC_dir_v(i,J) = -1 + endif ; enddo ; enddo + + ! Use the masked product of the 4 (or 2) thicknesses around a velocity-point interface for weights. + !$OMP parallel do default(shared) + do K=2,nz + do j=js-1,je+1 ; do I=is-1,ie + if (OBC_dir_u(I,j) == 0) then + h4_u(I,j,K) = G%mask2dCu(I,j) * ( (h(i,j,k)*h(i+1,j,k)) * (h(i,j,k-1)*h(i+1,j,k-1)) ) + elseif (OBC_dir_u(I,j) == 1) then ! OBC_DIRECTION_E + h4_u(I,j,K) = G%mask2dCu(I,j) * ( (h(i,j,k)**2) * (h(i,j,k-1)**2) ) + elseif (OBC_dir_u(I,j) == -1) then ! OBC_DIRECTION_W + h4_u(I,j,K) = G%mask2dCu(I,j) * ( (h(i+1,j,k)**2) * (h(i+1,j,k-1)**2) ) + endif + enddo ; enddo + do J=js-1,je ; do i=is-1,ie+1 + if (OBC_dir_v(i,J) == 0) then + h4_v(i,J,K) = G%mask2dCv(i,J) * ( (h(i,j,k)*h(i,j+1,k)) * (h(i,j,k-1)*h(i,j+1,k-1)) ) + elseif (OBC_dir_v(i,J) == 1) then ! OBC_DIRECTION_N + h4_v(i,J,K) = G%mask2dCv(i,J) * ( (h(i,j,k)**2) * (h(i,j,k-1)**2) ) + elseif (OBC_dir_v(i,J) == -1) then ! OBC_DIRECTION_S + h4_v(i,J,K) = G%mask2dCv(i,J) * ( (h(i,j+1,k)**2) * (h(i,j+1,k-1)**2) ) + endif + enddo ; enddo + enddo + else ! The land mask is sufficient and there are no special considerations taken at OBC points. + ! Use the masked product of the 4 thicknesses around a velocity-point interface for weights. + !$OMP parallel do default(shared) + do K=2,nz + do j=js-1,je+1 ; do I=is-1,ie + h4_u(I,j,K) = G%mask2dCu(I,j) * ( (h(i,j,k)*h(i+1,j,k)) * (h(i,j,k-1)*h(i+1,j,k-1)) ) + enddo ; enddo + do J=js-1,je ; do i=is-1,ie+1 + h4_v(i,J,K) = G%mask2dCv(i,J) * ( (h(i,j,k)*h(i,j+1,k)) * (h(i,j,k-1)*h(i,j+1,k-1)) ) + enddo ; enddo + enddo + endif ! To set the length scale based on the deformation radius, use wave_speed to ! calculate the first-mode gravity wave speed and then blend the equatorial @@ -734,10 +798,17 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C H_geom = sqrt( Hdn * Hup ) !H_geom = H_geom * sqrt(N2) ! WKB-ish !H_geom = H_geom * N2 ! WKB-ish - wSE = G%mask2dCv(i+1,J-1) * ( (h(i+1,j,k)*h(i+1,j-1,k)) * (h(i+1,j,k-1)*h(i+1,j-1,k-1)) ) - wNW = G%mask2dCv(i ,J ) * ( (h(i ,j,k)*h(i ,j+1,k)) * (h(i ,j,k-1)*h(i ,j+1,k-1)) ) - wNE = G%mask2dCv(i+1,J ) * ( (h(i+1,j,k)*h(i+1,j+1,k)) * (h(i+1,j,k-1)*h(i+1,j+1,k-1)) ) - wSW = G%mask2dCv(i ,J-1) * ( (h(i ,j,k)*h(i ,j-1,k)) * (h(i ,j,k-1)*h(i ,j-1,k-1)) ) + wSE = h4_v(i+1,J-1,K) + wNW = h4_v(i,J,K) + wNE = h4_v(i+1,J,K) + wSW = h4_v(i,J-1,K) + if (OBC_dir_u(I,j) == 1) then ! OBC_DIRECTION_E + wSE = 0.0 ; wNE = 0.0 + H_geom = sqrt( h(i,j,k) * h(i,j,k-1) ) + elseif (OBC_dir_u(I,j) == -1) then ! OBC_DIRECTION_W + wSW = 0.0 ; wNW = 0.0 + H_geom = sqrt( h(i+1,j,k) * h(i+1,j,k-1) ) + endif S2 = slope_x(I,j,K)**2 + & (((wNW*slope_y(i,J,K)**2) + (wSE*slope_y(i+1,J-1,K)**2)) + & ((wNE*slope_y(i+1,J,K)**2) + (wSW*slope_y(i,J-1,K)**2)) ) / & @@ -770,10 +841,17 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C H_geom = sqrt( Hdn * Hup ) !H_geom = H_geom * sqrt(N2) ! WKB-ish !H_geom = H_geom * N2 ! WKB-ish - wSE = G%mask2dCu(I,j) * ( (h(i,j ,k)*h(i+1,j ,k)) * (h(i,j ,k-1)*h(i+1,j ,k-1)) ) - wNW = G%mask2dCu(I-1,j+1) * ( (h(i,j+1,k)*h(i-1,j+1,k)) * (h(i,j+1,k-1)*h(i-1,j+1,k-1)) ) - wNE = G%mask2dCu(I,j+1) * ( (h(i,j+1,k)*h(i+1,j+1,k)) * (h(i,j+1,k-1)*h(i+1,j+1,k-1)) ) - wSW = G%mask2dCu(I-1,j) * ( (h(i,j ,k)*h(i-1,j ,k)) * (h(i,j ,k-1)*h(i-1,j ,k-1)) ) + wSE = h4_u(I,j,K) + wNW = h4_u(I-1,j+1,K) + wNE = h4_u(I,j+1,K) + wSW = h4_u(I-1,j,K) + if (OBC_dir_v(i,J) == 1) then ! OBC_DIRECTION_N + wNW = 0.0 ; wNE = 0.0 + H_geom = sqrt( h(i,j,k) * h(i,j,k-1) ) + elseif (OBC_dir_v(i,J) == -1) then ! OBC_DIRECTION_S + wSW = 0.0 ; wSE = 0.0 + H_geom = sqrt( h(i,j+1,k) * h(i,j+1,k-1) ) + endif S2 = slope_y(i,J,K)**2 + & (((wSE*slope_x(I,j,K)**2) + (wNW*slope_x(I-1,j+1,K)**2)) + & ((wNE*slope_x(I,j+1,K)**2) + (wSW*slope_x(I-1,j,K)**2)) ) / & @@ -804,6 +882,8 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C if (CS%debug) then call uvchksum("calc_Visbeck_coeffs_old slope_[xy]", slope_x, slope_y, G%HI, & unscale=US%Z_to_L, haloshift=1) + ! call uvchksum("calc_Visbeck_coeffs_old S2_[uv]", S2_u, S2_v, G%HI, & + ! unscale=US%Z_to_L**2, scalar_pair=.true.) call uvchksum("calc_Visbeck_coeffs_old N2_u, N2_v", N2_u, N2_v, G%HI, & unscale=US%L_to_Z**2*US%s_to_T**2, scalar_pair=.true.) call uvchksum("calc_Visbeck_coeffs_old SN_[uv]", CS%SN_u, CS%SN_v, G%HI, & @@ -1149,7 +1229,7 @@ subroutine calc_QG_slopes(h, tv, dt, G, GV, US, slope_x, slope_y, CS, OBC) call find_eta(h, tv, G, GV, US, e, halo_size=3) call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, & - slope_x, slope_y, halo=2, OBC=OBC) + slope_x, slope_y, halo=2, OBC=OBC, OBC_N2=CS%OBC_friendly) end subroutine calc_QG_slopes @@ -1337,10 +1417,14 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) ! mode wave speed as the starting point for iterations. real :: Stanley_coeff ! Coefficient relating the temperature gradient and sub-gridscale ! temperature variance [nondim] - logical :: om4_remap_via_sub_cells ! Use the OM4-era ramap_via_sub_cells for calculating the EBT structure + logical :: om4_remap_via_sub_cells ! Use the OM4-era remap_via_sub_cells for calculating the EBT structure + logical :: mixing_coefs_OBC_bug ! If false, use only interior data for thickness weighting in + ! lateral mixing coefficient calculations and to calculate stratification + ! and other fields at open boundary condition faces. ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "MOM_lateral_mixing_coeffs" ! This module's name. + integer :: number_of_OBC_segments integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -1471,6 +1555,14 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) if (Stanley_coeff < 0.0) call MOM_error(FATAL, & "STANLEY_COEFF must be set >= 0 if USE_STANLEY_ISO is true.") endif + call get_param(param_file, mdl, "OBC_NUMBER_OF_SEGMENTS", number_of_OBC_segments, & + default=0, do_not_log=.true.) + call get_param(param_file, mdl, "MIXING_COEFS_OBC_BUG", mixing_coefs_OBC_bug, & + "If false, use only interior data for thickness weighting in lateral mixing "//& + "coefficient calculations and to calculate stratification and other fields at "//& + "open boundary condition faces.", default=.true., do_not_log=(number_of_OBC_segments<=0)) + !### Change the default for MIXING_COEFS_OBC_BUG to false. + CS%OBC_friendly = .not. MIXING_COEFS_OBC_BUG if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct & .or. CS%BS_EBT_power>0. .or. CS%khtr_use_ebt_struct) then @@ -2027,7 +2119,7 @@ end subroutine VarMix_end !! \section section_vertical_structure_khth Vertical structure function for KhTh !! !! The thickness diffusivity can be prescribed a vertical distribution with the shape of the equivalent barotropic -!! velocity mode. The structure function is stored in the control structure for thie module (varmix_cs) but is +!! velocity mode. The structure function is stored in the control structure for this module (varmix_cs) but is !! calculated using subroutines in mom_wave_speed. !! !! | Symbol | Module parameter | diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index 53d6b36e4a..bb4d4040d7 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -57,12 +57,12 @@ module MOM_kappa_shear real :: lz_rescale !< A coefficient to rescale the distance to the nearest !! solid boundary. This adjustment is to account for !! regions where 3 dimensional turbulence prevents the - !! growth of shear instabilies [nondim]. + !! growth of shear instabilities [nondim]. real :: TKE_bg !< The background level of TKE [Z2 T-2 ~> m2 s-2]. real :: kappa_0 !< The background diapycnal diffusivity [H Z T-1 ~> m2 s-1 or Pa s] real :: kappa_seed !< A moderately large seed value of diapycnal diffusivity that !! is used as a starting turbulent diffusivity in the iterations - !! to findind an energetically constrained solution for the + !! to finding an energetically constrained solution for the !! shear-driven diffusivity [H Z T-1 ~> m2 s-1 or Pa s] real :: kappa_trunc !< Diffusivities smaller than this are rounded to 0 [H Z T-1 ~> m2 s-1 or Pa s] real :: kappa_tol_err !< The fractional error in kappa that is tolerated [nondim]. @@ -77,7 +77,7 @@ module MOM_kappa_shear !! to estimate the time-averaged diffusivity. logical :: dKdQ_iteration_bug !< If true. use an older, dimensionally inconsistent estimate of !! the derivative of diffusivity with energy in the Newton's method - !! iteration. The bug causes undercorrections when dz > 1m. + !! iteration. The bug causes under-corrections when dz > 1m. logical :: KS_at_vertex !< If true, do the calculations of the shear-driven mixing !! at the cell vertices (i.e., the vorticity points). logical :: eliminate_massless !< If true, massless layers are merged with neighboring @@ -103,6 +103,10 @@ module MOM_kappa_shear !! are some massless layers. logical :: VS_viscosity_bug !< If true, use a bug in the calculation of the viscosity that sets !! it to zero for all vertices that are on a coastline. + logical :: vertex_shear_OBC_bug !< If false, use extra masking when interpolating thicknesses to velocity + !! points for setting up the shear velocities at vertices to avoid using + !! external thicknesses at open boundaries. When OBCs are not in use, + !! this parameter does not change answers, but true is more efficient. logical :: VS_GeometricMean !< If true use geometric averaging for Kd from vertices to tracer points logical :: VS_ThicknessMean !< If true use thickness weighting when averaging Kd from vertices to !! tracer points @@ -214,8 +218,8 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & k0dt = dt*CS%kappa_0 dz_massless = 0.1*sqrt((US%Z_to_m*GV%m_to_H)*k0dt) - if (CS%id_N2_init>0) diag_N2_init(:,:,:) = 0.0 - if (CS%id_S2_init>0) diag_S2_init(:,:,:) = 0.0 + if ((CS%id_N2_init>0) .or. CS%debug) diag_N2_init(:,:,:) = 0.0 + if ((CS%id_S2_init>0) .or. CS%debug) diag_S2_init(:,:,:) = 0.0 if (CS%id_N2_mean>0) diag_N2_mean(:,:,:) = 0.0 if (CS%id_S2_mean>0) diag_S2_mean(:,:,:) = 0.0 @@ -340,10 +344,10 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & if (CS%id_S2_mean>0) then ; do K=1,nz+1 diag_S2_mean(i,j,K) = S2_mean(K) enddo ; endif - if (CS%id_N2_init>0) then ; do K=1,nz+1 + if ((CS%id_N2_init>0) .or. CS%debug) then ; do K=1,nz+1 diag_N2_init(i,j,K) = N2_init(K) enddo ; endif - if (CS%id_S2_init>0) then ; do K=1,nz+1 + if ((CS%id_S2_init>0) .or. CS%debug) then ; do K=1,nz+1 diag_S2_init(i,j,K) = S2_init(K) enddo ; endif else @@ -360,16 +364,16 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & if (kf(K) == 0.0) then if (CS%id_N2_mean>0) diag_N2_mean(i,j,K) = N2_mean(kc(K)) if (CS%id_S2_mean>0) diag_S2_mean(i,j,K) = S2_mean(kc(K)) - if (CS%id_N2_init>0) diag_N2_init(i,j,K) = N2_init(kc(K)) - if (CS%id_S2_init>0) diag_S2_init(i,j,K) = S2_init(kc(K)) + if ((CS%id_N2_init>0) .or. CS%debug) diag_N2_init(i,j,K) = N2_init(kc(K)) + if ((CS%id_S2_init>0) .or. CS%debug) diag_S2_init(i,j,K) = S2_init(kc(K)) else if (CS%id_N2_mean>0) & diag_N2_mean(i,j,K) = (1.0-kf(K)) * N2_mean(kc(K)) + kf(K) * N2_mean(kc(K)+1) if (CS%id_S2_mean>0) & diag_S2_mean(i,j,K) = (1.0-kf(K)) * S2_mean(kc(K)) + kf(K) * S2_mean(kc(K)+1) - if (CS%id_N2_init>0) & + if ((CS%id_N2_init>0) .or. CS%debug) & diag_N2_init(i,j,K) = (1.0-kf(K)) * N2_init(kc(K)) + kf(K) * N2_init(kc(K)+1) - if (CS%id_S2_init>0) & + if ((CS%id_S2_init>0) .or. CS%debug) & diag_S2_init(i,j,K) = (1.0-kf(K)) * S2_init(kc(K)) + kf(K) * S2_init(kc(K)+1) endif enddo @@ -391,6 +395,8 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & enddo ! end of j-loop if (CS%debug) then + call hchksum(diag_N2_init, "kappa_shear N2_init", G%HI, unscale=US%s_to_T**2) + call hchksum(diag_S2_init, "kappa_shear S2_init", G%HI, unscale=US%s_to_T**2) call hchksum(kappa_io, "kappa", G%HI, unscale=GV%HZ_T_to_m2_s) call hchksum(tke_io, "tke", G%HI, unscale=US%Z_to_m**2*US%s_to_T**2) endif @@ -453,8 +459,12 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ dz_3d ! Vertical distance between interface heights [Z ~> m]. real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1) :: & kappa_vertex ! Diffusivity at interfaces and vertices [H Z T-1 ~> m2 s-1 or Pa s] - real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1) :: & - h_vert ! Thicknesses interpolated to vertices [H Z T-1 ~> m2 s-1 or Pa s] + real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)) :: & + h_vert ! Thicknesses interpolated to vertices [H ~> m or kg m-2] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: & + h_at_u ! A mask-weighted thickness interpolated to u-points [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: & + h_at_v ! A mask-weighted thickness interpolated to v-points [H ~> m or kg m-2] real, dimension(SZIB_(G),SZK_(GV)) :: & h_2d, & ! A 2-D version of h interpolated to vertices [H ~> m or kg m-2]. dz_2d, & ! Vertical distance between interface heights [Z ~> m]. @@ -500,16 +510,18 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ real, dimension(SZK_(GV)+1) :: kf ! The fractional weight of interface kc+1 for ! interpolating back to the original index space [nondim]. real :: h_SW, h_SE, h_NW, h_NE ! Thicknesses at adjacent vertices [H ~> m or kg m-2] - real :: mks_to_HZ_T ! A factor used to restore dimensional scaling after the geomentric mean + real :: mks_to_HZ_T ! A factor used to restore dimensional scaling after the geometric mean ! diffusivity is taken using thickness weighted powers [H Z s m-2 T-1 ~> 1] ! or [H Z m s kg-1 T-1 ~> 1] + real :: H_tiny ! A sub-roundoff thickness to use in the denominator when calculating + ! thickness-weighted averages [H ~> m or kg m-2] integer :: IsB, IeB, JsB, JeB, i, j, k, nz, nzc ! Diagnostics that should be deleted? isB = G%isc-1 ; ieB = G%iecB ; jsB = G%jsc-1 ; jeB = G%jecB ; nz = GV%ke - if (CS%id_N2_init>0) diag_N2_init(:,:,:) = 0.0 - if (CS%id_S2_init>0) diag_S2_init(:,:,:) = 0.0 + if ((CS%id_N2_init>0) .or. CS%debug) diag_N2_init(:,:,:) = 0.0 + if ((CS%id_S2_init>0) .or. CS%debug) diag_S2_init(:,:,:) = 0.0 if (CS%id_N2_mean>0) diag_N2_mean(:,:,:) = 0.0 if (CS%id_S2_mean>0) diag_S2_mean(:,:,:) = 0.0 kappa_vertex(:,:,:) = 0.0 @@ -519,10 +531,39 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ k0dt = dt*CS%kappa_0 dz_massless = 0.1*sqrt((US%Z_to_m*GV%m_to_H)*k0dt) I_Prandtl = 0.0 ; if (CS%Prandtl_turb > 0.0) I_Prandtl = 1.0 / CS%Prandtl_turb + H_tiny = 0.5 * GV%H_subroundoff ! Convert layer thicknesses into geometric thickness in height units. call thickness_to_dz(h, tv, dz_3d, G, GV, US, halo_size=1) + if (CS%vertex_shear_OBC_bug) then + !$OMP parallel do default(shared) + do k=1,nz + do j=JsB,JeB+1 ; do I=IsB,IeB + h_at_u(I,j,k) = G%mask2dCu(I,j) * (h(i,j,k) + h(i+1,j,k)) * 0.5 + enddo ; enddo + do J=JsB,JeB ; do i=IsB,IeB+1 + h_at_v(i,J,k) = G%mask2dCv(i,J) * (h(i,j,k) + h(i,j+1,k)) * 0.5 + enddo ; enddo + enddo + else + ! Because G%mask2dCu(I,j) is zero if either G%mask2dT(i,j) or G%mask2dT(i+1,j) except at OBC + ! faces, the following form give equivalent answers to those above unless OBCs are in use, + ! although the former is clearly less complicated and costly. + !$OMP parallel do default(shared) + do k=1,nz + do j=JsB,JeB+1 ; do I=IsB,IeB + h_at_u(I,j,k) = G%mask2dCu(I,j) * (G%mask2dT(i,j) * h(i,j,k) + G%mask2dT(i+1,j) * h(i+1,j,k)) / & + (G%mask2dT(i,j) + G%mask2dT(i+1,j) + 1.0e-36) + enddo ; enddo + do J=JsB,JeB ; do i=IsB,IeB+1 + h_at_v(i,J,k) = G%mask2dCv(i,J) * (G%mask2dT(i,j) * h(i,j,k) + G%mask2dT(i,j+1) * h(i,j+1,k)) / & + (G%mask2dT(i,j) + G%mask2dT(i,j+1) + 1.0e-36) + enddo ; enddo + enddo + endif + + !$OMP parallel do default(private) shared(jsB,jeB,isB,ieB,nz,h,u_in,v_in,use_temperature,tv,G,GV,US,CS,kappa_io, & !$OMP dz_massless,k0dt,p_surf,dt,tke_io,kv_io,kappa_vertex,h_vert,I_Prandtl, & !$OMP diag_N2_init,diag_S2_init,diag_N2_mean,diag_S2_mean) @@ -530,14 +571,11 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ ! Interpolate the various quantities to the corners, using masks. do k=1,nz ; do I=IsB,IeB - u_2d(I,k) = (G%mask2dCu(I,j) * (u_in(I,j,k) * (h(i,j,k) + h(i+1,j,k))) + & - G%mask2dCu(I,j+1) * (u_in(I,j+1,k) * (h(i,j+1,k) + h(i+1,j+1,k))) ) / & - ((G%mask2dCu(I,j) * (h(i,j,k) + h(i+1,j,k)) + & - G%mask2dCu(I,j+1) * (h(i,j+1,k) + h(i+1,j+1,k))) + GV%H_subroundoff) - v_2d(I,k) = (G%mask2dCv(i,J) * (v_in(i,J,k) * (h(i,j,k) + h(i,j+1,k))) + & - G%mask2dCv(i+1,J) * (v_in(i+1,J,k) * (h(i+1,j,k) + h(i+1,j+1,k))) ) / & - ((G%mask2dCv(i,J) * (h(i,j,k) + h(i,j+1,k)) + & - G%mask2dCv(i+1,J) * (h(i+1,j,k) + h(i+1,j+1,k))) + GV%H_subroundoff) + u_2d(I,k) = ( (u_in(I,j,k) * h_at_u(I,j,k)) + (u_in(I,j+1,k) * h_at_u(I,j+1,k)) ) / & + ( (h_at_u(I,j,k) + h_at_u(I,j+1,k)) + H_tiny ) + v_2d(I,k) = ( (v_in(i,J,k) * h_at_v(i,J,k)) + (v_in(i+1,J,k) * h_at_v(i+1,J,k)) ) / & + ( (h_at_v(i,J,k) + h_at_v(i+1,J,k)) + H_tiny ) + I_hwt = 1.0 / (((G%mask2dT(i,j) * h(i,j,k) + G%mask2dT(i+1,j+1) * h(i+1,j+1,k)) + & (G%mask2dT(i+1,j) * h(i+1,j,k) + G%mask2dT(i,j+1) * h(i,j+1,k))) + & GV%H_subroundoff) @@ -668,22 +706,22 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ do K=1,nz+1 kappa_2d(I,K) = kappa_avg(K) if (CS%all_layer_TKE_bug) then - tke_2d(i,K) = tke(K) + tke_2d(I,K) = tke(K) else - tke_2d(i,K) = tke_avg(K) + tke_2d(I,K) = tke_avg(K) endif enddo if (CS%id_N2_mean>0) then ; do K=1,nz+1 - diag_N2_mean(i,j,K) = N2_mean(K) + diag_N2_mean(I,J,K) = N2_mean(K) enddo ; endif if (CS%id_S2_mean>0) then ; do K=1,nz+1 - diag_S2_mean(i,j,K) = S2_mean(K) + diag_S2_mean(I,J,K) = S2_mean(K) enddo ; endif - if (CS%id_N2_init>0) then ; do K=1,nz+1 - diag_N2_init(i,j,K) = N2_init(K) + if ((CS%id_N2_init>0) .or. CS%debug) then ; do K=1,nz+1 + diag_N2_init(I,J,K) = N2_init(K) enddo ; endif - if (CS%id_S2_init>0) then ; do K=1,nz+1 - diag_S2_init(i,j,K) = S2_init(K) + if ((CS%id_S2_init>0) .or. CS%debug) then ; do K=1,nz+1 + diag_S2_init(I,J,K) = S2_init(K) enddo ; endif else do K=1,nz+1 @@ -699,16 +737,16 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ if (kf(K) == 0.0) then if (CS%id_N2_mean>0) diag_N2_mean(I,J,K) = N2_mean(kc(K)) if (CS%id_S2_mean>0) diag_S2_mean(I,J,K) = S2_mean(kc(K)) - if (CS%id_N2_init>0) diag_N2_init(I,J,K) = N2_init(kc(K)) - if (CS%id_S2_init>0) diag_S2_init(I,J,K) = S2_init(kc(K)) + if ((CS%id_N2_init>0) .or. CS%debug) diag_N2_init(I,J,K) = N2_init(kc(K)) + if ((CS%id_S2_init>0) .or. CS%debug) diag_S2_init(I,J,K) = S2_init(kc(K)) else if (CS%id_N2_mean>0) & diag_N2_mean(I,J,K) = (1.0-kf(K)) * N2_mean(kc(K)) + kf(K) * N2_mean(kc(K)+1) if (CS%id_S2_mean>0) & diag_S2_mean(I,J,K) = (1.0-kf(K)) * S2_mean(kc(K)) + kf(K) * S2_mean(kc(K)+1) - if (CS%id_N2_init>0) & + if ((CS%id_N2_init>0) .or. CS%debug) & diag_N2_init(I,J,K) = (1.0-kf(K)) * N2_init(kc(K)) + kf(K) * N2_init(kc(K)+1) - if (CS%id_S2_init>0) & + if ((CS%id_S2_init>0) .or. CS%debug) & diag_S2_init(I,J,K) = (1.0-kf(K)) * S2_init(kc(K)) + kf(K) * S2_init(kc(K)+1) endif enddo @@ -749,40 +787,45 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ kappa_io(i,j,1) = 0.0 kappa_io(i,j,nz+1) = 0.0 enddo ; enddo - if (CS%VS_ThicknessMean) then - ! This conversion factor is required to allow for aribtrary fracional powers of the diffusivities. - if (CS%VS_GeometricMean) mks_to_HZ_T = 1.0 / GV%HZ_T_to_MKS + if (CS%VS_ThicknessMean .and. CS%VS_GeometricMean) then + ! This conversion factor is required to allow for arbitrary fractional powers of the diffusivities. + mks_to_HZ_T = 1.0 / GV%HZ_T_to_MKS !$OMP parallel do default(private) shared(nz,G,GV,CS,kappa_io,kappa_vertex,h_vert) do K=2,nz ; do j=G%jsc,G%jec ; do i=G%isc,G%iec h_SW = 0.5 * (h_vert(I-1,J-1,k) + h_vert(I-1,J-1,k-1)) h_NE = 0.5 * (h_vert(I,J,k) + h_vert(I,J,k-1)) h_NW = 0.5 * (h_vert(I-1,J,k) + h_vert(I-1,J,k-1)) h_SE = 0.5 * (h_vert(I,J-1,k) + h_vert(I,J-1,k-1)) - if (CS%VS_GeometricMean) then - if ((h_SW + h_NE) + (h_NW + h_SE) > 0.0) then - ! The geometric mean is zero if any component is zero, hence the need to use a floor - ! on the value of kappa_trunc in regions on boundaries of shear zones. - I_htot = 1.0 / ((h_SW + h_NE) + (h_NW + h_SE)) - kappa_io(i,j,K) = G%mask2dT(i,j) * mks_to_HZ_T * & - ( ((GV%HZ_T_to_MKS * max(kappa_vertex(I-1,J-1,K),CS%VS_GeoMean_Kdmin))**(h_SW*I_htot) * & - (GV%HZ_T_to_MKS * max(kappa_vertex(I,J,K),CS%VS_GeoMean_Kdmin))**(h_NE*I_htot)) * & - ((GV%HZ_T_to_MKS * max(kappa_vertex(I-1,J,K),CS%VS_GeoMean_Kdmin))**(h_NW*I_htot) * & - (GV%HZ_T_to_MKS * max(kappa_vertex(I,J-1,K),CS%VS_GeoMean_Kdmin))**(h_SE*I_htot)) ) - else - ! If all points have zero thickness, the thikncess-weighted geometric mean is undefined, so use - ! the non-thickness weighted geometric mean instead. - kappa_io(i,j,K) = G%mask2dT(i,j) * sqrt(sqrt( & - (max(kappa_vertex(I-1,J-1,K),CS%VS_GeoMean_Kdmin) * max(kappa_vertex(I,J,K),CS%VS_GeoMean_Kdmin)) * & - (max(kappa_vertex(I-1,J,K),CS%VS_GeoMean_Kdmin) * max(kappa_vertex(I,J-1,K),CS%VS_GeoMean_Kdmin)) )) - endif + if ((h_SW + h_NE) + (h_NW + h_SE) > 0.0) then + ! The geometric mean is zero if any component is zero, hence the need to use a floor + ! on the value of kappa_trunc in regions on boundaries of shear zones. + I_htot = 1.0 / ((h_SW + h_NE) + (h_NW + h_SE)) + kappa_io(i,j,K) = G%mask2dT(i,j) * mks_to_HZ_T * & + ( ((GV%HZ_T_to_MKS * max(kappa_vertex(I-1,J-1,K), CS%VS_GeoMean_Kdmin))**(h_SW*I_htot) * & + (GV%HZ_T_to_MKS * max(kappa_vertex(I,J,K), CS%VS_GeoMean_Kdmin))**(h_NE*I_htot)) * & + ((GV%HZ_T_to_MKS * max(kappa_vertex(I-1,J,K), CS%VS_GeoMean_Kdmin))**(h_NW*I_htot) * & + (GV%HZ_T_to_MKS * max(kappa_vertex(I,J-1,K), CS%VS_GeoMean_Kdmin))**(h_SE*I_htot)) ) else - ! The following expression is a thickness weighted arithmetic mean at tracer points: - I_htot = 1.0 / (((h_SW + h_NE) + (h_NW + h_SE)) + GV%H_subroundoff) - kappa_io(i,j,K) = G%mask2dT(i,j) * & - (((kappa_vertex(I-1,J-1,K)*h_SW) + (kappa_vertex(I,J,K)*h_NE)) + & - ((kappa_vertex(I-1,J,K)*h_NW) + (kappa_vertex(I,J-1,K)*h_SE))) * I_htot + ! If all points have zero thickness, the thickness-weighted geometric mean is undefined, so use + ! the non-thickness weighted geometric mean instead. + kappa_io(i,j,K) = G%mask2dT(i,j) * sqrt(sqrt( & + (max(kappa_vertex(I-1,J-1,K),CS%VS_GeoMean_Kdmin) * max(kappa_vertex(I,J,K),CS%VS_GeoMean_Kdmin)) * & + (max(kappa_vertex(I-1,J,K),CS%VS_GeoMean_Kdmin) * max(kappa_vertex(I,J-1,K),CS%VS_GeoMean_Kdmin)) )) endif enddo ; enddo ; enddo + elseif (CS%VS_ThicknessMean) then ! Use thickness-weighted arithmetic mean diffusivities. + !$OMP parallel do default(private) shared(nz,G,GV,CS,kappa_io,kappa_vertex,h_vert) + do K=2,nz ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + h_SW = 0.5 * (h_vert(I-1,J-1,k) + h_vert(I-1,J-1,k-1)) + h_NE = 0.5 * (h_vert(I,J,k) + h_vert(I,J,k-1)) + h_NW = 0.5 * (h_vert(I-1,J,k) + h_vert(I-1,J,k-1)) + h_SE = 0.5 * (h_vert(I,J-1,k) + h_vert(I,J-1,k-1)) + ! The following expression is a thickness weighted arithmetic mean at tracer points: + I_htot = 1.0 / (((h_SW + h_NE) + (h_NW + h_SE)) + GV%H_subroundoff) + kappa_io(i,j,K) = G%mask2dT(i,j) * & + (((kappa_vertex(I-1,J-1,K)*h_SW) + (kappa_vertex(I,J,K)*h_NE)) + & + ((kappa_vertex(I-1,J,K)*h_NW) + (kappa_vertex(I,J-1,K)*h_SE))) * I_htot + enddo ; enddo ; enddo elseif (CS%VS_GeometricMean) then ! The geometic mean diffusivities are not thickness weighted. !$OMP parallel do default(private) shared(nz,G,CS,kappa_io,kappa_vertex) do K=2,nz ; do j=G%jsc,G%jec ; do i=G%isc,G%iec @@ -800,6 +843,8 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ endif if (CS%debug) then + call Bchksum(diag_N2_init, "shear_vertex N2_init", G%HI, unscale=US%s_to_T**2) + call Bchksum(diag_S2_init, "shear_vertex S2_init", G%HI, unscale=US%s_to_T**2) call hchksum(kappa_io, "kappa", G%HI, unscale=GV%HZ_T_to_m2_s) call Bchksum(tke_io, "tke", G%HI, unscale=US%Z_to_m**2*US%s_to_T**2) endif @@ -1341,7 +1386,7 @@ subroutine calculate_projected_state(kappa, u0, v0, T0, S0, dt, nz, dz, I_dz_int real, dimension(nz), intent(in) :: S0 !< The initial salinity [S ~> ppt]. real, intent(in) :: dt !< The time step [T ~> s]. real, dimension(nz), intent(in) :: dz !< The layer thicknesses [H ~> m or kg m-2] - real, dimension(nz+1), intent(in) :: I_dz_int !< The inverse of the distance between succesive + real, dimension(nz+1), intent(in) :: I_dz_int !< The inverse of the distance between successive !! layer centers [Z-1 ~> m-1]. real, dimension(nz+1), intent(in) :: dbuoy_dT !< The partial derivative of buoyancy with !! temperature [Z T-2 C-1 ~> m s-2 degC-1]. @@ -2040,6 +2085,7 @@ function kappa_shear_init(Time, G, GV, US, param_file, diag, CS) ! for setting the default of KD_SMOOTH [Z2 T-1 ~> m2 s-1] real :: kappa_0_default ! The default value for KD_KAPPA_SHEAR_0 [Z2 T-1 ~> m2 s-1] logical :: merge_mixedlayer + integer :: number_of_OBC_segments logical :: debug_shear logical :: just_read ! If true, this module is not used, so only read the parameters. ! This include declares and sets the variable "version". @@ -2079,6 +2125,16 @@ function kappa_shear_init(Time, G, GV, US, param_file, diag, CS) "If true, use a bug in vertex shear that zeros out viscosities at "//& "vertices on coastlines.", & default=.true., do_not_log=just_read.or.(.not.CS%KS_at_vertex)) + call get_param(param_file, mdl, "OBC_NUMBER_OF_SEGMENTS", number_of_OBC_segments, & + default=0, do_not_log=.true.) + call get_param(param_file, mdl, "VERTEX_SHEAR_OBC_BUG", CS%vertex_shear_OBC_bug, & + "If false, use extra masking when interpolating thicknesses to velocity "//& + "points for setting up the shear velocities at vertices to avoid using "//& + "external thicknesses at open boundaries. When OBCs are not in use, "//& + "this parameter does not change answers, but true is more efficient.", & + default=.true., & + do_not_log=just_read.or.(.not.CS%KS_at_vertex).or.(number_of_OBC_segments<=0)) + !### Use OBC settings to set the default for VERTEX_SHEAR_OBC_BUG? call get_param(param_file, mdl, "VERTEX_SHEAR_GEOMETRIC_MEAN", CS%VS_GeometricMean, & "If true, use a geometric mean for moving diffusivity from "//& "vertices to tracer points. False uses algebraic mean.", & @@ -2092,7 +2148,7 @@ function kappa_shear_init(Time, G, GV, US, param_file, diag, CS) CS%VS_GeoMean_Kdmin, "If using the geometric mean in vertex shear, "//& "use this minimum value for Kd. This is an ad-hoc parameter, the "//& "diffusivities on the edge of shear regions are sensitive to the choice.",& - units="m2 s-1",default=0.0, scale=GV%m2_s_to_HZ_T, do_not_log=just_read) + units="m2 s-1", default=0.0, scale=GV%m2_s_to_HZ_T, do_not_log=just_read) endif call get_param(param_file, mdl, "RINO_CRIT", CS%RiNo_crit, & "The critical Richardson number for shear mixing.", &