From a4f9550f909ffd93c952751853b8d4f93e4668f0 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 19 Sep 2019 16:50:12 -0600 Subject: [PATCH 1/7] Added new equilibrium formula for MEKE * Follow equation 1 of Jansen et al. (2015), balancing the GEOMETRIC GM coefficient against bottom drag (Equations 3 and 12); * Added limited for SN in this formula, to avoid extremely large values. TODO: * Increase GEOMETRIC_ALPHA in this calculation * Use GEOMETRIC_EPSILON as a limiter for SN --- src/parameterizations/lateral/MOM_MEKE.F90 | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 3688c3dfea..f24d549970 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -45,6 +45,8 @@ module MOM_MEKE logical :: visc_drag !< If true use the vertvisc_type to calculate bottom drag. logical :: MEKE_GEOMETRIC !< If true, uses the GM coefficient formulation from the GEOMETRIC !! framework (Marshall et al., 2012) + real :: MEKE_GEOMETRIC_alpha !< The nondimensional coefficient governing the efficiency of the + !! GEOMETRIC thickness diffusion. logical :: MEKE_equilibrium_alt !< If true, use an alternative calculation for the !! equilibrium value of MEKE. logical :: GM_src_alt !< If true, use the GM energy conversion form S^2*N^2*kappa rather @@ -747,7 +749,13 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m EKE = 0. endif if (CS%MEKE_equilibrium_alt) then - MEKE%MEKE(i,j) = (US%Z_to_m*G%bathyT(i,j)*SN / (8*CS%cdrag))**2 + if (CS%MEKE_GEOMETRIC) then + ! Equation 1 of Jansen et al. (2015), balancing the GEOMETRIC GM coefficient against + ! bottom drag (Equations 3 and 12) + MEKE%MEKE(i,j) = (CS%MEKE_GEOMETRIC_alpha * MIN(SN,1.0e-7))**2 / ((I_H * CS%cdrag)**2 * (bottomFac2**3)) + else + MEKE%MEKE(i,j) = (US%Z_to_m*G%bathyT(i,j)*SN / (8*CS%cdrag))**2 + endif else MEKE%MEKE(i,j) = EKE endif @@ -978,6 +986,9 @@ logical function MEKE_init(Time, G, param_file, diag, CS, MEKE, restart_CS) call get_param(param_file, mdl, "MEKE_GEOMETRIC", CS%MEKE_GEOMETRIC, & "If MEKE_GEOMETRIC is true, uses the GM coefficient formulation "//& "from the GEOMETRIC framework (Marshall et al., 2012).", default=.false.) + call get_param(param_file, mdl, "MEKE_GEOMETRIC_ALPHA", CS%MEKE_GEOMETRIC_alpha, & + "The nondimensional coefficient governing the efficiency of the GEOMETRIC \n"//& + "thickness diffusion.", units="nondim", default=0.05) call get_param(param_file, mdl, "MEKE_EQUILIBRIUM_ALT", CS%MEKE_equilibrium_alt, & "If true, use an alternative formula for computing the (equilibrium)"//& "initial value of MEKE.", default=.false.) From bb46c38d8f7de49fcdcc8944dd34f8ede67049ff Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 23 Sep 2019 14:59:25 -0600 Subject: [PATCH 2/7] Calculates bottomFac2 IF CS%MEKE_GEOMETRIC=True In this commit, bottomFac2 is calculated when CS%MEKE_GEOMETRIC is set to true. Previously, bottomFac2 was calculated in MEKE_lengthScales_0d but something else on that subroutine was returning Nan so we decided to pull out just the bottomFac2 calculation from that. --- src/parameterizations/lateral/MOM_MEKE.F90 | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index f24d549970..645dcc5e8a 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -623,6 +623,8 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m real, parameter :: tolerance = 1.e-12 ! Width of EKE bracket [m2 s-2]. logical :: useSecant, debugIteration + real :: Lgrid, Ldeform, Lfrict + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec debugIteration = .false. @@ -750,9 +752,18 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m endif if (CS%MEKE_equilibrium_alt) then if (CS%MEKE_GEOMETRIC) then + Lgrid = sqrt(G%areaT(i,j)) ! Grid scale + Ldeform =Lgrid * MIN(1.0,MEKE%Rd_dx_h(i,j)) ! Deformation scale + Lfrict = (US%Z_to_m * G%bathyT(i,j)) / CS%cdrag ! Frictional arrest scale + ! gamma_b^2 is the ratio of bottom eddy energy to mean column eddy energy + ! used in calculating bottom drag + bottomFac2 = CS%MEKE_CD_SCALE**2 + if (Lfrict*CS%MEKE_Cb>0.) bottomFac2 = bottomFac2 + 1./( 1. + CS%MEKE_Cb*(Ldeform/Lfrict) )**0.8 + bottomFac2 = max(bottomFac2, CS%MEKE_min_gamma) ! Equation 1 of Jansen et al. (2015), balancing the GEOMETRIC GM coefficient against ! bottom drag (Equations 3 and 12) - MEKE%MEKE(i,j) = (CS%MEKE_GEOMETRIC_alpha * MIN(SN,1.0e-7))**2 / ((I_H * CS%cdrag)**2 * (bottomFac2**3)) + ! TODO: create a run time parameter for limitting SN. + MEKE%MEKE(i,j) = (CS%MEKE_GEOMETRIC_alpha * MIN(SN,1.e-5) * US%Z_to_m*G%bathyT(i,j))**2 / (CS%cdrag**2 * bottomFac2**3) else MEKE%MEKE(i,j) = (US%Z_to_m*G%bathyT(i,j)*SN / (8*CS%cdrag))**2 endif From 15c9f06463f1db23afd383b653f26dd10255d50a Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 2 Oct 2019 15:39:11 -0600 Subject: [PATCH 3/7] Hard-code min of SN to be 1.0e-7 --- src/parameterizations/lateral/MOM_MEKE.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index bb825cdd2d..a64a73ae4c 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -775,7 +775,7 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m ! Equation 1 of Jansen et al. (2015), balancing the GEOMETRIC GM coefficient against ! bottom drag (Equations 3 and 12) ! TODO: create a run time parameter for limitting SN. - MEKE%MEKE(i,j) = (CS%MEKE_GEOMETRIC_alpha * MIN(SN,1.e-5) * US%Z_to_m*G%bathyT(i,j))**2 / (CS%cdrag**2 * bottomFac2**3) + MEKE%MEKE(i,j) = (CS%MEKE_GEOMETRIC_alpha * MIN(SN,1.e-7) * US%Z_to_m*G%bathyT(i,j))**2 / (CS%cdrag**2 * bottomFac2**3) else MEKE%MEKE(i,j) = (US%Z_to_m*G%bathyT(i,j)*SN / (8*CS%cdrag))**2 endif From 839217d9d26f97580a6586697aed7de36ad440cc Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 15 Oct 2019 15:35:02 -0600 Subject: [PATCH 4/7] Rearranged MEKE_EQUILIBRIUM subroutine * Moved MEKE_equilibrium_alt toward top of subroutine to avoid unnecessary calculations. --- src/parameterizations/lateral/MOM_MEKE.F90 | 221 ++++++++++----------- 1 file changed, 103 insertions(+), 118 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index a64a73ae4c..d6ec7814ce 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -658,128 +658,113 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m ! This avoids extremes values in equilibrium solution due to bad values in SN_u, SN_v SN = min(SN_u(I,j), SN_u(I-1,j), SN_v(i,J), SN_v(i,J-1)) - FatH = 0.25*((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + & - (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))) ! Coriolis parameter at h points - - ! Since zero-bathymetry cells are masked, this avoids calculations on land - if (CS%MEKE_topographic_beta == 0. .or. G%bathyT(i,j) == 0.) then - beta_topo_x = 0. ; beta_topo_y = 0. + if (CS%MEKE_equilibrium_alt) then + MEKE%MEKE(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_m*G%bathyT(i,j))**2 / cd2 else - !### Consider different combinations of these estimates of topographic beta, and the use - ! of the water column thickness instead of the bathymetric depth. - beta_topo_x = CS%MEKE_topographic_beta * FatH * 0.5 * ( & - (G%bathyT(i+1,j)-G%bathyT(i,j)) * G%IdxCu(I,j) & - / max(G%bathyT(i+1,j),G%bathyT(i,j), GV%H_subroundoff) & - + (G%bathyT(i,j)-G%bathyT(i-1,j)) * G%IdxCu(I-1,j) & - / max(G%bathyT(i,j),G%bathyT(i-1,j), GV%H_subroundoff) ) - !### There is a bug in the 4th lne below, where IdxCu should be IdyCv. - beta_topo_y = CS%MEKE_topographic_beta * FatH * 0.5 * ( & - (G%bathyT(i,j+1)-G%bathyT(i,j)) * G%IdyCv(i,J) & - / max(G%bathyT(i,j+1),G%bathyT(i,j), GV%H_subroundoff) + & - (G%bathyT(i,j)-G%bathyT(i,j-1)) * G%IdxCu(i,J-1) & - / max(G%bathyT(i,j),G%bathyT(i,j-1), GV%H_subroundoff) ) - endif - beta = sqrt((G%dF_dx(i,j) - beta_topo_x)**2 + & - (G%dF_dy(i,j) - beta_topo_y)**2 ) - - I_H = US%L_to_m*GV%Rho0 * I_mass(i,j) - - if (KhCoeff*SN*I_H>0.) then - ! Solve resid(E) = 0, where resid = Kh(E) * (SN)^2 - damp_rate(E) E - EKEmin = 0. ! Use the trivial root as the left bracket - ResMin = 0. ! Need to detect direction of left residual - EKEmax = 0.01*US%m_s_to_L_T**2 ! First guess at right bracket - useSecant = .false. ! Start using a bisection method - - ! First find right bracket for which resid<0 - resid = 1.0*US%m_to_L**2*US%T_to_s**3 ; n1 = 0 - do while (resid>0.) - n1 = n1 + 1 - EKE = EKEmax - call MEKE_lengthScales_0d(CS, US, G%areaT(i,j), beta, G%bathyT(i,j), & - MEKE%Rd_dx_h(i,j), SN, EKE, & - bottomFac2, barotrFac2, LmixScale, LRhines, LEady) - ! TODO: Should include resolution function in Kh - Kh = (KhCoeff * sqrt(2.*barotrFac2*EKE) * LmixScale) - src = Kh * (SN * SN) - drag_rate = I_H * sqrt(drag_rate_visc(i,j)**2 + cd2 * ( 2.0*bottomFac2*EKE + Ubg2 ) ) - ldamping = CS%MEKE_damping + drag_rate * bottomFac2 - resid = src - ldamping * EKE - ! if (debugIteration) then - ! write(0,*) n1, 'EKE=',EKE,'resid=',resid - ! write(0,*) 'EKEmin=',EKEmin,'ResMin=',ResMin - ! write(0,*) 'src=',src,'ldamping=',ldamping - ! write(0,*) 'gamma-b=',bottomFac2,'gamma-t=',barotrFac2 - ! write(0,*) 'drag_visc=',drag_rate_visc(i,j),'Ubg2=',Ubg2 - ! endif - if (resid>0.) then ! EKE is to the left of the root - EKEmin = EKE ! so we move the left bracket here - EKEmax = 10. * EKE ! and guess again for the right bracket - if (resid 2.e17) then - if (debugIteration) stop 'Something has gone very wrong' - debugIteration = .true. - resid = 1. ; n1 = 0 - EKEmin = 0. ; ResMin = 0. - EKEmax = 0.01*US%m_s_to_L_T**2 - useSecant = .false. - endif - endif - enddo ! while(resid>0.) searching for right bracket - ResMax = resid - - ! Bisect the bracket - n2 = 0 ; EKEerr = EKEmax - EKEmin - do while (US%L_T_to_m_s**2*EKEerr>tolerance) - n2 = n2 + 1 - if (useSecant) then - EKE = EKEmin + (EKEmax - EKEmin) * (ResMin / (ResMin - ResMax)) - else - EKE = 0.5 * (EKEmin + EKEmax) - endif - EKEerr = min( EKE-EKEmin, EKEmax-EKE ) - ! TODO: Should include resolution function in Kh - Kh = (KhCoeff * sqrt(2.*barotrFac2*EKE) * LmixScale) - src = Kh * (SN * SN) - drag_rate = I_H * sqrt( drag_rate_visc(i,j)**2 + cd2 * ( 2.0*bottomFac2*EKE + Ubg2 ) ) - ldamping = CS%MEKE_damping + drag_rate * bottomFac2 - resid = src - ldamping * EKE - if (useSecant .and. resid>ResMin) useSecant = .false. - if (resid>0.) then ! EKE is to the left of the root - EKEmin = EKE ! so we move the left bracket here - if (resid EKE is exactly at the root - endif - if (n2>200) stop 'Failing to converge?' - enddo ! while(EKEmax-EKEmin>tolerance) - else - EKE = 0. - endif - if (CS%MEKE_equilibrium_alt) then - if (CS%MEKE_GEOMETRIC) then - Lgrid = sqrt(G%areaT(i,j)) ! Grid scale - Ldeform =Lgrid * MIN(1.0,MEKE%Rd_dx_h(i,j)) ! Deformation scale - Lfrict = (US%Z_to_m * G%bathyT(i,j)) / CS%cdrag ! Frictional arrest scale - ! gamma_b^2 is the ratio of bottom eddy energy to mean column eddy energy - ! used in calculating bottom drag - bottomFac2 = CS%MEKE_CD_SCALE**2 - if (Lfrict*CS%MEKE_Cb>0.) bottomFac2 = bottomFac2 + 1./( 1. + CS%MEKE_Cb*(Ldeform/Lfrict) )**0.8 - bottomFac2 = max(bottomFac2, CS%MEKE_min_gamma) - ! Equation 1 of Jansen et al. (2015), balancing the GEOMETRIC GM coefficient against - ! bottom drag (Equations 3 and 12) - ! TODO: create a run time parameter for limitting SN. - MEKE%MEKE(i,j) = (CS%MEKE_GEOMETRIC_alpha * MIN(SN,1.e-7) * US%Z_to_m*G%bathyT(i,j))**2 / (CS%cdrag**2 * bottomFac2**3) + FatH = 0.25*((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + & + (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))) ! Coriolis parameter at h points + + ! Since zero-bathymetry cells are masked, this avoids calculations on land + if (CS%MEKE_topographic_beta == 0. .or. G%bathyT(i,j) == 0.) then + beta_topo_x = 0. ; beta_topo_y = 0. + else + !### Consider different combinations of these estimates of topographic beta, and the use + ! of the water column thickness instead of the bathymetric depth. + beta_topo_x = CS%MEKE_topographic_beta * FatH * 0.5 * ( & + (G%bathyT(i+1,j)-G%bathyT(i,j)) * G%IdxCu(I,j) & + / max(G%bathyT(i+1,j),G%bathyT(i,j), GV%H_subroundoff) & + + (G%bathyT(i,j)-G%bathyT(i-1,j)) * G%IdxCu(I-1,j) & + / max(G%bathyT(i,j),G%bathyT(i-1,j), GV%H_subroundoff) ) + !### There is a bug in the 4th lne below, where IdxCu should be IdyCv. + beta_topo_y = CS%MEKE_topographic_beta * FatH * 0.5 * ( & + (G%bathyT(i,j+1)-G%bathyT(i,j)) * G%IdyCv(i,J) & + / max(G%bathyT(i,j+1),G%bathyT(i,j), GV%H_subroundoff) + & + (G%bathyT(i,j)-G%bathyT(i,j-1)) * G%IdxCu(i,J-1) & + / max(G%bathyT(i,j),G%bathyT(i,j-1), GV%H_subroundoff) ) + endif + beta = sqrt((G%dF_dx(i,j) - beta_topo_x)**2 + & + (G%dF_dy(i,j) - beta_topo_y)**2 ) + + I_H = US%L_to_m*GV%Rho0 * I_mass(i,j) + + if (KhCoeff*SN*I_H>0.) then + ! Solve resid(E) = 0, where resid = Kh(E) * (SN)^2 - damp_rate(E) E + EKEmin = 0. ! Use the trivial root as the left bracket + ResMin = 0. ! Need to detect direction of left residual + EKEmax = 0.01*US%m_s_to_L_T**2 ! First guess at right bracket + useSecant = .false. ! Start using a bisection method + + ! First find right bracket for which resid<0 + resid = 1.0*US%m_to_L**2*US%T_to_s**3 ; n1 = 0 + do while (resid>0.) + n1 = n1 + 1 + EKE = EKEmax + call MEKE_lengthScales_0d(CS, US, G%areaT(i,j), beta, G%bathyT(i,j), & + MEKE%Rd_dx_h(i,j), SN, EKE, & + bottomFac2, barotrFac2, LmixScale, LRhines, LEady) + ! TODO: Should include resolution function in Kh + Kh = (KhCoeff * sqrt(2.*barotrFac2*EKE) * LmixScale) + src = Kh * (SN * SN) + drag_rate = I_H * sqrt(drag_rate_visc(i,j)**2 + cd2 * ( 2.0*bottomFac2*EKE + Ubg2 ) ) + ldamping = CS%MEKE_damping + drag_rate * bottomFac2 + resid = src - ldamping * EKE + ! if (debugIteration) then + ! write(0,*) n1, 'EKE=',EKE,'resid=',resid + ! write(0,*) 'EKEmin=',EKEmin,'ResMin=',ResMin + ! write(0,*) 'src=',src,'ldamping=',ldamping + ! write(0,*) 'gamma-b=',bottomFac2,'gamma-t=',barotrFac2 + ! write(0,*) 'drag_visc=',drag_rate_visc(i,j),'Ubg2=',Ubg2 + ! endif + if (resid>0.) then ! EKE is to the left of the root + EKEmin = EKE ! so we move the left bracket here + EKEmax = 10. * EKE ! and guess again for the right bracket + if (resid 2.e17) then + if (debugIteration) stop 'Something has gone very wrong' + debugIteration = .true. + resid = 1. ; n1 = 0 + EKEmin = 0. ; ResMin = 0. + EKEmax = 0.01*US%m_s_to_L_T**2 + useSecant = .false. + endif + endif + enddo ! while(resid>0.) searching for right bracket + ResMax = resid + + ! Bisect the bracket + n2 = 0 ; EKEerr = EKEmax - EKEmin + do while (US%L_T_to_m_s**2*EKEerr>tolerance) + n2 = n2 + 1 + if (useSecant) then + EKE = EKEmin + (EKEmax - EKEmin) * (ResMin / (ResMin - ResMax)) + else + EKE = 0.5 * (EKEmin + EKEmax) + endif + EKEerr = min( EKE-EKEmin, EKEmax-EKE ) + ! TODO: Should include resolution function in Kh + Kh = (KhCoeff * sqrt(2.*barotrFac2*EKE) * LmixScale) + src = Kh * (SN * SN) + drag_rate = I_H * sqrt( drag_rate_visc(i,j)**2 + cd2 * ( 2.0*bottomFac2*EKE + Ubg2 ) ) + ldamping = CS%MEKE_damping + drag_rate * bottomFac2 + resid = src - ldamping * EKE + if (useSecant .and. resid>ResMin) useSecant = .false. + if (resid>0.) then ! EKE is to the left of the root + EKEmin = EKE ! so we move the left bracket here + if (resid EKE is exactly at the root + endif + if (n2>200) stop 'Failing to converge?' + enddo ! while(EKEmax-EKEmin>tolerance) else - MEKE%MEKE(i,j) = (US%Z_to_m*G%bathyT(i,j)*SN / (8*CS%cdrag))**2 + EKE = 0. endif - else MEKE%MEKE(i,j) = EKE endif enddo ; enddo From ebf5ee0d37f1e57dfaed0c56492369dfd0a1249d Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 15 Oct 2019 16:42:37 -0600 Subject: [PATCH 5/7] Adds MEKE_equilibrium_restoring This commit adds a subroutine that calculates a new equilibrium value for MEKE at each time step. This is not copied into MEKE%MEKE; rather, it is used as a restoring term to nudge MEKE%MEKE back to an equilibrium value. To select this option one needs to set MEKE_EQUILIBRIUM_RESTORING=True. The timescale for nudging is controlled via MEKE_RESTORING_TIMESCALE. --- src/parameterizations/lateral/MOM_MEKE.F90 | 65 ++++++++++++++++++++++ 1 file changed, 65 insertions(+) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index d6ec7814ce..853f3a8613 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -30,6 +30,8 @@ module MOM_MEKE !> Control structure that contains MEKE parameters and diagnostics handles type, public :: MEKE_CS ; private ! Parameters + real, dimension(:,:), pointer :: equilibrium_value => NULL() !< The equilbrium value + !! of MEKE to be calculated at each time step [L2 T-2 ~> m2 s-2] real :: MEKE_FrCoeff !< Efficiency of conversion of ME into MEKE [nondim] real :: MEKE_GMcoeff !< Efficiency of conversion of PE into MEKE [nondim] real :: MEKE_GMECoeff !< Efficiency of conversion of MEKE into ME by GME [nondim] @@ -47,6 +49,8 @@ module MOM_MEKE !! GEOMETRIC thickness diffusion. logical :: MEKE_equilibrium_alt !< If true, use an alternative calculation for the !! equilibrium value of MEKE. + logical :: MEKE_equilibrium_restoring !< If true, restore MEKE back to its equilibrium value, + !! which is calculated at each time step. logical :: GM_src_alt !< If true, use the GM energy conversion form S^2*N^2*kappa rather !! than the streamfunction for the MEKE GM source term. logical :: Rd_as_max_scale !< If true the length scale can not exceed the @@ -77,6 +81,8 @@ module MOM_MEKE real :: MEKE_advection_factor !< A scaling in front of the advection of MEKE [nondim] real :: MEKE_topographic_beta !< Weight for how much topographic beta is considered !! when computing beta in Rhines scale [nondim] + real :: MEKE_restoring_rate !< Inverse of the timescale used to nudge MEKE toward its equilibrium value [s-1]. + logical :: kh_flux_enabled !< If true, lateral diffusive MEKE flux is enabled. logical :: initialize !< If True, invokes a steady state solver to calculate MEKE. logical :: debug !< If true, write out checksums of data for debugging @@ -89,6 +95,7 @@ module MOM_MEKE integer :: id_KhMEKE_u = -1, id_KhMEKE_v = -1, id_Ku = -1, id_Au = -1 integer :: id_Le = -1, id_gamma_b = -1, id_gamma_t = -1 integer :: id_Lrhines = -1, id_Leady = -1 + integer :: id_MEKE_equilibrium = -1 !!@} ! Infrastructure @@ -325,6 +332,13 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h endif endif + if (CS%MEKE_equilibrium_restoring) then + call MEKE_equilibrium_restoring(CS, MEKE, G, GV, US, SN_u, SN_v) + do j=js,je ; do i=is,ie + src(i,j) = src(i,j) - CS%MEKE_restoring_rate*(MEKE%MEKE(i,j) - CS%equilibrium_value(i,j)) + enddo ; enddo + endif + ! Increase EKE by a full time-steps worth of source !$OMP parallel do default(shared) do j=js,je ; do i=is,ie @@ -772,6 +786,38 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m end subroutine MEKE_equilibrium +!< This subroutine calculates a new equilibrium value for MEKE at each time step. This is not copied into +!! MEKE%MEKE; rather, it is used as a restoring term to nudge MEKE%MEKE back to an equilibrium value +subroutine MEKE_equilibrium_restoring(CS, MEKE, G, GV, US, SN_u, SN_v) + type(ocean_grid_type), intent(inout) :: G !< Ocean grid. + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(MEKE_CS), pointer :: CS !< MEKE control structure. + type(MEKE_type), pointer :: MEKE !< A structure with MEKE data. + real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1]. + real, dimension(SZI_(G),SZJB_(G)), intent(in) :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1]. + ! Local variables + real :: SN ! The local Eady growth rate [T-1 ~> s-1] + integer :: i, j, is, ie, js, je, n1, n2 + real :: cd2 + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + cd2 = CS%cdrag**2 + +!$OMP do + do j=js,je ; do i=is,ie + ! SN = 0.25*max( (SN_u(I,j) + SN_u(I-1,j)) + (SN_v(i,J) + SN_v(i,J-1)), 0.) + ! This avoids extremes values in equilibrium solution due to bad values in SN_u, SN_v + SN = min(SN_u(I,j), SN_u(I-1,j), SN_v(i,J), SN_v(i,J-1)) + + CS%equilibrium_value(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_m*G%bathyT(i,j))**2 / cd2 + enddo ; enddo + + if (CS%id_MEKE_equilibrium>0) call post_data(CS%id_MEKE_equilibrium, CS%equilibrium_value, CS%diag) + +end subroutine MEKE_equilibrium_restoring + + !> Calculates the eddy mixing length scale and \f$\gamma_b\f$ and \f$\gamma_t\f$ !! functions that are ratios of either bottom or barotropic eddy energy to the !! column eddy energy, respectively. See \ref section_MEKE_equations. @@ -937,6 +983,7 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS) ! run to the representation in a restart file. real :: L_rescale ! A rescaling factor for length from the internal representation in this ! run to the representation in a restart file. + real :: MEKE_restoring_timescale ! The timescale used to nudge MEKE toward its equilibrium value. integer :: i, j, is, ie, js, je, isd, ied, jsd, jed logical :: laplacian, biharmonic, useVarMix, coldStart ! This include declares and sets the variable "version". @@ -1002,6 +1049,19 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS) call get_param(param_file, mdl, "MEKE_EQUILIBRIUM_ALT", CS%MEKE_equilibrium_alt, & "If true, use an alternative formula for computing the (equilibrium)"//& "initial value of MEKE.", default=.false.) + if (CS%MEKE_equilibrium_alt) then + call get_param(param_file, mdl, "MEKE_EQUILIBRIUM_RESTORING", CS%MEKE_equilibrium_restoring, & + "If true, restore MEKE back to its equilibrium value, which is calculated at"//& + "each time step.", default=.false.) + if (CS%MEKE_equilibrium_restoring) then + call get_param(param_file, mdl, "MEKE_RESTORING_TIMESCALE", MEKE_restoring_timescale, & + "The timescale used to nudge MEKE toward its equilibrium value.", units="s", & + default=1e6, scale=US%T_to_s) + allocate(CS%equilibrium_value(isd:ied,jsd:jed)) ; CS%equilibrium_value(:,:) = 0.0 + CS%MEKE_restoring_rate = 1.0 / MEKE_restoring_timescale + endif + + endif call get_param(param_file, mdl, "MEKE_FRCOEFF", CS%MEKE_FrCoeff, & "The efficiency of the conversion of mean energy into "//& "MEKE. If MEKE_FRCOEFF is negative, this conversion "//& @@ -1193,6 +1253,11 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS) 'Meridional diffusivity of MEKE', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T) endif + if (CS%MEKE_equilibrium_restoring) then + CS%id_MEKE_equilibrium = register_diag_field('ocean_model', 'MEKE_equilibrium', diag%axesT1, Time, & + 'Equilibrated Mesoscale Eddy Kinetic Energy', 'm2 s-2', conversion=US%L_T_to_m_s**2) + endif + CS%id_clock_pass = cpu_clock_id('(Ocean continuity halo updates)', grain=CLOCK_ROUTINE) ! Detect whether this instance of MEKE_init() is at the beginning of a run From 3f041d93fbb784b420456ca2e6b60df647133426 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 3 Oct 2019 18:24:14 -0400 Subject: [PATCH 6/7] MEKE diagnostic array fixes This patch fixes the following MEKE diagnostics: - MEKE_Ue, MEKE_Ub, MEKE_Ut The diagnostics were computed as inline operations inside post_data, e.g.: post_data(..., sqrt(0, max(0., MEKE*bottomFac2))) rather than computing the fields explicitly inside of array loops. This case causing floating point exceptions in Intel compilers, possibly likely due to evaluations inside of halos. We resolve these diagnostics by computing the values into a scratch array which is then passed to post_data. --- src/parameterizations/lateral/MOM_MEKE.F90 | 28 ++++++++++++++++++---- 1 file changed, 24 insertions(+), 4 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 853f3a8613..8d3e78262c 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -139,7 +139,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h del4MEKE, & ! Time-integrated MEKE tendency arising from the biharmonic of MEKE [L2 T-2 ~> m2 s-2]. LmixScale, & ! Eddy mixing length [L ~> m]. barotrFac2, & ! Ratio of EKE_barotropic / EKE [nondim] - bottomFac2 ! Ratio of EKE_bottom / EKE [nondim] + bottomFac2, & ! Ratio of EKE_bottom / EKE [nondim] + 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 @@ -593,10 +594,29 @@ 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) call post_data(CS%id_Ue, sqrt(max(0.,2.0*MEKE%MEKE)), CS%diag) - if (CS%id_Ub>0) call post_data(CS%id_Ub, sqrt(max(0.,2.0*MEKE%MEKE*bottomFac2)), CS%diag) - if (CS%id_Ut>0) call post_data(CS%id_Ut, sqrt(max(0.,2.0*MEKE%MEKE*barotrFac2)), CS%diag) + if (CS%id_Ue>0) then + do j=js,je ; do i=is,ie + tmp(i,j) = sqrt(max(0., 2. * MEKE%MEKE(i,j))) + enddo ; enddo + call post_data(CS%id_Ue, tmp, CS%diag) + endif + if (CS%id_Ub>0) then + do j=js,je ; do i=is,ie + tmp(i,j) = sqrt(max(0., 2. * MEKE%MEKE(i,j) * bottomFac2(i,j))) + enddo ; enddo + call post_data(CS%id_Ub, tmp, CS%diag) + endif + if (CS%id_Ut>0) then + do j=js,je ; do i=is,ie + tmp(i,j) = sqrt(max(0., 2. * MEKE%MEKE(i,j) * barotrFac2(i,j))) + enddo ; enddo + call post_data(CS%id_Ut, tmp, CS%diag) + endif if (CS%id_Kh>0) call post_data(CS%id_Kh, MEKE%Kh, CS%diag) if (CS%id_Ku>0) call post_data(CS%id_Ku, MEKE%Ku, CS%diag) if (CS%id_Au>0) call post_data(CS%id_Au, MEKE%Au, CS%diag) From 050aa31b38f78861a330cd871c1cdb3c11e3f689 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 22 Oct 2019 18:14:01 -0600 Subject: [PATCH 7/7] Moves allocation of CS%equilibrium_value inside subroutine MEKE_equilibrium_restoring * Also deletes unneeded variables from subroutine MEKE_equilibrium_restoring. --- src/parameterizations/lateral/MOM_MEKE.F90 | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 8d3e78262c..a009aea1f6 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -334,7 +334,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h endif if (CS%MEKE_equilibrium_restoring) then - call MEKE_equilibrium_restoring(CS, MEKE, G, GV, US, SN_u, SN_v) + call MEKE_equilibrium_restoring(CS, G, US, SN_u, SN_v) do j=js,je ; do i=is,ie src(i,j) = src(i,j) - CS%MEKE_restoring_rate*(MEKE%MEKE(i,j) - CS%equilibrium_value(i,j)) enddo ; enddo @@ -808,28 +808,28 @@ end subroutine MEKE_equilibrium !< This subroutine calculates a new equilibrium value for MEKE at each time step. This is not copied into !! MEKE%MEKE; rather, it is used as a restoring term to nudge MEKE%MEKE back to an equilibrium value -subroutine MEKE_equilibrium_restoring(CS, MEKE, G, GV, US, SN_u, SN_v) +subroutine MEKE_equilibrium_restoring(CS, G, US, SN_u, SN_v) type(ocean_grid_type), intent(inout) :: G !< Ocean grid. - type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure. - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type. type(MEKE_CS), pointer :: CS !< MEKE control structure. - type(MEKE_type), pointer :: MEKE !< A structure with MEKE data. real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1]. real, dimension(SZI_(G),SZJB_(G)), intent(in) :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1]. ! Local variables - real :: SN ! The local Eady growth rate [T-1 ~> s-1] - integer :: i, j, is, ie, js, je, n1, n2 - real :: cd2 + real :: SN ! The local Eady growth rate [T-1 ~> s-1] + integer :: i, j, is, ie, js, je ! local indices + real :: cd2 ! bottom drag is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec cd2 = CS%cdrag**2 + if (.not. associated(CS%equilibrium_value)) allocate(CS%equilibrium_value(SZI_(G),SZJ_(G))) + CS%equilibrium_value(:,:) = 0.0 + !$OMP do do j=js,je ; do i=is,ie ! SN = 0.25*max( (SN_u(I,j) + SN_u(I-1,j)) + (SN_v(i,J) + SN_v(i,J-1)), 0.) ! This avoids extremes values in equilibrium solution due to bad values in SN_u, SN_v SN = min(SN_u(I,j), SN_u(I-1,j), SN_v(i,J), SN_v(i,J-1)) - CS%equilibrium_value(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_m*G%bathyT(i,j))**2 / cd2 enddo ; enddo @@ -837,7 +837,6 @@ subroutine MEKE_equilibrium_restoring(CS, MEKE, G, GV, US, SN_u, SN_v) end subroutine MEKE_equilibrium_restoring - !> Calculates the eddy mixing length scale and \f$\gamma_b\f$ and \f$\gamma_t\f$ !! functions that are ratios of either bottom or barotropic eddy energy to the !! column eddy energy, respectively. See \ref section_MEKE_equations. @@ -1077,7 +1076,6 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS) call get_param(param_file, mdl, "MEKE_RESTORING_TIMESCALE", MEKE_restoring_timescale, & "The timescale used to nudge MEKE toward its equilibrium value.", units="s", & default=1e6, scale=US%T_to_s) - allocate(CS%equilibrium_value(isd:ied,jsd:jed)) ; CS%equilibrium_value(:,:) = 0.0 CS%MEKE_restoring_rate = 1.0 / MEKE_restoring_timescale endif