press_release/code/pressure_release_theta.F

C $Header: /u/gcmpack/MITgcm_contrib/verification_other/../pressure_release_thermodynamics.F,v 1.2 2016/12/12 11:32:17 ksnow Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
      SUBROUTINE PRESSURE_RELEASE_THETA( 
     U                       gT_arr,
     I                       iMin,iMax,jMin,jMax, k, bi,bj,
     I                       myTime, myIter, myThid )
C     *============================================================*
C     | SUBROUTINE PRESSURE_RELEASE_THETA
C     | o Transport theta with darcy flux
C     *============================================================*
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "SURFACE.h"
#include "FFIELDS.h"

C     === Routine arguments ===
C     myThid -  Number of this instance 

      _RL     gT_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      INTEGER k,bi,bj
      INTEGER iMin,iMax,jMin,jMax
      _RL     myTime
      INTEGER myIter
      INTEGER myThid
      INTEGER i,j,k_e,k_ce,k_s,k_cs,k_w,k_cw,k_n,k_cn
      _RL     T_trans_west,T_trans_east,T_trans_south,T_trans_north

CEndOfInterface

C Need to find if adjacent cells are deeper or shallower


C Mass flux into one side of cell is 
      DO j=jMin+1,jMax-1
        DO i=iMin+1,iMax-1

C calculate the k cells the tracers are transferred between in north,
C south east and west cells.
          IF (kLowC(i,j,bi,bj) .GE. kLowC(i+1,j,bi,bj)) THEN
             k_e = kLowC(i+1,j,bi,bj)
             k_ce = kSurfC(i,j,bi,bj)
          ELSE
             k_e = kSurfC(i+1,j,bi,bj)
             k_ce = kLowC(i,j,bi,bj)
          ENDIF

          IF (kLowC(i,j,bi,bj) .GE. kLowC(i-1,j,bi,bj)) THEN
             k_w = kLowC(i-1,j,bi,bj)
             k_cw = kSurfC(i,j,bi,bj)
          ELSE
             k_w = kSurfC(i-1,j,bi,bj)
             k_cw = kLowC(i,j,bi,bj)
          ENDIF

          IF (kLowC(i,j,bi,bj) .GE. kLowC(i,j+1,bi,bj)) THEN
             k_n = kLowC(i,j+1,bi,bj)
             k_cn = kSurfC(i,j,bi,bj)
          ELSE
             k_n = kSurfC(i,j+1,bi,bj)
             k_cn = kLowC(i,j,bi,bj)
          ENDIF

          IF (kLowC(i,j,bi,bj) .GE. kLowC(i,j-1,bi,bj)) THEN
             k_s = kLowC(i,j-1,bi,bj)
             k_cs = kSurfC(i,j,bi,bj)
          ELSE
             k_s = kSurfC(i,j-1,bi,bj)
             k_cs = kLowC(i,j,bi,bj)
          ENDIF

C calculate the net tracer flux through north, south east and west
C faces.
          T_trans_west =pReleaseTransX(i,j,bi,bj)*
     &                  (theta(i-1,j,k_w,bi,bj) -theta(i,j,k_cw,bi,bj))
C     &                  *rhoFacC(k)*mass2rUnit
C     &                  *_dyG(i,j,bi,bj)*recip_rA(i,j,bi,bj)*
     &                  *recip_dxG(i,j,bi,bj)
     &                  *recip_drF(k_cw)*_recip_hFacC(i,j,k_cw,bi,bj)

          T_trans_east =pReleaseTransX(i+1,j,bi,bj)*
     &                  (theta(i,j,k_ce,bi,bj) -theta(i+1,j,k_e,bi,bj))
     &                  *recip_dxG(i+1,j,bi,bj)
     &                  *recip_drF(k_ce)*_recip_hFacC(i+1,j,k_ce,bi,bj)

          T_trans_south =pReleaseTransY(i,j,bi,bj)*
     &                  (theta(i,j-1,k_s,bi,bj) -theta(i,j,k_cs,bi,bj))
     &                  *recip_dyG(i,j,bi,bj)
     &                  *recip_drF(k_cs)*_recip_hFacC(i,j,k_cs,bi,bj)

          T_trans_north =pReleaseTransY(i,j+1,bi,bj)*
     &                  (theta(i,j,k_cn,bi,bj) -theta(i,j+1,k_n,bi,bj))
     &                  *recip_dyG(i,j+1,bi,bj)
     &                  *recip_drF(k_cn)*_recip_hFacC(i,j+1,k_cn,bi,bj)

C Add to get total tracer tendency.
          gT_arr(i,j) = gT_arr(i,j) + T_trans_west - T_trans_east
     &                 + T_trans_south - T_trans_north       

         ENDDO
       ENDDO

      RETURN
      END
1	ksnow	1.1	C $Header: /u/gcmpack/MITgcm_contrib/verification_other/../pressure_release_thermodynamics.F,v 1.2 2016/12/12 11:32:17 ksnow Exp $
2			C $Name: $
3
4			#include "PACKAGES_CONFIG.h"
5			#include "CPP_OPTIONS.h"
6
7			CBOP
8			SUBROUTINE PRESSURE_RELEASE_THETA(
9			U gT_arr,
10			I iMin,iMax,jMin,jMax, k, bi,bj,
11			I myTime, myIter, myThid )
12			C ============================================================
13			C \| SUBROUTINE PRESSURE_RELEASE_THETA
14			C \| o Transport theta with darcy flux
15			C ============================================================
16			IMPLICIT NONE
17
18			C === Global variables ===
19			#include "SIZE.h"
20			#include "EEPARAMS.h"
21			#include "PARAMS.h"
22			#include "GRID.h"
23			#include "DYNVARS.h"
24			#include "SURFACE.h"
25			#include "FFIELDS.h"
26
27			C === Routine arguments ===
28			C myThid - Number of this instance
29
30			_RL gT_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31			INTEGER k,bi,bj
32			INTEGER iMin,iMax,jMin,jMax
33			_RL myTime
34			INTEGER myIter
35			INTEGER myThid
36			INTEGER i,j,k_e,k_ce,k_s,k_cs,k_w,k_cw,k_n,k_cn
37			_RL T_trans_west,T_trans_east,T_trans_south,T_trans_north
38
39			CEndOfInterface
40
41			C Need to find if adjacent cells are deeper or shallower
42
43
44			C Mass flux into one side of cell is
45			DO j=jMin+1,jMax-1
46			DO i=iMin+1,iMax-1
47
48			C calculate the k cells the tracers are transferred between in north,
49			C south east and west cells.
50			IF (kLowC(i,j,bi,bj) .GE. kLowC(i+1,j,bi,bj)) THEN
51			k_e = kLowC(i+1,j,bi,bj)
52			k_ce = kSurfC(i,j,bi,bj)
53			ELSE
54			k_e = kSurfC(i+1,j,bi,bj)
55			k_ce = kLowC(i,j,bi,bj)
56			ENDIF
57
58			IF (kLowC(i,j,bi,bj) .GE. kLowC(i-1,j,bi,bj)) THEN
59			k_w = kLowC(i-1,j,bi,bj)
60			k_cw = kSurfC(i,j,bi,bj)
61			ELSE
62			k_w = kSurfC(i-1,j,bi,bj)
63			k_cw = kLowC(i,j,bi,bj)
64			ENDIF
65
66			IF (kLowC(i,j,bi,bj) .GE. kLowC(i,j+1,bi,bj)) THEN
67			k_n = kLowC(i,j+1,bi,bj)
68			k_cn = kSurfC(i,j,bi,bj)
69			ELSE
70			k_n = kSurfC(i,j+1,bi,bj)
71			k_cn = kLowC(i,j,bi,bj)
72			ENDIF
73
74			IF (kLowC(i,j,bi,bj) .GE. kLowC(i,j-1,bi,bj)) THEN
75			k_s = kLowC(i,j-1,bi,bj)
76			k_cs = kSurfC(i,j,bi,bj)
77			ELSE
78			k_s = kSurfC(i,j-1,bi,bj)
79			k_cs = kLowC(i,j,bi,bj)
80			ENDIF
81
82			C calculate the net tracer flux through north, south east and west
83			C faces.
84			T_trans_west =pReleaseTransX(i,j,bi,bj)*
85			& (theta(i-1,j,k_w,bi,bj) -theta(i,j,k_cw,bi,bj))
86			C & rhoFacC(k)mass2rUnit
87			C & _dyG(i,j,bi,bj)recip_rA(i,j,bi,bj)*
88			& *recip_dxG(i,j,bi,bj)
89			& recip_drF(k_cw)_recip_hFacC(i,j,k_cw,bi,bj)
90
91			T_trans_east =pReleaseTransX(i+1,j,bi,bj)*
92			& (theta(i,j,k_ce,bi,bj) -theta(i+1,j,k_e,bi,bj))
93			& *recip_dxG(i+1,j,bi,bj)
94			& recip_drF(k_ce)_recip_hFacC(i+1,j,k_ce,bi,bj)
95
96			T_trans_south =pReleaseTransY(i,j,bi,bj)*
97			& (theta(i,j-1,k_s,bi,bj) -theta(i,j,k_cs,bi,bj))
98			& *recip_dyG(i,j,bi,bj)
99			& recip_drF(k_cs)_recip_hFacC(i,j,k_cs,bi,bj)
100
101			T_trans_north =pReleaseTransY(i,j+1,bi,bj)*
102			& (theta(i,j,k_cn,bi,bj) -theta(i,j+1,k_n,bi,bj))
103			& *recip_dyG(i,j+1,bi,bj)
104			& recip_drF(k_cn)_recip_hFacC(i,j+1,k_cn,bi,bj)
105
106			C Add to get total tracer tendency.
107			gT_arr(i,j) = gT_arr(i,j) + T_trans_west - T_trans_east
108			& + T_trans_south - T_trans_north
109
110			ENDDO
111			ENDDO
112
113			RETURN
114			END