high_res_cube/code-mods/calc_diffusivity.F

C $Header: /u/gcmpack/MITgcm_contrib/high_res_cube/code-mods/calc_diffusivity.F,v 1.1 2007/04/12 02:24:39 dimitri Exp $
C $Name:  $

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

CBOP
C     !ROUTINE: CALC_DIFFUSIVITY
C     !INTERFACE: 
      SUBROUTINE CALC_DIFFUSIVITY( 
     I        bi,bj,iMin,iMax,jMin,jMax,k,
     I        maskUp,
     O        KappaRT,KappaRS,
     I        myThid)

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE CALC_DIFFUSIVITY                               
C     | o Calculate net vertical diffusivity                               
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"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     bi, bj, :: tile indices
C     iMin,   :: Range of points for which calculation is performed.
C     iMax, 
C     jMin, 
C     jMax
C     maskUp  :: land/water mask for Wvel points (above tracer level)
C     myThid  :: Instance number for this innvocation of CALC_DIFFUSIVITY
C     KappaRT :: Net diffusivity for temperature
C     KappaRS :: Net diffusivity for salinity
      INTEGER bi,bj,iMin,iMax,jMin,jMax,K
      _RS maskUp(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL KappaRT(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL KappaRS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      INTEGER myThid

C     !LOCAL VARIABLES:
C     == Local variables ==
C     I, J :: Loop counters
      INTEGER i,j
      _RL KbryanLewis79, KbryanLewisHL, KbryanLewisEQ
CEOP

      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
        KappaRT(i,j) =
     &         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

      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
        KappaRS(i,j) =
     &         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

#ifdef ALLOW_GMREDI
      IF (useGMRedi) THEN 
         CALL GMREDI_CALC_DIFF(
     I        bi,bj,iMin,iMax,jMin,jMax,k,1,
     U        KappaRT,
     I        myThid)
         CALL GMREDI_CALC_DIFF(
     I        bi,bj,iMin,iMax,jMin,jMax,k,1,
     U        KappaRS,
     I        myThid)
      ENDIF
#endif

#ifdef ALLOW_KPP
      IF (useKPP) THEN
         CALL KPP_CALC_DIFF_T(
     I        bi,bj,iMin+1,iMax,jMin+1,jMax,k,1,
     U        KappaRT,
     I        myThid)
         CALL KPP_CALC_DIFF_S(
     I        bi,bj,iMin+1,iMax,jMin+1,jMax,k,1,
     U        KappaRS,
     I        myThid)
      ENDIF
#endif

#ifdef ALLOW_PP81
      IF (usePP81) THEN
         CALL PP81_CALC_DIFF(
     I        bi,bj,iMin+1,iMax,jMin+1,jMax,k,1,
     U        KappaRT,
     I        myThid)
         CALL PP81_CALC_DIFF(
     I        bi,bj,iMin+1,iMax,jMin+1,jMax,k,1,
     U        KappaRS,
     I        myThid)
      ENDIF
#endif

#ifdef ALLOW_MY82
      IF (useMY82) THEN
         CALL MY82_CALC_DIFF(
     I        bi,bj,iMin+1,iMax,jMin+1,jMax,k,1,
     U        KappaRT,
     I        myThid)
         CALL MY82_CALC_DIFF(
     I        bi,bj,iMin+1,iMax,jMin+1,jMax,k,1,
     U        KappaRS,
     I        myThid)
      ENDIF
#endif
      
#ifdef ALLOW_GGL90
      IF (useGGL90) THEN
         CALL GGL90_CALC_DIFF(
     I        bi,bj,iMin+1,iMax,jMin+1,jMax,k,1,
     O        KappaRT,
     I        myThid)
         CALL GGL90_CALC_DIFF(
     I        bi,bj,iMin+1,iMax,jMin+1,jMax,k,1,
     O        KappaRS,
     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.
      DO j = 1-Oly, sNy+Oly
       DO i = 1-Olx, sNx+Olx
        KappaRT(i,j) = maskUp(i,j)*KappaRT(i,j)
        KappaRS(i,j) = maskUp(i,j)*KappaRS(i,j)
       ENDDO
      ENDDO

      RETURN
      END
1	dimitri	1.2	C $Header: /u/gcmpack/MITgcm_contrib/high_res_cube/code-mods/calc_diffusivity.F,v 1.1 2007/04/12 02:24:39 dimitri Exp $
2	dimitri	1.1	C $Name: $
3
4			#include "PACKAGES_CONFIG.h"
5			#include "CPP_OPTIONS.h"
6
7			CBOP
8			C !ROUTINE: CALC_DIFFUSIVITY
9			C !INTERFACE:
10			SUBROUTINE CALC_DIFFUSIVITY(
11			I bi,bj,iMin,iMax,jMin,jMax,k,
12			I maskUp,
13			O KappaRT,KappaRS,
14			I myThid)
15
16			C !DESCRIPTION: \bv
17			C ==========================================================
18			C \| SUBROUTINE CALC_DIFFUSIVITY
19			C \| o Calculate net vertical diffusivity
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
35			C !INPUT/OUTPUT PARAMETERS:
36			C == Routine arguments ==
37			C bi, bj, :: tile indices
38			C iMin, :: Range of points for which calculation is performed.
39			C iMax,
40			C jMin,
41			C jMax
42			C maskUp :: land/water mask for Wvel points (above tracer level)
43			C myThid :: Instance number for this innvocation of CALC_DIFFUSIVITY
44			C KappaRT :: Net diffusivity for temperature
45			C KappaRS :: Net diffusivity for salinity
46			INTEGER bi,bj,iMin,iMax,jMin,jMax,K
47			_RS maskUp(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48			_RL KappaRT(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49			_RL KappaRS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50			INTEGER myThid
51
52			C !LOCAL VARIABLES:
53			C == Local variables ==
54			C I, J :: Loop counters
55			INTEGER i,j
56			_RL KbryanLewis79, KbryanLewisHL, KbryanLewisEQ
57			CEOP
58
59			KbryanLewisHL=diffKrBL79surf+(diffKrBL79deep-diffKrBL79surf)
60			& *( atan( -( rF(k)-diffKrBL79Ho )/diffKrBL79scl )/PI+0.5 _d 0)
61
62			KbryanLewisEQ=diffKrBLEQsurf+(diffKrBLEQdeep-diffKrBLEQsurf)
63			& *( atan( -( rF(k)-diffKrBLEQHo )/diffKrBLEQscl )/PI+0.5 _d 0)
64
65			DO j = 1-Oly, sNy+Oly
66			IF ( abs(YC(i,j,bi,bj)) .GT. KbryanLewisLatTransition ) THEN
67			KbryanLewis79 = KbryanLewisHL
68			ELSE
69			KbryanLewis79 = KbryanLewisHL - (KbryanLewisHL-KbryanLewisEQ) *
70			& (1+cos(YC(i,j,bi,bj)*pi/KbryanLewisLatTransition)) / 2
71			ENDIF
72			DO i = 1-Olx, sNx+Olx
73			KappaRT(i,j) =
74			& IVDConvCount(i,j,k,bi,bj)*ivdc_kappa
75			#if (defined ALLOW_3D_DIFFKR \|\| \
76			(defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_DIFFKR_CONTROL)))
77			& + diffKr(i,j,k,bi,bj)
78			#else
79			& + diffKrNrT(k)
80			#endif
81			& + KbryanLewis79
82			ENDDO
83			ENDDO
84
85			DO j = 1-Oly, sNy+Oly
86			IF ( abs(YC(i,j,bi,bj)) .GT. KbryanLewisLatTransition ) THEN
87			KbryanLewis79 = KbryanLewisHL
88			ELSE
89			KbryanLewis79 = KbryanLewisHL - (KbryanLewisHL-KbryanLewisEQ) *
90			& (1+cos(YC(i,j,bi,bj)*pi/KbryanLewisLatTransition)) / 2
91			ENDIF
92			DO i = 1-Olx, sNx+Olx
93			KappaRS(i,j) =
94			& IVDConvCount(i,j,k,bi,bj)*ivdc_kappa
95			#if (defined ALLOW_3D_DIFFKR \|\| \
96			(defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_DIFFKR_CONTROL)))
97			& + diffKr(i,j,k,bi,bj)
98			#else
99			& + diffKrNrS(k)
100			#endif
101			& + KbryanLewis79
102			ENDDO
103			ENDDO
104
105			#ifdef ALLOW_GMREDI
106			IF (useGMRedi) THEN
107			CALL GMREDI_CALC_DIFF(
108			I bi,bj,iMin,iMax,jMin,jMax,k,1,
109			U KappaRT,
110			I myThid)
111			CALL GMREDI_CALC_DIFF(
112			I bi,bj,iMin,iMax,jMin,jMax,k,1,
113			U KappaRS,
114			I myThid)
115			ENDIF
116			#endif
117
118			#ifdef ALLOW_KPP
119			IF (useKPP) THEN
120			CALL KPP_CALC_DIFF_T(
121			I bi,bj,iMin+1,iMax,jMin+1,jMax,k,1,
122			U KappaRT,
123			I myThid)
124			CALL KPP_CALC_DIFF_S(
125			I bi,bj,iMin+1,iMax,jMin+1,jMax,k,1,
126			U KappaRS,
127			I myThid)
128			ENDIF
129			#endif
130
131			#ifdef ALLOW_PP81
132			IF (usePP81) THEN
133			CALL PP81_CALC_DIFF(
134			I bi,bj,iMin+1,iMax,jMin+1,jMax,k,1,
135			U KappaRT,
136			I myThid)
137			CALL PP81_CALC_DIFF(
138			I bi,bj,iMin+1,iMax,jMin+1,jMax,k,1,
139			U KappaRS,
140			I myThid)
141			ENDIF
142			#endif
143
144			#ifdef ALLOW_MY82
145			IF (useMY82) THEN
146			CALL MY82_CALC_DIFF(
147			I bi,bj,iMin+1,iMax,jMin+1,jMax,k,1,
148			U KappaRT,
149			I myThid)
150			CALL MY82_CALC_DIFF(
151			I bi,bj,iMin+1,iMax,jMin+1,jMax,k,1,
152			U KappaRS,
153			I myThid)
154			ENDIF
155			#endif
156
157			#ifdef ALLOW_GGL90
158			IF (useGGL90) THEN
159			CALL GGL90_CALC_DIFF(
160			I bi,bj,iMin+1,iMax,jMin+1,jMax,k,1,
161			O KappaRT,
162			I myThid)
163			CALL GGL90_CALC_DIFF(
164			I bi,bj,iMin+1,iMax,jMin+1,jMax,k,1,
165			O KappaRS,
166			I myThid)
167			ENDIF
168			#endif
169
170			C- Apply mask to vertical diffusivity
171			C jmc: don't have the impression that masking is needed
172			C but could be removed later if it's the case.
173			DO j = 1-Oly, sNy+Oly
174			DO i = 1-Olx, sNx+Olx
175			KappaRT(i,j) = maskUp(i,j)*KappaRT(i,j)
176			KappaRS(i,j) = maskUp(i,j)*KappaRS(i,j)
177			ENDDO
178			ENDDO
179
180			RETURN
181			END