model/src/calc_phi_hyd.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_phi_hyd.F,v 1.8.2.1 2001/01/12 21:02:46 adcroft Exp $

#include "CPP_OPTIONS.h"

      SUBROUTINE CALC_PHI_HYD(
     I                         bi, bj, iMin, iMax, jMin, jMax, K,
     I                         theta, salt,
     U                         phiHyd, phiHydInterface,
     I                         myThid)
C     /==========================================================\
C     | SUBROUTINE CALC_PHI_HYD                                  |
C     | o Integrate the hydrostatic relation to find phiHyd.     |
C     |                                                          |
C     | On entry:                                                |
C     |   theta,salt    are the current thermodynamics quantities|
C     |                 (unchanged on exit)                      |
C     |   phiHyd(i,j,1:k-1) is the hydrostatic pressure/geopot.  |
C     |                 at cell centers (tracer points)          |
C     |                 - 1:k-1 layers are valid                 |
C     |                 - k:Nr layers are invalid                |
C     |   phiHydInterface(i,j) is the hydrostatic pressure/geop. |
C     |                 at cell the interface k (w point above)  |
C     | On exit:                                                 |
C     |   phiHyd(i,j,1:k) is the hydrostatic pressure/geopot.    |
C     |                 at cell centers (tracer points)          |
C     |                 - 1:k layers are valid                   |
C     |                 - k+1:Nr layers are invalid              |
C     |   phiHydInterface(i,j) is the hydrostatic pressure/geop. |
C     |                 at cell the interface k+1 (w point below)|
C     |                                                          |
C     \==========================================================/
      IMPLICIT NONE
C     == Global variables ==
#include "SIZE.h"
#include "GRID.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
C     == Routine arguments ==
      INTEGER bi,bj,iMin,iMax,jMin,jMax,K
      _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL salt(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL phiHydInterface(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      INTEGER myThid

#ifdef INCLUDE_PHIHYD_CALCULATION_CODE

C     == Local variables ==
      INTEGER i,j
      _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL dRloc,dRlocKp1
      _RL ddRm1, ddRp1, ddRm, ddRp
      _RL atm_cp, atm_kappa, atm_po

      IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN

        dRloc=drC(k)
        IF (k.EQ.1) dRloc=drF(1)
        IF (k.EQ.Nr) THEN
          dRlocKp1=0.
        ELSE
          dRlocKp1=drC(k+1)
        ENDIF

C--     If this is the top layer we impose the boundary condition
C       P(z=eta) = P(atmospheric_loading)
        IF (k.EQ.1) THEN
          DO j=jMin,jMax
            DO i=iMin,iMax
C             *NOTE* The loading should go here but has not been implemented yet
              phiHydInterface(i,j)=0.
              phiHyd(i,j,k)=0.
            ENDDO
          ENDDO
        ENDIF

C       Calculate density
        CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
     &                 theta, salt,
     &                 alphaRho, myThid)

C       Hydrostatic pressure at cell centers
        DO j=jMin,jMax
          DO i=iMin,iMax
#ifdef ALLOW_AUTODIFF_TAMC
          Is this directive correct or even necessary in this new code?
CADJ GENERAL
#endif
C           This discretization is the "finite volume" form
C           which has not been used to date since it does not
C           conserve KE+PE exactly even though it is more natural
C
C           phiHyd(i,j,k)=phiHydInterface(i,j)+
C    &          0.5*drF(K)*gravity*alphaRho(i,j)
C
C           This discretization is the "energy conserving" form
C           which has been used since at least Adcroft et al., MWR 1997
C
            phiHyd(i,j,k)=phiHyd(i,j,k)+
     &          0.5*dRloc*gravity*alphaRho(i,j)
            phiHyd(i,j,k+1)=phiHyd(i,j,k)+
     &          0.5*dRlocKp1*gravity*alphaRho(i,j)
          ENDDO
        ENDDO
        
C       Hydrostatic pressure at interface below
        DO j=jMin,jMax
          DO i=iMin,iMax
            phiHydInterface(i,j)=phiHydInterface(i,j)+
     &              drF(K)*gravity*alphaRho(i,j)
          ENDDO
        ENDDO

      ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN

        atm_cp=1004. _d 0
        atm_kappa=2. _d 0/7. _d 0
        atm_po=1. _d 5
        IF (K.EQ.1) THEN
C       Integrate d Phi / d R
          ddRp1=atm_cp*( ((rC(K)/atm_po)**atm_kappa)
     &                  -((rF(K)/atm_po)**atm_kappa) )
          DO j=jMin,jMax
           DO i=iMin,iMax
             ddRp=ddRp1
             IF (hFacC(I,J, K ,bi,bj).EQ.0.) ddRp=0.
             phiHyd(i,j,K)=0.
     &          -ddRp*(theta(I,J,K,bi,bj)-tRef(K))
           ENDDO
          ENDDO
        ELSE
C       Integrate d Phi / d R
          ddRp1=atm_cp*( ((rC( K )/atm_po)**atm_kappa)
     &                  -((rF( K )/atm_po)**atm_kappa) )
          ddRm1=atm_cp*( ((rF( K )/atm_po)**atm_kappa)
     &                  -((rC(K-1)/atm_po)**atm_kappa) )
          DO j=jMin,jMax
           DO i=iMin,iMax
             ddRp=ddRp1
             ddRm=ddRm1
             IF (hFacC(I,J, K ,bi,bj).EQ.0.) ddRp=0.
             IF (hFacC(I,J,K-1,bi,bj).EQ.0.) ddRm=0.
             phiHyd(i,j,K)=phiHyd(i,j,K-1)
     &          -( ddRm*(theta(I,J,K-1,bi,bj)-tRef(K-1))
     &            +ddRp*(theta(I,J, K ,bi,bj)-tRef( K )) )
           ENDDO
          ENDDO
        ENDIF
                                                                                        
      ELSE
        STOP 'CALC_PHI_HYD: We should never reach this point!'
      ENDIF

#endif

      return
      end
1	jmc	1.8.2.2	C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_phi_hyd.F,v 1.8.2.1 2001/01/12 21:02:46 adcroft Exp $
2	cnh	1.1
3	cnh	1.6	#include "CPP_OPTIONS.h"
4	cnh	1.1
5	adcroft	1.8.2.1	SUBROUTINE CALC_PHI_HYD(
6			I bi, bj, iMin, iMax, jMin, jMax, K,
7			I theta, salt,
8			U phiHyd, phiHydInterface,
9			I myThid)
10	cnh	1.1	C /==========================================================\
11			C \| SUBROUTINE CALC_PHI_HYD \|
12			C \| o Integrate the hydrostatic relation to find phiHyd. \|
13			C \| \|
14	adcroft	1.8.2.1	C \| On entry: \|
15			C \| theta,salt are the current thermodynamics quantities\|
16			C \| (unchanged on exit) \|
17			C \| phiHyd(i,j,1:k-1) is the hydrostatic pressure/geopot. \|
18			C \| at cell centers (tracer points) \|
19			C \| - 1:k-1 layers are valid \|
20			C \| - k:Nr layers are invalid \|
21			C \| phiHydInterface(i,j) is the hydrostatic pressure/geop. \|
22			C \| at cell the interface k (w point above) \|
23			C \| On exit: \|
24			C \| phiHyd(i,j,1:k) is the hydrostatic pressure/geopot. \|
25			C \| at cell centers (tracer points) \|
26			C \| - 1:k layers are valid \|
27			C \| - k+1:Nr layers are invalid \|
28			C \| phiHydInterface(i,j) is the hydrostatic pressure/geop. \|
29			C \| at cell the interface k+1 (w point below)\|
30			C \| \|
31	cnh	1.1	C \==========================================================/
32			IMPLICIT NONE
33			C == Global variables ==
34			#include "SIZE.h"
35			#include "GRID.h"
36			#include "EEPARAMS.h"
37			#include "PARAMS.h"
38			C == Routine arguments ==
39			INTEGER bi,bj,iMin,iMax,jMin,jMax,K
40	adcroft	1.8.2.1	_RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
41			_RL salt(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
42	cnh	1.2	_RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
43	adcroft	1.8.2.1	_RL phiHydInterface(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
44			INTEGER myThid
45	cnh	1.1
46	cnh	1.6	#ifdef INCLUDE_PHIHYD_CALCULATION_CODE
47
48	adcroft	1.8.2.1	C == Local variables ==
49			INTEGER i,j
50			_RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
51			_RL dRloc,dRlocKp1
52	jmc	1.8.2.2	_RL ddRm1, ddRp1, ddRm, ddRp
53			_RL atm_cp, atm_kappa, atm_po
54	adcroft	1.8.2.1
55			IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN
56
57			dRloc=drC(k)
58			IF (k.EQ.1) dRloc=drF(1)
59			IF (k.EQ.Nr) THEN
60			dRlocKp1=0.
61			ELSE
62			dRlocKp1=drC(k+1)
63			ENDIF
64
65			C-- If this is the top layer we impose the boundary condition
66			C P(z=eta) = P(atmospheric_loading)
67			IF (k.EQ.1) THEN
68			DO j=jMin,jMax
69			DO i=iMin,iMax
70			C NOTE The loading should go here but has not been implemented yet
71			phiHydInterface(i,j)=0.
72			phiHyd(i,j,k)=0.
73			ENDDO
74			ENDDO
75			ENDIF
76
77			C Calculate density
78			CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
79			& theta, salt,
80			& alphaRho, myThid)
81
82			C Hydrostatic pressure at cell centers
83			DO j=jMin,jMax
84			DO i=iMin,iMax
85	heimbach	1.7	#ifdef ALLOW_AUTODIFF_TAMC
86	adcroft	1.8.2.1	Is this directive correct or even necessary in this new code?
87	heimbach	1.7	CADJ GENERAL
88			#endif
89	adcroft	1.8.2.1	C This discretization is the "finite volume" form
90			C which has not been used to date since it does not
91			C conserve KE+PE exactly even though it is more natural
92			C
93			C phiHyd(i,j,k)=phiHydInterface(i,j)+
94			C & 0.5drF(K)gravity*alphaRho(i,j)
95			C
96			C This discretization is the "energy conserving" form
97			C which has been used since at least Adcroft et al., MWR 1997
98			C
99			phiHyd(i,j,k)=phiHyd(i,j,k)+
100			& 0.5dRlocgravity*alphaRho(i,j)
101			phiHyd(i,j,k+1)=phiHyd(i,j,k)+
102			& 0.5dRlocKp1gravity*alphaRho(i,j)
103			ENDDO
104			ENDDO
105
106			C Hydrostatic pressure at interface below
107			DO j=jMin,jMax
108			DO i=iMin,iMax
109			phiHydInterface(i,j)=phiHydInterface(i,j)+
110			& drF(K)gravityalphaRho(i,j)
111			ENDDO
112			ENDDO
113
114			ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
115
116	jmc	1.8.2.2	atm_cp=1004. _d 0
117			atm_kappa=2. _d 0/7. _d 0
118			atm_po=1. _d 5
119			IF (K.EQ.1) THEN
120			C Integrate d Phi / d R
121			ddRp1=atm_cp( ((rC(K)/atm_po)*atm_kappa)
122			& -((rF(K)/atm_po)**atm_kappa) )
123			DO j=jMin,jMax
124			DO i=iMin,iMax
125			ddRp=ddRp1
126			IF (hFacC(I,J, K ,bi,bj).EQ.0.) ddRp=0.
127			phiHyd(i,j,K)=0.
128			& -ddRp*(theta(I,J,K,bi,bj)-tRef(K))
129			ENDDO
130			ENDDO
131			ELSE
132			C Integrate d Phi / d R
133			ddRp1=atm_cp( ((rC( K )/atm_po)*atm_kappa)
134			& -((rF( K )/atm_po)**atm_kappa) )
135			ddRm1=atm_cp( ((rF( K )/atm_po)*atm_kappa)
136			& -((rC(K-1)/atm_po)**atm_kappa) )
137			DO j=jMin,jMax
138			DO i=iMin,iMax
139			ddRp=ddRp1
140			ddRm=ddRm1
141			IF (hFacC(I,J, K ,bi,bj).EQ.0.) ddRp=0.
142			IF (hFacC(I,J,K-1,bi,bj).EQ.0.) ddRm=0.
143			phiHyd(i,j,K)=phiHyd(i,j,K-1)
144			& -( ddRm*(theta(I,J,K-1,bi,bj)-tRef(K-1))
145			& +ddRp*(theta(I,J, K ,bi,bj)-tRef( K )) )
146			ENDDO
147			ENDDO
148			ENDIF
149
150	adcroft	1.8.2.1	ELSE
151			STOP 'CALC_PHI_HYD: We should never reach this point!'
152			ENDIF
153	cnh	1.1
154	cnh	1.6	#endif
155
156	cnh	1.1	return
157			end