diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 41deef61a7..b73590f8e6 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -384,7 +384,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, if (CS%use_GME) then ! GME tapers off above this depth H0 = 1000.0 - FWfrac = 0.1 + FWfrac = 1.0 ! initialize diag. array with zeros GME_coeff_h(:,:,:) = 0.0 GME_coeff_q(:,:,:) = 0.0 @@ -447,10 +447,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, ! Calculate horizontal tension do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 - sh_xx(i,j) = (CS%DY_dxT(i,j)*(G%IdyCu(I,j) * u(I,j,k) - & - G%IdyCu(I-1,j) * u(I-1,j,k)) - & - CS%DX_dyT(i,j)*(G%IdxCv(i,J) * v(i,J,k) - & - G%IdxCv(i,J-1)*v(i,J-1,k))) + dudx(i,j) = CS%DY_dxT(i,j)*(G%IdyCu(I,j) * u(I,j,k) - & + G%IdyCu(I-1,j) * u(I-1,j,k)) + dvdy(i,j) = CS%DX_dyT(i,j)*(G%IdxCv(i,J) * v(i,J,k) - & + G%IdxCv(i,J-1) * v(i,J-1,k)) + sh_xx(i,j) = dudx(i,j) - dvdy(i,j) enddo ; enddo ! Components for the shearing strain @@ -459,6 +460,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, dudy(I,J) = CS%DX_dyBu(I,J)*(u(I,j+1,k)*G%IdxCu(I,j+1) - u(I,j,k)*G%IdxCu(I,j)) enddo ; enddo + ! Interpolate the thicknesses to velocity points. ! The extra wide halos are to accommodate the cross-corner-point projections ! in OBCs, which are not ordinarily be necessary, and might not be necessary @@ -676,12 +678,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, (G%dxCv(i,J) * v(i,J,k) * (h(i,j,k)+h(i,j+1,k)) - & G%dxCv(i,J-1)*v(i,J-1,k)*(h(i,j,k)+h(i,j-1,k))))*G%IareaT(i,j)/ & (h(i,j,k) + GV%H_subroundoff) - dudx(i,j) = 0.5*(G%dyCu(I,j) * u(I,j,k) * (h(i+1,j,k)+h(i,j,k)) - & - G%dyCu(I-1,j) * u(I-1,j,k) * (h(i-1,j,k)+h(i,j,k)) )*G%IareaT(i,j)/ & - (h(i,j,k) + GV%H_subroundoff) * G%mask2dcu(I,j) - dvdy(i,j) = 0.5*(G%dxCv(i,J) * v(i,J,k) * (h(i,j,k)+h(i,j+1,k)) - & - G%dxCv(i,J-1)*v(i,J-1,k)*(h(i,j,k)+h(i,j-1,k)))*G%IareaT(i,j)/ & - (h(i,j,k) + GV%H_subroundoff) * G%mask2dcv(i,J) enddo ; enddo call pass_var(div_xx, G%Domain, complete=.true.) @@ -1031,17 +1027,21 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, if (find_FrictWork) then if (CS%biharmonic) call pass_vector(u0, v0, G%Domain) + call pass_var(dudx, G%Domain, complete=.true.) + call pass_var(dvdy, G%Domain, complete=.true.) + call pass_var(dvdx, G%Domain, position=CORNER, complete=.true.) + call pass_var(dudy, G%Domain, position=CORNER, complete=.true.) do j=js,je ; do i=is,ie ! Diagnose str_xx*d_x u - str_yy*d_y v + str_xy*(d_y u + d_x v) - FrictWork_diss(i,j,k) = Kh_h(i,j,k) * (dudx(i,j)**2 + dvdy(i,j)**2 + & + FrictWork_diss(i,j,k) = -Kh_h(i,j,k) * (dudx(i,j)**2 + dvdy(i,j)**2 + & (0.25*(dvdx(I,J)+dvdx(I-1,J)+dvdx(I,J-1)+dvdx(I-1,J-1)) )**2 + & - (0.25*(dudy(I,J)+dudy(I-1,J)+dudy(I,J-1)+dudy(I-1,J-1)) )**2) + & + (0.25*(dudy(I,J)+dudy(I-1,J)+dudy(I,J-1)+dudy(I-1,J-1)) )**2) - & Ah_h(i,j,k) * ((0.5*(u0(I,j) + u0(I-1,j)))**2 + & (0.5*(v0(i,J) + v0(i,J-1)))**2) if (associated(MEKE)) then ; if (associated(MEKE%mom_src)) then - FrictWorkMax(i,j,k) = MEKE%MEKE(i,j) * sqrt(dudx(i,j)**2 + dvdy(i,j)**2 + & + FrictWorkMax(i,j,k) = 2.0*MEKE%MEKE(i,j) * sqrt(dudx(i,j)**2 + dvdy(i,j)**2 + & (0.25*(dvdx(I,J)+dvdx(I-1,J)+dvdx(I,J-1)+dvdx(I-1,J-1)) )**2 + & (0.25*(dudy(I,J)+dudy(I-1,J)+dudy(I,J-1)+dudy(I-1,J-1)) )**2) @@ -1050,7 +1050,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, ! we need to add the FrictWork done by the dissipation operators, since this work ! is done only for numerical stability and is therefore spurious if (CS%use_GME) then - target_FrictWork_GME(i,j,k) = FWfrac * FrictWorkMax(i,j,k) - FrictWork(i,j,k) + target_FrictWork_GME(i,j,k) = FWfrac * FrictWorkMax(i,j,k) - FrictWork_diss(i,j,k) endif endif ; endif diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index b6e6caea8f..cabcf1f9f2 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -329,16 +329,6 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, MEKE, VarMix, CDp, CS do K=1,nz+1 ; do J=js-1,je ; do i=is,ie ; int_slope_v(i,J,K) = 0.0 ; enddo ; enddo ; enddo !$OMP end parallel -! if (associated(MEKE)) then -! do j=js,je ; do i=is,ie -! MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + 0.25*(Kh_u(I,j,k) + KH_u(I-1,j,k) + & -! KH_v(i,J,k) + KH_v(i,J-1,k)) -!! 0.25*MAX(VarMix%N2_u(I,j,k)+VarMix%N2_u(I-1,j,k)+VarMix%N2_v(i,J,k)+VarMix%N2_v(i,J-1,k),0.0) * & -!! ( (0.5*(VarMix%slope_x(I,j,k)+VarMix%slope_x(I-1,j,k)) )**2 + & -!! (0.5*(VarMix%slope_y(i,J,k)+VarMix%slope_y(i,J-1,k)) )**2 ) -! enddo; enddo -! endif - if (CS%detangle_interfaces) then call add_detangling_Kh(h, e, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, tv, dt, G, GV, & CS, int_slope_u, int_slope_v) @@ -362,8 +352,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, MEKE, VarMix, CDp, CS ! Calculate uhD, vhD from h, e, KH_u, KH_v, tv%T/S if (use_stored_slopes) then call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, MEKE, CS, & - int_slope_u, int_slope_v, VarMix%slope_x, VarMix%slope_y, & - VarMix%N2_u, VarMix%N2_v) + int_slope_u, int_slope_v, VarMix%slope_x, VarMix%slope_y) else call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, MEKE, CS, & int_slope_u, int_slope_v) @@ -457,7 +446,7 @@ end subroutine thickness_diffuse !! Fluxes are limited to give positive definite thicknesses. !! Called by thickness_diffuse(). subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, MEKE, & - CS, int_slope_u, int_slope_v, slope_x, slope_y, N2_x, N2_y) + CS, int_slope_u, int_slope_v, slope_x, slope_y) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness in H (m or kg/m2) @@ -483,10 +472,6 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV !! density gradients. real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), optional, intent(in) :: slope_x !< Isopycnal slope at u-points real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), optional, intent(in) :: slope_y !< Isopycnal slope at v-points - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), optional, intent(in) :: N2_x !< Brunt-Vaisala frequency at - !! u points (s-2) - real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), optional, intent(in) :: N2_y !< Brunt-Vaisala frequency at - !! v points (s-2) ! Local variables real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: & T, & ! The temperature (or density) in C, with the values in @@ -571,7 +556,6 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV real, dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: diag_sfn_y, diag_sfn_unlim_y ! Diagnostics logical :: present_int_slope_u, present_int_slope_v logical :: present_slope_x, present_slope_y, calc_derivatives - logical :: present_N2_x, present_N2_y integer :: is, ie, js, je, nz, IsdB integer :: i, j, k is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke ; IsdB = G%IsdB @@ -589,8 +573,6 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV present_int_slope_v = PRESENT(int_slope_v) present_slope_x = PRESENT(slope_x) present_slope_y = PRESENT(slope_y) - present_N2_x = PRESENT(N2_x) - present_N2_y = PRESENT(N2_y) nk_linear = max(GV%nkml, 1) @@ -1198,14 +1180,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ((Work_u(I-1,j) + Work_u(I,j)) + (Work_v(i,J-1) + Work_v(i,J))) if (associated(CS%GMwork)) CS%GMwork(i,j) = Work_h if (associated(MEKE)) then ; if (associated(MEKE%GM_src)) then - MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + Work_h -! if (present_N2_x .and. present_N2_y) & -! MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + 0.25*(Kh_u(I,j,k) + KH_u(I-1,j,k) + & -! KH_v(i,J,k) + KH_v(i,J-1,k)) * & -! 0.25*MAX(N2_x(I,j,k)+N2_x(I-1,j,k)+N2_y(i,J,k)+N2_y(i,J-1,k),0.0) * & -! ( (0.5*(slope_x(I,j,k)+slope_x(I-1,j,k)) )**2 + & -! (0.5*(slope_y(i,J,k)+slope_y(i,J-1,k)) )**2 ) - + MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + MIN(0.0,Work_h) endif ; endif enddo ; enddo ; endif