[MITgcm-devel] ptracers

[Patrick Heimbach]

Hi there,

I'm copying this message from Ralf.
He was looking into why the derivative code
of ptracers gave wrong gradient checks.

I'll hopefully check in the modif's tomorrow
together with AD-related modif.s for genmake2.

We need to discuss point 2.

-Patrick

--------------------------------------------------

Hi,

I have worked on the problem in the ptracer adjoint.
Thereby I found three problematic things in the code.
The third one causes the adjoint to be wrong.

Please find enclosed my changes to the code which fixes 1) and 3)
but not 2).

Ralf


1) Aliasing problem in ptracers_forcing.F

C !INTERFACE: ==========================================================
       SUBROUTINE PTRACERS_FORCING(
      I                            bi,bj,k,iTracer,
      U                            gPtracer,
      I                            myIter,myTime,myThid )

C !DESCRIPTION:
C     Adds sources and sinks of passive tracers to the tendancy arrays

C !USES: ===============================================================
       IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PTRACERS.h"
#include "PARAMS.h"
#include "FFIELDS.h"
#include "GRID.h"

It is not necessary to include PTRACERS.h and it even violates the Fortran standard.
Thats because the variable gPtr is passed as argument
and becomes dummy argument gPtracer.
Now there are local arrays gPtr and gPtracer being the same array!
This aliasing is not allowed if one of them is modified, which is the case here.
Fix: just remove #include PTRACERS.h  .


2) Similar in GAD_CALC_RHS, gTracer is dummy argument and is aliased to gS
by the call in calc_gs and aliased to gT by the call in calc_gt.
But here
   #include "DYNVARS.h"
cannot be removed since the arrays UVEL, VVEL, and WVEL are needed.
Fix: not as simple, split DYNVARS.h, pass UVEL, VVEL, and WVEL,
or use tracerIdentity to identify which array is to be used
instead of dummy argument ptracer.


3) In ptracers_integrate.F there is a local array rFlx which is only
initialized for k.EQ.Nr . My guess it that it is updated for all
k.lt.Nr . But in this case it must either be saved or declared in
thermodynamics and passed as argument to ptracers_integrate
(similar to fVerS in thermodynamics) otherwise its not initialised
and used if not k.EQ.Nr !

C  rFlx			:: vertical flux
       INTEGER i,j,iTracer
       INTEGER iMin,iMax,jMin,jMax
       INTEGER kUp,kDown,km1
       _RL rFlx(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num)
       LOGICAL calcAdvection
CEOP

C Loop over tracers
       DO iTracer=1,PTRACERS_numInUse

C Initialize vertical flux to zero and set no-flux across k=Nr+1
        IF (k.EQ.Nr) THEN
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           rFlx(i,j,1,iTracer)=0.
           rFlx(i,j,2,iTracer)=0.
          ENDDO
         ENDDO
        ENDIF


-- 
###################################################
Dr. Ralf Giering
NEW ADDRESS:
FastOpt, Schanzenstr. 36, 20357 Hamburg, Germany
Tel.: +49 40 48096347
Fax : +49 40 48096357
Email: Ralf.Giering at FastOpt.de
URL : http://www.FastOpt.de
###################################################