[MITgcm-devel] ptracers

Patrick Heimbach heimbach at MIT.EDU
Fri Oct 24 02:01:47 EDT 2003


Want to bet an ice cream that fixing 3 alone
doesn't do the trick?
-p.

Chris Hill wrote:
> You're right anything is possible! But I would try 3 first.
> 
> 
>>-----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 11:27 PM
>>To: MITgcm-devel at mitgcm.org
>>Subject: Re: [MITgcm-devel] ptracers
>>
>>
>>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
> 
> _______________________________________________
> 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




More information about the MITgcm-devel mailing list