[MITgcm-devel] bolus velocites

[Jean-Michel Campin]

Hi Jeff,

I took a quick look at your matlab script, sign and factor of 1/5 seems good to me.

However, I would prefer to do things differently:
 a) comute 3-D strean-function
 b) from 3-D stream-function, compute bolus velosity.
The reason is that we know what Boundary Condition to put on stream-function,
(and (b) is traitforward), and this way we are sure that it's divergence free.
Whereas the way you propose, i.e., computing a V_bollus at grid-cell center and 
then putting it back at vVel location, is more tricky regarding BC.
Also a little tricky on non-cartesian grid (for the same reasons).

So, I would suggest:
 for k=1:Nr
 psiY(i,j,k) = 0.25 * ( Kwy(i,j,k)*rAc(i,j) +  Kwy(i,j-1,k)*rAc(i,j-1) )/rAw(i,j) 
                    * maskS(i,j,k)    <-- this applies BC
 with PsiY(:,:,Nr+1) = 0 ; 
then:
 V_bolus(i,j,k) = ( psiY(i,j,k+1) - psiY(i,j,k) )/drF(k)/hFacS(i,j,k)

Cheers,
Jean-Michel

PS: I am cc MITgcm-devel so that we keep a record (and also if I am wrong other
  will complain).
On Tue, Sep 15, 2020 at 08:55:20PM +0000, Jeffery R Scott wrote:
> 
> Hi Jean-Michel,
> 
> This is relevant for the global tutorial, but in the short term I???m using cfc_example and I need to compute bolus velocity and add to Eul velocity.
> 
> I think below (from support list) is correct calculation, not completely straightforward how to get below but follows from what???s in manual.
> (better would be to give this explicitly in the manual, discuss skew vs. advective forms etc.)
> 
> If I want to add to v_bol to VVEL I need to average v_bol from RF points to RC points, and then also to the velocity points in the horizontal plane.
> This is why we push people to use ADV form, as we discuss in SO tutorial.
> 
> Confirm?
> 
> Thanks,
> Jeff
> 
> 
> 
>  %Eddy MOC from skew flux formulation
> >>> clear all
> >>> % load data
> >>> DXG = rdmds('DXG');
> >>> hFacS = rdmds('hFacS');
> >>> YG = rdmds('YG');
> >>> DRF = rdmds('DRF');
> >>> DRF = squeeze(DRF);
> >>> KWY = rdmds('GM_Kwy-T.0000172800');
> >>>
> >>> % compute mer. bolus velocity
> >>> KWY(:,:,16)=0.0;
> >>> for i=1:90
> >>> for j=1:40
> >>> for k=1:15
> >>> Vb(i,j,k)=KWY(i,j,k+1)*0.5-KWY(i,j,k)*0.5;
> >>> end
> >>> end
> >>> end
> >>> DX=repmat(DXG,[1,1,15]);
> >>> ddz=repmat(DRF,[1, 40, 90]);
> >>> ddz=permute(ddz,[3 2 1]);
> >>> n=find(hFacS > 0);
> >>> Vb(n)=Vb(n)./(ddz(n).*hFacS(n));
> >>>
> >>> % integrate MOC
> >>> Xint=sum(Vb.*DX.*hFacS);
> >>> Xint=squeeze(Xint);
> >>> dz=repmat(DRF,[1,40]);
> >>> dz=permute(dz,[2 1]);
> >>> MOC=cumsum(Xint.*dz./1e6, 2);
> >>>
> >>> %Plotting
> >>> [hcl cl]=contourf(YG(1,:),cumsum(DRF),MOC',[-30:5:-10 -10:2:10 10:5:30 ]);
> >>> set(gca,'YDIR','reverse')
> >>> clabel(hcl,cl)
> >>> caxis([-30,30])
> >>> colorbar
> >>> title('Skew Flux MOC [Sv]');
> 
>