--- MITgcm/pkg/seaice/seaice_tracer_phys.F	2011/06/09 19:37:01	1.2
+++ MITgcm/pkg/seaice/seaice_tracer_phys.F	2014/08/23 20:22:16	1.11
@@ -1,26 +1,30 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_tracer_phys.F,v 1.2 2011/06/09 19:37:01 gforget Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_tracer_phys.F,v 1.11 2014/08/23 20:22:16 torge 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     \=======================================================/
+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 "PARAMS.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:
@@ -41,8 +45,9 @@
       _RL HEFFprev, HEFFpost, growFact, meltPart, tmpscal1
       _RL SItrExpand  (1:sNx,1:sNy)
       _RL AREAprev, AREApost, expandFact
+      CHARACTER*8   diagName
 
-#ifdef ALLOW_SITRACER_DIAG
+#ifdef ALLOW_SITRACER_DEBUG_DIAG
       _RL DIAGarray     (1:sNx,1:sNy,Nr)
 #endif
 
@@ -54,37 +59,24 @@
 
       DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
-      DO iTr=1,SItrMaxNum
+      DO iTr=1,SItrNumInUse
 
-c 0) set ice-ocean and ice-snow exchange values       
+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
+        SItrFromOcean(i,j)=SItrFromOcean0(iTr)
+        SItrFromFlood(i,j)=SItrFromFlood0(iTr)
+        SItrExpand(i,j)=SItrExpand0(iTr)
        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: 
+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)=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
+         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
@@ -94,11 +86,11 @@
       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) = 
+#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 apply the sequence of thermodynamics increments to actual tracer
 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)
@@ -131,9 +123,9 @@
 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)
+#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
@@ -141,10 +133,10 @@
       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 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 simultaneous (which is sane), and the only generic effect
 c is due to expansion (new cover).
       DO J=1,sNy
        DO I=1,sNx
@@ -182,7 +174,6 @@
 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
@@ -196,55 +187,67 @@
           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
+#ifdef ALLOW_DIAGNOSTICS
+      IF ( useDiagnostics .AND. 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
-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
+        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
-c        saltPlumeFlux(I,J,bi,bj) = ...
-c       ENDDO
-c      ENDDO
-      elseif (SItrName(iTr).EQ.'one') then
-c "ice concentration" tracer: not passed to ocean
+#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_DIAG
-        diagArray(I,J,4+(iTr-1)*5) = - SItrBucket(i,j,bi,bj,iTr)
+#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
+c  but not for 'grease' (see seaice_growth.F)
+        if (SItrName(iTr).NE.'grease') 
+     &      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
+#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.'age') then
-          DIAGarray(I,J,3+(iTr-1)*5) = IceAgeTr(i,j,bi,bj,2)
-        elseif (SItrName(iTr).EQ.'salinity') then
+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)
+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
@@ -253,17 +256,13 @@
         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)
+#ifdef ALLOW_SITRACER_DEBUG_DIAG
+c     CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG1  ',0,Nr,3,bi,bj,myThid)
 #endif
       ENDDO
       ENDDO