/[MITgcm]/MITgcm/pkg/ggl90/ggl90_calc.F
ViewVC logotype

Diff of /MITgcm/pkg/ggl90/ggl90_calc.F

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

revision 1.19 by jmc, Thu Aug 19 23:52:37 2010 UTC revision 1.23 by gforget, Wed Aug 8 22:22:42 2012 UTC
# Line 8  C !ROUTINE: GGL90_CALC Line 8  C !ROUTINE: GGL90_CALC
8    
9  C !INTERFACE: ======================================================  C !INTERFACE: ======================================================
10        SUBROUTINE GGL90_CALC(        SUBROUTINE GGL90_CALC(
11       I     bi, bj, myTime, myIter, myThid )       I                 bi, bj, sigmaR, myTime, myIter, myThid )
12    
13    
14  C !DESCRIPTION: \bv  C !DESCRIPTION: \bv
15  C     *==========================================================*  C     *==========================================================*
# Line 41  C !USES: =============================== Line 42  C !USES: ===============================
42    
43  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
44  C Routine arguments  C Routine arguments
45  C     bi, bj :: array indices on which to apply calculations  C     bi, bj :: Current tile indices
46    C     sigmaR :: Vertical gradient of iso-neutral density
47  C     myTime :: Current time in simulation  C     myTime :: Current time in simulation
48  C     myIter :: Current time-step number  C     myIter :: Current time-step number
49  C     myThid :: My Thread Id number  C     myThid :: My Thread Id number
50        INTEGER bi, bj        INTEGER bi, bj
51          _RL     sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
52        _RL     myTime        _RL     myTime
53        INTEGER myIter        INTEGER myIter
54        INTEGER myThid        INTEGER myThid
# Line 72  C     GGL90mixingLength:: mixing length Line 75  C     GGL90mixingLength:: mixing length
75  C         rMixingLength:: inverse of mixing length  C         rMixingLength:: inverse of mixing length
76  C     totalDepth       :: thickness of water column (inverse of recip_Rcol)  C     totalDepth       :: thickness of water column (inverse of recip_Rcol)
77  C     TKEPrandtlNumber :: here, an empirical function of the Richardson number  C     TKEPrandtlNumber :: here, an empirical function of the Richardson number
 C     rhoK, rhoKm1     :: density at layer k and km1 (relative to k)  
78        INTEGER iMin ,iMax ,jMin ,jMax        INTEGER iMin ,iMax ,jMin ,jMax
79        INTEGER i, j, k, kp1, km1, kSurf, kBottom        INTEGER i, j, k, kp1, km1, kSurf, kBottom
80        _RL     explDissFac, implDissFac        _RL     explDissFac, implDissFac
# Line 93  c     _RL     SQRTTKE Line 95  c     _RL     SQRTTKE
95        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
96        _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
       _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL     rhoKm1           (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
98        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
100  C-    tri-diagonal matrix  C-    tri-diagonal matrix
101        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
102        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
103        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
104        INTEGER errCode        INTEGER errCode
105  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
106  C     xA, yA   :: area of lateral faces  C     xA, yA   :: area of lateral faces
# Line 133  C     explicit/implicit timestepping wei Line 133  C     explicit/implicit timestepping wei
133    
134  C     Initialize local fields  C     Initialize local fields
135        DO k = 1, Nr        DO k = 1, Nr
136         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
137          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
138           KappaE(i,j,k)            = 0. _d 0           KappaE(i,j,k)            = 0. _d 0
139           TKEPrandtlNumber(i,j,k)  = 1. _d 0           TKEPrandtlNumber(i,j,k)  = 1. _d 0
140           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
141           GGL90visctmp(i,j,k)      = 0. _d 0           GGL90visctmp(i,j,k)      = 0. _d 0
142    #ifndef SOLVE_DIAGONAL_LOWMEMORY
143             a3d(i,j,k) = 0. _d 0
144             b3d(i,j,k) = 1. _d 0
145             c3d(i,j,k) = 0. _d 0
146    #endif
147          ENDDO          ENDDO
148         ENDDO         ENDDO
149        ENDDO        ENDDO
150        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
151         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
         rhoK(i,j)          = 0. _d 0  
         rhoKm1(i,j)        = 0. _d 0  
152          totalDepth(i,j)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)          totalDepth(i,j)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
153          rMixingLength(i,j,1) = 0. _d 0          rMixingLength(i,j,1) = 0. _d 0
154          mxLength_Dn(i,j,1) = GGL90mixingLengthMin          mxLength_Dn(i,j,1) = GGL90mixingLengthMin
# Line 155  C     Initialize local fields Line 158  C     Initialize local fields
158    
159  C     start k-loop  C     start k-loop
160        DO k = 2, Nr        DO k = 2, Nr
161         km1 = k-1  c      km1 = k-1
162  c      kp1 = MIN(Nr,k+1)  c      kp1 = MIN(Nr,k+1)
        CALL FIND_RHO_2D(  
      I      iMin, iMax, jMin, jMax, k,  
      I      theta(1-OLx,1-OLy,km1,bi,bj), salt(1-OLx,1-OLy,km1,bi,bj),  
      O      rhoKm1,  
      I      km1, bi, bj, myThid )  
   
        CALL FIND_RHO_2D(  
      I      iMin, iMax, jMin, jMax, k,  
      I      theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),  
      O      rhoK,  
      I      k, bi, bj, myThid )  
163         DO j=jMin,jMax         DO j=jMin,jMax
164          DO i=iMin,iMax          DO i=iMin,iMax
165           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
166    
167  C     buoyancy frequency  C     buoyancy frequency
168           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(k)           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
169       &        * ( rhoKm1(i,j) - rhoK(i,j) )*maskC(i,j,k,bi,bj)       &                  * sigmaR(i,j,k)
170  cC     vertical shear term (dU/dz)^2+(dV/dz)^2  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
171  c         tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)  c         tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
172  c     &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )  c     &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
# Line 195  C     mixing length Line 187  C     mixing length
187         ENDDO         ENDDO
188        ENDDO        ENDDO
189    
190    C- ensure mixing between first and second level
191          IF (mxlSurfFlag) THEN
192           DO j=jMin,jMax
193            DO i=iMin,iMax
194             GGL90mixingLength(i,j,2)=drF(1)
195            ENDDO
196           ENDDO
197          ENDIF
198    
199  C- Impose upper and lower bound for mixing length  C- Impose upper and lower bound for mixing length
200        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) THEN
201  C-  
202         DO k=2,Nr         DO k=2,Nr
203          DO j=jMin,jMax          DO j=jMin,jMax
204           DO i=iMin,iMax           DO i=iMin,iMax
# Line 219  C- Line 220  C-
220         ENDDO         ENDDO
221    
222        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
223  C-  
224         DO k=2,Nr         DO k=2,Nr
225          DO j=jMin,jMax          DO j=jMin,jMax
226           DO i=iMin,iMax           DO i=iMin,iMax
# Line 242  c         MaxLength=MAX(MaxLength,20. _d Line 243  c         MaxLength=MAX(MaxLength,20. _d
243         ENDDO         ENDDO
244    
245        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
246  C-  
 cgf ensure mixing between first and second level  
 c      DO j=jMin,jMax  
 c        DO i=iMin,iMax  
 c         GGL90mixingLength(i,j,2)=drF(1)  
 c        ENDDO  
 c      ENDDO  
 cgf  
247         DO k=2,Nr         DO k=2,Nr
248          DO j=jMin,jMax          DO j=jMin,jMax
249           DO i=iMin,iMax           DO i=iMin,iMax
# Line 284  cgf Line 278  cgf
278         ENDDO         ENDDO
279    
280        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
281  C-  
282         DO k=2,Nr         DO k=2,Nr
283          DO j=jMin,jMax          DO j=jMin,jMax
284           DO i=iMin,iMax           DO i=iMin,iMax
# Line 343  c     ENDDO Line 337  c     ENDDO
337  C     horizontal diffusion of TKE (requires an exchange in  C     horizontal diffusion of TKE (requires an exchange in
338  C      do_fields_blocking_exchanges)  C      do_fields_blocking_exchanges)
339  C     common factors  C     common factors
340          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
341           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
342            xA(i,j) = _dyG(i,j,bi,bj)            xA(i,j) = _dyG(i,j,bi,bj)
343       &         *drF(k)*_hFacW(i,j,k,bi,bj)       &         *drF(k)*_hFacW(i,j,k,bi,bj)
344            yA(i,j) = _dxG(i,j,bi,bj)            yA(i,j) = _dxG(i,j,bi,bj)
# Line 353  C     common factors Line 347  C     common factors
347          ENDDO          ENDDO
348  C     Compute diffusive fluxes  C     Compute diffusive fluxes
349  C     ... across x-faces  C     ... across x-faces
350          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
351           dfx(1-Olx,j)=0. _d 0           dfx(1-OLx,j)=0. _d 0
352           DO i=1-Olx+1,sNx+Olx           DO i=1-OLx+1,sNx+OLx
353            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
354       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
355       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))
# Line 363  C     ... across x-faces Line 357  C     ... across x-faces
357           ENDDO           ENDDO
358          ENDDO          ENDDO
359  C     ... across y-faces  C     ... across y-faces
360          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
361           dfy(i,1-Oly)=0. _d 0           dfy(i,1-OLy)=0. _d 0
362          ENDDO          ENDDO
363          DO j=1-Oly+1,sNy+Oly          DO j=1-OLy+1,sNy+OLy
364           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
365            dfy(i,j) = -GGL90diffTKEh*yA(i,j)            dfy(i,j) = -GGL90diffTKEh*yA(i,j)
366       &      *_recip_dyC(i,j,bi,bj)       &      *_recip_dyC(i,j,bi,bj)
367       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))
# Line 377  C     ... across y-faces Line 371  C     ... across y-faces
371           ENDDO           ENDDO
372          ENDDO          ENDDO
373  C     Compute divergence of fluxes  C     Compute divergence of fluxes
374          DO j=1-Oly,sNy+Oly-1          DO j=1-OLy,sNy+OLy-1
375           DO i=1-Olx,sNx+Olx-1           DO i=1-OLx,sNx+OLx-1
376            gTKE(i,j) =            gTKE(i,j) =
377       &    -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)       &    -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)
378       &         *( (dfx(i+1,j)-dfx(i,j))       &         *( (dfx(i+1,j)-dfx(i,j))
# Line 454  C     set up matrix Line 448  C     set up matrix
448  C--   Lower diagonal  C--   Lower diagonal
449        DO j=jMin,jMax        DO j=jMin,jMax
450         DO i=iMin,iMax         DO i=iMin,iMax
451           a(i,j,1) = 0. _d 0           a3d(i,j,1) = 0. _d 0
452         ENDDO         ENDDO
453        ENDDO        ENDDO
454        DO k=2,Nr        DO k=2,Nr
# Line 464  C--   Lower diagonal Line 458  C--   Lower diagonal
458  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
459  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
460  C-    No need for maskC(k-1) with recip_hFacC(k-1)  C-    No need for maskC(k-1) with recip_hFacC(k-1)
461           a(i,j,k) = -deltaTggl90           a3d(i,j,k) = -deltaTggl90
462       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
463       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
464       &        *recip_drC(k)*maskC(i,j,k,bi,bj)       &        *recip_drC(k)*maskC(i,j,k,bi,bj)
# Line 474  C-    No need for maskC(k-1) with recip_ Line 468  C-    No need for maskC(k-1) with recip_
468  C--   Upper diagonal  C--   Upper diagonal
469        DO j=jMin,jMax        DO j=jMin,jMax
470         DO i=iMin,iMax         DO i=iMin,iMax
471           c(i,j,1)  = 0. _d 0           c3d(i,j,1)  = 0. _d 0
472         ENDDO         ENDDO
473        ENDDO        ENDDO
474        DO k=2,Nr        DO k=2,Nr
# Line 484  C--   Upper diagonal Line 478  C--   Upper diagonal
478  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
479  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
480  C-    No need for maskC(k) with recip_hFacC(k)  C-    No need for maskC(k) with recip_hFacC(k)
481            c(i,j,k) = -deltaTggl90            c3d(i,j,k) = -deltaTggl90
482       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
483       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
484       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
# Line 496  C--   Center diagonal Line 490  C--   Center diagonal
490         km1 = MAX(k-1,1)         km1 = MAX(k-1,1)
491         DO j=jMin,jMax         DO j=jMin,jMax
492          DO i=iMin,iMax          DO i=iMin,iMax
493            b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)            b3d(i,j,k) = 1. _d 0 - c3d(i,j,k) - a3d(i,j,k)
494       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
495       &        * rMixingLength(i,j,k)       &        * rMixingLength(i,j,k)
496       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
# Line 520  C     Dirichlet surface boundary conditi Line 514  C     Dirichlet surface boundary conditi
514          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
515       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
516          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
517       &               - a(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
518          a(i,j,kp1) = 0. _d 0          a3d(i,j,kp1) = 0. _d 0
519  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
520          kBottom   = MAX(kLowC(i,j,bi,bj),1)          kBottom   = MAX(kLowC(i,j,bi,bj),1)
521          GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)          GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
522       &                              - GGL90TKEbottom*c(i,j,kBottom)       &                              - GGL90TKEbottom*c3d(i,j,kBottom)
523          c(i,j,kBottom) = 0. _d 0          c3d(i,j,kBottom) = 0. _d 0
524         ENDDO         ENDDO
525        ENDDO        ENDDO
526    
527  C     solve tri-diagonal system  C     solve tri-diagonal system
528        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
529       I                        a, b, c,       I                        a3d, b3d, c3d,
530       U                        GGL90TKE,       U                        GGL90TKE,
531       O                        errCode,       O                        errCode,
532       I                        bi, bj, myThid )       I                        bi, bj, myThid )
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22