[MITgcm-devel] ptracers

[Patrick Heimbach]

1 and 2 may mess up the flow dependencies though,
depending on how TAF does them.
I'll try it incrementally.
-p.

Chris Hill wrote:
> 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
>>
> 
> 
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://dev.mitgcm.org/mailman/listinfo/mitgcm-devel


-- 
Patrick Heimbach ........................... M.I.T
FON: +1/617/253-5259 .......... EAPS, Room 54-1518
FAX: +1/617/253-4464 ..... 77 Massachusetts Avenue
mailto:heimbach at mit.edu ....... Cambridge MA 02139
http://www.mit.edu/~heimbach/ .............. U.S.A