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

Jean-Michel Campin jmc at ocean.mit.edu
Mon Jun 17 09:02:45 EDT 2013 

Hi Martin,

I sent a reply yesterday to the old message (Jan 2008) from Christopher Wolfe.

I also reply to the 2nd message from Lars, except that only the 1rst part of the reply
has been kept on the mitgcm.org/pipermail archive, so I put a copy here:
> Date: Fri, 14 Jun 2013 18:07:35 -0400
> From mitgcm-support-bounces at mitgcm.org Fri Jun 14 18: 8:05 2013
> From: Jean-Michel Campin <jmc at ocean.mit.edu>
> To: mitgcm-support at mitgcm.org
> Subject: Re: [MITgcm-support] GM: Skew flux vs Advective form
>
> 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

In short, I don't understand how the diagnose stream-function derived from KWY
can be non zero at the bottom since the Vb velocity is derived from d/dz(KWY)
and should therefore be non divergent. Do you agree with this or did I miss
something ?
My impression is that the divergent aspect (derived afterwards from KWY or 
with pkg/layers, but then it uses the same type of discretisation) is strange,
but cannot excluse a Pb in pkg/layers with skew-flux.

Cheers,
Jean-Michel

On Sat, Jun 15, 2013 at 08:55:46AM +0200, Martin Losch wrote:
> Hi Jean-Michel,
> 
> Lars and Carsten think it's a problem relatled to the skew-flux of gm. (originally they thought it is related to the layers package, but not any more).
> I don't know enough about the gm-implementation to be able to reply to this properly, but it smells like a problem:
> On June 6th Lars sent a message with a temperature difference plot for temperature with a range of 1.5deg, just because of switching between skew flux and advective form in gm, then there are the plots (email of June12) about the divergence of the skew flux-eddy MOC.
> 
> Carsten called me up (on the phone) and said that according to his view of the world, skew flux and advective form are the same mathematically, and even if there are numerical differences in the implementation, they shouldnot be very large, and certainly the eddy-MOC should be non-divergent. I think I can agree with this (but that does not say that Lars' analysis is correct). Carsten wanted to make sure that someone who understands the code has a look (and I don't understand the code, but I said I'd ask you).
> 
> Christopher's old email also points to a problem in the skew flux, but I don't think he reported any results with the advective form of gm. He just compared gm vs. nogm. 
> 
> Again, if I knew enough about the gm-implementation, I would have answered these questions myself. Maybe there is no problem at all, but we have had these occasional support-emails over the past years were people complained about unrealistic tracers (too cold temperatures) and I was always unable to explain these plausibly. As far as I remember these complains were always related to coarse model (probably with gm), but I'd have to look that up.
> I just thought that this is worth looking into, and if we find that there is no problem, all the better.
> 
> Martin
> 
> 
> On Jun 14, 2013, at 8:45 PM, Jean-Michel Campin wrote:
> 
> > Hi Martin,
> > 
> > two things:
> > The messages from Lars Czeschel are relative to the layers diagnostics.
> > Is it a pkg/layers problem (with the skew flux) or not just a diagnostics
> > problem but a more serious pkg/gmredi problem ?
> > And what was the specific point that Carsten mentioned to you ?
> > 
> > And the older message (from Christopher Wolfe, Jan 2008) that Lars pointed 
> > to: http://mitgcm.org/pipermail/mitgcm-support/2008-January/005180.html
> > is about difference between advective form and skew flux form 
> > plus the issue of monotonicity of the scheme.
> > I looked for a reply to this message but did not find one.
> > I can certainly answer something to this old message, if this is
> > useful for anyone.
> > 
> > Cheers,
> > Jean-Michel
> > 
> > On Fri, Jun 14, 2013 at 11:35:26AM +0200, Martin Losch wrote:
> >> Hi guys,
> >> 
> >> I just got a call from Carsten Eden, who wanted to raise awareness for the apparent problem with the skew flux formulation of GM. Here's the thread again: <http://mitgcm.org/pipermail/mitgcm-support/2013-June/008344.html>
> >> 
> >> Do you agree, that this is a serious issue? How can we fix that? May the adjoint stability even be affected in a positive way by this?
> >> 
> >> Martin
> >> 
> >> Begin forwarded message:
> >> 
> >>> From: Lars Czeschel <lars.czeschel at zmaw.de>
> >>> Subject: Re: [MITgcm-support] GM: Skew flux vs Advective form
> >>> Date: June 12, 2013 12:50:08 PM GMT+02:00
> >>> To: <mitgcm-support at mitgcm.org>
> >>> Reply-To: <mitgcm-support at mitgcm.org>
> >>> 
> >>> 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-devel mailing list
> >> MITgcm-devel at mitgcm.org
> >> http://mitgcm.org/mailman/listinfo/mitgcm-devel
> > 
> > 
> > _______________________________________________
> > MITgcm-devel mailing list
> > MITgcm-devel at mitgcm.org
> > http://mitgcm.org/mailman/listinfo/mitgcm-devel
> 
> 
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel

More information about the MITgcm-devel mailing list