/[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.5 by cnh, Wed Oct 28 03:11:36 1998 UTC revision 1.16 by jmc, Wed Sep 19 13:58:08 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6  C     /==========================================================\  C     /==========================================================\
7  C     | S/R CALC_DIV_GHAT                                        |  C     | S/R CALC_DIV_GHAT                                        |
8  C     | o Form the right hand-side of the surface pressure eqn.  |  C     | o Form the right hand-side of the surface pressure eqn.  |
9    C     |==========================================================|
10  C     \==========================================================/  C     \==========================================================/
11        SUBROUTINE CALC_DIV_GHAT(        SUBROUTINE CALC_DIV_GHAT(
12       I        bi,bj,iMin,iMax,jMin,jMax,       I        bi,bj,iMin,iMax,jMin,jMax,K,
      I        K,  
13       I        xA,yA,       I        xA,yA,
14         U        cg2d_b,
15       I        myThid)       I        myThid)
16    
17        IMPLICIT NONE        IMPLICIT NONE
# Line 21  C     == Global variables == Line 23  C     == Global variables ==
23  #include "EEPARAMS.h"  #include "EEPARAMS.h"
24  #include "PARAMS.h"  #include "PARAMS.h"
25  #include "GRID.h"  #include "GRID.h"
26  #include "CG2D.h"  #include "SOLVE_FOR_PRESSURE3D.h"
27    
28  C     == Routine arguments ==  C     == Routine arguments ==
 C     pH - Hydrostatic pressure  
29  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
30  C                                      results will be set.  C                                      results will be set.
31  C     kUp, kDown, kM1                - Index for upper and lower layers.  C     k                              - Index of layer.
32  C     myThid - Instance number for this innvocation of CALC_MOM_RHS  C     xA, yA                         - Cell face areas
33    C     cg2d_b - Conjugate Gradient 2-D solver : Right-hand side vector
34    C     myThid - Instance number for this call of CALC_DIV_GHAT
35        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
36        INTEGER K        INTEGER K
37        _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
38        _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
39          _RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
40        INTEGER myThid        INTEGER myThid
41    
42  C     == Local variables ==  C     == Local variables ==
# Line 44  C     Term is the vertical integral of t Line 48  C     Term is the vertical integral of t
48  C     time tendency terms along with a relaxation term that  C     time tendency terms along with a relaxation term that
49  C     pulls div(U) + dh/dt back toward zero.  C     pulls div(U) + dh/dt back toward zero.
50    
51        IF ( k .EQ. Nr ) THEN        IF (implicDiv2Dflow.EQ.1.) THEN
52  C      Initialise source term on first pass  C     Fully Implicit treatment of the Barotropic Flow Divergence
53         DO j=jMin,jMax          DO j=1,sNy
54          DO i=iMin,iMax           DO i=1,sNx+1
55  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
56  C    &    -freeSurfFac*_rA(i,j,bi,bj)*           ENDDO
57  C    &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTMom          ENDDO
58           cg2d_b(i,j,bi,bj) =        ELSEIF (exactConserv) THEN
59       &    -freeSurfFac*_rA(i,j,bi,bj)*horiVertRatio*  c     ELSEIF (nonlinFreeSurf.GT.0) THEN
60       &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTMom  C     Implicit treatment of the Barotropic Flow Divergence
61            DO j=1,sNy
62             DO i=1,sNx+1
63              pf(i,j) = implicDiv2Dflow
64         &             *xA(i,j)*gUNm1(i,j,k,bi,bj) / deltaTmom
65             ENDDO
66            ENDDO
67          ELSE
68    C     Explicit+Implicit part of the Barotropic Flow Divergence
69    C      => Filtering of uVel,vVel is necessary
70    C-- Now the filter are applied in the_correction_step().
71    C   We have left this code here to indicate where the filters used to be
72    C   in the algorithm before JMC moved them to after the pressure solver.
73    C-
74    C#ifdef ALLOW_ZONAL_FILT
75    C        IF (zonal_filt_lat.LT.90.) THEN
76    C          CALL ZONAL_FILTER(
77    C     &      uVel, hFacW, 1-1, sNy+1, k, k, bi, bj, 1, myThid)
78    C          CALL ZONAL_FILTER(
79    C     &      vVel, hFacS, 1-1, sNy+1, k, k, bi, bj, 2, myThid)
80    C        ENDIF
81    C#endif
82            DO j=1,sNy
83             DO i=1,sNx+1
84              pf(i,j) = ( implicDiv2Dflow * gUNm1(i,j,k,bi,bj)
85         &          + (1.-implicDiv2Dflow) * uVel(i,j,k,bi,bj)
86         &               ) * xA(i,j) / deltaTmom
87             ENDDO
88          ENDDO          ENDDO
        ENDDO  
89        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  
90        DO j=1,sNy        DO j=1,sNy
91         DO i=1,sNx         DO i=1,sNx
92          cg2d_b(i,j,bi,bj) =  cg2d_b(i,j,bi,bj) +          cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj) +
93       &   pf(i+1,j)-pf(i,j)       &   pf(i+1,j)-pf(i,j)
94         ENDDO         ENDDO
95        ENDDO        ENDDO
96    
97        DO j=jMin,jMax  #ifdef ALLOW_NONHYDROSTATIC
98         DO i=iMin,iMax        IF (nonHydrostatic) THEN
99          pf(i,j) = yA(i,j)*gVNm1(i,j,k,bi,bj) / deltatmom         DO j=1,sNy
100            DO i=1,sNx
101             cg3d_b(i,j,k,bi,bj) = pf(i+1,j)-pf(i,j)
102            ENDDO
103         ENDDO         ENDDO
104        ENDDO        ENDIF
105    #endif
106    
107          IF (implicDiv2Dflow.EQ.1.) THEN
108    C     Fully Implicit treatment of the Barotropic Flow Divergence
109            DO j=1,sNy+1
110             DO i=1,sNx
111              pf(i,j) = yA(i,j)*gVNm1(i,j,k,bi,bj) / deltatmom
112             ENDDO
113            ENDDO
114          ELSEIF (exactConserv) THEN
115    c     ELSEIF (nonlinFreeSurf.GT.0) THEN
116    C     Implicit treatment of the Barotropic Flow Divergence
117            DO j=1,sNy+1
118             DO i=1,sNx
119              pf(i,j) = implicDiv2Dflow
120         &             *yA(i,j)*gVNm1(i,j,k,bi,bj) / deltatmom
121             ENDDO
122            ENDDO
123          ELSE
124    C     Explicit+Implicit part of the Barotropic Flow Divergence
125            DO j=1,sNy+1
126             DO i=1,sNx
127              pf(i,j) = ( implicDiv2Dflow * gVNm1(i,j,k,bi,bj)
128         &          + (1.-implicDiv2Dflow) * vVel(i,j,k,bi,bj)
129         &               ) * yA(i,j) / deltaTmom
130             ENDDO
131            ENDDO
132          ENDIF
133        DO j=1,sNy        DO j=1,sNy
134         DO i=1,sNx         DO i=1,sNx
135          cg2d_b(i,j,bi,bj)   =  cg2d_b(i,j,bi,bj) +          cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj) +
136       &   pf(i,j+1)-pf(i,j)       &   pf(i,j+1)-pf(i,j)
137         ENDDO         ENDDO
138        ENDDO        ENDDO
139    
140    #ifdef ALLOW_NONHYDROSTATIC
141          IF (nonHydrostatic) THEN
142           DO j=1,sNy
143            DO i=1,sNx
144             cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj) +
145         &    pf(i,j+1)-pf(i,j)
146            ENDDO
147           ENDDO
148          ENDIF
149    #endif
150    
151        RETURN        RETURN
152        END        END
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22