source: trunk/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES/modcluster.F90 @ 5656

Last change on this file since 5656 was 5656, checked in by timgraham, 5 years ago

Merge of AGRIF branch (branches/2014/dev_r4765_CNRS_agrif) onto the trunk

  • Property svn:keywords set to Id
File size: 49.4 KB
Line 
1!
2! $Id$
3!
4!     AGRIF (Adaptive Grid Refinement In Fortran)
5!
6!     Copyright (C) 2003 Laurent Debreu (Laurent.Debreu@imag.fr)
7!                        Christophe Vouland (Christophe.Vouland@imag.fr)
8!
9!     This program is free software; you can redistribute it and/or modify
10!     it under the terms of the GNU General Public License as published by
11!     the Free Software Foundation; either version 2 of the License, or
12!     (at your option) any later version.
13!
14!     This program is distributed in the hope that it will be useful,
15!     but WITHOUT ANY WARRANTY; without even the implied warranty of
16!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17!     GNU General Public License for more details.
18!
19!     You should have received a copy of the GNU General Public License
20!     along with this program; if not, write to the Free Software
21!     Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
22!
23!
24!> Module Agrif_Clustering
25!>
26!> Contains subroutines to create and initialize the grid hierarchy from the
27!> AGRIF_FixedGrids.in file.
28!
29module Agrif_Clustering
30!
31    use Agrif_CurgridFunctions
32    use Agrif_Init_Vars
33    use Agrif_Save
34!
35    implicit none
36!
37    abstract interface
38        subroutine init_proc()
39        end subroutine init_proc
40    end interface
41!
42contains
43!
44!===================================================================================================
45!  subroutine Agrif_Cluster_All
46!
47!> Subroutine for the clustering. A temporary grid hierarchy, pointed by parent_rect, is created.
48!---------------------------------------------------------------------------------------------------
49recursive subroutine Agrif_Cluster_All ( g, parent_rect )
50!---------------------------------------------------------------------------------------------------
51    TYPE(Agrif_Grid)     , pointer   :: g        !< Pointer on the current grid
52    TYPE(Agrif_Rectangle), pointer   :: parent_rect
53!
54    TYPE(Agrif_LRectangle), pointer  :: parcours
55    TYPE(Agrif_Grid)      , pointer  :: newgrid
56    REAL                             :: g_eps
57    INTEGER                          :: i
58!
59    g_eps = huge(1.)
60    do i = 1,Agrif_Probdim
61        g_eps = min(g_eps, g % Agrif_dx(i))
62    enddo
63!
64    g_eps = g_eps / 100.
65!
66!   Necessary condition for clustering
67    do i = 1,Agrif_Probdim
68        if ( g%Agrif_dx(i)/Agrif_coeffref(i) < (Agrif_mind(i)-g_eps)) return
69    enddo
70!
71    nullify(parent_rect%childgrids)
72!
73    call Agrif_ClusterGridnD(g,parent_rect)
74!
75    parcours => parent_rect % childgrids
76!
77    do while ( associated(parcours) )
78!
79!       Newgrid is created. It is a copy of a fine grid created previously by clustering.
80        allocate(newgrid)
81!
82        do i = 1,Agrif_Probdim
83            newgrid % nb(i) = (parcours % r % imax(i) - parcours % r % imin(i)) * Agrif_Coeffref(i)
84            newgrid % Agrif_x(i)  = g % Agrif_x(i)  + (parcours % r % imin(i) -1) * g%Agrif_dx(i)
85            newgrid % Agrif_dx(i) = g % Agrif_dx(i) / Agrif_Coeffref(i)
86        enddo
87!
88        if ( Agrif_Probdim == 1 ) then
89            allocate(newgrid%tabpoint1D(newgrid%nb(1)+1))
90            newgrid%tabpoint1D = 0
91        endif
92!
93        if ( Agrif_Probdim == 2 ) then
94            allocate(newgrid%tabpoint2D(newgrid%nb(1)+1, newgrid%nb(2)+1))
95            newgrid%tabpoint2D = 0
96        endif
97!
98        if ( Agrif_Probdim == 3 ) then
99            allocate(newgrid%tabpoint3D(newgrid%nb(1)+1, newgrid%nb(2)+1, newgrid%nb(3)+1))
100            newgrid%tabpoint3D = 0
101        endif
102!
103!       Points detection on newgrid
104        call Agrif_TabpointsnD(Agrif_mygrid,newgrid)
105!
106!       recursive call to Agrif_Cluster_All
107        call Agrif_Cluster_All(newgrid, parcours % r)
108!
109        parcours => parcours % next
110!
111        if ( Agrif_Probdim == 1 ) deallocate(newgrid%tabpoint1D)
112        if ( Agrif_Probdim == 2 ) deallocate(newgrid%tabpoint2D)
113        if ( Agrif_Probdim == 3 ) deallocate(newgrid%tabpoint3D)
114!
115        deallocate(newgrid)
116!
117    enddo
118!---------------------------------------------------------------------------------------------------
119end subroutine Agrif_Cluster_All
120!===================================================================================================
121!
122!===================================================================================================
123!  subroutine Agrif_TabpointsnD
124!
125!> Copy on newgrid of points detected on the grid hierarchy pointed by g.
126!---------------------------------------------------------------------------------------------------
127recursive subroutine Agrif_TabpointsnD ( g, newgrid )
128!---------------------------------------------------------------------------------------------------
129    TYPE(Agrif_Grid), pointer   :: g, newgrid
130!
131    TYPE(Agrif_PGrid), pointer  :: parcours
132!
133    REAL                  :: g_eps, newgrid_eps, eps
134    REAL   , DIMENSION(3) :: newmin, newmax
135    REAL   , DIMENSION(3) :: gmin, gmax
136    REAL   , DIMENSION(3) :: xmin
137    INTEGER, DIMENSION(3) :: igmin, inewmin
138    INTEGER, DIMENSION(3) :: inewmax
139    INTEGER               :: i,  j,  k
140    INTEGER               :: i0, j0, k0
141!
142    parcours => g % child_list % first
143!
144    do while( associated(parcours) )
145        call Agrif_TabpointsnD(parcours%gr,newgrid)
146        parcours => parcours % next
147    enddo
148!
149    g_eps = 10.
150    newgrid_eps = 10.
151!
152    do i = 1,Agrif_probdim
153        g_eps = min( g_eps , g % Agrif_dx(i) )
154        newgrid_eps = min(newgrid_eps,newgrid%Agrif_dx(i))
155    enddo
156!
157    eps = min(g_eps,newgrid_eps)/100.
158!
159    do i = 1,Agrif_probdim
160!
161         if ( newgrid%Agrif_dx(i) < (g%Agrif_dx(i)-eps) ) return
162!
163         gmin(i) = g%Agrif_x(i)
164         gmax(i) = g%Agrif_x(i) + g%nb(i) * g%Agrif_dx(i)
165!
166         newmin(i) = newgrid%Agrif_x(i)
167         newmax(i) = newgrid%Agrif_x(i) + newgrid%nb(i) * newgrid%Agrif_dx(i)
168!
169         if (gmax(i) < newmin(i)) return
170         if (gmin(i) > newmax(i)) return
171!
172         inewmin(i) = 1 - floor(-(max(gmin(i),newmin(i))-newmin(i)) / newgrid%Agrif_dx(i))
173!
174         xmin(i) = newgrid%Agrif_x(i) + (inewmin(i)-1)*newgrid%Agrif_dx(i)
175!
176         igmin(i) = 1 + nint((xmin(i)-gmin(i))/g%Agrif_dx(i))
177!
178         inewmax(i) = 1 + int(   (min(gmax(i),newmax(i))-newmin(i)) / newgrid%Agrif_dx(i))
179!
180    enddo
181!
182    if ( Agrif_probdim == 1 ) then
183        i0 = igmin(1)
184        do i = inewmin(1),inewmax(1)
185            newgrid%tabpoint1D(i) = max( newgrid%tabpoint1D(i), g%tabpoint1D(i0) )
186        enddo
187        i0 = i0 + int(newgrid%Agrif_dx(1)/g%Agrif_dx(1))
188    endif
189!
190    if ( Agrif_probdim == 2 ) then
191        i0 = igmin(1)
192        do i = inewmin(1),inewmax(1)
193            j0 = igmin(2)
194            do j = inewmin(2),inewmax(2)
195                newgrid%tabpoint2D(i,j) = max( newgrid%tabpoint2D(i,j), g%tabpoint2D(i0,j0) )
196                j0 = j0 + int(newgrid%Agrif_dx(2)/g%Agrif_dx(2))
197            enddo
198            i0 = i0 + int(newgrid%Agrif_dx(1)/g%Agrif_dx(1))
199        enddo
200    endif
201!
202    if ( Agrif_probdim == 3 ) then
203        i0 = igmin(1)
204        do i = inewmin(1),inewmax(1)
205            j0 = igmin(2)
206            do j = inewmin(2),inewmax(2)
207                k0 = igmin(3)
208                do k = inewmin(3),inewmax(3)
209                    newgrid%tabpoint3D(i,j,k) = max( newgrid%tabpoint3D(i,j,k), g%tabpoint3D(i0,j0,k0))
210                    k0 = k0 + int(newgrid%Agrif_dx(3)/g%Agrif_dx(3))
211                enddo
212                j0 = j0 + int(newgrid%Agrif_dx(2)/g%Agrif_dx(2))
213            enddo
214            i0 = i0 + int(newgrid%Agrif_dx(1)/g%Agrif_dx(1))
215        enddo
216    endif
217!---------------------------------------------------------------------------------------------------
218end subroutine Agrif_TabpointsnD
219!===================================================================================================
220!
221!===================================================================================================
222!  subroutine Agrif_ClusterGridnD
223!
224!> Clustering on the grid pointed by g.
225!---------------------------------------------------------------------------------------------------
226subroutine Agrif_ClusterGridnD ( g, parent_rect )
227!---------------------------------------------------------------------------------------------------
228    TYPE(Agrif_Grid)     , pointer  :: g          !< Pointer on the current grid
229    TYPE(Agrif_Rectangle), pointer  :: parent_rect
230!
231    TYPE(Agrif_Rectangle) :: newrect
232    TYPE(Agrif_Variable_i)  :: newflag
233!
234    INTEGER               :: i,j,k
235    INTEGER, DIMENSION(3) :: sx
236    INTEGER               :: bufferwidth,flagpoints
237    INTEGER               :: n1,n2,m1,m2,o1,o2
238!
239    bufferwidth = int(Agrif_Minwidth/5.)
240!
241    do i = 1,Agrif_probdim
242        sx(i) = g % nb(i) + 1
243    enddo
244!
245    if ( Agrif_probdim == 1 ) then
246        allocate(newflag%iarray1(sx(1)))
247        newflag%iarray1 = 0
248    endif
249    if ( Agrif_probdim == 2 ) then
250        allocate(newflag%iarray2(sx(1),sx(2)))
251        newflag%iarray2 = 0
252    endif
253    if ( Agrif_probdim == 3 ) then
254        allocate(newflag%iarray3(sx(1),sx(2),sx(3)))
255        newflag%iarray3 = 0
256    endif
257!
258    flagpoints = 0
259!
260    if ( bufferwidth>0 ) then
261!
262        if ( Agrif_probdim == 1 ) then
263            do i = bufferwidth,sx(1)-bufferwidth+1
264                if (g % tabpoint1D(i) == 1) then
265                    m1 = i - bufferwidth + 1
266                    m2 = i + bufferwidth - 1
267                    flagpoints = flagpoints + 1
268                    newflag%iarray1(m1:m2) = 1
269                endif
270            enddo
271        endif
272!
273        if ( Agrif_probdim == 2 ) then
274            do i = bufferwidth,sx(1)-bufferwidth+1
275            do j = bufferwidth,sx(2)-bufferwidth+1
276                if (g % tabpoint2D(i,j) == 1) then
277                    n1 = j - bufferwidth + 1
278                    n2 = j + bufferwidth - 1
279                    m1 = i - bufferwidth + 1
280                    m2 = i + bufferwidth - 1
281                    flagpoints = flagpoints + 1
282                    newflag%iarray2(m1:m2,n1:n2) = 1
283                endif
284            enddo
285            enddo
286        endif
287!
288        if ( Agrif_probdim == 3 ) then
289            do i = bufferwidth,sx(1)-bufferwidth+1
290            do j = bufferwidth,sx(2)-bufferwidth+1
291            do k = bufferwidth,sx(3)-bufferwidth+1
292                if (g % tabpoint3D(i,j,k) == 1) then
293                    o1 = k - bufferwidth + 1
294                    o2 = k + bufferwidth - 1
295                    n1 = j - bufferwidth + 1
296                    n2 = j + bufferwidth - 1
297                    m1 = i - bufferwidth + 1
298                    m2 = i + bufferwidth - 1
299                    flagpoints = flagpoints + 1
300                    newflag%iarray3(m1:m2,n1:n2,o1:o2) = 1
301                endif
302            enddo
303            enddo
304            enddo
305        endif
306!
307    else
308!
309        flagpoints = 1
310        if ( Agrif_probdim == 1 ) newflag%iarray1 = g % tabpoint1D
311        if ( Agrif_probdim == 2 ) newflag%iarray2 = g % tabpoint2D
312        if ( Agrif_probdim == 3 ) newflag%iarray3 = g % tabpoint3D
313!
314    endif
315!
316    if (flagpoints == 0) then
317        if ( Agrif_probdim == 1 ) deallocate(newflag%iarray1)
318        if ( Agrif_probdim == 2 ) deallocate(newflag%iarray2)
319        if ( Agrif_probdim == 3 ) deallocate(newflag%iarray3)
320        return
321    endif
322!
323    do i = 1 , Agrif_probdim
324        newrect % imin(i) = 1
325        newrect % imax(i) = sx(i)
326    enddo
327!
328    call Agrif_Clusternd(newflag,parent_rect%childgrids,newrect)
329!
330    if ( Agrif_probdim == 1 ) deallocate(newflag%iarray1)
331    if ( Agrif_probdim == 2 ) deallocate(newflag%iarray2)
332    if ( Agrif_probdim == 3 ) deallocate(newflag%iarray3)
333!---------------------------------------------------------------------------------------------------
334end subroutine Agrif_ClusterGridnD
335!===================================================================================================
336!
337!===================================================================================================
338!  subroutine Agrif_ClusternD
339!
340!> Clustering on the grid pointed by oldB.
341!---------------------------------------------------------------------------------------------------
342recursive subroutine Agrif_Clusternd ( flag, boxlib, oldB )
343!---------------------------------------------------------------------------------------------------
344    TYPE(Agrif_Variable_i)            :: flag
345    TYPE(Agrif_LRectangle), pointer :: boxlib
346    TYPE(Agrif_Rectangle)           :: oldB
347!
348    TYPE(Agrif_LRectangle),pointer :: tempbox,parcbox,parcbox2
349    TYPE(Agrif_Rectangle) :: newB,newB2
350    INTEGER :: i,j,k
351    INTEGER, DIMENSION(:), allocatable :: i_sig, j_sig, k_sig
352    INTEGER, DIMENSION(3) :: ipu,ipl
353    INTEGER, DIMENSION(3) :: istr,islice
354    REAL :: cureff, neweff
355    INTEGER :: ValMax,ValSum,TailleTab
356    INTEGER :: nbpoints,nbpointsflag
357    LOGICAL :: test
358!
359                              allocate( i_sig(oldB%imin(1):oldB%imax(1)) )
360    if ( Agrif_probdim >= 2 ) allocate( j_sig(oldB%imin(2):oldB%imax(2)) )
361    if ( Agrif_probdim == 3 ) allocate( k_sig(oldB%imin(3):oldB%imax(3)) )
362!
363    test = .FALSE.
364    do i = 1,Agrif_probdim
365        test = test .OR. ( (oldB%imax(i)-oldB%imin(i)+1) <  Agrif_Minwidth)
366    enddo
367    if ( test ) return
368!
369    if ( Agrif_probdim == 1 ) i_sig = flag%iarray1
370    if ( Agrif_probdim == 2 ) then
371        do i = oldB%imin(1),oldB%imax(1)
372            i_sig(i) = SUM(flag%iarray2(i, oldB%imin(2):oldB%imax(2)))
373        enddo
374        do j = oldB%imin(2),oldB%imax(2)
375            j_sig(j) = SUM(flag%iarray2(oldB%imin(1):oldB%imax(1),j))
376        enddo
377    endif
378    if ( Agrif_probdim == 3 ) then
379        do i = oldB%imin(1),oldB%imax(1)
380            i_sig(i) = SUM(flag%iarray3(i,oldB%imin(2):oldB%imax(2),    &
381                                          oldB%imin(3):oldB%imax(3)))
382        enddo
383        do j = oldB%imin(2),oldB%imax(2)
384            j_sig(j) = SUM(flag%iarray3(  oldB%imin(1):oldB%imax(1), j, &
385                                          oldB%imin(3):oldB%imax(3)))
386        enddo
387        do k = oldB%imin(3),oldB%imax(3)
388            k_sig(k) = SUM(flag%iarray3(  oldB%imin(1):oldB%imax(1),    &
389                                          oldB%imin(2):oldB%imax(2), k) )
390        enddo
391    endif
392!
393    do i = 1,Agrif_probdim
394        ipl(i) = oldB%imin(i)
395        ipu(i) = oldB%imax(i)
396    enddo
397!
398                              call Agrif_Clusterprune(i_sig,ipl(1),ipu(1))
399    if ( Agrif_probdim >= 2 ) call Agrif_Clusterprune(j_sig,ipl(2),ipu(2))
400    if ( Agrif_probdim == 3 ) call Agrif_Clusterprune(k_sig,ipl(3),ipu(3))
401!
402    test = .TRUE.
403    do i = 1,Agrif_probdim
404        test = test .AND. (ipl(i) == oldB%imin(i))
405        test = test .AND. (ipu(i) == oldB%imax(i))
406    enddo
407
408    if (.NOT. test) then
409        do i = 1 , Agrif_probdim
410            newB%imin(i) = ipl(i)
411            newB%imax(i) = ipu(i)
412        enddo
413!
414        if ( Agrif_probdim == 1 ) nbpoints = SUM(flag%iarray1(newB%imin(1):newB%imax(1)))
415        if ( Agrif_probdim == 2 ) nbpoints = SUM(flag%iarray2(newB%imin(1):newB%imax(1),    &
416                                                              newB%imin(2):newB%imax(2)))
417        if ( Agrif_probdim == 3 ) nbpoints = SUM(flag%iarray3(newB%imin(1):newB%imax(1),    &
418                                                              newB%imin(2):newB%imax(2),    &
419                                                              newB%imin(3):newB%imax(3)))
420!
421        if ( Agrif_probdim == 1 ) TailleTab = (newB%imax(1)-newB%imin(1)+1)
422        if ( Agrif_probdim == 2 ) TailleTab = (newB%imax(1)-newB%imin(1)+1) * &
423                                              (newB%imax(2)-newB%imin(2)+1)
424        if ( Agrif_probdim == 3 ) TailleTab = (newB%imax(1)-newB%imin(1)+1) * &
425                                              (newB%imax(2)-newB%imin(2)+1) * &
426                                              (newB%imax(3)-newB%imin(3)+1)
427        neweff = REAL(nbpoints) / TailleTab
428!
429        if (nbpoints > 0) then
430!
431            if ((neweff > Agrif_efficiency)) then
432                call Agrif_Add_Rectangle(newB,boxlib)
433                return
434            else
435!
436                tempbox => boxlib
437                newB2 = newB
438                call Agrif_Clusternd(flag,boxlib,newB)
439!
440!               Compute new efficiency
441                cureff = neweff
442                parcbox2 => boxlib
443                nbpoints = 0
444                nbpointsflag = 0
445!
446                do while (associated(parcbox2))
447                    if (associated(parcbox2,tempbox)) exit
448                    newB = parcbox2%r
449!
450                    if ( Agrif_probdim == 1 ) Valsum = SUM(flag%iarray1(newB%imin(1):newB%imax(1)))
451                    if ( Agrif_probdim == 2 ) Valsum = SUM(flag%iarray2(newB%imin(1):newB%imax(1), &
452                                                                        newB%imin(2):newB%imax(2)))
453                    if ( Agrif_probdim == 3 ) Valsum = SUM(flag%iarray3(newB%imin(1):newB%imax(1), &
454                                                                        newB%imin(2):newB%imax(2), &
455                                                                        newB%imin(3):newB%imax(3)))
456                    nbpointsflag = nbpointsflag + ValSum
457!
458                    if ( Agrif_probdim == 1 ) TailleTab = (newB%imax(1)-newB%imin(1)+1)
459                    if ( Agrif_probdim == 2 ) TailleTab = (newB%imax(1)-newB%imin(1)+1) * &
460                                                          (newB%imax(2)-newB%imin(2)+1)
461                    if ( Agrif_probdim == 3 ) TailleTab = (newB%imax(1)-newB%imin(1)+1) * &
462                                                          (newB%imax(2)-newB%imin(2)+1) * &
463                                                          (newB%imax(3)-newB%imin(3)+1)
464                    nbpoints = nbpoints + TailleTab
465                    parcbox2 => parcbox2%next
466                enddo
467!
468! coefficient 1.05 avant 1.15 possibilite de laisser choix a l utilisateur
469                if ( REAL(nbpointsflag)/REAL(nbpoints) < (1.0001*cureff)) then
470                    parcbox2 => boxlib
471                    do while (associated(parcbox2))
472                        if (associated(parcbox2,tempbox)) exit
473                        deallocate(parcbox2%r)
474                        parcbox2 => parcbox2%next
475                    enddo
476                    boxlib => tempbox
477                    call Agrif_Add_Rectangle(newB2,boxlib)
478                    return
479                endif
480            endif
481        endif
482        return
483    endif
484!
485    do i = 1,Agrif_Probdim
486        istr(i) = oldB%imax(i)
487        islice(i) = oldB%imin(i)
488    enddo
489!
490                              call Agrif_Clusterslice(i_sig,islice(1),istr(1))
491    if ( Agrif_probdim >= 2 ) call Agrif_Clusterslice(j_sig,islice(2),istr(2))
492    if ( Agrif_probdim == 3 ) call Agrif_Clusterslice(k_sig,islice(3),istr(3))
493!
494    ValSum = 0
495    do i = 1,Agrif_Probdim
496        Valsum = valSum + islice(i)
497    enddo
498!
499    if ( Valsum == -Agrif_Probdim ) then
500        call Agrif_Add_Rectangle(oldB,boxlib)
501        return
502    endif
503!
504    nullify(tempbox)
505    tempbox => boxlib
506    if ( Agrif_probdim == 1 ) cureff = (oldB%imax(1)-oldB%imin(1)+1)
507    if ( Agrif_probdim == 2 ) cureff = (oldB%imax(1)-oldB%imin(1)+1) * &
508                                       (oldB%imax(2)-oldB%imin(2)+1)
509    if ( Agrif_probdim == 3 ) cureff = (oldB%imax(1)-oldB%imin(1)+1) * &
510                                       (oldB%imax(2)-oldB%imin(2)+1) * &
511                                       (oldB%imax(3)-oldB%imin(3)+1)
512    nullify(parcbox)
513!
514    do i = 1,Agrif_Probdim
515        newB%imax(i) = oldB%imax(i)
516        newB%imin(i) = oldB%imin(i)
517    enddo
518!
519    ValMax = 0
520    do i = 1 , Agrif_Probdim
521        ValMax = Max(ValMax,istr(i))
522    enddo
523!
524    if (istr(1) == ValMax ) then
525        newB%imax(1) = islice(1)
526        call Agrif_Add_Rectangle(newB,parcbox)
527        newB%imin(1) = islice(1)+1
528        newB%imax(1) = oldB%imax(1)
529        call Agrif_Add_Rectangle(newB,parcbox)
530    else if ( Agrif_probdim >= 2 ) then
531        if (istr(2) == ValMax ) then
532            newB%imax(2) = islice(2)
533            call Agrif_Add_Rectangle(newB,parcbox)
534            newB%imin(2) = islice(2)+1
535            newB%imax(2) = oldB%imax(2)
536            call Agrif_Add_Rectangle(newB,parcbox)
537        else if ( Agrif_probdim == 3 ) then
538            newB%imax(3) = islice(3)
539            call Agrif_Add_Rectangle(newB,parcbox)
540            newB%imin(3) = islice(3)+1
541            newB%imax(3) = oldB%imax(3)
542            call Agrif_Add_Rectangle(newB,parcbox)
543        endif
544    endif
545!
546    do while ( associated(parcbox) )
547        newB = parcbox%r
548!
549        if ( Agrif_probdim == 1 ) nbpoints = SUM(flag%iarray1(newB%imin(1):newB%imax(1)))
550        if ( Agrif_probdim == 2 ) nbpoints = SUM(flag%iarray2(newB%imin(1):newB%imax(1),    &
551                                                              newB%imin(2):newB%imax(2)))
552        if ( Agrif_probdim == 3 ) nbpoints = SUM(flag%iarray3(newB%imin(1):newB%imax(1),    &
553                                                              newB%imin(2):newB%imax(2),    &
554                                                              newB%imin(3):newB%imax(3)))
555!
556        if ( Agrif_probdim == 1 ) TailleTab = (newB%imax(1)-newB%imin(1)+1)
557        if ( Agrif_probdim == 2 ) TailleTab = (newB%imax(1)-newB%imin(1)+1) * &
558                                              (newB%imax(2)-newB%imin(2)+1)
559        if ( Agrif_probdim == 3 ) TailleTab = (newB%imax(1)-newB%imin(1)+1) * &
560                                              (newB%imax(2)-newB%imin(2)+1) * &
561                                              (newB%imax(3)-newB%imin(3)+1)
562        neweff = REAL(nbpoints) / TailleTab
563!
564        if ( nbpoints > 0 ) then
565!
566            if ( (neweff > Agrif_efficiency)) then
567                call Agrif_Add_Rectangle(newB,boxlib)
568            else
569                tempbox => boxlib
570                newB2 = newB
571                call Agrif_Clusternd(flag,boxlib,newB)
572!
573!               Compute new efficiency
574                cureff = neweff
575                parcbox2 => boxlib
576                nbpoints = 0
577                nbpointsflag = 0
578!
579                do while (associated(parcbox2))
580                    if (associated(parcbox2,tempbox)) exit
581                    newB = parcbox2%r
582!
583                    if ( Agrif_probdim == 1 ) ValSum = SUM(flag%iarray1(newB%imin(1):newB%imax(1)))
584                    if ( Agrif_probdim == 2 ) ValSum = SUM(flag%iarray2(newB%imin(1):newB%imax(1),  &
585                                                                        newB%imin(2):newB%imax(2)))
586                    if ( Agrif_probdim == 3 ) ValSum = SUM(flag%iarray3(newB%imin(1):newB%imax(1),  &
587                                                                        newB%imin(2):newB%imax(2),  &
588                                                                        newB%imin(3):newB%imax(3)))
589                    nbpointsflag = nbpointsflag + ValSum
590!
591                    if ( Agrif_probdim == 1 ) TailleTab = (newB%imax(1)-newB%imin(1)+1)
592                    if ( Agrif_probdim == 2 ) TailleTab = (newB%imax(1)-newB%imin(1)+1) * &
593                                                          (newB%imax(2)-newB%imin(2)+1)
594                    if ( Agrif_probdim == 3 ) TailleTab = (newB%imax(1)-newB%imin(1)+1) * &
595                                                          (newB%imax(2)-newB%imin(2)+1) * &
596                                                          (newB%imax(3)-newB%imin(3)+1)
597                    nbpoints = nbpoints + TailleTab
598                    parcbox2 => parcbox2%next
599                enddo
600!
601                if  ( REAL(nbpointsflag)/REAL(nbpoints) < (1.15*cureff)) then
602                    parcbox2 => boxlib
603                    do while (associated(parcbox2))
604                        if (associated(parcbox2,tempbox)) exit
605                        deallocate(parcbox2%r)
606                        parcbox2 => parcbox2%next
607                    enddo
608                    boxlib => tempbox
609                    call Agrif_Add_Rectangle(newB2,boxlib)
610                endif
611            endif
612        endif
613        parcbox => parcbox%next
614    enddo
615!---------------------------------------------------------------------------------------------------
616end subroutine Agrif_Clusternd
617!===================================================================================================
618!
619!===================================================================================================
620!  subroutine Agrif_Clusterslice
621!---------------------------------------------------------------------------------------------------
622subroutine Agrif_Clusterslice ( sig, slice, str )
623!---------------------------------------------------------------------------------------------------
624    INTEGER,                       intent(inout) :: slice
625    INTEGER,                       intent(inout) :: str
626    INTEGER, DIMENSION(slice:str), intent(in)    :: sig
627!
628    INTEGER                         :: ideb, ifin
629    INTEGER                         :: i, t, a, di, ts
630    INTEGER, DIMENSION(slice:str)   :: lap
631!
632    ideb = slice
633    ifin = str
634!
635    if (SIZE(sig) <= 2*Agrif_Minwidth) then
636        str = -1
637        slice = -1
638        return
639    endif
640!
641    t = -1
642    a = -1
643!
644    do i = ideb+Agrif_Minwidth-1,ifin-Agrif_Minwidth
645        if (sig(i) == 0) then
646            if ((i-ideb) < (ifin-i)) then
647                di = i - ideb
648            else
649                di = ifin - i
650            endif
651!
652            if (di > t) then
653                a = i
654                t = di
655            endif
656        endif
657    enddo
658!
659    if (a /= -1) then
660        slice = a
661        str = t
662        return
663    endif
664!
665    t = -1
666    a = -1
667!
668    do i = ideb+1,ifin-1
669        lap(i) = sig(i+1) + sig(i-1) - 2*sig(i)
670    enddo
671!
672    do i = ideb + Agrif_Minwidth-1,ifin-Agrif_Minwidth
673        if ((lap(i+1)*lap(i)) <= 0) then
674            ts = ABS(lap(i+1) - lap(i))
675            if (ts > t) then
676                t = ts
677                a = i
678            endif
679       endif
680    enddo
681!
682    if (a == (ideb + Agrif_Minwidth - 1)) then
683        a = -1
684        t = -1
685    endif
686!
687    slice = a
688    str = t
689!---------------------------------------------------------------------------------------------------
690end subroutine Agrif_Clusterslice
691!===================================================================================================
692!
693!===================================================================================================
694!  subroutine Agrif_Clusterprune
695!---------------------------------------------------------------------------------------------------
696subroutine Agrif_Clusterprune ( sig, pl, pu )
697!---------------------------------------------------------------------------------------------------
698    INTEGER,                   intent(inout)    :: pl, pu
699    INTEGER, DIMENSION(pl:pu), intent(in)       :: sig
700!
701    INTEGER :: ideb, ifin
702    INTEGER :: diff, addl, addu, udist, ldist
703!
704    ideb = pl
705    ifin = pu
706!
707    if (SIZE(sig) <= Agrif_Minwidth) return
708!
709    do while ((sig(pl) == 0) .AND. (pl < ifin))
710        pl = pl + 1
711    enddo
712!
713    do while ((sig(pu) == 0) .AND. (pu > ideb))
714        pu = pu - 1
715    enddo
716!
717    if ( (pu-pl) < Agrif_Minwidth ) then
718        diff = Agrif_Minwidth - (pu - pl + 1)
719        udist = ifin - pu
720        ldist = pl - ideb
721        addl = diff / 2
722        addu = diff - addl
723        if (addu > udist) then
724            addu = udist
725            addl = diff - addu
726        endif
727!
728        if (addl > ldist) then
729            addl = ldist
730            addu = diff - addl
731        endif
732!
733        pu = pu + addu
734        pl = pl - addl
735    endif
736!---------------------------------------------------------------------------------------------------
737end subroutine Agrif_Clusterprune
738!===================================================================================================
739!
740!===================================================================================================
741!  subroutine Agrif_Add_Rectangle
742!
743!> Adds the Agrif_Rectangle R in a list managed by LR.
744!---------------------------------------------------------------------------------------------------
745subroutine Agrif_Add_Rectangle ( R, LR )
746!---------------------------------------------------------------------------------------------------
747    TYPE(Agrif_Rectangle)           :: R
748    TYPE(Agrif_LRectangle), pointer :: LR
749!
750    TYPE(Agrif_LRectangle), pointer :: newrect
751!
752    integer                         :: i
753!
754    allocate(newrect)
755    allocate(newrect % r)
756!
757    newrect % r = R
758!
759    do i = 1,Agrif_Probdim
760        newrect % r % spaceref(i) = Agrif_Coeffref(i)
761        newrect % r % timeref(i)  = Agrif_Coeffreft(i)
762    enddo
763!
764    newrect % r % number = -1
765    newrect % next => LR
766    LR => newrect
767!---------------------------------------------------------------------------------------------------
768end subroutine Agrif_Add_Rectangle
769!===================================================================================================
770!
771!===================================================================================================
772!  subroutine Agrif_Copy_Rectangle
773!
774!> Creates and returns a copy of Agrif_Rectangle R.
775!---------------------------------------------------------------------------------------------------
776function Agrif_Copy_Rectangle ( R, expand ) result( copy )
777!---------------------------------------------------------------------------------------------------
778    TYPE(Agrif_Rectangle), pointer, intent(in) :: R
779    integer, optional,              intent(in) :: expand
780!
781    TYPE(Agrif_Rectangle), pointer :: copy
782!
783    allocate(copy)
784!
785    copy = R
786    if ( present(expand) ) then
787        copy % imin = copy % imin - expand
788        copy % imax = copy % imax + expand
789    endif
790    copy % childgrids => R % childgrids
791!---------------------------------------------------------------------------------------------------
792end function Agrif_Copy_Rectangle
793!===================================================================================================
794!
795!===================================================================================================
796!  subroutine Agrif_Read_Fix_Grd
797!
798!> Creates the grid hierarchy from the reading of the AGRIF_FixedGrids.in file.
799!---------------------------------------------------------------------------------------------------
800recursive subroutine Agrif_Read_Fix_Grd ( parent_rect, j, nunit )
801!---------------------------------------------------------------------------------------------------
802    TYPE(Agrif_Rectangle), pointer   :: parent_rect !< Pointer on the first grid of the grid hierarchy
803    INTEGER                          :: j          !< Number of the new grid
804    INTEGER                          :: nunit      !< unit associated with file
805!
806    TYPE(Agrif_Rectangle)            :: newrect    ! Pointer on a new grid
807    TYPE(Agrif_LRectangle), pointer  :: parcours   ! Pointer for the recursive procedure
808    TYPE(Agrif_LRectangle), pointer  :: newlrect
809    TYPE(Agrif_LRectangle), pointer  :: end_list
810    INTEGER                          :: i,n        ! for each child grid
811    INTEGER                          :: nb_grids   ! Number of child grids
812!
813!   Reading of the number of child grids
814    read(nunit,*,end=99) nb_grids
815!
816    allocate(end_list)
817!
818    parent_rect % childgrids => end_list
819!
820!   Reading of imin(1),imax(1),imin(2),imax(2),imin(3),imax(3), and space and
821!        time refinement factors for each child grid.
822!   Creation and addition of the new grid to the grid hierarchy.
823!
824    do n = 1,nb_grids
825!
826        allocate(newlrect)
827        newrect % number = j   ! Number of the grid
828!
829        if ( Agrif_USE_ONLY_FIXED_GRIDS == 0 ) then
830            if (Agrif_Probdim == 3) then
831                read(nunit,*) newrect % imin(1), newrect % imax(1), &
832                              newrect % imin(2), newrect % imax(2), &
833                              newrect % imin(3), newrect % imax(3), &
834                              newrect % spaceref(1), newrect % spaceref(2), newrect % spaceref(3), &
835                              newrect % timeref(1), newrect % timeref(2), newrect % timeref(3)
836            elseif (Agrif_Probdim == 2) then
837                read(nunit,*) newrect % imin(1), newrect % imax(1), &
838                              newrect % imin(2), newrect % imax(2), &
839                              newrect % spaceref(1), newrect % spaceref(2), &
840                              newrect % timeref(1),  newrect % timeref(2)
841            elseif (Agrif_Probdim == 1) then
842                read(nunit,*) newrect % imin(1), newrect % imax(1), &
843                              newrect % spaceref(1), &
844                              newrect % timeref(1)
845            endif
846        else
847            if (Agrif_Probdim == 3) then
848                read(nunit,*) newrect % imin(1), newrect % imax(1), &
849                              newrect % imin(2), newrect % imax(2), &
850                              newrect % imin(3), newrect % imax(3), &
851                              newrect % spaceref(1), newrect % spaceref(2), newrect % spaceref(3), &
852                              newrect % timeref(1)
853            elseif (Agrif_Probdim == 2) then
854                read(nunit,*) newrect % imin(1), newrect % imax(1), &
855                              newrect % imin(2), newrect % imax(2), &
856                              newrect % spaceref(1), newrect % spaceref(2), &
857                              newrect % timeref(1)
858            elseif (Agrif_Probdim == 1) then
859                read(nunit,*) newrect % imin(1), newrect % imax(1), &
860                              newrect % spaceref(1), &
861                              newrect % timeref(1)
862            endif
863!
864            if ( Agrif_probdim >= 2 ) then
865                do i = 2,Agrif_probdim
866                    newrect % timeref(i) = newrect % timeref(1)
867                enddo
868            endif
869!
870        endif
871!
872!       Addition to the grid hierarchy
873!
874        j = j + 1
875        allocate(newlrect%r)
876        newlrect % r = newrect
877        end_list % next => newlrect
878        end_list => end_list % next
879!
880    enddo
881!
882    parent_rect % childgrids => parent_rect % childgrids % next
883    parcours => parent_rect % childgrids
884!
885!   recursive operation to create the grid hierarchy branch by branch
886!
887    do while ( associated(parcours) )
888        call Agrif_Read_Fix_Grd(parcours % r, j, nunit)
889        parcours => parcours % next
890    enddo
89199  continue
892!---------------------------------------------------------------------------------------------------
893end subroutine Agrif_Read_Fix_Grd
894!===================================================================================================
895!
896!===================================================================================================
897!  subroutine Agrif_Create_Grids
898!
899!> Creates the grid hierarchy (g) from the one created with the #Agrif_Read_Fix_Grd or
900!! #Agrif_Cluster_All procedures (parent_rect).
901!---------------------------------------------------------------------------------------------------
902recursive subroutine Agrif_Create_Grids ( parent_grid, parent_rect )
903!---------------------------------------------------------------------------------------------------
904    TYPE(Agrif_Grid)     , pointer  :: parent_grid  !< Pointer on the root coarse grid
905    TYPE(Agrif_Rectangle), pointer  :: parent_rect  !< Pointer on the root coarse grid of the grid hierarchy
906                                                    !!   created with the #Agrif_Read_Fix_Grd subroutine
907!
908    TYPE(Agrif_Grid)      , pointer :: child_grid   ! Newly created child grid
909    TYPE(Agrif_PGrid)     , pointer :: child_grid_p
910    TYPE(Agrif_LRectangle), pointer :: child_rect_p
911    type(Agrif_Rectangle),  pointer :: child_rect
912!
913    INTEGER                         :: i
914    INTEGER, save                   :: moving_grid_id = 0
915!
916    child_rect_p => parent_rect % childgrids
917!
918!   Creation of the grid hierarchy from the one created by using the Agrif_Read_Fix_Grd subroutine
919!
920    do while ( associated(child_rect_p) )
921!
922        child_rect => child_rect_p % r
923!
924        allocate(child_grid)
925!
926!       Pointer on the parent grid
927        child_grid % parent => parent_grid
928        child_grid % rect_in_parent => Agrif_Copy_Rectangle(child_rect_p % r, expand=Agrif_Extra_Boundary_Cells)
929!
930        moving_grid_id = moving_grid_id+1
931        child_grid % grid_id = moving_grid_id
932!
933        do i = 1,Agrif_Probdim
934            child_grid % spaceref(i) = child_rect % spaceref(i)
935            child_grid % timeref(i)  = child_rect % timeref(i)
936            child_grid % nb(i) = (child_rect % imax(i) - child_rect % imin(i)) * child_rect % spaceref(i)
937            child_grid % ix(i) =  child_rect % imin(i)
938            child_grid % Agrif_dt(i) = parent_grid % Agrif_dt(i) / REAL(child_grid % timeref(i))
939            child_grid % Agrif_dx(i) = parent_grid % Agrif_dx(i) / REAL(child_grid % spaceref(i))
940            child_grid % Agrif_x(i)  = parent_grid % Agrif_x(i) + &
941                                            (child_rect % imin(i) - 1) * parent_grid % Agrif_dx(i)
942        enddo
943!
944!       Size of the grid in terms of cpu cost (nx*ny*timeref)
945        child_grid % size = product( child_grid % nb(1:Agrif_Probdim) ) * child_grid % timeref(1)
946!
947!       Level of the current grid
948        child_grid % level = child_grid % parent % level + 1
949        if (child_grid % level > Agrif_MaxLevelLoc) then
950          Agrif_MaxLevelLoc = child_grid%level
951        endif
952!
953!       Number of the grid pointed by child_grid
954        child_grid % fixedrank = child_rect % number
955!
956!       Grid pointed by child_grid is a fixed grid
957        child_grid % fixed = ( child_grid % fixedrank > 0 )
958!
959!       Update the total number of fixed grids
960        if (child_grid % fixed) then
961            Agrif_nbfixedgrids = Agrif_nbfixedgrids + 1
962        endif
963!
964!       Initialize integration counter
965        child_grid % ngridstep = 0
966!
967!       Test indicating if the current grid has a common border with the root
968!       coarse grid in the x direction
969        do i = 1 , Agrif_Probdim
970!
971            child_grid % NearRootBorder(i) = (                  &
972                (child_grid % parent % NearRootBorder(i)) .AND. &
973                (child_grid % ix(i) == 1) )
974!
975            child_grid % DistantRootBorder(i) = (                  &
976                (child_grid % parent % DistantRootBorder(i)) .AND. &
977                (child_grid % ix(i) + (child_grid%nb(i)/child_grid%spaceref(i))-1  == child_grid % parent % nb(i)) )
978        enddo
979!
980!       Writing in output files
981        child_grid % oldgrid = .FALSE.
982!
983#if defined AGRIF_MPI
984        child_grid % communicator = parent_grid % communicator
985#endif
986!
987!       Definition of the characteristics of the variable of the grid pointed by child_grid
988        call Agrif_Create_Var(child_grid)
989!
990!       Addition of this grid to the grid hierarchy
991        call Agrif_gl_append( parent_grid % child_list, child_grid )
992!
993        child_rect_p => child_rect_p % next
994!
995    enddo
996!
997!   Recursive call to the subroutine Agrif_Create_Fixed_Grids to create the grid hierarchy
998!
999    child_grid_p => parent_grid % child_list % first
1000    child_rect_p => parent_rect % childgrids
1001!
1002    do while ( associated(child_rect_p) )
1003        call Agrif_Create_Grids( child_grid_p % gr, child_rect_p % r )
1004        child_grid_p => child_grid_p % next
1005        child_rect_p => child_rect_p % next
1006    enddo
1007!---------------------------------------------------------------------------------------------------
1008end subroutine Agrif_Create_Grids
1009!===================================================================================================
1010!
1011!===================================================================================================
1012!  subroutine Agrif_Init_Hierarchy
1013!
1014!> Initializes all the grids except the root coarse grid (this one, pointed by Agrif_Types::Agrif_Mygrid, is
1015!! initialized by the subroutine Agrif_Util#Agrif_Init_Grids defined in the module Agrif_Util and
1016!! called in the main program).
1017!---------------------------------------------------------------------------------------------------
1018recursive subroutine Agrif_Init_Hierarchy ( g, procname )
1019!---------------------------------------------------------------------------------------------------
1020    use Agrif_seq
1021!
1022    type(Agrif_Grid), pointer       :: g            !< Pointer on the current grid
1023    procedure(init_proc), optional  :: procname     !< Initialisation subroutine (Default: Agrif_InitValues)
1024!
1025    TYPE(Agrif_PGrid), pointer :: parcours  ! Pointer for the recursive call
1026    LOGICAL                    :: Init_Hierarchy
1027!
1028! Initialise the grand mother grid (if any)
1029!
1030    if ( associated(g, Agrif_Mygrid) .and. agrif_coarse ) then
1031        call Agrif_Instance(Agrif_Coarsegrid)
1032        call Agrif_Allocation(Agrif_Coarsegrid)
1033        call Agrif_initialisations(Agrif_Coarsegrid)
1034        call Agrif_InitWorkspace()
1035!
1036!       Initialization by interpolation (this routine is written by the user)
1037        if (present(procname)) Then
1038            call procname()
1039        else
1040            call Agrif_InitValues()
1041        endif
1042        call Agrif_Instance(Agrif_Mygrid)
1043    endif
1044
1045    parcours => g % child_list % first
1046!
1047    do while ( associated(parcours) )
1048!
1049        Init_Hierarchy = .false.
1050        if ( Agrif_USE_FIXED_GRIDS == 1 .OR. Agrif_USE_ONLY_FIXED_GRIDS == 1 ) then
1051            if ( (parcours%gr%fixed) .AND. (Agrif_Mygrid%ngridstep == 0) ) then
1052                Init_Hierarchy = .true.
1053            endif
1054        endif
1055!
1056        if (.NOT. parcours % gr % fixed) Init_Hierarchy = .true.
1057        if (parcours % gr % oldgrid)     Init_Hierarchy = .false.
1058!
1059        if (Init_Hierarchy) then
1060!
1061!           Instanciation of the grid pointed by parcours%gr and its variables
1062            call Agrif_Instance(parcours % gr)
1063!
1064!           Allocation of the arrays containing values of the variables of the
1065!           grid pointed by parcours%gr
1066!
1067            call Agrif_Allocation(parcours % gr)
1068            call Agrif_initialisations(parcours % gr)
1069!
1070            if ( Agrif_USE_ONLY_FIXED_GRIDS == 0 ) then
1071!               Initialization by copy of the grids created by clustering
1072                call Agrif_Allocate_Restore(parcours % gr)
1073                call Agrif_CopyFromOld_All(parcours % gr, Agrif_oldmygrid)
1074            endif
1075!
1076!           Initialization by interpolation (this routine is written by the user)
1077            call Agrif_InitWorkSpace()
1078            if (present(procname)) Then
1079                call procname()
1080            else
1081                call Agrif_InitValues()
1082            endif
1083!
1084            if ( Agrif_USE_ONLY_FIXED_GRIDS == 0 ) then
1085                call Agrif_Free_Restore(parcours % gr)
1086            endif
1087!
1088        endif
1089!
1090        parcours => parcours % next
1091!
1092    enddo
1093!
1094    parcours => g % child_list % first
1095!
1096!   recursive operation to initialize all the grids
1097    do while ( associated(parcours) )
1098        call Agrif_Init_Hierarchy(parcours%gr,procname)
1099        parcours => parcours%next
1100    enddo
1101!---------------------------------------------------------------------------------------------------
1102end subroutine Agrif_Init_Hierarchy
1103!===================================================================================================
1104!
1105#if defined AGRIF_MPI
1106!===================================================================================================
1107!  subroutine Agrif_Init_Hierarchy_Parallel_1
1108!---------------------------------------------------------------------------------------------------
1109recursive subroutine Agrif_Init_Hierarchy_Parallel_1 ( g )
1110!---------------------------------------------------------------------------------------------------
1111    use Agrif_seq
1112!
1113    type(Agrif_Grid), pointer       :: g            !< Pointer on the current grid
1114!
1115    TYPE(Agrif_PGrid), pointer :: parcours  ! Pointer for the recursive call
1116    LOGICAL                    :: Init_Hierarchy
1117!
1118    parcours => g % child_list % first
1119!
1120    do while ( associated(parcours) )
1121!
1122        Init_Hierarchy = .false.
1123        if ( Agrif_USE_FIXED_GRIDS == 1 .OR. Agrif_USE_ONLY_FIXED_GRIDS == 1 ) then
1124            if ( (parcours%gr%fixed) .AND. (Agrif_Mygrid%ngridstep == 0) ) then
1125                Init_Hierarchy = .true.
1126            endif
1127        endif
1128!
1129        if (.NOT. parcours % gr % fixed) Init_Hierarchy = .true.
1130        if (parcours % gr % oldgrid)     Init_Hierarchy = .false.
1131!
1132        if (Init_Hierarchy) then
1133!
1134!           Instanciation of the grid pointed by parcours%gr and its variables
1135            call Agrif_Instance(parcours % gr)
1136!
1137!           Allocation of the arrays containing values of the variables of the
1138!           grid pointed by parcours%gr
1139!
1140            call Agrif_Allocation(parcours % gr)
1141            call Agrif_initialisations(parcours % gr)
1142!
1143        endif
1144!
1145        parcours => parcours % next
1146!
1147    enddo
1148!
1149    parcours => g % child_list % first
1150!
1151!   recursive operation to initialize all the grids
1152    do while ( associated(parcours) )
1153        call Agrif_Init_Hierarchy_Parallel_1(parcours%gr)
1154        parcours => parcours%next
1155    enddo
1156!
1157!---------------------------------------------------------------------------------------------------
1158end subroutine Agrif_Init_Hierarchy_Parallel_1
1159!===================================================================================================
1160!
1161!===================================================================================================
1162!  subroutine Agrif_Init_Hierarchy_Parallel_2
1163!---------------------------------------------------------------------------------------------------
1164recursive subroutine Agrif_Init_Hierarchy_Parallel_2 ( g, procname )
1165!---------------------------------------------------------------------------------------------------
1166    use Agrif_seq
1167!
1168    type(Agrif_Grid),     pointer   :: g            !< Pointer on the current grid
1169    procedure(init_proc), optional  :: procname     !< Initialisation subroutine (Default: Agrif_InitValues)
1170!
1171    type(Agrif_PGrid), pointer :: parcours  ! Pointer for the recursive call
1172    integer :: is
1173!
1174    call Agrif_Instance(g)
1175    call Agrif_seq_init_sequences( g )
1176!
1177    if ( .not. associated(g % child_seq) ) return
1178!
1179    do is = 1, g % child_seq % nb_seqs
1180!
1181        parcours => Agrif_seq_select_child(g,is)
1182!
1183!     Instanciation of the variables of the current grid
1184        call Agrif_Instance(parcours % gr)
1185!
1186!     Initialization by interpolation (this routine is written by the user)
1187        if (present(procname)) Then
1188            call procname()
1189        else
1190            call Agrif_InitValues()
1191        endif
1192!
1193        call Agrif_Init_ProcList(parcours % gr % proc_def_list, &
1194                                 parcours % gr % proc_def_in_parent_list % nitems)
1195!
1196        call Agrif_Init_Hierarchy_Parallel_2(parcours%gr,procname)
1197!
1198    enddo
1199!---------------------------------------------------------------------------------------------------
1200end subroutine Agrif_Init_Hierarchy_Parallel_2
1201!===================================================================================================
1202#endif
1203!
1204!===================================================================================================
1205!  subroutine Agrif_Allocate_Restore
1206!---------------------------------------------------------------------------------------------------
1207subroutine Agrif_Allocate_Restore ( Agrif_Gr )
1208!---------------------------------------------------------------------------------------------------
1209    TYPE(Agrif_Grid), pointer  :: Agrif_Gr   !< Pointer on the root coarse grid
1210!
1211    integer :: i
1212!
1213    do i = 1,Agrif_NbVariables(0)
1214!
1215        if ( Agrif_Mygrid%tabvars(i) % restore ) then
1216            if ( Agrif_Gr%tabvars(i) % nbdim == 1 ) then
1217                allocate( Agrif_Gr%tabvars(i)%Restore1D( &
1218                    lbound(Agrif_Gr%tabvars(i)%array1,1):&
1219                    ubound(Agrif_Gr%tabvars(i)%array1,1)))
1220                Agrif_Gr%tabvars(i)%Restore1D = 0
1221            endif
1222            if ( Agrif_Gr%tabvars(i) % nbdim == 2 ) then
1223                allocate( Agrif_Gr%tabvars(i)%Restore2D( &
1224                    lbound(Agrif_Gr%tabvars(i)%array2,1):&
1225                    ubound(Agrif_Gr%tabvars(i)%array2,1),&
1226                    lbound(Agrif_Gr%tabvars(i)%array2,2):&
1227                    ubound(Agrif_Gr%tabvars(i)%array2,2)))
1228                Agrif_Gr%tabvars(i)%Restore2D = 0
1229            endif
1230            if ( Agrif_Mygrid%tabvars(i) % nbdim == 3 ) then
1231                allocate( Agrif_Gr%tabvars(i)%Restore3D( &
1232                    lbound(Agrif_Gr%tabvars(i)%array3,1):&
1233                    ubound(Agrif_Gr%tabvars(i)%array3,1),&
1234                    lbound(Agrif_Gr%tabvars(i)%array3,2):&
1235                    ubound(Agrif_Gr%tabvars(i)%array3,2),&
1236                    lbound(Agrif_Gr%tabvars(i)%array3,3):&
1237                    ubound(Agrif_Gr%tabvars(i)%array3,3)))
1238                Agrif_Gr%tabvars(i)%Restore3D = 0
1239            endif
1240        endif
1241!
1242    enddo
1243!---------------------------------------------------------------------------------------------------
1244end subroutine Agrif_Allocate_Restore
1245!===================================================================================================
1246!
1247!===================================================================================================
1248!  subroutine Agrif_Free_Restore
1249!---------------------------------------------------------------------------------------------------
1250subroutine Agrif_Free_Restore ( Agrif_Gr )
1251!---------------------------------------------------------------------------------------------------
1252    TYPE(Agrif_Grid), pointer  :: Agrif_Gr   !< Pointer on the root coarse grid
1253!
1254    TYPE(Agrif_Variable), pointer :: var
1255    integer                       :: i
1256!
1257    do i = 1,Agrif_NbVariables(0)
1258!
1259        if ( Agrif_Mygrid % tabvars(i) % restore ) then
1260!
1261            var = Agrif_Gr % tabvars(i)
1262!
1263            if (associated(var%Restore1D)) deallocate(var%Restore1D)
1264            if (associated(var%Restore2D)) deallocate(var%Restore2D)
1265            if (associated(var%Restore3D)) deallocate(var%Restore3D)
1266            if (associated(var%Restore4D)) deallocate(var%Restore4D)
1267            if (associated(var%Restore5D)) deallocate(var%Restore5D)
1268            if (associated(var%Restore6D)) deallocate(var%Restore6D)
1269!
1270        endif
1271!
1272    enddo
1273!---------------------------------------------------------------------------------------------------
1274end subroutine Agrif_Free_Restore
1275!===================================================================================================
1276!
1277end module Agrif_Clustering
Note: See TracBrowser for help on using the repository browser.