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

[Jean-Michel Campin]

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