[MITgcm-support] GM: Skew flux vs Advective form

Tue Jun 18 04:54:25 EDT 2013

Hi Lars,

I am attaching two plots. I am also including a (python-)script (forgive me, I am trying to move away from matlab, but you'll understand the syntax, it's very similar to matlab), that created these plots, based on the last month-average of a ten year run of tutorial_global_oce_latlon. Basically the script does the same as your matlab code, but following Jean-Michel's instructions loosly, I compute the eddy-MOC from KWY on C-points with hFacCs, rather than V-(or S-)points with hFacS. This way I avoid extra averaging, but since we intergrate around the globe, it should not matter. Bottom line: with this integration, the divergence in the bottom layer is gone (see moc_bottom.png), for comparison I also include the wrong method, just to reproduce your results.

In moc.png you'll see that the eddy-MOCs are still quite different in the southern ocean. So I think we can say that the skew flux form is *NOT* divergent (that's good), but the advective and skew flux form *DO* give different results. I don't find that too surprising and it is explained in Jean-Michel's reply to the old question by Christopher: <http://mitgcm.org/pipermail/mitgcm-support/2013-June/008381.html>
I am not sure, what to suggest. Probably using the advective form with some sensible advection scheme and the GAD_SMOLARKIEWICZ_HACK to avoid overshoots in diffusion.

Martin

On Jun 17, 2013, at 5:59 PM, Jean-Michel Campin <jmc at ocean.mit.edu> wrote:

> Hi Lars,
> 
> Thanks to Martin, I think I see the problem illustrated in Skew_MOC_bottom.jpg:
> Kwy is defined at the same location as wVel (horizontally centered, grid-cell vertical
> interface), and would need to be 
> a) horizontally averaged above the V point
> b) and then masked.
> This is what pkg/layer does (pkg/layers/layers_fluxcalc.F, line 278 to 286).
> 
> The problem in this matlab script is that kwy(:,:,k) is not averaged (not the issue
> regarding non-zero stream-function at the bottom) but also not masked with 
> maskS(:,:,k), resulting in non-zero value at the edge of a bottom stair-case
> (where maskC(j) = 1 but maskC(j-1) = 0). After, the masking (using hFacS in
> Xint=sum(Vb.*DX.*hFacS); ) does not correct for the wrong bottom velocity and 
> it propagates in the stream-function down to k=Nr+1.
> 
> The other part of the matlab script dealing with advective-form case 
> has no problem since GM_PsiYtave is already at the right location (above vVel)
> and masked correctly.
> 
> So, my previous answer was wrong (sorry), but it would be great if you could fix you 
> matlab script and re-do this Skew_MOC_bottom.jpg figure to confirm. 
> 
> Cheers,
> Jean-Michel
> 
> On Fri, Jun 14, 2013 at 06:07:35PM -0400, Jean-Michel Campin wrote:
>> Hi Lars,
>> 
>> I have a problem with this plot: Skew_MOC_bottom.jpg
>> http://mitgcm.org/pipermail/mitgcm-support/attachments/20130612/322705a6/attachment-0008.jpg
>>> From the matlab code you provide:
>>> Vb(i,j,k)=KWY(i,j,k+1)*0.5-KWY(i,j,k)*0.5;
>> and since KWY(:,:,1)=0 (from MITgcm code) and KWY(:,:,16)=0 (as set earlier)
>> we should always get zero for the stream-function at the bottom,
>> independently of what is in the file GM_Kwy-T, as long as KWY(:,:,1)=0
>> 
>> Could you double check this part ?
>> 
>> Thanks,
>> Jean-Michel
>> 
>> On Wed, Jun 12, 2013 at 12:50:08PM +0200, Lars Czeschel wrote:
>>> Hi Ryan,
>>> it seems it is not a problem of the layer diagnostic.
>>> I calculated the eddy-MOC from the bolus velocities in z-levels 
>>> for the 'tutorial_global_oce_latlon' for year 10 with similar results.
>>> 
>>> I attached some figures showing the global eddy MOC for both cases and the differences.
>>> The last two figures show the eddy MOC in the deepest level to give an impression of the 
>>> divergence when using the skew flux formulation.
>>> 
>>> I also attached the matlab code in case I am missing something in my calculation.  
>>> 
>>> Thanks, Lars
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> %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]');
>>> 
>>> ----
>>> 
>>> %Eddy MOC from advective formulation
>>> clear all
>>> DXG = rdmds('DXG');
>>> hFacS = rdmds('hFacS');
>>> YG = rdmds('YG');
>>> DRF = rdmds('DRF');
>>> DRF = squeeze(DRF);
>>> PSIY = rdmds('GM_PsiYtave.0000172800');
>>> 
>>> % compute mer. bolus velocity
>>> PSIY(:,:,16)=0.0;
>>> for i=1:90
>>> for j=1:40
>>> for k=1:15
>>> Vb(i,j,k)=PSIY(i,j,k+1)-PSIY(i,j,k);
>>> 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));
>>> 
>>> % compute eddy 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('Advective form MOC [Sv]');
>>> 
>>> 
>>> 
>>> 
>>> On 06.06.2013, at 18:04, Ryan Abernathey wrote:
>>> 
>>>> Lars,
>>>> Unfortunately, I can't really comment on climatological differences between the skew flux vs adv form.
>>>> But I am curious whether using adv form fixed the problems you were having with the layers diagnostics.
>>>> -Ryan
>>>> 
>>> 
>>>> 
>>>> On Thu, Jun 6, 2013 at 2:24 AM, Lars Czeschel <lars.czeschel at zmaw.de> wrote:
>>>> Hi all,
>>>> 
>>>> i am running a 1deg global model with version MITgcm_63q. Using the layers pkg I noticed
>>>> that my global overturning in density space is divergent (up to 5 Sv in ACC latitudes).
>>>> Without bolus velocities or using  ' GM_AdvForm    = .TRUE.,' in data.gmredi
>>>> the overturning becomes non-divergent.
>>>> 
>>>> As we used the skew flux version frequently in the past (which is also the default) I hoped to find
>>>> a bug in the layers pkg but without success. As a further test I ran the tutorial example
>>>> 'tutorial_global_oce_latlon'  for 10 years with   GM_AdvForm    = .False. compared
>>>> to GM_AdvForm    = .TRUE. in data.gmredi. I attached the quite strong differences in theta (z=85m)
>>>> below.
>>>> 
>>>> To my understanding the differences should be very small even if the small slope assumption
>>>> is made using the skew flux version.
>>>> 
>>>> Did I miss something here or do we have an undiscovered bug in the skew flux code and
>>>> the advective form should be recommended ?
>>>> 
>>>>  Thanks , Lars
>>>> 
>>>> p.s: Several years ago Christopher Wolfe reported a related problem
>>>> http://mitgcm.org/pipermail/mitgcm-support/2008-January/005180.html
>>>> 
>>>> 
>>>> ----------------------------------------------
>>>> Lars Czeschel
>>>> Theoretical Oceanography
>>>> Institut für Meereskunde
>>>> KlimaCampus Universität Hamburg
>>>> Bundesstr. 53
>>>> 20146 Hamburg
>>>> 
>>>> 
>>>> _______________________________________________
>>>> MITgcm-support mailing list
>>>> MITgcm-support at mitgcm.org
>>>> http://mitgcm.org/mailman/listinfo/mitgcm-support
>>>> 
>>>> 
>>>> _______________________________________________
>>>> MITgcm-support mailing list
>>>> MITgcm-support at mitgcm.org
>>>> http://mitgcm.org/mailman/listinfo/mitgcm-support
>>> 
>> 
>>> _______________________________________________
>>> MITgcm-support mailing list
>>> MITgcm-support at mitgcm.org
>>> http://mitgcm.org/mailman/listinfo/mitgcm-support
>> 
>> 
>> _______________________________________________
>> MITgcm-support mailing list
>> MITgcm-support at mitgcm.org
>> http://mitgcm.org/mailman/listinfo/mitgcm-support
> 
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support
-------------- next part --------------
A non-text attachment was scrubbed...
Name: bottom_moc.png
Type: image/png
Size: 60874 bytes
Desc: not available
URL: <http://mitgcm.org/pipermail/mitgcm-support/attachments/20130618/d0a0c63a/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: moc.png
Type: image/png
Size: 65488 bytes
Desc: not available
URL: <http://mitgcm.org/pipermail/mitgcm-support/attachments/20130618/d0a0c63a/attachment-0003.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: moc.py
Type: text/x-python-script
Size: 2269 bytes
Desc: not available
URL: <http://mitgcm.org/pipermail/mitgcm-support/attachments/20130618/d0a0c63a/attachment-0001.bin>