From 8574c7c6a503b6dac7f50af2fe5ecd49d9f05eff Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 9 Oct 2022 23:14:38 -0400 Subject: [PATCH] *+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