model/src/calc_3d_diffusivity.F

C $Header: /u/gcmpack/MITgcm/model/src/calc_3d_diffusivity.F,v 1.6 2007/04/05 22:51:48 dimitri Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: CALC_3D_DIFFUSIVITY
C     !INTERFACE: 
      SUBROUTINE CALC_3D_DIFFUSIVITY( 
     I        bi,bj,iMin,iMax,jMin,jMax,
     I        trIdentity, trUseGMRedi, trUseKPP,
     O        KappaRTr,
     I        myThid)

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE CALC_3D_DIFFUSIVITY
C     | o Calculate net (3D) vertical diffusivity for 1 tracer
C     *==========================================================*
C     | Combines spatially varying diffusion coefficients from
C     | KPP and/or GM and/or convective stability test.
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     == GLobal variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
#ifdef ALLOW_GENERIC_ADVDIFF
#include "GAD.h"
#endif
#ifdef ALLOW_PTRACERS
#include "PTRACERS_SIZE.h"
#include "PTRACERS.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     bi, bj     :: tile indices
C     iMin,iMax  :: Range of points for which calculation is performed.
C     jMin,jMax  :: Range of points for which calculation is performed.
C     trIdentity :: tracer identifier
C     trUseGMRedi:: this tracer use GM-Redi
C     trUseKPP   :: this tracer use KPP
C     myThid     :: Instance number for this innvocation of CALC_3D_DIFFUSIVITY
C     KappaRTr   :: Net diffusivity for this tracer (trIdentity)
      INTEGER bi,bj,iMin,iMax,jMin,jMax
      INTEGER trIdentity
      LOGICAL trUseGMRedi, trUseKPP
      _RL KappaRTr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
      INTEGER myThid

#ifdef ALLOW_GENERIC_ADVDIFF
C     !LOCAL VARIABLES:
C     == Local variables ==
C     i, j, k    :: Loop counters
C     iTr        :: passive tracer index
C     msgBuf     :: message buffer
      INTEGER i,j,k
      _RL KbryanLewis79, KbryanLewisHL, KbryanLewisEQ
      CHARACTER*(MAX_LEN_MBUF) msgBuf
#ifdef ALLOW_PTRACERS
      INTEGER iTr
#endif
CEOP

      IF ( trIdentity.EQ.GAD_TEMPERATURE) THEN

       DO k = 1,Nr
        KbryanLewisHL=diffKrBL79surf+(diffKrBL79deep-diffKrBL79surf)
     &   *( atan( -( rF(k)-diffKrBL79Ho )/diffKrBL79scl )/PI+0.5 _d 0)
        KbryanLewisEQ=diffKrBLEQsurf+(diffKrBLEQdeep-diffKrBLEQsurf)
     &   *( atan( -( rF(k)-diffKrBLEQHo )/diffKrBLEQscl )/PI+0.5 _d 0)
        DO j = 1-Oly,sNy+Oly
         IF ( abs(YC(i,j,bi,bj)) .GT. KbryanLewisLatTransition ) THEN
          KbryanLewis79 = KbryanLewisHL
         ELSE
          KbryanLewis79 = KbryanLewisHL-(KbryanLewisHL-KbryanLewisEQ)*
     &         (1+cos(YC(i,j,bi,bj)*pi/KbryanLewisLatTransition)) / 2
         ENDIF
         DO i = 1-Olx,sNx+Olx
          KappaRTr(i,j,k) =
     &         IVDConvCount(i,j,k,bi,bj)*ivdc_kappa
#if (defined ALLOW_3D_DIFFKR || \
     (defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_DIFFKR_CONTROL)))
     &       + diffKr(i,j,k,bi,bj)
#else
     &       + diffKrNrT(k)
#endif
     &       + KbryanLewis79
         ENDDO
        ENDDO
       ENDDO

      ELSEIF ( trIdentity.EQ.GAD_SALINITY) THEN

       DO k = 1,Nr
        KbryanLewisHL=diffKrBL79surf+(diffKrBL79deep-diffKrBL79surf)
     &   *( atan( -( rF(k)-diffKrBL79Ho )/diffKrBL79scl )/PI+0.5 _d 0)
        KbryanLewisEQ=diffKrBLEQsurf+(diffKrBLEQdeep-diffKrBLEQsurf)
     &   *( atan( -( rF(k)-diffKrBLEQHo )/diffKrBLEQscl )/PI+0.5 _d 0)
        DO j = 1-Oly, sNy+Oly
         IF ( abs(YC(i,j,bi,bj)) .GT. KbryanLewisLatTransition ) THEN
          KbryanLewis79 = KbryanLewisHL
         ELSE
          KbryanLewis79 = KbryanLewisHL-(KbryanLewisHL-KbryanLewisEQ)*
     &         (1+cos(YC(i,j,bi,bj)*pi/KbryanLewisLatTransition)) / 2
         ENDIF
         DO i = 1-Olx, sNx+Olx
          KappaRTr(i,j,k) =
     &         IVDConvCount(i,j,k,bi,bj)*ivdc_kappa
#if (defined ALLOW_3D_DIFFKR || \
     (defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_DIFFKR_CONTROL)))
     &       + diffKr(i,j,k,bi,bj)
#else
     &       + diffKrNrS(k)
#endif
     &       + KbryanLewis79
         ENDDO
        ENDDO
       ENDDO

#ifdef ALLOW_PTRACERS
      ELSEIF ( trIdentity.GE.GAD_TR1
     &   .AND. trIdentity.LT.GAD_TR1+PTRACERS_numInUse) THEN

       iTr = trIdentity - GAD_TR1 + 1
       DO k = 1,Nr
        KbryanLewisHL=diffKrBL79surf+(diffKrBL79deep-diffKrBL79surf)
     &   *( atan( -( rF(k)-diffKrBL79Ho )/diffKrBL79scl )/PI+0.5 _d 0)
        KbryanLewisEQ=diffKrBLEQsurf+(diffKrBLEQdeep-diffKrBLEQsurf)
     &   *( atan( -( rF(k)-diffKrBLEQHo )/diffKrBLEQscl )/PI+0.5 _d 0)
        DO j = 1-Oly, sNy+Oly
         IF ( abs(YC(i,j,bi,bj)) .GT. KbryanLewisLatTransition ) THEN
          KbryanLewis79 = KbryanLewisHL
         ELSE
          KbryanLewis79 = KbryanLewisHL-(KbryanLewisHL-KbryanLewisEQ)*
     &         (1+cos(YC(i,j,bi,bj)*pi/KbryanLewisLatTransition)) / 2
         ENDIF
         DO i = 1-Olx, sNx+Olx
          KappaRTr(i,j,k) =
     &         IVDConvCount(i,j,k,bi,bj)*ivdc_kappa
#if (defined ALLOW_3D_DIFFKR || \
     (defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_DIFFKR_CONTROL)))
     &       + diffKr(i,j,k,bi,bj)
#else
     &       + PTRACERS_diffKrNr(k,iTr)
#endif
     &       + KbryanLewis79
         ENDDO
        ENDDO
       ENDDO
#endif /* ALLOW_PTRACERS */

      ELSE
        WRITE(msgBuf,'(A,I4)')
     &      ' CALC_3D_DIFFUSIVITY: Invalid tracer Id: ',trIdentity
        CALL PRINT_ERROR(msgBuf, myThid)
        STOP 'ABNORMAL END: S/R CALC_3D_DIFFUSIVITY'
      ENDIF

C--   Add physical pacakge contributions:

#ifdef ALLOW_GMREDI
      IF (trUseGMRedi) THEN 
         CALL GMREDI_CALC_DIFF(
     I        bi,bj,iMin,iMax,jMin,jMax,0,Nr,
     U        KappaRTr,
     I        myThid)
      ENDIF
#endif

#ifdef ALLOW_KPP
      IF (trUseKPP) THEN
       IF (trIdentity.EQ.GAD_TEMPERATURE) THEN
         CALL KPP_CALC_DIFF_T(
     I        bi,bj,iMin,iMax,jMin,jMax,0,Nr,
     U        KappaRTr,
     I        myThid)
       ELSE
         CALL KPP_CALC_DIFF_S(
     I        bi,bj,iMin,iMax,jMin,jMax,0,Nr,
     U        KappaRTr,
     I        myThid)
       ENDIF
      ENDIF
#endif

#ifdef ALLOW_PP81
      IF (usePP81) THEN
         CALL PP81_CALC_DIFF(
     I        bi,bj,iMin,iMax,jMin,jMax,0,Nr,
     U        KappaRTr,
     I        myThid)
      ENDIF
#endif

#ifdef ALLOW_MY82
      IF (useMY82) THEN
         CALL MY82_CALC_DIFF(
     I        bi,bj,iMin,iMax,jMin,jMax,0,Nr,
     U        KappaRTr,
     I        myThid)
      ENDIF
#endif
      
#ifdef ALLOW_GGL90
      IF (useGGL90) THEN
         CALL GGL90_CALC_DIFF(
     I        bi,bj,iMin,iMax,jMin,jMax,0,Nr,
     O        KappaRTr,
     I        myThid)
      ENDIF
#endif
      
C-    Apply mask to vertical diffusivity 
C jmc: don't have the impression that masking is needed 
C      but could be removed later if it's the case.
c     DO j = 1-Oly, sNy+Oly
c      DO i = 1-Olx, sNx+Olx
c       KappaRTr(i,j,k) = maskUp(i,j)*KappaRTr(i,j,k)
c      ENDDO
c     ENDDO

#endif /* ALLOW_GENERIC_ADVDIFF */

      RETURN
      END
1	dimitri	1.7	C $Header: /u/gcmpack/MITgcm/model/src/calc_3d_diffusivity.F,v 1.6 2007/04/05 22:51:48 dimitri Exp $
2	jmc	1.1	C $Name: $
3
4			#include "PACKAGES_CONFIG.h"
5			#include "CPP_OPTIONS.h"
6
7			CBOP
8			C !ROUTINE: CALC_3D_DIFFUSIVITY
9			C !INTERFACE:
10			SUBROUTINE CALC_3D_DIFFUSIVITY(
11			I bi,bj,iMin,iMax,jMin,jMax,
12			I trIdentity, trUseGMRedi, trUseKPP,
13			O KappaRTr,
14			I myThid)
15
16			C !DESCRIPTION: \bv
17			C ==========================================================
18			C \| SUBROUTINE CALC_3D_DIFFUSIVITY
19			C \| o Calculate net (3D) vertical diffusivity for 1 tracer
20			C ==========================================================
21			C \| Combines spatially varying diffusion coefficients from
22			C \| KPP and/or GM and/or convective stability test.
23			C ==========================================================
24			C \ev
25
26			C !USES:
27			IMPLICIT NONE
28			C == GLobal variables ==
29			#include "SIZE.h"
30			#include "EEPARAMS.h"
31			#include "PARAMS.h"
32			#include "DYNVARS.h"
33			#include "GRID.h"
34	jmc	1.3	#ifdef ALLOW_GENERIC_ADVDIFF
35	jmc	1.1	#include "GAD.h"
36	jmc	1.3	#endif
37	jmc	1.1	#ifdef ALLOW_PTRACERS
38			#include "PTRACERS_SIZE.h"
39			#include "PTRACERS.h"
40			#endif
41
42			C !INPUT/OUTPUT PARAMETERS:
43			C == Routine arguments ==
44			C bi, bj :: tile indices
45			C iMin,iMax :: Range of points for which calculation is performed.
46			C jMin,jMax :: Range of points for which calculation is performed.
47			C trIdentity :: tracer identifier
48			C trUseGMRedi:: this tracer use GM-Redi
49			C trUseKPP :: this tracer use KPP
50			C myThid :: Instance number for this innvocation of CALC_3D_DIFFUSIVITY
51			C KappaRTr :: Net diffusivity for this tracer (trIdentity)
52			INTEGER bi,bj,iMin,iMax,jMin,jMax
53			INTEGER trIdentity
54			LOGICAL trUseGMRedi, trUseKPP
55			_RL KappaRTr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
56			INTEGER myThid
57
58	jmc	1.3	#ifdef ALLOW_GENERIC_ADVDIFF
59	jmc	1.1	C !LOCAL VARIABLES:
60			C == Local variables ==
61			C i, j, k :: Loop counters
62			C iTr :: passive tracer index
63			C msgBuf :: message buffer
64			INTEGER i,j,k
65	dimitri	1.7	_RL KbryanLewis79, KbryanLewisHL, KbryanLewisEQ
66	jmc	1.1	CHARACTER*(MAX_LEN_MBUF) msgBuf
67	jmc	1.4	#ifdef ALLOW_PTRACERS
68			INTEGER iTr
69			#endif
70	jmc	1.1	CEOP
71
72			IF ( trIdentity.EQ.GAD_TEMPERATURE) THEN
73
74			DO k = 1,Nr
75	dimitri	1.7	KbryanLewisHL=diffKrBL79surf+(diffKrBL79deep-diffKrBL79surf)
76			& *( atan( -( rF(k)-diffKrBL79Ho )/diffKrBL79scl )/PI+0.5 _d 0)
77			KbryanLewisEQ=diffKrBLEQsurf+(diffKrBLEQdeep-diffKrBLEQsurf)
78			& *( atan( -( rF(k)-diffKrBLEQHo )/diffKrBLEQscl )/PI+0.5 _d 0)
79	jmc	1.1	DO j = 1-Oly,sNy+Oly
80	dimitri	1.7	IF ( abs(YC(i,j,bi,bj)) .GT. KbryanLewisLatTransition ) THEN
81			KbryanLewis79 = KbryanLewisHL
82			ELSE
83			KbryanLewis79 = KbryanLewisHL-(KbryanLewisHL-KbryanLewisEQ)*
84			& (1+cos(YC(i,j,bi,bj)*pi/KbryanLewisLatTransition)) / 2
85			ENDIF
86	jmc	1.1	DO i = 1-Olx,sNx+Olx
87			KappaRTr(i,j,k) =
88			& IVDConvCount(i,j,k,bi,bj)*ivdc_kappa
89	dimitri	1.6	#if (defined ALLOW_3D_DIFFKR \|\| \
90			(defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_DIFFKR_CONTROL)))
91	jmc	1.1	& + diffKr(i,j,k,bi,bj)
92			#else
93	dimitri	1.6	& + diffKrNrT(k)
94	jmc	1.1	#endif
95			& + KbryanLewis79
96			ENDDO
97			ENDDO
98			ENDDO
99
100			ELSEIF ( trIdentity.EQ.GAD_SALINITY) THEN
101
102			DO k = 1,Nr
103	dimitri	1.7	KbryanLewisHL=diffKrBL79surf+(diffKrBL79deep-diffKrBL79surf)
104			& *( atan( -( rF(k)-diffKrBL79Ho )/diffKrBL79scl )/PI+0.5 _d 0)
105			KbryanLewisEQ=diffKrBLEQsurf+(diffKrBLEQdeep-diffKrBLEQsurf)
106			& *( atan( -( rF(k)-diffKrBLEQHo )/diffKrBLEQscl )/PI+0.5 _d 0)
107	jmc	1.1	DO j = 1-Oly, sNy+Oly
108	dimitri	1.7	IF ( abs(YC(i,j,bi,bj)) .GT. KbryanLewisLatTransition ) THEN
109			KbryanLewis79 = KbryanLewisHL
110			ELSE
111			KbryanLewis79 = KbryanLewisHL-(KbryanLewisHL-KbryanLewisEQ)*
112			& (1+cos(YC(i,j,bi,bj)*pi/KbryanLewisLatTransition)) / 2
113			ENDIF
114	jmc	1.1	DO i = 1-Olx, sNx+Olx
115			KappaRTr(i,j,k) =
116			& IVDConvCount(i,j,k,bi,bj)*ivdc_kappa
117	dimitri	1.6	#if (defined ALLOW_3D_DIFFKR \|\| \
118			(defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_DIFFKR_CONTROL)))
119	jmc	1.1	& + diffKr(i,j,k,bi,bj)
120			#else
121	dimitri	1.6	& + diffKrNrS(k)
122	jmc	1.1	#endif
123			& + KbryanLewis79
124			ENDDO
125			ENDDO
126			ENDDO
127
128			#ifdef ALLOW_PTRACERS
129			ELSEIF ( trIdentity.GE.GAD_TR1
130			& .AND. trIdentity.LT.GAD_TR1+PTRACERS_numInUse) THEN
131
132			iTr = trIdentity - GAD_TR1 + 1
133			DO k = 1,Nr
134	dimitri	1.7	KbryanLewisHL=diffKrBL79surf+(diffKrBL79deep-diffKrBL79surf)
135			& *( atan( -( rF(k)-diffKrBL79Ho )/diffKrBL79scl )/PI+0.5 _d 0)
136			KbryanLewisEQ=diffKrBLEQsurf+(diffKrBLEQdeep-diffKrBLEQsurf)
137			& *( atan( -( rF(k)-diffKrBLEQHo )/diffKrBLEQscl )/PI+0.5 _d 0)
138	jmc	1.1	DO j = 1-Oly, sNy+Oly
139	dimitri	1.7	IF ( abs(YC(i,j,bi,bj)) .GT. KbryanLewisLatTransition ) THEN
140			KbryanLewis79 = KbryanLewisHL
141			ELSE
142			KbryanLewis79 = KbryanLewisHL-(KbryanLewisHL-KbryanLewisEQ)*
143			& (1+cos(YC(i,j,bi,bj)*pi/KbryanLewisLatTransition)) / 2
144			ENDIF
145	jmc	1.1	DO i = 1-Olx, sNx+Olx
146			KappaRTr(i,j,k) =
147			& IVDConvCount(i,j,k,bi,bj)*ivdc_kappa
148	dimitri	1.6	#if (defined ALLOW_3D_DIFFKR \|\| \
149			(defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_DIFFKR_CONTROL)))
150	jmc	1.1	& + diffKr(i,j,k,bi,bj)
151			#else
152			& + PTRACERS_diffKrNr(k,iTr)
153			#endif
154			& + KbryanLewis79
155			ENDDO
156			ENDDO
157			ENDDO
158			#endif /* ALLOW_PTRACERS */
159
160			ELSE
161			WRITE(msgBuf,'(A,I4)')
162			& ' CALC_3D_DIFFUSIVITY: Invalid tracer Id: ',trIdentity
163			CALL PRINT_ERROR(msgBuf, myThid)
164			STOP 'ABNORMAL END: S/R CALC_3D_DIFFUSIVITY'
165			ENDIF
166
167			C-- Add physical pacakge contributions:
168
169			#ifdef ALLOW_GMREDI
170			IF (trUseGMRedi) THEN
171			CALL GMREDI_CALC_DIFF(
172			I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
173			U KappaRTr,
174			I myThid)
175			ENDIF
176			#endif
177
178			#ifdef ALLOW_KPP
179			IF (trUseKPP) THEN
180			IF (trIdentity.EQ.GAD_TEMPERATURE) THEN
181			CALL KPP_CALC_DIFF_T(
182	jmc	1.2	I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
183	jmc	1.1	U KappaRTr,
184			I myThid)
185			ELSE
186			CALL KPP_CALC_DIFF_S(
187	jmc	1.2	I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
188	jmc	1.1	U KappaRTr,
189			I myThid)
190			ENDIF
191			ENDIF
192			#endif
193
194			#ifdef ALLOW_PP81
195			IF (usePP81) THEN
196			CALL PP81_CALC_DIFF(
197	jmc	1.2	I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
198	jmc	1.1	U KappaRTr,
199			I myThid)
200			ENDIF
201			#endif
202
203			#ifdef ALLOW_MY82
204			IF (useMY82) THEN
205			CALL MY82_CALC_DIFF(
206	jmc	1.2	I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
207	jmc	1.1	U KappaRTr,
208			I myThid)
209			ENDIF
210			#endif
211
212			#ifdef ALLOW_GGL90
213			IF (useGGL90) THEN
214			CALL GGL90_CALC_DIFF(
215	jmc	1.2	I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
216	jmc	1.1	O KappaRTr,
217			I myThid)
218			ENDIF
219			#endif
220
221			C- Apply mask to vertical diffusivity
222			C jmc: don't have the impression that masking is needed
223			C but could be removed later if it's the case.
224			c DO j = 1-Oly, sNy+Oly
225			c DO i = 1-Olx, sNx+Olx
226			c KappaRTr(i,j,k) = maskUp(i,j)*KappaRTr(i,j,k)
227			c ENDDO
228			c ENDDO
229
230	jmc	1.3	#endif /* ALLOW_GENERIC_ADVDIFF */
231
232	jmc	1.1	RETURN
233			END