/[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.9 by dfer, Mon Sep 29 13:47:23 2008 UTC revision 1.10 by dfer, Tue Oct 7 19:35:10 2008 UTC
# Line 65  C     RiNumber         - local Richardso Line 65  C     RiNumber         - local Richardso
65  C     KappaM           - (local) viscosity parameter (eq.10)  C     KappaM           - (local) viscosity parameter (eq.10)
66  C     KappaH           - (local) diffusivity parameter for temperature (eq.11)  C     KappaH           - (local) diffusivity parameter for temperature (eq.11)
67  C     KappaE           - (local) diffusivity parameter for TKE (eq.15)  C     KappaE           - (local) diffusivity parameter for TKE (eq.15)
 C     buoyFreq         - buoyancy freqency  
68  C     TKEdissipation   - dissipation of TKE  C     TKEdissipation   - dissipation of TKE
69  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse
70  C         rMixingLength- inverse of mixing length  C         rMixingLength- inverse of mixing length
# Line 138  C     Initialize local fields Line 137  C     Initialize local fields
137          rhoK    (I,J) = 0. _d 0          rhoK    (I,J) = 0. _d 0
138          rhoKm1  (I,J) = 0. _d 0          rhoKm1  (I,J) = 0. _d 0
139          totalDepth(I,J)   = 0. _d 0          totalDepth(I,J)   = 0. _d 0
140          IF ( recip_Rcol(I,J,bi,bj) .NE. 0. )          IF ( recip_Rcol(I,J,bi,bj) .NE. 0. _d 0 )
141       &       totalDepth(I,J) = 1./recip_Rcol(I,J,bi,bj)       &       totalDepth(I,J) = 1./recip_Rcol(I,J,bi,bj)
142         ENDDO         ENDDO
143        ENDDO        ENDDO
# Line 167  C Line 166  C
166           Nsquare = - gravity*recip_rhoConst*recip_drC(K)           Nsquare = - gravity*recip_rhoConst*recip_drC(K)
167       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)
168  C     vertical shear term (dU/dz)^2+(dV/dz)^2  C     vertical shear term (dU/dz)^2+(dV/dz)^2
169           tempu= .5*(  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)
170       &             - (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)) )
171       &        *recip_drC(K)       &        *recip_drC(K)
172           tempv= .5*(  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)
173       &             - (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)) )
174       &        *recip_drC(K)       &        *recip_drC(K)
175           verticalShear = tempU*tempU + tempV*tempV           verticalShear = tempU*tempU + tempV*tempV
176           RiNumber   = MAX(Nsquare,0. _d 0)/(verticalShear+GGL90eps)           RiNumber   = MAX(Nsquare,0. _d 0)/(verticalShear+GGL90eps)
177  C     compute Prandtl number (always greater than 0)  C     compute Prandtl number (always greater than 0)
178           prTemp = 1. _d 0           prTemp = 1. _d 0
179           IF ( RiNumber .GE. 0.2 ) prTemp = 5.0 * RiNumber           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
180           TKEPrandtlNumber(I,J,K) = MIN(10.0 _d 0,prTemp)           TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)
181  C     mixing length  C     mixing length
182           GGL90mixingLength(I,J,K) = SQRTTWO *           GGL90mixingLength(I,J,K) = SQRTTWO *
183       &        SQRTTKE/SQRT( MAX(Nsquare,GGL90eps) )       &        SQRTTKE/SQRT( MAX(Nsquare,GGL90eps) )
# Line 191  C     impose minimum mixing length (to a Line 190  C     impose minimum mixing length (to a
190           rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)           rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)
191  C     viscosity of last timestep  C     viscosity of last timestep
192           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE
193             KappaH = KappaM/TKEPrandtlNumber(I,J,K)
194    
195    C     Set a minium (= background) and maximum value
196             KappaM = MAX(KappaM,viscAr)
197             KappaH = MAX(KappaH,diffKrNrT(k))
198             KappaM = MIN(KappaM,GGL90viscMax)
199             KappaH = MIN(KappaH,GGL90diffMax)
200    
201    C     Mask land points and save
202           KappaE(I,J,K) = KappaM*GGL90alpha           KappaE(I,J,K) = KappaM*GGL90alpha
203             GGL90viscAr(I,J,K,bi,bj) = KappaM * maskC(I,J,K,bi,bj)
204             GGL90diffKr(I,J,K,bi,bj) = KappaH * maskC(I,J,K,bi,bj)
205    
206  C     dissipation term  C     dissipation term
207           TKEdissipation = ab05*GGL90ceps           TKEdissipation = ab05*GGL90ceps
208       &        *SQRTTKE*rMixingLength(I,J,K)       &        *SQRTTKE*rMixingLength(I,J,K)
# Line 200  C     sum up contributions to form the r Line 211  C     sum up contributions to form the r
211           gTKE(I,J,K) = GGL90TKE(I,J,K,bi,bj)           gTKE(I,J,K) = GGL90TKE(I,J,K,bi,bj)
212       &        + deltaTggl90*(       &        + deltaTggl90*(
213       &        + KappaM*verticalShear       &        + KappaM*verticalShear
214       &        - KappaM*Nsquare/TKEPrandtlNumber(I,J,K)       &        - KappaH*Nsquare
215    c    &        - KappaM*Nsquare/TKEPrandtlNumber(I,J,K)
216       &        - TKEdissipation       &        - TKEdissipation
217       &        )       &        )
218          ENDDO          ENDDO
# Line 223  C     common factors Line 235  C     common factors
235  C     Compute diffusive fluxes  C     Compute diffusive fluxes
236  C     ... across x-faces  C     ... across x-faces
237          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
238           dfx(1-Olx,j)=0.           dfx(1-Olx,j)=0. _d 0
239           DO i=1-Olx+1,sNx+Olx           DO i=1-Olx+1,sNx+Olx
240            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
241       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
# Line 233  C     ... across x-faces Line 245  C     ... across x-faces
245          ENDDO          ENDDO
246  C     ... across y-faces  C     ... across y-faces
247          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
248           dfy(i,1-Oly)=0.           dfy(i,1-Oly)=0. _d 0
249          ENDDO          ENDDO
250          DO j=1-Oly+1,sNy+Oly          DO j=1-Oly+1,sNy+Oly
251           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
# Line 278  C--   Lower diagonal Line 290  C--   Lower diagonal
290          DO i=iMin,iMax          DO i=iMin,iMax
291           a(i,j,k) = -deltaTggl90           a(i,j,k) = -deltaTggl90
292       &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)       &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)
293       &        *.5*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
294       &        *recip_drC(k)       &        *recip_drC(k)
295            IF (recip_hFacI(i,j,km1,bi,bj).EQ.0.) a(i,j,k)=0.            IF (recip_hFacI(i,j,km1,bi,bj).EQ.0. _d 0) a(i,j,k)=0. _d 0
296          ENDDO          ENDDO
297         ENDDO         ENDDO
298        ENDDO        ENDDO
# Line 298  CML      DO k=1,Nr-1 Line 310  CML      DO k=1,Nr-1
310          DO i=iMin,iMax          DO i=iMin,iMax
311            c(i,j,k) = -deltaTggl90            c(i,j,k) = -deltaTggl90
312       &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)       &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)
313       &               *.5*(KappaE(i,j,k)+KappaE(i,j,kp1))       &               *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
314       &        *recip_drC(k)       &        *recip_drC(k)
315            IF (recip_hFacI(i,j,kp1,bi,bj).EQ.0.) c(i,j,k)=0.            IF (recip_hFacI(i,j,kp1,bi,bj).EQ.0. _d 0) c(i,j,k)=0. _d 0
316          ENDDO          ENDDO
317         ENDDO         ENDDO
318        ENDDO        ENDDO
# Line 323  C Line 335  C
335         DO I=iMin,iMax         DO I=iMin,iMax
336  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
337          uStarSquare = SQRT(          uStarSquare = SQRT(
338       &         ( .5*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(I,  J,  bi,bj)
339       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2
340       &       + ( .5*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(I,  J,  bi,bj)
341       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2
342       &                     )       &                     )
343  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
# Line 361  C     end of time step Line 373  C     end of time step
373  C     ===============================  C     ===============================
374  C     compute viscosity coefficients  C     compute viscosity coefficients
375  C  C
376        DO K=2,Nr  c      DO K=2,Nr
377         DO J=jMin,jMax  c       DO J=jMin,jMax
378          DO I=iMin,iMax  c        DO I=iMin,iMax
379  C     Eq. (11), (18) and (21)  C     Eq. (11), (18) and (21)
380           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*  c         KappaM = GGL90ck*GGL90mixingLength(I,J,K)*
381       &                  SQRT( GGL90TKE(I,J,K,bi,bj) )  c     &                  SQRT( GGL90TKE(I,J,K,bi,bj) )
382           KappaH = KappaM/TKEPrandtlNumber(I,J,K)  c         KappaH = KappaM/TKEPrandtlNumber(I,J,K)
383  C     Set a minium (= background) value  C     Set a minium (= background) value
384           KappaM = MAX(KappaM,viscAr)  c         KappaM = MAX(KappaM,viscAr)
385           KappaH = MAX(KappaH,diffKrNrT(k))  c         KappaH = MAX(KappaH,diffKrNrT(k))
386  C     Set a maximum and mask land point  C     Set a maximum and mask land point
387           GGL90viscAr(I,J,K,bi,bj) = MIN(KappaM,GGL90viscMax)  c        GGL90viscAr(I,J,K,bi,bj) = MIN(KappaM,GGL90viscMax)
388       &        * maskC(I,J,K,bi,bj)  c    &        * maskC(I,J,K,bi,bj)
389           GGL90diffKr(I,J,K,bi,bj) = MIN(KappaH,GGL90diffMax)  c        GGL90diffKr(I,J,K,bi,bj) = MIN(KappaH,GGL90diffMax)
390       &        * maskC(I,J,K,bi,bj)  c    &        * maskC(I,J,K,bi,bj)
391          ENDDO  c        ENDDO
392         ENDDO  c       ENDDO
393  C     end third k-loop  C     end third k-loop
394        ENDDO  c      ENDDO
395    
396    #ifdef ALLOW_DIAGNOSTICS
397          IF ( useDiagnostics ) THEN
398             CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
399         &                          0,Nr, 1, bi, bj, myThid )
400             CALL DIAGNOSTICS_FILL( GGL90viscAr,'GGL90Ar ',
401         &                          0,Nr, 1, bi, bj, myThid )
402             CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
403         &                          0,Nr, 1, bi, bj, myThid )
404             CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
405         &                          0,Nr, 2, bi, bj, myThid )
406             CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
407         &                          0,Nr, 2, bi, bj, myThid )
408          ENDIF
409    #endif
410    
411  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
412    
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22