[MITgcm-devel] ptracers

[Chris Hill]

That's great. Lets hope that 3 fixes it. 

1 and 2 will be tidied up. Its very unlikely they could ever cause a
problem as is because all references within the routines and lower down
in the call trees are only to the dummy args and there are no refs to
the common members. However, it is messy and unecessary to have the
global refs visible so it should be cleaned up.

Chris

> -----Original Message-----
> From: mitgcm-devel-bounces at mitgcm.org 
> [mailto:mitgcm-devel-bounces at mitgcm.org] On Behalf Of Patrick Heimbach
> Sent: Thursday, October 23, 2003 10:51 PM
> To: MITgcm-devel at mitgcm.org
> Subject: [MITgcm-devel] ptracers
> 
> 
> 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 
> ###################################################
> 
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org 
> http://dev.mitgcm.org/mailman/listinfo/mitgcm-> devel
>