Some small bugs/questions about second-moment budgets

[BrodiePearson]

I've been going through the budget equations to match them with the Langmuir closure Overleaf document and had three of questions/comments:

• It looks like the first part of the pressure closure in the u'w' and v'w' budgets is multiplied by an extra factor of a half, as the KE is defined with the factor of half already built in (see below). The factor of half outside makes sense to me, but I'm not sure where the second 0.5 or 0.5_RKIND comes from
  

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

  Lines 572 to 573 in 4923b0e

  	uwtend2(k,iCell) = 0.5*((0.8 - 4.0*alpha1/3.0)*0.5*KE**2.0 + &
  	(alpha1 - alpha2)*u2(i1,k,iCell) + (alpha1 + &

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

  Lines 587 to 588 in 4923b0e

  	0.5_RKIND*((0.8_RKIND - 4.0*alpha1/3.0)*0.5_RKIND*KE**2.0_RKIND + &
  	(alpha1 - alpha2)*v2(i1,k,iCell) + (alpha1 + &

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

  Lines 498 to 499 in 4923b0e

  	KE = sqrt(0.5_RKIND*(u2(i1,k,iCell) + v2(i1,k,iCell) + areaFraction(k,iCell)* &
  	(1.0_RKIND - areaFraction(k,iCell))*wumd(k,iCell)**2.0))

• I don't follow the flux dissipation terms in the u'v', v'w', u'w', w't', w's'. (which contain kappa_FL or kappa_VAR; I'm also not sure why the latter is defined separate to kappa_FL). I assume these represent the viscous/diffusive terms that have the form $\partial^2(\overline{u'w'})/\partial z^2$, but the numerator looks like a first-order numerical derivative instead of second-order. In the u'w' budget these terms are:
  

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

  Lines 580 to 581 in 4923b0e

  	kappa_FL*(uw(i1,k-1,iCell) - uw(i1,k+1,iCell)) / &
  	(ze(k-1,iCell) - ze(k+1,iCell))**2.0

  
  But I think the simplest second-order derivative (assuming that vertical grid spacing is constant) would be kappa_FL*(uw(i1,k-1,iCell) - 2.0_RKIND*uw(i1,k,iCell) + uw(i1,k+1,iCell)) / ( 0.5_RKIND*(ze(k-1,iCell) - ze(k+1,iCell))**2.0_RKIND)

• The coefficients for the buoyancy pressure closure change between the different budgets (w'w' contains C2 while all the other budgets contain beta5). Are these connected in your parameter definitions (C2=1-beta5 for energy conservation?), and if so could they be combined for simplicity in the code? It also looks like B could be used to replace the buoyancy flux in the w'w' budget as in the other budgets.
  

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

  Lines 512 to 513 in 4923b0e

  	w2tend4(k,iCell) = (2.0_RKIND - 4.0_RKIND/3.0_RKIND*C_2)*Mc(k,iCell)* &
  	(grav*alphaT(k,iCell)*tumd(k,iCell) - grav*betaS(k,iCell)*sumd(k,iCell))

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

  Line 605 in 4923b0e

  u2tend4(k,iCell) = 2.0_RKIND/3.0_RKIND*(1.0_RKIND-beta5)*B

[BrodiePearson]

@amrapallig @vanroekel @qingli411

[vanroekel]

@BrodiePearson thanks for continuing to dig through. I'll take a shot at your equation

1. you're right on that extra half, that shouldn't be there. Previously I had defined the KE (L498-499) without the 0.5 but when I put it in I missed making that correction

2. Again right on. The second derivative is wrong. We never use this but would be important to fix. It is something CLUBB developers have suggested a need for. Testing by @amrapallig suggests it doesn't run with this, likely because of the incorrect derivative.

3. beta5 and c2 are identical 0.5, only b/c as you say c2 = 1-beta5. It would be good to unify that, it's just poorly mixed notation.

[qingli411]

Here is my branch of the code. I have changed the variable names so that the coefficients should be consistent with the Overleaf file. However, this branch was based on an earlier version of the code so a lot of the recent changes have not been merged in. So many of the issues we saw earlier haven't been fixed in this branch. But I just wanted to put it here for reference. It would be good to unify these branches as we move forward.

[BrodiePearson]

@qingli411 Thanks for the link to your code and for adding the new constant names, I also like the commenting that you have added for budget terms etc. I agree that unifying these branches would be a good idea, especially as it will only get harder the longer we wait. I hope that it is not too challenging to do so now.

@vanroekel Thanks for your responses, I'm glad the comments were sensible. It looks like Qing's code has actually fixed the bugs in points (1) and (3), although it looks like the viscous terms in (2) are removed in @qingli411's code - that's good to know as we merge the codes.