/[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.16 by jmc, Wed Sep 19 13:58:08 2001 UTC revision 1.27 by jmc, Fri Nov 9 22:37:05 2012 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6  C     /==========================================================\  CBOP
7  C     | S/R CALC_DIV_GHAT                                        |  C     !ROUTINE: CALC_DIV_GHAT
8  C     | o Form the right hand-side of the surface pressure eqn.  |  C     !INTERFACE:
9  C     |==========================================================|        SUBROUTINE CALC_DIV_GHAT(
10  C     \==========================================================/       I                bi,bj,k,
11        SUBROUTINE CALC_DIV_GHAT(       U                cg2d_b, cg3d_b,
12       I        bi,bj,iMin,iMax,jMin,jMax,K,       I                myThid)
13       I        xA,yA,  C     !DESCRIPTION: \bv
14       U        cg2d_b,  C     *==========================================================*
15       I        myThid)  C     | S/R CALC_DIV_GHAT
16    C     | o Form the right hand-side of the surface pressure eqn.
17    C     *==========================================================*
18    C     | Right hand side of pressure equation is divergence
19    C     | of veclocity tendency (GHAT) term along with a relaxation
20    C     | term equal to the barotropic flow field divergence.
21    C     *==========================================================*
22    C     \ev
23    
24    C     !USES:
25        IMPLICIT NONE        IMPLICIT NONE
   
26  C     == Global variables ==  C     == Global variables ==
27  #include "SIZE.h"  #include "SIZE.h"
 #include "DYNVARS.h"  
 #include "FFIELDS.h"  
28  #include "EEPARAMS.h"  #include "EEPARAMS.h"
29  #include "PARAMS.h"  #include "PARAMS.h"
30  #include "GRID.h"  #include "GRID.h"
31  #include "SOLVE_FOR_PRESSURE3D.h"  #include "DYNVARS.h"
32    #ifdef ALLOW_ADDFLUID
33    # include "FFIELDS.h"
34    #endif
35    
36    C     !INPUT/OUTPUT PARAMETERS:
37  C     == Routine arguments ==  C     == Routine arguments ==
38  C     bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation  C     bi, bj  :: tile indices
39  C                                      results will be set.  C     k       :: Index of layer.
40  C     k                              - Index of layer.  C     cg2d_b  :: Conjugate Gradient 2-D solver : Right-hand side vector
41  C     xA, yA                         - Cell face areas  C     cg3d_b  :: Conjugate Gradient 3-D solver : Right-hand side vector
42  C     cg2d_b - Conjugate Gradient 2-D solver : Right-hand side vector  C     myThid  :: Instance number for this call of CALC_DIV_GHAT
43  C     myThid - Instance number for this call of CALC_DIV_GHAT        INTEGER bi,bj
44        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER k
       INTEGER K  
       _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
45        _RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
46    #ifdef ALLOW_NONHYDROSTATIC
47          _RL cg3d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
48    #else
49          _RL cg3d_b(1)
50    #endif
51        INTEGER myThid        INTEGER myThid
52    
53    C     !LOCAL VARIABLES:
54  C     == Local variables ==  C     == Local variables ==
55    C     i,j    :: Loop counters
56    C     xA, yA :: Cell vertical face areas
57    C     pf     :: Intermediate array for building RHS source term.
58        INTEGER i,j        INTEGER i,j
59          _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
60          _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61        _RL pf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL pf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62    CEOP
63    
64    C     Calculate vertical face areas
65          DO j=1,sNy+1
66            DO i=1,sNx+1
67              xA(i,j) = _dyG(i,j,bi,bj)*deepFacC(k)
68         &              *drF(k)*_hFacW(i,j,k,bi,bj)*rhoFacC(k)
69              yA(i,j) = _dxG(i,j,bi,bj)*deepFacC(k)
70         &              *drF(k)*_hFacS(i,j,k,bi,bj)*rhoFacC(k)
71            ENDDO
72          ENDDO
73    
74  C--   Pressure equation source term  C--   Pressure equation source term
75  C     Term is the vertical integral of the divergence of the  C     Term is the vertical integral of the divergence of the
76  C     time tendency terms along with a relaxation term that  C     time tendency terms along with a relaxation term that
77  C     pulls div(U) + dh/dt back toward zero.  C     pulls div(U) + dh/dt back toward zero.
78    
# Line 52  C     pulls div(U) + dh/dt back toward z Line 80  C     pulls div(U) + dh/dt back toward z
80  C     Fully Implicit treatment of the Barotropic Flow Divergence  C     Fully Implicit treatment of the Barotropic Flow Divergence
81          DO j=1,sNy          DO j=1,sNy
82           DO i=1,sNx+1           DO i=1,sNx+1
83            pf(i,j) = xA(i,j)*gUNm1(i,j,k,bi,bj) / deltaTmom            pf(i,j) = xA(i,j)*gU(i,j,k,bi,bj) / deltaTMom
84           ENDDO           ENDDO
85          ENDDO          ENDDO
86        ELSEIF (exactConserv) THEN        ELSEIF (exactConserv) THEN
# Line 60  c     ELSEIF (nonlinFreeSurf.GT.0) THEN Line 88  c     ELSEIF (nonlinFreeSurf.GT.0) THEN
88  C     Implicit treatment of the Barotropic Flow Divergence  C     Implicit treatment of the Barotropic Flow Divergence
89          DO j=1,sNy          DO j=1,sNy
90           DO i=1,sNx+1           DO i=1,sNx+1
91            pf(i,j) = implicDiv2Dflow            pf(i,j) = implicDiv2Dflow
92       &             *xA(i,j)*gUNm1(i,j,k,bi,bj) / deltaTmom       &             *xA(i,j)*gU(i,j,k,bi,bj) / deltaTMom
93           ENDDO           ENDDO
94          ENDDO          ENDDO
95        ELSE        ELSE
# Line 69  C     Explicit+Implicit part of the Baro Line 97  C     Explicit+Implicit part of the Baro
97  C      => Filtering of uVel,vVel is necessary  C      => Filtering of uVel,vVel is necessary
98  C-- Now the filter are applied in the_correction_step().  C-- Now the filter are applied in the_correction_step().
99  C   We have left this code here to indicate where the filters used to be  C   We have left this code here to indicate where the filters used to be
100  C   in the algorithm before JMC moved them to after the pressure solver.  C   in the algorithm before JMC moved them to after the pressure solver.
101  C-  c#ifdef ALLOW_ZONAL_FILT
102  C#ifdef ALLOW_ZONAL_FILT  c        IF (zonal_filt_lat.LT.90.) THEN
103  C        IF (zonal_filt_lat.LT.90.) THEN  c          CALL ZONAL_FILTER(
104  C          CALL ZONAL_FILTER(  c    U                     uVel( 1-OLx,1-OLy,k,bi,bj),
105  C     &      uVel, hFacW, 1-1, sNy+1, k, k, bi, bj, 1, myThid)  c    I                     hFacW(1-OLx,1-OLy,k,bi,bj),
106  C          CALL ZONAL_FILTER(  c    I                     0, sNy+1, 1, bi, bj, 1, myThid )
107  C     &      vVel, hFacS, 1-1, sNy+1, k, k, bi, bj, 2, myThid)  c          CALL ZONAL_FILTER(
108  C        ENDIF  c    U                     vVel( 1-OLx,1-OLy,k,bi,bj),
109  C#endif  c    I                     hFacS(1-OLx,1-OLy,k,bi,bj),
110    c    I                     0, sNy+1, 1, bi, bj, 2, myThid )
111    c        ENDIF
112    c#endif
113          DO j=1,sNy          DO j=1,sNy
114           DO i=1,sNx+1           DO i=1,sNx+1
115            pf(i,j) = ( implicDiv2Dflow * gUNm1(i,j,k,bi,bj)            pf(i,j) = ( implicDiv2Dflow * gU(i,j,k,bi,bj)
116       &          + (1.-implicDiv2Dflow) * uVel(i,j,k,bi,bj)       &     + (1. _d 0-implicDiv2Dflow)* uVel(i,j,k,bi,bj)
117       &               ) * xA(i,j) / deltaTmom       &               ) * xA(i,j) / deltaTMom
118           ENDDO           ENDDO
119          ENDDO          ENDDO
120        ENDIF        ENDIF
# Line 95  C#endif Line 126  C#endif
126        ENDDO        ENDDO
127    
128  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
129        IF (nonHydrostatic) THEN        IF (use3Dsolver) THEN
130         DO j=1,sNy         DO j=1,sNy
131          DO i=1,sNx          DO i=1,sNx
132           cg3d_b(i,j,k,bi,bj) = pf(i+1,j)-pf(i,j)           cg3d_b(i,j,k,bi,bj) = ( pf(i+1,j)-pf(i,j) )
133          ENDDO          ENDDO
134         ENDDO         ENDDO
135        ENDIF        ENDIF
# Line 108  C#endif Line 139  C#endif
139  C     Fully Implicit treatment of the Barotropic Flow Divergence  C     Fully Implicit treatment of the Barotropic Flow Divergence
140          DO j=1,sNy+1          DO j=1,sNy+1
141           DO i=1,sNx           DO i=1,sNx
142            pf(i,j) = yA(i,j)*gVNm1(i,j,k,bi,bj) / deltatmom            pf(i,j) = yA(i,j)*gV(i,j,k,bi,bj) / deltatmom
143           ENDDO           ENDDO
144          ENDDO          ENDDO
145        ELSEIF (exactConserv) THEN        ELSEIF (exactConserv) THEN
# Line 117  C     Implicit treatment of the Barotrop Line 148  C     Implicit treatment of the Barotrop
148          DO j=1,sNy+1          DO j=1,sNy+1
149           DO i=1,sNx           DO i=1,sNx
150            pf(i,j) = implicDiv2Dflow            pf(i,j) = implicDiv2Dflow
151       &             *yA(i,j)*gVNm1(i,j,k,bi,bj) / deltatmom       &             *yA(i,j)*gV(i,j,k,bi,bj) / deltatmom
152           ENDDO           ENDDO
153          ENDDO          ENDDO
154        ELSE        ELSE
155  C     Explicit+Implicit part of the Barotropic Flow Divergence  C     Explicit+Implicit part of the Barotropic Flow Divergence
156          DO j=1,sNy+1          DO j=1,sNy+1
157           DO i=1,sNx           DO i=1,sNx
158            pf(i,j) = ( implicDiv2Dflow * gVNm1(i,j,k,bi,bj)            pf(i,j) = ( implicDiv2Dflow * gV(i,j,k,bi,bj)
159       &          + (1.-implicDiv2Dflow) * vVel(i,j,k,bi,bj)       &     + (1. _d 0-implicDiv2Dflow)* vVel(i,j,k,bi,bj)
160       &               ) * yA(i,j) / deltaTmom       &               ) * yA(i,j) / deltaTMom
161           ENDDO           ENDDO
162          ENDDO          ENDDO
163        ENDIF        ENDIF
# Line 138  C     Explicit+Implicit part of the Baro Line 169  C     Explicit+Implicit part of the Baro
169        ENDDO        ENDDO
170    
171  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
172        IF (nonHydrostatic) THEN        IF (use3Dsolver) THEN
173         DO j=1,sNy         DO j=1,sNy
174          DO i=1,sNx          DO i=1,sNx
175           cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj) +           cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
176       &    pf(i,j+1)-pf(i,j)       &                       + ( pf(i,j+1)-pf(i,j) )
177          ENDDO          ENDDO
178         ENDDO         ENDDO
179        ENDIF        ENDIF
180  #endif  #endif
181    
182    #ifdef ALLOW_ADDFLUID
183          IF ( selectAddFluid.GE.1 ) THEN
184            DO j=1,sNy
185             DO i=1,sNx
186              cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
187         &         - addMass(i,j,k,bi,bj)*mass2rUnit/deltaTMom
188             ENDDO
189            ENDDO
190    #ifdef ALLOW_NONHYDROSTATIC
191           IF (use3Dsolver) THEN
192            DO j=1,sNy
193             DO i=1,sNx
194              cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
195         &         - addMass(i,j,k,bi,bj)*mass2rUnit/deltaTMom
196             ENDDO
197            ENDDO
198           ENDIF
199    #endif
200          ENDIF
201    #endif /* ALLOW_ADDFLUID */
202    
203        RETURN        RETURN
204        END        END
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.22