Line data Source code
1 : module edge_mod
2 :
3 : use shr_kind_mod, only: r8=>shr_kind_r8, i8=>shr_kind_i8
4 : use dimensions_mod, only: max_neigh_edges, nelemd
5 : use perf_mod, only: t_startf, t_stopf, t_adj_detailf ! _EXTERNAL
6 : use thread_mod, only: max_num_threads, omp_get_num_threads, omp_get_thread_num
7 : use coordinate_systems_mod, only: cartesian3D_t
8 : use schedtype_mod, only: cycle_t, schedule_t, pgindex_t, schedule, HME_Ordinal,HME_Cardinal
9 : use cam_abortutils, only: endrun
10 : use cam_logfile, only: iulog
11 : use parallel_mod, only: parallel_t, &
12 : MAX_ACTIVE_MSG, HME_status_size, BNDRY_TAG_BASE, HME_BNDRY_A2A, HME_BNDRY_P2P, &
13 : HME_BNDRY_A2AO
14 : use edgetype_mod, only: edgedescriptor_t, edgebuffer_t, &
15 : Longedgebuffer_t, initedgebuffer_callid, Ghostbuffer3D_t
16 : use element_mod, only: element_t
17 : use gbarrier_mod, only: gbarrier_init, gbarrier_delete
18 : use spmd_utils, only: mpi_real8, mpi_integer, mpi_info_null, mpi_success
19 :
20 : implicit none
21 : private
22 : save
23 :
24 : ! 8-byte Integer routines
25 : public :: LongEdgeVpack, LongEdgeVunpackMIN
26 :
27 : ! 8-byte Real routines
28 : public :: zeroEdgeBuffer
29 :
30 : interface initEdgeBuffer
31 : module procedure initEdgeBuffer_r8
32 : module procedure initEdgeBuffer_i8
33 : end interface
34 : interface initEdgeSBuffer
35 : module procedure initEdgeSbuffer_r8
36 : end interface
37 : interface freeEdgeBuffer
38 : module procedure freeEdgeBuffer_r8
39 : module procedure freeEdgeBuffer_i8
40 : end interface
41 : interface freeGhostBuffer
42 : module procedure freeGhostBuffer_r8
43 : end interface
44 :
45 : public :: initEdgeBuffer
46 : public :: initEdgeSBuffer
47 : public :: freeEdgeBuffer
48 :
49 : public :: initGhostBuffer
50 : public :: ghostpack, ghostunpack
51 : public :: freeGhostBuffer
52 :
53 : !---------------------------------------------------------
54 : ! Pack/unpack routines that use the New format Edge buffer
55 : !---------------------------------------------------------
56 :
57 : public :: edgeVpack, edgeVunpack
58 : public :: edgeVunpackMIN, edgeVunpackMAX
59 : public :: edgeDGVpack, edgeDGVunpack
60 : public :: edgeVunpackVert
61 :
62 :
63 : public :: initGhostBuffer3D
64 : public :: FreeGhostBuffer3D
65 : public :: ghostVpack3D, ghostVunpack3D
66 :
67 : !----------------------------------------------------------------
68 : ! Pack/unpack routines that communicate a fixed number values
69 : ! per element. This is used to communicate MIN/MAX values from
70 : ! neighboring elemeents
71 : !----------------------------------------------------------------
72 : interface edgeSpack
73 : module procedure edgeSpack_r8
74 : end interface
75 : public :: edgeSpack
76 : public :: edgeSunpackMIN, edgeSunpackMAX
77 :
78 : logical, private :: threadsafe=.true.
79 :
80 : real(kind=r8), parameter, public :: edgeDefaultVal = 1.11e+100_r8
81 :
82 : ! NOTE ON ELEMENT ORIENTATION
83 : !
84 : ! Element orientation: index V(i,j)
85 : !
86 : ! (1,np) NWEST (np,np) NEAST
87 : !
88 : ! (1,1) SWEST (np,1) SEAST
89 : !
90 : !
91 : ! for the edge neighbors:
92 : ! we set the "reverse" flag if two elements who share an edge use a
93 : ! reverse orientation. The data is reversed during the *pack* stage
94 : ! For corner neighbors:
95 : ! for edge buffers, there is no orientation because two corner neighbors
96 : ! only share a single point.
97 : ! For ghost cell data, there is a again two posible orientations. For
98 : ! this case, we set the "reverse" flag if the corner element is using
99 : ! the reverse orientation. In this case, the data is reversed during the
100 : ! *unpack* stage (not sure why)
101 : !
102 : ! The edge orientation is set at startup. The corner orientation is computed
103 : ! at run time, via the call to compute_ghost_corner_orientation()
104 : ! This routine only works for meshes with at most 1 corner element. It's
105 : ! not called and the corner orientation flag is not set for unstructured meshes
106 :
107 : !
108 : ! Christoph Erath
109 : ! pack/unpack partial element of data of size (nx,nx) with user specifed halo size nh
110 : ! user specifies the sizes when creating the buffer
111 : ! buffer has 1 extra dimension (as compared to subroutines above) for multiple tracers
112 : ! input/output arrays are cartesian, and thus assume at most 1 element at each corner
113 : ! hence currently only supports cube-sphere grids.
114 : !
115 : !
116 : ! routines which including element edge data
117 : ! (used for FVM arrays where edge data is not shared by neighboring elements)
118 : ! these routines pack/unpack element data with user specified halo size
119 :
120 : ! Wrap pointer so we can make an array of them.
121 : type :: wrap_ptr
122 : real (kind=r8), dimension(:,:), pointer :: ptr => null()
123 : end type wrap_ptr
124 :
125 : type(wrap_ptr) :: edgebuff_ptrs(0:1)
126 :
127 : contains
128 :
129 0 : subroutine initEdgeSBuffer_r8(par,edge,elem,nlyr,bndry_type, nthreads)
130 : type (parallel_t), intent(in) :: par
131 : type (EdgeBuffer_t), target, intent(out) :: edge
132 : type (element_t), intent(in) :: elem(:)
133 : integer, intent(in) :: nlyr
134 : integer , optional, intent(in) :: bndry_type
135 : integer, optional, intent(in) :: nthreads
136 :
137 :
138 : call initEdgeBuffer(par,edge,elem,nlyr,bndry_type=bndry_type, &
139 0 : nthreads=nthreads,CardinalLength=1,OrdinalLength=1)
140 :
141 0 : end subroutine initEdgeSBuffer_r8
142 :
143 18432 : subroutine initGhostBuffer(par,edge,elem,nlyr,ndepth, npoints,bndry_type,nthreads)
144 :
145 : type (parallel_t), intent(in) :: par
146 : type (Edgebuffer_t), target, intent(out) :: edge
147 : type (element_t), intent(in) :: elem(:)
148 : integer,intent(in) :: nlyr,ndepth, npoints
149 : integer , optional, intent(in) :: bndry_type
150 : integer, optional, intent(in) :: nthreads
151 :
152 : call initEdgeBuffer(par,edge,elem,nlyr,bndry_type=bndry_type, &
153 18432 : nthreads=nthreads,CardinalLength=ndepth*npoints,OrdinalLength=ndepth*ndepth)
154 : ! set some parameters need to support deep halos
155 18432 : edge%ndepth = ndepth
156 18432 : edge%npoints = npoints
157 18432 : edge%lb = 1 - edge%ndepth
158 18432 : edge%ub = edge%npoints + edge%ndepth
159 :
160 18432 : end subroutine initGhostBuffer
161 :
162 :
163 :
164 0 : subroutine zeroEdgeBuffer(edge)
165 :
166 : type (EdgeBuffer_t), intent(inout) :: edge
167 : integer :: i
168 :
169 0 : do i=1,edge%nbuf
170 0 : edge%buf(i) = 0.0d0
171 0 : edge%receive(i) = 0.0d0
172 : enddo
173 :
174 0 : end subroutine zeroEdgeBuffer
175 :
176 : ! =========================================
177 : ! initEdgeBuffer:
178 : !
179 : ! create an Real based communication buffer
180 : ! =========================================
181 62976 : subroutine initEdgeBuffer_r8(par,edge,elem,nlyr, bndry_type,nthreads,CardinalLength,OrdinalLength)
182 : use dimensions_mod, only: np, nelemd, max_corner_elem
183 : use schedtype_mod, only: cycle_t, schedule_t, schedule
184 : use mpi, only: MPI_VERSION
185 :
186 : type (parallel_t), intent(in) :: par
187 : type (EdgeBuffer_t), target, intent(out) :: edge
188 : type (element_t), intent(in) :: elem(:)
189 : integer, intent(in) :: nlyr
190 : integer, optional, intent(in) :: bndry_type
191 : integer, optional, intent(in) :: nthreads
192 : integer, optional, intent(in) :: CardinalLength
193 : integer, optional, intent(in) :: OrdinalLength
194 :
195 : ! Notes about the buf_ptr/receive_ptr options:
196 : !
197 : ! If an EdgeBuffer_t object is initialized from pre-existing storage
198 : ! (i.e. buf_ptr is provided and not null), it must *not* be freed,
199 : ! and must not be used if the underlying storage has been deallocated.
200 : !
201 : ! All these restrictions also applied to the old newbuf and newreceive
202 : ! options.
203 :
204 : ! Workaround for NAG bug.
205 : ! NAG 5.3.1 dies if you use pointer bounds remapping to set
206 : ! a pointer that is also a component. So remap to temporary,
207 : ! then use that to set component pointer.
208 :
209 : ! Local variables
210 : integer :: nbuf,ith
211 : integer :: nSendCycles, nRecvCycles
212 : integer :: icycle, ierr
213 : integer :: ie, i
214 : integer :: edgeid,elemid
215 : integer :: ptr,llen,moveLength, mLen, tlen
216 : type (Cycle_t), pointer :: pCycle
217 : type (Schedule_t), pointer :: pSchedule
218 : integer :: dest, source, length, tag, iptr
219 : integer :: nlen, ithr
220 :
221 : integer :: len, lenP,lenS
222 : integer :: j,jj,il,mesgid, dst0,src0
223 : integer :: moveptr
224 : integer :: nbuf2,ilm1,iem1,lenm1
225 31488 : integer,allocatable :: putmap2(:,:),getmap2(:,:)
226 : integer,allocatable :: scounts(:), rcounts(:)
227 : integer,allocatable :: sdispls(:), rdispls(:)
228 : integer :: nInter, nIntra
229 : integer :: icInter, icIntra
230 :
231 : integer :: maxnsend
232 : integer :: tmpnMesg
233 : integer :: wintmpnMesg, wintmpDest, wintmpDisp
234 : integer(kind=i8) :: winsize
235 : integer :: win
236 : integer :: sizeofreal
237 : integer, allocatable :: tmpDest(:),tmpDisp(:)
238 : integer :: nFull
239 : integer :: disp, one
240 : integer :: errorcode,errorlen
241 : integer :: CardinalLen, OrdinalLen
242 : character(len=80) :: errorstring
243 : character(len=80), parameter :: subname='initedgeBuffer'
244 :
245 31488 : if(present(bndry_type)) then
246 6144 : if ( MPI_VERSION >= 3 ) then
247 6144 : edge%bndry_type = bndry_type
248 : else
249 : edge%bndry_type = HME_BNDRY_P2P
250 : endif
251 : else
252 25344 : edge%bndry_type = HME_BNDRY_P2P
253 : endif
254 :
255 : ! Set the length of the cardinal and ordinal message lengths
256 31488 : if(present(CardinalLength)) then
257 18432 : CardinalLen = CardinalLength
258 : else
259 13056 : CardinalLen = np
260 : endif
261 31488 : if(present(OrdinalLength)) then
262 18432 : OrdinalLen = OrdinalLength
263 : else
264 13056 : OrdinalLen = 1
265 : endif
266 :
267 : ! DO NOT REMOVE THIS NEXT BARRIER
268 : ! MT: This initial barrier fixes a long standing issue with Intel compilers on
269 : ! two different platforms. Without this barrier, edge buffers initialized from
270 : ! within the threaded region would not work in a reproducable way with certain
271 : ! thread combinations. I cant explain why, but this fixes that issue on Edison
272 : !$OMP BARRIER
273 :
274 31488 : if (nlyr==0) return ! tracer code might call initedgebuffer() with zero tracers
275 :
276 :
277 : !$OMP MASTER
278 : !
279 : ! Keep a counter of how many times initedgebuffer is called.
280 : ! This is used to assign a unique message ID for the boundary exchange
281 : !
282 31488 : initedgebuffer_callid=initedgebuffer_callid+1
283 31488 : edge%id = initedgebuffer_callid
284 31488 : edge%tag = BNDRY_TAG_BASE + MODULO(edge%id, MAX_ACTIVE_MSG)
285 :
286 125952 : allocate(edge%putmap(max_neigh_edges,nelemd))
287 94464 : allocate(edge%getmap(max_neigh_edges,nelemd))
288 94464 : allocate(edge%reverse(max_neigh_edges,nelemd))
289 :
290 2024088 : edge%putmap(:,:)=-1
291 2024088 : edge%getmap(:,:)=-1
292 :
293 125952 : allocate(putmap2(max_neigh_edges,nelemd))
294 94464 : allocate(getmap2(max_neigh_edges,nelemd))
295 2024088 : putmap2(:,:)=-1
296 2024088 : getmap2(:,:)=-1
297 252888 : do ie=1,nelemd
298 2024088 : do i=1,max_neigh_edges
299 1992600 : edge%reverse(i,ie) = elem(ie)%desc%reverse(i)
300 : enddo
301 : enddo
302 :
303 31488 : pSchedule => Schedule(1)
304 31488 : nSendCycles = pSchedule%nSendCycles
305 31488 : nRecvCycles = pSchedule%nRecvCycles
306 31488 : nInter = pSchedule%nInter
307 31488 : nIntra = pSchedule%nIntra
308 31488 : nFull = nInter+nIntra
309 :
310 31488 : edge%nInter=nInter
311 31488 : edge%nIntra=nIntra
312 :
313 31488 : if(nInter>0) then
314 125952 : allocate(edge%rcountsInter(nInter),edge%rdisplsInter(nInter))
315 94464 : allocate(edge%scountsInter(nInter),edge%sdisplsInter(nInter))
316 : endif
317 31488 : if(nIntra>0) then
318 0 : allocate(edge%rcountsIntra(nIntra),edge%rdisplsIntra(nIntra))
319 0 : allocate(edge%scountsIntra(nIntra),edge%sdisplsIntra(nIntra))
320 : endif
321 :
322 31488 : if (nSendCycles>0) then
323 125952 : allocate(edge%scountsFull(nSendCycles),edge%sdisplsFull(nSendCycles))
324 62976 : allocate(edge%Srequest(nSendCycles))
325 234274 : edge%scountsFull(:) = 0
326 : endif
327 : !
328 : ! Setup the data-structures for the sends
329 : !
330 31488 : j = 1
331 31488 : icycle = 1
332 31488 : dst0 = pSchedule%pIndx(j)%mesgid
333 31488 : il = pSchedule%pIndx(j)%edgeid
334 31488 : ie = pSchedule%pIndx(j)%elemid
335 31488 : len = CalcSegmentLength(pSchedule%pIndx(j),CardinalLen,OrdinalLen,nlyr)
336 31488 : edge%putmap(il,ie) = 0
337 31488 : if(nSendCycles>0) then
338 31488 : edge%sdisplsFull(icycle) = edge%putmap(il,ie)
339 31488 : edge%scountsFull(icycle) = len
340 : endif
341 31488 : ilm1 = il
342 31488 : iem1 = ie
343 31488 : lenm1 = len
344 :
345 1771200 : do j=2,SIZE(pSchedule%pIndx)
346 1739712 : il = pSchedule%pIndx(j)%edgeid
347 1739712 : ie = pSchedule%pIndx(j)%elemid
348 1739712 : mesgid = pSchedule%pIndx(j)%mesgid
349 1771200 : if(il>0 .and. ie >0) then
350 1738728 : len = CalcSegmentLength(pSchedule%pIndx(j),CardinalLen,OrdinalLen,nlyr)
351 1738728 : edge%putmap(il,ie) = edge%putmap(ilm1,iem1)+lenm1
352 1738728 : if(mesgid .ne. par%rank) then ! don't enter if this is a move cycle where (mesgid == par%rank)
353 948740 : if(mesgid .ne. dst0) then
354 171298 : icycle=icycle+1
355 171298 : if (nSendCycles>0) edge%sdisplsFull(icycle) = edge%putmap(il,ie)
356 : dst0=mesgid
357 : endif
358 948740 : if (nSendCycles>0) edge%scountsFull(icycle) = edge%scountsFull(icycle)+len
359 : endif
360 : ilm1=il
361 : iem1=ie
362 : lenm1=len
363 : endif
364 : enddo
365 :
366 : icInter=0
367 : icIntra=0
368 234274 : do icycle=1,nSendCycles
369 234274 : if(pSchedule%SendCycle(icycle)%onNode .eqv. .FALSE.) then
370 202786 : icInter=icInter+1
371 202786 : edge%sdisplsInter(icInter)=edge%sdisplsFull(icycle)
372 202786 : edge%scountsInter(icInter)=edge%scountsFull(icycle)
373 : else
374 0 : icIntra=icIntra+1
375 0 : edge%sdisplsIntra(icIntra)=edge%sdisplsFull(icycle)
376 0 : edge%scountsIntra(icIntra)=edge%scountsFull(icycle)
377 : endif
378 : enddo
379 :
380 31488 : if (nRecvCycles>0) then
381 125952 : allocate(edge%rcountsFull(nRecvCycles),edge%rdisplsFull(nRecvCycles))
382 94464 : allocate(edge%getDisplsFull(nRecvCycles),edge%putDisplsFull(nRecvCycles))
383 234274 : edge%rcountsFull(:) = 0
384 : ! allocate the MPI Send/Recv request handles
385 62976 : allocate(edge%Rrequest(nRecvCycles))
386 94464 : allocate(edge%status(HME_status_size,nRecvCycles))
387 : endif
388 :
389 : !
390 : ! Setup the data-structures for the receives
391 : !
392 31488 : j = 1
393 31488 : icycle = 1
394 31488 : src0 = pSchedule%gIndx(j)%mesgid
395 31488 : il = pSchedule%gIndx(j)%edgeid
396 31488 : ie = pSchedule%gIndx(j)%elemid
397 31488 : len = CalcSegmentLength(pSchedule%gIndx(j),CardinalLen,OrdinalLen,nlyr)
398 31488 : edge%getmap(il,ie) = 0
399 31488 : if (nRecvCycles>0) then
400 31488 : edge%rdisplsFull(icycle) = edge%getmap(il,ie)
401 31488 : edge%rcountsFull(icycle) = len
402 : endif
403 31488 : ilm1=il
404 31488 : iem1=ie
405 31488 : lenm1=len
406 :
407 1771200 : do j=2,SIZE(pSchedule%gIndx)
408 1739712 : il = pSchedule%gIndx(j)%edgeid
409 1739712 : ie = pSchedule%gIndx(j)%elemid
410 1739712 : mesgid = pSchedule%gIndx(j)%mesgid
411 1771200 : if(il>0 .and. ie >0) then
412 1738728 : len = CalcSegmentLength(pSchedule%gIndx(j),CardinalLen,OrdinalLen,nlyr)
413 1738728 : edge%getmap(il,ie) = edge%getmap(ilm1,iem1)+lenm1
414 1738728 : if(mesgid .ne. par%rank) then ! don't enter if this is a move cycle where (mesgid == par%rank)
415 948740 : if(mesgid .ne. src0) then
416 171298 : if (nRecvCycles>0) edge%rdisplsFull(icycle+1) = edge%getmap(il,ie)
417 171298 : icycle=icycle+1
418 171298 : src0=mesgid
419 : endif
420 948740 : if (nRecvCycles>0) edge%rcountsFull(icycle) = edge%rcountsFull(icycle)+len
421 : endif
422 : ilm1=il
423 : iem1=ie
424 : lenm1=len
425 : endif
426 : enddo
427 :
428 :
429 : !
430 : ! populate the Inter and Intra node communication data-structures
431 : !
432 : icInter=0
433 : icIntra=0
434 234274 : do icycle=1,nRecvCycles
435 234274 : if(pSchedule%RecvCycle(icycle)%onNode .eqv. .FALSE.) then
436 202786 : icInter=icInter+1
437 202786 : edge%rdisplsInter(icInter)=edge%rdisplsFull(icycle)
438 202786 : edge%rcountsInter(icInter)=edge%rcountsFull(icycle)
439 : else
440 0 : icIntra=icIntra+1
441 0 : edge%rdisplsIntra(icIntra)=edge%rdisplsFull(icycle)
442 0 : edge%rcountsIntra(icIntra)=edge%rcountsFull(icycle)
443 : endif
444 : enddo
445 :
446 :
447 : ! Setup the data-structures for the on process moves
448 : ! Note that this assumes that the data to move is at
449 : ! the end of the message buffer.
450 31488 : if(nRecvCycles>0) then
451 31488 : moveptr = edge%rdisplsFull(nRecvCycles)+edge%rcountsFull(nRecvCycles)+1
452 : else
453 : moveptr = 1
454 : endif
455 31488 : moveLength = 0
456 1802688 : do j=1,SIZE(pSchedule%gIndx)
457 1771200 : il = pSchedule%gIndx(j)%edgeid
458 1771200 : ie = pSchedule%gIndx(j)%elemid
459 1771200 : mesgid = pSchedule%gIndx(j)%mesgid
460 1802688 : if(mesgid == par%rank) then
461 789988 : len = CalcSegmentLength(pSchedule%gIndx(j),CardinalLen,OrdinalLen,nlyr)
462 789988 : moveLength = moveLength + len
463 : endif
464 : enddo
465 :
466 : ! decompose the move data between the available threads
467 31488 : if(max_num_threads<=0) then
468 0 : nlen = 1
469 : else
470 31488 : if(present(nthreads)) then
471 26112 : if (nthreads > 0) then
472 26112 : nlen = nthreads
473 : else
474 0 : nlen = max_num_threads
475 : end if
476 : else
477 5376 : nlen = 1
478 : end if
479 : end if
480 31488 : call gbarrier_init(edge%gbarrier, nlen)
481 :
482 94464 : allocate(edge%moveLength(nlen))
483 62976 : allocate(edge%movePtr(nlen))
484 :
485 31488 : if (nlen > 1) then
486 : ! the master thread performs no data movement because it is busy with the
487 : ! MPI messaging
488 0 : edge%moveLength(1) = -1
489 0 : edge%movePtr(1) = 0
490 :
491 : ! Calculate the length of the local copy in bndy_exchange
492 0 : llen = ceiling(real(moveLength,kind=r8)/real(nlen-1,kind=r8))
493 0 : iptr = moveptr
494 0 : mLen = 0
495 0 : do i=2,nlen
496 0 : if( (mLen+llen) <= moveLength) then
497 : tlen = llen
498 : else
499 0 : tlen = moveLength - mLen
500 : endif
501 0 : edge%moveLength(i) = tlen
502 0 : edge%movePtr(i) = iptr
503 0 : iptr = iptr + tlen
504 0 : mLen = mLen + tLen
505 : enddo
506 : else
507 31488 : edge%moveLength(1) = moveLength
508 31488 : edge%movePtr(1) = moveptr
509 : endif
510 :
511 : ! Set the maximum length of the message buffer
512 31488 : nbuf = movePtr+moveLength
513 :
514 31488 : edge%nlyr=nlyr
515 31488 : edge%nbuf=nbuf
516 :
517 94464 : allocate(edge%receive(nbuf))
518 62976 : allocate(edge%buf(nbuf))
519 :
520 : 21 format('RANK: ',i2, A,8(i6))
521 :
522 : !$OMP END MASTER
523 : ! MT: This next barrier is also needed - threads cannot start using edge()
524 : ! until MASTER is done initializing it
525 : !$OMP BARRIER
526 :
527 31488 : end subroutine initEdgeBuffer_r8
528 :
529 4330420 : integer function CalcSegmentLength(pgIndx,CardinalLength,OrdinalLength,nlyr) result(len)
530 :
531 : type(pgindex_t) :: pgIndx
532 : integer, intent(in) :: CardinalLength,OrdinalLength
533 : integer, intent(in) :: nlyr
534 :
535 : integer :: rem
536 : integer, parameter :: alignment=1 ! align on word boundaries
537 : ! integer, parameter :: alignment=2 ! align on 2 word boundaries
538 : ! integer, parameter :: alignment=8 ! align on 8 word boundaries
539 :
540 6589110 : select case(pgIndx%edgeType)
541 : CASE(HME_Cardinal)
542 2258690 : len = nlyr*CardinalLength
543 : CASE(HME_Ordinal)
544 4330420 : len = nlyr*OrdinalLength
545 : end select
546 :
547 4330420 : rem = MODULO(len,alignment)
548 : if(rem .ne. 0) then
549 : len = len + (alignment-rem)
550 : endif
551 :
552 4330420 : end function calcSegmentLength
553 :
554 : ! =========================================
555 : ! initEdgeBuffer:
556 : !
557 : ! create an Integer based communication buffer
558 : ! =========================================
559 1536 : subroutine initEdgeBuffer_i8(edge,nlyr)
560 : use dimensions_mod, only : np, nelemd, max_corner_elem
561 :
562 : integer, intent(in) :: nlyr
563 : type (LongEdgeBuffer_t), intent(out) :: edge
564 :
565 : ! Local variables
566 : integer :: nbuf
567 :
568 : ! sanity check for threading
569 1536 : if (omp_get_num_threads()>1) then
570 0 : call endrun('ERROR: initEdgeBuffer must be called before threaded reagion')
571 : endif
572 :
573 1536 : nbuf=4*(np+max_corner_elem)*nelemd
574 1536 : edge%nlyr=nlyr
575 1536 : edge%nbuf=nbuf
576 6144 : allocate(edge%buf(nlyr,nbuf))
577 433536 : edge%buf(:,:)=0
578 :
579 4608 : allocate(edge%receive(nlyr,nbuf))
580 433536 : edge%receive(:,:)=0
581 :
582 1536 : end subroutine initEdgeBuffer_i8
583 : ! =========================================
584 : ! edgeDGVpack:
585 : !
586 : ! Pack edges of v into buf for DG stencil
587 : ! =========================================
588 0 : subroutine edgeDGVpack(edge,v,vlyr,kptr,ielem)
589 : use dimensions_mod, only: np
590 :
591 : type (EdgeBuffer_t) :: edge
592 : integer, intent(in) :: vlyr
593 : real (kind=r8), intent(in) :: v(np,np,vlyr)
594 : integer, intent(in) :: kptr
595 : integer, intent(in) :: ielem
596 :
597 : ! =========================================
598 : ! This code is just a wrapper call the
599 : ! normal oldedgeVpack
600 : ! =========================================
601 0 : call edgeVpack(edge,v,vlyr,kptr,ielem)
602 :
603 0 : end subroutine edgeDGVpack
604 :
605 4608 : subroutine FreeGhostBuffer_r8(edge)
606 : type (EdgeBuffer_t), intent(inout) :: edge
607 4608 : call FreeEdgeBuffer_r8(edge)
608 4608 : end subroutine FreeGhostBuffer_r8
609 : ! ===========================================
610 : ! FreeEdgeBuffer:
611 : !
612 : ! Freed an edge communication buffer
613 : ! =========================================
614 8448 : subroutine FreeEdgeBuffer_r8(edge)
615 :
616 : type (EdgeBuffer_t),intent(inout) :: edge
617 :
618 : !$OMP BARRIER
619 : !$OMP MASTER
620 8448 : if(allocated(edge%buf)) deallocate(edge%buf)
621 8448 : if(allocated(edge%receive)) deallocate(edge%receive)
622 8448 : if(associated(edge%putmap)) deallocate(edge%putmap)
623 8448 : if(associated(edge%getmap)) deallocate(edge%getmap)
624 8448 : if(associated(edge%reverse)) deallocate(edge%reverse)
625 8448 : if(associated(edge%moveLength)) deallocate(edge%moveLength)
626 8448 : if(associated(edge%movePtr)) deallocate(edge%movePtr)
627 :
628 : ! All MPI communications
629 8448 : if(associated(edge%rcountsFull)) deallocate(edge%rcountsFull)
630 8448 : if(associated(edge%scountsFull)) deallocate(edge%scountsFull)
631 8448 : if(associated(edge%sdisplsFull)) deallocate(edge%sdisplsFull)
632 8448 : if(associated(edge%rdisplsFull)) deallocate(edge%rdisplsFull)
633 :
634 : ! Intra-node MPI Communication
635 8448 : if(edge%nIntra>0) then
636 0 : if(associated(edge%rcountsIntra)) deallocate(edge%rcountsIntra)
637 0 : if(associated(edge%scountsIntra)) deallocate(edge%scountsIntra)
638 0 : if(associated(edge%sdisplsIntra)) deallocate(edge%sdisplsIntra)
639 0 : if(associated(edge%rdisplsIntra)) deallocate(edge%rdisplsIntra)
640 : endif
641 :
642 : ! Inter-node MPI Communication
643 8448 : if(edge%nInter>0) then
644 8448 : if(associated(edge%rcountsInter)) deallocate(edge%rcountsInter)
645 8448 : if(associated(edge%scountsInter)) deallocate(edge%scountsInter)
646 8448 : if(associated(edge%sdisplsInter)) deallocate(edge%sdisplsInter)
647 8448 : if(associated(edge%rdisplsInter)) deallocate(edge%rdisplsInter)
648 : endif
649 8448 : if(allocated(edge%rRequest)) deallocate(edge%rRequest)
650 8448 : if(allocated(edge%sRequest)) deallocate(edge%sRequest)
651 8448 : if(allocated(edge%status)) deallocate(edge%status)
652 8448 : call gbarrier_delete(edge%gbarrier)
653 :
654 : !$OMP END MASTER
655 :
656 8448 : end subroutine FreeEdgeBuffer_r8
657 :
658 : ! ===========================================
659 : ! FreeEdgeBuffer:
660 : !
661 : ! Freed an edge communication buffer
662 : ! =========================================
663 1536 : subroutine FreeEdgeBuffer_i8(edge)
664 :
665 : type (LongEdgeBuffer_t),intent(inout) :: edge
666 :
667 1536 : edge%nbuf=0
668 1536 : edge%nlyr=0
669 1536 : deallocate(edge%buf)
670 1536 : deallocate(edge%receive)
671 :
672 1536 : end subroutine FreeEdgeBuffer_i8
673 :
674 : ! =========================================
675 : !
676 : !> @brief Pack edges of v into an edge buffer for boundary exchange.
677 : !
678 : !> This subroutine packs for one or more vertical layers into an edge
679 : !! buffer. If the buffer associated with edge is not large enough to
680 : !! hold all vertical layers you intent to pack, the method will
681 : !! halt the program with a call to endrum().
682 : !! @param[in] edge Edge Buffer into which the data will be packed.
683 : !! This buffer must be previously allocated with initEdgeBuffer().
684 : !! @param[in] v The data to be packed.
685 : !! @param[in] vlyr Number of vertical level coming into the subroutine
686 : !! for packing for input v.
687 : !! @param[in] kptr Vertical pointer to the place in the edge buffer where
688 : !! data will be located.
689 : ! =========================================
690 823548600 : subroutine edgeVpack(edge,v,vlyr,kptr,ielem)
691 : use dimensions_mod, only: np, max_corner_elem
692 : use control_mod, only: north, south, east, west, neast, nwest, seast, swest
693 :
694 : type (EdgeBuffer_t) :: edge
695 : integer, intent(in) :: vlyr
696 : real (kind=r8), intent(in) :: v(np,np,vlyr)
697 : integer, intent(in) :: kptr
698 : integer, intent(in) :: ielem
699 :
700 : ! Local variables
701 : integer :: i,k,ir,ll,iptr
702 : integer :: is,ie,in,iw,edgeptr
703 :
704 823548600 : is = edge%putmap(south,ielem)
705 823548600 : ie = edge%putmap(east,ielem)
706 823548600 : in = edge%putmap(north,ielem)
707 823548600 : iw = edge%putmap(west,ielem)
708 823548600 : if (edge%nlyr < (kptr+vlyr) ) then
709 0 : print *,'edge%nlyr = ',edge%nlyr
710 0 : print *,'kptr+vlyr = ',kptr+vlyr
711 0 : call endrun('edgeVpack: Buffer overflow: size of the vertical dimension must be increased!')
712 : endif
713 :
714 : !dir$ ivdep
715 68500355400 : do k=1,vlyr
716 67676806800 : iptr = np*(kptr+k-1)
717 >33920*10^7 : do i=1,np
718 >27070*10^7 : edge%buf(iptr+ie+i) = v(np ,i ,k) ! East
719 >27070*10^7 : edge%buf(iptr+is+i) = v(i ,1 ,k) ! South
720 >27070*10^7 : edge%buf(iptr+in+i) = v(i ,np,k) ! North
721 >33838*10^7 : edge%buf(iptr+iw+i) = v(1 ,i ,k) ! West
722 : enddo
723 : enddo
724 :
725 : ! This is really kludgy way to setup the index reversals
726 : ! But since it is so a rare event not real need to spend time optimizing
727 :
728 823548600 : if(edge%reverse(south,ielem)) then
729 1141672590 : do k=1,vlyr
730 1127946780 : iptr = np*(kptr+k-1)+is
731 5653459710 : do i=1,np
732 4511787120 : ir = np-i+1
733 5639733900 : edge%buf(iptr+ir)=v(i,1,k)
734 : enddo
735 : enddo
736 : endif
737 :
738 823548600 : if(edge%reverse(east,ielem)) then
739 380557530 : do k=1,vlyr
740 375982260 : iptr=np*(kptr+k-1)+ie
741 1884486570 : do i=1,np
742 1503929040 : ir = np-i+1
743 1879911300 : edge%buf(iptr+ir)=v(np,i,k)
744 : enddo
745 : enddo
746 : endif
747 :
748 823548600 : if(edge%reverse(north,ielem)) then
749 1141672590 : do k=1,vlyr
750 1127946780 : iptr=np*(kptr+k-1)+in
751 5653459710 : do i=1,np
752 4511787120 : ir = np-i+1
753 5639733900 : edge%buf(iptr+ir)=v(i,np,k)
754 : enddo
755 : enddo
756 : endif
757 :
758 823548600 : if(edge%reverse(west,ielem)) then
759 380557530 : do k=1,vlyr
760 375982260 : iptr=np*(kptr+k-1)+iw
761 1884486570 : do i=1,np
762 1503929040 : ir = np-i+1
763 1879911300 : edge%buf(iptr+ir)=v(1,i,k)
764 : enddo
765 : enddo
766 : endif
767 :
768 : ! SWEST
769 1647097200 : do ll=swest,swest+max_corner_elem-1
770 1647097200 : if (edge%putmap(ll,ielem) /= -1) then
771 822633546 : edgeptr = edge%putmap(ll,ielem)+1
772 68424243894 : do k=1,vlyr
773 67601610348 : iptr = (kptr+k-1)+edgeptr
774 67601610348 : if (iptr > size(edge%buf)) then
775 0 : write(6, *) 'ERROR SW: ',size(edge%buf),iptr,edge%putmap(ll,ielem)
776 0 : call endrun('pointer bounds ERROR SW')
777 : end if
778 68424243894 : edge%buf(iptr) = v(1, 1, k)
779 : end do
780 : end if
781 : end do
782 :
783 : ! SEAST
784 1647097200 : do ll=swest+max_corner_elem,swest+2*max_corner_elem-1
785 1647097200 : if (edge%putmap(ll,ielem) /= -1) then
786 822633546 : edgeptr = edge%putmap(ll,ielem)+1
787 68424243894 : do k=1,vlyr
788 67601610348 : iptr = (kptr+k-1)+edgeptr
789 67601610348 : if (iptr > size(edge%buf)) then
790 0 : write(6, *) 'ERROR SE: ',size(edge%buf),iptr,edge%putmap(ll,ielem)
791 0 : call endrun('pointer bounds ERROR SE')
792 : end if
793 68424243894 : edge%buf(iptr)=v(np, 1, k)
794 : end do
795 : end if
796 : end do
797 :
798 : ! NEAST
799 1647097200 : do ll=swest+3*max_corner_elem,swest+4*max_corner_elem-1
800 1647097200 : if (edge%putmap(ll,ielem) /= -1) then
801 822633546 : edgeptr = edge%putmap(ll,ielem)+1
802 68424243894 : do k=1,vlyr
803 67601610348 : iptr = (kptr+k-1)+edgeptr
804 67601610348 : if (iptr > size(edge%buf)) then
805 0 : write(6, *) 'ERROR NE: ',size(edge%buf),iptr,edge%putmap(ll,ielem)
806 0 : call endrun('pointer bounds ERROR NE')
807 : end if
808 68424243894 : edge%buf(iptr) = v(np, np, k)
809 : end do
810 : end if
811 : end do
812 :
813 : ! NWEST
814 1647097200 : do ll=swest+2*max_corner_elem,swest+3*max_corner_elem-1
815 1647097200 : if (edge%putmap(ll,ielem) /= -1) then
816 822633546 : edgeptr = edge%putmap(ll,ielem)+1
817 68424243894 : do k=1,vlyr
818 67601610348 : iptr = (kptr+k-1)+edgeptr
819 67601610348 : if (iptr > size(edge%buf)) then
820 0 : write(6, *) 'ERROR NW: ',size(edge%buf),iptr,edge%putmap(ll,ielem)
821 0 : call endrun('pointer bounds ERROR NW')
822 : end if
823 68424243894 : edge%buf(iptr) = v(1, np, k)
824 : end do
825 : end if
826 : end do
827 :
828 823548600 : end subroutine edgeVpack
829 :
830 0 : subroutine edgeSpack_r8(edge,v,vlyr,kptr,ielem)
831 : use dimensions_mod, only: np, max_corner_elem
832 : use control_mod, only: north, south, east, west, neast, nwest, seast, swest
833 :
834 : type (EdgeBuffer_t) :: edge
835 : integer, intent(in) :: vlyr
836 : real (kind=r8), intent(in) :: v(vlyr)
837 : integer, intent(in) :: kptr
838 : integer, intent(in) :: ielem
839 :
840 : ! Local variables
841 : integer :: i,k,ir,ll,iptr
842 : integer :: is,ie,in,iw,edgeptr
843 : real (kind=r8) :: tmp
844 :
845 0 : is = edge%putmap(south,ielem)
846 0 : ie = edge%putmap(east,ielem)
847 0 : in = edge%putmap(north,ielem)
848 0 : iw = edge%putmap(west,ielem)
849 0 : if (edge%nlyr < (kptr+vlyr) ) then
850 0 : call endrun('edgeSpack: Buffer overflow: size of the vertical dimension must be increased!')
851 : endif
852 :
853 0 : do k=1,vlyr
854 0 : iptr = kptr+k-1
855 0 : edge%buf(iptr+ie+1) = v(k) ! East
856 0 : edge%buf(iptr+is+1) = v(k) ! South
857 0 : edge%buf(iptr+in+1) = v(k) ! North
858 0 : edge%buf(iptr+iw+1) = v(k) ! West
859 : enddo
860 :
861 : ! SWEST
862 0 : do ll=swest,swest+max_corner_elem-1
863 0 : if (edge%putmap(ll,ielem) /= -1) then
864 0 : edgeptr=edge%putmap(ll,ielem)+1
865 0 : do k=1,vlyr
866 0 : iptr = (kptr+k-1)+edgeptr
867 0 : edge%buf(iptr)=v(k)
868 : end do
869 : end if
870 : end do
871 :
872 : ! SEAST
873 0 : do ll=swest+max_corner_elem,swest+2*max_corner_elem-1
874 0 : if (edge%putmap(ll,ielem) /= -1) then
875 0 : edgeptr=edge%putmap(ll,ielem)+1
876 0 : do k=1,vlyr
877 0 : iptr = (kptr+k-1)+edgeptr
878 0 : edge%buf(iptr)=v(k)
879 : end do
880 : end if
881 : end do
882 :
883 : ! NEAST
884 0 : do ll=swest+3*max_corner_elem,swest+4*max_corner_elem-1
885 0 : if (edge%putmap(ll,ielem) /= -1) then
886 0 : edgeptr=edge%putmap(ll,ielem)+1
887 0 : do k=1,vlyr
888 0 : iptr = (kptr+k-1)+edgeptr
889 0 : edge%buf(iptr)=v(k)
890 : end do
891 : end if
892 : end do
893 :
894 : ! NWEST
895 0 : do ll=swest+2*max_corner_elem,swest+3*max_corner_elem-1
896 0 : if (edge%putmap(ll,ielem) /= -1) then
897 0 : edgeptr=edge%putmap(ll,ielem)+1
898 0 : do k=1,vlyr
899 0 : iptr = (kptr+k-1)+edgeptr
900 0 : edge%buf(iptr)=v(k)
901 : end do
902 : end if
903 : end do
904 :
905 0 : end subroutine edgeSpack_r8
906 :
907 : ! =========================================
908 : ! LongEdgeVpack:
909 : !
910 : ! Pack edges of v into buf...
911 : ! =========================================
912 10800 : subroutine LongEdgeVpack(edge,v,vlyr,kptr,desc)
913 : use control_mod, only: north, south, east, west, neast, nwest, seast, swest
914 : use dimensions_mod, only: np, max_corner_elem
915 :
916 : type (LongEdgeBuffer_t) :: edge
917 : integer, intent(in) :: vlyr
918 : integer , intent(in) :: v(np,np,vlyr)
919 : integer, intent(in) :: kptr
920 : type (EdgeDescriptor_t), intent(in) :: desc
921 :
922 : ! Local variables
923 : logical, parameter :: UseUnroll = .TRUE.
924 : integer :: i,k,ir,l
925 : integer :: is,ie,in,iw
926 :
927 10800 : if(.not. threadsafe) then
928 : !$OMP BARRIER
929 0 : threadsafe=.true.
930 : end if
931 :
932 10800 : is = desc%putmapP(south)
933 10800 : ie = desc%putmapP(east)
934 10800 : in = desc%putmapP(north)
935 10800 : iw = desc%putmapP(west)
936 :
937 : if(MODULO(np,2) == 0 .and. UseUnroll) then
938 21600 : do k=1,vlyr
939 10800 : do i=1,np,2
940 21600 : edge%buf(kptr+k,is+i) = v(i ,1 ,k)
941 21600 : edge%buf(kptr+k,is+i+1) = v(i+1,1 ,k)
942 21600 : edge%buf(kptr+k,ie+i) = v(np ,i ,k)
943 21600 : edge%buf(kptr+k,ie+i+1) = v(np ,i+1 ,k)
944 21600 : edge%buf(kptr+k,in+i) = v(i ,np,k)
945 21600 : edge%buf(kptr+k,in+i+1) = v(i+1 ,np,k)
946 21600 : edge%buf(kptr+k,iw+i) = v(1 ,i ,k)
947 21600 : edge%buf(kptr+k,iw+i+1) = v(1 ,i+1 ,k)
948 :
949 : enddo
950 : end do
951 : else
952 : do k=1,vlyr
953 : do i=1,np
954 : edge%buf(kptr+k,is+i) = v(i ,1 ,k)
955 : edge%buf(kptr+k,ie+i) = v(np ,i ,k)
956 : edge%buf(kptr+k,in+i) = v(i ,np,k)
957 : edge%buf(kptr+k,iw+i) = v(1 ,i ,k)
958 : enddo
959 : end do
960 :
961 : endif
962 :
963 :
964 : ! This is really kludgy way to setup the index reversals
965 : ! But since it is so a rare event not real need to spend time optimizing
966 :
967 10800 : if(desc%reverse(south)) then
968 180 : is = desc%putmapP(south)
969 360 : do k=1,vlyr
970 1080 : do i=1,np
971 720 : ir = np-i+1
972 900 : edge%buf(kptr+k,is+ir)=v(i,1,k)
973 : enddo
974 : enddo
975 : endif
976 :
977 10800 : if(desc%reverse(east)) then
978 60 : ie = desc%putmapP(east)
979 120 : do k=1,vlyr
980 360 : do i=1,np
981 240 : ir = np-i+1
982 300 : edge%buf(kptr+k,ie+ir)=v(np,i,k)
983 : enddo
984 : enddo
985 : endif
986 :
987 10800 : if(desc%reverse(north)) then
988 180 : in = desc%putmapP(north)
989 360 : do k=1,vlyr
990 1080 : do i=1,np
991 720 : ir = np-i+1
992 900 : edge%buf(kptr+k,in+ir)=v(i,np,k)
993 : enddo
994 : enddo
995 : endif
996 :
997 10800 : if(desc%reverse(west)) then
998 60 : iw = desc%putmapP(west)
999 120 : do k=1,vlyr
1000 360 : do i=1,np
1001 240 : ir = np-i+1
1002 300 : edge%buf(kptr+k,iw+ir)=v(1,i,k)
1003 : enddo
1004 : enddo
1005 : endif
1006 :
1007 : ! SWEST
1008 21600 : do l=swest,swest+max_corner_elem-1
1009 21600 : if (desc%putmapP(l) /= -1) then
1010 21576 : do k=1,vlyr
1011 21576 : edge%buf(kptr+k,desc%putmapP(l)+1)=v(1 ,1 ,k)
1012 : end do
1013 : end if
1014 : end do
1015 :
1016 : ! SEAST
1017 21600 : do l=swest+max_corner_elem,swest+2*max_corner_elem-1
1018 21600 : if (desc%putmapP(l) /= -1) then
1019 21576 : do k=1,vlyr
1020 21576 : edge%buf(kptr+k,desc%putmapP(l)+1)=v(np ,1 ,k)
1021 : end do
1022 : end if
1023 : end do
1024 :
1025 : ! NEAST
1026 21600 : do l=swest+3*max_corner_elem,swest+4*max_corner_elem-1
1027 21600 : if (desc%putmapP(l) /= -1) then
1028 21576 : do k=1,vlyr
1029 21576 : edge%buf(kptr+k,desc%putmapP(l)+1)=v(np ,np,k)
1030 : end do
1031 : end if
1032 : end do
1033 :
1034 : ! NWEST
1035 21600 : do l=swest+2*max_corner_elem,swest+3*max_corner_elem-1
1036 21600 : if (desc%putmapP(l) /= -1) then
1037 21576 : do k=1,vlyr
1038 21576 : edge%buf(kptr+k,desc%putmapP(l)+1)=v(1 ,np,k)
1039 : end do
1040 : end if
1041 : end do
1042 :
1043 10800 : end subroutine LongEdgeVpack
1044 :
1045 823548600 : subroutine edgeVunpack(edge,v,vlyr,kptr,ielem,rank)
1046 : use dimensions_mod, only: np, max_corner_elem
1047 : use control_mod, only: north, south, east, west, neast, nwest, seast, swest
1048 :
1049 : type (EdgeBuffer_t), intent(in) :: edge
1050 : integer, intent(in) :: vlyr
1051 : real (kind=r8), intent(inout) :: v(np,np,vlyr)
1052 : integer, intent(in) :: kptr
1053 : integer, intent(in) :: ielem
1054 : integer, optional, intent(in) :: rank
1055 :
1056 : ! Local
1057 : integer :: i,k,ll,iptr
1058 : integer :: is,ie,in,iw,edgeptr
1059 : integer :: ise,isw,ine,inw
1060 : integer :: ks,ke,kblock
1061 : logical :: done
1062 :
1063 823548600 : is=edge%getmap(south,ielem)
1064 823548600 : ie=edge%getmap(east,ielem)
1065 823548600 : in=edge%getmap(north,ielem)
1066 823548600 : iw=edge%getmap(west,ielem)
1067 823548600 : isw=edge%getmap(swest,ielem)
1068 823548600 : ise=edge%getmap(seast,ielem)
1069 823548600 : inw=edge%getmap(nwest,ielem)
1070 823548600 : ine=edge%getmap(neast,ielem)
1071 :
1072 : !DIR$ IVDEP
1073 68500355400 : do k=1,vlyr
1074 67676806800 : iptr=np*(kptr+k-1)
1075 >33920*10^7 : do i=1,np
1076 >27070*10^7 : v(np ,i ,k) = v(np ,i ,k)+edge%receive(iptr+i+ie) ! East
1077 >27070*10^7 : v(i ,1 ,k) = v(i ,1 ,k)+edge%receive(iptr+i+is) ! South
1078 >27070*10^7 : v(i ,np ,k) = v(i ,np ,k)+edge%receive(iptr+i+in) ! North
1079 >33838*10^7 : v(1 ,i ,k) = v(1 ,i ,k)+edge%receive(iptr+i+iw) ! West
1080 : enddo
1081 : enddo
1082 :
1083 : ! SWEST
1084 1647097200 : do ll=swest,swest+max_corner_elem-1
1085 1647097200 : if(edge%getmap(ll,ielem) /= -1) then
1086 822633546 : edgeptr=edge%getmap(ll,ielem)+1
1087 68424243894 : do k=1,vlyr
1088 67601610348 : iptr = (kptr+k-1)+edgeptr
1089 68424243894 : v(1 ,1 ,k)=v(1 ,1 ,k)+edge%receive(iptr)
1090 : enddo
1091 : endif
1092 : end do
1093 :
1094 : ! SEAST
1095 1647097200 : do ll=swest+max_corner_elem,swest+2*max_corner_elem-1
1096 1647097200 : if(edge%getmap(ll,ielem) /= -1) then
1097 822633546 : edgeptr=edge%getmap(ll,ielem)+1
1098 68424243894 : do k=1,vlyr
1099 67601610348 : iptr = (kptr+k-1)+edgeptr
1100 68424243894 : v(np ,1 ,k)=v(np,1 ,k)+edge%receive(iptr)
1101 : enddo
1102 : endif
1103 : end do
1104 :
1105 : ! NEAST
1106 1647097200 : do ll=swest+3*max_corner_elem,swest+4*max_corner_elem-1
1107 1647097200 : if(edge%getmap(ll,ielem) /= -1) then
1108 822633546 : edgeptr=edge%getmap(ll,ielem)+1
1109 68424243894 : do k=1,vlyr
1110 67601610348 : iptr = (kptr+k-1)+edgeptr
1111 68424243894 : v(np ,np,k)=v(np,np,k)+edge%receive(iptr)
1112 : enddo
1113 : endif
1114 : end do
1115 :
1116 : ! NWEST
1117 1647097200 : do ll=swest+2*max_corner_elem,swest+3*max_corner_elem-1
1118 1647097200 : if(edge%getmap(ll,ielem) /= -1) then
1119 822633546 : edgeptr=edge%getmap(ll,ielem)+1
1120 68424243894 : do k=1,vlyr
1121 67601610348 : iptr = (kptr+k-1)+edgeptr
1122 68424243894 : v(1 ,np,k)=v(1 ,np,k)+edge%receive(iptr)
1123 : enddo
1124 : endif
1125 : end do
1126 :
1127 :
1128 823548600 : end subroutine edgeVunpack
1129 : !
1130 0 : subroutine edgeVunpackVert(edge,v,ielem)
1131 : use control_mod, only: north, south, east, west, neast, nwest, seast, swest
1132 : use dimensions_mod, only: np, max_corner_elem, ne
1133 : use coordinate_systems_mod, only: cartesian3D_t
1134 :
1135 : type (EdgeBuffer_t), intent(inout) :: edge
1136 : type (cartesian3D_t), intent(inout) :: v(:,:,:)
1137 : integer, intent(in) :: ielem
1138 :
1139 : ! Local
1140 : logical, parameter :: UseUnroll = .TRUE.
1141 : integer :: i,k,l, nce
1142 : integer :: is,ie,in,iw,ine,inw,isw,ise
1143 :
1144 0 : threadsafe=.false.
1145 :
1146 0 : if (max_corner_elem.ne.1 .and. ne==0) then
1147 : ! MNL: this is used to construct the dual grid on the cube,
1148 : ! currently only supported for the uniform grid. If
1149 : ! this is desired on a refined grid, a little bit of
1150 : ! work will be required.
1151 0 : call endrun("edgeVunpackVert should not be called with unstructured meshes")
1152 : end if
1153 :
1154 0 : is=edge%getmap(south,ielem)
1155 0 : ie=edge%getmap(east,ielem)
1156 0 : in=edge%getmap(north,ielem)
1157 0 : iw=edge%getmap(west,ielem)
1158 :
1159 :
1160 : ! N+S
1161 0 : do i=1,np/2
1162 : ! North
1163 0 : v(3,i ,np)%x = edge%receive(in+i)
1164 0 : v(3,i ,np)%y = edge%receive(np+in+i)
1165 0 : v(3,i ,np)%z = edge%receive(2*np+in+i)
1166 :
1167 : ! South
1168 0 : v(2,i ,1)%x = edge%receive(is+i)
1169 0 : v(2,i ,1)%y = edge%receive(np+is+i)
1170 0 : v(2,i ,1)%z = edge%receive(2*np+is+i)
1171 : enddo
1172 :
1173 0 : do i=np/2+1,np
1174 : ! North
1175 0 : v(4,i ,np)%x = edge%receive(in+i)
1176 0 : v(4,i ,np)%y = edge%receive(np+in+i)
1177 0 : v(4,i ,np)%z = edge%receive(2*np+in+i)
1178 : ! South
1179 0 : v(1,i ,1)%x = edge%receive(is+i)
1180 0 : v(1,i ,1)%y = edge%receive(np+is+i)
1181 0 : v(1,i ,1)%z = edge%receive(2*np+is+i)
1182 : enddo
1183 :
1184 0 : do i=1,np/2
1185 : ! East
1186 0 : v(3,np,i)%x = edge%receive(ie+i)
1187 0 : v(3,np,i)%y = edge%receive(np+ie+i)
1188 0 : v(3,np,i)%z = edge%receive(2*np+ie+i)
1189 : ! West
1190 0 : v(4,1,i)%x = edge%receive(iw+i)
1191 0 : v(4,1,i)%y = edge%receive(np+iw+i)
1192 0 : v(4,1,i)%z = edge%receive(2*np+iw+i)
1193 : end do
1194 :
1195 0 : do i=np/2+1,np
1196 : ! East
1197 0 : v(2,np,i)%x = edge%receive(ie+i)
1198 0 : v(2,np,i)%y = edge%receive(np+ie+i)
1199 0 : v(2,np,i)%z = edge%receive(2*np+ie+i)
1200 : ! West
1201 0 : v(1,1,i)%x = edge%receive(iw+i)
1202 0 : v(1,1,i)%y = edge%receive(np+iw+i)
1203 0 : v(1,1,i)%z = edge%receive(2*np+iw+i)
1204 : end do
1205 :
1206 : ! SWEST
1207 0 : nce = max_corner_elem
1208 0 : do l=swest,swest+max_corner_elem-1
1209 : ! find the one active corner, then exist
1210 0 : isw=edge%getmap(l,ielem)
1211 0 : if(isw /= -1) then
1212 0 : v(1,1,1)%x=edge%receive(isw+1)
1213 0 : v(1,1,1)%y=edge%receive(nce+isw+1)
1214 0 : v(1,1,1)%z=edge%receive(2*nce+isw+1)
1215 0 : exit
1216 : else
1217 0 : v(1,1,1)%x=0_r8
1218 0 : v(1,1,1)%y=0_r8
1219 0 : v(1,1,1)%z=0_r8
1220 : endif
1221 : end do
1222 :
1223 : ! SEAST
1224 0 : do l=swest+max_corner_elem,swest+2*max_corner_elem-1
1225 : ! find the one active corner, then exist
1226 0 : ise=edge%getmap(l,ielem)
1227 0 : if(ise /= -1) then
1228 0 : v(2,np,1)%x=edge%receive(ise+1)
1229 0 : v(2,np,1)%y=edge%receive(nce+ise+1)
1230 0 : v(2,np,1)%z=edge%receive(2*nce+ise+1)
1231 0 : exit
1232 : else
1233 0 : v(2,np,1)%x=0_r8
1234 0 : v(2,np,1)%y=0_r8
1235 0 : v(2,np,1)%z=0_r8
1236 : endif
1237 : end do
1238 :
1239 : ! NEAST
1240 0 : do l=swest+3*max_corner_elem,swest+4*max_corner_elem-1
1241 : ! find the one active corner, then exist
1242 0 : ine=edge%getmap(l,ielem)
1243 0 : if(ine /= -1) then
1244 0 : v(3,np,np)%x=edge%receive(ine+1)
1245 0 : v(3,np,np)%y=edge%receive(nce+ine+1)
1246 0 : v(3,np,np)%z=edge%receive(2*nce+ine+1)
1247 0 : exit
1248 : else
1249 0 : v(3,np,np)%x=0_r8
1250 0 : v(3,np,np)%y=0_r8
1251 0 : v(3,np,np)%z=0_r8
1252 : endif
1253 : end do
1254 :
1255 : ! NWEST
1256 0 : do l=swest+2*max_corner_elem,swest+3*max_corner_elem-1
1257 : ! find the one active corner, then exist
1258 0 : inw = edge%getmap(l,ielem)
1259 0 : if(inw/= -1) then
1260 0 : v(4,1,np)%x=edge%receive(inw+1)
1261 0 : v(4,1,np)%y=edge%receive(nce+inw+1)
1262 0 : v(4,1,np)%z=edge%receive(2*nce+inw+1)
1263 0 : exit
1264 : else
1265 0 : v(4,1,np)%x=0_r8
1266 0 : v(4,1,np)%y=0_r8
1267 0 : v(4,1,np)%z=0_r8
1268 : endif
1269 : end do
1270 :
1271 : ! Fill the missing vertex info
1272 :
1273 0 : do i=2,np/2
1274 : ! North
1275 0 : v(4,i ,np)%x = v(3,i-1 ,np)%x
1276 0 : v(4,i ,np)%y = v(3,i-1 ,np)%y
1277 0 : v(4,i ,np)%z = v(3,i-1 ,np)%z
1278 : ! South
1279 0 : v(1,i ,1)%x = v(2,i-1 ,1)%x
1280 0 : v(1,i ,1)%y = v(2,i-1 ,1)%y
1281 0 : v(1,i ,1)%z = v(2,i-1 ,1)%z
1282 : enddo
1283 :
1284 0 : do i=np/2+1,np-1
1285 : ! North
1286 0 : v(3,i ,np)%x = v(4,i+1 ,np)%x
1287 0 : v(3,i ,np)%y = v(4,i+1 ,np)%y
1288 0 : v(3,i ,np)%z = v(4,i+1 ,np)%z
1289 : ! South
1290 0 : v(2,i ,1)%x = v(1,i+1 ,1)%x
1291 0 : v(2,i ,1)%y = v(1,i+1 ,1)%y
1292 0 : v(2,i ,1)%z = v(1,i+1 ,1)%z
1293 : enddo
1294 :
1295 0 : do i=2,np/2
1296 : ! East
1297 0 : v(2,np,i)%x = v(3,np,i-1)%x
1298 0 : v(2,np,i)%y = v(3,np,i-1)%y
1299 0 : v(2,np,i)%z = v(3,np,i-1)%z
1300 : ! West
1301 0 : v(1,1,i)%x = v(4,1,i-1)%x
1302 0 : v(1,1,i)%y = v(4,1,i-1)%y
1303 0 : v(1,1,i)%z = v(4,1,i-1)%z
1304 : end do
1305 :
1306 0 : do i=np/2+1,np-1
1307 : ! East
1308 0 : v(3,np,i)%x = v(2,np,i+1)%x
1309 0 : v(3,np,i)%y = v(2,np,i+1)%y
1310 0 : v(3,np,i)%z = v(2,np,i+1)%z
1311 : ! West
1312 0 : v(4,1,i)%x = v(1,1,i+1)%x
1313 0 : v(4,1,i)%y = v(1,1,i+1)%y
1314 0 : v(4,1,i)%z = v(1,1,i+1)%z
1315 : end do
1316 :
1317 0 : end subroutine edgeVunpackVert
1318 : ! ========================================
1319 : ! edgeDGVunpack:
1320 : !
1321 : ! Unpack edges from edge buffer into v...
1322 : ! ========================================
1323 :
1324 124675200 : subroutine edgeDGVunpack(edge,v,vlyr,kptr,ielem)
1325 : use dimensions_mod, only: np, max_corner_elem
1326 : use control_mod, only: north, south, east, west, neast, nwest, seast, swest
1327 :
1328 : type (EdgeBuffer_t), intent(in) :: edge
1329 : integer, intent(in) :: vlyr
1330 : real (kind=r8), intent(inout) :: v(0:np+1,0:np+1,vlyr)
1331 : integer, intent(in) :: kptr
1332 : integer, intent(in) :: ielem
1333 :
1334 : ! Local
1335 : integer :: i,k,iptr
1336 : integer :: is,ie,in,iw
1337 :
1338 124675200 : threadsafe=.false.
1339 :
1340 124675200 : is=edge%getmap(south,ielem)
1341 124675200 : ie=edge%getmap(east,ielem)
1342 124675200 : in=edge%getmap(north,ielem)
1343 124675200 : iw=edge%getmap(west,ielem)
1344 7558434000 : do k=1,vlyr
1345 7433758800 : iptr=np*(kptr+k-1)
1346 37293469200 : do i=1,np
1347 29735035200 : v(i ,0 ,k)=edge%receive(iptr+is+i)
1348 29735035200 : v(np+1,i ,k)=edge%receive(iptr+ie+i)
1349 29735035200 : v(i ,np+1,k)=edge%receive(iptr+in+i)
1350 37168794000 : v(0 ,i ,k)=edge%receive(iptr+iw+i)
1351 : end do
1352 : end do
1353 :
1354 124675200 : i = swest
1355 124675200 : if(edge%getmap(i,ielem) /= -1) then
1356 7550035740 : do k=1,vlyr
1357 7425499068 : iptr=(kptr+k-1)
1358 7550035740 : v(0,0,k) = edge%receive(iptr+edge%getmap(i,ielem)+1)
1359 : end do
1360 : end if
1361 124675200 : i = swest+max_corner_elem
1362 124675200 : if(edge%getmap(i,ielem) /= -1) then
1363 7550035740 : do k=1,vlyr
1364 7425499068 : iptr=(kptr+k-1)
1365 7550035740 : v(np+1,0,k) = edge%receive(iptr+edge%getmap(i,ielem)+1)
1366 : end do
1367 : end if
1368 124675200 : i = swest+3*max_corner_elem
1369 124675200 : if(edge%getmap(i,ielem) /= -1) then
1370 7550035740 : do k=1,vlyr
1371 7425499068 : iptr=(kptr+k-1)
1372 7550035740 : v(np+1,np+1,k) = edge%receive(iptr+edge%getmap(i,ielem)+1)
1373 : end do
1374 : end if
1375 124675200 : i = swest+2*max_corner_elem
1376 124675200 : if(edge%getmap(i,ielem) /= -1) then
1377 7550035740 : do k=1,vlyr
1378 7425499068 : iptr=(kptr+k-1)
1379 7550035740 : v(0,np+1,k) = edge%receive(iptr+edge%getmap(i,ielem)+1)
1380 : end do
1381 : end if
1382 :
1383 124675200 : end subroutine edgeDGVunpack
1384 :
1385 : ! ========================================
1386 : ! edgeVunpackMIN/MAX:
1387 : !
1388 : ! Finds the Min/Max edges from edge buffer into v...
1389 : ! ========================================
1390 0 : subroutine edgeVunpackMAX(edge,v,vlyr,kptr,ielem)
1391 : use dimensions_mod, only: np, max_corner_elem
1392 : use control_mod, only: north, south, east, west, neast, nwest, seast, swest
1393 :
1394 : type (EdgeBuffer_t), intent(in) :: edge
1395 : integer, intent(in) :: vlyr
1396 : real (kind=r8), intent(inout) :: v(np,np,vlyr)
1397 : integer, intent(in) :: kptr
1398 : integer, intent(in) :: ielem
1399 :
1400 : ! Local
1401 : integer :: i,k,l,iptr
1402 : integer :: is,ie,in,iw
1403 :
1404 0 : threadsafe=.false.
1405 :
1406 0 : is=edge%getmap(south,ielem)
1407 0 : ie=edge%getmap(east,ielem)
1408 0 : in=edge%getmap(north,ielem)
1409 0 : iw=edge%getmap(west,ielem)
1410 0 : do k=1,vlyr
1411 0 : iptr=np*(kptr+k-1)
1412 0 : do i=1,np
1413 0 : v(np ,i ,k) = MAX(v(np ,i ,k),edge%receive(iptr+ie+i ))
1414 0 : v(i ,1 ,k) = MAX(v(i ,1 ,k),edge%receive(iptr+is+i ))
1415 0 : v(i ,np ,k) = MAX(v(i ,np ,k),edge%receive(iptr+in+i ))
1416 0 : v(1 ,i ,k) = MAX(v(1 ,i ,k),edge%receive(iptr+iw+i ))
1417 : end do
1418 : end do
1419 :
1420 : ! SWEST
1421 0 : do l=swest,swest+max_corner_elem-1
1422 0 : if(edge%getmap(l,ielem) /= -1) then
1423 0 : do k=1,vlyr
1424 0 : v(1 ,1 ,k)=MAX(v(1 ,1 ,k),edge%receive((kptr+k-1)+edge%getmap(l,ielem)+1))
1425 : enddo
1426 : endif
1427 : end do
1428 :
1429 : ! SEAST
1430 0 : do l=swest+max_corner_elem,swest+2*max_corner_elem-1
1431 0 : if(edge%getmap(l,ielem) /= -1) then
1432 0 : do k=1,vlyr
1433 0 : v(np ,1 ,k)=MAX(v(np,1 ,k),edge%receive((kptr+k-1)+edge%getmap(l,ielem)+1))
1434 : enddo
1435 : endif
1436 : end do
1437 :
1438 : ! NEAST
1439 0 : do l=swest+3*max_corner_elem,swest+4*max_corner_elem-1
1440 0 : if(edge%getmap(l,ielem) /= -1) then
1441 0 : do k=1,vlyr
1442 0 : v(np ,np,k)=MAX(v(np,np,k),edge%receive((kptr+k-1)+edge%getmap(l,ielem)+1))
1443 : enddo
1444 : endif
1445 : end do
1446 :
1447 : ! NWEST
1448 0 : do l=swest+2*max_corner_elem,swest+3*max_corner_elem-1
1449 0 : if(edge%getmap(l,ielem) /= -1) then
1450 0 : do k=1,vlyr
1451 0 : v(1 ,np,k)=MAX(v(1 ,np,k),edge%receive((kptr+k-1)+edge%getmap(l,ielem)+1))
1452 : enddo
1453 : endif
1454 : end do
1455 :
1456 0 : end subroutine edgeVunpackMAX
1457 :
1458 0 : subroutine edgeSunpackMAX(edge,v,vlyr,kptr,ielem)
1459 : use dimensions_mod, only: np, max_corner_elem
1460 : use control_mod, only: north, south, east, west, neast, nwest, seast, swest
1461 :
1462 : type (EdgeBuffer_t), intent(in) :: edge
1463 : integer, intent(in) :: vlyr
1464 : real (kind=r8), intent(inout) :: v(vlyr)
1465 : integer, intent(in) :: kptr
1466 : integer, intent(in) :: ielem
1467 :
1468 : ! Local
1469 : integer :: i,k,l,iptr
1470 : integer :: is,ie,in,iw,edgeptr
1471 :
1472 0 : threadsafe=.false.
1473 :
1474 0 : is=edge%getmap(south,ielem)
1475 0 : ie=edge%getmap(east,ielem)
1476 0 : in=edge%getmap(north,ielem)
1477 0 : iw=edge%getmap(west,ielem)
1478 0 : do k=1,vlyr
1479 0 : iptr=(kptr+k-1)
1480 0 : v(k) = MAX(v(k),edge%receive(iptr+is+1),edge%receive(iptr+ie+1),edge%receive(iptr+in+1),edge%receive(iptr+iw+1))
1481 : end do
1482 :
1483 : ! SWEST
1484 0 : do l=swest,swest+max_corner_elem-1
1485 0 : if(edge%getmap(l,ielem) /= -1) then
1486 0 : edgeptr = edge%getmap(l,ielem)+1
1487 0 : do k=1,vlyr
1488 0 : iptr = (kptr+k-1)+edgeptr
1489 0 : v(k)=MAX(v(k),edge%receive(iptr))
1490 : enddo
1491 : endif
1492 : end do
1493 :
1494 : ! SEAST
1495 0 : do l=swest+max_corner_elem,swest+2*max_corner_elem-1
1496 0 : if(edge%getmap(l,ielem) /= -1) then
1497 0 : edgeptr = edge%getmap(l,ielem)+1
1498 0 : do k=1,vlyr
1499 0 : iptr = (kptr+k-1)+edgeptr
1500 0 : v(k)=MAX(v(k),edge%receive(iptr))
1501 : enddo
1502 : endif
1503 : end do
1504 :
1505 : ! NEAST
1506 0 : do l=swest+3*max_corner_elem,swest+4*max_corner_elem-1
1507 0 : if(edge%getmap(l,ielem) /= -1) then
1508 0 : edgeptr = edge%getmap(l,ielem)+1
1509 0 : do k=1,vlyr
1510 0 : iptr = (kptr+k-1)+edgeptr
1511 0 : v(k)=MAX(v(k),edge%receive(iptr))
1512 : enddo
1513 : endif
1514 : end do
1515 :
1516 : ! NWEST
1517 0 : do l=swest+2*max_corner_elem,swest+3*max_corner_elem-1
1518 0 : if(edge%getmap(l,ielem) /= -1) then
1519 0 : edgeptr = edge%getmap(l,ielem)+1
1520 0 : do k=1,vlyr
1521 0 : iptr = (kptr+k-1)+edgeptr
1522 0 : v(k)=MAX(v(k),edge%receive(iptr))
1523 : enddo
1524 : endif
1525 : end do
1526 :
1527 0 : end subroutine edgeSunpackMAX
1528 :
1529 0 : subroutine edgeSunpackMIN(edge,v,vlyr,kptr,ielem)
1530 : use dimensions_mod, only: np, max_corner_elem
1531 : use control_mod, only: north, south, east, west, neast, nwest, seast, swest
1532 :
1533 : type (EdgeBuffer_t), intent(in) :: edge
1534 : integer, intent(in) :: vlyr
1535 : real (kind=r8), intent(inout) :: v(vlyr)
1536 : integer, intent(in) :: kptr
1537 : integer, intent(in) :: ielem
1538 :
1539 : ! Local
1540 : integer :: i,k,l,iptr
1541 : integer :: is,ie,in,iw,edgeptr
1542 :
1543 0 : threadsafe=.false.
1544 :
1545 0 : is=edge%getmap(south,ielem)
1546 0 : ie=edge%getmap(east,ielem)
1547 0 : in=edge%getmap(north,ielem)
1548 0 : iw=edge%getmap(west,ielem)
1549 0 : do k=1,vlyr
1550 0 : iptr=(kptr+k-1)
1551 0 : v(k) = MIN(v(k),edge%receive(iptr+is+1),edge%receive(iptr+ie+1),edge%receive(iptr+in+1),edge%receive(iptr+iw+1))
1552 : end do
1553 :
1554 : ! SWEST
1555 0 : do l=swest,swest+max_corner_elem-1
1556 0 : if(edge%getmap(l,ielem) /= -1) then
1557 0 : edgeptr = edge%getmap(l,ielem)+1
1558 0 : do k=1,vlyr
1559 0 : iptr = (kptr+k-1)+edgeptr
1560 0 : v(k)=MiN(v(k),edge%receive(iptr))
1561 : enddo
1562 : endif
1563 : end do
1564 :
1565 : ! SEAST
1566 0 : do l=swest+max_corner_elem,swest+2*max_corner_elem-1
1567 0 : if(edge%getmap(l,ielem) /= -1) then
1568 0 : edgeptr = edge%getmap(l,ielem)+1
1569 0 : do k=1,vlyr
1570 0 : iptr = (kptr+k-1)+edgeptr
1571 0 : v(k)=MIN(v(k),edge%receive(iptr))
1572 : enddo
1573 : endif
1574 : end do
1575 :
1576 : ! NEAST
1577 0 : do l=swest+3*max_corner_elem,swest+4*max_corner_elem-1
1578 0 : if(edge%getmap(l,ielem) /= -1) then
1579 0 : edgeptr = edge%getmap(l,ielem)+1
1580 0 : do k=1,vlyr
1581 0 : iptr = (kptr+k-1)+edgeptr
1582 0 : v(k)=MIN(v(k),edge%receive(iptr))
1583 : enddo
1584 : endif
1585 : end do
1586 :
1587 : ! NWEST
1588 0 : do l=swest+2*max_corner_elem,swest+3*max_corner_elem-1
1589 0 : if(edge%getmap(l,ielem) /= -1) then
1590 0 : edgeptr = edge%getmap(l,ielem)+1
1591 0 : do k=1,vlyr
1592 0 : iptr = (kptr+k-1)+edgeptr
1593 0 : v(k)=MIN(v(k),edge%receive(iptr))
1594 : enddo
1595 : endif
1596 : end do
1597 :
1598 0 : end subroutine edgeSunpackMIN
1599 :
1600 0 : subroutine edgeVunpackMIN(edge,v,vlyr,kptr,ielem)
1601 : use dimensions_mod, only: np, max_corner_elem
1602 : use control_mod, only: north, south, east, west, neast, nwest, seast, swest
1603 :
1604 : type (EdgeBuffer_t), intent(in) :: edge
1605 : integer, intent(in) :: vlyr
1606 : real (kind=r8), intent(inout) :: v(np,np,vlyr)
1607 : integer, intent(in) :: kptr
1608 : integer, intent(in) :: ielem
1609 :
1610 : ! Local
1611 : integer :: i,k,l,iptr
1612 : integer :: is,ie,in,iw,edgeptr
1613 :
1614 0 : threadsafe=.false.
1615 :
1616 0 : is=edge%getmap(south,ielem)
1617 0 : ie=edge%getmap(east,ielem)
1618 0 : in=edge%getmap(north,ielem)
1619 0 : iw=edge%getmap(west,ielem)
1620 0 : do k=1,vlyr
1621 0 : iptr = np*(kptr+k-1)
1622 0 : do i=1,np
1623 0 : v(np ,i ,k) = MIN(v(np ,i ,k),edge%receive(iptr+ie+i ))
1624 0 : v(i ,1 ,k) = MIN(v(i ,1 ,k),edge%receive(iptr+is+i ))
1625 0 : v(i ,np ,k) = MIN(v(i ,np ,k),edge%receive(iptr+in+i ))
1626 0 : v(1 ,i ,k) = MIN(v(1 ,i ,k),edge%receive(iptr+iw+i ))
1627 : end do
1628 : end do
1629 :
1630 : ! SWEST
1631 0 : do l=swest,swest+max_corner_elem-1
1632 0 : if(edge%getmap(l,ielem) /= -1) then
1633 0 : edgeptr=edge%getmap(l,ielem)+1
1634 0 : do k=1,vlyr
1635 0 : iptr=(kptr+k-1)+edgeptr
1636 0 : v(1 ,1 ,k)=MIN(v(1 ,1 ,k),edge%receive(iptr))
1637 : enddo
1638 : endif
1639 : end do
1640 :
1641 : ! SEAST
1642 0 : do l=swest+max_corner_elem,swest+2*max_corner_elem-1
1643 0 : if(edge%getmap(l,ielem) /= -1) then
1644 0 : edgeptr=edge%getmap(l,ielem)+1
1645 0 : do k=1,vlyr
1646 0 : iptr=(kptr+k-1)+edgeptr
1647 0 : v(np ,1 ,k)=MIN(v(np,1 ,k),edge%receive(iptr))
1648 : enddo
1649 : endif
1650 : end do
1651 :
1652 : ! NEAST
1653 0 : do l=swest+3*max_corner_elem,swest+4*max_corner_elem-1
1654 0 : if(edge%getmap(l,ielem) /= -1) then
1655 0 : edgeptr=edge%getmap(l,ielem)+1
1656 0 : do k=1,vlyr
1657 0 : iptr=(kptr+k-1)+edgeptr
1658 0 : v(np ,np,k)=MIN(v(np,np,k),edge%receive(iptr))
1659 : enddo
1660 : endif
1661 : end do
1662 :
1663 : ! NWEST
1664 0 : do l=swest+2*max_corner_elem,swest+3*max_corner_elem-1
1665 0 : if(edge%getmap(l,ielem) /= -1) then
1666 0 : edgeptr=edge%getmap(l,ielem)+1
1667 0 : do k=1,vlyr
1668 0 : iptr=(kptr+k-1)+edgeptr
1669 0 : v(1 ,np,k)=MIN(v(1 ,np,k),edge%receive(iptr))
1670 : enddo
1671 : endif
1672 : end do
1673 :
1674 0 : end subroutine edgeVunpackMIN
1675 :
1676 : ! ========================================
1677 : ! LongEdgeVunpackMIN:
1678 : !
1679 : ! Finds the Min edges from edge buffer into v...
1680 : ! ========================================
1681 10800 : subroutine LongEdgeVunpackMIN(edge,v,vlyr,kptr,desc)
1682 : use control_mod, only: north, south, east, west, neast, nwest, seast, swest
1683 : use dimensions_mod, only: np, max_corner_elem
1684 :
1685 : type (LongEdgeBuffer_t), intent(in) :: edge
1686 : integer, intent(in) :: vlyr
1687 : integer , intent(inout) :: v(np,np,vlyr)
1688 : integer, intent(in) :: kptr
1689 : type (EdgeDescriptor_t), intent(in) :: desc
1690 :
1691 : ! Local
1692 :
1693 : integer :: i,k,l
1694 : integer :: is,ie,in,iw
1695 :
1696 10800 : threadsafe=.false.
1697 :
1698 10800 : is=desc%getmapP(south)
1699 10800 : ie=desc%getmapP(east)
1700 10800 : in=desc%getmapP(north)
1701 10800 : iw=desc%getmapP(west)
1702 21600 : do k=1,vlyr
1703 64800 : do i=1,np
1704 43200 : v(i ,1 ,k) = MIN(v(i ,1 ,k),edge%buf(kptr+k,is+i ))
1705 43200 : v(np ,i ,k) = MIN(v(np ,i ,k),edge%buf(kptr+k,ie+i ))
1706 43200 : v(i ,np ,k) = MIN(v(i ,np ,k),edge%buf(kptr+k,in+i ))
1707 54000 : v(1 ,i ,k) = MIN(v(1 ,i ,k),edge%buf(kptr+k,iw+i ))
1708 : end do
1709 : end do
1710 :
1711 : ! SWEST
1712 21600 : do l=swest,swest+max_corner_elem-1
1713 21600 : if(desc%getmapP(l) /= -1) then
1714 21576 : do k=1,vlyr
1715 21576 : v(1 ,1 ,k)=MIN(v(1 ,1 ,k),edge%buf(kptr+k,desc%getmapP(l)+1))
1716 : enddo
1717 : endif
1718 : end do
1719 :
1720 : ! SEAST
1721 21600 : do l=swest+max_corner_elem,swest+2*max_corner_elem-1
1722 21600 : if(desc%getmapP(l) /= -1) then
1723 21576 : do k=1,vlyr
1724 21576 : v(np ,1 ,k)=MIN(v(np,1 ,k),edge%buf(kptr+k,desc%getmapP(l)+1))
1725 : enddo
1726 : endif
1727 : end do
1728 :
1729 : ! NEAST
1730 21600 : do l=swest+3*max_corner_elem,swest+4*max_corner_elem-1
1731 21600 : if(desc%getmapP(l) /= -1) then
1732 21576 : do k=1,vlyr
1733 21576 : v(np ,np,k)=MIN(v(np,np,k),edge%buf(kptr+k,desc%getmapP(l)+1))
1734 : enddo
1735 : endif
1736 : end do
1737 :
1738 : ! NWEST
1739 21600 : do l=swest+2*max_corner_elem,swest+3*max_corner_elem-1
1740 21600 : if(desc%getmapP(l) /= -1) then
1741 21576 : do k=1,vlyr
1742 21576 : v(1 ,np,k)=MIN(v(1 ,np,k),edge%buf(kptr+k,desc%getmapP(l)+1))
1743 : enddo
1744 : endif
1745 : end do
1746 :
1747 10800 : end subroutine LongEdgeVunpackMIN
1748 :
1749 :
1750 138018600 : subroutine ghostpack(edge,v,vlyr,kptr,ielem)
1751 :
1752 : use dimensions_mod, only : max_corner_elem
1753 : use control_mod, only : north, south, east, west, neast, nwest, seast, swest
1754 : use edgetype_mod, only : EdgeDescriptor_t
1755 :
1756 : implicit none
1757 :
1758 : type (Edgebuffer_t) :: edge
1759 : integer, intent(in) :: vlyr
1760 : integer, intent(in) :: kptr
1761 :
1762 : real (kind=r8),intent(in) :: v(edge%lb:edge%ub,edge%lb:edge%ub,vlyr)
1763 : integer, intent(in) :: ielem
1764 :
1765 : ! Local variables
1766 : integer :: i,j,k,ir,l,itr,ktmp
1767 :
1768 : integer :: is,ie,in,iw,isw,ise,inw,ine
1769 : integer :: nhc, npoints
1770 : integer :: edgeptr,iptr
1771 :
1772 138018600 : is = edge%putmap(south,ielem)
1773 138018600 : ie = edge%putmap(east,ielem)
1774 138018600 : in = edge%putmap(north,ielem)
1775 138018600 : iw = edge%putmap(west,ielem)
1776 138018600 : if (edge%nlyr < (kptr+vlyr) ) then
1777 0 : print *,'edge%nlyr = ',edge%nlyr
1778 0 : print *,'kptr+vlyr = ',kptr+vlyr
1779 0 : call endrun('ghostpack: Buffer overflow: size of the vertical dimension must be increased!')
1780 : endif
1781 :
1782 :
1783 138018600 : nhc = edge%ndepth
1784 138018600 : npoints = edge%npoints
1785 :
1786 : !DIR$ IVDEP
1787 17289007200 : do k=1,vlyr
1788 17150988600 : ktmp = nhc*(kptr+k-1)
1789 53282237400 : do j=1,nhc
1790 35993230200 : iptr = npoints*(ktmp + j - 1)
1791 >16112*10^7 : do i=1,npoints
1792 >10797*10^7 : edge%buf(iptr+is+i) = v(i ,j ,k)
1793 >10797*10^7 : edge%buf(iptr+ie+i) = v(npoints-j+1 ,i ,k)
1794 >10797*10^7 : edge%buf(iptr+in+i) = v(i ,npoints-j+1,k)
1795 >14397*10^7 : edge%buf(iptr+iw+i) = v(j ,i ,k)
1796 : enddo
1797 : end do
1798 : end do
1799 :
1800 :
1801 : ! This is really kludgy way to setup the index reversals
1802 : ! But since it is so a rare event not real need to spend time optimizing
1803 : ! Check if the edge orientation of the recieving element is different
1804 : ! if it is, swap the order of data in the edge
1805 138018600 : if(edge%reverse(south,ielem)) then
1806 : !DIR$ IVDEP
1807 288150120 : do k=1,vlyr
1808 285849810 : ktmp = nhc*(kptr+k-1)
1809 888037290 : do j=1,nhc
1810 599887170 : iptr = npoints*(ktmp + j - 1)
1811 2685398130 : do i=1,npoints
1812 1799661150 : ir = npoints-i+1
1813 2399548320 : edge%buf(iptr+is+i)=v(ir,j,k)
1814 : enddo
1815 : enddo
1816 : enddo
1817 : endif
1818 :
1819 138018600 : if(edge%reverse(east,ielem)) then
1820 : !DIR$ IVDEP
1821 96050040 : do k=1,vlyr
1822 95283270 : ktmp = nhc*(kptr+k-1)
1823 296012430 : do j=1,nhc
1824 199962390 : iptr = npoints*(ktmp + j - 1)
1825 895132710 : do i=1,npoints
1826 599887050 : ir = npoints-i+1
1827 799849440 : edge%buf(iptr+ie+i)=v(npoints-j+1,ir,k)
1828 : enddo
1829 : enddo
1830 : enddo
1831 : endif
1832 :
1833 138018600 : if(edge%reverse(north,ielem)) then
1834 : !DIR$ IVDEP
1835 288150120 : do k=1,vlyr
1836 285849810 : ktmp = nhc*(kptr+k-1)
1837 888037290 : do j=1,nhc
1838 599887170 : iptr = npoints*(ktmp + j - 1)
1839 2685398130 : do i=1,npoints
1840 1799661150 : ir = npoints-i+1
1841 2399548320 : edge%buf(iptr+in+i)=v(ir,npoints-j+1,k)
1842 : enddo
1843 : enddo
1844 : enddo
1845 : endif
1846 :
1847 138018600 : if(edge%reverse(west,ielem)) then
1848 : !DIR$ IVDEP
1849 96050040 : do k=1,vlyr
1850 95283270 : ktmp = nhc*(kptr+k-1)
1851 296012430 : do j=1,nhc
1852 199962390 : iptr = npoints*(ktmp + j - 1)
1853 895132710 : do i=1,npoints
1854 599887050 : ir = npoints-i+1
1855 799849440 : edge%buf(iptr+iw+i)=v(j,ir,k)
1856 : enddo
1857 : enddo
1858 : enddo
1859 : endif
1860 :
1861 :
1862 : ! corners. this is difficult because we dont know the orientaton
1863 : ! of the corners, and this which (i,j) dimension maps to which dimension
1864 : ! SWEST
1865 276037200 : do l=swest,swest+max_corner_elem-1
1866 276037200 : if (edge%putmap(l,ielem) /= -1) then
1867 : isw = edge%putmap(l,ielem)
1868 : !DIR$ IVDEP
1869 17269797192 : do k=1,vlyr
1870 17131931946 : ktmp = nhc*(kptr+k-1)
1871 53223034914 : do j=1,nhc
1872 35953237722 : iptr = nhc*(ktmp + j - 1)
1873 >14550*10^7 : do i=1,nhc
1874 >12837*10^7 : edge%buf(iptr+isw+i)=v(i ,j ,k)
1875 : enddo
1876 : end do
1877 : end do
1878 : end if
1879 : end do
1880 :
1881 : ! SEAST
1882 276037200 : do l=swest+max_corner_elem,swest+2*max_corner_elem-1
1883 276037200 : if (edge%putmap(l,ielem) /= -1) then
1884 : ise = edge%putmap(l,ielem)
1885 : !DIR$ IVDEP
1886 17269797192 : do k=1,vlyr
1887 17131931946 : ktmp = nhc*(kptr+k-1)
1888 53223034914 : do j=1,nhc
1889 35953237722 : iptr = nhc*(ktmp + j - 1)
1890 >14550*10^7 : do i=1,nhc
1891 >12837*10^7 : edge%buf(iptr+ise+i)=v(npoints-i+1 ,j ,k)
1892 : enddo
1893 : end do
1894 : end do
1895 : end if
1896 : end do
1897 :
1898 : ! NEAST
1899 276037200 : do l=swest+3*max_corner_elem,swest+4*max_corner_elem-1
1900 276037200 : if (edge%putmap(l,ielem) /= -1) then
1901 : ine = edge%putmap(l,ielem)
1902 : !DIR$ IVDEP
1903 17269797192 : do k=1,vlyr
1904 17131931946 : ktmp = nhc*(kptr+k-1)
1905 53223034914 : do j=1,nhc
1906 35953237722 : iptr = nhc*(ktmp + j - 1)
1907 >14550*10^7 : do i=1,nhc
1908 >12837*10^7 : edge%buf(iptr+ine+i)=v(npoints-i+1,npoints-j+1,k)
1909 : enddo
1910 : enddo
1911 : end do
1912 : end if
1913 : end do
1914 :
1915 : ! NWEST
1916 276037200 : do l=swest+2*max_corner_elem,swest+3*max_corner_elem-1
1917 276037200 : if (edge%putmap(l,ielem) /= -1) then
1918 : inw = edge%putmap(l,ielem)
1919 : !DIR$ IVDEP
1920 17269797192 : do k=1,vlyr
1921 17131931946 : ktmp = nhc*(kptr+k-1)
1922 53223034914 : do j=1,nhc
1923 35953237722 : iptr = nhc*(ktmp + j - 1)
1924 >14550*10^7 : do i=1,nhc
1925 >12837*10^7 : edge%buf(iptr+inw+i)=v(i ,npoints-j+1,k)
1926 : enddo
1927 : end do
1928 : end do
1929 : end if
1930 : end do
1931 :
1932 138018600 : end subroutine ghostpack
1933 :
1934 138018600 : subroutine ghostunpack(edge,v,vlyr,kptr,ielem)
1935 : use dimensions_mod, only : max_corner_elem
1936 : use control_mod, only : north, south, east, west, neast, nwest, seast, swest
1937 : type (Edgebuffer_t), intent(in) :: edge
1938 :
1939 : integer, intent(in) :: vlyr
1940 : integer, intent(in) :: kptr
1941 : integer, intent(in) :: ielem
1942 :
1943 : real (kind=r8), intent(inout) :: v(edge%lb:edge%ub,edge%lb:edge%ub,vlyr)
1944 :
1945 :
1946 : ! Local
1947 : logical, parameter :: UseUnroll = .TRUE.
1948 : integer :: i,j,k,l,itr, ktmp
1949 : integer :: is,ie,in,iw,isw,ise,inw,ine
1950 : integer :: nhc,npoints,iptr
1951 : logical :: reverse
1952 :
1953 138018600 : threadsafe=.false.
1954 :
1955 138018600 : is=edge%getmap(south,ielem)
1956 138018600 : ie=edge%getmap(east,ielem)
1957 138018600 : in=edge%getmap(north,ielem)
1958 138018600 : iw=edge%getmap(west,ielem)
1959 :
1960 138018600 : nhc = edge%ndepth
1961 138018600 : npoints = edge%npoints
1962 :
1963 : ! example for north buffer
1964 : ! first row ('edge') goes in v(:,np+1,k)
1965 : ! 2nd row ('edge') goes in v(:,np+2,k)
1966 : ! etc...
1967 : !DIR$ IVDEP
1968 17289007200 : do k=1,vlyr
1969 17150988600 : ktmp = nhc*(kptr+k-1)
1970 53282237400 : do j=1,nhc
1971 35993230200 : iptr = npoints*(ktmp + j - 1)
1972 >16112*10^7 : do i=1,npoints
1973 >10797*10^7 : v(i ,1-j ,k) = edge%receive(iptr+is+i) ! South
1974 >10797*10^7 : v(npoints+j ,i ,k) = edge%receive(iptr+ie+i) ! East
1975 >10797*10^7 : v(i ,npoints+j ,k) = edge%receive(iptr+in+i) ! North
1976 >14397*10^7 : v(1-j ,i ,k) = edge%receive(iptr+iw+i) ! West
1977 : end do
1978 : end do
1979 : end do
1980 :
1981 :
1982 : ! SWEST
1983 276037200 : do l=swest,swest+max_corner_elem-1
1984 138018600 : isw = edge%getmap(l,ielem)
1985 276037200 : if(isw /= -1) then
1986 : ! note the following is the the correct meaning of reverse in this code.
1987 : ! It is best described as a transponse operation
1988 137865246 : if (edge%reverse(l,ielem)) then
1989 371393024 : do k=1,vlyr
1990 368428412 : ktmp = nhc*(kptr+k-1)
1991 1144580468 : do j=1,nhc
1992 773187444 : iptr = nhc*(ktmp + j - 1)
1993 3129080396 : do i=1,nhc
1994 2760651984 : v(1-j,1-i,k)=edge%receive(iptr+isw+i)
1995 : enddo
1996 : enddo
1997 : enddo
1998 : else
1999 16898404168 : do k=1,vlyr
2000 16763503534 : ktmp = nhc*(kptr+k-1)
2001 52078454446 : do i=1,nhc
2002 35180050278 : iptr = nhc*(ktmp + i - 1)
2003 >14237*10^7 : do j=1,nhc
2004 >12560*10^7 : v(1-j,1-i,k)=edge%receive(iptr+isw+j)
2005 : enddo
2006 : enddo
2007 : enddo
2008 : endif
2009 : else
2010 19210008 : do k=1,vlyr
2011 59202486 : do j=1,nhc
2012 161849070 : do i=1,nhc
2013 142792416 : v(1-i,1-j,k)=edgeDefaultVal
2014 : enddo
2015 : enddo
2016 : enddo
2017 : endif
2018 : end do
2019 :
2020 : ! SEAST
2021 276037200 : do l=swest+max_corner_elem,swest+2*max_corner_elem-1
2022 138018600 : ise = edge%getmap(l,ielem)
2023 276037200 : if(ise /= -1) then
2024 137865246 : if (edge%reverse(l,ielem)) then
2025 371393024 : do k=1,vlyr
2026 368428412 : ktmp = nhc*(kptr+k-1)
2027 1144580468 : do i=1,nhc
2028 773187444 : iptr = nhc*(ktmp + i - 1)
2029 3129080396 : do j=1,nhc
2030 2760651984 : v(npoints+i,1-j,k)=edge%receive(iptr+ise+j)
2031 : enddo
2032 : enddo
2033 : enddo
2034 : else
2035 16898404168 : do k=1,vlyr
2036 16763503534 : ktmp = nhc*(kptr+k-1)
2037 52078454446 : do j=1,nhc
2038 35180050278 : iptr = nhc*(ktmp + j - 1)
2039 >14237*10^7 : do i=1,nhc
2040 >12560*10^7 : v(npoints+i ,1-j ,k)=edge%receive(iptr+ise+i)
2041 : enddo
2042 : enddo
2043 : enddo
2044 : endif
2045 : else
2046 19210008 : do k=1,vlyr
2047 59202486 : do j=1,nhc
2048 161849070 : do i=1,nhc
2049 142792416 : v(npoints+i,1-j,k)=edgeDefaultVal
2050 : enddo
2051 : enddo
2052 : enddo
2053 : endif
2054 : end do
2055 :
2056 : ! NEAST
2057 276037200 : do l=swest+3*max_corner_elem,swest+4*max_corner_elem-1
2058 138018600 : ine = edge%getmap(l,ielem)
2059 276037200 : if(ine /= -1) then
2060 137865246 : if (edge%reverse(l,ielem)) then
2061 371393024 : do k=1,vlyr
2062 368428412 : ktmp = nhc*(kptr+k-1)
2063 1144580468 : do j=1,nhc
2064 3129080396 : do i=1,nhc
2065 1987464540 : iptr = nhc*(ktmp + i - 1)
2066 2760651984 : v(npoints+i ,npoints+j,k)=edge%receive(iptr+ine+j)
2067 : enddo
2068 : enddo
2069 : enddo
2070 : else
2071 16898404168 : do k=1,vlyr
2072 16763503534 : ktmp = nhc*(kptr+k-1)
2073 52078454446 : do j=1,nhc
2074 35180050278 : iptr = nhc*(ktmp + j - 1)
2075 >14237*10^7 : do i=1,nhc
2076 >12560*10^7 : v(npoints+i ,npoints+j,k)=edge%receive(iptr+ine+i)
2077 : enddo
2078 : enddo
2079 : enddo
2080 : endif
2081 : else
2082 19210008 : do k=1,vlyr
2083 59202486 : do j=1,nhc
2084 161849070 : do i=1,nhc
2085 142792416 : v(npoints+i,npoints+j,k)=edgeDefaultVal
2086 : enddo
2087 : enddo
2088 : enddo
2089 : endif
2090 : end do
2091 :
2092 : ! NWEST
2093 276037200 : do l=swest+2*max_corner_elem,swest+3*max_corner_elem-1
2094 138018600 : inw = edge%getmap(l,ielem)
2095 276037200 : if(inw /= -1) then
2096 137865246 : if (edge%reverse(l,ielem)) then
2097 371393024 : do k=1,vlyr
2098 368428412 : ktmp = nhc*(kptr+k-1)
2099 1144580468 : do i=1,nhc
2100 773187444 : iptr = nhc*(ktmp + i - 1)
2101 3129080396 : do j=1,nhc
2102 2760651984 : v(1-i ,npoints+j,k)=edge%receive(iptr+inw+j)
2103 : enddo
2104 : enddo
2105 : enddo
2106 : else
2107 16898404168 : do k=1,vlyr
2108 16763503534 : ktmp = nhc*(kptr+k-1)
2109 52078454446 : do j=1,nhc
2110 35180050278 : iptr = nhc*(ktmp + j - 1)
2111 >14237*10^7 : do i=1,nhc
2112 >12560*10^7 : v(1-i ,npoints+j,k)=edge%receive(iptr+inw+i)
2113 : enddo
2114 : enddo
2115 : enddo
2116 : endif
2117 : else
2118 19210008 : do k=1,vlyr
2119 59202486 : do j=1,nhc
2120 161849070 : do i=1,nhc
2121 142792416 : v(1-i,npoints+j,k)=edgeDefaultVal
2122 : enddo
2123 : enddo
2124 : enddo
2125 : endif
2126 : end do
2127 :
2128 138018600 : end subroutine ghostunpack
2129 :
2130 : ! =========================================
2131 : ! initGhostBuffer3d:
2132 : ! Author: James Overfelt
2133 : ! create an Real based communication buffer
2134 : ! npoints is the number of points on one side
2135 : ! nhc is the deep of the ghost/halo zone
2136 : ! =========================================
2137 0 : subroutine initGhostBuffer3d(ghost,nlyr,np,nhc_in)
2138 :
2139 : implicit none
2140 : integer,intent(in) :: nlyr, np
2141 : integer,intent(in),optional :: nhc_in
2142 : type (Ghostbuffer3d_t),intent(out) :: ghost
2143 :
2144 : ! Local variables
2145 :
2146 : integer :: nbuf,nhc,i
2147 :
2148 : ! sanity check for threading
2149 0 : if (omp_get_num_threads()>1) then
2150 0 : call endrun('ERROR: initGhostBuffer must be called before threaded region')
2151 : endif
2152 :
2153 0 : if (present(nhc_in)) then
2154 0 : nhc=nhc_in
2155 : else
2156 0 : nhc = np-1
2157 : endif
2158 :
2159 0 : nbuf=max_neigh_edges*nelemd
2160 :
2161 0 : ghost%nlyr = nlyr
2162 0 : ghost%nhc = nhc
2163 0 : ghost%np = np
2164 0 : ghost%nbuf = nbuf
2165 0 : ghost%elem_size = np*(nhc+1)
2166 0 : allocate(ghost%buf (np,(nhc+1),nlyr,nbuf))
2167 0 : allocate(ghost%receive(np,(nhc+1),nlyr,nbuf))
2168 0 : ghost%buf=0
2169 0 : ghost%receive=0
2170 :
2171 0 : end subroutine initGhostBuffer3d
2172 :
2173 : ! =================================================================================
2174 : ! GHOSTVPACK3D
2175 : ! AUTHOR: James Overfelt (from a subroutine of Christoph Erath, ghostvpack2D)
2176 : ! Pack edges of v into an ghost buffer for boundary exchange.
2177 : !
2178 : ! This subroutine packs for many vertical layers into an ghost
2179 : ! buffer.
2180 : ! If the buffer associated with edge is not large enough to
2181 : ! hold all vertical layers you intent to pack, the method will
2182 : ! halt the program with a call to endrun().
2183 : ! INPUT:
2184 : ! - ghost Buffer into which the data will be packed.
2185 : ! This buffer must be previously allocated with initGhostBuffer().
2186 : ! - v The data to be packed.
2187 : ! - nhc deep of ghost/halo zone
2188 : ! - npoints number of points on on side
2189 : ! - kptr Vertical pointer to the place in the edge buffer where
2190 : ! data will be located.
2191 : ! =================================================================================
2192 0 : subroutine ghostVpack3d(ghost, v, vlyr, kptr, desc)
2193 : use dimensions_mod, only : max_corner_elem
2194 : use control_mod, only : north, south, east, west, neast, nwest, seast, swest
2195 : use edgetype_mod, only : edgedescriptor_t, ghostbuffer3d_t
2196 : implicit none
2197 :
2198 : type (Ghostbuffer3d_t) :: ghost
2199 : integer, intent(in) :: kptr,vlyr
2200 : real (kind=r8),intent(in) :: v(ghost%np, ghost%np, vlyr)
2201 : type (EdgeDescriptor_t),intent(in) :: desc
2202 :
2203 : integer :: nhc, np
2204 :
2205 : ! Local variables
2206 : integer :: i,j,k,ir,l,e
2207 :
2208 : integer :: is,ie,in,iw
2209 :
2210 0 : if(.not. threadsafe) then
2211 : !$OMP BARRIER
2212 0 : threadsafe=.true.
2213 : end if
2214 : ! Example convenction for buffer to the north:
2215 : ! buf(:,,:,i,e)
2216 : ! each "edge" is a row of data (i=1,np) in the element
2217 : ! north most row of data goes into e=1
2218 : ! next row of data goes into e=2
2219 : ! ....
2220 : ! south most row of data goes into e=np
2221 : ! We need to pack this way to preserve the orientation
2222 : ! so the data can be unpacked correctly
2223 :
2224 : ! note: we think of buf as dimensioned buf(k,is,i,e)
2225 : ! but this array is flatted to: buf(k,is+(i-1)+(e-1)*np)
2226 : !
2227 0 : nhc = ghost%nhc
2228 0 : np = ghost%np
2229 0 : is = desc%putmapP_ghost(south)
2230 0 : ie = desc%putmapP_ghost(east)
2231 0 : in = desc%putmapP_ghost(north)
2232 0 : iw = desc%putmapP_ghost(west)
2233 :
2234 0 : do k=1,vlyr
2235 0 : do j=1,nhc
2236 0 : do i=1,np
2237 0 : ghost%buf(i,j,kptr+k,is) = v(i, j+1 , k)
2238 0 : ghost%buf(i,j,kptr+k,ie) = v(np-j, i , k)
2239 0 : ghost%buf(i,j,kptr+k,in) = v(i, np-j , k)
2240 0 : ghost%buf(i,j,kptr+k,iw) = v(j+1, i , k)
2241 : enddo
2242 : end do
2243 : end do
2244 : ! This is really kludgy way to setup the index reversals
2245 : ! But since it is so a rare event not real need to spend time optimizing
2246 : ! Check if the edge orientation of the recieving element is different
2247 : ! if it is, swap the order of data in the edge
2248 0 : if(desc%reverse(south)) then
2249 0 : do k=1,vlyr
2250 0 : do j=1,nhc
2251 0 : do i=1,np
2252 0 : ir = np-i+1
2253 0 : ghost%buf(ir, j, kptr+k, is)=v(i, j+1, k)
2254 : enddo
2255 : enddo
2256 : enddo
2257 : endif
2258 :
2259 0 : if(desc%reverse(east)) then
2260 0 : do k=1,vlyr
2261 0 : do j=1,nhc
2262 0 : do i=1,np
2263 0 : ir = np-i+1
2264 0 : ghost%buf(ir, j, kptr+k, ie)=v(np-j, i, k)
2265 : enddo
2266 : enddo
2267 : enddo
2268 : endif
2269 :
2270 0 : if(desc%reverse(north)) then
2271 0 : do k=1,vlyr
2272 0 : do j=1,nhc
2273 0 : do i=1,np
2274 0 : ir = np-i+1
2275 0 : ghost%buf(ir, j, kptr+k, in)=v(i, np-j, k)
2276 : enddo
2277 : enddo
2278 : enddo
2279 : endif
2280 :
2281 0 : if(desc%reverse(west)) then
2282 0 : do k=1,vlyr
2283 0 : do j=1,nhc
2284 0 : do i=1,np
2285 0 : ir = np-i+1
2286 0 : ghost%buf(ir, j, kptr+k, iw)=v(j+1, i, k)
2287 : enddo
2288 : enddo
2289 : enddo
2290 : endif
2291 :
2292 : ! corners. this is difficult because we dont know the orientaton
2293 : ! of the corners, and this which (i,j) dimension maps to which dimension
2294 : ! SWEST
2295 0 : do l=swest, swest+max_corner_elem-1
2296 0 : if (desc%putmapP_ghost(l) /= -1) then
2297 0 : do k=1,vlyr
2298 0 : do j=1,nhc+1
2299 0 : do i=1,nhc+1
2300 0 : ghost%buf(i, j, kptr+k, desc%putmapP_ghost(l))=v(i, j, k)
2301 : enddo
2302 : enddo
2303 : enddo
2304 : end if
2305 : end do
2306 :
2307 : ! SEAST
2308 0 : do l=swest+max_corner_elem,swest+2*max_corner_elem-1
2309 0 : if (desc%putmapP_ghost(l) /= -1) then
2310 0 : do k=1,vlyr
2311 0 : do j=1,nhc+1
2312 0 : do i=1,nhc+1
2313 0 : ghost%buf(i, j, kptr+k, desc%putmapP_ghost(l))=v(np-i+1, j, k)
2314 : enddo
2315 : enddo
2316 : enddo
2317 : end if
2318 : end do
2319 :
2320 : ! NEAST
2321 0 : do l=swest+3*max_corner_elem,swest+4*max_corner_elem-1
2322 0 : if (desc%putmapP_ghost(l) /= -1) then
2323 0 : do k=1,vlyr
2324 0 : do j=1,nhc+1
2325 0 : do i=1,nhc+1
2326 0 : ghost%buf(i, j, kptr+k,desc%putmapP_ghost(l))=v(np-i+1, np-j+1, k)
2327 : enddo
2328 : enddo
2329 : enddo
2330 : end if
2331 : end do
2332 :
2333 : ! NWEST
2334 0 : do l=swest+2*max_corner_elem,swest+3*max_corner_elem-1
2335 0 : if (desc%putmapP_ghost(l) /= -1) then
2336 0 : do k=1,vlyr
2337 0 : do j=1,nhc+1
2338 0 : do i=1,nhc+1
2339 0 : ghost%buf(i, j, kptr+k,desc%putmapP_ghost(l))=v(i, np-j+1, k)
2340 : enddo
2341 : enddo
2342 : enddo
2343 : end if
2344 : end do
2345 0 : end subroutine ghostVpack3d
2346 :
2347 : ! =================================================================================
2348 : ! GHOSTVUNPACK3D
2349 : ! AUTHOR: James Overfelt (from a subroutine of Christoph Erath,
2350 : ! ghostVunpack2d)
2351 : ! Unpack ghost points from ghost buffer into v...
2352 : ! It is for cartesian points (v is only two dimensional).
2353 : ! INPUT SAME arguments as for GHOSTVPACK
2354 : ! =================================================================================
2355 :
2356 0 : subroutine ghostVunpack3d(g, v, vlyr, kptr, desc, sw, se, nw, ne, mult)
2357 : use dimensions_mod, only : max_corner_elem
2358 : use control_mod, only : north, south, east, west, neast, nwest, seast, swest
2359 : use edgetype_mod, only : edgedescriptor_t, ghostbuffer3d_t
2360 : implicit none
2361 : type (Ghostbuffer3d_t), intent(in) :: g
2362 :
2363 : integer, intent(in) :: kptr,vlyr
2364 : real (kind=r8), intent(inout) :: v (1-g%nhc : g%np+g%nhc, 1-g%nhc : g%np+g%nhc, vlyr)
2365 : integer, intent(out) :: mult(5:8)
2366 : real (kind=r8), intent(out) :: sw(1-g%nhc : 1, 1-g%nhc : 1, vlyr, max_corner_elem-1)
2367 : real (kind=r8), intent(out) :: se( g%np : g%np+g%nhc, 1-g%nhc : 1, vlyr, max_corner_elem-1)
2368 : real (kind=r8), intent(out) :: ne( g%np : g%np+g%nhc, g%np : g%np+g%nhc, vlyr, max_corner_elem-1)
2369 : real (kind=r8), intent(out) :: nw(1-g%nhc : 1, g%np : g%np+g%nhc, vlyr, max_corner_elem-1)
2370 : type (EdgeDescriptor_t) :: desc
2371 :
2372 : integer :: nhc, np
2373 :
2374 : ! Local
2375 : logical, parameter :: UseUnroll = .TRUE.
2376 : integer :: i,j,k,l
2377 : integer :: is,ie,in,iw,ic
2378 : logical :: reverse
2379 :
2380 0 : threadsafe=.false.
2381 :
2382 0 : nhc = g%nhc
2383 0 : np = g%np
2384 :
2385 0 : is=desc%getmapP_ghost(south)
2386 0 : ie=desc%getmapP_ghost(east)
2387 0 : in=desc%getmapP_ghost(north)
2388 0 : iw=desc%getmapP_ghost(west)
2389 :
2390 : ! fill in optional values with edgeDefaultVal
2391 0 : do k=1,vlyr
2392 0 : do j=1,nhc
2393 0 : do i=1,nhc
2394 0 : v(1-i, 1-j, k)=edgeDefaultVal
2395 0 : v(np+i , 1-j, k)=edgeDefaultVal
2396 0 : v(np+i, np+j, k)=edgeDefaultVal
2397 0 : v(1-i , np+j, k)=edgeDefaultVal
2398 : enddo
2399 : enddo
2400 : enddo
2401 :
2402 : ! example for north buffer
2403 : ! first row ('edge') goes in v(:,np+1)
2404 : ! 2nd row ('edge') goes in v(:,np+2)
2405 : ! etc...
2406 :
2407 0 : do k=1,vlyr
2408 0 : do j=1,nhc
2409 0 : do i=1,np
2410 0 : v(i , 1-j , k) = g%buf(i,j,kptr+k,is )
2411 0 : v(np+j , i , k) = g%buf(i,j,kptr+k,ie )
2412 0 : v(i , np+j, k) = g%buf(i,j,kptr+k,in )
2413 0 : v(1-j , i , k) = g%buf(i,j,kptr+k,iw )
2414 : end do
2415 : end do
2416 : end do
2417 :
2418 : ! four sides are always just one
2419 0 : mult(swest) = 0
2420 0 : mult(seast) = 0
2421 0 : mult(neast) = 0
2422 0 : mult(nwest) = 0
2423 :
2424 :
2425 :
2426 : ! SWEST
2427 0 : do l=swest, swest+max_corner_elem-1
2428 0 : ic = desc%getmapP_ghost(l)
2429 0 : if(ic /= -1) then
2430 0 : reverse=desc%reverse(l)
2431 0 : if (mult(swest) .eq. 0) then
2432 0 : if (reverse) then
2433 0 : do k=1,vlyr
2434 0 : do j=1,nhc
2435 0 : do i=1,nhc
2436 0 : v(1-i, 1-j, k)=g%buf(j+1, i+1, kptr+k, ic)
2437 : enddo
2438 : enddo
2439 : enddo
2440 : else
2441 0 : do k=1,vlyr
2442 0 : do j=1,nhc
2443 0 : do i=1,nhc
2444 0 : v(1-i,1-j,k)=g%buf(i+1,j+1,kptr+k,ic)
2445 : enddo
2446 : enddo
2447 : enddo
2448 : endif
2449 : else
2450 0 : if (reverse) then
2451 0 : do k=1,vlyr
2452 0 : do j=0,nhc
2453 0 : do i=0,nhc
2454 0 : sw(1-i,1-j,k,mult(swest))=g%buf(j+1,i+1,kptr+k,ic)
2455 : enddo
2456 : enddo
2457 : enddo
2458 : else
2459 0 : do k=1,vlyr
2460 0 : do j=0,nhc
2461 0 : do i=0,nhc
2462 0 : sw(1-i,1-j,k,mult(swest))=g%buf(i+1,j+1,kptr+k,ic)
2463 : enddo
2464 : enddo
2465 : enddo
2466 : endif
2467 : endif
2468 0 : mult(swest) = mult(swest) + 1
2469 : endif
2470 : end do
2471 :
2472 : ! SEAST
2473 0 : do l=swest+max_corner_elem,swest+2*max_corner_elem-1
2474 0 : ic = desc%getmapP_ghost(l)
2475 0 : if(ic /= -1) then
2476 0 : reverse=desc%reverse(l)
2477 0 : if (mult(seast) .eq. 0) then
2478 0 : if (reverse) then
2479 0 : do k=1,vlyr
2480 0 : do j=1,nhc
2481 0 : do i=1,nhc
2482 0 : v(np+i,1-j,k)=g%buf(j+1,i+1,kptr+k,ic)
2483 : enddo
2484 : enddo
2485 : enddo
2486 : else
2487 0 : do k=1,vlyr
2488 0 : do j=1,nhc
2489 0 : do i=1,nhc
2490 0 : v(np+i ,1-j,k)=g%buf(i+1,j+1,kptr+k,ic)
2491 : enddo
2492 : enddo
2493 : enddo
2494 : endif
2495 : else
2496 0 : if (reverse) then
2497 0 : do k=1,vlyr
2498 0 : do j=0,nhc
2499 0 : do i=0,nhc
2500 0 : se(np+i,1-j,k,mult(seast))=g%buf(j+1,i+1,kptr+k,ic)
2501 : enddo
2502 : enddo
2503 : enddo
2504 : else
2505 0 : do k=1,vlyr
2506 0 : do j=0,nhc
2507 0 : do i=0,nhc
2508 0 : se(np+i ,1-j,k,mult(seast))=g%buf(i+1,j+1,kptr+k,ic)
2509 : enddo
2510 : enddo
2511 : enddo
2512 : endif
2513 : endif
2514 0 : mult(seast) = mult(seast) + 1
2515 : endif
2516 : end do
2517 :
2518 :
2519 : ! NEAST
2520 0 : do l=swest+3*max_corner_elem,swest+4*max_corner_elem-1
2521 0 : ic = desc%getmapP_ghost(l)
2522 0 : if(ic /= -1) then
2523 0 : reverse=desc%reverse(l)
2524 0 : if (mult(neast) .eq. 0) then
2525 0 : if (reverse) then
2526 0 : do k=1,vlyr
2527 0 : do j=1,nhc
2528 0 : do i=1,nhc
2529 0 : v(np+i ,np+j,k)=g%buf(j+1,i+1,kptr+k,ic)
2530 : enddo
2531 : enddo
2532 : enddo
2533 : else
2534 0 : do k=1,vlyr
2535 0 : do j=1,nhc
2536 0 : do i=1,nhc
2537 0 : v(np+i ,np+j,k)=g%buf(i+1,j+1,kptr+k,ic)
2538 : enddo
2539 : enddo
2540 : enddo
2541 : endif
2542 : else
2543 0 : if (reverse) then
2544 0 : do k=1,vlyr
2545 0 : do j=0,nhc
2546 0 : do i=0,nhc
2547 0 : ne(np+i ,np+j,k,mult(neast))=g%buf(j+1,i+1,kptr+k,ic)
2548 : enddo
2549 : enddo
2550 : enddo
2551 : else
2552 0 : do k=1,vlyr
2553 0 : do j=0,nhc
2554 0 : do i=0,nhc
2555 0 : ne(np+i ,np+j,k,mult(neast))=g%buf(i+1,j+1,kptr+k,ic)
2556 : enddo
2557 : enddo
2558 : enddo
2559 : endif
2560 : endif
2561 0 : mult(neast) = mult(neast) + 1
2562 : endif
2563 : end do
2564 :
2565 : ! NWEST
2566 0 : do l=swest+2*max_corner_elem,swest+3*max_corner_elem-1
2567 0 : ic = desc%getmapP_ghost(l)
2568 0 : if(ic /= -1) then
2569 0 : reverse=desc%reverse(l)
2570 0 : if (mult(nwest) .eq. 0) then
2571 0 : if (reverse) then
2572 0 : do k=1,vlyr
2573 0 : do j=1,nhc
2574 0 : do i=1,nhc
2575 0 : v(1-i ,np+j,k)=g%buf(j+1,i+1,kptr+k,ic)
2576 : enddo
2577 : enddo
2578 : enddo
2579 : else
2580 0 : do k=1,vlyr
2581 0 : do j=1,nhc
2582 0 : do i=1,nhc
2583 0 : v(1-i ,np+j,k)=g%buf(i+1,j+1,kptr+k,ic)
2584 : enddo
2585 : enddo
2586 : enddo
2587 : endif
2588 : else
2589 0 : if (reverse) then
2590 0 : do k=1,vlyr
2591 0 : do j=0,nhc
2592 0 : do i=0,nhc
2593 0 : nw(1-i ,np+j,k,mult(nwest))=g%buf(j+1,i+1,kptr+k,ic)
2594 : enddo
2595 : enddo
2596 : enddo
2597 : else
2598 0 : do k=1,vlyr
2599 0 : do j=0,nhc
2600 0 : do i=0,nhc
2601 0 : nw(1-i ,np+j,k,mult(nwest))=g%buf(i+1,j+1,kptr+k,ic)
2602 : enddo
2603 : enddo
2604 : enddo
2605 : endif
2606 : endif
2607 0 : mult(nwest) = mult(nwest) + 1
2608 : endif
2609 : end do
2610 :
2611 0 : end subroutine ghostVunpack3d
2612 :
2613 0 : subroutine FreeGhostBuffer3D(buffer)
2614 : use edgetype_mod, only : ghostbuffer3d_t
2615 : implicit none
2616 : type (Ghostbuffer3d_t),intent(inout) :: buffer
2617 :
2618 : !$OMP BARRIER
2619 : !$OMP MASTER
2620 0 : buffer%nbuf=0
2621 0 : buffer%nlyr=0
2622 0 : deallocate(buffer%buf)
2623 0 : deallocate(buffer%receive)
2624 : !$OMP END MASTER
2625 :
2626 0 : end subroutine FreeGhostBuffer3D
2627 :
2628 :
2629 0 : End module edge_mod
|
