Line data Source code
1 : !-----------------------------------------------------------------------
2 : !BOP
3 : ! !ROUTINE: benergy --- Calculate the total energy (based on GFDL)
4 : !
5 : ! !INTERFACE:
6 :
7 0 : subroutine benergy(grid, u, v, t3, delp, &
8 0 : qqq, pe, peln, phis, &
9 0 : r_vir, cp3v, r3v, tte, te0 )
10 :
11 : ! !USES:
12 :
13 : use shr_kind_mod, only: r8 => shr_kind_r8
14 : use dynamics_vars, only: T_FVDYCORE_GRID
15 : use cam_logfile, only: iulog
16 : use par_vecsum_mod,only: par_vecsum
17 :
18 : #if defined( SPMD )
19 : use mod_comm, only: mp_send3d, mp_recv3d
20 : #endif
21 : implicit none
22 :
23 : ! !INPUT PARAMETERS:
24 : type (T_FVDYCORE_GRID), intent(in) :: grid ! YZ decomposition
25 :
26 : ! U-winds
27 : real(r8), intent(in) :: u(grid%ifirstxy:grid%ilastxy, &
28 : grid%jfirstxy:grid%jlastxy, &
29 : grid%km)
30 : ! V-winds
31 : real(r8), intent(in) :: v(grid%ifirstxy:grid%ilastxy, &
32 : grid%jfirstxy:grid%jlastxy, &
33 : grid%km)
34 :
35 : ! Temperature (K)
36 : real(r8), intent(in) :: t3(grid%ifirstxy:grid%ilastxy, &
37 : grid%jfirstxy:grid%jlastxy, &
38 : grid%km)
39 :
40 : ! Delta pressure
41 : real(r8), intent(in) :: delp(grid%ifirstxy:grid%ilastxy, &
42 : grid%jfirstxy:grid%jlastxy, &
43 : grid%km)
44 :
45 : ! Specific humidity
46 : real(r8), intent(in) :: qqq(grid%ifirstxy:grid%ilastxy, &
47 : grid%jfirstxy:grid%jlastxy, &
48 : grid%km)
49 :
50 : ! Edge pressure
51 : real(r8), intent(in) :: pe(grid%ifirstxy:grid%ilastxy, &
52 : grid%km+1, &
53 : grid%jfirstxy:grid%jlastxy)
54 :
55 : ! Edge pressure
56 : real(r8), intent(in) :: peln(grid%ifirstxy:grid%ilastxy, &
57 : grid%km+1, &
58 : grid%jfirstxy:grid%jlastxy)
59 :
60 : ! Surface heights
61 : real(r8), intent(in) :: phis(grid%ifirstxy:grid%ilastxy, &
62 : grid%jfirstxy:grid%jlastxy)
63 :
64 : real(r8), intent(in) :: r_vir ! Virtual effect constant ( rwv/rg-1 )
65 :
66 : ! C_p
67 : real(r8), intent(in) :: cp3v(grid%ifirstxy:grid%ilastxy, &
68 : grid%jfirstxy:grid%jlastxy, &
69 : grid%km)
70 : ! R (gas "constant")
71 : real(r8), intent(in) :: r3v(grid%ifirstxy:grid%ilastxy, &
72 : grid%jfirstxy:grid%jlastxy, &
73 : grid%km)
74 :
75 : ! !OUTPUT PARAMETERS:
76 :
77 : ! column integrated Total Energy
78 : real(r8), intent(out) :: tte(grid%jm)
79 : ! globally integrated total energy
80 : real(r8), intent(out) :: te0
81 :
82 : ! !DESCRIPTION:
83 : ! Determines the column and globally integrated total energy
84 : !
85 : ! !REVISION HISTORY:
86 : !
87 : ! SJL 99.04.13 : Delivered as release 0.9.8
88 : ! WS 99.05.18 : Added im, jm, km, te, dz as arguments
89 : ! WS 99.05.25 : Replaced IMR by IM, JMR by JM-1; removed fvcore.h
90 : ! WS 99.10.11 : Ghosted U, now fully limited to jfirst:jlast
91 : ! WS 99.11.23 : Pruned te, additional cleaning
92 : ! WS 00.05.14 : Renamed ghost indices as per Kevin's definitions
93 : ! WS 00.07.13 : Changed PILGRIM API
94 : ! WS 00.08.28 : Cosmetic changes
95 : ! AAM 00.08.28 : Added kfirst,klast
96 : ! WS 00.12.01 : Replaced MPI_ON with SPMD; hs now distributed
97 : ! AAM 01.06.15 : Changes for zero diff
98 : ! WS 01.12.10 : Ghosted PT
99 : ! WS 01.12.31 : Ghosted U,V
100 : ! WS 02.07.04 : Fixed 2D decomposition bug dest/src for mp_send3d
101 : ! WS 03.10.22 : pmgrid removed (now spmd_dyn)
102 : ! WS 03.12.03 : added grid as input argument
103 : ! WS 04.10.07 : Removed dependency on spmd_dyn; info now in GRID
104 : ! WS 06.05.02 : Rewritten for XY decomposition based on GFDL-code
105 : ! WS 06.06.21 : Extensive debugging of revised version
106 : !
107 : !EOP
108 : !---------------------------------------------------------------------
109 : !BOC
110 :
111 : ! Local
112 : real (r8), parameter :: D0_0 = 0.0_r8
113 : real (r8), parameter :: D0_25 = 0.25_r8
114 : real (r8), parameter :: D0_5 = 0.5_r8
115 : real (r8), parameter :: D1_0 = 1.0_r8
116 :
117 : integer :: im, jm, km, ifirstxy, ilastxy, jfirstxy, jlastxy
118 : integer :: iam, myidxy_x, myidxy_y, nprxy_x, nprxy_y, dest, src ! SPMD related
119 : integer :: i, j, k, js1g0, js2g0, jn1g0, jn1g1, jn2g0, ktot, jtot, itot
120 :
121 0 : real (r8) :: u2(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy+1)
122 0 : real (r8) :: v2(grid%ifirstxy:grid%ilastxy+1,grid%jfirstxy:grid%jlastxy)
123 :
124 0 : real (r8) :: tm(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy)
125 0 : real (r8) :: bte(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy)
126 0 : real (r8) :: te_sp(grid%ifirstxy:grid%ilastxy,grid%km)
127 0 : real (r8) :: te_np(grid%ifirstxy:grid%ilastxy,grid%km)
128 0 : real (r8) :: gztop(grid%ifirstxy:grid%ilastxy)
129 0 : real (r8) :: xsum(grid%jfirstxy:grid%jlastxy)
130 0 : real (r8) :: sp_sum(grid%km), np_sum(grid%km)
131 0 : real (r8) :: tm_sp(grid%km), tm_np(grid%km)
132 : real (r8) :: tmp
133 :
134 : real (r8) :: te(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy, &
135 0 : grid%km)
136 : real (r8) :: dz(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy, &
137 0 : grid%km)
138 0 : real(r8) :: veast(grid%jfirstxy:grid%jlastxy,grid%km) ! East halo
139 0 : real(r8) :: unorth(grid%ifirstxy:grid%ilastxy,grid%km) ! North halo
140 :
141 0 : im = grid%im
142 0 : jm = grid%jm
143 0 : km = grid%km
144 :
145 0 : ifirstxy = grid%ifirstxy
146 0 : ilastxy = grid%ilastxy
147 0 : jfirstxy = grid%jfirstxy
148 0 : jlastxy = grid%jlastxy
149 :
150 0 : iam = grid%iam
151 0 : myidxy_x = grid%myidxy_x
152 0 : myidxy_y = grid%myidxy_y
153 0 : nprxy_x = grid%nprxy_x
154 0 : nprxy_y = grid%nprxy_y
155 :
156 0 : js1g0 = max(1,jfirstxy)
157 0 : js2g0 = max(2,jfirstxy)
158 0 : jn2g0 = min(jm-1,jlastxy)
159 0 : jn1g0 = min(jm,jlastxy)
160 0 : jn1g1 = min(jm,jlastxy+1)
161 :
162 0 : itot = ilastxy - ifirstxy + 1
163 0 : jtot = jlastxy - jfirstxy + 1
164 :
165 : #if defined(SPMD)
166 : call mp_send3d( grid%commxy, iam-nprxy_x, iam+nprxy_x, im, jm, km, &
167 : ifirstxy, ilastxy, jfirstxy, jlastxy, 1, km, &
168 0 : ifirstxy, ilastxy, jfirstxy, jfirstxy, 1, km, u )
169 : call mp_recv3d( grid%commxy, iam+nprxy_x, im, jm, km, &
170 : ifirstxy, ilastxy, jlastxy+1, jlastxy+1, 1, km, &
171 0 : ifirstxy, ilastxy, jlastxy+1, jlastxy+1, 1, km, unorth )
172 :
173 0 : if (itot .ne. im) then
174 0 : dest = myidxy_y*nprxy_x + MOD(iam+nprxy_x-1,nprxy_x)
175 0 : src = myidxy_y*nprxy_x + MOD(iam+1,nprxy_x)
176 : call mp_send3d( grid%commxy, dest, src, im, jm, km, &
177 : ifirstxy, ilastxy, jfirstxy, jlastxy, 1, km, &
178 0 : ifirstxy, ifirstxy, jfirstxy, jlastxy, 1, km, v )
179 : call mp_recv3d( grid%commxy, src, im, jm, km, &
180 : ilastxy+1, ilastxy+1, jfirstxy, jlastxy, 1, km, &
181 0 : ilastxy+1, ilastxy+1, jfirstxy, jlastxy, 1, km, veast )
182 : else
183 : !$omp parallel do private(j, k)
184 0 : do k = 1,km
185 0 : do j=jfirstxy,jlastxy
186 0 : veast(j,k) = v(1,j,k)
187 : enddo
188 : enddo
189 : endif
190 : #else
191 : !$omp parallel do private(j, k)
192 : do k = 1,km
193 : do j=1,jm
194 : veast(j,k) = v(1,j,k)
195 : enddo
196 : enddo
197 : #endif
198 :
199 :
200 : !-----------------------------------------------------------------------------------------------
201 :
202 :
203 : !$omp parallel do private(i, j, k, u2, v2, tm)
204 0 : do k=1,km
205 : !
206 : ! Check the poles for consistent values
207 :
208 0 : do j=js2g0,jlastxy
209 0 : do i=ifirstxy,ilastxy
210 0 : u2(i,j) = grid%cose(j) * u(i,j,k)**2
211 : enddo
212 : enddo
213 :
214 0 : if ( jlastxy /= jm ) then ! Pull information out of northern halo
215 0 : do i=ifirstxy,ilastxy
216 0 : u2(i,jlastxy+1) = grid%cose(jlastxy+1) * unorth(i,k)**2
217 : enddo
218 : endif
219 :
220 0 : do j=js2g0,jn2g0
221 0 : do i=ifirstxy,ilastxy
222 0 : v2(i,j) = v(i,j,k)**2
223 : enddo
224 0 : v2(ilastxy+1,j) = veast(j,k)**2 ! eastern halo
225 : enddo
226 :
227 0 : do j=js2g0,jn2g0
228 0 : do i=ifirstxy,ilastxy
229 0 : te(i,j,k) = D0_25*((u2(i,j)+u2(i,j+1))*grid%acosu(j) + &
230 0 : v2(i,j) + v2(i+1,j))
231 : enddo
232 : enddo
233 :
234 0 : do j=jfirstxy,jlastxy
235 0 : do i=ifirstxy,ilastxy
236 0 : tm(i,j) = t3(i,j,k)*(D1_0+r_vir*qqq(i,j,k))
237 : enddo
238 : enddo
239 :
240 0 : do j=js2g0,jn2g0
241 0 : do i=ifirstxy, ilastxy
242 0 : te(i,j,k) = delp(i,j,k) * ( te(i,j,k) + cp3v(i,j,k)*tm(i,j) )
243 : enddo
244 : enddo
245 :
246 0 : if ( jfirstxy == 1 ) then
247 0 : do i=ifirstxy,ilastxy
248 0 : te_sp(i,k) = D0_5*u2(i,2)/grid%cose(2)
249 : enddo
250 0 : tm_sp(k) = tm(ifirstxy,1) ! All tm(:,1) should be the same
251 : endif
252 :
253 0 : if ( jlastxy == jm ) then
254 0 : do i=ifirstxy,ilastxy
255 0 : te_np(i,k)= D0_5*u2(i,jm)/grid%cose(jm)
256 : enddo
257 0 : tm_np(k) = tm(ifirstxy,jm) ! All tm(:,jm) should be the same
258 : endif
259 :
260 0 : do j=jfirstxy,jlastxy
261 0 : do i=ifirstxy,ilastxy
262 0 : dz(i,j,k) = r3v(i,j,k)*tm(i,j)
263 : enddo
264 : enddo
265 : enddo
266 :
267 :
268 0 : if ( jfirstxy == 1 ) then
269 0 : call par_xsum( grid, te_sp, km, sp_sum )
270 : !$omp parallel do private(i, k, tmp)
271 0 : do k=1,km
272 0 : tmp = delp(ifirstxy,1,k) * (D0_5*sp_sum(k)/real(im,r8) + &
273 0 : cp3v(ifirstxy,1,k)*tm_sp(k))
274 0 : do i=ifirstxy,ilastxy
275 0 : te(i,1,k) = tmp
276 : enddo
277 : enddo
278 : endif
279 0 : if ( jlastxy == jm ) then
280 0 : call par_xsum( grid, te_np, km, np_sum )
281 : !$omp parallel do private(i, k, tmp)
282 0 : do k=1,km
283 0 : tmp = delp(ifirstxy,jm,k) * (D0_5*np_sum(k)/real(im,r8) +&
284 0 : cp3v(ifirstxy,jm,k)*tm_np(k))
285 0 : do i=ifirstxy,ilastxy
286 0 : te(i,jm,k) = tmp
287 : enddo
288 : enddo
289 : endif
290 :
291 0 : bte = D0_0
292 : !$omp parallel do private(i,j,k,gztop)
293 0 : do j=jfirstxy,jlastxy
294 : ! Perform vertical integration
295 0 : do i=ifirstxy,ilastxy
296 0 : gztop(i) = phis(i,j)
297 0 : do k=1,km
298 0 : gztop(i) = gztop(i) + dz(i,j,k)*(peln(i,k+1,j)-peln(i,k,j))
299 : enddo
300 : enddo
301 :
302 0 : if (j == 1) then
303 : ! gztop(:) should all have identical values WS 2006.06.22: this checks out
304 : ! SP
305 0 : tte(1) = pe(ifirstxy,km+1,1)*phis(ifirstxy,1) - pe(ifirstxy,1,1)*gztop(ifirstxy)
306 0 : do k=1,km
307 0 : tte(1) = tte(1) + te(ifirstxy,1,k)
308 : enddo
309 0 : tte(1) = grid%acap * tte(1)
310 0 : elseif (j == jm) then
311 : ! gztop(:) should all have identical values WS 2006.06.22: this checks out
312 : ! NP
313 0 : tte(jm) = pe(ifirstxy,km+1,jm)*phis(ifirstxy,jm) - pe(ifirstxy,1,jm)*gztop(ifirstxy)
314 0 : do k=1,km
315 0 : tte(jm) = tte(jm) + te(ifirstxy,jm,k)
316 : enddo
317 0 : tte(jm) = grid%acap * tte(jm)
318 : else
319 : ! Interior
320 :
321 0 : do i=ifirstxy,ilastxy
322 0 : bte(i,j) = pe(i,km+1,j)*phis(i,j) - pe(i,1,j)*gztop(i)
323 : enddo
324 :
325 0 : do k=1,km
326 0 : do i=ifirstxy,ilastxy
327 0 : bte(i,j) = bte(i,j) + te(i,j,k)
328 : enddo
329 : enddo
330 : endif
331 : enddo
332 :
333 0 : call par_xsum(grid, bte, jtot, xsum)
334 :
335 : !$omp parallel do private(j)
336 0 : do j=js2g0,jn2g0
337 0 : tte(j) = xsum(j)*grid%cosp(j)
338 : enddo
339 :
340 0 : call par_vecsum(jm, jfirstxy, jlastxy, tte, te0, grid%commxy_y, grid%nprxy_y)
341 :
342 0 : write(iulog,*) "myidxy_x/y:", myidxy_x, myidxy_y, "The total energy is", te0
343 :
344 : !EOC
345 0 : end subroutine benergy
346 : !-----------------------------------------------------------------------
|
