fix porous flux tracers

dngoldberg · Apr 9, 2021 · 8c1a266 · 8c1a266
1 parent 34f1fdc
commit 8c1a266
Show file tree

Hide file tree

Showing 3 changed files with 179 additions and 22 deletions.
diff --git a/model/src/integrate_for_w.F b/model/src/integrate_for_w.F
@@ -23,6 +23,10 @@ SUBROUTINE INTEGRATE_FOR_W(
 #include "PARAMS.h"
 #include "GRID.h"
 #include "SURFACE.h"
+#ifdef ALLOW_PRESSURE_RELEASE_CODE
+#include "DYNVARS.h"
+#endif
+
 
 C !INPUT/OUTPUT PARAMETERS:
 C == Routine arguments ==
@@ -54,6 +58,10 @@ SUBROUTINE INTEGRATE_FOR_W(
  _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
  _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
  _RL conv2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+#ifdef ALLOW_PRESSURE_RELEASE_CODE
+! _RL rTmp
+ INTEGER kl, ks
+#endif
 CEOP
 
 C-- Calculate velocity field "volume transports" through
@@ -84,6 +92,81 @@ SUBROUTINE INTEGRATE_FOR_W(
  ENDDO
  ENDIF
 #endif /* ALLOW_ADDFLUID */
+#ifdef ALLOW_PRESSURE_RELEASE_CODE
+ DO j=1,sNy
+ DO i=1,sNx
+ kl = kLowC(i,j,bi,bj)
+ ks = kSurfC(i,j,bi,bj)
+
+ if (ks .gt. 0) then
+
+ if ((k.eq.ks).and.(kl.eq.ks)) then
+
+ conv2d(i,j) = conv2d(i,j) +
+ & _dyG(i,j,bi,bj)*
+ & pReleaseTransX(i,j,bi,bj)
+ conv2d(i,j) = conv2d(i,j) -
+ & _dyG(i+1,j,bi,bj)*
+ & pReleaseTransX(i+1,j,bi,bj)
+ conv2d(i,j) = conv2d(i,j) +
+ & _dxG(i,j,bi,bj)*
+ & pReleaseTransY(i,j,bi,bj)
+ conv2d(i,j) = conv2d(i,j) -
+ & _dxG(i,j+1,bi,bj)*
+ & pReleaseTransY(i+1,j,bi,bj)
+
+ elseif (k.eq.kl) then
+
+ if (k.le.kSurfC(i-1,j,bi,bj)) then
+ conv2d(i,j) = conv2d(i,j) +
+ & _dyG(i,j,bi,bj)*
+ & pReleaseTransX(i,j,bi,bj)
+ endif
+ if (k.le.kSurfC(i+1,j,bi,bj)) then
+ conv2d(i,j) = conv2d(i,j) -
+ & _dyG(i+1,j,bi,bj)*
+ & pReleaseTransX(i+1,j,bi,bj)
+ endif
+ if (k.le.kSurfC(i,j-1,bi,bj)) then
+ conv2d(i,j) = conv2d(i,j) +
+ & _dxG(i,j,bi,bj)*
+ & pReleaseTransY(i,j,bi,bj)
+ endif
+ if (k.le.kSurfC(i,j+1,bi,bj)) then
+ conv2d(i,j) = conv2d(i,j) -
+ & _dxG(i,j+1,bi,bj)*
+ & pReleaseTransY(i,j+1,bi,bj)
+ endif
+
+ elseif (k.eq.ks) then
+
+ if (k.ge.kLowC(i-1,j,bi,bj)) then
+ conv2d(i,j) = conv2d(i,j) +
+ & _dyG(i,j,bi,bj)*
+ & pReleaseTransX(i,j,bi,bj)
+ endif
+ if (k.ge.kLowC(i+1,j,bi,bj)) then
+ conv2d(i,j) = conv2d(i,j) -
+ & _dyG(i+1,j,bi,bj)*
+ & pReleaseTransX(i+1,j,bi,bj)
+ endif
+ if (k.ge.kLowC(i,j-1,bi,bj)) then
+ conv2d(i,j) = conv2d(i,j) +
+ & _dxG(i,j,bi,bj)*
+ & pReleaseTransY(i,j,bi,bj)
+ endif
+ if (k.ge.kLowC(i,j+1,bi,bj)) then
+ conv2d(i,j) = conv2d(i,j) -
+ & _dxG(i,j+1,bi,bj)*
+ & pReleaseTransY(i,j+1,bi,bj)
+ endif
+
+ endif
+ endif
+
+ ENDDO
+ ENDDO
+#endif
 
 C-- Calculate vertical "volume transport" through face k
 C between tracer cell k-1 & k

diff --git a/model/src/pressure_release_salt.F b/model/src/pressure_release_salt.F
@@ -38,6 +38,8 @@ SUBROUTINE PRESSURE_RELEASE_SALT(
 C === Local Variables ===
  INTEGER i,j,k_e,k_ce,k_s,k_cs,k_w,k_cw,k_n,k_cn
  _RL S_trans_west,S_trans_east,S_trans_south,S_trans_north
+ _RL vol_trans_west,vol_trans_east,
+ & vol_trans_south,vol_trans_north
 
  DO j=jMin+1,jMax-1
  DO i=iMin+1,iMax-1
@@ -46,35 +48,39 @@ SUBROUTINE PRESSURE_RELEASE_SALT(
  S_trans_north = 0.0
  S_trans_east = 0.0 
  S_trans_south = 0.0
+ vol_trans_west = 0.0 
+ vol_trans_north = 0.0
+ vol_trans_east = 0.0 
+ vol_trans_south = 0.0
 
 C calculate the k cells the tracers are transferred between in north,
 C south east and west cells.
 C Need to find if adjacent cells are deeper or shallower
- IF (kLowC(i,j,bi,bj) .GE. kLowC(i+1,j,bi,bj)) THEN
+ IF (kLowC(i,j,bi,bj) .GT. ksurfC(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
+ IF (kLowC(i,j,bi,bj) .GT. ksurfC(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
+ IF (kLowC(i,j,bi,bj) .GT. ksurfC(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
+ IF (kLowC(i,j,bi,bj) .GT. ksurfC(i,j-1,bi,bj)) THEN
  k_s = kLowC(i,j-1,bi,bj)
  k_cs = kSurfC(i,j,bi,bj)
  ELSE
@@ -85,14 +91,19 @@ SUBROUTINE PRESSURE_RELEASE_SALT(
 C calculate the net tracer flux through north, south east and west
 C faces.
 
+
  IF (k .EQ. k_cw) THEN
  if (k_cw.gt.0 .and. k_w.gt.0) then
  S_trans_west =0.5 _d 0 * 
  & ( pReleaseTransX(i,j,bi,bj) *
  & (salt(i-1,j,k_w,bi,bj)+salt(i,j,k_cw,bi,bj))
  & +abs(pReleaseTransX(i,j,bi,bj)) * 
  & (salt(i-1,j,k_w,bi,bj)-salt(i,j,k_cw,bi,bj)))
- & *recip_rA(i,j,bi,bj)
+ & *recip_rA(i,j,bi,bj)*_dyG(i,j,bi,bj)
+ & *recip_deepFac2C(k_cw)*recip_rhoFacC(k_cw)
+ & *recip_drF(k_cw)*_recip_hFacC(i,j,k_cw,bi,bj)
+ vol_trans_west = pReleaseTransX(i,j,bi,bj)
+ & *recip_rA(i,j,bi,bj)*_dyG(i,j,bi,bj)
  & *recip_deepFac2C(k_cw)*recip_rhoFacC(k_cw)
  & *recip_drF(k_cw)*_recip_hFacC(i,j,k_cw,bi,bj)
  endif
@@ -104,7 +115,11 @@ SUBROUTINE PRESSURE_RELEASE_SALT(
  & (salt(i,j,k_ce,bi,bj)+salt(i+1,j,k_e,bi,bj))
  & +abs(pReleaseTransX(i+1,j,bi,bj)) *
  & (salt(i,j,k_ce,bi,bj)-salt(i+1,j,k_e,bi,bj)))
- & *recip_rA(i,j,bi,bj)
+ & *recip_rA(i,j,bi,bj)*_dyG(i+1,j,bi,bj)
+ & *recip_deepFac2C(k_ce)*recip_rhoFacC(k_ce)
+ & *recip_drF(k_ce)*_recip_hFacC(i,j,k_ce,bi,bj)
+ vol_trans_east = pReleaseTransX(i+1,j,bi,bj)
+ & *recip_rA(i,j,bi,bj)*_dyG(i+1,j,bi,bj)
  & *recip_deepFac2C(k_ce)*recip_rhoFacC(k_ce)
  & *recip_drF(k_ce)*_recip_hFacC(i,j,k_ce,bi,bj)
 ! S_trans_east =pReleaseTransX(i+1,j,bi,bj)*
@@ -120,7 +135,11 @@ SUBROUTINE PRESSURE_RELEASE_SALT(
  & (salt(i,j-1,k_s,bi,bj)+salt(i,j,k_cs,bi,bj))
  & +abs(pReleaseTransY(i,j,bi,bj)) *
  & (salt(i,j-1,k_s,bi,bj)-salt(i,j,k_cs,bi,bj)))
- & *recip_rA(i,j,bi,bj)
+ & *recip_rA(i,j,bi,bj)*_dxG(i,j,bi,bj)
+ & *recip_deepFac2C(k_cs)*recip_rhoFacC(k_cs)
+ & *recip_drF(k_cs)*_recip_hFacC(i,j,k_cs,bi,bj)
+ vol_trans_south = pReleaseTransY(i,j,bi,bj)
+ & *recip_rA(i,j,bi,bj)*_dxG(i,j,bi,bj)
  & *recip_deepFac2C(k_cs)*recip_rhoFacC(k_cs)
  & *recip_drF(k_cs)*_recip_hFacC(i,j,k_cs,bi,bj)
 ! S_trans_south =pReleaseTransY(i,j,bi,bj)*
@@ -136,7 +155,11 @@ SUBROUTINE PRESSURE_RELEASE_SALT(
  & (salt(i,j,k_cn,bi,bj)+salt(i,j+1,k_n,bi,bj))
  & +abs(pReleaseTransY(i,j+1,bi,bj)) *
  & (salt(i,j,k_cn,bi,bj)-salt(i,j+1,k_n,bi,bj)))
- & *recip_rA(i,j,bi,bj)
+ & *recip_rA(i,j,bi,bj)*_dxG(i,j+1,bi,bj)
+ & *recip_deepFac2C(k_cn)*recip_rhoFacC(k_cn)
+ & *recip_drF(k_cn)*_recip_hFacC(i,j,k_cn,bi,bj)
+ vol_trans_north = pReleaseTransY(i,j+1,bi,bj)
+ & *recip_rA(i,j,bi,bj)*_dxG(i,j+1,bi,bj)
  & *recip_deepFac2C(k_cn)*recip_rhoFacC(k_cn)
  & *recip_drF(k_cn)*_recip_hFacC(i,j,k_cn,bi,bj)
 ! S_trans_north =pReleaseTransY(i,j+1,bi,bj)*
@@ -147,16 +170,29 @@ SUBROUTINE PRESSURE_RELEASE_SALT(
  ENDIF
 
 
-
 C Add to get total tracer tendency.
  gS_arr(i,j) = gS_arr(i,j) + S_trans_west - S_trans_east
- & + S_trans_south - S_trans_north 
-
-
-
+ & + S_trans_south - S_trans_north 
+ if (k.eq.k_cn) then
+ gS_arr(i,j) = gS_arr(i,j) 
+ & - salt(i,j,k_cn,bi,bj) * (-vol_trans_north)
+ endif
+ if (k.eq.k_cs) then
+ gS_arr(i,j) = gS_arr(i,j) 
+ & - salt(i,j,k_cs,bi,bj) * (vol_trans_south)
+ endif
+ if (k.eq.k_ce) then
+ gS_arr(i,j) = gS_arr(i,j) 
+ & - salt(i,j,k_ce,bi,bj) * (-vol_trans_east)
+ endif
+ if (k.eq.k_cw) then
+ gS_arr(i,j) = gS_arr(i,j) 
+ & - salt(i,j,k_cw,bi,bj) * (vol_trans_west)
+ endif
  ENDDO
  ENDDO
 
+
 #endif /* ALLOW_PRESSURE_RELEASE_CODE */
 
  RETURN

diff --git a/model/src/pressure_release_theta.F b/model/src/pressure_release_theta.F
@@ -37,6 +37,8 @@ SUBROUTINE PRESSURE_RELEASE_THETA(
 C === Local Variables ===
  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
+ _RL vol_trans_west,vol_trans_east,
+ & vol_trans_south,vol_trans_north
 
  DO j=jMin+1,jMax-1
  DO i=iMin+1,iMax-1
@@ -45,42 +47,47 @@ SUBROUTINE PRESSURE_RELEASE_THETA(
  T_trans_north = 0.0
  T_trans_east = 0.0 
  T_trans_south = 0.0
+ vol_trans_west = 0.0 
+ vol_trans_north = 0.0
+ vol_trans_east = 0.0 
+ vol_trans_south = 0.0
 
 C calculate the k cells the tracers are transferred between in north,
 C south east and west cells.
 C Need to find if adjacent cells are deeper or shallower
- IF (kLowC(i,j,bi,bj) .GE. kLowC(i+1,j,bi,bj)) THEN
+ IF (kLowC(i,j,bi,bj) .GT. ksurfC(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
+ IF (kLowC(i,j,bi,bj) .GT. ksurfC(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
+ IF (kLowC(i,j,bi,bj) .GT. ksurfC(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
+ IF (kLowC(i,j,bi,bj) .GT. ksurfC(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.
 
@@ -92,7 +99,11 @@ SUBROUTINE PRESSURE_RELEASE_THETA(
  & (theta(i-1,j,k_w,bi,bj)+theta(i,j,k_cw,bi,bj))
  & +abs(pReleaseTransX(i,j,bi,bj)) * 
  & (theta(i-1,j,k_w,bi,bj)-theta(i,j,k_cw,bi,bj)))
- & *recip_rA(i,j,bi,bj)
+ & *recip_rA(i,j,bi,bj)*_dyG(i,j,bi,bj)
+ & *recip_deepFac2C(k_cw)*recip_rhoFacC(k_cw)
+ & *recip_drF(k_cw)*_recip_hFacC(i,j,k_cw,bi,bj)
+ vol_trans_west = pReleaseTransX(i,j,bi,bj)
+ & *recip_rA(i,j,bi,bj)*_dyG(i,j,bi,bj)
  & *recip_deepFac2C(k_cw)*recip_rhoFacC(k_cw)
  & *recip_drF(k_cw)*_recip_hFacC(i,j,k_cw,bi,bj)
  endif
@@ -104,7 +115,11 @@ SUBROUTINE PRESSURE_RELEASE_THETA(
  & (theta(i,j,k_ce,bi,bj)+theta(i+1,j,k_e,bi,bj))
  & +abs(pReleaseTransX(i+1,j,bi,bj)) *
  & (theta(i,j,k_ce,bi,bj)-theta(i+1,j,k_e,bi,bj)))
- & *recip_rA(i,j,bi,bj)
+ & *recip_rA(i,j,bi,bj)*_dyG(i+1,j,bi,bj)
+ & *recip_deepFac2C(k_ce)*recip_rhoFacC(k_ce)
+ & *recip_drF(k_ce)*_recip_hFacC(i,j,k_ce,bi,bj)
+ vol_trans_east = pReleaseTransX(i+1,j,bi,bj)
+ & *recip_rA(i,j,bi,bj)*_dyG(i+1,j,bi,bj)
  & *recip_deepFac2C(k_ce)*recip_rhoFacC(k_ce)
  & *recip_drF(k_ce)*_recip_hFacC(i,j,k_ce,bi,bj)
  endif
@@ -116,7 +131,11 @@ SUBROUTINE PRESSURE_RELEASE_THETA(
  & (theta(i,j-1,k_s,bi,bj)+theta(i,j,k_cs,bi,bj))
  & +abs(pReleaseTransY(i,j,bi,bj)) *
  & (theta(i,j-1,k_s,bi,bj)-theta(i,j,k_cs,bi,bj)))
- & *recip_rA(i,j,bi,bj)
+ & *recip_rA(i,j,bi,bj)*_dxG(i,j,bi,bj)
+ & *recip_deepFac2C(k_cs)*recip_rhoFacC(k_cs)
+ & *recip_drF(k_cs)*_recip_hFacC(i,j,k_cs,bi,bj)
+ vol_trans_south = pReleaseTransY(i,j,bi,bj)
+ & *recip_rA(i,j,bi,bj)*_dxG(i,j,bi,bj)
  & *recip_deepFac2C(k_cs)*recip_rhoFacC(k_cs)
  & *recip_drF(k_cs)*_recip_hFacC(i,j,k_cs,bi,bj)
  endif
@@ -128,7 +147,11 @@ SUBROUTINE PRESSURE_RELEASE_THETA(
  & (theta(i,j,k_cn,bi,bj)+theta(i,j+1,k_n,bi,bj))
  & +abs(pReleaseTransY(i,j+1,bi,bj)) *
  & (theta(i,j,k_cn,bi,bj)-theta(i,j+1,k_n,bi,bj)))
- & *recip_rA(i,j,bi,bj)
+ & *recip_rA(i,j,bi,bj)*_dxG(i,j+1,bi,bj)
+ & *recip_deepFac2C(k_cn)*recip_rhoFacC(k_cn)
+ & *recip_drF(k_cn)*_recip_hFacC(i,j,k_cn,bi,bj)
+ vol_trans_north = pReleaseTransY(i,j+1,bi,bj)
+ & *recip_rA(i,j,bi,bj)*_dxG(i,j+1,bi,bj)
  & *recip_deepFac2C(k_cn)*recip_rhoFacC(k_cn)
  & *recip_drF(k_cn)*_recip_hFacC(i,j,k_cn,bi,bj)
  endif
@@ -173,7 +196,22 @@ SUBROUTINE PRESSURE_RELEASE_THETA(
 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 
-
+ if (k.eq.k_cn) then
+ gT_arr(i,j) = gT_arr(i,j) 
+ & - theta(i,j,k_cn,bi,bj) * (-vol_trans_north)
+ endif
+ if (k.eq.k_cs) then
+ gT_arr(i,j) = gT_arr(i,j) 
+ & - theta(i,j,k_cs,bi,bj) * (vol_trans_south)
+ endif
+ if (k.eq.k_ce) then
+ gT_arr(i,j) = gT_arr(i,j) 
+ & - theta(i,j,k_ce,bi,bj) * (-vol_trans_east)
+ endif
+ if (k.eq.k_cw) then
+ gT_arr(i,j) = gT_arr(i,j) 
+ & - theta(i,j,k_cw,bi,bj) * (vol_trans_west)
+ endif
  ENDDO
  ENDDO