/[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.27 by mlosch, Thu Feb 19 15:44:12 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          _RL     IDEMIX_RiNumber
89        _RL     TKEdissipation        _RL     TKEdissipation
90        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
91        _RL     MaxLength, tmpmlx, tmpVisc        _RL     MaxLength, tmpmlx, tmpVisc
# Line 96  c     _RL     SQRTTKE Line 96  c     _RL     SQRTTKE
96        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
99    #ifdef ALLOW_DIAGNOSTICS
100          _RL     surf_flx_tke     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101    #endif /* ALLOW_DIAGNOSTICS */
102  C-    tri-diagonal matrix  C-    tri-diagonal matrix
103        _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
104        _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 118  C     gTKE     :: right hand side of dif
118  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
119  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
120        _RL p4, p8, p16        _RL p4, p8, p16
121    CEOP
122        p4=0.25 _d 0        p4=0.25 _d 0
123        p8=0.125 _d 0        p8=0.125 _d 0
124        p16=0.0625 _d 0        p16=0.0625 _d 0
# Line 157  C     Initialize local fields Line 161  C     Initialize local fields
161         ENDDO         ENDDO
162        ENDDO        ENDDO
163    
164    #ifdef ALLOW_GGL90_IDEMIX
165          IF ( useIDEMIX) CALL GGL90_IDEMIX(
166         &   bi, bj, sigmaR, myTime, myIter, myThid )
167    #endif /* ALLOW_GGL90_IDEMIX */
168    
169  C     start k-loop  C     start k-loop
170        DO k = 2, Nr        DO k = 2, Nr
171  c      km1 = k-1  c      km1 = k-1
# Line 378  C     ... across y-faces Line 387  C     ... across y-faces
387  C     Compute divergence of fluxes  C     Compute divergence of fluxes
388          DO j=1-OLy,sNy+OLy-1          DO j=1-OLy,sNy+OLy-1
389           DO i=1-OLx,sNx+OLx-1           DO i=1-OLx,sNx+OLx-1
390            hFac = min(.5 _d 0,_hFacC(i,j,k-1,bi,bj) ) +  #ifdef ALLOW_GGL90_IDEMIX
391       &          min(.5 _d 0,_hFacC(i,j,k  ,bi,bj) )            gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)
392         &         *recip_hFacI(i,j,k,bi,bj)
393    #else
394              hFac = MIN(.5 _d 0,_hFacC(i,j,k-1,bi,bj) ) +
395         &           MIN(.5 _d 0,_hFacC(i,j,k  ,bi,bj) )
396            gTKE(i,j) = 0.0            gTKE(i,j) = 0.0
397            if ( hFac .ne. 0.0 )            IF ( hFac .ne. 0.0 )
398       &      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
399    #endif
400       &         *((dfx(i+1,j)-dfx(i,j))       &         *((dfx(i+1,j)-dfx(i,j))
401       &          +(dfy(i,j+1)-dfy(i,j)) )       &          +(dfy(i,j+1)-dfy(i,j)) )
402           ENDDO           ENDDO
# Line 397  C     vertical shear term (dU/dz)^2+(dV/ Line 411  C     vertical shear term (dU/dz)^2+(dV/
411           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)
412       &                 -( 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)) )
413       &        *recip_drC(k)       &        *recip_drC(k)
414    #ifdef ALLOW_GGL90_IDEMIX
415         &        *recip_hFacI(i,j,k,bi,bj)
416    #endif
417           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)
418       &                 -( 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)) )
419       &        *recip_drC(k)       &        *recip_drC(k)
420    #ifdef ALLOW_GGL90_IDEMIX
421         &        *recip_hFacI(i,j,k,bi,bj)
422    #endif
423           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  
424    
425  C     viscosity and diffusivity  C     viscosity and diffusivity
426           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 429  C     viscosity and diffusivity
429  c        note: storing GGL90visctmp like this, and using it later to compute  c        note: storing GGL90visctmp like this, and using it later to compute
430  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)
431           KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)           KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)
432    
433    C     compute Prandtl number (always greater than 0)
434             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
435    CML         IDEMIX_RiNumber = 1./GGL90eps
436    #ifdef ALLOW_GGL90_IDEMIX
437             IDEMIX_RiNumber = MAX( KappaM*Nsquare(i,j,k), 0. _d 0)/
438         &    (GGL90eps+IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2)
439             prTemp         = MIN(5.*RiNumber, 6.6*IDEMIX_RiNumber)
440    #else
441             IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
442    #endif /* ALLOW_GGL90_IDEMIX */
443             TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
444             TKEPrandtlNumber(i,j,k) = MAX( 1. _d 0,TKEPrandtlNumber(i,j,k))
445    
446    c        diffusivity
447           KappaH = KappaM/TKEPrandtlNumber(i,j,k)           KappaH = KappaM/TKEPrandtlNumber(i,j,k)
448           KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)           KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)
449    
# Line 428  C     partial update with sum of explici Line 457  C     partial update with sum of explici
457       &        + KappaM*verticalShear       &        + KappaM*verticalShear
458       &        - KappaH*Nsquare(i,j,k)       &        - KappaH*Nsquare(i,j,k)
459       &        - TKEdissipation       &        - TKEdissipation
460    #ifdef ALLOW_GGL90_IDEMIX
461         &        + IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2
462    #endif
463       &        )       &        )
464          ENDDO          ENDDO
465         ENDDO         ENDDO
# Line 469  C-    No need for maskC(k-1) with recip_ Line 501  C-    No need for maskC(k-1) with recip_
501       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
502       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
503       &        *recip_drC(k)*maskC(i,j,k,bi,bj)       &        *recip_drC(k)*maskC(i,j,k,bi,bj)
504    #ifdef ALLOW_GGL90_IDEMIX
505         &        *recip_hFacI(i,j,k,bi,bj)
506    #endif
507          ENDDO          ENDDO
508         ENDDO         ENDDO
509        ENDDO        ENDDO
# Line 489  C-    No need for maskC(k) with recip_hF Line 524  C-    No need for maskC(k) with recip_hF
524       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
525       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
526       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
527    #ifdef ALLOW_GGL90_IDEMIX
528         &        *recip_hFacI(i,j,k,bi,bj)
529    #endif
530          ENDDO          ENDDO
531         ENDDO         ENDDO
532        ENDDO        ENDDO
533    
534          IF (.NOT.GGL90_dirichlet) THEN
535    C      Neumann bottom boundary condition for TKE: no flux from bottom
536           DO j=jMin,jMax
537            DO i=iMin,iMax
538             kBottom   = MAX(kLowC(i,j,bi,bj),1)
539             c3d(i,j,kBottom) = 0. _d 0
540            ENDDO
541           ENDDO
542          ENDIF
543    
544  C--   Center diagonal  C--   Center diagonal
545        DO k=1,Nr        DO k=1,Nr
546         km1 = MAX(k-1,1)         km1 = MAX(k-1,1)
# Line 523  C     Dirichlet surface boundary conditi Line 572  C     Dirichlet surface boundary conditi
572          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
573       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
574          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  
575         ENDDO         ENDDO
576        ENDDO        ENDDO
577    
578          IF (GGL90_dirichlet) THEN
579    C      Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
580           DO j=jMin,jMax
581            DO i=iMin,iMax
582             kBottom   = MAX(kLowC(i,j,bi,bj),1)
583             GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
584         &                              - GGL90TKEbottom*c3d(i,j,kBottom)
585             c3d(i,j,kBottom) = 0. _d 0
586            ENDDO
587           ENDDO
588          ENDIF
589    
590  C     solve tri-diagonal system  C     solve tri-diagonal system
591        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
592       I                        a3d, b3d, c3d,       I                        a3d, b3d, c3d,
# Line 667  C     =============================== Line 723  C     ===============================
723       &                          0,Nr, 2, bi, bj, myThid )       &                          0,Nr, 2, bi, bj, myThid )
724           CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',           CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
725       &                          0,Nr, 2, bi, bj, myThid )       &                          0,Nr, 2, bi, bj, myThid )
726    
727          kp1 = MIN(Nr,kSurf+1)
728          DO j=jMin,jMax
729           DO i=iMin,iMax
730    c        diagnose surface flux of TKE  
731             surf_flx_tke(i,j) =(GGL90TKE(i,j,kSurf,bi,bj)-
732         &                       GGL90TKE(i,j,kp1,bi,bj))
733         &        *recip_drF(kSurf)*recip_hFacC(i,j,kSurf,bi,bj)
734         &        *KappaE(i,j,kp1)
735    
736           ENDDO
737          ENDDO
738          CALL DIAGNOSTICS_FILL(surf_flx_tke,'GGL90flx',
739         &                      0,1,1,bi,bj,myThid)
740    
741          k=kSurf
742          DO j=jMin,jMax
743           DO i=iMin,iMax
744    c        diagnose work done by the wind
745             surf_flx_tke(i,j) =
746         &     .5 _d 0*( surfaceForcingU(i,  j,bi,bj)*uvel(i  ,j,k,bi,bj)
747         &              +surfaceForcingU(i+1,j,bi,bj)*uvel(i+1,j,k,bi,bj))
748         &  +  .5 _d 0*( surfaceForcingV(i,j,  bi,bj)*vvel(i,j  ,k,bi,bj)
749         &              +surfaceForcingV(i,j+1,bi,bj)*vvel(i,j+1,k,bi,bj))
750           ENDDO
751          ENDDO
752          CALL DIAGNOSTICS_FILL(surf_flx_tke,'GGL90tau',
753         &                      0,1,1,bi,bj,myThid)
754    
755    
756        ENDIF        ENDIF
757  #endif  #endif /* ALLOW_DIAGNOSTICS */
758    
759  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
760    
Legend:
Removed from v.1.26  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.22