pkg/seaice/seaice_tracer_phys.F

C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_tracer_phys.F,v 1.1 2011/06/07 03:58:23 gforget Exp $
C $Name:  $

#include "SEAICE_OPTIONS.h"

CStartOfInterface
      SUBROUTINE SEAICE_TRACER_PHYS( myTime, myIter, myThid )
C     /=======================================================\
C     | SUBROUTINE seaice_tracer_phys                         |
C     | o Time step SItr/SItrEFF as a result of               |
C     |   seaice thermodynamics and specific tracer physics   |
C     \=======================================================/
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "FFIELDS.h"
#include "DYNVARS.h"
#include "SEAICE_SIZE.h"
#include "SEAICE.h"
#include "SEAICE_PARAMS.h"
#include "SEAICE_TRACER.h"

C     === Routine arguments ===
C     INPUT:
C     myTime  :: Simulation time
C     myIter  :: Simulation timestep number
C     myThid  :: Thread no. that called this routine.
C     OUTPUT:
      _RL myTime
      INTEGER myIter, myThid
CEndOfInterface

C     === Local variables ===
#ifdef ALLOW_SITRACER

      INTEGER iTr, jTh, I, J, bi, bj, ks
      _RL SItrFromOcean  (1:sNx,1:sNy)
      _RL SItrFromFlood   (1:sNx,1:sNy)
      _RL HEFFprev, HEFFpost, growFact, meltPart, tmpscal1
      _RL SItrExpand  (1:sNx,1:sNy)
      _RL AREAprev, AREApost, expandFact

#ifdef ALLOW_SITRACER_DIAG
      _RL DIAGarray     (1:sNx,1:sNy,Nr)
#endif

cgf for now I do not fully account for ocean-ice fluxes of tracer
cgf -> I just prescribe it consistent with age tracer
cgf eventually I will need to handle them as function params

      ks=1

      DO bj=myByLo(myThid),myByHi(myThid)
      DO bi=myBxLo(myThid),myBxHi(myThid)
      DO iTr=1,SItrMaxNum

c 0) set ice-ocean and ice-snow exchange values       
c =============================================
      DO J=1,sNy
       DO I=1,sNx
        SItrFromOcean(i,j)=0. _d 0
        SItrFromFlood(i,j)=0. _d 0
        SItrExpand(i,j)=0. _d 0
       ENDDO
      ENDDO
      if (SItrName(iTr).EQ.'age') then
c age tracer: no age in ocean, or effect from ice cover changes
      elseif (SItrName(iTr).EQ.'salinity') then
c salinity tracer: 
       DO J=1,sNy
        DO I=1,sNx
         SItrFromOcean(i,j)=SIsal0
#ifdef SEAICE_VARIABLE_SALINITY
         if (SIsalFRAC.GT.0.) 
     &   SItrFromOcean(i,j)=SIsalFRAC*salt(I,j,ks,bi,bj)
#endif
c as of now, flooding implies no salt extraction from ocean
        ENDDO
       ENDDO
      elseif (SItrName(iTr).EQ.'one') then
c "ice concentration" tracer that should remain .EQ.1.
       DO J=1,sNy
        DO I=1,sNx
         SItrFromOcean(i,j)=1. _d 0
         SItrFromFlood(i,j)=1. _d 0
        ENDDO
       ENDDO
      endif
c 1) seaice thermodynamics processes
c ==================================
      if (SItrMate(iTr).EQ.'HEFF') then
      DO J=1,sNy
       DO I=1,sNx
        HEFFprev=SItrHEFF(i,j,bi,bj,1)
#ifdef ALLOW_SITRACER_DIAG
        diagArray(I,J,5+(iTr-1)*5) = 
     &    HEFFprev*SItracer(i,j,bi,bj,iTr) + SItrBucket(i,j,bi,bj,iTr)
#endif
c apply the sequence of thermodynamics increments to actual traceur
c (see seaice_growth.F)
c (jTh=1 tendency due to ice-ocean interaction)
c (jTh=2 tendency due to the atmosphere, over ice covered part)
c (jTh=3 tendency due to the atmosphere, over open water part)
c (jTh=4 tendency due to flooding)
        DO jTh=1,3
         HEFFprev=SItrHEFF(i,j,bi,bj,jTh)
         HEFFpost=SItrHEFF(i,j,bi,bj,jTh+1)
c compute ratio in [0. 1.] range for either growth or melt
         growFact=1. _d 0
         meltPart=0. _d 0
         if (HEFFpost.GT.HEFFprev) growFact=HEFFprev/HEFFpost
         if (HEFFpost.LT.HEFFprev) meltPart=HEFFprev-HEFFpost
c update SItr accordingly
         SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*growFact
     &                      +SItrFromOcean(i,j)*(1. _d 0 - growFact)
         SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
     &             -HEFFpost*SItrFromOcean(i,j)*(1. _d 0 - growFact)
         SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
     &                            +meltPart*SItracer(i,j,bi,bj,iTr)
        ENDDO
c apply flooding term
        growFact=1. _d 0
        HEFFprev=SItrHEFF(i,j,bi,bj,4)
        HEFFpost=SItrHEFF(i,j,bi,bj,5)
        if (HEFFpost.GT.HEFFprev) growFact=HEFFprev/HEFFpost
        SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*growFact
     &                     +SItrFromFlood(i,j) *(1. _d 0 - growFact)
c rk: flooding can only imply an ocean-ice tracer exchange, as long
c as we dont have snow tracers, so it goes through SItrBucket.
          SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
     &             -HEFFpost*SItrFromFlood(i,j)*(1. _d 0 - growFact)
#ifdef ALLOW_SITRACER_DIAG
        diagArray(I,J,5+(iTr-1)*5) = HEFFpost*SItracer(i,j,bi,bj,iTr) 
     &  +SItrBucket(i,j,bi,bj,iTr)-diagArray(I,J,5+(iTr-1)*5)
#endif
       ENDDO
      ENDDO
c TAF?      if (SItrMate(iTr).EQ.'AREA') then
      else
c 1) or seaice cover expansion
c ============================
c this is much simpler than for ice volume/mass tracers, because 
c properties of the ice surface are not be conserved across the 
c ocean-ice system, the contraction/expansion terms are all
c simultaneous (which is sane), and the only generic effect 
c is due to expansion (new cover).
      DO J=1,sNy
       DO I=1,sNx
c apply expansion
        AREAprev=SItrAREA(i,j,bi,bj,2)
        AREApost=SItrAREA(i,j,bi,bj,3)
c compute ratio in [0. 1.] range for expansion/contraction
        expandFact=1. _d 0
        if (AREApost.GT.AREAprev) expandFact=AREAprev/AREApost
c update SItr accordingly
         SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*expandFact
     &                      +SItrExpand(i,j)*(1. _d 0 - expandFact)
       ENDDO
      ENDDO
      endif
c 2) very ice tracer processes
c ============================
      if (SItrName(iTr).EQ.'age') then
c age tracer: grow old as time passes by
       DO J=1,sNy
        DO I=1,sNx
          if (( (SItrHEFF(i,j,bi,bj,5).GT.0. _d 0).AND.(SItrMate(iTr)
     &     .EQ.'HEFF') ).OR.( (SItrAREA(i,j,bi,bj,3).GT.0. _d 0).AND.
     &     (SItrMate(iTr).EQ.'AREA') )) then
            SItracer(i,j,bi,bj,iTr)=
     &      SItracer(i,j,bi,bj,iTr)+SEAICE_deltaTtherm
          else
            SItracer(i,j,bi,bj,iTr)=0. _d 0
          endif
        ENDDO
       ENDDO
      elseif (SItrName(iTr).EQ.'salinity') then
c salinity tracer: no specific process
      elseif (SItrName(iTr).EQ.'one') then
c "ice concentration" tracer: no specific process
      elseif (SItrName(iTr).EQ.'ridge') then
c simple, made up, ice surface roughness index prototype
#ifndef SEAICE_GROWTH_LEGACY
       DO J=1,sNy
        DO I=1,sNx
c ridging increases roughness
          SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)+
     &    MAX(0. _d 0, SItrAREA(i,j,bi,bj,1)-SItrAREA(i,j,bi,bj,2))
c ice melt reduces ridges/roughness
          HEFFprev=SItrHEFF(i,j,bi,bj,1)
          HEFFpost=SItrHEFF(i,j,bi,bj,4)
          tmpscal1=1. _d 0
          if (HEFFprev.GT.HEFFpost) tmpscal1=HEFFpost/HEFFprev
          SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*tmpscal1
        ENDDO
       ENDDO
#endif
      endif
c 3) ice-ocean tracer exchange/mapping to external variables
c ==========================================================
      if (SItrName(iTr).EQ.'age') then
c age tracer: not passed to ocean
      elseif (SItrName(iTr).EQ.'salinity') then
c salinity tracer: salt flux
c      DO J=1,sNy
c       DO I=1,sNx
c        saltFlux(I,J,bi,bj) = - SItrBucket(i,j,bi,bj,iTr)
c    &     *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce
c note: at this point of the time step, that is the correct sign
c        saltPlumeFlux(I,J,bi,bj) = ...
c       ENDDO
c      ENDDO
      elseif (SItrName(iTr).EQ.'one') then
c "ice concentration" tracer: not passed to ocean
      endif
      DO J=1,sNy
       DO I=1,sNx
#ifdef ALLOW_SITRACER_DIAG
        diagArray(I,J,4+(iTr-1)*5) = - SItrBucket(i,j,bi,bj,iTr)
     &  *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce
#endif
c empty bucket
        SItrBucket(i,j,bi,bj,iTr)=0. _d 0
       ENDDO
      ENDDO
c TAF? elseif (SItrMate(iTr).EQ.'AREA') then
c 4) diagnostics
c ==============
#ifdef ALLOW_SITRACER_DIAG
      if (SItrMate(iTr).EQ.'HEFF') then
      DO J=1,sNy
       DO I=1,sNx
        HEFFpost=SItrHEFF(i,j,bi,bj,5)
        DIAGarray(I,J,1+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)
        DIAGarray(I,J,2+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)*HEFFpost
c diagArray(:,:,3) is the term of comparison for diagArray(:,:,2)
        if (SItrName(iTr).EQ.'age') then
          DIAGarray(I,J,3+(iTr-1)*5) = IceAgeTr(i,j,bi,bj,2)
        elseif (SItrName(iTr).EQ.'salinity') then
          DIAGarray(I,J,3+(iTr-1)*5) = HSALT(i,j,bi,bj)/SEAICE_rhoIce
        elseif (SItrName(iTr).EQ.'one') then
          DIAGarray(I,J,3+(iTr-1)*5) = HEFFpost
        endif
c diagArray(:,:,4) allows check of conservation : del(SItrBucket)+del(SItr*HEFF)=0. over do_phys
c diagArray(:,:,5) is the tracer flux from the ocean (<0 incr. ocean tracer)
       ENDDO
      ENDDO
      else
      DO J=1,sNy
       DO I=1,sNx
        AREApost=SItrAREA(i,j,bi,bj,3)
        DIAGarray(I,J,1+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)
        DIAGarray(I,J,2+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)*AREApost
c diagArray(:,:,3) is the term of comparison for diagArray(:,:,2)
        if (SItrName(iTr).EQ.'age') then
          DIAGarray(I,J,3+(iTr-1)*5) = IceAgeTr(i,j,bi,bj,1)
        endif
       ENDDO
      ENDDO
      endif
#endif
      ENDDO
#ifdef ALLOW_SITRACER_DIAG
      CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1  ',0,Nr,3,bi,bj,myThid)
#endif
      ENDDO
      ENDDO

#endif /* ALLOW_SITRACER */

      RETURN
      END
1	C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_tracer_phys.F,v 1.1 2011/06/07 03:58:23 gforget Exp $
2	C $Name: $
3
4	#include "SEAICE_OPTIONS.h"
5
6	CStartOfInterface
7	SUBROUTINE SEAICE_TRACER_PHYS( myTime, myIter, myThid )
8	C /=======================================================\
9	C \| SUBROUTINE seaice_tracer_phys \|
10	C \| o Time step SItr/SItrEFF as a result of \|
11	C \| seaice thermodynamics and specific tracer physics \|
12	C \=======================================================/
13	IMPLICIT NONE
14
15	C === Global variables ===
16	#include "SIZE.h"
17	#include "EEPARAMS.h"
18	#include "FFIELDS.h"
19	#include "DYNVARS.h"
20	#include "SEAICE_SIZE.h"
21	#include "SEAICE.h"
22	#include "SEAICE_PARAMS.h"
23	#include "SEAICE_TRACER.h"
24
25	C === Routine arguments ===
26	C INPUT:
27	C myTime :: Simulation time
28	C myIter :: Simulation timestep number
29	C myThid :: Thread no. that called this routine.
30	C OUTPUT:
31	_RL myTime
32	INTEGER myIter, myThid
33	CEndOfInterface
34
35	C === Local variables ===
36	#ifdef ALLOW_SITRACER
37
38	INTEGER iTr, jTh, I, J, bi, bj, ks
39	_RL SItrFromOcean (1:sNx,1:sNy)
40	_RL SItrFromFlood (1:sNx,1:sNy)
41	_RL HEFFprev, HEFFpost, growFact, meltPart, tmpscal1
42	_RL SItrExpand (1:sNx,1:sNy)
43	_RL AREAprev, AREApost, expandFact
44
45	#ifdef ALLOW_SITRACER_DIAG
46	_RL DIAGarray (1:sNx,1:sNy,Nr)
47	#endif
48
49	cgf for now I do not fully account for ocean-ice fluxes of tracer
50	cgf -> I just prescribe it consistent with age tracer
51	cgf eventually I will need to handle them as function params
52
53	ks=1
54
55	DO bj=myByLo(myThid),myByHi(myThid)
56	DO bi=myBxLo(myThid),myBxHi(myThid)
57	DO iTr=1,SItrMaxNum
58
59	c 0) set ice-ocean and ice-snow exchange values
60	c =============================================
61	DO J=1,sNy
62	DO I=1,sNx
63	SItrFromOcean(i,j)=0. _d 0
64	SItrFromFlood(i,j)=0. _d 0
65	SItrExpand(i,j)=0. _d 0
66	ENDDO
67	ENDDO
68	if (SItrName(iTr).EQ.'age') then
69	c age tracer: no age in ocean, or effect from ice cover changes
70	elseif (SItrName(iTr).EQ.'salinity') then
71	c salinity tracer:
72	DO J=1,sNy
73	DO I=1,sNx
74	SItrFromOcean(i,j)=SIsal0
75	#ifdef SEAICE_VARIABLE_SALINITY
76	if (SIsalFRAC.GT.0.)
77	& SItrFromOcean(i,j)=SIsalFRAC*salt(I,j,ks,bi,bj)
78	#endif
79	c as of now, flooding implies no salt extraction from ocean
80	ENDDO
81	ENDDO
82	elseif (SItrName(iTr).EQ.'one') then
83	c "ice concentration" tracer that should remain .EQ.1.
84	DO J=1,sNy
85	DO I=1,sNx
86	SItrFromOcean(i,j)=1. _d 0
87	SItrFromFlood(i,j)=1. _d 0
88	ENDDO
89	ENDDO
90	endif
91	c 1) seaice thermodynamics processes
92	c ==================================
93	if (SItrMate(iTr).EQ.'HEFF') then
94	DO J=1,sNy
95	DO I=1,sNx
96	HEFFprev=SItrHEFF(i,j,bi,bj,1)
97	#ifdef ALLOW_SITRACER_DIAG
98	diagArray(I,J,5+(iTr-1)*5) =
99	& HEFFprev*SItracer(i,j,bi,bj,iTr) + SItrBucket(i,j,bi,bj,iTr)
100	#endif
101	c apply the sequence of thermodynamics increments to actual traceur
102	c (see seaice_growth.F)
103	c (jTh=1 tendency due to ice-ocean interaction)
104	c (jTh=2 tendency due to the atmosphere, over ice covered part)
105	c (jTh=3 tendency due to the atmosphere, over open water part)
106	c (jTh=4 tendency due to flooding)
107	DO jTh=1,3
108	HEFFprev=SItrHEFF(i,j,bi,bj,jTh)
109	HEFFpost=SItrHEFF(i,j,bi,bj,jTh+1)
110	c compute ratio in [0. 1.] range for either growth or melt
111	growFact=1. _d 0
112	meltPart=0. _d 0
113	if (HEFFpost.GT.HEFFprev) growFact=HEFFprev/HEFFpost
114	if (HEFFpost.LT.HEFFprev) meltPart=HEFFprev-HEFFpost
115	c update SItr accordingly
116	SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*growFact
117	& +SItrFromOcean(i,j)*(1. _d 0 - growFact)
118	SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
119	& -HEFFpostSItrFromOcean(i,j)(1. _d 0 - growFact)
120	SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
121	& +meltPart*SItracer(i,j,bi,bj,iTr)
122	ENDDO
123	c apply flooding term
124	growFact=1. _d 0
125	HEFFprev=SItrHEFF(i,j,bi,bj,4)
126	HEFFpost=SItrHEFF(i,j,bi,bj,5)
127	if (HEFFpost.GT.HEFFprev) growFact=HEFFprev/HEFFpost
128	SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*growFact
129	& +SItrFromFlood(i,j) *(1. _d 0 - growFact)
130	c rk: flooding can only imply an ocean-ice tracer exchange, as long
131	c as we dont have snow tracers, so it goes through SItrBucket.
132	SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
133	& -HEFFpostSItrFromFlood(i,j)(1. _d 0 - growFact)
134	#ifdef ALLOW_SITRACER_DIAG
135	diagArray(I,J,5+(iTr-1)5) = HEFFpostSItracer(i,j,bi,bj,iTr)
136	& +SItrBucket(i,j,bi,bj,iTr)-diagArray(I,J,5+(iTr-1)*5)
137	#endif
138	ENDDO
139	ENDDO
140	c TAF? if (SItrMate(iTr).EQ.'AREA') then
141	else
142	c 1) or seaice cover expansion
143	c ============================
144	c this is much simpler than for ice volume/mass tracers, because
145	c properties of the ice surface are not be conserved across the
146	c ocean-ice system, the contraction/expansion terms are all
147	c simultaneous (which is sane), and the only generic effect
148	c is due to expansion (new cover).
149	DO J=1,sNy
150	DO I=1,sNx
151	c apply expansion
152	AREAprev=SItrAREA(i,j,bi,bj,2)
153	AREApost=SItrAREA(i,j,bi,bj,3)
154	c compute ratio in [0. 1.] range for expansion/contraction
155	expandFact=1. _d 0
156	if (AREApost.GT.AREAprev) expandFact=AREAprev/AREApost
157	c update SItr accordingly
158	SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*expandFact
159	& +SItrExpand(i,j)*(1. _d 0 - expandFact)
160	ENDDO
161	ENDDO
162	endif
163	c 2) very ice tracer processes
164	c ============================
165	if (SItrName(iTr).EQ.'age') then
166	c age tracer: grow old as time passes by
167	DO J=1,sNy
168	DO I=1,sNx
169	if (( (SItrHEFF(i,j,bi,bj,5).GT.0. _d 0).AND.(SItrMate(iTr)
170	& .EQ.'HEFF') ).OR.( (SItrAREA(i,j,bi,bj,3).GT.0. _d 0).AND.
171	& (SItrMate(iTr).EQ.'AREA') )) then
172	SItracer(i,j,bi,bj,iTr)=
173	& SItracer(i,j,bi,bj,iTr)+SEAICE_deltaTtherm
174	else
175	SItracer(i,j,bi,bj,iTr)=0. _d 0
176	endif
177	ENDDO
178	ENDDO
179	elseif (SItrName(iTr).EQ.'salinity') then
180	c salinity tracer: no specific process
181	elseif (SItrName(iTr).EQ.'one') then
182	c "ice concentration" tracer: no specific process
183	elseif (SItrName(iTr).EQ.'ridge') then
184	c simple, made up, ice surface roughness index prototype
185	#ifndef SEAICE_GROWTH_LEGACY
186	DO J=1,sNy
187	DO I=1,sNx
188	c ridging increases roughness
189	SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)+
190	& MAX(0. _d 0, SItrAREA(i,j,bi,bj,1)-SItrAREA(i,j,bi,bj,2))
191	c ice melt reduces ridges/roughness
192	HEFFprev=SItrHEFF(i,j,bi,bj,1)
193	HEFFpost=SItrHEFF(i,j,bi,bj,4)
194	tmpscal1=1. _d 0
195	if (HEFFprev.GT.HEFFpost) tmpscal1=HEFFpost/HEFFprev
196	SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)*tmpscal1
197	ENDDO
198	ENDDO
199	#endif
200	endif
201	c 3) ice-ocean tracer exchange/mapping to external variables
202	c ==========================================================
203	if (SItrName(iTr).EQ.'age') then
204	c age tracer: not passed to ocean
205	elseif (SItrName(iTr).EQ.'salinity') then
206	c salinity tracer: salt flux
207	c DO J=1,sNy
208	c DO I=1,sNx
209	c saltFlux(I,J,bi,bj) = - SItrBucket(i,j,bi,bj,iTr)
210	c & HEFFM(I,J,bi,bj)/SEAICE_deltaTthermSEAICE_rhoIce
211	c note: at this point of the time step, that is the correct sign
212	c saltPlumeFlux(I,J,bi,bj) = ...
213	c ENDDO
214	c ENDDO
215	elseif (SItrName(iTr).EQ.'one') then
216	c "ice concentration" tracer: not passed to ocean
217	endif
218	DO J=1,sNy
219	DO I=1,sNx
220	#ifdef ALLOW_SITRACER_DIAG
221	diagArray(I,J,4+(iTr-1)*5) = - SItrBucket(i,j,bi,bj,iTr)
222	& HEFFM(I,J,bi,bj)/SEAICE_deltaTthermSEAICE_rhoIce
223	#endif
224	c empty bucket
225	SItrBucket(i,j,bi,bj,iTr)=0. _d 0
226	ENDDO
227	ENDDO
228	c TAF? elseif (SItrMate(iTr).EQ.'AREA') then
229	c 4) diagnostics
230	c ==============
231	#ifdef ALLOW_SITRACER_DIAG
232	if (SItrMate(iTr).EQ.'HEFF') then
233	DO J=1,sNy
234	DO I=1,sNx
235	HEFFpost=SItrHEFF(i,j,bi,bj,5)
236	DIAGarray(I,J,1+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)
237	DIAGarray(I,J,2+(iTr-1)5) = SItracer(i,j,bi,bj,iTr)HEFFpost
238	c diagArray(:,:,3) is the term of comparison for diagArray(:,:,2)
239	if (SItrName(iTr).EQ.'age') then
240	DIAGarray(I,J,3+(iTr-1)*5) = IceAgeTr(i,j,bi,bj,2)
241	elseif (SItrName(iTr).EQ.'salinity') then
242	DIAGarray(I,J,3+(iTr-1)*5) = HSALT(i,j,bi,bj)/SEAICE_rhoIce
243	elseif (SItrName(iTr).EQ.'one') then
244	DIAGarray(I,J,3+(iTr-1)*5) = HEFFpost
245	endif
246	c diagArray(:,:,4) allows check of conservation : del(SItrBucket)+del(SItr*HEFF)=0. over do_phys
247	c diagArray(:,:,5) is the tracer flux from the ocean (<0 incr. ocean tracer)
248	ENDDO
249	ENDDO
250	else
251	DO J=1,sNy
252	DO I=1,sNx
253	AREApost=SItrAREA(i,j,bi,bj,3)
254	DIAGarray(I,J,1+(iTr-1)*5) = SItracer(i,j,bi,bj,iTr)
255	DIAGarray(I,J,2+(iTr-1)5) = SItracer(i,j,bi,bj,iTr)AREApost
256	c diagArray(:,:,3) is the term of comparison for diagArray(:,:,2)
257	if (SItrName(iTr).EQ.'age') then
258	DIAGarray(I,J,3+(iTr-1)*5) = IceAgeTr(i,j,bi,bj,1)
259	endif
260	ENDDO
261	ENDDO
262	endif
263	#endif
264	ENDDO
265	#ifdef ALLOW_SITRACER_DIAG
266	CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1 ',0,Nr,3,bi,bj,myThid)
267	#endif
268	ENDDO
269	ENDDO
270
271	#endif /* ALLOW_SITRACER */
272
273	RETURN
274	END