C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/Attic/salt_plume_frac.F,v 1.2 2014/04/20 03:49:35 atn dead $
C $Name:  $

#include "SALT_PLUME_OPTIONS.h"

CBOP
C     !ROUTINE: SALT_PLUME_FRAC
C     !INTERFACE:
      SUBROUTINE SALT_PLUME_FRAC(
     I                  imax, fact,SPDepth,
     U                  plumek,
     I                  myTime, myIter, myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE SALT_PLUME_FRAC
C     | o Compute saltplume penetration.
C     *==========================================================*
C     | Compute fraction of saltplume (flux) penetrating to
C     | specified depth, plumek, due to rejected salt
C     | during freezing.
C     | For example, if surface value is Saltplume0,
C     | and each level gets equal fraction 1/5 down to SPDepth=5,
C     | SALT_PLUME_FRAC will report plumek = 4/5 on output if the input
C     | plumek = 1. Else, output plumek = 0.                
C     | Reference : Duffy et al, (GRL 1999)
C     |
C     | =====
C     | Written by   : ATN (based on SWFRAC)
C     | Date         : Sep 13, 2007
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
#include "SIZE.h"
#include "SALT_PLUME.h"
#include "EEPARAMS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     input arguments
C     imax    :: number of vertical grid points
C     fact    :: scale  factor to apply to depth array
C     SPDpeth :: corresponding SaltPlumeDepth(i,j) at this grid point
C     myTime  :: Current time in simulation
C     myIter  :: Current iteration number in simulation
C     myThid  :: My Thread Id. number
      INTEGER imax
      _RL     fact
      _RL     myTime
      INTEGER myIter
      INTEGER myThid
C     input/output arguments
C     SPDEPTH(2):: identical value of salt_plume_depth
C     plumek :: on input: vertical depth for desired plume fraction
C               (fact*plumek) is negative distance (m) from surface
C     plumek :: on output: saltplume contribution fraction
      _RL     plumek(imax), SPDepth(imax)
CEOP

#ifdef ALLOW_SALT_PLUME

C     !LOCAL VARIABLES:
      _RL facz, dd, dd20
      INTEGER i
#ifndef TARGET_NEC_SX
      INTEGER kk
#endif
      _RL     three, S, So, tiny
      parameter( three = 3. _d 0, tiny = 1. _d -5 )
C     This is an abbreviation of 1./(exp(1.)-1.)
      _RL     recip_expOneM1
      parameter( recip_expOneM1 = 0.581976706869326343 )
      
      IF(plumek(1).LE.SPDepth(1) .and. plumek(2).GE.SPDepth(1)) THEN
         plumek(2) = SPDepth(1)
      ENDIF

      DO i = 1,imax
       facz = abs(fact*plumek(i))
       IF (SPDepth(i).GE.facz .and. SPDepth(i) .GT. zeroRL) THEN

C     Default: uniform distribution, PlumeMethod=1, Npower=0
        IF (PlumeMethod .EQ. 1) THEN
         dd20 = (abs(SPDepth(i)))
#ifdef TARGET_NEC_SX
         IF ( dd20 .GT. 0. _d 0 ) THEN
          S   = (facz/dd20)
C     crazy attempt to make the code faster and raise tempN to
C     (Npower+1)
          IF (Npower .GT. 0) S = S*S**Npower
         ELSE
          S = 0. _d 0
         ENDIF
         plumek(i) = max(zeroRL,S)
#else
         S  = oneRL                  !input depth temp
         So = oneRL
         DO kk=1,Npower+1
            S  = facz*S            !raise to the Npower+1
            So = dd20*So
         ENDDO
         plumek(i) = max(zeroRL,S/So)
#endif /* TARGET_NEC_SX */

        ELSEIF (PlumeMethod .EQ. 2) THEN  !exponential distribution
         dd = abs(SPDepth(i))
         So=exp(oneRL)-oneRL
         S =exp(facz/dd)-oneRL
         plumek(i) = max(zeroRL,S/So)

C     PlumeMethod = 3, distribute salt LINEARLY between SPDepth and SPDepth/SPovershoot
C     (1-SPovershoot)% has already been taken into account in SPDepth calculation,
C     i.e., SPDepth = SPovershoot*SPDepth.
        ELSEIF (PlumeMethod .EQ. 3) THEN  !overshoot 20%
         dd20 = (abs(SPDepth(i)))
         dd   = dd20/SPovershoot
         So=dd20-dd
         S =facz-dd
         IF( (facz.GE.dd).AND.(facz.LT.dd20) ) THEN
          plumek(i) = max(zeroRL,S/So)
         ELSE
          plumek(i) = zeroRL
         ENDIF
C     PlumeMethod = 5, dumping all salt at the top layer
        ELSEIF (PlumeMethod .EQ. 5) THEN
         IF( (facz.GE.0.0).AND.(facz.LT.1.0) ) THEN
          plumek(i) = zeroRL 
         ELSE
          plumek(i) = oneRL
         ENDIF
        ELSEIF (PlumeMethod .EQ. 6) THEN
C     PLumeMethod = 6, currently only works for Npower = 1 and 2.
         dd20 = (abs(SPDepth(i)))
#ifdef TARGET_NEC_SX
         IF ( dd20 .GT. 0. _d 0 ) THEN
          S  = (facz/dd20)
C     crazy attempt to make the code faster and raise tempN to
C     (Npower+1)
          IF (Npower .GT. 0) S = S*S**Npower
          So = 1. _d 0/dd20
         ELSE
          S  = 0. _d 0
          So = 0. _d 0
         ENDIF
         IF(Npower.EQ.1) THEN   !Npower=1
          plumek(i) = max(zeroRL,twoRL*So*facz-S)
         ELSE                   !Npower=2
          plumek(i) = max(zeroRL,
     &         three*So*facz - three*So*So*facz*facz
     &         + tempN)
         ENDIF
#else
          S    = oneRL               !input depth temp
          So   = oneRL 
          DO kk=1,Npower+1
             S  = facz*S             !raise to the Npower+1
             So = dd20*So
          ENDDO
          IF(Npower.EQ.1) THEN           !Npower=1
             plumek(i) = max(zeroRL,twoRL/dd20*facz-S/So)
          ELSE                           !Npower=2
             plumek(i) = max(zeroRL,
     &            three/dd20*facz - three/(dd20*dd20)*facz*facz
     &            + S/So)
          ENDIF
#endif /* TARGET_NEC_SX */
          
        ELSEIF (PlumeMethod .EQ. 7) THEN 
C     PLumeMethod = 7, dump an offset parabolla with more salt at surface
C        tapered to zero at depth SPDepth/2, then increased until SPDepth
C        need input SPDepth, NPower = percentage of SPDepth
C        Ex: if Npower = 10 -> (10/2)=5% of SPDepth
C        NPower can be negative integer here.
C        0 -> parabola centered at SPDepth/2;
C        + -> parabola offset, salt @ surface < @ SPDepth
C        - -> parabola offset, salt @ surface > @ SPDepth
C        S and So are dummy variables
         dd   = (abs(SPDepth(i)))
         dd20 = dd*(oneRL/twoRL-Npower/200. _d 0)
         So   = (dd*dd*dd/three)
     &            -(dd*dd)      *dd20
     &            +(dd)         *dd20*dd20
         S    = (facz*facz *facz/three)
     &             - (facz*facz)*dd20
     &             + (facz)     *dd20*dd20
         plumek(i) = max(zeroRL,(S/So))

#ifndef TARGET_NEC_SX
        ELSE
         WRITE(*,*) 'salt_plume_frac: PLumeMethod =', PLumeMethod,
     &        'not implemented'
         STOP 'ABNORMAL in S/R SALT_PLUME_FRAC'
#endif /* not TARGET_NEC_SX */
        ENDIF
       ELSE
        plumek(i) = zeroRL 
       ENDIF
      ENDDO
      
#endif /* ALLOW_SALT_PLUME */

      RETURN
      END