[MITgcm-devel] changes in tracer vertical implicit scheme

[Patrick Heimbach]

Hi Jean-Michel,

it appears that the broken impl.vert.adv. adjoint is not related to your changes below,
but happened much earlier (Oct 13, 2011), and we didn't notice since it isn't tested yet
(the culprit left the country just in time ;)

Reverting to solve_pentadiagonal.F revision 1.4
fixes the problems with the adjoint.

The changes to vectorize this one need to be revisited.
We should also start testing this in global_ocean.cs32x15

-Patrick

On Feb 1, 2012, at 5:20 PM, Jean-Michel Campin wrote:

> Hi Patrick,
> 
> The 2 majors changes (on Dec 01, 2011) are:
> 1) new argument "recip_hFacNew" computed in thermodynamics.F
> and passed to gad_implicit_r.F , and then passed to the selected
> gad_*_impl_r.F S/R ;
> 2) the term: tracer* del_k(wtrans)* deltaT has been moved from
> gad_implicit_r.F (in version 1.14i of gad_implicit_r.F, lines 210 - 213):
>           gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
>     &      + deltaTLev(1)*recip_rA(i,j,bi,bj)
>     &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
>     &       *tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))*rkSign
> to gad_calc_rhs.F :
> in gad_calc_rhs.F, now we have line 164:
>    IF (implicitAdvection) rAdvFac = rkSign
> but before it was:
>    IF (implicitAdvection) rAdvFac = 0. _d 0
> 
> What could be tried is to have only 1 of the 2 changes, and see 
> how this goes.
> 
> Cheers,
> Jean-Michel
> 
> On Wed, Feb 01, 2012 at 12:01:16PM -0500, Patrick Heimbach wrote:
>> 
>> Jean-Michel,
>> 
>> for the record, are you aware that your changes to tracer vertical implicit scheme
>> break the adjoint (which we now use in several production calculations)?
>> 
>> Easy way to test  (I had thought it's already tested by default) is in
>> verification/global_ocean.cs32x15/input_ad/ 
>> to set
>> tempImplVertAdv=.TRUE.,
>> saltImplVertAdv=.TRUE.,
>> tempVertAdvScheme=3,
>> saltVertAdvScheme=3,
>> 
>> If you point me to likely culprit routines, I'll take a look.
>> 
>> Thanks
>> -Patrick
>> 
>> On Dec 1, 2011, at 9:43 AM, Jean-Michel Campin wrote:
>> 
>>> Hi,
>>> 
>>> Some fix/changes related to tracer vertical direction implicit method:
>>> 
>>> 1) implicit vertical advection conservation with Adams-Bashforth 
>>> and/or non-lin free-surf:
>>> the term Tr*d/dz(w) was not going through AB and neither was
>>> rescaled (NonLin-FreeSurf) and this was breaking tracer conservation.
>>> I remove Tr*d/dz(w) in gad_implicit_r.F and add it in gad_calc_rhs.F;
>>> this fix conservation in both cases.
>>> 2) implicit vertical diffusion or advection with non-lin free-surf:
>>> Was using current recip_hFacC instead of future recip_hFacC.
>>> This lead to small error when using NonLin-FreeSurf and also break 
>>> conservation with NonLin-FreeSurf in r-coordinate (but not with r*).
>>> The right recip_hFac is now computed in thermodynamics.F and passed
>>> to all vertical implicit S/R.
>>> 
>>> I also added diagnostics for gtNm1 & gsNm1, so that it's easier
>>> to check tracer conservation with NonLin-FreeSurf + AB (the quantity
>>> being conserved is Tr + (0.5+abEps)*deltaT*gtNm1).
>>> And will add a secondary test to advect_xz to test implicit vertical
>>> advection (+ NonLin-FreeSurf).
>>> 
>>> I don't know why those changes affect bottom_ctrl_5x5 (adjoint gradient)
>>> and break agreement with forward gradient.
>>> I removed all the implicit vertical code in thermodynamics.F (checked-in),
>>> and it's better, but still not as good than what it used to be.
>>> If someone has time to take a look at this experiment.
>>> 
>>> Cheers,
>>> Jean-Michel
>>> 
>>> _______________________________________________
>>> MITgcm-devel mailing list
>>> MITgcm-devel at mitgcm.org
>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>> 
>> ---
>> Patrick Heimbach | heimbach at mit.edu | http://www.mit.edu/~heimbach
>> MIT | EAPS 54-1420 | 77 Massachusetts Ave | Cambridge MA 02139 USA
>> FON +1-617-253-5259 | FAX +1-617-253-4464 | SKYPE patrick.heimbach
>> 
>> 
>> 
>> _______________________________________________
>> 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

---
Patrick Heimbach | heimbach at mit.edu | http://www.mit.edu/~heimbach
MIT | EAPS 54-1420 | 77 Massachusetts Ave | Cambridge MA 02139 USA
FON +1-617-253-5259 | FAX +1-617-253-4464 | SKYPE patrick.heimbach