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

Lars Czeschel lars.czeschel at zmaw.de
Tue Jun 18 08:45:29 EDT 2013 

Hi Martin and Jean-Michel,

thanks for your efforts. I am relieved there is no bug in the skew flux
formulation. It saves me more than 75000 CPUh.

The differences between both formulations are also strong where steep slopes are close to
boundaries (e.g. Labrador Sea). In these regions the differences (e.g. theta (z=85m) 
are of the same order as between GM_Redi / no GM_Redi. If I understand Jean-Michel's 
reply correctly the differences will be reduced if there are more grid
points between steep slopes and boundaries, i.e. with higher resolution.

Anyway, I think I will stick with the skew-flux formulation as the slopes are the same 
for isopycnic and thickness diffusion (sounds more physical to me). 

The bug in the layers pkg before 1.5 shows me again that I should keep my code updated.

Thanks, Lars

 
On 18.06.2013, at 10:54, Martin Losch wrote:

> 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
> <bottom_moc.png><moc.png><moc.py>_______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support

More information about the MITgcm-support mailing list