/[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.1 by mlosch, Thu Sep 16 11:27:18 2004 UTC revision 1.2 by mlosch, Mon Sep 27 08:02:04 2004 UTC
# Line 95  C     tri-diagonal matrix Line 95  C     tri-diagonal matrix
95        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
96        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
98    #ifdef ALLOW_GGL90_HORIZDIFF
99    C     xA, yA   - area of lateral faces
100    C     dfx, dfy - diffusive flux across lateral faces
101          _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102          _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103          _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104          _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105    #endif /* ALLOW_GGL90_HORIZDIFF */
106  CEOP  CEOP
107        iMin = 2-OLx        iMin = 2-OLx
108        iMax = sNx+OLx-1        iMax = sNx+OLx-1
# Line 192  C     sum up contributions to form the r Line 200  C     sum up contributions to form the r
200       &        )       &        )
201          ENDDO            ENDDO  
202         ENDDO         ENDDO
203        ENDDO            ENDDO
204  C      C     horizontal diffusion of TKE (requires an exchange in
205  C     Implicit time step to update TKE for k=1,Nr; TKE(Nr+1)=0 by default  C       do_fields_blocking_exchanges)
206  C      #ifdef ALLOW_GGL90_HORIZDIFF
207          IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
208           DO K=2,Nr
209    C     common factors
210            DO j=1-Oly,sNy+Oly
211             DO i=1-Olx,sNx+Olx
212              xA(i,j) = _dyG(i,j,bi,bj)
213         &         *drF(k)*_hFacW(i,j,k,bi,bj)
214              yA(i,j) = _dxG(i,j,bi,bj)
215         &         *drF(k)*_hFacS(i,j,k,bi,bj)
216             ENDDO
217            ENDDO  
218    C     Compute diffusive fluxes
219    C     ... across x-faces
220            DO j=1-Oly,sNy+Oly
221             dfx(1-Olx,j)=0.
222             DO i=1-Olx+1,sNx+Olx
223              dfx(i,j) = -GGL90diffTKEh*xA(i,j)
224         &      *_recip_dxC(i,j,bi,bj)
225         &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))
226         &      *CosFacU(j,bi,bj)
227             ENDDO
228            ENDDO
229    C     ... across y-faces
230            DO i=1-Olx,sNx+Olx
231             dfy(i,1-Oly)=0.
232            ENDDO
233            DO j=1-Oly+1,sNy+Oly
234             DO i=1-Olx,sNx+Olx
235              dfy(i,j) = -GGL90diffTKEh*yA(i,j)
236         &      *_recip_dyC(i,j,bi,bj)
237         &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))
238    #ifdef ISOTROPIC_COS_SCALING
239         &      *CosFacV(j,bi,bj)
240    #endif /* ISOTROPIC_COS_SCALING */
241             ENDDO
242            ENDDO  
243    C     Compute divergence of fluxes
244            DO j=1-Oly,sNy+Oly-1
245             DO i=1-Olx,sNx+Olx-1
246              gTKE(i,j,k)=gTKE(i,j,k)
247         &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)
248         &         *( (dfx(i+1,j)-dfx(i,j))
249         &           +(dfy(i,j+1)-dfy(i,j))
250         &           )
251             ENDDO  
252            ENDDO
253    C       end of k-loop
254           ENDDO
255    C     end if GGL90diffTKEh .eq. 0.
256          ENDIF
257    #endif /* ALLOW_GGL90_HORIZDIFF */
258    
259    C     ============================================
260    C     Implicit time step to update TKE for k=1,Nr;
261    C     TKE(Nr+1)=0 by default
262    C     ============================================
263  C     set up matrix  C     set up matrix
264  C--   Old aLower  C--   Lower diagonal
265        DO j=jMin,jMax        DO j=jMin,jMax
266         DO i=iMin,iMax         DO i=iMin,iMax
267           a(i,j,1) = 0. _d 0           a(i,j,1) = 0. _d 0
# Line 215  C--   Old aLower Line 279  C--   Old aLower
279          ENDDO          ENDDO
280         ENDDO         ENDDO
281        ENDDO        ENDDO
282    C--   Upper diagonal
 C--   Old aUpper  
283        DO j=jMin,jMax        DO j=jMin,jMax
284         DO i=iMin,iMax         DO i=iMin,iMax
285           c(i,j,1)  = 0. _d 0           c(i,j,1)  = 0. _d 0
# Line 224  C--   Old aUpper Line 287  C--   Old aUpper
287         ENDDO         ENDDO
288        ENDDO        ENDDO
289  CML      DO k=1,Nr-1  CML      DO k=1,Nr-1
290        DO k=2,Nr        DO k=2,Nr-1
291         kp1=min(Nr,k+1)         kp1=min(Nr,k+1)
292         DO j=jMin,jMax         DO j=jMin,jMax
293          DO i=iMin,iMax          DO i=iMin,iMax
# Line 232  CML      DO k=1,Nr-1 Line 295  CML      DO k=1,Nr-1
295       &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)       &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)
296       &               *.5*(KappaE(i,j,k)+KappaE(i,j,kp1))       &               *.5*(KappaE(i,j,k)+KappaE(i,j,kp1))
297       &        *recip_drC(k)       &        *recip_drC(k)
298  C          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.) c(i,j,k)=0.
299          ENDDO          ENDDO
300         ENDDO         ENDDO
301        ENDDO        ENDDO
302    C--   Center diagonal
 C--   Old aCenter  
303        DO k=1,Nr        DO k=1,Nr
304         DO j=jMin,jMax         DO j=jMin,jMax
305          DO i=iMin,iMax          DO i=iMin,iMax
# Line 286  C Line 348  C
348  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
349           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )
350       &        * maskC(I,J,K,bi,bj)       &        * maskC(I,J,K,bi,bj)
 C  
 C     end of time step  
 C  
351          ENDDO          ENDDO
352         ENDDO         ENDDO
353        ENDDO            ENDDO    
354  C  C
355    C     end of time step
356    C     ===============================
357  C     compute viscosity coefficients  C     compute viscosity coefficients
358  C      C    
359        DO K=2,Nr        DO K=2,Nr
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22