From 58f6ffaef79758be4998a774fbf393befd5a2431 Mon Sep 17 00:00:00 2001 From: Scott Bachman Date: Mon, 5 Jun 2017 13:15:56 -0600 Subject: [PATCH 1/6] This adds the 2D Leith viscosity. Viscosity is proportional to the magnitude of the vertical vorticity gradient. This commit adds this viscosity to MOM_hor_visc.F90 by mimicking the Smagorinsky code. --- .../lateral/MOM_hor_visc.F90 | 198 ++++++++++++++---- 1 file changed, 161 insertions(+), 37 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 767f09ec37..073d35f544 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -131,6 +131,10 @@ module MOM_hor_visc ! viscosity. KH is the background value. logical :: Smagorinsky_Ah ! If true, use a biharmonic form of Smagorinsky ! nonlinear eddy viscosity. AH is the background. + logical :: Leith_Kh ! If true, use 2D Leith nonlinear eddy + ! viscosity. KH is the background value. + logical :: Leith_Ah ! If true, use a biharmonic form of 2D Leith + ! nonlinear eddy viscosity. AH is the background. logical :: bound_Coriolis ! If true & SMAGORINSKY_AH is used, the biharmonic ! viscosity is modified to include a term that ! scales quadratically with the velocity shears. @@ -141,13 +145,13 @@ module MOM_hor_visc real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & Kh_bg_xx, &! The background Laplacian viscosity at h points, in units ! of m2 s-1. The actual viscosity may be the larger of this - ! viscosity and the Smagorinsky viscosity. + ! viscosity and the Smagorinsky and Leith viscosities. Kh_bg_2d, &! The background Laplacian viscosity at h points, in units ! of m2 s-1. The actual viscosity may be the larger of this - ! viscosity and the Smagorinsky viscosity. + ! viscosity and the Smagorinsky and Leith viscosities. Ah_bg_xx, &! The background biharmonic viscosity at h points, in units ! of m4 s-1. The actual viscosity may be the larger of this - ! viscosity and the Smagorinsky viscosity. + ! viscosity and the Smagorinsky and Leith viscosities. Kh_Max_xx, &! The maximum permitted Laplacian viscosity, m2 s-1. Ah_Max_xx, &! The maximum permitted biharmonic viscosity, m4 s-1. Biharm_Const2_xx,&! A constant relating the biharmonic viscosity to the @@ -161,10 +165,10 @@ module MOM_hor_visc real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: & Kh_bg_xy, &! The background Laplacian viscosity at q points, in units ! of m2 s-1. The actual viscosity may be the larger of this - ! viscosity and the Smagorinsky viscosity. + ! viscosity and the Smagorinsky and Leith viscosities. Ah_bg_xy, &! The background biharmonic viscosity at q points, in units ! of m4 s-1. The actual viscosity may be the larger of this - ! viscosity and the Smagorinsky viscosity. + ! viscosity and the Smagorinsky and Leith viscosities. Kh_Max_xy, &! The maximum permitted Laplacian viscosity, m2 s-1. Ah_Max_xy, &! The maximum permitted biharmonic viscosity, m4 s-1. Biharm_Const2_xy,&! A constant relating the biharmonic viscosity to the @@ -189,11 +193,16 @@ module MOM_hor_visc ! The following variables are precalculated time-invariant combinations of ! parameters and metric terms. real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & - Laplac_Const_xx, & ! Laplacian metric-dependent constants (nondim) - Biharm_Const_xx ! Biharmonic metric-dependent constants (nondim) + Laplac_Const_xx, & ! Laplacian metric-dependent constants (nondim) + Biharm_Const_xx, & ! Biharmonic metric-dependent constants (nondim) + Laplac3_Const_xx, & ! Laplacian metric-dependent constants (nondim) + Biharm5_Const_xx ! Biharmonic metric-dependent constants (nondim) + real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: & - Laplac_Const_xy, & ! Laplacian metric-dependent constants (nondim) - Biharm_Const_xy ! Biharmonic metric-dependent constants (nondim) + Laplac_Const_xy, & ! Laplacian metric-dependent constants (nondim) + Biharm_Const_xy, & ! Biharmonic metric-dependent constants (nondim) + Laplac3_Const_xy, & ! Laplacian metric-dependent constants (nondim) + Biharm5_Const_xy ! Biharmonic metric-dependent constants (nondim) type(diag_ctrl), pointer :: diag ! structure to regulate diagnostic timing @@ -269,6 +278,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, str_xy, & ! str_xy is the cross term in the stress tensor (H m2 s-2) bhstr_xy ! A copy of str_xy that only contains the biharmonic contribution (H m2 s-2) + real, dimension(SZIB_(G),SZJB_(G)) :: & + vort_xy, & ! vertical vorticity (dv/dx - du/dy) (1/sec) including metric terms + vort_xy_dx, & ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) (m-1 sec-1) including metric terms + vort_xy_dy ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) (m-1 sec-1) including metric terms + real, dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: & Ah_q, & ! biharmonic viscosity at corner points (m4/s) Kh_q ! Laplacian viscosity at corner points (m2/s) @@ -282,7 +296,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, real :: Kh ! Laplacian viscosity (m2/s) real :: AhSm ! Smagorinsky biharmonic viscosity (m4/s) real :: KhSm ! Smagorinsky Laplacian viscosity (m2/s) + real :: AhLth ! 2D Leith biharmonic viscosity (m4/s) + real :: KhLth ! 2D Leith Laplacian viscosity (m2/s) real :: Shear_mag ! magnitude of the shear (1/s) + real :: Vort_mag ! magnitude of the vorticity (1/s) real :: h2uq, h2vq ! temporary variables in units of H^2 (i.e. m2 or kg2 m-4). real :: hu, hv ! Thicknesses interpolated by arithmetic means to corner ! points; these are first interpolated to u or v velocity @@ -480,6 +497,22 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, enddo ; enddo endif + if (CS%no_slip) then + do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 + vort_xy(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx(I,J) - dudy(I,J) ) + enddo ; enddo + else + do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 + vort_xy(I,J) = G%mask2dBu(I,J) * ( dvdx(I,J) - dudy(I,J) ) + enddo ; enddo + endif + +! Vorticity gradient + do J=js-1,Jeq+1 ; do I=is-1,Ieq+1 + vort_xy_dx(I,J) = CS%DY_dxBu(I,J)*(vort_xy(I,J)*G%IdyCu(I,J) - vort_xy(I-1,J)*G%IdyCu(I-1,J)) + vort_xy_dy(I,J) = CS%DX_dyBu(I,J)*(vort_xy(I,J)*G%IdxCv(I,J) - vort_xy(I,J-1)*G%IdxCv(I,J-1)) + enddo ; enddo + ! Evaluate u0 = x.Div(Grad u) and v0 = y.Div( Grad u) if (CS%biharmonic) then do j=js-1,Jeq+1 ; do I=Isq-1,Ieq+1 @@ -507,10 +540,18 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, endif do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) & + if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then Shear_mag = sqrt(sh_xx(i,j)*sh_xx(i,j) + & 0.25*((sh_xy(I-1,J-1)*sh_xy(I-1,J-1) + sh_xy(I,J)*sh_xy(I,J)) + & (sh_xy(I-1,J)*sh_xy(I-1,J) + sh_xy(I,J-1)*sh_xy(I,J-1)))) +! Vort_mag = sqrt( & +! 0.5*((vort_xy_dx(I,J-1)*vort_xy_dx(I,J-1) + vort_xy_dx(I,J)*vort_xy_dx(I,J)) + & +! (vort_xy_dy(I-1,J)*vort_xy_dy(I-1,J) + vort_xy_dy(I,J)*vort_xy_dy(I,J)))) + endif + if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) & + Vort_mag = sqrt( & + 0.5*((vort_xy_dx(I,J-1)*vort_xy_dx(I,J-1) + vort_xy_dx(I,J)*vort_xy_dx(I,J)) + & + (vort_xy_dy(I-1,J)*vort_xy_dy(I-1,J) + vort_xy_dy(I,J)*vort_xy_dy(I,J)))) if (CS%better_bound_Ah .or. CS%better_bound_Kh) then hrat_min = min(1.0, min(h_u(I,j), h_u(I-1,j), h_v(i,J), h_v(i,J-1)) / & (h(i,j,k) + h_neglect) ) @@ -521,14 +562,19 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, ! Determine the Laplacian viscosity at h points, using the ! largest value from several parameterizations. Kh_scale = 1.0 ; if (rescale_Kh) Kh_scale = VarMix%Res_fn_h(i,j) - if (CS%Smagorinsky_Kh) then - KhSm = CS%LAPLAC_CONST_xx(i,j) * Shear_mag - Kh = Kh_scale * MAX(CS%Kh_bg_xx(i,j), KhSm) + KhSm = 0.0; KhLth = 0.0 + if ((CS%Smagorinsky_Kh) .or. (CS%Leith_Kh)) then + if (CS%Smagorinsky_Kh) & + KhSm = CS%LAPLAC_CONST_xx(i,j) * Shear_mag + if (CS%Leith_Kh) & + KhLth = CS%LAPLAC3_CONST_xx(i,j) * Vort_mag + Kh = Kh_scale * MAX(KhLth, MAX(CS%Kh_bg_xx(i,j), KhSm)) if (CS%bound_Kh .and. .not.CS%better_bound_Kh) & Kh = MIN(Kh, CS%Kh_Max_xx(i,j)) else Kh = Kh_scale * CS%Kh_bg_xx(i,j) endif + if (use_MEKE_Ku) then Kh = Kh + MEKE%Ku(i,j) endif @@ -552,19 +598,25 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, if (CS%biharmonic) then ! Determine the biharmonic viscosity at h points, using the ! largest value from several parameterizations. - if (CS%Smagorinsky_Ah) then - if (CS%bound_Coriolis) then - AhSm = Shear_mag * (CS%BIHARM_CONST_xx(i,j) + & + AhSm = 0.0; AhLth = 0.0 + if ((CS%Smagorinsky_Ah) .or. (CS%Leith_Ah)) then + if (CS%Smagorinsky_Ah) then + if (CS%bound_Coriolis) then + AhSm = Shear_mag * (CS%BIHARM_CONST_xx(i,j) + & CS%Biharm_Const2_xx(i,j)*Shear_mag) - else - AhSm = CS%BIHARM_CONST_xx(i,j) * Shear_mag + else + AhSm = CS%BIHARM_CONST_xx(i,j) * Shear_mag + endif endif - Ah = MAX(CS%Ah_bg_xx(i,j), AhSm) + if (CS%Leith_Ah) & + AhLth = Vort_mag * (CS%BIHARM_CONST_xx(i,j)) + Ah = MAX(MAX(CS%Ah_bg_xx(i,j), AhSm),AhLth) if (CS%bound_Ah .and. .not.CS%better_bound_Ah) & Ah = MIN(Ah, CS%Ah_Max_xx(i,j)) - else - Ah = CS%Ah_bg_xx(i,j) - endif ! Smagorinsky_Ah + else + Ah = CS%Ah_bg_xx(i,j) + endif ! Smagorinsky_Ah or Leith_Ah + if (CS%better_bound_Ah) then Ah = MIN(Ah, visc_bound_rem*hrat_min*CS%Ah_Max_xx(i,j)) endif @@ -618,10 +670,18 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, endif do J=js-1,Jeq ; do I=is-1,Ieq - if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) & + if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then Shear_mag = sqrt(sh_xy(I,J)*sh_xy(I,J) + & 0.25*((sh_xx(i,j)*sh_xx(i,j) + sh_xx(i+1,j+1)*sh_xx(i+1,j+1)) + & (sh_xx(i,j+1)*sh_xx(i,j+1) + sh_xx(i+1,j)*sh_xx(i+1,j)))) +! Vort_mag = sqrt( & +! 0.5*((vort_xy_dx(I,J)*vort_xy_dx(I,J) + vort_xy_dx(I+1,J)*vort_xy_dx(I+1,J)) + & +! (vort_xy_dy(I,J)*vort_xy_dy(I,J) + vort_xy_dy(I,J+1)*vort_xy_dy(I,J+1)))) + endif + if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) & + Vort_mag = sqrt( & + 0.5*((vort_xy_dx(I,J)*vort_xy_dx(I,J) + vort_xy_dx(I+1,J)*vort_xy_dx(I+1,J)) + & + (vort_xy_dy(I,J)*vort_xy_dy(I,J) + vort_xy_dy(I,J+1)*vort_xy_dy(I,J+1)))) h2uq = 4.0 * h_u(I,j) * h_u(I,j+1) h2vq = 4.0 * h_v(i,J) * h_v(i+1,J) @@ -660,9 +720,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, ! Determine the Laplacian viscosity at q points, using the ! largest value from several parameterizations. Kh_scale = 1.0 ; if (rescale_Kh) Kh_scale = VarMix%Res_fn_q(I,J) - if (CS%Smagorinsky_Kh) then - KhSm = CS%LAPLAC_CONST_xy(I,J) * Shear_mag - Kh = Kh_scale * MAX(CS%Kh_bg_xy(I,J), KhSm) + KhSm = 0.0; KhLth = 0.0 + if ((CS%Smagorinsky_Kh) .or. (CS%Leith_Kh)) then + if (CS%Smagorinsky_Kh) & + KhSm = CS%LAPLAC_CONST_xy(I,J) * Shear_mag + if (CS%Leith_Kh) & + KhLth = CS%LAPLAC3_CONST_xy(I,J) * Vort_mag + Kh = Kh_scale * MAX(MAX(CS%Kh_bg_xy(I,J), KhSm), KhLth) if (CS%bound_Kh .and. .not.CS%better_bound_Kh) & Kh = MIN(Kh, CS%Kh_Max_xy(I,J)) else @@ -695,19 +759,24 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, if (CS%biharmonic) then ! Determine the biharmonic viscosity at q points, using the ! largest value from several parameterizations. - if (CS%Smagorinsky_Ah) then - if (CS%bound_Coriolis) then - AhSm = Shear_mag * (CS%BIHARM_CONST_xy(I,J) + & + AhSm = 0.0; AhLth = 0.0 + if (CS%Smagorinsky_Ah .or. CS%Leith_Ah) then + if (CS%Smagorinsky_Ah) then + if (CS%bound_Coriolis) then + AhSm = Shear_mag * (CS%BIHARM_CONST_xy(I,J) + & CS%Biharm_Const2_xy(I,J)*Shear_mag) - else - AhSm = CS%BIHARM_CONST_xy(I,J) * Shear_mag + else + AhSm = CS%BIHARM_CONST_xy(I,J) * Shear_mag + endif endif - Ah = MAX(CS%Ah_bg_xy(I,J), AhSm) + if (CS%Leith_Ah) & + AhLth = Vort_mag * (CS%BIHARM5_CONST_xy(I,J)) + Ah = MAX(MAX(CS%Ah_bg_xy(I,J), AhSm),AhLth) if (CS%bound_Ah .and. .not.CS%better_bound_Ah) & Ah = MIN(Ah, CS%Ah_Max_xy(I,J)) else Ah = CS%Ah_bg_xy(I,J) - endif ! Smagorinsky_Ah + endif ! Smagorinsky_Ah or Leith_Ah if (CS%better_bound_Ah) then Ah = MIN(Ah, visc_bound_rem*hrat_min*CS%Ah_Max_xy(I,J)) endif @@ -886,7 +955,9 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) ! u0v is the Laplacian sensitivities to the v velocities ! at u points, in m-2, with u0u, v0u, and v0v defined similarly. real :: grid_sp_h2 ! Harmonic mean of the squares of the grid + real :: grid_sp_h3 ! Harmonic mean of the squares of the grid^(3/2) real :: grid_sp_q2 ! spacings at h and q points (m2) + real :: grid_sp_q3 ! spacings at h and q points^(3/2) (m3) real :: Kh_Limit ! A coefficient (1/s) used, along with the ! grid spacing, to limit Laplacian viscosity. real :: fmax ! maximum absolute value of f at the four @@ -900,6 +971,8 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) real :: Ah_vel_scale ! this speed (m/s) times grid spacing cubed gives bih visc real :: Smag_Lap_const ! nondimensional Laplacian Smagorinsky constant real :: Smag_bi_const ! nondimensional biharmonic Smagorinsky constant + real :: Leith_Lap_const ! nondimensional Laplacian Leith constant + real :: Leith_bi_const ! nondimensional biharmonic Leith constant real :: dt ! dynamics time step (sec) real :: Idt ! inverse of dt (1/s) real :: denom ! work variable; the denominator of a fraction @@ -939,8 +1012,8 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) ! It is not clear whether these initialization lines are needed for the ! cases where the corresponding parameters are not read. - CS%bound_Kh = .false. ; CS%better_bound_Kh = .false. ; CS%Smagorinsky_Kh = .false. - CS%bound_Ah = .false. ; CS%better_bound_Ah = .false. ; CS%Smagorinsky_Ah = .false. + CS%bound_Kh = .false. ; CS%better_bound_Kh = .false. ; CS%Smagorinsky_Kh = .false. ; CS%Leith_Kh = .false. + CS%bound_Ah = .false. ; CS%better_bound_Ah = .false. ; CS%Smagorinsky_Ah = .false. ; CS%Leith_Ah = .false. CS%bound_Coriolis = .false. Kh = 0.0 ; Ah = 0.0 @@ -963,7 +1036,7 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) "The velocity scale which is multiplied by the grid \n"//& "spacing to calculate the Laplacian viscosity. \n"//& "The final viscosity is the largest of this scaled \n"//& - "viscosity, the Smagorinsky viscosity and KH.", & + "viscosity, the Smagorinsky and Leith viscosities, and KH.", & units="m s-1", default=0.0) call get_param(param_file, mod, "SMAGORINSKY_KH", CS%Smagorinsky_Kh, & @@ -975,6 +1048,15 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) "often 0.15.", units="nondim", default=0.0, & fail_if_missing = CS%Smagorinsky_Kh) + call get_param(param_file, mod, "LEITH_KH", CS%Leith_Kh, & + "If true, use a Leith nonlinear eddy viscosity.", & + default=.false.) + if (CS%Leith_Kh .or. get_all) & + call get_param(param_file, mod, "LEITH_LAP_CONST", Leith_Lap_const, & + "The nondimensional Laplacian Leith constant, \n"//& + "often ??", units="nondim", default=0.0, & + fail_if_missing = CS%Leith_Kh) + call get_param(param_file, mod, "BOUND_KH", CS%bound_Kh, & "If true, the Laplacian coefficient is locally limited \n"//& "to be stable.", default=.true.) @@ -996,11 +1078,14 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) "The velocity scale which is multiplied by the cube of \n"//& "the grid spacing to calculate the biharmonic viscosity. \n"//& "The final viscosity is the largest of this scaled \n"//& - "viscosity, the Smagorinsky viscosity and AH.", & + "viscosity, the Smagorinsky and Leith viscosities, and AH.", & units="m s-1", default=0.0) call get_param(param_file, mod, "SMAGORINSKY_AH", CS%Smagorinsky_Ah, & "If true, use a biharmonic Smagorinsky nonlinear eddy \n"//& "viscosity.", default=.false.) + call get_param(param_file, mod, "LEITH_AH", CS%Leith_Ah, & + "If true, use a biharmonic Leith nonlinear eddy \n"//& + "viscosity.", default=.false.) call get_param(param_file, mod, "BOUND_AH", CS%bound_Ah, & "If true, the biharmonic coefficient is locally limited \n"//& @@ -1034,6 +1119,14 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) units="m s-1", default=maxvel) endif endif + + if (CS%Leith_Ah .or. get_all) then + call get_param(param_file, mod, "LEITH_BI_CONST",Leith_bi_const, & + "The nondimensional biharmonic Leith constant, \n"//& + "typically ??", units="nondim", default=0.0, & + fail_if_missing = CS%Leith_Ah) + endif + endif if (CS%better_bound_Ah .or. CS%better_bound_Kh .or. get_all) & @@ -1090,6 +1183,11 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) ALLOC_(CS%Laplac_Const_xx(isd:ied,jsd:jed)) ; CS%Laplac_Const_xx(:,:) = 0.0 ALLOC_(CS%Laplac_Const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Laplac_Const_xy(:,:) = 0.0 endif + if (CS%Leith_Kh) then + ALLOC_(CS%Laplac3_Const_xx(isd:ied,jsd:jed)) ; CS%Laplac3_Const_xx(:,:) = 0.0 + ALLOC_(CS%Laplac3_Const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Laplac3_Const_xy(:,:) = 0.0 + endif + endif ALLOC_(CS%reduction_xx(isd:ied,jsd:jed)) ; CS%reduction_xx(:,:) = 0.0 ALLOC_(CS%reduction_xy(IsdB:IedB,JsdB:JedB)) ; CS%reduction_xy(:,:) = 0.0 @@ -1126,6 +1224,10 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) ALLOC_(CS%Biharm_Const2_xy(IsdB:IedB,JsdB:JedB)) ; CS%Biharm_Const2_xy(:,:) = 0.0 endif endif + if (CS%Leith_Ah) then + ALLOC_(CS%Biharm5_Const_xx(isd:ied,jsd:jed)) ; CS%Biharm5_Const_xx(:,:) = 0.0 + ALLOC_(CS%Biharm5_Const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Biharm5_Const_xy(:,:) = 0.0 + endif endif do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 @@ -1175,7 +1277,9 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) if (CS%bound_Kh .or. CS%bound_Ah) Kh_Limit = 0.3 / (dt*4.0) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 grid_sp_h2 = (2.0*CS%DX2h(i,j)*CS%DY2h(i,j)) / (CS%DX2h(i,j) + CS%DY2h(i,j)) + grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2) if (CS%Smagorinsky_Kh) CS%LAPLAC_CONST_xx(i,j) = Smag_Lap_const * grid_sp_h2 + if (CS%Leith_Kh) CS%LAPLAC3_CONST_xx(i,j) = Leith_Lap_const * grid_sp_h3 CS%Kh_bg_xx(i,j) = MAX(Kh, Kh_vel_scale * sqrt(grid_sp_h2)) @@ -1188,7 +1292,9 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) enddo ; enddo do J=js-1,Jeq ; do I=is-1,Ieq grid_sp_q2 = (2.0*CS%DX2q(I,J)*CS%DY2q(I,J)) / (CS%DX2q(I,J) + CS%DY2q(I,J)) + grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2) if (CS%Smagorinsky_Kh) CS%LAPLAC_CONST_xy(I,J) = Smag_Lap_const * grid_sp_q2 + if (CS%Leith_Kh) CS%LAPLAC3_CONST_xy(i,j) = Leith_Lap_const * grid_sp_q3 CS%Kh_bg_xy(I,J) = MAX(Kh, Kh_vel_scale * sqrt(grid_sp_q2)) @@ -1220,6 +1326,7 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) BoundCorConst = 1.0 / (5.0*(bound_Cor_vel*bound_Cor_vel)) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 grid_sp_h2 = (2.0*CS%DX2h(i,j)*CS%DY2h(i,j)) / (CS%DX2h(i,j)+CS%DY2h(i,j)) + grid_sp_h3 = grid_sp_h2**1.5 if (CS%Smagorinsky_Ah) then CS%BIHARM_CONST_xx(i,j) = Smag_bi_const * (grid_sp_h2 * grid_sp_h2) @@ -1230,6 +1337,10 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) (fmax * BoundCorConst) endif endif + if (CS%Leith_Ah) then + CS%BIHARM5_CONST_xx(i,j) = Leith_bi_const * (grid_sp_h2 * grid_sp_h3) + endif + CS%Ah_bg_xx(i,j) = MAX(Ah, Ah_vel_scale * grid_sp_h2 * sqrt(grid_sp_h2)) if (CS%bound_Ah .and. .not.CS%better_bound_Ah) then CS%Ah_Max_xx(i,j) = Ah_Limit * (grid_sp_h2 * grid_sp_h2) @@ -1238,6 +1349,8 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) enddo ; enddo do J=js-1,Jeq ; do I=is-1,Ieq grid_sp_q2 = (2.0*CS%DX2q(I,J)*CS%DY2q(I,J)) / (CS%DX2q(I,J)+CS%DY2q(I,J)) + grid_sp_q3 = grid_sp_q2**1.5 + if (CS%Smagorinsky_Ah) then CS%BIHARM_CONST_xy(I,J) = Smag_bi_const * (grid_sp_q2 * grid_sp_q2) if (CS%bound_Coriolis) then @@ -1245,6 +1358,11 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) (abs(G%CoriolisBu(I,J)) * BoundCorConst) endif endif + + if (CS%Leith_Ah) then + CS%BIHARM5_CONST_xy(I,J) = Leith_bi_const * (grid_sp_q2 * grid_sp_q3) + endif + CS%Ah_bg_xy(I,J) = MAX(Ah, Ah_vel_scale * grid_sp_q2 * sqrt(grid_sp_q2)) if (CS%bound_Ah .and. .not.CS%better_bound_Ah) then CS%Ah_Max_xy(I,J) = Ah_Limit * (grid_sp_q2 * grid_sp_q2) @@ -1399,6 +1517,9 @@ subroutine hor_visc_end(CS) if (CS%Smagorinsky_Kh) then DEALLOC_(CS%Laplac_Const_xx) ; DEALLOC_(CS%Laplac_Const_xy) endif + if (CS%Leith_Kh) then + DEALLOC_(CS%Laplac3_Const_xx) ; DEALLOC_(CS%Laplac3_Const_xy) + endif endif if (CS%biharmonic) then @@ -1414,6 +1535,9 @@ subroutine hor_visc_end(CS) DEALLOC_(CS%Biharm_Const2_xx) ; DEALLOC_(CS%Biharm_Const2_xy) endif endif + if (CS%Leith_Ah) then + DEALLOC_(CS%Biharm5_Const_xx) ; DEALLOC_(CS%Biharm5_Const_xy) + endif endif deallocate(CS) From 8e49c02e7295a54a273a95fae3a41987fd8e41b3 Mon Sep 17 00:00:00 2001 From: Scott Bachman Date: Mon, 5 Jun 2017 13:55:37 -0600 Subject: [PATCH 2/6] Fixed whitespace. No details. --- src/parameterizations/lateral/MOM_hor_visc.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 073d35f544..cfaf3a94c2 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -299,7 +299,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, real :: AhLth ! 2D Leith biharmonic viscosity (m4/s) real :: KhLth ! 2D Leith Laplacian viscosity (m2/s) real :: Shear_mag ! magnitude of the shear (1/s) - real :: Vort_mag ! magnitude of the vorticity (1/s) + real :: Vort_mag ! magnitude of the vorticity (1/s) real :: h2uq, h2vq ! temporary variables in units of H^2 (i.e. m2 or kg2 m-4). real :: hu, hv ! Thicknesses interpolated by arithmetic means to corner ! points; these are first interpolated to u or v velocity @@ -551,7 +551,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) & Vort_mag = sqrt( & 0.5*((vort_xy_dx(I,J-1)*vort_xy_dx(I,J-1) + vort_xy_dx(I,J)*vort_xy_dx(I,J)) + & - (vort_xy_dy(I-1,J)*vort_xy_dy(I-1,J) + vort_xy_dy(I,J)*vort_xy_dy(I,J)))) + (vort_xy_dy(I-1,J)*vort_xy_dy(I-1,J) + vort_xy_dy(I,J)*vort_xy_dy(I,J)))) if (CS%better_bound_Ah .or. CS%better_bound_Kh) then hrat_min = min(1.0, min(h_u(I,j), h_u(I-1,j), h_v(i,J), h_v(i,J-1)) / & (h(i,j,k) + h_neglect) ) @@ -564,7 +564,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, Kh_scale = 1.0 ; if (rescale_Kh) Kh_scale = VarMix%Res_fn_h(i,j) KhSm = 0.0; KhLth = 0.0 if ((CS%Smagorinsky_Kh) .or. (CS%Leith_Kh)) then - if (CS%Smagorinsky_Kh) & + if (CS%Smagorinsky_Kh) & KhSm = CS%LAPLAC_CONST_xx(i,j) * Shear_mag if (CS%Leith_Kh) & KhLth = CS%LAPLAC3_CONST_xx(i,j) * Vort_mag @@ -1125,7 +1125,7 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) "The nondimensional biharmonic Leith constant, \n"//& "typically ??", units="nondim", default=0.0, & fail_if_missing = CS%Leith_Ah) - endif + endif endif From e10afceb23bb26a37c8044c9915098d8457275d3 Mon Sep 17 00:00:00 2001 From: Scott Bachman Date: Mon, 5 Jun 2017 15:58:55 -0600 Subject: [PATCH 3/6] Removed commented lines from adding Leith. Deleted some commented lines, per request. --- src/parameterizations/lateral/MOM_hor_visc.F90 | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index cfaf3a94c2..c48d5bc874 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -544,9 +544,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, Shear_mag = sqrt(sh_xx(i,j)*sh_xx(i,j) + & 0.25*((sh_xy(I-1,J-1)*sh_xy(I-1,J-1) + sh_xy(I,J)*sh_xy(I,J)) + & (sh_xy(I-1,J)*sh_xy(I-1,J) + sh_xy(I,J-1)*sh_xy(I,J-1)))) -! Vort_mag = sqrt( & -! 0.5*((vort_xy_dx(I,J-1)*vort_xy_dx(I,J-1) + vort_xy_dx(I,J)*vort_xy_dx(I,J)) + & -! (vort_xy_dy(I-1,J)*vort_xy_dy(I-1,J) + vort_xy_dy(I,J)*vort_xy_dy(I,J)))) endif if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) & Vort_mag = sqrt( & @@ -674,9 +671,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, Shear_mag = sqrt(sh_xy(I,J)*sh_xy(I,J) + & 0.25*((sh_xx(i,j)*sh_xx(i,j) + sh_xx(i+1,j+1)*sh_xx(i+1,j+1)) + & (sh_xx(i,j+1)*sh_xx(i,j+1) + sh_xx(i+1,j)*sh_xx(i+1,j)))) -! Vort_mag = sqrt( & -! 0.5*((vort_xy_dx(I,J)*vort_xy_dx(I,J) + vort_xy_dx(I+1,J)*vort_xy_dx(I+1,J)) + & -! (vort_xy_dy(I,J)*vort_xy_dy(I,J) + vort_xy_dy(I,J+1)*vort_xy_dy(I,J+1)))) endif if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) & Vort_mag = sqrt( & @@ -1123,7 +1117,7 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) if (CS%Leith_Ah .or. get_all) then call get_param(param_file, mod, "LEITH_BI_CONST",Leith_bi_const, & "The nondimensional biharmonic Leith constant, \n"//& - "typically ??", units="nondim", default=0.0, & + "typical values are thus far undetermined", units="nondim", default=0.0, & fail_if_missing = CS%Leith_Ah) endif From 133490ad6b7aa67c9eda2436e8b59d517927fef3 Mon Sep 17 00:00:00 2001 From: Scott Bachman Date: Mon, 5 Jun 2017 16:47:37 -0600 Subject: [PATCH 4/6] Changed indices and array sizes for Leith terms. Changed all indices pertaining to Leith to adhere to capitalization soft convention. Changed the array sizes for vort_xy_** arrays. Changed **1.5 exponents to sqrt()*(). --- .../lateral/MOM_hor_visc.F90 | 30 ++++++++++--------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index c48d5bc874..aa86ae964d 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -276,12 +276,14 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, dvdx, dudy, & ! components in the shearing strain (s-1) sh_xy, & ! horizontal shearing strain (du/dy + dv/dx) (1/sec) including metric terms str_xy, & ! str_xy is the cross term in the stress tensor (H m2 s-2) - bhstr_xy ! A copy of str_xy that only contains the biharmonic contribution (H m2 s-2) + bhstr_xy, & ! A copy of str_xy that only contains the biharmonic contribution (H m2 s-2) + vort_xy ! vertical vorticity (dv/dx - du/dy) (1/sec) including metric terms + + real, dimension(SZI_(G),SZJB_(G)) :: & + vort_xy_dx ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) (m-1 sec-1) including metric terms - real, dimension(SZIB_(G),SZJB_(G)) :: & - vort_xy, & ! vertical vorticity (dv/dx - du/dy) (1/sec) including metric terms - vort_xy_dx, & ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) (m-1 sec-1) including metric terms - vort_xy_dy ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) (m-1 sec-1) including metric terms + real, dimension(SZIB_(G),SZJ_(G)) :: & + vort_xy_dy ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) (m-1 sec-1) including metric terms real, dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: & Ah_q, & ! biharmonic viscosity at corner points (m4/s) @@ -509,8 +511,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, ! Vorticity gradient do J=js-1,Jeq+1 ; do I=is-1,Ieq+1 - vort_xy_dx(I,J) = CS%DY_dxBu(I,J)*(vort_xy(I,J)*G%IdyCu(I,J) - vort_xy(I-1,J)*G%IdyCu(I-1,J)) - vort_xy_dy(I,J) = CS%DX_dyBu(I,J)*(vort_xy(I,J)*G%IdxCv(I,J) - vort_xy(I,J-1)*G%IdxCv(I,J-1)) + vort_xy_dx(i,J) = CS%DY_dxBu(I,J)*(vort_xy(I,J)*G%IdyCu(I,j) - vort_xy(I-1,J)*G%IdyCu(I-1,j)) + vort_xy_dy(I,j) = CS%DX_dyBu(I,J)*(vort_xy(I,J)*G%IdxCv(i,J) - vort_xy(I,J-1)*G%IdxCv(i,J-1)) enddo ; enddo ! Evaluate u0 = x.Div(Grad u) and v0 = y.Div( Grad u) @@ -547,8 +549,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, endif if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) & Vort_mag = sqrt( & - 0.5*((vort_xy_dx(I,J-1)*vort_xy_dx(I,J-1) + vort_xy_dx(I,J)*vort_xy_dx(I,J)) + & - (vort_xy_dy(I-1,J)*vort_xy_dy(I-1,J) + vort_xy_dy(I,J)*vort_xy_dy(I,J)))) + 0.5*((vort_xy_dx(i,J-1)*vort_xy_dx(i,J-1) + vort_xy_dx(i,J)*vort_xy_dx(i,J)) + & + (vort_xy_dy(I-1,j)*vort_xy_dy(I-1,j) + vort_xy_dy(I,j)*vort_xy_dy(I,j)))) if (CS%better_bound_Ah .or. CS%better_bound_Kh) then hrat_min = min(1.0, min(h_u(I,j), h_u(I-1,j), h_v(i,J), h_v(i,J-1)) / & (h(i,j,k) + h_neglect) ) @@ -674,8 +676,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, endif if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) & Vort_mag = sqrt( & - 0.5*((vort_xy_dx(I,J)*vort_xy_dx(I,J) + vort_xy_dx(I+1,J)*vort_xy_dx(I+1,J)) + & - (vort_xy_dy(I,J)*vort_xy_dy(I,J) + vort_xy_dy(I,J+1)*vort_xy_dy(I,J+1)))) + 0.5*((vort_xy_dx(i,J)*vort_xy_dx(i,J) + vort_xy_dx(i+1,J)*vort_xy_dx(i+1,J)) + & + (vort_xy_dy(I,j)*vort_xy_dy(I,j) + vort_xy_dy(I,j+1)*vort_xy_dy(I,j+1)))) h2uq = 4.0 * h_u(I,j) * h_u(I,j+1) h2vq = 4.0 * h_v(i,J) * h_v(i+1,J) @@ -1288,7 +1290,7 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) grid_sp_q2 = (2.0*CS%DX2q(I,J)*CS%DY2q(I,J)) / (CS%DX2q(I,J) + CS%DY2q(I,J)) grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2) if (CS%Smagorinsky_Kh) CS%LAPLAC_CONST_xy(I,J) = Smag_Lap_const * grid_sp_q2 - if (CS%Leith_Kh) CS%LAPLAC3_CONST_xy(i,j) = Leith_Lap_const * grid_sp_q3 + if (CS%Leith_Kh) CS%LAPLAC3_CONST_xy(I,J) = Leith_Lap_const * grid_sp_q3 CS%Kh_bg_xy(I,J) = MAX(Kh, Kh_vel_scale * sqrt(grid_sp_q2)) @@ -1320,7 +1322,7 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) BoundCorConst = 1.0 / (5.0*(bound_Cor_vel*bound_Cor_vel)) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 grid_sp_h2 = (2.0*CS%DX2h(i,j)*CS%DY2h(i,j)) / (CS%DX2h(i,j)+CS%DY2h(i,j)) - grid_sp_h3 = grid_sp_h2**1.5 + grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2) if (CS%Smagorinsky_Ah) then CS%BIHARM_CONST_xx(i,j) = Smag_bi_const * (grid_sp_h2 * grid_sp_h2) @@ -1343,7 +1345,7 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) enddo ; enddo do J=js-1,Jeq ; do I=is-1,Ieq grid_sp_q2 = (2.0*CS%DX2q(I,J)*CS%DY2q(I,J)) / (CS%DX2q(I,J)+CS%DY2q(I,J)) - grid_sp_q3 = grid_sp_q2**1.5 + grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2) if (CS%Smagorinsky_Ah) then CS%BIHARM_CONST_xy(I,J) = Smag_bi_const * (grid_sp_q2 * grid_sp_q2) From fe5c10f2e4b299e73dd23d92bafc6483ef8e2223 Mon Sep 17 00:00:00 2001 From: Scott Bachman Date: Tue, 6 Jun 2017 10:52:23 -0600 Subject: [PATCH 5/6] Fixed loop ranges for vorticity gradient. Fixed loop ranges for vort_xy_dx and vort_xy_dy. --- src/parameterizations/lateral/MOM_hor_visc.F90 | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index aa86ae964d..7ff6ac4018 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -278,7 +278,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, str_xy, & ! str_xy is the cross term in the stress tensor (H m2 s-2) bhstr_xy, & ! A copy of str_xy that only contains the biharmonic contribution (H m2 s-2) vort_xy ! vertical vorticity (dv/dx - du/dy) (1/sec) including metric terms - + real, dimension(SZI_(G),SZJB_(G)) :: & vort_xy_dx ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) (m-1 sec-1) including metric terms @@ -510,10 +510,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, endif ! Vorticity gradient - do J=js-1,Jeq+1 ; do I=is-1,Ieq+1 + do J=Jsq-1,Jeq+1 ; do I=Isq,Ieq+1 vort_xy_dx(i,J) = CS%DY_dxBu(I,J)*(vort_xy(I,J)*G%IdyCu(I,j) - vort_xy(I-1,J)*G%IdyCu(I-1,j)) + enddo ; enddo + + do J=Jsq,Jeq+1 ; do I=Isq-1,Ieq+1 vort_xy_dy(I,j) = CS%DX_dyBu(I,J)*(vort_xy(I,J)*G%IdxCv(i,J) - vort_xy(I,J-1)*G%IdxCv(i,J-1)) - enddo ; enddo + enddo ; enddo ! Evaluate u0 = x.Div(Grad u) and v0 = y.Div( Grad u) if (CS%biharmonic) then From 03df3e460a81f611557ac863921e37ef0bd24c6d Mon Sep 17 00:00:00 2001 From: Scott Bachman Date: Tue, 6 Jun 2017 15:33:31 -0600 Subject: [PATCH 6/6] Fixed loop ranges for vorticity gradients (again). Fixed an inconsistency with loop range when calculating vorticity gradients in non-symmetric memory. --- src/parameterizations/lateral/MOM_hor_visc.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 7ff6ac4018..3b0014f92a 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -510,11 +510,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, endif ! Vorticity gradient - do J=Jsq-1,Jeq+1 ; do I=Isq,Ieq+1 + do J=js-2,Jeq+1 ; do I=is-1,Ieq+1 vort_xy_dx(i,J) = CS%DY_dxBu(I,J)*(vort_xy(I,J)*G%IdyCu(I,j) - vort_xy(I-1,J)*G%IdyCu(I-1,j)) enddo ; enddo - do J=Jsq,Jeq+1 ; do I=Isq-1,Ieq+1 + do J=js-1,Jeq+1 ; do I=is-2,Ieq+1 vort_xy_dy(I,j) = CS%DX_dyBu(I,J)*(vort_xy(I,J)*G%IdxCv(i,J) - vort_xy(I,J-1)*G%IdxCv(i,J-1)) enddo ; enddo