model/src/calc_3d_diffusivity.F

C $Header: /u/gcmpack/MITgcm/model/src/calc_3d_diffusivity.F,v 1.14 2008/02/12 23:04:44 dfer 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_PARAMS.h"
#endif
#ifdef ALLOW_LONGSTEP
#include "LONGSTEP.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
#ifdef ALLOW_BL79_LAT_VARY
      _RL KbryanLewisEQ
#endif
      CHARACTER*(MAX_LEN_MBUF) msgBuf
#ifdef ALLOW_PTRACERS
      INTEGER iTr
#endif
CEOP

      IF ( .NOT. trUseKPP ) THEN
       DO k = 1,Nr
        KbryanLewis79=diffKrBL79surf+(diffKrBL79deep-diffKrBL79surf)
     &       *(atan(-(rF(k)-diffKrBL79Ho)/diffKrBL79scl)/PI+0.5 _d 0)
#ifdef ALLOW_BL79_LAT_VARY
        KbryanLewisEQ=diffKrBLEQsurf+(diffKrBLEQdeep-diffKrBLEQsurf)
     &       *(atan(-(rF(k)-diffKrBLEQHo)/diffKrBLEQscl)/PI+0.5 _d 0)
#endif
        DO j = 1-Oly,sNy+Oly
         DO i = 1-Olx,sNx+Olx
#ifdef ALLOW_LONGSTEP
          IF ( trIdentity .GE. GAD_TR1) THEN
           KappaRTr(i,j,k) =
     &         LS_IVDConvCount(i,j,k,bi,bj)*ivdc_kappa
     &         + KbryanLewis79
#ifdef ALLOW_BL79_LAT_VARY
     &         + (KbryanLewisEQ-KbryanLewis79)*BL79LatArray(i,j,bi,bj)
#endif
          ELSE
#else
          IF ( .TRUE. ) THEN
#endif /* ALLOW_LONGSTEP */
           KappaRTr(i,j,k) =
     &         IVDConvCount(i,j,k,bi,bj)*ivdc_kappa
     &         + KbryanLewis79
#ifdef ALLOW_BL79_LAT_VARY
     &         + (KbryanLewisEQ-KbryanLewis79)*BL79LatArray(i,j,bi,bj)
#endif
          ENDIF
         ENDDO
        ENDDO
       ENDDO
       IF ( trIdentity.EQ.GAD_TEMPERATURE ) THEN
        DO k = 1,Nr
         DO j = 1-Oly,sNy+Oly
          DO i = 1-Olx,sNx+Olx
           KappaRTr(i,j,k) = KappaRTr(i,j,k)
#if (defined ALLOW_3D_DIFFKR || \
     (defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_DIFFKR_CONTROL)))
     &          + diffKr(i,j,k,bi,bj)
#else
     &          + diffKrNrT(k)
#endif
          ENDDO
         ENDDO
        ENDDO
       ELSEIF ( trIdentity.EQ.GAD_SALINITY) THEN
        DO k = 1,Nr
         DO j = 1-Oly, sNy+Oly
          DO i = 1-Olx, sNx+Olx
           KappaRTr(i,j,k) = KappaRTr(i,j,k)
#if (defined ALLOW_3D_DIFFKR || \
     (defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_DIFFKR_CONTROL)))
     &          + diffKr(i,j,k,bi,bj)
#else
     &          + diffKrNrS(k)
#endif
          ENDDO
         ENDDO
        ENDDO
#ifdef ALLOW_PTRACERS
       ELSEIF ( trIdentity.GE.GAD_TR1) THEN

        iTr = trIdentity - GAD_TR1 + 1
        DO k = 1,Nr
         DO j = 1-Oly, sNy+Oly
          DO i = 1-Olx, sNx+Olx
           KappaRTr(i,j,k) = KappaRTr(i,j,k)
#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
          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
      ENDIF

C--   Add physical pacakge contributions:

#ifdef ALLOW_KPP
      IF (trUseKPP) THEN
C--   Set vertical diffusivity contribution from KPP
       IF (trIdentity.EQ.GAD_TEMPERATURE) THEN
         CALL KPP_CALC_DIFF_T(
     I        bi,bj,iMin,iMax,jMin,jMax,0,Nr,
     O        KappaRTr,
     I        myThid)
       ELSEIF (trIdentity.EQ.GAD_SALINITY) THEN
         CALL KPP_CALC_DIFF_S(
     I        bi,bj,iMin,iMax,jMin,jMax,0,Nr,
     O        KappaRTr,
     I        myThid)
       ELSEIF ( trIdentity.GE.GAD_TR1) THEN
         CALL KPP_CALC_DIFF_Ptr(
     I        bi,bj,iMin,iMax,jMin,jMax,0,Nr,
     O        KappaRTr,
     I        myThid)
       ENDIF
#if (defined ALLOW_PTRACERS && ! (defined ALLOW_3D_DIFFKR || \
     (defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_DIFFKR_CONTROL))))
       IF ( trIdentity.GE.GAD_TR1) THEN
        iTr = trIdentity - GAD_TR1 + 1
        DO k = 1,Nr
         DO j = 1-Oly, sNy+Oly
          DO i = 1-Olx, sNx+Olx
           KappaRTr(i,j,k) = KappaRTr(i,j,k)
     &          - diffKrNrS(k) + PTRACERS_diffKrNr(k,iTr)
          ENDDO
         ENDDO
        ENDDO
       ENDIF
#endif
      ENDIF
#endif /* ALLOW_KPP */

#ifdef ALLOW_GMREDI
      IF (trUseGMRedi) THEN
         CALL GMREDI_CALC_DIFF(
     I        bi,bj,iMin,iMax,jMin,jMax,0,Nr,
     U        KappaRTr,
     I        trIdentity,myThid)
      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	jahn	1.15	C $Header: /u/gcmpack/MITgcm/model/src/calc_3d_diffusivity.F,v 1.14 2008/02/12 23:04:44 dfer 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	jmc	1.11	C !INTERFACE:
10			SUBROUTINE CALC_3D_DIFFUSIVITY(
11	jmc	1.1	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	jmc	1.13	#include "PTRACERS_PARAMS.h"
40	jmc	1.1	#endif
41	jahn	1.15	#ifdef ALLOW_LONGSTEP
42			#include "LONGSTEP.h"
43			#endif
44	jmc	1.1
45			C !INPUT/OUTPUT PARAMETERS:
46			C == Routine arguments ==
47			C bi, bj :: tile indices
48			C iMin,iMax :: Range of points for which calculation is performed.
49			C jMin,jMax :: Range of points for which calculation is performed.
50			C trIdentity :: tracer identifier
51			C trUseGMRedi:: this tracer use GM-Redi
52			C trUseKPP :: this tracer use KPP
53			C myThid :: Instance number for this innvocation of CALC_3D_DIFFUSIVITY
54			C KappaRTr :: Net diffusivity for this tracer (trIdentity)
55			INTEGER bi,bj,iMin,iMax,jMin,jMax
56			INTEGER trIdentity
57			LOGICAL trUseGMRedi, trUseKPP
58			_RL KappaRTr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
59			INTEGER myThid
60
61	jmc	1.3	#ifdef ALLOW_GENERIC_ADVDIFF
62	jmc	1.1	C !LOCAL VARIABLES:
63			C == Local variables ==
64			C i, j, k :: Loop counters
65			C iTr :: passive tracer index
66			C msgBuf :: message buffer
67			INTEGER i,j,k
68	jmc	1.11	_RL KbryanLewis79
69			#ifdef ALLOW_BL79_LAT_VARY
70			_RL KbryanLewisEQ
71			#endif
72	jmc	1.1	CHARACTER*(MAX_LEN_MBUF) msgBuf
73	jmc	1.4	#ifdef ALLOW_PTRACERS
74			INTEGER iTr
75			#endif
76	jmc	1.1	CEOP
77
78	dimitri	1.10	IF ( .NOT. trUseKPP ) THEN
79	jmc	1.1	DO k = 1,Nr
80	jmc	1.8	KbryanLewis79=diffKrBL79surf+(diffKrBL79deep-diffKrBL79surf)
81	dimitri	1.10	& *(atan(-(rF(k)-diffKrBL79Ho)/diffKrBL79scl)/PI+0.5 _d 0)
82	dimitri	1.9	#ifdef ALLOW_BL79_LAT_VARY
83			KbryanLewisEQ=diffKrBLEQsurf+(diffKrBLEQdeep-diffKrBLEQsurf)
84	dimitri	1.10	& *(atan(-(rF(k)-diffKrBLEQHo)/diffKrBLEQscl)/PI+0.5 _d 0)
85	dimitri	1.9	#endif
86	jmc	1.1	DO j = 1-Oly,sNy+Oly
87			DO i = 1-Olx,sNx+Olx
88	jahn	1.15	#ifdef ALLOW_LONGSTEP
89			IF ( trIdentity .GE. GAD_TR1) THEN
90			KappaRTr(i,j,k) =
91			& LS_IVDConvCount(i,j,k,bi,bj)*ivdc_kappa
92			& + KbryanLewis79
93			#ifdef ALLOW_BL79_LAT_VARY
94			& + (KbryanLewisEQ-KbryanLewis79)*BL79LatArray(i,j,bi,bj)
95			#endif
96			ELSE
97			#else
98			IF ( .TRUE. ) THEN
99			#endif /* ALLOW_LONGSTEP */
100			KappaRTr(i,j,k) =
101	jmc	1.1	& IVDConvCount(i,j,k,bi,bj)*ivdc_kappa
102	dimitri	1.10	& + KbryanLewis79
103			#ifdef ALLOW_BL79_LAT_VARY
104			& + (KbryanLewisEQ-KbryanLewis79)*BL79LatArray(i,j,bi,bj)
105			#endif
106	jahn	1.15	ENDIF
107	dimitri	1.10	ENDDO
108			ENDDO
109			ENDDO
110			IF ( trIdentity.EQ.GAD_TEMPERATURE ) THEN
111			DO k = 1,Nr
112			DO j = 1-Oly,sNy+Oly
113			DO i = 1-Olx,sNx+Olx
114	jmc	1.11	KappaRTr(i,j,k) = KappaRTr(i,j,k)
115	dimitri	1.6	#if (defined ALLOW_3D_DIFFKR \|\| \
116			(defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_DIFFKR_CONTROL)))
117	dimitri	1.10	& + diffKr(i,j,k,bi,bj)
118	jmc	1.1	#else
119	dimitri	1.10	& + diffKrNrT(k)
120	dimitri	1.9	#endif
121	dimitri	1.10	ENDDO
122	jmc	1.1	ENDDO
123			ENDDO
124	dimitri	1.10	ELSEIF ( trIdentity.EQ.GAD_SALINITY) THEN
125			DO k = 1,Nr
126			DO j = 1-Oly, sNy+Oly
127			DO i = 1-Olx, sNx+Olx
128			KappaRTr(i,j,k) = KappaRTr(i,j,k)
129	dimitri	1.6	#if (defined ALLOW_3D_DIFFKR \|\| \
130			(defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_DIFFKR_CONTROL)))
131	dimitri	1.10	& + diffKr(i,j,k,bi,bj)
132	jmc	1.1	#else
133	dimitri	1.10	& + diffKrNrS(k)
134	dimitri	1.9	#endif
135	dimitri	1.10	ENDDO
136	jmc	1.1	ENDDO
137			ENDDO
138			#ifdef ALLOW_PTRACERS
139	jahn	1.15	ELSEIF ( trIdentity.GE.GAD_TR1) THEN
140	jmc	1.1
141	dimitri	1.10	iTr = trIdentity - GAD_TR1 + 1
142			DO k = 1,Nr
143			DO j = 1-Oly, sNy+Oly
144			DO i = 1-Olx, sNx+Olx
145			KappaRTr(i,j,k) = KappaRTr(i,j,k)
146	dimitri	1.6	#if (defined ALLOW_3D_DIFFKR \|\| \
147			(defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_DIFFKR_CONTROL)))
148	dimitri	1.10	& + diffKr(i,j,k,bi,bj)
149	jmc	1.1	#else
150	dimitri	1.10	& + PTRACERS_diffKrNr(k,iTr)
151	dimitri	1.9	#endif
152	dimitri	1.10	ENDDO
153	jmc	1.1	ENDDO
154			ENDDO
155			#endif /* ALLOW_PTRACERS */
156	dimitri	1.10	ELSE
157	jmc	1.1	WRITE(msgBuf,'(A,I4)')
158	dimitri	1.10	& ' CALC_3D_DIFFUSIVITY: Invalid tracer Id: ',trIdentity
159	jmc	1.1	CALL PRINT_ERROR(msgBuf, myThid)
160			STOP 'ABNORMAL END: S/R CALC_3D_DIFFUSIVITY'
161	dimitri	1.10	ENDIF
162	jmc	1.1	ENDIF
163
164			C-- Add physical pacakge contributions:
165
166			#ifdef ALLOW_KPP
167			IF (trUseKPP) THEN
168	dimitri	1.12	C-- Set vertical diffusivity contribution from KPP
169	jmc	1.1	IF (trIdentity.EQ.GAD_TEMPERATURE) THEN
170			CALL KPP_CALC_DIFF_T(
171	jmc	1.2	I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
172	dimitri	1.12	O KappaRTr,
173	jmc	1.1	I myThid)
174	jahn	1.15	ELSEIF (trIdentity.EQ.GAD_SALINITY) THEN
175	jmc	1.1	CALL KPP_CALC_DIFF_S(
176	jmc	1.2	I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
177	dimitri	1.12	O KappaRTr,
178	jmc	1.1	I myThid)
179	jahn	1.15	ELSEIF ( trIdentity.GE.GAD_TR1) THEN
180			CALL KPP_CALC_DIFF_Ptr(
181			I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
182			O KappaRTr,
183			I myThid)
184	jmc	1.1	ENDIF
185	dimitri	1.10	#if (defined ALLOW_PTRACERS && ! (defined ALLOW_3D_DIFFKR \|\| \
186			(defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_DIFFKR_CONTROL))))
187	jahn	1.15	IF ( trIdentity.GE.GAD_TR1) THEN
188	dimitri	1.10	iTr = trIdentity - GAD_TR1 + 1
189			DO k = 1,Nr
190			DO j = 1-Oly, sNy+Oly
191			DO i = 1-Olx, sNx+Olx
192			KappaRTr(i,j,k) = KappaRTr(i,j,k)
193			& - diffKrNrS(k) + PTRACERS_diffKrNr(k,iTr)
194			ENDDO
195			ENDDO
196			ENDDO
197			ENDIF
198			#endif
199			ENDIF
200			#endif /* ALLOW_KPP */
201
202			#ifdef ALLOW_GMREDI
203	jmc	1.11	IF (trUseGMRedi) THEN
204	dimitri	1.10	CALL GMREDI_CALC_DIFF(
205			I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
206			U KappaRTr,
207	dfer	1.14	I trIdentity,myThid)
208	jmc	1.1	ENDIF
209			#endif
210
211			#ifdef ALLOW_PP81
212			IF (usePP81) THEN
213			CALL PP81_CALC_DIFF(
214	jmc	1.2	I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
215	jmc	1.1	U KappaRTr,
216			I myThid)
217			ENDIF
218			#endif
219
220			#ifdef ALLOW_MY82
221			IF (useMY82) THEN
222			CALL MY82_CALC_DIFF(
223	jmc	1.2	I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
224	jmc	1.1	U KappaRTr,
225			I myThid)
226			ENDIF
227			#endif
228	jmc	1.11
229	jmc	1.1	#ifdef ALLOW_GGL90
230			IF (useGGL90) THEN
231			CALL GGL90_CALC_DIFF(
232	jmc	1.2	I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
233	jmc	1.1	O KappaRTr,
234			I myThid)
235			ENDIF
236			#endif
237	jmc	1.11
238			C- Apply mask to vertical diffusivity
239			C jmc: don't have the impression that masking is needed
240	jmc	1.1	C but could be removed later if it's the case.
241			c DO j = 1-Oly, sNy+Oly
242			c DO i = 1-Olx, sNx+Olx
243			c KappaRTr(i,j,k) = maskUp(i,j)*KappaRTr(i,j,k)
244			c ENDDO
245			c ENDDO
246
247	jmc	1.3	#endif /* ALLOW_GENERIC_ADVDIFF */
248
249	jmc	1.1	RETURN
250			END