pkg/seaice/seaice_tracer_phys.F

C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_tracer_phys.F,v 1.7 2012/02/16 03:32:39 jmc 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"
#ifdef ALLOW_SALT_PLUME
# include "SALT_PLUME.h"
#endif

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
      CHARACTER*8   diagName

#ifdef ALLOW_SITRACER_DEBUG_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,SItrNumInUse

c 0) set ice-ocean and ice-snow exchange values
c =============================================
      DO J=1,sNy
       DO I=1,sNx
        SItrFromOcean(i,j)=SItrFromOcean0(iTr)
        SItrFromFlood(i,j)=SItrFromFlood0(iTr)
        SItrExpand(i,j)=SItrExpand0(iTr)
       ENDDO
      ENDDO
c salinity tracer:
      if ( (SItrName(iTr).EQ.'salinity').AND.
     &      (SItrFromOceanFrac(iTr).GT.ZERO) ) then
       DO J=1,sNy
        DO I=1,sNx
         SItrFromOcean(i,j)=SItrFromOceanFrac(iTr)*salt(I,j,ks,bi,bj)
         SItrFromFlood(i,j)=SItrFromFloodFrac(iTr)*salt(I,j,ks,bi,bj)
        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_DEBUG_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_DEBUG_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 ==========================================================
#ifdef ALLOW_DIAGNOSTICS
      if (SItrMate(iTr).EQ.'HEFF') then
        WRITE(diagName,'(A4,I2.2,A2)') 'SItr',iTr,'Fx'
        tmpscal1=-ONE/SEAICE_deltaTtherm*SEAICE_rhoIce
        CALL DIAGNOSTICS_SCALE_FILL(SItrBucket(1-OLx,1-OLy,bi,bj,iTr),
     &   tmpscal1, 1, diagName,0,1,2,bi,bj,myThid)
      endif
#endif

      if ( (SItrName(iTr).EQ.'salinity').AND.
     &     (SEAICE_salinityTracer) ) then
c salinity tracer: salt flux
        DO J=1,sNy
         DO I=1,sNx
          saltFlux(I,J,bi,bj) = - SItrBucket(i,j,bi,bj,iTr)
     &      *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce
c note: at this point of the time step, that is the correct sign
#ifdef ALLOW_SALT_PLUME
c should work for both constant and variable ice salinity -- to be tested
          saltPlumeFlux(I,J,bi,bj) = MAX(ZERO,saltFlux(I,J,bi,bj))
     &      *SPsalFRAC*(salt(I,j,ks,bi,bj)-SItrFromOcean(i,j))
#endif
         ENDDO
        ENDDO
      endif

      DO J=1,sNy
       DO I=1,sNx
#ifdef ALLOW_SITRACER_DEBUG_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_DEBUG_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.'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
       ENDDO
      ENDDO
      endif
#endif
      ENDDO
#ifdef ALLOW_SITRACER_DEBUG_DIAG
c     CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1  ',0,Nr,3,bi,bj,myThid)
#endif
      ENDDO
      ENDDO

#endif /* ALLOW_SITRACER */

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