File: N:\mfix\model\des\desgrid_mod.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Module: desgrid                                                     !
4     !  Author: Pradeep G                                                   !
5     !                                                                      !
6     !  Purpose: Defines the desgrid and sets the indices; also sets the    !
7     !  communication between the desgrid for ghost cell exchange.          !
8     !                                                                      !
9     !  Comments: The variables are named similar to the fluid grid. More   !
10     !  descriptions about naming can be found in compar_mod under          !
11     !  dmp_modules folder.                                                 !
12     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
13           MODULE DESGRID
14     
15     ! variables related to global block structure
16     !---------------------------------------------------------------------//
17           integer  :: dg_imin1, dg_imax1, &
18                       dg_imin2, dg_imax2
19     
20           integer  :: dg_jmin1, dg_jmax1, &
21                       dg_jmin2, dg_jmax2
22     
23           integer  :: dg_kmin1, dg_kmax1, &
24                       dg_kmin2, dg_kmax2
25     
26     ! variables contain information of local indices of all processors
27     !---------------------------------------------------------------------//
28           integer,dimension (:),allocatable ::      &
29                       dg_istart1_all, dg_iend1_all, &
30                       dg_istart2_all, dg_iend2_all, &
31                       dg_isize_all
32     
33           integer,dimension (:),allocatable ::      &
34                       dg_jstart1_all, dg_jend1_all, &
35                       dg_jstart2_all, dg_jend2_all, &
36                       dg_jsize_all
37     
38           integer,dimension (:),allocatable ::      &
39                       dg_kstart1_all, dg_kend1_all, &
40                       dg_kstart2_all, dg_kend2_all, &
41                       dg_ksize_all
42     
43     ! variables relate to processor specific details
44     !---------------------------------------------------------------------//
45           integer  :: dg_istart,  dg_iend,  &
46                       dg_istart1, dg_iend1, &
47                       dg_istart2, dg_iend2
48     
49           integer  :: dg_jstart,  dg_jend,  &
50                       dg_jstart1, dg_jend1, &
51                       dg_jstart2, dg_jend2
52     
53           integer  :: dg_kstart,  dg_kend,  &
54                       dg_kstart1, dg_kend1, &
55                       dg_kstart2, dg_kend2
56     
57           integer  :: dg_ijkstart2, dg_ijkend2, dg_ijksize2
58     
59     
60     ! variables related to processor domain size
61     !---------------------------------------------------------------------//
62           double precision :: dg_xstart, dg_xend, dg_dxinv
63           double precision :: dg_ystart, dg_yend, dg_dyinv
64           double precision :: dg_zstart, dg_zend, dg_dzinv
65     
66     ! Variables related to cyclic boundary  used in desgrid_functions
67           integer,dimension(:,:),allocatable :: dg_cycoffset, icycoffset
68     
69           integer :: dg_pic_max_init = 25
70     
71     ! constants required for functions computing local and global ijkvalues
72           integer dg_c1_gl, dg_c2_gl, dg_c3_gl  ! global
73           integer dg_c1_lo, dg_c2_lo, dg_c3_lo  ! local
74     
75           integer, dimension(:), allocatable :: dg_c1_all
76           integer, dimension(:), allocatable :: dg_c2_all
77           integer, dimension(:), allocatable :: dg_c3_all
78     
79           contains
80     
81     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
82     !  Function: PROCIJJK                                                  !
83     !                                                                      !
84     !  Purpose: Calculate the process that owns the cell with the given    !
85     !  I, J, K indicies.                                                   !
86     !......................................................................!
87           integer function procijk(fi,fj,fk)
88           use compar, only: nodesi, nodesj
89           implicit none
90           integer fi,fj,fk
91           procijk =fi+fj*nodesi+fk*nodesi*nodesj
92           end function procijk
93     
94     
95     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
96     !  Function: IofPROC                                                   !
97     !                                                                      !
98     !  Purpose: Return the I index of the current rank given an IJK value. !
99     !......................................................................!
100           integer function iofproc(fijk)
101           use compar, only: nodesi, nodesj
102           implicit none
103           integer fijk
104           iofproc = mod(mod(fijk,nodesi*nodesj),nodesi)
105           end function iofproc
106     
107     
108     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
109     !  Function: JofPROC                                                   !
110     !                                                                      !
111     !  Purpose: Return the J index of the current rank given an IJK value. !
112     !......................................................................!
113           integer function jofproc(fijk)
114           use compar, only: nodesi, nodesj
115           implicit none
116           integer fijk
117           jofproc = mod(fijk - iofproc(fijk),nodesi*nodesj)/nodesi
118           end function jofproc
119     
120     
121     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
122     !  Function: KofPROC                                                   !
123     !                                                                      !
124     !  Purpose: Return the K index of the current rank given an IJK value. !
125     !......................................................................!
126           integer function kofproc(fijk)
127           use compar, only: nodesi, nodesj
128           implicit none
129           integer fijk
130           kofproc = (fijk - iofproc(fijk) - jofproc(fijk)*nodesi) / &
131              (nodesi*nodesj)
132           end function kofproc
133     
134     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
135     !  Function: DG_FUNIJK                                                 !
136     !                                                                      !
137     !  Purpose: Calculate the desgrid IJK given I, J, K components. The    !
138     !  result is local to the rank calculating the IJK.                    !
139     !......................................................................!
140           integer function dg_funijk(fi,fj,fk)
141           implicit none
142           integer fi,fj,fk
143           dg_funijk = fj+fi*dg_c1_lo+fk*dg_c2_lo+dg_c3_lo
144           end function dg_funijk
145     
146     
147     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
148     !  Function: DG_FUNIJK_GL                                              !
149     !                                                                      !
150     !  Purpose: Calculate the global desgrid IJK given I, J, K components. !
151     !  result should be consistent across all ranks.                       !
152     !......................................................................!
153           integer function dg_funijk_gl(fi,fj,fk)
154           implicit none
155           integer fi,fj,fk
156           dg_funijk_gl = fj+fi*dg_c1_gl+fk*dg_c2_gl+dg_c3_gl
157           end function dg_funijk_gl
158     
159     
160     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
161     !  Function: DG_FUNIJK_PROC                                            !
162     !                                                                      !
163     !  Purpose: Calculate the local desgrid IJK given I, J, K components   !
164     !  for the specified processor.                                        !
165     !......................................................................!
166           integer function dg_funijk_proc(fi,fj,fk,fproc)
167           implicit none
168           integer fi,fj,fk,fproc
169           dg_funijk_proc = fj + fi*dg_c1_all(fproc) + &
170              fk*dg_c2_all(fproc) + dg_c3_all(fproc)
171           end function dg_funijk_proc
172     
173     
174     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
175     !  Function: DG_FUNIM                                                  !
176     !                                                                      !
177     !  Purpose: Return the local IJK value for the cell to the west.       !
178     !  This is equivalent to calculating: dg_funijk(I-1,J,K)               !
179     !......................................................................!
180           integer function dg_funim(fijk)
181           implicit none
182           integer fijk
183           dg_funim = fijk - dg_c1_lo
184           end function dg_funim
185     
186     
187     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
188     !  Function: DG_FUNIP                                                  !
189     !                                                                      !
190     !  Purpose: Return the local IJK value for the cell to the east.       !
191     !  This is equivalent to calculating: dg_funijk(I+1,J,K)               !
192     !......................................................................!
193           integer function dg_funip(fijk)
194           implicit none
195           integer fijk
196           dg_funip = fijk + dg_c1_lo
197           end function dg_funip
198     
199     
200     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
201     !  Function: DG_FUNJM                                                  !
202     !                                                                      !
203     !  Purpose: Return the local IJK value for the cell to the south.      !
204     !  This is equivalent to calculating: dg_funijk(I,J-1,K)               !
205     !......................................................................!
206           integer function dg_funjm(fijk)
207           implicit none
208           integer fijk
209           dg_funjm = fijk - 1
210           end function dg_funjm
211     
212     
213     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
214     !  Function: DG_FUNJP                                                  !
215     !                                                                      !
216     !  Purpose: Return the local IJK value for the cell to the north.      !
217     !  This is equivalent to calculating: dg_funijk(I,J+1,K)               !
218     !......................................................................!
219           integer function dg_funjp(fijk)
220           implicit none
221           integer fijk
222           dg_funjp = fijk + 1
223           end function dg_funjp
224     
225     
226     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
227     !  Function: DG_FUNKM                                                  !
228     !                                                                      !
229     !  Purpose: Return the local IJK value for the cell to the bottom.     !
230     !  This is equivalent to calculating: dg_funijk(I,J,K-1)               !
231     !......................................................................!
232           integer function dg_funkm(fijk)
233           implicit none
234           integer fijk
235           dg_funkm = fijk - dg_c2_lo
236           end function dg_funkm
237     
238     
239     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
240     !  Function: DG_FUNJM                                                  !
241     !                                                                      !
242     !  Purpose: Return the local IJK value for the cell to the top.        !
243     !  This is equivalent to calculating: dg_funijk(I,J,K+1)               !
244     !......................................................................!
245           integer function dg_funkp(fijk)
246           implicit none
247           integer fijk
248           dg_funkp = fijk + dg_c2_lo
249           end function dg_funkp
250     
251     
252     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
253     !  Function: DG_JOF_GL                                                 !
254     !                                                                      !
255     !  Purpose: Calculate the J index from a global IJK value.             !
256     !......................................................................!
257           integer function dg_jof_gl(fijk)
258           implicit none
259           integer fijk
260           dg_jof_gl = mod(mod(fijk-1,dg_c2_gl),dg_c1_gl)+dg_jmin2
261           end function dg_jof_gl
262     
263     
264     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
265     !  Function: DG_IOF_GL                                                 !
266     !                                                                      !
267     !  Purpose: Calculate the I index from a global IJK value.             !
268     !......................................................................!
269           integer function dg_iof_gl(fijk)
270           implicit none
271           integer fijk
272           dg_iof_gl = (mod(fijk-dg_jof_gl(fijk)+dg_jmin2-1,dg_c2_gl)) /    &
273              dg_c1_gl + dg_imin2
274           end function dg_iof_gl
275     
276     
277     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
278     !  Function: DG_KOF_GL                                                 !
279     !                                                                      !
280     !  Purpose: Calculate the K index from a global IJK value.             !
281     !......................................................................!
282           integer function dg_kof_gl(fijk)
283           implicit none
284           integer fijk
285           dg_kof_gl = (fijk-dg_c3_gl-dg_iof_gl(fijk)*                      &
286              dg_c1_gl-dg_jof_gl(fijk))/dg_c2_gl
287           end function dg_kof_gl
288     
289     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
290     !  Function: DG_JOF_LO                                                 !
291     !                                                                      !
292     !  Purpose: Calculate the J index from a local IJK value.              !
293     !......................................................................!
294           integer function dg_jof_lo(fijk)
295           implicit none
296           integer fijk
297           dg_jof_lo = mod(mod(fijk-1,dg_c2_lo),dg_c1_lo)+dg_jstart2
298           end function dg_jof_lo
299     
300     
301     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
302     !  Function: DG_IOF_LO                                                 !
303     !                                                                      !
304     !  Purpose: Calculate the I index from a local IJK value.              !
305     !......................................................................!
306           integer function dg_iof_lo(fijk)
307           implicit none
308           integer fijk
309           dg_iof_lo = (mod(fijk-dg_jof_lo(fijk)+dg_jstart2-1,dg_c2_lo)) /  &
310              dg_c1_lo + dg_istart2
311           end function dg_iof_lo
312     
313     
314     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
315     !  Function: DG_KOF_LO                                                 !
316     !                                                                      !
317     !  Purpose: Calculate the K index from a local IJK value.              !
318     !......................................................................!
319           integer function dg_kof_lo(fijk)
320           implicit none
321           integer fijk
322           dg_kof_lo = (fijk-dg_c3_lo-dg_iof_lo(fijk)*                      &
323              dg_c1_lo-dg_jof_lo(fijk))/dg_c2_lo
324           end function dg_kof_lo
325     
326     
327     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
328     !  Function: DG_IJKCONV                                                !
329     !                                                                      !
330     !  Purpose: Convert the IJK from once process to another.              !
331     !                                                                      !
332     !......................................................................!
333           integer function dg_ijkconv(fijk,fface,fto_proc)
334           implicit none
335           integer fijk,fto_proc,fface
336           dg_ijkconv = dg_funijk_proc(dg_iof_lo(fijk)+                    &
337             dg_cycoffset(fface,1), dg_jof_lo(fijk)+dg_cycoffset(fface,2), &
338             dg_kof_lo(fijk)+dg_cycoffset(fface,3), fto_proc)
339           end function dg_ijkconv
340     
341     
342     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
343     !  Function: IofPOS                                                    !
344     !                                                                      !
345     !  Purpose: Calculate the cell I index containing the given position.  !
346     !......................................................................!
347           integer function iofpos(fpos)
348           implicit none
349           double precision fpos
350           iofpos = floor((fpos-dg_xstart)*dg_dxinv) + dg_istart1
351           end function iofpos
352     
353     
354     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
355     !  Function: JofPOS                                                    !
356     !                                                                      !
357     !  Purpose: Calculate the cell J index containing the given position.  !
358     !......................................................................!
359           integer function jofpos(fpos)
360           implicit none
361           double precision fpos
362           jofpos = floor((fpos-dg_ystart)*dg_dyinv) + dg_jstart1
363           end function jofpos
364     
365     
366     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
367     !  Function: KofPOS                                                    !
368     !                                                                      !
369     !  Purpose: Calculate the cell K index containing the given position.  !
370     !......................................................................!
371           integer function kofpos(fpos)
372           implicit none
373           double precision fpos
374           kofpos = floor((fpos-dg_zstart)*dg_dzinv) + dg_kstart1
375           end function kofpos
376     
377     
378     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
379     !  Function: dg_is_ON_myPE_OWNS                                        !
380     !                                                                      !
381     !  Purpose: Determine if the current rank owns this cell based on the  !
382     !  I, J, K values.                                                     !
383     !......................................................................!
384           logical function dg_is_ON_myPE_OWNS(lI, lJ, lK)
385           implicit none
386           integer, intent(in) :: lI, lJ, lK
387     
388           dg_is_ON_myPE_OWNS = (&
389              (dg_istart <= lI) .AND. (lI <= dg_iend) .AND. &
390              (dg_jstart <= lJ) .AND. (lJ <= dg_jend) .AND. &
391              (dg_kstart <= lK) .AND. (lK <= dg_kend))
392     
393           end function
394     
395     
396     !''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
397     !  Function: dg_is_ON_myPE_plus1layers                                 !
398     !                                                                      !
399     !  Purpose: Determine if the current rank contains the current cell    !
400     !  as its own or in its single layer of desgrid ghost cells.           !
401     !......................................................................!
402           logical function dg_is_ON_myPE_plus1layers(lI, lJ, lK)
403           implicit none
404           integer, intent(in) :: lI, lJ, lK
405     
406           dg_is_ON_myPE_plus1layers = (&
407              (dg_istart2 <= lI) .AND. (lI <= dg_iend2) .AND. &
408              (dg_jstart2 <= lJ) .AND. (lJ <= dg_jend2) .AND. &
409              (dg_kstart2 <= lK) .AND. (lK <= dg_kend2))
410     
411           end function
412     
413     
414     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
415     !                                                                      !
416     !  Subroutine: DESGRID_INIT                                            !
417     !  Author: Pradeep G                                                   !
418     !                                                                      !
419     !  Purpose: Sets indices for the DESGRID and defines constants         !
420     !  required for the DESGRID functions. This is needed for              !
421     !  communication between the DESGRID for ghost cell exchanges.         !
422     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
423           subroutine desgrid_init()
424     
425           use funits
426           use compar
427           use discretelement, only: DES_PERIODIC_WALLS_X
428           use discretelement, only: DES_PERIODIC_WALLS_Y
429           use discretelement, only: DES_PERIODIC_WALLS_Z
430           use discretelement, only: DIMN
431           use discretelement, only: XE, YN, ZT
432           use constant
433           use desmpi_wrapper
434     
435     ! Global Variables:
436     !---------------------------------------------------------------------//
437     ! Domain partition for DEM background mesh.
438           use discretelement, only: DESGRIDSEARCH_IMAX
439           use discretelement, only: DESGRIDSEARCH_JMAX
440           use discretelement, only: DESGRIDSEARCH_KMAX
441     ! Domain size specifed by the user.
442           use geometry, only: XLENGTH, YLENGTH, ZLENGTH, NO_K
443           use derived_types, only: dg_pic
444     
445     ! Use the error manager for posting error messages.
446     !---------------------------------------------------------------------//
447           use error_manager
448     
449           implicit none
450     !-----------------------------------------------
451     ! Local varibles
452     !-----------------------------------------------
453           double precision :: ltempdx,ltempdy,ltempdz
454           integer :: lijkproc,liproc,ljproc,lkproc
455           integer :: lijk
456     !......................................................................!
457     
458     ! Initialize the error manager.
459           CALL INIT_ERR_MSG("DESGRID_INIT")
460     
461     ! set indices for all processors
462          allocate (dg_istart1_all(0:numpes-1), dg_iend1_all(0:numpes-1))
463          allocate (dg_istart2_all(0:numpes-1), dg_iend2_all(0:numpes-1))
464          allocate (dg_isize_all(0:nodesi-1))
465     
466          allocate (dg_jstart1_all(0:numpes-1), dg_jend1_all(0:numpes-1))
467          allocate (dg_jstart2_all(0:numpes-1), dg_jend2_all(0:numpes-1))
468          allocate (dg_jsize_all(0:nodesj-1))
469     
470          allocate (dg_kstart1_all(0:numpes-1), dg_kend1_all(0:numpes-1))
471          allocate (dg_kstart2_all(0:numpes-1), dg_kend2_all(0:numpes-1))
472          allocate (dg_ksize_all(0:nodesk-1))
473     
474     
475           dg_istart1_all=0; dg_iend1_all=0
476           dg_istart2_all=0; dg_iend2_all=0
477           dg_isize_all=0
478     
479           dg_jstart1_all=0; dg_jend1_all=0
480           dg_jstart2_all=0; dg_jend2_all=0
481           dg_jsize_all=0
482     
483           dg_kstart1_all=0; dg_kend1_all=0
484           dg_kstart2_all=0; dg_kend2_all=0
485           dg_ksize_all=0
486     
487     
488     ! set grid size based on user input desgridsearch_<ijk>max
489           ltempdx = xlength/desgridsearch_imax
490           ltempdy = ylength/desgridsearch_jmax
491           if(do_k) ltempdz = zlength/desgridsearch_kmax
492           dg_ksize_all(:) = 1
493     
494           lijkproc = 0
495           do lkproc=0,nodesk-1
496              do ljproc=0,nodesj-1
497                 do liproc=0,nodesi-1
498                    dg_isize_all(liproc) = NINT((xe(iend1_all(lijkproc))-xe(istart1_all(lijkproc)-1))/ltempdx)
499                    dg_jsize_all(ljproc) = NINT((yn(jend1_all(lijkproc))-yn(jstart1_all(lijkproc)-1))/ltempdy)
500                    if(do_k) dg_ksize_all(lkproc) = NINT((zt(kend1_all(lijkproc))-zt(kstart1_all(lijkproc)-1))/ltempdz)
501                    dg_istart1_all(lijkproc) = sum(dg_isize_all(0:liproc-1)) + 2
502                    dg_jstart1_all(lijkproc) = sum(dg_jsize_all(0:ljproc-1)) + 2
503                    dg_kstart1_all(lijkproc) = sum(dg_ksize_all(0:lkproc-1)) + 2
504                    dg_iend1_all(lijkproc) = dg_isize_all(liproc)+dg_istart1_all(lijkproc)-1
505                    dg_jend1_all(lijkproc) = dg_jsize_all(ljproc)+dg_jstart1_all(lijkproc)-1
506                    dg_kend1_all(lijkproc) = dg_ksize_all(lkproc)+dg_kstart1_all(lijkproc)-1
507                    dg_istart2_all(lijkproc) = dg_istart1_all(lijkproc)-1
508                    dg_jstart2_all(lijkproc) = dg_jstart1_all(lijkproc)-1
509                    dg_kstart2_all(lijkproc) = dg_kstart1_all(lijkproc)-1
510                    dg_iend2_all(lijkproc) =  dg_iend1_all(lijkproc)+1
511                    dg_jend2_all(lijkproc) = dg_jend1_all(lijkproc)+1
512                    dg_kend2_all(lijkproc) = dg_kend1_all(lijkproc)+1
513                    lijkproc = lijkproc + 1
514                 end do
515              end do
516           end do
517           if(no_k) then
518              dg_kstart1_all(:) = 1
519              dg_kend1_all(:) = 1
520              dg_kstart2_all(:) = 1
521              dg_kend2_all(:) = 1
522           end if
523     
524     ! set indices for global block
525           dg_imin2 = 1
526           dg_imin1 = 2
527           dg_imax1 = dg_imin1+sum(dg_isize_all(0:nodesi-1))-1
528           dg_imax2 = dg_imax1+1
529           dg_jmin2 = 1
530           dg_jmin1 = 2
531           dg_jmax1 = dg_jmin1+sum(dg_jsize_all(0:nodesj-1))-1
532           dg_jmax2 = dg_jmax1+1
533           if (no_k) then
534              dg_kmin2 = 1
535              dg_kmin1 = 1
536              dg_kmax1 = 1
537              dg_kmax2 = 1
538           else
539              dg_kmin2 = 1
540              dg_kmin1 = 2
541              dg_kmax1 = dg_kmin1+sum(dg_ksize_all(0:nodesk-1))-1
542              dg_kmax2 = dg_kmax1+1
543           end if
544     
545     ! set offset indices for periodic boundaries
546           lijkproc = mype
547           liproc = iofproc(lijkproc)
548           ljproc = jofproc(lijkproc)
549           lkproc = kofproc(lijkproc)
550           allocate(dg_cycoffset(2*dimn,3),icycoffset(2*dimn,3))
551           dg_cycoffset(:,:) = 0; icycoffset(:,:) = 0
552           if (des_periodic_walls_x) then
553              if(liproc.eq.0) then
554                 dg_cycoffset(2,1)= (dg_imax2-dg_imin1)
555                 icycoffset(2,1)= (imax2-imin1)
556              end if
557              if(liproc.eq.nodesi-1) then
558                 dg_cycoffset(1,1)=-(dg_imax2-dg_imin1)
559                 icycoffset(1,1)= -(imax2-imin1)
560              end if
561           end if
562           if (des_periodic_walls_y) then
563              if(ljproc.eq.0) then
564                 dg_cycoffset(4,2)= (dg_jmax2-dg_jmin1)
565                 icycoffset(4,2)= (jmax2-jmin1)
566              end if
567              if(ljproc.eq.nodesj-1) then
568                 dg_cycoffset(3,2)=-(dg_jmax2-dg_jmin1)
569                 icycoffset(3,2)= -(jmax2-jmin1)
570              end if
571           end if
572           if (des_periodic_walls_z) then
573              if(lkproc.eq.0) then
574                 dg_cycoffset(6,3)=(dg_kmax2-dg_kmin1)
575                 icycoffset(6,3)= (kmax2-kmin1)
576              end if
577              if(lkproc.eq.nodesk-1) then
578                 dg_cycoffset(5,3)=-(dg_kmax2-dg_kmin1)
579                 icycoffset(5,3)= -(kmax2-kmin1)
580              end if
581           end if
582     
583     
584     ! set the indices and variables for the current processor
585           dg_istart2 = dg_istart2_all(mype)
586           dg_istart1 = dg_istart1_all(mype)
587           dg_iend1 = dg_iend1_all(mype)
588           dg_iend2 = dg_iend2_all(mype)
589           dg_jstart2 = dg_jstart2_all(mype)
590           dg_jstart1 = dg_jstart1_all(mype)
591           dg_jend1 = dg_jend1_all(mype)
592           dg_jend2 = dg_jend2_all(mype)
593           dg_kstart2 = dg_kstart2_all(mype)
594           dg_kstart1 = dg_kstart1_all(mype)
595           dg_kend1 = dg_kend1_all(mype)
596           dg_kend2 = dg_kend2_all(mype)
597     
598     !       Defining new set of varaibles to define upper and lower bound of the
599     ! indices to include actual physical boundaries of the problem.
600           dg_istart = dg_istart1;   dg_iend = dg_iend1
601           dg_jstart = dg_jstart1;   dg_jend = dg_jend1
602           dg_kstart = dg_kstart1;   dg_kend = dg_kend1
603     
604           if(dg_istart .eq. dg_imin1) dg_istart = dg_imin2
605           if(dg_iend   .eq. dg_imax1) dg_iend   = dg_imax2
606     
607           if(dg_jstart .eq. dg_jmin1) dg_jstart = dg_jmin2
608           if(dg_jend   .eq. dg_jmax1) dg_jend   = dg_jmax2
609     
610           if(dg_kstart .eq. dg_kmin1) dg_kstart = dg_kmin2
611           if(dg_kend   .eq. dg_kmax1) dg_kend   = dg_kmax2
612     
613     ! set constants required for dg_funijk dg_funijk_gl for all processors
614           allocate(dg_c1_all(0:numpes-1),dg_c2_all(0:numpes-1),dg_c3_all(0:numpes-1))
615     
616           dg_c1_all=0;dg_c2_all=0;dg_c3_all=0
617           do lijkproc = 0,numpes-1
618              dg_c1_all(lijkproc) = (dg_jend2_all(lijkproc)-dg_jstart2_all(lijkproc)+1)
619              dg_c2_all(lijkproc) = dg_c1_all(lijkproc)*(dg_iend2_all(lijkproc)-dg_istart2_all(lijkproc)+1)
620              dg_c3_all(lijkproc) = -dg_c1_all(lijkproc)*dg_istart2_all(lijkproc) &
621                                       -dg_c2_all(lijkproc)*dg_kstart2_all(lijkproc) &
622                                       -dg_jstart2_all(lijkproc)+1
623           end do
624     ! set global constants
625           dg_c1_gl = (dg_jmax2-dg_jmin2+1)
626           dg_c2_gl = dg_c1_gl*(dg_imax2-dg_imin2+1)
627           dg_c3_gl = -dg_c1_gl*imin2-dg_c2_gl*kmin2-dg_jmin2+1
628     
629     ! set local constants
630           dg_c1_lo = dg_c1_all(mype)
631           dg_c2_lo = dg_c2_all(mype)
632           dg_c3_lo = dg_c3_all(mype)
633     
634           dg_ijksize2 = (dg_iend2-dg_istart2+1)* &
635                         (dg_jend2-dg_jstart2+1)* &
636                         (dg_kend2-dg_kstart2+1)
637           dg_ijkstart2 = dg_funijk(dg_istart2,dg_jstart2,dg_kstart2)
638           dg_ijkend2 = dg_funijk(dg_iend2,dg_jend2,dg_kend2)
639     
640     ! Confirmation checks
641           IF (DG_IJKSTART2.NE.1) THEN
642              WRITE(ERR_MSG,1100)'DG_IJKStart2'
643              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
644           ENDIF
645     
646           IF(DG_IJKEND2 /= DG_IJKSIZE2) THEN
647              WRITE(ERR_MSG,1100)'DG_IJKEnd2'
648              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
649           ENDIF
650     
651      1100 FORMAT('Error 1100: Invalid DG_IJKStart2. FATAL')
652     
653           IF(DG_IMIN1 > DG_IMAX1) THEN
654              WRITE(ERR_MSG,1105) 'DG_IMIN1 > DG_IMAX1'
655              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
656           ELSEIF(DG_JMIN1 > DG_JMAX1) THEN
657              WRITE(ERR_MSG,1105) 'DG_JMIN1 > DG_JMAX1'
658              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
659           ELSEIF(DG_KMIN1 > DG_KMAX1) THEN
660              WRITE(ERR_MSG,1105) 'DG_KMIN1 > DG_KMAX1'
661              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
662           ENDIF
663     
664      1105 FORMAT('Error 1105: Invalid DES grid indices: ',A,/'This is ',   &
665              'likely the result of automated grid calculations based',/    &
666              'on the maximum particle size. Specify the number of ',       &
667              'partitions',/'for the DES grid in the mfix.dat file ',       &
668              '(e.g., DESGRIDSEARCH_IMAX)')
669     
670     
671     ! set the domain length and dx,dy and dz values used in particle_in_cell
672     ! to bin the particles
673           dg_xstart = xe(istart2)
674           dg_xend = xe(iend1)
675           dg_ystart = yn(jstart2)
676           dg_yend = yn(jend1)
677           dg_dxinv = (dg_iend1-dg_istart1+1)/(dg_xend-dg_xstart)
678           dg_dyinv = (dg_jend1-dg_jstart1+1)/(dg_yend-dg_ystart)
679           if(no_k) then
680              dg_zstart = 1
681              dg_zend = 1
682              dg_dzinv = 1
683           else
684              dg_zstart = zt(kstart2)
685              dg_zend = zt(kend1)
686              dg_dzinv = (dg_kend1-dg_kstart1+1)/(dg_zend-dg_zstart)
687           end if
688     
689     ! allocate the desgridsearch related variables
690           allocate(dg_pic(dg_ijksize2))
691           dg_pic(:)%isize = 0
692           do lijk = 1,dg_ijksize2
693              allocate(dg_pic(lijk)%p(dg_pic_max_init))
694           end do
695     
696     !      call des_dbggrid
697           CALL FINL_ERR_MSG
698           end subroutine desgrid_init
699     
700     !------------------------------------------------------------------------
701     ! Subroutine       : desgrid_pic
702     ! Purpose          : it updates the particle in cell information
703     !                    and also located the particle
704     !                    moving across processor boundary
705     ! Parameters       : plocate - Locate the particles
706     ! Comments         : plocate should be set to true only when particle has to
707     !                    be located; if it is false only PIC information will be
708     !                    updated by this routine
709     !------------------------------------------------------------------------
710           subroutine desgrid_pic(plocate)
711     
712           use derived_types, only: dg_pic
713           use discretelement
714           use geometry, only: no_k
715     
716           implicit none
717     !-----------------------------------------------
718     ! Dummy arguments
719     !-----------------------------------------------
720           logical, INTENT(IN) :: plocate
721     !-----------------------------------------------
722     ! Local variables
723     !-----------------------------------------------
724           integer, dimension(:), allocatable :: lpic,lindx
725           integer li,lj,lk,lijk,lijk_count,lcurpar,lparcount,lcurpic
726           logical, save :: first_pass = .true.
727     !-----------------------------------------------
728     
729     ! locate the particles including ghost cells
730           max_isize = 0
731           lparcount = 1
732           allocate(lpic(dg_ijksize2))
733           allocate(lindx(dg_ijksize2))
734           lpic(:) = 0
735           if (plocate) then
736              dg_pijkprv(:)= dg_pijk(:)
737              do lcurpar = 1,max_pip
738                 if(lparcount.gt.pip) exit
739                 if(is_nonexistent(lcurpar)) cycle
740     
741                 lparcount = lparcount + 1
742                 li = min(dg_iend2,max(dg_istart2,iofpos(des_pos_new(lcurpar,1))))
743                 lj = min(dg_jend2,max(dg_jstart2,jofpos(des_pos_new(lcurpar,2))))
744                 if(no_k) then
745                    lk = 1
746                 else
747                    lk = min(dg_kend2,max(dg_kstart2,kofpos(des_pos_new(lcurpar,3))))
748                 end if
749                 dg_pijk(lcurpar) = dg_funijk(li,lj,lk)
750                 lijk = dg_pijk(lcurpar)
751                 lpic(lijk) = lpic(lijk) + 1
752              end do
753           else
754              do lcurpar = 1,max_pip
755                 if(lparcount.gt.pip) exit
756                 if(is_nonexistent(lcurpar)) cycle
757                 lparcount = lparcount + 1
758                 lijk = dg_pijk(lcurpar)
759                 lpic(lijk) = lpic(lijk) + 1
760              end do
761           end if
762     
763           if (first_pass) then
764              dg_pijkprv(:) = dg_pijk(:)
765              first_pass = .false.
766           end if
767     
768     ! redfine the array of dg_pic
769     !$omp parallel do default(shared)                               &
770     !$omp private(lijk,lcurpic) schedule (guided,50) reduction(max:max_isize)
771           do lijk = dg_ijkstart2,dg_ijkend2
772              lcurpic = lpic(lijk)
773              if(lcurpic > size(dg_pic(lijk)%p)) then
774                 deallocate(dg_pic(lijk)%p)
775                 allocate(dg_pic(lijk)%p(2*lcurpic))
776              end if
777              dg_pic(lijk)%isize = lcurpic
778              max_isize = max(max_isize,dg_pic(lijk)%isize)
779           end do
780     !$omp end parallel do
781     
782     ! assign the particle info in pic array
783           lindx(:) = 1
784           lparcount = 1
785     
786     #if ( defined(__GFORTRAN__) &&  ( __GNUC__ < 4 || ( __GNUC__ == 4 && __GNUC_MINOR__ < 7 ) ) ) \
787           ||  ( defined(__INTEL_COMPILER) && (__INTEL_COMPILER < 1400) )
788     
789           do lcurpar = 1, max_pip
790              if(lparcount.gt.pip) exit
791              if(is_nonexistent(lcurpar)) cycle
792              lparcount = lparcount + 1
793              lijk = dg_pijk(lcurpar)
794              dg_pic(lijk)%p(lindx(lijk)) = lcurpar
795              lindx(lijk) = lindx(lijk) +  1
796           enddo
797     #else
798     
799     !$omp parallel default(none) private(lcurpar,lijk,lijk_count) shared(max_pip,pip,lparcount,dg_pijk,dg_pic,lindx)
800     !$omp do
801           do lcurpar = 1, max_pip
802              if(lparcount.gt.pip) cycle
803              if(is_nonexistent(lcurpar)) cycle
804              lijk = dg_pijk(lcurpar)
805     
806     !$omp atomic capture
807              lijk_count = lindx(lijk)
808              lindx(lijk) = lindx(lijk) +  1
809     !$omp end atomic
810     
811              dg_pic(lijk)%p(lijk_count) = lcurpar
812     
813     !$omp atomic
814              lparcount = lparcount + 1
815           enddo
816     !$omp end do
817     !$omp end parallel
818     
819     #endif
820     
821     !      open (unit=100,file='desgrid.txt',status='unknown')
822     !      write(100,*)lpic
823     !      close(100)
824     
825           deallocate(lpic)
826           deallocate(lindx)
827     
828         contains
829     
830           include 'functions.inc'
831     
832         end subroutine desgrid_pic
833     
834     !------------------------------------------------------------------------
835     ! subroutine       : desgrid_neigh_build ()
836     ! Purpose          : This particles build the neigh list for the particles
837     !                    currently active in the system
838     !------------------------------------------------------------------------
839           subroutine desgrid_neigh_build ()
840     
841     !-----------------------------------------------
842     ! Modules
843     !-----------------------------------------------
844           Use des_thermo
845           use compar
846           use constant
847           use des_allocate, only: add_pair
848           use desmpi_wrapper
849           use derived_types, only: dg_pic
850           use discretelement
851           use funits
852           use geometry
853           use param1
854     
855           implicit none
856     !-----------------------------------------------
857     ! Local variables
858     !-----------------------------------------------
859           integer lcurpar,lkoffset
860           integer lijk,lic,ljc,lkc,li,lj,lk,ltotpic,lpicloc,lneigh,cc
861           double precision lsearch_rad,ldistsquared
862           double precision :: ldistvec(3)
863           double precision :: lcurpar_pos(3)
864           double precision :: lcur_off
865           integer il_off,iu_off,jl_off,ju_off,kl_off,ku_off
866           integer, dimension(:), allocatable :: tmp_neigh
867     !$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: int_tmp
868     !$    INTEGER :: PAIR_NUM_SMP,PAIR_MAX_SMP, tt, curr_tt, diff, mm, lSIZE2, dd
869     !$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: PAIRS_SMP
870     !$    INTEGER omp_get_num_threads
871     !$    INTEGER omp_get_thread_num
872     
873     !-----------------------------------------------
874     
875     ! loop through neighbours and build the contact particles list for particles
876     ! present in the system
877           lkoffset = dimn-2
878     
879     !$omp parallel default(none) private(lcurpar,lijk,lic,ljc,lkc,cc,curr_tt,diff, &
880     !$omp    il_off,iu_off,jl_off,ju_off,kl_off,ku_off,lcurpar_pos,lcur_off,lSIZE2,tmp_neigh,   &
881     !$omp    ltotpic, lneigh,lsearch_rad,ldistvec,ldistsquared, pair_num_smp, pair_max_smp, pairs_smp, int_tmp) &
882     !$omp    shared(max_pip,neighbors,neighbor_index,dg_pijk,NO_K,des_pos_new,dg_pic, factor_RLM,dd,  &
883     !$omp           des_radius, dg_xstart,dg_ystart,dg_zstart,dg_dxinv,dg_dyinv,dg_dzinv,dg_ijkstart2,dg_ijkend2, max_isize)
884     
885           allocate(tmp_neigh(max_isize))
886     
887     !$      PAIR_NUM_SMP = 0
888     !$      PAIR_MAX_SMP = 1024
889     !$      Allocate(  PAIRS_SMP(2,PAIR_MAX_SMP) )
890     
891     !$omp do
892     
893           do lcurpar =1,max_pip
894     
895     !$  if (.false.) then
896                 if (NEIGHBOR_INDEX(lcurpar) .eq. 0) then
897                     if (lcurpar .eq. 1) then
898                         NEIGHBOR_INDEX(lcurpar) = 1
899                     else
900                         NEIGHBOR_INDEX(lcurpar) = NEIGHBOR_INDEX(lcurpar-1)
901                     endif
902                 endif
903     !$  endif
904     
905              if (is_nonexistent(lcurpar) .or.is_entering(lcurpar) .or. is_entering_ghost(lcurpar) .or. is_ghost(lcurpar) .or. is_exiting_ghost(lcurpar)) cycle
906              lijk = dg_pijk(lcurpar)
907              lic = dg_iof_lo(lijk)
908              ljc = dg_jof_lo(lijk)
909              lkc = dg_kof_lo(lijk)
910     
911              il_off = 1
912              iu_off = 1
913              jl_off = 1
914              ju_off = 1
915              kl_off = 1
916              ku_off = 1
917     
918              lcurpar_pos(:) = des_pos_new(lcurpar,:)
919     !   The desgrid size should not be less than 2*dia*rlm_factor
920              lcur_off = (lcurpar_pos(1)-dg_xstart)*dg_dxinv - &
921                 floor((lcurpar_pos(1)-dg_xstart)*dg_dxinv)
922              if(lcur_off .ge. 0.5) then
923                 il_off = 0
924              else
925                 iu_off = 0
926              endif
927     
928              lcur_off = (lcurpar_pos(2)-dg_ystart)*dg_dyinv - &
929                 floor((lcurpar_pos(2)-dg_ystart)*dg_dyinv)
930              if(lcur_off .ge. 0.5) then
931                 jl_off = 0
932              else
933                 ju_off = 0
934              endif
935     
936              if(NO_K)then
937                 kl_off = 0
938                 ku_off = 0
939              else
940                lcur_off = (lcurpar_pos(3)-dg_zstart)*dg_dzinv - &
941                   floor((lcurpar_pos(3)-dg_zstart)*dg_dzinv)
942                if(lcur_off .ge. 0.5) then
943                   kl_off = 0
944                else
945                   ku_off = 0
946                endif
947             endif
948     
949              do lk = lkc-kl_off,lkc+ku_off
950              do lj = ljc-jl_off,ljc+ju_off
951              do li = lic-il_off,lic+iu_off
952                 lijk = dg_funijk(li,lj,lk)
953                 ltotpic =dg_pic(lijk)%isize
954                 do lpicloc = 1,ltotpic
955                    lneigh = dg_pic(lijk)%p(lpicloc)
956     
957                    tmp_neigh(lpicloc) = 0
958     ! Only skip real particles otherwise collisions with ghost, entering,
959     ! and exiting particles are missed.
960                    if (lneigh .eq. lcurpar) cycle
961                    if (lneigh < lcurpar .and.is_normal(lneigh)) cycle
962                    if (is_nonexistent(lneigh)) THEN
963                       cycle
964                    endif
965     
966                    lsearch_rad = factor_RLM*(des_radius(lcurpar)+des_radius(lneigh))
967                    ldistvec(1) = lcurpar_pos(1)-des_pos_new(lneigh,1)
968                    ldistvec(2) = lcurpar_pos(2)-des_pos_new(lneigh,2)
969                    ldistvec(3) = lcurpar_pos(3)-des_pos_new(lneigh,3)
970                    ldistsquared = dot_product(ldistvec,ldistvec)
971                    if (ldistsquared.gt.lsearch_rad*lsearch_rad) cycle
972                    tmp_neigh(lpicloc) = lneigh
973                 end do
974     
975                 do lpicloc = 1,ltotpic
976                    lneigh = tmp_neigh(lpicloc)
977                    if (lneigh .ne. 0) then
978     !$  if (.true.) then
979     !!!!! SMP version
980     
981     !$      pair_num_smp = pair_num_smp + 1
982     ! Reallocate to double the size of the arrays if needed.
983     !$      IF(PAIR_NUM_SMP > PAIR_MAX_SMP) THEN
984     !$         PAIR_MAX_SMP = 2*PAIR_MAX_SMP
985     !$         lSIZE2 = size(pairs_smp,2)
986     
987     !$         allocate(int_tmp(2,PAIR_MAX_SMP))
988     !$         int_tmp(:,1:lSIZE2) = pairs_smp(:,1:lSIZE2)
989     !$         call move_alloc(int_tmp,pairs_smp)
990     !$      ENDIF
991     
992     !$      pairs_smp(1,pair_num_smp) = lcurpar
993     !$      pairs_smp(2,pair_num_smp) = lneigh
994     
995     !$  else
996     !!!!! Serial version
997                       cc = add_pair(lcurpar, lneigh)
998     !$  endif
999                    endif
1000                 end do
1001              end do
1002              end do
1003              end do
1004           end do
1005     !$omp end do
1006     
1007     !$  curr_tt = omp_get_thread_num()+1  ! add one because thread numbering starts at zero
1008     
1009     !$omp single
1010     !$  NEIGHBOR_INDEX(1) = 1
1011     !$  dd = 1
1012     !$omp end single
1013     
1014     !$omp do ordered schedule(static,1)
1015     !$    do tt = 1, omp_get_num_threads()
1016     !$omp ordered
1017     !$        do MM = 1, PAIR_NUM_SMP
1018     !$            lcurpar = PAIRS_SMP(1,MM)
1019     !$            do while (dd .lt. lcurpar)
1020     !$                dd = dd + 1
1021     !$                NEIGHBOR_INDEX(dd) = NEIGHBOR_INDEX(dd-1)
1022     !$            enddo
1023     !$            cc = add_pair(lcurpar, PAIRS_SMP(2,MM))
1024     !$        enddo
1025     !$omp end ordered
1026     
1027     !$     end do
1028     !$omp  end do
1029     
1030     !$    deallocate( PAIRS_SMP )
1031     
1032           deallocate(tmp_neigh)
1033     
1034     !$omp end parallel
1035     
1036         contains
1037     
1038           include 'functions.inc'
1039     
1040         end subroutine desgrid_neigh_build
1041     
1042     !------------------------------------------------------------------------
1043     ! subroutine       : des_dbggrid
1044     ! Purpose          : For printing the indices used for desgrid
1045     !------------------------------------------------------------------------
1046           subroutine des_dbggrid()
1047     
1048           use param1
1049           use funits
1050           use geometry
1051           use compar
1052           use discretelement
1053           use constant
1054           use desmpi_wrapper
1055     
1056           implicit none
1057     !-----------------------------------------------
1058     ! Local variables
1059     !-----------------------------------------------
1060           integer lproc,liproc,ljproc,lkproc
1061           character (255) filename
1062     !-----------------------------------------------
1063     
1064           write(filename,'("dbg_desgridn",I4.4,".dat")') mype
1065           open(44,file=filename,convert='big_endian')
1066           do lproc = 0,numpes-1
1067              write(44,*) "Information for Proc =", lproc
1068              liproc= iofproc(lproc)
1069              ljproc= jofproc(lproc)
1070              lkproc= kofproc(lproc)
1071              write(44,*) "   "
1072              write(44,*) "i,j,k location of proc",liproc,ljproc,lkproc
1073              write(44,*) "i,j,k size of proc", dg_isize_all(liproc),dg_jsize_all(ljproc),dg_ksize_all(lkproc)
1074              write(44,*) "-------------------------------------------------"
1075              write(44,*) "indices        start     end"
1076              write(44,*) "-------------------------------------------------"
1077              write(44,*) "for i1:      ",dg_istart1_all(lproc),dg_iend1_all(lproc)
1078              write(44,*) "for i2:      ",dg_istart2_all(lproc),dg_iend2_all(lproc)
1079              write(44,*) "for j1:      ",dg_jstart1_all(lproc),dg_jend1_all(lproc)
1080              write(44,*) "for j2:      ",dg_jstart2_all(lproc),dg_jend2_all(lproc)
1081              write(44,*) "for k1:      ",dg_kstart1_all(lproc),dg_kend1_all(lproc)
1082              write(44,*) "for k2:      ",dg_kstart2_all(lproc),dg_kend2_all(lproc)
1083           end do
1084           write(44,*) "   "
1085           write(44,*) "-------------------------------------------------"
1086           write(44,*) "Local Start and end"
1087           write(44,*) "-------------------------------------------------"
1088           write(44,*) "for i :      ",dg_istart, dg_iend
1089           write(44,*) "for i1:      ",dg_istart1,dg_iend1
1090           write(44,*) "for i2:      ",dg_istart2,dg_iend2
1091           write(44,*) "   "
1092           write(44,*) "for j :      ",dg_jstart, dg_jend
1093           write(44,*) "for j1:      ",dg_jstart1,dg_jend1
1094           write(44,*) "for j2:      ",dg_jstart2,dg_jend2
1095           write(44,*) "   "
1096           write(44,*) "for k :      ",dg_kstart, dg_kend
1097           write(44,*) "for k1:      ",dg_kstart1,dg_kend1
1098           write(44,*) "for k2:      ",dg_kstart2,dg_kend2
1099           write(44,*) "   "
1100           write(44,*) "-------------------------------------------------"
1101           write(44,*) "global Start and end"
1102           write(44,*) "-------------------------------------------------"
1103           write(44,*) "for i1:      ",dg_imin1,dg_imax1
1104           write(44,*) "for i2:      ",dg_imin2,dg_imax2
1105           write(44,*) "for j1:      ",dg_jmin1,dg_jmax1
1106           write(44,*) "for j2:      ",dg_jmin2,dg_jmax2
1107           write(44,*) "for k1:      ",dg_kmin1,dg_kmax1
1108           write(44,*) "for k2:      ",dg_kmin2,dg_kmax2
1109           write(44,*) "   "
1110           write(44,*) "-------------------------------------------------"
1111           write(44,*) "dg_xstart:   ",dg_xstart
1112           write(44,*) "dg_xend:     ",dg_xend
1113           write(44,*) "dg_dxinv:    ",dg_dxinv
1114           write(44,*) "   "
1115           write(44,*) "dg_ystart:   ",dg_ystart
1116           write(44,*) "dg_yend:     ",dg_yend
1117           write(44,*) "dg_dyinv:    ",dg_dyinv
1118           write(44,*) "   "
1119           write(44,*) "dg_zstart:   ",dg_zstart
1120           write(44,*) "dg_zend:     ",dg_zend
1121           write(44,*) "dg_dzinv:    ",dg_dzinv
1122           write(44,*) "   "
1123           write(44,*) "dg_c1_lo/gl: ",dg_c1_lo, dg_c1_gl
1124           write(44,*) "dg_c2_lo/gl: ",dg_c2_lo, dg_c2_gl
1125           write(44,*) "dg_c3_lo/gl: ",dg_c3_lo, dg_c3_gl
1126           write(44,*) "   "
1127           write(44,*) "-------------------------------------------------"
1128     
1129     
1130     
1131           close(44)
1132           end subroutine des_dbggrid
1133     
1134           end module
1135