/[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.28 by mlosch, Sat Feb 21 01:46:39 2015 UTC revision 1.29 by jmc, Sat Feb 21 17:13:20 2015 UTC
# Line 440  C     viscosity and diffusivity Line 440  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)
441           GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))           GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))
442       &                            * maskC(i,j,k,bi,bj)       &                            * maskC(i,j,k,bi,bj)
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)  C     compute Prandtl number (always greater than 0)
448           RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)           RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
449  #ifdef ALLOW_GGL90_IDEMIX  #ifdef ALLOW_GGL90_IDEMIX
450  CML         IDEMIX_RiNumber = 1./GGL90eps  CML         IDEMIX_RiNumber = 1./GGL90eps
451           IDEMIX_RiNumber = MAX( KappaM*Nsquare(i,j,k), 0. _d 0)/           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)       &    (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)           prTemp         = MIN(5.*RiNumber, 6.6*IDEMIX_RiNumber)
454  #else  #else
455           prTemp = 1. _d 0                   prTemp = 1. _d 0
456           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
457  #endif /* ALLOW_GGL90_IDEMIX */  #endif /* ALLOW_GGL90_IDEMIX */
458           TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)           TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
459           TKEPrandtlNumber(i,j,k) = MAX( 1. _d 0,TKEPrandtlNumber(i,j,k))           TKEPrandtlNumber(i,j,k) = MAX( 1. _d 0,TKEPrandtlNumber(i,j,k))
460    
461  c        diffusivity  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 726  C     =============================== Line 726  C     ===============================
726    
727  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
728        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
729           CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',          CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
730       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
731           CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',          CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
732       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
733           CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',          CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
734       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
735           CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',          CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
736       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
737           CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',          CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
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)          kp1 = MIN(Nr,kSurf+1)
743        DO j=jMin,jMax          DO j=jMin,jMax
744         DO i=iMin,iMax           DO i=iMin,iMax
745  c        diagnose surface flux of TKE    C     diagnose surface flux of TKE
746           surf_flx_tke(i,j) =(GGL90TKE(i,j,kSurf,bi,bj)-            surf_flx_tke(i,j) =(GGL90TKE(i,j,kSurf,bi,bj)-
747       &                       GGL90TKE(i,j,kp1,bi,bj))       &                        GGL90TKE(i,j,kp1,bi,bj))
748       &        *recip_drF(kSurf)*recip_hFacC(i,j,kSurf,bi,bj)       &        *recip_drF(kSurf)*recip_hFacC(i,j,kSurf,bi,bj)
749       &        *KappaE(i,j,kp1)       &        *KappaE(i,j,kp1)
750             ENDDO
751            ENDDO
752            CALL DIAGNOSTICS_FILL( surf_flx_tke,'GGL90flx',
753         &                         0, 1, 2, bi, bj, myThid )
754    
755         ENDDO          k=kSurf
756        ENDDO          DO j=jMin,jMax
757        CALL DIAGNOSTICS_FILL(surf_flx_tke,'GGL90flx',           DO i=iMin,iMax
758       &                      0,1,1,bi,bj,myThid)  C     diagnose work done by the wind
759              surf_flx_tke(i,j) =
760        k=kSurf       &      halfRL*( surfaceForcingU(i,  j,bi,bj)*uVel(i  ,j,k,bi,bj)
761        DO j=jMin,jMax       &              +surfaceForcingU(i+1,j,bi,bj)*uVel(i+1,j,k,bi,bj))
762         DO i=iMin,iMax       &    + halfRL*( surfaceForcingV(i,j,  bi,bj)*vVel(i,j  ,k,bi,bj)
763  c        diagnose work done by the wind       &              +surfaceForcingV(i,j+1,bi,bj)*vVel(i,j+1,k,bi,bj))
764           surf_flx_tke(i,j) =           ENDDO
765       &     .5 _d 0*( surfaceForcingU(i,  j,bi,bj)*uvel(i  ,j,k,bi,bj)          ENDDO
766       &              +surfaceForcingU(i+1,j,bi,bj)*uvel(i+1,j,k,bi,bj))          CALL DIAGNOSTICS_FILL( surf_flx_tke,'GGL90tau',
767       &  +  .5 _d 0*( surfaceForcingV(i,j,  bi,bj)*vvel(i,j  ,k,bi,bj)       &                         0, 1, 2, bi, bj, myThid )
      &              +surfaceForcingV(i,j+1,bi,bj)*vvel(i,j+1,k,bi,bj))  
        ENDDO  
       ENDDO  
       CALL DIAGNOSTICS_FILL(surf_flx_tke,'GGL90tau',  
      &                      0,1,1,bi,bj,myThid)  
   
768    
769        ENDIF        ENDIF
770  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
Legend:
Removed from v.1.28  
changed lines
  Added in v.1.29
  ViewVC Help
Powered by ViewVC 1.1.22