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

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

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

revision 1.17 by jmc, Tue Mar 6 16:57:10 2001 UTC revision 1.22 by adcroft, Tue May 29 14:01:37 2001 UTC
# Line 26  C     == Global variables Line 26  C     == Global variables
26  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
27  #include "OBCS.h"  #include "OBCS.h"
28  #endif  #endif
29    #include "SOLVE_FOR_PRESSURE.h"
30    
31  C     == Routine arguments ==  C     == Routine arguments ==
32  C     myThid - Number of this instance of SOLVE_FOR_PRESSURE  C     myThid - Number of this instance of SOLVE_FOR_PRESSURE
33        INTEGER myThid        INTEGER myThid
34  CEndOfInterface  CEndOfInterface
35    
36  C     Local variables  C     == Local variables ==
 C     cg2d_x - Conjugate Gradient 2-D solver : Solution vector  
 C     cg2d_b - Conjugate Gradient 2-D solver : Right-hand side vector  
37        INTEGER i,j,k,bi,bj        INTEGER i,j,k,bi,bj
38        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
39        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
40        _RL cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL firstResidual,lastResidual
41        _RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        INTEGER numIters
42    
43  C--   Save previous solution & Initialise Vector solution and source term :  C--   Save previous solution & Initialise Vector solution and source term :
44        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
45         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
46          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
47           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
 #ifdef INCLUDE_CD_CODE  
48            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)            etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
49  #endif            cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
           cg2d_x(i,j,bi,bj) = etaN(i,j,bi,bj)  
50            cg2d_b(i,j,bi,bj) = 0.            cg2d_b(i,j,bi,bj) = 0.
51  #ifdef USE_NATURAL_BCS  #ifdef USE_NATURAL_BCS
52       &     + freeSurfFac*_rA(i,j,bi,bj)*horiVertRatio*       &     + freeSurfFac*_rA(i,j,bi,bj)*
53       &       EmPmR(I,J,bi,bj)/deltaTMom       &       EmPmR(I,J,bi,bj)/deltaTMom
54  #endif  #endif
55           ENDDO           ENDDO
# Line 87  C--   Add source term arising from w=d/d Line 84  C--   Add source term arising from w=d/d
84          DO j=1,sNy          DO j=1,sNy
85           DO i=1,sNx           DO i=1,sNx
86            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
87       &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(       &      -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
88       &          -cg2d_x(I,J,bi,bj)       &        *( etaN(i,j,bi,bj)
89       &          -cg3d_x(I,J,1,bi,bj)       &          +cg3d_x(i,j,1,bi,bj)*horiVertRatio/gravity )
      &        )/deltaTMom/deltaTMom  
90            cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)            cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
91       &      +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(       &      -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
92       &         -cg2d_x(I,J,bi,bj)       &        *( etaN(i,j,bi,bj)
93       &         -cg3d_x(I,J,1,bi,bj)       &          +cg3d_x(i,j,1,bi,bj)*horiVertRatio/gravity )
94       &       )/deltaTMom/deltaTMom  C-jmc
95    c    &      -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)
96    c    &        *( cg2d_x(i,j,bi,bj) + cg3d_x(i,j,1,bi,bj) )
97    c    &        /deltaTMom/deltaTMom
98    C-jmc
99           ENDDO           ENDDO
100          ENDDO          ENDDO
101  #else  #else
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       &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(       &      -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTMom
106       &          -cg2d_x(I,J,bi,bj)       &        * etaN(i,j,bi,bj)
      &        )/deltaTMom/deltaTMom  
107           ENDDO           ENDDO
108          ENDDO          ENDDO
109  #endif  #endif
# Line 139  C Western boundary Line 138  C Western boundary
138    
139  C--   Find the surface pressure using a two-dimensional conjugate  C--   Find the surface pressure using a two-dimensional conjugate
140  C--   gradient solver.  C--   gradient solver.
141  C     see CG2D_INTERNAL.h for the interface to this routine.  C     see CG2D.h for the interface to this routine.
142          firstResidual=0.
143          lastResidual=0.
144          numIters=cg2dMaxIters
145        CALL CG2D(        CALL CG2D(
146       I           cg2d_b,       U           cg2d_b,
147       U           cg2d_x,       U           cg2d_x,
148         O           firstResidual,
149         O           lastResidual,
150         U           numIters,
151       I           myThid )       I           myThid )
   
152        _EXCH_XY_R8(cg2d_x, myThid )        _EXCH_XY_R8(cg2d_x, myThid )
153    
154          _BEGIN_MASTER( myThid )
155           WRITE(*,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
156         &                               0, firstResidual
157           WRITE(*,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
158         &                         numIters, lastResidual
159          _END_MASTER( )
160    
161  C--   Transfert the 2D-solution to "etaN" :  C--   Transfert the 2D-solution to "etaN" :
162        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
163         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
164          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
165           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
166            etaN(i,j,bi,bj) = cg2d_x(i,j,bi,bj)            etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
167           ENDDO           ENDDO
168          ENDDO          ENDDO
169         ENDDO         ENDDO
# Line 167  C     see CG3D.h for the interface to th Line 178  C     see CG3D.h for the interface to th
178          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
179           DO j=1,sNy+1           DO j=1,sNy+1
180            DO i=1,sNx+1            DO i=1,sNx+1
181             uf(i,j)=-gBaro*_recip_dxC(i,j,bi,bj)*             uf(i,j)=-_recip_dxC(i,j,bi,bj)*
182       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
183             vf(i,j)=-gBaro*_recip_dyC(i,j,bi,bj)*             vf(i,j)=-_recip_dyC(i,j,bi,bj)*
184       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
185            ENDDO            ENDDO
186           ENDDO           ENDDO
# Line 207  C Western boundary Line 218  C Western boundary
218       &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)       &       -dRF(K)*dYG( i ,j,bi,bj)*hFacW( i ,j,k,bi,bj)*uf( i ,j)
219       &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)       &       +dRF(K)*dXG(i,j+1,bi,bj)*hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
220       &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )       &       -dRF(K)*dXG(i, j ,bi,bj)*hFacS(i, j ,k,bi,bj)*vf(i, j )
221       &       +(       &       +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
222       &         -wVel(i,j,k+1,bi,bj)       &          -wVel(i,j,k+1,bi,bj)
223       &        )*_rA(i,j,bi,bj)/deltaTmom       &        )*_rA(i,j,bi,bj)/deltaTmom
      &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(  
      &          +cg2d_x(I,J,bi,bj)  
      &        )/deltaTMom/deltaTMom  
224            ENDDO            ENDDO
225           ENDDO           ENDDO
226           DO K=2,Nr-1           DO K=2,Nr-1
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22