[MITgcm-devel] bolus velocites

[Jeffery R Scott]

Thanks Jean-Michel

I looked at your suggested approach and definitely seems preferable to the old/alternate method, particularly
when the ultimate objective is to add to VVEL and compute a residual MOC.

I compared the two for the 2.8deg global experiment and results were pretty similar;
your way produced a bit stronger cancellation of the Deacon Cell in the upper layers.

A note on your matlab code: one needs to zero out the NaNs in V_bolus field (to put in a cumsum) or
rewrite w/o for loop using something like the n=find(hFacS > 0) Vb(n)= trick in the old approach code.

For now we’ll keep this as posted to devel and allow others a chance to chime in;
hopefully at some point I can update the pkg/gmredi doc and provide a link to this algorithm.
As I mentioned, in the new SO tutorial I skirt the issue by strongly recommending using the advective form :)
(which I think is the best approach for the global setups makeover too)

Jeff


> On Sep 15, 2020, at 6:43 PM, Jean-Michel Campin <jmc at mit.edu> wrote:
> 
> 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]');
>> 
>>