atnguyen/code_21Dec2012_saltplume/salt_plume_volfrac.F

C $Header: /u/gcmpack/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/salt_plume_volfrac.F,v 1.4 2014/04/29 06:49:40 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     SPkBottom    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     dplumek
      _RL locz, zo, z20
      INTEGER i,j,k,kk,Nlev,Nr,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
         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
        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
        SPbrineVolFlux(i,j,bi,bj)        = 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
        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*SPbrineVolFlux(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)
        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

C#ifdef ALLOW_DIAGNOSTICS
C      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 )
C       CALL DIAGNOSTICS_FILL(
C     &      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 )
C      ENDIF
C#endif /* ALLOW_DIAGNOSTICS */

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

      RETURN
      END
1	atn	1.5	C $Header: /u/gcmpack/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/salt_plume_volfrac.F,v 1.4 2014/04/29 06:49:40 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.5	_RL SPkBottom (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66	atn	1.1	_RL dplumek
67			_RL locz, zo, z20
68	atn	1.5	INTEGER i,j,k,kk,Nlev,Nr,Nrp1
69			_RL three, S, So
70			parameter( three = 3. _d 0 )
71	atn	1.1	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.5	Nrp1=Nr+1
77	atn	1.1	dplumek = 0. _d 0
78			DO k=1,Nr
79			DO j=1-OLy,sNy+OLy
80			DO i=1-OLx,sNx+OLx
81	atn	1.2	dSPvolSurf2kLev(i,j,k,bi,bj) = 0. _d 0
82			dSPvolBelow2kLev(i,j,k,bi,bj) = 0. _d 0
83			dSPvolkLev2Above(i,j,k,bi,bj) = 0. _d 0
84			SPplumek(i,j,k,bi,bj) = 1. _d 0
85	atn	1.1	ENDDO
86			ENDDO
87			ENDDO
88			DO j=1-OLy,sNy+OLy
89			DO i=1-OLx,sNx+OLx
90	atn	1.5	dSPvolBelow2kLev(i,j,Nrp1,bi,bj) = 0. _d 0
91			SPplumek(i,j,Nrp1,bi,bj) = 1. _d 0
92			C
93			C SPbrineSalt(i,j,bi,bj) = 0. _d 0
94			SPbrineVolFlux(i,j,bi,bj) = 0. _d 0
95			dMbdt(i,j) = 0. _d 0
96			SPkBottom(i,j) = 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.5	& SaltPlumeDepth(i,j,bi,bj).GT.zeroRL ) THEN
107	atn	1.1
108	atn	1.5	SPkBottom(i,j)=min(k,Nr)
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.5	S = oneRL !input depth temp
113			So = oneRL
114	atn	1.1	DO kk=1,Npower+1
115			S = locz*S !raise to the Npower+1
116			So = zo*So
117			ENDDO
118	atn	1.5	SPplumek(i,j,k,bi,bj) = max(zeroRL,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.5	S = exp(locz/z20)-oneRL
123	atn	1.1	So = recip_expOneM1 !So = exp(one)-one
124	atn	1.5	SPplumek(i,j,k,bi,bj) = max(zeroRL,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.5	SPplumek(i,j,k,bi,bj) = max(zeroRL,S/So)
137	atn	1.1	ELSE
138	atn	1.5	SPplumek(i,j,k,bi,bj) = zeroRL
139	atn	1.1	ENDIF
140
141			C PlumeMethod = 5, dumping all salt at the top layer
142			ELSEIF (PlumeMethod .EQ. 5) THEN
143	atn	1.5	z20 = oneRL
144			zo = zeroRL
145	atn	1.1	IF( (locz.GE.zo).AND.(locz.LT.z20) ) THEN
146	atn	1.5	SPplumek(i,j,k,bi,bj) = zeroRL
147	atn	1.1	ELSE
148	atn	1.5	SPplumek(i,j,k,bi,bj) = oneRL
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.5	S = oneRL !input depth temp
154			So = oneRL
155	atn	1.1	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.5	SPplumek(i,j,k,bi,bj) = max(zeroRL,twoRL/z20*locz-S/So)
161	atn	1.1	ELSE !Npower=2
162	atn	1.5	SPplumek(i,j,k,bi,bj) = max(zeroRL,
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.5	zo = z20*(oneRL/twoRL-Npower/200. _d 0)
178	atn	1.1	So = (z20z20z20/three)
179			& - (z20z20)zo
180			& + (z20) zozo
181			S = (loczloczlocz/three)
182			& - (loczlocz)zo
183			& + (locz) zozo
184	atn	1.5	SPplumek(i,j,k,bi,bj) = max(zeroRL,(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.5	SPplumek(i,j,k,bi,bj) = oneRL
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	atn	1.4	C Assuming we go with (B) here:
220	atn	1.1	DO j=1-OLy,sNy+OLy
221			DO i=1-OLx,sNx+OLx
222	atn	1.5	C brine salinity at surface:
223	atn	1.4	C SPbrineSalt(i,j,bi,bj)= saltPlumeFlux(i,j,bi,bj)/SPalpha
224			C & mass2rUnitrecip_drF(1)*dTtracerLev(1)
225	atn	1.5	C SPbrineSalt(i,j,bi,bj) = SPbrineSconst
226			C SPbrineSalt(i,j,bi,bj) = min( SPbrineSalt(i,j,bi,bj)
227			C & ,SPbrineSaltmax )
228			C brine mass and volume at surface:
229			dMbdt(i,j)=saltPlumeFlux(i,j,bi,bj)/SPbrineSconst
230			SPbrineVolFlux(i,j,bi,bj)=dMbdt(i,j)*mass2rUnit
231	atn	1.1	ENDDO
232			ENDDO
233
234			C Distributing down: this is always from level 1 to depth
235			DO k=Nr,1,-1
236			DO j=1-OLy,sNy+OLy
237			DO i=1-OLx,sNx+OLx
238	atn	1.2	dplumek=SPplumek(i,j,k+1,bi,bj)-SPplumek(i,j,k,bi,bj)
239	atn	1.5	dSPvolSurf2kLev(i,j,k,bi,bj)=dplumek*SPbrineVolFlux(i,j)
240	atn	1.1	ENDDO
241			ENDDO
242			ENDDO
243
244			C Now volume up: need to scan from bottom of SPDepth
245			DO j=1-OLy,sNy+OLy
246			DO i=1-OLx,sNx+OLx
247	atn	1.5	Nlev=SPkBottom(i,j)
248			IF(Nlev.GE.1 .AND. Nlev.LE.Nr) THEN
249	atn	1.1	DO k=Nlev,1,-1
250	atn	1.5	kp1=k+1
251			C dSPvolBelow2kLev(i,j,k,bi,bj)=-dSPvolkLev2Above(i,j,k+1,bi,bj)
252			C dSPvolkLev2Above(i,j,k,bi,bj)=-(dSPvolBelow2kLev(i,j,k,bi,bj)
253			C & + dSPvolSurf2kLev(i,j,k,bi,bj))
254			dSPvolkLev2Above(i,j,k,bi,bj)=dSPvolkLev2Above(i,j,kp1,bi,bj)
255			& - dSPvolSurf2kLev(i,j,k,bi,bj)
256	atn	1.1	ENDDO
257			ENDIF
258			ENDDO
259			ENDDO
260	atn	1.5
261	atn	1.4	C#ifdef ALLOW_DIAGNOSTICS
262			C IF ( useDiagnostics ) THEN
263			C CALL DIAGNOSTICS_FILL(
264			C & SPkBottom,'SPkbottm',0,1,1,bi,bj,myThid )
265			C CALL DIAGNOSTICS_FILL(
266			C & SPbrineSalt,'SPbrineS',0,1,1,bi,bj,myThid )
267			C CALL DIAGNOSTICS_FILL(
268			C & SPplumek,'PLUMEKB ',0,Nr,1,bi,bj,myThid )
269			C CALL DIAGNOSTICS_FILL(
270			C & dSPvolSurf2kLev,'SPVs2k ',0,Nr,1,bi,bj,myThid )
271			C CALL DIAGNOSTICS_FILL(
272			C & dSPvolBelow2kLev,'SPVp2k ',0,Nr,1,bi,bj,myThid )
273			C CALL DIAGNOSTICS_FILL(
274			C & dSPvolkLev2Above,'SPVk2m ',0,Nr,1,bi,bj,myThid )
275			C ENDIF
276			C#endif /* ALLOW_DIAGNOSTICS */
277	atn	1.2
278			#endif /* SALT_PLUME_VOLUME */
279	atn	1.1	#endif /* ALLOW_SALT_PLUME */
280
281			RETURN
282			END