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

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

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

revision 1.42 by jmc, Tue Dec 18 01:16:53 2012 UTC revision 1.43 by jmc, Fri Jan 4 21:07:07 2013 UTC
# Line 65  C     myThid     :: thread number for th Line 65  C     myThid     :: thread number for th
65        INTEGER bi,bj,iMin,iMax,jMin,jMax,k        INTEGER bi,bj,iMin,iMax,jMin,jMax,k
66        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
67        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
 c     _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
68        _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69        _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70        _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 77  c     _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy Line 76  c     _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy
76    
77  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
78  C     == Local variables ==  C     == Local variables ==
79    C     phiHydU    :: hydrostatic potential anomaly at interface k+1 (Upper k)
80    C     pKappaF    :: (p/Po)^kappa at interface k
81    C     pKappaU    :: (p/Po)^kappa at interface k+1 (Upper k)
82    C     pKappaC    :: (p/Po)^kappa at cell center k
83        INTEGER i,j        INTEGER i,j
84        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85    #ifndef DISABLE_SIGMA_CODE
86          _RL phiHydU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87          _RL pKappaF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88          _RL pKappaU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89          _RL pKappaC, locDepth, fullDepth
90    #endif /* DISABLE_SIGMA_CODE */
91        _RL dRlocM,dRlocP, ddRloc, locAlpha        _RL dRlocM,dRlocP, ddRloc, locAlpha
92        _RL ddPIm, ddPIp, rec_dRm, rec_dRp        _RL ddPIm, ddPIp, rec_dRm, rec_dRp
93        _RL surfPhiFac        _RL surfPhiFac
94        LOGICAL useDiagPhiRlow, addSurfPhiAnom        LOGICAL useDiagPhiRlow, addSurfPhiAnom
95          LOGICAL useFVgradPhi
96  CEOP  CEOP
97        useDiagPhiRlow = .TRUE.        useDiagPhiRlow = .TRUE.
98        addSurfPhiAnom = select_rStar.EQ.0 .AND. nonlinFreeSurf.GE.4        addSurfPhiAnom = select_rStar.EQ.0 .AND. nonlinFreeSurf.GE.4
99          useFVgradPhi   = selectSigmaCoord.NE.0
100    
101        surfPhiFac = 0.        surfPhiFac = 0.
102        IF (addSurfPhiAnom) surfPhiFac = 1.        IF (addSurfPhiAnom) surfPhiFac = 1.
103    
# Line 430  C--     Quasi-hydrostatic terms are adde Line 442  C--     Quasi-hydrostatic terms are adde
442    
443  C---    Integrate d Phi / d pi  C---    Integrate d Phi / d pi
444    
445    #ifndef DISABLE_SIGMA_CODE
446    C  --  Finite Volume Form, integrated to r-level (cell center & upper interface)
447           IF ( useFVgradPhi ) THEN
448    
449             fullDepth = rF(1)-rF(Nr+1)
450             DO j=jMin,jMax
451              DO i=iMin,iMax
452               locDepth = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
453    #ifdef NONLIN_FRSURF
454               locDepth = locDepth + surfPhiFac*etaH(i,j,bi,bj)
455    #endif
456               pKappaF(i,j) = (
457         &           ( R_low(i,j,bi,bj) + aHybSigmF( k )*fullDepth
458         &                              + bHybSigmF( k )*locDepth
459         &           )/atm_Po )**atm_kappa
460               pKappaC      = (
461         &           ( R_low(i,j,bi,bj) + aHybSigmC( k )*fullDepth
462         &                              + bHybSigmC( k )*locDepth
463         &           )/atm_Po )**atm_kappa
464               pKappaU(i,j) = (
465         &           ( R_low(i,j,bi,bj)+ aHybSigmF(k+1)*fullDepth
466         &                             + bHybSigmF(k+1)*locDepth
467         &           )/atm_Po )**atm_kappa
468               phiHydC(i,j) = phiHydF(i,j)
469         &        + atm_Cp*( pKappaF(i,j) - pKappaC      )*alphaRho(i,j)
470               phiHydU(i,j) = phiHydF(i,j)
471         &        + atm_Cp*( pKappaF(i,j) - pKappaU(i,j) )*alphaRho(i,j)
472              ENDDO
473             ENDDO
474    C end: Finite Volume Form, integrated to r-level.
475    
476           ELSEIF (integr_GeoPot.EQ.0) THEN
477    #else /* DISABLE_SIGMA_CODE */
478         IF (integr_GeoPot.EQ.0) THEN         IF (integr_GeoPot.EQ.0) THEN
479    #endif /* DISABLE_SIGMA_CODE */
480  C  --  Energy Conserving Form, accurate with Full cell topo  --  C  --  Energy Conserving Form, accurate with Full cell topo  --
481  C------------ The integration for the first level phi(k=1) is the same  C------------ The integration for the first level phi(k=1) is the same
482  C             for both the "finite volume" and energy conserving methods.  C             for both the "finite volume" and energy conserving methods.
# Line 548  C---+----1----+----2----+----3----+----4 Line 594  C---+----1----+----2----+----3----+----4
594          STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'          STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
595        ENDIF        ENDIF
596    
597          IF ( .NOT. useFVgradPhi ) THEN
598    C--   r-coordinate and r*-coordinate cases:
599    
600           IF ( momPressureForcing ) THEN
601            CALL CALC_GRAD_PHI_HYD(
602         I                         k, bi, bj, iMin,iMax, jMin,jMax,
603         I                         phiHydC, alphaRho, tFld, sFld,
604         O                         dPhiHydX, dPhiHydY,
605         I                         myTime, myIter, myThid)
606           ENDIF
607    
608    #ifndef DISABLE_SIGMA_CODE
609          ELSE
610    C--   else (SigmaCoords part)
611    
612           IF ( fluidIsWater ) THEN
613            STOP 'CALC_PHI_HYD: missing code for SigmaCoord'
614           ENDIF
615           IF ( momPressureForcing ) THEN
616            CALL CALC_GRAD_PHI_FV(
617         I                         k, bi, bj, iMin,iMax, jMin,jMax,
618         I                         phiHydF, phiHydU, pKappaF, pKappaU,
619         O                         dPhiHydX, dPhiHydY,
620         I                         myTime, myIter, myThid)
621           ENDIF
622           DO j=jMin,jMax
623             DO i=iMin,iMax
624               phiHydF(i,j) = phiHydU(i,j)
625             ENDDO
626           ENDDO
627    
628    #endif /* DISABLE_SIGMA_CODE */
629    C--   end if-not/else useFVgradPhi
630          ENDIF
631    
632  C---   Diagnose Phi at boundary r=R_low :  C---   Diagnose Phi at boundary r=R_low :
633  C       = Ocean bottom pressure (Ocean, Z-coord.)  C       = Ocean bottom pressure (Ocean, Z-coord.)
634  C       = Sea-surface height    (Ocean, P-coord.)  C       = Sea-surface height    (Ocean, P-coord.)
# Line 565  C---   Diagnose Full Hydrostatic Potenti Line 646  C---   Diagnose Full Hydrostatic Potenti
646       I                      phiHydC,       I                      phiHydC,
647       I                      myTime, myIter, myThid)       I                      myTime, myIter, myThid)
648    
       IF (momPressureForcing) THEN  
         CALL CALC_GRAD_PHI_HYD(  
      I                         k, bi, bj, iMin,iMax, jMin,jMax,  
      I                         phiHydC, alphaRho, tFld, sFld,  
      O                         dPhiHydX, dPhiHydY,  
      I                         myTime, myIter, myThid)  
       ENDIF  
   
649  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
650    
651        RETURN        RETURN
Legend:
Removed from v.1.42  
changed lines
  Added in v.1.43
  ViewVC Help
Powered by ViewVC 1.1.22