atnguyen/code_21Dec2012_saltplume/salt_plume_calc_depth.F

C $Header: /u/gcmpack/MITgcm/pkg/salt_plume/salt_plume_calc_depth.F,v 1.7 2012/04/28 22:18:23 heimbach Exp $
C $Name:  $

#include "SALT_PLUME_OPTIONS.h"

CBOP
C     !ROUTINE: SALT_PLUME_CALC_DEPTH
C     !INTERFACE:
      SUBROUTINE SALT_PLUME_CALC_DEPTH(
     I                       rhoSurf, sigmaR,
     I                       bi, bj, myTime, myIter, myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R SALT_PLUME_CALC_DEPTH
C     | o Compute depth of penetration of salt plumes rejected
C     |   during sea ice growth
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
#include "DYNVARS.h"
#include "SALT_PLUME.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
# include "tamc_keys.h"
#endif /* ALLOW_AUTODIFF_TAMC */

C     !INPUT/OUTPUT PARAMETERS:
C     rhoSurf   :: Surface density anomaly
C     sigmaR    :: Vertical gradient of potential density
C     bi,bj     :: tile indices
C     myTime    :: Current time in simulation
C     myIter    :: Current iteration number in simulation
C     myThid    :: my Thread Id number
      _RL     rhoSurf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      INTEGER bi, bj
      _RL     myTime
      INTEGER myIter
      INTEGER myThid
CEOP

#ifdef ALLOW_SALT_PLUME

C     !LOCAL VARIABLES:
C     i,j :: Loop counters
      INTEGER i,j,k
      _RL     rhoBigNb, tmpFac
      _RL     rhoMxL(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     rhoKm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     rhoLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     GG, GGm1
      _RL     SPIND (1-OLx:sNx+OLx,1-OLy:sNy+OLy)

#ifdef ALLOW_AUTODIFF_TAMC
        act1 = bi - myBxLo(myThid)
        max1 = myBxHi(myThid) - myBxLo(myThid) + 1
        act2 = bj - myByLo(myThid)
        max2 = myByHi(myThid) - myByLo(myThid) + 1
        act3 = myThid - 1
        max3 = nTx*nTy
        act4 = ikey_dynamics - 1
        ikey = (act1 + 1) + act2*max1
     &                    + act3*max1*max2
     &                    + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */

C Initializing the saltplume depth to bottom topography
      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
        SaltPlumeDepth(i,j,bi,bj) = rF(1)-R_low(I,J,bi,bj)
        SPIND(i,j)  = 0. _d 0
        rhoKm1(i,j) = 0. _d 0
        rhoMxL(i,j) = 0. _d 0
       ENDDO
      ENDDO

C CriterionType 1 = use delta_rho to determine salt plume depth
      IF (CriterionType.EQ.1) THEN

       rhoBigNb  = rhoConst*1. _d 10
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
         SaltPlumeDepth(i,j,bi,bj) = rF(1)-R_low(I,J,bi,bj)
         rhoKm1(i,j) = rhoSurf(i,j)
         rhoMxL(i,j) = rhoSurf(i,j) + SaltPlumeCriterion
        ENDDO
       ENDDO

       DO k = 2,Nr
#ifdef ALLOW_AUTODIFF_TAMC
          kkey = (ikey-1)*Nr + k
CADJ STORE rhoKm1(:,:) = comlev1_bibj_k, 
CADJ &     key=kkey, byte=isbyte, kind = isbyte
CADJ STORE rhoMxL(:,:) = comlev1_bibj_k, 
CADJ &     key=kkey, byte=isbyte, kind = isbyte
CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, 
CADJ &     key=kkey, byte=isbyte, kind = isbyte
CADJ STORE salt(:,:,k,bi,bj) = comlev1_bibj_k, 
CADJ &     key=kkey, byte=isbyte, kind = isbyte
#endif
C-     potential density (reference level = surface level)
        CALL FIND_RHO_2D(
     I       1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1,
     I       theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj),
     O       rhoLoc,
     I       K, bi, bj, myThid )

        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          IF ( k.LE.klowC(i,j,bi,bj) .AND.
     &         rhoLoc(i,j).GE.rhoMxL(i,j) ) THEN
           IF ( rhoLoc(i,j).GT.rhoKm1(i,j) ) THEN
              tmpFac = ( rhoMxL(i,j) - rhoKm1(i,j) )
     &             / ( rhoLoc(i,j) - rhoKm1(i,j) )
           ELSE
              tmpFac = 0.
           ENDIF
           SaltPlumeDepth(i,j,bi,bj) = rF(1)-rC(k-1)+tmpFac*drC(k)
           rhoMxL(i,j) = rhoBigNb
          ELSE
            rhoKm1(i,j) = rhoLoc(i,j)
          ENDIF
         ENDDO
        ENDDO
       ENDDO

      ELSEIF ( CriterionType.EQ.2 ) THEN

        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          SaltPlumeDepth(i,j,bi,bj) = rF(1)-R_low(I,J,bi,bj)
          SPIND(i,j)  = 0. _d 0
         ENDDO
        ENDDO

        DO k=2,Nr
#ifdef ALLOW_AUTODIFF_TAMC
          kkey = (ikey-1)*Nr + k
CADJ STORE SPIND(:,:) = comlev1_bibj_k, 
CADJ &     key=kkey, byte=isbyte, kind = isbyte
#endif
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           GG  =-1.0*sigmaR(i,j,k)
           GGm1=-1.0*sigmaR(i,j,k-1)
           IF ( k.LE.klowC(i,j,bi,bj) .AND.
     &          GG.GE.SaltPlumeCriterion ) THEN
            IF (GGm1.LE.SaltPlumeCriterion) THEN
             tmpFac = (SaltPlumeCriterion - GGm1)
     &              / (GG                 - GGm1)
             IF(SPIND(i,j) .LT. 0.5) THEN
               SaltPlumeDepth(i,j,bi,bj) = rF(1)-rC(k-1)+tmpFac*drC(k)
               SPIND(i,j)=1.
             ENDIF
            ELSE
             tmpFac = 0.
            ENDIF
           ENDIF
          ENDDO
         ENDDO
        ENDDO

        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          SaltPlumeDepth(i,j,bi,bj) =
     &       min( SaltPlumeDepth(i,j,bi,bj)*SPovershoot,
     &            rF(1)-R_low(i,j,bi,bj) )
         ENDDO
        ENDDO
      ENDIF

C Make sure that the deepest SaltPlumeDepth is bottom topography:
      DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          SaltPlumeDepth(i,j,bi,bj) =
     &       min( SaltPlumeDepth(i,j,bi,bj),
     &            rF(1)-R_low(i,j,bi,bj) )
         ENDDO
      ENDDO

C#ifdef ALLOW_DIAGNOSTICS
C      IF ( useDiagnostics )
C     &      CALL SALT_PLUME_DIAGNOSTICS_FILL(bi,bj,myThid)
C#endif /* ALLOW_DIAGNOSTICS */

#endif /* ALLOW_SALT_PLUME */

      RETURN
      END
1	atn	1.3	C $Header: /u/gcmpack/MITgcm/pkg/salt_plume/salt_plume_calc_depth.F,v 1.7 2012/04/28 22:18:23 heimbach Exp $
2	atn	1.1	C $Name: $
3
4			#include "SALT_PLUME_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: SALT_PLUME_CALC_DEPTH
8			C !INTERFACE:
9			SUBROUTINE SALT_PLUME_CALC_DEPTH(
10			I rhoSurf, sigmaR,
11			I bi, bj, myTime, myIter, myThid )
12
13			C !DESCRIPTION: \bv
14			C ==========================================================
15			C \| S/R SALT_PLUME_CALC_DEPTH
16			C \| o Compute depth of penetration of salt plumes rejected
17			C \| during sea ice growth
18			C ==========================================================
19			C \ev
20
21			C !USES:
22			IMPLICIT NONE
23			C == Global variables ==
24			#include "SIZE.h"
25			#include "EEPARAMS.h"
26			#include "PARAMS.h"
27			#include "GRID.h"
28			#include "SURFACE.h"
29			#include "DYNVARS.h"
30			#include "SALT_PLUME.h"
31			#ifdef ALLOW_AUTODIFF_TAMC
32			# include "tamc.h"
33			# include "tamc_keys.h"
34			#endif /* ALLOW_AUTODIFF_TAMC */
35
36			C !INPUT/OUTPUT PARAMETERS:
37			C rhoSurf :: Surface density anomaly
38			C sigmaR :: Vertical gradient of potential density
39			C bi,bj :: tile indices
40			C myTime :: Current time in simulation
41			C myIter :: Current iteration number in simulation
42			C myThid :: my Thread Id number
43			_RL rhoSurf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
44			_RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
45			INTEGER bi, bj
46			_RL myTime
47			INTEGER myIter
48			INTEGER myThid
49			CEOP
50
51			#ifdef ALLOW_SALT_PLUME
52
53			C !LOCAL VARIABLES:
54			C i,j :: Loop counters
55			INTEGER i,j,k
56			_RL rhoBigNb, tmpFac
57			_RL rhoMxL(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58			_RL rhoKm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59			_RL rhoLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
60			_RL GG, GGm1
61			_RL SPIND (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62
63			#ifdef ALLOW_AUTODIFF_TAMC
64			act1 = bi - myBxLo(myThid)
65			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
66			act2 = bj - myByLo(myThid)
67			max2 = myByHi(myThid) - myByLo(myThid) + 1
68			act3 = myThid - 1
69			max3 = nTx*nTy
70			act4 = ikey_dynamics - 1
71			ikey = (act1 + 1) + act2*max1
72			& + act3max1max2
73			& + act4max1max2*max3
74			#endif /* ALLOW_AUTODIFF_TAMC */
75
76			C Initializing the saltplume depth to bottom topography
77			DO j=1-Oly,sNy+Oly
78			DO i=1-Olx,sNx+Olx
79			SaltPlumeDepth(i,j,bi,bj) = rF(1)-R_low(I,J,bi,bj)
80			SPIND(i,j) = 0. _d 0
81			rhoKm1(i,j) = 0. _d 0
82			rhoMxL(i,j) = 0. _d 0
83			ENDDO
84			ENDDO
85
86			C CriterionType 1 = use delta_rho to determine salt plume depth
87			IF (CriterionType.EQ.1) THEN
88
89			rhoBigNb = rhoConst*1. _d 10
90			DO j=1-Oly,sNy+Oly
91			DO i=1-Olx,sNx+Olx
92			SaltPlumeDepth(i,j,bi,bj) = rF(1)-R_low(I,J,bi,bj)
93			rhoKm1(i,j) = rhoSurf(i,j)
94			rhoMxL(i,j) = rhoSurf(i,j) + SaltPlumeCriterion
95			ENDDO
96			ENDDO
97
98			DO k = 2,Nr
99			#ifdef ALLOW_AUTODIFF_TAMC
100			kkey = (ikey-1)*Nr + k
101			CADJ STORE rhoKm1(:,:) = comlev1_bibj_k,
102			CADJ & key=kkey, byte=isbyte, kind = isbyte
103			CADJ STORE rhoMxL(:,:) = comlev1_bibj_k,
104			CADJ & key=kkey, byte=isbyte, kind = isbyte
105			CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k,
106			CADJ & key=kkey, byte=isbyte, kind = isbyte
107			CADJ STORE salt(:,:,k,bi,bj) = comlev1_bibj_k,
108			CADJ & key=kkey, byte=isbyte, kind = isbyte
109			#endif
110			C- potential density (reference level = surface level)
111			CALL FIND_RHO_2D(
112			I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1,
113			I theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj),
114			O rhoLoc,
115			I K, bi, bj, myThid )
116
117			DO j=1-Oly,sNy+Oly
118			DO i=1-Olx,sNx+Olx
119			IF ( k.LE.klowC(i,j,bi,bj) .AND.
120			& rhoLoc(i,j).GE.rhoMxL(i,j) ) THEN
121			IF ( rhoLoc(i,j).GT.rhoKm1(i,j) ) THEN
122			tmpFac = ( rhoMxL(i,j) - rhoKm1(i,j) )
123			& / ( rhoLoc(i,j) - rhoKm1(i,j) )
124			ELSE
125			tmpFac = 0.
126			ENDIF
127			SaltPlumeDepth(i,j,bi,bj) = rF(1)-rC(k-1)+tmpFac*drC(k)
128			rhoMxL(i,j) = rhoBigNb
129			ELSE
130			rhoKm1(i,j) = rhoLoc(i,j)
131			ENDIF
132			ENDDO
133			ENDDO
134			ENDDO
135
136			ELSEIF ( CriterionType.EQ.2 ) THEN
137
138			DO j=1-Oly,sNy+Oly
139			DO i=1-Olx,sNx+Olx
140			SaltPlumeDepth(i,j,bi,bj) = rF(1)-R_low(I,J,bi,bj)
141			SPIND(i,j) = 0. _d 0
142			ENDDO
143			ENDDO
144
145			DO k=2,Nr
146			#ifdef ALLOW_AUTODIFF_TAMC
147			kkey = (ikey-1)*Nr + k
148			CADJ STORE SPIND(:,:) = comlev1_bibj_k,
149			CADJ & key=kkey, byte=isbyte, kind = isbyte
150			#endif
151			DO j=1-Oly,sNy+Oly
152			DO i=1-Olx,sNx+Olx
153			GG =-1.0*sigmaR(i,j,k)
154			GGm1=-1.0*sigmaR(i,j,k-1)
155			IF ( k.LE.klowC(i,j,bi,bj) .AND.
156			& GG.GE.SaltPlumeCriterion ) THEN
157			IF (GGm1.LE.SaltPlumeCriterion) THEN
158			tmpFac = (SaltPlumeCriterion - GGm1)
159			& / (GG - GGm1)
160			IF(SPIND(i,j) .LT. 0.5) THEN
161			SaltPlumeDepth(i,j,bi,bj) = rF(1)-rC(k-1)+tmpFac*drC(k)
162			SPIND(i,j)=1.
163			ENDIF
164			ELSE
165			tmpFac = 0.
166			ENDIF
167			ENDIF
168			ENDDO
169			ENDDO
170			ENDDO
171
172			DO j=1-Oly,sNy+Oly
173			DO i=1-Olx,sNx+Olx
174			SaltPlumeDepth(i,j,bi,bj) =
175			& min( SaltPlumeDepth(i,j,bi,bj)*SPovershoot,
176			& rF(1)-R_low(i,j,bi,bj) )
177			ENDDO
178			ENDDO
179			ENDIF
180
181			C Make sure that the deepest SaltPlumeDepth is bottom topography:
182			DO j=1-Oly,sNy+Oly
183			DO i=1-Olx,sNx+Olx
184			SaltPlumeDepth(i,j,bi,bj) =
185			& min( SaltPlumeDepth(i,j,bi,bj),
186			& rF(1)-R_low(i,j,bi,bj) )
187			ENDDO
188			ENDDO
189
190	atn	1.3	C#ifdef ALLOW_DIAGNOSTICS
191			C IF ( useDiagnostics )
192			C & CALL SALT_PLUME_DIAGNOSTICS_FILL(bi,bj,myThid)
193			C#endif /* ALLOW_DIAGNOSTICS */
194	atn	1.1
195			#endif /* ALLOW_SALT_PLUME */
196
197			RETURN
198			END