/[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.48 by jmc, Tue Nov 8 02:14:10 2005 UTC revision 1.49 by jmc, Mon Dec 12 15:50:51 2005 UTC
# Line 63  C     == Local variables == Line 63  C     == Local variables ==
63        LOGICAL putPmEinXvector        LOGICAL putPmEinXvector
64        INTEGER numIters        INTEGER numIters
65        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
66    #ifdef ALLOW_NONHYDROSTATIC
67          INTEGER ks
68          LOGICAL zeroPsNH
69    #endif
70  CEOP  CEOP
71    
72  #ifdef TIME_PER_TIMESTEP  #ifdef TIME_PER_TIMESTEP
# Line 73  C     == Timing variables == Line 77  C     == Timing variables ==
77        COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold        COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold
78  #endif  #endif
79    
80    #ifdef ALLOW_NONHYDROSTATIC
81    c       zeroPsNH = .FALSE.
82            zeroPsNH = exactConserv
83    #endif
84    
85  C--   Initialise the Vector solution with etaN + deltaT*Global_mean_PmE  C--   Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
86  C     instead of simply etaN ; This can speed-up the solver convergence in  C     instead of simply etaN ; This can speed-up the solver convergence in
87  C     the case where |Global_mean_PmE| is large.  C     the case where |Global_mean_PmE| is large.
# Line 154  C--   Add source term arising from w=d/d Line 163  C--   Add source term arising from w=d/d
163        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
164         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
165  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
166          IF ( nonHydrostatic ) THEN          IF ( nonHydrostatic .AND. zeroPsNH ) THEN
167             DO j=1,sNy
168              DO i=1,sNx
169               cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
170         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
171         &         * etaH(i,j,bi,bj)
172               cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
173         &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
174         &         * etaH(i,j,bi,bj)
175              ENDDO
176             ENDDO
177            ELSEIF ( nonHydrostatic ) THEN
178           DO j=1,sNy           DO j=1,sNy
179            DO i=1,sNx            DO i=1,sNx
180             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
# Line 215  C Western boundary Line 235  C Western boundary
235            ENDIF            ENDIF
236           ENDDO           ENDDO
237          ENDIF          ENDIF
238  #endif  #endif /* ALLOW_OBCS */
239    C-    end bi,bj loops
240         ENDDO         ENDDO
241        ENDDO        ENDDO
242    
# Line 313  C Western boundary Line 334  C Western boundary
334            ENDIF            ENDIF
335            ENDDO            ENDDO
336           ENDIF           ENDIF
337  #endif  #endif /* ALLOW_OBCS */
338    
339           K=1           K=1
340           DO j=1,sNy           DO j=1,sNy
341            DO i=1,sNx            DO i=1,sNx
342             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)
343       &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)       &       +drF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
344       &       -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)
345       &       +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)
346       &       -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 )
347       &       +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom       &       +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
348       &          -wVel(i,j,k+1,bi,bj)       &          -wVel(i,j,k+1,bi,bj)
349       &        )*_rA(i,j,bi,bj)/deltaTmom       &        )*_rA(i,j,bi,bj)/deltaTmom
# Line 332  C Western boundary Line 353  C Western boundary
353            DO j=1,sNy            DO j=1,sNy
354             DO i=1,sNx             DO i=1,sNx
355              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)
356       &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)       &       +drF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
357       &       -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)
358       &       +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)
359       &       -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 )
360       &       +( wVel(i,j,k  ,bi,bj)       &       +( wVel(i,j,k  ,bi,bj)
361       &         -wVel(i,j,k+1,bi,bj)       &         -wVel(i,j,k+1,bi,bj)
362       &        )*_rA(i,j,bi,bj)/deltaTmom       &        )*_rA(i,j,bi,bj)/deltaTmom
# Line 347  C Western boundary Line 368  C Western boundary
368           DO j=1,sNy           DO j=1,sNy
369            DO i=1,sNx            DO i=1,sNx
370              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)
371       &       +dRF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)       &       +drF(K)*dYG(i+1,j,bi,bj)*hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
372       &       -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)
373       &       +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)
374       &       -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 )
375       &       +( wVel(i,j,k  ,bi,bj)       &       +( wVel(i,j,k  ,bi,bj)
376       &        )*_rA(i,j,bi,bj)/deltaTmom       &        )*_rA(i,j,bi,bj)/deltaTmom
377    
# Line 382  C Western boundary Line 403  C Western boundary
403            ENDDO            ENDDO
404            ENDDO            ENDDO
405           ENDIF           ENDIF
406  #endif  #endif /* ALLOW_OBCS */
407    C-    end bi,bj loops
408          ENDDO ! bi          ENDDO
409         ENDDO ! bj         ENDDO
410    
411        firstResidual=0.        firstResidual=0.
412        lastResidual=0.        lastResidual=0.
413        numIters=cg2dMaxIters        numIters=cg3dMaxIters
414        CALL CG3D(        CALL CG3D(
415       U           cg3d_b,       U           cg3d_b,
416       U           phi_nh,       U           phi_nh,
# Line 413  C Western boundary Line 434  C Western boundary
434         ENDIF         ENDIF
435        ENDIF        ENDIF
436    
437    C--   Update surface pressure (account for NH-p @ surface level) and NH pressure:
438          IF ( zeroPsNH ) THEN
439           DO bj=myByLo(myThid),myByHi(myThid)
440            DO bi=myBxLo(myThid),myBxHi(myThid)
441    
442             IF ( usingZCoords ) THEN
443    C-       Z coordinate: assume surface @ level k=1
444              DO k=2,Nr
445               DO j=1-OLy,sNy+OLy
446                DO i=1-OLx,sNx+OLx
447                 phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
448         &                           - phi_nh(i,j,1,bi,bj)
449                ENDDO
450               ENDDO
451              ENDDO
452              DO j=1-OLy,sNy+OLy
453               DO i=1-OLx,sNx+OLx
454                 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
455         &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,1,bi,bj))
456                 phi_nh(i,j,1,bi,bj) = 0.
457               ENDDO
458              ENDDO
459             ELSE
460    C-       Other than Z coordinate: no assumption on surface level index
461              DO j=1-OLy,sNy+OLy
462               DO i=1-OLx,sNx+OLx
463                ks = ksurfC(i,j,bi,bj)
464                IF ( ks.LE.Nr ) THEN
465                 etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
466         &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,ks,bi,bj))
467                 DO k=Nr,1,-1
468                  phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
469         &                            - phi_nh(i,j,ks,bi,bj)
470                 ENDDO
471                ENDIF
472               ENDDO
473              ENDDO
474             ENDIF
475    
476            ENDDO
477           ENDDO
478        ENDIF        ENDIF
479  #endif  
480          ENDIF
481    #endif /* ALLOW_NONHYDROSTATIC */
482    
483  #ifdef TIME_PER_TIMESTEP  #ifdef TIME_PER_TIMESTEP
484  CCE107 Time per timestep information  CCE107 Time per timestep information
Legend:
Removed from v.1.48  
changed lines
  Added in v.1.49
  ViewVC Help
Powered by ViewVC 1.1.22