/[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.35 by jmc, Tue Dec 5 05:25:08 2006 UTC revision 1.36 by jmc, Mon Mar 12 23:48:24 2007 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    #define CALC_GW_NEW_THICK
6    
7  CBOP  CBOP
8  C     !ROUTINE: CALC_GW  C     !ROUTINE: CALC_GW
# Line 56  C     iMin,iMax Line 57  C     iMin,iMax
57  C     jMin,jMax  C     jMin,jMax
58  C     xA            :: W-Cell face area normal to X  C     xA            :: W-Cell face area normal to X
59  C     yA            :: W-Cell face area normal to Y  C     yA            :: W-Cell face area normal to Y
60    C     rMinW,rMaxW   :: column boundaries (r-units) at Western  Edge
61    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 70  C     i,j,k         :: Loop counters Line 73  C     i,j,k         :: Loop counters
73        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
74        _RS    xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS    xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75        _RS    yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS    yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76          _RS    rMinW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77          _RS    rMaxW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78          _RS    rMinS (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79          _RS    rMaxS (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80        _RL    rThickC_W    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    rThickC_W    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81        _RL    rThickC_S    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    rThickC_S    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82        _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 98  C     i,j,k         :: Loop counters
98        _RS  halfRS, zeroRS        _RS  halfRS, zeroRS
99        PARAMETER( halfRL = 0.5D0 )        PARAMETER( halfRL = 0.5D0 )
100        PARAMETER( halfRS = 0.5 , zeroRS = 0. )        PARAMETER( halfRS = 0.5 , zeroRS = 0. )
101          PARAMETER( iMin = 1 , iMax = sNx )
102          PARAMETER( jMin = 1 , jMax = sNy )
103  CEOP  CEOP
104    
105  C Catch barotropic mode  C Catch barotropic mode
106        IF ( Nr .LT. 2 ) RETURN        IF ( Nr .LT. 2 ) RETURN
107    
       iMin = 1  
       iMax = sNx  
       jMin = 1  
       jMax = sNy  
   
108  C--   Initialise gW to zero  C--   Initialise gW to zero
109        DO k=1,Nr        DO k=1,Nr
110          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
# Line 123  C--   Boundaries condition at top Line 127  C--   Boundaries condition at top
127           flxDisUp(i,j) = 0.           flxDisUp(i,j) = 0.
128         ENDDO         ENDDO
129        ENDDO        ENDDO
130    C--   column boundaries :
131          IF (momViscosity) THEN
132            DO j=1-Oly,sNy+Oly
133             DO i=1-Olx+1,sNx+Olx
134               rMaxW(i,j) = MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
135               rMinW(i,j) = MAX(   R_low(i-1,j,bi,bj),   R_low(i,j,bi,bj) )
136             ENDDO
137            ENDDO
138            DO j=1-Oly+1,sNy+Oly
139             DO i=1-Olx,sNx+Olx
140               rMaxS(i,j) = MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
141               rMinS(i,j) = MAX(   R_low(i,j-1,bi,bj),   R_low(i,j,bi,bj) )
142             ENDDO
143            ENDDO
144          ENDIF
145    
146  C---  Sweep down column  C---  Sweep down column
147        DO k=2,Nr        DO k=2,Nr
# Line 133  C---  Sweep down column Line 152  C---  Sweep down column
152            wOverRide=0.            wOverRide=0.
153          ENDIF          ENDIF
154  C--   Compute grid factor arround a W-point:  C--   Compute grid factor arround a W-point:
155    #ifdef CALC_GW_NEW_THICK
156            DO j=1-Oly,sNy+Oly
157             DO i=1-Olx,sNx+Olx
158               IF ( maskC(i,j,k-1,bi,bj).EQ.0. .OR.
159         &          maskC(i,j, k ,bi,bj).EQ.0. ) THEN
160                 recip_rThickC(i,j) = 0.
161               ELSE
162    C-    valid in z & p coord.; also accurate if Interface @ middle between 2 centers
163                 recip_rThickC(i,j) = 1. _d 0 /
164         &        (  MIN( Ro_surf(i,j,bi,bj),rC(k-1) )
165         &         - MAX( R_low(i,j,bi,bj),  rC(k)   )
166         &        )
167               ENDIF
168             ENDDO
169            ENDDO
170            IF (momViscosity) THEN
171             DO j=1-Oly,sNy+Oly
172              DO i=1-Olx,sNx+Olx
173               rThickC_C(i,j) = MAX( zeroRS,
174         &                           MIN( Ro_surf(i,j,bi,bj), rC(k-1) )
175         &                          -MAX(   R_low(i,j,bi,bj),  rC(k)  )
176         &                         )
177              ENDDO
178             ENDDO
179             DO j=1-Oly,sNy+Oly
180              DO i=1-Olx+1,sNx+Olx
181               rThickC_W(i,j) = MAX( zeroRS,
182         &                           MIN( rMaxW(i,j), rC(k-1) )
183         &                          -MAX( rMinW(i,j), rC(k)   )
184         &                         )
185    C     W-Cell Western face area:
186               xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
187    c    &                              *deepFacF(k)
188              ENDDO
189             ENDDO
190             DO j=1-Oly+1,sNy+Oly
191              DO i=1-Olx,sNx+Olx
192               rThickC_S(i,j) = MAX( zeroRS,
193         &                           MIN( rMaxS(i,j), rC(k-1) )
194         &                          -MAX( rMinS(i,j), rC(k)   )
195         &                         )
196    C     W-Cell Southern face area:
197               yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)
198    c    &                              *deepFacF(k)
199    C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC.
200    C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)
201              ENDDO
202             ENDDO
203            ENDIF
204    #else /* CALC_GW_NEW_THICK */
205          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
206           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
207  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 !
# Line 171  C this gives deepFacF*recip_deepFacF => Line 240  C this gives deepFacF*recip_deepFacF =>
240  c          ENDIF  c          ENDIF
241           ENDDO           ENDDO
242          ENDDO          ENDDO
243    #endif /* CALC_GW_NEW_THICK */
244    
245  C--   horizontal bi-harmonic dissipation  C--   horizontal bi-harmonic dissipation
246          IF (momViscosity .AND. viscA4W.NE.0. ) THEN          IF (momViscosity .AND. viscA4W.NE.0. ) THEN
# Line 227  CML ************************************ Line 297  CML ************************************
297  CML   No-slip Boundary conditions for bi-harmonic dissipation  CML   No-slip Boundary conditions for bi-harmonic dissipation
298  CML   need to be implemented here!  CML   need to be implemented here!
299  CML ************************************************************  CML ************************************************************
300          ELSE          ELSEIF (momViscosity) THEN
301  C-    Initialize del2w to zero:  C-    Initialize del2w to zero:
302            DO j=1-Oly,sNy+Oly            DO j=1-Oly,sNy+Oly
303             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx
# Line 322  C--        prepare for next level (k+1) Line 392  C--        prepare for next level (k+1)
392            ENDDO            ENDDO
393          ENDIF          ENDIF
394    
395          IF (no_slip_sides) THEN          IF ( momViscosity .AND. no_slip_sides ) THEN
396  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
397            CALL MOM_W_SIDEDRAG(            CALL MOM_W_SIDEDRAG(
398       I               bi,bj,k,       I               bi,bj,k,
# Line 374  cOld &                )*halfRL Line 444  cOld &                )*halfRL
444  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)
445            DO j=jMin,jMax            DO j=jMin,jMax
446             DO i=iMin,iMax             DO i=iMin,iMax
447               tmp_WbarZ = halfRL*( wVel(i,j, k ,bi,bj)  C     NH in p-coord.: advect wSpeed [m/s] with rTrans
448       &                           +wVel(i,j,kp1,bi,bj)*wOverRide )               tmp_WbarZ = halfRL*
449         &              ( wVel(i,j, k ,bi,bj)*rVel2wUnit(k)
450         &               +wVel(i,j,kp1,bi,bj)*rVel2wUnit(kp1)*wOverRide )
451  C     transport through Lower face area:  C     transport through Lower face area:
452               rTrans = halfRL*               rTrans = halfRL*
453       &              ( wVel(i,j, k ,bi,bj)*deepFac2F( k )*rhoFacF( k )       &              ( wVel(i,j, k ,bi,bj)*deepFac2F( k )*rhoFacF( k )
# Line 393  C     anelastic: all transports & advect Line 465  C     anelastic: all transports & advect
465               gW(i,j,k,bi,bj) =               gW(i,j,k,bi,bj) =
466       &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )       &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
467       &           + ( flx_NS(i,j+1)-flx_NS(i,j) )       &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
468       &           + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign       &           + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign*wUnit2rVel(k)
469       &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)       &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
470       &          *recip_deepFac2F(k)*recip_rhoFacF(k)       &          *recip_deepFac2F(k)*recip_rhoFacF(k)
471  cOld         gW(i,j,k,bi,bj) =  cOld         gW(i,j,k,bi,bj) =
Legend:
Removed from v.1.35  
changed lines
  Added in v.1.36
  ViewVC Help
Powered by ViewVC 1.1.22