From aac71f58c53c186fad7a0937841ffea1141010ca Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 2 Dec 2019 18:59:26 +0000 Subject: [PATCH 1/5] Fix units for terms used in advection of MEKE - Code was passing dimensional testing because a conversion factor was included to obtain the right units for the advective flux terms. The origin of the unit problem was in the rate of inter-column exchange (baroHu and baroHv) which were in units of [H L2] and should have been in units of [R Z L2]: - [H L2 ~> m3 or kg] so differs between Boussinesq and non-Boussinesq; - [R Z L2 ~> kg] for both modes. - This commit moves the conversion to the summation for the barotropic exchange, removes the previous scaling factors, and updates the relevant comments. - Closes #1036 --- src/parameterizations/lateral/MOM_MEKE.F90 | 35 ++++++++++------------ 1 file changed, 15 insertions(+), 20 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 9513937c9d..9f99017882 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -131,17 +131,17 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h tmp ! Temporary variable for diagnostic computation real, dimension(SZIB_(G),SZJ_(G)) :: & - MEKE_uflux, & ! The zonal advective and diffusive flux of MEKE with different units in different - ! places of [L2 T-2 ~> m2 s-2] or [m L4 T-3 ~> m5 s-3] or [kg m-2 L4 T-3 ~> kg m-2 s-3]. + MEKE_uflux, & ! The zonal advective and diffusive flux of MEKE with units of [R Z L4 T-3 ~> kg m-2 s-3]. + ! In one place, MEKE_uflux is used as temporary work space with units of [L2 T-2 ~> m2 s-2]. Kh_u, & ! The zonal diffusivity that is actually used [L2 T-1 ~> m2 s-1]. - baroHu, & ! Depth integrated accumulated zonal mass flux [H L2 ~> m3 or kg]. + baroHu, & ! Depth integrated accumulated zonal mass flux [R Z L2 ~> kg]. drag_vel_u ! A (vertical) viscosity associated with bottom drag at ! u-points [Z T-1 ~> m s-1]. real, dimension(SZI_(G),SZJB_(G)) :: & - MEKE_vflux, & ! The meridional advective and diffusive flux of MEKE with different units in different - ! places of [L2 T-2 ~> m2 s-2] or [m L4 T-3 ~> m5 s-3] or [kg m-2 L4 T-3 ~> kg m-2 s-3]. + MEKE_vflux, & ! The meridional advective and diffusive flux of MEKE with units of [R Z L4 T-3 ~> kg m-2 s-3]. + ! In one place, MEKE_vflux is used as temporary work space with units of [L2 T-2 ~> m2 s-2]. Kh_v, & ! The meridional diffusivity that is actually used [L2 T-1 ~> m2 s-1]. - baroHv, & ! Depth integrated accumulated meridional mass flux [H L2 ~> m3 or kg]. + baroHv, & ! Depth integrated accumulated meridional mass flux [R Z L2 ~> kg]. drag_vel_v ! A (vertical) viscosity associated with bottom drag at ! v-points [Z T-1 ~> m s-1]. real :: Kh_here ! The local horizontal viscosity [L2 T-1 ~> m2 s-1] @@ -149,8 +149,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h real :: K4_here ! The local horizontal biharmonic viscosity [L4 T-1 ~> m4 s-1] real :: Inv_K4_max ! The inverse of the local horizontal biharmonic viscosity [T L-4 ~> s m-4] real :: cdrag2 - real :: advFac ! The product of the advection scaling factor and some unit conversion - ! factors divided by the timestep [m H-1 T-1 ~> s-1 or m3 kg-1 s-1] + real :: advFac ! The product of the advection scaling factor and 1/dt [T-1 ~> s-1] real :: mass_neglect ! A negligible mass [R Z ~> kg m-2]. real :: ldamping ! The MEKE damping rate [T-1 ~> s-1]. real :: Rho0 ! A density used to convert mass to distance [R ~> kg m-3]. @@ -201,14 +200,14 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h ! advisable to use Strang splitting between the damping and diffusion. sdt_damp = sdt ; if (CS%MEKE_KH >= 0.0 .or. CS%MEKE_K4 >= 0.) sdt_damp = 0.5*sdt - ! Calculate depth integrated mass flux if doing advection + ! Calculate depth integrated mass exchange if doing advection [R Z L2 ~> kg] if (CS%MEKE_advection_factor>0.) then do j=js,je ; do I=is-1,ie baroHu(I,j) = 0. enddo ; enddo do k=1,nz do j=js,je ; do I=is-1,ie - baroHu(I,j) = hu(I,j,k) + baroHu(I,j) = hu(I,j,k) * GV%H_to_RZ enddo ; enddo enddo do J=js-1,je ; do i=is,ie @@ -216,7 +215,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo ; enddo do k=1,nz do J=js-1,je ; do i=is,ie - baroHv(i,J) = hv(i,J,k) + baroHv(i,J) = hv(i,J,k) * GV%H_to_RZ enddo ; enddo enddo endif @@ -358,7 +357,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h ! Calculate Laplacian of MEKE !$OMP parallel do default(shared) do j=js-1,je+1 ; do I=is-2,ie+1 - ! Here the units of MEKE_uflux are [L2 T-2 ~> m2 s-2]. + ! MEKE_uflux is used here as workspace with units of [L2 T-2 ~> m2 s-2]. MEKE_uflux(I,j) = ((G%dy_Cu(I,j)*G%IdxCu(I,j)) * G%mask2dCu(I,j)) * & (MEKE%MEKE(i+1,j) - MEKE%MEKE(i,j)) ! This would have units of [R Z L2 T-2 ~> kg s-2] @@ -368,7 +367,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo ; enddo !$OMP parallel do default(shared) do J=js-2,je+1 ; do i=is-1,ie+1 - ! Here the units of MEKE_vflux are [L2 T-2 ~> m2 s-2]. + ! MEKE_vflux is used here as workspace with units of [L2 T-2 ~> m2 s-2]. MEKE_vflux(i,J) = ((G%dx_Cv(i,J)*G%IdyCv(i,J)) * G%mask2dCv(i,J)) * & (MEKE%MEKE(i,j+1) - MEKE%MEKE(i,j)) ! This would have units of [R Z L2 T-2 ~> kg s-2] @@ -457,12 +456,10 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h (MEKE%MEKE(i,j) - MEKE%MEKE(i,j+1)) enddo ; enddo if (CS%MEKE_advection_factor>0.) then - !### I think that for dimensional consistency, this should be: - ! advFac = GV%H_to_RZ * CS%MEKE_advection_factor / sdt - advFac = US%kg_m3_to_R*GV%H_to_Z * CS%MEKE_advection_factor / dt + advFac = CS%MEKE_advection_factor / sdt ! [T-1] !$OMP parallel do default(shared) do j=js,je ; do I=is-1,ie - ! Here the units of the quantities added to MEKE_uflux and MEKE_vflux are [m L4 T-3]. + ! Here the units of the quantities added to MEKE_uflux are [R Z L4 T-3 ~> kg m2 s-3]. if (baroHu(I,j)>0.) then MEKE_uflux(I,j) = MEKE_uflux(I,j) + baroHu(I,j)*MEKE%MEKE(i,j)*advFac elseif (baroHu(I,j)<0.) then @@ -471,7 +468,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo ; enddo !$OMP parallel do default(shared) do J=js-1,je ; do i=is,ie - ! Here the units of the quantities added to MEKE_uflux and MEKE_vflux are [m L4 T-3]. + ! Here the units of the quantities added to MEKE_vflux are [R Z L4 T-3 ~> kg m2 s-3]. if (baroHv(i,J)>0.) then MEKE_vflux(i,J) = MEKE_vflux(i,J) + baroHv(i,J)*MEKE%MEKE(i,j)*advFac elseif (baroHv(i,J)<0.) then @@ -480,10 +477,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo ; enddo endif - !$OMP parallel do default(shared) do j=js,je ; do i=is,ie - ! This expression is correct if the units of MEKE_uflux and MEKE_vflux are [kg m-2 L4 T-3]. MEKE%MEKE(i,j) = MEKE%MEKE(i,j) + (sdt*(G%IareaT(i,j)*I_mass(i,j))) * & ((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + & (MEKE_vflux(i,J-1) - MEKE_vflux(i,J))) From b7b7fbce954f5541cfd568811e6686b469c2032a Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Wed, 27 Nov 2019 21:19:16 +0000 Subject: [PATCH 2/5] Removed out-of-date commented code - Removed commented code snippet that was out-of-date w.r.t. to the current code and completely misleading. --- src/parameterizations/lateral/MOM_MEKE.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 9f99017882..5632f63956 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -380,8 +380,6 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h do j=js-1,je+1 ; do i=is-1,ie+1 del2MEKE(i,j) = G%IareaT(i,j) * & ((MEKE_uflux(I,j) - MEKE_uflux(I-1,j)) + (MEKE_vflux(i,J) - MEKE_vflux(i,J-1))) - ! del2MEKE(i,j) = (G%IareaT(i,j)*I_mass(i,j)) * & - ! ((MEKE_uflux(I,j) - MEKE_uflux(I-1,j)) + (MEKE_vflux(i,J) - MEKE_vflux(i,J-1))) enddo ; enddo ! Bi-harmonic diffusion of MEKE From 9c4d22d236838460634a686719fe630f425b3f2e Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Wed, 27 Nov 2019 21:21:54 +0000 Subject: [PATCH 3/5] Updated/added comments for units of variables in MOM_MEKE.F90 - Cleaned up/added comments annotating units. --- src/parameterizations/lateral/MOM_MEKE.F90 | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 5632f63956..eed783db98 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -265,7 +265,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo ; enddo do i=is-1,ie+1 I_mass(i,j) = 0.0 - if (mass(i,j) > 0.0) I_mass(i,j) = 1.0 / mass(i,j) + if (mass(i,j) > 0.0) I_mass(i,j) = 1.0 / mass(i,j) ! [m2 kg-1] enddo enddo @@ -354,7 +354,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h endif if (CS%MEKE_K4 >= 0.0) then - ! Calculate Laplacian of MEKE + ! Calculate Laplacian of MEKE using MEKE_uflux and MEKE_vflux as temporary work space. !$OMP parallel do default(shared) do j=js-1,je+1 ; do I=is-2,ie+1 ! MEKE_uflux is used here as workspace with units of [L2 T-2 ~> m2 s-2]. @@ -377,7 +377,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo ; enddo !$OMP parallel do default(shared) - do j=js-1,je+1 ; do i=is-1,ie+1 + do j=js-1,je+1 ; do i=is-1,ie+1 ! del2MEKE has units [T-2 ~> s-2]. del2MEKE(i,j) = G%IareaT(i,j) * & ((MEKE_uflux(I,j) - MEKE_uflux(I-1,j)) + (MEKE_vflux(i,J) - MEKE_vflux(i,J-1))) enddo ; enddo @@ -385,7 +385,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h ! Bi-harmonic diffusion of MEKE !$OMP parallel do default(shared) private(K4_here,Inv_K4_max) do j=js,je ; do I=is-1,ie - K4_here = CS%MEKE_K4 + K4_here = CS%MEKE_K4 ! [L4 T-1 ~> m4 s-1] ! Limit Kh to avoid CFL violations. Inv_K4_max = 64.0 * sdt * ((G%dy_Cu(I,j)*G%IdxCu(I,j)) * & max(G%IareaT(i,j), G%IareaT(i+1,j)))**2 @@ -398,15 +398,16 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo ; enddo !$OMP parallel do default(shared) private(K4_here,Inv_K4_max) do J=js-1,je ; do i=is,ie - K4_here = CS%MEKE_K4 + K4_here = CS%MEKE_K4 ! [L4 T-1 ~> m4 s-1] Inv_K4_max = 64.0 * sdt * ((G%dx_Cv(i,J)*G%IdyCv(i,J)) * max(G%IareaT(i,j), G%IareaT(i,j+1)))**2 if (K4_here*Inv_K4_max > 0.3) K4_here = 0.3 / Inv_K4_max + ! Here the units of MEKE_vflux are [kg m-2 L4 T-3]. MEKE_vflux(i,J) = ((K4_here * (G%dx_Cv(i,J)*G%IdyCv(i,J))) * & ((2.0*mass(i,j)*mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * & (del2MEKE(i,j+1) - del2MEKE(i,j)) enddo ; enddo - ! Store tendency arising from the bi-harmonic in del4MEKE + ! Store change in MEKE arising from the bi-harmonic in del4MEKE [L2 T-2]. !$OMP parallel do default(shared) do j=js,je ; do i=is,ie del4MEKE(i,j) = (sdt*(G%IareaT(i,j)*I_mass(i,j))) * & @@ -415,7 +416,6 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo ; enddo endif ! - if (CS%kh_flux_enabled) then ! Lateral diffusion of MEKE Kh_here = max(0., CS%MEKE_Kh) @@ -570,10 +570,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h endif ! Offer fields for averaging. - if (any([CS%id_Ue, CS%id_Ub, CS%id_Ut] > 0)) & tmp(:,:) = 0. - if (CS%id_MEKE>0) call post_data(CS%id_MEKE, MEKE%MEKE, CS%diag) if (CS%id_Ue>0) then do j=js,je ; do i=is,ie From 285273ac7a2758fe54913e69f0813fb863da48c6 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 2 Dec 2019 20:17:57 +0000 Subject: [PATCH 4/5] Corrected/added converted MKS units in comments --- src/parameterizations/lateral/MOM_MEKE.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index eed783db98..5e15b66ab6 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -261,11 +261,11 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h do j=js-1,je+1 do i=is-1,ie+1 ; mass(i,j) = 0.0 ; enddo do k=1,nz ; do i=is-1,ie+1 - mass(i,j) = mass(i,j) + G%mask2dT(i,j) * (GV%H_to_RZ * h(i,j,k)) + mass(i,j) = mass(i,j) + G%mask2dT(i,j) * (GV%H_to_RZ * h(i,j,k)) ! [R Z ~> kg m-2] enddo ; enddo do i=is-1,ie+1 I_mass(i,j) = 0.0 - if (mass(i,j) > 0.0) I_mass(i,j) = 1.0 / mass(i,j) ! [m2 kg-1] + if (mass(i,j) > 0.0) I_mass(i,j) = 1.0 / mass(i,j) ! [R-1 Z-1 ~> m2 kg-1] enddo enddo @@ -402,12 +402,12 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h Inv_K4_max = 64.0 * sdt * ((G%dx_Cv(i,J)*G%IdyCv(i,J)) * max(G%IareaT(i,j), G%IareaT(i,j+1)))**2 if (K4_here*Inv_K4_max > 0.3) K4_here = 0.3 / Inv_K4_max - ! Here the units of MEKE_vflux are [kg m-2 L4 T-3]. + ! Here the units of MEKE_vflux are [R Z L4 T-3]. MEKE_vflux(i,J) = ((K4_here * (G%dx_Cv(i,J)*G%IdyCv(i,J))) * & ((2.0*mass(i,j)*mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * & (del2MEKE(i,j+1) - del2MEKE(i,j)) enddo ; enddo - ! Store change in MEKE arising from the bi-harmonic in del4MEKE [L2 T-2]. + ! Store change in MEKE arising from the bi-harmonic in del4MEKE [L2 T-2 ~> m2 s-2]. !$OMP parallel do default(shared) do j=js,je ; do i=is,ie del4MEKE(i,j) = (sdt*(G%IareaT(i,j)*I_mass(i,j))) * & @@ -454,7 +454,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h (MEKE%MEKE(i,j) - MEKE%MEKE(i,j+1)) enddo ; enddo if (CS%MEKE_advection_factor>0.) then - advFac = CS%MEKE_advection_factor / sdt ! [T-1] + advFac = CS%MEKE_advection_factor / sdt ! [T-1 ~> s-1] !$OMP parallel do default(shared) do j=js,je ; do I=is-1,ie ! Here the units of the quantities added to MEKE_uflux are [R Z L4 T-3 ~> kg m2 s-3]. From 7a404fe27321bbff7d3253dbdd67722b2a9773e9 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 2 Dec 2019 21:17:06 +0000 Subject: [PATCH 5/5] Added more converted MKS units in comments --- src/parameterizations/lateral/MOM_MEKE.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 5e15b66ab6..38bf24ee60 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -391,7 +391,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h max(G%IareaT(i,j), G%IareaT(i+1,j)))**2 if (K4_here*Inv_K4_max > 0.3) K4_here = 0.3 / Inv_K4_max - ! Here the units of MEKE_uflux are [R Z L4 T-3]. + ! Here the units of MEKE_uflux are [R Z L4 T-3 ~> kg m2 s-3]. MEKE_uflux(I,j) = ((K4_here * (G%dy_Cu(I,j)*G%IdxCu(I,j))) * & ((2.0*mass(i,j)*mass(i+1,j)) / ((mass(i,j)+mass(i+1,j)) + mass_neglect)) ) * & (del2MEKE(i+1,j) - del2MEKE(i,j)) @@ -402,7 +402,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h Inv_K4_max = 64.0 * sdt * ((G%dx_Cv(i,J)*G%IdyCv(i,J)) * max(G%IareaT(i,j), G%IareaT(i,j+1)))**2 if (K4_here*Inv_K4_max > 0.3) K4_here = 0.3 / Inv_K4_max - ! Here the units of MEKE_vflux are [R Z L4 T-3]. + ! Here the units of MEKE_vflux are [R Z L4 T-3 ~> kg m2 s-3]. MEKE_vflux(i,J) = ((K4_here * (G%dx_Cv(i,J)*G%IdyCv(i,J))) * & ((2.0*mass(i,j)*mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * & (del2MEKE(i,j+1) - del2MEKE(i,j))