Line data Source code
1 : module te_map_mod
2 :
3 : implicit none
4 : save
5 : private
6 : public :: te_map
7 :
8 : contains
9 : !-----------------------------------------------------------------------
10 : !BOP
11 : ! !ROUTINE: te_map --- Map vertical Lagrangian coordinates to normal grid
12 : !
13 : ! !INTERFACE:
14 :
15 16128 : subroutine te_map(grid, consv, convt, ps, omga, &
16 16128 : pe, delp, pkz, pk, mdt, &
17 16128 : nx, u, v, pt, tracer, &
18 16128 : hs, cp3v, cap3v, kord, peln, &
19 16128 : te0, te, dz, mfx, mfy, &
20 16128 : uc, vc, du_s, du_w, &
21 : am_geom_crrct, am_diag_lbl)
22 : !
23 : ! !USES:
24 :
25 : use shr_kind_mod, only : r8 => shr_kind_r8
26 : use spmd_utils, only : masterproc
27 : use dynamics_vars, only : T_FVDYCORE_GRID
28 : use mapz_module, only : map1_cubic_te, map1_ppm, mapn_ppm_tracer
29 : use cam_logfile, only : iulog
30 : use cam_abortutils,only : endrun
31 :
32 : #if defined( SPMD )
33 : use mod_comm, only : mp_send3d, mp_recv3d
34 : #endif
35 : use phys_control, only: waccmx_is !WACCM-X runtime switch
36 : use cam_thermo, only: cam_thermo_calc_kappav
37 : use par_vecsum_mod,only: par_vecsum
38 :
39 : implicit none
40 :
41 : #if defined( SPMD )
42 : #define CPP_PRT_PREFIX if(grid%iam==0)
43 : #else
44 : #define CPP_PRT_PREFIX
45 : #endif
46 :
47 : ! !INPUT PARAMETERS:
48 : type (T_FVDYCORE_GRID), intent(inout) :: grid ! grid for XY decomp
49 : logical consv ! flag to force TE conservation
50 : logical convt ! flag to control pt output (see below)
51 : integer mdt ! mapping time step (same as phys)
52 : integer nx ! number of SMP "decomposition" in x
53 : real(r8) hs(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy) ! surface geopotential
54 : real(r8) te0
55 :
56 : ! !INPUT/OUTPUT PARAMETERS:
57 : real(r8) pk(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km+1) ! pe to the kappa
58 : real(r8) u(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km) ! u-wind (m/s)
59 : real(r8) v(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km) ! v-wind (m/s)
60 : ! tracers including specific humidity
61 : !!! real(r8) q(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km,grid%ntotq)
62 :
63 : real(r8), intent(inout) :: &
64 : cp3v(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km) ! C_p
65 : real(r8), intent(inout) :: &
66 : cap3v(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km) ! cappa
67 :
68 : real(r8), intent(inout) :: &
69 : tracer(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km,grid%ntotq) ! Tracer array
70 : real(r8) pe(grid%ifirstxy:grid%ilastxy,grid%km+1,grid%jfirstxy:grid%jlastxy) ! pressure at layer edges
71 : real(r8) ps(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy) ! surface pressure
72 : real(r8) pt(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km) ! virtual potential temperature as input
73 : ! Output: virtual temperature if convt is true
74 : ! false: output is (virtual) potential temperature
75 : real(r8) te(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km) ! Work array (cache performance)
76 : real(r8) dz(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km) ! Work array (cache performance)
77 : real(r8) mfx(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km)
78 : real(r8) mfy(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km)
79 :
80 : ! !OUTPUT PARAMETERS:
81 : real(r8) delp(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km) ! pressure thickness
82 : real(r8) omga(grid%ifirstxy:grid%ilastxy,grid%km,grid%jfirstxy:grid%jlastxy) ! vertical press. velocity (pascal/sec)
83 : real(r8) peln(grid%ifirstxy:grid%ilastxy,grid%km+1,grid%jfirstxy:grid%jlastxy) ! log(pe)
84 : real(r8) pkz(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km) ! layer-mean pk for converting t to pt
85 :
86 : ! AM conservation mods
87 : logical, intent(in) :: am_geom_crrct ! logical switch for AM correction
88 : logical, intent(in) :: am_diag_lbl ! input
89 :
90 : real(r8), intent(in) :: du_s(grid%km)
91 : real(r8), intent(inout), allocatable :: du_w(:,:,:)
92 :
93 : real(r8), intent(inout), allocatable :: uc(:,:,:)
94 : real(r8), intent(inout), allocatable :: vc(:,:,:)
95 :
96 : ! !DESCRIPTION:
97 : !
98 : ! !REVISION HISTORY:
99 : !
100 : ! WS 99.05.19 : Replaced IMR, JMR, JNP and NL with IM, JM-1, JM and KM
101 : ! WS 99.05.25 : Revised conversions with IMR and JMR; removed fvcore.h
102 : ! WS 99.07.29 : Reduced computation region to grid%jfirstxy:jlast
103 : ! WS 99.07.30 : Tested concept with message concatenation of te_map calls
104 : ! WS 99.10.01 : Documentation; indentation; cleaning
105 : ! SJL 99.12.31: SMP "decomposition" in-E-W direction
106 : ! WS 00.05.14 : Renamed ghost indices as per Kevin's definitions
107 : ! WS 00.07.13 : Changed PILGRIM API
108 : ! AM 00.08.29 : Variables in this routine will ultimately be decomposed in (i,j).
109 : ! AM 01.06.13 : 2-D decomposition; reordering summation causes roundoff difference.
110 : ! WS 01.06.10 : Removed "if(first)" section in favor of a variable module
111 : ! AM 01.06.27 : Merged yz decomposition technology into ccm code.
112 : ! WS 02.01.14 : Upgraded to mod_comm
113 : ! WS 02.04.22 : Use mapz_module from FVGCM
114 : ! WS 02.04.25 : New mod_comm interfaces
115 : ! WS 03.08.12 : Introduced unorth
116 : ! WS 03.11.19 : Merged in CAM changes by Mirin
117 : ! WS 03.12.03 : Added GRID as argument, dynamics_vars removed
118 : ! WS 04.08.25 : Simplified interface by using GRID
119 : ! WS 04.10.07 : Removed dependency on spmd_dyn; info now in GRID
120 : ! WS 05.03.25 : Changed tracer to type T_TRACERS
121 : ! WS 05.04.12 : Call mapn_ppm_tracer instead of mapn_ppm
122 : ! AT 05.05.11 : Merged with the version Cerebus (unique pole issues)
123 : ! WS 05.05.18 : Merged CAM and GEOS5 versions (mostly GEOS5)
124 : ! LT 05.11.14 : Call map1_cubic_te for Cubic Interpolation of Total Energy
125 : ! WP 06.01.18 : Added calls to map1_ppm for horizontal mass fluxes
126 : ! LT 06.02.08 : Implement code for partial remapping option
127 : ! WS 06.11.29 : Merge CAM/GEOS5; magic numbers isolated
128 : ! CC 07.01.29 : Additions for proper calculation of OMGA
129 : !
130 : !EOP
131 : !-----------------------------------------------------------------------
132 : !BOC
133 : ! !LOCAL VARIABLES:
134 :
135 : ! Magic numbers used in this module
136 : real(r8), parameter :: D0_0 = 0.0_r8
137 : real(r8), parameter :: D0_25 = 0.25_r8
138 : real(r8), parameter :: D0_5 = 0.5_r8
139 : real(r8), parameter :: D1_0 = 1.0_r8
140 : real(r8), parameter :: D2_0 = 2.0_r8
141 : real(r8), parameter :: D10_0 = 10.0_r8
142 : real(r8), parameter :: D1E4 = 1.0e4_r8
143 :
144 : integer :: im, jm, km ! x, y, z dimensions
145 : integer :: nq ! number of tracers to be advected
146 : integer :: ifirst, ilast ! starting & ending longitude index
147 : integer :: jfirst, jlast ! starting & ending latitude index
148 : integer :: myidxy_y, iam
149 : integer :: nprxy_x, nprxy_y
150 : integer :: kk
151 :
152 : ! Local variables for Partial Remapping
153 : ! -------------------------------------
154 32256 : real(r8) :: pref(grid%km+1)
155 32256 : real(r8) :: zz(grid%km+1)
156 : real(r8) :: z1,z2
157 :
158 : real(r8), parameter :: alf = 0.042_r8
159 : real(r8), parameter :: pa = 1.0_r8
160 : real(r8), parameter :: pb = 500.0_r8
161 : real(r8), parameter :: psurf = 100001.0_r8
162 : real(r8), parameter :: bet = D2_0*alf/(D1_0+alf)
163 :
164 : real(r8), parameter :: lagrangianlevcrit = 1.0e-11_r8 ! Criteria for "Lagrangian levels are crossing" error
165 :
166 : ! Local arrays:
167 : ! -------------
168 32256 : real(r8) rmin(nx*grid%jm), rmax(nx*grid%jm)
169 32256 : real(r8) tte(grid%jm)
170 : ! x-y
171 32256 : real(r8) u2(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy+1)
172 32256 : real(r8) v2(grid%ifirstxy:grid%ilastxy+1,grid%jfirstxy:grid%jlastxy)
173 32256 : real(r8) t2(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy)
174 32256 : real(r8) veast(grid%jfirstxy:grid%jlastxy,grid%km)
175 : ! y-z
176 32256 : real(r8) pe0(grid%ifirstxy:grid%ilastxy,grid%km+1)
177 32256 : real(r8) pe1(grid%ifirstxy:grid%ilastxy,grid%km+1)
178 32256 : real(r8) pe2(grid%ifirstxy:grid%ilastxy,grid%km+1)
179 32256 : real(r8) pe3(grid%ifirstxy:grid%ilastxy,grid%km+1)
180 32256 : real(r8) u2_sp(grid%ifirstxy:grid%ilastxy,grid%km)
181 32256 : real(r8) v2_sp(grid%ifirstxy:grid%ilastxy,grid%km)
182 32256 : real(r8) t2_sp(grid%ifirstxy:grid%ilastxy,grid%km)
183 32256 : real(r8) u2_np(grid%ifirstxy:grid%ilastxy,grid%km)
184 32256 : real(r8) v2_np(grid%ifirstxy:grid%ilastxy,grid%km)
185 32256 : real(r8) t2_np(grid%ifirstxy:grid%ilastxy,grid%km)
186 :
187 : ! Log of pe1/pe2.
188 32256 : real(r8) log_pe1(grid%ifirstxy:grid%ilastxy,grid%km+1)
189 32256 : real(r8) log_pe2(grid%ifirstxy:grid%ilastxy,grid%km+1)
190 :
191 : ! x
192 : real(r8) gz(grid%ifirstxy:grid%ilastxy)
193 32256 : real(r8) ratio(grid%ifirstxy:grid%ilastxy)
194 32256 : real(r8) bte(grid%ifirstxy:grid%ilastxy)
195 : ! z
196 32256 : real(r8) pe1w(grid%km+1)
197 32256 : real(r8) pe2w(grid%km+1)
198 :
199 : ! AM correction
200 : ! variable for zonal momentum
201 32256 : real(r8) :: dum(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy)
202 :
203 32256 : real(r8) cap3vi(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km+1) ! cappa interface
204 :
205 : integer i1w, nxu
206 : integer i, j, k, js2g0, jn2g0, jn1g1
207 : integer kord
208 : integer krd
209 :
210 : real(r8) dak, bkh, qmax, qmin
211 32256 : real(r8) te_sp(grid%km), te_np(grid%km)
212 32256 : real(r8) xysum(grid%jfirstxy:grid%jlastxy,2)
213 32256 : real(r8) tmpik(grid%ifirstxy:grid%ilastxy,grid%km)
214 32256 : real(r8) tmpij(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,2)
215 32256 : real(r8) omga_ik(grid%ifirstxy:grid%ilastxy,grid%km) ! vertical press. velocity (tmp 2-d array)
216 : real(r8) dtmp
217 : real(r8) sum
218 : real(r8) te1
219 : real(r8) dlnp
220 :
221 : integer ixj, jp, it, i1, i2
222 :
223 : #if defined( SPMD )
224 : integer :: dest, src
225 32256 : real(r8) unorth(grid%ifirstxy:grid%ilastxy,grid%km)
226 32256 : real(r8) pewest(grid%km+1,grid%jfirstxy:grid%jlastxy)
227 16128 : real(r8), allocatable :: pesouth(:,:)
228 : #endif
229 :
230 : integer comm_use, npry_use, itot
231 :
232 : logical diag
233 : logical :: high_alt
234 :
235 : data diag /.false./
236 :
237 16128 : z1 = log(pa/psurf)
238 16128 : z2 = log(pb/psurf)
239 :
240 16128 : high_alt = grid%high_alt
241 :
242 16128 : im = grid%im
243 16128 : jm = grid%jm
244 16128 : km = grid%km
245 16128 : nq = grid%nq
246 :
247 16128 : ifirst = grid%ifirstxy
248 16128 : ilast = grid%ilastxy
249 16128 : jfirst = grid%jfirstxy
250 16128 : jlast = grid%jlastxy
251 :
252 16128 : iam = grid%iam
253 16128 : myidxy_y = grid%myidxy_y
254 16128 : nprxy_x = grid%nprxy_x
255 16128 : nprxy_y = grid%nprxy_y
256 :
257 : ! Intialize PREF for Partial Remapping (above 100-mb)
258 : ! -----------------------------------------------------------
259 2128896 : do k=1,km+1
260 2128896 : pref(k) = grid%ak(k) + grid%bk(k)*psurf
261 : enddo
262 2128896 : zz = log( pref/psurf )
263 2128896 : zz = D10_0*(zz-z2)/z1
264 2128896 : zz = (D1_0-bet)*tanh(zz)
265 :
266 : ! WS 99.07.27 : Set loop limits appropriately
267 : ! --------------------------------------------
268 16128 : js2g0 = max(2,jfirst)
269 16128 : jn1g1 = min(jm,jlast+1)
270 16128 : jn2g0 = min(jm-1,jlast)
271 64512 : do j=jfirst,jlast
272 48384 : xysum(j,1) = D0_0
273 64512 : xysum(j,2) = D0_0
274 : enddo
275 64512 : do j=jfirst,jlast
276 645120 : do i=ifirst,ilast
277 580608 : tmpij(i,j,1) = D0_0
278 628992 : tmpij(i,j,2) = D0_0
279 : enddo
280 : enddo
281 2112768 : do k=1,km
282 27272448 : do i=ifirst,ilast
283 27256320 : tmpik(i,k) = D0_0
284 : enddo
285 : enddo
286 :
287 16128 : itot = ilast-ifirst+1
288 16128 : nxu = 1
289 16128 : if (itot == im) nxu = nx
290 :
291 : #if defined( SPMD )
292 16128 : comm_use = grid%commxy_y
293 16128 : npry_use = nprxy_y
294 :
295 : call mp_send3d( grid%commxy, iam-nprxy_x, iam+nprxy_x, im, jm, km, &
296 : ifirst, ilast, jfirst, jlast, 1, km, &
297 16128 : ifirst, ilast, jfirst, jfirst, 1, km, u )
298 : ! Nontrivial x decomposition
299 16128 : if (itot /= im) then
300 16128 : dest = myidxy_y*nprxy_x + MOD(iam+nprxy_x-1,nprxy_x)
301 16128 : src = myidxy_y*nprxy_x + MOD(iam+1,nprxy_x)
302 : call mp_send3d( grid%commxy, dest, src, im, jm, km, &
303 : ifirst, ilast, jfirst, jlast, 1,km, &
304 16128 : ifirst, ifirst, jfirst, jlast, 1, km, v )
305 : endif
306 : #endif
307 : call pkez(nxu, im, km, jfirst, jlast, 1, km, ifirst, ilast, &
308 16128 : pe, pk, cap3v, grid%ks, peln, pkz, .false., high_alt)
309 :
310 : ! Single subdomain case (periodic)
311 2112768 : do k=1,km
312 8402688 : do j=jfirst,jlast
313 8386560 : veast(j,k) = v(ifirst,j,k)
314 : enddo
315 : enddo
316 : #if defined( SPMD )
317 : call mp_recv3d( grid%commxy, iam+nprxy_x, im, jm, km, &
318 : ifirst, ilast, jlast+1, jlast+1, 1, km, &
319 16128 : ifirst, ilast, jlast+1, jlast+1, 1, km, unorth )
320 : ! Nontrivial x decomposition
321 16128 : if (itot /= im) then
322 : call mp_recv3d( grid%commxy, src, im, jm, km, &
323 : ilast+1, ilast+1, jfirst, jlast, 1, km, &
324 16128 : ilast+1, ilast+1, jfirst, jlast, 1, km, veast )
325 16128 : dest = myidxy_y*nprxy_x + MOD(iam+1,nprxy_x)
326 16128 : src = myidxy_y*nprxy_x + MOD(iam+nprxy_x-1,nprxy_x)
327 : call mp_send3d( grid%commxy, dest, src, im, km+1, jm, &
328 : ifirst, ilast, 1, km+1, jfirst, jlast, &
329 16128 : ilast, ilast, 1, km+1, jfirst, jlast, pe )
330 : endif
331 : call mp_send3d( grid%commxy, iam+nprxy_x, iam-nprxy_x, im, km+1,jm, &
332 : ifirst, ilast, 1, km+1, jfirst, jlast, &
333 16128 : ifirst, ilast, 1, km+1, jlast, jlast, pe )
334 : #endif
335 :
336 16128 : if (high_alt) then
337 16128 : call cam_thermo_calc_kappav( tracer, cap3v, cpv=cp3v )
338 : endif
339 :
340 : !$omp parallel do &
341 : !$omp default(shared) &
342 : !$omp private(i,j, k, u2, v2, t2)
343 :
344 : ! Compute cp*T + KE
345 :
346 2112768 : do 1000 k=1,km
347 :
348 8321040 : do j=js2g0,jlast
349 83013840 : do i=ifirst,ilast
350 80917200 : u2(i,j) = u(i,j,k)**2
351 : enddo
352 : enddo
353 : #if defined( SPMD )
354 2096640 : if ( jlast < jm ) then
355 26404560 : do i=ifirst,ilast
356 26404560 : u2(i,jlast+1) = unorth(i,k)**2 ! fill ghost zone
357 : enddo
358 : endif
359 : #endif
360 :
361 8255520 : do j=js2g0,jn2g0
362 80065440 : do i=ifirst,ilast
363 80065440 : v2(i,j) = v(i,j,k)**2
364 : enddo
365 8255520 : v2(ilast+1,j) = veast(j,k)**2
366 : enddo
367 :
368 8386560 : do j=jfirst,jlast
369 83865600 : do i=ifirst,ilast
370 : ! convert to Cv*T
371 81768960 : t2(i,j) = (cp3v(i,j,k)-cap3v(i,j,k)*cp3v(i,j,k))*pt(i,j,k)
372 : enddo
373 : enddo
374 :
375 8255520 : do j=js2g0,jn2g0
376 82162080 : do i=ifirst,ilast
377 295626240 : te(i,j,k) = D0_25 * ( u2(i,j) + u2(i,j+1) + &
378 147813120 : v2(i,j) + v2(i+1,j) ) + &
379 449598240 : t2(i,j)*pkz(i,j,k)
380 : enddo
381 : enddo
382 :
383 : ! WS 99.07.29 : Restructuring creates small round-off. Not clear why...
384 :
385 : ! Do collective Mpisum (in i) for te_sp and te_np below (AAM)
386 : !
387 2096640 : if ( jfirst == 1 ) then
388 : ! South pole
389 851760 : do i=ifirst,ilast
390 786240 : u2_sp(i,k) = u2(i,2)
391 786240 : v2_sp(i,k) = v2(i,2)
392 851760 : t2_sp(i,k) = t2(i,1)
393 : enddo
394 : endif
395 :
396 2096640 : if ( jlast == jm ) then
397 : ! North pole
398 851760 : do i=ifirst,ilast
399 786240 : u2_np(i,k) = u2(i,jm)
400 786240 : v2_np(i,k) = v2(i,jm-1)
401 851760 : t2_np(i,k) = t2(i,jm)
402 : enddo
403 : endif
404 :
405 : ! Compute dz; geo-potential increments
406 8386560 : do j=jfirst,jlast
407 83865600 : do i=ifirst,ilast
408 81768960 : dz(i,j,k) = t2(i,j)*(pk(i,j,k+1)-pk(i,j,k))
409 : enddo
410 : enddo
411 16128 : 1000 continue
412 :
413 : #if defined( SPMD )
414 64512 : allocate( pesouth(ifirst:ilast,km+1) )
415 16128 : if (itot /= im) then
416 : call mp_recv3d( grid%commxy, src, im, km+1, jm, &
417 : ifirst-1, ifirst-1, 1, km+1, jfirst, jlast, &
418 16128 : ifirst-1, ifirst-1, 1, km+1, jfirst, jlast, pewest )
419 : endif
420 : call mp_recv3d( grid%commxy, iam-nprxy_x, im, km+1, jm, &
421 : ifirst, ilast, 1, km+1, jfirst-1, jfirst-1, &
422 16128 : ifirst, ilast, 1, km+1, jfirst-1, jfirst-1, pesouth )
423 : #endif
424 :
425 16128 : if ( jfirst == 1 ) then
426 :
427 : !$omp parallel do &
428 : !$omp default(shared) &
429 : !$omp private(i, k)
430 :
431 66024 : do k = 1, km
432 852264 : do i=ifirst,ilast
433 851760 : tmpik(i,k) = D0_5*( u2_sp(i,k) + v2_sp(i,k) ) + t2_sp(i,k)*pkz(i,1,k)
434 : enddo
435 : enddo
436 :
437 504 : call par_xsum( grid, tmpik, km, te_sp)
438 :
439 : !$omp parallel do &
440 : !$omp default(shared) &
441 : !$omp private(i, k)
442 :
443 66024 : do k = 1, km
444 65520 : te_sp(k) = te_sp(k)/real(im,r8)
445 852264 : do i=ifirst,ilast
446 851760 : te(i, 1,k) = te_sp(k)
447 : enddo
448 : enddo
449 : endif
450 :
451 16128 : if ( jlast == jm ) then
452 :
453 : !$omp parallel do &
454 : !$omp default(shared) &
455 : !$omp private(i, k)
456 :
457 66024 : do k = 1, km
458 852264 : do i=ifirst,ilast
459 851760 : tmpik(i,k) = D0_5*( u2_np(i,k) + v2_np(i,k) ) + t2_np(i,k)*pkz(i,jm,k)
460 : enddo
461 : enddo
462 :
463 504 : call par_xsum( grid, tmpik, km, te_np)
464 :
465 : !$omp parallel do &
466 : !$omp default(shared) &
467 : !$omp private(i, k)
468 :
469 66024 : do k = 1, km
470 65520 : te_np(k) = te_np(k)/real(im,r8)
471 852264 : do i=ifirst,ilast
472 851760 : te(i,jm,k) = te_np(k)
473 : enddo
474 : enddo
475 : endif
476 :
477 : ! Converting pt to t
478 209664 : do i=ifirst,ilast
479 790272 : do j=jfirst,jlast
480 76253184 : do k=1,km
481 76059648 : pt(i,j,k) = pt(i,j,k)*pkz(i,j,k)
482 : enddo
483 : enddo
484 : enddo
485 :
486 16128 : it = itot / nxu
487 16128 : jp = nxu * ( jlast - jfirst + 1 )
488 :
489 : !$omp parallel do &
490 : !$omp default(shared) &
491 : !$omp private(i,j,k,i1w,pe0,pe1,pe2,pe3,log_pe1,log_pe2,ratio) &
492 : !$omp private(dak,bkh,krd, ixj,i1,i2) &
493 : !$omp private(pe1w, pe2w, omga_ik )
494 :
495 : ! do 2000 j=jfirst,jlast
496 64512 : do 2000 ixj=1,jp
497 :
498 48384 : j = jfirst + (ixj-1) / nxu
499 48384 : i1 = ifirst + it * mod(ixj-1, nxu)
500 48384 : i2 = i1 + it - 1
501 :
502 : ! Copy data to local 2D arrays.
503 48384 : i1w = i1-1
504 48384 : if (i1 == 1) i1w = im
505 6386688 : do k=1,km+1
506 82397952 : do i=i1,i2
507 76059648 : pe1(i,k) = pe(i,k,j)
508 82397952 : if (k>1) then
509 75479040 : if (pe1(i,k)-pe1(i,k-1)<lagrangianlevcrit) then
510 0 : write(iulog,*) "Lagrangian levels are crossing", lagrangianlevcrit
511 0 : write(iulog,*) "Run will ABORT!"
512 0 : write(iulog,*) "Suggest to increase FV_NSPLTVRM"
513 0 : do kk=1,km
514 0 : write(iulog,'(A21,I5,A1,3f16.12)') "k,dp(unit=hPa),u,v: ",&
515 0 : kk," ",(pe(i,kk,j)-pe(i,kk-1,j))/100.0_r8,u(i,j,kk),v(i,j,kk)
516 : end do
517 0 : call endrun('te_map: Lagrangian levels are crossing')
518 : endif
519 : endif
520 : enddo
521 6386688 : if( itot == im ) then
522 0 : pe1w(k) = pe(i1w,k,j)
523 : #if defined( SPMD )
524 : else
525 6338304 : pe1w(k) = pewest(k,j)
526 : #endif
527 : endif
528 : enddo
529 :
530 5515776 : do k=1,grid%ks+1
531 71124480 : do i=i1,i2
532 65608704 : pe0(i,k) = grid%ak(k)
533 65608704 : pe2(i,k) = grid%ak(k)
534 71076096 : pe3(i,k) = grid%ak(k)
535 : enddo
536 : enddo
537 :
538 870912 : do k=grid%ks+2,km
539 10741248 : do i=i1,i2
540 9870336 : pe0(i,k) = grid%ak(k) + grid%bk(k)* ps(i,j) ! Remapped PLE based on Old PS
541 10692864 : pe2(i,k) = grid%ak(k) + grid%bk(k)*pe1(i,km+1) ! Remapped PLE based on Updated PS
542 : enddo
543 : enddo
544 :
545 628992 : do i=i1,i2
546 580608 : pe0(i,km+1) = ps(i,j)
547 628992 : pe2(i,km+1) = pe1(i,km+1)
548 : enddo
549 :
550 : ! Ghosting for v mapping
551 870912 : do k=grid%ks+2,km
552 870912 : pe2w(k) = grid%ak(k) + grid%bk(k)*pe1w(km+1)
553 : enddo
554 48384 : pe2w(km+1) = pe1w(km+1)
555 :
556 : ! update ps
557 : ! ---------
558 628992 : do i=i1,i2
559 628992 : ps(i,j) = pe1(i,km+1)
560 : enddo
561 :
562 : ! #######################################################################
563 : ! # ReMap Temperature over log(P)
564 : ! #######################################################################
565 6386688 : do k=1,km+1
566 82446336 : do i=i1,i2
567 76059648 : log_pe1(i,k) = log(pe1(i,k))
568 82397952 : log_pe2(i,k) = log(pe2(i,k))
569 : end do
570 : end do
571 : call map1_ppm ( km, log_pe1, pt, &
572 : km, log_pe2, pt, 0, 0, &
573 : itot, i1-ifirst+1, i2-ifirst+1, &
574 48384 : j, jfirst, jlast, 1, kord)
575 :
576 : ! Update Delta-Pressure (from final remapped pressures)
577 : ! -----------------------------------------------------
578 6338304 : do k=1,km
579 81817344 : do i=i1,i2
580 81768960 : delp(i,j,k) = pe2(i,k+1) - pe2(i,k)
581 : enddo
582 : enddo
583 :
584 : ! Compute omga (dp/dt)
585 : ! --------------------
586 6338304 : do k=2,km+1
587 81817344 : do i=i1,i2
588 81768960 : pe0(i,k) = pe1(i,k) - pe0(i,k) ! Delta-P: PLE(After CD_Core) minus PLE(Remapped based on old PS)
589 : enddo
590 : enddo
591 :
592 : ! C.-C. Chen
593 : ! Map omga
594 6338304 : do k=1,km
595 81817344 : do i=i1,i2
596 75479040 : omga_ik(i,k) = omga(i,k,j)
597 :
598 81768960 : if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then
599 75479040 : if (k == 1) omga_ik(i,k) = D0_0
600 : endif
601 :
602 : end do
603 : end do
604 : call map1_ppm ( km, pe1, omga_ik, &
605 : km, pe2, omga_ik, 0, 0, &
606 : itot, i1-ifirst+1, i2-ifirst+1, &
607 48384 : 1, 1, 1, 1, kord)
608 6338304 : do k=1,km
609 81817344 : do i=i1,i2
610 81768960 : omga(i,k,j) = omga_ik(i,k)
611 : end do
612 : end do
613 :
614 : ! #######################################################################
615 : ! # ReMap Constituents
616 : ! #######################################################################
617 :
618 48384 : if( nq /= 0 ) then
619 48384 : if(kord == 8) then
620 0 : krd = 8
621 : else
622 48384 : krd = 7
623 : endif
624 :
625 : call mapn_ppm_tracer ( km, pe1, tracer, nq, &
626 : km, pe2, i1, i2, &
627 48384 : j, ifirst, ilast, jfirst, jlast, 0, krd)
628 : endif
629 :
630 : ! #######################################################################
631 : ! # ReMap U-Wind
632 : ! #######################################################################
633 :
634 48384 : if(j /= 1) then
635 :
636 47880 : if (am_geom_crrct) then
637 :
638 : ! WS 99.07.29 : protect j==jfirst case
639 0 : if (j > jfirst) then
640 0 : do k=2,km+1
641 0 : do i=i1,i2
642 : ! extensive integral weight -> use cosines
643 0 : pe0(i,k) = (pe1(i,k)*grid%cosp(j) + pe(i,k,j-1)*grid%cosp(j-1)) &
644 0 : / (grid%cosp(j) + grid%cosp(j-1))
645 : enddo
646 : enddo
647 0 : do k=grid%ks+2,km+1
648 0 : bkh = D0_5*grid%bk(k)
649 0 : do i=i1,i2
650 0 : pe3(i,k) = grid%ak(k) + grid%bk(k)*(pe1(i,km+1)*grid%cosp(j) + &
651 0 : pe(i,km+1,j-1)*grid%cosp(j-1)) / &
652 0 : (grid%cosp(j) + grid%cosp(j-1))
653 : enddo
654 : enddo
655 :
656 : #if defined( SPMD )
657 : else
658 : ! WS 99.10.01 : Read in pe(:,:,jfirst-1) from the pesouth buffer
659 0 : do k=2,km+1
660 0 : do i=i1,i2
661 0 : pe0(i,k) = (pe1(i,k)*grid%cosp(j) + pesouth(i,k)*grid%cosp(j-1)) &
662 0 : / (grid%cosp(j) + grid%cosp(j-1))
663 : enddo
664 : enddo
665 0 : do k=grid%ks+2,km+1
666 0 : bkh = D0_5*grid%bk(k)
667 0 : do i=i1,i2
668 0 : pe3(i,k) = grid%ak(k) + grid%bk(k)*(pe1(i,km+1)*grid%cosp(j) + &
669 0 : pesouth(i,km+1)*grid%cosp(j-1)) / &
670 0 : (grid%cosp(j) + grid%cosp(j-1))
671 : enddo
672 : enddo
673 : #endif
674 : endif ! (j > jfirst)
675 :
676 : ! total zonal momentum
677 0 : do i=i1,i2
678 0 : dum(i,j)=0._r8
679 : enddo
680 0 : do k=1,km
681 0 : do i=i1,i2
682 0 : dum(i,j)=dum(i,j)-u(i,j,k)*(pe0(i,k+1)-pe0(i,k))
683 : enddo
684 : enddo
685 :
686 : else ! not am_geom_crrct
687 :
688 : ! WS 99.07.29 : protect j==jfirst case
689 47880 : if (j > jfirst) then
690 4225536 : do k=2,km+1
691 54544896 : do i=i1,i2
692 54512640 : pe0(i,k) = D0_5*(pe1(i,k)+pe(i,k,j-1))
693 : enddo
694 : enddo
695 612864 : do k=grid%ks+2,km+1
696 580608 : bkh = D0_5*grid%bk(k)
697 7580160 : do i=i1,i2
698 7547904 : pe3(i,k) = grid%ak(k) + bkh*(pe1(i,km+1)+pe(i,km+1,j-1))
699 : enddo
700 : enddo
701 : #if defined( SPMD )
702 : else
703 : ! WS 99.10.01 : Read in pe(:,:,jfirst-1) from the pesouth buffer
704 2046744 : do k=2,km+1
705 26420184 : do i=i1,i2
706 26404560 : pe0(i,k) = D0_5*(pe1(i,k)+pesouth(i,k))
707 : enddo
708 : enddo
709 296856 : do k=grid%ks+2,km+1
710 281232 : bkh = D0_5*grid%bk(k)
711 3671640 : do i=i1,i2
712 3656016 : pe3(i,k) = grid%ak(k) + bkh*(pe1(i,km+1)+pesouth(i,km+1))
713 : enddo
714 : enddo
715 : #endif
716 : endif ! (j > jfirst)
717 :
718 : endif ! (am_geom_crrct)
719 :
720 : !-------------------------------
721 :
722 : ! ReMap U-Wind (D-Grid Location)
723 : ! ------------------------------
724 : call map1_ppm ( km, pe0, u, km, pe3, u, &
725 : 0, 0, itot, i1-ifirst+1, i2-ifirst+1, &
726 47880 : j, jfirst, jlast, -1, kord)
727 :
728 47880 : if (am_geom_crrct) then
729 :
730 : ! compute zonal momentum difference due to remapping
731 0 : do k=1,km
732 0 : do i=i1,i2
733 0 : dum(i,j)=dum(i,j)+u(i,j,k)*(pe3(i,k+1)-pe3(i,k))
734 : enddo
735 : enddo
736 :
737 : ! correct zonal wind to preserve momentum
738 0 : do k=1,km
739 0 : do i=i1,i2
740 0 : u(i,j,k)=u(i,j,k)-dum(i,j)/(pe3(i,km+1)-pe3(i,1))
741 : enddo
742 : enddo
743 : endif
744 :
745 47880 : if (am_diag_lbl) then
746 :
747 : ! Remap advective wind increment uc
748 :
749 : call map1_ppm ( km, pe0, uc, km, pe3, uc, &
750 : 0, 0, itot, i1-ifirst+1, i2-ifirst+1, &
751 0 : j, jfirst, jlast, -1, kord)
752 :
753 0 : do k=1,km
754 0 : do i=i1,i2
755 0 : du_w(i,j,k)=du_s(k)
756 : enddo
757 : enddo
758 : call map1_ppm ( km, pe0, du_w, km, pe3, du_w, &
759 : 0, 0, itot, i1-ifirst+1, i2-ifirst+1, &
760 0 : j, jfirst, jlast, -1, kord)
761 : endif
762 :
763 : ! ReMap Y-Mass Flux (C-Grid Location)
764 : ! -----------------------------------
765 6272280 : do k=1,km
766 80965080 : do i=i1,i2
767 80917200 : mfy(i,j,k) = mfy(i,j,k)/(pe0(i,k+1)-pe0(i,k))
768 : enddo
769 : enddo
770 : call map1_ppm ( km, pe0, mfy, km, pe3, mfy, &
771 : 0, 0, itot, i1-ifirst+1, i2-ifirst+1, &
772 47880 : j, jfirst, jlast, -1, kord)
773 6272280 : do k=1,km
774 80965080 : do i=i1,i2
775 80917200 : mfy(i,j,k) = mfy(i,j,k)*(pe3(i,k+1)-pe3(i,k))
776 : enddo
777 : enddo
778 : endif
779 :
780 : ! #######################################################################
781 : ! # ReMap V-Wind
782 : ! #######################################################################
783 :
784 48384 : if(j /= 1 .and. j /= jm) then
785 6206256 : do k=2,km+1
786 : ! pe1(i1-1,1:km+1) must be ghosted
787 6158880 : pe0(i1,k) = D0_5*(pe1(i1,k)+pe1w(k))
788 73953936 : do i=i1+1,i2
789 73906560 : pe0(i ,k) = D0_5*(pe1(i,k)+pe1(i-1,k))
790 : enddo
791 : enddo
792 :
793 900144 : do k=grid%ks+2,km+1
794 : ! pe2(i1-1,grid%ks+2:km+1) must be ghosted
795 852768 : pe3(i1,k) = D0_5*(pe2(i1,k)+pe2w(k))
796 10280592 : do i=i1+1,i2
797 10233216 : pe3(i,k) = D0_5*(pe2(i,k)+pe2(i-1,k))
798 : enddo
799 : enddo
800 :
801 : ! ReMap V-Wind (D-Grid Location)
802 : ! ------------------------------
803 : call map1_ppm ( km, pe0, v, km, pe3, v, &
804 : 0, 0, itot, i1-ifirst+1, i2-ifirst+1, &
805 47376 : j, jfirst, jlast, -1, kord)
806 :
807 :
808 47376 : if (am_diag_lbl) then
809 : call map1_ppm ( km, pe0, vc, km, pe3, vc, &
810 : 0, 0, itot, i1-ifirst+1, i2-ifirst+1, &
811 0 : j, jfirst, jlast, -1, kord)
812 : end if
813 :
814 : ! ReMap X-Mass Flux (C-Grid Location)
815 : ! -----------------------------------
816 6206256 : do k=1,km
817 80112816 : do i=i1,i2
818 80065440 : mfx(i,j,k) = mfx(i,j,k)/(pe0(i,k+1)-pe0(i,k))
819 : enddo
820 : enddo
821 : call map1_ppm ( km, pe0, mfx, km, pe3, mfx, &
822 : 0, 0, itot, i1-ifirst+1, i2-ifirst+1, &
823 47376 : j, jfirst, jlast, -1, kord)
824 6206256 : do k=1,km
825 80112816 : do i=i1,i2
826 80065440 : mfx(i,j,k) = mfx(i,j,k)*(pe3(i,k+1)-pe3(i,k))
827 : enddo
828 : enddo
829 : endif
830 :
831 : ! Save new PE to temp storage peln
832 : ! --------------------------------
833 6289920 : do k=2,km
834 81188352 : do i=i1,i2
835 81139968 : peln(i,k,j) = pe2(i,k)
836 : enddo
837 : enddo
838 :
839 : ! Check deformation.
840 48384 : if( diag ) then
841 0 : rmax(ixj) = D0_0
842 0 : rmin(ixj) = D1_0
843 0 : do k=1,km
844 0 : do i=i1,i2
845 0 : ratio(i) = (pe1(i,k+1)-pe1(i,k)) / (pe2(i,k+1)-pe2(i,k))
846 : enddo
847 :
848 0 : do i=i1,i2
849 0 : if(ratio(i) > rmax(ixj)) then
850 0 : rmax(ixj) = ratio(i)
851 0 : elseif(ratio(i) < rmin(ixj)) then
852 0 : rmin(ixj) = ratio(i)
853 : endif
854 : enddo
855 : enddo
856 : endif
857 16128 : 2000 continue
858 :
859 16128 : if (high_alt) then
860 16128 : call cam_thermo_calc_kappav( tracer, cap3v, cpv=cp3v )
861 : !$omp parallel do private(i,j,k)
862 2096640 : do k=2,km
863 8338176 : do j=jfirst,jlast
864 83220480 : do i=ifirst,ilast
865 81139968 : cap3vi(i,j,k) = 0.5_r8*(cap3v(i,j,k-1)+cap3v(i,j,k))
866 : enddo
867 : enddo
868 : enddo
869 645120 : cap3vi(:,:,1) = 1.5_r8 * cap3v(:,:,1) - 0.5_r8 * cap3v(:,:,2)
870 645120 : cap3vi(:,:,km+1) = 1.5_r8 * cap3v(:,:,km) - 0.5_r8 * cap3v(:,:,km-1)
871 : else
872 0 : cap3vi(:,:,:) = cap3v(grid%ifirstxy,grid%jfirstxy,1)
873 : endif
874 :
875 : #if defined( SPMD )
876 16128 : deallocate( pesouth )
877 :
878 : ! Send u southward
879 : call mp_send3d( grid%commxy, iam-nprxy_x, iam+nprxy_x, im, jm, km,&
880 : ifirst, ilast, jfirst, jlast, 1, km, &
881 16128 : ifirst, ilast, jfirst, jfirst, 1, km, u )
882 16128 : if (itot /= im) then
883 16128 : dest = myidxy_y*nprxy_x + MOD(iam+nprxy_x-1,nprxy_x)
884 16128 : src = myidxy_y*nprxy_x + MOD(iam+1,nprxy_x)
885 : call mp_send3d( grid%commxy, dest, src, im, jm, km, &
886 : ifirst, ilast, jfirst, jlast, 1, km, &
887 16128 : ifirst, ifirst, jfirst, jlast, 1, km, v )
888 : endif
889 : #endif
890 :
891 16128 : if( diag ) then
892 0 : qmin = rmin(1)
893 0 : do ixj=2, jp
894 0 : if(rmin(ixj) < qmin) then
895 0 : qmin = rmin(ixj)
896 : endif
897 : enddo
898 0 : CPP_PRT_PREFIX write(iulog,*) 'rmin=', qmin
899 :
900 0 : qmax = rmax(1)
901 0 : do ixj=2, jp
902 0 : if(rmax(ixj) > qmax) then
903 0 : qmax = rmax(ixj)
904 : endif
905 : enddo
906 0 : CPP_PRT_PREFIX write(iulog,*) 'rmax=', qmax
907 : endif
908 :
909 : ! Recover Final Edge-Pressures and Compute Mid-Level PKZ
910 : ! ------------------------------------------------------
911 :
912 : !$omp parallel do &
913 : !$omp default(shared) &
914 : !$omp private(i,j,k)
915 :
916 64512 : do j=jfirst,jlast
917 6306048 : do k=2,km
918 81188352 : do i=ifirst,ilast
919 81139968 : pe(i,k,j) = peln(i,k,j)
920 : enddo
921 : enddo
922 : enddo
923 :
924 2128896 : do k=1,km+1
925 8467200 : do j=jfirst,jlast
926 84510720 : do i=ifirst,ilast
927 82397952 : pk(i,j,k) = pe(i,k,j)**cap3vi(i,j,k)
928 : enddo
929 : enddo
930 : enddo
931 : call pkez(nxu, im, km, jfirst, jlast, 1, km, ifirst, ilast, &
932 16128 : pe, pk, cap3v, grid%ks, peln, pkz, .false., high_alt)
933 :
934 : ! Single x-subdomain case (periodic)
935 2112768 : do k = 1, km
936 8402688 : do j = jfirst, jlast
937 8386560 : veast(j,k) = v(ifirst,j,k)
938 : enddo
939 : enddo
940 :
941 : #if defined( SPMD )
942 : ! Recv u from north
943 : call mp_recv3d( grid%commxy, iam+nprxy_x, im, jm, km, &
944 : ifirst, ilast, jlast+1, jlast+1, 1, km, &
945 16128 : ifirst, ilast, jlast+1, jlast+1, 1, km, unorth )
946 16128 : if (itot /= im) then
947 : call mp_recv3d( grid%commxy, src, im, jm, km, &
948 : ilast+1, ilast+1, jfirst, jlast, 1, km, &
949 16128 : ilast+1, ilast+1, jfirst, jlast, 1, km, veast )
950 : endif
951 : #endif
952 :
953 : ! ((((((((((((((((( compute globally integrated TE )))))))))))))))))
954 :
955 16128 : if( consv ) then
956 :
957 : !$omp parallel do &
958 : !$omp default(shared) &
959 : !$omp private(i,j,k)
960 :
961 0 : do k=1,km
962 0 : do j=jfirst,jlast
963 0 : do i=ifirst,ilast
964 0 : dz(i,j,k) = te(i,j,k) * delp(i,j,k)
965 : enddo
966 : enddo
967 : enddo
968 :
969 : !$omp parallel do &
970 : !$omp default(shared) &
971 : !$omp private(i,j,k,bte)
972 :
973 : ! Perform vertical integration
974 :
975 0 : do 4000 j=jfirst,jlast
976 :
977 0 : if ( j == 1 ) then
978 : ! SP
979 0 : tte(1) = D0_0
980 :
981 0 : do k=1,km
982 0 : tte(1) = tte(1) + dz(ifirst,1,k)
983 : enddo
984 :
985 0 : elseif ( j == jm) then
986 : ! NP
987 0 : tte(jm) = D0_0
988 :
989 0 : do k=1,km
990 0 : tte(jm) = tte(jm) + dz(ifirst,jm,k)
991 : enddo
992 :
993 : else
994 : ! Interior
995 0 : do i=ifirst,ilast
996 0 : bte(i) = D0_0
997 : enddo
998 :
999 0 : do k=1,km
1000 0 : do i=ifirst,ilast
1001 0 : bte(i) = bte(i) + dz(i,j,k)
1002 : enddo
1003 : enddo
1004 :
1005 0 : do i=ifirst,ilast
1006 0 : tmpij(i,j,1) = bte(i)
1007 : enddo
1008 :
1009 : endif
1010 0 : 4000 continue
1011 :
1012 0 : call par_xsum( grid, tmpij, jlast-jfirst+1, xysum)
1013 :
1014 : !$omp parallel do &
1015 : !$omp default(shared) &
1016 : !$omp private(j)
1017 :
1018 0 : do j = max(jfirst,2), min(jlast,jm-1)
1019 0 : tte(j) = xysum(j,1)*grid%cosp(j)
1020 : enddo
1021 :
1022 0 : if ( jfirst == 1 ) tte(1) = grid%acap * tte(1)
1023 0 : if ( jlast == jm ) tte(jm) = grid%acap * tte(jm)
1024 :
1025 0 : te1 = D0_0
1026 0 : call par_vecsum(jm, jfirst, jlast, tte, te1, comm_use, npry_use)
1027 :
1028 : endif ! consv
1029 :
1030 16128 : if( consv ) then
1031 :
1032 : !$omp parallel do &
1033 : !$omp& default(shared) &
1034 : !$omp& private(i,j)
1035 :
1036 0 : do j=js2g0, jn2g0
1037 0 : do i=ifirst,ilast
1038 0 : tmpij(i,j,1) = ps(i,j)
1039 0 : tmpij(i,j,2) = peln(i,km+1,j)
1040 : enddo
1041 : enddo
1042 :
1043 0 : call par_xsum( grid, tmpij, 2*(jlast-jfirst+1), xysum)
1044 :
1045 : !$omp parallel do &
1046 : !$omp default(shared) &
1047 : !$omp private(j)
1048 :
1049 0 : do j=js2g0, jn2g0
1050 0 : tte(j) = cp3v(ifirst,j,1)*grid%cosp(j)*(xysum(j,1) - grid%ptop*real(im,r8) - &
1051 0 : cap3v(ifirst,j,1)*grid%ptop*(xysum(j,2) - peln(ifirst,1,j)*real(im,r8)) )
1052 : ! peln(i,1,j) should be independent of i (AAM)
1053 : enddo
1054 :
1055 0 : if ( jfirst == 1 ) tte(1) = grid%acap*cp3v(ifirst,1,km) * (ps(ifirst,1) - grid%ptop - &
1056 0 : cap3v(ifirst,1,km)*grid%ptop*(peln(ifirst,km+1,1) - peln(ifirst,1,1) ) )
1057 0 : if ( jlast == jm ) tte(jm)= grid%acap*cp3v(ifirst,jm,km) * (ps(ifirst,jm) - grid%ptop - &
1058 0 : cap3v(ifirst,jm,km)*grid%ptop*(peln(ifirst,km+1,jm) - peln(ifirst,1,jm) ) )
1059 : endif ! consv
1060 :
1061 16128 : if (consv) then
1062 :
1063 0 : sum=D0_0
1064 0 : call par_vecsum(jm, jfirst, jlast, tte, sum, comm_use, npry_use)
1065 :
1066 0 : dtmp = (te0 - te1) / sum
1067 0 : if( diag ) then
1068 0 : CPP_PRT_PREFIX write(iulog,*) 'te=',te0, ' Energy deficit in T = ', dtmp
1069 : endif
1070 :
1071 : endif ! end consv check
1072 :
1073 : !$omp parallel do &
1074 : !$omp default(shared) &
1075 : !$omp private(i,j,k, u2, v2)
1076 :
1077 : ! --------------------------------------------------------------------
1078 : ! --- Recover Tv from remapped Total Energy and its components ---
1079 : ! --------------------------------------------------------------------
1080 :
1081 2112768 : do 8000 k=1,km
1082 :
1083 : ! Intialize Kinetic Energy
1084 : ! ------------------------
1085 8321040 : do j=js2g0,jlast
1086 83013840 : do i=ifirst,ilast
1087 80917200 : u2(i,j) = u(i,j,k)**2
1088 : enddo
1089 : enddo
1090 : #if defined( SPMD )
1091 2096640 : if ( jlast < jm ) then
1092 26404560 : do i=ifirst,ilast
1093 26404560 : u2(i,jlast+1) = unorth(i,k)**2 ! fill ghost zone
1094 : enddo
1095 : endif
1096 : #endif
1097 :
1098 8255520 : do j=js2g0,jn2g0
1099 80065440 : do i=ifirst,ilast
1100 80065440 : v2(i,j) = v(i,j,k)**2
1101 : enddo
1102 8255520 : v2(ilast+1,j) = veast(j,k)**2
1103 : enddo
1104 :
1105 : ! Subtract Kinetic Energy from Total Energy (Leaving Internal + Potential)
1106 : ! ------------------------------------------------------------------------
1107 8255520 : do j=js2g0,jn2g0
1108 82162080 : do i=ifirst,ilast
1109 295626240 : te(i,j,k) = D0_25 * ( u2(i,j) + u2(i,j+1) &
1110 375691680 : +v2(i,j) + v2(i+1,j) )
1111 : enddo
1112 : enddo
1113 :
1114 : ! South pole
1115 : ! ----------
1116 2096640 : if ( jfirst == 1 ) then
1117 851760 : do i=ifirst,ilast
1118 786240 : u2_sp(i,k) = u2(i,2)
1119 851760 : v2_sp(i,k) = v2(i,2)
1120 : enddo
1121 : endif
1122 :
1123 : ! North pole
1124 : ! ----------
1125 2096640 : if ( jlast == jm ) then
1126 851760 : do i=ifirst,ilast
1127 786240 : u2_np(i,k) = u2(i,jm)
1128 851760 : v2_np(i,k) = v2(i,jm-1)
1129 : enddo
1130 : endif
1131 :
1132 16128 : 8000 continue
1133 :
1134 : ! South pole
1135 : ! ----------
1136 16128 : if ( jfirst == 1 ) then
1137 :
1138 : !$omp parallel do &
1139 : !$omp default(shared) &
1140 : !$omp private(i, k)
1141 :
1142 66024 : do k = 1, km
1143 852264 : do i=ifirst,ilast
1144 851760 : tmpik(i,k) = D0_5*( u2_sp(i,k) + v2_sp(i,k) )
1145 : enddo
1146 : enddo
1147 :
1148 504 : call par_xsum( grid, tmpik, km, te_sp)
1149 :
1150 : !$omp parallel do &
1151 : !$omp default(shared) &
1152 : !$omp private(i, k)
1153 :
1154 66024 : do k = 1, km
1155 65520 : te_sp(k) = te_sp(k)/real(im,r8)
1156 852264 : do i=ifirst,ilast
1157 851760 : te(i,1,k) = te_sp(k)
1158 : enddo
1159 : enddo
1160 : endif
1161 :
1162 : ! North pole
1163 : ! ----------
1164 16128 : if ( jlast == jm ) then
1165 :
1166 : !$omp parallel do &
1167 : !$omp default(shared) &
1168 : !$omp private(i, k)
1169 :
1170 66024 : do k = 1, km
1171 852264 : do i=ifirst,ilast
1172 851760 : tmpik(i,k) = D0_5*( u2_np(i,k) + v2_np(i,k) )
1173 : enddo
1174 : enddo
1175 :
1176 504 : call par_xsum( grid, tmpik, km, te_np)
1177 :
1178 : !$omp parallel do &
1179 : !$omp default(shared) &
1180 : !$omp private(i, k)
1181 :
1182 66024 : do k = 1, km
1183 65520 : te_np(k) = te_np(k)/real(im,r8)
1184 852264 : do i=ifirst,ilast
1185 851760 : te(i,jm,k) = te_np(k)
1186 : enddo
1187 : enddo
1188 : endif
1189 :
1190 :
1191 16128 : if( .not. convt ) then
1192 104832 : do i=ifirst,ilast
1193 395136 : do j=jfirst,jlast
1194 38126592 : do k=1,km
1195 38029824 : pt(i,j,k) = pt(i,j,k)/pkz(i,j,k) ! Scaled Virtual Potential Temperature
1196 : enddo
1197 : enddo
1198 : enddo
1199 : endif
1200 :
1201 16128 : return
1202 : !EOC
1203 16128 : end subroutine te_map
1204 : !-----------------------------------------------------------------------
1205 : end module te_map_mod
|
