/[MITgcm]/MITgcm_contrib/ksnow/press_release/code/pressure_release_theta.F
ViewVC logotype

Contents of /MITgcm_contrib/ksnow/press_release/code/pressure_release_theta.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.4 - (show annotations) (download)
Fri Apr 7 13:49:05 2017 UTC (8 years, 8 months ago) by dgoldberg
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +5 -5 lines
correct k index for cell cross-transfer

1 C $Header: /u/gcmpack/MITgcm_contrib/ksnow/press_release/code/pressure_release_theta.F,v 1.3 2017/03/14 15:57:18 dgoldberg 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
37 CEndOfInterface
38
39 #ifdef ALLOW_PRESSURE_RELEASE_CODE
40 C === Local Variables ===
41 INTEGER i,j,k_e,k_ce,k_s,k_cs,k_w,k_cw,k_n,k_cn
42 _RL T_trans_west,T_trans_east,T_trans_south,T_trans_north
43
44 DO j=jMin+1,jMax-1
45 DO i=iMin+1,iMax-1
46
47 T_trans_west = 0.0
48 T_trans_north = 0.0
49 T_trans_east = 0.0
50 T_trans_south = 0.0
51
52 C calculate the k cells the tracers are transferred between in north,
53 C south east and west cells.
54 C Need to find if adjacent cells are deeper or shallower
55 IF (kLowC(i,j,bi,bj) .GE. kLowC(i+1,j,bi,bj)) THEN
56 k_e = kLowC(i+1,j,bi,bj)
57 k_ce = kSurfC(i,j,bi,bj)
58 ELSE
59 k_e = kSurfC(i+1,j,bi,bj)
60 k_ce = kLowC(i,j,bi,bj)
61 ENDIF
62
63 IF (kLowC(i,j,bi,bj) .GE. kLowC(i-1,j,bi,bj)) THEN
64 k_w = kLowC(i-1,j,bi,bj)
65 k_cw = kSurfC(i,j,bi,bj)
66 ELSE
67 k_w = kSurfC(i-1,j,bi,bj)
68 k_cw = kLowC(i,j,bi,bj)
69 ENDIF
70
71 IF (kLowC(i,j,bi,bj) .GE. kLowC(i,j+1,bi,bj)) THEN
72 k_n = kLowC(i,j+1,bi,bj)
73 k_cn = kSurfC(i,j,bi,bj)
74 ELSE
75 k_n = kSurfC(i,j+1,bi,bj)
76 k_cn = kLowC(i,j,bi,bj)
77 ENDIF
78
79 IF (kLowC(i,j,bi,bj) .GE. kLowC(i,j-1,bi,bj)) THEN
80 k_s = kLowC(i,j-1,bi,bj)
81 k_cs = kSurfC(i,j,bi,bj)
82 ELSE
83 k_s = kSurfC(i,j-1,bi,bj)
84 k_cs = kLowC(i,j,bi,bj)
85 ENDIF
86
87 C calculate the net tracer flux through north, south east and west
88 C faces.
89
90
91 IF (k .EQ. k_cw) THEN
92 if (k_cw.gt.0 .and. k_w.gt.0) then
93 T_trans_west =0.5 _d 0 *
94 & ( pReleaseTransX(i,j,bi,bj) *
95 & (theta(i-1,j,k_w,bi,bj)+theta(i,j,k_cw,bi,bj))
96 & +abs(pReleaseTransX(i,j,bi,bj)) *
97 & (theta(i-1,j,k_w,bi,bj)-theta(i,j,k_cw,bi,bj)))
98 & *recip_rA(i,j,bi,bj)
99 & *recip_deepFac2C(k_cw)*recip_rhoFacC(k_cw)
100 & *recip_drF(k_cw)*_recip_hFacC(i,j,k_cw,bi,bj)
101 endif
102 ENDIF
103 IF (k .EQ. k_ce) THEN
104 if (k_ce.gt.0 .and. k_e.gt.0) then
105 T_trans_east =0.5 _d 0 *
106 & ( pReleaseTransX(i+1,j,bi,bj) *
107 & (theta(i,j,k_ce,bi,bj)+theta(i+1,j,k_e,bi,bj))
108 & +abs(pReleaseTransX(i+1,j,bi,bj)) *
109 & (theta(i,j,k_ce,bi,bj)-theta(i+1,j,k_e,bi,bj)))
110 & *recip_rA(i,j,bi,bj)
111 & *recip_deepFac2C(k_ce)*recip_rhoFacC(k_ce)
112 & *recip_drF(k_ce)*_recip_hFacC(i,j,k_ce,bi,bj)
113 endif
114 ENDIF
115 IF (k .EQ. k_cs) THEN
116 if (k_cs.gt.0 .and. k_s.gt.0) then
117 T_trans_south =0.5 _d 0 *
118 & ( pReleaseTransY(i,j,bi,bj) *
119 & (theta(i,j-1,k_s,bi,bj)+theta(i,j,k_cs,bi,bj))
120 & +abs(pReleaseTransY(i,j,bi,bj)) *
121 & (theta(i,j-1,k_s,bi,bj)-theta(i,j,k_cs,bi,bj)))
122 & *recip_rA(i,j,bi,bj)
123 & *recip_deepFac2C(k_cs)*recip_rhoFacC(k_cs)
124 & *recip_drF(k_cs)*_recip_hFacC(i,j,k_cs,bi,bj)
125 endif
126 ENDIF
127 IF (k .EQ. k_cn) THEN
128 if (k_cn.gt.0 .and. k_n.gt.0) then
129 T_trans_north =0.5 _d 0 *
130 & ( pReleaseTransY(i,j+1,bi,bj) *
131 & (theta(i,j,k_cn,bi,bj)+theta(i,j+1,k_n,bi,bj))
132 & +abs(pReleaseTransY(i,j+1,bi,bj)) *
133 & (theta(i,j,k_cn,bi,bj)-theta(i,j+1,k_n,bi,bj)))
134 & *recip_rA(i,j,bi,bj)
135 & *recip_deepFac2C(k_cn)*recip_rhoFacC(k_cn)
136 & *recip_drF(k_cn)*_recip_hFacC(i,j,k_cn,bi,bj)
137 endif
138 ENDIF
139
140
141 ! IF (k .EQ. k_cw) THEN
142 ! if (k_w.gt.0 .and. k_cw.gt.0) then
143 ! T_trans_west =pReleaseTransX(i,j,bi,bj)*
144 ! & (theta(i-1,j,k_w,bi,bj) -theta(i,j,k_cw,bi,bj))
145 C & *rhoFacC(k)*mass2rUnit
146 C & *_dyG(i,j,bi,bj)*recip_rA(i,j,bi,bj)*
147 ! & *recip_dxG(i,j,bi,bj)
148 ! & *recip_drF(k_cw)*_recip_hFacC(i,j,k_cw,bi,bj)
149 ! endif
150 ! ENDIF
151 ! IF (k .EQ. k_ce) THEN
152 ! if (k_e.gt.0 .and. k_ce.gt.0) then
153 ! T_trans_east =pReleaseTransX(i+1,j,bi,bj)*
154 ! & (theta(i,j,k_ce,bi,bj) -theta(i+1,j,k_e,bi,bj))
155 ! & *recip_dxG(i,j,bi,bj)
156 ! & *recip_drF(k_ce)*_recip_hFacC(i,j,k_ce,bi,bj)
157 ! endif
158 ! ENDIF
159 ! IF (k .EQ. k_cs) THEN
160 ! if (k_s.gt.0 .and. k_cs.gt.0) then
161 ! T_trans_south =pReleaseTransY(i,j,bi,bj)*
162 ! & (theta(i,j-1,k_s,bi,bj) -theta(i,j,k_cs,bi,bj))
163 ! & *recip_dyG(i,j,bi,bj)
164 ! & *recip_drF(k_cs)*_recip_hFacC(i,j,k_cs,bi,bj)
165 ! endif
166 ! ENDIF
167 ! IF (k .EQ. k_cn) THEN
168 ! if (k_n.gt.0 .and. k_cn.gt.0) then
169 ! T_trans_north =pReleaseTransY(i,j+1,bi,bj)*
170 ! & (theta(i,j,k_cn,bi,bj) -theta(i,j+1,k_n,bi,bj))
171 ! & *recip_dyG(i,j,bi,bj)
172 ! & *recip_drF(k_cn)*_recip_hFacC(i,j,k_cn,bi,bj)
173 ! endif
174 ! ENDIF
175
176 C Add to get total tracer tendency.
177 gT_arr(i,j) = gT_arr(i,j) + T_trans_west - T_trans_east
178 & + T_trans_south - T_trans_north
179
180 ENDDO
181 ENDDO
182
183 #endif /* ALLOW_PRESSURE_RELEASE_CODE */
184
185 RETURN
186 END

  ViewVC Help
Powered by ViewVC 1.1.22