/[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.70 by jmc, Wed Nov 25 20:56:15 2009 UTC revision 1.71 by jmc, Sun Nov 29 03:17:05 2009 UTC
# Line 7  C $Name$ Line 7  C $Name$
7  CBOP  CBOP
8  C     !ROUTINE: SOLVE_FOR_PRESSURE  C     !ROUTINE: SOLVE_FOR_PRESSURE
9  C     !INTERFACE:  C     !INTERFACE:
10        SUBROUTINE SOLVE_FOR_PRESSURE(myTime, myIter, myThid)        SUBROUTINE SOLVE_FOR_PRESSURE( myTime, myIter, myThid )
11    
12  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
13  C     *==========================================================*  C     *==========================================================*
# Line 59  C     == Local variables == Line 59  C     == Local variables ==
59        _RL tmpFac        _RL tmpFac
60        _RL sumEmP, tileEmP(nSx,nSy)        _RL sumEmP, tileEmP(nSx,nSy)
61        LOGICAL putPmEinXvector        LOGICAL putPmEinXvector
62        INTEGER numIters, ks        INTEGER numIters, ks, ioUnit
63        CHARACTER*10 sufx        CHARACTER*10 sufx
64        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
65  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
66        INTEGER kp1        INTEGER kp1
67        _RL     wFacKm, wFacKp        _RL     wFacKm, wFacKp
68        LOGICAL zeroPsNH        LOGICAL zeroPsNH, zeroMeanPnh, oldFreeSurfTerm
69          _RL     tmpVar(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
70        _RL     uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL     uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
71        _RL     vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL     vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
72  #else  #else
# Line 75  CEOP Line 76  CEOP
76    
77  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
78          zeroPsNH = .FALSE.          zeroPsNH = .FALSE.
79  c       zeroPsNH = exactConserv  c       zeroPsNH = use3Dsolver .AND. exactConserv
80    c    &                         .AND. select_rStar.EQ.0
81            zeroMeanPnh = .FALSE.
82    c       zeroMeanPnh = use3Dsolver .AND. select_rStar.NE.0
83            oldFreeSurfTerm = use3Dsolver .AND. select_rStar.EQ.0
84         &                                .AND. .NOT.zeroPsNH
85    c       oldFreeSurfTerm = use3Dsolver .AND. .NOT.exactConserv
86  #else  #else
87          cg3d_b(1) = 0.          cg3d_b(1) = 0.
88  #endif  #endif
# Line 95  C     the case where |Global_mean_PmE| i Line 102  C     the case where |Global_mean_PmE| i
102        putPmEinXvector = .FALSE.        putPmEinXvector = .FALSE.
103  c     putPmEinXvector = useRealFreshWaterFlux.AND.fluidIsWater  c     putPmEinXvector = useRealFreshWaterFlux.AND.fluidIsWater
104    
105          IF ( myIter.EQ.1+nIter0 .AND. debugLevel .GE. debLevA ) THEN
106            _BEGIN_MASTER( myThid )
107            ioUnit = standardMessageUnit
108            WRITE(msgBuf,'(2A,L5)') 'SOLVE_FOR_PRESSURE:',
109         &       ' putPmEinXvector =', putPmEinXvector
110            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
111    #ifdef ALLOW_NONHYDROSTATIC
112            WRITE(msgBuf,'(A,2(A,L5))') 'SOLVE_FOR_PRESSURE:',
113         &       ' zeroPsNH=', zeroPsNH, ' , zeroMeanPnh=', zeroMeanPnh
114            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
115            WRITE(msgBuf,'(2A,L5)') 'SOLVE_FOR_PRESSURE:',
116         &       ' oldFreeSurfTerm =', oldFreeSurfTerm
117            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
118    #endif
119            _END_MASTER( myThid )
120          ENDIF
121    
122  C--   Save previous solution & Initialise Vector solution and source term :  C--   Save previous solution & Initialise Vector solution and source term :
123        sumEmP = 0.        sumEmP = 0.
124        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
# Line 179  C      (has been done for cg2d_b ; and a Line 203  C      (has been done for cg2d_b ; and a
203             ENDDO             ENDDO
204            ENDDO            ENDDO
205          ENDIF          ENDIF
206          IF ( use3Dsolver .AND. zeroPsNH ) THEN  
207           DO j=1,sNy          IF ( oldFreeSurfTerm ) THEN
           DO i=1,sNx  
            ks = ksurfC(i,j,bi,bj)  
            IF ( ks.LE.Nr ) THEN  
             cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)  
      &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)  
      &         /deltaTMom/deltaTfreesurf  
      &         * etaH(i,j,bi,bj)  
             cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)  
      &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)  
      &         /deltaTMom/deltaTfreesurf  
      &         * etaH(i,j,bi,bj)  
            ENDIF  
           ENDDO  
          ENDDO  
         ELSEIF ( use3Dsolver ) THEN  
208           DO j=1,sNy           DO j=1,sNy
209            DO i=1,sNx            DO i=1,sNx
210             ks = ksurfC(i,j,bi,bj)             ks = ksurfC(i,j,bi,bj)
# Line 369  C--   Solve for a three-dimensional pres Line 378  C--   Solve for a three-dimensional pres
378  C     see CG3D.h for the interface to this routine.  C     see CG3D.h for the interface to this routine.
379         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
380          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
381    
382    C--   Update or Add free-surface contribution to cg3d_b:
383    c        IF ( select_rStar.EQ.0 .AND. exactConserv ) THEN
384             IF ( select_rStar.EQ.0 .AND. .NOT.oldFreeSurfTerm ) THEN
385               tmpFac = 0.
386               DO j=1,sNy
387                DO i=1,sNx
388                  ks = ksurfC(i,j,bi,bj)
389                  IF ( ks.LE.Nr ) THEN
390                   cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
391         &                  +freeSurfFac*(etaN(i,j,bi,bj)-etaH(i,j,bi,bj))
392         &                              *_rA(i,j,bi,bj)*deepFac2F(ks)
393         &                              /deltaTMom/deltaTfreesurf
394                  ENDIF
395                ENDDO
396               ENDDO
397    #ifdef NONLIN_FRSURF
398             ELSEIF ( select_rStar.NE.0 ) THEN
399               tmpFac = 0.
400               DO j=1,sNy
401                DO i=1,sNx
402                  ks = ksurfC(i,j,bi,bj)
403                  tmpVar(i,j) = freeSurfFac
404         &                    *( etaN(i,j,bi,bj) - etaH(i,j,bi,bj) )
405         &                    *_rA(i,j,bi,bj)*deepFac2F(ks)
406         &                    /deltaTMom/deltaTfreesurf
407         &                    *recip_Rcol(i,j,bi,bj)
408                ENDDO
409               ENDDO
410               DO k=1,Nr
411                DO j=1,sNy
412                 DO i=1,sNx
413                  cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
414         &          + tmpVar(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
415                 ENDDO
416                ENDDO
417               ENDDO
418    #endif /* NONLIN_FRSURF */
419             ELSEIF ( usingZCoords ) THEN
420    C-       Z coordinate: assume surface @ level k=1
421               tmpFac = freeSurfFac*deepFac2F(1)
422             ELSE
423    C-       Other than Z coordinate: no assumption on surface level index
424               tmpFac = 0.
425               DO j=1,sNy
426                DO i=1,sNx
427                  ks = ksurfC(i,j,bi,bj)
428                  IF ( ks.LE.Nr ) THEN
429                   cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
430         &              +freeSurfFac*etaN(i,j,bi,bj)/deltaTfreesurf
431         &                  *_rA(i,j,bi,bj)*deepFac2F(ks)/deltaTmom
432                  ENDIF
433                ENDDO
434               ENDDO
435             ENDIF
436    
437    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
438    
439    C--   Finish updating cg3d_b: 1) increment in horiz velocity due to new cg2d_x
440    C                             2) add vertical velocity contribution.
441           DO j=1,sNy+1           DO j=1,sNy+1
442            DO i=1,sNx+1            DO i=1,sNx+1
443             uf(i,j)=-_recip_dxC(i,j,bi,bj)*             uf(i,j) = -_recip_dxC(i,j,bi,bj)
444       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))       &             * implicSurfPress*implicDiv2DFlow
445             vf(i,j)=-_recip_dyC(i,j,bi,bj)*       &             *(cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
446       &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))             vf(i,j) = -_recip_dyC(i,j,bi,bj)
447         &             * implicSurfPress*implicDiv2DFlow
448         &             *(cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
449            ENDDO            ENDDO
450           ENDDO           ENDDO
451    
# Line 382  C     see CG3D.h for the interface to th Line 453  C     see CG3D.h for the interface to th
453           IF (useOBCS) THEN           IF (useOBCS) THEN
454            DO i=1,sNx+1            DO i=1,sNx+1
455  C Northern boundary  C Northern boundary
456            IF (OB_Jn(i,bi,bj).NE.0) THEN             IF (OB_Jn(i,bi,bj).NE.0)
457             vf(i,OB_Jn(i,bi,bj))=0.       &      vf(i,OB_Jn(i,bi,bj)) = 0.
           ENDIF  
458  C Southern boundary  C Southern boundary
459            IF (OB_Js(i,bi,bj).NE.0) THEN             IF (OB_Js(i,bi,bj).NE.0)
460             vf(i,OB_Js(i,bi,bj)+1)=0.       &      vf(i,OB_Js(i,bi,bj)+1) = 0.
           ENDIF  
461            ENDDO            ENDDO
462            DO j=1,sNy+1            DO j=1,sNy+1
463  C Eastern boundary  C Eastern boundary
464            IF (OB_Ie(j,bi,bj).NE.0) THEN             IF (OB_Ie(j,bi,bj).NE.0)
465             uf(OB_Ie(j,bi,bj),j)=0.       &      uf(OB_Ie(j,bi,bj),j) = 0.
           ENDIF  
466  C Western boundary  C Western boundary
467            IF (OB_Iw(j,bi,bj).NE.0) THEN             IF (OB_Iw(j,bi,bj).NE.0)
468             uf(OB_Iw(j,bi,bj)+1,J)=0.       &      uf(OB_Iw(j,bi,bj)+1,j) = 0.
           ENDIF  
469            ENDDO            ENDDO
470           ENDIF           ENDIF
471  #endif /* ALLOW_OBCS */  #endif /* ALLOW_OBCS */
472    
473           IF ( usingZCoords ) THEN  C Note: with implicDiv2DFlow < 1, wVel contribution to cg3d_b is similar to
474  C-       Z coordinate: assume surface @ level k=1  C       uVel,vVel contribution to cg2d_b when exactConserv=T, since wVel is
475             tmpFac = freeSurfFac*deepFac2F(1)  C       always recomputed from continuity eq (like eta when exactConserv=T)
          ELSE  
 C-       Other than Z coordinate: no assumption on surface level index  
            tmpFac = 0.  
            DO j=1,sNy  
             DO i=1,sNx  
               ks = ksurfC(i,j,bi,bj)  
               IF ( ks.LE.Nr ) THEN  
                cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)  
      &              +freeSurfFac*etaN(i,j,bi,bj)/deltaTfreesurf  
      &                  *_rA(i,j,bi,bj)*deepFac2F(ks)/deltaTmom  
               ENDIF  
             ENDDO  
            ENDDO  
          ENDIF  
476           k=1           k=1
477           kp1 = MIN(k+1,Nr)           kp1 = MIN(k+1,Nr)
478           wFacKp = deepFac2F(kp1)*rhoFacF(kp1)           wFacKp = implicDiv2DFlow*deepFac2F(kp1)*rhoFacF(kp1)
479           IF (k.GE.Nr) wFacKp = 0.           IF (k.GE.Nr) wFacKp = 0.
480           DO j=1,sNy           DO j=1,sNy
481            DO i=1,sNx            DO i=1,sNx
# Line 440  C-       Other than Z coordinate: no ass Line 493  C-       Other than Z coordinate: no ass
493            kp1 = MIN(k+1,Nr)            kp1 = MIN(k+1,Nr)
494  C-       deepFac & rhoFac cancel with the ones in uf[=del_i(Phi)/dx],vf ;  C-       deepFac & rhoFac cancel with the ones in uf[=del_i(Phi)/dx],vf ;
495  C        both appear in wVel term, but at 2 different levels  C        both appear in wVel term, but at 2 different levels
496            wFacKm = deepFac2F( k )*rhoFacF( k )            wFacKm = implicDiv2DFlow*deepFac2F( k )*rhoFacF( k )
497            wFacKp = deepFac2F(kp1)*rhoFacF(kp1)            wFacKp = implicDiv2DFlow*deepFac2F(kp1)*rhoFacF(kp1)
498            IF (k.GE.Nr) wFacKp = 0.            IF (k.GE.Nr) wFacKp = 0.
499            DO j=1,sNy            DO j=1,sNy
500             DO i=1,sNx             DO i=1,sNx
# Line 461  C        both appear in wVel term, but a Line 514  C        both appear in wVel term, but a
514  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
515           IF (useOBCS) THEN           IF (useOBCS) THEN
516            DO k=1,Nr            DO k=1,Nr
517            DO i=1,sNx             DO i=1,sNx
518  C Northern boundary  C Northern boundary
519            IF (OB_Jn(i,bi,bj).NE.0) THEN              IF (OB_Jn(i,bi,bj).NE.0)
520             cg3d_b(i,OB_Jn(i,bi,bj),k,bi,bj)=0.       &       cg3d_b(i,OB_Jn(i,bi,bj),k,bi,bj) = 0.
           ENDIF  
521  C Southern boundary  C Southern boundary
522            IF (OB_Js(i,bi,bj).NE.0) THEN              IF (OB_Js(i,bi,bj).NE.0)
523             cg3d_b(i,OB_Js(i,bi,bj),k,bi,bj)=0.       &       cg3d_b(i,OB_Js(i,bi,bj),k,bi,bj) = 0.
524            ENDIF             ENDDO
525            ENDDO             DO j=1,sNy
           DO j=1,sNy  
526  C Eastern boundary  C Eastern boundary
527            IF (OB_Ie(j,bi,bj).NE.0) THEN              IF (OB_Ie(j,bi,bj).NE.0)
528             cg3d_b(OB_Ie(j,bi,bj),j,k,bi,bj)=0.       &       cg3d_b(OB_Ie(j,bi,bj),j,k,bi,bj) = 0.
           ENDIF  
529  C Western boundary  C Western boundary
530            IF (OB_Iw(j,bi,bj).NE.0) THEN              IF (OB_Iw(j,bi,bj).NE.0)
531             cg3d_b(OB_Iw(j,bi,bj),j,k,bi,bj)=0.       &       cg3d_b(OB_Iw(j,bi,bj),j,k,bi,bj) = 0.
532            ENDIF             ENDDO
           ENDDO  
533            ENDDO            ENDDO
534           ENDIF           ENDIF
535  #endif /* ALLOW_OBCS */  #endif /* ALLOW_OBCS */
# Line 509  C-    end bi,bj loops Line 558  C-    end bi,bj loops
558       O           firstResidual,       O           firstResidual,
559       O           lastResidual,       O           lastResidual,
560       U           numIters,       U           numIters,
561       I           myThid )       I           myIter, myThid )
562        _EXCH_XYZ_RL( phi_nh, myThid )        _EXCH_XYZ_RL( phi_nh, myThid )
563        CALL TIMER_STOP ('CG3D   [SOLVE_FOR_PRESSURE]',myThid)        CALL TIMER_STOP ('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
564    
# Line 528  C-    end bi,bj loops Line 577  C-    end bi,bj loops
577        ENDIF        ENDIF
578    
579  C--   Update surface pressure (account for NH-p @ surface level) and NH pressure:  C--   Update surface pressure (account for NH-p @ surface level) and NH pressure:
580        IF ( zeroPsNH ) THEN        IF ( zeroPsNH .OR. zeroMeanPnh ) THEN
581         IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN         IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
582          WRITE(sufx,'(I10.10)') myIter          WRITE(sufx,'(I10.10)') myIter
583          CALL WRITE_FLD_XYZ_RL( 'cg3d_x.',sufx,phi_nh, myIter, myThid )          CALL WRITE_FLD_XYZ_RL( 'cg3d_x.',sufx, phi_nh, myIter, myThid )
584         ENDIF         ENDIF
585         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
586          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
587    
588           IF ( usingZCoords ) THEN           IF ( zeroPsNH .AND. usingZCoords ) THEN
589  C-       Z coordinate: assume surface @ level k=1  C-       Z coordinate: assume surface @ level k=1
           DO k=2,Nr  
            DO j=1-OLy,sNy+OLy  
             DO i=1-OLx,sNx+OLx  
              phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)  
      &                           - phi_nh(i,j,1,bi,bj)  
             ENDDO  
            ENDDO  
           ENDDO  
590            DO j=1-OLy,sNy+OLy            DO j=1-OLy,sNy+OLy
591             DO i=1-OLx,sNx+OLx             DO i=1-OLx,sNx+OLx
592               etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)               tmpVar(i,j) = phi_nh(i,j,1,bi,bj)
      &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,1,bi,bj))  
              phi_nh(i,j,1,bi,bj) = 0.  
593             ENDDO             ENDDO
594            ENDDO            ENDDO
595           ELSE           ELSEIF ( zeroPsNH ) THEN
596  C-       Other than Z coordinate: no assumption on surface level index  C-       Other than Z coordinate: no assumption on surface level index
597            DO j=1-OLy,sNy+OLy            DO j=1-OLy,sNy+OLy
598             DO i=1-OLx,sNx+OLx             DO i=1-OLx,sNx+OLx
599              ks = ksurfC(i,j,bi,bj)              ks = ksurfC(i,j,bi,bj)
600              IF ( ks.LE.Nr ) THEN              IF ( ks.LE.Nr ) THEN
601               etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)               tmpVar(i,j) = phi_nh(i,j,ks,bi,bj)
602       &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,ks,bi,bj))              ELSE
603               DO k=Nr,1,-1               tmpVar(i,j) = 0.
               phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)  
      &                            - phi_nh(i,j,ks,bi,bj)  
              ENDDO  
604              ENDIF              ENDIF
605             ENDDO             ENDDO
606            ENDDO            ENDDO
607    #ifdef NONLIN_FRSURF
608             ELSE
609    C        zeroMeanPnh : transfert vertical average of P_NH to EtaN
610              DO j=1-OLy,sNy+OLy
611               DO i=1-OLx,sNx+OLx
612                 tmpVar(i,j) = 0.
613               ENDDO
614              ENDDO
615              DO k=1,Nr
616               DO j=1-OLy,sNy+OLy
617                DO i=1-OLx,sNx+OLx
618                 tmpVar(i,j) = tmpVar(i,j)
619         &         + phi_nh(i,j,k,bi,bj)*drF(k)*h0FacC(i,j,k,bi,bj)
620                ENDDO
621               ENDDO
622              ENDDO
623              DO j=1-OLy,sNy+OLy
624               DO i=1-OLx,sNx+OLx
625                 tmpVar(i,j) = tmpVar(i,j)*recip_Rcol(i,j,bi,bj)
626               ENDDO
627              ENDDO
628    #endif /* NONLIN_FRSURF */
629           ENDIF           ENDIF
630             DO k=1,Nr
631              DO j=1-OLy,sNy+OLy
632               DO i=1-OLx,sNx+OLx
633                phi_nh(i,j,k,bi,bj) = ( phi_nh(i,j,k,bi,bj)
634         &                            - tmpVar(i,j)
635         &                            )*maskC(i,j,k,bi,bj)
636               ENDDO
637              ENDDO
638             ENDDO
639             DO j=1-OLy,sNy+OLy
640              DO i=1-OLx,sNx+OLx
641                etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
642         &                      *( cg2d_x(i,j,bi,bj) + tmpVar(i,j) )
643              ENDDO
644             ENDDO
645    
646          ENDDO          ENDDO
647         ENDDO         ENDDO
Legend:
Removed from v.1.70  
changed lines
  Added in v.1.71
  ViewVC Help
Powered by ViewVC 1.1.22