C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/salt_plume_tendency_apply_s.F,v 1.3 2014/04/20 09:00:36 atn Exp $
C $Name:  $

#include "SALT_PLUME_OPTIONS.h"

CBOP 0
C     !ROUTINE: SALT_PLUME_TENDENCY_APPLY_S
C     !INTERFACE:
      SUBROUTINE SALT_PLUME_TENDENCY_APPLY_S(
     &     iMin, iMax, jMin, jMax,
     &     bi,bj,kLev,myTime,myThid)

C     !DESCRIPTION:
C     Add salt_plume tendency terms to S tendency.
C     Routine works for one level at a time.
C     SaltPlume is the amount of salt rejected by ice while freezing;
C     it is here redistributed to multiple vertical levels as per
C     Duffy et al. (GRL 1999).

C     !INPUT PARAMETERS:
      IMPLICIT NONE
#include "SIZE.h"
#include "GRID.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "SALT_PLUME.h"

C     !INPUT PARAMETERS:
      integer iMin, iMax, jMin, jMax, kLev, bi, bj, myThid
      _RL myTime
CEOP

#ifdef ALLOW_SALT_PLUME

C     !LOCAL VARIABLES:
      integer i, j
      _RL minusone
      parameter(minusone = -1.)
      _RL plumefrac(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
Catn: unit plumeStend [kg/s/m2 psu] = [g/s/m2]; same as saltPlumeFlux
      _RL plumeStend(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#ifdef SALT_PLUME_VOLUME
      _RL gS_Surf2kLev, gS_Below2kLev, gS_kLev2Above
      integer kp1
#endif

#ifdef TARGET_NEC_SX
      integer imt
      parameter( imt=(sNx+2*OLx)*(sNy+2*OLy) )
      _RL plumekb2D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#else
      integer two2
      parameter(two2 = 2)
      _RL plumekb(two2), SPdepth(two2)
#endif

#ifndef SALT_PLUME_VOLUME
      DO j=jMin,jMax
       DO i=iMin,iMax
C Penetrating saltplume fraction:cumulativeSP(klev+1)-cumulativeSP(klev)
        IF ( SaltPlumeDepth(i,j,bi,bj) .GT. abs(rF(kLev)) ) THEN
         plumefrac(I,J) = 0. _d 0
         plumekb(1)=abs(rF(klev))
         plumekb(2)=abs(rF(klev+1))
         SPdepth(1)=SaltPlumeDepth(i,j,bi,bj)
         SPdepth(2)=SaltPlumeDepth(i,j,bi,bj)
         CALL SALT_PLUME_FRAC(
     I                   two2,minusone,SPdepth,
     U                   plumekb,
     I                   myTime, 1, myThid )
         plumefrac(I,J) = (plumekb(2)-plumekb(1))*maskC(i,j,klev,bi,bj)
         plumeStend(I,J) = saltPlumeFlux(i,j,bi,bj)*plumefrac(I,J)
         gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)+plumeStend(I,J)
     &        *recip_drF(kLev)*mass2rUnit*_recip_hFacC(i,j,kLev,bi,bj)
        ELSE
         plumefrac(I,J) = 0. _d 0
         plumeStend(I,J)= 0. _d 0
        ENDIF
       ENDDO
      ENDDO
#else /* SALT_PLUME_VOLUME */
      kp1 = kLev + 1
      if(kLev.EQ.Nr) kp1 = kLev
      DO j=jMin,jMax
       DO i=iMin,iMax
Catn: [m/s * psu * kg/m3] = [kg/s/m2 psu] = unit of saltPlumeFlux
        gS_Surf2kLev  = dSPvolSurf2kLev(i,j,kLev,bi,bj)
     &                * SPbrineSalt(i,j,bi,bj) * rhoConst
        gS_Below2kLev = dSPvolBelow2kLev(i,j,kLev,bi,bj)
                      * salt(i,j,kp1,bi,bj) * rhoConst
Catn: gS_kLev2Above works even for kLev=1 because this is how much
C volume of original salinity was replaced by same volume of brine.
C Note: by design, dSPvolkLev2Above already is negative
        gS_kLev2Above = dSPvolkLev2Above(i,j,kLev,bi,bj)
                      * salt(i,j,kLev,bi,bj) * rhoConst
        plumeStend(i,j)=gS_Surf2kLev + gS_Below2kLev + gS_kLev2Above
        gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)+plumeStend(I,J)
     &        *recip_drF(kLev)*mass2rUnit*_recip_hFacC(i,j,kLev,bi,bj)
       ENDDO
      ENDDO
#endif /* SALT_PLUME_VOLUME */

#ifdef ALLOW_DIAGNOSTICS
#ifndef SALT_PLUME_VOLUME
      IF ( useDiagnostics ) THEN
       CALL DIAGNOSTICS_FILL (
     &      plumefrac,'PLUMEKB ',kLev,1,2,bi,bj,myThid )
#endif /* SALT_PLUME_VOLUME */
       CALL DIAGNOSTICS_FILL (
     &      plumeStend,'oceSPtnd',kLev,1,2,bi,bj,myThid )
      ENDIF
#endif /* ALLOW_DIAGNOSTICS */
      
#endif /* ALLOW_SALT_PLUME */
      
      RETURN
      END