Skip to content

Commit 5dc8d6e

Browse files
vanroekelvanroekel
committed
Adds splat parameterization
1 parent a4211df commit 5dc8d6e

File tree

3 files changed

+190
-62
lines changed

3 files changed

+190
-62
lines changed

‎src/core_ocean/adc_mixing/Registry_adc_mixing_fields_opts.xml

Lines changed: 48 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -106,7 +106,39 @@
106106
<nml_option name="config_adc_kappaW3" type="real" default_value="0" units="m^2/s"
107107
description="background diffusion coefficient for w3"
108108
possible_values="small positive real numbers"
109-
/>
109+
/>
110+
<nml_option name="config_adc_bc_wstar" type="real" default_value="0.3" units="unitless"
111+
description="scaling factor in front of the surface boundary condition on wstar"
112+
possible_values="small positive real numbers"
113+
/>
114+
<nml_option name="config_adc_frictionVelocityMin" type="real" default_value="1.0e-5" units="m^2/s^2"
115+
description="minimum allowable surface friction velocity squared"
116+
possible_values="very small positive real numbers"
117+
/>
118+
<nml_option name="config_adc_bc_const" type="real" default_value="1.8" units="unitless"
119+
description="constant multiplying surface boundary condition, taken from CLUBB"
120+
possible_values="small positive real numbers"
121+
/>
122+
<nml_option name="config_adc_bc_const_wp2" type="real" default_value="1.8" units="unitless"
123+
description="identical to previous parameter but for wp2 to allow disabling"
124+
possible_values="small positive real numbers or zero"
125+
/>
126+
<nml_option name="config_adc_use_splat_parameterization" type="logical" default_value=".true." units="unitless"
127+
description="flag to use the splat parameterization"
128+
possible_values=".true. or .false."
129+
/>
130+
<nml_option name="config_adc_splat_tend_max" type="real" default_value="1.0e-5" units="m^2/s^3"
131+
description="maximum tendency of the splat term"
132+
possible_values="very small positive values"
133+
/>
134+
<nml_option name="config_adc_splat_wp2_val" type="real" default_value="2.0" units="unitless"
135+
description="factor multiplying the splat term in wp2 and wp3"
136+
possible_values="small positive real numbers"
137+
/>
138+
<nml_option name="config_adc_up2_vp2_factor" type="real" default_value="4.0" units="unitless"
139+
description="factor multiplying the boundary condition on up2 and vp2"
140+
possible_values="small positive real numbers"
141+
/>
110142
</nml_record>
111143
<packages>
112144
<package name="adcMixingPKG" description="This package includes variables required for the assumed distribution closure model."/>
@@ -241,6 +273,9 @@
241273
<var name="w2tend5" type="real" dimensions="nVertLevelsP1 nCells Time" units="m^2/s^3"
242274
description="shear production term"
243275
/>
276+
<var name="w2tend6" type="real" dimensions="nVertLevelsP1 nCells Time" units="m^2/s^3"
277+
description="conversion of w2 due to splat effects"
278+
/>
244279
<var name="w3tend1" type="real" dimensions="nVertLevels nCells Time" units="m^2/s^3"
245280
description="entrainment production term -- similar to dissipation"
246281
/>
@@ -255,7 +290,10 @@
255290
/>
256291
<var name="w3tend5" type="real" dimensions="nVertLevels nCells Time" units="m^2/s^3"
257292
description="third order moment of buoyancy, turb transport of buoyancy fluxes"
258-
/>
293+
/>
294+
<var name="w3tend6" type="real" dimensions="nVertLevels nCells Time" units="m^2/s^3"
295+
description="destruction due to splat effects"
296+
/>
259297
<var name="wttend1" type="real" dimensions="nVertLevelsP1 nCells Time" units="m^2/s^3"
260298
description="divergence of turbulent transport of flux"
261299
/>
@@ -336,7 +374,10 @@
336374
/>
337375
<var name="u2tend5" type="real" dimensions="nVertLevelsP1 nCells Time" units="m^2/s^3"
338376
description="return to isotropy"
339-
/>
377+
/>
378+
<var name="u2tend6" type="real" dimensions="nVertLevelsP1 nCells Time" units="m^2/s^3"
379+
description="production due to splatting of w'2"
380+
/>
340381
<var name="v2tend1" type="real" dimensions="nVertLevelsP1 nCells Time" units="m^2/s^3"
341382
description="divergence of turbulent transport of flux"
342383
/>
@@ -351,7 +392,10 @@
351392
/>
352393
<var name="v2tend5" type="real" dimensions="nVertLevelsP1 nCells Time" units="m^2/s^3"
353394
description="return to isotropy"
354-
/>
395+
/>
396+
<var name="v2tend6" type="real" dimensions="nVertLevelsP1 nCells Time" units="m^2/s^3"
397+
description="production due to splatting of w'2"
398+
/>
355399
<var name="u2cliptend" type="real" dimensions="nVertLevelsP1 nCells Time" units="m^2/s^3"
356400
description="if less than zero, it is reset to zero, this stores when that is active"
357401
/>

‎src/core_ocean/shared/mpas_ocn_adcReconstruct.F

Lines changed: 93 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -229,15 +229,17 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
229229
real :: KEm1, KEp1, tauUP, tauDN, tomUP, tomDN
230230
real :: tauM1, tau, tauP1, tauAV, utemp, vtemp
231231
real :: B, Cval, diff, wtav, dzmid, Ksps, Sz, Tz, w4k, w4kp1, w2k, w2kp1
232-
real :: lareaFraction, wstar, Q, w3av, tempMoment
232+
real :: lareaFraction, wstar, Q, w3av, tempMoment, frictionVelocity
233233
real :: sfcFrictionVelocitySquared, wtSumUp, wtSumDn, wsSumUp, wsSumDn
234234

235235
real,dimension(nVertLevels,nCells) :: Swumd
236236
real,dimension(nVertLevels,nCells) :: tauw3, tauTemp, tauSalt, tauVel, tauvVel
237237
real,dimension(nVertLevels,nCells) :: areaFractionMid, tumdMid, McMid, wumdMid, sumdMid
238238

239-
real :: Swk
239+
real :: Swk, tau_sfc, d_sqrt_wp2_dz, tp2, sp2, min_wp2_sfc_val, wp2_splat_sfc_correction
240+
real :: min_wps_sfc_val
240241

242+
min_wps_sfc_val = 1.0E-10_RKIND
241243
dt_small = config_adc_timestep
242244
niter = dt / dt_small
243245

@@ -246,6 +248,7 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
246248

247249
stopflag = .false.
248250

251+
do iIter=1,niter
249252
!on further examination build_diagnostics array can live outside the iter loop
250253
do iCell=1,nCells
251254
Q = grav*(alphaT(1,iCell)*wtsfc(iCell) - betaS(1,iCell)*wssfc(iCell))* &
@@ -263,10 +266,13 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
263266
w2t(1,iCell) = -0.3_RKIND*wstar * wtsfc(iCell)
264267
w2s(1,iCell) = 0.3_RKIND*wstar * wssfc(iCell)
265268

266-
sfcFrictionVelocitySquared = sqrt(uwsfc(iCell)**2 + vwsfc(iCell)**2)
269+
sfcFrictionVelocitySquared = uwsfc(iCell) + vwsfc(iCell)
270+
frictionVelocity = sqrt(sfcFrictionVelocitySquared + config_adc_bc_wstar * wstar * wstar)
271+
frictionVelocity = max( config_adc_frictionVelocityMin, frictionVelocity )
267272
do k=1,2
268-
u2(k,1,iCell) = uwsfc(iCell)! sfcFrictionVelocitySquared !+ 0.3*wstar**2.0
269-
v2(k,1,iCell) = vwsfc(iCell)! sfcFrictionVelocitySquared !+ 0.3*wstar**2.0
273+
u2(k,1,iCell) = config_adc_up2_vp2_factor*config_adc_bc_const*frictionVelocity**2.0
274+
v2(k,1,iCell) = config_adc_up2_vp2_factor*config_adc_bc_const*frictionVelocity**2.0
275+
w2(k,1,iCell) = config_adc_bc_const_wp2*frictionVelocity**2.0
270276
uw(k,1,iCell) = -uwsfc(iCell)
271277
vw(k,1,iCell) = -vwsfc(iCell)
272278
wt(k,1,iCell) = wtsfc(iCell)
@@ -276,45 +282,87 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
276282
enddo
277283
enddo
278284

279-
do iIter=1,niter
280-
do iCell=1,nCells-1000000
281-
do k=2,nVertLevels
282-
w3av = 0.5_RKIND*(w2(i1,k-1,iCell) + w2(i1,k,iCell))
283-
284-
Sw = w3(i1,k-1,iCell) / (w3av**1.5_RKIND + 1.0E-15_RKIND)
285-
lareaFraction = 0.5_RKIND + 0.5_RKIND*Sw / sqrt(4.0_RKIND + Sw**2)
285+
! compute the splat effect, adds w3tend and w2tend
286+
if(config_adc_use_splat_parameterization) then
287+
! do k=1 separately for performance reasons, for k=1 we use one sided derivatives
288+
! note that the splat_factor forces the term to be smaller than factor * dt
289+
k=1
290+
do iCell=1,nCells
291+
d_sqrt_wp2_dz = (sqrt(w2(i1,1,iCell)) - sqrt(0.5_RKIND*(w2(i1,1,iCell)+ &
292+
w2(i1,2,iCell)))) / (ze(1,iCell) - zm(1,iCell))
293+
tau_sfc = length(1,iCell) / sqrt(0.5_RKIND*(u2(i1,k,iCell) + v2(i1,k,iCell)))
294+
w2tend6(1,iCell) = min(max(-config_adc_splat_tend_max, -w2(i1,1,iCell)* &
295+
config_adc_splat_wp2_val*tau_sfc*d_sqrt_wp2_dz**2), &
296+
config_adc_splat_tend_max)
297+
tau_sfc = 0.5_RKIND*(length(1,iCell) + length(2,iCell)) / sqrt(0.5_RKIND*( &
298+
0.5_RKIND*(u2(i1,1,iCell) + u2(i1,2,iCell)) + 0.5_RKIND*( &
299+
v2(i1,1,iCell) + v2(i1,2,iCell)) + 0.5_RKIND*(w2(i1,1,iCell) + &
300+
w2(i1,2,iCell))))
301+
d_sqrt_wp2_dz = (sqrt(w2(i1,1,iCell)) - sqrt(w2(i1,2,iCell))) / &
302+
(ze(1,iCell) - ze(2,iCell))
303+
w3tend6(1,iCell) = min(max(-config_adc_splat_tend_max, -w3(i1,1,iCell)* &
304+
config_adc_splat_wp2_val*tau_sfc*d_sqrt_wp2_dz**2), &
305+
config_adc_splat_tend_max)
306+
end do
286307

287-
!This clips w3 to be consistent with the clipped areaFraction
288-
if(lareaFraction < 0.01_RKIND) then
289-
w3(i1,k-1,iCell) = -9.85_RKIND*w3av**1.5
290-
end if
308+
do iCell=1,nCells
309+
do k=2,nVertLevels
310+
d_sqrt_wp2_dz = (sqrt(w2(i1,k-1,iCell)) - sqrt(w2(i1,k+1,iCell))) / &
311+
(ze(k-1,iCell) - ze(k+1,iCell))
312+
tau_sfc = length(k,iCell) / sqrt(0.5_RKIND*(u2(i1,k,iCell) + v2(i1,k,iCell) + &
313+
w2(i1,k,iCell)))
314+
w2tend6(k,iCell) = min(max(-config_adc_splat_tend_max, -w2(i1,k,iCell)* &
315+
config_adc_splat_wp2_val*tau_sfc*d_sqrt_wp2_dz**2), &
316+
config_adc_splat_tend_max)
317+
tau_sfc = 0.5_RKIND*(length(k,iCell) + length(k+1,iCell)) / sqrt(0.5_RKIND*( &
318+
0.5_RKIND*(u2(i1,k,iCell) + u2(i1,k+1,iCell)) + 0.5_RKIND*( &
319+
v2(i1,k,iCell) + v2(i1,k+1,iCell)) + 0.5_RKIND*(w2(i1,k,iCell) + &
320+
w2(i1,k+1,iCell))))
321+
d_sqrt_wp2_dz = (sqrt(w2(i1,k,iCell)) - sqrt(w2(i1,k+1,iCell))) / &
322+
(ze(k,iCell) - ze(k+1,iCell))
323+
w3tend6(1,iCell) = min(max(-config_adc_splat_tend_max, -w3(i1,k,iCell)* &
324+
config_adc_splat_wp2_val*tau_sfc*d_sqrt_wp2_dz**2), &
325+
config_adc_splat_tend_max)
326+
end do
327+
end do
291328

292-
if(lareaFraction > 0.99_RKIND) then
293-
w3(i1,k-1,iCell) = 9.85_RKIND*w3av**1.5
329+
do iCell=1,nCells
330+
tp2 = 0.4_RKIND * config_adc_bc_const * (wt(i1,1,iCell) / frictionVelocity)**2
331+
sp2 = 0.4_RKIND * config_adc_bc_const * (ws(i1,1,iCell) / frictionVelocity)**2
332+
min_wp2_sfc_val = max(1.0E-10_RKIND, wt(i1,1,iCell)**2 / (tp2 * 0.99_RKIND**2 + 1.0E-15_RKIND), &
333+
ws(i1,1,iCell)**2 / (sp2 * 0.99_RKIND**2 + 1.0E-15_RKIND))
334+
tau_sfc = length(1,iCell) / sqrt(0.5_RKIND*(u2(i1,1,iCell) + v2(i1,1,iCell)))
335+
336+
if(w2(k,1,iCell) + tau_sfc * w2tend6(1,iCell) < min_wps_sfc_val) then
337+
wp2_splat_sfc_correction = -w2(i1,1,iCell) + min_wp2_sfc_val
338+
w2(i1,1,iCell) = min_wp2_sfc_val
339+
else
340+
wp2_splat_sfc_correction = tau_sfc * w2tend6(1,iCell)
341+
w2(i1,1,iCell) = w2(i1,1,iCell) + wp2_splat_sfc_correction
294342
end if
295-
296-
Sw = w3(i1,k-1,iCell) / (w3av**1.5 + 1.0E-15_RKIND)
297-
lareaFraction = 0.5_RKIND + 0.5_RKIND*Sw / sqrt(4.0_RKIND + Sw**2)
298-
299-
areaFractionMid(k-1,iCell) = lareaFraction
300-
wumdMid(k-1,iCell) = sqrt(w3av / (areaFractionMid(k-1,iCell) * &
301-
(1.0_RKIND - areaFractionMid(k-1,iCell))))
302-
McMid(k-1,iCell) = areaFractionMid(k-1,iCell)*(1.0_RKIND - &
303-
areaFractionMid(k-1,iCell)) * wumdMid(k-1,iCell)
304-
tumdMid(k-1,iCell) = (0.5*(wt(i1,k-1,iCell) + wt(i1,k,iCell))) / &
305-
(1.0E-12_RKIND + McMid(k-1,iCell))
306-
sumdMid(k-1,iCell) = (0.5*(ws(i1,k-1,iCell) + ws(i1,k,iCell))) / &
307-
(1.0E-12_RKIND + McMid(k-1,iCell))
308-
if(w3av <= epsilon + 1.0e-9) then
309-
areaFractionMid(k-1,iCell) = 0.5_RKIND
310-
! wumdMid(k-1,iCell) = 0.0_RKIND
311-
tumdMid(k-1,iCell) = 0.0_RKIND
312-
sumdMid(k-1,iCell) = 0.0_RKIND
313-
! McMid(k-1,iCell) = 0.0_RKIND
314-
endif
315-
enddo
316-
enddo
317-
343+
u2(i1,1,iCell) = u2(i1,1,iCell) - 0.5_RKIND * wp2_splat_sfc_correction
344+
v2(i1,1,iCell) = v2(i1,1,iCell) - 0.5_RKIND * wp2_splat_sfc_correction
345+
end do
346+
347+
!build up splat tend for u2 and v2
348+
do iCell=1,nCells
349+
do k=2,nVertLevels
350+
u2tend6(k,iCell) = -0.5_RKIND*w2tend6(k,iCell)
351+
v2tend6(k,iCell) = -0.5_RKIND*w2tend6(k,iCell)
352+
end do
353+
end do
354+
else
355+
do iCell=1,nCells
356+
do k=1,nVertLevels
357+
w3tend6(k,iCell) = 0.0_RKIND
358+
w2tend6(k,iCell) = 0.0_RKIND
359+
u2tend6(k,iCell) = 0.0_RKIND
360+
v2tend6(k,iCell) = 0.0_RKIND
361+
end do
362+
end do
363+
end if ! end use splat correction
364+
365+
318366
!Kernel 1 inlined versions of the base arrays, needed for later to make them collapsible
319367
do iCell=1,nCells
320368
do k=2,nVertLevels
@@ -459,7 +507,7 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
459507
betaS(k,iCell)*w2s(k,iCell))
460508

461509
w3tend(i3_f,k,iCell) = w3tend1(k,iCell) + w3tend2(k,iCell) + w3tend3(k,iCell) + &
462-
w3tend4(k,iCell) + w3tend5(k,iCell)
510+
w3tend4(k,iCell) + w3tend5(k,iCell) + w3tend6(k,iCell)
463511

464512
if(k>1 .and. k < nVertLevels .and. kappa_w3 > 0.0) then
465513
w3tend(i3_f,k,iCell) = w3tend(i3_f,k,iCell) + kappa_w3*(w3(i1,k-1,iCell) &
@@ -527,7 +575,7 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
527575
Mc(k,iCell)*(Swumd(k-1,iCell) + Swumd(k,iCell))
528576

529577
w2tend(i3_f,k,iCell) = w2tend1(k,iCell) + w2tend2(k,iCell) + &
530-
w2tend3(k,iCell) + w2tend4(k,iCell) + w2tend5(k,iCell)
578+
w2tend3(k,iCell) + w2tend4(k,iCell) + w2tend5(k,iCell) + w2tend6(k,iCell)
531579

532580
wttend1(k,iCell) = -1.0_RKIND*(Entrainment(k,iCell) + Detrainment(k,iCell)) * &
533581
wumd(k,iCell)*tumd(k,iCell)
@@ -621,7 +669,7 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
621669
1.0_RKIND - areaFraction(k,iCell))*wumd(k,iCell)**2)/3.0_RKIND
622670

623671
u2tend(i3_f,k,iCell) = u2tend1(k,iCell) + u2tend2(k,iCell) + &
624-
u2tend3(k,iCell) + u2tend4(k,iCell) + u2tend5(k,iCell)
672+
u2tend3(k,iCell) + u2tend4(k,iCell) + u2tend5(k,iCell) + u2tend6(k,iCell)
625673

626674
v2tend1(k,iCell) = -(v2w(k-1,iCell) - v2w(k,iCell)) / dzmid
627675
v2tend2(k,iCell) = (1.0_RKIND/3.0_RKIND*alpha1 + alpha2 - &
@@ -635,7 +683,7 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
635683
1.0_RKIND - areaFraction(k,iCell))*wumd(k,iCell)**2)/3.0_RKIND
636684

637685
v2tend(i3_f,k,iCell) = v2tend1(k,iCell) + v2tend2(k,iCell) + &
638-
v2tend3(k,iCell) + v2tend4(k,iCell) + v2tend5(k,iCell)
686+
v2tend3(k,iCell) + v2tend4(k,iCell) + v2tend5(k,iCell) + v2tend6(k,iCell)
639687

640688
uttend(i3_f,k,iCell) = (-(uwt(k-1,iCell) - uwt(k,iCell))/dz - &
641689
uw(i1,k,iCell)*Tz - (1.0_RKIND - alpha3)*wt(i1,k,iCell) &

0 commit comments

Comments
 (0)