/[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.42 by jmc, Fri Dec 11 13:53:07 2009 UTC revision 1.43 by jmc, Sat Jan 23 00:04:03 2010 UTC
# Line 89  C     i,j,k         :: Loop counters Line 89  C     i,j,k         :: Loop counters
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)        _RL    wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92        INTEGER i,j,k, kp1        INTEGER i,j,k, km1, kp1
93        _RL  wOverride, surfFac        _RL  mskM1, mskP1
94        _RL  tmp_WbarZ        _RL  tmp_WbarZ
95        _RL  uTrans, vTrans, rTrans        _RL  uTrans, vTrans, rTrans
96        _RL  viscLoc        _RL  viscLoc
97        _RL  halfRL        _RL  halfRL
98        _RS  halfRS, zeroRS        _RS  halfRS, zeroRS
99        PARAMETER( halfRL = 0.5D0 )        PARAMETER( halfRL = 0.5 _d 0 )
100        PARAMETER( halfRS = 0.5 , zeroRS = 0. )        PARAMETER( halfRS = 0.5 , zeroRS = 0. )
101        PARAMETER( iMin = 1 , iMax = sNx )        PARAMETER( iMin = 1 , iMax = sNx )
102        PARAMETER( jMin = 1 , jMax = sNy )        PARAMETER( jMin = 1 , jMax = sNy )
# Line 107  CEOP Line 107  CEOP
107        EXTERNAL DIAGNOSTICS_IS_ON        EXTERNAL DIAGNOSTICS_IS_ON
108  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
109    
110  C--   Catch barotropic mode  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
       IF ( Nr .LT. 2 ) RETURN  
111    
112  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
113        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
# Line 144  C-    Initialize del2w to zero: Line 143  C-    Initialize del2w to zero:
143        ENDIF        ENDIF
144    
145  C--   Boundaries condition at top (vertical advection of vertical momentum):  C--   Boundaries condition at top (vertical advection of vertical momentum):
146  C     Setting surfFac=0 ignores surface wVel (as if wVel_k=1 = 0)        DO j=1-OLy,sNy+OLy
147        surfFac = 1. _d 0         DO i=1-OLx,sNx+OLx
148  c     surfFac = 0. _d 0           flxAdvUp(i,j) = 0.
149    c        flxDisUp(i,j) = 0.
150           ENDDO
151          ENDDO
152    
153    
154  C---  Sweep down column  C---  Sweep down column
155        DO k=2,Nr        DO k=1,Nr
156          kp1=k+1          km1 = MAX( k-1, 1 )
157          wOverRide=1.          kp1 = MIN( k+1,Nr )
158          IF (k.EQ.Nr) THEN          mskM1 = 1.
159            kp1=Nr          mskP1 = 1.
160            wOverRide=0.          IF ( k.EQ. 1 ) mskM1 = 0.
161          ENDIF          IF ( k.EQ.Nr ) mskP1 = 0.
162            IF ( k.GT.1 ) THEN
163  C--   Compute grid factor arround a W-point:  C--   Compute grid factor arround a W-point:
164  #ifdef CALC_GW_NEW_THICK  #ifdef CALC_GW_NEW_THICK
165          DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
166           DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
167             IF ( maskC(i,j,k-1,bi,bj).EQ.0. .OR.             IF ( maskC(i,j,k-1,bi,bj).EQ.0. .OR.
168       &          maskC(i,j, k ,bi,bj).EQ.0. ) THEN       &          maskC(i,j, k ,bi,bj).EQ.0. ) THEN
169               recip_rThickC(i,j) = 0.               recip_rThickC(i,j) = 0.
# Line 170  C-    valid in z & p coord.; also accura Line 174  C-    valid in z & p coord.; also accura
174       &         - MAX( R_low(i,j,bi,bj),  rC(k)   )       &         - MAX( R_low(i,j,bi,bj),  rC(k)   )
175       &        )       &        )
176             ENDIF             ENDIF
177              ENDDO
178           ENDDO           ENDDO
179          ENDDO           IF (momViscosity) THEN
         IF (momViscosity) THEN  
180           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
181            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
182             rThickC_C(i,j) = MAX( zeroRS,             rThickC_C(i,j) = MAX( zeroRS,
# Line 205  C deep-model: xA,yA is only used for vis Line 209  C deep-model: xA,yA is only used for vis
209  C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)  C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)
210            ENDDO            ENDDO
211           ENDDO           ENDDO
212          ENDIF           ENDIF
213  #else /* CALC_GW_NEW_THICK */  #else /* CALC_GW_NEW_THICK */
214          DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
215           DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
216  C-    note: assume fluid @ smaller k than bottom: does not work in p-coordinate !  C-    note: assume fluid @ smaller k than bottom: does not work in p-coordinate !
217             IF ( maskC(i,j,k,bi,bj).EQ.0. ) THEN             IF ( maskC(i,j,k,bi,bj).EQ.0. ) THEN
218               recip_rThickC(i,j) = 0.               recip_rThickC(i,j) = 0.
# Line 243  c    &                              *dee Line 247  c    &                              *dee
247  C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC.  C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC.
248  C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)  C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)
249  c          ENDIF  c          ENDIF
250              ENDDO
251           ENDDO           ENDDO
         ENDDO  
252  #endif /* CALC_GW_NEW_THICK */  #endif /* CALC_GW_NEW_THICK */
253            ELSEIF ( selectNHfreeSurf.GE.1 ) THEN
254             DO j=1-Oly,sNy+Oly
255              DO i=1-Olx,sNx+Olx
256               recip_rThickC(i,j) = recip_drC(k)
257    c          rThickC_C(i,j) = drC(k)
258    c          rThickC_W(i,j) = drC(k)
259    c          rThickC_S(i,j) = drC(k)
260    c          xA(i,j) = _dyG(i,j,bi,bj)*drC(k)
261    c          yA(i,j) = _dxG(i,j,bi,bj)*drC(k)
262              ENDDO
263             ENDDO
264            ENDIF
265    
266    C--   local copy of wVel:
267            DO j=1-Oly,sNy+Oly
268              DO i=1-Olx,sNx+Olx
269                wFld(i,j) = wVel(i,j,k,bi,bj)
270              ENDDO
271            ENDDO
272    
273  C--   horizontal bi-harmonic dissipation  C--   horizontal bi-harmonic dissipation
274          IF (momViscosity .AND. viscA4W.NE.0. ) THEN          IF ( momViscosity .AND. k.GT.1 .AND. viscA4W.NE.0. ) THEN
275    
 C-    local copy of wVel:  
           DO j=1-Oly,sNy+Oly  
            DO i=1-Olx,sNx+Olx  
              wFld(i,j) = wVel(i,j,k,bi,bj)  
            ENDDO  
           ENDDO  
276  C-    calculate the horizontal Laplacian of vertical flow  C-    calculate the horizontal Laplacian of vertical flow
277  C     Zonal flux d/dx W  C     Zonal flux d/dx W
278            IF ( useCubedSphereExchange ) THEN            IF ( useCubedSphereExchange ) THEN
# Line 310  C     Divergence of horizontal fluxes Line 327  C     Divergence of horizontal fluxes
327  C end if biharmonic viscosity  C end if biharmonic viscosity
328          ENDIF          ENDIF
329    
330          IF (momViscosity) THEN          IF ( momViscosity .AND. k.GT.1 ) THEN
331  C Viscous Flux on Western face  C Viscous Flux on Western face
332            DO j=jMin,jMax            DO j=jMin,jMax
333             DO i=iMin,iMax+1             DO i=iMin,iMax+1
# Line 361  C     Interpolate vert viscosity to cent Line 378  C     Interpolate vert viscosity to cent
378       &                  +KappaRV(i,j,kp1)+KappaRV(i,j+1,kp1)       &                  +KappaRV(i,j,kp1)+KappaRV(i,j+1,kp1)
379       &                 )*0.125 _d 0       &                 )*0.125 _d 0
380               flx_Dn(i,j) =               flx_Dn(i,j) =
381       &          - viscLoc*( wVel(i,j,kp1,bi,bj)*wOverRide       &          - viscLoc*( wVel(i,j,kp1,bi,bj)*mskP1
382       &                     -wVel(i,j, k ,bi,bj) )*rkSign       &                     -wVel(i,j, k ,bi,bj) )*rkSign
383       &                   *recip_drF(k)*rA(i,j,bi,bj)       &                   *recip_drF(k)*rA(i,j,bi,bj)
384       &                   *deepFac2C(k)*rhoFacC(k)       &                   *deepFac2C(k)*rhoFacC(k)
# Line 405  C--        prepare for next level (k+1) Line 422  C--        prepare for next level (k+1)
422            ENDDO            ENDDO
423          ENDIF          ENDIF
424    
425          IF ( momViscosity .AND. no_slip_sides ) THEN          IF ( momViscosity .AND. k.GT.1 .AND. no_slip_sides ) THEN
426  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
427            CALL MOM_W_SIDEDRAG(            CALL MOM_W_SIDEDRAG(
428       I               bi,bj,k,       I               bi,bj,k,
# Line 424  C-     No-slip BCs impose a drag at wall Line 441  C-     No-slip BCs impose a drag at wall
441  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
442    
443          IF ( momAdvection ) THEN          IF ( momAdvection ) THEN
444    
445             IF ( k.GT.1 .OR. selectNHfreeSurf.GE.1 ) THEN
446  C Advective Flux on Western face  C Advective Flux on Western face
447            DO j=jMin,jMax            DO j=jMin,jMax
448             DO i=iMin,iMax+1             DO i=iMin,iMax+1
449  C     transport through Western face area:  C     transport through Western face area:
450               uTrans = (               uTrans = (
451       &          drF(k-1)*_hFacW(i,j,k-1,bi,bj)*uVel(i,j,k-1,bi,bj)       &          drF(km1)*_hFacW(i,j,km1,bi,bj)*uVel(i,j,km1,bi,bj)
452       &                  *rhoFacC(k-1)       &                  *rhoFacC(km1)*mskM1
453       &        + drF( k )*_hFacW(i,j, k ,bi,bj)*uVel(i,j, k ,bi,bj)       &        + drF( k )*_hFacW(i,j, k ,bi,bj)*uVel(i,j, k ,bi,bj)
454       &                  *rhoFacC(k)       &                  *rhoFacC(k)
455       &                )*halfRL*_dyG(i,j,bi,bj)*deepFacF(k)       &                )*halfRL*_dyG(i,j,bi,bj)*deepFacF(k)
456               flx_EW(i,j)=               flx_EW(i,j) = uTrans*(wFld(i,j)+wFld(i-1,j))*halfRL
457       &         uTrans*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))*halfRL  c            flx_EW(i,j)=
458    c    &         uTrans*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))*halfRL
459             ENDDO             ENDDO
460            ENDDO            ENDDO
461  C Advective Flux on Southern face  C Advective Flux on Southern face
# Line 443  C Advective Flux on Southern face Line 463  C Advective Flux on Southern face
463             DO i=iMin,iMax             DO i=iMin,iMax
464  C     transport through Southern face area:  C     transport through Southern face area:
465               vTrans = (               vTrans = (
466       &          drF(k-1)*_hFacS(i,j,k-1,bi,bj)*vVel(i,j,k-1,bi,bj)       &          drF(km1)*_hFacS(i,j,km1,bi,bj)*vVel(i,j,km1,bi,bj)
467       &                  *rhoFacC(k-1)       &                  *rhoFacC(km1)*mskM1
468       &         +drF( k )*_hFacS(i,j, k ,bi,bj)*vVel(i,j, k ,bi,bj)       &         +drF( k )*_hFacS(i,j, k ,bi,bj)*vVel(i,j, k ,bi,bj)
469       &                  *rhoFacC(k)       &                  *rhoFacC(k)
470       &                )*halfRL*_dxG(i,j,bi,bj)*deepFacF(k)       &                )*halfRL*_dxG(i,j,bi,bj)*deepFacF(k)
471               flx_NS(i,j)=               flx_NS(i,j) = vTrans*(wFld(i,j)+wFld(i,j-1))*halfRL
472       &         vTrans*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))*halfRL  c            flx_NS(i,j)=
473    c    &         vTrans*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))*halfRL
474             ENDDO             ENDDO
475            ENDDO            ENDDO
476             ENDIF
477  C Advective Flux on Lower face of W-Cell (= at tracer-cell center, level k)  C Advective Flux on Lower face of W-Cell (= at tracer-cell center, level k)
478    c        IF (.TRUE.) THEN
479            DO j=jMin,jMax            DO j=jMin,jMax
480             DO i=iMin,iMax             DO i=iMin,iMax
481  C     NH in p-coord.: advect wSpeed [m/s] with rTrans  C     NH in p-coord.: advect wSpeed [m/s] with rTrans
482               tmp_WbarZ = halfRL*               tmp_WbarZ = halfRL*
483       &              ( wVel(i,j, k ,bi,bj)*rVel2wUnit(k)       &              ( wVel(i,j, k ,bi,bj)*rVel2wUnit( k )
484       &               +wVel(i,j,kp1,bi,bj)*rVel2wUnit(kp1)*wOverRide )       &               +wVel(i,j,kp1,bi,bj)*rVel2wUnit(kp1)*mskP1 )
485  C     transport through Lower face area:  C     transport through Lower face area:
486               rTrans = halfRL*               rTrans = halfRL*
487       &              ( wVel(i,j, k ,bi,bj)*deepFac2F( k )*rhoFacF( k )       &              ( wVel(i,j, k ,bi,bj)*deepFac2F( k )*rhoFacF( k )
488       &               +wVel(i,j,kp1,bi,bj)*deepFac2F(kp1)*rhoFacF(kp1)       &               +wVel(i,j,kp1,bi,bj)*deepFac2F(kp1)*rhoFacF(kp1)
489       &                                   *wOverRide       &                                   *mskP1
490       &              )*rA(i,j,bi,bj)       &              )*rA(i,j,bi,bj)
491               flx_Dn(i,j) = rTrans*tmp_WbarZ               flx_Dn(i,j) = rTrans*tmp_WbarZ
492             ENDDO             ENDDO
493            ENDDO            ENDDO
494            IF ( k.EQ.2 ) THEN  c        ENDIF
495  C Advective Flux on Upper face of W-Cell (= at tracer-cell center, level k-1)           IF ( k.EQ.1 .AND. selectNHfreeSurf.GE.1 ) THEN
496    C Advective Flux on Upper face of W-Cell (= at surface)
497             DO j=jMin,jMax             DO j=jMin,jMax
498              DO i=iMin,iMax              DO i=iMin,iMax
499               tmp_WbarZ = halfRL*               tmp_WbarZ = wVel(i,j,k,bi,bj)*rVel2wUnit(k)
500       &              ( wVel(i,j,k-1,bi,bj)*rVel2wUnit(k-1)*surfFac               rTrans = wVel(i,j,k,bi,bj)*deepFac2F(k)*rhoFacF(k)
501       &               +wVel(i,j, k, bi,bj)*rVel2wUnit( k ) )       &               *rA(i,j,bi,bj)
              rTrans = halfRL*  
      &              ( wVel(i,j,k-1,bi,bj)*deepFac2F(k-1)*rhoFacF(k-1)  
      &                                   *surfFac  
      &               +wVel(i,j, k ,bi,bj)*deepFac2F( k )*rhoFacF( k )  
      &              )*rA(i,j,bi,bj)  
502               flxAdvUp(i,j) = rTrans*tmp_WbarZ               flxAdvUp(i,j) = rTrans*tmp_WbarZ
 C to recover old (before 2009/11/30) results (since flxAdvUp(k=2) was zero)  
503  c            flxAdvUp(i,j) = 0.  c            flxAdvUp(i,j) = 0.
504              ENDDO              ENDDO
505             ENDDO             ENDDO
506            ENDIF           ENDIF
507             IF ( k.GT.1 .OR. selectNHfreeSurf.GE.1 ) THEN
508  C     Tendency is minus divergence of advective fluxes:  C     Tendency is minus divergence of advective fluxes:
509  C     anelastic: all transports & advect. fluxes are scaled by rhoFac  C     anelastic: all transports & advect. fluxes are scaled by rhoFac
510            DO j=jMin,jMax            DO j=jMin,jMax
# Line 496  C     anelastic: all transports & advect Line 515  C     anelastic: all transports & advect
515       &           + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign*wUnit2rVel(k)       &           + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign*wUnit2rVel(k)
516       &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)       &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
517       &          *recip_deepFac2F(k)*recip_rhoFacF(k)       &          *recip_deepFac2F(k)*recip_rhoFacF(k)
518               ENDDO
519              ENDDO
520             ENDIF
521    
522             DO j=jMin,jMax
523               DO i=iMin,iMax
524  C--          prepare for next level (k+1)  C--          prepare for next level (k+1)
525               flxAdvUp(i,j)=flx_Dn(i,j)               flxAdvUp(i,j)=flx_Dn(i,j)
526             ENDDO             ENDDO
527            ENDDO           ENDDO
528    
529    c       ELSE
530    C-    if momAdvection / else
531    c         DO j=1-OLy,sNy+OLy
532    c          DO i=1-OLx,sNx+OLx
533    c            gW(i,j,k,bi,bj) = 0. _d 0
534    c          ENDDO
535    c         ENDDO
536    
537    C-    endif momAdvection.
538          ENDIF          ENDIF
539    
540          IF ( useNHMTerms ) THEN  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
541    
542            IF ( useNHMTerms .AND. k.GT.1 ) THEN
543            CALL MOM_W_METRIC_NH(            CALL MOM_W_METRIC_NH(
544       I               bi,bj,k,       I               bi,bj,k,
545       I               uVel, vVel,       I               uVel, vVel,
# Line 514  C--          prepare for next level (k+1 Line 551  C--          prepare for next level (k+1
551             ENDDO             ENDDO
552            ENDDO            ENDDO
553          ENDIF          ENDIF
554          IF ( use3dCoriolis ) THEN          IF ( use3dCoriolis .AND. k.GT.1 ) THEN
555            CALL MOM_W_CORIOLIS_NH(            CALL MOM_W_CORIOLIS_NH(
556       I               bi,bj,k,       I               bi,bj,k,
557       I               uVel, vVel,       I               uVel, vVel,
# Line 535  C---+----1----+----2----+----3----+----4 Line 572  C---+----1----+----2----+----3----+----4
572       &                           k, 1, 2, bi,bj, myThid )       &                           k, 1, 2, bi,bj, myThid )
573  C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL  C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
574  C        does it only if k=1 (never the case here)  C        does it only if k=1 (never the case here)
575            IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Diss ',bi,bj,myThid)  c         IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Diss ',bi,bj,myThid)
576          ENDIF          ENDIF
577          IF ( diagAdvec ) THEN          IF ( diagAdvec ) THEN
578            CALL DIAGNOSTICS_FILL( gW, 'Wm_Advec',            CALL DIAGNOSTICS_FILL( gW, 'Wm_Advec',
579       &                           k,Nr, 1, bi,bj, myThid )       &                           k,Nr, 1, bi,bj, myThid )
580            IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Advec',bi,bj,myThid)  c         IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Advec',bi,bj,myThid)
581          ENDIF          ENDIF
582  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
583    
Legend:
Removed from v.1.42  
changed lines
  Added in v.1.43
  ViewVC Help
Powered by ViewVC 1.1.22