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