C $Header: /u/gcmpack/MITgcm/pkg/salt_plume/salt_plume_apply.F,v 1.2 2014/05/22 22:55:07 jmc Exp $
C $Name:  $

#include "SALT_PLUME_OPTIONS.h"

CBOP
C     !ROUTINE: SALT_PLUME_APPLY
C     !INTERFACE:
      SUBROUTINE SALT_PLUME_APPLY(
     I            trIdentity, bi, bj,
     I            recip_hFac_arg,
     I            tracer,trApplyFlag,
     I            myTime, myIter, myThid )

C     U            trStar,

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE SALT_PLUME_APPLY
C     | o Apply the salt_pume-transport to tracer field
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "GRID.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "SALT_PLUME.h"
#ifdef ALLOW_GENERIC_ADVDIFF
# include "GAD.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     trIdentity :: tracer identification number
C     bi,bj      :: Tile indices
C     recip_drF  :: Reciprol of cell thickness
C     recip_hFac_arg :: Reciprol of cell open-depth factor
C     tracer     :: tracer field at current time (input)
C     trApplyFlag:: [0]=update saltplume forcing T/S terms
C                :: [1]=update gTr tendency
C     trStar     :: future tracer field (modified)
C     myTime     :: Current time in simulation
C     myIter     :: Current time-step number
C     myThid     :: my Thread Id. number

      INTEGER trIdentity, trApplyFlag
      INTEGER bi, bj
      _RS recip_hFac_arg(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL tracer        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
C      _RL trStar        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)

      _RL     myTime
      INTEGER myIter, myThid

#ifdef ALLOW_SALT_PLUME
#ifdef SALT_PLUME_VOLUME

C     !LOCAL VARIABLES:
C     === Local variables ===
C     msgBuf     :: Informational/error message buffer
C     plumetend  :: forcing terms [W/m2 or kg/m2/s*psu]
C     work       :: working array
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      INTEGER i, j, k
      INTEGER upward
      LOGICAL onOffFlag
      _RL gTr_Surf2kLev, gTr_Below2kLev, gTr_kLev2Above,
     &    dSPvolBelow2kLev, gTr_totSurf2Below,
     &    SurfVal, SEAICE_Tfrz, ConvertFac, recip_ConvertFac
      integer kp1, Nrp1
      _RL         plumetend(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL         work(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)

#ifdef ALLOW_DIAGNOSTICS
      CHARACTER*8 diagName
      CHARACTER*4 diagSufx
      LOGICAL     doDiagSPtend
C-    Functions:
      LOGICAL  DIAGNOSTICS_IS_ON
      EXTERNAL 
#ifdef ALLOW_GENERIC_ADVDIFF
      CHARACTER*5 GAD_DIAG_SUFX
      EXTERNAL    
#endif /* ALLOW_GENERIC_ADVDIFF */
#endif /* ALLOW_DIAGNOSTICS */

CEOP

C      WRITE(msgBuf,'(A,2i4)') 'trApplyFlag,trIdentity: ',
C     & trApplyFlag,trIdentity
C         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
C     &                      SQUEEZE_RIGHT, myThid )

      IF ( trApplyFlag.LT.0 .OR. trApplyFlag.GT.1) THEN
       STOP 'S/R SALT_PLUME_APPLY: incorrect setting of trApplyFlag!'
      ELSE

       SEAICE_Tfrz = -1.96 _d 0

       onOffFlag = .FALSE.
#ifdef ALLOW_GENERIC_ADVDIFF
       IF ( trIdentity.EQ.GAD_TEMPERATURE ) onOffFlag = .TRUE.
       IF ( trIdentity.EQ.GAD_SALINITY    ) onOffFlag = .TRUE.
#endif
       IF ( onOffFlag ) THEN
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

c       upward = rkSign*NINT(-gravitySign)
        upward = 1
        IF (usingZCoords) upward = -1

        IF ( trIdentity.EQ.GAD_TEMPERATURE ) THEN
         SurfVal          = SEAICE_Tfrz
         ConvertFac       = HeatCapacity_Cp*rhoConst
         recip_ConvertFac = mass2rUnit/HeatCapacity_Cp
#ifdef ALLOW_DIAGNOSTICS
         IF ( useDiagnostics ) diagSufx = '_TH '
#endif /* ALLOW_DIAGNOSTICS */
        ENDIF
        IF ( trIdentity.EQ.GAD_SALINITY ) THEN
         SurfVal          = SPbrineSconst
         ConvertFac       = rhoConst
         recip_ConvertFac = mass2rUnit
#ifdef ALLOW_DIAGNOSTICS
         IF ( useDiagnostics ) diagSufx = '_SLT'
#endif /* ALLOW_DIAGNOSTICS */
        ENDIF

#ifdef ALLOW_DIAGNOSTICS
        doDiagSPtend = .FALSE.
        diagName = 'SPtd'
        IF ( useDiagnostics ) THEN
C--   Set diagnostic suffix for the current tracer
#ifdef ALLOW_GENERIC_ADVDIFF
         diagSufx = GAD_DIAG_SUFX( trIdentity, myThid )
#endif /* ALLOW_GENERIC_ADVDIFF */
         diagName = 'SPtd'//diagSufx
         doDiagSPtend = DIAGNOSTICS_IS_ON(diagName,myThid)
        ENDIF
#endif /* ALLOW_DIAGNOSTICS */

C-- initializing:
        Nrp1=Nr+1
        DO k=1,Nr
         DO j=1,OLy
          DO i=1,OLx
           plumetend(i,j,k) = 0. _d 0
           work(i,j,k) = tracer(i,j,k,bi,bj)
          ENDDO
         ENDDO
        ENDDO
        DO j=1,OLy
         DO i=1,OLx
           work(i,j,Nrp1) = 0. _d 0
         ENDDO
        ENDDO
C-----------------

Catn: After discussion with JM, it is cleaner to remove the negative
c     buoyancy associated with saltplumeflux from the surface lev
c     here instead of inside salt_plume_forcing_surf.F and
c     kpp_forcing_surf.F:

        IF(trApplyFlag.LT.1) THEN
C ======= if trApplyFlag==0: calc SPforcing[T,S] =================
         IF(trIdentity.EQ.GAD_TEMPERATURE) THEN
          DO k=1,Nr
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             SPforcingT(i,j,k,bi,bj)=0. _d 0
            ENDDO
           ENDDO
          ENDDO
         ENDIF
         IF(trIdentity.EQ.GAD_SALINITY) THEN
          DO k=1,Nr
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             SPforcingS(i,j,k,bi,bj)=0. _d 0
            ENDDO
           ENDDO
          ENDDO
         ENDIF

         DO k=Nr,1,-1
          kp1=k+1
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
C           IF(trIdentity.EQ.GAD_SALINITY) SurfVal=SPbrineSalt(i,j,bi,bj)
           IF(trIdentity.EQ.GAD_TEMPERATURE) SurfVal=tracer(i,j,1,bi,bj)
Catn: m/s*[degC,psu]
            gTr_totSurf2Below = SPbrineVolFlux(i,j,bi,bj)*SurfVal
            gTr_Surf2kLev = dSPvolSurf2kLev(i,j,k,bi,bj) * SurfVal
            dSPvolBelow2kLev = -dSPvolkLev2Above(i,j,kp1,bi,bj)
            gTr_Below2kLev= dSPvolBelow2kLev * work(i,j,kp1)
Catn: gTr_kLev2Above works even for kLev=1 because this is how much
C volume of original [salinity,heat] associated with [SSS,SST]
C was replaced by same volume of brine [salt,heat(from SEAICE_Tfrz)].
C Note: by design, dSPvolkLev2Above already is negative
            gTr_kLev2Above= dSPvolkLev2Above(i,j,k,bi,bj) * work(i,j,k)

C salt: [m/s * psu * kg/m3] = [kg/s/m2 psu] = unit of saltPlumeFlux
C theta:[m/s * kg/m3 * J/kg/degC * degC] = [W/m2]
            plumetend(i,j,k) = ConvertFac *
     &          ( gTr_Surf2kLev + gTr_Below2kLev + gTr_kLev2Above )
            IF(k.eq.1) plumetend(i,j,k) = ConvertFac *
     &          ( gTr_Surf2kLev + gTr_Below2kLev)
C remove negative buoyancy from surface,
C used to be in salt_plume_forcing_surf.F
C           IF(k.EQ.1) THEN
C            plumetend(i,j,k) = plumetend(i,j,k)
C     &                         - ConvertFac*gTr_totSurf2Below
C           ENDIF

Catn: report T/S SPforcing[T,S] related to saltplumeflux for kpp
C     and return zero to do_oceanic_phys.F; unit: [g/m2/s or W/m2]
C            trStar(i,j,k,bi,bj) = 0. _d 0
            IF (trIdentity.EQ.GAD_SALINITY) THEN
             SPforcingS(i,j,k,bi,bj)=plumetend(i,j,k)
            ENDIF
            IF (trIdentity.EQ.GAD_TEMPERATURE) THEN
             SPforcingT(i,j,k,bi,bj)=plumetend(i,j,k)
            ENDIF
           ENDDO
          ENDDO
         ENDDO
#ifdef ALLOW_DIAGNOSTICS
        IF ( doDiagSPtend )
     &   CALL DIAGNOSTICS_FILL(plumetend, diagName, 0,Nr,2,bi,bj,myThid)
#endif /* ALLOW_DIAGNOSTICS */

        ELSE
C ============ updating gTr =====================================
Catn: updating tendency gTr (gS,gT); unit: [psu/s or degC/s]
C         WRITE(msgBuf,'(a,2e24.17)') 'b4SPap: [Tr,gTr](1,1,1,bi,bj):',
C     &     tracer(1,1,1,bi,bj),trStar(1,1,1,bi,bj)
C         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
C     &                      SQUEEZE_RIGHT, myThid )
         DO k=Nr,1,-1
          kp1=k+1
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            IF (trIdentity.EQ.GAD_TEMPERATURE)
     &       plumetend(i,j,k)=SPforcingT(i,j,k,bi,bj)
            IF (trIdentity.EQ.GAD_SALINITY) THEN
             plumetend(i,j,k)=SPforcingS(i,j,k,bi,bj)
            ENDIF

C            trStar(i,j,k,bi,bj)=trStar(i,j,k,bi,bj)+plumetend(i,j,k)
C     &      *recip_drF(k)*recip_hFac_arg(i,j,k)
C     &      *recip_ConvertFac
           ENDDO
          ENDDO
         ENDDO

        ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C--    end if on-off-flag
       ENDIF
C--   end trApplyFlag
      ENDIF

#endif /* SALT_PLUME_VOLUME */
#endif /* ALLOW_SALT_PLUME */

      RETURN
      END