/[MITgcm]/MITgcm/model/src/calc_gw.F
ViewVC logotype

Diff of /MITgcm/model/src/calc_gw.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.39 by jmc, Tue Apr 22 22:20:31 2008 UTC revision 1.40 by jmc, Fri Feb 27 01:47:41 2009 UTC
# Line 59  C     iMin,iMax Line 59  C     iMin,iMax
59  C     jMin,jMax  C     jMin,jMax
60  C     xA            :: W-Cell face area normal to X  C     xA            :: W-Cell face area normal to X
61  C     yA            :: W-Cell face area normal to Y  C     yA            :: W-Cell face area normal to Y
 C     rMinW,rMaxW   :: column boundaries (r-units) at Western  Edge  
 C     rMinS,rMaxS   :: column boundaries (r-units) at Southern Edge  
62  C     rThickC_W     :: thickness (in r-units) of W-Cell at Western Edge  C     rThickC_W     :: thickness (in r-units) of W-Cell at Western Edge
63  C     rThickC_S     :: thickness (in r-units) of W-Cell at Southern Edge  C     rThickC_S     :: thickness (in r-units) of W-Cell at Southern Edge
64  C     rThickC_C     :: thickness (in r-units) of W-Cell (centered on W pt)  C     rThickC_C     :: thickness (in r-units) of W-Cell (centered on W pt)
# Line 71  C     flxAdvUp      :: vertical mom. adv Line 69  C     flxAdvUp      :: vertical mom. adv
69  C     flxDisUp      :: vertical mom. dissipation flux, vertical direction (@ level k-1)  C     flxDisUp      :: vertical mom. dissipation flux, vertical direction (@ level k-1)
70  C     flx_Dn        :: vertical momentum flux, vertical direction (@ level k)  C     flx_Dn        :: vertical momentum flux, vertical direction (@ level k)
71  C     gwDiss        :: vertical momentum dissipation tendency  C     gwDiss        :: vertical momentum dissipation tendency
72    C     gwAdd         :: other tendencies (Coriolis, Metric-terms)
73    C     del2w         :: laplacian of wVel
74    C     wFld          :: local copy of wVel
75  C     i,j,k         :: Loop counters  C     i,j,k         :: Loop counters
76        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
77        _RS    xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS    xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78        _RS    yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS    yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RS    rMinW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS    rMaxW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS    rMinS (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS    rMaxS (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
79        _RL    rThickC_W    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    rThickC_W    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80        _RL    rThickC_S    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    rThickC_S    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81        _RL    rThickC_C    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    rThickC_C    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 91  C     i,j,k         :: Loop counters Line 88  C     i,j,k         :: Loop counters
88        _RL    gwDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    gwDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89        _RL    gwAdd (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    gwAdd (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90        _RL    del2w (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    del2w (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91          _RL    wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92        INTEGER i,j,k, kp1        INTEGER i,j,k, kp1
93        _RL  wOverride        _RL  wOverride
94        _RL  tmp_WbarZ        _RL  tmp_WbarZ
# Line 136  C-    Initialise gwDiss to zero Line 134  C-    Initialise gwDiss to zero
134           gwDiss(i,j) = 0.           gwDiss(i,j) = 0.
135         ENDDO         ENDDO
136        ENDDO        ENDDO
137          IF (momViscosity) THEN
138    C-    Initialize del2w to zero:
139            DO j=1-Oly,sNy+Oly
140             DO i=1-Olx,sNx+Olx
141               del2w(i,j) = 0. _d 0
142             ENDDO
143            ENDDO
144          ENDIF
145    
146  C--   Boundaries condition at top  C--   Boundaries condition at top
147        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
# Line 144  C--   Boundaries condition at top Line 150  C--   Boundaries condition at top
150           flxDisUp(i,j) = 0.           flxDisUp(i,j) = 0.
151         ENDDO         ENDDO
152        ENDDO        ENDDO
 C--   column boundaries :  
       IF (momViscosity) THEN  
         DO j=1-Oly,sNy+Oly  
          DO i=1-Olx+1,sNx+Olx  
            rMaxW(i,j) = MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )  
            rMinW(i,j) = MAX(   R_low(i-1,j,bi,bj),   R_low(i,j,bi,bj) )  
          ENDDO  
         ENDDO  
         DO j=1-Oly+1,sNy+Oly  
          DO i=1-Olx,sNx+Olx  
            rMaxS(i,j) = MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )  
            rMinS(i,j) = MAX(   R_low(i,j-1,bi,bj),   R_low(i,j,bi,bj) )  
          ENDDO  
         ENDDO  
       ENDIF  
153    
154  C---  Sweep down column  C---  Sweep down column
155        DO k=2,Nr        DO k=2,Nr
# Line 196  C-    valid in z & p coord.; also accura Line 187  C-    valid in z & p coord.; also accura
187           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
188            DO i=1-Olx+1,sNx+Olx            DO i=1-Olx+1,sNx+Olx
189             rThickC_W(i,j) = MAX( zeroRS,             rThickC_W(i,j) = MAX( zeroRS,
190       &                           MIN( rMaxW(i,j), rC(k-1) )       &                           MIN( rSurfW(i,j,bi,bj), rC(k-1) )
191       &                          -MAX( rMinW(i,j), rC(k)   )       &                          -MAX(  rLowW(i,j,bi,bj), rC(k)   )
192       &                         )       &                         )
193  C     W-Cell Western face area:  C     W-Cell Western face area:
194             xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)             xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
# Line 207  c    &                              *dee Line 198  c    &                              *dee
198           DO j=1-Oly+1,sNy+Oly           DO j=1-Oly+1,sNy+Oly
199            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
200             rThickC_S(i,j) = MAX( zeroRS,             rThickC_S(i,j) = MAX( zeroRS,
201       &                           MIN( rMaxS(i,j), rC(k-1) )       &                           MIN( rSurfS(i,j,bi,bj), rC(k-1) )
202       &                          -MAX( rMinS(i,j), rC(k)   )       &                          -MAX(  rLowS(i,j,bi,bj), rC(k)   )
203       &                         )       &                         )
204  C     W-Cell Southern face area:  C     W-Cell Southern face area:
205             yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)             yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)
# Line 261  c          ENDIF Line 252  c          ENDIF
252    
253  C--   horizontal bi-harmonic dissipation  C--   horizontal bi-harmonic dissipation
254          IF (momViscosity .AND. viscA4W.NE.0. ) THEN          IF (momViscosity .AND. viscA4W.NE.0. ) THEN
255    
256    C-    local copy of wVel:
257              DO j=1-Oly,sNy+Oly
258               DO i=1-Olx,sNx+Olx
259                 wFld(i,j) = wVel(i,j,k,bi,bj)
260               ENDDO
261              ENDDO
262  C-    calculate the horizontal Laplacian of vertical flow  C-    calculate the horizontal Laplacian of vertical flow
263  C     Zonal flux d/dx W  C     Zonal flux d/dx W
264              IF ( useCubedSphereExchange ) THEN
265    C     to compute d/dx(W), fill corners with appropriate values:
266                CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
267         &                                 wFld, bi,bj, myThid )
268              ENDIF
269            DO j=1-Oly,sNy+Oly            DO j=1-Oly,sNy+Oly
270             flx_EW(1-Olx,j)=0.             flx_EW(1-Olx,j)=0.
271             DO i=1-Olx+1,sNx+Olx             DO i=1-Olx+1,sNx+Olx
272              flx_EW(i,j) =              flx_EW(i,j) =
273       &              (wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))       &               ( wFld(i,j) - wFld(i-1,j) )
274       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
275  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
276       &              *sqCosFacU(j,bi,bj)       &              *sqCosFacU(j,bi,bj)
277  #endif  #endif
278             ENDDO             ENDDO
279            ENDDO            ENDDO
280    
281  C     Meridional flux d/dy W  C     Meridional flux d/dy W
282              IF ( useCubedSphereExchange ) THEN
283    C     to compute d/dy(W), fill corners with appropriate values:
284                CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
285         &                                 wFld, bi,bj, myThid )
286              ENDIF
287            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
288             flx_NS(i,1-Oly)=0.             flx_NS(i,1-Oly)=0.
289            ENDDO            ENDDO
290            DO j=1-Oly+1,sNy+Oly            DO j=1-Oly+1,sNy+Oly
291             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx
292              flx_NS(i,j) =              flx_NS(i,j) =
293       &              (wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))       &               ( wFld(i,j) - wFld(i,j-1) )
294       &              *_recip_dyC(i,j,bi,bj)*yA(i,j)       &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
295  #ifdef ISOTROPIC_COS_SCALING  #ifdef ISOTROPIC_COS_SCALING
296  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
# Line 292  C     Meridional flux d/dy W Line 301  C     Meridional flux d/dy W
301            ENDDO            ENDDO
302    
303  C     del^2 W  C     del^2 W
304  C     Difference of zonal fluxes ...  C     Divergence of horizontal fluxes
           DO j=1-Oly,sNy+Oly  
            DO i=1-Olx,sNx+Olx-1  
             del2w(i,j)=flx_EW(i+1,j)-flx_EW(i,j)  
            ENDDO  
            del2w(sNx+Olx,j)=0.  
           ENDDO  
   
 C     ... add difference of meridional fluxes and divide by volume  
305            DO j=1-Oly,sNy+Oly-1            DO j=1-Oly,sNy+Oly-1
306             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx-1
307              del2w(i,j) = ( del2w(i,j)              del2w(i,j) = ( ( flx_EW(i+1,j)-flx_EW(i,j) )
308       &                    +(flx_NS(i,j+1)-flx_NS(i,j))       &                    +( flx_NS(i,j+1)-flx_NS(i,j) )
309       &                   )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)       &                   )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
310       &                    *recip_deepFac2F(k)       &                    *recip_deepFac2F(k)
311             ENDDO             ENDDO
312            ENDDO            ENDDO
313  C-- No-slip BCs impose a drag at walls...  C end if biharmonic viscosity
 CML ************************************************************  
 CML   No-slip Boundary conditions for bi-harmonic dissipation  
 CML   need to be implemented here!  
 CML ************************************************************  
         ELSEIF (momViscosity) THEN  
 C-    Initialize del2w to zero:  
           DO j=1-Oly,sNy+Oly  
            DO i=1-Olx,sNx+Olx  
             del2w(i,j) = 0. _d 0  
            ENDDO  
           ENDDO  
314          ENDIF          ENDIF
315    
316          IF (momViscosity) THEN          IF (momViscosity) THEN
# Line 504  C--          prepare for next level (k+1 Line 494  C--          prepare for next level (k+1
494  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
495    
496  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
497        IF ( diagDiss )  CALL DIAGNOSTICS_FILL( gwDiss, 'Wm_Diss ',          IF ( diagDiss )  THEN
498       &                                        k, 1, 2, bi,bj, myThid )            CALL DIAGNOSTICS_FILL( gwDiss, 'Wm_Diss ',
499        IF ( diagAdvec ) CALL DIAGNOSTICS_FILL( gW, 'Wm_Advec',       &                           k, 1, 2, bi,bj, myThid )
500       &                                        k,Nr, 1, bi,bj, myThid )  C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
501    C        does it only if k=1 (never the case here)
502              IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Diss ',bi,bj,myThid)
503            ENDIF
504            IF ( diagAdvec ) THEN
505              CALL DIAGNOSTICS_FILL( gW, 'Wm_Advec',
506         &                           k,Nr, 1, bi,bj, myThid )
507              IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Advec',bi,bj,myThid)
508            ENDIF
509  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
510    
511  C--   Dissipation term inside the Adams-Bashforth:  C--   Dissipation term inside the Adams-Bashforth:
Legend:
Removed from v.1.39  
changed lines
  Added in v.1.40
  ViewVC Help
Powered by ViewVC 1.1.22