/[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.7 by mlosch, Tue Jun 6 22:15:19 2006 UTC revision 1.8 by jmc, Mon Aug 11 22:28:06 2008 UTC
# Line 60  C     K, Kp1, km1, kSurf, kBottom  - ver Line 60  C     K, Kp1, km1, kSurf, kBottom  - ver
60  C     ab15, ab05       - weights for implicit timestepping  C     ab15, ab05       - weights for implicit timestepping
61  C     uStarSquare      - square of friction velocity  C     uStarSquare      - square of friction velocity
62  C     verticalShear    - (squared) vertical shear of horizontal velocity  C     verticalShear    - (squared) vertical shear of horizontal velocity
63  C     Nsquare          - squared buoyancy freqency  C     Nsquare          - squared buoyancy freqency
64  C     RiNumber         - local Richardson number  C     RiNumber         - local Richardson number
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)
# Line 113  CEOP Line 113  CEOP
113    
114  C     set separate time step (should be deltaTtracer)  C     set separate time step (should be deltaTtracer)
115        deltaTggl90 = dTtracerLev(1)        deltaTggl90 = dTtracerLev(1)
116  C      C
117        kSurf = 1        kSurf = 1
118  C     implicit timestepping weights for dissipation  C     implicit timestepping weights for dissipation
119        ab15 =  1.5 _d 0        ab15 =  1.5 _d 0
# Line 131  C     Initialize local fields Line 131  C     Initialize local fields
131           GGL90mixingLength(I,J,K) = GGL90mixingLengthMin           GGL90mixingLength(I,J,K) = GGL90mixingLengthMin
132               rMixingLength(I,J,K) = 0. _d 0               rMixingLength(I,J,K) = 0. _d 0
133          ENDDO          ENDDO
134         ENDDO             ENDDO
135        ENDDO        ENDDO
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          rhoK    (I,J) = 0. _d 0          rhoK    (I,J) = 0. _d 0
139          rhoKm1  (I,J) = 0. _d 0          rhoKm1  (I,J) = 0. _d 0
140          totalDepth(I,J)   = 0. _d 0          totalDepth(I,J)   = 0. _d 0
141          IF ( recip_Rcol(I,J,bi,bj) .NE. 0. )          IF ( recip_Rcol(I,J,bi,bj) .NE. 0. )
142       &       totalDepth(I,J) = 1./recip_Rcol(I,J,bi,bj)       &       totalDepth(I,J) = 1./recip_Rcol(I,J,bi,bj)
143         ENDDO         ENDDO
144        ENDDO        ENDDO
# Line 147  C     start k-loop Line 147  C     start k-loop
147        DO K = 2, Nr        DO K = 2, Nr
148         Km1 = K-1         Km1 = K-1
149         Kp1 = MIN(Nr,K+1)         Kp1 = MIN(Nr,K+1)
150         CALL FIND_RHO(         CALL FIND_RHO_2D(
151       I      bi, bj, iMin, iMax, jMin, jMax, Km1, K,       I      iMin, iMax, jMin, jMax, K,
152       I      theta, salt,       I      theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),
153       O      rhoKm1,       O      rhoKm1,
154       I      myThid )       I      Km1, bi, bj, myThid )
155         CALL FIND_RHO(  
156       I      bi, bj, iMin, iMax, jMin, jMax, K, K,         CALL FIND_RHO_2D(
157       I      theta, salt,       I      iMin, iMax, jMin, jMax, K,
158         I      theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj),
159       O      rhoK,       O      rhoK,
160       I      myThid )       I      K, bi, bj, myThid )
161         DO J=jMin,jMax         DO J=jMin,jMax
162          DO I=iMin,iMax          DO I=iMin,iMax
163           SQRTTKE=SQRT( GGL90TKE(I,J,K,bi,bj) )           SQRTTKE=SQRT( GGL90TKE(I,J,K,bi,bj) )
# Line 179  C     compute Prandtl number (always gre Line 180  C     compute Prandtl number (always gre
180           IF ( RiNumber .GE. 0.2 ) prTemp = 5.0 * RiNumber           IF ( RiNumber .GE. 0.2 ) prTemp = 5.0 * RiNumber
181           TKEPrandtlNumber(I,J,K) = MIN(10.0 _d 0,prTemp)           TKEPrandtlNumber(I,J,K) = MIN(10.0 _d 0,prTemp)
182  C     mixing length  C     mixing length
183           GGL90mixingLength(I,J,K) =           GGL90mixingLength(I,J,K) =
184       &        SQRTTKE/SQRT( MAX(Nsquare,GGL90eps) )       &        SQRTTKE/SQRT( MAX(Nsquare,GGL90eps) )
185  C     impose upper bound for mixing length (total depth)  C     impose upper bound for mixing length (total depth)
186           GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),           GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
187       &        totalDepth(I,J))       &        totalDepth(I,J))
188  C     impose minimum mixing length (to avoid division by zero)  C     impose minimum mixing length (to avoid division by zero)
189           GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),           GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),
190       &        GGL90mixingLengthMin)       &        GGL90mixingLengthMin)
191           rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)           rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)
192  C     viscosity of last timestep  C     viscosity of last timestep
# Line 194  C     viscosity of last timestep Line 195  C     viscosity of last timestep
195  C     dissipation term  C     dissipation term
196           TKEdissipation = ab05*GGL90ceps           TKEdissipation = ab05*GGL90ceps
197       &        *SQRTTKE*rMixingLength(I,J,K)       &        *SQRTTKE*rMixingLength(I,J,K)
198       &        *GGL90TKE(I,J,K,bi,bj)             &        *GGL90TKE(I,J,K,bi,bj)
199  C     sum up contributions to form the right hand side  C     sum up contributions to form the right hand side
200           gTKE(I,J,K) = GGL90TKE(I,J,K,bi,bj)           gTKE(I,J,K) = GGL90TKE(I,J,K,bi,bj)
201       &        + deltaTggl90*(       &        + deltaTggl90*(
202       &        + KappaM*verticalShear       &        + KappaM*verticalShear
203       &        - KappaM*Nsquare/TKEPrandtlNumber(I,J,K)       &        - KappaM*Nsquare/TKEPrandtlNumber(I,J,K)
204       &        - TKEdissipation       &        - TKEdissipation
205       &        )       &        )
206          ENDDO            ENDDO
207         ENDDO         ENDDO
208        ENDDO        ENDDO
209  C     horizontal diffusion of TKE (requires an exchange in  C     horizontal diffusion of TKE (requires an exchange in
210  C       do_fields_blocking_exchanges)  C      do_fields_blocking_exchanges)
211  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
212        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
213         DO K=2,Nr         DO K=2,Nr
# Line 218  C     common factors Line 219  C     common factors
219            yA(i,j) = _dxG(i,j,bi,bj)            yA(i,j) = _dxG(i,j,bi,bj)
220       &         *drF(k)*_hFacS(i,j,k,bi,bj)       &         *drF(k)*_hFacS(i,j,k,bi,bj)
221           ENDDO           ENDDO
222          ENDDO            ENDDO
223  C     Compute diffusive fluxes  C     Compute diffusive fluxes
224  C     ... across x-faces  C     ... across x-faces
225          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
# Line 243  C     ... across y-faces Line 244  C     ... across y-faces
244       &      *CosFacV(j,bi,bj)       &      *CosFacV(j,bi,bj)
245  #endif /* ISOTROPIC_COS_SCALING */  #endif /* ISOTROPIC_COS_SCALING */
246           ENDDO           ENDDO
247          ENDDO            ENDDO
248  C     Compute divergence of fluxes  C     Compute divergence of fluxes
249          DO j=1-Oly,sNy+Oly-1          DO j=1-Oly,sNy+Oly-1
250           DO i=1-Olx,sNx+Olx-1           DO i=1-Olx,sNx+Olx-1
251            gTKE(i,j,k)=gTKE(i,j,k)            gTKE(i,j,k)=gTKE(i,j,k)
252       &   -_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)
253       &         *( (dfx(i+1,j)-dfx(i,j))       &         *( (dfx(i+1,j)-dfx(i,j))
254       &           +(dfy(i,j+1)-dfy(i,j))       &           +(dfy(i,j+1)-dfy(i,j))
255       &           )       &           )
256           ENDDO             ENDDO
257          ENDDO          ENDDO
258  C       end of k-loop  C       end of k-loop
259         ENDDO         ENDDO
260  C     end if GGL90diffTKEh .eq. 0.  C     end if GGL90diffTKEh .eq. 0.
261        ENDIF        ENDIF
# Line 268  C     set up matrix Line 269  C     set up matrix
269  C--   Lower diagonal  C--   Lower diagonal
270        DO j=jMin,jMax        DO j=jMin,jMax
271         DO i=iMin,iMax         DO i=iMin,iMax
272           a(i,j,1) = 0. _d 0           a(i,j,1) = 0. _d 0
273         ENDDO         ENDDO
274        ENDDO        ENDDO
275        DO k=2,Nr        DO k=2,Nr
# Line 317  C     end set up matrix Line 318  C     end set up matrix
318    
319  C  C
320  C     Apply boundary condition  C     Apply boundary condition
321  C      C
322        DO J=jMin,jMax        DO J=jMin,jMax
323         DO I=iMin,iMax         DO I=iMin,iMax
324  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
325          uStarSquare = SQRT(          uStarSquare = SQRT(
326       &         ( .5*( surfaceForcingU(I,  J,  bi,bj)       &         ( .5*( surfaceForcingU(I,  J,  bi,bj)
327       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2
328       &       + ( .5*( surfaceForcingV(I,  J,  bi,bj)       &       + ( .5*( surfaceForcingV(I,  J,  bi,bj)
329       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2
330       &                     )       &                     )
331  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
# Line 332  C     Dirichlet surface boundary conditi Line 333  C     Dirichlet surface boundary conditi
333       &                     *maskC(I,J,kSurf,bi,bj)       &                     *maskC(I,J,kSurf,bi,bj)
334  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
335          kBottom   = MIN(MAX(kLowC(I,J,bi,bj),1),Nr)          kBottom   = MIN(MAX(kLowC(I,J,bi,bj),1),Nr)
336          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)
337       &       - GGL90TKEbottom*c(I,J,kBottom)       &       - GGL90TKEbottom*c(I,J,kBottom)
338         ENDDO         ENDDO
339        ENDDO            ENDDO
340  C  C
341  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)
342  C  C
# Line 345  C Line 346  C
346       I     myThid )       I     myThid )
347  C  C
348  C     now update TKE  C     now update TKE
349  C      C
350        DO K=1,Nr        DO K=1,Nr
351         DO J=jMin,jMax         DO J=jMin,jMax
352          DO I=iMin,iMax          DO I=iMin,iMax
353  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
354           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )
355       &        * maskC(I,J,K,bi,bj)       &        * maskC(I,J,K,bi,bj)
356          ENDDO          ENDDO
357         ENDDO         ENDDO
358        ENDDO            ENDDO
359  C  C
360  C     end of time step  C     end of time step
361  C     ===============================  C     ===============================
362  C     compute viscosity coefficients  C     compute viscosity coefficients
363  C      C
364        DO K=2,Nr        DO K=2,Nr
365         DO J=jMin,jMax         DO J=jMin,jMax
366          DO I=iMin,iMax          DO I=iMin,iMax
# Line 376  C     Set a maximum and mask land point Line 377  C     Set a maximum and mask land point
377           GGL90diffKr(I,J,K,bi,bj) = MIN(KappaH,GGL90diffMax)           GGL90diffKr(I,J,K,bi,bj) = MIN(KappaH,GGL90diffMax)
378       &        * maskC(I,J,K,bi,bj)       &        * maskC(I,J,K,bi,bj)
379          ENDDO          ENDDO
380         ENDDO         ENDDO
381  C     end third k-loop  C     end third k-loop
382        ENDDO            ENDDO
383    
384  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
385    
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22