File: /nfs/home/0/users/jenkins/mfix.git/model/dmp_modules/gridmap_mod.f

1     !--------------------------------------------------------------------
2     ! Purpose:
3     ! Contains following subroutines:
4     !    partition, gridmap_init
5     !--------------------------------------------------------------------
6            module gridmap
7     
8     !-----------------------------------------------
9     ! Modules
10     !-----------------------------------------------
11             use mpi_utility
12             use parallel_mpi
13             use geometry
14             use sendrecv
15             use compar
16             use run
17             use indices
18     
19             use error_manager
20     
21             implicit none
22     
23             contains
24     
25     !----------------------------------------------------------------------!
26     ! Purpose: Routine to partition the grid. It works for 1-d, 2-d        !
27     ! decomposition in the current implementation                          !
28     !----------------------------------------------------------------------!
29           SUBROUTINE PARTITION(CYC_XLL, CYC_YLL, CYC_ZLL)
30     
31           implicit none
32     
33     ! DUMMY Arguments
34     !---------------------------------------------------------------------//
35     ! Local flags for cyclic boundarys
36            LOGICAL :: CYC_XLL, CYC_YLL, CYC_ZLL
37     
38     ! Local variables
39     !---------------------------------------------------------------------//
40           INTEGER, dimension(0:nodesi-1) :: isize1_all
41           INTEGER, dimension(0:nodesj-1) :: jsize1_all
42           INTEGER, dimension(0:nodesk-1) :: ksize1_all
43     
44           INTEGER :: ip, iproc, isize, iremain
45           INTEGER :: jp, jproc, jsize, jremain
46           INTEGER :: kp, kproc, ksize, kremain
47     
48           INTEGER :: ijkproc
49     
50           LOGICAL :: PRESENT
51     
52     !......................................................................!
53     
54     
55     ! Initialize the error manager.
56           CALL INIT_ERR_MSG("PARTITION")
57     
58     ! Set the number of layers for BICGSTAB
59           IF(NODESI .NE. 1 .AND. CYC_XLL) nlayers_bicgs = 2
60           IF(NODESJ .NE. 1 .AND. CYC_YLL) nlayers_bicgs = 2
61           IF(NODESK .NE. 1 .AND. CYC_ZLL) nlayers_bicgs = 2
62     
63     ! Flag that the current setup may not be efficient.
64           IF(NODESJ .NE. 1) THEN
65              WRITE(ERR_MSG,1000)
66              CALL FLUSH_ERR_MSG
67           ENDIF
68     
69     ! Get Domain size from ADJUST_IJK_SIZE Subroutine
70           IF(DOMAIN_SIZE_ADJUSTED) THEN
71              isize1_all = isize_all
72              jsize1_all = jsize_all
73              ksize1_all = ksize_all
74           ELSE
75     ! Determine the size in i direction and add the remainder sequentially
76              isize = (imax1-imin1+1)/nodesi
77              isize1_all(0:nodesi-1) = isize
78              iremain = (imax1-imin1+1) - nodesi*isize
79              IF (iremain.ge.1) isize1_all( 0:(iremain-1) ) = isize + 1
80     
81     ! Determine the size in j direction and add the remainder sequentially
82              jsize = (jmax1-jmin1+1)/nodesj
83              jsize1_all(0:nodesj-1) = jsize
84              jremain = (jmax1-jmin1+1) - nodesj*jsize
85              IF (jremain.ge.1) jsize1_all( 0:(jremain-1) ) = jsize + 1
86     
87     ! Determine the size in k direction and add the remainder sequentially
88              ksize = (kmax1-kmin1+1)/nodesk
89              ksize1_all(0:nodesk-1) = ksize
90              kremain = (kmax1-kmin1+1) - nodesk*ksize
91              IF (kremain.ge.1) ksize1_all( 0:(kremain-1) ) = ksize + 1
92           ENDIF
93     
94     
95     ! Get Domain size from gridmap.dat
96     ! This works only in the j-direction   <-------------------------
97           IF(.NOT.DOMAIN_SIZE_ADJUSTED) THEN
98              INQUIRE(FILE='gridmap.dat',EXIST=PRESENT)
99              IF(PRESENT) THEN
100                 IF(MyPE == PE_IO) THEN
101                    WRITE(*,*)'Reading gridmap from grimap.dat...'
102                    OPEN(UNIT=777, FILE='gridmap.dat', STATUS='OLD')
103     
104                    READ (777, *) NODESI,NODESJ,NODESK
105                    DO IPROC = 0,NODESI-1
106                       READ(777,*) jPROC,Isize1_all(IPROC)
107                    ENDDO
108                    DO IPROC = 0,NODESJ-1
109                       READ(777,*) jPROC,Jsize1_all(IPROC)
110                    ENDDO
111                    DO IPROC = 0,NODESK-1
112                       READ(777,*) jPROC,Ksize1_all(IPROC)
113                    ENDDO
114     
115                    CLOSE(777)
116                 ENDIF
117                 CALL BCAST(ISIZE1_ALL)
118                 CALL BCAST(JSIZE1_ALL)
119                 CALL BCAST(KSIZE1_ALL)
120                 allocate( ISIZE_ALL(0:NODESI-1))
121                 allocate( JSIZE_ALL(0:NODESJ-1))
122                 allocate( KSIZE_ALL(0:NODESK-1))
123                 isize_all = isize1_all
124                 jsize_all = jsize1_all
125                 ksize_all = ksize1_all
126              ENDIF
127           ENDIF
128     
129     ! The following is general for 1-d or 2-d or 3-d decompostion
130     ! Determining  istart, jstart and kstart for all the processors
131           ijkproc = 0
132           kp = kmin1
133           do kproc=0,nodesk-1
134              jp = jmin1
135              do jproc=0,nodesj-1
136                 ip = imin1
137                 do iproc=0,nodesi-1
138     
139                    istart1_all(ijkproc) = ip + sum(isize1_all(0:iproc-1))
140                    iend1_all(ijkproc) = istart1_all(ijkproc) + isize1_all(iproc)-1
141     
142                    jstart1_all(ijkproc) = jp + sum(jsize1_all(0:jproc-1))
143                    jend1_all(ijkproc) = jstart1_all(ijkproc) + jsize1_all(jproc)-1
144     
145                    kstart1_all(ijkproc) = kp + sum(ksize1_all(0:kproc-1))
146                    kend1_all(ijkproc) = kstart1_all(ijkproc) + ksize1_all(kproc)-1
147     
148                    ijkproc = ijkproc+1
149     
150                 ENDDO
151              ENDDO
152           ENDDO
153     
154           CALL FINL_ERR_MSG
155     
156           RETURN
157     
158      1000 FORMAT('WARNING 1000: The preconditioner for the linear solver',/&
159              'MIGHT NOT be very efficient with DMP partitions in the y-',  &
160              'axis.')
161     
162           END SUBROUTINE PARTITION
163     
164     
165     !----------------------------------------------------------------------!
166     ! Purpose: Initializing all the variables from the information         !
167     ! obtained in the above routine.                                       !
168     !----------------------------------------------------------------------!
169             SUBROUTINE GRIDMAP_INIT
170     
171             use functions
172             use toleranc
173     
174             implicit none
175     
176     ! Local variables
177     !---------------------------------------------------------------------//
178     ! Loop indicies
179           integer :: iproc, ijk, ii, jj, kk, idebug
180     ! Local flags for cyclic boundarys
181           LOGICAL :: CYC_XL, CYC_YL, CYC_ZL
182     ! Communicator (MPI_COMM_WORLD)
183           INTEGER :: COMM
184     ! Amount of load imbalance
185           INTEGER :: IMBALANCE
186     ! Theoritical speedup (based on load imbalance)
187           CHARACTER(len=32) :: AMDAHL_SPEEDUP
188     
189     !......................................................................!
190     
191     ! Initialize the error manager.
192           CALL INIT_ERR_MSG("GRIDMAP_INIT")
193     
194     ! Set local flags for cyclic boundaries.
195           CYC_XL = (CYCLIC_X .OR. CYCLIC_X_PD)
196           CYC_YL = (CYCLIC_Y .OR. CYCLIC_Y_PD)
197           CYC_ZL = (CYCLIC_Z .OR. CYCLIC_Z_PD)
198           IF(DO_K .AND. COMPARE(ZLENGTH,8.D0*ATAN(ONE))  .AND. &
199              (COORDINATES == 'CYLINDRICAL')) CYC_ZL = .TRUE.
200     
201     
202           IF(.NOT.ALLOCATED(ijksize3_all))   allocate( ijksize3_all(0:numPEs-1) )
203           IF(.NOT.ALLOCATED(ijkstart3_all))  allocate( ijkstart3_all(0:numPEs-1) )
204           IF(.NOT.ALLOCATED(ijkend3_all))    allocate( ijkend3_all(0:numPEs-1) )
205     
206           IF(.NOT.ALLOCATED(ijksize4_all))   allocate( ijksize4_all(0:numPEs-1) )
207           IF(.NOT.ALLOCATED(ijkstart4_all))  allocate( ijkstart4_all(0:numPEs-1) )
208           IF(.NOT.ALLOCATED(ijkend4_all))    allocate( ijkend4_all(0:numPEs-1) )
209     
210           IF(.NOT.ALLOCATED(istart_all))     allocate( istart_all(0:numPEs-1) )
211           IF(.NOT.ALLOCATED(jstart_all))     allocate( jstart_all(0:numPEs-1) )
212           IF(.NOT.ALLOCATED(kstart_all))     allocate( kstart_all(0:numPEs-1) )
213     
214           IF(.NOT.ALLOCATED(istart1_all))    allocate( istart1_all(0:numPEs-1) )
215           IF(.NOT.ALLOCATED(jstart1_all))    allocate( jstart1_all(0:numPEs-1) )
216           IF(.NOT.ALLOCATED(kstart1_all))    allocate( kstart1_all(0:numPEs-1) )
217     
218           IF(.NOT.ALLOCATED(istart2_all))    allocate( istart2_all(0:numPEs-1) )
219           IF(.NOT.ALLOCATED(jstart2_all))    allocate( jstart2_all(0:numPEs-1) )
220           IF(.NOT.ALLOCATED(kstart2_all))    allocate( kstart2_all(0:numPEs-1) )
221     
222           IF(.NOT.ALLOCATED(istart3_all))    allocate( istart3_all(0:numPEs-1) )
223           IF(.NOT.ALLOCATED(jstart3_all))    allocate( jstart3_all(0:numPEs-1) )
224           IF(.NOT.ALLOCATED(kstart3_all))    allocate( kstart3_all(0:numPEs-1) )
225     
226           IF(.NOT.ALLOCATED(istart4_all))    allocate( istart4_all(0:numPEs-1) )
227           IF(.NOT.ALLOCATED(jstart4_all))    allocate( jstart4_all(0:numPEs-1) )
228           IF(.NOT.ALLOCATED(kstart4_all))    allocate( kstart4_all(0:numPEs-1) )
229     
230           IF(.NOT.ALLOCATED(iend_all))       allocate( iend_all(0:numPEs-1) )
231           IF(.NOT.ALLOCATED(jend_all))       allocate( jend_all(0:numPEs-1) )
232           IF(.NOT.ALLOCATED(kend_all))       allocate( kend_all(0:numPEs-1) )
233     
234           IF(.NOT.ALLOCATED(iend1_all))      allocate( iend1_all(0:numPEs-1) )
235           IF(.NOT.ALLOCATED(jend1_all))      allocate( jend1_all(0:numPEs-1) )
236           IF(.NOT.ALLOCATED(kend1_all))      allocate( kend1_all(0:numPEs-1) )
237     
238           IF(.NOT.ALLOCATED(iend2_all))      allocate( iend2_all(0:numPEs-1) )
239           IF(.NOT.ALLOCATED(jend2_all))      allocate( jend2_all(0:numPEs-1) )
240           IF(.NOT.ALLOCATED(kend2_all))      allocate( kend2_all(0:numPEs-1) )
241     
242           IF(.NOT.ALLOCATED(iend3_all))      allocate( iend3_all(0:numPEs-1) )
243           IF(.NOT.ALLOCATED(jend3_all))      allocate( jend3_all(0:numPEs-1) )
244           IF(.NOT.ALLOCATED(kend3_all))      allocate( kend3_all(0:numPEs-1) )
245     
246           IF(.NOT.ALLOCATED(iend4_all))      allocate( iend4_all(0:numPEs-1) )
247           IF(.NOT.ALLOCATED(jend4_all))      allocate( jend4_all(0:numPEs-1) )
248           IF(.NOT.ALLOCATED(kend4_all))      allocate( kend4_all(0:numPEs-1) )
249     
250           IF(.NOT.ALLOCATED(displs))         allocate( displs(0:numPEs-1) )
251     
252     
253           CALL PARTITION(CYC_XL, CYC_YL, CYC_ZL)
254     
255     ! The upper and lower bounds are prescribed such that two ghost
256     ! layers are allowed at the physical boundaries - this is consistent
257     ! with our present approach - need to be generalized if only one ghost
258     ! layer is needed
259           do iproc=0,numPEs-1
260              istart2_all(iproc) = max(imin1-1,min(imax1+1,istart1_all(iproc)-1))
261              if(nodesi.ne.1) then
262                 istart3_all(iproc) = max(imin1-2,min(imax1+2,istart2_all(iproc)-1))
263                 istart4_all(iproc) = max(imin1-3,min(imax1+3,istart3_all(iproc)-1))
264              else
265                 istart3_all(iproc) = istart2_all(iproc)
266                 istart4_all(iproc) = istart3_all(iproc)
267              endif
268     
269              jstart2_all(iproc) = max(jmin1-1,min(jmax1+1,jstart1_all(iproc)-1))
270              if(nodesj.ne.1) then
271                 jstart3_all(iproc) = max(jmin1-2,min(jmax1+2,jstart2_all(iproc)-1))
272                 jstart4_all(iproc) = max(jmin1-3,min(jmax1+3,jstart3_all(iproc)-1))
273              else
274                 jstart3_all(iproc) = jstart2_all(iproc)
275                 jstart4_all(iproc) = jstart3_all(iproc)
276              endif
277     
278              if(no_k) then
279                 kstart2_all(iproc) = kstart1_all(iproc)
280                 kstart3_all(iproc) = kstart1_all(iproc)
281                 kstart4_all(iproc) = kstart1_all(iproc)
282              else
283                 kstart2_all(iproc) = max(kmin1-1,min(kmax1+1,kstart1_all(iproc)-1))
284                 if(nodesk.ne.1) then
285                    kstart3_all(iproc) = max(kmin1-2,min(kmax1+2,kstart2_all(iproc)-1))
286                    kstart4_all(iproc) = max(kmin1-3,min(kmax1+3,kstart3_all(iproc)-1))
287                 else
288                    kstart3_all(iproc) =  kstart2_all(iproc)
289                    kstart4_all(iproc) =  kstart3_all(iproc)
290                 endif
291              endif
292     
293              iend2_all(iproc) = max(imin1-1,min(imax1+1,iend1_all(iproc)+1))
294              if(nodesi.ne.1) then
295                 iend3_all(iproc) = max(imin1-2,min(imax1+2,iend2_all(iproc)+1))
296                 iend4_all(iproc) = max(imin1-3,min(imax1+3,iend3_all(iproc)+1))
297              else
298                 iend3_all(iproc) = iend2_all(iproc)
299                 iend4_all(iproc) = iend3_all(iproc)
300              endif
301     
302              jend2_all(iproc) = max(jmin1-1,min(jmax1+1,jend1_all(iproc)+1))
303              if(nodesj.ne.1) then
304                 jend3_all(iproc) = max(jmin1-2,min(jmax1+2,jend2_all(iproc)+1))
305                 jend4_all(iproc) = max(jmin1-3,min(jmax1+3,jend3_all(iproc)+1))
306              else
307                 jend3_all(iproc) = jend2_all(iproc)
308                 jend4_all(iproc) = jend3_all(iproc)
309              endif
310     
311              if(no_k) then
312                 kend2_all(iproc) = kend1_all(iproc)
313                 kend3_all(iproc) = kend1_all(iproc)
314                 kend4_all(iproc) = kend1_all(iproc)
315              else
316                 kend2_all(iproc) = max(kmin1-1,min(kmax1+1,kend1_all(iproc)+1))
317                 if(nodesk.ne.1) then
318                    kend3_all(iproc) = max(kmin1-2,min(kmax1+2,kend2_all(iproc)+1))
319                    kend4_all(iproc) = max(kmin1-3,min(kmax1+3,kend3_all(iproc)+1))
320                 else
321                    kend3_all(iproc) = kend2_all(iproc)
322                    kend4_all(iproc) = kend3_all(iproc)
323                 endif
324              endif
325           enddo
326     
327     ! for higher order methods
328           if(.not.fpfoi) then
329              do iproc=0,numPEs-1
330                 istart4_all(iproc) = istart3_all(iproc)
331                 jstart4_all(iproc) = jstart3_all(iproc)
332                 kstart4_all(iproc) = kstart3_all(iproc)
333                 iend4_all(iproc)   = iend3_all(iproc)
334                 jend4_all(iproc)   = jend3_all(iproc)
335                 kend4_all(iproc)   = kend3_all(iproc)
336              enddo
337           endif
338     
339           do iproc=0,numPEs-1
340              ijkstart3_all(iproc) = 1
341              ijkend3_all(iproc) =  1 + (iend3_all(iproc) - istart3_all(iproc)) &
342                + (jend3_all(iproc)-jstart3_all(iproc))*(iend3_all(iproc)-istart3_all(iproc)+1) &
343                + (kend3_all(iproc)-kstart3_all(iproc))*(jend3_all(iproc)-jstart3_all(iproc)+1)* &
344                  (iend3_all(iproc)-istart3_all(iproc)+1)
345     
346              ijkstart4_all(iproc) = 1
347              ijkend4_all(iproc) =  1 + (iend4_all(iproc) - istart4_all(iproc)) &
348                 + (jend4_all(iproc)-jstart4_all(iproc))*(iend4_all(iproc)-istart4_all(iproc)+1) &
349                 + (kend4_all(iproc)-kstart4_all(iproc))*(jend4_all(iproc)-jstart4_all(iproc)+1)* &
350                   (iend4_all(iproc)-istart4_all(iproc)+1)
351           enddo
352     
353           do iproc=0,numPEs-1
354              ijksize3_all(iproc) = ijkend3_all(iproc) - ijkstart3_all(iproc) + 1
355              ijksize4_all(iproc) = ijkend4_all(iproc) - ijkstart4_all(iproc) + 1
356           enddo
357     
358           displs(0) = 0
359           do iproc=1,numPEs-1
360              displs(iproc) = displs(iproc-1)+ijksize3_all(iproc-1)
361     !       write(*,*) 'displ',displs(iproc),iproc, ijksize3_all(iproc)
362           enddo
363     
364     
365           ijkstart3 = ijkstart3_all(myPE)
366           ijkend3   = ijkend3_all(myPE)
367           ijksize3  = ijksize3_all(myPE)
368     
369           ijkstart4 = ijkstart4_all(myPE)
370           ijkend4   = ijkend4_all(myPE)
371           ijksize4  = ijksize4_all(myPE)
372     
373           istart1   = istart1_all(myPE)
374           iend1     = iend1_all(myPE)
375           jstart1   = jstart1_all(myPE)
376           jend1     = jend1_all(myPE)
377           kstart1   = kstart1_all(myPE)
378           kend1     = kend1_all(myPE)
379     
380           istart2   = istart2_all(myPE)
381           iend2     = iend2_all(myPE)
382           jstart2   = jstart2_all(myPE)
383           jend2     = jend2_all(myPE)
384           kstart2   = kstart2_all(myPE)
385           kend2     = kend2_all(myPE)
386     
387           istart3   = istart3_all(myPE)
388           iend3     = iend3_all(myPE)
389           jstart3   = jstart3_all(myPE)
390           jend3     = jend3_all(myPE)
391           kstart3   = kstart3_all(myPE)
392           kend3     = kend3_all(myPE)
393     
394           istart4   = istart4_all(myPE)
395           iend4     = iend4_all(myPE)
396           jstart4   = jstart4_all(myPE)
397           jend4     = jend4_all(myPE)
398           kstart4   = kstart4_all(myPE)
399           kend4     = kend4_all(myPE)
400     
401           IF(.not.allocated(NCPP_UNIFORM)) allocate( NCPP_UNIFORM(0:NumPEs-1))
402     
403           IF(.NOT.NCPP_UNIFORM_BACKED_UP) THEN
404              NCPP_UNIFORM(MyPE) = ijksize3_all(MyPE)
405           ENDIF
406           NCPP_UNIFORM_BACKED_UP = .TRUE.
407     
408           IF(SHORT_GRIDMAP_INIT) THEN
409     !        do iproc=0,numPEs-1
410     !           NCPP_UNIFORM(iproc) = ijksize3_all(iproc)
411     !        enddo
412              RETURN
413           ENDIF
414     
415     ! Setup mapping to take care of cyclic boundary conditions
416     ! ---------------------------------------------------------------->>>
417     ! consider cyclic boundary condition using the imap(:),jmap(:),kmap(:)
418     ! indirection arrays
419     
420           allocate( imap( imin4:imax4 ) )
421           allocate( jmap( jmin4:jmax4 ) )
422           allocate( kmap( kmin4:kmax4 ) )
423     
424           allocate( imap_c( imin4:imax4 ) )
425           allocate( jmap_c( jmin4:jmax4 ) )
426           allocate( kmap_c( kmin4:kmax4 ) )
427     
428           do kk=kmin4,kmax4
429             kmap(kk) = kk
430           enddo
431     
432           do jj=jmin4,jmax4
433             jmap(jj) = jj
434           enddo
435     
436           do ii=imin4,imax4
437             imap(ii) = ii
438           enddo
439     
440           if (CYC_ZL) then
441              kmap( kmax2 ) = kmin1
442              kmap( kmin2 ) = kmax1
443              if (kmax3.gt.kmax2) kmap(kmax3) = kmap(kmax2)+1
444              if (kmin3.lt.kmin2) kmap(kmin3) = kmap(kmin2)-1
445              if (kmax4.gt.kmax3) kmap(kmax4) = kmap(kmax3)+1
446              if (kmin4.lt.kmin3) kmap(kmin4) = kmap(kmin3)-1
447           endif
448     
449           if (CYC_YL) then
450              jmap( jmax2 ) = jmin1
451              jmap( jmin2 ) = jmax1
452              if (jmax3.gt.jmax2) jmap(jmax3) = jmap(jmax2)+1
453              if (jmin3.lt.jmin2) jmap(jmin3) = jmap(jmin2)-1
454              if (jmax4.gt.jmax3) jmap(jmax4) = jmap(jmax3)+1
455              if (jmin4.lt.jmin3) jmap(jmin4) = jmap(jmin3)-1
456           endif
457     
458           if (CYC_XL) then
459              imap( imax2 ) = imin1
460              imap( imin2 ) = imax1
461              if (imax3.gt.imax2) imap(imax3) = imap(imax2)+1
462              if (imin3.lt.imin2) imap(imin3) = imap(imin2)-1
463              if (imax4.gt.imax3) imap(imax4) = imap(imax3)+1
464              if (imin4.lt.imin3) imap(imin4) = imap(imin3)-1
465           endif
466     
467           do kk=kmin4,kmax4
468             kmap_c(kk) = kk
469           enddo
470     
471           do jj=jmin4,jmax4
472             jmap_c(jj) = jj
473           enddo
474     
475           do ii=imin4,imax4
476             imap_c(ii) = ii
477           enddo
478     
479           if (CYC_ZL.and.nodesk.eq.1) then
480              kmap_c( kmax2 ) = kmin1
481              kmap_c( kmin2 ) = kmax1
482              if (kmax3.gt.kmax2) kmap_c(kmax3) = kmap_c(kmax2)+1
483              if (kmin3.lt.kmin2) kmap_c(kmin3) = kmap_c(kmin2)-1
484              if (kmax4.gt.kmax3) kmap_c(kmax4) = kmap_c(kmax3)+1
485              if (kmin4.lt.kmin3) kmap_c(kmin4) = kmap_c(kmin3)-1
486           endif
487     
488           if (CYC_YL.and.nodesj.eq.1) then
489              jmap_c( jmax2 ) = jmin1
490              jmap_c( jmin2 ) = jmax1
491              if (jmax3.gt.jmax2) jmap_c(jmax3) = jmap_c(jmax2)+1
492              if (jmin3.lt.jmin2) jmap_c(jmin3) = jmap_c(jmin2)-1
493              if (jmax4.gt.jmax3) jmap_c(jmax4) = jmap_c(jmax3)+1
494              if (jmin4.lt.jmin3) jmap_c(jmin4) = jmap_c(jmin3)-1
495           endif
496     
497           if (CYC_XL.and.nodesi.eq.1) then
498              imap_c( imax2 ) = imin1
499              imap_c( imin2 ) = imax1
500              if (imax3.gt.imax2) imap_c(imax3) = imap_c(imax2)+1
501              if (imin3.lt.imin2) imap_c(imin3) = imap_c(imin2)-1
502              if (imax4.gt.imax3) imap_c(imax4) = imap_c(imax3)+1
503              if (imin4.lt.imin3) imap_c(imin4) = imap_c(imin3)-1
504           endif
505     ! End setup mapping to take care of cyclic boundary conditions
506     ! ----------------------------------------------------------------<<<
507     
508     
509     ! Defining new set of varaibles to define upper and lower bound of the
510     ! indices to include actual physical boundaries of the problem
511           do iproc = 0, numPEs-1
512              istart = istart1
513              iend = iend1
514              jstart = jstart1
515              jend = jend1
516              kstart = kstart1
517              kend = kend1
518     
519              if(istart3.eq.imin3) istart = istart2
520              if(iend3.eq.imax3) iend = iend2
521              if(jstart3.eq.jmin3) jstart = jstart2
522              if(jend3.eq.jmax3) jend = jend2
523              if(kstart3.eq.kmin3) kstart = kstart2
524              if(kend3.eq.kmax3) kend = kend2
525     
526              istart_all(iproc) = istart
527              iend_all(iproc)   = iend
528              jstart_all(iproc) = jstart
529              jend_all(iproc)   = jend
530              kstart_all(iproc) = kstart
531              kend_all(iproc)   = kend
532           enddo
533     
534           IF(numPEs .GT. 1) THEN
535     ! Calculate any load imbalance.
536              IMBALANCE = DBLE(maxval(ijksize3_all)- minval(ijksize3_all)) /    &
537                 minval(ijksize3_all)*100.0
538     ! Calculate potential speedup based on Amdahl's Law
539              IF(IMBALANCE == 0)THEN
540                 AMDAHL_SPEEDUP='+Inf'
541              ELSE
542                 AMDAHL_SPEEDUP=''
543                 WRITE(AMDAHL_SPEEDUP,*)1.0/dble(IMBALANCE)
544              ENDIF
545     ! Construct a message for the user telling them the grid partition info.
546              WRITE(ERR_MSG,1000)maxval(ijksize3_all), maxloc(ijksize3_all),&
547                 minval(ijksize3_all), minloc(ijksize3_all),                &
548                 sum(ijksize3_all)/numPEs, trim(AMDAHL_SPEEDUP)
549              CALL FLUSH_ERR_MSG
550           ENDIF
551     
552     ! Setup coefficients of FUINIJK
553             c0 = 1 - jstart3_all(myPE)
554             c1 = (jend3_all(myPE)-jstart3_all(myPE)+1)
555             c2 = (jend3_all(myPE)-jstart3_all(myPE)+1)* (iend3_all(myPE)-istart3_all(myPE)+1)
556             c0 =  c0  - c1*istart3_all(myPE) - c2*kstart3_all(myPE)
557     
558     ! Setup coefficients of FUINIJK3
559             c0_3 = 1 - jstart4_all(myPE)
560             c1_3 = (jend4_all(myPE)-jstart4_all(myPE)+1)
561             c2_3 = (jend4_all(myPE)-jstart4_all(myPE)+1)* (iend4_all(myPE)-istart4_all(myPE)+1)
562             c0_3 =  c0_3  - c1_3*istart4_all(myPE) - c2_3*kstart4_all(myPE)
563     
564     !   Initialize Array mapping (I,J,K) to IJK
565             INCREMENT_ARRAYS_ALLOCATED = .FALSE.
566     
567     ! These arrays could already be allocated in post_mfix
568     ! when interpolating old data to new grid
569     
570             if(allocated(IJK_ARRAY_OF)) deallocate(IJK_ARRAY_OF)
571             if(allocated(FUNIJK_MAP_C)) deallocate(FUNIJK_MAP_C)
572             if(allocated(DEAD_CELL_AT)) deallocate(DEAD_CELL_AT)
573     
574             ! Must extend range such that neighbors (IM,JP etc...) stay in bound
575             allocate(IJK_ARRAY_OF(istart3-1:iend3+1,jstart3-1:jend3+1,kstart3-1:kend3+1))
576             allocate(FUNIJK_MAP_C(istart3-1:iend3+1,jstart3-1:jend3+1,kstart3-1:kend3+1))
577             allocate(DEAD_CELL_AT(imin3-1:imax3+1,jmin3-1:jmax3+1,kmin3-1:kmax3+1))
578     
579             DEAD_CELL_AT = .FALSE.
580     
581     ! Save IJK value of (I,J,K) cell in an array
582     ! IJK_ARRAY_OF(I,J,K) will replace the use of FUNIJK(I,J,K)
583             DO ii = istart3,iend3
584                DO jj = jstart3,jend3
585                   DO kk = kstart3,kend3
586                      IJK_ARRAY_OF(ii,jj,kk)=FUNIJK_0(ii,jj,kk)
587                   ENDDO
588                ENDDO
589             ENDDO
590     
591             DO ii = istart3,iend3
592                DO jj = jstart3,jend3
593                   DO kk = kstart3,kend3
594                      FUNIJK_MAP_C(ii,jj,kk)=IJK_ARRAY_OF(IMAP_C(ii),JMAP_C(jj),KMAP_C(kk))
595                   ENDDO
596                ENDDO
597             ENDDO
598     
599     
600     ! Call to sendrecv_init to set all the communication pattern
601           COMM = MPI_COMM_WORLD
602           CALL SENDRECV_INIT(COMM, CYC_XL, CYC_YL, CYC_ZL, idebug=0)
603     
604           CALL FINL_ERR_MSG
605           RETURN
606     
607      1000 FORMAT('Parallel load balancing statistics:',2/,13x,'Comp. cells',&
608           4X,'Processor',/3X,'maximum   ',I11,4X,I9,/3X,'minimum   ',I11,&
609           4X,I9,/3X,'average   ',I11,6X,'-N/A-',2/,3X,'Maximum speedup ',&
610           '(Amdahls Law) = ',A)
611     
612           END SUBROUTINE GRIDMAP_INIT
613     
614           END MODULE GRIDMAP
615