File: RELATIVE:/../../../mfix.git/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 discretelement, 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 discretelement
713           use geometry, only: no_k
714     
715           implicit none
716     !-----------------------------------------------
717     ! Dummy arguments
718     !-----------------------------------------------
719           logical, INTENT(IN) :: plocate
720     !-----------------------------------------------
721     ! Local variables
722     !-----------------------------------------------
723           integer, dimension(dg_ijksize2) :: lpic,lindx
724           integer li,lj,lk,lijk,lijk_count,lcurpar,lparcount,lcurpic
725           logical, save :: first_pass = .true.
726     !-----------------------------------------------
727     
728     ! locate the particles including ghost cells
729           max_isize = 0
730           lparcount = 1
731           lpic(:) = 0
732           if (plocate) then
733              dg_pijkprv(:)= dg_pijk(:)
734              do lcurpar = 1,max_pip
735                 if(lparcount.gt.pip) exit
736                 if(is_nonexistent(lcurpar)) cycle
737     
738                 lparcount = lparcount + 1
739                 li = min(dg_iend2,max(dg_istart2,iofpos(des_pos_new(1,lcurpar))))
740                 lj = min(dg_jend2,max(dg_jstart2,jofpos(des_pos_new(2,lcurpar))))
741                 if(no_k) then
742                    lk = 1
743                 else
744                    lk = min(dg_kend2,max(dg_kstart2,kofpos(des_pos_new(3,lcurpar))))
745                 end if
746                 dg_pijk(lcurpar) = dg_funijk(li,lj,lk)
747                 lijk = dg_pijk(lcurpar)
748                 lpic(lijk) = lpic(lijk) + 1
749              end do
750           else
751              do lcurpar = 1,max_pip
752                 if(lparcount.gt.pip) exit
753                 if(is_nonexistent(lcurpar)) cycle
754                 lparcount = lparcount + 1
755                 lijk = dg_pijk(lcurpar)
756                 lpic(lijk) = lpic(lijk) + 1
757              end do
758           end if
759     
760           if (first_pass) then
761              dg_pijkprv(:) = dg_pijk(:)
762              first_pass = .false.
763           end if
764     
765     ! redfine the array of dg_pic
766     !!$omp parallel do default(shared)                               &
767     !!$omp private(lijk,lcurpic) schedule (guided,50)
768           do lijk = dg_ijkstart2,dg_ijkend2
769              lcurpic = lpic(lijk)
770              if(lcurpic > size(dg_pic(lijk)%p)) then
771                 deallocate(dg_pic(lijk)%p)
772                 allocate(dg_pic(lijk)%p(2*lcurpic))
773              end if
774              dg_pic(lijk)%isize = lcurpic
775              max_isize = max(max_isize,dg_pic(lijk)%isize)
776           end do
777     !!$omp end parallel do
778     
779     ! assign the particle info in pic array
780           lindx(:) = 1
781           lparcount = 1
782     
783     #if ( defined(__GFORTRAN__) &&  ( __GNUC__ < 4 || ( __GNUC__ == 4 && __GNUC_MINOR__ < 7 ) ) ) \
784           ||  ( defined(__INTEL_COMPILER) && (__INTEL_COMPILER < 1400) )
785     
786           do lcurpar = 1, max_pip
787              if(lparcount.gt.pip) exit
788              if(is_nonexistent(lcurpar)) cycle
789              lparcount = lparcount + 1
790              lijk = dg_pijk(lcurpar)
791              dg_pic(lijk)%p(lindx(lijk)) = lcurpar
792              lindx(lijk) = lindx(lijk) +  1
793           enddo
794     #else
795     
796     !$omp parallel default(none) private(lcurpar,lijk,lijk_count) shared(max_pip,pip,lparcount,dg_pijk,dg_pic,lindx)
797     !$omp do
798           do lcurpar = 1, max_pip
799              if(lparcount.gt.pip) cycle
800              if(is_nonexistent(lcurpar)) cycle
801              lijk = dg_pijk(lcurpar)
802     
803              !$omp atomic capture
804              lijk_count = lindx(lijk)
805              lindx(lijk) = lindx(lijk) +  1
806              !$omp end atomic
807     
808              dg_pic(lijk)%p(lijk_count) = lcurpar
809     
810              !$omp atomic
811              lparcount = lparcount + 1
812           enddo
813     !$omp end do
814     !$omp end parallel
815     
816     #endif
817     
818     !      open (unit=100,file='desgrid.txt',status='unknown')
819     !      write(100,*)lpic
820     !      close(100)
821     
822         contains
823     
824           include 'functions.inc'
825     
826         end subroutine desgrid_pic
827     
828     !------------------------------------------------------------------------
829     ! subroutine       : desgrid_neigh_build ()
830     ! Purpose          : This particles build the neigh list for the particles
831     !                    currently active in the system
832     !------------------------------------------------------------------------
833           subroutine desgrid_neigh_build ()
834     
835     !-----------------------------------------------
836     ! Modules
837     !-----------------------------------------------
838           Use des_thermo
839           use compar
840           use constant
841           use des_allocate, only: add_pair
842           use desmpi_wrapper
843           use discretelement
844           use funits
845           use geometry
846           use param1
847     
848           implicit none
849     !-----------------------------------------------
850     ! Local variables
851     !-----------------------------------------------
852           integer lcurpar,lkoffset
853           integer lijk,lic,ljc,lkc,li,lj,lk,ltotpic,lpicloc,lneigh,cc
854           double precision lsearch_rad,ldistsquared
855           double precision :: ldistvec(3)
856           double precision :: lcurpar_pos(3)
857           double precision :: lcur_off
858           integer il_off,iu_off,jl_off,ju_off,kl_off,ku_off
859           integer, dimension(:), allocatable :: tmp_neigh
860     !$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: int_tmp
861     !$    INTEGER :: PAIR_NUM_SMP,PAIR_MAX_SMP, tt, curr_tt, diff, mm, lSIZE2, dd
862     !$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: PAIRS_SMP
863     !$    INTEGER omp_get_num_threads
864     !$    INTEGER omp_get_thread_num
865     
866     !-----------------------------------------------
867     
868     ! loop through neighbours and build the contact particles list for particles
869     ! present in the system
870           lkoffset = dimn-2
871     
872     !$omp parallel default(none) private(lcurpar,lijk,lic,ljc,lkc,cc,curr_tt,diff, &
873     !$omp    il_off,iu_off,jl_off,ju_off,kl_off,ku_off,lcurpar_pos,lcur_off,lSIZE2,tmp_neigh,   &
874     !$omp    ltotpic, lneigh,lsearch_rad,ldistvec,ldistsquared, pair_num_smp, pair_max_smp, pairs_smp, int_tmp) &
875     !$omp    shared(max_pip,neighbors,neighbor_index,neigh_max,dg_pijk,NO_K,des_pos_new,dg_pic, factor_RLM,dd,  &
876     !$omp           des_radius, dg_xstart,dg_ystart,dg_zstart,dg_dxinv,dg_dyinv,dg_dzinv,dg_ijkstart2,dg_ijkend2, max_isize)
877     
878           allocate(tmp_neigh(max_isize))
879     
880     !$      PAIR_NUM_SMP = 0
881     !$      PAIR_MAX_SMP = 1024
882     !$      Allocate(  PAIRS_SMP(2,PAIR_MAX_SMP) )
883     
884     !$omp do
885     
886           do lcurpar =1,max_pip
887     
888     !$  if (.false.) then
889                 if (NEIGHBOR_INDEX(lcurpar) .eq. 0) then
890                     if (lcurpar .eq. 1) then
891                         NEIGHBOR_INDEX(lcurpar) = 1
892                     else
893                         NEIGHBOR_INDEX(lcurpar) = NEIGHBOR_INDEX(lcurpar-1)
894                     endif
895                 endif
896     !$  endif
897     
898              if (is_nonexistent(lcurpar) .or.is_entering(lcurpar) .or. is_entering_ghost(lcurpar) .or. is_ghost(lcurpar) .or. is_exiting_ghost(lcurpar)) cycle
899              lijk = dg_pijk(lcurpar)
900              lic = dg_iof_lo(lijk)
901              ljc = dg_jof_lo(lijk)
902              lkc = dg_kof_lo(lijk)
903     
904              il_off = 1
905              iu_off = 1
906              jl_off = 1
907              ju_off = 1
908              kl_off = 1
909              ku_off = 1
910     
911              lcurpar_pos(:) = des_pos_new(:,lcurpar)
912     !   The desgrid size should not be less than 2*dia*rlm_factor
913              lcur_off = (lcurpar_pos(1)-dg_xstart)*dg_dxinv - &
914                 floor((lcurpar_pos(1)-dg_xstart)*dg_dxinv)
915              if(lcur_off .ge. 0.5) then
916                 il_off = 0
917              else
918                 iu_off = 0
919              endif
920     
921              lcur_off = (lcurpar_pos(2)-dg_ystart)*dg_dyinv - &
922                 floor((lcurpar_pos(2)-dg_ystart)*dg_dyinv)
923              if(lcur_off .ge. 0.5) then
924                 jl_off = 0
925              else
926                 ju_off = 0
927              endif
928     
929              if(NO_K)then
930                 kl_off = 0
931                 ku_off = 0
932              else
933                lcur_off = (lcurpar_pos(3)-dg_zstart)*dg_dzinv - &
934                   floor((lcurpar_pos(3)-dg_zstart)*dg_dzinv)
935                if(lcur_off .ge. 0.5) then
936                   kl_off = 0
937                else
938                   ku_off = 0
939                endif
940             endif
941     
942              do lk = lkc-kl_off,lkc+ku_off
943              do lj = ljc-jl_off,ljc+ju_off
944              do li = lic-il_off,lic+iu_off
945                 lijk = dg_funijk(li,lj,lk)
946                 ltotpic =dg_pic(lijk)%isize
947                 do lpicloc = 1,ltotpic
948                    lneigh = dg_pic(lijk)%p(lpicloc)
949     
950                    tmp_neigh(lpicloc) = 0
951     ! Only skip real particles otherwise collisions with ghost, entering,
952     ! and exiting particles are missed.
953                    if (lneigh .eq. lcurpar) cycle
954                    if (lneigh < lcurpar .and.is_normal(lneigh)) cycle
955                    if (is_nonexistent(lneigh)) THEN
956                       cycle
957                    endif
958     
959                    lsearch_rad = factor_RLM*(des_radius(lcurpar)+des_radius(lneigh))
960                    ldistvec(1) = lcurpar_pos(1)-des_pos_new(1,lneigh)
961                    ldistvec(2) = lcurpar_pos(2)-des_pos_new(2,lneigh)
962                    ldistvec(3) = lcurpar_pos(3)-des_pos_new(3,lneigh)
963                    ldistsquared = dot_product(ldistvec,ldistvec)
964                    if (ldistsquared.gt.lsearch_rad*lsearch_rad) cycle
965                    tmp_neigh(lpicloc) = lneigh
966                 end do
967     
968                 do lpicloc = 1,ltotpic
969                    lneigh = tmp_neigh(lpicloc)
970                    if (lneigh .ne. 0) then
971     !$  if (.true.) then
972     !!!!! SMP version
973     
974     !$      pair_num_smp = pair_num_smp + 1
975     ! Reallocate to double the size of the arrays if needed.
976     !$      IF(PAIR_NUM_SMP > PAIR_MAX_SMP) THEN
977     !$         PAIR_MAX_SMP = 2*PAIR_MAX_SMP
978     !$         lSIZE2 = size(pairs_smp,2)
979     
980     !$         allocate(int_tmp(2,PAIR_MAX_SMP))
981     !$         int_tmp(:,1:lSIZE2) = pairs_smp(:,1:lSIZE2)
982     !$         call move_alloc(int_tmp,pairs_smp)
983     !$      ENDIF
984     
985     !$      pairs_smp(1,pair_num_smp) = lcurpar
986     !$      pairs_smp(2,pair_num_smp) = lneigh
987     
988     !$  else
989     !!!!! Serial version
990                       cc = add_pair(lcurpar, lneigh)
991     !$  endif
992                    endif
993                 end do
994              end do
995              end do
996              end do
997           end do
998     !$omp end do
999     
1000     !$  curr_tt = omp_get_thread_num()+1  ! add one because thread numbering starts at zero
1001     
1002     !$omp single
1003     !$  NEIGHBOR_INDEX(1) = 1
1004     !$  dd = 1
1005     !$omp end single
1006     
1007     !$omp do ordered schedule(static,1)
1008     !$    do tt = 1, omp_get_num_threads()
1009     !$omp ordered
1010     !$        do MM = 1, PAIR_NUM_SMP
1011     !$            lcurpar = PAIRS_SMP(1,MM)
1012     !$            do while (dd .lt. lcurpar)
1013     !$                dd = dd + 1
1014     !$                NEIGHBOR_INDEX(dd) = NEIGHBOR_INDEX(dd-1)
1015     !$            enddo
1016     !$            cc = add_pair(lcurpar, PAIRS_SMP(2,MM))
1017     !$        enddo
1018     !$omp end ordered
1019     
1020     !$     end do
1021     !$omp  end do
1022     
1023     !$    deallocate( PAIRS_SMP )
1024     
1025           deallocate(tmp_neigh)
1026     
1027     !$omp end parallel
1028     
1029         contains
1030     
1031           include 'functions.inc'
1032     
1033         end subroutine desgrid_neigh_build
1034     
1035     !------------------------------------------------------------------------
1036     ! subroutine       : des_dbggrid
1037     ! Purpose          : For printing the indices used for desgrid
1038     !------------------------------------------------------------------------
1039           subroutine des_dbggrid()
1040     
1041           use param1
1042           use funits
1043           use geometry
1044           use compar
1045           use discretelement
1046           use constant
1047           use desmpi_wrapper
1048     
1049           implicit none
1050     !-----------------------------------------------
1051     ! Local variables
1052     !-----------------------------------------------
1053           integer lproc,liproc,ljproc,lkproc
1054           character (255) filename
1055     !-----------------------------------------------
1056     
1057           write(filename,'("dbg_desgridn",I4.4,".dat")') mype
1058           open(44,file=filename,convert='big_endian')
1059           do lproc = 0,numpes-1
1060              write(44,*) "Information for Proc =", lproc
1061              liproc= iofproc(lproc)
1062              ljproc= jofproc(lproc)
1063              lkproc= kofproc(lproc)
1064              write(44,*) "   "
1065              write(44,*) "i,j,k location of proc",liproc,ljproc,lkproc
1066              write(44,*) "i,j,k size of proc", dg_isize_all(liproc),dg_jsize_all(ljproc),dg_ksize_all(lkproc)
1067              write(44,*) "-------------------------------------------------"
1068              write(44,*) "indices        start     end"
1069              write(44,*) "-------------------------------------------------"
1070              write(44,*) "for i1:      ",dg_istart1_all(lproc),dg_iend1_all(lproc)
1071              write(44,*) "for i2:      ",dg_istart2_all(lproc),dg_iend2_all(lproc)
1072              write(44,*) "for j1:      ",dg_jstart1_all(lproc),dg_jend1_all(lproc)
1073              write(44,*) "for j2:      ",dg_jstart2_all(lproc),dg_jend2_all(lproc)
1074              write(44,*) "for k1:      ",dg_kstart1_all(lproc),dg_kend1_all(lproc)
1075              write(44,*) "for k2:      ",dg_kstart2_all(lproc),dg_kend2_all(lproc)
1076           end do
1077           write(44,*) "   "
1078           write(44,*) "-------------------------------------------------"
1079           write(44,*) "Local Start and end"
1080           write(44,*) "-------------------------------------------------"
1081           write(44,*) "for i :      ",dg_istart, dg_iend
1082           write(44,*) "for i1:      ",dg_istart1,dg_iend1
1083           write(44,*) "for i2:      ",dg_istart2,dg_iend2
1084           write(44,*) "   "
1085           write(44,*) "for j :      ",dg_jstart, dg_jend
1086           write(44,*) "for j1:      ",dg_jstart1,dg_jend1
1087           write(44,*) "for j2:      ",dg_jstart2,dg_jend2
1088           write(44,*) "   "
1089           write(44,*) "for k :      ",dg_kstart, dg_kend
1090           write(44,*) "for k1:      ",dg_kstart1,dg_kend1
1091           write(44,*) "for k2:      ",dg_kstart2,dg_kend2
1092           write(44,*) "   "
1093           write(44,*) "-------------------------------------------------"
1094           write(44,*) "global Start and end"
1095           write(44,*) "-------------------------------------------------"
1096           write(44,*) "for i1:      ",dg_imin1,dg_imax1
1097           write(44,*) "for i2:      ",dg_imin2,dg_imax2
1098           write(44,*) "for j1:      ",dg_jmin1,dg_jmax1
1099           write(44,*) "for j2:      ",dg_jmin2,dg_jmax2
1100           write(44,*) "for k1:      ",dg_kmin1,dg_kmax1
1101           write(44,*) "for k2:      ",dg_kmin2,dg_kmax2
1102           write(44,*) "   "
1103           write(44,*) "-------------------------------------------------"
1104           write(44,*) "dg_xstart:   ",dg_xstart
1105           write(44,*) "dg_xend:     ",dg_xend
1106           write(44,*) "dg_dxinv:    ",dg_dxinv
1107           write(44,*) "   "
1108           write(44,*) "dg_ystart:   ",dg_ystart
1109           write(44,*) "dg_yend:     ",dg_yend
1110           write(44,*) "dg_dyinv:    ",dg_dyinv
1111           write(44,*) "   "
1112           write(44,*) "dg_zstart:   ",dg_zstart
1113           write(44,*) "dg_zend:     ",dg_zend
1114           write(44,*) "dg_dzinv:    ",dg_dzinv
1115           write(44,*) "   "
1116           write(44,*) "dg_c1_lo/gl: ",dg_c1_lo, dg_c1_gl
1117           write(44,*) "dg_c2_lo/gl: ",dg_c2_lo, dg_c2_gl
1118           write(44,*) "dg_c3_lo/gl: ",dg_c3_lo, dg_c3_gl
1119           write(44,*) "   "
1120           write(44,*) "-------------------------------------------------"
1121     
1122     
1123     
1124           close(44)
1125           end subroutine des_dbggrid
1126     
1127           end module
1128