File: /nfs/home/0/users/jenkins/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, nodesk
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     ! Maximum particle size.
444           use discretelement, only: MAX_RADIUS
445     
446           use discretelement, only: dg_pic
447     
448     ! Global Parameters:
449     !---------------------------------------------------------------------//
450           use param1, only: UNDEFINED_I
451     
452     ! Use the error manager for posting error messages.
453     !---------------------------------------------------------------------//
454           use error_manager
455     
456     
457           implicit none
458     !-----------------------------------------------
459     ! Local varibles
460     !-----------------------------------------------
461           double precision :: ltempdx,ltempdy,ltempdz
462           integer :: lijkproc,liproc,ljproc,lkproc
463           integer :: lijk
464     !......................................................................!
465     
466     ! Initialize the error manager.
467           CALL INIT_ERR_MSG("DESGRID_INIT")
468     
469     ! set indices for all processors
470          allocate (dg_istart1_all(0:numpes-1), dg_iend1_all(0:numpes-1))
471          allocate (dg_istart2_all(0:numpes-1), dg_iend2_all(0:numpes-1))
472          allocate (dg_isize_all(0:nodesi-1))
473     
474          allocate (dg_jstart1_all(0:numpes-1), dg_jend1_all(0:numpes-1))
475          allocate (dg_jstart2_all(0:numpes-1), dg_jend2_all(0:numpes-1))
476          allocate (dg_jsize_all(0:nodesj-1))
477     
478          allocate (dg_kstart1_all(0:numpes-1), dg_kend1_all(0:numpes-1))
479          allocate (dg_kstart2_all(0:numpes-1), dg_kend2_all(0:numpes-1))
480          allocate (dg_ksize_all(0:nodesk-1))
481     
482     
483           dg_istart1_all=0; dg_iend1_all=0
484           dg_istart2_all=0; dg_iend2_all=0
485           dg_isize_all=0
486     
487           dg_jstart1_all=0; dg_jend1_all=0
488           dg_jstart2_all=0; dg_jend2_all=0
489           dg_jsize_all=0
490     
491           dg_kstart1_all=0; dg_kend1_all=0
492           dg_kstart2_all=0; dg_kend2_all=0
493           dg_ksize_all=0
494     
495     
496     ! set grid size based on user input desgridsearch_<ijk>max
497           ltempdx = xlength/desgridsearch_imax
498           ltempdy = ylength/desgridsearch_jmax
499           if(do_k) ltempdz = zlength/desgridsearch_kmax
500           dg_ksize_all(:) = 1
501     
502           lijkproc = 0
503           do lkproc=0,nodesk-1
504              do ljproc=0,nodesj-1
505                 do liproc=0,nodesi-1
506                    dg_isize_all(liproc) = (xe(iend1_all(lijkproc))-xe(istart1_all(lijkproc)-1))/ltempdx
507                    dg_jsize_all(ljproc) = (yn(jend1_all(lijkproc))-yn(jstart1_all(lijkproc)-1))/ltempdy
508                    if(do_k) dg_ksize_all(lkproc) = (zt(kend1_all(lijkproc))-zt(kstart1_all(lijkproc)-1))/ltempdz
509                    dg_istart1_all(lijkproc) = sum(dg_isize_all(0:liproc-1)) + 2
510                    dg_jstart1_all(lijkproc) = sum(dg_jsize_all(0:ljproc-1)) + 2
511                    dg_kstart1_all(lijkproc) = sum(dg_ksize_all(0:lkproc-1)) + 2
512                    dg_iend1_all(lijkproc) = dg_isize_all(liproc)+dg_istart1_all(lijkproc)-1
513                    dg_jend1_all(lijkproc) = dg_jsize_all(ljproc)+dg_jstart1_all(lijkproc)-1
514                    dg_kend1_all(lijkproc) = dg_ksize_all(lkproc)+dg_kstart1_all(lijkproc)-1
515                    dg_istart2_all(lijkproc) = dg_istart1_all(lijkproc)-1
516                    dg_jstart2_all(lijkproc) = dg_jstart1_all(lijkproc)-1
517                    dg_kstart2_all(lijkproc) = dg_kstart1_all(lijkproc)-1
518                    dg_iend2_all(lijkproc) =  dg_iend1_all(lijkproc)+1
519                    dg_jend2_all(lijkproc) = dg_jend1_all(lijkproc)+1
520                    dg_kend2_all(lijkproc) = dg_kend1_all(lijkproc)+1
521                    lijkproc = lijkproc + 1
522                 end do
523              end do
524           end do
525           if(no_k) then
526              dg_kstart1_all(:) = 1
527              dg_kend1_all(:) = 1
528              dg_kstart2_all(:) = 1
529              dg_kend2_all(:) = 1
530           end if
531     
532     ! set indices for global block
533           dg_imin2 = 1
534           dg_imin1 = 2
535           dg_imax1 = dg_imin1+sum(dg_isize_all(0:nodesi-1))-1
536           dg_imax2 = dg_imax1+1
537           dg_jmin2 = 1
538           dg_jmin1 = 2
539           dg_jmax1 = dg_jmin1+sum(dg_jsize_all(0:nodesj-1))-1
540           dg_jmax2 = dg_jmax1+1
541           if (no_k) then
542              dg_kmin2 = 1
543              dg_kmin1 = 1
544              dg_kmax1 = 1
545              dg_kmax2 = 1
546           else
547              dg_kmin2 = 1
548              dg_kmin1 = 2
549              dg_kmax1 = dg_kmin1+sum(dg_ksize_all(0:nodesk-1))-1
550              dg_kmax2 = dg_kmax1+1
551           end if
552     
553     ! set offset indices for periodic boundaries
554           lijkproc = mype
555           liproc = iofproc(lijkproc)
556           ljproc = jofproc(lijkproc)
557           lkproc = kofproc(lijkproc)
558           allocate(dg_cycoffset(2*dimn,3),icycoffset(2*dimn,3))
559           dg_cycoffset(:,:) = 0; icycoffset(:,:) = 0
560           if (des_periodic_walls_x) then
561              if(liproc.eq.0) then
562                 dg_cycoffset(2,1)= (dg_imax2-dg_imin1)
563                 icycoffset(2,1)= (imax2-imin1)
564              end if
565              if(liproc.eq.nodesi-1) then
566                 dg_cycoffset(1,1)=-(dg_imax2-dg_imin1)
567                 icycoffset(1,1)= -(imax2-imin1)
568              end if
569           end if
570           if (des_periodic_walls_y) then
571              if(ljproc.eq.0) then
572                 dg_cycoffset(4,2)= (dg_jmax2-dg_jmin1)
573                 icycoffset(4,2)= (jmax2-jmin1)
574              end if
575              if(ljproc.eq.nodesj-1) then
576                 dg_cycoffset(3,2)=-(dg_jmax2-dg_jmin1)
577                 icycoffset(3,2)= -(jmax2-jmin1)
578              end if
579           end if
580           if (des_periodic_walls_z) then
581              if(lkproc.eq.0) then
582                 dg_cycoffset(6,3)=(dg_kmax2-dg_kmin1)
583                 icycoffset(6,3)= (kmax2-kmin1)
584              end if
585              if(lkproc.eq.nodesk-1) then
586                 dg_cycoffset(5,3)=-(dg_kmax2-dg_kmin1)
587                 icycoffset(5,3)= -(kmax2-kmin1)
588              end if
589           end if
590     
591     
592     ! set the indices and variables for the current processor
593           dg_istart2 = dg_istart2_all(mype)
594           dg_istart1 = dg_istart1_all(mype)
595           dg_iend1 = dg_iend1_all(mype)
596           dg_iend2 = dg_iend2_all(mype)
597           dg_jstart2 = dg_jstart2_all(mype)
598           dg_jstart1 = dg_jstart1_all(mype)
599           dg_jend1 = dg_jend1_all(mype)
600           dg_jend2 = dg_jend2_all(mype)
601           dg_kstart2 = dg_kstart2_all(mype)
602           dg_kstart1 = dg_kstart1_all(mype)
603           dg_kend1 = dg_kend1_all(mype)
604           dg_kend2 = dg_kend2_all(mype)
605     
606     !       Defining new set of varaibles to define upper and lower bound of the
607     ! indices to include actual physical boundaries of the problem.
608           dg_istart = dg_istart1;   dg_iend = dg_iend1
609           dg_jstart = dg_jstart1;   dg_jend = dg_jend1
610           dg_kstart = dg_kstart1;   dg_kend = dg_kend1
611     
612           if(dg_istart .eq. dg_imin1) dg_istart = dg_imin2
613           if(dg_iend   .eq. dg_imax1) dg_iend   = dg_imax2
614     
615           if(dg_jstart .eq. dg_jmin1) dg_jstart = dg_jmin2
616           if(dg_jend   .eq. dg_jmax1) dg_jend   = dg_jmax2
617     
618           if(dg_kstart .eq. dg_kmin1) dg_kstart = dg_kmin2
619           if(dg_kend   .eq. dg_kmax1) dg_kend   = dg_kmax2
620     
621     ! set constants required for dg_funijk dg_funijk_gl for all processors
622           allocate(dg_c1_all(0:numpes-1),dg_c2_all(0:numpes-1),dg_c3_all(0:numpes-1))
623     
624           dg_c1_all=0;dg_c2_all=0;dg_c3_all=0
625           do lijkproc = 0,numpes-1
626              dg_c1_all(lijkproc) = (dg_jend2_all(lijkproc)-dg_jstart2_all(lijkproc)+1)
627              dg_c2_all(lijkproc) = dg_c1_all(lijkproc)*(dg_iend2_all(lijkproc)-dg_istart2_all(lijkproc)+1)
628              dg_c3_all(lijkproc) = -dg_c1_all(lijkproc)*dg_istart2_all(lijkproc) &
629                                       -dg_c2_all(lijkproc)*dg_kstart2_all(lijkproc) &
630                                       -dg_jstart2_all(lijkproc)+1
631           end do
632     ! set global constants
633           dg_c1_gl = (dg_jmax2-dg_jmin2+1)
634           dg_c2_gl = dg_c1_gl*(dg_imax2-dg_imin2+1)
635           dg_c3_gl = -dg_c1_gl*imin2-dg_c2_gl*kmin2-dg_jmin2+1
636     
637     ! set local constants
638           dg_c1_lo = dg_c1_all(mype)
639           dg_c2_lo = dg_c2_all(mype)
640           dg_c3_lo = dg_c3_all(mype)
641     
642           dg_ijksize2 = (dg_iend2-dg_istart2+1)* &
643                         (dg_jend2-dg_jstart2+1)* &
644                         (dg_kend2-dg_kstart2+1)
645           dg_ijkstart2 = dg_funijk(dg_istart2,dg_jstart2,dg_kstart2)
646           dg_ijkend2 = dg_funijk(dg_iend2,dg_jend2,dg_kend2)
647     
648     ! Confirmation checks
649           IF (DG_IJKSTART2.NE.1) THEN
650              WRITE(ERR_MSG,1100)'DG_IJKStart2'
651              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
652           ENDIF
653     
654           IF(DG_IJKEND2 /= DG_IJKSIZE2) THEN
655              WRITE(ERR_MSG,1100)'DG_IJKEnd2'
656              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
657           ENDIF
658     
659      1100 FORMAT('Error 1100: Invalid DG_IJKStart2. FATAL')
660     
661           IF(DG_IMIN1 > DG_IMAX1) THEN
662              WRITE(ERR_MSG,1105) 'DG_IMIN1 > DG_IMAX1'
663              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
664           ELSEIF(DG_JMIN1 > DG_JMAX1) THEN
665              WRITE(ERR_MSG,1105) 'DG_JMIN1 > DG_JMAX1'
666              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
667           ELSEIF(DG_KMIN1 > DG_KMAX1) THEN
668              WRITE(ERR_MSG,1105) 'DG_KMIN1 > DG_KMAX1'
669              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
670           ENDIF
671     
672      1105 FORMAT('Error 1105: Invalid DES grid indices: ',A,/'This is ',   &
673              'likely the result of automated grid calculations based',/    &
674              'on the maximum particle size. Specify the number of ',       &
675              'partitions',/'for the DES grid in the mfix.dat file ',       &
676              '(e.g., DESGRIDSEARCH_IMAX)')
677     
678     
679     ! set the domain length and dx,dy and dz values used in particle_in_cell
680     ! to bin the particles
681           dg_xstart = xe(istart2)
682           dg_xend = xe(iend1)
683           dg_ystart = yn(jstart2)
684           dg_yend = yn(jend1)
685           dg_dxinv = (dg_iend1-dg_istart1+1)/(dg_xend-dg_xstart)
686           dg_dyinv = (dg_jend1-dg_jstart1+1)/(dg_yend-dg_ystart)
687           if(no_k) then
688              dg_zstart = 1
689              dg_zend = 1
690              dg_dzinv = 1
691           else
692              dg_zstart = zt(kstart2)
693              dg_zend = zt(kend1)
694              dg_dzinv = (dg_kend1-dg_kstart1+1)/(dg_zend-dg_zstart)
695           end if
696     
697     ! allocate the desgridsearch related variables
698           allocate(dg_pic(dg_ijksize2))
699           dg_pic(:)%isize = 0
700           do lijk = 1,dg_ijksize2
701              allocate(dg_pic(lijk)%p(dg_pic_max_init))
702           end do
703     
704     !      call des_dbggrid
705           CALL FINL_ERR_MSG
706           end subroutine desgrid_init
707     
708     !------------------------------------------------------------------------
709     ! Subroutine       : desgrid_pic
710     ! Purpose          : it updates the particle in cell information
711     !                    and also located the particle
712     !                    moving across processor boundary
713     ! Parameters       : plocate - Locate the particles
714     ! Comments         : plocate should be set to true only when particle has to
715     !                    be located; if it is false only PIC information will be
716     !                    updated by this routine
717     !------------------------------------------------------------------------
718           subroutine desgrid_pic(plocate)
719     
720           use param1
721           use funits
722           use geometry
723           use compar
724           use discretelement
725           use constant
726           use desmpi_wrapper
727     
728     
729           implicit none
730     !-----------------------------------------------
731     ! Dummy arguments
732     !-----------------------------------------------
733           logical, INTENT(IN) :: plocate
734     !-----------------------------------------------
735     ! Local variables
736     !-----------------------------------------------
737           integer, dimension(dg_ijksize2) :: lpic,lindx
738           integer li,lj,lk,lijk,lijk_count,lcurpar,lparcount,lcurpic
739           logical, save :: first_pass = .true.
740     !-----------------------------------------------
741     
742     ! locate the particles including ghost cells
743           lparcount = 1
744           lpic(:) = 0
745           if (plocate) then
746              dg_pijkprv(:)= dg_pijk(:)
747              do lcurpar = 1,max_pip
748                 if(lparcount.gt.pip) exit
749                 if(.not.pea(lcurpar,1)) cycle
750     
751     
752                 lparcount = lparcount + 1
753                 li = min(dg_iend2,max(dg_istart2,iofpos(des_pos_new(1,lcurpar))))
754                 lj = min(dg_jend2,max(dg_jstart2,jofpos(des_pos_new(2,lcurpar))))
755                 if(no_k) then
756                    lk = 1
757                 else
758                    lk = min(dg_kend2,max(dg_kstart2,kofpos(des_pos_new(3,lcurpar))))
759                 end if
760                 dg_pijk(lcurpar) = dg_funijk(li,lj,lk)
761                 lijk = dg_pijk(lcurpar)
762                 lpic(lijk) = lpic(lijk) + 1
763              end do
764           else
765              do lcurpar = 1,max_pip
766                 if(lparcount.gt.pip) exit
767                 if(.not.pea(lcurpar,1)) cycle
768                 lparcount = lparcount + 1
769                 lijk = dg_pijk(lcurpar)
770                 lpic(lijk) = lpic(lijk) + 1
771              end do
772           end if
773     
774           if (first_pass) then
775              dg_pijkprv(:) = dg_pijk(:)
776              first_pass = .false.
777           end if
778     
779     ! redfine the array of dg_pic
780     !!$omp parallel do default(shared)                               &
781     !!$omp private(lijk,lcurpic) schedule (guided,50)
782           do lijk = dg_ijkstart2,dg_ijkend2
783              lcurpic = lpic(lijk)
784              if(lcurpic > size(dg_pic(lijk)%p)) then
785                 deallocate(dg_pic(lijk)%p)
786                 allocate(dg_pic(lijk)%p(2*lcurpic))
787              end if
788              dg_pic(lijk)%isize = lcurpic
789           end do
790     !!$omp end parallel do
791     
792     
793     ! assign the particle info in pic array
794           lindx(:) = 1
795           lparcount = 1
796     
797     #if ( defined(__GFORTRAN__) &&  ( __GNUC__ < 4 || ( __GNUC__ == 4 && __GNUC_MINOR__ < 7 ) ) ) \
798           ||  ( defined(__INTEL_COMPILER) && (__INTEL_COMPILER < 1400) )
799     
800           do lcurpar = 1, max_pip
801              if(lparcount.gt.pip) exit
802              if(.not.pea(lcurpar,1)) cycle
803              lparcount = lparcount + 1
804              lijk = dg_pijk(lcurpar)
805              dg_pic(lijk)%p(lindx(lijk)) = lcurpar
806              lindx(lijk) = lindx(lijk) +  1
807           enddo
808     
809     #else
810     
811     !$omp parallel default(none) private(lcurpar,lijk,lijk_count) shared(max_pip,pip,lparcount,pea,dg_pijk,dg_pic,lindx)
812     !$omp do
813           do lcurpar = 1, max_pip
814              if(lparcount.gt.pip) cycle
815              if(.not.pea(lcurpar,1)) cycle
816              lijk = dg_pijk(lcurpar)
817     
818              !$omp atomic capture
819              lijk_count = lindx(lijk)
820              lindx(lijk) = lindx(lijk) +  1
821              !$omp end atomic
822     
823              dg_pic(lijk)%p(lijk_count) = lcurpar
824     
825              !$omp atomic
826              lparcount = lparcount + 1
827           enddo
828     !$omp end do
829     !$omp end parallel
830     
831     #endif
832     
833     !      open (unit=100,file='desgrid.txt',status='unknown')
834     !      write(100,*)lpic
835     !      close(100)
836     
837           end subroutine desgrid_pic
838     
839     !------------------------------------------------------------------------
840     ! subroutine       : desgrid_neigh_build ()
841     ! Purpose          : This particles build the neigh list for the particles
842     !                    currently active in the system
843     !------------------------------------------------------------------------
844           subroutine desgrid_neigh_build ()
845     
846     !-----------------------------------------------
847     ! Modules
848     !-----------------------------------------------
849           Use des_thermo
850           use compar
851           use constant
852           use des_allocate, only: add_pair
853           use desmpi_wrapper
854           use discretelement
855           use funits
856           use geometry
857           use param1
858     
859           implicit none
860     !-----------------------------------------------
861     ! Local variables
862     !-----------------------------------------------
863           integer lcurpar,lkoffset
864           integer lijk,lic,ljc,lkc,li,lj,lk,ltotpic,lpicloc,lneigh,lneighcnt
865           double precision lsearch_rad,ldistsquared
866           double precision :: ldistvec(3)
867           double precision :: lcurpar_pos(3)
868           double precision :: lcur_off
869           integer il_off,iu_off,jl_off,ju_off,kl_off,ku_off,mm,lSIZE2,lcurpar_index,lijk2
870     !$      INTEGER, DIMENSION(:,:), ALLOCATABLE :: int_tmp
871     !$    INTEGER :: PAIR_NUM_SMP,PAIR_MAX_SMP
872     !$    INTEGER, DIMENSION(:,:), ALLOCATABLE :: PAIRS_SMP
873     !$    integer num_threads
874     !-----------------------------------------------
875     
876     ! loop through neighbours and build the contact particles list for particles
877     ! present in the system
878           lkoffset = dimn-2
879     
880     !$omp parallel default(none) private(lcurpar,lijk,lic,ljc,lkc,lneighcnt,   &
881     !$omp    il_off,iu_off,jl_off,ju_off,kl_off,ku_off,lcurpar_pos,lcur_off,   &
882     !$omp    ltotpic, lneigh,lsearch_rad,ldistvec,ldistsquared, pair_num_smp, pair_max_smp, pairs_smp, lSIZE2, int_tmp) &
883     !$omp    shared(max_pip,pea,dg_pijk,NO_K,des_pos_new,dg_pic, factor_RLM,   &
884     !$omp           num_threads, des_radius, dg_xstart,dg_ystart,dg_zstart,dg_dxinv,dg_dyinv,dg_dzinv,dg_ijkstart2,dg_ijkend2)
885     
886     !$      PAIR_NUM_SMP = 0
887     !$      PAIR_MAX_SMP = 1024
888     !$      Allocate(  PAIRS_SMP(2,PAIR_MAX_SMP) )
889     
890     !$omp do
891     
892     
893           do lijk2 = dg_ijkstart2,dg_ijkend2
894     
895              do lcurpar_index = 1, dg_pic(lijk2)%isize
896                 lcurpar = dg_pic(lijk2)%p(lcurpar_index)
897     
898     !      do lcurpar =1,max_pip
899              if (.not. pea(lcurpar,1)) cycle
900              if (pea(lcurpar,2)) cycle
901              if (pea(lcurpar,4)) cycle
902              lneighcnt = 0
903     !         lijk = dg_pijk(lcurpar)
904              lic = dg_iof_lo(lijk2)
905              ljc = dg_jof_lo(lijk2)
906              lkc = dg_kof_lo(lijk2)
907     
908              il_off = 1
909              iu_off = 1
910              jl_off = 1
911              ju_off = 1
912              kl_off = 1
913              ku_off = 1
914     
915              lcurpar_pos(:) = des_pos_new(:,lcurpar)
916     !   The desgrid size should not be less than 2*dia*rlm_factor
917              lcur_off = (lcurpar_pos(1)-dg_xstart)*dg_dxinv - &
918                 floor((lcurpar_pos(1)-dg_xstart)*dg_dxinv)
919              if(lcur_off .ge. 0.5) then
920                 il_off = 0
921              else
922                 iu_off = 0
923              endif
924     
925              lcur_off = (lcurpar_pos(2)-dg_ystart)*dg_dyinv - &
926                 floor((lcurpar_pos(2)-dg_ystart)*dg_dyinv)
927              if(lcur_off .ge. 0.5) then
928                 jl_off = 0
929              else
930                 ju_off = 0
931              endif
932     
933              if(NO_K)then
934                 kl_off = 0
935                 ku_off = 0
936              else
937                lcur_off = (lcurpar_pos(3)-dg_zstart)*dg_dzinv - &
938                   floor((lcurpar_pos(3)-dg_zstart)*dg_dzinv)
939                if(lcur_off .ge. 0.5) then
940                   kl_off = 0
941                else
942                   ku_off = 0
943                endif
944             endif
945     
946              do lk = lkc-kl_off,lkc+ku_off
947              do lj = ljc-jl_off,ljc+ju_off
948              do li = lic-il_off,lic+iu_off
949                 lijk = dg_funijk(li,lj,lk)
950                 ltotpic =dg_pic(lijk)%isize
951                 do lpicloc = 1,ltotpic
952                    lneigh = dg_pic(lijk)%p(lpicloc)
953     
954     ! Only skip real particles otherwise collisions with ghost, entering,
955     ! and exiting particles are missed.
956                    if (lneigh <= lcurpar) then
957                       if(.not.any(pea(lneigh,2:4))) cycle
958                    endif
959     
960                    lsearch_rad = factor_RLM*(des_radius(lcurpar)+des_radius(lneigh))
961                    ldistvec = lcurpar_pos(:)-des_pos_new(:,lneigh)
962                    ldistsquared = dot_product(ldistvec,ldistvec)
963                    if (ldistsquared.gt.lsearch_rad*lsearch_rad) cycle
964                    lneighcnt = lneighcnt + 1
965                    if (pea(lcurpar,1) .and. .not.pea(lcurpar,4) .and. pea(lneigh,1)) THEN
966     !$  if (.true.) then
967     !!!!! SMP version
968     
969     !$      pair_num_smp = pair_num_smp + 1
970     ! Reallocate to double the size of the arrays if needed.
971     !$      IF(PAIR_NUM_SMP > PAIR_MAX_SMP) THEN
972     !$         PAIR_MAX_SMP = 2*PAIR_MAX_SMP
973     !$         lSIZE2 = size(pairs_smp,2)
974     
975     !$         allocate(int_tmp(2,PAIR_MAX_SMP))
976     !$         int_tmp(:,1:lSIZE2) = pairs_smp(:,1:lSIZE2)
977     !$         call move_alloc(int_tmp,pairs_smp)
978     !$      ENDIF
979     
980     !$      pairs_smp(1,pair_num_smp) = lcurpar
981     !$      pairs_smp(2,pair_num_smp) = lneigh
982     
983     !$  else
984     !!!!! Serial version
985                       call add_pair(lcurpar, lneigh)
986     !$  endif
987                    endif
988                 end do
989              end do
990              end do
991              end do
992           end do
993           end do
994     !$omp end do
995     
996     !$omp critical
997     !$     do MM = 1, PAIR_NUM_SMP
998     !$        call add_pair(PAIRS_SMP(1,MM), PAIRS_SMP(2,MM))
999     !$     enddo
1000     !$omp end critical
1001     
1002     !$    deallocate( PAIRS_SMP )
1003     
1004     !$omp end parallel
1005     
1006           end subroutine desgrid_neigh_build
1007     
1008     !------------------------------------------------------------------------
1009     ! subroutine       : des_dbggrid
1010     ! Purpose          : For printing the indices used for desgrid
1011     !------------------------------------------------------------------------
1012           subroutine des_dbggrid()
1013     
1014           use param1
1015           use funits
1016           use geometry
1017           use compar
1018           use discretelement
1019           use constant
1020           use desmpi_wrapper
1021     
1022           implicit none
1023     !-----------------------------------------------
1024     ! Local variables
1025     !-----------------------------------------------
1026           integer lproc,liproc,ljproc,lkproc
1027           character (30) filename
1028     !-----------------------------------------------
1029     
1030           write(filename,'("dbg_desgridn",I4.4,".dat")') mype
1031           open(44,file=filename)
1032           do lproc = 0,numpes-1
1033              write(44,*) "Information for Proc =", lproc
1034              liproc= iofproc(lproc)
1035              ljproc= jofproc(lproc)
1036              lkproc= kofproc(lproc)
1037              write(44,*) "   "
1038              write(44,*) "i,j,k location of proc",liproc,ljproc,lkproc
1039              write(44,*) "i,j,k size of proc", dg_isize_all(liproc),dg_jsize_all(ljproc),dg_ksize_all(lkproc)
1040              write(44,*) "-------------------------------------------------"
1041              write(44,*) "indices        start     end"
1042              write(44,*) "-------------------------------------------------"
1043              write(44,*) "for i1:      ",dg_istart1_all(lproc),dg_iend1_all(lproc)
1044              write(44,*) "for i2:      ",dg_istart2_all(lproc),dg_iend2_all(lproc)
1045              write(44,*) "for j1:      ",dg_jstart1_all(lproc),dg_jend1_all(lproc)
1046              write(44,*) "for j2:      ",dg_jstart2_all(lproc),dg_jend2_all(lproc)
1047              write(44,*) "for k1:      ",dg_kstart1_all(lproc),dg_kend1_all(lproc)
1048              write(44,*) "for k2:      ",dg_kstart2_all(lproc),dg_kend2_all(lproc)
1049           end do
1050           write(44,*) "   "
1051           write(44,*) "-------------------------------------------------"
1052           write(44,*) "Local Start and end"
1053           write(44,*) "-------------------------------------------------"
1054           write(44,*) "for i :      ",dg_istart, dg_iend
1055           write(44,*) "for i1:      ",dg_istart1,dg_iend1
1056           write(44,*) "for i2:      ",dg_istart2,dg_iend2
1057           write(44,*) "   "
1058           write(44,*) "for j :      ",dg_jstart, dg_jend
1059           write(44,*) "for j1:      ",dg_jstart1,dg_jend1
1060           write(44,*) "for j2:      ",dg_jstart2,dg_jend2
1061           write(44,*) "   "
1062           write(44,*) "for k :      ",dg_kstart, dg_kend
1063           write(44,*) "for k1:      ",dg_kstart1,dg_kend1
1064           write(44,*) "for k2:      ",dg_kstart2,dg_kend2
1065           write(44,*) "   "
1066           write(44,*) "-------------------------------------------------"
1067           write(44,*) "global Start and end"
1068           write(44,*) "-------------------------------------------------"
1069           write(44,*) "for i1:      ",dg_imin1,dg_imax1
1070           write(44,*) "for i2:      ",dg_imin2,dg_imax2
1071           write(44,*) "for j1:      ",dg_jmin1,dg_jmax1
1072           write(44,*) "for j2:      ",dg_jmin2,dg_jmax2
1073           write(44,*) "for k1:      ",dg_kmin1,dg_kmax1
1074           write(44,*) "for k2:      ",dg_kmin2,dg_kmax2
1075           write(44,*) "   "
1076           write(44,*) "-------------------------------------------------"
1077           write(44,*) "dg_xstart:   ",dg_xstart
1078           write(44,*) "dg_xend:     ",dg_xend
1079           write(44,*) "dg_dxinv:    ",dg_dxinv
1080           write(44,*) "   "
1081           write(44,*) "dg_ystart:   ",dg_ystart
1082           write(44,*) "dg_yend:     ",dg_yend
1083           write(44,*) "dg_dyinv:    ",dg_dyinv
1084           write(44,*) "   "
1085           write(44,*) "dg_zstart:   ",dg_zstart
1086           write(44,*) "dg_zend:     ",dg_zend
1087           write(44,*) "dg_dzinv:    ",dg_dzinv
1088           write(44,*) "   "
1089           write(44,*) "dg_c1_lo/gl: ",dg_c1_lo, dg_c1_gl
1090           write(44,*) "dg_c2_lo/gl: ",dg_c2_lo, dg_c2_gl
1091           write(44,*) "dg_c3_lo/gl: ",dg_c3_lo, dg_c3_gl
1092           write(44,*) "   "
1093           write(44,*) "-------------------------------------------------"
1094     
1095     
1096     
1097           close(44)
1098           end subroutine des_dbggrid
1099     
1100           end module
1101