/[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.37 by jmc, Fri Oct 19 14:41:39 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 29  C     == Global variables == Line 30  C     == Global variables ==
30  #include "EEPARAMS.h"  #include "EEPARAMS.h"
31  #include "PARAMS.h"  #include "PARAMS.h"
32  #include "GRID.h"  #include "GRID.h"
33    #include "RESTART.h"
34  #include "SURFACE.h"  #include "SURFACE.h"
35  #include "DYNVARS.h"  #include "DYNVARS.h"
36  #include "NH_VARS.h"  #include "NH_VARS.h"
# Line 56  C     iMin,iMax Line 58  C     iMin,iMax
58  C     jMin,jMax  C     jMin,jMax
59  C     xA            :: W-Cell face area normal to X  C     xA            :: W-Cell face area normal to X
60  C     yA            :: W-Cell face area normal to Y  C     yA            :: W-Cell face area normal to Y
61    C     rMinW,rMaxW   :: column boundaries (r-units) at Western  Edge
62    C     rMinS,rMaxS   :: column boundaries (r-units) at Southern Edge
63  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
64  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
65  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 74  C     i,j,k         :: Loop counters
74        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
75        _RS    xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS    xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76        _RS    yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS    yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77          _RS    rMinW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78          _RS    rMaxW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79          _RS    rMinS (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80          _RS    rMaxS (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81        _RL    rThickC_W    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    rThickC_W    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82        _RL    rThickC_S    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    rThickC_S    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83        _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 99  C     i,j,k         :: Loop counters
99        _RS  halfRS, zeroRS        _RS  halfRS, zeroRS
100        PARAMETER( halfRL = 0.5D0 )        PARAMETER( halfRL = 0.5D0 )
101        PARAMETER( halfRS = 0.5 , zeroRS = 0. )        PARAMETER( halfRS = 0.5 , zeroRS = 0. )
102          PARAMETER( iMin = 1 , iMax = sNx )
103          PARAMETER( jMin = 1 , jMax = sNy )
104  CEOP  CEOP
105    
106  C Catch barotropic mode  C Catch barotropic mode
107        IF ( Nr .LT. 2 ) RETURN        IF ( Nr .LT. 2 ) RETURN
108    
       iMin = 1  
       iMax = sNx  
       jMin = 1  
       jMax = sNy  
   
109  C--   Initialise gW to zero  C--   Initialise gW to zero
110        DO k=1,Nr        DO k=1,Nr
111          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
# Line 123  C--   Boundaries condition at top Line 128  C--   Boundaries condition at top
128           flxDisUp(i,j) = 0.           flxDisUp(i,j) = 0.
129         ENDDO         ENDDO
130        ENDDO        ENDDO
131    C--   column boundaries :
132          IF (momViscosity) THEN
133            DO j=1-Oly,sNy+Oly
134             DO i=1-Olx+1,sNx+Olx
135               rMaxW(i,j) = MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
136               rMinW(i,j) = MAX(   R_low(i-1,j,bi,bj),   R_low(i,j,bi,bj) )
137             ENDDO
138            ENDDO
139            DO j=1-Oly+1,sNy+Oly
140             DO i=1-Olx,sNx+Olx
141               rMaxS(i,j) = MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
142               rMinS(i,j) = MAX(   R_low(i,j-1,bi,bj),   R_low(i,j,bi,bj) )
143             ENDDO
144            ENDDO
145          ENDIF
146    
147  C---  Sweep down column  C---  Sweep down column
148        DO k=2,Nr        DO k=2,Nr
# Line 133  C---  Sweep down column Line 153  C---  Sweep down column
153            wOverRide=0.            wOverRide=0.
154          ENDIF          ENDIF
155  C--   Compute grid factor arround a W-point:  C--   Compute grid factor arround a W-point:
156    #ifdef CALC_GW_NEW_THICK
157            DO j=1-Oly,sNy+Oly
158             DO i=1-Olx,sNx+Olx
159               IF ( maskC(i,j,k-1,bi,bj).EQ.0. .OR.
160         &          maskC(i,j, k ,bi,bj).EQ.0. ) THEN
161                 recip_rThickC(i,j) = 0.
162               ELSE
163    C-    valid in z & p coord.; also accurate if Interface @ middle between 2 centers
164                 recip_rThickC(i,j) = 1. _d 0 /
165         &        (  MIN( Ro_surf(i,j,bi,bj),rC(k-1) )
166         &         - MAX( R_low(i,j,bi,bj),  rC(k)   )
167         &        )
168               ENDIF
169             ENDDO
170            ENDDO
171            IF (momViscosity) THEN
172             DO j=1-Oly,sNy+Oly
173              DO i=1-Olx,sNx+Olx
174               rThickC_C(i,j) = MAX( zeroRS,
175         &                           MIN( Ro_surf(i,j,bi,bj), rC(k-1) )
176         &                          -MAX(   R_low(i,j,bi,bj),  rC(k)  )
177         &                         )
178              ENDDO
179             ENDDO
180             DO j=1-Oly,sNy+Oly
181              DO i=1-Olx+1,sNx+Olx
182               rThickC_W(i,j) = MAX( zeroRS,
183         &                           MIN( rMaxW(i,j), rC(k-1) )
184         &                          -MAX( rMinW(i,j), rC(k)   )
185         &                         )
186    C     W-Cell Western face area:
187               xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
188    c    &                              *deepFacF(k)
189              ENDDO
190             ENDDO
191             DO j=1-Oly+1,sNy+Oly
192              DO i=1-Olx,sNx+Olx
193               rThickC_S(i,j) = MAX( zeroRS,
194         &                           MIN( rMaxS(i,j), rC(k-1) )
195         &                          -MAX( rMinS(i,j), rC(k)   )
196         &                         )
197    C     W-Cell Southern face area:
198               yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)
199    c    &                              *deepFacF(k)
200    C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC.
201    C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)
202              ENDDO
203             ENDDO
204            ENDIF
205    #else /* CALC_GW_NEW_THICK */
206          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
207           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
208  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 241  C this gives deepFacF*recip_deepFacF =>
241  c          ENDIF  c          ENDIF
242           ENDDO           ENDDO
243          ENDDO          ENDDO
244    #endif /* CALC_GW_NEW_THICK */
245    
246  C--   horizontal bi-harmonic dissipation  C--   horizontal bi-harmonic dissipation
247          IF (momViscosity .AND. viscA4W.NE.0. ) THEN          IF (momViscosity .AND. viscA4W.NE.0. ) THEN
# Line 227  CML ************************************ Line 298  CML ************************************
298  CML   No-slip Boundary conditions for bi-harmonic dissipation  CML   No-slip Boundary conditions for bi-harmonic dissipation
299  CML   need to be implemented here!  CML   need to be implemented here!
300  CML ************************************************************  CML ************************************************************
301          ELSE          ELSEIF (momViscosity) THEN
302  C-    Initialize del2w to zero:  C-    Initialize del2w to zero:
303            DO j=1-Oly,sNy+Oly            DO j=1-Oly,sNy+Oly
304             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx
# Line 322  C--        prepare for next level (k+1) Line 393  C--        prepare for next level (k+1)
393            ENDDO            ENDDO
394          ENDIF          ENDIF
395    
396          IF (no_slip_sides) THEN          IF ( momViscosity .AND. no_slip_sides ) THEN
397  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
398            CALL MOM_W_SIDEDRAG(            CALL MOM_W_SIDEDRAG(
399       I               bi,bj,k,       I               bi,bj,k,
# Line 374  cOld &                )*halfRL Line 445  cOld &                )*halfRL
445  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)
446            DO j=jMin,jMax            DO j=jMin,jMax
447             DO i=iMin,iMax             DO i=iMin,iMax
448               tmp_WbarZ = halfRL*( wVel(i,j, k ,bi,bj)  C     NH in p-coord.: advect wSpeed [m/s] with rTrans
449       &                           +wVel(i,j,kp1,bi,bj)*wOverRide )               tmp_WbarZ = halfRL*
450         &              ( wVel(i,j, k ,bi,bj)*rVel2wUnit(k)
451         &               +wVel(i,j,kp1,bi,bj)*rVel2wUnit(kp1)*wOverRide )
452  C     transport through Lower face area:  C     transport through Lower face area:
453               rTrans = halfRL*               rTrans = halfRL*
454       &              ( 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 466  C     anelastic: all transports & advect
466               gW(i,j,k,bi,bj) =               gW(i,j,k,bi,bj) =
467       &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )       &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
468       &           + ( flx_NS(i,j+1)-flx_NS(i,j) )       &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
469       &           + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign       &           + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign*wUnit2rVel(k)
470       &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)       &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
471       &          *recip_deepFac2F(k)*recip_rhoFacF(k)       &          *recip_deepFac2F(k)*recip_rhoFacF(k)
472  cOld         gW(i,j,k,bi,bj) =  cOld         gW(i,j,k,bi,bj) =
# Line 451  c#ifdef ALLOW_ADAMSBASHFORTH_3 Line 524  c#ifdef ALLOW_ADAMSBASHFORTH_3
524  c       CALL ADAMS_BASHFORTH3(  c       CALL ADAMS_BASHFORTH3(
525  c    I                         bi, bj, k,  c    I                         bi, bj, k,
526  c    U                         gW, gwNm,  c    U                         gW, gwNm,
527  c    I                         momStartAB, myIter, myThid )  c    I                         nHydStartAB, myIter, myThid )
528  c#else /* ALLOW_ADAMSBASHFORTH_3 */  c#else /* ALLOW_ADAMSBASHFORTH_3 */
529          CALL ADAMS_BASHFORTH2(          CALL ADAMS_BASHFORTH2(
530       I                         bi, bj, k,       I                         bi, bj, k,
531       U                         gW, gwNm1,       U                         gW, gwNm1,
532       I                         myIter, myThid )       I                         nHydStartAB, myIter, myThid )
533  c#endif /* ALLOW_ADAMSBASHFORTH_3 */  c#endif /* ALLOW_ADAMSBASHFORTH_3 */
534    
535  C--   Dissipation term outside the Adams-Bashforth:  C--   Dissipation term outside the Adams-Bashforth:
Legend:
Removed from v.1.35  
changed lines
  Added in v.1.37
  ViewVC Help
Powered by ViewVC 1.1.22