[MITgcm-devel] your 2004 Ocean Model. paper

[Jean-Michel Campin]

Hi Samar,

This is a technical point:
the multi-dimensional advection with "cubed-sphere grid" topology
(useCubedSphereExchange=T) is different from a normal lat-lon or
cartesian grid because there are in fact 3 directions (instead of 2 dir,
x and y) and therefore 3 steps in multi-dim advection. This is 
represented on fig 2.18 in the manual.

Now, let's assume the advection scheme needs "adv" pts on each side 
(e.g., adv=2 with DST-3), if I consider any tile on face 3 that is touching 
the Western edge of the face, and thinking in reverse order (from npass = 3 to 1):
To get a valid updated value in the SW corner of this tile after npass=3 (Y update),
I need valid value in the Southern halo region (i=1:sNx, j=1-adv:0) of 
this tile at the end of npass=2.
This requires (npass=2, X update) to have valid value in the SW corner 
of the halo region of this tile (i=1-adv:0, j=1-adv:0) at the end of npass=1.
And this later requirement can only be achieved (npass=1, Y update in Western halo
at face Western edge) if I start with valid value in the region:
 i=1-adv:0 and j=1-2*adv:0+adv.

Note that this requirement is not directly related to a face "wet-corner"
(I mean, face corner that is in the ocean) which does not happen
with the lat-lon-cap, ecco-4 set-up, but to the fact that
some tile-corners are wet (and this happen, I think, in the lat-lon-cap
set-up: e.g, face 5 is a large one and has the same issue as face-3,
the small one, as pointed above).

In the MITgcm, we try to stop the run if using CubedSphereExchange
and multi-dim advection with Olx,Oly < 2*adv 

Cheers,
Jean-Michel

PS: I cc to MITgcm-devel so that we keep a record of this (since I always
forgot why we need to double the minimum size of Olx,Oly )

On Sun, Aug 31, 2014 at 01:16:08PM -0400, Samar Khatiwala wrote:
> Hi Jean-Michel,
> 
> Thanks so much for taking the time to reply in such detail. Sorry for my late reply. I've been doing 
> various tests to make sure I understood what was going on. Thanks to your explanation things are 
> making a lot more sense now.
> 
> I asked this question because I've been trying to figure out a small problem with the transport matrix 
> I extracted from the new ECCO v4 state estimate. The way this works is that I initialize a ptracer with 
> a "1" at some point and "0" everywhere else, take one time step, and compute the tendency from 
> advection and horizontal diffusion. (Implicit vertical mixing is treated separately.) In practice I initialize 
> with a pattern of "1"'s and "0"'s rather than a "1" at a single location. As long as the "1"'s are separated 
> by a sufficient distance that depends on the stencil of the advection operator all should be well. This has 
> worked without fail in the past. I usually use DST3 which has a halo of 2 points on either side, that is in a 
> single time step tracer will only spread within a 5 x 5 grid point square centered on the point at which I 
> put a "1".
> 
> But I just discovered that for certain points on the boundary between the Arctic face and Atlantic face of 
> the LLC90 grid the halo is actually 6 x 6 grid points. That turns out to be the source of my problem and it 
> took me a very long time to track it down. I could easily solve my problem this by increasing the halo 
> size in the matrix extraction (it just means more computations), but I would like to understand what is 
> going on. In all these years of using the DST3 routine I've never had this happen. The halo was always 
> 5 x 5. 
> 
> Can you think of an explanation? Is there something different about the LLC90 grid (this is the first time 
> I'm using it). I thought there might be some other package like salt_plume that does some extra transport 
> but that doesn't seem to act on ptracers.
> 
> Again, thanks very much for your help!
> 
> Best,
> 
> Samar
> 
> On Aug 22, 2014, at 3:58 PM, Jean-Michel Campin <jmc at ocean.mit.edu> wrote:
> 
> > Hi Samar,
> > 
> > Ok, I will try to clarify this. 
> > First, when using 
> >> forward-in-time advection schemes like DST3.
> > the Adams-Bashforth is off, which makes it simpler, 
> > so I will continue with this simpler case (no AB in tracer Eq).
> > 
> >> 1) Where in your paper does the FREESURF_RESCALE_G step come in? And what exactly is the rescaling factor? 
> > This is not explicitly written in the OM paper, but in the MITgcm manual,
> > it corresponds to the scaling h^n/h^{n+1} (or V^n/V^{n+1} if you prefer)
> > that can be found in equation (without number) at the end of 2.10.2.3, just after:
> >> Then, in a second step, the thickness variation (expansion/reduction) is taken into account:
> > 
> >> From what I can make out at a given (i,j) point the same term is applied to all points in the vertical.
> > This all depend on which coordinate is used: i
> > in z-coords, it's just applied to the surface tendency:
> >>     IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
> > and with z* coords, the scaling is the same for all pts of the same column.
> > (but with Hybrid sigma-z coords, it also vary with depth).
> > 
> >> 2) I've gone over S/R GAD_ADVECTION, the end result of which seems to be gTracer = (T^{n+1}-T^{n})/dt. 
> > Regarding the multi-dim advection (S/R GAD_ADVECTION), there are 2 versions
> > (within the same file), the original method which corresponds to 
> > #undef GAD_MULTIDIM_COMPRESSIBLE and the alternative method that I call 
> > "compressible flow method" that corresponds to 
> > #define GAD_MULTIDIM_COMPRESSIBLE.
> > This later one was added in March 2013 into gad_advection.F but was already
> > used in gad_som_advect.F (2nd order moment advection scheme, from Prather 86).
> > 
> > The original method follows equation 2.201 - 2.203 of the manual,
> > (except may be a sign error ?) and it match also what you wrote here:
> >> From what I can make out this is gTracer = -1/V^{n} [delta(u A T)^{n} - T^{n} delta(u A)^{n}], 
> >> where delta(x) = x_i+1 - x_i (plus the j and k terms), 
> >> A are the face areas, and V is the cell volume.
> > 
> > The reason to add the term "- T^{n} delta(u A)" is to (almost) recover 
> > the advective form: u.grad(T) = delta(u.T) - T.delta(u) 
> > (this is important because, even if the 3-D flow is non-divergent,
> > the flow in 1 direction is not, and it matters if we want to split
> > the advection by direction)
> > but instead of using t^{n+1/3} in 2.202 we keep using t^{n} so that
> > when adding these 3 contributions (and because of continuity eq), we are 
> > just left with "-w_s T_1^n" (in OM paper, eq 14 & 15 & 16).
> > This way we get a perfect conservation (but a less precise time
> > discretisation).
> > 
> > The "compressible flow method" (#define GAD_MULTIDIM_COMPRESSIBLE) 
> > also conserve perfectly and is also more accurate (time-discretisation).
> > It updates 3 times (one per direction) the volume and the tracer:
> >     V^{n+1/3} = V^n     - dt * delta_i(u A)
> > (V*T)^{n+1/3} = (V*T)^n - dt * delta_i(u A T)
> > It's little bit more complicated than the original method regarding
> > the adjoint, and no attempt have been made to pass it through TAF. 
> > 
> > I hope this will help, but if you have a question, don't hesitate.
> > 
> > Cheers,
> > Jean-Michel
> > 
> > On Thu, Aug 21, 2014 at 05:06:15PM -0400, Samar Khatiwala wrote:
> >> Hi Jean-Michel,
> >> 
> >> I realize you're very busy but I've been stuck trying to understand the time-stepping scheme/algorithm for passive tracer 
> >> advection/diffusion with the nonlinear free surface and exact conservation turned on. I've read both Sec. 2.10.x of the 
> >> documentation and your 2004 Ocean Modelling paper. Of course I've also stared at the code! But I'm still confused 
> >> at the correspondence between your paper and what is implemented in MITgcm. I was wondering if you could please 
> >> clarify the precise algorithm/implementation. I'm particularly interested in the forward-in-time advection schemes 
> >> like DST3.
> >> 
> >> Specifically:
> >> 
> >> 1) Where in your paper does the FREESURF_RESCALE_G step come in? And what exactly is the rescaling factor? From 
> >> what I can make out at a given (i,j) point the same term is applied to all points in the vertical.
> >> 
> >> 2) I've gone over S/R GAD_ADVECTION, the end result of which seems to be gTracer = (T^{n+1}-T^{n})/dt. From what I can 
> >> make out this is gTracer = -1/V^{n} [delta(u A T)^{n} - T^{n} delta(u A)^{n}], where delta(x) = x_i+1 - x_i (plus the j and k terms), 
> >> A are the face areas, and V is the cell volume.
> >> 
> >> If I were to start with tracer conservation for each cell: d(V T)/dt = -delta(u A T) and discretize I *almost* but not quite get to the 
> >> above and I'm trying to understand what exactly this routine is doing and how it connects to the subsequent rescaling step.
> >> 
> >> If it is easier for you I'm happy to give you a call or Skype at your convenience. I would even visit you but I'm now on the other 
> >> side of the Atlantic!
> >> 
> >> Thanks very much!
> >> 
> >> Best,
> >> 
> >> Samar
> >> 
>