atnguyen/code_21Dec2012_saltplume/salt_plume_apply.F

C $Header: /u/gcmpack/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/salt_plume_apply.F,v 1.2 2014/04/22 10:32:25 atn Exp $
C $Name:  $

#include "SALT_PLUME_OPTIONS.h"

CBOP
C     !ROUTINE: SALT_PLUME_APPLY
C     !INTERFACE:
      SUBROUTINE SALT_PLUME_APPLY(
     I            trIdentity, bi, bj,
     I            recip_hFac_arg,
     I            tracer,
     U            trStar,
     I            myTime, myIter, myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE SALT_PLUME_APPLY
C     | o Apply the salt_pume-transport to tracer field
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "GRID.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "SALT_PLUME.h"

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     trIdentity :: tracer identification number
C     bi,bj      :: Tile indices
C     recip_drF  :: Reciprol of cell thickness
C     recip_hFac_arg :: Reciprol of cell open-depth factor
C     tracer     :: tracer field at current time (input)
C     trStar     :: future tracer field (modified)
C     myTime     :: Current time in simulation
C     myIter     :: Current time-step number
C     myThid     :: my Thread Id. number

      INTEGER trIdentity
      INTEGER bi, bj
      _RS recip_hFac_arg(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL tracer        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL trStar        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)

      _RL     myTime
      INTEGER myIter, myThid

#ifdef ALLOW_SALT_PLUME
#ifdef SALT_PLUME_VOLUME

C     !LOCAL VARIABLES:
C     === Local variables ===
C     msgBuf     :: Informational/error message buffer
C      CHARACTER*(MAX_LEN_MBUF) msgBuf
      INTEGER i, j, k
      INTEGER upward
      _RL gTr_Surf2kLev, gTr_Below2kLev, gTr_kLev2Above
      _RL SurfVal, SEAICE_Tfrz, Const, recip_Const
      integer kp1

#ifdef ALLOW_DIAGNOSTICS
      CHARACTER*8 diagName
      CHARACTER*5 diagSufx
      LOGICAL     doDiagSPtend
      _RL         plumetend(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
C-    Functions:
      LOGICAL  DIAGNOSTICS_IS_ON
      EXTERNAL DIAGNOSTICS_IS_ON
#endif /* ALLOW_DIAGNOSTICS */

CEOP

      SEAICE_Tfrz = -1.96 _d 0

      IF ( trIdentity.EQ.1 .OR. trIdentity.EQ.2 ) THEN
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

c      upward = rkSign*NINT(-gravitySign)
       upward = 1
       IF (usingZCoords) upward = -1

       IF ( trIdentity.EQ.1 ) THEN
        SurfVal    = SEAICE_Tfrz
        Const      = HeatCapacity_Cp
        recip_Const=recip_Cp
#ifdef ALLOW_DIAGNOSTICS
        IF ( useDiagnostics ) diagSufx = 'EPtnd'
#endif /* ALLOW_DIAGNOSTICS */
       ENDIF
       IF ( trIdentity.EQ.2 ) THEN
        SurfVal    = 0. _d 0
        Const      = 1. _d 0
        recip_Const= 1. _d 0
#ifdef ALLOW_DIAGNOSTICS
        IF ( useDiagnostics ) diagSufx = 'SPtnd'
#endif /* ALLOW_DIAGNOSTICS */
       ENDIF
       
#ifdef ALLOW_DIAGNOSTICS
       doDiagSPtend = .FALSE.
       diagName = 'oce'
       IF ( useDiagnostics ) THEN
C--   Set diagnostic suffix for the current tracer
        diagName = 'oce'//diagSufx
        doDiagSPtend = DIAGNOSTICS_IS_ON(diagName,myThid)
        IF ( doDiagSPtend ) THEN
         DO k=1,Nr
          DO j=1,OLy
           DO i=1,OLx
            plumetend(i,j,k) = 0. _d 0
           ENDDO
          ENDDO
         ENDDO
        ENDIF
       ELSE
        doDiagSPtend = .FALSE.
       ENDIF
#endif /* ALLOW_DIAGNOSTICS */

       DO k=Nr,1,-1
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          IF(trIdentity.EQ.2) SurfVal = SPbrineSalt(i,j,bi,bj)
Catn:
C salt: [m/s * psu * kg/m3] = [kg/s/m2 psu] = unit of saltPlumeFlux
C theta:[m/s * kg/m3 * J/kg/degC * degC] = [W/m2]
          gTr_Surf2kLev = dSPvolSurf2kLev(i,j,k,bi,bj)
     &                  * SurfVal * rhoConst * Const
          gTr_Below2kLev= dSPvolBelow2kLev(i,j,k,bi,bj)
     &                  * tracer(i,j,k,bi,bj) * rhoConst * Const
Catn: gTr_kLev2Above works even for kLev=1 because this is how much
C volume of original [salinity,heat] associated with [SSS,SST] 
C was replaced by same volume of brine [salt,heat(from SEAICE_Tfrz)].
C Note: by design, dSPvolkLev2Above already is negative
          gTr_kLev2Above= dSPvolkLev2Above(i,j,k,bi,bj)
     &                  * tracer(i,j,k,bi,bj) * rhoConst * Const
          plumetend(i,j,k)=
     &          gTr_Surf2kLev + gTr_Below2kLev + gTr_kLev2Above
          trStar(i,j,k,bi,bj)=trStar(i,j,k,bi,bj)+plumetend(i,j,k)
     &         *recip_drF(k)*mass2rUnit*recip_hFac_arg(i,j,k,bi,bj)
     &         *recip_Const
         ENDDO
        ENDDO
       ENDDO

#ifdef ALLOW_DIAGNOSTICS
       IF ( doDiagSPtend )
     &  CALL DIAGNOSTICS_FILL(plumetend, diagName, 0,Nr,2,bi,bj,myThid)
#endif /* ALLOW_DIAGNOSTICS */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C--   end if trIdentity=[1,2]
      ENDIF

#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_apply.F,v 1.2 2014/04/22 10:32:25 atn Exp $
2	atn	1.1	C $Name: $
3
4			#include "SALT_PLUME_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: SALT_PLUME_APPLY
8			C !INTERFACE:
9			SUBROUTINE SALT_PLUME_APPLY(
10			I trIdentity, bi, bj,
11			I recip_hFac_arg,
12			I tracer,
13			U trStar,
14			I myTime, myIter, myThid )
15
16			C !DESCRIPTION: \bv
17			C ==========================================================
18			C \| SUBROUTINE SALT_PLUME_APPLY
19			C \| o Apply the salt_pume-transport to tracer field
20			C ==========================================================
21			C \ev
22
23			C !USES:
24			IMPLICIT NONE
25
26			C === Global variables ===
27			#include "SIZE.h"
28			#include "GRID.h"
29			#include "EEPARAMS.h"
30			#include "PARAMS.h"
31			#include "DYNVARS.h"
32			#include "SALT_PLUME.h"
33
34			C !INPUT/OUTPUT PARAMETERS:
35			C === Routine arguments ===
36			C trIdentity :: tracer identification number
37			C bi,bj :: Tile indices
38			C recip_drF :: Reciprol of cell thickness
39			C recip_hFac_arg :: Reciprol of cell open-depth factor
40			C tracer :: tracer field at current time (input)
41			C trStar :: future tracer field (modified)
42			C myTime :: Current time in simulation
43			C myIter :: Current time-step number
44			C myThid :: my Thread Id. number
45
46			INTEGER trIdentity
47			INTEGER bi, bj
48			_RS recip_hFac_arg(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
49			_RL tracer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
50			_RL trStar (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
51
52			_RL myTime
53			INTEGER myIter, myThid
54
55			#ifdef ALLOW_SALT_PLUME
56			#ifdef SALT_PLUME_VOLUME
57
58			C !LOCAL VARIABLES:
59			C === Local variables ===
60			C msgBuf :: Informational/error message buffer
61			C CHARACTER*(MAX_LEN_MBUF) msgBuf
62			INTEGER i, j, k
63			INTEGER upward
64			_RL gTr_Surf2kLev, gTr_Below2kLev, gTr_kLev2Above
65			_RL SurfVal, SEAICE_Tfrz, Const, recip_Const
66			integer kp1
67
68			#ifdef ALLOW_DIAGNOSTICS
69			CHARACTER*8 diagName
70	atn	1.3	CHARACTER*5 diagSufx
71	atn	1.1	LOGICAL doDiagSPtend
72			_RL plumetend(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
73			C- Functions:
74			LOGICAL DIAGNOSTICS_IS_ON
75			EXTERNAL DIAGNOSTICS_IS_ON
76			#endif /* ALLOW_DIAGNOSTICS */
77
78			CEOP
79
80			SEAICE_Tfrz = -1.96 _d 0
81
82			IF ( trIdentity.EQ.1 .OR. trIdentity.EQ.2 ) THEN
83			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
84
85			c upward = rkSign*NINT(-gravitySign)
86			upward = 1
87			IF (usingZCoords) upward = -1
88
89			IF ( trIdentity.EQ.1 ) THEN
90			SurfVal = SEAICE_Tfrz
91			Const = HeatCapacity_Cp
92			recip_Const=recip_Cp
93			#ifdef ALLOW_DIAGNOSTICS
94			IF ( useDiagnostics ) diagSufx = 'EPtnd'
95			#endif /* ALLOW_DIAGNOSTICS */
96			ENDIF
97	atn	1.3	IF ( trIdentity.EQ.2 ) THEN
98	atn	1.1	SurfVal = 0. _d 0
99			Const = 1. _d 0
100			recip_Const= 1. _d 0
101			#ifdef ALLOW_DIAGNOSTICS
102			IF ( useDiagnostics ) diagSufx = 'SPtnd'
103			#endif /* ALLOW_DIAGNOSTICS */
104			ENDIF
105
106			#ifdef ALLOW_DIAGNOSTICS
107	atn	1.3	doDiagSPtend = .FALSE.
108			diagName = 'oce'
109	atn	1.1	IF ( useDiagnostics ) THEN
110			C-- Set diagnostic suffix for the current tracer
111			diagName = 'oce'//diagSufx
112			doDiagSPtend = DIAGNOSTICS_IS_ON(diagName,myThid)
113			IF ( doDiagSPtend ) THEN
114			DO k=1,Nr
115			DO j=1,OLy
116			DO i=1,OLx
117			plumetend(i,j,k) = 0. _d 0
118			ENDDO
119			ENDDO
120			ENDDO
121			ENDIF
122			ELSE
123			doDiagSPtend = .FALSE.
124			ENDIF
125			#endif /* ALLOW_DIAGNOSTICS */
126
127			DO k=Nr,1,-1
128			DO j=1-OLy,sNy+OLy
129			DO i=1-OLx,sNx+OLx
130			IF(trIdentity.EQ.2) SurfVal = SPbrineSalt(i,j,bi,bj)
131			Catn:
132			C salt: [m/s * psu * kg/m3] = [kg/s/m2 psu] = unit of saltPlumeFlux
133			C theta:[m/s * kg/m3 * J/kg/degC * degC] = [W/m2]
134			gTr_Surf2kLev = dSPvolSurf2kLev(i,j,k,bi,bj)
135			& * SurfVal * rhoConst * Const
136			gTr_Below2kLev= dSPvolBelow2kLev(i,j,k,bi,bj)
137	atn	1.2	& * tracer(i,j,k,bi,bj) * rhoConst * Const
138	atn	1.1	Catn: gTr_kLev2Above works even for kLev=1 because this is how much
139			C volume of original [salinity,heat] associated with [SSS,SST]
140			C was replaced by same volume of brine [salt,heat(from SEAICE_Tfrz)].
141			C Note: by design, dSPvolkLev2Above already is negative
142			gTr_kLev2Above= dSPvolkLev2Above(i,j,k,bi,bj)
143	atn	1.2	& * tracer(i,j,k,bi,bj) * rhoConst * Const
144	atn	1.1	plumetend(i,j,k)=
145			& gTr_Surf2kLev + gTr_Below2kLev + gTr_kLev2Above
146	atn	1.2	trStar(i,j,k,bi,bj)=trStar(i,j,k,bi,bj)+plumetend(i,j,k)
147			& recip_drF(k)mass2rUnit*recip_hFac_arg(i,j,k,bi,bj)
148	atn	1.1	& *recip_Const
149			ENDDO
150			ENDDO
151			ENDDO
152
153			#ifdef ALLOW_DIAGNOSTICS
154			IF ( doDiagSPtend )
155			& CALL DIAGNOSTICS_FILL(plumetend, diagName, 0,Nr,2,bi,bj,myThid)
156			#endif /* ALLOW_DIAGNOSTICS */
157
158			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
159			C-- end if trIdentity=[1,2]
160			ENDIF
161
162			#endif /* SALT_PLUME_VOLUME */
163			#endif /* ALLOW_SALT_PLUME */
164
165			RETURN
166			END