atnguyen/code_21Dec2012_saltplume/salt_plume_frac.F

C $Header: /u/gcmpack/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/salt_plume_frac.F,v 1.1 2012/12/21 10:00:29 atn Exp $
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
1	atn	1.2	C $Header: /u/gcmpack/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/salt_plume_frac.F,v 1.1 2012/12/21 10:00:29 atn Exp $
2	atn	1.1	C $Name: $
3
4			#include "SALT_PLUME_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: SALT_PLUME_FRAC
8			C !INTERFACE:
9			SUBROUTINE SALT_PLUME_FRAC(
10			I imax, fact,SPDepth,
11			U plumek,
12			I myTime, myIter, myThid )
13
14			C !DESCRIPTION: \bv
15			C ==========================================================
16			C \| SUBROUTINE SALT_PLUME_FRAC
17			C \| o Compute saltplume penetration.
18			C ==========================================================
19			C \| Compute fraction of saltplume (flux) penetrating to
20			C \| specified depth, plumek, due to rejected salt
21			C \| during freezing.
22			C \| For example, if surface value is Saltplume0,
23			C \| and each level gets equal fraction 1/5 down to SPDepth=5,
24			C \| SALT_PLUME_FRAC will report plumek = 4/5 on output if the input
25			C \| plumek = 1. Else, output plumek = 0.
26			C \| Reference : Duffy et al, (GRL 1999)
27			C \|
28			C \| =====
29			C \| Written by : ATN (based on SWFRAC)
30			C \| Date : Sep 13, 2007
31			C ==========================================================
32			C \ev
33
34			C !USES:
35			IMPLICIT NONE
36			#include "SIZE.h"
37			#include "SALT_PLUME.h"
38			#include "EEPARAMS.h"
39
40			C !INPUT/OUTPUT PARAMETERS:
41			C input arguments
42			C imax :: number of vertical grid points
43			C fact :: scale factor to apply to depth array
44			C SPDpeth :: corresponding SaltPlumeDepth(i,j) at this grid point
45			C myTime :: Current time in simulation
46			C myIter :: Current iteration number in simulation
47			C myThid :: My Thread Id. number
48			INTEGER imax
49			_RL fact
50			_RL myTime
51			INTEGER myIter
52			INTEGER myThid
53			C input/output arguments
54			C SPDEPTH(2):: identical value of salt_plume_depth
55			C plumek :: on input: vertical depth for desired plume fraction
56			C (fact*plumek) is negative distance (m) from surface
57			C plumek :: on output: saltplume contribution fraction
58			_RL plumek(imax), SPDepth(imax)
59			CEOP
60
61			#ifdef ALLOW_SALT_PLUME
62
63			C !LOCAL VARIABLES:
64			_RL facz, dd, dd20
65			INTEGER i
66			#ifndef TARGET_NEC_SX
67			INTEGER kk
68			#endif
69			_RL three, S, So, tiny
70			parameter( three = 3. _d 0, tiny = 1. _d -5 )
71			C This is an abbreviation of 1./(exp(1.)-1.)
72			_RL recip_expOneM1
73			parameter( recip_expOneM1 = 0.581976706869326343 )
74
75			IF(plumek(1).LE.SPDepth(1) .and. plumek(2).GE.SPDepth(1)) THEN
76			plumek(2) = SPDepth(1)
77			ENDIF
78
79			DO i = 1,imax
80			facz = abs(fact*plumek(i))
81			IF (SPDepth(i).GE.facz .and. SPDepth(i) .GT. zeroRL) THEN
82
83			C Default: uniform distribution, PlumeMethod=1, Npower=0
84			IF (PlumeMethod .EQ. 1) THEN
85			dd20 = (abs(SPDepth(i)))
86			#ifdef TARGET_NEC_SX
87			IF ( dd20 .GT. 0. _d 0 ) THEN
88			S = (facz/dd20)
89			C crazy attempt to make the code faster and raise tempN to
90			C (Npower+1)
91			IF (Npower .GT. 0) S = SS*Npower
92			ELSE
93			S = 0. _d 0
94			ENDIF
95			plumek(i) = max(zeroRL,S)
96			#else
97			S = oneRL !input depth temp
98			So = oneRL
99			DO kk=1,Npower+1
100			S = facz*S !raise to the Npower+1
101			So = dd20*So
102			ENDDO
103			plumek(i) = max(zeroRL,S/So)
104			#endif /* TARGET_NEC_SX */
105
106			ELSEIF (PlumeMethod .EQ. 2) THEN !exponential distribution
107			dd = abs(SPDepth(i))
108			So=exp(oneRL)-oneRL
109			S =exp(facz/dd)-oneRL
110			plumek(i) = max(zeroRL,S/So)
111
112			C PlumeMethod = 3, distribute salt LINEARLY between SPDepth and SPDepth/SPovershoot
113			C (1-SPovershoot)% has already been taken into account in SPDepth calculation,
114			C i.e., SPDepth = SPovershoot*SPDepth.
115			ELSEIF (PlumeMethod .EQ. 3) THEN !overshoot 20%
116			dd20 = (abs(SPDepth(i)))
117			dd = dd20/SPovershoot
118			So=dd20-dd
119			S =facz-dd
120			IF( (facz.GE.dd).AND.(facz.LT.dd20) ) THEN
121			plumek(i) = max(zeroRL,S/So)
122			ELSE
123			plumek(i) = zeroRL
124			ENDIF
125			C PlumeMethod = 5, dumping all salt at the top layer
126			ELSEIF (PlumeMethod .EQ. 5) THEN
127			IF( (facz.GE.0.0).AND.(facz.LT.1.0) ) THEN
128			plumek(i) = zeroRL
129			ELSE
130			plumek(i) = oneRL
131			ENDIF
132			ELSEIF (PlumeMethod .EQ. 6) THEN
133			C PLumeMethod = 6, currently only works for Npower = 1 and 2.
134			dd20 = (abs(SPDepth(i)))
135			#ifdef TARGET_NEC_SX
136			IF ( dd20 .GT. 0. _d 0 ) THEN
137			S = (facz/dd20)
138			C crazy attempt to make the code faster and raise tempN to
139			C (Npower+1)
140			IF (Npower .GT. 0) S = SS*Npower
141			So = 1. _d 0/dd20
142			ELSE
143			S = 0. _d 0
144			So = 0. _d 0
145			ENDIF
146			IF(Npower.EQ.1) THEN !Npower=1
147			plumek(i) = max(zeroRL,twoRLSofacz-S)
148			ELSE !Npower=2
149			plumek(i) = max(zeroRL,
150			& threeSofacz - threeSoSofaczfacz
151			& + tempN)
152			ENDIF
153			#else
154			S = oneRL !input depth temp
155			So = oneRL
156			DO kk=1,Npower+1
157			S = facz*S !raise to the Npower+1
158			So = dd20*So
159			ENDDO
160			IF(Npower.EQ.1) THEN !Npower=1
161			plumek(i) = max(zeroRL,twoRL/dd20*facz-S/So)
162			ELSE !Npower=2
163			plumek(i) = max(zeroRL,
164			& three/dd20facz - three/(dd20dd20)faczfacz
165			& + S/So)
166			ENDIF
167			#endif /* TARGET_NEC_SX */
168
169			ELSEIF (PlumeMethod .EQ. 7) THEN
170			C PLumeMethod = 7, dump an offset parabolla with more salt at surface
171			C tapered to zero at depth SPDepth/2, then increased until SPDepth
172			C need input SPDepth, NPower = percentage of SPDepth
173			C Ex: if Npower = 10 -> (10/2)=5% of SPDepth
174			C NPower can be negative integer here.
175			C 0 -> parabola centered at SPDepth/2;
176			C + -> parabola offset, salt @ surface < @ SPDepth
177			C - -> parabola offset, salt @ surface > @ SPDepth
178			C S and So are dummy variables
179			dd = (abs(SPDepth(i)))
180			dd20 = dd*(oneRL/twoRL-Npower/200. _d 0)
181			So = (dddddd/three)
182			& -(dddd) dd20
183			& +(dd) dd20dd20
184			S = (faczfacz facz/three)
185			& - (faczfacz)dd20
186			& + (facz) dd20dd20
187			plumek(i) = max(zeroRL,(S/So))
188
189			#ifndef TARGET_NEC_SX
190			ELSE
191			WRITE(,) 'salt_plume_frac: PLumeMethod =', PLumeMethod,
192			& 'not implemented'
193			STOP 'ABNORMAL in S/R SALT_PLUME_FRAC'
194			#endif /* not TARGET_NEC_SX */
195			ENDIF
196			ELSE
197			plumek(i) = zeroRL
198			ENDIF
199			ENDDO
200
201			#endif /* ALLOW_SALT_PLUME */
202
203			RETURN
204			END