/[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.26 by jmc, Thu Aug 14 16:42:35 2014 UTC revision 1.28 by mlosch, Sat Feb 21 01:46:39 2015 UTC
# Line 51  C     myThid :: My Thread Id number Line 51  C     myThid :: My Thread Id number
51        _RL     myTime        _RL     myTime
52        INTEGER myIter        INTEGER myIter
53        INTEGER myThid        INTEGER myThid
 CEOP  
54    
55  #ifdef ALLOW_GGL90  #ifdef ALLOW_GGL90
56    
# Line 86  c     _RL     Nsquare Line 85  c     _RL     Nsquare
85  c     _RL     SQRTTKE  c     _RL     SQRTTKE
86        _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
87        _RL     RiNumber        _RL     RiNumber
88    #ifdef ALLOW_GGL90_IDEMIX
89          _RL     IDEMIX_RiNumber
90    #endif
91        _RL     TKEdissipation        _RL     TKEdissipation
92        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
93        _RL     MaxLength, tmpmlx, tmpVisc        _RL     MaxLength, tmpmlx, tmpVisc
# Line 96  c     _RL     SQRTTKE Line 98  c     _RL     SQRTTKE
98        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
99        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
101    #ifdef ALLOW_DIAGNOSTICS
102          _RL     surf_flx_tke     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103    #endif /* ALLOW_DIAGNOSTICS */
104  C-    tri-diagonal matrix  C-    tri-diagonal matrix
105        _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
106        _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 115  C     gTKE     :: right hand side of dif Line 120  C     gTKE     :: right hand side of dif
120  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
121  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
122        _RL p4, p8, p16        _RL p4, p8, p16
123    CEOP
124        p4=0.25 _d 0        p4=0.25 _d 0
125        p8=0.125 _d 0        p8=0.125 _d 0
126        p16=0.0625 _d 0        p16=0.0625 _d 0
# Line 136  C     Initialize local fields Line 142  C     Initialize local fields
142        DO k = 1, Nr        DO k = 1, Nr
143         DO j=1-OLy,sNy+OLy         DO j=1-OLy,sNy+OLy
144          DO i=1-OLx,sNx+OLx          DO i=1-OLx,sNx+OLx
145             rMixingLength(i,j,k)     = 0. _d 0
146             mxLength_Dn(i,j,k)       = 0. _d 0
147             GGL90visctmp(i,j,k)      = 0. _d 0
148           KappaE(i,j,k)            = 0. _d 0           KappaE(i,j,k)            = 0. _d 0
149           TKEPrandtlNumber(i,j,k)  = 1. _d 0           TKEPrandtlNumber(i,j,k)  = 1. _d 0
150           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
# Line 145  C     Initialize local fields Line 154  C     Initialize local fields
154           b3d(i,j,k) = 1. _d 0           b3d(i,j,k) = 1. _d 0
155           c3d(i,j,k) = 0. _d 0           c3d(i,j,k) = 0. _d 0
156  #endif  #endif
157             Nsquare(i,j,k) = 0. _d 0
158             SQRTTKE(i,j,k) = 0. _d 0
159          ENDDO          ENDDO
160         ENDDO         ENDDO
161        ENDDO        ENDDO
# Line 154  C     Initialize local fields Line 165  C     Initialize local fields
165          rMixingLength(i,j,1) = 0. _d 0          rMixingLength(i,j,1) = 0. _d 0
166          mxLength_Dn(i,j,1) = GGL90mixingLengthMin          mxLength_Dn(i,j,1) = GGL90mixingLengthMin
167          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
168    #ifdef ALLOW_GGL90_HORIZDIFF
169            xA(i,j)  = 0. _d 0
170            yA(i,j)  = 0. _d 0
171            dfx(i,j) = 0. _d 0
172            dfy(i,j) = 0. _d 0
173            gTKE(i,j) = 0. _d 0
174    #endif /* ALLOW_GGL90_HORIZDIFF */
175         ENDDO         ENDDO
176        ENDDO        ENDDO
177    
178    #ifdef ALLOW_GGL90_IDEMIX
179          IF ( useIDEMIX) CALL GGL90_IDEMIX(
180         &   bi, bj, sigmaR, myTime, myIter, myThid )
181    #endif /* ALLOW_GGL90_IDEMIX */
182    
183  C     start k-loop  C     start k-loop
184        DO k = 2, Nr        DO k = 2, Nr
185  c      km1 = k-1  c      km1 = k-1
# Line 378  C     ... across y-faces Line 401  C     ... across y-faces
401  C     Compute divergence of fluxes  C     Compute divergence of fluxes
402          DO j=1-OLy,sNy+OLy-1          DO j=1-OLy,sNy+OLy-1
403           DO i=1-OLx,sNx+OLx-1           DO i=1-OLx,sNx+OLx-1
404            hFac = min(.5 _d 0,_hFacC(i,j,k-1,bi,bj) ) +  #ifdef ALLOW_GGL90_IDEMIX
405       &          min(.5 _d 0,_hFacC(i,j,k  ,bi,bj) )            gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)
406         &         *recip_hFacI(i,j,k,bi,bj)
407    #else
408              hFac = MIN(.5 _d 0,_hFacC(i,j,k-1,bi,bj) ) +
409         &           MIN(.5 _d 0,_hFacC(i,j,k  ,bi,bj) )
410            gTKE(i,j) = 0.0            gTKE(i,j) = 0.0
411            if ( hFac .ne. 0.0 )            IF ( hFac .ne. 0.0 )
412       &      gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)/hFac       &      gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)/hFac
413    #endif
414       &         *((dfx(i+1,j)-dfx(i,j))       &         *((dfx(i+1,j)-dfx(i,j))
415       &          +(dfy(i,j+1)-dfy(i,j)) )       &          +(dfy(i,j+1)-dfy(i,j)) )
416           ENDDO           ENDDO
# Line 397  C     vertical shear term (dU/dz)^2+(dV/ Line 425  C     vertical shear term (dU/dz)^2+(dV/
425           tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)           tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
426       &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )       &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
427       &        *recip_drC(k)       &        *recip_drC(k)
428    #ifdef ALLOW_GGL90_IDEMIX
429         &        *recip_hFacI(i,j,k,bi,bj)
430    #endif
431           tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)           tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
432       &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )       &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
433       &        *recip_drC(k)       &        *recip_drC(k)
434    #ifdef ALLOW_GGL90_IDEMIX
435         &        *recip_hFacI(i,j,k,bi,bj)
436    #endif
437           verticalShear = tempU*tempU + tempV*tempV           verticalShear = tempU*tempU + tempV*tempV
          RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)  
 C     compute Prandtl number (always greater than 0)  
          prTemp = 1. _d 0  
          IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber  
          TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)  
 c         TKEPrandtlNumber(i,j,k) = 1. _d 0  
438    
439  C     viscosity and diffusivity  C     viscosity and diffusivity
440           KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)           KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
# Line 415  C     viscosity and diffusivity Line 443  C     viscosity and diffusivity
443  c        note: storing GGL90visctmp like this, and using it later to compute  c        note: storing GGL90visctmp like this, and using it later to compute
444  c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)  c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
445           KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)           KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)
446    
447    C     compute Prandtl number (always greater than 0)
448             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
449    #ifdef ALLOW_GGL90_IDEMIX
450    CML         IDEMIX_RiNumber = 1./GGL90eps
451             IDEMIX_RiNumber = MAX( KappaM*Nsquare(i,j,k), 0. _d 0)/
452         &    (GGL90eps+IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2)
453             prTemp         = MIN(5.*RiNumber, 6.6*IDEMIX_RiNumber)
454    #else
455             prTemp = 1. _d 0        
456             IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
457    #endif /* ALLOW_GGL90_IDEMIX */
458             TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
459             TKEPrandtlNumber(i,j,k) = MAX( 1. _d 0,TKEPrandtlNumber(i,j,k))
460    
461    c        diffusivity
462           KappaH = KappaM/TKEPrandtlNumber(i,j,k)           KappaH = KappaM/TKEPrandtlNumber(i,j,k)
463           KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)           KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)
464    
# Line 428  C     partial update with sum of explici Line 472  C     partial update with sum of explici
472       &        + KappaM*verticalShear       &        + KappaM*verticalShear
473       &        - KappaH*Nsquare(i,j,k)       &        - KappaH*Nsquare(i,j,k)
474       &        - TKEdissipation       &        - TKEdissipation
475    #ifdef ALLOW_GGL90_IDEMIX
476         &        + IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2
477    #endif
478       &        )       &        )
479          ENDDO          ENDDO
480         ENDDO         ENDDO
# Line 469  C-    No need for maskC(k-1) with recip_ Line 516  C-    No need for maskC(k-1) with recip_
516       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
517       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
518       &        *recip_drC(k)*maskC(i,j,k,bi,bj)       &        *recip_drC(k)*maskC(i,j,k,bi,bj)
519    #ifdef ALLOW_GGL90_IDEMIX
520         &        *recip_hFacI(i,j,k,bi,bj)
521    #endif
522          ENDDO          ENDDO
523         ENDDO         ENDDO
524        ENDDO        ENDDO
# Line 489  C-    No need for maskC(k) with recip_hF Line 539  C-    No need for maskC(k) with recip_hF
539       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
540       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
541       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
542    #ifdef ALLOW_GGL90_IDEMIX
543         &        *recip_hFacI(i,j,k,bi,bj)
544    #endif
545          ENDDO          ENDDO
546         ENDDO         ENDDO
547        ENDDO        ENDDO
548    
549          IF (.NOT.GGL90_dirichlet) THEN
550    C      Neumann bottom boundary condition for TKE: no flux from bottom
551           DO j=jMin,jMax
552            DO i=iMin,iMax
553             kBottom   = MAX(kLowC(i,j,bi,bj),1)
554             c3d(i,j,kBottom) = 0. _d 0
555            ENDDO
556           ENDDO
557          ENDIF
558    
559  C--   Center diagonal  C--   Center diagonal
560        DO k=1,Nr        DO k=1,Nr
561         km1 = MAX(k-1,1)         km1 = MAX(k-1,1)
# Line 523  C     Dirichlet surface boundary conditi Line 587  C     Dirichlet surface boundary conditi
587          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
588       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
589          a3d(i,j,kp1) = 0. _d 0          a3d(i,j,kp1) = 0. _d 0
 C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  
         kBottom   = MAX(kLowC(i,j,bi,bj),1)  
         GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)  
      &                              - GGL90TKEbottom*c3d(i,j,kBottom)  
         c3d(i,j,kBottom) = 0. _d 0  
590         ENDDO         ENDDO
591        ENDDO        ENDDO
592    
593          IF (GGL90_dirichlet) THEN
594    C      Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
595           DO j=jMin,jMax
596            DO i=iMin,iMax
597             kBottom   = MAX(kLowC(i,j,bi,bj),1)
598             GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
599         &                              - GGL90TKEbottom*c3d(i,j,kBottom)
600             c3d(i,j,kBottom) = 0. _d 0
601            ENDDO
602           ENDDO
603          ENDIF
604    
605  C     solve tri-diagonal system  C     solve tri-diagonal system
606        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
607       I                        a3d, b3d, c3d,       I                        a3d, b3d, c3d,
# Line 667  C     =============================== Line 738  C     ===============================
738       &                          0,Nr, 2, bi, bj, myThid )       &                          0,Nr, 2, bi, bj, myThid )
739           CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',           CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
740       &                          0,Nr, 2, bi, bj, myThid )       &                          0,Nr, 2, bi, bj, myThid )
741    
742          kp1 = MIN(Nr,kSurf+1)
743          DO j=jMin,jMax
744           DO i=iMin,iMax
745    c        diagnose surface flux of TKE  
746             surf_flx_tke(i,j) =(GGL90TKE(i,j,kSurf,bi,bj)-
747         &                       GGL90TKE(i,j,kp1,bi,bj))
748         &        *recip_drF(kSurf)*recip_hFacC(i,j,kSurf,bi,bj)
749         &        *KappaE(i,j,kp1)
750    
751           ENDDO
752          ENDDO
753          CALL DIAGNOSTICS_FILL(surf_flx_tke,'GGL90flx',
754         &                      0,1,1,bi,bj,myThid)
755    
756          k=kSurf
757          DO j=jMin,jMax
758           DO i=iMin,iMax
759    c        diagnose work done by the wind
760             surf_flx_tke(i,j) =
761         &     .5 _d 0*( surfaceForcingU(i,  j,bi,bj)*uvel(i  ,j,k,bi,bj)
762         &              +surfaceForcingU(i+1,j,bi,bj)*uvel(i+1,j,k,bi,bj))
763         &  +  .5 _d 0*( surfaceForcingV(i,j,  bi,bj)*vvel(i,j  ,k,bi,bj)
764         &              +surfaceForcingV(i,j+1,bi,bj)*vvel(i,j+1,k,bi,bj))
765           ENDDO
766          ENDDO
767          CALL DIAGNOSTICS_FILL(surf_flx_tke,'GGL90tau',
768         &                      0,1,1,bi,bj,myThid)
769    
770    
771        ENDIF        ENDIF
772  #endif  #endif /* ALLOW_DIAGNOSTICS */
773    
774  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
775    
Legend:
Removed from v.1.26  
changed lines
  Added in v.1.28
  ViewVC Help
Powered by ViewVC 1.1.22