/[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.44 by jmc, Fri Apr 4 20:54:11 2014 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "PACKAGES_CONFIG.h"  #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    #ifdef ALLOW_AUTODIFF
7    # include "AUTODIFF_OPTIONS.h"
8    #endif
9    
10  CBOP  CBOP
11  C     !ROUTINE: CALC_PHI_HYD  C     !ROUTINE: CALC_PHI_HYD
# Line 42  C     == Global variables == Line 45  C     == Global variables ==
45  #include "GRID.h"  #include "GRID.h"
46  #include "EEPARAMS.h"  #include "EEPARAMS.h"
47  #include "PARAMS.h"  #include "PARAMS.h"
48  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF
49  #include "tamc.h"  #include "tamc.h"
50  #include "tamc_keys.h"  #include "tamc_keys.h"
51  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF */
52  #include "SURFACE.h"  #include "SURFACE.h"
53  #include "DYNVARS.h"  #include "DYNVARS.h"
54    
# Line 65  C     myThid     :: thread number for th Line 68  C     myThid     :: thread number for th
68        INTEGER bi,bj,iMin,iMax,jMin,jMax,k        INTEGER bi,bj,iMin,iMax,jMin,jMax,k
69        _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)
70        _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)  
71        _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72        _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73        _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 79  c     _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy
79    
80  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
81  C     == Local variables ==  C     == Local variables ==
82    C     phiHydU    :: hydrostatic potential anomaly at interface k+1 (Upper k)
83    C     pKappaF    :: (p/Po)^kappa at interface k
84    C     pKappaU    :: (p/Po)^kappa at interface k+1 (Upper k)
85    C     pKappaC    :: (p/Po)^kappa at cell center k
86        INTEGER i,j        INTEGER i,j
87        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88    #ifndef DISABLE_SIGMA_CODE
89          _RL phiHydU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90          _RL pKappaF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91          _RL pKappaU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92          _RL pKappaC, locDepth, fullDepth
93    #endif /* DISABLE_SIGMA_CODE */
94        _RL dRlocM,dRlocP, ddRloc, locAlpha        _RL dRlocM,dRlocP, ddRloc, locAlpha
95        _RL ddPIm, ddPIp, rec_dRm, rec_dRp        _RL ddPIm, ddPIp, rec_dRm, rec_dRp
96        _RL surfPhiFac        _RL surfPhiFac
97        LOGICAL useDiagPhiRlow, addSurfPhiAnom        LOGICAL useDiagPhiRlow, addSurfPhiAnom
98          LOGICAL useFVgradPhi
99  CEOP  CEOP
100        useDiagPhiRlow = .TRUE.        useDiagPhiRlow = .TRUE.
101        addSurfPhiAnom = select_rStar.EQ.0 .AND. nonlinFreeSurf.GE.4        addSurfPhiAnom = select_rStar.EQ.0 .AND. nonlinFreeSurf.GE.4
102          useFVgradPhi   = selectSigmaCoord.NE.0
103    
104        surfPhiFac = 0.        surfPhiFac = 0.
105        IF (addSurfPhiAnom) surfPhiFac = 1.        IF (addSurfPhiAnom) surfPhiFac = 1.
106    
# Line 430  C--     Quasi-hydrostatic terms are adde Line 445  C--     Quasi-hydrostatic terms are adde
445    
446  C---    Integrate d Phi / d pi  C---    Integrate d Phi / d pi
447    
448    #ifndef DISABLE_SIGMA_CODE
449    C  --  Finite Volume Form, integrated to r-level (cell center & upper interface)
450           IF ( useFVgradPhi ) THEN
451    
452             fullDepth = rF(1)-rF(Nr+1)
453             DO j=jMin,jMax
454              DO i=iMin,iMax
455               locDepth = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
456    #ifdef NONLIN_FRSURF
457               locDepth = locDepth + surfPhiFac*etaH(i,j,bi,bj)
458    #endif
459               pKappaF(i,j) = (
460         &           ( R_low(i,j,bi,bj) + aHybSigmF( k )*fullDepth
461         &                              + bHybSigmF( k )*locDepth
462         &           )/atm_Po )**atm_kappa
463               pKappaC      = (
464         &           ( R_low(i,j,bi,bj) + aHybSigmC( k )*fullDepth
465         &                              + bHybSigmC( k )*locDepth
466         &           )/atm_Po )**atm_kappa
467               pKappaU(i,j) = (
468         &           ( R_low(i,j,bi,bj)+ aHybSigmF(k+1)*fullDepth
469         &                             + bHybSigmF(k+1)*locDepth
470         &           )/atm_Po )**atm_kappa
471               phiHydC(i,j) = phiHydF(i,j)
472         &        + atm_Cp*( pKappaF(i,j) - pKappaC      )*alphaRho(i,j)
473               phiHydU(i,j) = phiHydF(i,j)
474         &        + atm_Cp*( pKappaF(i,j) - pKappaU(i,j) )*alphaRho(i,j)
475              ENDDO
476             ENDDO
477    C end: Finite Volume Form, integrated to r-level.
478    
479           ELSEIF (integr_GeoPot.EQ.0) THEN
480    #else /* DISABLE_SIGMA_CODE */
481         IF (integr_GeoPot.EQ.0) THEN         IF (integr_GeoPot.EQ.0) THEN
482    #endif /* DISABLE_SIGMA_CODE */
483  C  --  Energy Conserving Form, accurate with Full cell topo  --  C  --  Energy Conserving Form, accurate with Full cell topo  --
484  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
485  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 597  C---+----1----+----2----+----3----+----4
597          STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'          STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
598        ENDIF        ENDIF
599    
600          IF ( .NOT. useFVgradPhi ) THEN
601    C--   r-coordinate and r*-coordinate cases:
602    
603           IF ( momPressureForcing ) THEN
604            CALL CALC_GRAD_PHI_HYD(
605         I                         k, bi, bj, iMin,iMax, jMin,jMax,
606         I                         phiHydC, alphaRho, tFld, sFld,
607         O                         dPhiHydX, dPhiHydY,
608         I                         myTime, myIter, myThid)
609           ENDIF
610    
611    #ifndef DISABLE_SIGMA_CODE
612          ELSE
613    C--   else (SigmaCoords part)
614    
615           IF ( fluidIsWater ) THEN
616            STOP 'CALC_PHI_HYD: missing code for SigmaCoord'
617           ENDIF
618           IF ( momPressureForcing ) THEN
619            CALL CALC_GRAD_PHI_FV(
620         I                         k, bi, bj, iMin,iMax, jMin,jMax,
621         I                         phiHydF, phiHydU, pKappaF, pKappaU,
622         O                         dPhiHydX, dPhiHydY,
623         I                         myTime, myIter, myThid)
624           ENDIF
625           DO j=jMin,jMax
626             DO i=iMin,iMax
627               phiHydF(i,j) = phiHydU(i,j)
628             ENDDO
629           ENDDO
630    
631    #endif /* DISABLE_SIGMA_CODE */
632    C--   end if-not/else useFVgradPhi
633          ENDIF
634    
635  C---   Diagnose Phi at boundary r=R_low :  C---   Diagnose Phi at boundary r=R_low :
636  C       = Ocean bottom pressure (Ocean, Z-coord.)  C       = Ocean bottom pressure (Ocean, Z-coord.)
637  C       = Sea-surface height    (Ocean, P-coord.)  C       = Sea-surface height    (Ocean, P-coord.)
# Line 565  C---   Diagnose Full Hydrostatic Potenti Line 649  C---   Diagnose Full Hydrostatic Potenti
649       I                      phiHydC,       I                      phiHydC,
650       I                      myTime, myIter, myThid)       I                      myTime, myIter, myThid)
651    
       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  
   
652  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
653    
654        RETURN        RETURN
Legend:
Removed from v.1.42  
changed lines
  Added in v.1.44
  ViewVC Help
Powered by ViewVC 1.1.22