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

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

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

revision 1.7 by adcroft, Tue Dec 15 00:20:34 1998 UTC revision 1.14 by ljmc, Mon Jun 25 20:38:15 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    #undef ALLOW_ZONAL_FILT
6    #undef ALLOW_SHAP_FILT
7    
8  C     /==========================================================\  C     /==========================================================\
9  C     | S/R CALC_DIV_GHAT                                        |  C     | S/R CALC_DIV_GHAT                                        |
10  C     | o Form the right hand-side of the surface pressure eqn.  |  C     | o Form the right hand-side of the surface pressure eqn.  |
11    C     |==========================================================|
12  C     \==========================================================/  C     \==========================================================/
13        SUBROUTINE CALC_DIV_GHAT(        SUBROUTINE CALC_DIV_GHAT(
14       I        bi,bj,iMin,iMax,jMin,jMax,       I        bi,bj,iMin,iMax,jMin,jMax,K,
      I        K,  
15       I        xA,yA,       I        xA,yA,
16         U        cg2d_b,
17       I        myThid)       I        myThid)
18    
19        IMPLICIT NONE        IMPLICIT NONE
# Line 21  C     == Global variables == Line 25  C     == Global variables ==
25  #include "EEPARAMS.h"  #include "EEPARAMS.h"
26  #include "PARAMS.h"  #include "PARAMS.h"
27  #include "GRID.h"  #include "GRID.h"
28  #include "CG2D.h"  #ifdef ALLOW_NONHYDROSTATIC
29    #include "CG3D.h"
30    #endif
31    
32  C     == Routine arguments ==  C     == Routine arguments ==
 C     pH - Hydrostatic pressure  
33  C     bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation  C     bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation
34  C                                      results will be set.  C                                      results will be set.
35  C     kUp, kDown, kM1                - Index for upper and lower layers.  C     k                              - Index of layer.
36  C     myThid - Instance number for this innvocation of CALC_MOM_RHS  C     xA, yA                         - Cell face areas
37    C     cg2d_b - Conjugate Gradient 2-D solver : Right-hand side vector
38    C     myThid - Instance number for this call of CALC_DIV_GHAT
39        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
40        INTEGER K        INTEGER K
41        _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
42        _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
43          _RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
44        INTEGER myThid        INTEGER myThid
45    
46  C     == Local variables ==  C     == Local variables ==
# Line 44  C     Term is the vertical integral of t Line 52  C     Term is the vertical integral of t
52  C     time tendency terms along with a relaxation term that  C     time tendency terms along with a relaxation term that
53  C     pulls div(U) + dh/dt back toward zero.  C     pulls div(U) + dh/dt back toward zero.
54    
55        IF ( k .EQ. Nr ) THEN        IF (implicDiv2Dflow.EQ.1.) then
56  C      Initialise source term on first pass  C     Fully Implicit treatment of the Barotropic Flow Divergence
57         DO j=jMin,jMax          DO j=1,sNy
58          DO i=iMin,iMax           DO i=1,sNx+1
59  C        cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)            pf(i,j) = xA(i,j)*gUNm1(i,j,k,bi,bj) / deltaTmom
60  C    &    -freeSurfFac*_rA(i,j,bi,bj)*           ENDDO
61  C    &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTMom          ENDDO
62           cg2d_b(i,j,bi,bj) =        ELSE
63       &     freeSurfFac*_rA(i,j,bi,bj)*horiVertRatio*(  C     Explicit+Implicit part of the Barotropic Flow Divergence
64       &      -cg2d_x(I,J,bi,bj)/deltaTMom/deltaTMom  C      => Filtering of uVel,vVel is necessary
65  #ifdef USE_NATURAL_BCS  #ifdef ALLOW_ZONAL_FILT
66       &      +EmPmR(I,J,bi,bj)/deltaTMom          IF (zonal_filt_lat.LT.90.) THEN
67              CALL ZONAL_FILTER(
68         &      uVel, hFacW, 1-1, sNy+1, k, k, bi, bj, 1, myThid)
69              CALL ZONAL_FILTER(
70         &      vVel, hFacS, 1-1, sNy+1, k, k, bi, bj, 2, myThid)
71            ENDIF
72  #endif  #endif
73       &     )          DO j=1,sNy
74             DO i=1,sNx+1
75              pf(i,j) = ( implicDiv2Dflow * gUNm1(i,j,k,bi,bj)
76         &          + (1.-implicDiv2Dflow) * uVel(i,j,k,bi,bj)
77         &               ) * xA(i,j) / deltaTmom
78             ENDDO
79          ENDDO          ENDDO
        ENDDO  
80        ENDIF        ENDIF
   
       DO j=jMin,jMax  
        DO i=iMin,iMax  
         pf(i,j) = xA(i,j)*gUNm1(i,j,k,bi,bj) / deltaTmom  
        ENDDO  
       ENDDO  
81        DO j=1,sNy        DO j=1,sNy
82         DO i=1,sNx         DO i=1,sNx
83          cg2d_b(i,j,bi,bj) =  cg2d_b(i,j,bi,bj) +          cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj) +
84       &   pf(i+1,j)-pf(i,j)       &   pf(i+1,j)-pf(i,j)
85         ENDDO         ENDDO
86        ENDDO        ENDDO
87    
88        DO j=jMin,jMax  #ifdef ALLOW_NONHYDROSTATIC
89         DO i=iMin,iMax        IF (nonHydrostatic) THEN
90          pf(i,j) = yA(i,j)*gVNm1(i,j,k,bi,bj) / deltatmom         DO j=1,sNy
91            DO i=1,sNx
92             cg3d_b(i,j,k,bi,bj) = pf(i+1,j)-pf(i,j)
93            ENDDO
94         ENDDO         ENDDO
95        ENDDO        ENDIF
96    #endif
97    
98          IF (implicDiv2Dflow.EQ.1.) then
99    C     Fully Implicit treatment of the Barotropic Flow Divergence
100            DO j=1,sNy+1
101             DO i=1,sNx
102              pf(i,j) = yA(i,j)*gVNm1(i,j,k,bi,bj) / deltatmom
103             ENDDO
104            ENDDO
105          ELSE
106    C     Explicit+Implicit part of the Barotropic Flow Divergence
107            DO j=1,sNy+1
108             DO i=1,sNx
109              pf(i,j) = ( implicDiv2Dflow * gVNm1(i,j,k,bi,bj)
110         &          + (1.-implicDiv2Dflow) * vVel(i,j,k,bi,bj)
111         &               ) * yA(i,j) / deltaTmom
112             ENDDO
113            ENDDO
114          ENDIF
115        DO j=1,sNy        DO j=1,sNy
116         DO i=1,sNx         DO i=1,sNx
117          cg2d_b(i,j,bi,bj)   =  cg2d_b(i,j,bi,bj) +          cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj) +
118       &   pf(i,j+1)-pf(i,j)       &   pf(i,j+1)-pf(i,j)
119         ENDDO         ENDDO
120        ENDDO        ENDDO
121    
122    #ifdef ALLOW_NONHYDROSTATIC
123          IF (nonHydrostatic) THEN
124           DO j=1,sNy
125            DO i=1,sNx
126             cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj) +
127         &    pf(i,j+1)-pf(i,j)
128            ENDDO
129           ENDDO
130          ENDIF
131    #endif
132    
133        RETURN        RETURN
134        END        END
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.22