Line data Source code
1 : !-----------------------------------------------------------------------
2 : !BOP
3 : ! !ROUTINE: geopk --- Calculate geopotential to the kappa
4 : !
5 : !-----------------------------------------------------------------------
6 : ! There are three versions of geopk below. The first is the standard
7 : ! version and is typically used with transposes between yz and xy
8 : ! space. The second (called geopk16) operates in yz space and performs
9 : ! semi-global communication in the z direction (to avoid transposes).
10 : ! It also can use 16-byte reals to preserve accuracy through round-off;
11 : ! this is accomplished by toggling DSIZE to 16 immediately below.
12 : ! The third version (called geopk_d) also operates in yz space
13 : ! and implements a ring-pipeline algorithm in the z direction.
14 : ! Numerics are identical with the first version without requiring
15 : ! 16-byte arithmetic. While less parallel, communication costs are
16 : ! smaller, and this is often the fastest option.
17 : !
18 : ! Note that the interfaces to the first, second, and third versions are
19 : ! slightly different. Also, geopk (the standard version with transposes)
20 : ! is called for the D-grid during the last two small timesteps in cd_core.
21 : ! Geopk16 uses mod_comm communication calls; one can activate the old
22 : ! Pilgrim calls (for debugging) by activating PaREXCH immediately below.
23 :
24 : !#define PAREXCH
25 : !#define DSIZE 16
26 : #define DSIZE 8
27 :
28 : #if (DSIZE == 16)
29 : # define DTWO 2
30 : #else
31 : # define DTWO 1
32 : #endif
33 : !-----------------------------------------------------------------------
34 : !
35 : ! !INTERFACE:
36 258048 : subroutine geopk(grid, pe, delp, pk, wz, hs, pt, cp3v, cap3v, nx)
37 :
38 : use shr_kind_mod, only: r8 => shr_kind_r8
39 : use dynamics_vars, only: T_FVDYCORE_GRID
40 :
41 : implicit none
42 :
43 : ! !INPUT PARAMETERS:
44 :
45 : type (T_FVDYCORE_GRID), intent(in) :: grid
46 : integer nx ! # of pieces in longitude direction
47 : real(r8) cp3v(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km)
48 : real(r8) cap3v(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km)
49 : real(r8) hs(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy)
50 : real(r8) pt(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km)
51 : real(r8) delp(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km)
52 :
53 : ! !OUTPUT PARAMETERS:
54 : real(r8) wz(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km+1) ! space N*1 S*1
55 : real(r8) pk(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km+1) ! space N*1 S*1
56 : real(r8) pe(grid%ifirstxy:grid%ilastxy,grid%km+1,grid%jfirstxy:grid%jlastxy)
57 :
58 : ! !DESCRIPTION:
59 : ! Calculates geopotential and pressure to the kappa. This is an expensive
60 : ! operation and several out arrays are kept around for future use.
61 : !
62 : ! !REVISION HISTORY:
63 : !
64 : ! WS 99.10.22: MPIed SJ's original SMP version
65 : ! SJL 00.01.01: Merged C-core and D-core computation
66 : ! SMP "decmposition" in E-W by combining i and j loops
67 : ! WS 00.12.01: Replaced MPI_ON with SPMD; hs now distributed
68 : ! AAM 01.06.27: Generalize for 2D decomposition
69 : ! AAM 01.07.24: Removed dpcheck
70 : ! WS 04.10.07: Simplified interface using Grid as input argument
71 : ! WS 05.05.25: Merged CAM and GEOS5 versions (mostly CAM)
72 : !
73 : !EOP
74 : !---------------------------------------------------------------------
75 : !BOC
76 :
77 : ! Local:
78 : real(r8), parameter :: D0_0 = 0.0_r8
79 : integer :: im, jm, km, jfirst, jlast, ifirst, ilast
80 : real(r8) :: ptop
81 :
82 : integer i, j, k
83 : integer ixj, jp, it, i1, i2, nxu, itot
84 516096 : real(r8) delpp(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km)
85 516096 : real(r8) cap3vi(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km+1)
86 516096 : real(r8) peln(grid%ifirstxy:grid%ilastxy,grid%km+1,grid%jfirstxy:grid%jlastxy)
87 :
88 : logical :: high_alt
89 258048 : high_alt = grid%high_alt
90 :
91 258048 : ptop = grid%ptop
92 258048 : im = grid%im
93 258048 : jm = grid%jm
94 258048 : km = grid%km
95 258048 : ifirst = grid%ifirstxy
96 258048 : ilast = grid%ilastxy
97 258048 : jfirst = grid%jfirstxy
98 258048 : jlast = grid%jlastxy
99 :
100 258048 : itot = ilast - ifirst + 1
101 : ! nxu = nx
102 258048 : nxu = 1
103 258048 : it = itot / nxu
104 258048 : jp = nxu * ( jlast - jfirst + 1 )
105 :
106 258048 : if (grid%high_alt) then
107 : !$omp parallel do private(i,j,k)
108 0 : do k=2,km
109 0 : do j=grid%jfirstxy,grid%jlastxy
110 0 : do i=grid%ifirstxy,grid%ilastxy
111 0 : cap3vi(i,j,k) = 0.5_r8*(cap3v(i,j,k-1)+cap3v(i,j,k))
112 : enddo
113 : enddo
114 : enddo
115 0 : cap3vi(:,:,1) = 1.5_r8 * cap3v(:,:,1) - 0.5_r8 * cap3v(:,:,2)
116 0 : cap3vi(:,:,km+1) = 1.5_r8 * cap3v(:,:,km) - 0.5_r8 * cap3v(:,:,km-1)
117 : else
118 647442432 : cap3vi(:,:,:) = cap3v(grid%ifirstxy,grid%jfirstxy,1)
119 : endif
120 :
121 : !$omp parallel do &
122 : !$omp default(shared) &
123 : !$omp private(i1, i2, ixj, i, j, k )
124 :
125 : ! do 2000 j=jfirst,jlast
126 1032192 : do 2000 ixj=1, jp
127 :
128 774144 : j = jfirst + (ixj-1)/nxu
129 774144 : i1 = ifirst + it * mod(ixj-1, nxu)
130 774144 : i2 = i1 + it - 1
131 :
132 19353600 : do i=i1,i2
133 18579456 : pe(i,1,j) = D0_0
134 19353600 : wz(i,j,km+1) = D0_0
135 : enddo
136 :
137 : ! Top down
138 25546752 : do k=2,km+1
139 620089344 : do i=i1,i2
140 619315200 : pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
141 : enddo
142 : enddo
143 26320896 : do k=1,km+1
144 639442944 : do i=i1,i2
145 613122048 : pe(i,k,j) = pe(i,k,j) + ptop
146 613122048 : peln(i,k,j) = log(pe(i,k,j))
147 638668800 : pk(i,j,k) = pe(i,k,j)**cap3vi(i,j,k)
148 : enddo
149 : enddo
150 :
151 : ! Bottom up
152 774144 : if (high_alt) then
153 0 : do k=1,km
154 0 : do i=i1,i2
155 0 : delpp(i,j,k) = cp3v(i,j,k)*cap3v(i,j,k)*pt(i,j,k)*pk(i,j,k)*(peln(i,k+1,j)-peln(i,k,j))
156 : enddo
157 : enddo
158 : else
159 25546752 : do k=1,km
160 620089344 : do i=i1,i2
161 619315200 : delpp(i,j,k) = cp3v(i,j,k)*pt(i,j,k)*(pk(i,j,k+1)-pk(i,j,k))
162 : enddo
163 : enddo
164 : endif
165 :
166 25546752 : do k=km,1,-1
167 620089344 : do i=i1,i2
168 619315200 : wz(i,j,k) = wz(i,j,k+1)+delpp(i,j,k)
169 : enddo
170 : enddo
171 26320896 : do k=1,km+1
172 639442944 : do i=i1,i2
173 638668800 : wz(i,j,k) = wz(i,j,k)+hs(i,j)
174 : enddo
175 : enddo
176 258048 : 2000 continue
177 :
178 258048 : return
179 : !EOC
180 : end subroutine geopk
181 : !-----------------------------------------------------------------------
182 :
183 :
184 : !-----------------------------------------------------------------------
185 : !BOP
186 : ! !ROUTINE: geopk16 --- Calculate geopotential to the kappa
187 : !
188 : ! !INTERFACE:
189 0 : subroutine geopk16(grid, pe, delp, pk, wz, hs, pt, ng, cp, akap )
190 :
191 : use shr_kind_mod, only : r8 => shr_kind_r8, i8 => shr_kind_i8
192 : use decompmodule, only : decomptype
193 : use dynamics_vars, only : T_FVDYCORE_GRID
194 :
195 : #if defined( SPMD )
196 : use parutilitiesmodule, only : parexchangevector
197 : use mod_comm, only : blockdescriptor, get_partneroffset, &
198 : mp_sendirr, mp_recvirr, max_nparcels
199 : #endif
200 :
201 : implicit none
202 :
203 : #if defined ( SPMD )
204 : #include "mpif.h"
205 : #endif
206 :
207 : ! !INPUT PARAMETERS:
208 :
209 : type (T_FVDYCORE_GRID), intent(in) :: grid
210 : integer, intent(in) :: ng ! Halo size (not always = ng_d)
211 :
212 : real(r8) akap, cp
213 : real(r8) hs(1:grid%im,grid%jfirst:grid%jlast)
214 :
215 : ! !INPUT PARAMETERS:
216 : real(r8) pt(1:grid%im,grid%jfirst-ng:grid%jlast+ng,grid%kfirst:grid%klast)
217 : real(r8) delp(1:grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast)
218 :
219 : ! !OUTPUT PARAMETERS:
220 : real(r8) wz(1:grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast+1) ! space N*1 S*1
221 : real(r8) pk(1:grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast+1) ! space N*1 S*1
222 : real(r8) pe(1:grid%im,grid%kfirst:grid%klast+1,grid%jfirst:grid%jlast) ! temporary variable
223 :
224 : ! !DESCRIPTION:
225 : ! Calculates geopotential and pressure to the kappa. This is an expensive
226 : ! operation and several out arrays are kept around for future use.
227 : ! To preserve accuracy through round-off, 16-byte reals are used
228 : ! for some intermediate data.
229 : !
230 : ! !REVISION HISTORY:
231 : !
232 : ! AAM 00.12.18: Original version
233 : ! AAM 03.01.21: Use mod_comm
234 : ! WS 03.11.19: Merged latest CAM version (by AAM)
235 : ! WS 04.10.07: Simplified interface using Grid as input argument
236 : ! WS 05.05.17: Merged CAM and GEOS5 versions
237 : !
238 : !EOP
239 : !---------------------------------------------------------------------
240 : !BOC
241 :
242 : #ifndef NO_CRAY_POINTERS
243 :
244 : ! Local:
245 : integer :: i, j, k, nk, ijtot, ierror, ione
246 :
247 : integer :: im,jm,km, ifirst, ilast, jfirst, jlast, kfirst, klast
248 : real(r8):: ptop
249 :
250 : integer :: npr_y, npr_z, npes_yz, myid_y, myid_z
251 : integer :: twod_decomp, mod_geopk
252 :
253 : #if (DSIZE == 16)
254 : #ifdef NO_R16
255 : integer,parameter :: r16= selected_real_kind(12) ! 8 byte real
256 : #else
257 : integer,parameter :: r16= selected_real_kind(24) ! 16 byte real
258 : #endif
259 : real(r16), parameter :: DP0_0 = 0.0_r16
260 : real(r16) delp16(1:grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast)
261 : real(r16) pe16(1:grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast+1)
262 : real(r16) inbuf(1:grid%im,grid%jfirst:grid%jlast,0:grid%npr_z-1)
263 : real(r16) outbuf(1:grid%im,grid%jfirst:grid%jlast,0:grid%npr_z-1)
264 : #else
265 : real (r8), parameter :: DP0_0 = 0.0_r8
266 0 : real (r8) delp16(1:grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast)
267 0 : real (r8) pe16(1:grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast+1)
268 : real (r8) inbuf(1:grid%im,grid%jfirst:grid%jlast,0:grid%npr_z-1)
269 0 : real (r8) outbuf(1:grid%im,grid%jfirst:grid%jlast,0:grid%npr_z-1)
270 : #endif
271 : integer sendcount(0:grid%npr_z-1), recvcount(0:grid%npr_z-1)
272 :
273 : #if defined(SPMD)
274 : !
275 : ! data structures for mp_sendirr, mp_recvirr
276 : !
277 : type (blockdescriptor), allocatable, save :: sendbl1(:), recvbl1(:)
278 : type (blockdescriptor), allocatable, save :: sendbl2(:), recvbl2(:)
279 :
280 : #endif
281 :
282 : integer first_time_through
283 : data first_time_through / 0 /
284 :
285 : ! Arrays inbuf8 and outbuf8 are created to fool the compiler
286 : ! into accepting them as calling arguments for parexchangevector.
287 : ! The trickery below equivalences them to inbuf and outbuf
288 : real (r8) inbuf8(1), outbuf8(1)
289 : pointer (ptr_inbuf8, inbuf8)
290 : pointer (ptr_outbuf8, outbuf8)
291 : integer (i8) locinbuf, locoutbuf
292 :
293 : !
294 : ! Initialize variables from Grid
295 : !
296 0 : ptop = grid%ptop
297 :
298 0 : im = grid%im
299 0 : jm = grid%jm
300 0 : km = grid%km
301 :
302 0 : ifirst = 1 ! 2004.10.04 (WS): Now hardwired for 1..im
303 0 : ilast = grid%im ! Code was always used in this mode
304 0 : jfirst = grid%jfirst
305 0 : jlast = grid%jlast
306 0 : kfirst = grid%kfirst
307 0 : klast = grid%klast
308 :
309 0 : myid_y = grid%myid_y
310 0 : myid_z = grid%myid_z
311 :
312 0 : npr_y = grid%npr_y
313 0 : npr_z = grid%npr_z
314 0 : npes_yz = grid%npes_yz
315 :
316 0 : twod_decomp = grid%twod_decomp
317 0 : mod_geopk = grid%mod_geopk
318 :
319 0 : ijtot = (jlast-jfirst+1) * (ilast-ifirst+1)
320 :
321 : #if defined (SPMD)
322 0 : if (first_time_through .eq. 0) then
323 0 : first_time_through = 1
324 0 : ione = 1
325 0 : if (npr_z .gt. 1) then
326 0 : allocate( sendbl1(0:npes_yz-1) )
327 0 : allocate( recvbl1(0:npes_yz-1) )
328 0 : allocate( sendbl2(0:npes_yz-1) )
329 0 : allocate( recvbl2(0:npes_yz-1) )
330 :
331 0 : do nk = 0,npes_yz-1
332 :
333 0 : sendbl1(nk)%method = mod_geopk
334 0 : sendbl2(nk)%method = mod_geopk
335 0 : recvbl1(nk)%method = mod_geopk
336 0 : recvbl2(nk)%method = mod_geopk
337 :
338 : ! allocate for either method (safety)
339 0 : allocate( sendbl1(nk)%blocksizes(1) )
340 0 : allocate( sendbl1(nk)%displacements(1) )
341 0 : allocate( recvbl1(nk)%blocksizes(1) )
342 0 : allocate( recvbl1(nk)%displacements(1) )
343 0 : allocate( sendbl2(nk)%blocksizes(1) )
344 0 : allocate( sendbl2(nk)%displacements(1) )
345 0 : allocate( recvbl2(nk)%blocksizes(1) )
346 0 : allocate( recvbl2(nk)%displacements(1) )
347 :
348 0 : sendbl1(nk)%type = MPI_DATATYPE_NULL
349 :
350 0 : if ( (nk/npr_y) > myid_z .and. mod(nk,npr_y) == myid_y ) then
351 :
352 0 : if (mod_geopk .ne. 0) then
353 : call MPI_TYPE_INDEXED(ione, DTWO*ijtot, &
354 : DTWO*ijtot*(klast-kfirst+1), MPI_DOUBLE_PRECISION, &
355 0 : sendbl1(nk)%type, ierror)
356 0 : call MPI_TYPE_COMMIT(sendbl1(nk)%type, ierror)
357 : endif
358 :
359 0 : sendbl1(nk)%blocksizes(1) = DTWO*ijtot
360 0 : sendbl1(nk)%displacements(1) = DTWO*ijtot*(klast-kfirst+1)
361 0 : sendbl1(nk)%partneroffset = myid_z * ijtot * DTWO
362 :
363 : else
364 :
365 0 : sendbl1(nk)%blocksizes(1) = 0
366 0 : sendbl1(nk)%displacements(1) = 0
367 0 : sendbl1(nk)%partneroffset = 0
368 :
369 : endif
370 0 : sendbl1(nk)%nparcels = size(sendbl1(nk)%displacements)
371 0 : sendbl1(nk)%tot_size = sum(sendbl1(nk)%blocksizes)
372 0 : max_nparcels = max(max_nparcels, sendbl1(nk)%nparcels)
373 :
374 0 : recvbl1(nk)%type = MPI_DATATYPE_NULL
375 :
376 0 : if ( (nk/npr_y) < myid_z .and. mod(nk,npr_y) == myid_y ) then
377 :
378 0 : if (mod_geopk .ne. 0) then
379 : call MPI_TYPE_INDEXED(ione, DTWO*ijtot, &
380 : nk/npr_y * ijtot * DTWO, MPI_DOUBLE_PRECISION, &
381 0 : recvbl1(nk)%type, ierror)
382 0 : call MPI_TYPE_COMMIT(recvbl1(nk)%type, ierror)
383 : endif
384 :
385 0 : recvbl1(nk)%blocksizes(1) = DTWO*ijtot
386 0 : recvbl1(nk)%displacements(1) = nk/npr_y * ijtot * DTWO
387 0 : recvbl1(nk)%partneroffset = 0
388 :
389 : else
390 :
391 0 : recvbl1(nk)%blocksizes(1) = 0
392 0 : recvbl1(nk)%displacements(1) = 0
393 0 : recvbl1(nk)%partneroffset = 0
394 :
395 : endif
396 0 : recvbl1(nk)%nparcels = size(recvbl1(nk)%displacements)
397 0 : recvbl1(nk)%tot_size = sum(recvbl1(nk)%blocksizes)
398 0 : max_nparcels = max(max_nparcels, recvbl1(nk)%nparcels)
399 :
400 0 : if ( (nk/npr_y) < myid_z .and. mod(nk,npr_y) == myid_y ) then
401 :
402 : call MPI_TYPE_INDEXED(ione, DTWO*ijtot, &
403 : 0, MPI_DOUBLE_PRECISION, &
404 0 : sendbl2(nk)%type, ierror)
405 0 : call MPI_TYPE_COMMIT(sendbl2(nk)%type, ierror)
406 :
407 0 : sendbl2(nk)%blocksizes(1) = DTWO*ijtot
408 0 : sendbl2(nk)%displacements(1) = 0
409 0 : sendbl2(nk)%partneroffset = (myid_z-nk/npr_y-1) * ijtot * DTWO
410 :
411 : else
412 :
413 0 : sendbl2(nk)%type = MPI_DATATYPE_NULL
414 :
415 0 : sendbl2(nk)%blocksizes(1) = 0
416 0 : sendbl2(nk)%displacements(1) = 0
417 0 : sendbl2(nk)%partneroffset = 0
418 :
419 : endif
420 0 : sendbl2(nk)%nparcels = size(sendbl2(nk)%displacements)
421 0 : sendbl2(nk)%tot_size = sum(sendbl2(nk)%blocksizes)
422 0 : max_nparcels = max(max_nparcels, sendbl2(nk)%nparcels)
423 :
424 0 : if ( (nk/npr_y) > myid_z .and. mod(nk,npr_y) == myid_y ) then
425 :
426 : call MPI_TYPE_INDEXED(ione, DTWO*ijtot, &
427 : nk/npr_y * ijtot * DTWO, MPI_DOUBLE_PRECISION, &
428 0 : recvbl2(nk)%type, ierror)
429 0 : call MPI_TYPE_COMMIT(recvbl2(nk)%type, ierror)
430 :
431 0 : recvbl2(nk)%blocksizes(1) = DTWO*ijtot
432 0 : recvbl2(nk)%displacements(1) = nk/npr_y * ijtot * DTWO
433 0 : recvbl2(nk)%partneroffset = 0
434 :
435 : else
436 :
437 0 : recvbl2(nk)%type = MPI_DATATYPE_NULL
438 :
439 0 : recvbl2(nk)%blocksizes(1) = 0
440 0 : recvbl2(nk)%displacements(1) = 0
441 0 : recvbl2(nk)%partneroffset = 0
442 :
443 : endif
444 0 : recvbl2(nk)%nparcels = size(recvbl2(nk)%displacements)
445 0 : recvbl2(nk)%tot_size = sum(recvbl2(nk)%blocksizes)
446 0 : max_nparcels = max(max_nparcels, recvbl2(nk)%nparcels)
447 : enddo
448 :
449 0 : call get_partneroffset(grid%commyz, sendbl1, recvbl1)
450 0 : call get_partneroffset(grid%commyz, sendbl2, recvbl2)
451 :
452 : endif
453 : endif
454 :
455 : #if (!defined PAREXCH)
456 0 : locinbuf = loc(pe16)
457 : #else
458 : locinbuf = loc(inbuf)
459 : #endif
460 0 : locoutbuf = loc(outbuf)
461 0 : ptr_inbuf8 = locinbuf
462 0 : ptr_outbuf8 = locoutbuf
463 : #endif
464 :
465 : ! Top down
466 :
467 : #if (DSIZE == 16)
468 : !$omp parallel do &
469 : !$omp default(shared) &
470 : !$omp private(i, j, k)
471 : do k = kfirst,klast
472 : do j = jfirst,jlast
473 : do i = ifirst,ilast
474 : delp16(i,j,k) = delp(i,j,k)
475 : enddo
476 : enddo
477 : enddo
478 : #endif
479 :
480 : !$omp parallel do &
481 : !$omp default(shared) &
482 : !$omp private(i, j)
483 0 : do j = jfirst,jlast
484 0 : do i = ifirst,ilast
485 0 : pe16(i,j,kfirst) = DP0_0
486 : enddo
487 : enddo
488 :
489 : ! compute partial sums
490 :
491 : !$omp parallel do &
492 : !$omp default(shared) &
493 : !$omp private(i, j, k)
494 0 : do j = jfirst,jlast
495 0 : do k = kfirst+1,klast+1
496 0 : do i = ifirst,ilast
497 : #if (DSIZE == 16)
498 : pe16(i,j,k) = pe16(i,j,k-1) + delp16(i,j,k-1)
499 : #else
500 0 : pe16(i,j,k) = pe16(i,j,k-1) + delp(i,j,k-1)
501 : #endif
502 : enddo
503 : enddo
504 : enddo
505 :
506 : #if defined( SPMD )
507 0 : if (npr_z .gt. 1) then
508 :
509 : ! communicate upward
510 :
511 : # if !defined (PAREXCH)
512 : call mp_sendirr(grid%commyz, sendbl1, recvbl1, inbuf8, outbuf8, &
513 0 : modc=grid%modc_cdcore )
514 : call mp_recvirr(grid%commyz, sendbl1, recvbl1, inbuf8, outbuf8, &
515 0 : modc=grid%modc_cdcore )
516 : # else
517 :
518 : do nk = 0, npr_z-1
519 : sendcount(nk) = 0
520 : recvcount(nk) = 0
521 : enddo
522 :
523 : !$omp parallel do &
524 : !$omp default(shared) &
525 : !$omp private(i, j, nk)
526 : do nk = myid_z+1, npr_z-1
527 : do j = jfirst,jlast
528 : do i = ifirst,ilast
529 : inbuf(i,j,nk-myid_z-1) = pe16(i,j,klast+1)
530 : enddo
531 : enddo
532 : ! Double sendcount since quantities are 16-bytes long
533 : sendcount(nk) = DTWO*ijtot
534 : enddo
535 :
536 : call parexchangevector(grid%comm_z, sendcount, inbuf8, recvcount, outbuf8)
537 :
538 : # endif
539 :
540 : !$omp parallel do &
541 : !$omp default(shared) &
542 : !$omp private(i, j, k, nk)
543 0 : do k = kfirst,klast+1
544 0 : do nk = 0, myid_z-1
545 0 : do j = jfirst,jlast
546 0 : do i = ifirst,ilast
547 0 : pe16(i,j,k) = pe16(i,j,k) + outbuf(i,j,nk)
548 : enddo
549 : enddo
550 : enddo
551 : enddo
552 :
553 : endif
554 : #endif
555 :
556 : !$omp parallel do &
557 : !$omp default(shared) &
558 : !$omp private(i, j, k)
559 0 : do k = kfirst,klast+1
560 0 : do j = jfirst,jlast
561 0 : do i = ifirst,ilast
562 0 : pe(i,k,j) = pe16(i,j,k) + ptop
563 0 : pk(i,j,k) = pe(i,k,j) ** akap
564 : enddo
565 : enddo
566 : enddo
567 :
568 : ! Bottom up
569 :
570 : !$omp parallel do &
571 : !$omp default(shared) &
572 : !$omp private(i, j, k)
573 0 : do k = kfirst,klast
574 0 : do j = jfirst,jlast
575 0 : do i = ifirst,ilast
576 0 : delp16(i,j,k) = cp*pt(i,j,k)*(pk(i,j,k+1)-pk(i,j,k))
577 : enddo
578 : enddo
579 : enddo
580 :
581 : !$omp parallel do &
582 : !$omp default(shared) &
583 : !$omp private(i, j)
584 0 : do j = jfirst,jlast
585 0 : do i = ifirst,ilast
586 0 : pe16(i,j,klast+1) = DP0_0
587 : enddo
588 : enddo
589 :
590 : ! compute partial sums
591 :
592 : !$omp parallel do &
593 : !$omp default(shared) &
594 : !$omp private(i, j, k)
595 0 : do j = jfirst,jlast
596 0 : do k = klast,kfirst,-1
597 0 : do i = ifirst,ilast
598 0 : pe16(i,j,k) = pe16(i,j,k+1) + delp16(i,j,k)
599 : enddo
600 : enddo
601 : enddo
602 :
603 : #if defined( SPMD )
604 0 : if (npr_z .gt. 1) then
605 :
606 : ! communicate downward
607 :
608 : # if !defined (PAREXCH)
609 : call mp_sendirr(grid%commyz, sendbl2, recvbl2, inbuf8, outbuf8, &
610 0 : modc=grid%modc_cdcore )
611 : call mp_recvirr(grid%commyz, sendbl2, recvbl2, inbuf8, outbuf8, &
612 0 : modc=grid%modc_cdcore )
613 : # else
614 :
615 : do nk = 0, npr_z-1
616 : sendcount(nk) = 0
617 : recvcount(nk) = 0
618 : enddo
619 :
620 : !$omp parallel do &
621 : !$omp default(shared) &
622 : !$omp private(i, j, nk)
623 : do nk = 0, myid_z-1
624 : do j = jfirst,jlast
625 : do i = ifirst,ilast
626 : inbuf(i,j,nk) = pe16(i,j,kfirst)
627 : enddo
628 : enddo
629 : ! Double sendcount since quantities are 16-bytes long
630 : sendcount(nk) = DTWO*ijtot
631 : enddo
632 :
633 : call parexchangevector(grid%comm_z, sendcount, inbuf8, recvcount, outbuf8)
634 :
635 : # endif
636 :
637 : !$omp parallel do &
638 : !$omp default(shared) &
639 : !$omp private(i, j, k, nk)
640 0 : do k = kfirst,klast+1
641 0 : do nk = myid_z+1, npr_z-1
642 0 : do j = jfirst,jlast
643 0 : do i = ifirst,ilast
644 : # if !defined (PAREXCH)
645 0 : pe16(i,j,k) = pe16(i,j,k) + outbuf(i,j,nk)
646 : # else
647 : pe16(i,j,k) = pe16(i,j,k) + outbuf(i,j,nk-myid_z-1)
648 : # endif
649 : enddo
650 : enddo
651 : enddo
652 : enddo
653 :
654 : endif
655 : #endif
656 :
657 : !$omp parallel do &
658 : !$omp default(shared) &
659 : !$omp private(i, j, k)
660 0 : do k = kfirst,klast+1
661 0 : do j = jfirst,jlast
662 0 : do i = ifirst,ilast
663 0 : wz(i,j,k) = pe16(i,j,k) + hs(i,j)
664 : enddo
665 : enddo
666 : enddo
667 :
668 0 : return
669 : ! endif for NO_CRAY_POINTERS
670 : #endif
671 : !EOC
672 : end subroutine geopk16
673 : !-----------------------------------------------------------------------
674 :
675 : !-----------------------------------------------------------------------
676 : !BOP
677 : ! !ROUTINE: geopk_d --- Calculate geopotential to the kappa
678 : !
679 : ! !INTERFACE:
680 0 : subroutine geopk_d( grid, pe, delp, pk, wz, hs, pt, ng, cp, akap )
681 :
682 : use shr_kind_mod, only : r8 => shr_kind_r8, i8 => shr_kind_i8
683 : use dynamics_vars, only : T_FVDYCORE_GRID
684 :
685 : implicit none
686 :
687 : #if defined ( SPMD )
688 : #include "mpif.h"
689 : #endif
690 :
691 : ! !INPUT PARAMETERS:
692 :
693 : type (T_FVDYCORE_GRID), intent(in) :: grid
694 : integer, intent(in) :: ng ! Halo size (not always = ng_d)
695 :
696 : real(r8) akap, cp
697 : real(r8) hs(1:grid%im,grid%jfirst:grid%jlast)
698 :
699 : ! !INPUT PARAMETERS:
700 : real(r8) pt(1:grid%im,grid%jfirst-ng:grid%jlast+ng,grid%kfirst:grid%klast)
701 : real(r8) delp(1:grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast)
702 :
703 : ! !OUTPUT PARAMETERS:
704 : real(r8) wz(1:grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast+1) ! space N*1 S*1
705 : real(r8) pk(1:grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast+1) ! space N*1 S*1
706 : real(r8) pe(1:grid%im,grid%kfirst:grid%klast+1,grid%jfirst:grid%jlast) ! temporary variable
707 :
708 : ! !DESCRIPTION:
709 : ! Calculates geopotential and pressure to the kappa. This is an expensive
710 : ! operation and several out arrays are kept around for future use.
711 : ! To preserve reproducibility, ordering of transposed-based geopk algorithm
712 : ! is preserved at the cost of a serialization of computation in the Z-direction.
713 : !
714 : ! !REVISION HISTORY:
715 : !
716 : ! PW 08.06.27: Original: simple ring r8 version of geopk16 -
717 : ! serialized Z-direction but minimized communication overhead
718 : !
719 : !EOP
720 : !---------------------------------------------------------------------
721 : !BOC
722 :
723 : ! Local:
724 : real(r8):: ptop
725 : integer :: km, ifirst, ilast, jfirst, jlast, kfirst, klast
726 : integer :: npr_z
727 :
728 : integer :: itot, jtot
729 :
730 : integer :: n_blocks
731 : logical :: sendd
732 :
733 : integer :: i, il, ilmax, ib, iblksiz
734 : integer :: j, jl, jlmax, jb, jblksiz
735 : integer :: k, block, ierror
736 :
737 0 : integer :: klimits(2), klimits_all(2,0:grid%npr_z-1)
738 : integer, save :: k_succ_pid, k_pred_pid
739 :
740 0 : integer, allocatable :: rcvreq(:), sndreq(:)
741 :
742 : real (r8), parameter :: DP0_0 = 0.0_r8
743 0 : real (r8) l_pe(1:grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast+1)
744 0 : real (r8) l_delp(1:grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast)
745 :
746 : real (r8) inbuf(1:grid%im,grid%jfirst:grid%jlast,0:grid%npr_z-1)
747 : real (r8) outbuf(1:grid%im,grid%jfirst:grid%jlast,0:grid%npr_z-1)
748 : integer sendcount(0:grid%npr_z-1), recvcount(0:grid%npr_z-1)
749 :
750 : integer first_time_through
751 : data first_time_through / 0 /
752 :
753 : #if defined ( SPMD )
754 : integer status (MPI_STATUS_SIZE) ! Status of message
755 : #endif
756 :
757 : !
758 : ! Initialize variables from Grid
759 : !
760 0 : ptop = grid%ptop
761 :
762 0 : km = grid%km
763 :
764 0 : ifirst = 1 ! 2004.10.04 (WS): Now hardwired for 1..im
765 0 : ilast = grid%im ! Code was always used in this mode
766 0 : jfirst = grid%jfirst
767 0 : jlast = grid%jlast
768 0 : kfirst = grid%kfirst
769 0 : klast = grid%klast
770 :
771 0 : npr_z = grid%npr_z
772 :
773 0 : itot = (ilast-ifirst+1)
774 0 : jtot = (jlast-jfirst+1)
775 :
776 0 : if (grid%modc_cdcore(3) .eq. 1) then
777 : sendd = .true.
778 : else
779 0 : sendd = .false.
780 : endif
781 :
782 0 : n_blocks = max(1,grid%geopkblocks)
783 :
784 0 : if (n_blocks < jtot) then
785 0 : jblksiz = ceiling(float(jtot)/float(n_blocks))
786 : iblksiz = itot
787 : else
788 0 : jblksiz = 1
789 0 : iblksiz = ceiling(float(itot*jtot)/float(n_blocks))
790 : endif
791 :
792 0 : block = 0
793 0 : do j=jfirst,jlast,jblksiz
794 0 : do i=ifirst,ilast,iblksiz
795 0 : block = block + 1
796 : enddo
797 : enddo
798 :
799 0 : allocate( sndreq(block) )
800 0 : allocate( rcvreq(block) )
801 :
802 0 : if (first_time_through .eq. 0) then
803 0 : first_time_through = 1
804 0 : k_pred_pid = -1
805 0 : k_succ_pid = -1
806 : #if defined (SPMD)
807 0 : klimits(1) = kfirst
808 0 : klimits(2) = klast
809 : call mpi_allgather (klimits, 2, mpi_integer, &
810 : klimits_all, 2, mpi_integer, &
811 0 : grid%comm_z, ierror)
812 0 : do i=0,npr_z-1
813 0 : if (klimits_all(2,i) == kfirst-1) k_pred_pid = i
814 0 : if (klimits_all(1,i) == klast+1) k_succ_pid = i
815 : enddo
816 : #endif
817 : endif
818 :
819 : ! Top down
820 :
821 : ! prepost first set of receive requests
822 : #if defined (SPMD)
823 0 : if (k_pred_pid /= -1) then
824 0 : block = 0
825 0 : do j=jfirst,jlast,jblksiz
826 0 : if (j+jblksiz > jlast) then
827 0 : jb = jlast-j+1
828 : else
829 : jb = jblksiz
830 : endif
831 :
832 0 : do i=ifirst,ilast,iblksiz
833 0 : if (i+iblksiz > ilast) then
834 0 : ib = ilast-i+1
835 : else
836 : ib = iblksiz
837 : endif
838 :
839 0 : block = block + 1
840 0 : call mpi_irecv (l_pe(i,j,kfirst), jb*ib, &
841 : mpi_real8, k_pred_pid, block, &
842 0 : grid%comm_z, rcvreq(block), ierror)
843 : enddo
844 : enddo
845 : endif
846 : #endif
847 :
848 0 : block = 0
849 0 : do j=jfirst,jlast,jblksiz
850 :
851 0 : if (j+jblksiz > jlast) then
852 0 : jb = jlast-j+1
853 : else
854 : jb = jblksiz
855 : endif
856 0 : jlmax = j+jb-1
857 :
858 0 : do i=ifirst,ilast,iblksiz
859 0 : if (i+iblksiz > ilast) then
860 0 : ib = ilast-i+1
861 : else
862 : ib = iblksiz
863 : endif
864 0 : ilmax = i+ib-1
865 :
866 0 : block = block + 1
867 :
868 : ! get data from k predecessor
869 0 : if (k_pred_pid /= -1) then
870 : #if defined (SPMD)
871 0 : call mpi_wait (rcvreq(block), status, ierror)
872 : #endif
873 : else
874 0 : do jl=j,jlmax
875 0 : do il = i,ilmax
876 0 : l_pe(il,jl,kfirst) = DP0_0
877 : enddo
878 : enddo
879 : endif
880 :
881 : ! compute partial sums (note that can not thread over k-loop)
882 0 : do k = kfirst+1,klast+1
883 0 : do jl=j,jlmax
884 0 : do il = i,ilmax
885 0 : l_pe(il,jl,k) = l_pe(il,jl,k-1) + delp(il,jl,k-1)
886 : enddo
887 : enddo
888 : enddo
889 :
890 : ! send results to k successor
891 : #if defined (SPMD)
892 0 : if (k_succ_pid /= -1) then
893 0 : if (sendd) then
894 0 : call mpi_send (l_pe(i,j,klast+1), jb*ib, mpi_real8, &
895 : k_succ_pid, block, grid%comm_z, &
896 0 : ierror)
897 : else
898 0 : call mpi_isend (l_pe(i,j,klast+1), jb*ib, mpi_real8, &
899 : k_succ_pid, block, grid%comm_z, &
900 0 : sndreq(block), ierror)
901 : endif
902 : endif
903 : #endif
904 : !$omp parallel do &
905 : !$omp default(shared) &
906 : !$omp private(il, jl, k)
907 0 : do k = kfirst,klast+1
908 0 : do jl = j,jlmax
909 0 : do il = i,ilmax
910 0 : pe(il,k,jl) = l_pe(il,jl,k) + ptop
911 0 : pk(il,jl,k) = pe(il,k,jl) ** akap
912 : enddo
913 : enddo
914 : enddo
915 :
916 : #if defined (SPMD)
917 0 : if (k_succ_pid /= -1) then
918 0 : if (.not. sendd) then
919 0 : call mpi_wait (sndreq(block), status, ierror)
920 : endif
921 : endif
922 : #endif
923 : enddo
924 : enddo
925 :
926 : ! Bottom up
927 :
928 : ! prepost second set of receive requests
929 : #if defined (SPMD)
930 0 : if (k_succ_pid /= -1) then
931 0 : block = 0
932 0 : do j=jfirst,jlast,jblksiz
933 0 : if (j+jblksiz > jlast) then
934 0 : jb = jlast-j+1
935 : else
936 : jb = jblksiz
937 : endif
938 :
939 0 : do i=ifirst,ilast,iblksiz
940 0 : if (i+iblksiz > ilast) then
941 0 : ib = ilast-i+1
942 : else
943 : ib = iblksiz
944 : endif
945 :
946 0 : block = block + 1
947 0 : call mpi_irecv (l_pe(i,j,klast+1), jb*ib, &
948 : mpi_real8, k_succ_pid, block, &
949 0 : grid%comm_z, rcvreq(block), ierror)
950 : enddo
951 : enddo
952 : endif
953 : #endif
954 :
955 0 : block = 0
956 0 : do j=jfirst,jlast,jblksiz
957 :
958 0 : if (j+jblksiz > jlast) then
959 0 : jb = jlast-j+1
960 : else
961 : jb = jblksiz
962 : endif
963 0 : jlmax = j+jb-1
964 :
965 0 : do i=ifirst,ilast,iblksiz
966 0 : if (i+iblksiz > ilast) then
967 0 : ib = ilast-i+1
968 : else
969 : ib = iblksiz
970 : endif
971 0 : ilmax = i+ib-1
972 :
973 0 : block = block + 1
974 :
975 : !$omp parallel do &
976 : !$omp default(shared) &
977 : !$omp private(il, jl, k)
978 0 : do k = kfirst,klast
979 0 : do jl=j,jlmax
980 0 : do il = i,ilmax
981 0 : l_delp(il,jl,k) = &
982 0 : cp*pt(il,jl,k)*(pk(il,jl,k+1)-pk(il,jl,k))
983 : enddo
984 : enddo
985 : enddo
986 :
987 : ! get data from k predecessor
988 0 : if (k_succ_pid /= -1) then
989 : #if defined (SPMD)
990 0 : call mpi_wait (rcvreq(block), status, ierror)
991 : #endif
992 : else
993 0 : do jl=j,jlmax
994 0 : do il = i,ilmax
995 0 : l_pe(il,jl,klast+1) = DP0_0
996 : enddo
997 : enddo
998 : endif
999 :
1000 : ! compute partial sums (note that can not thread over k-loop)
1001 0 : do k = klast,kfirst,-1
1002 0 : do jl=j,jlmax
1003 0 : do il = i,ilmax
1004 0 : l_pe(il,jl,k) = l_pe(il,jl,k+1) + l_delp(il,jl,k)
1005 : enddo
1006 : enddo
1007 : enddo
1008 :
1009 : ! send results to k predecessor
1010 : #if defined (SPMD)
1011 0 : if (k_pred_pid /= -1) then
1012 0 : if (sendd) then
1013 0 : call mpi_send (l_pe(i,j,kfirst), jb*ib, mpi_real8, &
1014 : k_pred_pid, block, &
1015 0 : grid%comm_z, ierror)
1016 : else
1017 0 : call mpi_isend (l_pe(i,j,kfirst), jb*ib, mpi_real8, &
1018 : k_pred_pid, block, &
1019 0 : grid%comm_z, sndreq(block), ierror)
1020 : endif
1021 : endif
1022 : #endif
1023 :
1024 : !$omp parallel do &
1025 : !$omp default(shared) &
1026 : !$omp private(il, jl, k)
1027 0 : do k = kfirst,klast+1
1028 0 : do jl=j,jlmax
1029 0 : do il = i,ilmax
1030 0 : wz(il,jl,k) = l_pe(il,jl,k) + hs(il,jl)
1031 : enddo
1032 : enddo
1033 : enddo
1034 :
1035 : #if defined (SPMD)
1036 0 : if (k_pred_pid /= -1) then
1037 0 : if (.not. sendd) then
1038 0 : call mpi_wait (sndreq(block), status, ierror)
1039 : endif
1040 : endif
1041 : #endif
1042 :
1043 : enddo
1044 :
1045 : enddo
1046 :
1047 0 : deallocate( sndreq )
1048 0 : deallocate( rcvreq )
1049 :
1050 0 : return
1051 : !EOC
1052 : end subroutine geopk_d
1053 : !-----------------------------------------------------------------------
|
