Skip to content

adds reduced formulation (down gradient) option for horiz momentum flux terms #25

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 3 commits into
base: ocean/addADCMixing
Choose a base branch
from
Open

adds reduced formulation (down gradient) option for horiz momentum flux terms #25

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
ef8692c
adds reduced formulation (down gradiet) option for horiz momentum flu…
katsmith133 Oct 25, 2021
50045e8
fixes sign of horiz momentum flux terms
katsmith133 Oct 25, 2021
a149486
reduces unnecessary computations if using the horix down gradient mom…
katsmith133 Oct 26, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 10 additions & 2 deletions src/core_ocean/adc_mixing/Registry_adc_mixing_fields_opts.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,10 @@
description="flag to use buoyancy length scale from previous closure, temporary"
possible_values=".true. or .false."
/>
<nml_option name="config_use_grad_diff_horiz_mom_flux" type="logical" default_value=".false." units="unitless"
description="flag to use down gradient diffusion closure for horizontal momentum fluxes"
possible_values=".true. or .false."
/>
<nml_option name="config_adc_truncate_tend" type="logical" default_value=".false." units="unitless"
description="flag to truncate tendency terms"
possible_values=".true. or .false."
Expand Down Expand Up @@ -66,15 +70,19 @@
<nml_option name="config_adc_alpha_tracer1" type="real" default_value="0.6" units="unitless"
description="alpha_scalar1 first rapid pressure coefficient in tracer flux budgets (Mironov 2001)"
possible_values="suggested values: 0.6 (Canuto 2001) or equal to alpha_scalar2, >0, <1 (Harcourt 2013)"
/>
/>
<nml_option name="config_adc_alpha_tracer2" type="real" default_value="1.0" units="unitless"
description="alpha_scalar2 rapid pressure coefficient in tracer flux budgets (Mironov 2001)"
possible_values="suggested values: 1.0 (Canuto 2001) or equal to alpha_scalar2, >0, <1 (Harcourt 2013)"
/>
/>
<nml_option name="config_adc_c11" type="real" default_value="0.1" units="unitless"
description="c_11 is the buoyancy contribution in the pressure correlation in w3 eqn"
possible_values="values near 0.1"
/>
<nml_option name="config_adc_Cmom_flux" type="real" default_value="0.1" units="unitless"
description="horizontal momentum flux down gradient coefficient"
possible_values="values near 0.1"
/>
<nml_option name="config_adc_Cmom" type="real" default_value="0.0" units="unitless"
description="third order down gradient coefficient, from lappen and randall 2001b"
possible_values="values near 0.5"
Expand Down
128 changes: 69 additions & 59 deletions src/core_ocean/shared/mpas_ocn_adcReconstruct.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -488,11 +488,13 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
lenav = 0.5_RKIND*(length(k,iCell) + length(k+1,iCell))
diff = C_mom * sqrt(0.5_RKIND*(KE + KEp1)) * lenav
dz = ze(k,iCell) - ze(k+1,iCell)
uw2(k,iCell) = -diff*(uw(i1,k,iCell) - uw(i1,k+1,iCell)) / dz
vw2(k,iCell) = -diff*(vw(i1,k,iCell) - vw(i1,k+1,iCell)) / dz
u2w(k,iCell) = -diff*(u2(i1,k,iCell) - u2(i1,k+1,iCell)) / dz
v2w(k,iCell) = -diff*(v2(i1,k,iCell) - v2(i1,k+1,iCell)) / dz
uvw(k,iCell) = -diff*(uv(i1,k,iCell) - uv(i1,k+1,iCell)) / dz
if (.not. config_use_grad_diff_horiz_mom_flux) then
uw2(k,iCell) = -diff*(uw(i1,k,iCell) - uw(i1,k+1,iCell)) / dz
vw2(k,iCell) = -diff*(vw(i1,k,iCell) - vw(i1,k+1,iCell)) / dz
uvw(k,iCell) = -diff*(uv(i1,k,iCell) - uv(i1,k+1,iCell)) / dz
endif

diff = C_therm*sqrt(0.5_RKIND*(KE + KEp1)) * lenav
uwt(k,iCell) = -diff*(ut(i1,k,iCell) - ut(i1,k+1,iCell)) / dz
Expand Down Expand Up @@ -675,6 +677,7 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
endif


if (.not. config_use_grad_diff_horiz_mom_flux) then
! CALCULATE UW TENDENCY TERMS
! The uw tendency is the sum of the following terms (which are saved individually)
! (1) turbulent transport of uw [ -d(uww)/dz ]
Expand All @@ -683,30 +686,30 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
! (4) buoyancy production of uw [ u*buoyancy]
! (5) dissipation of uw [ -kappa_FL*[d(uw)/dz]^2 ]
! The tendency term is truncated if config_adc_truncate_tend flag it switched on
uwtend1(k,iCell) = -(uw2(k-1,iCell) - uw2(k,iCell)) / dzmid
uwtend1(k,iCell) = -(uw2(k-1,iCell) - uw2(k,iCell)) / dzmid

uwtend2(k,iCell) = -tauVel(k,iCell)*uw(i1,k,iCell) + &
0.5_RKIND*((alpha_0 - 4.0_RKIND*alpha_1/3.0_RKIND)*KE**2.0_RKIND + &
(alpha_1 - alpha_2)*u2(i1,k,iCell) + (alpha_1 + &
alpha_2)*w2(i1,k,iCell))*Uz + 0.5_RKIND*(alpha_1 - alpha_2)* &
uv(i1,k,iCell)*Vz - C_b*grav*(alphaT(k,iCell)* &
ut(i1,k,iCell) - betaS(k,iCell)*us(i1,k,iCell))
uwtend2(k,iCell) = -tauVel(k,iCell)*uw(i1,k,iCell) &
+ 0.5_RKIND*((alpha_0 - 4.0_RKIND*alpha_1/3.0_RKIND)*KE**2.0_RKIND + &
(alpha_1 - alpha_2)*u2(i1,k,iCell) + (alpha_1 + &
alpha_2)*w2(i1,k,iCell))*Uz + 0.5_RKIND*(alpha_1 - alpha_2)* &
uv(i1,k,iCell)*Vz - C_b*grav*(alphaT(k,iCell)* &
ut(i1,k,iCell) - betaS(k,iCell)*us(i1,k,iCell))

uwtend3(k,iCell) = - w2(i1,k,iCell)*Uz
uwtend3(k,iCell) = - w2(i1,k,iCell)*Uz

uwtend4(k,iCell) = grav*(alphaT(k,iCell)* &
ut(i1,k,iCell) - betaS(k,iCell)*us(i1,k,iCell))
uwtend4(k,iCell) = grav*(alphaT(k,iCell)* &
ut(i1,k,iCell) - betaS(k,iCell)*us(i1,k,iCell))

uwtend5(k,iCell) = kappa_FL*(uw(i1,k-1,iCell) - uw(i1,k+1,iCell)) / &
(ze(k-1,iCell) - ze(k+1,iCell))**2.0_RKIND
uwtend5(k,iCell) = kappa_FL*(uw(i1,k-1,iCell) - uw(i1,k+1,iCell)) / &
(ze(k-1,iCell) - ze(k+1,iCell))**2.0_RKIND
! Sum the above tendency terms to find the total uw tendency.
uwtend(i3_f,k,iCell) = uwtend1(k,iCell) + uwtend2(k,iCell) + &
uwtend3(k,iCell) + uwtend4(k,iCell) + uwtend5(k,iCell)
uwtend(i3_f,k,iCell) = uwtend1(k,iCell) + uwtend2(k,iCell) + &
uwtend3(k,iCell) + uwtend4(k,iCell) + uwtend5(k,iCell)

if (config_adc_truncate_tend) then
uwtend(i3_f,k,iCell) = FLOAT (INT(uwtend(i3_f,k,iCell) * adcRound &
+ 0.5_RKIND)) / adcRound
endif
if (config_adc_truncate_tend) then
uwtend(i3_f,k,iCell) = FLOAT (INT(uwtend(i3_f,k,iCell) * adcRound &
+ 0.5_RKIND)) / adcRound
endif


! CALCULATE VW TENDENCY TERMS
Expand All @@ -717,30 +720,30 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
! (4) buoyancy production of vw [ v*buoyancy]
! (5) dissipation of vw [ -kappa_FL*[d(vw)/dz]^2 ]
! The tendency term is truncated if config_adc_truncate_tend flag it switched on
vwtend1(k,iCell) = -(vw2(k-1,iCell) - vw2(k,iCell)) / dzmid
vwtend1(k,iCell) = -(vw2(k-1,iCell) - vw2(k,iCell)) / dzmid

vwtend2(k,iCell) = -tauVel(k,iCell)*vw(i1,k,iCell) + &
0.5_RKIND*((alpha_0 - 4.0_RKIND*alpha_1/3.0_RKIND)*KE**2.0_RKIND + &
(alpha_1 - alpha_2)*v2(i1,k,iCell) + (alpha_1 + &
alpha_2)*w2(i1,k,iCell))*Vz + 0.5_RKIND*(alpha_1 - alpha_2)* &
uv(i1,k,iCell)*Uz - C_b*grav*(alphaT(k,iCell)* &
vt(i1,k,iCell) - betaS(k,iCell)*vs(i1,k,iCell))
vwtend2(k,iCell) = -tauVel(k,iCell)*vw(i1,k,iCell) &
+ 0.5_RKIND*((alpha_0 - 4.0_RKIND*alpha_1/3.0_RKIND)*KE**2.0_RKIND + &
(alpha_1 - alpha_2)*v2(i1,k,iCell) + (alpha_1 + &
alpha_2)*w2(i1,k,iCell))*Vz + 0.5_RKIND*(alpha_1 - alpha_2)* &
uv(i1,k,iCell)*Uz - C_b*grav*(alphaT(k,iCell)* &
vt(i1,k,iCell) - betaS(k,iCell)*vs(i1,k,iCell))

vwtend3(k,iCell) = - w2(i1,k,iCell)*Vz
vwtend3(k,iCell) = - w2(i1,k,iCell)*Vz

vwtend4(k,iCell) = grav*(alphaT(k,iCell)* &
vt(i1,k,iCell) - betaS(k,iCell)*vs(i1,k,iCell))
vwtend4(k,iCell) = grav*(alphaT(k,iCell)* &
vt(i1,k,iCell) - betaS(k,iCell)*vs(i1,k,iCell))

vwtend5(k,iCell) = kappa_FL*(vw(i1,k-1,iCell) - vw(i1,k+1,iCell)) / &
(ze(k-1,iCell) - ze(k+1,iCell))**2.0_RKIND
vwtend5(k,iCell) = kappa_FL*(vw(i1,k-1,iCell) - vw(i1,k+1,iCell)) / &
(ze(k-1,iCell) - ze(k+1,iCell))**2.0_RKIND
! Sum the above tendency terms to find the total vw tendency.
vwtend(i3_f,k,iCell) = vwtend1(k,iCell) + vwtend2(k,iCell) + &
vwtend3(k,iCell) + vwtend4(k,iCell) + vwtend5(k,iCell)
vwtend(i3_f,k,iCell) = vwtend1(k,iCell) + vwtend2(k,iCell) + &
vwtend3(k,iCell) + vwtend4(k,iCell) + vwtend5(k,iCell)

if (config_adc_truncate_tend) then
vwtend(i3_f,k,iCell) = FLOAT (INT(vwtend(i3_f,k,iCell) * adcRound &
+ 0.5_RKIND)) / adcRound
endif
if (config_adc_truncate_tend) then
vwtend(i3_f,k,iCell) = FLOAT (INT(vwtend(i3_f,k,iCell) * adcRound &
+ 0.5_RKIND)) / adcRound
endif


! CALCULATE UV TENDENCY TERMS (these are not saved individually)
Expand All @@ -750,20 +753,20 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
! (3) shear production of uv [uw*dV/dz + vw*dU/dz]
! (4) dissipation of uv [ -kappa_FL*[d(uv)/dz]^2 ]
! The tendency term is truncated if config_adc_truncate_tend flag it switched on
uvtend(i3_f,k,iCell) = -(uvw(k-1,iCell) - uvw(k,iCell)) / dz & ! transport
- tauVel(k,iCell)*uv(i1,k,iCell) & ! Rotta
+ 0.5_RKIND*(alpha_1+alpha_2)* &
(uw(i1,k,iCell)*Vz + vw(i1,k,iCell)*Uz) & ! rapid
- (uw(i1,k,iCell)*Vz + vw(i1,k,iCell)*Uz) & ! shear production
+ kappa_VAR*(uv(i1,k-1,iCell) - uv(i1,k+1,iCell)) / &
(ze(k-1,iCell) - ze(k+1,iCell))**2.0_RKIND ! dissipation
uvtend(i3_f,k,iCell) = -(uvw(k-1,iCell) - uvw(k,iCell)) / dz & ! transport
- tauVel(k,iCell)*uv(i1,k,iCell) & ! Rotta
+ 0.5_RKIND*(alpha_1+alpha_2)* &
(uw(i1,k,iCell)*Vz + vw(i1,k,iCell)*Uz) & ! rapid
- (uw(i1,k,iCell)*Vz + vw(i1,k,iCell)*Uz) & ! shear production
+ kappa_VAR*(uv(i1,k-1,iCell) - uv(i1,k+1,iCell)) / &
(ze(k-1,iCell) - ze(k+1,iCell))**2.0_RKIND ! dissipation

! Truncate the uv tendency term if the truncation flag is switched on
if (config_adc_truncate_tend) then
uvtend(i3_f,k,iCell) = FLOAT (INT(uvtend(i3_f,k,iCell) * adcRound &
+ 0.5_RKIND)) / adcRound
endif

if (config_adc_truncate_tend) then
uvtend(i3_f,k,iCell) = FLOAT (INT(uvtend(i3_f,k,iCell) * adcRound &
+ 0.5_RKIND)) / adcRound
endif
endif

! CALCULATE U2 TENDENCY TERMS
! The u2 tendency is the sum of the following terms (which are saved individually)
Expand Down Expand Up @@ -1064,19 +1067,26 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
v2(i2,k,iCell) = epsilon
endif

uw(i2,k,iCell) = uw(i1,k,iCell) + dt_small*(Cw3 * uwtend(i3_f,k,iCell) + &
Cw2 * uwtend(i2_f,k,iCell) + Cw1 * uwtend(i1_f,k,iCell))
vw(i2,k,iCell) = vw(i1,k,iCell) + dt_small*(Cw3 * vwtend(i3_f,k,iCell) + &
Cw2 * vwtend(i2_f,k,iCell) + Cw1 * vwtend(i1_f,k,iCell))
uv(i2,k,iCell) = uv(i1,k,iCell) + dt_small*(Cw3 * uvtend(i3_f,k,iCell) + &
Cw2 * uvtend(i2_f,k,iCell) + Cw1 * uvtend(i1_f,k,iCell))
if (config_use_grad_diff_horiz_mom_flux) then
Uz = (uvel(k-1,iCell) - uvel(k,iCell)) / dzmid
Vz = (vvel(k-1,iCell) - vvel(k,iCell)) / dzmid
uw(i2,k,iCell) = -C_mom_flux*sqrt(0.5_RKIND*(u2(i1,k,iCell) + v2(i1,k,iCell) + w2(i1,k,iCell)))*length(k,iCell)*Uz
vw(i2,k,iCell) = -C_mom_flux*sqrt(0.5_RKIND*(u2(i1,k,iCell) + v2(i1,k,iCell) + w2(i1,k,iCell)))*length(k,iCell)*Vz
else
uw(i2,k,iCell) = uw(i1,k,iCell) + dt_small*(Cw3 * uwtend(i3_f,k,iCell) + &
Cw2 * uwtend(i2_f,k,iCell) + Cw1 * uwtend(i1_f,k,iCell))
vw(i2,k,iCell) = vw(i1,k,iCell) + dt_small*(Cw3 * vwtend(i3_f,k,iCell) + &
Cw2 * vwtend(i2_f,k,iCell) + Cw1 * vwtend(i1_f,k,iCell))
uv(i2,k,iCell) = uv(i1,k,iCell) + dt_small*(Cw3 * uvtend(i3_f,k,iCell) + &
Cw2 * uvtend(i2_f,k,iCell) + Cw1 * uvtend(i1_f,k,iCell))
end if
ut(i2,k,iCell) = ut(i1,k,iCell) + dt_small*(Cw3 * uttend(i3_f,k,iCell) + &
Cw2 * uttend(i2_f,k,iCell) + Cw1 * uttend(i1_f,k,iCell))
vt(i2,k,iCell) = vt(i1,k,iCell) + dt_small*(Cw3 * vttend(i3_f,k,iCell) + &
Cw2 * vttend(i2_f,k,iCell) + Cw1 * vttend(i1_f,k,iCell))
wt(i2,k,iCell) = (wt(i1,k,iCell) + dt_small*(Cw3 * wttend(i3_f,k,iCell) + &
Cw2 * wttend(i2_f,k,iCell) + Cw1 * wttend(i1_f,k,iCell))) / &
(1.0_RKIND + dt_small*tau_tracer(k,iCell))
vt(i2,k,iCell) = vt(i1,k,iCell) + dt_small*(Cw3 * vttend(i3_f,k,iCell) + &
Cw2 * vttend(i2_f,k,iCell) + Cw1 * vttend(i1_f,k,iCell))
us(i2,k,iCell) = us(i1,k,iCell) + dt_small*(Cw3 * ustend(i3_f,k,iCell) + &
Cw2 * ustend(i2_f,k,iCell) + Cw1 * ustend(i1_f,k,iCell))
vs(i2,k,iCell) = vs(i1,k,iCell) + dt_small*(Cw3 * vstend(i3_f,k,iCell) + &
Expand Down
7 changes: 6 additions & 1 deletion src/core_ocean/shared/mpas_ocn_turbulence.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ module ocn_turbulence

real(kind=RKIND),public :: epsilon, &! small value below which w2, v2, u2 cannot go
sigmat, Ko, c_b, c_b_tracer, alpha_tracer1, alpha_tracer2, c11, &
alpha_0, alpha_1, alpha_2, B1, Kt, grav, c_mom, &
alpha_0, alpha_1, alpha_2, B1, Kt, grav, c_mom, c_mom_flux, &
c_therm, c_mom_w3, c_epsilon, c_slow, c_slow_tracer, slow_w_factor, &
kappa_FL, kappa_w3, kappa_VAR, Cww_D, Cww_E, adcRound

Expand Down Expand Up @@ -162,6 +162,7 @@ subroutine ocn_turbulenceCreate(domain)
Kt = 0.4_RKIND
c_b = config_adc_c_b
c11 = config_adc_c11
c_mom_flux = config_adc_Cmom_flux
c_mom = config_adc_Cmom
c_therm = config_adc_Ctherm
c_mom_w3 = config_adc_Cmom_w3
Expand Down Expand Up @@ -560,6 +561,7 @@ subroutine ocn_turbulenceCreate(domain)
!$acc Kt, &
!$acc c_1, &
!$acc c_2, &
!$acc c_mom_flux, &
!$acc c_mom, &
!$acc c_therm, &
!$acc c_mom_w3, &
Expand Down Expand Up @@ -766,6 +768,7 @@ subroutine ocn_turbulenceUpdateFields(domain)
!$acc Kt, &
!$acc c_1, &
!$acc c_2, &
!$acc c_mom_flux, &
!$acc c_mom, &
!$acc c_therm, &
!$acc c_mom_w3, &
Expand Down Expand Up @@ -1186,6 +1189,7 @@ subroutine ocn_turbulenceDestroy(ierr)
alpha_tracer2 = 0.0_RKIND
B1 = 0.0_RKIND
Kt = 0.0_RKIND
c_mom_flux = 0.0_RKIND
c_mom = 0.0_RKIND
c_therm = 0.0_RKIND
c_mom_w3 = 0.0_RKIND
Expand Down Expand Up @@ -1220,6 +1224,7 @@ subroutine ocn_turbulenceDestroy(ierr)
!$acc Kt, &
!$acc c_1, &
!$acc c_2, &
!$acc c_mom_flux, &
!$acc c_mom, &
!$acc c_therm, &
!$acc c_mom_w3, &
Expand Down