model/src/diags_phi_hyd.F

C $Header: /u/gcmpack/MITgcm/model/src/diags_phi_hyd.F,v 1.3 2011/06/10 20:54:54 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: DIAGS_PHI_HYD
C     !INTERFACE:
      SUBROUTINE DIAGS_PHI_HYD(
     I                       k, bi, bj, iMin,iMax, jMin,jMax,
     I                       phiHydC,
     I                       myTime, myIter, myThid)
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R DIAGS_PHI_HYD
C     | o Diagnose full hydrostatic Potential at cell center ;
C     |   used for output & with EOS funct. of P
C     *==========================================================*
C     | NOTE: For now, only contains the (total) Potential anomaly
C     |  since phiRef (for Atmos) is not available (not in common)
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
#include "DYNVARS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine Arguments ==
C     k, bi,bj      :: level & tile indices
C     iMin,iMax,jMin,jMax :: Loop counters
C     phiHydC    :: hydrostatic potential anomaly at cell center
C                  (atmos: =Geopotential ; ocean-z: =Pressure/rho)
C     myTime :: Current time
C     myIter :: Current iteration number
C     myThid :: Instance number for this call of the routine.
      INTEGER k, bi,bj, iMin,iMax, jMin,jMax
      _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL myTime
      INTEGER myIter, myThid

#ifdef INCLUDE_PHIHYD_CALCULATION_CODE

C     !LOCAL VARIABLES:
C     == Local variables ==
C     i,j        :: Loop counters
C     phiHydCstR :: total hydrostatic Potential (anomaly, for now),
C                   at fixed r-position, cell center level location.
      INTEGER i,j
#ifdef NONLIN_FRSURF
      _RL facP, dPhiRef
      _RL phiHydCstR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#endif /* NONLIN_FRSURF */
CEOP

        DO j=jMin,jMax
         DO i=iMin,iMax
           totPhiHyd(i,j,k,bi,bj) = phiHydC(i,j)
     &            + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
     &            + phi0surf(i,j,bi,bj)
#ifdef NONLIN_FRSURF
           phiHydCstR(i,j) = totPhiHyd(i,j,k,bi,bj)
#endif /* NONLIN_FRSURF */
         ENDDO
        ENDDO

#ifdef NONLIN_FRSURF
c     IF (select_rStar.GE.2 .AND. nonlinFreeSurf.GE.4 ) THEN
      IF (select_rStar.GE.1 .AND. nonlinFreeSurf.GE.4 ) THEN
c# ifndef DISABLE_RSTAR_CODE
C-    Integral of b.dr = rStarFac * Integral of b.dr* :
       IF ( fluidIsAir ) THEN
C-     Consistent with Phi'= Integr[ theta'.dPi ] :
        DO j=jMin,jMax
         DO i=iMin,iMax
           facP = rStarFacC(i,j,bi,bj)**atm_kappa
           dPhiRef = phiRef(2*k) - gravity*topoZ(i,j,bi,bj)
     &                           - phi0surf(i,j,bi,bj)
           totPhiHyd(i,j,k,bi,bj) =
     &              phiHydC(i,j)*facP
     &            + MAX( dPhiRef, 0. _d 0 )*( facP - 1. _d 0 )
     &            + phi0surf(i,j,bi,bj)
c          phiHydCstR(i,j) = phiHydCstR(i,j)
c    &            + phiHydC(i,j)*( facP - 1. _d 0 )
         ENDDO
        ENDDO
       ELSEIF ( usingPCoords ) THEN
        DO j=jMin,jMax
         DO i=iMin,iMax
c    &     dPhiRef = phiRef(2*k) - gravity*topoZ(i,j,bi,bj)
c    &                           - phi0surf(i,j,bi,bj)
C--    assume PhiRef is just (ps0 - p)/rhoConst :
           dPhiRef =( Ro_surf(i,j,bi,bj)-rC(k) )*recip_rhoConst
           totPhiHyd(i,j,k,bi,bj) =
     &              phiHydC(i,j)*rStarFacC(i,j,bi,bj)
     &            + MAX( dPhiRef, 0. _d 0 )
     &                *( rStarFacC(i,j,bi,bj) - 1. _d 0 )
     &            + phi0surf(i,j,bi,bj)
c          totPhiHyd(i,j,k,bi,bj) = phiHydCstR(i,j)
         ENDDO
        ENDDO
       ELSE
        DO j=jMin,jMax
         DO i=iMin,iMax
           dPhiRef =( Ro_surf(i,j,bi,bj)-rC(k) )*gravity
           totPhiHyd(i,j,k,bi,bj) =
     &              phiHydC(i,j)*rStarFacC(i,j,bi,bj)
     &            + MAX( dPhiRef, 0. _d 0 )
     &                *( rStarFacC(i,j,bi,bj) - 1. _d 0 )
     &            + phi0surf(i,j,bi,bj)
c          totPhiHyd(i,j,k,bi,bj) = phiHydCstR(i,j)
         ENDDO
        ENDDO
       ENDIF
       CALL DIAGNOSTICS_FILL(phiHydCstR,'PHIHYDcR',k,1,2,bi,bj,myThid)
c# endif /* DISABLE_RSTAR_CODE */
      ENDIF
#endif /* NONLIN_FRSURF */

#endif /* INCLUDE_PHIHYD_CALCULATION_CODE */

      RETURN
      END
1	jmc	1.4	C $Header: /u/gcmpack/MITgcm/model/src/diags_phi_hyd.F,v 1.3 2011/06/10 20:54:54 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "CPP_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: DIAGS_PHI_HYD
8			C !INTERFACE:
9	jmc	1.3	SUBROUTINE DIAGS_PHI_HYD(
10	jmc	1.1	I k, bi, bj, iMin,iMax, jMin,jMax,
11			I phiHydC,
12			I myTime, myIter, myThid)
13			C !DESCRIPTION: \bv
14			C ==========================================================
15	jmc	1.3	C \| S/R DIAGS_PHI_HYD
16			C \| o Diagnose full hydrostatic Potential at cell center ;
17	jmc	1.1	C \| used for output & with EOS funct. of P
18			C ==========================================================
19			C \| NOTE: For now, only contains the (total) Potential anomaly
20			C \| since phiRef (for Atmos) is not available (not in common)
21			C ==========================================================
22			C \ev
23
24			C !USES:
25			IMPLICIT NONE
26			C == Global variables ==
27			#include "SIZE.h"
28			#include "EEPARAMS.h"
29			#include "PARAMS.h"
30			#include "GRID.h"
31			#include "SURFACE.h"
32			#include "DYNVARS.h"
33
34			C !INPUT/OUTPUT PARAMETERS:
35			C == Routine Arguments ==
36			C k, bi,bj :: level & tile indices
37			C iMin,iMax,jMin,jMax :: Loop counters
38			C phiHydC :: hydrostatic potential anomaly at cell center
39			C (atmos: =Geopotential ; ocean-z: =Pressure/rho)
40			C myTime :: Current time
41			C myIter :: Current iteration number
42			C myThid :: Instance number for this call of the routine.
43			INTEGER k, bi,bj, iMin,iMax, jMin,jMax
44			_RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
45			_RL myTime
46			INTEGER myIter, myThid
47
48			#ifdef INCLUDE_PHIHYD_CALCULATION_CODE
49
50			C !LOCAL VARIABLES:
51			C == Local variables ==
52	jmc	1.4	C i,j :: Loop counters
53			C phiHydCstR :: total hydrostatic Potential (anomaly, for now),
54			C at fixed r-position, cell center level location.
55	jmc	1.1	INTEGER i,j
56	jmc	1.4	#ifdef NONLIN_FRSURF
57			_RL facP, dPhiRef
58			_RL phiHydCstR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59			#endif /* NONLIN_FRSURF */
60	jmc	1.1	CEOP
61	jmc	1.3
62	jmc	1.4	DO j=jMin,jMax
63			DO i=iMin,iMax
64			totPhiHyd(i,j,k,bi,bj) = phiHydC(i,j)
65			& + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
66			& + phi0surf(i,j,bi,bj)
67			#ifdef NONLIN_FRSURF
68			phiHydCstR(i,j) = totPhiHyd(i,j,k,bi,bj)
69			#endif /* NONLIN_FRSURF */
70			ENDDO
71			ENDDO
72
73	jmc	1.3	#ifdef NONLIN_FRSURF
74			c IF (select_rStar.GE.2 .AND. nonlinFreeSurf.GE.4 ) THEN
75			IF (select_rStar.GE.1 .AND. nonlinFreeSurf.GE.4 ) THEN
76			c# ifndef DISABLE_RSTAR_CODE
77			C- Integral of b.dr = rStarFac * Integral of b.dr* :
78	jmc	1.4	IF ( fluidIsAir ) THEN
79	jmc	1.3	C- Consistent with Phi'= Integr[ theta'.dPi ] :
80			DO j=jMin,jMax
81			DO i=iMin,iMax
82	jmc	1.4	facP = rStarFacC(i,j,bi,bj)**atm_kappa
83			dPhiRef = phiRef(2k) - gravitytopoZ(i,j,bi,bj)
84			& - phi0surf(i,j,bi,bj)
85	jmc	1.3	totPhiHyd(i,j,k,bi,bj) =
86	jmc	1.4	& phiHydC(i,j)*facP
87			& + MAX( dPhiRef, 0. _d 0 )*( facP - 1. _d 0 )
88	jmc	1.3	& + phi0surf(i,j,bi,bj)
89	jmc	1.4	c phiHydCstR(i,j) = phiHydCstR(i,j)
90			c & + phiHydC(i,j)*( facP - 1. _d 0 )
91	jmc	1.3	ENDDO
92			ENDDO
93	jmc	1.4	ELSEIF ( usingPCoords ) THEN
94	jmc	1.3	DO j=jMin,jMax
95			DO i=iMin,iMax
96	jmc	1.4	c & dPhiRef = phiRef(2k) - gravitytopoZ(i,j,bi,bj)
97			c & - phi0surf(i,j,bi,bj)
98			C-- assume PhiRef is just (ps0 - p)/rhoConst :
99			dPhiRef =( Ro_surf(i,j,bi,bj)-rC(k) )*recip_rhoConst
100	jmc	1.3	totPhiHyd(i,j,k,bi,bj) =
101			& phiHydC(i,j)*rStarFacC(i,j,bi,bj)
102	jmc	1.4	& + MAX( dPhiRef, 0. _d 0 )
103			& *( rStarFacC(i,j,bi,bj) - 1. _d 0 )
104	jmc	1.3	& + phi0surf(i,j,bi,bj)
105	jmc	1.4	c totPhiHyd(i,j,k,bi,bj) = phiHydCstR(i,j)
106	jmc	1.3	ENDDO
107			ENDDO
108	jmc	1.4	ELSE
109	jmc	1.1	DO j=jMin,jMax
110			DO i=iMin,iMax
111	jmc	1.4	dPhiRef =( Ro_surf(i,j,bi,bj)-rC(k) )*gravity
112			totPhiHyd(i,j,k,bi,bj) =
113			& phiHydC(i,j)*rStarFacC(i,j,bi,bj)
114			& + MAX( dPhiRef, 0. _d 0 )
115			& *( rStarFacC(i,j,bi,bj) - 1. _d 0 )
116	jmc	1.3	& + phi0surf(i,j,bi,bj)
117	jmc	1.4	c totPhiHyd(i,j,k,bi,bj) = phiHydCstR(i,j)
118	jmc	1.1	ENDDO
119			ENDDO
120	jmc	1.4	ENDIF
121			CALL DIAGNOSTICS_FILL(phiHydCstR,'PHIHYDcR',k,1,2,bi,bj,myThid)
122			c# endif /* DISABLE_RSTAR_CODE */
123	jmc	1.3	ENDIF
124	jmc	1.4	#endif /* NONLIN_FRSURF */
125	jmc	1.1
126			#endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
127
128			RETURN
129			END