atnguyen/code_21Dec2012_saltplume/salt_plume_volfrac.F

C $Header: /u/gcmpack/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/salt_plume_volfrac.F,v 1.5 2014/05/01 08:09:30 atn Exp $
C $Name:  $

#include "SALT_PLUME_OPTIONS.h"

CBOP
C     !ROUTINE: SALT_PLUME_VOLFRAC
C     !INTERFACE:
      SUBROUTINE SALT_PLUME_VOLFRAC(
     I                       bi, bj, myTime, myIter, myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE SALT_PLUME_VOLFRAC
C     | o Compute saltplume penetration.
C     *==========================================================*
C     | Compute fraction of volume flux associated with saltplume 
C     | flux penetrating through the entire water columns due to 
C     | rejected salt during freezing.
C     | 
C     | For example, if surface value is Saltplume0,
C     | and each level gets equal fraction 1/5 down to SPDepth=5,
C     | SALT_PLUME_VOLFRAC will report 
C     | dSPvolkLev2Above[2to1,3to2,4to3,5to4,6to5] = [4/5,3/5,2/5,1/5,  0]
C     | dSPvolSurf2kLev [1to1,1to2,1to3,1to4,1to5] = [1/5,1/5,1/5,1/5,1/5]
C     | sum [into5] = 1to5 + 6to5 - 5to4 = 1/5 +   0 - 1/5 = 0
C     |     [into4] = 1to4 + 5to4 - 4to3 = 1/5 + 1/5 - 2/5 = 0
C     |     [into3] = 1to3 + 4to3 - 3to2 = 1/5 + 2/5 - 3/5 = 0
C     |     [into2] = 1to2 + 3to2 - 2to1 = 1/5 + 3/5 - 4/5 = 0
C     |     [into1] = 1to1 + 2to1 - 1to[1,2,3,4,5] = 1/5 + 4/5 - 5/5 = 0
C     | NOTE: volume will always be conserved.
C     | =====
C     | Written by   : ATN (based on SALT_PLUME_FRAC)
C     | Date         : Apr 14, 2014
C     *==========================================================*
C     \ev

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

C     !INPUT/OUTPUT PARAMETERS:
C     input arguments
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 bi,bj
      _RL     myTime
      INTEGER myIter
      INTEGER myThid
C     input/output arguments
C      CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP

#ifdef ALLOW_SALT_PLUME
#ifdef SALT_PLUME_VOLUME

C     !LOCAL VARIABLES:
      _RL     dMbdt        (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     dVbdt        (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     SPkBottom    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     dplumek
      _RL locz, zo, z20
      INTEGER i,j,k,kp1,kk,Nlev,Nrp1
      _RL     three, S, So
      parameter( three = 3. _d 0 )
C     This is an abbreviation of 1./(exp(1.)-1.)
      _RL     recip_expOneM1
      parameter( recip_expOneM1 = 0.581976706869326343 )

C initialize at every time step
      Nrp1=Nr+1
      dplumek = 0. _d 0
      DO k=1,Nr
       DO j=1-OLy,sNy+OLy
        DO i=1-OLx,sNx+OLx
         dSPvolSurf2kLev(i,j,k,bi,bj)  = 0. _d 0
C         dSPvolBelow2kLev(i,j,k,bi,bj) = 0. _d 0
         dSPvolkLev2Above(i,j,k,bi,bj) = 0. _d 0
         SPplumek(i,j,k,bi,bj)         = 1. _d 0
        ENDDO
       ENDDO
      ENDDO
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
C        dSPvolBelow2kLev(i,j,Nrp1,bi,bj) = 0. _d 0
        SPplumek(i,j,Nrp1,bi,bj)         = 1. _d 0
C
C       SPbrineSalt(i,j,bi,bj)           = 0. _d 0
C        SPbrineVolFlux(i,j,bi,bj)        = 0. _d 0
        dVbdt(i,j)                       = 0. _d 0
        dMbdt(i,j)                       = 0. _d 0
        SPkBottom(i,j)                   = 0
       ENDDO
      ENDDO

      DO k = 1,Nr+1
       DO j=1-OLy,sNy+OLy
        DO i=1-OLx,sNx+OLx
         locz = abs(rF(k))

         IF ( SaltPlumeDepth(i,j,bi,bj).GT.locz .and.
     &        SaltPlumeDepth(i,j,bi,bj).GT.zeroRL ) THEN

          SPkBottom(i,j)=min(k,Nr)
C     Default: uniform distribution, PlumeMethod=1, Npower=0
          IF (PlumeMethod .EQ. 1) THEN
           zo = abs(SaltPlumeDepth(i,j,bi,bj))
           S  = oneRL                !input depth temp
           So = oneRL
           DO kk=1,Npower+1
            S  = locz*S              !raise to the Npower+1
            So = zo*So
           ENDDO
           SPplumek(i,j,k,bi,bj) = max(zeroRL,S/So)

          ELSEIF (PlumeMethod .EQ. 2) THEN !exponential distribution
           z20 = abs(SaltPlumeDepth(i,j,bi,bj))
           S   = exp(locz/z20)-oneRL
           So  = recip_expOneM1      !So = exp(one)-one
           SPplumek(i,j,k,bi,bj) = max(zeroRL,S*So)

C     PlumeMethod = 3, distribute salt LINEARLY between SPDepth and 
C     SPDepth/SPovershoot
C     (1-SPovershoot)percent has already been taken into account in
C     SPDepth calculation, i.e., SPDepth = SPovershoot*SPDepth.
          ELSEIF (PlumeMethod .EQ. 3) THEN !overshoot 20%
           z20 = abs(SaltPlumeDepth(i,j,bi,bj))
           zo  = z20/SPovershoot
           So  = z20-zo
           S   = locz-zo
           IF( (locz.GE.zo).AND.(locz.LT.z20) ) THEN
            SPplumek(i,j,k,bi,bj) = max(zeroRL,S/So)
           ELSE
            SPplumek(i,j,k,bi,bj) = zeroRL
           ENDIF
 
C     PlumeMethod = 5, dumping all salt at the top layer
          ELSEIF (PlumeMethod .EQ. 5) THEN
           z20  = oneRL
           zo   = zeroRL
           IF( (locz.GE.zo).AND.(locz.LT.z20) ) THEN
            SPplumek(i,j,k,bi,bj) = zeroRL 
           ELSE
            SPplumek(i,j,k,bi,bj) = oneRL
           ENDIF
          ELSEIF (PlumeMethod .EQ. 6) THEN
C     PLumeMethod = 6, currently only works for Npower = 1 and 2.
           z20 = abs(SaltPlumeDepth(i,j,bi,bj))
           S   = oneRL               !input depth temp
           So  = oneRL
           DO kk=1,Npower+1
            S  = locz*S              !raise to the Npower+1
            So = z20*So
           ENDDO
           IF(Npower.EQ.1) THEN      !Npower=1
            SPplumek(i,j,k,bi,bj) = max(zeroRL,twoRL/z20*locz-S/So)
           ELSE                      !Npower=2
            SPplumek(i,j,k,bi,bj) = max(zeroRL,
     &           three/z20*locz - three/(z20*z20)*locz*locz + S/So)
           ENDIF
           
          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
           z20 = abs(SaltPlumeDepth(i,j,bi,bj))
           zo  = z20*(oneRL/twoRL-Npower/200. _d 0)
           So  = (z20*z20*z20/three)
     &              - (z20*z20)*zo
     &              + (z20)    *zo*zo
           S    = (locz*locz*locz/three)
     &               - (locz*locz)*zo
     &               + (locz)     *zo*zo
           SPplumek(i,j,k,bi,bj) = max(zeroRL,(S/So))

          ELSE
           WRITE(*,*) 'salt_plume_frac: PLumeMethod =', PLumeMethod,
     &          'not implemented'
           STOP 'ABNORMAL in S/R SALT_PLUME_FRAC'
          ENDIF
         ELSE
          SPplumek(i,j,k,bi,bj) = oneRL
         ENDIF
        ENDDO
       ENDDO
      ENDDO

C Now calculating dplumek, dSPvolumeUp, dSPvolSurf2kLev
C units:
C Sbrine=dsb/dt*dt/(rhoConst*SPalpha*drF)[psu kg/m2/s*s/(kg/m3*m)]=[psu]
C SPplumek : fraction : unitless
C SaltPlumeFlux: dsb/dt [psu.kg/m^2/s = g/m^2/s]
C brine_mass_flux dMb/dt = dsb/dt / Sbrine [kg/m2/s]
C                        = dsb/dt / (dsb/dt*dt/(rhoConst*SPalpha*drF))
C                        = rhoConst*SPalpha*drF/dt [kg/m3 m/s]=[kg/m2/s]
C dVbrine/dt = dMb/dt 1/rhoConst [m/s]

C has 2 ways to define brine properties: either provide 
C (A) SPalpha: vol frac or (B) SPbrineSalt: brine salinity.
C (A) SPalpha:  can calc SPbrineSalt as fxn of dhice/dt, 
C     constrained by SPbrineSaltmax:
C     SPbrineSalt=SaltPlumeFlux/rhoConst/SPalpha/drF(1)*dt
C     SPbrineSalt=min(SPbrineSalt,SPbrineSaltmax)
C     dMbdt = saltPlumeFlux / SPbrineSalt
C           = rhoConst*SPalpha*drF(1)/dt <-- a function of SPalpha
C (B) SPbrinesalt provided 
C     dMbdt = saltPlumeFlux / SPbrineSalt <-- fxn of SPbrineSalt

C Assuming we go with (B) here:
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
C brine salinity at surface:
C        SPbrineSalt(i,j,bi,bj)= saltPlumeFlux(i,j,bi,bj)/SPalpha
C     &                         *mass2rUnit*recip_drF(1)*dTtracerLev(1)
C       SPbrineSalt(i,j,bi,bj) = SPbrineSconst
C        SPbrineSalt(i,j,bi,bj) = min( SPbrineSalt(i,j,bi,bj)
C     &                               ,SPbrineSaltmax )
C brine mass and volume at surface:
        dMbdt(i,j)=saltPlumeFlux(i,j,bi,bj)/SPbrineSconst
        dVbdt(i,j)=dMbdt(i,j)*mass2rUnit
C        SPbrineVolFlux(i,j,bi,bj)=dMbdt(i,j)*mass2rUnit
       ENDDO
      ENDDO

C Distributing down: this is always from level 1 to depth
      DO k=Nr,1,-1
       DO j=1-OLy,sNy+OLy
        DO i=1-OLx,sNx+OLx
         dplumek=SPplumek(i,j,k+1,bi,bj)-SPplumek(i,j,k,bi,bj)
         dSPvolSurf2kLev(i,j,k,bi,bj)=dplumek*dVbdt(i,j)
C         dSPvolSurf2kLev(i,j,k,bi,bj)=dplumek*SPbrineVolFlux(i,j,bi,bj)
        ENDDO
       ENDDO
      ENDDO

C Now volume up: need to scan from bottom of SPDepth
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
        Nlev=SPkBottom(i,j)
        IF(Nlev.GE.1 .AND. Nlev.LE.Nr) THEN
         DO k=Nlev,1,-1
          kp1=k+1
C          dSPvolBelow2kLev(i,j,k,bi,bj)=-dSPvolkLev2Above(i,j,k+1,bi,bj)
C          dSPvolkLev2Above(i,j,k,bi,bj)=-(dSPvolBelow2kLev(i,j,k,bi,bj)
C     &                                  + dSPvolSurf2kLev(i,j,k,bi,bj))
          dSPvolkLev2Above(i,j,k,bi,bj)=dSPvolkLev2Above(i,j,kp1,bi,bj)
     &                                  - dSPvolSurf2kLev(i,j,k,bi,bj)
         ENDDO
        ENDIF
       ENDDO
      ENDDO

#ifdef ALLOW_DIAGNOSTICS
      IF ( useDiagnostics ) THEN
C       CALL DIAGNOSTICS_FILL(
C     &      SPkBottom,'SPkbottm',0,1,1,bi,bj,myThid )
C       CALL DIAGNOSTICS_FILL(
C     &      SPbrineSalt,'SPbrineS',0,1,1,bi,bj,myThid )
       CALL DIAGNOSTICS_FILL(
     &      SPplumek,'PLUMEKB ',0,Nr,1,bi,bj,myThid )
C       CALL DIAGNOSTICS_FILL(
C     &      dSPvolSurf2kLev,'SPVs2k  ',0,Nr,1,bi,bj,myThid )
C       CALL DIAGNOSTICS_FILL(
C     &      dSPvolBelow2kLev,'SPVp2k  ',0,Nr,1,bi,bj,myThid )
C       CALL DIAGNOSTICS_FILL(
C     &      dSPvolkLev2Above,'SPVk2m  ',0,Nr,1,bi,bj,myThid )
      ENDIF
#endif /* ALLOW_DIAGNOSTICS */

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

      RETURN
      END
1	atn	1.6	C $Header: /u/gcmpack/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/salt_plume_volfrac.F,v 1.5 2014/05/01 08:09:30 atn Exp $
2	atn	1.1	C $Name: $
3
4			#include "SALT_PLUME_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: SALT_PLUME_VOLFRAC
8			C !INTERFACE:
9			SUBROUTINE SALT_PLUME_VOLFRAC(
10	atn	1.2	I bi, bj, myTime, myIter, myThid )
11	atn	1.1
12			C !DESCRIPTION: \bv
13			C ==========================================================
14			C \| SUBROUTINE SALT_PLUME_VOLFRAC
15			C \| o Compute saltplume penetration.
16			C ==========================================================
17			C \| Compute fraction of volume flux associated with saltplume
18			C \| flux penetrating through the entire water columns due to
19			C \| rejected salt during freezing.
20			C \|
21			C \| For example, if surface value is Saltplume0,
22			C \| and each level gets equal fraction 1/5 down to SPDepth=5,
23			C \| SALT_PLUME_VOLFRAC will report
24			C \| dSPvolkLev2Above[2to1,3to2,4to3,5to4,6to5] = [4/5,3/5,2/5,1/5, 0]
25			C \| dSPvolSurf2kLev [1to1,1to2,1to3,1to4,1to5] = [1/5,1/5,1/5,1/5,1/5]
26			C \| sum [into5] = 1to5 + 6to5 - 5to4 = 1/5 + 0 - 1/5 = 0
27			C \| [into4] = 1to4 + 5to4 - 4to3 = 1/5 + 1/5 - 2/5 = 0
28			C \| [into3] = 1to3 + 4to3 - 3to2 = 1/5 + 2/5 - 3/5 = 0
29			C \| [into2] = 1to2 + 3to2 - 2to1 = 1/5 + 3/5 - 4/5 = 0
30			C \| [into1] = 1to1 + 2to1 - 1to[1,2,3,4,5] = 1/5 + 4/5 - 5/5 = 0
31			C \| NOTE: volume will always be conserved.
32			C \| =====
33			C \| Written by : ATN (based on SALT_PLUME_FRAC)
34			C \| Date : Apr 14, 2014
35			C ==========================================================
36			C \ev
37
38			C !USES:
39			IMPLICIT NONE
40			#include "SIZE.h"
41			#include "GRID.h"
42			#include "SALT_PLUME.h"
43	atn	1.3	#include "EEPARAMS.h"
44			#include "PARAMS.h"
45	atn	1.1
46			C !INPUT/OUTPUT PARAMETERS:
47			C input arguments
48			C SPDpeth :: corresponding SaltPlumeDepth(i,j) at this grid point
49			C myTime :: Current time in simulation
50			C myIter :: Current iteration number in simulation
51			C myThid :: My Thread Id. number
52	atn	1.3	INTEGER bi,bj
53	atn	1.1	_RL myTime
54			INTEGER myIter
55			INTEGER myThid
56			C input/output arguments
57	atn	1.4	C CHARACTER*(MAX_LEN_MBUF) msgBuf
58	atn	1.1	CEOP
59
60			#ifdef ALLOW_SALT_PLUME
61	atn	1.2	#ifdef SALT_PLUME_VOLUME
62	atn	1.1
63			C !LOCAL VARIABLES:
64	atn	1.2	_RL dMbdt (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65	atn	1.6	_RL dVbdt (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66	atn	1.5	_RL SPkBottom (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67	atn	1.1	_RL dplumek
68			_RL locz, zo, z20
69	atn	1.6	INTEGER i,j,k,kp1,kk,Nlev,Nrp1
70	atn	1.5	_RL three, S, So
71			parameter( three = 3. _d 0 )
72	atn	1.1	C This is an abbreviation of 1./(exp(1.)-1.)
73			_RL recip_expOneM1
74			parameter( recip_expOneM1 = 0.581976706869326343 )
75
76	atn	1.2	C initialize at every time step
77	atn	1.5	Nrp1=Nr+1
78	atn	1.1	dplumek = 0. _d 0
79			DO k=1,Nr
80			DO j=1-OLy,sNy+OLy
81			DO i=1-OLx,sNx+OLx
82	atn	1.2	dSPvolSurf2kLev(i,j,k,bi,bj) = 0. _d 0
83	atn	1.6	C dSPvolBelow2kLev(i,j,k,bi,bj) = 0. _d 0
84	atn	1.2	dSPvolkLev2Above(i,j,k,bi,bj) = 0. _d 0
85			SPplumek(i,j,k,bi,bj) = 1. _d 0
86	atn	1.1	ENDDO
87			ENDDO
88			ENDDO
89			DO j=1-OLy,sNy+OLy
90			DO i=1-OLx,sNx+OLx
91	atn	1.6	C dSPvolBelow2kLev(i,j,Nrp1,bi,bj) = 0. _d 0
92	atn	1.5	SPplumek(i,j,Nrp1,bi,bj) = 1. _d 0
93			C
94			C SPbrineSalt(i,j,bi,bj) = 0. _d 0
95	atn	1.6	C SPbrineVolFlux(i,j,bi,bj) = 0. _d 0
96			dVbdt(i,j) = 0. _d 0
97	atn	1.5	dMbdt(i,j) = 0. _d 0
98			SPkBottom(i,j) = 0
99	atn	1.1	ENDDO
100			ENDDO
101
102			DO k = 1,Nr+1
103			DO j=1-OLy,sNy+OLy
104			DO i=1-OLx,sNx+OLx
105	atn	1.3	locz = abs(rF(k))
106	atn	1.1
107	atn	1.3	IF ( SaltPlumeDepth(i,j,bi,bj).GT.locz .and.
108	atn	1.5	& SaltPlumeDepth(i,j,bi,bj).GT.zeroRL ) THEN
109	atn	1.1
110	atn	1.5	SPkBottom(i,j)=min(k,Nr)
111	atn	1.1	C Default: uniform distribution, PlumeMethod=1, Npower=0
112			IF (PlumeMethod .EQ. 1) THEN
113	atn	1.2	zo = abs(SaltPlumeDepth(i,j,bi,bj))
114	atn	1.5	S = oneRL !input depth temp
115			So = oneRL
116	atn	1.1	DO kk=1,Npower+1
117			S = locz*S !raise to the Npower+1
118			So = zo*So
119			ENDDO
120	atn	1.5	SPplumek(i,j,k,bi,bj) = max(zeroRL,S/So)
121	atn	1.1
122			ELSEIF (PlumeMethod .EQ. 2) THEN !exponential distribution
123	atn	1.2	z20 = abs(SaltPlumeDepth(i,j,bi,bj))
124	atn	1.5	S = exp(locz/z20)-oneRL
125	atn	1.1	So = recip_expOneM1 !So = exp(one)-one
126	atn	1.5	SPplumek(i,j,k,bi,bj) = max(zeroRL,S*So)
127	atn	1.1
128			C PlumeMethod = 3, distribute salt LINEARLY between SPDepth and
129			C SPDepth/SPovershoot
130			C (1-SPovershoot)percent has already been taken into account in
131			C SPDepth calculation, i.e., SPDepth = SPovershoot*SPDepth.
132			ELSEIF (PlumeMethod .EQ. 3) THEN !overshoot 20%
133	atn	1.2	z20 = abs(SaltPlumeDepth(i,j,bi,bj))
134	atn	1.1	zo = z20/SPovershoot
135			So = z20-zo
136			S = locz-zo
137			IF( (locz.GE.zo).AND.(locz.LT.z20) ) THEN
138	atn	1.5	SPplumek(i,j,k,bi,bj) = max(zeroRL,S/So)
139	atn	1.1	ELSE
140	atn	1.5	SPplumek(i,j,k,bi,bj) = zeroRL
141	atn	1.1	ENDIF
142
143			C PlumeMethod = 5, dumping all salt at the top layer
144			ELSEIF (PlumeMethod .EQ. 5) THEN
145	atn	1.5	z20 = oneRL
146			zo = zeroRL
147	atn	1.1	IF( (locz.GE.zo).AND.(locz.LT.z20) ) THEN
148	atn	1.5	SPplumek(i,j,k,bi,bj) = zeroRL
149	atn	1.1	ELSE
150	atn	1.5	SPplumek(i,j,k,bi,bj) = oneRL
151	atn	1.1	ENDIF
152			ELSEIF (PlumeMethod .EQ. 6) THEN
153			C PLumeMethod = 6, currently only works for Npower = 1 and 2.
154	atn	1.2	z20 = abs(SaltPlumeDepth(i,j,bi,bj))
155	atn	1.5	S = oneRL !input depth temp
156			So = oneRL
157	atn	1.1	DO kk=1,Npower+1
158			S = locz*S !raise to the Npower+1
159			So = z20*So
160			ENDDO
161			IF(Npower.EQ.1) THEN !Npower=1
162	atn	1.5	SPplumek(i,j,k,bi,bj) = max(zeroRL,twoRL/z20*locz-S/So)
163	atn	1.1	ELSE !Npower=2
164	atn	1.5	SPplumek(i,j,k,bi,bj) = max(zeroRL,
165	atn	1.1	& three/z20locz - three/(z20z20)loczlocz + S/So)
166			ENDIF
167
168			ELSEIF (PlumeMethod .EQ. 7) THEN
169			C PLumeMethod = 7, dump an offset parabolla with more salt at surface
170			C tapered to zero at depth SPDepth/2, then increased until SPDepth
171			C need input SPDepth, NPower = percentage of SPDepth
172			C Ex: if Npower = 10 -> (10/2)=5% of SPDepth
173			C NPower can be negative integer here.
174			C 0 -> parabola centered at SPDepth/2;
175			C + -> parabola offset, salt @ surface < @ SPDepth
176			C - -> parabola offset, salt @ surface > @ SPDepth
177			C S and So are dummy variables
178	atn	1.2	z20 = abs(SaltPlumeDepth(i,j,bi,bj))
179	atn	1.5	zo = z20*(oneRL/twoRL-Npower/200. _d 0)
180	atn	1.1	So = (z20z20z20/three)
181			& - (z20z20)zo
182			& + (z20) zozo
183			S = (loczloczlocz/three)
184			& - (loczlocz)zo
185			& + (locz) zozo
186	atn	1.5	SPplumek(i,j,k,bi,bj) = max(zeroRL,(S/So))
187	atn	1.1
188			ELSE
189			WRITE(,) 'salt_plume_frac: PLumeMethod =', PLumeMethod,
190			& 'not implemented'
191			STOP 'ABNORMAL in S/R SALT_PLUME_FRAC'
192			ENDIF
193			ELSE
194	atn	1.5	SPplumek(i,j,k,bi,bj) = oneRL
195	atn	1.1	ENDIF
196			ENDDO
197			ENDDO
198			ENDDO
199
200			C Now calculating dplumek, dSPvolumeUp, dSPvolSurf2kLev
201			C units:
202			C Sbrine=dsb/dtdt/(rhoConstSPalphadrF)[psu kg/m2/ss/(kg/m3*m)]=[psu]
203	atn	1.2	C SPplumek : fraction : unitless
204	atn	1.1	C SaltPlumeFlux: dsb/dt [psu.kg/m^2/s = g/m^2/s]
205			C brine_mass_flux dMb/dt = dsb/dt / Sbrine [kg/m2/s]
206			C = dsb/dt / (dsb/dtdt/(rhoConstSPalpha*drF))
207			C = rhoConstSPalphadrF/dt [kg/m3 m/s]=[kg/m2/s]
208			C dVbrine/dt = dMb/dt 1/rhoConst [m/s]
209
210			C has 2 ways to define brine properties: either provide
211			C (A) SPalpha: vol frac or (B) SPbrineSalt: brine salinity.
212			C (A) SPalpha: can calc SPbrineSalt as fxn of dhice/dt,
213			C constrained by SPbrineSaltmax:
214	atn	1.2	C SPbrineSalt=SaltPlumeFlux/rhoConst/SPalpha/drF(1)*dt
215	atn	1.1	C SPbrineSalt=min(SPbrineSalt,SPbrineSaltmax)
216			C dMbdt = saltPlumeFlux / SPbrineSalt
217	atn	1.2	C = rhoConstSPalphadrF(1)/dt <-- a function of SPalpha
218	atn	1.1	C (B) SPbrinesalt provided
219			C dMbdt = saltPlumeFlux / SPbrineSalt <-- fxn of SPbrineSalt
220
221	atn	1.4	C Assuming we go with (B) here:
222	atn	1.1	DO j=1-OLy,sNy+OLy
223			DO i=1-OLx,sNx+OLx
224	atn	1.5	C brine salinity at surface:
225	atn	1.4	C SPbrineSalt(i,j,bi,bj)= saltPlumeFlux(i,j,bi,bj)/SPalpha
226			C & mass2rUnitrecip_drF(1)*dTtracerLev(1)
227	atn	1.5	C SPbrineSalt(i,j,bi,bj) = SPbrineSconst
228			C SPbrineSalt(i,j,bi,bj) = min( SPbrineSalt(i,j,bi,bj)
229			C & ,SPbrineSaltmax )
230			C brine mass and volume at surface:
231			dMbdt(i,j)=saltPlumeFlux(i,j,bi,bj)/SPbrineSconst
232	atn	1.6	dVbdt(i,j)=dMbdt(i,j)*mass2rUnit
233			C SPbrineVolFlux(i,j,bi,bj)=dMbdt(i,j)*mass2rUnit
234	atn	1.1	ENDDO
235			ENDDO
236
237			C Distributing down: this is always from level 1 to depth
238			DO k=Nr,1,-1
239			DO j=1-OLy,sNy+OLy
240			DO i=1-OLx,sNx+OLx
241	atn	1.2	dplumek=SPplumek(i,j,k+1,bi,bj)-SPplumek(i,j,k,bi,bj)
242	atn	1.6	dSPvolSurf2kLev(i,j,k,bi,bj)=dplumek*dVbdt(i,j)
243			C dSPvolSurf2kLev(i,j,k,bi,bj)=dplumek*SPbrineVolFlux(i,j,bi,bj)
244	atn	1.1	ENDDO
245			ENDDO
246			ENDDO
247
248			C Now volume up: need to scan from bottom of SPDepth
249			DO j=1-OLy,sNy+OLy
250			DO i=1-OLx,sNx+OLx
251	atn	1.5	Nlev=SPkBottom(i,j)
252			IF(Nlev.GE.1 .AND. Nlev.LE.Nr) THEN
253	atn	1.1	DO k=Nlev,1,-1
254	atn	1.5	kp1=k+1
255			C dSPvolBelow2kLev(i,j,k,bi,bj)=-dSPvolkLev2Above(i,j,k+1,bi,bj)
256			C dSPvolkLev2Above(i,j,k,bi,bj)=-(dSPvolBelow2kLev(i,j,k,bi,bj)
257			C & + dSPvolSurf2kLev(i,j,k,bi,bj))
258			dSPvolkLev2Above(i,j,k,bi,bj)=dSPvolkLev2Above(i,j,kp1,bi,bj)
259			& - dSPvolSurf2kLev(i,j,k,bi,bj)
260	atn	1.1	ENDDO
261			ENDIF
262			ENDDO
263			ENDDO
264	atn	1.5
265	atn	1.6	#ifdef ALLOW_DIAGNOSTICS
266			IF ( useDiagnostics ) THEN
267	atn	1.4	C CALL DIAGNOSTICS_FILL(
268			C & SPkBottom,'SPkbottm',0,1,1,bi,bj,myThid )
269			C CALL DIAGNOSTICS_FILL(
270			C & SPbrineSalt,'SPbrineS',0,1,1,bi,bj,myThid )
271	atn	1.6	CALL DIAGNOSTICS_FILL(
272			& SPplumek,'PLUMEKB ',0,Nr,1,bi,bj,myThid )
273	atn	1.4	C CALL DIAGNOSTICS_FILL(
274			C & dSPvolSurf2kLev,'SPVs2k ',0,Nr,1,bi,bj,myThid )
275			C CALL DIAGNOSTICS_FILL(
276			C & dSPvolBelow2kLev,'SPVp2k ',0,Nr,1,bi,bj,myThid )
277			C CALL DIAGNOSTICS_FILL(
278			C & dSPvolkLev2Above,'SPVk2m ',0,Nr,1,bi,bj,myThid )
279	atn	1.6	ENDIF
280			#endif /* ALLOW_DIAGNOSTICS */
281	atn	1.2
282			#endif /* SALT_PLUME_VOLUME */
283	atn	1.1	#endif /* ALLOW_SALT_PLUME */
284
285			RETURN
286			END