atnguyen/code_21Dec2012_saltplume/salt_plume_volfrac.F

C $Header: /u/gcmpack/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/salt_plume_volfrac.F,v 1.2 2014/04/22 04:32:26 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
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     dplumek
      _RL locz, zo, z20
      INTEGER i,j,k,kk,Nlev
      _RL     one, two, three, S, So, zero
      parameter( one = 1. _d 0, two = 2. _d 0, three = 3. _d 0 )
      parameter( zero = 0. _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
      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
         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
         IF(k.EQ.Nr) THEN
          SPplumek(i,j,k+1,bi,bj)         = 1. _d 0
          dSPvolBelow2kLev(i,j,k+1,bi,bj) = 0. _d 0
         ENDIF
        ENDDO
       ENDDO
      ENDDO
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
        SPkBottom(i,j,bi,bj)   = 0
        SPbrineSalt(i,j,bi,bj) = 0. _d 0
        dMbdt(i,j)             = 0. _d 0
        dVbdt(i,j)             = 0. _d 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.zero ) THEN

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

          ELSEIF (PlumeMethod .EQ. 2) THEN !exponential distribution
           z20 = abs(SaltPlumeDepth(i,j,bi,bj))
           S   = exp(locz/z20)-one
           So  = recip_expOneM1      !So = exp(one)-one
           SPplumek(i,j,k,bi,bj) = max(zero,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(zero,S/So)
           ELSE
            SPplumek(i,j,k,bi,bj) = zero
           ENDIF
 
C     PlumeMethod = 5, dumping all salt at the top layer
          ELSEIF (PlumeMethod .EQ. 5) THEN
           z20  = one
           zo   = zero
           IF( (locz.GE.zo).AND.(locz.LT.z20) ) THEN
            SPplumek(i,j,k,bi,bj) = zero 
           ELSE
            SPplumek(i,j,k,bi,bj) = one
           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   = one                 !input depth temp
           So  = one
           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(zero,two/z20*locz-S/So)
           ELSE                      !Npower=2
            SPplumek(i,j,k,bi,bj) = max(zero,
     &           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*(one/two-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(zero,(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) = one
         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 (A) here:
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
        SPbrineSalt(i,j,bi,bj)= saltPlumeFlux(i,j,bi,bj)/SPalpha
     &                         *mass2rUnit*recip_drF(1)*dTtracerLev(1)
        SPbrineSalt(i,j,bi,bj)= min( SPbrineSalt(i,j,bi,bj)
     &                              ,SPbrineSaltmax )
        dMbdt(i,j)=saltPlumeFlux(i,j,bi,bj)/SPbrineSalt(i,j,bi,bj)
        dVbdt(i,j)=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)
        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,bi,bj)
        IF(Nlev.GE.1) THEN
         DO k=Nlev,1,-1
          dSPvolBelow2kLev(i,j,k,bi,bj)=-dSPvolkLev2Above(i,j,k+1,bi,bj)
          dSPvolkLev2Above(i,j,k,bi,bj)=-(dSPvolBelow2kLev(i,j,k,bi,bj)
     &                                  + dSPvolSurf2kLev(i,j,k,bi,bj))
         ENDDO
        ENDIF
       ENDDO
      ENDDO
      
#ifdef ALLOW_DIAGNOSTICS
      IF ( useDiagnostics ) THEN
       CALL DIAGNOSTICS_FILL(
     &      SPkBottom,'SPkbottm',0,1,1,bi,bj,myThid )
       CALL DIAGNOSTICS_FILL(
     &      SPbrineSalt,'SPbrineS',0,1,1,bi,bj,myThid )
       CALL DIAGNOSTICS_FILL(
     &      SPplumek,'PLUMEKB ',0,Nr,1,bi,bj,myThid )
       CALL DIAGNOSTICS_FILL(
     &      dSPvolSurf2kLev,'SPVs2k  ',0,Nr,1,bi,bj,myThid )
       CALL DIAGNOSTICS_FILL(
     &      dSPvolBelow2kLev,'SPVp2k  ',0,Nr,1,bi,bj,myThid )
       CALL DIAGNOSTICS_FILL(
     &      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.3	C $Header: /u/gcmpack/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/salt_plume_volfrac.F,v 1.2 2014/04/22 04:32:26 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			CEOP
58
59			#ifdef ALLOW_SALT_PLUME
60	atn	1.2	#ifdef SALT_PLUME_VOLUME
61	atn	1.1
62			C !LOCAL VARIABLES:
63	atn	1.2	_RL dMbdt (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64			_RL dVbdt (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65	atn	1.1	_RL dplumek
66			_RL locz, zo, z20
67	atn	1.3	INTEGER i,j,k,kk,Nlev
68	atn	1.1	_RL one, two, three, S, So, zero
69			parameter( one = 1. _d 0, two = 2. _d 0, three = 3. _d 0 )
70			parameter( zero = 0. _d 0 )
71			C This is an abbreviation of 1./(exp(1.)-1.)
72			_RL recip_expOneM1
73			parameter( recip_expOneM1 = 0.581976706869326343 )
74
75	atn	1.2	C initialize at every time step
76	atn	1.1	dplumek = 0. _d 0
77			DO k=1,Nr
78			DO j=1-OLy,sNy+OLy
79			DO i=1-OLx,sNx+OLx
80	atn	1.2	dSPvolSurf2kLev(i,j,k,bi,bj) = 0. _d 0
81			dSPvolBelow2kLev(i,j,k,bi,bj) = 0. _d 0
82			dSPvolkLev2Above(i,j,k,bi,bj) = 0. _d 0
83			SPplumek(i,j,k,bi,bj) = 1. _d 0
84	atn	1.1	IF(k.EQ.Nr) THEN
85	atn	1.2	SPplumek(i,j,k+1,bi,bj) = 1. _d 0
86			dSPvolBelow2kLev(i,j,k+1,bi,bj) = 0. _d 0
87	atn	1.1	ENDIF
88			ENDDO
89			ENDDO
90			ENDDO
91			DO j=1-OLy,sNy+OLy
92			DO i=1-OLx,sNx+OLx
93	atn	1.2	SPkBottom(i,j,bi,bj) = 0
94			SPbrineSalt(i,j,bi,bj) = 0. _d 0
95			dMbdt(i,j) = 0. _d 0
96			dVbdt(i,j) = 0. _d 0
97	atn	1.1	ENDDO
98			ENDDO
99
100			DO k = 1,Nr+1
101			DO j=1-OLy,sNy+OLy
102			DO i=1-OLx,sNx+OLx
103	atn	1.3	locz = abs(rF(k))
104	atn	1.1
105	atn	1.3	IF ( SaltPlumeDepth(i,j,bi,bj).GT.locz .and.
106	atn	1.2	& SaltPlumeDepth(i,j,bi,bj).GT.zero ) THEN
107	atn	1.1
108	atn	1.2	SPkBottom(i,j,bi,bj)=k
109	atn	1.1	C Default: uniform distribution, PlumeMethod=1, Npower=0
110			IF (PlumeMethod .EQ. 1) THEN
111	atn	1.2	zo = abs(SaltPlumeDepth(i,j,bi,bj))
112	atn	1.1	S = one !input depth temp
113			So = one
114			DO kk=1,Npower+1
115			S = locz*S !raise to the Npower+1
116			So = zo*So
117			ENDDO
118	atn	1.2	SPplumek(i,j,k,bi,bj) = max(zero,S/So)
119	atn	1.1
120			ELSEIF (PlumeMethod .EQ. 2) THEN !exponential distribution
121	atn	1.2	z20 = abs(SaltPlumeDepth(i,j,bi,bj))
122	atn	1.1	S = exp(locz/z20)-one
123			So = recip_expOneM1 !So = exp(one)-one
124	atn	1.2	SPplumek(i,j,k,bi,bj) = max(zero,S*So)
125	atn	1.1
126			C PlumeMethod = 3, distribute salt LINEARLY between SPDepth and
127			C SPDepth/SPovershoot
128			C (1-SPovershoot)percent has already been taken into account in
129			C SPDepth calculation, i.e., SPDepth = SPovershoot*SPDepth.
130			ELSEIF (PlumeMethod .EQ. 3) THEN !overshoot 20%
131	atn	1.2	z20 = abs(SaltPlumeDepth(i,j,bi,bj))
132	atn	1.1	zo = z20/SPovershoot
133			So = z20-zo
134			S = locz-zo
135			IF( (locz.GE.zo).AND.(locz.LT.z20) ) THEN
136	atn	1.2	SPplumek(i,j,k,bi,bj) = max(zero,S/So)
137	atn	1.1	ELSE
138	atn	1.2	SPplumek(i,j,k,bi,bj) = zero
139	atn	1.1	ENDIF
140
141			C PlumeMethod = 5, dumping all salt at the top layer
142			ELSEIF (PlumeMethod .EQ. 5) THEN
143			z20 = one
144			zo = zero
145			IF( (locz.GE.zo).AND.(locz.LT.z20) ) THEN
146	atn	1.2	SPplumek(i,j,k,bi,bj) = zero
147	atn	1.1	ELSE
148	atn	1.2	SPplumek(i,j,k,bi,bj) = one
149	atn	1.1	ENDIF
150			ELSEIF (PlumeMethod .EQ. 6) THEN
151			C PLumeMethod = 6, currently only works for Npower = 1 and 2.
152	atn	1.2	z20 = abs(SaltPlumeDepth(i,j,bi,bj))
153	atn	1.1	S = one !input depth temp
154			So = one
155			DO kk=1,Npower+1
156			S = locz*S !raise to the Npower+1
157			So = z20*So
158			ENDDO
159			IF(Npower.EQ.1) THEN !Npower=1
160	atn	1.2	SPplumek(i,j,k,bi,bj) = max(zero,two/z20*locz-S/So)
161	atn	1.1	ELSE !Npower=2
162	atn	1.2	SPplumek(i,j,k,bi,bj) = max(zero,
163	atn	1.1	& three/z20locz - three/(z20z20)loczlocz + S/So)
164			ENDIF
165
166			ELSEIF (PlumeMethod .EQ. 7) THEN
167			C PLumeMethod = 7, dump an offset parabolla with more salt at surface
168			C tapered to zero at depth SPDepth/2, then increased until SPDepth
169			C need input SPDepth, NPower = percentage of SPDepth
170			C Ex: if Npower = 10 -> (10/2)=5% of SPDepth
171			C NPower can be negative integer here.
172			C 0 -> parabola centered at SPDepth/2;
173			C + -> parabola offset, salt @ surface < @ SPDepth
174			C - -> parabola offset, salt @ surface > @ SPDepth
175			C S and So are dummy variables
176	atn	1.2	z20 = abs(SaltPlumeDepth(i,j,bi,bj))
177	atn	1.1	zo = z20*(one/two-Npower/200. _d 0)
178			So = (z20z20z20/three)
179			& - (z20z20)zo
180			& + (z20) zozo
181			S = (loczloczlocz/three)
182			& - (loczlocz)zo
183			& + (locz) zozo
184	atn	1.2	SPplumek(i,j,k,bi,bj) = max(zero,(S/So))
185	atn	1.1
186			ELSE
187			WRITE(,) 'salt_plume_frac: PLumeMethod =', PLumeMethod,
188			& 'not implemented'
189			STOP 'ABNORMAL in S/R SALT_PLUME_FRAC'
190			ENDIF
191			ELSE
192	atn	1.2	SPplumek(i,j,k,bi,bj) = one
193	atn	1.1	ENDIF
194			ENDDO
195			ENDDO
196			ENDDO
197
198			C Now calculating dplumek, dSPvolumeUp, dSPvolSurf2kLev
199			C units:
200			C Sbrine=dsb/dtdt/(rhoConstSPalphadrF)[psu kg/m2/ss/(kg/m3*m)]=[psu]
201	atn	1.2	C SPplumek : fraction : unitless
202	atn	1.1	C SaltPlumeFlux: dsb/dt [psu.kg/m^2/s = g/m^2/s]
203			C brine_mass_flux dMb/dt = dsb/dt / Sbrine [kg/m2/s]
204			C = dsb/dt / (dsb/dtdt/(rhoConstSPalpha*drF))
205			C = rhoConstSPalphadrF/dt [kg/m3 m/s]=[kg/m2/s]
206			C dVbrine/dt = dMb/dt 1/rhoConst [m/s]
207
208			C has 2 ways to define brine properties: either provide
209			C (A) SPalpha: vol frac or (B) SPbrineSalt: brine salinity.
210			C (A) SPalpha: can calc SPbrineSalt as fxn of dhice/dt,
211			C constrained by SPbrineSaltmax:
212	atn	1.2	C SPbrineSalt=SaltPlumeFlux/rhoConst/SPalpha/drF(1)*dt
213	atn	1.1	C SPbrineSalt=min(SPbrineSalt,SPbrineSaltmax)
214			C dMbdt = saltPlumeFlux / SPbrineSalt
215	atn	1.2	C = rhoConstSPalphadrF(1)/dt <-- a function of SPalpha
216	atn	1.1	C (B) SPbrinesalt provided
217			C dMbdt = saltPlumeFlux / SPbrineSalt <-- fxn of SPbrineSalt
218
219			C Assuming we go with (A) here:
220			DO j=1-OLy,sNy+OLy
221			DO i=1-OLx,sNx+OLx
222	atn	1.2	SPbrineSalt(i,j,bi,bj)= saltPlumeFlux(i,j,bi,bj)/SPalpha
223	atn	1.3	& mass2rUnitrecip_drF(1)*dTtracerLev(1)
224	atn	1.2	SPbrineSalt(i,j,bi,bj)= min( SPbrineSalt(i,j,bi,bj)
225			& ,SPbrineSaltmax )
226			dMbdt(i,j)=saltPlumeFlux(i,j,bi,bj)/SPbrineSalt(i,j,bi,bj)
227			dVbdt(i,j)=dMbdt(i,j)*mass2rUnit
228	atn	1.1	ENDDO
229			ENDDO
230
231			C Distributing down: this is always from level 1 to depth
232			DO k=Nr,1,-1
233			DO j=1-OLy,sNy+OLy
234			DO i=1-OLx,sNx+OLx
235	atn	1.2	dplumek=SPplumek(i,j,k+1,bi,bj)-SPplumek(i,j,k,bi,bj)
236			dSPvolSurf2kLev(i,j,k,bi,bj)=dplumek*dVbdt(i,j)
237	atn	1.1	ENDDO
238			ENDDO
239			ENDDO
240
241			C Now volume up: need to scan from bottom of SPDepth
242			DO j=1-OLy,sNy+OLy
243			DO i=1-OLx,sNx+OLx
244	atn	1.2	Nlev=SPkBottom(i,j,bi,bj)
245	atn	1.1	IF(Nlev.GE.1) THEN
246			DO k=Nlev,1,-1
247	atn	1.2	dSPvolBelow2kLev(i,j,k,bi,bj)=-dSPvolkLev2Above(i,j,k+1,bi,bj)
248			dSPvolkLev2Above(i,j,k,bi,bj)=-(dSPvolBelow2kLev(i,j,k,bi,bj)
249			& + dSPvolSurf2kLev(i,j,k,bi,bj))
250	atn	1.1	ENDDO
251			ENDIF
252			ENDDO
253			ENDDO
254
255	atn	1.2	#ifdef ALLOW_DIAGNOSTICS
256			IF ( useDiagnostics ) THEN
257			CALL DIAGNOSTICS_FILL(
258			& SPkBottom,'SPkbottm',0,1,1,bi,bj,myThid )
259			CALL DIAGNOSTICS_FILL(
260			& SPbrineSalt,'SPbrineS',0,1,1,bi,bj,myThid )
261			CALL DIAGNOSTICS_FILL(
262			& SPplumek,'PLUMEKB ',0,Nr,1,bi,bj,myThid )
263			CALL DIAGNOSTICS_FILL(
264			& dSPvolSurf2kLev,'SPVs2k ',0,Nr,1,bi,bj,myThid )
265			CALL DIAGNOSTICS_FILL(
266			& dSPvolBelow2kLev,'SPVp2k ',0,Nr,1,bi,bj,myThid )
267			CALL DIAGNOSTICS_FILL(
268			& dSPvolkLev2Above,'SPVk2m ',0,Nr,1,bi,bj,myThid )
269			ENDIF
270			#endif /* ALLOW_DIAGNOSTICS */
271
272			#endif /* SALT_PLUME_VOLUME */
273	atn	1.1	#endif /* ALLOW_SALT_PLUME */
274
275			RETURN
276			END