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

Annotation 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 - (hide 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 dgoldberg 1.4 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 ksnow 1.1 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 ksnow 1.2 #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 ksnow 1.1
44     DO j=jMin+1,jMax-1
45     DO i=iMin+1,iMax-1
46    
47 ksnow 1.2 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 ksnow 1.1 C calculate the k cells the tracers are transferred between in north,
53     C south east and west cells.
54 ksnow 1.2 C Need to find if adjacent cells are deeper or shallower
55 ksnow 1.1 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 ksnow 1.2
90 dgoldberg 1.3
91 ksnow 1.2 IF (k .EQ. k_cw) THEN
92 dgoldberg 1.3 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 ksnow 1.1 & *recip_drF(k_cw)*_recip_hFacC(i,j,k_cw,bi,bj)
101 dgoldberg 1.3 endif
102     ENDIF
103 ksnow 1.2 IF (k .EQ. k_ce) THEN
104 dgoldberg 1.3 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 dgoldberg 1.4 & (theta(i,j,k_ce,bi,bj)+theta(i+1,j,k_e,bi,bj))
108 dgoldberg 1.3 & +abs(pReleaseTransX(i+1,j,bi,bj)) *
109 dgoldberg 1.4 & (theta(i,j,k_ce,bi,bj)-theta(i+1,j,k_e,bi,bj)))
110 dgoldberg 1.3 & *recip_rA(i,j,bi,bj)
111     & *recip_deepFac2C(k_ce)*recip_rhoFacC(k_ce)
112 ksnow 1.2 & *recip_drF(k_ce)*_recip_hFacC(i,j,k_ce,bi,bj)
113 dgoldberg 1.3 endif
114 ksnow 1.2 ENDIF
115     IF (k .EQ. k_cs) THEN
116 dgoldberg 1.3 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 ksnow 1.1 & *recip_drF(k_cs)*_recip_hFacC(i,j,k_cs,bi,bj)
125 dgoldberg 1.3 endif
126 ksnow 1.2 ENDIF
127     IF (k .EQ. k_cn) THEN
128 dgoldberg 1.3 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 dgoldberg 1.4 & (theta(i,j,k_cn,bi,bj)+theta(i,j+1,k_n,bi,bj))
132 dgoldberg 1.3 & +abs(pReleaseTransY(i,j+1,bi,bj)) *
133 dgoldberg 1.4 & (theta(i,j,k_cn,bi,bj)-theta(i,j+1,k_n,bi,bj)))
134 dgoldberg 1.3 & *recip_rA(i,j,bi,bj)
135     & *recip_deepFac2C(k_cn)*recip_rhoFacC(k_cn)
136 ksnow 1.2 & *recip_drF(k_cn)*_recip_hFacC(i,j,k_cn,bi,bj)
137 dgoldberg 1.3 endif
138 ksnow 1.2 ENDIF
139 ksnow 1.1
140 dgoldberg 1.3
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 ksnow 1.1 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 ksnow 1.2 #endif /* ALLOW_PRESSURE_RELEASE_CODE */
184    
185 ksnow 1.1 RETURN
186     END

  ViewVC Help
Powered by ViewVC 1.1.22