Patchwork of changes to ADC implementation

[katsmith133]

This PR makes several changes:

1. Changes the formulation of the transport terms in w3tend to that in Gryanik & Hartmann 2001. This was suggested by Brodie in issue Issues with transport terms in $\overline{w'w'w'}$ budget #7 .

2. Removes extra 0.5 factor from certain KE terms, as was mentioned as point 1 in issue Some small bugs/questions about second-moment budgets #9 .

3. Turns on u2 and v2 boundary conditions and changes their formulations.

Will add figures to demonstrate effect of each change shortly...

[katsmith133]

Results for change 1. Plots of various terms after 4 days of simulation for "Cooling2" case (-100 W heat flux). Red line = results using code prior to this PR, green line = with Gryanik & Hartmann 2001 changes (commit a98b72f), black line = LES results.

Same plots but after only 1 hour of simulation (crashed) for "Wind1" case (tau = 0.1).

[katsmith133]

Results for change 2. For both cases, "Cooling2" and "Wind1", there were no discernable differences made by this change.

[BrodiePearson]

Thanks for putting these figures together Kat, it's great to see that the change in w4 seems to improve the w3 profiles significantly! The plots you added are really interesting, and I have a few questions/points that we could discuss later:

• Just to clarify, are you running the ADC code and LES separately but with identical forcing and initial conditions?
• Did the ADC code crash for the Wind1 case, or was it the LES that crashed?
• I assume there is no wind stress in Cooling2, so u'w' and v'w' are just noise in the LES and there is no mean shear in U or V. This means that the only depth-integrated source of u2 and v2 in this setup are the pressure terms moving w2 into u2 and v2.
• Based on the above point, the u2 and v2 profiles for the Cooling2 case are concerning, particularly at the surface where they are way below the LES values. It seems like they tend to zero at the surface, but perhaps in the surface cell we should take the w2 production produced by buoyancy flux and split it between u2 and v2 (since w2==0 at the surface). It looks like something like that is in the code already but is now commented - I guess this is related to your third point above (while I'm writing this I realized that you are doing each point separately so you may already be working on this).
• The w2 profiles at mid-depth are not great either...my first guess would be that this is related to the low u2 and v2 near the surface. If this surface production were present, the near surface high u2 and v2 could be mixed down into the water column, and the balance at mid-depths would be w2 > u2 = v2, where u2 and v2 are supplied by a mix of conversion from w2 and non-local transport of u2 and v2 from the surface. Without that production, the total TKE of the mixed layer is much lower, and the interior balance becomes w2 = u2 = v2 (lower TKE means a shorter timescale for return-to-isotropy which speeds up this equilibration of velocity components).

[katsmith133]

Results from change 3. For the "Cooling2" case after 4 days. Red line = with the u2 and v2 BC changes (Cmom=Ctherm=4), green line = same previous G&H 2001 lines (code at commit a51d872), black line = LES.


Same plots as before but for "Wind1" case, this time it ran to 4 days (G&H2001 comparison has been omitted since it crashed at 1 hour, as mentioned previously).

Plots with varying Cmom and Ctherm to come...

[katsmith133]

Hi Brody,

1. Yes, the ADC and LES were run separately and to the best of my knowledge (since I did not run these LES), they have the same forcing and initial conditions.
2. ADC code crashed. As you can see from the post I just made, after adding in the u2/v2 boundary conditions, it no longer crashes.
3. Correct, there is no wind stress in this case or mean shear in U or V. For the sake of testing I am just focusing on the two extremes right now, an isolated convective case ("Cooling2") and an isolated shear case ("Wind1"). I will test combinations of them after we get the isolated cases correct. Let me know if you'd like to see test results for a combination.
4. Are you referring to the wstar term in the u2/v2 top boundary condition? I have left this commented out per Luke's suggestion, but I can easily test the effects of adding that term, just let me know if that's what you are referring to.
5. Same response as above. Happy to test the addition of wstar, if that is what you are referring to, as it seems even in the change 3 results this is still an issue (probably because, as you said, u2 and v2 are still low at the surface).

Thanks for putting these figures together Kat, it's great to see that the change in w4 seems to improve the w3 profiles significantly! The plots you added are really interesting, and I have a few questions/points that we could discuss later:

• Just to clarify, are you running the ADC code and LES separately but with identical forcing and initial conditions?
• Did the ADC code crash for the Wind1 case, or was it the LES that crashed?
• I assume there is no wind stress in Cooling2, so u'w' and v'w' are just noise in the LES and there is no mean shear in U or V. This means that the only depth-integrated source of u2 and v2 in this setup are the pressure terms moving w2 into u2 and v2.
• Based on the above point, the u2 and v2 profiles for the Cooling2 case are concerning, particularly at the surface where they are way below the LES values. It seems like they tend to zero at the surface, but perhaps in the surface cell we should take the w2 production produced by buoyancy flux and split it between u2 and v2 (since w2==0 at the surface). It looks like something like that is in the code already but is now commented - I guess this is related to your third point above (while I'm writing this I realized that you are doing each point separately so you may already be working on this).
• The w2 profiles at mid-depth are not great either...my first guess would be that this is related to the low u2 and v2 near the surface. If this surface production were present, the near surface high u2 and v2 could be mixed down into the water column, and the balance at mid-depths would be w2 > u2 = v2, where u2 and v2 are supplied by a mix of conversion from w2 and non-local transport of u2 and v2 from the surface. Without that production, the total TKE of the mixed layer is much lower, and the interior balance becomes w2 = u2 = v2 (lower TKE means a shorter timescale for return-to-isotropy which speeds up this equilibration of velocity components).

[BrodiePearson]

Cool, I agree that these test cases are the most useful ones for identifying issues. I don't see the need for mixed ones yet.

You're right, I'm referring to the wstar term in the u2 and v2 boundary conditions. It would be interesting to see whether u2 and v2 improve in the convective case when this term is included, perhaps with 0.3 (or 0.5) coefficients out front.

Some of the results after you added the uwsfc boundary condition in your above cases makes sense to me (u2 is non-zero at the surface in the Wind simulation), but other changes don't. For example, why do u2 and v2 profiles change in the Cooling simulation where uwsfc is zero so shouldn't change anything. Did you adjust other parameters? I'm not sure what Cmom and Ctherm are, and whether they were also changed.

[katsmith133]

I've got runs going right now that test the inclusion of the wstar term (with 0.5 coefficient out front) for u2 and v2. Should have those results soon.

Yes, oops! I forgot to mention that by turning on C_mom and C_therm in the namelist file (they were =0 before), it turns on the 3rd order moments here:

MPAS-Model/src/core_ocean/shared/mpas_ocn_adcReconstruct.F

Lines 462 to 477 in a96c470

	KE = 0.5_RKIND*(u2(i1,k,iCell) + v2(i1,k,iCell) + w2(i1,k,iCell))
	KEp1 = 0.5_RKIND*(u2(i1,k+1,iCell) + v2(i1,k+1,iCell) + w2(i1,k+1,iCell))
	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
	diff = C_therm*sqrt(0.5*(KE + KEp1)) * lenav
	uwt(k,iCell) = -diff*(ut(i1,k,iCell) - ut(i1,k+1,iCell)) / dz
	vwt(k,iCell) = -diff*(vt(i1,k,iCell) - vt(i1,k+1,iCell)) / dz
	uws(k,iCell) = -diff*(us(i1,k,iCell) - us(i1,k+1,iCell)) / dz
	vws(k,iCell) = -diff*(vs(i1,k,iCell) - vs(i1,k+1,iCell)) / dz

I believe C_mom is equivalent to the S_q in Eq (21) of your ADC w/CL manuscript...

If I turn C_mom=C_therm=0, but leave the u2 and v2 terms non-zero, then the cooling case results go back to the G&H2001 results, no change, and the wind case it again crashes at 1hour and the only differences from the G&H2001 results are shown here (removed LES results as they make it hard to see):

Hopefully that clarifies things, let me know if it does not.

[katsmith133]

Note, @BrodiePearson and @vanroekel, this has been updated: Results with 0.5*wstar^2 added to u2 and v2 boundary condition
Cooling case crashes after 62 hours:

Wind case appears that there is no effects of adding the wstar term.

[vanroekel]

approving based on visual inspection and testing by @katsmith133