C $Header: /u/gcmpack/MITgcm/pkg/salt_plume/salt_plume_frac.F,v 1.10 2015/02/21 20:49:45 jmc Exp $
C $Name: $
#include "SALT_PLUME_OPTIONS.h"
CBOP
C !ROUTINE: SALT_PLUME_FRAC
C !INTERFACE:
SUBROUTINE SALT_PLUME_FRAC(
I imax, fact,SPDepth,
#ifdef SALT_PLUME_SPLIT_BASIN
I lon,lat,
#endif
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 = 1/5,2/5,3/5,4/5,5/5 on
C | output if input plumek = 1,2,3,4,5. 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 | Modified : Mar 16, 2014 by atn to improve/clean up
C | -> replace 1-[cummulative plumefrac] code which was based
C | on swfrac with cleaner [cummulative plumefrac] on output
C | in order to avoid 1-[1-[cummulative_plumefrac]] whenever
C | SALT_PLUME_FRAC is called and used from outside.
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
#include "SIZE.h"
#include "SALT_PLUME.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 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)
#ifdef SALT_PLUME_SPLIT_BASIN
_RL lon(imax), lat(imax)
#endif
CEOP
#ifdef ALLOW_SALT_PLUME
C !LOCAL VARIABLES:
_RL facz, dd, dd20
INTEGER i, Npowerloc
#ifndef TARGET_NEC_SX
INTEGER kk
#endif
_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 )
DO i = 1,imax
facz = abs(fact*plumek(i))
#ifdef SALT_PLUME_SPLIT_BASIN
IF(SaltPlumeSplitBasin) THEN
Npowerloc = Npower(2)
IF(lat(imax).LT. 85.0 .AND. lat(imax).GT. 71.0
& .AND. lon(imax) .LT. -90.0) Npowerloc = Npower(1)
ELSE
Npowerloc = Npower(1)
ENDIF
#else
Npowerloc = Npower
#endif
IF (SPDepth(i).GE.facz .and. SPDepth(i) .GT. zero) 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. zero ) THEN
S = (facz/dd20)
C crazy attempt to make the code faster and raise S to (Npower+1)
IF (Npowerloc .GT. 0) S = S*S**Npowerloc
ELSE
S = zero
ENDIF
plumek(i) = max(zero,S)
#else
S = one !input depth temp
So = one
DO kk=1,Npowerloc+1
S = facz*S !raise to the Npowerloc+1
So = dd20*So
ENDDO
plumek(i) = max(zero,S/So)
#endif /* TARGET_NEC_SX */
ELSEIF (PlumeMethod .EQ. 2) THEN !exponential distribution
dd = abs(SPDepth(i))
S = exp(facz/dd)-one
So = recip_expOneM1 !So = exp(one)-one
plumek(i) = 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%
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(zero,S/So)
ELSE
plumek(i) = zero
ENDIF
C PlumeMethod = 5, dumping all salt at the top layer
ELSEIF (PlumeMethod .EQ. 5) THEN
dd = zero
dd20 = one
IF( (facz.GE.dd).AND.(facz.LT.dd20) ) THEN
plumek(i) = zero
ELSE
plumek(i) = one
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. zero ) THEN
S = (facz/dd20)
C crazy attempt to make the code faster and raise S to (Npower+1)
IF (Npowerloc .GT. 0) S = S*S**Npowerloc
So = 1. _d 0/dd20
ELSE
S = zero
So = zero
ENDIF
IF(Npowerloc.EQ.1) THEN !Npower=1
plumek(i) = max(zero,two*So*facz-S)
ELSE !Npower=2
plumek(i) = max(zero,
& three*So*facz - three*So*So*facz*facz + S)
ENDIF
#else
S = one !input depth temp
So = one
DO kk=1,Npowerloc+1
S = facz*S !raise to the Npower+1
So = dd20*So
ENDDO
IF(Npowerloc.EQ.1) THEN !Npower=1
plumek(i) = max(zero,two/dd20*facz-S/So)
ELSE !Npower=2
plumek(i) = max(zero,
& three/dd20*facz - three/(dd20*dd20)*facz*facz + S/So)
ENDIF
#endif /* TARGET_NEC_SX */
catn: 15.Mar.2014
catn: this is a new method by atn. After fixing adjoint compiling error,
catn: will switch this on. Currently commenting out for purpose of
catn: comparing old (1-abc) vs new (abc) way of coding
c ELSEIF (PlumeMethod .EQ. 7) THEN
cC PLumeMethod = 7, dump an offset parabolla with more salt at surface
cC tapered to zero at depth SPDepth/2, then increased until SPDepth
cC need input SPDepth, NPower = percentage of SPDepth
cC Ex: if Npower = 10 -> (10/2)=5% of SPDepth
cC NPower can be negative integer here.
cC 0 -> parabola centered at SPDepth/2;
cC + -> parabola offset, salt @ surface < @ SPDepth
cC - -> parabola offset, salt @ surface > @ SPDepth
cC S and So are dummy variables
c dd = (abs(SPDepth(i)))
c dd20 = dd*(one/two-Npower/200. _d 0)
c So = (dd*dd*dd/three)
c & -(dd*dd) *dd20
c & +(dd) *dd20*dd20
c S = (facz*facz *facz/three)
c & - (facz*facz)*dd20
c & + (facz) *dd20*dd20
c plumek(i) = max(zero,(S/So))
c
ELSE
plumek(i) = one
#ifndef TARGET_NEC_SX
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) = one
ENDIF
ENDDO
#endif /* ALLOW_SALT_PLUME */
RETURN
END