New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
modinterp.F90 in vendors/AGRIF/CMEMS_2020/AGRIF_FILES – NEMO

source: vendors/AGRIF/CMEMS_2020/AGRIF_FILES/modinterp.F90 @ 10087

Last change on this file since 10087 was 10087, checked in by rblod, 6 years ago

update AGRIF library

  • Property svn:keywords set to Id
File size: 90.8 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 to initialize a fine grid from its parent grid, by using a space interpolation
25!
26module Agrif_Interpolation
27!
28  use Agrif_Init
29  use Agrif_Arrays
30  use Agrif_InterpBasic
31  use Agrif_User_Functions
32
33#if defined AGRIF_MPI
34    use Agrif_Mpp
35#endif
36 
37    use Agrif_Mask
38!
39    implicit none
40!
41    logical, private:: precomputedone(7) = .FALSE.
42!
43    private :: Agrif_Parentbounds
44    private :: Agrif_Interp_1D_recursive, Agrif_Interp_2D_recursive, Agrif_Interp_3D_recursive
45    private :: Agrif_Interp_4D_recursive, Agrif_Interp_5D_recursive, Agrif_Interp_6D_recursive
46    private :: Agrif_InterpBase
47    private :: Agrif_Find_list_interp, Agrif_AddTo_list_interp
48!
49contains
50!
51!===================================================================================================
52!  subroutine Agrif_InterpVariable
53!
54!> Sets some arguments of subroutine Agrif_InterpnD, n being the dimension of the grid variable
55!---------------------------------------------------------------------------------------------------
56subroutine Agrif_InterpVariable ( parent, child, torestore, procname )
57!---------------------------------------------------------------------------------------------------
58    type(Agrif_Variable), pointer       :: parent       !< Variable on the parent grid
59    type(Agrif_Variable), pointer       :: child        !< Variable on the child grid
60    logical,              intent(in)    :: torestore    !< .false. indicates that the results of the
61                                                        !! interpolation are applied on the whole current grid
62    procedure()                         :: procname     !< Data recovery procedure
63!---------------------------------------------------------------------------------------------------
64    logical               :: memberin
65    integer               :: nbdim        ! Number of dimensions of the current grid
66    integer, dimension(6) :: type_interp  ! Type of interpolation (linear,spline,...)
67    integer, dimension(6) :: nb_child
68    integer, dimension(6) :: lb_child
69    integer, dimension(6) :: ub_child
70    integer, dimension(6) :: lb_parent
71    real(kind=8)   , dimension(6) :: s_child,   s_parent
72    real(kind=8)   , dimension(6) :: ds_child, ds_parent
73    integer, dimension(child % root_var % nbdim,2,2) :: childarray
74!
75    nbdim       = child % root_var % nbdim
76    type_interp = child % root_var % type_interp
77!
78    call PreProcessToInterpOrUpdate( parent,   child,       &
79                                     nb_child, ub_child,    &
80                                     lb_child, lb_parent,   &
81                                      s_child,  s_parent,   &
82                                     ds_child, ds_parent, nbdim, interp=.true.)
83!
84!   Call to a procedure of interpolation against the number of dimensions of the grid variable
85!
86    call Agrif_InterpnD(type_interp, parent, child, &
87                        lb_child, ub_child,         &
88                        lb_child, lb_parent,        &
89                        s_child,   s_parent,        &
90                        ds_child, ds_parent,        &
91                        child, torestore, nbdim,    &
92                        childarray, memberin,       &
93                        .false., procname, 0, 0)
94!---------------------------------------------------------------------------------------------------
95end subroutine Agrif_InterpVariable
96!===================================================================================================
97!
98!===================================================================================================
99!  subroutine Agrif_InterpnD
100!
101!> Interpolates a nD grid variable from its parent grid, by using a space interpolation
102!---------------------------------------------------------------------------------------------------
103subroutine Agrif_InterpnD ( type_interp, parent, child, pttab, petab, pttab_Child, pttab_Parent, &
104                            s_Child, s_Parent, ds_Child, ds_Parent, restore, torestore,          &
105                            nbdim, childarray, memberin, in_bc, procname, nb, ndir )
106!---------------------------------------------------------------------------------------------------
107#if defined AGRIF_MPI
108    include 'mpif.h'
109#endif
110!
111    INTEGER, DIMENSION(6),     INTENT(in)   :: type_interp  !< Type of interpolation !    (linear,...)
112    TYPE(Agrif_Variable),      pointer      :: parent       !< Variable of the parent grid
113    TYPE(Agrif_Variable),      pointer      :: child        !< Variable of the child grid
114    INTEGER, DIMENSION(nbdim), INTENT(in)   :: pttab        !< Index of the first point inside the domain
115    INTEGER, DIMENSION(nbdim), INTENT(in)   :: petab        !< Index of the first point inside the domain
116    INTEGER, DIMENSION(nbdim), INTENT(in)   :: pttab_Child  !< Index of the first point inside the domain
117                                                            !<    for the child grid variable
118    INTEGER, DIMENSION(nbdim), INTENT(in)   :: pttab_Parent !< Index of the first point inside the domain
119                                                            !<    for the parent grid variable
120    REAL(kind=8),    DIMENSION(nbdim), INTENT(in)   :: s_Child,s_Parent   !< Positions of the parent and child grids
121    REAL(kind=8),    DIMENSION(nbdim), INTENT(in)   :: ds_Child,ds_Parent !< Space steps of the parent and child grids
122    TYPE(Agrif_Variable),      pointer      :: restore            !< Indicates points where interpolation
123    LOGICAL,                   INTENT(in)   :: torestore          !< Indicates if the array restore is used
124    INTEGER,                   INTENT(in)   :: nbdim
125    LOGICAL,                   INTENT(out)  :: memberin
126    LOGICAL,                   INTENT(in)   :: in_bc              !< .true. if called from Agrif_CorrectVariable \n
127                                                                  !! .false. if called from Agrif_InterpVariable
128    procedure()                             :: procname           !< Data recovery procedure
129    INTEGER,                   INTENT(in)   :: nb, ndir
130!
131    INTEGER                       :: i,j,k,l,m,n
132    INTEGER, DIMENSION(nbdim)     :: pttruetab,cetruetab
133    INTEGER, DIMENSION(nbdim)     :: indmin,     indmax, indmin_required_p, indmax_required_p
134    INTEGER, DIMENSION(nbdim)     :: indminglob, indmaxglob
135#if defined AGRIF_MPI
136    INTEGER, DIMENSION(nbdim)     :: indminglob2,indmaxglob2
137    INTEGER, DIMENSION(nbdim)     :: indminglob3,indmaxglob3
138#endif
139    LOGICAL, DIMENSION(nbdim)     :: noraftab
140    REAL(kind=8)   , DIMENSION(nbdim)     :: s_Child_temp,s_Parent_temp,s_Parent_temp_p
141    INTEGER, DIMENSION(nbdim)     :: lowerbound, upperbound, coords
142    INTEGER, DIMENSION(nbdim,2,2), INTENT(OUT) :: childarray
143    INTEGER, DIMENSION(nbdim,2,2)              :: parentarray
144    LOGICAL :: member
145    LOGICAL :: find_list_interp
146!
147#if defined AGRIF_MPI
148!
149    INTEGER, PARAMETER          :: etiquette = 100
150    INTEGER                     :: code, local_proc
151    INTEGER, DIMENSION(nbdim,4)                     :: tab3
152    INTEGER, DIMENSION(nbdim,4,0:Agrif_Nbprocs-1)   :: tab4
153    INTEGER, DIMENSION(nbdim,0:Agrif_Nbprocs-1,8)   :: tab4t
154    INTEGER,DIMENSION(nbdim,2) :: tab5
155    INTEGER,DIMENSION(nbdim,2,0:Agrif_Nbprocs-1) :: tab6
156    INTEGER,DIMENSION(nbdim,0:Agrif_Nbprocs-1,2) :: tab5t
157    LOGICAL, DIMENSION(0:Agrif_Nbprocs-1)           :: memberinall
158    LOGICAL, DIMENSION(0:Agrif_Nbprocs-1)           :: sendtoproc1, recvfromproc1
159    LOGICAL, DIMENSION(1)                           :: memberin1
160    LOGICAL                                         :: memberout
161!
162#endif
163!
164    type(Agrif_Variable), pointer, save :: tempC => NULL()        ! Temporary child grid variable
165    type(Agrif_Variable), pointer, save :: tempP => NULL()        ! Temporary parent grid variable
166    type(Agrif_Variable), pointer, save :: tempPextend => NULL()  ! Temporary parent grid variable
167    type(Agrif_Variable), pointer, save :: parentvalues => NULL()
168!
169    coords = child % root_var % coords
170!
171!   Boundaries of the current grid where interpolation is done
172    find_list_interp =                                              &
173        Agrif_Find_list_interp(                                     &
174            child % list_interp,                                    &
175            pttab, petab, pttab_Child, pttab_Parent, nbdim,         &
176            indmin, indmax, indmin_required_p, indmax_required_p,   &
177            indminglob, indmaxglob,                                 &
178            pttruetab, cetruetab, memberin                          &
179#if defined AGRIF_MPI
180           ,indminglob2, indmaxglob2, parentarray,                  &
181            member, tab4t,memberinall,  sendtoproc1, recvfromproc1  &
182#endif
183        )
184
185!
186    if (.not.find_list_interp) then
187!
188! output : lowerbound and upperbound are the (local) lower and upper bounds of the child arrays
189
190        call Agrif_get_var_bounds_array(child, lowerbound, upperbound, nbdim)
191
192! input : pttab, petab : global indexes where the interpolation is required
193! output : pttruetab, cetruetab : global indexes restricted to the bounds of the current processor
194! output : memberin is false if the current processor is not involved in the interpolation
195
196        call Agrif_Childbounds(nbdim, lowerbound, upperbound,               &
197                               pttab, petab, Agrif_Procrank, coords,        &
198                               pttruetab, cetruetab, memberin)
199         
200
201
202! output : indminglob, indmaxglob : global indexes required on the parent grid for the interpolation
203! output : s_Parent_temp, s_Child_temp : local s_Parent, s_Child relatively to indmin ou pttab
204        call Agrif_Parentbounds(type_interp,nbdim,indminglob,indmaxglob,    &
205                                indmin_required_p, indmax_required_p,           &
206                                s_Parent_temp,s_Child_temp,                 &
207                                s_Child,ds_Child,                           &
208                                s_Parent,ds_Parent,                         &
209                                pttab,petab,                                &
210                                pttab_Child,pttab_Parent,                   &
211                                child%root_var % posvar, coords)
212#if defined AGRIF_MPI
213        if (memberin) then
214
215! output : indmin, indmax : global indexes required on the parent grid for the interpolation on the current processor (i.e. on pttruetab, cetruetab)
216! output : s_Parent_temp, s_Child_temp : local s_Parent, s_Child relatively to indmin ou pttruetab
217            call Agrif_Parentbounds(type_interp,nbdim,indmin,indmax,        &
218                                    indmin_required_p, indmax_required_p,       &
219                                    s_Parent_temp,s_Child_temp,             &
220                                    s_Child,ds_Child,                       &
221                                    s_Parent,ds_Parent,                     &
222                                    pttruetab,cetruetab,                    &
223                                    pttab_Child,pttab_Parent,               &
224                                    child%root_var % posvar, coords)
225        endif
226
227        local_proc = Agrif_Procrank
228
229! output : lowerbound and upperbound are the (local) lower and upper bounds of the parent arrays
230        call Agrif_get_var_bounds_array(parent,lowerbound,upperbound,nbdim)
231        call Agrif_ChildGrid_to_ParentGrid()
232
233! input : indminglob,indmaxglob : global indexes where data are required for the interpolation
234! output : indminglob2,indmaxglob2 : global indexes restricted to the bounds of the current processor
235! output : member is false if the current processor does not need to send data
236! output : indminglob3,indmaxglob3 : global bounds on the current processor for the parent array
237
238        call Agrif_Childbounds(nbdim,lowerbound,upperbound,                 &
239                               indminglob,indmaxglob, local_proc, coords,   &
240                               indminglob2,indmaxglob2,member,              &
241                               indminglob3,indmaxglob3)
242!
243        if (member) then
244
245! output : parentarray
246! output : parentarray (:,:,2) : indminglob2, indmaxglob2 in term of local indexes on current processor
247! output : parentarray (:,:,1) : indminglob2, indmaxglob2 restricted to the current processor (different from indminglob2 ???)
248! output : member is .false. is the current processor has not data to send
249
250            call Agrif_GlobalToLocalBounds(parentarray,                     &
251                                           lowerbound,  upperbound,         &
252                                           indminglob2, indmaxglob2, coords,&
253                                           nbdim, local_proc, member)
254        endif
255
256        call Agrif_ParentGrid_to_ChildGrid()
257#else
258
259! In the sequentiel case, the following lines ensure that the bounds needed on the parent grid in the interpolation
260! do not exceed lower and upper bounds of the parent array
261
262! output : lowerbound and upperbound are the (local) lower and upper bounds of the parent arrays
263        call Agrif_get_var_bounds_array(parent,lowerbound,upperbound,nbdim)
264        call Agrif_ChildGrid_to_ParentGrid()
265
266! input : indminglob,indmaxglob : global indexes where data are required for the interpolation
267! output : indminglob2,indmaxglob2 : global indexes restricted to the bounds of the current processor
268! output : member is false if the current processor does not need to send data
269
270        call Agrif_Childbounds(nbdim,lowerbound,upperbound,                 &
271                               indminglob,indmaxglob, Agrif_Procrank, coords,   &
272                               indmin,indmax,member)
273
274        call Agrif_ParentGrid_to_ChildGrid()
275
276        indminglob = indmin
277        indmaxglob = indmax
278
279        parentarray(:,1,1) = indminglob
280        parentarray(:,2,1) = indmaxglob
281        parentarray(:,1,2) = indminglob
282        parentarray(:,2,2) = indmaxglob
283 
284!       indmin = indminglob
285!        indmax = indmaxglob
286
287        member = .TRUE.
288        s_Parent_temp = s_Parent + (indminglob - pttab_Parent) * ds_Parent
289
290#endif
291!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
292! Correct for non refined directions
293!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
294        do i=1,nbdim
295          if (coords(i) == 0) then
296             indmin(i) = indminglob(i)
297             indmax(i) = indmaxglob(i)
298             pttruetab(i) = indminglob(i)
299             cetruetab(i) = indmaxglob(i)
300          endif
301        enddo
302
303    else
304
305#if defined AGRIF_MPI
306        s_Parent_temp = s_Parent + (indmin - pttab_Parent) * ds_Parent
307        s_Child_temp  = s_Child + (pttruetab - pttab_Child) * ds_Child
308#else
309
310! In the sequentiel case, the following lines ensure that the bounds needed on the parent grid in the interpolation
311! do not exceed lower and upper bounds of the parent array
312
313! output : lowerbound and upperbound are the (local) lower and upper bounds of the parent arrays
314        call Agrif_get_var_bounds_array(parent,lowerbound,upperbound,nbdim)
315        call Agrif_ChildGrid_to_ParentGrid()
316
317! input : indminglob,indmaxglob : global indexes where data are required for the interpolation
318! output : indminglob2,indmaxglob2 : global indexes restricted to the bounds of the current processor
319! output : member is false if the current processor does not need to send data
320
321        call Agrif_Childbounds(nbdim,lowerbound,upperbound,                 &
322                               indminglob,indmaxglob, Agrif_Procrank, coords,   &
323                               indmin,indmax,member)
324
325        call Agrif_ParentGrid_to_ChildGrid()
326
327        indminglob = indmin
328        indmaxglob = indmax
329
330        parentarray(:,1,1) = indminglob
331        parentarray(:,2,1) = indmaxglob
332        parentarray(:,1,2) = indminglob
333        parentarray(:,2,2) = indmaxglob
334 !       indmin = indminglob
335 !       indmax = indmaxglob
336        member = .TRUE.
337        s_Parent_temp = s_Parent + (indminglob - pttab_Parent) * ds_Parent
338        s_Child_temp  = s_Child + (pttab - pttab_Child) * ds_Child
339#endif
340    endif
341!
342    if (member) then
343        if (.not.associated(tempP)) allocate(tempP)
344!
345
346        call Agrif_array_allocate(tempP,parentarray(:,1,1),parentarray(:,2,1),nbdim)
347        call Agrif_var_set_array_tozero(tempP,nbdim)
348
349        call Agrif_ChildGrid_to_ParentGrid()
350!
351        select case (nbdim)
352        case(1)
353            call procname(tempP%array1,                         &
354                      parentarray(1,1,2),parentarray(1,2,2),.TRUE.,nb,ndir)
355        case(2)
356            call procname(tempP%array2,                         &
357                      parentarray(1,1,2),parentarray(1,2,2),    &
358                      parentarray(2,1,2),parentarray(2,2,2),.TRUE.,nb,ndir)
359        case(3)
360            call procname(tempP%array3,                         &
361                      parentarray(1,1,2),parentarray(1,2,2),    &
362                      parentarray(2,1,2),parentarray(2,2,2),    &
363                      parentarray(3,1,2),parentarray(3,2,2),.TRUE.,nb,ndir)
364        case(4)
365            call procname(tempP%array4,                         &
366                      parentarray(1,1,2),parentarray(1,2,2),    &
367                      parentarray(2,1,2),parentarray(2,2,2),    &
368                      parentarray(3,1,2),parentarray(3,2,2),    &
369                      parentarray(4,1,2),parentarray(4,2,2),.TRUE.,nb,ndir)
370        case(5)
371            call procname(tempP%array5,                         &
372                      parentarray(1,1,2),parentarray(1,2,2),    &
373                      parentarray(2,1,2),parentarray(2,2,2),    &
374                      parentarray(3,1,2),parentarray(3,2,2),    &
375                      parentarray(4,1,2),parentarray(4,2,2),    &
376                      parentarray(5,1,2),parentarray(5,2,2),.TRUE.,nb,ndir)
377        case(6)
378            call procname(tempP%array6,                         &
379                      parentarray(1,1,2),parentarray(1,2,2),    &
380                      parentarray(2,1,2),parentarray(2,2,2),    &
381                      parentarray(3,1,2),parentarray(3,2,2),    &
382                      parentarray(4,1,2),parentarray(4,2,2),    &
383                      parentarray(5,1,2),parentarray(5,2,2),    &
384                      parentarray(6,1,2),parentarray(6,2,2),.TRUE.,nb,ndir)
385        end select
386
387!
388        call Agrif_ParentGrid_to_ChildGrid()
389!
390    endif
391
392#if defined AGRIF_MPI
393    if (.not.find_list_interp) then
394!
395        tab3(:,1) = indminglob2(:)
396        tab3(:,2) = indmaxglob2(:)
397        tab3(:,3) = indmin(:)
398        tab3(:,4) = indmax(:)
399        tab5(:,1) = indminglob3(:)
400        tab5(:,2) = indmaxglob3(:)
401!
402        call MPI_ALLGATHER(tab3,4*nbdim,MPI_INTEGER,tab4,4*nbdim,MPI_INTEGER,Agrif_mpi_comm,code)
403        call MPI_ALLGATHER(tab5,2*nbdim,MPI_INTEGER,tab6,2*nbdim,MPI_INTEGER,Agrif_mpi_comm,code)
404        if (.not.associated(tempPextend))   allocate(tempPextend)
405
406        do k=0,Agrif_Nbprocs-1
407            do j=1,4
408                do i=1,nbdim
409                    tab4t(i,k,j) = tab4(i,j,k)
410                enddo
411            enddo
412        enddo
413
414        do k=0,Agrif_Nbprocs-1
415          do j=1,2
416            do i=1,nbdim
417               tab5t(i,k,j) = tab6(i,j,k)
418            enddo
419          enddo
420        enddo
421     
422        memberin1(1) = memberin
423        call MPI_ALLGATHER(memberin1,1,MPI_LOGICAL,memberinall,1,MPI_LOGICAL,Agrif_mpi_comm,code)
424
425        call Get_External_Data_first(tab4t(:,:,1),tab4t(:,:,2),         &
426                                     tab4t(:,:,3),tab4t(:,:,4),         &
427                                     nbdim,memberinall, coords,         &
428                                     sendtoproc1,recvfromproc1,         &
429                                     tab4t(:,:,5),tab4t(:,:,6),         &
430                                     tab4t(:,:,7),tab4t(:,:,8),         &
431                                     tab5t(:,:,1),tab5t(:,:,2))
432    endif
433
434    call ExchangeSameLevel(sendtoproc1,recvfromproc1,nbdim,         &
435            tab4t(:,:,3),tab4t(:,:,4),tab4t(:,:,5),tab4t(:,:,6),    &
436            tab4t(:,:,7),tab4t(:,:,8),memberin,tempP,tempPextend)
437#else
438    tempPextend => tempP
439#endif
440
441    if (.not.find_list_interp) then
442        call Agrif_Addto_list_interp(                           &
443                child%list_interp,pttab,petab,                  &
444                pttab_Child,pttab_Parent,indmin,indmax,         &
445                indmin_required_p, indmax_required_p,           &
446                indminglob,indmaxglob,                          &
447                pttruetab,cetruetab,                            &
448                memberin,nbdim                                  &
449#if defined AGRIF_MPI
450               ,indminglob2,indmaxglob2,                        &
451                parentarray,                                    &
452                member,                                         &
453                tab4t,memberinall,sendtoproc1,recvfromproc1     &
454#endif
455        )
456    endif
457!
458
459    if (memberin) then
460!
461        if (.not.associated(tempC)) allocate(tempC)
462!
463
464        call Agrif_array_allocate(tempC,pttruetab,cetruetab,nbdim)
465
466!
467!       Special values on the parent grid
468        if (Agrif_UseSpecialValue) then
469!
470            noraftab(1:nbdim) = child % root_var % interptab(1:nbdim) == 'N'
471!
472            if (.not.associated(parentvalues))  allocate(parentvalues)
473!t
474
475            call Agrif_array_allocate(parentvalues,indmin,indmax,nbdim)
476            call Agrif_var_full_copy_array(parentvalues,tempPextend,nbdim)
477!
478            call Agrif_CheckMasknD(tempPextend,parentvalues,                &
479                    indmin(1:nbdim),indmax(1:nbdim),                        &
480                    indmin(1:nbdim),indmax(1:nbdim),                        &
481                    indmin_required_p(1:nbdim),indmax_required_p(1:nbdim),  &
482                    noraftab(1:nbdim),nbdim)
483!
484            call Agrif_array_deallocate(parentvalues,nbdim)
485!
486        endif
487!
488!       Interpolation of the current grid
489!
490        if ( memberin ) then
491            select case(nbdim)
492            case(1)
493                call Agrif_Interp_1D_recursive( type_interp(1),                         &
494                                                tempPextend%array1,                     &
495                                                tempC%array1,                           &
496                                                indmin(1), indmax(1),                   &
497                                                pttruetab(1),    cetruetab(1),          &
498                                                s_Child_temp(1), s_Parent_temp(1),      &
499                                                ds_Child(1),     ds_Parent(1) )
500            case(2)
501                call Agrif_Interp_2D_recursive( type_interp(1:2),                       &
502                                                tempPextend % array2,                   &
503                                                tempC       % array2,                   &
504                                                indmin(1:2), indmax(1:2),               &
505                                                pttruetab(1:2),    cetruetab(1:2),      &
506                                                s_Child_temp(1:2), s_Parent_temp(1:2),  &
507                                                ds_Child(1:2),    ds_Parent(1:2) )
508            case(3)
509                s_Parent_temp_p = s_Parent + (indmin_required_p - pttab_Parent) * ds_Parent
510                call Agrif_Interp_3D_recursive( type_interp(1:3),                                 &
511                                                tempPextend % array3(                             &
512                                                indmin_required_p(1):indmax_required_p(1),        &
513                                                indmin_required_p(2):indmax_required_p(2),        &
514                                                indmin_required_p(3):indmax_required_p(3)),       &
515                                                tempC       % array3,                             &
516                                                indmin_required_p(1:3), indmax_required_p(1:3),   &
517                                                pttruetab(1:3),    cetruetab(1:3),                &
518                                                s_Child_temp(1:3), s_Parent_temp_p(1:3),          &
519                                                ds_Child(1:3),    ds_Parent(1:3) )
520
521            case(4)
522                s_Parent_temp_p = s_Parent + (indmin_required_p - pttab_Parent) * ds_Parent
523                call Agrif_Interp_4D_recursive( type_interp(1:4),                                 &
524                                                tempPextend % array4(                             &
525                                                indmin_required_p(1):indmax_required_p(1),        &
526                                                indmin_required_p(2):indmax_required_p(2),        &
527                                                indmin_required_p(3):indmax_required_p(3),        &
528                                                indmin_required_p(4):indmax_required_p(4)),       &
529                                                tempC       % array4,                             &
530                                                indmin_required_p(1:4), indmax_required_p(1:4),   &
531                                                pttruetab(1:4),    cetruetab(1:4),                &
532                                                s_Child_temp(1:4), s_Parent_temp_p(1:4),          &
533                                                ds_Child(1:4),    ds_Parent(1:4) )
534            case(5)
535                call Agrif_Interp_5D_recursive( type_interp(1:5),                       &
536                                                tempPextend % array5,                   &
537                                                tempC       % array5,                   &
538                                                indmin(1:5), indmax(1:5),               &
539                                                pttruetab(1:5),    cetruetab(1:5),      &
540                                                s_Child_temp(1:5), s_Parent_temp(1:5),  &
541                                                ds_Child(1:5),    ds_Parent(1:5) )
542            case(6)
543                call Agrif_Interp_6D_recursive( type_interp(1:6),                       &
544                                                tempPextend % array6,                   &
545                                                tempC       % array6,                   &
546                                                indmin(1:6), indmax(1:6),               &
547                                                pttruetab(1:6),    cetruetab(1:6),      &
548                                                s_Child_temp(1:6), s_Parent_temp(1:6),  &
549                                                ds_Child(1:6),    ds_Parent(1:6) )
550            end select
551!
552            call Agrif_get_var_bounds_array(child,lowerbound,upperbound,nbdim)
553
554#if defined AGRIF_MPI
555            call Agrif_GlobalToLocalBounds(childarray, lowerbound, upperbound,  &
556                                            pttruetab, cetruetab, coords,       &
557                                            nbdim, Agrif_Procrank, memberout)
558#else
559            childarray(:,1,1) = pttruetab
560            childarray(:,2,1) = cetruetab
561            childarray(:,1,2) = pttruetab
562            childarray(:,2,2) = cetruetab
563!cccccccccccccc       memberout = .TRUE.
564#endif
565!
566!           Special values on the child grid
567            if (Agrif_UseSpecialValueFineGrid) then
568                call GiveAgrif_SpecialValueToTab_mpi( child, tempC, childarray, Agrif_SpecialValueFineGrid,nbdim )
569            endif
570!
571        endif   ! ( memberin )
572!
573        if (torestore) then
574!
575#if defined AGRIF_MPI
576!
577            SELECT CASE (nbdim)
578            CASE (1)
579                do i = pttruetab(1),cetruetab(1)
580!hildarrayAModifier     if (restore%restore1D(i) == 0)      &
581!hildarrayAModifier         child%array1(childarray(i,1,2)) = tempC%array1(i)
582                enddo
583            CASE (2)
584                do i = pttruetab(1),cetruetab(1)
585                do j = pttruetab(2),cetruetab(2)
586!hildarrayAModifier     if (restore%restore2D(i,j) == 0)    &
587!hildarrayAModifier         child%array2(childarray(i,1,2), &
588!hildarrayAModifier                          childarray(j,2,2)) = tempC%array2(i,j)
589                enddo
590                enddo
591            CASE (3)
592                do i = pttruetab(1),cetruetab(1)
593                do j = pttruetab(2),cetruetab(2)
594                do k = pttruetab(3),cetruetab(3)
595!hildarrayAModifier     if (restore%restore3D(i,j,k) == 0)  &
596!hildarrayAModifier         child%array3(childarray(i,1,2), &
597!hildarrayAModifier                          childarray(j,2,2), &
598!hildarrayAModifier                          childarray(k,3,2)) = tempC%array3(i,j,k)
599                enddo
600                enddo
601                enddo
602            CASE (4)
603                do i = pttruetab(1),cetruetab(1)
604                do j = pttruetab(2),cetruetab(2)
605                do k = pttruetab(3),cetruetab(3)
606                do l = pttruetab(4),cetruetab(4)
607!hildarrayAModifier     if (restore%restore4D(i,j,k,l) == 0)    &
608!hildarrayAModifier         child%array4(childarray(i,1,2),     &
609!hildarrayAModifier                          childarray(j,2,2),     &
610!hildarrayAModifier                          childarray(k,3,2),     &
611!hildarrayAModifier                          childarray(l,4,2)) = tempC%array4(i,j,k,l)
612                enddo
613                enddo
614                enddo
615                enddo
616            CASE (5)
617                do i = pttruetab(1),cetruetab(1)
618                do j = pttruetab(2),cetruetab(2)
619                do k = pttruetab(3),cetruetab(3)
620                do l = pttruetab(4),cetruetab(4)
621                do m = pttruetab(5),cetruetab(5)
622!hildarrayAModifier     if (restore%restore5D(i,j,k,l,m) == 0)  &
623!hildarrayAModifier         child%array5(childarray(i,1,2),     &
624!hildarrayAModifier                          childarray(j,2,2),     &
625!hildarrayAModifier                          childarray(k,3,2),     &
626!hildarrayAModifier                          childarray(l,4,2),     &
627!hildarrayAModifier                          childarray(m,5,2)) = tempC%array5(i,j,k,l,m)
628                enddo
629                enddo
630                enddo
631                enddo
632                enddo
633            CASE (6)
634                do i = pttruetab(1),cetruetab(1)
635                do j = pttruetab(2),cetruetab(2)
636                do k = pttruetab(3),cetruetab(3)
637                do l = pttruetab(4),cetruetab(4)
638                do m = pttruetab(5),cetruetab(5)
639                do n = pttruetab(6),cetruetab(6)
640!hildarrayAModifier     if (restore%restore6D(i,j,k,l,m,n) == 0)    &
641!hildarrayAModifier         child%array6(childarray(i,1,2),         &
642!hildarrayAModifier                          childarray(j,2,2),         &
643!hildarrayAModifier                          childarray(k,3,2),         &
644!hildarrayAModifier                          childarray(l,4,2),         &
645!hildarrayAModifier                          childarray(m,5,2),         &
646!hildarrayAModifier                          childarray(n,6,2)) = tempC%array6(i,j,k,l,m,n)
647                enddo
648                enddo
649                enddo
650                enddo
651                enddo
652                enddo
653            END SELECT
654!
655#else
656            select case (nbdim)
657            case (1)
658                do i = pttruetab(1),cetruetab(1)
659                    if (restore%restore1D(i) == 0)          &
660                        parray1(i) = tempC % array1(i)
661                enddo
662            case (2)
663                do j = pttruetab(2),cetruetab(2)
664                do i = pttruetab(1),cetruetab(1)
665                    if (restore%restore2D(i,j) == 0)        &
666                        parray2(i,j) = tempC % array2(i,j)
667                enddo
668                enddo
669            case (3)
670                do k = pttruetab(3),cetruetab(3)
671                do j = pttruetab(2),cetruetab(2)
672                do i = pttruetab(1),cetruetab(1)
673                    if (restore%restore3D(i,j,k) == 0)      &
674                        parray3(i,j,k) = tempC % array3(i,j,k)
675                enddo
676                enddo
677                enddo
678            case (4)
679                do l = pttruetab(4),cetruetab(4)
680                do k = pttruetab(3),cetruetab(3)
681                do j = pttruetab(2),cetruetab(2)
682                do i = pttruetab(1),cetruetab(1)
683                    if (restore%restore4D(i,j,k,l) == 0)    &
684                        parray4(i,j,k,l) = tempC % array4(i,j,k,l)
685                enddo
686                enddo
687                enddo
688                enddo
689            case (5)
690                do m = pttruetab(5),cetruetab(5)
691                do l = pttruetab(4),cetruetab(4)
692                do k = pttruetab(3),cetruetab(3)
693                do j = pttruetab(2),cetruetab(2)
694                do i = pttruetab(1),cetruetab(1)
695                    if (restore%restore5D(i,j,k,l,m) == 0)  &
696                        parray5(i,j,k,l,m) = tempC % array5(i,j,k,l,m)
697                enddo
698                enddo
699                enddo
700                enddo
701                enddo
702            case (6)
703                do n = pttruetab(6),cetruetab(6)
704                do m = pttruetab(5),cetruetab(5)
705                do l = pttruetab(4),cetruetab(4)
706                do k = pttruetab(3),cetruetab(3)
707                do j = pttruetab(2),cetruetab(2)
708                do i = pttruetab(1),cetruetab(1)
709                    if (restore%restore6D(i,j,k,l,m,n) == 0)    &
710                        parray6(i,j,k,l,m,n) = tempC % array6(i,j,k,l,m,n)
711                enddo
712                enddo
713                enddo
714                enddo
715                enddo
716                enddo
717            end select
718!
719#endif
720!
721        else    ! .not.to_restore
722!
723
724            if (memberin) then
725    !
726                if ( .not.in_bc ) then
727                    select case(nbdim)
728                    case(1)
729                        call procname(tempC % array1(            &
730                                childarray(1,1,1):childarray(1,2,1)), &
731                                childarray(1,1,2),childarray(1,2,2),.FALSE.,nb,ndir)
732                    case(2)
733                        call procname(                                &
734                                tempC % array2(            &
735                                childarray(1,1,1):childarray(1,2,1),  &
736                                childarray(2,1,1):childarray(2,2,1)), &
737                                childarray(1,1,2),childarray(1,2,2),  &
738                                childarray(2,1,2),childarray(2,2,2),.FALSE.,nb,ndir)
739                    case(3)
740                        call procname(                                &
741                                tempC % array3(            &
742                                childarray(1,1,1):childarray(1,2,1),  &
743                                childarray(2,1,1):childarray(2,2,1),  &
744                                childarray(3,1,1):childarray(3,2,1)), &
745                                childarray(1,1,2),childarray(1,2,2),  &
746                                childarray(2,1,2),childarray(2,2,2),  &
747                                childarray(3,1,2),childarray(3,2,2),.FALSE.,nb,ndir)
748                    case(4)
749                        call procname(                                &
750                                tempC % array4(            &
751                                childarray(1,1,1):childarray(1,2,1),  &
752                                childarray(2,1,1):childarray(2,2,1),  &
753                                childarray(3,1,1):childarray(3,2,1),  &
754                                childarray(4,1,1):childarray(4,2,1)), &
755                                childarray(1,1,2),childarray(1,2,2),  &
756                                childarray(2,1,2),childarray(2,2,2),  &
757                                childarray(3,1,2),childarray(3,2,2),  &
758                                childarray(4,1,2),childarray(4,2,2),.FALSE.,nb,ndir)
759                    case(5)
760                        call procname(                                &
761                                tempC % array5(            &
762                                childarray(1,1,1):childarray(1,2,1),  &
763                                childarray(2,1,1):childarray(2,2,1),  &
764                                childarray(3,1,1):childarray(3,2,1),  &
765                                childarray(4,1,1):childarray(4,2,1),  &
766                                childarray(5,1,1):childarray(5,2,1)), &
767                                childarray(1,1,2),childarray(1,2,2),  &
768                                childarray(2,1,2),childarray(2,2,2),  &
769                                childarray(3,1,2),childarray(3,2,2),  &
770                                childarray(4,1,2),childarray(4,2,2),  &
771                                childarray(5,1,2),childarray(5,2,2),.FALSE.,nb,ndir)
772                    case(6)
773                        call procname(                                &
774                                tempC % array6(            &
775                                childarray(1,1,1):childarray(1,2,1),  &
776                                childarray(2,1,1):childarray(2,2,1),  &
777                                childarray(3,1,1):childarray(3,2,1),  &
778                                childarray(4,1,1):childarray(4,2,1),  &
779                                childarray(5,1,1):childarray(5,2,1),  &
780                                childarray(6,1,1):childarray(6,2,1)), &
781                                childarray(1,1,2),childarray(1,2,2),  &
782                                childarray(2,1,2),childarray(2,2,2),  &
783                                childarray(3,1,2),childarray(3,2,2),  &
784                                childarray(4,1,2),childarray(4,2,2),  &
785                                childarray(5,1,2),childarray(5,2,2),  &
786                                childarray(6,1,2),childarray(6,2,2),.FALSE.,nb,ndir)
787                    end select
788                else    ! we are in_bc
789                    select case (nbdim)
790                    case (1)
791                        parray1(childarray(1,1,2):childarray(1,2,2)) =    &
792                         tempC%array1(childarray(1,1,1):childarray(1,2,1))
793                    case (2)
794                        parray2(childarray(1,1,2):childarray(1,2,2),      &
795                                      childarray(2,1,2):childarray(2,2,2)) =    &
796                         tempC%array2(childarray(1,1,1):childarray(1,2,1),      &
797                                      childarray(2,1,1):childarray(2,2,1))
798                    case (3)
799                        parray3(childarray(1,1,2):childarray(1,2,2),      &
800                                      childarray(2,1,2):childarray(2,2,2),      &
801                                      childarray(3,1,2):childarray(3,2,2)) =    &
802                         tempC%array3(childarray(1,1,1):childarray(1,2,1),      &
803                                      childarray(2,1,1):childarray(2,2,1),      &
804                                      childarray(3,1,1):childarray(3,2,1))
805                    case (4)
806                        parray4(childarray(1,1,2):childarray(1,2,2),      &
807                                      childarray(2,1,2):childarray(2,2,2),      &
808                                      childarray(3,1,2):childarray(3,2,2),      &
809                                      childarray(4,1,2):childarray(4,2,2)) =    &
810                         tempC%array4(childarray(1,1,1):childarray(1,2,1),      &
811                                      childarray(2,1,1):childarray(2,2,1),      &
812                                      childarray(3,1,1):childarray(3,2,1),      &
813                                      childarray(4,1,1):childarray(4,2,1))
814                    case (5)
815                        parray5(childarray(1,1,2):childarray(1,2,2),      &
816                                      childarray(2,1,2):childarray(2,2,2),      &
817                                      childarray(3,1,2):childarray(3,2,2),      &
818                                      childarray(4,1,2):childarray(4,2,2),      &
819                                      childarray(5,1,2):childarray(5,2,2)) =    &
820                         tempC%array5(childarray(1,1,1):childarray(1,2,1),      &
821                                      childarray(2,1,1):childarray(2,2,1),      &
822                                      childarray(3,1,1):childarray(3,2,1),      &
823                                      childarray(4,1,1):childarray(4,2,1),      &
824                                      childarray(5,1,1):childarray(5,2,1))
825                    case (6)
826                        parray6(childarray(1,1,2):childarray(1,2,2),      &
827                                      childarray(2,1,2):childarray(2,2,2),      &
828                                      childarray(3,1,2):childarray(3,2,2),      &
829                                      childarray(4,1,2):childarray(4,2,2),      &
830                                      childarray(5,1,2):childarray(5,2,2),      &
831                                      childarray(6,1,2):childarray(6,2,2)) =    &
832                         tempC%array6(childarray(1,1,1):childarray(1,2,1),      &
833                                      childarray(2,1,1):childarray(2,2,1),      &
834                                      childarray(3,1,1):childarray(3,2,1),      &
835                                      childarray(4,1,1):childarray(4,2,1),      &
836                                      childarray(5,1,1):childarray(5,2,1),      &
837                                      childarray(6,1,1):childarray(6,2,1))
838                    end select
839                endif  ! < (.not.in_bc)
840            endif  ! < memberin
841!
842        endif
843
844
845        call Agrif_array_deallocate(tempPextend,nbdim)
846        call Agrif_array_deallocate(tempC,nbdim)
847
848    endif
849!
850!   Deallocations
851#if defined AGRIF_MPI
852    if (member) then
853        call Agrif_array_deallocate(tempP,nbdim)
854    endif
855#endif
856!---------------------------------------------------------------------------------------------------
857end subroutine Agrif_InterpnD
858!===================================================================================================
859!
860!===================================================================================================
861!  subroutine Agrif_Parentbounds
862!
863!> Calculates the bounds of the parent grid for the interpolation of the child grid
864!---------------------------------------------------------------------------------------------------
865subroutine Agrif_Parentbounds ( type_interp, nbdim, indmin, indmax, &
866                                indmin_required,indmax_required,    &
867                                s_Parent_temp, s_Child_temp,        &
868                                s_Child, ds_Child,                  &
869                                s_Parent,ds_Parent,                 &
870                                pttruetab, cetruetab,               &
871                                pttab_Child, pttab_Parent, posvar, coords )
872!---------------------------------------------------------------------------------------------------
873    INTEGER, DIMENSION(6),     intent(in)  :: type_interp
874    INTEGER,                   intent(in)  :: nbdim
875    INTEGER, DIMENSION(nbdim), intent(out) :: indmin, indmax
876    INTEGER, DIMENSION(nbdim), intent(out) :: indmin_required, indmax_required
877    REAL(kind=8),    DIMENSION(nbdim), intent(out) :: s_Parent_temp, s_child_temp
878    REAL(kind=8),    DIMENSION(nbdim), intent(in)  :: s_Child, ds_child
879    REAL(kind=8),    DIMENSION(nbdim), intent(in)  :: s_Parent,ds_Parent
880    INTEGER, DIMENSION(nbdim), intent(in)  :: pttruetab, cetruetab
881    INTEGER, DIMENSION(nbdim), intent(in)  :: pttab_Child, pttab_Parent
882    INTEGER, DIMENSION(nbdim), intent(in)  :: posvar
883    INTEGER, DIMENSION(nbdim), intent(in)  :: coords
884!
885    REAL(kind=8) :: xpmin, xpmax
886    INTEGER :: coeffraf
887    INTEGER :: i
888    REAL(kind=8),DIMENSION(nbdim) :: dim_newmin, dim_newmax
889!
890    dim_newmin = s_Child + (pttruetab - pttab_Child) * ds_Child
891    dim_newmax = s_Child + (cetruetab - pttab_Child) * ds_Child
892!
893    do i = 1,nbdim
894!
895        indmin(i) = pttab_Parent(i) + agrif_int((dim_newmin(i)-s_Parent(i))/ds_Parent(i))
896        indmax(i) = pttab_Parent(i) + agrif_ceiling((dim_newmax(i)-s_Parent(i))/ds_Parent(i))
897
898        coeffraf = nint(ds_Parent(i)/ds_Child(i))
899       
900        indmin_required(i) = indmin(i)
901        indmax_required(i) = indmax(i)
902!
903!       Necessary for the Quadratic interpolation
904!
905
906        if ( (pttruetab(i) == cetruetab(i)) .and. (posvar(i) == 1) ) then
907            if (Agrif_UseSpecialValue) then
908               indmin(i) = indmin(i)-MaxSearch
909               indmax(i) = indmax(i)+MaxSearch
910            endif
911        elseif ( coords(i) == 0 ) then  ! (interptab == 'N')
912        elseif ( (type_interp(i) == Agrif_ppm)     .or.     &
913                 (type_interp(i) == Agrif_eno)     .or.     &
914                 (type_interp(i) == Agrif_ppm_lim) .or.     &
915                 (type_interp(i) == Agrif_weno) ) then
916                 
917            if ((mod(coeffraf,2) == 0).AND.(posvar(i)==2)) then
918           
919              xpmax = s_Parent(i)+(indmax(i)-pttab_Parent(i))*ds_Parent(i)
920              if (xpmax > dim_newmax(i)+ds_Child(i)) then
921                indmax(i) = indmax(i) + 1
922              else
923                indmax(i) = indmax(i) + 2
924              endif
925             
926              xpmin = s_Parent(i)+(indmin(i)-pttab_Parent(i))*ds_Parent(i)
927              if (xpmin < dim_newmin(i)-ds_Child(i)) then
928                indmin(i) = indmin(i) - 1
929              else
930                indmin(i) = indmin(i) - 2
931              endif
932             
933            else
934              indmin(i) = indmin(i) - 2
935              indmax(i) = indmax(i) + 2
936            endif
937
938            indmin_required(i) = indmin(i)
939            indmax_required(i) = indmax(i)
940       
941            if (Agrif_UseSpecialValue) then
942               indmin(i) = indmin(i)-MaxSearch
943               indmax(i) = indmax(i)+MaxSearch
944            endif
945        elseif (type_interp(i) == Agrif_linearconservlim) then
946       
947            if ((mod(coeffraf,2) == 0).AND.(posvar(i)==2)) then
948           
949              xpmax = s_Parent(i)+(indmax(i)-pttab_Parent(i))*ds_Parent(i)
950              if (xpmax > dim_newmax(i)+ds_Child(i)) then
951                indmax(i) = indmax(i)
952              else
953                indmax(i) = indmax(i) + 1
954              endif
955             
956              xpmin = s_Parent(i)+(indmin(i)-pttab_Parent(i))*ds_Parent(i)
957              if (xpmin < dim_newmin(i)-ds_Child(i)) then
958                indmin(i) = indmin(i)
959              else
960                indmin(i) = indmin(i) - 1
961              endif
962             
963            else
964              indmin(i) = indmin(i) - 1
965              indmax(i) = indmax(i) + 1
966            endif
967
968            indmin_required(i) = indmin(i)
969            indmax_required(i) = indmax(i)
970       
971            if (Agrif_UseSpecialValue) then
972               indmin(i) = indmin(i)-MaxSearch
973               indmax(i) = indmax(i)+MaxSearch
974            endif
975           
976        elseif ( (type_interp(i) /= Agrif_constant) .and.   &
977                 (type_interp(i) /= Agrif_linear) ) then
978            indmin(i) = indmin(i) - 1
979            indmax(i) = indmax(i) + 1
980           
981            indmin_required(i) = indmin(i)
982            indmax_required(i) = indmax(i)
983
984            if (Agrif_UseSpecialValue) then
985               indmin(i) = indmin(i)-MaxSearch
986               indmax(i) = indmax(i)+MaxSearch
987            endif
988        elseif ( (type_interp(i) == Agrif_constant) .or.   &
989                 (type_interp(i) == Agrif_linear) ) then
990            indmin_required(i) = indmin(i)
991            indmax_required(i) = indmax(i)
992            if (Agrif_UseSpecialValue) then
993               indmin(i) = indmin(i)-MaxSearch
994               indmax(i) = indmax(i)+MaxSearch
995            endif
996        endif
997
998!
999    enddo
1000!
1001    s_Parent_temp = s_Parent + (indmin - pttab_Parent) * ds_Parent
1002    s_Child_temp  = s_Child + (pttruetab - pttab_Child) * ds_Child
1003!---------------------------------------------------------------------------------------------------
1004end subroutine Agrif_Parentbounds
1005!===================================================================================================
1006!
1007!===================================================================================================
1008!  subroutine Agrif_Interp_1D_Recursive
1009!
1010!> Subroutine for the interpolation of a 1D grid variable.
1011!> It calls Agrif_InterpBase.
1012!---------------------------------------------------------------------------------------------------
1013subroutine Agrif_Interp_1D_recursive ( type_interp, tabin, tabout,  &
1014                                       indmin, indmax,              &
1015                                       pttab_child, petab_child,    &
1016                                       s_child,  s_parent,          &
1017                                       ds_child, ds_parent )
1018!---------------------------------------------------------------------------------------------------
1019    integer,            intent(in)  :: type_interp
1020    integer,            intent(in)  :: indmin, indmax
1021    integer,            intent(in)  :: pttab_child, petab_child
1022    real(kind=8),               intent(in)  :: s_child, s_parent
1023    real(kind=8),               intent(in)  :: ds_child, ds_parent
1024    real, dimension(            &
1025        indmin:indmax           &
1026    ),                  intent(in)  :: tabin
1027    real, dimension(            &
1028        pttab_child:petab_child &
1029    ),                  intent(out) :: tabout
1030!---------------------------------------------------------------------------------------------------
1031    call Agrif_InterpBase(type_interp,                      &
1032                          tabin(indmin:indmax),             &
1033                          tabout(pttab_child:petab_child),  &
1034                          indmin, indmax,                   &
1035                          pttab_child, petab_child,         &
1036                          s_parent,    s_child,             &
1037                         ds_parent,   ds_child)
1038!---------------------------------------------------------------------------------------------------
1039end subroutine Agrif_Interp_1D_recursive
1040!===================================================================================================
1041!
1042!===================================================================================================
1043!  subroutine Agrif_Interp_2D_Recursive
1044!
1045!> Subroutine for the interpolation of a 2D grid variable.
1046!> It calls Agrif_Interp_1D_recursive and Agrif_InterpBase.
1047!---------------------------------------------------------------------------------------------------
1048subroutine Agrif_Interp_2D_recursive ( type_interp, tabin, tabout,  &
1049                                       indmin, indmax,              &
1050                                       pttab_child, petab_child,    &
1051                                       s_child,  s_parent,          &
1052                                       ds_child, ds_parent )
1053!---------------------------------------------------------------------------------------------------
1054    integer, dimension(2),              intent(in)  :: type_interp
1055    integer, dimension(2),              intent(in)  :: indmin, indmax
1056    integer, dimension(2),              intent(in)  :: pttab_child, petab_child
1057    real(kind=8),    dimension(2),              intent(in)  :: s_child, s_parent
1058    real(kind=8),    dimension(2),              intent(in)  :: ds_child, ds_parent
1059    real,    dimension(                 &
1060        indmin(1):indmax(1),            &
1061        indmin(2):indmax(2)),           intent(in)  :: tabin
1062    real,    dimension(                 &
1063        pttab_child(1):petab_child(1),  &
1064        pttab_child(2):petab_child(2)), intent(out) :: tabout
1065!---------------------------------------------------------------------------------------------------
1066    real, dimension(                    &
1067        pttab_child(1):petab_child(1),  &
1068        indmin(2):indmax(2))            :: tabtemp
1069    real, dimension(                    &
1070        pttab_child(2):petab_child(2),  &
1071        pttab_child(1):petab_child(1))  :: tabout_trsp
1072    real, dimension(                    &
1073        indmin(2):indmax(2),            &
1074        pttab_child(1):petab_child(1))  :: tabtemp_trsp
1075    integer                             :: i, j, coeffraf, locind_child_left, ideb
1076!---------------------------------------------------------------------------------------------------
1077!
1078    coeffraf = nint ( ds_parent(1) / ds_child(1) )
1079!
1080    if ((type_interp(1) == Agrif_Linear) .and. (coeffraf /= 1)) then
1081!---CDIR NEXPAND
1082        if(.NOT. precomputedone(1))     &
1083            call Linear1dPrecompute2d(                  &
1084                    indmax(2)-indmin(2)+1,              &
1085                    indmax(1)-indmin(1)+1,              &
1086                    petab_child(1)-pttab_child(1)+1,    &
1087                    s_parent(1),s_child(1),ds_parent(1),ds_child(1),1)
1088!---CDIR NEXPAND
1089        call Linear1dAfterCompute(tabin,tabtemp,size(tabin),size(tabtemp),1)
1090!
1091    elseif ((type_interp(1) == Agrif_PPM) .and. (coeffraf /= 1)) then
1092!---CDIR NEXPAND
1093        if(.NOT. precomputedone(1))     &
1094            call PPM1dPrecompute2d(                     &
1095                    indmax(2)-indmin(2)+1,              &
1096                    indmax(1)-indmin(1)+1,              &
1097                    petab_child(1)-pttab_child(1)+1,        &
1098                    s_parent(1),s_child(1),ds_parent(1),ds_child(1),1)
1099!---CDIR NEXPAND
1100        call PPM1dAfterCompute(tabin,tabtemp,size(tabin),size(tabtemp),1,indchildppm(:,1),tabppm(:,:,1))
1101    else if (coeffraf == 1) then
1102        locind_child_left = indmin(1) + agrif_int((s_child(1)-s_parent(1))/ds_parent(1))
1103       
1104            do j = indmin(2), indmax(2)
1105            ideb = locind_child_left
1106            do i = pttab_child(1), petab_child(1)
1107                tabtemp(i,j) = tabin(ideb,j)
1108                ideb = ideb + 1
1109            enddo
1110            enddo
1111
1112    else
1113        do j = indmin(2),indmax(2)
1114!
1115!---CDIR NEXPAND
1116            call Agrif_Interp_1D_recursive(type_interp(1),                  &
1117                    tabin(indmin(1):indmax(1),j),               &
1118                    tabtemp(pttab_child(1):petab_child(1),j),   &
1119                    indmin(1),indmax(1),                    &
1120                    pttab_child(1),petab_child(1),          &
1121                    s_child(1), s_parent(1),                &
1122                    ds_child(1),ds_parent(1))
1123!
1124        enddo
1125    endif
1126
1127    coeffraf = nint(ds_parent(2)/ds_child(2))
1128    tabtemp_trsp = TRANSPOSE(tabtemp)
1129
1130    if ((type_interp(2) == Agrif_Linear) .and. (coeffraf /= 1)) then
1131!---CDIR NEXPAND
1132        if(.NOT. precomputedone(2))     &
1133            call Linear1dPrecompute2d(                  &
1134                    petab_child(1)-pttab_child(1)+1,    &
1135                    indmax(2)-indmin(2)+1,              &
1136                    petab_child(2)-pttab_child(2)+1,    &
1137                    s_parent(2),s_child(2),ds_parent(2),ds_child(2),2)
1138!---CDIR NEXPAND
1139        call Linear1dAfterCompute(tabtemp_trsp,tabout_trsp, &
1140                size(tabtemp_trsp),size(tabout_trsp),2)
1141
1142    elseif ((type_interp(2) == Agrif_PPM) .and. (coeffraf /= 1)) then
1143!---CDIR NEXPAND
1144        if(.NOT. precomputedone(2))     &
1145            call PPM1dPrecompute2d(                     &
1146                    petab_child(1)-pttab_child(1)+1,    &
1147                    indmax(2)-indmin(2)+1,              &
1148                    petab_child(2)-pttab_child(2)+1,    &
1149                    s_parent(2),s_child(2),ds_parent(2),ds_child(2),2)
1150!---CDIR NEXPAND
1151        call PPM1dAfterCompute(tabtemp_trsp, tabout_trsp,    &
1152                               size(tabtemp_trsp), size(tabout_trsp), 2, &
1153                              indchildppm(:,2),tabppm(:,:,2))
1154    else
1155        do i = pttab_child(1), petab_child(1)
1156!
1157!---CDIR NEXPAND
1158            call Agrif_InterpBase(type_interp(2),                                   &
1159                                  tabtemp_trsp(indmin(2):indmax(2), i),             &
1160                                  tabout_trsp(pttab_child(2):petab_child(2), i),    &
1161                                  indmin(2), indmax(2),                             &
1162                                  pttab_child(2), petab_child(2),                   &
1163                                  s_parent(2),    s_child(2),                       &
1164                                 ds_parent(2),   ds_child(2) )
1165        enddo
1166    endif
1167!
1168    tabout = TRANSPOSE(tabout_trsp)
1169!---------------------------------------------------------------------------------------------------
1170end subroutine Agrif_Interp_2D_recursive
1171!===================================================================================================
1172!
1173!===================================================================================================
1174!  subroutine Agrif_Interp_3D_Recursive
1175!
1176!> Subroutine for the interpolation of a 3D grid variable.
1177!> It calls #Agrif_Interp_2D_recursive and #Agrif_InterpBase.
1178!---------------------------------------------------------------------------------------------------
1179subroutine Agrif_Interp_3D_recursive ( type_interp, tabin, tabout,  &
1180                                       indmin, indmax,              &
1181                                       pttab_child, petab_child,    &
1182                                       s_child,  s_parent,          &
1183                                       ds_child, ds_parent )
1184!---------------------------------------------------------------------------------------------------
1185    integer, dimension(3),              intent(in)  :: type_interp
1186    integer, dimension(3),              intent(in)  :: indmin, indmax
1187    integer, dimension(3),              intent(in)  :: pttab_child, petab_child
1188    real(kind=8),    dimension(3),              intent(in)  :: s_child, s_parent
1189    real(kind=8),    dimension(3),              intent(in)  :: ds_child, ds_parent
1190    real,    dimension(                 &
1191        indmin(1):indmax(1),            &
1192        indmin(2):indmax(2),            &
1193        indmin(3):indmax(3)),           intent(in)  :: tabin
1194    real,    dimension(                 &
1195        pttab_child(1):petab_child(1),  &
1196        pttab_child(2):petab_child(2),  &
1197        pttab_child(3):petab_child(3)), intent(out) :: tabout
1198!---------------------------------------------------------------------------------------------------
1199    real, dimension(                    &
1200        pttab_child(1):petab_child(1),  &
1201        pttab_child(2):petab_child(2),  &
1202        indmin(3):indmax(3))            :: tabtemp
1203    integer                             :: i, j, k, coeffraf,kp,kp1,kp2,kp3,kp4,kref
1204    integer                             :: locind_child_left, kdeb
1205    real(kind=8)    :: ypos,globind_parent_left
1206    real(kind=8)    :: deltax, invdsparent
1207    real    :: t2,t3,t4,t5,t6,t7,t8
1208    integer :: locind_parent_left
1209
1210!
1211    coeffraf = nint ( ds_parent(1) / ds_child(1) )
1212    if ( (type_interp(1) == Agrif_Linear) .and. (coeffraf/=1) ) then
1213        call Linear1dPrecompute2d(indmax(2)-indmin(2)+1,            &
1214                                  indmax(1)-indmin(1)+1,            &
1215                                  petab_child(1)-pttab_child(1)+1,  &
1216                                  s_parent(1),s_child(1),ds_parent(1),ds_child(1),1)
1217        precomputedone(1) = .TRUE.
1218    elseif ( (type_interp(1) == Agrif_PPM) .and. (coeffraf/=1) ) then
1219        call PPM1dPrecompute2d(indmax(2)-indmin(2)+1,           &
1220                               indmax(1)-indmin(1)+1,           &
1221                               petab_child(1)-pttab_child(1)+1, &
1222                               s_parent(1),s_child(1),ds_parent(1),ds_child(1),1)
1223        precomputedone(1) = .TRUE.
1224    endif
1225
1226    coeffraf = nint ( ds_parent(2) / ds_child(2) )
1227    if ( (type_interp(2) == Agrif_Linear) .and. (coeffraf/=1) ) then
1228        call Linear1dPrecompute2d(petab_child(1)-pttab_child(1)+1,  &
1229                                  indmax(2)-indmin(2)+1,            &
1230                                  petab_child(2)-pttab_child(2)+1,  &
1231                                  s_parent(2),s_child(2),ds_parent(2),ds_child(2),2)
1232        precomputedone(2) = .TRUE.
1233    elseif ( (type_interp(2) == Agrif_PPM) .and. (coeffraf/=1) ) then
1234        call PPM1dPrecompute2d(petab_child(1)-pttab_child(1)+1, &
1235                               indmax(2)-indmin(2)+1,           &
1236                               petab_child(2)-pttab_child(2)+1, &
1237                               s_parent(2),s_child(2),ds_parent(2),ds_child(2),2)
1238        precomputedone(2) = .TRUE.
1239    endif
1240!
1241    do k = indmin(3), indmax(3)
1242        call Agrif_Interp_2D_recursive(type_interp(1:2),                            &
1243                                       tabin(indmin(1):indmax(1),                   &
1244                                             indmin(2):indmax(2), k),               &
1245                                       tabtemp(pttab_child(1):petab_child(1),       &
1246                                               pttab_child(2):petab_child(2), k),   &
1247                                       indmin(1:2), indmax(1:2),                    &
1248                                       pttab_child(1:2), petab_child(1:2),          &
1249                                       s_child(1:2),     s_parent(1:2),             &
1250                                       ds_child(1:2),   ds_parent(1:2) )
1251    enddo
1252!
1253    precomputedone(1) = .FALSE.
1254    precomputedone(2) = .FALSE.
1255    coeffraf = nint(ds_parent(3)/ds_child(3))
1256!
1257    if ( coeffraf == 1 ) then
1258        locind_child_left = 1 + agrif_int((s_child(3)-s_parent(3))/ds_parent(3))
1259        kdeb = indmin(3)+locind_child_left-2
1260        do k = pttab_child(3),petab_child(3)
1261            kdeb = kdeb + 1
1262            do j = pttab_child(2), petab_child(2)
1263            do i = pttab_child(1), petab_child(1)
1264                tabout(i,j,k) = tabtemp(i,j,kdeb)
1265            enddo
1266            enddo
1267        enddo
1268    else if (type_interp(3) == Agrif_LAGRANGE) then
1269      invdsparent = 1./ds_parent(3)
1270      ypos = s_child(3)
1271      do k=pttab_child(3), petab_child(3)
1272        locind_parent_left = indmin(3)+agrif_int((ypos - s_parent(3))/ds_parent(3))
1273        globind_parent_left = s_parent(3) + (locind_parent_left - indmin(3))*ds_parent(3)
1274        deltax = invdsparent*(ypos-globind_parent_left)
1275        deltax = nint(coeffraf*deltax)/real(coeffraf)
1276        ypos = ypos + ds_child(3)
1277
1278        if (abs(deltax) <= 0.0001) then
1279          do j = pttab_child(2), petab_child(2)
1280          do i = pttab_child(1), petab_child(1)
1281             tabout(i,j,k) = tabtemp(i,j,locind_parent_left)
1282          enddo
1283          enddo
1284        else
1285         t2 = deltax - 2.
1286        t3 = deltax - 1.
1287        t4 = deltax + 1.
1288
1289        t5 = -(1./6.)*deltax*t2*t3
1290        t6 = 0.5*t2*t3*t4
1291        t7 = -0.5*deltax*t2*t4
1292        t8 = (1./6.)*deltax*t3*t4     
1293          do j = pttab_child(2), petab_child(2)
1294          do i = pttab_child(1), petab_child(1)
1295             tabout(i,j,k) = t5*tabtemp(i,j,locind_parent_left-1) + t6*tabtemp(i,j,locind_parent_left)    &
1296              +t7*tabtemp(i,j,locind_parent_left+1) + t8*tabtemp(i,j,locind_parent_left+2)
1297          enddo
1298          enddo
1299
1300        endif
1301
1302      enddo
1303    else if (type_interp(3) == Agrif_PPM) then
1304     call PPM1dPrecompute2d(1, &
1305                               indmax(3)-indmin(3)+1,           &
1306                               petab_child(3)-pttab_child(3)+1, &
1307                               s_parent(3),s_child(3),ds_parent(3),ds_child(3),1)
1308
1309     do k=pttab_child(3),petab_child(3)
1310        kref = k-pttab_child(3)+1
1311        kp=indmin(3)+indparentppm(kref,1)-1
1312        kp1 = kp + 1
1313        kp2 = kp1 + 1
1314        kp3 = kp2 + 1
1315        kp4 = kp3 + 1
1316        do j = pttab_child(2), petab_child(2)
1317        do i = pttab_child(1), petab_child(1)
1318         tabout(i,j,k) = tabppm(1,indchildppm(kref,1),1)*tabtemp(i,j,kp)   + &
1319                         tabppm(2,indchildppm(kref,1),1)*tabtemp(i,j,kp1)  + &
1320                         tabppm(3,indchildppm(kref,1),1)*tabtemp(i,j,kp2)  + &
1321                         tabppm(4,indchildppm(kref,1),1)*tabtemp(i,j,kp3)  + &
1322                         tabppm(5,indchildppm(kref,1),1)*tabtemp(i,j,kp4)
1323        enddo
1324        enddo
1325     enddo
1326
1327    else
1328
1329        do j = pttab_child(2), petab_child(2)
1330        do i = pttab_child(1), petab_child(1)
1331            call Agrif_InterpBase(type_interp(3),                                   &
1332                                  tabtemp(i,j,indmin(3):indmax(3)),                 &
1333                                  tabout(i,j,pttab_child(3):petab_child(3)),        &
1334                                  indmin(3), indmax(3),                             &
1335                                  pttab_child(3), petab_child(3),                   &
1336                                  s_parent(3),    s_child(3),                       &
1337                                 ds_parent(3),   ds_child(3) )
1338        enddo
1339        enddo
1340
1341    endif
1342!---------------------------------------------------------------------------------------------------
1343end subroutine Agrif_Interp_3D_recursive
1344!===================================================================================================
1345!
1346!===================================================================================================
1347!  subroutine Agrif_Interp_4D_Recursive
1348!
1349!> Subroutine for the interpolation of a 4D grid variable.
1350!> It calls #Agrif_Interp_3D_recursive and #Agrif_InterpBase.
1351!---------------------------------------------------------------------------------------------------
1352subroutine Agrif_Interp_4D_recursive ( type_interp, tabin, tabout,  &
1353                                       indmin, indmax,              &
1354                                       pttab_child, petab_child,    &
1355                                       s_child,  s_parent,          &
1356                                       ds_child, ds_parent )
1357!---------------------------------------------------------------------------------------------------
1358    integer, dimension(4),              intent(in)  :: type_interp
1359    integer, dimension(4),              intent(in)  :: indmin, indmax
1360    integer, dimension(4),              intent(in)  :: pttab_child, petab_child
1361    real(kind=8),    dimension(4),              intent(in)  :: s_child, s_parent
1362    real(kind=8),    dimension(4),              intent(in)  :: ds_child, ds_parent
1363    real,    dimension(                 &
1364        indmin(1):indmax(1),            &
1365        indmin(2):indmax(2),            &
1366        indmin(3):indmax(3),            &
1367        indmin(4):indmax(4)),           intent(in)  :: tabin
1368    real,    dimension(                 &
1369        pttab_child(1):petab_child(1),  &
1370        pttab_child(2):petab_child(2),  &
1371        pttab_child(3):petab_child(3),  &
1372        pttab_child(4):petab_child(4)), intent(out) :: tabout
1373!---------------------------------------------------------------------------------------------------
1374    real, dimension(                    &
1375        pttab_child(1):petab_child(1),  &
1376        pttab_child(2):petab_child(2),  &
1377        pttab_child(3):petab_child(3),  &
1378        indmin(4):indmax(4))            :: tabtemp
1379    integer                             :: i, j, k, l
1380
1381    real(kind=8)    :: ypos,globind_parent_left
1382    real(kind=8)    :: deltax, invdsparent
1383    real    :: t2,t3,t4,t5,t6,t7,t8
1384    integer :: locind_parent_left, coeffraf
1385!
1386    do l = indmin(4), indmax(4)
1387        call Agrif_Interp_3D_recursive(type_interp(1:3),                            &
1388                                       tabin(indmin(1):indmax(1),                   &
1389                                             indmin(2):indmax(2),                   &
1390                                             indmin(3):indmax(3), l),               &
1391                                       tabtemp(pttab_child(1):petab_child(1),       &
1392                                               pttab_child(2):petab_child(2),       &
1393                                               pttab_child(3):petab_child(3), l),   &
1394                                       indmin(1:3), indmax(1:3),                    &
1395                                       pttab_child(1:3), petab_child(1:3),          &
1396                                       s_child(1:3),    s_parent(1:3),              &
1397                                       ds_child(1:3),  ds_parent(1:3) )
1398    enddo
1399!
1400     if (type_interp(4) == Agrif_LAGRANGE) then
1401      coeffraf = nint(ds_parent(4)/ds_child(4))
1402      invdsparent = 1./ds_parent(4)
1403      ypos = s_child(4)
1404      do l=pttab_child(4), petab_child(4)
1405        locind_parent_left = indmin(4)+agrif_int((ypos - s_parent(4))/ds_parent(4))
1406        globind_parent_left = s_parent(4) + (locind_parent_left - indmin(4))*ds_parent(4)
1407        deltax = invdsparent*(ypos-globind_parent_left)
1408        deltax = nint(coeffraf*deltax)/real(coeffraf)
1409        ypos = ypos + ds_child(4)
1410
1411        if (abs(deltax) <= 0.0001) then
1412          do k = pttab_child(3), petab_child(3)
1413          do j = pttab_child(2), petab_child(2)
1414          do i = pttab_child(1), petab_child(1)
1415             tabout(i,j,k,l) = tabtemp(i,j,k,locind_parent_left)
1416          enddo
1417          enddo
1418          enddo
1419        else
1420         t2 = deltax - 2.
1421        t3 = deltax - 1.
1422        t4 = deltax + 1.
1423
1424        t5 = -(1./6.)*deltax*t2*t3
1425        t6 = 0.5*t2*t3*t4
1426        t7 = -0.5*deltax*t2*t4
1427        t8 = (1./6.)*deltax*t3*t4     
1428          do k = pttab_child(3), petab_child(3)
1429          do j = pttab_child(2), petab_child(2)
1430          do i = pttab_child(1), petab_child(1)
1431             tabout(i,j,k,l) = t5*tabtemp(i,j,k,locind_parent_left-1) + t6*tabtemp(i,j,k,locind_parent_left)    &
1432              +t7*tabtemp(i,j,k,locind_parent_left+1) + t8*tabtemp(i,j,k,locind_parent_left+2)
1433          enddo
1434          enddo
1435          enddo
1436        endif
1437
1438      enddo
1439    else
1440    do k = pttab_child(3), petab_child(3)
1441    do j = pttab_child(2), petab_child(2)
1442    do i = pttab_child(1), petab_child(1)
1443        call Agrif_InterpBase(type_interp(4),                                       &
1444                              tabtemp(i,j,k,indmin(4):indmax(4)),                   &
1445                              tabout(i,j,k,pttab_child(4):petab_child(4)),          &
1446                              indmin(4), indmax(4),                                 &
1447                              pttab_child(4), petab_child(4),                       &
1448                              s_parent(4),    s_child(4),                           &
1449                             ds_parent(4),   ds_child(4) )
1450    enddo
1451    enddo
1452    enddo
1453    endif
1454!---------------------------------------------------------------------------------------------------
1455end subroutine Agrif_Interp_4D_recursive
1456!===================================================================================================
1457!
1458!===================================================================================================
1459!  subroutine Agrif_Interp_5D_Recursive
1460!
1461!> Subroutine for the interpolation of a 5D grid variable.
1462!> It calls #Agrif_Interp_4D_recursive and #Agrif_InterpBase.
1463!---------------------------------------------------------------------------------------------------
1464subroutine Agrif_Interp_5D_recursive ( type_interp, tabin, tabout,  &
1465                                       indmin, indmax,              &
1466                                       pttab_child, petab_child,    &
1467                                       s_child,  s_parent,          &
1468                                       ds_child, ds_parent )
1469!---------------------------------------------------------------------------------------------------
1470    integer, dimension(5),              intent(in)  :: type_interp
1471    integer, dimension(5),              intent(in)  :: indmin, indmax
1472    integer, dimension(5),              intent(in)  :: pttab_child, petab_child
1473    real(kind=8),    dimension(5),              intent(in)  :: s_child, s_parent
1474    real(kind=8),    dimension(5),              intent(in)  :: ds_child, ds_parent
1475    real,    dimension(                 &
1476        indmin(1):indmax(1),            &
1477        indmin(2):indmax(2),            &
1478        indmin(3):indmax(3),            &
1479        indmin(4):indmax(4),            &
1480        indmin(5):indmax(5)),           intent(in)  :: tabin
1481    real,    dimension(                 &
1482        pttab_child(1):petab_child(1),  &
1483        pttab_child(2):petab_child(2),  &
1484        pttab_child(3):petab_child(3),  &
1485        pttab_child(4):petab_child(4),  &
1486        pttab_child(5):petab_child(5)), intent(out) :: tabout
1487!---------------------------------------------------------------------------------------------------
1488    real, dimension(                    &
1489        pttab_child(1):petab_child(1),  &
1490        pttab_child(2):petab_child(2),  &
1491        pttab_child(3):petab_child(3),  &
1492        pttab_child(4):petab_child(4),  &
1493        indmin(5):indmax(5))            :: tabtemp
1494    integer                             :: i, j, k, l, m
1495!
1496    do m = indmin(5), indmax(5)
1497        call Agrif_Interp_4D_recursive(type_interp(1:4),                            &
1498                                       tabin(indmin(1):indmax(1),                   &
1499                                             indmin(2):indmax(2),                   &
1500                                             indmin(3):indmax(3),                   &
1501                                             indmin(4):indmax(4),m),                &
1502                                       tabtemp(pttab_child(1):petab_child(1),       &
1503                                               pttab_child(2):petab_child(2),       &
1504                                               pttab_child(3):petab_child(3),       &
1505                                               pttab_child(4):petab_child(4), m),   &
1506                                       indmin(1:4),indmax(1:4),                     &
1507                                       pttab_child(1:4), petab_child(1:4),          &
1508                                       s_child(1:4),     s_parent(1:4),             &
1509                                       ds_child(1:4),   ds_parent(1:4) )
1510    enddo
1511!
1512    do l = pttab_child(4), petab_child(4)
1513    do k = pttab_child(3), petab_child(3)
1514    do j = pttab_child(2), petab_child(2)
1515    do i = pttab_child(1), petab_child(1)
1516        call Agrif_InterpBase(type_interp(5),                                       &
1517                              tabtemp(i,j,k,l,indmin(5):indmax(5)),                 &
1518                              tabout(i,j,k,l,pttab_child(5):petab_child(5)),        &
1519                              indmin(5), indmax(5),                                 &
1520                              pttab_child(5), petab_child(5),                       &
1521                              s_parent(5),   s_child(5),                            &
1522                              ds_parent(5), ds_child(5) )
1523    enddo
1524    enddo
1525    enddo
1526    enddo
1527!---------------------------------------------------------------------------------------------------
1528end subroutine Agrif_Interp_5D_recursive
1529!===================================================================================================
1530!
1531!===================================================================================================
1532!  subroutine Agrif_Interp_6D_Recursive
1533!
1534!> Subroutine for the interpolation of a 6D grid variable.
1535!> It calls #Agrif_Interp_5D_recursive and Agrif_InterpBase.
1536!---------------------------------------------------------------------------------------------------
1537subroutine Agrif_Interp_6D_recursive ( type_interp, tabin, tabout,  &
1538                                       indmin, indmax,              &
1539                                       pttab_child, petab_child,    &
1540                                       s_child,  s_parent,          &
1541                                       ds_child, ds_parent )
1542!---------------------------------------------------------------------------------------------------
1543    integer, dimension(6),                  intent(in)  :: type_interp
1544    integer, dimension(6),                  intent(in)  :: indmin, indmax
1545    integer, dimension(6),                  intent(in)  :: pttab_child, petab_child
1546    real(kind=8),    dimension(6),                  intent(in)  :: s_child, s_parent
1547    real(kind=8),    dimension(6),                  intent(in)  :: ds_child, ds_parent
1548    real,    dimension(                 &
1549        indmin(1):indmax(1),            &
1550        indmin(2):indmax(2),            &
1551        indmin(3):indmax(3),            &
1552        indmin(4):indmax(4),            &
1553        indmin(5):indmax(5),            &
1554        indmin(6):indmax(6)),               intent(in)  :: tabin
1555    real,    dimension(                 &
1556        pttab_child(1):petab_child(1),  &
1557        pttab_child(2):petab_child(2),  &
1558        pttab_child(3):petab_child(3),  &
1559        pttab_child(4):petab_child(4),  &
1560        pttab_child(5):petab_child(5),  &
1561        pttab_child(6):petab_child(6)),     intent(out) :: tabout
1562!---------------------------------------------------------------------------------------------------
1563    real, dimension(                    &
1564        pttab_child(1):petab_child(1),  &
1565        pttab_child(2):petab_child(2),  &
1566        pttab_child(3):petab_child(3),  &
1567        pttab_child(4):petab_child(4),  &
1568        pttab_child(5):petab_child(5),  &
1569        indmin(6):indmax(6))            :: tabtemp
1570    integer                             :: i, j, k, l, m, n
1571!
1572    do n = indmin(6), indmax(6)
1573        call Agrif_Interp_5D_recursive(type_interp(1:5),                            &
1574                                       tabin(indmin(1):indmax(1),                   &
1575                                             indmin(2):indmax(2),                   &
1576                                             indmin(3):indmax(3),                   &
1577                                             indmin(4):indmax(4),                   &
1578                                             indmin(5):indmax(5), n),               &
1579                                        tabtemp(pttab_child(1):petab_child(1),      &
1580                                                pttab_child(2):petab_child(2),      &
1581                                                pttab_child(3):petab_child(3),      &
1582                                                pttab_child(4):petab_child(4),      &
1583                                                pttab_child(5):petab_child(5), n),  &
1584                                        indmin(1:5),indmax(1:5),                    &
1585                                        pttab_child(1:5), petab_child(1:5),         &
1586                                        s_child(1:5), s_parent(1:5),                &
1587                                        ds_child(1:5),ds_parent(1:5) )
1588    enddo
1589!
1590    do m = pttab_child(5), petab_child(5)
1591    do l = pttab_child(4), petab_child(4)
1592    do k = pttab_child(3), petab_child(3)
1593    do j = pttab_child(2), petab_child(2)
1594    do i = pttab_child(1), petab_child(1)
1595        call Agrif_InterpBase(type_interp(6),                                       &
1596                              tabtemp(i,j,k,l,m,indmin(6):indmax(6)),               &
1597                              tabout(i,j,k,l,m,pttab_child(6):petab_child(6)),      &
1598                              indmin(6), indmax(6),                                 &
1599                              pttab_child(6), petab_child(6),                       &
1600                              s_parent(6),   s_child(6),                            &
1601                              ds_parent(6), ds_child(6) )
1602    enddo
1603    enddo
1604    enddo
1605    enddo
1606    enddo
1607!---------------------------------------------------------------------------------------------------
1608end subroutine Agrif_Interp_6D_recursive
1609!===================================================================================================
1610!
1611!===================================================================================================
1612!  subroutine Agrif_InterpBase
1613!
1614!> Calls the interpolation method chosen by the user (linear, lagrange, spline, etc.).
1615!---------------------------------------------------------------------------------------------------
1616subroutine Agrif_InterpBase ( type_interp, parenttab, childtab, indmin, indmax, &
1617                              pttab_child, petab_child,                         &
1618                              s_parent, s_child, ds_parent, ds_child )
1619!---------------------------------------------------------------------------------------------------
1620    INTEGER                                                 :: type_interp
1621    INTEGER                                                 :: indmin, indmax
1622    INTEGER                                                 :: pttab_child, petab_child
1623    REAL, DIMENSION(indmin:indmax),           INTENT(IN)    :: parenttab
1624    REAL, DIMENSION(pttab_child:petab_child), INTENT(OUT)   :: childtab
1625    REAL(kind=8)                                            :: s_parent, s_child
1626    REAL(kind=8)                                            :: ds_parent,ds_child
1627!
1628    if ( (indmin == indmax) .and. (pttab_child == petab_child) ) then
1629!
1630        childtab(pttab_child) = parenttab(indmin)
1631!
1632    elseif (type_interp == Agrif_LINEAR) then    !       Linear interpolation
1633!
1634        call Agrif_basicinterp_linear1D(parenttab,childtab,                   &
1635                indmax-indmin+1,petab_child-pttab_child+1,  &
1636                s_parent,s_child,ds_parent,ds_child)
1637!
1638    elseif ( type_interp == Agrif_PPM ) then     !       PPM interpolation
1639
1640        call PPM1d(parenttab,childtab,                      &
1641                indmax-indmin+1,petab_child-pttab_child+1,  &
1642                s_parent,s_child,ds_parent,ds_child)
1643!
1644    elseif ( type_interp == Agrif_PPM_LIM ) then     !       PPM interpolation
1645
1646        call PPM1d_lim(parenttab,childtab,                      &
1647                indmax-indmin+1,petab_child-pttab_child+1,  &
1648                s_parent,s_child,ds_parent,ds_child)
1649!
1650    elseif (type_interp == Agrif_LAGRANGE) then  !       Lagrange interpolation
1651!
1652        call lagrange1D(parenttab,childtab,                 &
1653                indmax-indmin+1,petab_child-pttab_child+1,  &
1654                s_parent,s_child,ds_parent,ds_child)
1655!
1656    elseif (type_interp == Agrif_ENO) then       !       Eno interpolation
1657!
1658        call ENO1d(parenttab,childtab,                      &
1659                indmax-indmin+1,petab_child-pttab_child+1,  &
1660                s_parent,s_child,ds_parent,ds_child)
1661!
1662    elseif (type_interp == Agrif_WENO) then      !       Weno interpolation
1663!
1664        call WENO1d(parenttab,childtab,                     &
1665                indmax-indmin+1,petab_child-pttab_child+1,  &
1666                s_parent,s_child,ds_parent,ds_child)
1667!
1668    elseif (type_interp == Agrif_LINEARCONSERV) then !   Linear conservative interpolation
1669!
1670        call Linear1dConserv(parenttab,childtab,            &
1671                indmax-indmin+1,petab_child-pttab_child+1,  &
1672                s_parent,s_child,ds_parent,ds_child)
1673!
1674    elseif (type_interp == Agrif_LINEARCONSERVLIM) then !Linear conservative interpolation
1675!
1676        call Linear1dConservLim(parenttab,childtab,         &
1677                indmax-indmin+1,petab_child-pttab_child+1,  &
1678                s_parent,s_child,ds_parent,ds_child)
1679!
1680    elseif (type_interp == Agrif_CONSTANT) then
1681!
1682        call Constant1d(parenttab,childtab,                 &
1683                indmax-indmin+1,petab_child-pttab_child+1,  &
1684                s_parent,s_child,ds_parent,ds_child)
1685!
1686    endif
1687!---------------------------------------------------------------------------------------------------
1688end subroutine Agrif_InterpBase
1689!===================================================================================================
1690!
1691!===================================================================================================
1692!  subroutine Agrif_Find_list_interp
1693!---------------------------------------------------------------------------------------------------
1694function Agrif_Find_list_interp ( list_interp, pttab, petab, pttab_Child, pttab_Parent,     &
1695                                    nbdim, indmin, indmax, indmin_required_p, indmax_required_p, &
1696                                    indminglob,  indmaxglob,         &
1697                                    pttruetab, cetruetab, memberin                          &
1698#if defined AGRIF_MPI
1699                                   ,indminglob2, indmaxglob2, parentarray,                  &
1700                                    member, tab4t, memberinall, sendtoproc1, recvfromproc1  &
1701#endif
1702    ) result(find_list_interp)
1703!---------------------------------------------------------------------------------------------------
1704    type(Agrif_List_Interp_Loc),   pointer     :: list_interp
1705    integer,                       intent(in)  :: nbdim
1706    integer, dimension(nbdim),     intent(in)  :: pttab, petab, pttab_Child, pttab_Parent
1707    integer, dimension(nbdim),     intent(out) :: indmin, indmax, indmin_required_p, indmax_required_p
1708    integer, dimension(nbdim),     intent(out) :: indminglob, indmaxglob
1709    integer, dimension(nbdim),     intent(out) :: pttruetab, cetruetab
1710    logical,                       intent(out) :: memberin
1711#if defined AGRIF_MPI
1712    integer, dimension(nbdim),     intent(out) :: indminglob2, indmaxglob2
1713    integer, dimension(nbdim,2,2), intent(out) :: parentarray
1714    logical,                       intent(out) :: member
1715    integer, dimension(nbdim,0:Agrif_Nbprocs-1,8), intent(out) :: tab4t
1716    logical, dimension(0:Agrif_Nbprocs-1),         intent(out) :: memberinall
1717    logical, dimension(0:Agrif_Nbprocs-1),         intent(out) :: sendtoproc1, recvfromproc1
1718#endif
1719    logical :: find_list_interp
1720!
1721    integer :: i
1722    type(Agrif_List_Interp_Loc), pointer :: parcours
1723    type(Agrif_Interp_Loc),      pointer :: pil
1724
1725    find_list_interp = .false.
1726
1727    if ( .not. associated(list_interp) )    return
1728
1729    parcours => list_interp
1730    find_loop : do while ( associated(parcours) )
1731
1732        pil => parcours % interp_loc
1733
1734        do i = 1,nbdim
1735            if ( (pttab(i) /= pil % pttab(i)) .or. &
1736                 (petab(i) /= pil % petab(i)) .or. &
1737                 (pttab_child(i)  /= pil % pttab_child(i)) .or. &
1738                 (pttab_parent(i) /= pil % pttab_parent(i)) ) then
1739                parcours => parcours % suiv
1740                cycle find_loop
1741            endif
1742        enddo
1743
1744        indmin = pil % indmin(1:nbdim)
1745        indmax = pil % indmax(1:nbdim)
1746        indmin_required_p = pil % indmin_required_p(1:nbdim)
1747        indmax_required_p = pil % indmax_required_p(1:nbdim)
1748
1749        pttruetab = pil % pttruetab(1:nbdim)
1750        cetruetab = pil % cetruetab(1:nbdim)
1751
1752#if !defined AGRIF_MPI
1753        indminglob  = pil % indminglob(1:nbdim)
1754        indmaxglob  = pil % indmaxglob(1:nbdim)
1755#else
1756        indminglob  = pil % indminglob2(1:nbdim)
1757        indmaxglob  = pil % indmaxglob2(1:nbdim)
1758        indminglob2 = pil % indminglob2(1:nbdim)
1759        indmaxglob2 = pil % indmaxglob2(1:nbdim)
1760        parentarray = pil % parentarray(1:nbdim,:,:)
1761        member      = pil % member
1762        tab4t         = pil % tab4t(1:nbdim, 0:Agrif_Nbprocs-1, 1:8)
1763        memberinall   = pil % memberinall(0:Agrif_Nbprocs-1)
1764        sendtoproc1   = pil % sendtoproc1(0:Agrif_Nbprocs-1)
1765        recvfromproc1 = pil % recvfromproc1(0:Agrif_Nbprocs-1)
1766#endif
1767        memberin = pil % memberin
1768        find_list_interp = .true.
1769        exit find_loop
1770    enddo find_loop
1771!---------------------------------------------------------------------------------------------------
1772end function Agrif_Find_list_interp
1773!===================================================================================================
1774!
1775!===================================================================================================
1776!  subroutine Agrif_AddTo_list_interp
1777!---------------------------------------------------------------------------------------------------
1778subroutine Agrif_AddTo_list_interp ( list_interp, pttab, petab, pttab_Child, pttab_Parent,  &
1779                                     indmin, indmax, indmin_required_p, indmax_required_p,  &
1780                                     indminglob, indmaxglob,                &
1781                                     pttruetab, cetruetab,                                  &
1782                                     memberin, nbdim                                        &
1783#if defined AGRIF_MPI
1784                                    ,indminglob2, indmaxglob2,                              &
1785                                     parentarray,                                           &
1786                                     member,                                                &
1787                                     tab4t, memberinall, sendtoproc1, recvfromproc1         &
1788#endif
1789    )
1790!---------------------------------------------------------------------------------------------------
1791    type(Agrif_List_Interp_Loc), pointer    :: list_interp
1792    integer                                 :: nbdim
1793    integer, dimension(nbdim)               :: pttab, petab, pttab_Child, pttab_Parent
1794    integer, dimension(nbdim)               :: indmin,indmax, indmin_required_p, indmax_required_p
1795    integer, dimension(nbdim)               :: indminglob, indmaxglob
1796    integer, dimension(nbdim)               :: pttruetab, cetruetab
1797    logical                                 :: memberin
1798#if defined AGRIF_MPI
1799    integer, dimension(nbdim,2,2)           :: parentarray
1800    logical                                 :: member
1801    integer, dimension(nbdim)                       :: indminglob2,indmaxglob2
1802    integer, dimension(nbdim,0:Agrif_Nbprocs-1,8)   :: tab4t
1803    logical, dimension(0:Agrif_Nbprocs-1)           :: memberinall
1804    logical, dimension(0:Agrif_Nbprocs-1)           :: sendtoproc1
1805    logical, dimension(0:Agrif_Nbprocs-1)           :: recvfromproc1
1806#endif
1807!
1808    type(Agrif_List_Interp_Loc), pointer    :: parcours
1809    type(Agrif_Interp_Loc),      pointer    :: pil
1810!
1811    allocate(parcours)
1812    allocate(parcours % interp_loc)
1813
1814    pil => parcours % interp_loc
1815
1816    pil % pttab(1:nbdim) = pttab(1:nbdim)
1817    pil % petab(1:nbdim) = petab(1:nbdim)
1818    pil % pttab_child(1:nbdim) = pttab_child(1:nbdim)
1819    pil % pttab_parent(1:nbdim) = pttab_parent(1:nbdim)
1820
1821    pil % indmin(1:nbdim) = indmin(1:nbdim)
1822    pil % indmax(1:nbdim) = indmax(1:nbdim)
1823
1824    pil % indmin_required_p(1:nbdim) = indmin_required_p(1:nbdim)
1825    pil % indmax_required_p(1:nbdim) = indmax_required_p(1:nbdim)
1826
1827    pil % memberin = memberin
1828#if !defined AGRIF_MPI
1829    pil % indminglob(1:nbdim) = indminglob(1:nbdim)
1830    pil % indmaxglob(1:nbdim) = indmaxglob(1:nbdim)
1831#else
1832    pil % indminglob2(1:nbdim) = indminglob2(1:nbdim)
1833    pil % indmaxglob2(1:nbdim) = indmaxglob2(1:nbdim)
1834    pil % parentarray(1:nbdim,:,:) = parentarray(1:nbdim,:,:)
1835    pil % member = member
1836    allocate(pil % tab4t(nbdim, 0:Agrif_Nbprocs-1, 8))
1837    allocate(pil % memberinall(0:Agrif_Nbprocs-1))
1838    allocate(pil % sendtoproc1(0:Agrif_Nbprocs-1))
1839    allocate(pil % recvfromproc1(0:Agrif_Nbprocs-1))
1840    pil % tab4t         = tab4t
1841    pil % memberinall   = memberinall
1842    pil % sendtoproc1   = sendtoproc1
1843    pil % recvfromproc1 = recvfromproc1
1844#endif
1845
1846    pil % pttruetab(1:nbdim) = pttruetab(1:nbdim)
1847    pil % cetruetab(1:nbdim) = cetruetab(1:nbdim)
1848
1849    parcours % suiv => list_interp
1850    list_interp => parcours
1851!---------------------------------------------------------------------------------------------------
1852end subroutine Agrif_Addto_list_interp
1853!===================================================================================================
1854!
1855end module Agrif_Interpolation
Note: See TracBrowser for help on using the repository browser.