From fd4faf566d336c77dfd800959e3885e24358ea65 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 4 Dec 2017 15:03:37 -0500 Subject: [PATCH 1/9] Converted local scalar Vort_mag into local 2d array - In preparation for modularizing the Leith calculations so that we can add QG-Leith (Bachman et al., 2017), we have made the local scalar that held the vorticity magnitude a 2d array. --- .../lateral/MOM_hor_visc.F90 | 25 +++++++++++-------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index bba797523b..46d9ce52fc 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -281,14 +281,16 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, str_xx,& ! str_xx is the diagonal term in the stress tensor (H m2 s-2) bhstr_xx,& ! A copy of str_xx that only contains the biharmonic contribution (H m2 s-2) div_xx, & ! horizontal divergence (du/dx + dv/dy) (1/sec) including metric terms - FrictWorkIntz ! depth integrated energy dissipated by lateral friction (W/m2) + FrictWorkIntz, & ! depth integrated energy dissipated by lateral friction (W/m2) + vert_vort_mag_h ! Magnitude of vertical vorticity at h-points |dv/dx - du/dy| (1/sec) real, dimension(SZIB_(G),SZJB_(G)) :: & 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) - vort_xy ! vertical vorticity (dv/dx - du/dy) (1/sec) including metric terms + vort_xy, & ! Vertical vorticity (dv/dx - du/dy) (1/sec) + vert_vort_mag_q ! Magnitude of vertical vorticity at q-points |dv/dx - du/dy| (1/sec) 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 @@ -317,7 +319,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, ! viscosity. Here set equal to nondimensional Laplacian Leith constant. ! This is set equal to zero if modified Leith is not used. 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 @@ -379,11 +381,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,CS,G,GV,u,v,is,js,ie,je,h, & !$OMP rescale_Kh,VarMix,h_neglect,h_neglect3, & !$OMP Kh_h,Ah_h,Kh_q,Ah_q,diffu,apply_OBC,OBC,diffv, & -!$OMP find_FrictWork,FrictWork,use_MEKE_Ku,MEKE) & +!$OMP find_FrictWork,FrictWork,use_MEKE_Ku,MEKE, & +!$OMP vert_vort_mag_h, vert_vort_mag_q) & !$OMP private(u0, v0, sh_xx, str_xx, visc_bound_rem, & !$OMP sh_xy, str_xy, Ah, Kh, AhSm, KhSm, dvdx, dudy, & !$OMP bhstr_xx, bhstr_xy,FatH,RoScl, hu, hv, h_u, h_v, & -!$OMP vort_xy,vort_xy_dx,vort_xy_dy,Vort_mag,AhLth,KhLth, & +!$OMP vort_xy,vort_xy_dx,vort_xy_dy,AhLth,KhLth, & !$OMP div_xx, div_xx_dx, div_xx_dy, mod_Leith, & !$OMP Shear_mag, h2uq, h2vq, hq, Kh_scale, hrat_min) do k=1,nz @@ -600,7 +603,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, (sh_xy(I-1,J)*sh_xy(I-1,J) + sh_xy(I,J-1)*sh_xy(I,J-1)))) endif if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then - Vort_mag = sqrt( & + vert_vort_mag_h(i,j) = 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))) + & mod_Leith*0.5*((div_xx_dx(I,j)*div_xx_dx(I,j) + div_xx_dx(I-1,j)*div_xx_dx(I-1,j)) + & @@ -621,7 +624,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, 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 + KhLth = CS%LAPLAC3_CONST_xx(i,j) * vert_vort_mag_h(i,j) 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)) @@ -663,7 +666,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, endif endif if (CS%Leith_Ah) & - AhLth = Vort_mag * (CS%BIHARM_CONST_xx(i,j)) + AhLth = vert_vort_mag_h(i,j) * (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)) @@ -730,7 +733,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, (sh_xx(i,j+1)*sh_xx(i,j+1) + sh_xx(i+1,j)*sh_xx(i+1,j)))) endif if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) & - Vort_mag = sqrt( & + vert_vort_mag_q(I,J) = 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))) + & mod_Leith*0.5*((div_xx_dx(I,j)*div_xx_dx(I,j) + div_xx_dx(I,j+1)*div_xx_dx(I,j+1)) + & @@ -778,7 +781,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, 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 + KhLth = CS%LAPLAC3_CONST_xy(I,J) * vert_vort_mag_q(I,J) 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)) @@ -823,7 +826,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, endif endif if (CS%Leith_Ah) & - AhLth = Vort_mag * (CS%BIHARM5_CONST_xy(I,J)) + AhLth = vert_vort_mag_q(I,J) * (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)) From e5d4bdfde0da7fe5c2a5ceccfe7510209e2cc178 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 4 Dec 2017 15:36:43 -0500 Subject: [PATCH 2/9] Consolidated Leith related calculations - Grouped all the Leith related calculations in preparations for turning them into a single subroutine call. - Also duplicated the calculation of dvdx and dudy expecting that we might soon change the discretization for vorticity. --- .../lateral/MOM_hor_visc.F90 | 78 +++++++++++-------- 1 file changed, 47 insertions(+), 31 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 46d9ce52fc..451c366409 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -535,25 +535,6 @@ 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-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=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 - ! Divergence gradient do j=Jsq-1,Jeq+2 ; do I=is-2,Ieq+1 div_xx_dx(I,j) = G%IdxCu(I,j)*(div_xx(i+1,j) - div_xx(i,j)) @@ -596,18 +577,59 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, endif; endif endif - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - 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)))) + if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then + ! Components for the vertical vorticity + ! Note this a simple re-calculation of shearing components using the same discretization. + ! We will consider using a circulation based calculation of vorticity later. + ! Also note this will need OBC boundary conditions re-applied... + do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 + dvdx(I,J) = CS%DY_dxBu(I,J)*(v(i+1,J,k)*G%IdyCv(i+1,J) - v(i,J,k)*G%IdyCv(i,J)) + 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 + ! Vorticity + 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 - if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then + + ! Vorticity gradient + 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=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 + + ! Magnitude of vorticity at h-points + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 vert_vort_mag_h(i,j) = 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))) + & mod_Leith*0.5*((div_xx_dx(I,j)*div_xx_dx(I,j) + div_xx_dx(I-1,j)*div_xx_dx(I-1,j)) + & (div_xx_dy(i,J)*div_xx_dy(i,J) + div_xx_dy(i,J-1)*div_xx_dy(i,J-1)))) + enddo ; enddo + + ! Magnitude of vorticity at q-points + do J=js-1,Jeq ; do I=is-1,Ieq + vert_vort_mag_q(I,J) = 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))) + & + mod_Leith*0.5*((div_xx_dx(I,j)*div_xx_dx(I,j) + div_xx_dx(I,j+1)*div_xx_dx(I,j+1)) + & + (div_xx_dy(i,J)*div_xx_dy(i,J) + div_xx_dy(i+1,J)*div_xx_dy(i+1,J)))) + enddo ; enddo + endif + + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + 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)))) endif 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)) / & @@ -732,12 +754,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, 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)))) endif - if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) & - vert_vort_mag_q(I,J) = 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))) + & - mod_Leith*0.5*((div_xx_dx(I,j)*div_xx_dx(I,j) + div_xx_dx(I,j+1)*div_xx_dx(I,j+1)) + & - (div_xx_dy(i,J)*div_xx_dy(i,J) + div_xx_dy(i+1,J)*div_xx_dy(i+1,J)))) h2uq = 4.0 * h_u(I,j) * h_u(I,j+1) h2vq = 4.0 * h_v(i,J) * h_v(i+1,J) From 273745c0b00d2c71d266d5dcf08de16aeda968c5 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 4 Dec 2017 16:17:12 -0500 Subject: [PATCH 3/9] Moved Leith vorticity magnitude calc into MOM_lateral_mixing_coeffs.F90 - Block of calculations for vorticity magnitude used in Leith are now in a subroutine in MOM_lateral_mixing_coeffs.F90 --- .../lateral/MOM_hor_visc.F90 | 75 +------------ .../lateral/MOM_lateral_mixing_coeffs.F90 | 106 +++++++++++++++++- 2 files changed, 109 insertions(+), 72 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 451c366409..a12721a2ea 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -76,6 +76,7 @@ module MOM_hor_visc use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type use MOM_lateral_mixing_coeffs, only : VarMix_CS +use MOM_lateral_mixing_coeffs, only : calc_vert_vort_mag use MOM_MEKE_types, only : MEKE_type use MOM_open_boundary, only : ocean_OBC_type, OBC_DIRECTION_E, OBC_DIRECTION_W use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S, OBC_NONE @@ -280,7 +281,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, sh_xx, & ! horizontal tension (du/dx - dv/dy) (1/sec) including metric terms str_xx,& ! str_xx is the diagonal term in the stress tensor (H m2 s-2) bhstr_xx,& ! A copy of str_xx that only contains the biharmonic contribution (H m2 s-2) - div_xx, & ! horizontal divergence (du/dx + dv/dy) (1/sec) including metric terms FrictWorkIntz, & ! depth integrated energy dissipated by lateral friction (W/m2) vert_vort_mag_h ! Magnitude of vertical vorticity at h-points |dv/dx - du/dy| (1/sec) @@ -289,17 +289,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, 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) - vort_xy, & ! Vertical vorticity (dv/dx - du/dy) (1/sec) vert_vort_mag_q ! Magnitude of vertical vorticity at q-points |dv/dx - du/dy| (1/sec) - 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 - div_xx_dy ! y-derivative of horizontal divergence (d/dy(du/dx + dv/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 - div_xx_dx ! x-derivative of horizontal divergence (d/dx(du/dx + dv/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) @@ -386,8 +377,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, !$OMP private(u0, v0, sh_xx, str_xx, visc_bound_rem, & !$OMP sh_xy, str_xy, Ah, Kh, AhSm, KhSm, dvdx, dudy, & !$OMP bhstr_xx, bhstr_xy,FatH,RoScl, hu, hv, h_u, h_v, & -!$OMP vort_xy,vort_xy_dx,vort_xy_dy,AhLth,KhLth, & -!$OMP div_xx, div_xx_dx, div_xx_dy, mod_Leith, & +!$OMP AhLth,KhLth, & +!$OMP mod_Leith, & !$OMP Shear_mag, h2uq, h2vq, hq, Kh_scale, hrat_min) do k=1,nz @@ -409,11 +400,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, 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))) - div_xx(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%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) + h_neglect) enddo ; enddo ! Components for the shearing strain do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 @@ -535,15 +521,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, enddo ; enddo endif -! Divergence gradient - do j=Jsq-1,Jeq+2 ; do I=is-2,Ieq+1 - div_xx_dx(I,j) = G%IdxCu(I,j)*(div_xx(i+1,j) - div_xx(i,j)) - enddo ; enddo - - do J=js-2,Jeq+1 ; do i=Isq-1,Ieq+2 - div_xx_dy(i,J) = G%IdyCv(i,J)*(div_xx(i,j+1) - div_xx(i,j)) - enddo ; enddo - ! Coefficient for modified Leith if (CS%Modified_Leith) then mod_Leith = 1.0 @@ -578,51 +555,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, endif if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then - ! Components for the vertical vorticity - ! Note this a simple re-calculation of shearing components using the same discretization. - ! We will consider using a circulation based calculation of vorticity later. - ! Also note this will need OBC boundary conditions re-applied... - do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 - dvdx(I,J) = CS%DY_dxBu(I,J)*(v(i+1,J,k)*G%IdyCv(i+1,J) - v(i,J,k)*G%IdyCv(i,J)) - 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 - ! Vorticity - 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-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=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 - - ! Magnitude of vorticity at h-points - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - vert_vort_mag_h(i,j) = 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))) + & - mod_Leith*0.5*((div_xx_dx(I,j)*div_xx_dx(I,j) + div_xx_dx(I-1,j)*div_xx_dx(I-1,j)) + & - (div_xx_dy(i,J)*div_xx_dy(i,J) + div_xx_dy(i,J-1)*div_xx_dy(i,J-1)))) - enddo ; enddo - - ! Magnitude of vorticity at q-points - do J=js-1,Jeq ; do I=is-1,Ieq - vert_vort_mag_q(I,J) = 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))) + & - mod_Leith*0.5*((div_xx_dx(I,j)*div_xx_dx(I,j) + div_xx_dx(I,j+1)*div_xx_dx(I,j+1)) + & - (div_xx_dy(i,J)*div_xx_dy(i,J) + div_xx_dy(i+1,J)*div_xx_dy(i+1,J)))) - enddo ; enddo + call calc_vert_vort_mag(G, GV, u, v, h, k, CS%no_slip, mod_Leith, vert_vort_mag_h, vert_vort_mag_q) endif do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 7f981f7c04..1e5e7fdbea 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -125,6 +125,7 @@ module MOM_lateral_mixing_coeffs end type VarMix_CS public VarMix_init, calc_slope_functions, calc_resoln_function +public calc_vert_vort_mag contains @@ -719,7 +720,110 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, CS, e, calculate_slopes) end subroutine calc_slope_functions_using_just_e -!> Initializes the variables mixing coefficients container +!> Calculates the magnitude of the vertical component of vorticity for use in the Leith-like schemes +subroutine calc_vert_vort_mag(G, GV, u, v, h, k, no_slip, mod_Leith, vert_vort_mag_h, vert_vort_mag_q) + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal flow (m s-1) + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional flow (m s-1) + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: h !< Layer thickness (m or kg m-2) + integer, intent(in) :: k !< Layer for which to calculate vorticity magnitude + logical, intent(in) :: no_slip !< True if vorticity should have no-slip BCs + real, intent(in) :: mod_Leith !< Non-dimensional coefficient multiplying the + !! divergence contribution to the Leith viscosity. + !! Set to zero for conventional Leith. + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: vert_vort_mag_h !< Magnitude of vertical component + !! of vorticity at h-ponts (s-1) + real, dimension(SZIB_(G),SZJB_(G)), intent(out) :: vert_vort_mag_q !< Magnitude of vertical component + !! of vorticity at q-ponts (s-1) + ! Local variables + real, dimension(SZIB_(G),SZJB_(G)) :: vort_xy, & ! Vertical vorticity (dv/dx - du/dy) (s-1) + dudy, & ! Meridional shear of zonal velocity (s-1) + dvdx ! Zonal shear of meridional velocity (s-1) + real, dimension(SZI_(G),SZJB_(G)) :: & + vort_xy_dx, & ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) (m-1 s-1) + div_xx_dy ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) (m-1 s-1) + + real, dimension(SZIB_(G),SZJ_(G)) :: & + vort_xy_dy, & ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) (m-1 s-1) + div_xx_dx ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) (m-1 s-1) + real, dimension(SZI_(G),SZJ_(G)) :: div_xx ! Estimate of horizontal divergence at h-points (s-1) + real :: DY_dxBu, DX_dyBu + integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + + ! Divergence + do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 + div_xx(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%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) + enddo ; enddo + + ! Divergence gradient + do j=Jsq-1,Jeq+2 ; do I=is-2,Ieq+1 + div_xx_dx(I,j) = G%IdxCu(I,j)*(div_xx(i+1,j) - div_xx(i,j)) + enddo ; enddo + + do J=js-2,Jeq+1 ; do i=Isq-1,Ieq+2 + div_xx_dy(i,J) = G%IdyCv(i,J)*(div_xx(i,j+1) - div_xx(i,j)) + enddo ; enddo + + ! Components for the vertical vorticity + ! Note this a simple re-calculation of shearing components using the same discretization. + ! We will consider using a circulation based calculation of vorticity later. + ! Also note this will need OBC boundary conditions re-applied... + do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 + DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J) + dvdx(I,J) = DY_dxBu * (v(i+1,J,k) * G%IdyCv(i+1,J) - v(i,J,k) * G%IdyCv(i,J)) + DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J) + dudy(I,J) = DX_dyBu * (u(I,j+1,k) * G%IdxCu(I,j+1) - u(I,j,k) * G%IdxCu(I,j)) + enddo ; enddo + + ! Vorticity + if (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-2,Jeq+1 ; do I=is-1,Ieq+1 + DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J) + vort_xy_dx(i,J) = DY_dxBu * (vort_xy(I,J) * G%IdyCu(I,j) - vort_xy(I-1,J) * G%IdyCu(I-1,j)) + enddo ; enddo + + do J=js-1,Jeq+1 ; do I=is-2,Ieq+1 + DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J) + vort_xy_dy(I,j) = DX_dyBu * (vort_xy(I,J) * G%IdxCv(i,J) - vort_xy(I,J-1) * G%IdxCv(i,J-1)) + enddo ; enddo + + ! Magnitude of vorticity at h-points + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + vert_vort_mag_h(i,j) = 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))) + & + mod_Leith*0.5*((div_xx_dx(I,j)*div_xx_dx(I,j) + div_xx_dx(I-1,j)*div_xx_dx(I-1,j)) + & + (div_xx_dy(i,J)*div_xx_dy(i,J) + div_xx_dy(i,J-1)*div_xx_dy(i,J-1)))) + enddo ; enddo + + ! Magnitude of vorticity at q-points + do J=js-1,Jeq ; do I=is-1,Ieq + vert_vort_mag_q(I,J) = 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))) + & + mod_Leith*0.5*((div_xx_dx(I,j)*div_xx_dx(I,j) + div_xx_dx(I,j+1)*div_xx_dx(I,j+1)) + & + (div_xx_dy(i,J)*div_xx_dy(i,J) + div_xx_dy(i+1,J)*div_xx_dy(i+1,J)))) + enddo ; enddo + +end subroutine calc_vert_vort_mag + subroutine VarMix_init(Time, G, param_file, diag, CS) type(time_type), intent(in) :: Time !< Current model time type(ocean_grid_type), intent(in) :: G !< Ocean grid structure From 7dc2510d124e703052e4961380b1f6501859e60f Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 4 Dec 2017 16:43:16 -0500 Subject: [PATCH 4/9] Added beta to Leith vorticity magnitude - Added beta from grid type to vorticity magnitude --- src/parameterizations/lateral/MOM_hor_visc.F90 | 2 +- .../lateral/MOM_lateral_mixing_coeffs.F90 | 13 ++++++++++++- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index a12721a2ea..1959a200f8 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -555,7 +555,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, endif if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then - call calc_vert_vort_mag(G, GV, u, v, h, k, CS%no_slip, mod_Leith, vert_vort_mag_h, vert_vort_mag_q) + call calc_vert_vort_mag(G, GV, u, v, h, k, CS%no_slip, mod_Leith, .false., vert_vort_mag_h, vert_vort_mag_q) endif do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 1e5e7fdbea..bc4b06fa8b 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -721,7 +721,7 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, CS, e, calculate_slopes) end subroutine calc_slope_functions_using_just_e !> Calculates the magnitude of the vertical component of vorticity for use in the Leith-like schemes -subroutine calc_vert_vort_mag(G, GV, u, v, h, k, no_slip, mod_Leith, vert_vort_mag_h, vert_vort_mag_q) +subroutine calc_vert_vort_mag(G, GV, u, v, h, k, no_slip, mod_Leith, QG_Leith, vert_vort_mag_h, vert_vort_mag_q) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal flow (m s-1) @@ -732,6 +732,7 @@ subroutine calc_vert_vort_mag(G, GV, u, v, h, k, no_slip, mod_Leith, vert_vort_m real, intent(in) :: mod_Leith !< Non-dimensional coefficient multiplying the !! divergence contribution to the Leith viscosity. !! Set to zero for conventional Leith. + logical, intent(in) :: QG_Leith !< True if using QG Leith scheme real, dimension(SZI_(G),SZJ_(G)), intent(out) :: vert_vort_mag_h !< Magnitude of vertical component !! of vorticity at h-ponts (s-1) real, dimension(SZIB_(G),SZJB_(G)), intent(out) :: vert_vort_mag_q !< Magnitude of vertical component @@ -804,6 +805,16 @@ subroutine calc_vert_vort_mag(G, GV, u, v, h, k, no_slip, mod_Leith, vert_vort_m vort_xy_dy(I,j) = DX_dyBu * (vort_xy(I,J) * G%IdxCv(i,J) - vort_xy(I,J-1) * G%IdxCv(i,J-1)) enddo ; enddo + ! Add in beta for QG Leith + if (QG_Leith) then + do J=js-2,Jeq+1 ; do I=is-1,Ieq+1 + vort_xy_dx(i,J) = vort_xy_dx(i,J) + 0.5 * ( G%dF_dx(i,j) + G%dF_dx(i,j+1) ) + enddo ; enddo + do J=js-1,Jeq+1 ; do I=is-2,Ieq+1 + vort_xy_dy(I,j) = vort_xy_dy(I,j) + 0.5 * ( G%dF_dy(i,j) + G%dF_dy(i+1,j) ) + enddo ; enddo + endif + ! Magnitude of vorticity at h-points do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 vert_vort_mag_h(i,j) = sqrt( & From c59f87e1cb2a0e7fa203785cd9e521d387f21405 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 5 Dec 2017 11:23:21 -0500 Subject: [PATCH 5/9] Read Leith parameters in MOM_lateral_coeffs.F90 - Now reading parameters for Leith in MOM_lateral_coeffs.F90 (in additional to MOM_hor_visc.F90 which we'll remove later). - Removed parameters from dummy argument list. --- .../lateral/MOM_hor_visc.F90 | 12 +--- .../lateral/MOM_lateral_mixing_coeffs.F90 | 64 ++++++++++++++++--- 2 files changed, 56 insertions(+), 20 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 1959a200f8..f9da5ab15e 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -306,9 +306,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, real :: KhSm ! Smagorinsky Laplacian viscosity (m2/s) real :: AhLth ! 2D Leith biharmonic viscosity (m4/s) real :: KhLth ! 2D Leith Laplacian viscosity (m2/s) - real :: mod_Leith ! nondimensional coefficient for divergence part of modified Leith - ! viscosity. Here set equal to nondimensional Laplacian Leith constant. - ! This is set equal to zero if modified Leith is not used. 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). @@ -521,13 +518,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, enddo ; enddo endif -! Coefficient for modified Leith - if (CS%Modified_Leith) then - mod_Leith = 1.0 - else - mod_Leith = 0.0 - endif - ! 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 @@ -555,7 +545,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, endif if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then - call calc_vert_vort_mag(G, GV, u, v, h, k, CS%no_slip, mod_Leith, .false., vert_vort_mag_h, vert_vort_mag_q) + call calc_vert_vort_mag(VarMix, G, GV, u, v, h, k, vert_vort_mag_h, vert_vort_mag_q) endif do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index bc4b06fa8b..e5275c78ce 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -109,6 +109,21 @@ module MOM_lateral_mixing_coeffs !! and especially 2 are coded to be more efficient. real :: Visbeck_S_max !< Upper bound on slope used in Eady growth rate (nondim). + ! Leith parameters + logical :: use_QG_Leith !< If true, enables the QG Leith scheme + logical :: Leith_Kh !< If true, enables the Leith scheme + logical :: modified_Leith !< if true, include the divergence contribution to Leith viscosity + real :: Leith_Lap_const !< The non-dimensional coefficient in the Leith viscosity + logical :: Leith_Ah !< If true, enables the bi-harmonic Leith scheme + real :: Leith_bi_const !< The non-dimensional coefficient in the bi-harmonic Leith viscosity + logical :: no_slip !< If true, no slip boundary conditions are used. + !! Otherwise free slip boundary conditions are assumed. + !! The implementation of the free slip boundary + !! conditions on a C-grid is much cleaner than the + !! no slip boundary conditions. The use of free slip + !! b.c.s is strongly encouraged. The no slip b.c.s + !! are not implemented with the biharmonic viscosity. + ! Diagnostics !>@{ !! Diagnostic identifier @@ -721,18 +736,14 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, CS, e, calculate_slopes) end subroutine calc_slope_functions_using_just_e !> Calculates the magnitude of the vertical component of vorticity for use in the Leith-like schemes -subroutine calc_vert_vort_mag(G, GV, u, v, h, k, no_slip, mod_Leith, QG_Leith, vert_vort_mag_h, vert_vort_mag_q) +subroutine calc_vert_vort_mag(CS, G, GV, u, v, h, k, vert_vort_mag_h, vert_vort_mag_q) + type(VarMix_CS), pointer :: CS !< Variable mixing coefficients type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal flow (m s-1) real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional flow (m s-1) real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: h !< Layer thickness (m or kg m-2) integer, intent(in) :: k !< Layer for which to calculate vorticity magnitude - logical, intent(in) :: no_slip !< True if vorticity should have no-slip BCs - real, intent(in) :: mod_Leith !< Non-dimensional coefficient multiplying the - !! divergence contribution to the Leith viscosity. - !! Set to zero for conventional Leith. - logical, intent(in) :: QG_Leith !< True if using QG Leith scheme real, dimension(SZI_(G),SZJ_(G)), intent(out) :: vert_vort_mag_h !< Magnitude of vertical component !! of vorticity at h-ponts (s-1) real, dimension(SZIB_(G),SZJB_(G)), intent(out) :: vert_vort_mag_q !< Magnitude of vertical component @@ -749,7 +760,7 @@ subroutine calc_vert_vort_mag(G, GV, u, v, h, k, no_slip, mod_Leith, QG_Leith, v vort_xy_dy, & ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) (m-1 s-1) div_xx_dx ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) (m-1 s-1) real, dimension(SZI_(G),SZJ_(G)) :: div_xx ! Estimate of horizontal divergence at h-points (s-1) - real :: DY_dxBu, DX_dyBu + real :: mod_Leith, DY_dxBu, DX_dyBu integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -784,7 +795,7 @@ subroutine calc_vert_vort_mag(G, GV, u, v, h, k, no_slip, mod_Leith, QG_Leith, v enddo ; enddo ! Vorticity - if (no_slip) then + 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 @@ -806,7 +817,7 @@ subroutine calc_vert_vort_mag(G, GV, u, v, h, k, no_slip, mod_Leith, QG_Leith, v enddo ; enddo ! Add in beta for QG Leith - if (QG_Leith) then + if (CS%use_QG_Leith) then do J=js-2,Jeq+1 ; do I=is-1,Ieq+1 vort_xy_dx(i,J) = vort_xy_dx(i,J) + 0.5 * ( G%dF_dx(i,j) + G%dF_dx(i,j+1) ) enddo ; enddo @@ -815,6 +826,8 @@ subroutine calc_vert_vort_mag(G, GV, u, v, h, k, no_slip, mod_Leith, QG_Leith, v enddo ; enddo endif + mod_Leith = 0.; if (CS%modified_Leith) mod_Leith = 1.0 + ! Magnitude of vorticity at h-points do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 vert_vort_mag_h(i,j) = sqrt( & @@ -1131,6 +1144,39 @@ subroutine VarMix_init(Time, G, param_file, diag, CS) call wave_speed_init(CS%wave_speed_CSp, use_ebt_mode=CS%Resoln_use_ebt, mono_N2_depth=N2_filter_depth) endif + ! Leith parameters + call get_param(param_file, mdl, "LEITH_KH", CS%Leith_Kh, & + "If true, use a Leith nonlinear eddy viscosity.", & + default=CS%use_QG_Leith) + call get_param(param_file, mdl, "LEITH_AH", CS%Leith_Ah, & + "If true, use a biharmonic Leith nonlinear eddy \n"//& + "viscosity.", default=.false.) + call get_param(param_file, mdl, "USE_QG_LEITH", CS%use_QG_Leith, & + "If true, use the QG Leith nonlinear eddy viscosity.", & + default=.false.) + call get_param(param_file, mdl, "MODIFIED_LEITH", CS%modified_Leith, & + "If true, add a term to Leith viscosity which is \n"//& + "proportional to the gradient of divergence.", & + default=.false.) + call get_param(param_file, mdl, "LEITH_LAP_CONST", CS%Leith_Lap_const, & + "The nondimensional Laplacian Leith constant, \n"//& + "often set to 1.0", units="nondim", default=0.0, & + fail_if_missing = CS%Leith_Kh) + call get_param(param_file, mdl, "LEITH_BI_CONST", CS%Leith_bi_const, & + "The nondimensional biharmonic Leith constant, \n"//& + "typical values are thus far undetermined.", units="nondim", default=0.0, & + fail_if_missing = CS%Leith_Ah) + if (CS%Leith_Kh .or. CS%Leith_Ah) then + in_use = .true. + call get_param(param_file, mdl, "NOSLIP", CS%no_slip, & + "If true, no slip boundary conditions are used; otherwise \n"//& + "free slip boundary conditions are assumed. The \n"//& + "implementation of the free slip BCs on a C-grid is much \n"//& + "cleaner than the no slip BCs. The use of free slip BCs \n"//& + "is strongly encouraged, and no slip BCs are not used with \n"//& + "the biharmonic viscosity.", default=.false.) + endif + ! If nothing is being stored in this class then deallocate if (in_use) then CS%use_variable_mixing = .true. From c546d3ddc3276a0a1cda73d4852374d843f441e4 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 5 Dec 2017 11:25:13 -0500 Subject: [PATCH 6/9] Replaced calc_vert_vort_mag() with calc_Leith_viscosity() - The subroutine that was calculating terms used in the Leith viscosities now returns the viscosities themselves. --- .../lateral/MOM_hor_visc.F90 | 27 +++---- .../lateral/MOM_lateral_mixing_coeffs.F90 | 81 +++++++++++++++---- 2 files changed, 76 insertions(+), 32 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index f9da5ab15e..cc3359a070 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -76,7 +76,7 @@ module MOM_hor_visc use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type use MOM_lateral_mixing_coeffs, only : VarMix_CS -use MOM_lateral_mixing_coeffs, only : calc_vert_vort_mag +use MOM_lateral_mixing_coeffs, only : calc_Leith_viscosity use MOM_MEKE_types, only : MEKE_type use MOM_open_boundary, only : ocean_OBC_type, OBC_DIRECTION_E, OBC_DIRECTION_W use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S, OBC_NONE @@ -282,14 +282,16 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, str_xx,& ! str_xx is the diagonal term in the stress tensor (H m2 s-2) bhstr_xx,& ! A copy of str_xx that only contains the biharmonic contribution (H m2 s-2) FrictWorkIntz, & ! depth integrated energy dissipated by lateral friction (W/m2) - vert_vort_mag_h ! Magnitude of vertical vorticity at h-points |dv/dx - du/dy| (1/sec) + Leith_Kh_h, & ! Leith Laplacian viscosity at h-points (m2 s-1) + Leith_Ah_h ! Leith bi-harmonic viscosity at h-points (m4 s-1) real, dimension(SZIB_(G),SZJB_(G)) :: & 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) - vert_vort_mag_q ! Magnitude of vertical vorticity at q-points |dv/dx - du/dy| (1/sec) + Leith_Kh_q, & ! Leith Laplacian viscosity at q-points (m2 s-1) + Leith_Ah_q ! Leith bi-harmonic viscosity at q-points (m4 s-1) real, dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: & Ah_q, & ! biharmonic viscosity at corner points (m4/s) @@ -307,7 +309,6 @@ 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 :: 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 @@ -370,12 +371,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, !$OMP rescale_Kh,VarMix,h_neglect,h_neglect3, & !$OMP Kh_h,Ah_h,Kh_q,Ah_q,diffu,apply_OBC,OBC,diffv, & !$OMP find_FrictWork,FrictWork,use_MEKE_Ku,MEKE, & -!$OMP vert_vort_mag_h, vert_vort_mag_q) & +!$OMP Leith_Kh_h, Leith_Kh_q, Leith_Ah_h, Leith_Ah_q) & !$OMP private(u0, v0, sh_xx, str_xx, visc_bound_rem, & !$OMP sh_xy, str_xy, Ah, Kh, AhSm, KhSm, dvdx, dudy, & !$OMP bhstr_xx, bhstr_xy,FatH,RoScl, hu, hv, h_u, h_v, & !$OMP AhLth,KhLth, & -!$OMP mod_Leith, & !$OMP Shear_mag, h2uq, h2vq, hq, Kh_scale, hrat_min) do k=1,nz @@ -545,7 +545,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, endif if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then - call calc_vert_vort_mag(VarMix, G, GV, u, v, h, k, vert_vort_mag_h, vert_vort_mag_q) + call calc_Leith_viscosity(VarMix, G, GV, u, v, h, k, Leith_Kh_h, Leith_Kh_q, Leith_Ah_h, Leith_Ah_q) endif do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -568,8 +568,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, 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) * vert_vort_mag_h(i,j) + if (CS%Leith_Kh) KhLth = Leith_Kh_h(i,j) + ! Note: move Leith outside of resolution function 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)) @@ -610,8 +610,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, AhSm = CS%BIHARM_CONST_xx(i,j) * Shear_mag endif endif - if (CS%Leith_Ah) & - AhLth = vert_vort_mag_h(i,j) * (CS%BIHARM_CONST_xx(i,j)) + if (CS%Leith_Ah) AhLth = Leith_Ah_h(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)) @@ -719,8 +718,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, 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) * vert_vort_mag_q(I,J) + if (CS%Leith_Kh) KhLth = Leith_Kh_q(I,J) 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)) @@ -764,8 +762,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, AhSm = CS%BIHARM_CONST_xy(I,J) * Shear_mag endif endif - if (CS%Leith_Ah) & - AhLth = vert_vort_mag_q(I,J) * (CS%BIHARM5_CONST_xy(I,J)) + if (CS%Leith_Ah) AhLth = Leith_Ah_q(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)) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index e5275c78ce..764e393188 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -113,9 +113,7 @@ module MOM_lateral_mixing_coeffs logical :: use_QG_Leith !< If true, enables the QG Leith scheme logical :: Leith_Kh !< If true, enables the Leith scheme logical :: modified_Leith !< if true, include the divergence contribution to Leith viscosity - real :: Leith_Lap_const !< The non-dimensional coefficient in the Leith viscosity logical :: Leith_Ah !< If true, enables the bi-harmonic Leith scheme - real :: Leith_bi_const !< The non-dimensional coefficient in the bi-harmonic Leith viscosity logical :: no_slip !< If true, no slip boundary conditions are used. !! Otherwise free slip boundary conditions are assumed. !! The implementation of the free slip boundary @@ -123,6 +121,12 @@ module MOM_lateral_mixing_coeffs !! no slip boundary conditions. The use of free slip !! b.c.s is strongly encouraged. The no slip b.c.s !! are not implemented with the biharmonic viscosity. + real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & + Laplac3_const_xx, & ! Laplacian metric-dependent constants (nondim) + biharm5_const_xx ! Biharmonic metric-dependent constants (nondim) + real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: & + Laplac3_const_xy, & ! Laplacian metric-dependent constants (nondim) + biharm5_const_xy ! Biharmonic metric-dependent constants (nondim) ! Diagnostics !>@{ @@ -140,7 +144,7 @@ module MOM_lateral_mixing_coeffs end type VarMix_CS public VarMix_init, calc_slope_functions, calc_resoln_function -public calc_vert_vort_mag +public calc_Leith_viscosity contains @@ -735,8 +739,8 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, CS, e, calculate_slopes) end subroutine calc_slope_functions_using_just_e -!> Calculates the magnitude of the vertical component of vorticity for use in the Leith-like schemes -subroutine calc_vert_vort_mag(CS, G, GV, u, v, h, k, vert_vort_mag_h, vert_vort_mag_q) +!> Calculates the Leith Laplacian and bi-harmonic viscosity coefficients +subroutine calc_Leith_viscosity(CS, G, GV, u, v, h, k, Leith_Kh_h, Leith_Kh_q, Leith_Ah_h, Leith_Ah_q) type(VarMix_CS), pointer :: CS !< Variable mixing coefficients type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. @@ -744,10 +748,10 @@ subroutine calc_vert_vort_mag(CS, G, GV, u, v, h, k, vert_vort_mag_h, vert_vort_ real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional flow (m s-1) real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: h !< Layer thickness (m or kg m-2) integer, intent(in) :: k !< Layer for which to calculate vorticity magnitude - real, dimension(SZI_(G),SZJ_(G)), intent(out) :: vert_vort_mag_h !< Magnitude of vertical component - !! of vorticity at h-ponts (s-1) - real, dimension(SZIB_(G),SZJB_(G)), intent(out) :: vert_vort_mag_q !< Magnitude of vertical component - !! of vorticity at q-ponts (s-1) + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: Leith_Kh_h !< Leith Laplacian viscosity at h-points (m2 s-1) + real, dimension(SZIB_(G),SZJB_(G)), intent(out) :: Leith_Kh_q !< Leith Laplacian viscosity at q-points (m2 s-1) + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: Leith_Ah_h !< Leith bi-harmonic viscosity at h-points (m4 s-1) + real, dimension(SZIB_(G),SZJB_(G)), intent(out) :: Leith_Ah_q !< Leith bi-harmonic viscosity at q-points (m4 s-1) ! Local variables real, dimension(SZIB_(G),SZJB_(G)) :: vort_xy, & ! Vertical vorticity (dv/dx - du/dy) (s-1) dudy, & ! Meridional shear of zonal velocity (s-1) @@ -760,7 +764,7 @@ subroutine calc_vert_vort_mag(CS, G, GV, u, v, h, k, vert_vort_mag_h, vert_vort_ vort_xy_dy, & ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) (m-1 s-1) div_xx_dx ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) (m-1 s-1) real, dimension(SZI_(G),SZJ_(G)) :: div_xx ! Estimate of horizontal divergence at h-points (s-1) - real :: mod_Leith, DY_dxBu, DX_dyBu + real :: mod_Leith, DY_dxBu, DX_dyBu, vert_vort_mag integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -828,25 +832,29 @@ subroutine calc_vert_vort_mag(CS, G, GV, u, v, h, k, vert_vort_mag_h, vert_vort_ mod_Leith = 0.; if (CS%modified_Leith) mod_Leith = 1.0 - ! Magnitude of vorticity at h-points + ! h-point viscosities do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - vert_vort_mag_h(i,j) = sqrt( & + vert_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))) + & mod_Leith*0.5*((div_xx_dx(I,j)*div_xx_dx(I,j) + div_xx_dx(I-1,j)*div_xx_dx(I-1,j)) + & (div_xx_dy(i,J)*div_xx_dy(i,J) + div_xx_dy(i,J-1)*div_xx_dy(i,J-1)))) + if (CS%Leith_Kh) Leith_Kh_h(i,j) = CS%Laplac3_const_xx(i,j) * vert_vort_mag + if (CS%Leith_Ah) Leith_Ah_h(i,j) = CS%biharm5_const_xx(i,j) * vert_vort_mag enddo ; enddo - ! Magnitude of vorticity at q-points + ! q-point viscosities do J=js-1,Jeq ; do I=is-1,Ieq - vert_vort_mag_q(I,J) = sqrt( & + vert_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))) + & mod_Leith*0.5*((div_xx_dx(I,j)*div_xx_dx(I,j) + div_xx_dx(I,j+1)*div_xx_dx(I,j+1)) + & (div_xx_dy(i,J)*div_xx_dy(i,J) + div_xx_dy(i+1,J)*div_xx_dy(i+1,J)))) + if (CS%Leith_Kh) Leith_Kh_q(I,J) = CS%Laplac3_const_xy(I,J) * vert_vort_mag + if (CS%Leith_Ah) Leith_Ah_q(I,J) = CS%biharm5_const_xx(I,J) * vert_vort_mag enddo ; enddo -end subroutine calc_vert_vort_mag +end subroutine calc_Leith_viscosity subroutine VarMix_init(Time, G, param_file, diag, CS) type(time_type), intent(in) :: Time !< Current model time @@ -862,6 +870,9 @@ subroutine VarMix_init(Time, G, param_file, diag, CS) ! value is roughly (pi / (the age of the universe) )^2. logical :: Gill_equatorial_Ld, use_FGNV_streamfn, use_MEKE, in_use real :: MLE_front_length + real :: Leith_Lap_const ! The non-dimensional coefficient in the Leith viscosity + real :: Leith_bi_const ! The non-dimensional coefficient in the bi-harmonic Leith viscosity + real :: DX2, DY2, grid_sp_2, grid_sp_3 ! Intermediate quantities for Leith metrics ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mdl = "MOM_lateral_mixing_coeffs" ! This module's name. @@ -1158,11 +1169,11 @@ subroutine VarMix_init(Time, G, param_file, diag, CS) "If true, add a term to Leith viscosity which is \n"//& "proportional to the gradient of divergence.", & default=.false.) - call get_param(param_file, mdl, "LEITH_LAP_CONST", CS%Leith_Lap_const, & + call get_param(param_file, mdl, "LEITH_LAP_CONST", Leith_Lap_const, & "The nondimensional Laplacian Leith constant, \n"//& "often set to 1.0", units="nondim", default=0.0, & fail_if_missing = CS%Leith_Kh) - call get_param(param_file, mdl, "LEITH_BI_CONST", CS%Leith_bi_const, & + call get_param(param_file, mdl, "LEITH_BI_CONST", Leith_bi_const, & "The nondimensional biharmonic Leith constant, \n"//& "typical values are thus far undetermined.", units="nondim", default=0.0, & fail_if_missing = CS%Leith_Ah) @@ -1176,6 +1187,42 @@ subroutine VarMix_init(Time, G, param_file, diag, CS) "is strongly encouraged, and no slip BCs are not used with \n"//& "the biharmonic viscosity.", default=.false.) endif + if (CS%Leith_Kh) then + allocate(CS%Laplac3_const_xx(isd:ied,jsd:jed)) ; CS%Laplac3_const_xx(:,:) = 0.0 + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + DX2 = G%dxT(i,j)*G%dxT(i,j) + DY2 = G%dyT(i,j)*G%dyT(i,j) + grid_sp_2 = (2.0*DX2*DY2) / (DX2 + DY2) + grid_sp_3 = grid_sp_2*sqrt(grid_sp_2) + CS%Laplac3_const_xx(i,j) = Leith_Lap_const * grid_sp_3 + enddo ; enddo + allocate(CS%Laplac3_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Laplac3_const_xy(:,:) = 0.0 + do J=js-1,Jeq ; do I=is-1,Ieq + DX2 = G%dxBu(I,J)*G%dxBu(I,J) + DY2 = G%dyBu(I,J)*G%dyBu(I,J) + grid_sp_2 = (2.0*DX2*DY2) / (DX2 + DY2) + grid_sp_3 = grid_sp_2*sqrt(grid_sp_2) + CS%Laplac3_const_xy(I,J) = Leith_Lap_const * grid_sp_3 + enddo ; enddo + endif + if (CS%Leith_Ah) then + allocate(CS%biharm5_const_xx(isd:ied,jsd:jed)) ; CS%biharm5_const_xx(:,:) = 0.0 + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + DX2 = G%dxT(i,j)*G%dxT(i,j) + DY2 = G%dyT(i,j)*G%dyT(i,j) + grid_sp_2 = (2.0*DX2*DY2) / (DX2+DY2) + grid_sp_3 = grid_sp_2*sqrt(grid_sp_2) + CS%biharm5_const_xx(i,j) = Leith_bi_const * (grid_sp_2 * grid_sp_3) + enddo ; enddo + allocate(CS%biharm5_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%biharm5_const_xy(:,:) = 0.0 + do J=js-1,Jeq ; do I=is-1,Ieq + DX2 = G%dxBu(I,J)*G%dxBu(I,J) + DY2 = G%dyBu(I,J)*G%dyBu(I,J) + grid_sp_2 = (2.0*DX2*DY2) / (DX2+DY2) + grid_sp_3 = grid_sp_2*sqrt(grid_sp_2) + CS%biharm5_const_xy(I,J) = Leith_bi_const * (grid_sp_2 * grid_sp_3) + enddo ; enddo + endif ! If nothing is being stored in this class then deallocate if (in_use) then From b39edff245773cdbd61e7afa02d99321e29f8e0e Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 5 Dec 2017 11:36:07 -0500 Subject: [PATCH 7/9] Cleaned up left over arrays from moving Leith out of MOM_hor_visc.F90 - Removed allocatable arrays and a left over parameter. --- .../lateral/MOM_hor_visc.F90 | 52 +------------------ 1 file changed, 2 insertions(+), 50 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index cc3359a070..22dd8f4b22 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -117,8 +117,6 @@ module MOM_hor_visc ! nonlinear eddy viscosity. AH is the background. logical :: Leith_Kh ! If true, use 2D Leith nonlinear eddy ! viscosity. KH is the background value. - logical :: Modified_Leith ! If true, use extra component of Leith viscosity - ! to damp divergent flow. To use, still set Leith_Kh=.TRUE. 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 @@ -184,15 +182,11 @@ module MOM_hor_visc ! 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) - Laplac3_Const_xx, & ! Laplacian metric-dependent constants (nondim) - Biharm5_Const_xx ! Biharmonic metric-dependent constants (nondim) + Biharm_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) - Laplac3_Const_xy, & ! Laplacian metric-dependent constants (nondim) - Biharm5_Const_xy ! Biharmonic metric-dependent constants (nondim) + Biharm_Const_xy ! Biharmonic metric-dependent constants (nondim) type(diag_ctrl), pointer :: diag ! structure to regulate diagnostic timing @@ -1012,7 +1006,6 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) 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. - CS%Modified_Leith = .false. Kh = 0.0 ; Ah = 0.0 @@ -1050,17 +1043,6 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) "If true, use a Leith nonlinear eddy viscosity.", & default=.false.) - call get_param(param_file, mdl, "MODIFIED_LEITH", CS%Modified_Leith, & - "If true, add a term to Leith viscosity which is \n"//& - "proportional to the gradient of divergence.", & - default=.false.) - - if (CS%Leith_Kh .or. get_all) & - call get_param(param_file, mdl, "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, mdl, "BOUND_KH", CS%bound_Kh, & "If true, the Laplacian coefficient is locally limited \n"//& "to be stable.", default=.true.) @@ -1124,13 +1106,6 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) endif endif - if (CS%Leith_Ah .or. get_all) then - call get_param(param_file, mdl, "LEITH_BI_CONST",Leith_bi_const, & - "The nondimensional biharmonic Leith constant, \n"//& - "typical values are thus far undetermined", units="nondim", default=0.0, & - fail_if_missing = CS%Leith_Ah) - endif - endif call get_param(param_file, mdl, "USE_LAND_MASK_FOR_HVISC", CS%use_land_mask, & @@ -1198,10 +1173,6 @@ 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 @@ -1239,10 +1210,6 @@ 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 @@ -1294,7 +1261,6 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) 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)) @@ -1309,7 +1275,6 @@ 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 CS%Kh_bg_xy(I,J) = MAX(Kh, Kh_vel_scale * sqrt(grid_sp_q2)) @@ -1352,9 +1317,6 @@ 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 @@ -1374,10 +1336,6 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS) 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) @@ -1533,9 +1491,6 @@ 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 @@ -1551,9 +1506,6 @@ 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 f08592f7727fcd19dbade006124970482e393c57 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Wed, 6 Dec 2017 13:59:51 -0500 Subject: [PATCH 8/9] Added independent flag for beta-term in Leith - USE_BETA_IN_LEITH now controls the beta-term in the Leith viscosity. --- .../lateral/MOM_lateral_mixing_coeffs.F90 | 20 ++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 764e393188..4298963ac7 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -111,6 +111,7 @@ module MOM_lateral_mixing_coeffs ! Leith parameters logical :: use_QG_Leith !< If true, enables the QG Leith scheme + logical :: use_beta_in_Leith !< If true, includes the beta term in the Leith viscosity logical :: Leith_Kh !< If true, enables the Leith scheme logical :: modified_Leith !< if true, include the divergence contribution to Leith viscosity logical :: Leith_Ah !< If true, enables the bi-harmonic Leith scheme @@ -820,8 +821,8 @@ subroutine calc_Leith_viscosity(CS, G, GV, u, v, h, k, Leith_Kh_h, Leith_Kh_q, L vort_xy_dy(I,j) = DX_dyBu * (vort_xy(I,J) * G%IdxCv(i,J) - vort_xy(I,J-1) * G%IdxCv(i,J-1)) enddo ; enddo - ! Add in beta for QG Leith - if (CS%use_QG_Leith) then + ! Add in beta for the Leith viscosity + if (CS%use_beta_in_Leith) then do J=js-2,Jeq+1 ; do I=is-1,Ieq+1 vort_xy_dx(i,J) = vort_xy_dx(i,J) + 0.5 * ( G%dF_dx(i,j) + G%dF_dx(i,j+1) ) enddo ; enddo @@ -1156,15 +1157,21 @@ subroutine VarMix_init(Time, G, param_file, diag, CS) endif ! Leith parameters + call get_param(param_file, mdl, "USE_QG_LEITH", CS%use_QG_Leith, & + "If true, use the QG Leith nonlinear eddy viscosity.", & + default=.false.) call get_param(param_file, mdl, "LEITH_KH", CS%Leith_Kh, & "If true, use a Leith nonlinear eddy viscosity.", & default=CS%use_QG_Leith) + if (CS%use_QG_Leith .and. .not. CS%Leith_Kh) call MOM_error(FATAL, & + "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//& + "LEITH_KH must be True when USE_QG_LEITH=True.") + call get_param(param_file, mdl, "USE_BETA_IN_LEITH", CS%use_beta_in_Leith, & + "If true, include the beta term in the QG Leith nonlinear eddy viscosity.", & + default=CS%use_QG_Leith) call get_param(param_file, mdl, "LEITH_AH", CS%Leith_Ah, & "If true, use a biharmonic Leith nonlinear eddy \n"//& "viscosity.", default=.false.) - call get_param(param_file, mdl, "USE_QG_LEITH", CS%use_QG_Leith, & - "If true, use the QG Leith nonlinear eddy viscosity.", & - default=.false.) call get_param(param_file, mdl, "MODIFIED_LEITH", CS%modified_Leith, & "If true, add a term to Leith viscosity which is \n"//& "proportional to the gradient of divergence.", & @@ -1186,6 +1193,9 @@ subroutine VarMix_init(Time, G, param_file, diag, CS) "cleaner than the no slip BCs. The use of free slip BCs \n"//& "is strongly encouraged, and no slip BCs are not used with \n"//& "the biharmonic viscosity.", default=.false.) + if (.not. CS%use_stored_slopes) call MOM_error(FATAL, & + "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//& + "USE_STORED_SLOPES must be True when ing Leith.") endif if (CS%Leith_Kh) then allocate(CS%Laplac3_const_xx(isd:ied,jsd:jed)) ; CS%Laplac3_const_xx(:,:) = 0.0 From 69f2c755e7acbed9c6f5ff016b4817c34fdf90f9 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Wed, 6 Dec 2017 15:46:36 -0500 Subject: [PATCH 9/9] Added stretching term for QG Leith viscosity - Nasty expressions for harmonic mean h's at u/v points for vertical derivatives of slopes... Eek! --- .../lateral/MOM_lateral_mixing_coeffs.F90 | 63 ++++++++++++++++--- 1 file changed, 55 insertions(+), 8 deletions(-) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 4298963ac7..2820ea7e45 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -413,12 +413,14 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, CS) if (.not. ASSOCIATED(CS)) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions:"//& "Module must be initialized before it is used.") - if (CS%calculate_Eady_growth_rate) then + if (CS%calculate_Eady_growth_rate .or. CS%use_stored_slopes) then call find_eta(h, tv, GV%g_Earth, G, GV, e, halo_size=2) if (CS%use_stored_slopes) then call calc_isoneutral_slopes(G, GV, h, e, tv, dt*CS%kappa_smooth, & CS%slope_x, CS%slope_y, N2_u, N2_v, 1) - call calc_Visbeck_coeffs(h, e, CS%slope_x, CS%slope_y, N2_u, N2_v, G, GV, CS) + if (CS%calculate_Eady_growth_rate) then + call calc_Visbeck_coeffs(h, e, CS%slope_x, CS%slope_y, N2_u, N2_v, G, GV, CS) + endif ! call calc_slope_functions_using_just_e(h, G, CS, e, .false.) else !call calc_isoneutral_slopes(G, GV, h, e, tv, dt*CS%kappa_smooth, CS%slope_x, CS%slope_y) @@ -759,13 +761,18 @@ subroutine calc_Leith_viscosity(CS, G, GV, u, v, h, k, Leith_Kh_h, Leith_Kh_q, L dvdx ! Zonal shear of meridional velocity (s-1) real, dimension(SZI_(G),SZJB_(G)) :: & vort_xy_dx, & ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) (m-1 s-1) - div_xx_dy ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) (m-1 s-1) + div_xx_dy, & ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) (m-1 s-1) + dslopey_dz, & ! z-derivative of y-slope at v-points (m-1) + h_at_v ! Thickness at v-points (m or kg m-2) real, dimension(SZIB_(G),SZJ_(G)) :: & vort_xy_dy, & ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) (m-1 s-1) - div_xx_dx ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) (m-1 s-1) + div_xx_dx, & ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) (m-1 s-1) + dslopex_dz, & ! z-derivative of x-slope at u-points (m-1) + h_at_u ! Thickness at u-points (m or kg m-2) real, dimension(SZI_(G),SZJ_(G)) :: div_xx ! Estimate of horizontal divergence at h-points (s-1) real :: mod_Leith, DY_dxBu, DX_dyBu, vert_vort_mag + real :: h_at_slope_above, h_at_slope_below, Ih, f integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -811,26 +818,66 @@ subroutine calc_Leith_viscosity(CS, G, GV, u, v, h, k, Leith_Kh_h, Leith_Kh_q, L endif ! Vorticity gradient - do J=js-2,Jeq+1 ; do I=is-1,Ieq+1 + do J=js-2,Jeq+1 ; do i=is-1,Ieq+1 DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J) vort_xy_dx(i,J) = DY_dxBu * (vort_xy(I,J) * G%IdyCu(I,j) - vort_xy(I-1,J) * G%IdyCu(I-1,j)) enddo ; enddo - do J=js-1,Jeq+1 ; do I=is-2,Ieq+1 + do j=js-1,Jeq+1 ; do I=is-2,Ieq+1 DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J) vort_xy_dy(I,j) = DX_dyBu * (vort_xy(I,J) * G%IdxCv(i,J) - vort_xy(I,J-1) * G%IdxCv(i,J-1)) enddo ; enddo ! Add in beta for the Leith viscosity if (CS%use_beta_in_Leith) then - do J=js-2,Jeq+1 ; do I=is-1,Ieq+1 + do J=js-2,Jeq+1 ; do i=is-1,Ieq+1 vort_xy_dx(i,J) = vort_xy_dx(i,J) + 0.5 * ( G%dF_dx(i,j) + G%dF_dx(i,j+1) ) enddo ; enddo - do J=js-1,Jeq+1 ; do I=is-2,Ieq+1 + do j=js-1,Jeq+1 ; do I=is-2,Ieq+1 vort_xy_dy(I,j) = vort_xy_dy(I,j) + 0.5 * ( G%dF_dy(i,j) + G%dF_dy(i+1,j) ) enddo ; enddo endif + ! Add in stretching term for the QG Leith vsicosity + if (CS%use_QG_Leith) then + do j=js-1,Jeq+1 ; do I=is-2,Ieq+1 + h_at_slope_above = 2. * ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) * h(i+1,j,k) ) / & + ( ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) + h(i+1,j,k) ) & + + ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k-1) + h(i+1,j,k-1) ) + GV%H_subroundoff ) + h_at_slope_below = 2. * ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) * h(i+1,j,k+1) ) / & + ( ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) + h(i+1,j,k+1) ) & + + ( h(i,j,k+1) * h(i+1,j,k+1) ) * ( h(i,j,k) + h(i+1,j,k) ) + GV%H_subroundoff ) + Ih = 1./ ( ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff ) * GV%H_to_m ) + dslopex_dz(I,j) = 2. * ( CS%slope_x(i,j,k) - CS%slope_x(i,j,k+1) ) * Ih + h_at_u(I,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih + enddo ; enddo + do J=js-2,Jeq+1 ; do i=is-1,Ieq+1 + h_at_slope_above = 2. * ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) * h(i,j+1,k) ) / & + ( ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) + h(i,j+1,k) ) & + + ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k-1) + h(i,j+1,k-1) ) + GV%H_subroundoff ) + h_at_slope_below = 2. * ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) * h(i,j+1,k+1) ) / & + ( ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) + h(i,j+1,k+1) ) & + + ( h(i,j,k+1) * h(i,j+1,k+1) ) * ( h(i,j,k) + h(i,j+1,k) ) + GV%H_subroundoff ) + Ih = 1./ ( ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff ) * GV%H_to_m ) + dslopey_dz(i,J) = 2. * ( CS%slope_y(i,j,k) - CS%slope_y(i,j,k+1) ) * Ih + h_at_v(i,J) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih + enddo ; enddo + do J=js-2,Jeq+1 ; do i=is-1,Ieq+1 + f = 0.5 * ( G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J) ) + vort_xy_dx(i,J) = vort_xy_dx(i,J) - f * & + ( ( h_at_u(I,j) * dslopex_dz(I,j) + h_at_u(I-1,j+1) * dslopex_dz(I-1,j+1) ) & + + ( h_at_u(I-1,j) * dslopex_dz(I-1,j) + h_at_u(I,j+1) * dslopex_dz(I,j+1) ) ) / & + ( ( h_at_u(I,j) + h_at_u(I-1,j+1) ) + ( h_at_u(I-1,j) + h_at_u(I,j+1) ) ) + enddo ; enddo + do j=js-1,Jeq+1 ; do I=is-2,Ieq+1 + f = 0.5 * ( G%CoriolisBu(I,J) + G%CoriolisBu(I,J-1) ) + vort_xy_dy(I,j) = vort_xy_dx(I,j) - f * & + ( ( h_at_v(i,J) * dslopey_dz(i,J) + h_at_v(i+1,J-1) * dslopey_dz(i+1,J-1) ) & + + ( h_at_v(i,J-1) * dslopey_dz(i,J-1) + h_at_v(i+1,J) * dslopey_dz(i+1,J) ) ) / & + ( ( h_at_v(i,J) + h_at_v(i+1,J-1) ) + ( h_at_v(i,J-1) + h_at_v(i+1,J) ) ) + enddo ; enddo + endif + mod_Leith = 0.; if (CS%modified_Leith) mod_Leith = 1.0 ! h-point viscosities