/[MITgcm]/MITgcm/model/src/calc_3d_diffusivity.F
ViewVC logotype

Diff of /MITgcm/model/src/calc_3d_diffusivity.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.7 by dimitri, Wed Apr 11 00:00:47 2007 UTC revision 1.19 by jmc, Thu Nov 2 18:00:52 2017 UTC
# Line 6  C $Name$ Line 6  C $Name$
6    
7  CBOP  CBOP
8  C     !ROUTINE: CALC_3D_DIFFUSIVITY  C     !ROUTINE: CALC_3D_DIFFUSIVITY
9  C     !INTERFACE:  C     !INTERFACE:
10        SUBROUTINE CALC_3D_DIFFUSIVITY(        SUBROUTINE CALC_3D_DIFFUSIVITY(
11       I        bi,bj,iMin,iMax,jMin,jMax,       I        bi, bj, iMin,iMax,jMin,jMax,
12       I        trIdentity, trUseGMRedi, trUseKPP,       I        trIdentity, trUseGMRedi, trUseKPP,
13       O        KappaRTr,       O        KappaRTr,
14       I        myThid)       I        myThid )
15    
16  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
17  C     *==========================================================*  C     *==========================================================*
# Line 36  C     == GLobal variables == Line 36  C     == GLobal variables ==
36  #endif  #endif
37  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
38  #include "PTRACERS_SIZE.h"  #include "PTRACERS_SIZE.h"
39  #include "PTRACERS.h"  #include "PTRACERS_PARAMS.h"
40    #endif
41    #ifdef ALLOW_LONGSTEP
42    #include "LONGSTEP.h"
43  #endif  #endif
44    
45  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
# Line 52  C     KappaRTr   :: Net diffusivity for Line 55  C     KappaRTr   :: Net diffusivity for
55        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
56        INTEGER trIdentity        INTEGER trIdentity
57        LOGICAL trUseGMRedi, trUseKPP        LOGICAL trUseGMRedi, trUseKPP
58        _RL KappaRTr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRTr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
59        INTEGER myThid        INTEGER myThid
60    
61  #ifdef ALLOW_GENERIC_ADVDIFF  #ifdef ALLOW_GENERIC_ADVDIFF
# Line 62  C     i, j, k    :: Loop counters Line 65  C     i, j, k    :: Loop counters
65  C     iTr        :: passive tracer index  C     iTr        :: passive tracer index
66  C     msgBuf     :: message buffer  C     msgBuf     :: message buffer
67        INTEGER i,j,k        INTEGER i,j,k
68        _RL KbryanLewis79, KbryanLewisHL, KbryanLewisEQ        _RL KbryanLewis79
69    #ifdef ALLOW_BL79_LAT_VARY
70          _RL KbryanLewisEQ
71    #endif
72        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
73  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
74        INTEGER iTr        INTEGER iTr
75  #endif  #endif
76    #ifndef EXCLUDE_PCELL_MIX_CODE
77          INTEGER km, mixSurf, mixBott
78          _RL pC_kFac
79          _RL tmpFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80    #endif
81  CEOP  CEOP
82    
83        IF ( trIdentity.EQ.GAD_TEMPERATURE) THEN        IF ( .NOT. trUseKPP ) THEN
   
84         DO k = 1,Nr         DO k = 1,Nr
85          KbryanLewisHL=diffKrBL79surf+(diffKrBL79deep-diffKrBL79surf)          KbryanLewis79=diffKrBL79surf+(diffKrBL79deep-diffKrBL79surf)
86       &   *( atan( -( rF(k)-diffKrBL79Ho )/diffKrBL79scl )/PI+0.5 _d 0)       &       *(atan(-(rF(k)-diffKrBL79Ho)/diffKrBL79scl)/PI+0.5 _d 0)
87    #ifdef ALLOW_BL79_LAT_VARY
88          KbryanLewisEQ=diffKrBLEQsurf+(diffKrBLEQdeep-diffKrBLEQsurf)          KbryanLewisEQ=diffKrBLEQsurf+(diffKrBLEQdeep-diffKrBLEQsurf)
89       &   *( atan( -( rF(k)-diffKrBLEQHo )/diffKrBLEQscl )/PI+0.5 _d 0)       &       *(atan(-(rF(k)-diffKrBLEQHo)/diffKrBLEQscl)/PI+0.5 _d 0)
90          DO j = 1-Oly,sNy+Oly  #endif
91           IF ( abs(YC(i,j,bi,bj)) .GT. KbryanLewisLatTransition ) THEN          DO j = 1-OLy,sNy+OLy
92            KbryanLewis79 = KbryanLewisHL           DO i = 1-OLx,sNx+OLx
93           ELSE  #ifdef ALLOW_LONGSTEP
94            KbryanLewis79 = KbryanLewisHL-(KbryanLewisHL-KbryanLewisEQ)*            IF ( trIdentity .GE. GAD_TR1) THEN
95       &         (1+cos(YC(i,j,bi,bj)*pi/KbryanLewisLatTransition)) / 2             KappaRTr(i,j,k) =
96           ENDIF       &         LS_IVDConvCount(i,j,k,bi,bj)*ivdc_kappa
97           DO i = 1-Olx,sNx+Olx       &         + KbryanLewis79
98            KappaRTr(i,j,k) =  #ifdef ALLOW_BL79_LAT_VARY
99       &         IVDConvCount(i,j,k,bi,bj)*ivdc_kappa       &         + (KbryanLewisEQ-KbryanLewis79)*BL79LatArray(i,j,bi,bj)
100  #if (defined ALLOW_3D_DIFFKR || \  #endif
101       (defined (ALLOW_AUTODIFF_TAMC) && defined (ALLOW_DIFFKR_CONTROL)))            ELSE
      &       + diffKr(i,j,k,bi,bj)  
102  #else  #else
103       &       + diffKrNrT(k)            IF ( .TRUE. ) THEN
104    #endif /* ALLOW_LONGSTEP */
105               KappaRTr(i,j,k) =
106         &         IVDConvCount(i,j,k,bi,bj)*ivdc_kappa
107         &         + KbryanLewis79
108    #ifdef ALLOW_BL79_LAT_VARY
109         &         + (KbryanLewisEQ-KbryanLewis79)*BL79LatArray(i,j,bi,bj)
110  #endif  #endif
111       &       + KbryanLewis79            ENDIF
112           ENDDO           ENDDO
113          ENDDO          ENDDO
114         ENDDO         ENDDO
115           IF ( trIdentity.EQ.GAD_TEMPERATURE ) THEN
116        ELSEIF ( trIdentity.EQ.GAD_SALINITY) THEN          DO k = 1,Nr
117             DO j = 1-OLy,sNy+OLy
118         DO k = 1,Nr            DO i = 1-OLx,sNx+OLx
119          KbryanLewisHL=diffKrBL79surf+(diffKrBL79deep-diffKrBL79surf)             KappaRTr(i,j,k) = KappaRTr(i,j,k)
120       &   *( atan( -( rF(k)-diffKrBL79Ho )/diffKrBL79scl )/PI+0.5 _d 0)  #ifdef ALLOW_3D_DIFFKR
121          KbryanLewisEQ=diffKrBLEQsurf+(diffKrBLEQdeep-diffKrBLEQsurf)       &          + diffKr(i,j,k,bi,bj)
      &   *( 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)  
122  #else  #else
123       &       + diffKrNrS(k)       &          + diffKrNrT(k)
124  #endif  #endif
125       &       + KbryanLewis79            ENDDO
126             ENDDO
127            ENDDO
128           ELSEIF ( trIdentity.EQ.GAD_SALINITY) THEN
129            DO k = 1,Nr
130             DO j = 1-OLy, sNy+OLy
131              DO i = 1-OLx, sNx+OLx
132               KappaRTr(i,j,k) = KappaRTr(i,j,k)
133    #ifdef ALLOW_3D_DIFFKR
134         &          + diffKr(i,j,k,bi,bj)
135    #else
136         &          + diffKrNrS(k)
137    #endif
138              ENDDO
139           ENDDO           ENDDO
140          ENDDO          ENDDO
        ENDDO  
   
141  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
142        ELSEIF ( trIdentity.GE.GAD_TR1         ELSEIF ( trIdentity.GE.GAD_TR1) THEN
      &   .AND. trIdentity.LT.GAD_TR1+PTRACERS_numInUse) THEN  
143    
144         iTr = trIdentity - GAD_TR1 + 1          iTr = trIdentity - GAD_TR1 + 1
145         DO k = 1,Nr          DO k = 1,Nr
146          KbryanLewisHL=diffKrBL79surf+(diffKrBL79deep-diffKrBL79surf)           DO j = 1-OLy, sNy+OLy
147       &   *( atan( -( rF(k)-diffKrBL79Ho )/diffKrBL79scl )/PI+0.5 _d 0)            DO i = 1-OLx, sNx+OLx
148          KbryanLewisEQ=diffKrBLEQsurf+(diffKrBLEQdeep-diffKrBLEQsurf)             KappaRTr(i,j,k) = KappaRTr(i,j,k)
149       &   *( atan( -( rF(k)-diffKrBLEQHo )/diffKrBLEQscl )/PI+0.5 _d 0)  #ifdef ALLOW_3D_DIFFKR
150          DO j = 1-Oly, sNy+Oly       &          + diffKr(i,j,k,bi,bj)
          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)  
151  #else  #else
152       &       + PTRACERS_diffKrNr(k,iTr)       &          + PTRACERS_diffKrNr(k,iTr)
153  #endif  #endif
154       &       + KbryanLewis79            ENDDO
155           ENDDO           ENDDO
156          ENDDO          ENDDO
        ENDDO  
157  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
158           ELSE
       ELSE  
159          WRITE(msgBuf,'(A,I4)')          WRITE(msgBuf,'(A,I4)')
160       &      ' CALC_3D_DIFFUSIVITY: Invalid tracer Id: ',trIdentity       &       ' CALC_3D_DIFFUSIVITY: Invalid tracer Id: ',trIdentity
161          CALL PRINT_ERROR(msgBuf, myThid)          CALL PRINT_ERROR(msgBuf, myThid)
162          STOP 'ABNORMAL END: S/R CALC_3D_DIFFUSIVITY'          STOP 'ABNORMAL END: S/R CALC_3D_DIFFUSIVITY'
163           ENDIF
164        ENDIF        ENDIF
165    
166  C--   Add physical pacakge contributions:  C--   Add physical pacakge contributions:
167    
 #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  
   
168  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
169        IF (trUseKPP) THEN        IF (trUseKPP) THEN
170    C--   Set vertical diffusivity contribution from KPP
171         IF (trIdentity.EQ.GAD_TEMPERATURE) THEN         IF (trIdentity.EQ.GAD_TEMPERATURE) THEN
172           CALL KPP_CALC_DIFF_T(           CALL KPP_CALC_DIFF_T(
173       I        bi,bj,iMin,iMax,jMin,jMax,0,Nr,       I        bi,bj,iMin,iMax,jMin,jMax,0,Nr,
174       U        KappaRTr,       O        KappaRTr,
175       I        myThid)       I        myThid)
176         ELSE         ELSEIF (trIdentity.EQ.GAD_SALINITY) THEN
177           CALL KPP_CALC_DIFF_S(           CALL KPP_CALC_DIFF_S(
178       I        bi,bj,iMin,iMax,jMin,jMax,0,Nr,       I        bi,bj,iMin,iMax,jMin,jMax,0,Nr,
179       U        KappaRTr,       O        KappaRTr,
180       I        myThid)       I        myThid)
181    #ifdef ALLOW_PTRACERS
182           ELSEIF ( trIdentity.GE.GAD_TR1) THEN
183             iTr = trIdentity - GAD_TR1 + 1
184             CALL KPP_CALC_DIFF_Ptr(
185         I        bi,bj,iMin,iMax,jMin,jMax,0,Nr,
186         O        KappaRTr,
187         I        iTr, myThid )
188    #endif /* ALLOW_PTRACERS */
189           ELSE
190            WRITE(msgBuf,'(A,I4)')
191         &       ' CALC_3D_DIFFUSIVITY: Invalid tracer Id: ',trIdentity
192            CALL PRINT_ERROR( msgBuf, myThid )
193            STOP 'ABNORMAL END: S/R CALC_3D_DIFFUSIVITY'
194         ENDIF         ENDIF
195        ENDIF        ENDIF
196    #endif /* ALLOW_KPP */
197    
198    #ifdef ALLOW_GMREDI
199          IF (trUseGMRedi) THEN
200             CALL GMREDI_CALC_DIFF(
201         I        bi,bj,iMin,iMax,jMin,jMax,0,Nr,
202         U        KappaRTr,
203         I        trIdentity,myThid)
204          ENDIF
205  #endif  #endif
206    
207  #ifdef ALLOW_PP81  #ifdef ALLOW_PP81
# Line 200  C--   Add physical pacakge contributions Line 213  C--   Add physical pacakge contributions
213        ENDIF        ENDIF
214  #endif  #endif
215    
216    #ifdef ALLOW_KL10
217          IF (useKL10) THEN
218             CALL KL10_CALC_DIFF(
219         I        bi,bj,iMin,iMax,jMin,jMax,0,Nr,
220         U        KappaRTr,
221         I        myThid)
222          ENDIF
223    #endif
224    
225  #ifdef ALLOW_MY82  #ifdef ALLOW_MY82
226        IF (useMY82) THEN        IF (useMY82) THEN
227           CALL MY82_CALC_DIFF(           CALL MY82_CALC_DIFF(
# Line 208  C--   Add physical pacakge contributions Line 230  C--   Add physical pacakge contributions
230       I        myThid)       I        myThid)
231        ENDIF        ENDIF
232  #endif  #endif
233          
234  #ifdef ALLOW_GGL90  #ifdef ALLOW_GGL90
235        IF (useGGL90) THEN        IF (useGGL90) THEN
236           CALL GGL90_CALC_DIFF(           CALL GGL90_CALC_DIFF(
# Line 217  C--   Add physical pacakge contributions Line 239  C--   Add physical pacakge contributions
239       I        myThid)       I        myThid)
240        ENDIF        ENDIF
241  #endif  #endif
242          
243  C-    Apply mask to vertical diffusivity  #ifndef EXCLUDE_PCELL_MIX_CODE
244  C jmc: don't have the impression that masking is needed        IF ( interDiffKr_pCell ) THEN
245  C      but could be removed later if it's the case.  C--   This is a hack: alter vertical diffusivity (instead of changing many S/R)
246  c     DO j = 1-Oly, sNy+Oly  C     in order to account for missing hFac in diffusion term
247  c      DO i = 1-Olx, sNx+Olx         DO k = 2,Nr
248             km = k - 1
249    C-    account for true distance (including hFac) in vertical gradient
250             DO j = 2-OLy, sNy+OLy
251              DO i = 2-OLx, sNx+OLx
252               IF ( k.GT.kSurfC(i,j,bi,bj) .AND.
253         &          k.LE.kLowC(i,j,bi,bj) ) THEN
254                 KappaRTr(i,j,k) = KappaRTr(i,j,k)
255         &                *twoRL/(hFacC(i,j,km,bi,bj)+hFacC(i,j,k,bi,bj))
256               ENDIF
257              ENDDO
258             ENDDO
259           ENDDO
260          ENDIF
261    
262          IF ( pCellMix_select.GT.0 ) THEN
263    C--   This is a hack: alter vertical diffusivity (instead of changing many S/R)
264    C     in order to to increase mixing for too thin surface/bottom partial cell
265           mixSurf = pCellMix_select/10
266           mixBott = MOD(pCellMix_select,10)
267           DO k = 2,Nr
268             km = k - 1
269             pC_kFac = 1.
270             IF ( pCellMix_delR.LT.drF(k) )
271         &     pC_kFac = pCellMix_delR*recip_drF(k)
272    
273    C-    Increase KappaRTr above bottom level:
274             IF ( mixBott.GE.1 ) THEN
275              DO j = 2-OLy, sNy+OLy
276               DO i = 2-OLx, sNx+OLx
277                 tmpFac(i,j) = 0. _d 0
278                 IF ( k.EQ.kLowC(i,j,bi,bj) .AND.
279         &            k.GT.kSurfC(i,j,bi,bj) ) THEN
280                   tmpFac(i,j) = pC_kFac*_recip_hFacC(i,j,k,bi,bj)
281                 ENDIF
282               ENDDO
283              ENDDO
284             ENDIF
285             IF ( mixBott.EQ.2 ) THEN
286              DO j = 2-OLy, sNy+OLy
287               DO i = 2-OLx, sNx+OLx
288                 tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)
289               ENDDO
290              ENDDO
291             ELSEIF ( mixBott.EQ.3 ) THEN
292              DO j = 2-OLy, sNy+OLy
293               DO i = 2-OLx, sNx+OLx
294                 tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)*tmpFac(i,j)
295               ENDDO
296              ENDDO
297             ELSEIF ( mixBott.EQ.4 ) THEN
298              DO j = 2-OLy, sNy+OLy
299               DO i = 2-OLx, sNx+OLx
300                 tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)
301         &                   * tmpFac(i,j)*tmpFac(i,j)
302               ENDDO
303              ENDDO
304             ENDIF
305             IF ( mixBott.GE.1 ) THEN
306    C-    increase mixing above bottom (by ~(1/hFac)^mixBott) if too thin p-cell
307              DO j = 2-OLy, sNy+OLy
308               DO i = 2-OLx, sNx+OLx
309                 tmpFac(i,j) = MIN( tmpFac(i,j), pCellMix_maxFac )
310                 KappaRTr(i,j,k) = MAX( KappaRTr(i,j,k),
311         &                              pCellMix_diffKr(k)*tmpFac(i,j) )
312               ENDDO
313              ENDDO
314             ENDIF
315    
316             pC_kFac = 1.
317             IF ( pCellMix_delR.LT.drF(km) )
318         &     pC_kFac = pCellMix_delR*recip_drF(km)
319    
320    C-    Increase KappaRTr below surface level:
321             IF ( mixSurf.GE.1 ) THEN
322              DO j = 2-OLy, sNy+OLy
323               DO i = 2-OLx, sNx+OLx
324                 tmpFac(i,j) = 0. _d 0
325                 IF ( km.EQ.kSurfC(i,j,bi,bj) .AND.
326         &            km.LT.kLowC(i,j,bi,bj) ) THEN
327                   tmpFac(i,j) = pC_kFac*_recip_hFacW(i,j,km,bi,bj)
328                 ENDIF
329               ENDDO
330              ENDDO
331             ENDIF
332             IF ( mixSurf.EQ.2 ) THEN
333              DO j = 2-OLy, sNy+OLy
334               DO i = 2-OLx, sNx+OLx
335                 tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)
336               ENDDO
337              ENDDO
338             ELSEIF ( mixSurf.EQ.3 ) THEN
339              DO j = 2-OLy, sNy+OLy
340               DO i = 2-OLx, sNx+OLx
341                 tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)*tmpFac(i,j)
342               ENDDO
343              ENDDO
344             ELSEIF ( mixSurf.EQ.4 ) THEN
345              DO j = 2-OLy, sNy+OLy
346               DO i = 2-OLx, sNx+OLx
347                 tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)
348         &                   * tmpFac(i,j)*tmpFac(i,j)
349               ENDDO
350              ENDDO
351             ENDIF
352             IF ( mixSurf.GE.1 ) THEN
353    C-    increase mixing below surface (by ~(1/hFac)^mixSurf) if too thin p-cell
354              DO j = 2-OLy, sNy+OLy
355               DO i = 2-OLx, sNx+OLx
356                 tmpFac(i,j) = MIN( tmpFac(i,j), pCellMix_maxFac )
357                 KappaRTr(i,j,k) = MAX( KappaRTr(i,j,k),
358         &                              pCellMix_diffKr(k)*tmpFac(i,j) )
359               ENDDO
360              ENDDO
361             ENDIF
362    
363    C--   end of k loop
364           ENDDO
365          ENDIF
366    #endif /* ndef EXCLUDE_PCELL_MIX_CODE */
367    
368    C-    Apply mask to vertical diffusivity
369    C jmc: do not have the impression that masking is needed
370    C      but could be removed later if it is the case.
371    c     DO j = 1-OLy, sNy+OLy
372    c      DO i = 1-OLx, sNx+OLx
373  c       KappaRTr(i,j,k) = maskUp(i,j)*KappaRTr(i,j,k)  c       KappaRTr(i,j,k) = maskUp(i,j)*KappaRTr(i,j,k)
374  c      ENDDO  c      ENDDO
375  c     ENDDO  c     ENDDO
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.22