atnguyen/code_21Dec2012_saltplume/salt_plume_volfrac.F

C $Header: /u/gcmpack/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/salt_plume_volfrac.F,v 1.1 2014/04/20 04:03:07 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"

C     !INPUT/OUTPUT PARAMETERS:
C     input arguments
C     imax    :: number of vertical grid points
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     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
      _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
         zloc = abs(rF(k))

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