From ca0cdb36895b4c96f3d1db90941b2d11a89c3dbc Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 19 Sep 2022 18:28:00 -0400 Subject: [PATCH 1/4] *+Corrected some oblique OBC restarts Added separate restart variables for some of the oblique OBC restart variables at north-south or east-west faces. Before this fix, some of the full 3-d arrays for restarts were being overwritten for oblique OBC segments that join at adjacent corners on the north-east side of tracer points, as is noted in github.com/NOAA-GFDL/MOM6/issues/208. The discussion surrounding this issue confirms that there are cases using oblique OBCs (like some North-West Atlantic cases) that do not past the restart reproducibility tests. Some halo updates were also corrected, in accordance with the introduction of these new variables and their proper staggering locations. In a number of places, the proper case-sensitive horizontal indexing convention, as described in github.com/NOAA-GFDL/MOM6/wiki/Code-style-guide#soft-case-convention, is now being used. This commit also corrected the comments describing a number of the variables in the OBC_segment_type to make it clearer where they are discretized. This changes the names of oblique OBC-related variables in the restart files, but since the previous version did not reproduce across restarts, there does not seem to be any point in retaining those incorrect answers. All answers in the MOM6-examples test suite are bitwise identical, but this will change (fix) solutions with oblique OBCs and OBC_RAD_VEL_WT < 1. --- src/core/MOM_open_boundary.F90 | 242 +++++++++++++++++---------------- 1 file changed, 123 insertions(+), 119 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index cabf9272ca..a187ac5b2f 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -116,7 +116,8 @@ module MOM_open_boundary !! Not sure who should lock it or when... end type segment_tracer_registry_type -!> Open boundary segment data structure. +!> Open boundary segment data structure. Unless otherwise noted, 2-d and 3-d arrays are discretized +!! at the same position as normal velocity points in the middle of the OBC segments. type, public :: OBC_segment_type logical :: Flather !< If true, applies Flather + Chapman radiation of barotropic gravity waves. logical :: radiation !< If true, 1D Orlanksi radiation boundary conditions are applied. @@ -178,10 +179,10 @@ module MOM_open_boundary real, allocatable :: h(:,:,:) !< The cell thickness [H ~> m or kg m-2] at OBC-points. real, allocatable :: normal_vel(:,:,:) !< The layer velocity normal to the OB !! segment [L T-1 ~> m s-1]. - real, allocatable :: tangential_vel(:,:,:) !< The layer velocity tangential to the - !! OB segment [L T-1 ~> m s-1]. - real, allocatable :: tangential_grad(:,:,:) !< The gradient of the velocity tangential - !! to the OB segment [T-1 ~> s-1]. + real, allocatable :: tangential_vel(:,:,:) !< The layer velocity tangential to the OB segment + !! [L T-1 ~> m s-1], discretized at the corner points. + real, allocatable :: tangential_grad(:,:,:) !< The gradient of the velocity tangential to the OB + !! segment [T-1 ~> s-1], discretized at the corner points. real, allocatable :: normal_trans(:,:,:) !< The layer transport normal to the OB !! segment [H L2 T-1 ~> m3 s-1]. real, allocatable :: normal_vel_bt(:,:) !< The barotropic velocity normal to @@ -189,25 +190,38 @@ module MOM_open_boundary real, allocatable :: eta(:,:) !< The sea-surface elevation along the !! segment [H ~> m or kg m-2]. real, allocatable :: grad_normal(:,:,:) !< The gradient of the normal flow along the - !! segment times the grid spacing [L T-1 ~> m s-1] + !! segment times the grid spacing [L T-1 ~> m s-1], + !! with the first index being the corner-point index + !! along the segment, and the second index being 1 (for + !! values one point into the domain) or 2 (for values + !! along the OBC itself) real, allocatable :: grad_tan(:,:,:) !< The gradient of the tangential flow along the - !! segment times the grid spacing [L T-1 ~> m s-1] - real, allocatable :: grad_gradient(:,:,:) !< The gradient of the gradient of tangential flow along - !! the segment times the grid spacing [T-1 ~> s-1] + !! segment times the grid spacing [L T-1 ~> m s-1], with the + !! first index being the velocity/tracer point index along the + !! segment, and the second being 1 for the value 1.5 points + !! inside the domain and 2 for the value half a point + !! inside the domain. + real, allocatable :: grad_gradient(:,:,:) !< The gradient normal to the segment of the gradient + !! tangetial to the segment of tangential flow along the segment + !! times the grid spacing [T-1 ~> s-1], with the first + !! index being the velocity/tracer point index along the segment, + !! and the second being 1 for the value 2 points into the domain + !! and 2 for the value 1 point into the domain. real, allocatable :: rx_norm_rad(:,:,:) !< The previous normal phase speed use for EW radiation !! OBC, in grid points per timestep [nondim] real, allocatable :: ry_norm_rad(:,:,:) !< The previous normal phase speed use for NS radiation !! OBC, in grid points per timestep [nondim] - real, allocatable :: rx_norm_obl(:,:,:) !< The previous normal radiation coefficient for EW - !! oblique OBCs [L2 T-2 ~> m2 s-2] - real, allocatable :: ry_norm_obl(:,:,:) !< The previous normal radiation coefficient for NS - !! oblique OBCs [L2 T-2 ~> m2 s-2] - real, allocatable :: cff_normal(:,:,:) !< The denominator for oblique radiation - !! for normal velocity [L2 T-2 ~> m2 s-2] + real, allocatable :: rx_norm_obl(:,:,:) !< The previous x-direction normalized radiation coefficient + !! for either EW or NS oblique OBCs [L2 T-2 ~> m2 s-2] + real, allocatable :: ry_norm_obl(:,:,:) !< The previous y-direction normalized radiation coefficient + !! for either EW or NS oblique OBCs [L2 T-2 ~> m2 s-2] + real, allocatable :: cff_normal(:,:,:) !< The denominator for oblique radiation of the normal + !! velocity [L2 T-2 ~> m2 s-2] real, allocatable :: nudged_normal_vel(:,:,:) !< The layer velocity normal to the OB segment !! that values should be nudged towards [L T-1 ~> m s-1]. real, allocatable :: nudged_tangential_vel(:,:,:) !< The layer velocity tangential to the OB segment - !! that values should be nudged towards [L T-1 ~> m s-1]. + !! that values should be nudged towards [L T-1 ~> m s-1], + !! discretized at the corner (PV) points. real, allocatable :: nudged_tangential_grad(:,:,:) !< The layer dvdx or dudy towards which nudging !! can occur [T-1 ~> s-1]. type(segment_tracer_registry_type), pointer :: tr_Reg=> NULL()!< A pointer to the tracer registry for the segment. @@ -304,9 +318,18 @@ module MOM_open_boundary !! grid points per timestep [nondim] real, allocatable :: ry_normal(:,:,:) !< Array storage for normal phase speed for NS radiation OBCs in units of !! grid points per timestep [nondim] - real, allocatable :: rx_oblique(:,:,:) !< Array storage for oblique boundary condition restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: ry_oblique(:,:,:) !< Array storage for oblique boundary condition restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: cff_normal(:,:,:) !< Array storage for oblique boundary condition restarts [L2 T-2 ~> m2 s-2] + real, allocatable :: rx_oblique_u(:,:,:) !< X-direction oblique boundary condition radiation speeds squared + !! at u points for restarts [L2 T-2 ~> m2 s-2] + real, allocatable :: ry_oblique_u(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared + !! at u points for restarts [L2 T-2 ~> m2 s-2] + real, allocatable :: rx_oblique_v(:,:,:) !< X-direction oblique boundary condition radiation speeds squared + !! at v points for restarts [L2 T-2 ~> m2 s-2] + real, allocatable :: ry_oblique_v(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared + !! at v points for restarts [L2 T-2 ~> m2 s-2] + real, allocatable :: cff_normal_u(:,:,:) !< Denominator for normalizing EW oblique boundary condition radiation + !! rates at u points for restarts [L2 T-2 ~> m2 s-2] + real, allocatable :: cff_normal_v(:,:,:) !< Denominator for normalizing NS oblique boundary condition radiation + !! rates at v points for restarts [L2 T-2 ~> m2 s-2] real, allocatable :: tres_x(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc] real, allocatable :: tres_y(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc] logical :: debug !< If true, write verbose checksums for debugging purposes. @@ -1794,9 +1817,11 @@ subroutine open_boundary_init(G, GV, US, param_file, OBC, restart_CS) id_clock_pass = cpu_clock_id('(Ocean OBC halo updates)', grain=CLOCK_ROUTINE) if (OBC%radiation_BCs_exist_globally) call pass_vector(OBC%rx_normal, OBC%ry_normal, G%Domain, & To_All+Scalar_Pair) - if (OBC%oblique_BCs_exist_globally) call pass_vector(OBC%rx_oblique, OBC%ry_oblique, G%Domain, & - To_All+Scalar_Pair) - if (allocated(OBC%cff_normal)) call pass_var(OBC%cff_normal, G%Domain, position=CORNER) + if (OBC%oblique_BCs_exist_globally) then + call pass_vector(OBC%rx_oblique_u, OBC%ry_oblique_v, G%Domain, To_All+Scalar_Pair) + call pass_vector(OBC%ry_oblique_u, OBC%rx_oblique_v, G%Domain, To_All+Scalar_Pair) + call pass_vector(OBC%cff_normal_u, OBC%cff_normal_v, G%Domain, To_All+Scalar_Pair) + endif if (allocated(OBC%tres_x) .and. allocated(OBC%tres_y)) then do m=1,OBC%ntr call pass_vector(OBC%tres_x(:,:,:,m), OBC%tres_y(:,:,:,m), G%Domain, To_All+Scalar_Pair) @@ -1811,45 +1836,6 @@ subroutine open_boundary_init(G, GV, US, param_file, OBC, restart_CS) enddo endif - ! The rx_normal and ry_normal arrays used with radiation OBCs are currently in units of grid - ! points per timestep, but if this were to be corrected to [L T-1 ~> m s-1] or [T-1 ~> s-1] to - ! permit timesteps to change between calls to the OBC code, the following would be needed: -! if ( OBC%radiation_BCs_exist_globally .and. (US%s_to_T_restart * US%m_to_L_restart /= 0.0) .and. & -! (US%s_to_T_restart /= US%m_to_L_restart) ) then -! vel_rescale = US%s_to_T_restart / US%m_to_L_restart -! if (query_initialized(OBC%rx_normal, "rx_normal", restart_CS)) then -! do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB -! OBC%rx_normal(I,j,k) = vel_rescale * OBC%rx_normal(I,j,k) -! enddo ; enddo ; enddo -! endif -! if (query_initialized(OBC%ry_normal, "ry_normal", restart_CS)) then -! do k=1,nz ; do J=JsdB,JedB ; do i=isd,ied -! OBC%ry_normal(i,J,k) = vel_rescale * OBC%ry_normal(i,J,k) -! enddo ; enddo ; enddo -! endif -! endif - - ! The oblique boundary condition terms have units of [L2 T-2 ~> m2 s-2] and may need to be rescaled. - if ( OBC%oblique_BCs_exist_globally .and. (US%s_to_T_restart * US%m_to_L_restart /= 0.0) .and. & - (US%s_to_T_restart /= US%m_to_L_restart) ) then - vel2_rescale = US%s_to_T_restart**2 / US%m_to_L_restart**2 - if (query_initialized(OBC%rx_oblique, "rx_oblique", restart_CS)) then - do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB - OBC%rx_oblique(I,j,k) = vel2_rescale * OBC%rx_oblique(I,j,k) - enddo ; enddo ; enddo - endif - if (query_initialized(OBC%ry_oblique, "ry_oblique", restart_CS)) then - do k=1,nz ; do J=JsdB,JedB ; do i=isd,ied - OBC%ry_oblique(i,J,k) = vel2_rescale * OBC%ry_oblique(i,J,k) - enddo ; enddo ; enddo - endif - if (query_initialized(OBC%cff_normal, "cff_normal", restart_CS)) then - do k=1,nz ; do J=JsdB,JedB ; do I=IsdB,IedB - OBC%cff_normal(I,J,k) = vel2_rescale * OBC%cff_normal(I,J,k) - enddo ; enddo ; enddo - endif - endif - end subroutine open_boundary_init logical function open_boundary_query(OBC, apply_open_OBC, apply_specified_OBC, apply_Flather_OBC, & @@ -1891,9 +1877,12 @@ subroutine open_boundary_dealloc(OBC) if (allocated(OBC%segnum_v)) deallocate(OBC%segnum_v) if (allocated(OBC%rx_normal)) deallocate(OBC%rx_normal) if (allocated(OBC%ry_normal)) deallocate(OBC%ry_normal) - if (allocated(OBC%rx_oblique)) deallocate(OBC%rx_oblique) - if (allocated(OBC%ry_oblique)) deallocate(OBC%ry_oblique) - if (allocated(OBC%cff_normal)) deallocate(OBC%cff_normal) + if (allocated(OBC%rx_oblique_u)) deallocate(OBC%rx_oblique_u) + if (allocated(OBC%ry_oblique_u)) deallocate(OBC%ry_oblique_u) + if (allocated(OBC%rx_oblique_v)) deallocate(OBC%rx_oblique_v) + if (allocated(OBC%ry_oblique_v)) deallocate(OBC%ry_oblique_v) + if (allocated(OBC%cff_normal_u)) deallocate(OBC%cff_normal_u) + if (allocated(OBC%cff_normal_v)) deallocate(OBC%cff_normal_v) if (allocated(OBC%tres_x)) deallocate(OBC%tres_x) if (allocated(OBC%tres_y)) deallocate(OBC%tres_y) deallocate(OBC) @@ -2129,12 +2118,17 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, real :: cff_new, cff_avg ! denominator in oblique [L2 T-2 ~> m2 s-2] real, allocatable, dimension(:,:,:) :: & rx_tang_rad, & ! The phase speed at u-points for tangential oblique OBCs - ! in units of grid points per timestep [nondim] + ! in units of grid points per timestep [nondim], + ! discretized at the corner (PV) points. ry_tang_rad, & ! The phase speed at v-points for tangential oblique OBCs - ! in units of grid points per timestep [nondim] - rx_tang_obl, & ! The x-coefficient for tangential oblique OBCs [L2 T-2 ~> m2 s-2] - ry_tang_obl, & ! The y-coefficient for tangential oblique OBCs [L2 T-2 ~> m2 s-2] - cff_tangential ! The denominator for tangential oblique OBCs [L2 T-2 ~> m2 s-2] + ! in units of grid points per timestep [nondim], + ! discretized at the corner (PV) points. + rx_tang_obl, & ! The x-coefficient for tangential oblique OBCs [L2 T-2 ~> m2 s-2], + ! discretized at the corner (PV) points. + ry_tang_obl, & ! The y-coefficient for tangential oblique OBCs [L2 T-2 ~> m2 s-2], + ! discretized at the corner (PV) points. + cff_tangential ! The denominator for tangential oblique OBCs [L2 T-2 ~> m2 s-2], + ! discretized at the corner (PV) points. real :: eps ! A small velocity squared [L2 T-2 ~> m2 s-2] type(OBC_segment_type), pointer :: segment => NULL() integer :: i, j, k, is, ie, js, je, m, nz, n @@ -2175,18 +2169,18 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, do k=1,GV%ke I=segment%HI%IsdB do j=segment%HI%jsd,segment%HI%jed - segment%rx_norm_obl(I,j,k) = OBC%rx_oblique(I,j,k) - segment%ry_norm_obl(I,j,k) = OBC%ry_oblique(I,j,k) - segment%cff_normal(I,j,k) = OBC%cff_normal(I,j,k) + segment%rx_norm_obl(I,j,k) = OBC%rx_oblique_u(I,j,k) + segment%ry_norm_obl(I,j,k) = OBC%ry_oblique_u(I,j,k) + segment%cff_normal(I,j,k) = OBC%cff_normal_u(I,j,k) enddo enddo elseif (segment%is_N_or_S .and. segment%oblique) then do k=1,GV%ke J=segment%HI%JsdB do i=segment%HI%isd,segment%HI%ied - segment%rx_norm_obl(i,J,k) = OBC%rx_oblique(i,J,k) - segment%ry_norm_obl(i,J,k) = OBC%ry_oblique(i,J,k) - segment%cff_normal(i,J,k) = OBC%cff_normal(i,J,k) + segment%rx_norm_obl(i,J,k) = OBC%rx_oblique_v(i,J,k) + segment%ry_norm_obl(i,J,k) = OBC%ry_oblique_v(i,J,k) + segment%cff_normal(i,J,k) = OBC%cff_normal_v(i,J,k) enddo enddo endif @@ -2269,16 +2263,16 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, ry_new = min(cff_new,max(dhdt*dhdy,-cff_new)) if (gamma_u < 1.0) then rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(I,j,k) + gamma_u*rx_new - ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,J,k) + gamma_u*ry_new - cff_avg = (1.0-gamma_u)*segment%cff_normal(i,J,k) + gamma_u*cff_new + ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(I,j,k) + gamma_u*ry_new + cff_avg = (1.0-gamma_u)*segment%cff_normal(I,j,k) + gamma_u*cff_new else rx_avg = rx_new ry_avg = ry_new cff_avg = cff_new endif segment%rx_norm_obl(I,j,k) = rx_avg - segment%ry_norm_obl(i,J,k) = ry_avg - segment%cff_normal(i,J,k) = cff_avg + segment%ry_norm_obl(I,j,k) = ry_avg + segment%cff_normal(I,j,k) = cff_avg segment%normal_vel(I,j,k) = ((cff_avg*u_new(I,j,k) + rx_avg*u_new(I-1,j,k)) - & (max(ry_avg,0.0)*segment%grad_normal(J-1,2,k) + & min(ry_avg,0.0)*segment%grad_normal(J,2,k))) / & @@ -2286,9 +2280,9 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, if (gamma_u < 1.0) then ! Copy restart fields into 3-d arrays. This is an inefficient and temporary ! implementation as a work-around to limitations in restart capability - OBC%rx_oblique(I,j,k) = segment%rx_norm_obl(I,j,k) - OBC%ry_oblique(i,J,k) = segment%ry_norm_obl(i,J,k) - OBC%cff_normal(I,j,k) = segment%cff_normal(I,j,k) + OBC%rx_oblique_u(I,j,k) = segment%rx_norm_obl(I,j,k) + OBC%ry_oblique_u(I,j,k) = segment%ry_norm_obl(I,j,k) + OBC%cff_normal_u(I,j,k) = segment%cff_normal(I,j,k) endif elseif (segment%gradient) then segment%normal_vel(I,j,k) = u_new(I-1,j,k) @@ -2409,9 +2403,9 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, cff_new = max(dhdx*dhdx + dhdy*dhdy, eps) rx_new = min(dhdt*dhdx, cff_new*rx_max) ry_new = min(cff_new,max(dhdt*dhdy,-cff_new)) - rx_tang_obl(I,j,k) = rx_new - ry_tang_obl(i,J,k) = ry_new - cff_tangential(i,J,k) = cff_new + rx_tang_obl(I,J,k) = rx_new + ry_tang_obl(I,J,k) = ry_new + cff_tangential(I,J,k) = cff_new enddo endif enddo @@ -2514,7 +2508,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, ry_new = min(cff_new,max(dhdt*dhdy,-cff_new)) if (gamma_u < 1.0) then rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(I,j,k) + gamma_u*rx_new - ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,J,k) + gamma_u*ry_new + ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(I,j,k) + gamma_u*ry_new cff_avg = (1.0-gamma_u)*segment%cff_normal(I,j,k) + gamma_u*cff_new else rx_avg = rx_new @@ -2522,8 +2516,8 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, cff_avg = cff_new endif segment%rx_norm_obl(I,j,k) = rx_avg - segment%ry_norm_obl(i,J,k) = ry_avg - segment%cff_normal(i,J,k) = cff_avg + segment%ry_norm_obl(I,j,k) = ry_avg + segment%cff_normal(I,j,k) = cff_avg segment%normal_vel(I,j,k) = ((cff_avg*u_new(I,j,k) + rx_avg*u_new(I+1,j,k)) - & (max(ry_avg,0.0)*segment%grad_normal(J-1,2,k) + & min(ry_avg,0.0)*segment%grad_normal(J,2,k))) / & @@ -2531,9 +2525,9 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, if (gamma_u < 1.0) then ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues ! implemented as a work-around to limitations in restart capability - OBC%rx_oblique(I,j,k) = segment%rx_norm_obl(I,j,k) - OBC%ry_oblique(i,J,k) = segment%ry_norm_obl(i,J,k) - OBC%cff_normal(I,j,k) = segment%cff_normal(I,j,k) + OBC%rx_oblique_u(I,j,k) = segment%rx_norm_obl(I,j,k) + OBC%ry_oblique_u(I,j,k) = segment%ry_norm_obl(I,j,k) + OBC%cff_normal_u(I,j,k) = segment%cff_normal(I,j,k) endif elseif (segment%gradient) then segment%normal_vel(I,j,k) = u_new(I+1,j,k) @@ -2654,9 +2648,9 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, cff_new = max(dhdx*dhdx + dhdy*dhdy, eps) rx_new = min(dhdt*dhdx, cff_new*rx_max) ry_new = min(cff_new,max(dhdt*dhdy,-cff_new)) - rx_tang_obl(I,j,k) = rx_new - ry_tang_obl(i,J,k) = ry_new - cff_tangential(i,J,k) = cff_new + rx_tang_obl(I,J,k) = rx_new + ry_tang_obl(I,J,k) = ry_new + cff_tangential(I,J,k) = cff_new enddo endif enddo @@ -2765,7 +2759,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, ry_avg = ry_new cff_avg = cff_new endif - segment%rx_norm_obl(I,j,k) = rx_avg + segment%rx_norm_obl(i,J,k) = rx_avg segment%ry_norm_obl(i,J,k) = ry_avg segment%cff_normal(i,J,k) = cff_avg segment%normal_vel(i,J,k) = ((cff_avg*v_new(i,J,k) + ry_avg*v_new(i,J-1,k)) - & @@ -2775,9 +2769,9 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, if (gamma_u < 1.0) then ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues ! implemented as a work-around to limitations in restart capability - OBC%rx_oblique(I,j,k) = segment%rx_norm_obl(I,j,k) - OBC%ry_oblique(i,J,k) = segment%ry_norm_obl(i,J,k) - OBC%cff_normal(i,J,k) = segment%cff_normal(i,J,k) + OBC%rx_oblique_v(i,J,k) = segment%rx_norm_obl(i,J,k) + OBC%ry_oblique_v(i,J,k) = segment%ry_norm_obl(i,J,k) + OBC%cff_normal_v(i,J,k) = segment%cff_normal(i,J,k) endif elseif (segment%gradient) then segment%normal_vel(i,J,k) = v_new(i,J-1,k) @@ -2898,9 +2892,9 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, cff_new = max(dhdx*dhdx + dhdy*dhdy, eps) ry_new = min(dhdt*dhdy, cff_new*ry_max) rx_new = min(cff_new,max(dhdt*dhdx,-cff_new)) - rx_tang_obl(I,j,k) = rx_new - ry_tang_obl(i,J,k) = ry_new - cff_tangential(i,J,k) = cff_new + rx_tang_obl(I,J,k) = rx_new + ry_tang_obl(I,J,k) = ry_new + cff_tangential(I,J,k) = cff_new enddo endif enddo @@ -3002,7 +2996,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, ry_new = min(dhdt*dhdy, cff_new*ry_max) rx_new = min(cff_new,max(dhdt*dhdx,-cff_new)) if (gamma_u < 1.0) then - rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(I,j,k) + gamma_u*rx_new + rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(i,J,k) + gamma_u*rx_new ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,J,k) + gamma_u*ry_new cff_avg = (1.0-gamma_u)*segment%cff_normal(i,J,k) + gamma_u*cff_new else @@ -3010,7 +3004,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, ry_avg = ry_new cff_avg = cff_new endif - segment%rx_norm_obl(I,j,k) = rx_avg + segment%rx_norm_obl(i,J,k) = rx_avg segment%ry_norm_obl(i,J,k) = ry_avg segment%cff_normal(i,J,k) = cff_avg segment%normal_vel(i,J,k) = ((cff_avg*v_new(i,J,k) + ry_avg*v_new(i,J+1,k)) - & @@ -3020,9 +3014,9 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, if (gamma_u < 1.0) then ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues ! implemented as a work-around to limitations in restart capability - OBC%rx_oblique(I,j,k) = segment%rx_norm_obl(I,j,k) - OBC%ry_oblique(i,J,k) = segment%ry_norm_obl(i,J,k) - OBC%cff_normal(i,J,k) = segment%cff_normal(i,J,k) + OBC%rx_oblique_v(i,J,k) = segment%rx_norm_obl(i,J,k) + OBC%ry_oblique_v(i,J,k) = segment%ry_norm_obl(i,J,k) + OBC%cff_normal_v(i,J,k) = segment%cff_normal(i,J,k) endif elseif (segment%gradient) then segment%normal_vel(i,J,k) = v_new(i,J+1,k) @@ -3143,9 +3137,9 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, cff_new = max(dhdx*dhdx + dhdy*dhdy, eps) ry_new = min(dhdt*dhdy, cff_new*ry_max) rx_new = min(cff_new,max(dhdt*dhdx,-cff_new)) - rx_tang_obl(I,j,k) = rx_new - ry_tang_obl(i,J,k) = ry_new - cff_tangential(i,J,k) = cff_new + rx_tang_obl(I,J,k) = rx_new + ry_tang_obl(I,J,k) = ry_new + cff_tangential(I,J,k) = cff_new enddo endif enddo @@ -4991,18 +4985,28 @@ subroutine open_boundary_register_restarts(HI, GV, US, OBC, Reg, param_file, res endif if (OBC%oblique_BCs_exist_globally) then - allocate(OBC%rx_oblique(HI%isdB:HI%iedB,HI%jsd:HI%jed,GV%ke), source=0.0) - allocate(OBC%ry_oblique(HI%isd:HI%ied,HI%jsdB:HI%jedB,GV%ke), source=0.0) - - vd(1) = var_desc("rx_oblique", "m2 s-2", "Radiation Speed Squared for EW oblique OBCs", 'u', 'L') - vd(2) = var_desc("ry_oblique", "m2 s-2", "Radiation Speed Squared for NS oblique OBCs", 'v', 'L') - call register_restart_pair(OBC%rx_oblique, OBC%ry_oblique, vd(1), vd(2), .false., & + allocate(OBC%rx_oblique_u(HI%isdB:HI%iedB,HI%jsd:HI%jed,GV%ke), source=0.0) + allocate(OBC%ry_oblique_u(HI%isdB:HI%iedB,HI%jsd:HI%jed,GV%ke), source=0.0) + allocate(OBC%cff_normal_u(HI%IsdB:HI%IedB,HI%jsd:HI%jed,GV%ke), source=0.0) + allocate(OBC%rx_oblique_v(HI%isd:HI%ied,HI%jsdB:HI%jedB,GV%ke), source=0.0) + allocate(OBC%ry_oblique_v(HI%isd:HI%ied,HI%jsdB:HI%jedB,GV%ke), source=0.0) + allocate(OBC%cff_normal_v(HI%isd:HI%ied,HI%jsdB:HI%jedB,GV%ke), source=0.0) + + vd(1) = var_desc("rx_oblique_u", "m2 s-2", "X-Direction Radiation Speed Squared for EW oblique OBCs", 'u', 'L') + vd(2) = var_desc("ry_oblique_v", "m2 s-2", "Y-Direction Radiation Speed Squared for NS oblique OBCs", 'v', 'L') + call register_restart_pair(OBC%rx_oblique_u, OBC%ry_oblique_v, vd(1), vd(2), .false., & + restart_CS, conversion=US%L_T_to_m_s**2) + vd(1) = var_desc("ry_oblique_u", "m2 s-2", "Y-Direction Radiation Speed Squared for EW oblique OBCs", 'u', 'L') + vd(2) = var_desc("rx_oblique_v", "m2 s-2", "X-Direction Radiation Speed Squared for NS oblique OBCs", 'v', 'L') + call register_restart_pair(OBC%ry_oblique_u, OBC%rx_oblique_v, vd(1), vd(2), .false., & restart_CS, conversion=US%L_T_to_m_s**2) - allocate(OBC%cff_normal(HI%IsdB:HI%IedB,HI%jsdB:HI%jedB,GV%ke), source=0.0) - call register_restart_field(OBC%cff_normal, "cff_normal", .false., restart_CS, & - longname="denominator for oblique OBCs", & - units="m2 s-2", conversion=US%L_T_to_m_s**2, hor_grid="q") + vd(1) = var_desc("norm_oblique_u", "m2 s-2", "Denominator for normalizing EW oblique OBC radiation rates", & + 'u', 'L') + vd(2) = var_desc("norm_oblique_v", "m2 s-2", "Denominator for normalizing NS oblique OBC radiation rates", & + 'v', 'L') + call register_restart_pair(OBC%cff_normal_u, OBC%cff_normal_v, vd(1), vd(2), .false., & + restart_CS, conversion=US%L_T_to_m_s**2) endif if (Reg%ntr == 0) return From 5fe02e018165caf894261af356205b11d6623717 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Tue, 27 Sep 2022 16:42:48 -0800 Subject: [PATCH 2/4] +Attempt to fix an imperfect restart issue. --- src/core/MOM_barotropic.F90 | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 67f07db3f3..f93be0ef47 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -1061,6 +1061,30 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, enddo ; enddo enddo + if (apply_OBCs) then + do n=1,OBC%number_of_segments + if (.not. OBC%segment(n)%on_pe) cycle + I = OBC%segment(n)%HI%IsdB ; J = OBC%segment(n)%HI%JsdB + if (OBC%segment(n)%is_N_or_S .and. (J >= Jsq-1) .and. (J <= Jeq+1)) then + do i = max(Isq-1,OBC%segment(n)%HI%isd), min(Ieq+2,OBC%segment(n)%HI%ied) + if (OBC%segment(n)%direction == OBC_DIRECTION_N) then + gtot_S(i,j+1) = gtot_S(i,j) + else ! (OBC%segment(n)%direction == OBC_DIRECTION_S) + gtot_N(i,j) = gtot_N(i,j+1) + endif + enddo + elseif (OBC%segment(n)%is_E_or_W .and. (I >= Isq-1) .and. (I <= Ieq+1)) then + do j = max(Jsq-1,OBC%segment(n)%HI%jsd), min(Jeq+2,OBC%segment(n)%HI%jed) + if (OBC%segment(n)%direction == OBC_DIRECTION_E) then + gtot_W(i+1,j) = gtot_W(i,j) + else ! (OBC%segment(n)%direction == OBC_DIRECTION_W) + gtot_E(i,j) = gtot_E(i+1,j) + endif + enddo + endif + enddo + endif + if (CS%tides) then call tidal_forcing_sensitivity(G, CS%tides_CSp, det_de) if (CS%tidal_sal_bug) then From 2bd92ba89c2a7d423b61563823a960eced415c04 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 9 Oct 2022 23:14:38 -0400 Subject: [PATCH 3/4] *+Added G%OBCmaskCu and G%OBCmaskCv Added two new arrays (OBCmaskCu and OBCmaskCv) to the ocean_grid_type and dyn_horgrid_type to mask out values over land or at open boundary condition points. Without open boundary conditions, these arrays are identical to mask2dCu and mask2dCv. These arrays are used in some of the lateral parameterization modules to zero out certain gradient-dependent fluxes at open boundary points. With these changes, the Bering test case solutions no longer exhibit any dependence on whether DEBUG_OBC is true or false or the value of OBC_SILLY_THICK, so this commit should help to address some of the issues discussed in github.com/NOAA-GFDL/MOM6/issues/208. All answers are bitwise identical in the MOM6-examples test cases, but they change for some other tests that use open boundary conditions more extensively, and there are new arrays in some transparent types. --- src/core/MOM_grid.F90 | 11 ++- src/core/MOM_open_boundary.F90 | 36 +++++-- src/core/MOM_transcribe_grid.F90 | 6 ++ src/framework/MOM_dyn_horgrid.F90 | 9 +- src/initialization/MOM_grid_initialize.F90 | 6 ++ src/parameterizations/lateral/MOM_MEKE.F90 | 4 +- .../lateral/MOM_interface_filter.F90 | 4 +- .../lateral/MOM_lateral_mixing_coeffs.F90 | 99 ++----------------- .../lateral/MOM_mixed_layer_restrat.F90 | 12 +-- .../lateral/MOM_thickness_diffuse.F90 | 28 +++--- 10 files changed, 88 insertions(+), 127 deletions(-) diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index 384e4a6d2b..f3a48f3ded 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -90,6 +90,7 @@ module MOM_grid real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: & mask2dCu, & !< 0 for boundary points and 1 for ocean points on the u grid [nondim]. + OBCmaskCu, & !< 0 for boundary or OBC points and 1 for ocean points on the u grid [nondim]. geoLatCu, & !< The geographic latitude at u points in degrees of latitude or m. geoLonCu, & !< The geographic longitude at u points in degrees of longitude or m. dxCu, & !< dxCu is delta x at u points [L ~> m]. @@ -102,6 +103,7 @@ module MOM_grid real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: & mask2dCv, & !< 0 for boundary points and 1 for ocean points on the v grid [nondim]. + OBCmaskCv, & !< 0 for boundary or OBC points and 1 for ocean points on the v grid [nondim]. geoLatCv, & !< The geographic latitude at v points in degrees of latitude or m. geoLonCv, & !< The geographic longitude at v points in degrees of longitude or m. dxCv, & !< dxCv is delta x at v points [L ~> m]. @@ -573,7 +575,9 @@ subroutine allocate_metrics(G) ALLOC_(G%mask2dT(isd:ied,jsd:jed)) ; G%mask2dT(:,:) = 0.0 ALLOC_(G%mask2dCu(IsdB:IedB,jsd:jed)) ; G%mask2dCu(:,:) = 0.0 + ALLOC_(G%OBCmaskCu(IsdB:IedB,jsd:jed)) ; G%OBCmaskCu(:,:) = 0.0 ALLOC_(G%mask2dCv(isd:ied,JsdB:JedB)) ; G%mask2dCv(:,:) = 0.0 + ALLOC_(G%OBCmaskCv(isd:ied,JsdB:JedB)) ; G%OBCmaskCv(:,:) = 0.0 ALLOC_(G%mask2dBu(IsdB:IedB,JsdB:JedB)) ; G%mask2dBu(:,:) = 0.0 ALLOC_(G%geoLatT(isd:ied,jsd:jed)) ; G%geoLatT(:,:) = 0.0 ALLOC_(G%geoLatCu(IsdB:IedB,jsd:jed)) ; G%geoLatCu(:,:) = 0.0 @@ -637,8 +641,8 @@ subroutine MOM_grid_end(G) DEALLOC_(G%areaCu) ; DEALLOC_(G%IareaCu) DEALLOC_(G%areaCv) ; DEALLOC_(G%IareaCv) - DEALLOC_(G%mask2dT) ; DEALLOC_(G%mask2dCu) - DEALLOC_(G%mask2dCv) ; DEALLOC_(G%mask2dBu) + DEALLOC_(G%mask2dT) ; DEALLOC_(G%mask2dCu) ; DEALLOC_(G%OBCmaskCu) + DEALLOC_(G%mask2dCv) ; DEALLOC_(G%OBCmaskCv) ; DEALLOC_(G%mask2dBu) DEALLOC_(G%geoLatT) ; DEALLOC_(G%geoLatCu) DEALLOC_(G%geoLatCv) ; DEALLOC_(G%geoLatBu) @@ -686,6 +690,7 @@ end subroutine MOM_grid_end !! !! Each location also has a 2D mask indicating whether the entire column is land or ocean. !! `mask2dT` is 1 if the column is wet or 0 if the T-cell is land. -!! `mask2dCu` is 1 if both neighboring column are ocean, and 0 if either is land. +!! `mask2dCu` is 1 if both neighboring columns are ocean, and 0 if either is land. +!! `OBCmasku` is 1 if both neighboring columns are ocean, and 0 if either is land of if this is OBC point. end module MOM_grid diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index a187ac5b2f..1cc8505d17 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -137,7 +137,8 @@ module MOM_open_boundary logical :: specified !< Boundary normal velocity fixed to external value. logical :: specified_tan !< Boundary tangential velocity fixed to external value. logical :: specified_grad !< Boundary gradient of tangential velocity fixed to external value. - logical :: open !< Boundary is open for continuity solver. + logical :: open !< Boundary is open for continuity solver, and there are no other + !! parameterized mass fluxes at the open boundary. logical :: gradient !< Zero gradient at boundary. logical :: values_needed !< Whether or not any external OBC fields are needed. logical :: u_values_needed !< Whether or not external u OBC fields are needed. @@ -1963,16 +1964,16 @@ subroutine open_boundary_impose_land_mask(OBC, G, areaCu, areaCv, US) do j=segment%HI%jsd,segment%HI%jed if (G%mask2dCu(I,j) == 0) OBC%segnum_u(I,j) = OBC_NONE if (segment%direction == OBC_DIRECTION_W) then - G%mask2dT(i,j) = 0 + G%mask2dT(i,j) = 0.0 else - G%mask2dT(i+1,j) = 0 + G%mask2dT(i+1,j) = 0.0 endif enddo do J=segment%HI%JsdB+1,segment%HI%JedB-1 if (segment%direction == OBC_DIRECTION_W) then - G%mask2dCv(i,J) = 0 + G%mask2dCv(i,J) = 0 ; G%OBCmaskCv(i,J) = 0.0 else - G%mask2dCv(i+1,J) = 0 + G%mask2dCv(i+1,J) = 0.0 ; G%OBCmaskCv(i+1,J) = 0.0 endif enddo else @@ -1981,21 +1982,38 @@ subroutine open_boundary_impose_land_mask(OBC, G, areaCu, areaCv, US) do i=segment%HI%isd,segment%HI%ied if (G%mask2dCv(i,J) == 0) OBC%segnum_v(i,J) = OBC_NONE if (segment%direction == OBC_DIRECTION_S) then - G%mask2dT(i,j) = 0 + G%mask2dT(i,j) = 0.0 else - G%mask2dT(i,j+1) = 0 + G%mask2dT(i,j+1) = 0.0 endif enddo do I=segment%HI%IsdB+1,segment%HI%IedB-1 if (segment%direction == OBC_DIRECTION_S) then - G%mask2dCu(I,j) = 0 + G%mask2dCu(I,j) = 0.0 ; G%OBCmaskCu(I,j) = 0.0 else - G%mask2dCu(I,j+1) = 0 + G%mask2dCu(I,j+1) = 0.0 ; G%OBCmaskCu(I,j+1) = 0.0 endif enddo endif enddo + do n=1,OBC%number_of_segments + segment=>OBC%segment(n) + if (.not. (segment%on_pe .and. segment%open)) cycle + ! Set the OBCmask values to help eliminate certain terms at u- or v- OBC points. + if (segment%is_E_or_W) then + I=segment%HI%IsdB + do j=segment%HI%jsd,segment%HI%jed + G%OBCmaskCu(I,j) = 0.0 + enddo + else + J=segment%HI%JsdB + do i=segment%HI%isd,segment%HI%ied + G%OBCmaskCv(i,J) = 0.0 + enddo + endif + enddo + do n=1,OBC%number_of_segments segment=>OBC%segment(n) if (.not. segment%on_pe .or. .not. segment%specified) cycle diff --git a/src/core/MOM_transcribe_grid.F90 b/src/core/MOM_transcribe_grid.F90 index 7ab15d542e..8f8da21ef3 100644 --- a/src/core/MOM_transcribe_grid.F90 +++ b/src/core/MOM_transcribe_grid.F90 @@ -76,6 +76,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US) oG%porous_DavgU(I,j) = dG%porous_DavgU(I+ido,j+jdo) - oG%Z_ref oG%mask2dCu(I,j) = dG%mask2dCu(I+ido,j+jdo) + oG%OBCmaskCu(I,j) = dG%OBCmaskCu(I+ido,j+jdo) oG%areaCu(I,j) = dG%areaCu(I+ido,j+jdo) oG%IareaCu(I,j) = dG%IareaCu(I+ido,j+jdo) enddo ; enddo @@ -92,6 +93,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US) oG%porous_DavgV(i,J) = dG%porous_DavgV(i+ido,J+jdo) - oG%Z_ref oG%mask2dCv(i,J) = dG%mask2dCv(i+ido,J+jdo) + oG%OBCmaskCv(i,J) = dG%OBCmaskCv(i+ido,J+jdo) oG%areaCv(i,J) = dG%areaCv(i+ido,J+jdo) oG%IareaCv(i,J) = dG%IareaCv(i+ido,J+jdo) enddo ; enddo @@ -152,6 +154,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US) call pass_vector(oG%dxCu, oG%dyCv, oG%Domain, To_All+Scalar_Pair, CGRID_NE) call pass_vector(oG%dy_Cu, oG%dx_Cv, oG%Domain, To_All+Scalar_Pair, CGRID_NE) call pass_vector(oG%mask2dCu, oG%mask2dCv, oG%Domain, To_All+Scalar_Pair, CGRID_NE) + call pass_vector(oG%OBCmaskCu, oG%OBCmaskCv, oG%Domain, To_All+Scalar_Pair, CGRID_NE) call pass_vector(oG%IareaCu, oG%IareaCv, oG%Domain, To_All+Scalar_Pair, CGRID_NE) call pass_vector(oG%IareaCu, oG%IareaCv, oG%Domain, To_All+Scalar_Pair, CGRID_NE) call pass_vector(oG%geoLatCu, oG%geoLatCv, oG%Domain, To_All+Scalar_Pair, CGRID_NE) @@ -230,6 +233,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) dG%porous_DavgU(I,j) = oG%porous_DavgU(I+ido,j+jdo) + oG%Z_ref dG%mask2dCu(I,j) = oG%mask2dCu(I+ido,j+jdo) + dG%OBCmaskCu(I,j) = oG%OBCmaskCu(I+ido,j+jdo) dG%areaCu(I,j) = oG%areaCu(I+ido,j+jdo) dG%IareaCu(I,j) = oG%IareaCu(I+ido,j+jdo) enddo ; enddo @@ -246,6 +250,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) dG%porous_DavgV(i,J) = oG%porous_DavgU(i+ido,J+jdo) + oG%Z_ref dG%mask2dCv(i,J) = oG%mask2dCv(i+ido,J+jdo) + dG%OBCmaskCv(i,J) = oG%OBCmaskCv(i+ido,J+jdo) dG%areaCv(i,J) = oG%areaCv(i+ido,J+jdo) dG%IareaCv(i,J) = oG%IareaCv(i+ido,J+jdo) enddo ; enddo @@ -307,6 +312,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) call pass_vector(dG%dxCu, dG%dyCv, dG%Domain, To_All+Scalar_Pair, CGRID_NE) call pass_vector(dG%dy_Cu, dG%dx_Cv, dG%Domain, To_All+Scalar_Pair, CGRID_NE) call pass_vector(dG%mask2dCu, dG%mask2dCv, dG%Domain, To_All+Scalar_Pair, CGRID_NE) + call pass_vector(dG%OBCmaskCu, dG%OBCmaskCv, dG%Domain, To_All+Scalar_Pair, CGRID_NE) call pass_vector(dG%IareaCu, dG%IareaCv, dG%Domain, To_All+Scalar_Pair, CGRID_NE) call pass_vector(dG%IareaCu, dG%IareaCv, dG%Domain, To_All+Scalar_Pair, CGRID_NE) call pass_vector(dG%geoLatCu, dG%geoLatCv, dG%Domain, To_All+Scalar_Pair, CGRID_NE) diff --git a/src/framework/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90 index e97c8d981b..60c30d8e94 100644 --- a/src/framework/MOM_dyn_horgrid.F90 +++ b/src/framework/MOM_dyn_horgrid.F90 @@ -87,6 +87,7 @@ module MOM_dyn_horgrid real, allocatable, dimension(:,:) :: & mask2dCu, & !< 0 for boundary points and 1 for ocean points on the u grid [nondim]. + OBCmaskCu, & !< 0 for boundary or OBC points and 1 for ocean points on the u grid [nondim]. geoLatCu, & !< The geographic latitude at u points [degrees of latitude] or [m]. geoLonCu, & !< The geographic longitude at u points [degrees of longitude] or [m]. dxCu, & !< dxCu is delta x at u points [L ~> m]. @@ -99,6 +100,7 @@ module MOM_dyn_horgrid real, allocatable, dimension(:,:) :: & mask2dCv, & !< 0 for boundary points and 1 for ocean points on the v grid [nondim]. + OBCmaskCv, & !< 0 for boundary or OBC points and 1 for ocean points on the v grid [nondim]. geoLatCv, & !< The geographic latitude at v points [degrees of latitude] or [m]. geoLonCv, & !< The geographic longitude at v points [degrees of longitude] or [m]. dxCv, & !< dxCv is delta x at v points [L ~> m]. @@ -250,6 +252,8 @@ subroutine create_dyn_horgrid(G, HI, bathymetry_at_vel) allocate(G%mask2dCu(IsdB:IedB,jsd:jed), source=0.0) allocate(G%mask2dCv(isd:ied,JsdB:JedB), source=0.0) allocate(G%mask2dBu(IsdB:IedB,JsdB:JedB), source=0.0) + allocate(G%OBCmaskCu(IsdB:IedB,jsd:jed), source=0.0) + allocate(G%OBCmaskCv(isd:ied,JsdB:JedB), source=0.0) allocate(G%geoLatT(isd:ied,jsd:jed), source=0.0) allocate(G%geoLatCu(IsdB:IedB,jsd:jed), source=0.0) allocate(G%geoLatCv(isd:ied,JsdB:JedB), source=0.0) @@ -331,6 +335,7 @@ subroutine rotate_dyn_horgrid(G_in, G, US, turns) call rotate_array_pair(G_in%dx_Cv, G_in%dy_Cu, turns, G%dx_Cv, G%dy_Cu) call rotate_array_pair(G_in%mask2dCu, G_in%mask2dCv, turns, G%mask2dCu, G%mask2dCv) + call rotate_array_pair(G_in%OBCmaskCu, G_in%OBCmaskCv, turns, G%OBCmaskCu, G%OBCmaskCv) call rotate_array_pair(G_in%areaCu, G_in%areaCv, turns, G%areaCu, G%areaCv) call rotate_array_pair(G_in%IareaCu, G_in%IareaCv, turns, G%IareaCu, G%IareaCv) @@ -501,8 +506,8 @@ subroutine destroy_dyn_horgrid(G) deallocate(G%areaCu) ; deallocate(G%IareaCu) deallocate(G%areaCv) ; deallocate(G%IareaCv) - deallocate(G%mask2dT) ; deallocate(G%mask2dCu) - deallocate(G%mask2dCv) ; deallocate(G%mask2dBu) + deallocate(G%mask2dT) ; deallocate(G%mask2dCu) ; deallocate(G%OBCmaskCu) + deallocate(G%mask2dCv) ; deallocate(G%OBCmaskCv) ; deallocate(G%mask2dBu) deallocate(G%geoLatT) ; deallocate(G%geoLatCu) deallocate(G%geoLatCv) ; deallocate(G%geoLatBu) diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index bc004daa95..d84d2275e4 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -1206,6 +1206,8 @@ subroutine initialize_masks(G, PF, US) else G%mask2dCu(I,j) = 1.0 endif + ! This mask may be revised later after the open boundary positions are specified. + G%OBCmaskCu(I,j) = G%mask2dCu(I,j) enddo ; enddo do J=G%jsd,G%jed-1 ; do i=G%isd,G%ied @@ -1214,6 +1216,8 @@ subroutine initialize_masks(G, PF, US) else G%mask2dCv(i,J) = 1.0 endif + ! This mask may be revised later after the open boundary positions are specified. + G%OBCmaskCv(i,J) = G%mask2dCv(i,J) enddo ; enddo do J=G%jsd,G%jed-1 ; do I=G%isd,G%ied-1 @@ -1229,12 +1233,14 @@ subroutine initialize_masks(G, PF, US) call pass_vector(G%mask2dCu, G%mask2dCv, G%Domain, To_All+Scalar_Pair, CGRID_NE) do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB + ! This open face length may be revised later. G%dy_Cu(I,j) = G%mask2dCu(I,j) * G%dyCu(I,j) G%areaCu(I,j) = G%dxCu(I,j) * G%dy_Cu(I,j) G%IareaCu(I,j) = G%mask2dCu(I,j) * Adcroft_reciprocal(G%areaCu(I,j)) enddo ; enddo do J=G%JsdB,G%JedB ; do i=G%isd,G%ied + ! This open face length may be revised later. G%dx_Cv(i,J) = G%mask2dCv(i,J) * G%dxCv(i,J) G%areaCv(i,J) = G%dyCv(i,J) * G%dx_Cv(i,J) G%IareaCv(i,J) = G%mask2dCv(i,J) * Adcroft_reciprocal(G%areaCv(i,J)) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 9b024e62b0..ead6086346 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -462,7 +462,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h !$OMP parallel do default(shared) do j=js-1,je+1 ; do I=is-2,ie+1 ! MEKE_uflux is used here as workspace with units of [L2 T-2 ~> m2 s-2]. - MEKE_uflux(I,j) = ((G%dy_Cu(I,j)*G%IdxCu(I,j)) * G%mask2dCu(I,j)) * & + MEKE_uflux(I,j) = ((G%dy_Cu(I,j)*G%IdxCu(I,j)) * G%OBCmaskCu(I,j)) * & (MEKE%MEKE(i+1,j) - MEKE%MEKE(i,j)) ! This would have units of [R Z L2 T-2 ~> kg s-2] ! MEKE_uflux(I,j) = ((G%dy_Cu(I,j)*G%IdxCu(I,j)) * & @@ -472,7 +472,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h !$OMP parallel do default(shared) do J=js-2,je+1 ; do i=is-1,ie+1 ! MEKE_vflux is used here as workspace with units of [L2 T-2 ~> m2 s-2]. - MEKE_vflux(i,J) = ((G%dx_Cv(i,J)*G%IdyCv(i,J)) * G%mask2dCv(i,J)) * & + MEKE_vflux(i,J) = ((G%dx_Cv(i,J)*G%IdyCv(i,J)) * G%OBCmaskCv(i,J)) * & (MEKE%MEKE(i,j+1) - MEKE%MEKE(i,j)) ! This would have units of [R Z L2 T-2 ~> kg s-2] ! MEKE_vflux(i,J) = ((G%dx_Cv(i,J)*G%IdyCv(i,J)) * & diff --git a/src/parameterizations/lateral/MOM_interface_filter.F90 b/src/parameterizations/lateral/MOM_interface_filter.F90 index feed4bd930..dd082f1558 100644 --- a/src/parameterizations/lateral/MOM_interface_filter.F90 +++ b/src/parameterizations/lateral/MOM_interface_filter.F90 @@ -284,7 +284,7 @@ subroutine filter_interface(h, e, Lsm2_u, Lsm2_v, uhD, vhD, G, GV, US, halo_size do I=is-1,ie ; uhtot(I,j) = 0.0 ; enddo do K=nz,2,-1 do I=is-1,ie - Slope = ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) * G%mask2dCu(I,j) + Slope = ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) * G%OBCmaskCu(I,j) Sfn_est = (Lsm2_u(I,j)*G%dy_Cu(I,j)) * (GV%Z_to_H * Slope) @@ -316,7 +316,7 @@ subroutine filter_interface(h, e, Lsm2_u, Lsm2_v, uhD, vhD, G, GV, US, halo_size do i=is,ie ; vhtot(i,J) = 0.0 ; enddo do K=nz,2,-1 do i=is,ie - Slope = ((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) * G%mask2dCv(i,J) + Slope = ((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) * G%OBCmaskCv(i,J) Sfn_est = (Lsm2_v(i,J)*G%dx_Cv(i,J)) * (GV%Z_to_H * Slope) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index dc23042916..e6dd57fc2e 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -533,7 +533,6 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C real :: H_u(SZIB_(G)), H_v(SZI_(G)) real :: S2_u(SZIB_(G), SZJ_(G)) real :: S2_v(SZI_(G), SZJB_(G)) - logical :: local_open_u_BC, local_open_v_BC if (.not. CS%initialized) call MOM_error(FATAL, "calc_Visbeck_coeffs_old: "// & "Module must be initialized before it is used.") @@ -546,13 +545,6 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - local_open_u_BC = .false. - local_open_v_BC = .false. - if (associated(OBC)) then - local_open_u_BC = OBC%open_u_BCs_exist_globally - local_open_v_BC = OBC%open_v_BCs_exist_globally - endif - S2max = CS%Visbeck_S_max**2 !$OMP parallel do default(shared) @@ -593,20 +585,11 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C enddo ; enddo do I=is-1,ie if (H_u(I)>0.) then - CS%SN_u(I,j) = G%mask2dCu(I,j) * CS%SN_u(I,j) / H_u(I) - S2_u(I,j) = G%mask2dCu(I,j) * S2_u(I,j) / H_u(I) + CS%SN_u(I,j) = G%OBCmaskCu(I,j) * CS%SN_u(I,j) / H_u(I) + S2_u(I,j) = G%OBCmaskCu(I,j) * S2_u(I,j) / H_u(I) else CS%SN_u(I,j) = 0. endif - if (local_open_u_BC) then - l_seg = OBC%segnum_u(I,j) - - if (l_seg /= OBC_NONE) then - if (OBC%segment(l_seg)%open) then - CS%SN_u(i,J) = 0. - endif - endif - endif enddo enddo @@ -638,20 +621,11 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C enddo ; enddo do i=is,ie if (H_v(i)>0.) then - CS%SN_v(i,J) = G%mask2dCv(i,J) * CS%SN_v(i,J) / H_v(i) - S2_v(i,J) = G%mask2dCv(i,J) * S2_v(i,J) / H_v(i) + CS%SN_v(i,J) = G%OBCmaskCv(i,J) * CS%SN_v(i,J) / H_v(i) + S2_v(i,J) = G%OBCmaskCv(i,J) * S2_v(i,J) / H_v(i) else CS%SN_v(i,J) = 0. endif - if (local_open_v_BC) then - l_seg = OBC%segnum_v(i,J) - - if (l_seg /= OBC_NONE) then - if (OBC%segment(OBC%segnum_v(i,J))%open) then - CS%SN_v(i,J) = 0. - endif - endif - endif enddo enddo @@ -699,7 +673,7 @@ subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, OBC, h, e, dzu, dzv, dzSxN, d real :: r_crp_dist ! The inverse of the distance over which to scale the cropping [Z-1 ~> m-1] real :: dB, dT ! Elevation variables used when cropping [Z ~> m] integer :: i, j, k, l_seg - logical :: local_open_u_BC, local_open_v_BC, crop + logical :: crop dz_neglect = GV%H_subroundoff * GV%H_to_Z D_scale = CS%Eady_GR_D_scale @@ -707,13 +681,6 @@ subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, OBC, h, e, dzu, dzv, dzSxN, d r_crp_dist = 1. / max( dz_neglect, CS%cropping_distance ) crop = CS%cropping_distance>=0. ! Only filter out in-/out-cropped interface is parameter if non-negative - local_open_u_BC = .false. - local_open_v_BC = .false. - if (associated(OBC)) then - local_open_u_BC = OBC%open_u_BCs_exist_globally - local_open_v_BC = OBC%open_v_BCs_exist_globally - endif - if (CS%debug) then call uvchksum("calc_Eady_growth_rate_2D dz[uv]", dzu, dzv, G%HI, scale=US%Z_to_m, scalar_pair=.true.) call uvchksum("calc_Eady_growth_rate_2D dzS2N2[uv]", dzSxN, dzSyN, G%HI, & @@ -764,19 +731,9 @@ subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, OBC, h, e, dzu, dzv, dzSxN, d enddo ; enddo endif do I=G%isc-1,G%iec - CS%SN_u(I,j) = G%mask2dCu(I,j) * ( vint_SN(I) / sum_dz(I) ) - SN_cpy(I,j) = G%mask2dCu(I,j) * ( vint_SN(I) / sum_dz(I) ) + CS%SN_u(I,j) = G%OBCmaskCu(I,j) * ( vint_SN(I) / sum_dz(I) ) + SN_cpy(I,j) = G%OBCmaskCu(I,j) * ( vint_SN(I) / sum_dz(I) ) enddo - if (local_open_u_BC) then - do I=G%isc-1,G%iec - l_seg = OBC%segnum_u(I,j) - if (l_seg /= OBC_NONE) then - if (OBC%segment(l_seg)%open) then - CS%SN_u(i,J) = 0. - endif - endif - enddo - endif enddo !$OMP parallel do default(shared) private(dnew,dz,weight,l_seg) @@ -817,18 +774,8 @@ subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, OBC, h, e, dzu, dzv, dzSxN, d enddo ; enddo endif do i=G%isc-1,G%iec+1 - CS%SN_v(i,J) = G%mask2dCv(i,J) * ( vint_SN(i) / sum_dz(i) ) + CS%SN_v(i,J) = G%OBCmaskCv(i,J) * ( vint_SN(i) / sum_dz(i) ) enddo - if (local_open_v_BC) then - do i=G%isc-1,G%iec+1 - l_seg = OBC%segnum_v(i,J) - if (l_seg /= OBC_NONE) then - if (OBC%segment(l_seg)%open) then - CS%SN_v(i,J) = 0. - endif - endif - enddo - endif enddo do j = G%jsc,G%jec @@ -881,7 +828,6 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop integer :: l_seg real :: S2N2_u_local(SZIB_(G), SZJ_(G),SZK_(GV)) real :: S2N2_v_local(SZI_(G), SZJB_(G),SZK_(GV)) - logical :: local_open_u_BC, local_open_v_BC if (.not. CS%initialized) call MOM_error(FATAL, "calc_slope_functions_using_just_e: "// & "Module must be initialized before it is used.") @@ -894,13 +840,6 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - local_open_u_BC = .false. - local_open_v_BC = .false. - if (associated(OBC)) then - local_open_u_BC = OBC%open_u_BCs_exist_globally - local_open_v_BC = OBC%open_v_BCs_exist_globally - endif - one_meter = 1.0 * GV%m_to_H h_neglect = GV%H_subroundoff H_cutoff = real(2*nz) * (GV%Angstrom_H + h_neglect) @@ -972,20 +911,11 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop !SN_u(I,j) = sqrt( SN_u(I,j) / ( max(G%bathyT(i,j), G%bathyT(i+1,j)) + (G%Z_ref + GV%Angstrom_Z) ) ) !The code below behaves better than the line above. Not sure why? AJA if ( min(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref > H_cutoff*GV%H_to_Z ) then - CS%SN_u(I,j) = G%mask2dCu(I,j) * sqrt( CS%SN_u(I,j) / & + CS%SN_u(I,j) = G%OBCmaskCu(I,j) * sqrt( CS%SN_u(I,j) / & (max(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref) ) else CS%SN_u(I,j) = 0.0 endif - if (local_open_u_BC) then - l_seg = OBC%segnum_u(I,j) - - if (l_seg /= OBC_NONE) then - if (OBC%segment(l_seg)%open) then - CS%SN_u(I,j) = 0. - endif - endif - endif enddo enddo !$OMP parallel do default(shared) @@ -999,20 +929,11 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop !SN_v(i,J) = sqrt( SN_v(i,J) / ( max(G%bathyT(i,J), G%bathyT(i,J+1)) + (G%Z_ref + GV%Angstrom_Z) ) ) !The code below behaves better than the line above. Not sure why? AJA if ( min(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref > H_cutoff*GV%H_to_Z ) then - CS%SN_v(i,J) = G%mask2dCv(i,J) * sqrt( CS%SN_v(i,J) / & + CS%SN_v(i,J) = G%OBCmaskCv(i,J) * sqrt( CS%SN_v(i,J) / & (max(G%bathyT(i,j), G%bathyT(i,j+1)) + G%Z_ref) ) else CS%SN_v(i,J) = 0.0 endif - if (local_open_v_BC) then - l_seg = OBC%segnum_v(i,J) - - if (l_seg /= OBC_NONE) then - if (OBC%segment(OBC%segnum_v(i,J))%open) then - CS%SN_v(i,J) = 0. - endif - endif - endif enddo enddo diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index a8efa4cf12..0aef33ddc6 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -393,7 +393,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2) timescale = timescale * CS%ml_restrat_coef if (res_upscale) timescale = timescale * res_scaling_fac - uDml(I) = timescale * G%mask2dCu(I,j)*G%dyCu(I,j)*G%IdxCu(I,j) * & + uDml(I) = timescale * G%OBCmaskCu(I,j)*G%dyCu(I,j)*G%IdxCu(I,j) * & (Rml_av_fast(i+1,j)-Rml_av_fast(i,j)) * (h_vel**2 * GV%Z_to_H) ! As above but using the slow filtered MLD h_vel = 0.5*((htot_slow(i,j) + htot_slow(i+1,j)) + h_neglect) * GV%H_to_Z @@ -402,7 +402,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2) timescale = timescale * CS%ml_restrat_coef2 if (res_upscale) timescale = timescale * res_scaling_fac - uDml_slow(I) = timescale * G%mask2dCu(I,j)*G%dyCu(I,j)*G%IdxCu(I,j) * & + uDml_slow(I) = timescale * G%OBCmaskCu(I,j)*G%dyCu(I,j)*G%IdxCu(I,j) * & (Rml_av_slow(i+1,j)-Rml_av_slow(i,j)) * (h_vel**2 * GV%Z_to_H) if (uDml(I) + uDml_slow(I) == 0.) then @@ -468,7 +468,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2) timescale = timescale * CS%ml_restrat_coef if (res_upscale) timescale = timescale * res_scaling_fac - vDml(i) = timescale * G%mask2dCv(i,J)*G%dxCv(i,J)*G%IdyCv(i,J) * & + vDml(i) = timescale * G%OBCmaskCv(i,J)*G%dxCv(i,J)*G%IdyCv(i,J) * & (Rml_av_fast(i,j+1)-Rml_av_fast(i,j)) * (h_vel**2 * GV%Z_to_H) ! As above but using the slow filtered MLD h_vel = 0.5*((htot_slow(i,j) + htot_slow(i,j+1)) + h_neglect) * GV%H_to_Z @@ -477,7 +477,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2) timescale = timescale * CS%ml_restrat_coef2 if (res_upscale) timescale = timescale * res_scaling_fac - vDml_slow(i) = timescale * G%mask2dCv(i,J)*G%dxCv(i,J)*G%IdyCv(i,J) * & + vDml_slow(i) = timescale * G%OBCmaskCv(i,J)*G%dxCv(i,J)*G%IdyCv(i,J) * & (Rml_av_slow(i,j+1)-Rml_av_slow(i,j)) * (h_vel**2 * GV%Z_to_H) if (vDml(i) + vDml_slow(i) == 0.) then @@ -716,7 +716,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) timescale = timescale * CS%ml_restrat_coef ! timescale = timescale*(2?)*(L_def/L_MLI) * min(EKE/MKE,1.0 + (G%dyCv(i,j)/L_def)**2) - uDml(I) = timescale * G%mask2dCu(I,j)*G%dyCu(I,j)*G%IdxCu(I,j) * & + uDml(I) = timescale * G%OBCmaskCu(I,j)*G%dyCu(I,j)*G%IdxCu(I,j) * & (Rml_av(i+1,j)-Rml_av(i,j)) * (h_vel**2 * GV%Z_to_H) if (uDml(I) == 0) then @@ -762,7 +762,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) timescale = timescale * CS%ml_restrat_coef ! timescale = timescale*(2?)*(L_def/L_MLI) * min(EKE/MKE,1.0 + (G%dyCv(i,j)/L_def)**2) - vDml(i) = timescale * G%mask2dCv(i,J)*G%dxCv(i,J)*G%IdyCv(i,J) * & + vDml(i) = timescale * G%OBCmaskCv(i,J)*G%dxCv(i,J)*G%IdyCv(i,J) * & (Rml_av(i,j+1)-Rml_av(i,j)) * (h_vel**2 * GV%Z_to_H) if (vDml(i) == 0) then do k=1,nkml ; vhml(i,J,k) = 0.0 ; enddo diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 3cab1030da..c7310e1560 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -233,7 +233,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp if (CS%MEKE_GEOMETRIC) then !$OMP do do j=js,je ; do I=is-1,ie - Khth_loc_u(I,j) = Khth_loc_u(I,j) + G%mask2dCu(I,j) * CS%MEKE_GEOMETRIC_alpha * & + Khth_loc_u(I,j) = Khth_loc_u(I,j) + G%OBCmaskCu(I,j) * CS%MEKE_GEOMETRIC_alpha * & 0.5*(MEKE%MEKE(i,j)+MEKE%MEKE(i+1,j)) / & (VarMix%SN_u(I,j) + CS%MEKE_GEOMETRIC_epsilon) enddo ; enddo @@ -319,7 +319,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp if (CS%MEKE_GEOMETRIC) then !$OMP do do J=js-1,je ; do i=is,ie - Khth_loc_v(i,J) = Khth_loc_v(i,J) + G%mask2dCv(i,J) * CS%MEKE_GEOMETRIC_alpha * & + Khth_loc_v(i,J) = Khth_loc_v(i,J) + G%OBCmaskCv(i,J) * CS%MEKE_GEOMETRIC_alpha * & 0.5*(MEKE%MEKE(i,j)+MEKE%MEKE(i,j+1)) / & (VarMix%SN_v(i,J) + CS%MEKE_GEOMETRIC_epsilon) enddo ; enddo @@ -956,7 +956,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (present_slope_x) then Slope = slope_x(I,j,k) else - Slope = ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) * G%mask2dCu(I,j) + Slope = ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) * G%OBCmaskCu(I,j) endif if (CS%id_slope_x > 0) CS%diagSlopeX(I,j,k) = Slope Sfn_unlim_u(I,K) = ((KH_u(I,j,K)*G%dy_Cu(I,j))*Slope) @@ -971,7 +971,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV enddo ! k-loop if (CS%use_FGNV_streamfn) then - do k=1,nz ; do I=is-1,ie ; if (G%mask2dCu(I,j)>0.) then + do k=1,nz ; do I=is-1,ie ; if (G%OBCmaskCu(I,j)>0.) then h_harm = max( h_neglect, & 2. * h(i,j,k) * h(i+1,j,k) / ( ( h(i,j,k) + h(i+1,j,k) ) + h_neglect ) ) c2_h_u(I,k) = CS%FGNV_scale * & @@ -980,7 +980,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! Solve an elliptic equation for the streamfunction following Ferrari et al., 2010. do I=is-1,ie - if (G%mask2dCu(I,j)>0.) then + if (G%OBCmaskCu(I,j)>0.) then do K=2,nz Sfn_unlim_u(I,K) = (1. + CS%FGNV_scale) * Sfn_unlim_u(I,K) enddo @@ -1238,7 +1238,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (present_slope_y) then Slope = slope_y(i,J,k) else - Slope = ((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) * G%mask2dCv(i,J) + Slope = ((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) * G%OBCmaskCv(i,J) endif if (CS%id_slope_y > 0) CS%diagSlopeY(I,j,k) = Slope Sfn_unlim_v(i,K) = ((KH_v(i,J,K)*G%dx_Cv(i,J))*Slope) @@ -1253,7 +1253,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV enddo ! k-loop if (CS%use_FGNV_streamfn) then - do k=1,nz ; do i=is,ie ; if (G%mask2dCv(i,J)>0.) then + do k=1,nz ; do i=is,ie ; if (G%OBCmaskCv(i,J)>0.) then h_harm = max( h_neglect, & 2. * h(i,j,k) * h(i,j+1,k) / ( ( h(i,j,k) + h(i,j+1,k) ) + h_neglect ) ) c2_h_v(i,k) = CS%FGNV_scale * & @@ -1262,7 +1262,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! Solve an elliptic equation for the streamfunction following Ferrari et al., 2010. do i=is,ie - if (G%mask2dCv(i,J)>0.) then + if (G%OBCmaskCv(i,J)>0.) then do K=2,nz Sfn_unlim_v(i,K) = (1. + CS%FGNV_scale) * Sfn_unlim_v(i,K) enddo @@ -1651,7 +1651,7 @@ subroutine add_detangling_Kh(h, e, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, tv, dt, G, GV de_bot(i,j) = de_bot(i,j) + h(i,j,k+1) enddo ; enddo - do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0.0) then + do j=js,je ; do I=is-1,ie ; if (G%OBCmaskCu(I,j) > 0.0) then if (h(i,j,k) > h(i+1,j,k)) then h2 = h(i,j,k) h1 = max( h(i+1,j,k), h2 - min(de_bot(i+1,j), de_top(i+1,j,k)) ) @@ -1663,7 +1663,7 @@ subroutine add_detangling_Kh(h, e, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, tv, dt, G, GV KH_lay_u(I,j,k) = (Kh_scale * KH_u_CFL(I,j)) * jag_Rat**2 endif ; enddo ; enddo - do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.0) then + do J=js-1,je ; do i=is,ie ; if (G%OBCmaskCv(i,J) > 0.0) then if (h(i,j,k) > h(i,j+1,k)) then h2 = h(i,j,k) h1 = max( h(i,j+1,k), h2 - min(de_bot(i,j+1), de_top(i,j+1,k)) ) @@ -1689,7 +1689,7 @@ subroutine add_detangling_Kh(h, e, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, tv, dt, G, GV ! First, populate the diffusivities if (n==1) then ! This is a u-column. do i=ish,ie - do_i(I) = (G%mask2dCu(I,j) > 0.0) + do_i(I) = (G%OBCmaskCu(I,j) > 0.0) Kh_Max_max(I) = KH_u_CFL(I,j) enddo do K=1,nz+1 ; do i=ish,ie @@ -1699,7 +1699,7 @@ subroutine add_detangling_Kh(h, e, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, tv, dt, G, GV enddo ; enddo else ! This is a v-column. do i=ish,ie - do_i(i) = (G%mask2dCv(i,J) > 0.0) ; Kh_Max_max(I) = KH_v_CFL(i,J) + do_i(i) = (G%OBCmaskCv(i,J) > 0.0) ; Kh_Max_max(I) = KH_v_CFL(i,J) enddo do K=1,nz+1 ; do i=ish,ie Kh_bg(I,K) = KH_v(I,j,K) ; Kh(I,K) = Kh_bg(I,K) @@ -2003,11 +2003,11 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) allocate(CS%Kh_eta_v(G%isd:G%ied, G%JsdB:G%JedB), source=0.) do j=G%jsc,G%jec ; do I=G%isc-1,G%iec grid_sp = sqrt((2.0*G%dxCu(I,j)**2 * G%dyCu(I,j)**2) / (G%dxCu(I,j)**2 + G%dyCu(I,j)**2)) - CS%Kh_eta_u(I,j) = G%mask2dCu(I,j) * MAX(0.0, CS%Kh_eta_bg + CS%Kh_eta_vel * grid_sp) + CS%Kh_eta_u(I,j) = G%OBCmaskCu(I,j) * MAX(0.0, CS%Kh_eta_bg + CS%Kh_eta_vel * grid_sp) enddo ; enddo do J=G%jsc-1,G%jec ; do i=G%isc,G%iec grid_sp = sqrt((2.0*G%dxCv(i,J)**2 * G%dyCv(i,J)**2) / (G%dxCv(i,J)**2 + G%dyCv(i,J)**2)) - CS%Kh_eta_v(i,J) = G%mask2dCv(i,J) * MAX(0.0, CS%Kh_eta_bg + CS%Kh_eta_vel * grid_sp) + CS%Kh_eta_v(i,J) = G%OBCmaskCv(i,J) * MAX(0.0, CS%Kh_eta_bg + CS%Kh_eta_vel * grid_sp) enddo ; enddo endif From 70a9e1521419461aafa12e95f6b845ab36bc29b7 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 13 Oct 2022 06:10:51 -0400 Subject: [PATCH 4/4] +Eliminate unused OBC arguments Eliminated OBC arguments that are no longer used by three internal routines in MOM_lateral_mixing_coeffs. All answers are bitwise identical. --- .../lateral/MOM_lateral_mixing_coeffs.F90 | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index e6dd57fc2e..87562a9c83 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -17,7 +17,7 @@ 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, OBC_NONE +use MOM_open_boundary, only : ocean_OBC_type implicit none ; private @@ -475,16 +475,16 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC) 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) - call calc_Eady_growth_rate_2D(CS, G, GV, US, OBC, h, e, dzu, dzv, dzSxN, dzSyN, CS%SN_u, CS%SN_v) + call calc_Eady_growth_rate_2D(CS, G, GV, US, h, e, dzu, dzv, dzSxN, dzSyN, CS%SN_u, CS%SN_v) else call find_eta(h, tv, G, GV, US, e, halo_size=2) if (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, OBC) + call calc_Visbeck_coeffs_old(h, CS%slope_x, CS%slope_y, N2_u, N2_v, G, GV, US, CS) else !call calc_isoneutral_slopes(G, GV, h, e, tv, dt*CS%kappa_smooth, CS%slope_x, CS%slope_y) - call calc_slope_functions_using_just_e(h, G, GV, US, CS, e, .true., OBC) + call calc_slope_functions_using_just_e(h, G, GV, US, CS, e, .true.) endif endif endif @@ -507,7 +507,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, OBC) +subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, CS) 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] @@ -519,7 +519,6 @@ 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 struct - type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. ! Local variables real :: S2 ! Interface slope squared [nondim] @@ -647,12 +646,11 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C end subroutine calc_Visbeck_coeffs_old !> Calculates the Eady growth rate (2D fields) for use in MEKE and the Visbeck schemes -subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, OBC, h, e, dzu, dzv, dzSxN, dzSyN, SN_u, SN_v) +subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, h, e, dzu, dzv, dzSxN, dzSyN, SN_u, SN_v) type(VarMix_CS), intent(inout) :: CS !< Variable mixing coefficients 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 - type(ocean_OBC_type), pointer, intent(in) :: OBC !< Open boundaries control structure. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Interface height [Z ~> m] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface height [Z ~> m] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: dzu !< dz at u-points [Z ~> m] @@ -802,7 +800,7 @@ end subroutine calc_Eady_growth_rate_2D !> The original calc_slope_function() that calculated slopes using !! interface positions only, not accounting for density variations. -subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slopes, OBC) +subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slopes) 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(inout) :: h !< Layer thickness [H ~> m or kg m-2] @@ -811,7 +809,6 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface position [Z ~> m] logical, intent(in) :: calculate_slopes !< If true, calculate slopes !! internally otherwise use slopes stored in CS - type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. ! Local variables real :: E_x(SZIB_(G), SZJ_(G)) ! X-slope of interface at u points [nondim] (for diagnostics) real :: E_y(SZI_(G), SZJB_(G)) ! Y-slope of interface at v points [nondim] (for diagnostics)