/[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.18 by jmc, Thu Apr 17 13:40:06 2003 UTC revision 1.23 by jmc, Sun Aug 24 21:38:18 2008 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    
6  CBOP  CBOP
7  C     !ROUTINE: CALC_DIV_GHAT  C     !ROUTINE: CALC_DIV_GHAT
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE CALC_DIV_GHAT(        SUBROUTINE CALC_DIV_GHAT(
10       I        bi,bj,iMin,iMax,jMin,jMax,K,       I                bi,bj,k,
11       I        xA,yA,       U                cg2d_b, cg3d_b,
12       U        cg2d_b,       I                myThid)
      I        myThid)  
13  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
14  C     *==========================================================*  C     *==========================================================*
15  C     | S/R CALC_DIV_GHAT                                          C     | S/R CALC_DIV_GHAT
16  C     | o Form the right hand-side of the surface pressure eqn.    C     | o Form the right hand-side of the surface pressure eqn.
17  C     *==========================================================*  C     *==========================================================*
18  C     | Right hand side of pressure equation is divergence  C     | Right hand side of pressure equation is divergence
19  C     | of veclocity tendency (GHAT) term along with a relaxation  C     | of veclocity tendency (GHAT) term along with a relaxation
# Line 25  C     !USES: Line 25  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    
33  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
34  C     == Routine arguments ==  C     == Routine arguments ==
35  C     bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation  C     bi, bj  :: tile indices
36  C                                      results will be set.  C     k       :: Index of layer.
37  C     k                              - Index of layer.  C     cg2d_b  :: Conjugate Gradient 2-D solver : Right-hand side vector
38  C     xA, yA                         - Cell face areas  C     cg3d_b  :: Conjugate Gradient 3-D solver : Right-hand side vector
39  C     cg2d_b - Conjugate Gradient 2-D solver : Right-hand side vector  C     myThid  :: Instance number for this call of CALC_DIV_GHAT
40  C     myThid - Instance number for this call of CALC_DIV_GHAT        INTEGER bi,bj
41        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)  
42        _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)
43    c#ifdef ALLOW_NONHYDROSTATIC
44          _RL cg3d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
45    c#else
46    c     _RL cg3d_b(1)
47    c#endif
48        INTEGER myThid        INTEGER myThid
49    
50  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
51  C     == Local variables ==  C     == Local variables ==
52  C     i,j :: Loop counters  C     i,j    :: Loop counters
53  C     pf  :: Intermediate array for building RHS source term.  C     xA, yA :: Cell vertical face areas
54    C     pf     :: Intermediate array for building RHS source term.
55        INTEGER i,j        INTEGER i,j
56          _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57          _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58        _RL pf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL pf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59  CEOP  CEOP
60    
61    C     Calculate vertical face areas
62          DO j=1,sNy+1
63            DO i=1,sNx+1
64              xA(i,j) = _dyG(i,j,bi,bj)*deepFacC(k)
65         &              *drF(k)*_hFacW(i,j,k,bi,bj)*rhoFacC(k)
66              yA(i,j) = _dxG(i,j,bi,bj)*deepFacC(k)
67         &              *drF(k)*_hFacS(i,j,k,bi,bj)*rhoFacC(k)
68            ENDDO
69          ENDDO
70    
71  C--   Pressure equation source term  C--   Pressure equation source term
72  C     Term is the vertical integral of the divergence of the  C     Term is the vertical integral of the divergence of the
73  C     time tendency terms along with a relaxation term that  C     time tendency terms along with a relaxation term that
74  C     pulls div(U) + dh/dt back toward zero.  C     pulls div(U) + dh/dt back toward zero.
75    
# Line 72  c     ELSEIF (nonlinFreeSurf.GT.0) THEN Line 85  c     ELSEIF (nonlinFreeSurf.GT.0) THEN
85  C     Implicit treatment of the Barotropic Flow Divergence  C     Implicit treatment of the Barotropic Flow Divergence
86          DO j=1,sNy          DO j=1,sNy
87           DO i=1,sNx+1           DO i=1,sNx+1
88            pf(i,j) = implicDiv2Dflow            pf(i,j) = implicDiv2Dflow
89       &             *xA(i,j)*gU(i,j,k,bi,bj) / deltaTmom       &             *xA(i,j)*gU(i,j,k,bi,bj) / deltaTmom
90           ENDDO           ENDDO
91          ENDDO          ENDDO
# Line 81  C     Explicit+Implicit part of the Baro Line 94  C     Explicit+Implicit part of the Baro
94  C      => Filtering of uVel,vVel is necessary  C      => Filtering of uVel,vVel is necessary
95  C-- Now the filter are applied in the_correction_step().  C-- Now the filter are applied in the_correction_step().
96  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
97  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.
98  C-  C-
99  C#ifdef ALLOW_ZONAL_FILT  C#ifdef ALLOW_ZONAL_FILT
100  C        IF (zonal_filt_lat.LT.90.) THEN  C        IF (zonal_filt_lat.LT.90.) THEN
101  C          CALL ZONAL_FILTER(  C          CALL ZONAL_FILTER(
102  C     &      uVel, hFacW, 1-1, sNy+1, k, k, bi, bj, 1, myThid)  C     &      uVel, hFacW, 1-1, sNy+1, k, k, bi, bj, 1, myThid)
103  C          CALL ZONAL_FILTER(  C          CALL ZONAL_FILTER(
# Line 94  C#endif Line 107  C#endif
107          DO j=1,sNy          DO j=1,sNy
108           DO i=1,sNx+1           DO i=1,sNx+1
109            pf(i,j) = ( implicDiv2Dflow * gU(i,j,k,bi,bj)            pf(i,j) = ( implicDiv2Dflow * gU(i,j,k,bi,bj)
110       &          + (1.-implicDiv2Dflow) * uVel(i,j,k,bi,bj)       &     + (1. _d 0-implicDiv2Dflow)* uVel(i,j,k,bi,bj)
111       &               ) * xA(i,j) / deltaTmom       &               ) * xA(i,j) / deltaTmom
112           ENDDO           ENDDO
113          ENDDO          ENDDO
# Line 107  C#endif Line 120  C#endif
120        ENDDO        ENDDO
121    
122  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
123        IF (nonHydrostatic) THEN        IF (use3Dsolver) THEN
124         DO j=1,sNy         DO j=1,sNy
125          DO i=1,sNx          DO i=1,sNx
126           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)
# Line 137  C     Explicit+Implicit part of the Baro Line 150  C     Explicit+Implicit part of the Baro
150          DO j=1,sNy+1          DO j=1,sNy+1
151           DO i=1,sNx           DO i=1,sNx
152            pf(i,j) = ( implicDiv2Dflow * gV(i,j,k,bi,bj)            pf(i,j) = ( implicDiv2Dflow * gV(i,j,k,bi,bj)
153       &          + (1.-implicDiv2Dflow) * vVel(i,j,k,bi,bj)       &     + (1. _d 0-implicDiv2Dflow)* vVel(i,j,k,bi,bj)
154       &               ) * yA(i,j) / deltaTmom       &               ) * yA(i,j) / deltaTmom
155           ENDDO           ENDDO
156          ENDDO          ENDDO
# Line 150  C     Explicit+Implicit part of the Baro Line 163  C     Explicit+Implicit part of the Baro
163        ENDDO        ENDDO
164    
165  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
166        IF (nonHydrostatic) THEN        IF (use3Dsolver) THEN
167         DO j=1,sNy         DO j=1,sNy
168          DO i=1,sNx          DO i=1,sNx
169           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) +
# Line 160  C     Explicit+Implicit part of the Baro Line 173  C     Explicit+Implicit part of the Baro
173        ENDIF        ENDIF
174  #endif  #endif
175    
176    #ifdef ALLOW_ADDFLUID
177          IF ( selectAddFluid.GE.1 ) THEN
178            DO j=1,sNy
179             DO i=1,sNx
180              cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
181         &         - addMass(i,j,k,bi,bj)*mass2rUnit/deltaTmom
182             ENDDO
183            ENDDO
184    #ifdef ALLOW_NONHYDROSTATIC
185           IF (use3Dsolver) THEN
186            DO j=1,sNy
187             DO i=1,sNx
188              cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
189         &         - addMass(i,j,k,bi,bj)*mass2rUnit/deltaTmom
190             ENDDO
191            ENDDO
192           ENDIF
193    #endif
194          ENDIF
195    #endif /* ALLOW_ADDFLUID */
196    
197        RETURN        RETURN
198        END        END
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22