From 6b537abbbddc03ed6798d5d9c87ce71a8332f309 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Lagaert <Jean-Baptiste.Lagaert@math.u-psud.fr>
Date: Wed, 12 Oct 2011 13:07:43 +0000
Subject: [PATCH] =?UTF-8?q?enfin=20un=20commit=20correct=20du=20code=20du?=
=?UTF-8?q?=20LEGI.=20Bon,=20j'arr=C3=AAte=20de=20commiter=20toutes=20les?=
=?UTF-8?q?=202=20secondes,=20et=20je=20vais=20valider=20et=20ajouter=20l'?=
=?UTF-8?q?ordre=204?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
LEGI/AutresRemaillage/advec_O2SM_x.f90 | 536 +++++++++++++++++++++
LEGI/AutresRemaillage/advec_O2SM_y.f90 | 336 +++++++++++++
LEGI/AutresRemaillage/advec_O2SM_z.f90 | 335 +++++++++++++
LEGI/AutresRemaillage/advec_O4_x.f90 | 541 +++++++++++++++++++++
LEGI/AutresRemaillage/advec_O4_y.f90 | 340 +++++++++++++
LEGI/AutresRemaillage/advec_O4_z.f90 | 339 +++++++++++++
LEGI/AutresRemaillage/advec_O4v2_x.f90 | 640 +++++++++++++++++++++++++
LEGI/AutresRemaillage/advec_O4v2_y.f90 | 445 +++++++++++++++++
LEGI/AutresRemaillage/advec_O4v2_z.f90 | 447 +++++++++++++++++
LEGI/advec.F90 | 286 +++++++++++
LEGI/advecVAR.F90 | 36 ++
LEGI/advecX.F90 | 527 ++++++++++++++++++++
LEGI/advecY.F90 | 529 ++++++++++++++++++++
LEGI/advecZ.F90 | 514 ++++++++++++++++++++
14 files changed, 5851 insertions(+)
create mode 100644 LEGI/AutresRemaillage/advec_O2SM_x.f90
create mode 100644 LEGI/AutresRemaillage/advec_O2SM_y.f90
create mode 100644 LEGI/AutresRemaillage/advec_O2SM_z.f90
create mode 100644 LEGI/AutresRemaillage/advec_O4_x.f90
create mode 100644 LEGI/AutresRemaillage/advec_O4_y.f90
create mode 100644 LEGI/AutresRemaillage/advec_O4_z.f90
create mode 100644 LEGI/AutresRemaillage/advec_O4v2_x.f90
create mode 100644 LEGI/AutresRemaillage/advec_O4v2_y.f90
create mode 100644 LEGI/AutresRemaillage/advec_O4v2_z.f90
create mode 100644 LEGI/advec.F90
create mode 100644 LEGI/advecVAR.F90
create mode 100644 LEGI/advecX.F90
create mode 100644 LEGI/advecY.F90
create mode 100644 LEGI/advecZ.F90
diff --git a/LEGI/AutresRemaillage/advec_O2SM_x.f90 b/LEGI/AutresRemaillage/advec_O2SM_x.f90
new file mode 100644
index 000000000..634703862
--- /dev/null
+++ b/LEGI/AutresRemaillage/advec_O2SM_x.f90
@@ -0,0 +1,536 @@
+module x_advec_m
+
+ use solver
+ use data
+
+ implicit none
+
+ integer :: npart, ntag_total, npg
+ real(WP) :: xmin, ymin, zmin, xmax, ymax, zmax, cfl
+ integer, dimension(:), allocatable :: itype
+ integer, dimension(:), allocatable :: itype_aux
+ integer, dimension(:), allocatable :: itag
+ integer, dimension(:), allocatable :: icfl
+ real(WP), dimension(:), allocatable :: up
+ real(WP), dimension(:), allocatable :: up_aux
+ real(WP), dimension(:), allocatable :: xp
+ real(WP), dimension(:), allocatable :: xp0
+ real(WP), dimension(:), allocatable :: xp_aux
+ real(WP), dimension(:), allocatable :: vx
+ real(WP), dimension(:), allocatable :: vx_aux
+
+ real(WP) :: circlim
+
+end module x_advec_m
+
+
+
+subroutine x_advec_init
+
+ use x_advec_m
+
+ implicit none
+
+ npg = nxsc !100000
+
+ allocate(up(npg))
+ allocate(up_aux(npg))
+ allocate(xp(npg)) !!-> code original npg = 256 mais ca marche pas
+ allocate(xp0(npg))
+ allocate(xp_aux(npg))
+ allocate(vx(npg))
+ allocate(vx_aux(npg))
+ allocate(itype(npg))
+ allocate(itype_aux(npg))
+ allocate(itag(npg))
+ allocate(icfl(npg))
+
+ up = 0.0
+ up_aux = 0.0
+ xp = 0.0
+ xp0 = 0.0
+ xp_aux = 0.0
+ vx = 0.0
+ vx_aux = 0.0
+ cfl = 0.0
+
+ itype = 0
+ itype_aux = 0
+ itag = 0
+ icfl = 0
+
+ npart = 0
+
+ xmin = 0.
+ ymin = 0.
+ zmin = 0.
+ xmax = L
+ ymax = L
+ zmax = L
+
+ circlim = -1 !10**(-5)
+
+end subroutine x_advec_init
+
+
+
+subroutine x_advec
+
+ use x_advec_m
+
+ implicit none
+
+ real(WP) :: dx, yy, zz
+ integer :: i,j,k,np,n, jr, kr, np_aux, ini, ntag
+
+ dx = L/nxsc
+ cfl = dt/dx
+
+ npart = 0
+ ntag_total = 0
+ ntag = 0
+
+ do k=1,nxsc
+ zz=xmin+float(k-1)*dx
+ do j=1,nxsc
+ np=0
+ yy=xmin+float(j-1)*dx
+ do i=1,nxsc
+ if (abs(SC(i,j,k)).gt.circlim) then
+ np=np+1
+ up(np)=SC(i,j,k)
+ xp(np)=xmin+float(i-1)*dx
+ end if
+ end do
+ if (np.ne.0) then
+ call velox_x(np,j,k)
+ do n=1,np
+ xp0(n)=xp(n)
+ xp(n)=xp(n)+0.5*dt*vx(n)
+ if (xp(n).gt.xmax) xp(n)=xp(n)-xmax+xmin
+ if (xp(n).lt.xmin) xp(n)=xp(n)+xmax-xmin
+ end do
+ call velox_x(np,j,k)
+
+!! debut modif
+
+
+ call tag_particles(np,np_aux,ntag)
+ ntag_total = ntag_total + ntag
+ ! test
+ !itype(1:np_aux) = 1
+
+! np integer: deja defini (input) - local
+! np_aux integer: deja defini (output) - local
+! ntag integer: deja defini (output) - local
+! np_bl=1 --> supprimer de l'appel
+! icfl tablo integer: deja defini - common --> supprimer de l'appel
+! itype tablo integer: deja defini - common --> supprimer de l'appel
+! itype_aux tablo integer: deja defini - common --> supprimer de l'appel
+! itag tablo integer: deja defini - common --> supprimer de l'appel
+! xp_aux tablo real integer: deja defini - common --> supprimer de l'appel
+! up_aux tablo real integer: deja defini - common --> supprimer de l'appel
+! vx_aux tablo real integer: deja defini - common --> supprimer de l'appel
+
+! ---- > trouver les dimension max des tablo
+! ---- > trouver si les tablo doivent etre declarer en local ou en common
+
+ if (ntag.ne.0) then
+ do n=1,ntag
+ ini=itag(n)
+ xp(ini)=xp0(ini)+dt*vx(ini)
+ if (xp(ini).gt.xmax) xp(ini)=xp(ini)-xmax+xmin
+ if (xp(ini).lt.xmin) xp(ini)=xp(ini)+xmax-xmin !! TO CHECK WITH GH
+ end do
+ end if
+
+ do n=1,np_aux
+ xp_aux(n)=xp_aux(n)+dt*vx_aux(n)
+ if (xp_aux(n).gt.xmax) xp_aux(n)=xp_aux(n)-xmax+xmin
+ if (xp_aux(n).lt.xmin) xp_aux(n)=xp_aux(n)+xmax-xmin
+ end do
+ jr=j
+ kr=k
+ call remeshx(np_aux,jr,kr)
+ if (ntag.ne.0) call remeshx_tag(ntag,jr,kr)
+
+!! fin modif
+
+ npart=npart+np
+ end if
+ end do
+ end do
+
+ print*, 'NPART, NTAG ', npart,ntag_total
+ print*,'ntag x_advec :', ntag_total
+
+
+end subroutine x_advec
+
+
+
+subroutine velox_x(np,j,k)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: np, j, k, i, nx3, ny3, nz3, jp1, jp2, kp1, kp2, ip1, ip2
+ real(WP) :: dx3, dy3, dz3, dxinv, yy, zz, x0, y0, z0, dx2
+ real(WP) :: yy1, yy2, zz1, zz2, b1, b2, c1, c2, a1, a2
+ real(WP) :: x, xx1, xx2
+
+ do i=1,np
+ vx(i)=0.
+ end do
+
+ nx3=nxsc
+ ny3=nx3
+ nz3=nx3
+ dx3=xmax/float(nx3)
+ dy3=dx3
+ dz3=dx3
+
+ dxinv=1./dx3
+
+ x0=xmin
+ y0=ymin
+ z0=zmin
+
+ dx2 = L/nxsc
+
+ yy=+float(j-1)*dx2
+ zz=+float(k-1)*dx2
+
+ jp1 = int((yy)/dx3)
+ jp2 = jp1 + 1
+ kp1 = int((zz)/dx3)
+ kp2 = kp1 + 1
+ yy1 = (yy-float(jp1)*dx3)/dx3
+ yy2 = 1.0 - yy1
+ zz1 = (zz-float(kp1)*dx3)/dx3
+ zz2 = 1.0 - zz1
+ b1=yy2
+ b2=yy1
+ c1=zz2
+ c2=zz1
+
+ jp1=mod(jp1+nx3,nx3) +1
+ jp2=mod(jp2+nx3,nx3) +1
+ kp1=mod(kp1+nx3,nx3) +1
+ kp2=mod(kp2+nx3,nx3) +1
+
+ do i = 1,np
+
+ x = xp(i)
+
+ ip1 = int((x-x0)*dxinv)
+ ip2 = ip1 + 1
+
+ xx1 = (x - float(ip1)*dx3-x0)*dxinv
+ xx2=1-xx1
+
+ ip1=mod(ip1+nx3,nx3) +1
+ ip2=mod(ip2+nx3,nx3) +1
+
+ a1 = xx2
+ a2 = xx1
+
+ vx(i)= vx(i) + Usc(ip1,jp1,kp1)*a1*b1*c1
+ vx(i)= vx(i) + Usc(ip2,jp1,kp1)*a2*b1*c1
+ vx(i)= vx(i) + Usc(ip1,jp2,kp1)*a1*b2*c1
+ vx(i)= vx(i) + Usc(ip2,jp2,kp1)*a2*b2*c1
+ vx(i)= vx(i) + Usc(ip1,jp1,kp2)*a1*b1*c2
+ vx(i)= vx(i) + Usc(ip2,jp1,kp2)*a2*b1*c2
+ vx(i)= vx(i) + Usc(ip1,jp2,kp2)*a1*b2*c2
+ vx(i)= vx(i) + Usc(ip2,jp2,kp2)*a2*b2*c2
+
+ end do
+
+end subroutine velox_x
+
+subroutine remeshx(np1,jj,kk)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: np1, jj, kk
+ integer :: i, nx2, n, ip0, ip1, ip2, k, j, ip3
+ real(WP) :: dxinv, dx2, x0, g1, x, xx0, xx1, xx2, a0, a1, a2, a3, xx3
+
+ nx2 = nxsc
+ dx2 = L/nx2
+
+ dxinv=1./dx2
+
+ do i=1,nx2
+ SC(i,jj,kk)=0.
+ enddo
+
+ x0=xmin
+
+ do n = 1,np1
+ g1 = up_aux(n)
+ x = xp_aux(n)
+
+ if (itype_aux(n).eq.0) then
+ ip1 = int((x-x0)*dxinv)
+ ip0 = ip1 - 1
+ ip2 = ip1 + 1
+
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=1-xx1
+
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+
+ a0=-0.5*xx1*xx2
+ a1=1.-xx1**2
+ a2=0.5*xx0*xx1
+
+ SC(ip0,jj,kk) = SC(ip0,jj,kk) + g1*a0
+ SC(ip1,jj,kk) = SC(ip1,jj,kk) + g1*a1
+ SC(ip2,jj,kk) = SC(ip2,jj,kk) + g1*a2
+
+ else
+
+ ip1 = int((x-x0)*dxinv)
+ ip0 = ip1 - 1
+ ip2 = ip1 + 1
+ ip3 =ip1+2
+
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=1-xx1
+ xx3=2.-xx1
+
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ ip3=mod(ip3+nx2,nx2) +1
+
+ a0 = .5*((2.-xx0)**2)*(1.-xx0)
+ a1 = 1.-2.5*xx1*xx1 + 1.5*xx1*xx1*xx1
+ a2 = 1.-2.5*xx2*xx2 + 1.5*xx2*xx2*xx2
+ a3 = .5*((2.-xx3)**2)*(1.-xx3)
+
+
+ SC(ip0,jj,kk) = SC(ip0,jj,kk) + g1*a0
+ SC(ip1,jj,kk) = SC(ip1,jj,kk) + g1*a1
+ SC(ip2,jj,kk) = SC(ip2,jj,kk) + g1*a2
+ SC(ip3,jj,kk) = SC(ip3,jj,kk) + g1*a3
+ end if
+ end do
+
+!!! go to 222 !! -> ca veut dire quoi goto ?
+!!! do k=3,nx2-2
+!!! do j=1,nx2
+!!! do i=1,nx2
+!!! if (SC(i,j,k).gt.1.2) then
+!!! print*,'BOUM', i,j,k,SC(i,j,k)
+!!! goto 222
+!!! end if
+!!! end do
+!!! end do
+!!! end do
+
+end subroutine remeshx
+
+subroutine tag_particles(npt,npart_aux,ntag)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: npt, npart_aux, ntag
+
+ integer :: ntype(npg),ncfl(npg),npart_bl(npg),i_nbl(npg)
+ real(WP) :: amin_lambda(npg), cflm
+
+ real(WP) :: x0, dx, dx_bl, dx_bl_inv
+ integer :: nblock, m, nbl, i, j, ini, jj, jc
+
+ x0=xmin
+
+ cflm = cfl
+
+
+ m = 2 !np_bl+1
+ nblock=nxsc/(m)
+ dx = L/nxsc
+ dx_bl=float(m)*dx
+ dx_bl_inv=1./dx_bl
+
+ do nbl=1,nblock
+ amin_lambda(nbl)=111.
+ npart_bl(nbl)=0
+ i_nbl(nbl)=0
+ enddo
+
+ do i=1,npt
+ nbl=1+int((xp0(i)-x0+0.00001)*dx_bl_inv)
+ amin_lambda(nbl)=amin1(amin_lambda(nbl),vx(i)*cflm)
+ npart_bl(nbl)=npart_bl(nbl)+1
+ i_nbl(nbl)=i
+ enddo
+
+ do nbl=1,nblock-1
+ if (i_nbl(nbl).ne.0) then
+ ini=i_nbl(nbl)
+ if ((ini.lt.npt).and.(xp0(ini+1).lt.xp0(ini)+1.5*dx)) then
+ amin_lambda(nbl)=amin1(amin_lambda(nbl),vx(ini+1)*cflm)
+ end if
+ end if
+ end do
+
+ nbl=nblock
+ if (i_nbl(nbl).ne.0) then
+ if ((xp0(npt).ge.xmax-1.5*dx).and.(xp0(1).lt.xmin+0.5*dx)) then
+ amin_lambda(nbl)=amin1(amin_lambda(nbl),vx(1)*cflm)
+ end if
+ end if
+
+ do nbl=1,nblock
+ if (amin_lambda(nbl).lt.nint(amin_lambda(nbl))) then
+ ntype(nbl)=1
+ ncfl(nbl)=nint(amin_lambda(nbl))
+ else
+ ntype(nbl)=0
+ ncfl(nbl)=int(amin_lambda(nbl))
+ if (amin_lambda(nbl).lt.0) ncfl(nbl)=int(amin_lambda(nbl))-1
+ endif
+ end do
+
+ do i=1,npt
+ nbl=1+int((xp0(i)-x0+0.00001)*dx_bl_inv)
+ itype(i)=ntype(nbl)
+ icfl(i)=ncfl(nbl)
+ end do
+
+ ntag=0
+ npart_aux=0
+ j=2
+ jc=0
+ do i=2,npt-1
+ j=j+jc
+ if (j.lt.npt) then
+ jj=j+1
+ if ((icfl(j).ne.icfl(jj)).and.(itype(j).ne.itype(jj)).and.(xp0(jj).le.xp0(j)+1.5*dx).and.1.eq.2) then !!! TO AVOID TO TAG PARTICLES
+ ntag=ntag+1
+ itag(ntag)=j
+ ntag=ntag+1
+ itag(ntag)=j+1
+ jc=2
+ else
+ npart_aux=npart_aux+1
+ xp_aux(npart_aux)=xp0(j)
+ up_aux(npart_aux)=up(j)
+ vx_aux(npart_aux)=vx(j)
+ itype_aux(npart_aux)=itype(j)
+ jc=1
+ end if
+ end if
+ end do
+
+ if (npt.ge.1) then
+
+ if ((icfl(1).ne.icfl(npt)).and.(itype(1).ne.itype(npt)) &
+ & .and.(xp0(npt).ge.xp0(1)+(float(nxsc)-1.5)*dx) &
+ & .and.(itag(ntag).ne.npt).and.1.eq.2) then !!! TO AVOID TO TAG PARTICLES
+ ntag=ntag+1
+ itag(ntag)=npt
+ ntag=ntag+1
+ itag(ntag)=1
+ else
+ npart_aux=npart_aux+1
+ xp_aux(npart_aux)=xp0(1)
+ up_aux(npart_aux)=up(1)
+ vx_aux(npart_aux)=vx(1)
+ itype_aux(npart_aux)=itype(1)
+ if (npt.gt.1) then
+ npart_aux=npart_aux+1
+ xp_aux(npart_aux)=xp0(npt)
+ up_aux(npart_aux)=up(npt)
+ vx_aux(npart_aux)=vx(npt)
+ itype_aux(npart_aux)=itype(npt)
+ end if
+ end if
+ end if
+
+end subroutine tag_particles
+
+subroutine remeshx_tag(ntag,jj,kk)
+
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: ntag, jj, kk
+ integer :: i, nx2, n, ip0, ip1, ip2, k, j, ini, jp1, ip3, ip4
+ real(WP) :: dxinv, dx2, x0, g1, x, xx0, xx1, xx2, a0, a1, a2, b1, b2, y, u1, u2, yy0, yy1, yy2
+
+ nx2 = nxsc
+ dx2 = L/nx2
+
+ dxinv=1./dx2
+ x0=xmin
+
+ do n=1,ntag-1,2
+ i=itag(n)
+ ini=itag(n+1)
+ x=xp(i)
+ y=xp(ini)
+ u1=up(i)
+ u2=up(ini)
+ if (itype(i).eq.0) then
+ ip1 = int((x-x0)*dxinv)
+ jp1 = nint((y-x0)*dxinv)
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ yy1 = (y - float(jp1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=1-xx1
+ yy0=yy1+1
+ yy2=1-yy1
+ ip0=ip1-1
+ ip2=ip1+1
+ ip3=ip1+2
+ ip4=ip1+3
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ ip3=mod(ip3+nx2,nx2) +1
+ ip4=mod(ip4+nx2,nx2) +1
+ a0=-xx1*xx2/2.
+ a1=xx0*xx2
+ b1=yy0*yy2
+ b2=yy0*yy1/2.
+ SC(ip0,jj,kk)=SC(ip0,jj,kk)+a0*u1
+ SC(ip1,jj,kk)=SC(ip1,jj,kk)+a1*u1+(1.+yy1-b1-b2)*u2
+ SC(ip2,jj,kk)=SC(ip2,jj,kk)+xx1*u1-yy1*u2
+ SC(ip3,jj,kk)=SC(ip3,jj,kk)+(1.-a0-a1-xx1)*u1+b1*u2
+ SC(ip4,jj,kk)=SC(ip4,jj,kk)+b2*u2
+ else
+ ip1 = nint((x-x0)*dxinv)
+ jp1 = int((y-x0)*dxinv)
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ yy1 = (y - float(jp1)*dx2-x0)*dxinv
+ xx2=1-xx1
+ yy0=yy1+1
+ ip0=ip1-1
+ ip2=ip1+1
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ a0=-0.5*xx1*xx2
+ b2=0.5*yy0*yy1
+ SC(ip0,jj,kk)=SC(ip0,jj,kk)+a0*u1
+ SC(ip1,jj,kk)=SC(ip1,jj,kk)+(1.-a0)*u1+(1.-b2)*u2
+ SC(ip2,jj,kk)=SC(ip2,jj,kk)+b2*u2
+ endif
+ enddo
+
+end subroutine remeshx_tag
diff --git a/LEGI/AutresRemaillage/advec_O2SM_y.f90 b/LEGI/AutresRemaillage/advec_O2SM_y.f90
new file mode 100644
index 000000000..78e74f93d
--- /dev/null
+++ b/LEGI/AutresRemaillage/advec_O2SM_y.f90
@@ -0,0 +1,336 @@
+subroutine y_advec
+
+ use x_advec_m
+
+ implicit none
+
+ real(WP) :: dx, yy, zz
+ integer :: i,j,k,np,n, ir, kr, ini, np_aux, ntag
+
+
+ ntag = 0
+ ntag_total = 0
+ np_aux = 0
+
+ dx = L/nxsc
+ cfl = dt/dx
+
+ do k=1,nxsc
+ zz=xmin+float(k-1)*dx
+ do i=1,nxsc
+ np=0
+ yy=xmin+float(i-1)*dx
+ do j=1,nxsc
+ if (abs(SC(i,j,k)).gt.circlim) then
+ np=np+1
+ up(np)=SC(i,j,k)
+ xp(np)=xmin+float(j-1)*dx
+ endif
+ enddo
+ if (np.ne.0) then
+ call velox_y(np,i,k)
+ do n=1,np
+ xp0(n)=xp(n)
+ xp(n)=xp(n)+0.5*dt*vx(n)
+ if (xp(n).gt.xmax) xp(n)=xp(n)-xmax+xmin
+ if (xp(n).lt.xmin) xp(n)=xp(n)+xmax-xmin
+ enddo
+ call velox_y(np,i,k)
+
+!! debut modif
+
+ call tag_particles(np,np_aux,ntag)
+ ntag_total = ntag_total + ntag
+ ! test
+ !itype(1:np_aux) = 1
+
+! np integer: deja defini (input) - local
+! np_aux integer: deja defini (output) - local
+! ntag integer: deja defini (output) - local
+! np_bl=1 --> supprimer de l'appel
+! icfl tablo integer: deja defini - common --> supprimer de l'appel
+! itype tablo integer: deja defini - common --> supprimer de l'appel
+! itype_aux tablo integer: deja defini - common --> supprimer de l'appel
+! itag tablo integer: deja defini - common --> supprimer de l'appel
+! xp_aux tablo real integer: deja defini - common --> supprimer de l'appel
+! up_aux tablo real integer: deja defini - common --> supprimer de l'appel
+! vx_aux tablo real integer: deja defini - common --> supprimer de l'appel
+
+! ---- > trouver les dimension max des tablo
+! ---- > trouver si les tablo doivent etre declarer en local ou en common
+
+ if (ntag.ne.0) then
+ do n=1,ntag
+ ini=itag(n)
+ xp(ini)=xp0(ini)+dt*vx(ini)
+ if (xp(ini).gt.xmax) xp(ini)=xp(ini)-xmax+xmin
+ if (xp(ini).lt.xmin) xp(ini)=xp(ini)+xmax-xmin !! TO CHECK WITH GH
+ end do
+ end if
+
+ do n=1,np_aux
+ xp_aux(n)=xp_aux(n)+dt*vx_aux(n)
+ if (xp_aux(n).gt.xmax) xp_aux(n)=xp_aux(n)-xmax+xmin
+ if (xp_aux(n).lt.xmin) xp_aux(n)=xp_aux(n)+xmax-xmin
+ end do
+ ir=i
+ kr=k
+ call remeshy(np_aux,ir,kr)
+ if (ntag.ne.0) call remeshy_tag(ntag,ir,kr)
+
+!! fin modif
+
+ endif
+ enddo
+ enddo
+
+ print*,'ntag y_advec :', ntag_total
+
+
+end subroutine y_advec
+
+
+
+subroutine velox_y(np,j,k)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: np, j, k, i, nx3, ny3, nz3, jp1, jp2, kp1, kp2, ip1, ip2
+ real(WP) :: dx3, dy3, dz3, dxinv, yy, zz, x0, y0, z0, dx2
+ real(WP) :: yy1, yy2, zz1, zz2, b1, b2, c1, c2, a1, a2
+ real(WP) :: x, xx1, xx2
+
+ do i=1,np
+ vx(i)=0.
+ end do
+
+ nx3=nxsc
+ ny3=nx3
+ nz3=nx3
+ dx3=xmax/float(nx3)
+ dy3=dx3
+ dz3=dx3
+
+ dxinv=1./dx3
+
+ x0=xmin
+ y0=ymin
+ z0=zmin
+
+ dx2 = L/nxsc
+
+ yy=+float(j-1)*dx2
+ zz=+float(k-1)*dx2
+
+ jp1 = int((yy)/dx3)
+ jp2 = jp1 + 1
+ kp1 = int((zz)/dx3)
+ kp2 = kp1 + 1
+ yy1 = (yy-float(jp1)*dx3)/dx3
+ yy2 = 1.0 - yy1
+ zz1 = (zz-float(kp1)*dx3)/dx3
+ zz2 = 1.0 - zz1
+ b1=yy2
+ b2=yy1
+ c1=zz2
+ c2=zz1
+
+ jp1=mod(jp1+nx3,nx3) +1
+ jp2=mod(jp2+nx3,nx3) +1
+ kp1=mod(kp1+nx3,nx3) +1
+ kp2=mod(kp2+nx3,nx3) +1
+
+ do i = 1,np
+
+ x = xp(i)
+
+ ip1 = int((x-x0)*dxinv)
+ ip2 = ip1 + 1
+
+ xx1 = (x - float(ip1)*dx3-x0)*dxinv
+ xx2=1-xx1
+
+ ip1=mod(ip1+nx3,nx3) +1
+ ip2=mod(ip2+nx3,nx3) +1
+
+ a1 = xx2
+ a2 = xx1
+
+ vx(i)= vx(i) + Vsc(jp1,ip1,kp1)*a1*b1*c1
+ vx(i)= vx(i) + Vsc(jp1,ip2,kp1)*a2*b1*c1
+ vx(i)= vx(i) + Vsc(jp2,ip1,kp1)*a1*b2*c1
+ vx(i)= vx(i) + Vsc(jp2,ip2,kp1)*a2*b2*c1
+ vx(i)= vx(i) + Vsc(jp1,ip1,kp2)*a1*b1*c2
+ vx(i)= vx(i) + Vsc(jp1,ip2,kp2)*a2*b1*c2
+ vx(i)= vx(i) + Vsc(jp2,ip1,kp2)*a1*b2*c2
+ vx(i)= vx(i) + Vsc(jp2,ip2,kp2)*a2*b2*c2
+
+ end do
+
+end subroutine velox_y
+
+subroutine remeshy(np1,jj,kk)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: np1, jj, kk
+ integer :: i, nx2, n, ip0, ip1, ip2, k, j, ip3
+ real(WP) :: dxinv, dx2, x0, g1, x, xx0, xx1, xx2, a0, a1, a2, a3, xx3
+
+ nx2 = nxsc
+ dx2 = L/nx2
+
+ dxinv=1./dx2
+
+ do i=1,nx2
+ SC(jj,i,kk)=0.
+ enddo
+
+ x0=xmin
+
+ do n = 1,np1
+ g1 = up_aux(n)
+ x = xp_aux(n)
+
+ if (itype_aux(n).eq.0) then
+ ip1 = int((x-x0)*dxinv)
+ ip0 = ip1 - 1
+ ip2 = ip1 + 1
+
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=1-xx1
+
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+
+ a0=-0.5*xx1*xx2
+ a1=1.-xx1**2
+ a2=0.5*xx0*xx1
+
+ SC(jj,ip0,kk) = SC(jj,ip0,kk) + g1*a0
+ SC(jj,ip1,kk) = SC(jj,ip1,kk) + g1*a1
+ SC(jj,ip2,kk) = SC(jj,ip2,kk) + g1*a2
+
+ else
+
+ ip1 = int((x-x0)*dxinv)
+ ip0 = ip1 - 1
+ ip2 = ip1 + 1
+ ip3 =ip1+2
+
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=1-xx1
+ xx3=2.-xx1
+
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ ip3=mod(ip3+nx2,nx2) +1
+
+ a0 = .5*((2.-xx0)**2)*(1.-xx0)
+ a1 = 1.-2.5*xx1*xx1 + 1.5*xx1*xx1*xx1
+ a2 = 1.-2.5*xx2*xx2 + 1.5*xx2*xx2*xx2
+ a3 = .5*((2.-xx3)**2)*(1.-xx3)
+
+
+ SC(jj,ip0,kk) = SC(jj,ip0,kk) + g1*a0
+ SC(jj,ip1,kk) = SC(jj,ip1,kk) + g1*a1
+ SC(jj,ip2,kk) = SC(jj,ip2,kk) + g1*a2
+ SC(jj,ip3,kk) = SC(jj,ip3,kk) + g1*a3
+
+ end if
+ end do
+
+!!! go to 222 !! -> ca veut dire quoi goto ?
+!!! do k=3,nx2-2
+!!! do j=1,nx2
+!!! do i=1,nx2
+!!! if (SC(i,j,k).gt.1.2) then
+!!! print*,'BOUM', i,j,k,SC(i,j,k)
+!!! goto 222
+!!! end if
+!!! end do
+!!! end do
+!!! end do
+
+end subroutine remeshy
+
+subroutine remeshy_tag(ntag,jj,kk)
+
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: jj, kk, ntag
+ integer :: i, nx2, n, ip0, ip1, ip2, k, j, ini, jp1, ip3, ip4
+ real(WP) :: dxinv, dx2, x0, g1, x, xx0, xx1, xx2, a0, a1, a2, u1, u2, y, yy1, yy0, yy2, b1, b2
+
+ nx2 = nxsc
+ dx2 = L/nx2
+
+ dxinv=1./dx2
+ x0=xmin
+
+ do n=1,ntag,2
+ i=itag(n)
+ ini=itag(n+1)
+ x=xp(i)
+ y=xp(ini)
+ u1=up(i)
+ u2=up(ini)
+ if (itype(i).eq.0) then
+ ip1 = int((x-x0)*dxinv)
+ jp1 = nint((y-x0)*dxinv)
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ yy1 = (y - float(jp1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=1-xx1
+ yy0=yy1+1
+ yy2=1-yy1
+ ip0=ip1-1
+ ip2=ip1+1
+ ip3=ip1+2
+ ip4=ip1+3
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ ip3=mod(ip3+nx2,nx2) +1
+ ip4=mod(ip4+nx2,nx2) +1
+ a0=-xx1*xx2/2.
+ a1=xx0*xx2
+ b1=yy0*yy2
+ b2=yy0*yy1/2.
+ SC(jj,ip0,kk)=SC(jj,ip0,kk)+a0*u1
+ SC(jj,ip1,kk)=SC(jj,ip1,kk)+a1*u1+(1.+yy1-b1-b2)*u2
+ SC(jj,ip2,kk)=SC(jj,ip2,kk)+xx1*u1-yy1*u2
+ SC(jj,ip3,kk)=SC(jj,ip3,kk)+(1.-a0-a1-xx1)*u1+b1*u2
+ SC(jj,ip4,kk)=SC(jj,ip4,kk)+b2*u2
+ else
+ ip1 = nint((x-x0)*dxinv)
+ jp1 = int((y-x0)*dxinv)
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ yy1 = (y - float(jp1)*dx2-x0)*dxinv
+ xx2=1-xx1
+ yy0=yy1+1
+ ip0=ip1-1
+ ip2=ip1+1
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ a0=-0.5*xx1*xx2
+ b2=0.5*yy0*yy1
+ SC(jj,ip0,kk)=SC(jj,ip0,kk)+a0*u1
+ SC(jj,ip1,kk)=SC(jj,ip1,kk)+(1.-a0)*u1+(1.-b2)*u2
+ SC(jj,ip2,kk)=SC(jj,ip2,kk)+b2*u2
+ endif
+ enddo
+
+end subroutine remeshy_tag
diff --git a/LEGI/AutresRemaillage/advec_O2SM_z.f90 b/LEGI/AutresRemaillage/advec_O2SM_z.f90
new file mode 100644
index 000000000..e59525a59
--- /dev/null
+++ b/LEGI/AutresRemaillage/advec_O2SM_z.f90
@@ -0,0 +1,335 @@
+subroutine z_advec
+
+ use x_advec_m
+
+ implicit none
+
+ real(WP) :: dx, yy, zz
+ integer :: i,j,k,np,n, ir, jr, ini, ntag, np_aux
+
+
+ dx = L/nxsc
+ cfl = dt/dx
+
+ ntag = 0
+ ntag_total = 0
+ np_aux = 0
+
+ do j=1,nxsc
+ zz=xmin+float(k-1)*dx
+ do i=1,nxsc
+ np=0
+ yy=xmin+float(i-1)*dx
+ do k=1,nxsc
+ if (abs(SC(i,j,k)).gt.circlim) then
+ np=np+1
+ up(np)=SC(i,j,k)
+ xp(np)=xmin+float(k-1)*dx
+ endif
+ enddo
+ if (np.ne.0) then
+ call velox_z(np,i,j)
+ do n=1,np
+ xp0(n)=xp(n)
+ xp(n)=xp(n)+0.5*dt*vx(n)
+ if (xp(n).gt.xmax) xp(n)=xp(n)-xmax+xmin
+ if (xp(n).lt.xmin) xp(n)=xp(n)+xmax-xmin
+ enddo
+ call velox_z(np,i,j)
+
+!! debut modif
+
+ call tag_particles(np,np_aux,ntag)
+ ntag_total = ntag_total + ntag
+ ! test
+ !itype(1:np_aux) = 1
+
+! np integer: deja defini (input) - local
+! np_aux integer: deja defini (output) - local
+! ntag integer: deja defini (output) - local
+! np_bl=1 --> supprimer de l'appel
+! icfl tablo integer: deja defini - common --> supprimer de l'appel
+! itype tablo integer: deja defini - common --> supprimer de l'appel
+! itype_aux tablo integer: deja defini - common --> supprimer de l'appel
+! itag tablo integer: deja defini - common --> supprimer de l'appel
+! xp_aux tablo real integer: deja defini - common --> supprimer de l'appel
+! up_aux tablo real integer: deja defini - common --> supprimer de l'appel
+! vx_aux tablo real integer: deja defini - common --> supprimer de l'appel
+
+! ---- > trouver les dimension max des tablo
+! ---- > trouver si les tablo doivent etre declarer en local ou en common
+
+ if (ntag.ne.0) then
+ do n=1,ntag
+ ini=itag(n)
+ xp(ini)=xp0(ini)+dt*vx(ini)
+ if (xp(ini).gt.xmax) xp(ini)=xp(ini)-xmax+xmin
+ if (xp(ini).lt.xmin) xp(ini)=xp(ini)+xmax-xmin !! TO CHECK WITH GH
+ end do
+ end if
+
+ do n=1,np_aux
+ xp_aux(n)=xp_aux(n)+dt*vx_aux(n)
+ if (xp_aux(n).gt.xmax) xp_aux(n)=xp_aux(n)-xmax+xmin
+ if (xp_aux(n).lt.xmin) xp_aux(n)=xp_aux(n)+xmax-xmin
+ end do
+ ir=i
+ jr=j
+ call remeshz(np_aux,ir,jr)
+ if (ntag.ne.0) call remeshz_tag(ntag,ir,jr)
+
+!! fin modif
+
+ endif
+ enddo
+ enddo
+
+ print*,'ntag z_advec :', ntag_total
+
+end subroutine z_advec
+
+
+
+subroutine velox_z(np,j,k)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: np, j, k, i, nx3, ny3, nz3, jp1, jp2, kp1, kp2, ip1, ip2
+ real(WP) :: dx3, dy3, dz3, dxinv, yy, zz, x0, y0, z0, dx2
+ real(WP) :: yy1, yy2, zz1, zz2, b1, b2, c1, c2, a1, a2
+ real(WP) :: x, xx1, xx2
+
+ do i=1,np
+ vx(i)=0.
+ end do
+
+ nx3=nxsc
+ ny3=nx3
+ nz3=nx3
+ dx3=xmax/float(nx3)
+ dy3=dx3
+ dz3=dx3
+
+ dxinv=1./dx3
+
+ x0=xmin
+ y0=ymin
+ z0=zmin
+
+ dx2 = L/nxsc
+
+ yy=+float(j-1)*dx2
+ zz=+float(k-1)*dx2
+
+ jp1 = int((yy)/dx3)
+ jp2 = jp1 + 1
+ kp1 = int((zz)/dx3)
+ kp2 = kp1 + 1
+ yy1 = (yy-float(jp1)*dx3)/dx3
+ yy2 = 1.0 - yy1
+ zz1 = (zz-float(kp1)*dx3)/dx3
+ zz2 = 1.0 - zz1
+ b1=yy2
+ b2=yy1
+ c1=zz2
+ c2=zz1
+
+ jp1=mod(jp1+nx3,nx3) +1
+ jp2=mod(jp2+nx3,nx3) +1
+ kp1=mod(kp1+nx3,nx3) +1
+ kp2=mod(kp2+nx3,nx3) +1
+
+ do i = 1,np
+
+ x = xp(i)
+
+ ip1 = int((x-x0)*dxinv)
+ ip2 = ip1 + 1
+
+ xx1 = (x - float(ip1)*dx3-x0)*dxinv
+ xx2=1-xx1
+
+ ip1=mod(ip1+nx3,nx3) +1
+ ip2=mod(ip2+nx3,nx3) +1
+
+ a1 = xx2
+ a2 = xx1
+
+ vx(i)= vx(i) + Wsc(jp1,kp1,ip1)*a1*b1*c1
+ vx(i)= vx(i) + Wsc(jp1,kp1,ip2)*a2*b1*c1
+ vx(i)= vx(i) + Wsc(jp2,kp1,ip1)*a1*b2*c1
+ vx(i)= vx(i) + Wsc(jp2,kp1,ip2)*a2*b2*c1
+ vx(i)= vx(i) + Wsc(jp1,kp2,ip1)*a1*b1*c2
+ vx(i)= vx(i) + Wsc(jp1,kp2,ip2)*a2*b1*c2
+ vx(i)= vx(i) + Wsc(jp2,kp2,ip1)*a1*b2*c2
+ vx(i)= vx(i) + Wsc(jp2,kp2,ip2)*a2*b2*c2
+
+ end do
+
+end subroutine velox_z
+
+subroutine remeshz(np1,jj,kk)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: np1, jj, kk
+ integer :: i, nx2, n, ip0, ip1, ip2, k, j, ip3
+ real(WP) :: dxinv, dx2, x0, g1, x, xx0, xx1, xx2, a0, a1, a2, a3, xx3
+
+ nx2 = nxsc
+ dx2 = L/nx2
+
+ dxinv=1./dx2
+
+ do i=1,nx2
+ SC(jj,kk,i)=0.
+ enddo
+
+ x0=xmin
+
+ do n = 1,np1
+ g1 = up_aux(n)
+ x = xp_aux(n)
+
+ if (itype_aux(n).eq.0) then
+ ip1 = int((x-x0)*dxinv)
+ ip0 = ip1 - 1
+ ip2 = ip1 + 1
+
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=1-xx1
+
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+
+ a0=-0.5*xx1*xx2
+ a1=1.-xx1**2
+ a2=0.5*xx0*xx1
+
+ SC(jj,kk,ip0) = SC(jj,kk,ip0) + g1*a0
+ SC(jj,kk,ip1) = SC(jj,kk,ip1) + g1*a1
+ SC(jj,kk,ip2) = SC(jj,kk,ip2) + g1*a2
+
+ else
+
+ ip1 = int((x-x0)*dxinv)
+ ip0 = ip1 - 1
+ ip2 = ip1 + 1
+ ip3 =ip1+2
+
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=1-xx1
+ xx3=2.-xx1
+
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ ip3=mod(ip3+nx2,nx2) +1
+
+ a0 = .5*((2.-xx0)**2)*(1.-xx0)
+ a1 = 1.-2.5*xx1*xx1 + 1.5*xx1*xx1*xx1
+ a2 = 1.-2.5*xx2*xx2 + 1.5*xx2*xx2*xx2
+ a3 = .5*((2.-xx3)**2)*(1.-xx3)
+
+
+ SC(jj,kk,ip0) = SC(jj,kk,ip0) + g1*a0
+ SC(jj,kk,ip1) = SC(jj,kk,ip1) + g1*a1
+ SC(jj,kk,ip2) = SC(jj,kk,ip2) + g1*a2
+ SC(jj,kk,ip3) = SC(jj,kk,ip3) + g1*a3
+
+ end if
+ end do
+
+!!! go to 222 !! -> ca veut dire quoi goto ?
+!!! do k=3,nx2-2
+!!! do j=1,nx2
+!!! do i=1,nx2
+!!! if (SC(i,j,k).gt.1.2) then
+!!! print*,'BOUM', i,j,k,SC(i,j,k)
+!!! goto 222
+!!! end if
+!!! end do
+!!! end do
+!!! end do
+
+end subroutine remeshz
+
+subroutine remeshz_tag(ntag,kk,jj)
+
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: ntag, jj, kk
+ integer :: i, nx2, n, ip0, ip1, ip2, k, j, ini, jp1, ip3, ip4
+ real(WP) :: dxinv, dx2, x0, g1, x, xx0, xx1, xx2, a0, a1, a2, b1, b2, y, u1, u2, yy0, yy1, yy2
+
+ nx2 = nxsc
+ dx2 = L/nx2
+
+ dxinv=1./dx2
+ x0=xmin
+
+ do n=1,ntag,2
+ i=itag(n)
+ ini=itag(n+1)
+ x=xp(i)
+ y=xp(ini)
+ u1=up(i)
+ u2=up(ini)
+ if (itype(i).eq.0) then
+ ip1 = int((x-x0)*dxinv)
+ jp1 = nint((y-x0)*dxinv)
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ yy1 = (y - float(jp1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=1-xx1
+ yy0=yy1+1
+ yy2=1-yy1
+ ip0=ip1-1
+ ip2=ip1+1
+ ip3=ip1+2
+ ip4=ip1+3
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ ip3=mod(ip3+nx2,nx2) +1
+ ip4=mod(ip4+nx2,nx2) +1
+ a0=-xx1*xx2/2.
+ a1=xx0*xx2
+ b1=yy0*yy2
+ b2=yy0*yy1/2.
+ SC(kk,jj,ip0)=SC(kk,jj,ip0)+a0*u1
+ SC(kk,jj,ip1)=SC(kk,jj,ip1)+a1*u1+(1.+yy1-b1-b2)*u2
+ SC(kk,jj,ip2)=SC(kk,jj,ip2)+xx1*u1-yy1*u2
+ SC(kk,jj,ip3)=SC(kk,jj,ip3)+(1.-a0-a1-xx1)*u1+b1*u2
+ SC(kk,jj,ip4)=SC(kk,jj,ip4)+b2*u2
+ else
+ ip1 = nint((x-x0)*dxinv)
+ jp1 = int((y-x0)*dxinv)
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ yy1 = (y - float(jp1)*dx2-x0)*dxinv
+ xx2=1-xx1
+ yy0=yy1+1
+ ip0=ip1-1
+ ip2=ip1+1
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ a0=-0.5*xx1*xx2
+ b2=0.5*yy0*yy1
+ SC(kk,jj,ip0)=SC(kk,jj,ip0)+a0*u1
+ SC(kk,jj,ip1)=SC(kk,jj,ip1)+(1.-a0)*u1+(1.-b2)*u2
+ SC(kk,jj,ip2)=SC(kk,jj,ip2)+b2*u2
+ endif
+ enddo
+
+end subroutine remeshz_tag
diff --git a/LEGI/AutresRemaillage/advec_O4_x.f90 b/LEGI/AutresRemaillage/advec_O4_x.f90
new file mode 100644
index 000000000..f5cd0c6db
--- /dev/null
+++ b/LEGI/AutresRemaillage/advec_O4_x.f90
@@ -0,0 +1,541 @@
+module x_advec_m
+
+ use solver
+ use data
+
+ implicit none
+
+ integer :: npart, ntag_total, npg
+ real(WP) :: xmin, ymin, zmin, xmax, ymax, zmax, cfl
+ integer, dimension(:), allocatable :: itype
+ integer, dimension(:), allocatable :: itype_aux
+ integer, dimension(:), allocatable :: itag
+ integer, dimension(:), allocatable :: icfl
+ real(WP), dimension(:), allocatable :: up
+ real(WP), dimension(:), allocatable :: up_aux
+ real(WP), dimension(:), allocatable :: xp
+ real(WP), dimension(:), allocatable :: xp0
+ real(WP), dimension(:), allocatable :: xp_aux
+ real(WP), dimension(:), allocatable :: vx
+ real(WP), dimension(:), allocatable :: vx_aux
+
+ real(WP) :: circlim
+
+end module x_advec_m
+
+
+
+subroutine x_advec_init
+
+ use x_advec_m
+
+ implicit none
+
+ npg = nxsc !100000
+
+ allocate(up(npg))
+ allocate(up_aux(npg))
+ allocate(xp(npg)) !!-> code original npg = 256 mais ca marche pas
+ allocate(xp0(npg))
+ allocate(xp_aux(npg))
+ allocate(vx(npg))
+ allocate(vx_aux(npg))
+ allocate(itype(npg))
+ allocate(itype_aux(npg))
+ allocate(itag(npg))
+ allocate(icfl(npg))
+
+ up = 0.0
+ up_aux = 0.0
+ xp = 0.0
+ xp0 = 0.0
+ xp_aux = 0.0
+ vx = 0.0
+ vx_aux = 0.0
+ cfl = 0.0
+
+ itype = 0
+ itype_aux = 0
+ itag = 0
+ icfl = 0
+
+ npart = 0
+
+ xmin = 0.
+ ymin = 0.
+ zmin = 0.
+ xmax = L
+ ymax = L
+ zmax = L
+
+ circlim = -1 !10**(-5)
+
+end subroutine x_advec_init
+
+
+
+subroutine x_advec
+
+ use x_advec_m
+
+ implicit none
+
+ real(WP) :: dx, yy, zz
+ integer :: i,j,k,np,n, jr, kr, np_aux, ini, ntag
+
+ dx = L/nxsc
+ cfl = dt/dx
+
+ npart = 0
+ ntag_total = 0
+ ntag = 0
+
+ do k=1,nxsc
+ zz=xmin+float(k-1)*dx
+ do j=1,nxsc
+ np=0
+ yy=xmin+float(j-1)*dx
+ do i=1,nxsc
+ if (abs(SC(i,j,k)).gt.circlim) then
+ np=np+1
+ up(np)=SC(i,j,k)
+ xp(np)=xmin+float(i-1)*dx
+ end if
+ end do
+ if (np.ne.0) then
+ call velox_x(np,j,k)
+ do n=1,np
+ xp0(n)=xp(n)
+ xp(n)=xp(n)+0.5*dt*vx(n)
+ if (xp(n).gt.xmax) xp(n)=xp(n)-xmax+xmin
+ if (xp(n).lt.xmin) xp(n)=xp(n)+xmax-xmin
+ end do
+ call velox_x(np,j,k)
+
+!! debut modif
+
+
+ call tag_particles(np,np_aux,ntag)
+ ntag_total = ntag_total + ntag
+ ! test
+ !itype(1:np_aux) = 1
+
+! np integer: deja defini (input) - local
+! np_aux integer: deja defini (output) - local
+! ntag integer: deja defini (output) - local
+! np_bl=1 --> supprimer de l'appel
+! icfl tablo integer: deja defini - common --> supprimer de l'appel
+! itype tablo integer: deja defini - common --> supprimer de l'appel
+! itype_aux tablo integer: deja defini - common --> supprimer de l'appel
+! itag tablo integer: deja defini - common --> supprimer de l'appel
+! xp_aux tablo real integer: deja defini - common --> supprimer de l'appel
+! up_aux tablo real integer: deja defini - common --> supprimer de l'appel
+! vx_aux tablo real integer: deja defini - common --> supprimer de l'appel
+
+! ---- > trouver les dimension max des tablo
+! ---- > trouver si les tablo doivent etre declarer en local ou en common
+
+ if (ntag.ne.0) then
+ do n=1,ntag
+ ini=itag(n)
+ xp(ini)=xp0(ini)+dt*vx(ini)
+ if (xp(ini).gt.xmax) xp(ini)=xp(ini)-xmax+xmin
+ if (xp(ini).lt.xmin) xp(ini)=xp(ini)+xmax-xmin !! TO CHECK WITH GH
+ end do
+ end if
+
+ do n=1,np_aux
+ xp_aux(n)=xp_aux(n)+dt*vx_aux(n)
+ if (xp_aux(n).gt.xmax) xp_aux(n)=xp_aux(n)-xmax+xmin
+ if (xp_aux(n).lt.xmin) xp_aux(n)=xp_aux(n)+xmax-xmin
+ end do
+ jr=j
+ kr=k
+ call remeshx(np_aux,jr,kr)
+ if (ntag.ne.0) call remeshx_tag(ntag,jr,kr)
+
+!! fin modif
+
+ npart=npart+np
+ end if
+ end do
+ end do
+
+ print*, 'NPART, NTAG ', npart,ntag_total
+ print*,'ntag x_advec :', ntag_total
+
+
+end subroutine x_advec
+
+
+
+subroutine velox_x(np,j,k)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: np, j, k, i, nx3, ny3, nz3, jp1, jp2, kp1, kp2, ip1, ip2
+ real(WP) :: dx3, dy3, dz3, dxinv, yy, zz, x0, y0, z0, dx2
+ real(WP) :: yy1, yy2, zz1, zz2, b1, b2, c1, c2, a1, a2
+ real(WP) :: x, xx1, xx2
+
+ do i=1,np
+ vx(i)=0.
+ end do
+
+ nx3=nxsc
+ ny3=nx3
+ nz3=nx3
+ dx3=xmax/float(nx3)
+ dy3=dx3
+ dz3=dx3
+
+ dxinv=1./dx3
+
+ x0=xmin
+ y0=ymin
+ z0=zmin
+
+ dx2 = L/nxsc
+
+ yy=+float(j-1)*dx2
+ zz=+float(k-1)*dx2
+
+ jp1 = int((yy)/dx3)
+ jp2 = jp1 + 1
+ kp1 = int((zz)/dx3)
+ kp2 = kp1 + 1
+ yy1 = (yy-float(jp1)*dx3)/dx3
+ yy2 = 1.0 - yy1
+ zz1 = (zz-float(kp1)*dx3)/dx3
+ zz2 = 1.0 - zz1
+ b1=yy2
+ b2=yy1
+ c1=zz2
+ c2=zz1
+
+ jp1=mod(jp1+nx3,nx3) +1
+ jp2=mod(jp2+nx3,nx3) +1
+ kp1=mod(kp1+nx3,nx3) +1
+ kp2=mod(kp2+nx3,nx3) +1
+
+ do i = 1,np
+
+ x = xp(i)
+
+ ip1 = int((x-x0)*dxinv)
+ ip2 = ip1 + 1
+
+ xx1 = (x - float(ip1)*dx3-x0)*dxinv
+ xx2=1-xx1
+
+ ip1=mod(ip1+nx3,nx3) +1
+ ip2=mod(ip2+nx3,nx3) +1
+
+ a1 = xx2
+ a2 = xx1
+
+ vx(i)= vx(i) + Usc(ip1,jp1,kp1)*a1*b1*c1
+ vx(i)= vx(i) + Usc(ip2,jp1,kp1)*a2*b1*c1
+ vx(i)= vx(i) + Usc(ip1,jp2,kp1)*a1*b2*c1
+ vx(i)= vx(i) + Usc(ip2,jp2,kp1)*a2*b2*c1
+ vx(i)= vx(i) + Usc(ip1,jp1,kp2)*a1*b1*c2
+ vx(i)= vx(i) + Usc(ip2,jp1,kp2)*a2*b1*c2
+ vx(i)= vx(i) + Usc(ip1,jp2,kp2)*a1*b2*c2
+ vx(i)= vx(i) + Usc(ip2,jp2,kp2)*a2*b2*c2
+
+ end do
+
+end subroutine velox_x
+
+subroutine remeshx(np1,jj,kk)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: np1, jj, kk
+ integer :: i, nx2, n, ip0, ip1, ip2, k, j, ip3, ip4
+ real(WP) :: dxinv, dx2, x0, g1, x, xx0, xx1, xx2, a0, a1, a2, xx3, xx4, a3, a4
+
+ nx2 = nxsc
+ dx2 = L/nx2
+
+ dxinv=1./dx2
+
+ do i=1,nx2
+ SC(i,jj,kk)=0.
+ enddo
+
+ x0=xmin
+
+ do n = 1,np1
+ g1 = up_aux(n)
+ x = xp_aux(n)
+
+ if (itype_aux(n).eq.0) then
+ ip1 = int((x-x0)*dxinv)
+ ip0 = ip1 - 1
+ ip2 = ip1 + 1
+
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=1-xx1
+
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+
+ a0=-0.5*xx1*xx2
+ a1=1.-xx1**2
+ a2=0.5*xx0*xx1
+
+ SC(ip0,jj,kk) = SC(ip0,jj,kk) + g1*a0
+ SC(ip1,jj,kk) = SC(ip1,jj,kk) + g1*a1
+ SC(ip2,jj,kk) = SC(ip2,jj,kk) + g1*a2
+
+ else
+
+ ip1 = nint((x-x0)*dxinv)
+ ip0 = ip1 - 1
+ ip2 = ip1 + 1
+ ip3 =ip1-2
+ ip4=ip1+2
+
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=xx1-1
+ xx3=xx1+2
+ xx4=xx1-2
+
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ ip3=mod(ip3+nx2,nx2) +1
+ ip4=mod(ip4+nx2,nx2) +1
+
+ a0=-xx1*xx2*xx3*xx4/6.
+ a1=xx0*xx2*xx3*xx4/4.
+ a2=-xx0*xx1*xx3*xx4/6.
+ a3=xx0*xx1*xx2*xx4/24.
+ a4=xx0*xx1*xx2*xx3/24.
+
+ SC(ip0,jj,kk) = SC(ip0,jj,kk) + g1*a0
+ SC(ip1,jj,kk) = SC(ip1,jj,kk) + g1*a1
+ SC(ip2,jj,kk) = SC(ip2,jj,kk) + g1*a2
+ SC(ip3,jj,kk) = SC(ip3,jj,kk) + g1*a3
+ SC(ip4,jj,kk) = SC(ip4,jj,kk) + g1*a4
+
+ end if
+ end do
+
+!!! go to 222 !! -> ca veut dire quoi goto ?
+!!! do k=3,nx2-2
+!!! do j=1,nx2
+!!! do i=1,nx2
+!!! if (SC(i,j,k).gt.1.2) then
+!!! print*,'BOUM', i,j,k,SC(i,j,k)
+!!! goto 222
+!!! end if
+!!! end do
+!!! end do
+!!! end do
+
+end subroutine remeshx
+
+subroutine tag_particles(npt,npart_aux,ntag)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: npt, npart_aux, ntag
+
+ integer :: ntype(npg),ncfl(npg),npart_bl(npg),i_nbl(npg)
+ real(WP) :: amin_lambda(npg), cflm
+
+ real(WP) :: x0, dx, dx_bl, dx_bl_inv
+ integer :: nblock, m, nbl, i, j, ini, jj, jc
+
+ x0=xmin
+
+ cflm = cfl
+
+
+ m = 2 !np_bl+1
+ nblock=nxsc/(m)
+ dx = L/nxsc
+ dx_bl=float(m)*dx
+ dx_bl_inv=1./dx_bl
+
+ do nbl=1,nblock
+ amin_lambda(nbl)=111.
+ npart_bl(nbl)=0
+ i_nbl(nbl)=0
+ enddo
+
+ do i=1,npt
+ nbl=1+int((xp0(i)-x0+0.00001)*dx_bl_inv)
+ amin_lambda(nbl)=amin1(amin_lambda(nbl),vx(i)*cflm)
+ npart_bl(nbl)=npart_bl(nbl)+1
+ i_nbl(nbl)=i
+ enddo
+
+ do nbl=1,nblock-1
+ if (i_nbl(nbl).ne.0) then
+ ini=i_nbl(nbl)
+ if ((ini.lt.npt).and.(xp0(ini+1).lt.xp0(ini)+1.5*dx)) then
+ amin_lambda(nbl)=amin1(amin_lambda(nbl),vx(ini+1)*cflm)
+ end if
+ end if
+ end do
+
+ nbl=nblock
+ if (i_nbl(nbl).ne.0) then
+ if ((xp0(npt).ge.xmax-1.5*dx).and.(xp0(1).lt.xmin+0.5*dx)) then
+ amin_lambda(nbl)=amin1(amin_lambda(nbl),vx(1)*cflm)
+ end if
+ end if
+
+ do nbl=1,nblock
+ if (amin_lambda(nbl).lt.nint(amin_lambda(nbl))) then
+ ntype(nbl)=1
+ ncfl(nbl)=nint(amin_lambda(nbl))
+ else
+ ntype(nbl)=0
+ ncfl(nbl)=int(amin_lambda(nbl))
+ if (amin_lambda(nbl).lt.0) ncfl(nbl)=int(amin_lambda(nbl))-1
+ endif
+ end do
+
+ do i=1,npt
+ nbl=1+int((xp0(i)-x0+0.00001)*dx_bl_inv)
+ itype(i)=ntype(nbl)
+ icfl(i)=ncfl(nbl)
+ end do
+
+ ntag=0
+ npart_aux=0
+ j=2
+ jc=0
+ do i=2,npt-1
+ j=j+jc
+ if (j.lt.npt) then
+ jj=j+1
+ if ((icfl(j).ne.icfl(jj)).and.(itype(j).ne.itype(jj)).and.(xp0(jj).le.xp0(j)+1.5*dx)) then
+ ntag=ntag+1
+ itag(ntag)=j
+ ntag=ntag+1
+ itag(ntag)=j+1
+ jc=2
+ else
+ npart_aux=npart_aux+1
+ xp_aux(npart_aux)=xp0(j)
+ up_aux(npart_aux)=up(j)
+ vx_aux(npart_aux)=vx(j)
+ itype_aux(npart_aux)=itype(j)
+ jc=1
+ end if
+ end if
+ end do
+
+ if (npt.ge.1) then
+
+ if ((icfl(1).ne.icfl(npt)).and.(itype(1).ne.itype(npt)) &
+ & .and.(xp0(npt).ge.xp0(1)+(float(nxsc)-1.5)*dx) &
+ & .and.(itag(ntag).ne.npt)) then
+ ntag=ntag+1
+ itag(ntag)=npt
+ ntag=ntag+1
+ itag(ntag)=1
+ else
+ npart_aux=npart_aux+1
+ xp_aux(npart_aux)=xp0(1)
+ up_aux(npart_aux)=up(1)
+ vx_aux(npart_aux)=vx(1)
+ itype_aux(npart_aux)=itype(1)
+ if (npt.gt.1) then
+ npart_aux=npart_aux+1
+ xp_aux(npart_aux)=xp0(npt)
+ up_aux(npart_aux)=up(npt)
+ vx_aux(npart_aux)=vx(npt)
+ itype_aux(npart_aux)=itype(npt)
+ end if
+ end if
+ end if
+
+end subroutine tag_particles
+
+subroutine remeshx_tag(ntag,jj,kk)
+
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: ntag, jj, kk
+ integer :: i, nx2, n, ip0, ip1, ip2, k, j, ini, jp1, ip3, ip4
+ real(WP) :: dxinv, dx2, x0, g1, x, xx0, xx1, xx2, a0, a1, a2, b1, b2, y, u1, u2, yy0, yy1, yy2
+
+ nx2 = nxsc
+ dx2 = L/nx2
+
+ dxinv=1./dx2
+ x0=xmin
+
+ do n=1,ntag-1,2
+ i=itag(n)
+ ini=itag(n+1)
+ x=xp(i)
+ y=xp(ini)
+ u1=up(i)
+ u2=up(ini)
+ if (itype(i).eq.0) then
+ ip1 = int((x-x0)*dxinv)
+ jp1 = nint((y-x0)*dxinv)
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ yy1 = (y - float(jp1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=1-xx1
+ yy0=yy1+1
+ yy2=1-yy1
+ ip0=ip1-1
+ ip2=ip1+1
+ ip3=ip1+2
+ ip4=ip1+3
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ ip3=mod(ip3+nx2,nx2) +1
+ ip4=mod(ip4+nx2,nx2) +1
+ a0=-xx1*xx2/2.
+ a1=xx0*xx2
+ b1=yy0*yy2
+ b2=yy0*yy1/2.
+ SC(ip0,jj,kk)=SC(ip0,jj,kk)+a0*u1
+ SC(ip1,jj,kk)=SC(ip1,jj,kk)+a1*u1+(1.+yy1-b1-b2)*u2
+ SC(ip2,jj,kk)=SC(ip2,jj,kk)+xx1*u1-yy1*u2
+ SC(ip3,jj,kk)=SC(ip3,jj,kk)+(1.-a0-a1-xx1)*u1+b1*u2
+ SC(ip4,jj,kk)=SC(ip4,jj,kk)+b2*u2
+ else
+ ip1 = nint((x-x0)*dxinv)
+ jp1 = int((y-x0)*dxinv)
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ yy1 = (y - float(jp1)*dx2-x0)*dxinv
+ xx2=1-xx1
+ yy0=yy1+1
+ ip0=ip1-1
+ ip2=ip1+1
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ a0=-0.5*xx1*xx2
+ b2=0.5*yy0*yy1
+ SC(ip0,jj,kk)=SC(ip0,jj,kk)+a0*u1
+ SC(ip1,jj,kk)=SC(ip1,jj,kk)+(1.-a0)*u1+(1.-b2)*u2
+ SC(ip2,jj,kk)=SC(ip2,jj,kk)+b2*u2
+ endif
+ enddo
+
+end subroutine remeshx_tag
diff --git a/LEGI/AutresRemaillage/advec_O4_y.f90 b/LEGI/AutresRemaillage/advec_O4_y.f90
new file mode 100644
index 000000000..cfc81c2d0
--- /dev/null
+++ b/LEGI/AutresRemaillage/advec_O4_y.f90
@@ -0,0 +1,340 @@
+subroutine y_advec
+
+ use x_advec_m
+
+ implicit none
+
+ real(WP) :: dx, yy, zz
+ integer :: i,j,k,np,n, ir, kr, ini, np_aux, ntag
+
+
+ ntag = 0
+ ntag_total = 0
+ np_aux = 0
+
+ dx = L/nxsc
+ cfl = dt/dx
+
+ do k=1,nxsc
+ zz=xmin+float(k-1)*dx
+ do i=1,nxsc
+ np=0
+ yy=xmin+float(i-1)*dx
+ do j=1,nxsc
+ if (abs(SC(i,j,k)).gt.circlim) then
+ np=np+1
+ up(np)=SC(i,j,k)
+ xp(np)=xmin+float(j-1)*dx
+ endif
+ enddo
+ if (np.ne.0) then
+ call velox_y(np,i,k)
+ do n=1,np
+ xp0(n)=xp(n)
+ xp(n)=xp(n)+0.5*dt*vx(n)
+ if (xp(n).gt.xmax) xp(n)=xp(n)-xmax+xmin
+ if (xp(n).lt.xmin) xp(n)=xp(n)+xmax-xmin
+ enddo
+ call velox_y(np,i,k)
+
+!! debut modif
+
+ call tag_particles(np,np_aux,ntag)
+ ntag_total = ntag_total + ntag
+ ! test
+ !itype(1:np_aux) = 1
+
+! np integer: deja defini (input) - local
+! np_aux integer: deja defini (output) - local
+! ntag integer: deja defini (output) - local
+! np_bl=1 --> supprimer de l'appel
+! icfl tablo integer: deja defini - common --> supprimer de l'appel
+! itype tablo integer: deja defini - common --> supprimer de l'appel
+! itype_aux tablo integer: deja defini - common --> supprimer de l'appel
+! itag tablo integer: deja defini - common --> supprimer de l'appel
+! xp_aux tablo real integer: deja defini - common --> supprimer de l'appel
+! up_aux tablo real integer: deja defini - common --> supprimer de l'appel
+! vx_aux tablo real integer: deja defini - common --> supprimer de l'appel
+
+! ---- > trouver les dimension max des tablo
+! ---- > trouver si les tablo doivent etre declarer en local ou en common
+
+ if (ntag.ne.0) then
+ do n=1,ntag
+ ini=itag(n)
+ xp(ini)=xp0(ini)+dt*vx(ini)
+ if (xp(ini).gt.xmax) xp(ini)=xp(ini)-xmax+xmin
+ if (xp(ini).lt.xmin) xp(ini)=xp(ini)+xmax-xmin !! TO CHECK WITH GH
+ end do
+ end if
+
+ do n=1,np_aux
+ xp_aux(n)=xp_aux(n)+dt*vx_aux(n)
+ if (xp_aux(n).gt.xmax) xp_aux(n)=xp_aux(n)-xmax+xmin
+ if (xp_aux(n).lt.xmin) xp_aux(n)=xp_aux(n)+xmax-xmin
+ end do
+ ir=i
+ kr=k
+ call remeshy(np_aux,ir,kr)
+ if (ntag.ne.0) call remeshy_tag(ntag,ir,kr)
+
+!! fin modif
+
+ endif
+ enddo
+ enddo
+
+ print*,'ntag y_advec :', ntag_total
+
+
+end subroutine y_advec
+
+
+
+subroutine velox_y(np,j,k)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: np, j, k, i, nx3, ny3, nz3, jp1, jp2, kp1, kp2, ip1, ip2
+ real(WP) :: dx3, dy3, dz3, dxinv, yy, zz, x0, y0, z0, dx2
+ real(WP) :: yy1, yy2, zz1, zz2, b1, b2, c1, c2, a1, a2
+ real(WP) :: x, xx1, xx2
+
+ do i=1,np
+ vx(i)=0.
+ end do
+
+ nx3=nxsc
+ ny3=nx3
+ nz3=nx3
+ dx3=xmax/float(nx3)
+ dy3=dx3
+ dz3=dx3
+
+ dxinv=1./dx3
+
+ x0=xmin
+ y0=ymin
+ z0=zmin
+
+ dx2 = L/nxsc
+
+ yy=+float(j-1)*dx2
+ zz=+float(k-1)*dx2
+
+ jp1 = int((yy)/dx3)
+ jp2 = jp1 + 1
+ kp1 = int((zz)/dx3)
+ kp2 = kp1 + 1
+ yy1 = (yy-float(jp1)*dx3)/dx3
+ yy2 = 1.0 - yy1
+ zz1 = (zz-float(kp1)*dx3)/dx3
+ zz2 = 1.0 - zz1
+ b1=yy2
+ b2=yy1
+ c1=zz2
+ c2=zz1
+
+ jp1=mod(jp1+nx3,nx3) +1
+ jp2=mod(jp2+nx3,nx3) +1
+ kp1=mod(kp1+nx3,nx3) +1
+ kp2=mod(kp2+nx3,nx3) +1
+
+ do i = 1,np
+
+ x = xp(i)
+
+ ip1 = int((x-x0)*dxinv)
+ ip2 = ip1 + 1
+
+ xx1 = (x - float(ip1)*dx3-x0)*dxinv
+ xx2=1-xx1
+
+ ip1=mod(ip1+nx3,nx3) +1
+ ip2=mod(ip2+nx3,nx3) +1
+
+ a1 = xx2
+ a2 = xx1
+
+ vx(i)= vx(i) + Vsc(jp1,ip1,kp1)*a1*b1*c1
+ vx(i)= vx(i) + Vsc(jp1,ip2,kp1)*a2*b1*c1
+ vx(i)= vx(i) + Vsc(jp2,ip1,kp1)*a1*b2*c1
+ vx(i)= vx(i) + Vsc(jp2,ip2,kp1)*a2*b2*c1
+ vx(i)= vx(i) + Vsc(jp1,ip1,kp2)*a1*b1*c2
+ vx(i)= vx(i) + Vsc(jp1,ip2,kp2)*a2*b1*c2
+ vx(i)= vx(i) + Vsc(jp2,ip1,kp2)*a1*b2*c2
+ vx(i)= vx(i) + Vsc(jp2,ip2,kp2)*a2*b2*c2
+
+ end do
+
+end subroutine velox_y
+
+subroutine remeshy(np1,jj,kk)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: np1, jj, kk
+ integer :: i, nx2, n, ip0, ip1, ip2, k, j, ip3, ip4
+ real(WP) :: dxinv, dx2, x0, g1, x, xx0, xx1, xx2, a0, a1, a2, xx3, xx4, a3, a4
+
+ nx2 = nxsc
+ dx2 = L/nx2
+
+ dxinv=1./dx2
+
+ do i=1,nx2
+ SC(jj,i,kk)=0.
+ enddo
+
+ x0=xmin
+
+ do n = 1,np1
+ g1 = up_aux(n)
+ x = xp_aux(n)
+
+ if (itype_aux(n).eq.0) then
+ ip1 = int((x-x0)*dxinv)
+ ip0 = ip1 - 1
+ ip2 = ip1 + 1
+
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=1-xx1
+
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+
+ a0=-0.5*xx1*xx2
+ a1=1.-xx1**2
+ a2=0.5*xx0*xx1
+
+ SC(jj,ip0,kk) = SC(jj,ip0,kk) + g1*a0
+ SC(jj,ip1,kk) = SC(jj,ip1,kk) + g1*a1
+ SC(jj,ip2,kk) = SC(jj,ip2,kk) + g1*a2
+
+ else
+
+ ip1 = nint((x-x0)*dxinv)
+ ip0 = ip1 - 1
+ ip2 = ip1 + 1
+ ip3 =ip1-2
+ ip4=ip1+2
+
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=xx1-1
+ xx3=xx1+2
+ xx4=xx1-2
+
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ ip3=mod(ip3+nx2,nx2) +1
+ ip4=mod(ip4+nx2,nx2) +1
+
+ a0=-xx1*xx2*xx3*xx4/6.
+ a1=xx0*xx2*xx3*xx4/4.
+ a2=-xx0*xx1*xx3*xx4/6.
+ a3=xx0*xx1*xx2*xx4/24.
+ a4=xx0*xx1*xx2*xx3/24.
+
+ SC(jj,ip0,kk) = SC(jj,ip0,kk) + g1*a0
+ SC(jj,ip1,kk) = SC(jj,ip1,kk) + g1*a1
+ SC(jj,ip2,kk) = SC(jj,ip2,kk) + g1*a2
+ SC(jj,ip3,kk) = SC(jj,ip3,kk) + g1*a3
+ SC(jj,ip4,kk) = SC(jj,ip4,kk) + g1*a4
+
+ end if
+ end do
+
+!!! go to 222 !! -> ca veut dire quoi goto ?
+!!! do k=3,nx2-2
+!!! do j=1,nx2
+!!! do i=1,nx2
+!!! if (SC(i,j,k).gt.1.2) then
+!!! print*,'BOUM', i,j,k,SC(i,j,k)
+!!! goto 222
+!!! end if
+!!! end do
+!!! end do
+!!! end do
+
+end subroutine remeshy
+
+subroutine remeshy_tag(ntag,jj,kk)
+
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: jj, kk, ntag
+ integer :: i, nx2, n, ip0, ip1, ip2, k, j, ini, jp1, ip3, ip4
+ real(WP) :: dxinv, dx2, x0, g1, x, xx0, xx1, xx2, a0, a1, a2, u1, u2, y, yy1, yy0, yy2, b1, b2
+
+ nx2 = nxsc
+ dx2 = L/nx2
+
+ dxinv=1./dx2
+ x0=xmin
+
+ do n=1,ntag,2
+ i=itag(n)
+ ini=itag(n+1)
+ x=xp(i)
+ y=xp(ini)
+ u1=up(i)
+ u2=up(ini)
+ if (itype(i).eq.0) then
+ ip1 = int((x-x0)*dxinv)
+ jp1 = nint((y-x0)*dxinv)
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ yy1 = (y - float(jp1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=1-xx1
+ yy0=yy1+1
+ yy2=1-yy1
+ ip0=ip1-1
+ ip2=ip1+1
+ ip3=ip1+2
+ ip4=ip1+3
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ ip3=mod(ip3+nx2,nx2) +1
+ ip4=mod(ip4+nx2,nx2) +1
+ a0=-xx1*xx2/2.
+ a1=xx0*xx2
+ b1=yy0*yy2
+ b2=yy0*yy1/2.
+ SC(jj,ip0,kk)=SC(jj,ip0,kk)+a0*u1
+ SC(jj,ip1,kk)=SC(jj,ip1,kk)+a1*u1+(1.+yy1-b1-b2)*u2
+ SC(jj,ip2,kk)=SC(jj,ip2,kk)+xx1*u1-yy1*u2
+ SC(jj,ip3,kk)=SC(jj,ip3,kk)+(1.-a0-a1-xx1)*u1+b1*u2
+ SC(jj,ip4,kk)=SC(jj,ip4,kk)+b2*u2
+ else
+ ip1 = nint((x-x0)*dxinv)
+ jp1 = int((y-x0)*dxinv)
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ yy1 = (y - float(jp1)*dx2-x0)*dxinv
+ xx2=1-xx1
+ yy0=yy1+1
+ ip0=ip1-1
+ ip2=ip1+1
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ a0=-0.5*xx1*xx2
+ b2=0.5*yy0*yy1
+ SC(jj,ip0,kk)=SC(jj,ip0,kk)+a0*u1
+ SC(jj,ip1,kk)=SC(jj,ip1,kk)+(1.-a0)*u1+(1.-b2)*u2
+ SC(jj,ip2,kk)=SC(jj,ip2,kk)+b2*u2
+ endif
+ enddo
+
+end subroutine remeshy_tag
diff --git a/LEGI/AutresRemaillage/advec_O4_z.f90 b/LEGI/AutresRemaillage/advec_O4_z.f90
new file mode 100644
index 000000000..504d655c0
--- /dev/null
+++ b/LEGI/AutresRemaillage/advec_O4_z.f90
@@ -0,0 +1,339 @@
+subroutine z_advec
+
+ use x_advec_m
+
+ implicit none
+
+ real(WP) :: dx, yy, zz
+ integer :: i,j,k,np,n, ir, jr, ini, ntag, np_aux
+
+
+ dx = L/nxsc
+ cfl = dt/dx
+
+ ntag = 0
+ ntag_total = 0
+ np_aux = 0
+
+ do j=1,nxsc
+ zz=xmin+float(k-1)*dx
+ do i=1,nxsc
+ np=0
+ yy=xmin+float(i-1)*dx
+ do k=1,nxsc
+ if (abs(SC(i,j,k)).gt.circlim) then
+ np=np+1
+ up(np)=SC(i,j,k)
+ xp(np)=xmin+float(k-1)*dx
+ endif
+ enddo
+ if (np.ne.0) then
+ call velox_z(np,i,j)
+ do n=1,np
+ xp0(n)=xp(n)
+ xp(n)=xp(n)+0.5*dt*vx(n)
+ if (xp(n).gt.xmax) xp(n)=xp(n)-xmax+xmin
+ if (xp(n).lt.xmin) xp(n)=xp(n)+xmax-xmin
+ enddo
+ call velox_z(np,i,j)
+
+!! debut modif
+
+ call tag_particles(np,np_aux,ntag)
+ ntag_total = ntag_total + ntag
+ ! test
+ !itype(1:np_aux) = 1
+
+! np integer: deja defini (input) - local
+! np_aux integer: deja defini (output) - local
+! ntag integer: deja defini (output) - local
+! np_bl=1 --> supprimer de l'appel
+! icfl tablo integer: deja defini - common --> supprimer de l'appel
+! itype tablo integer: deja defini - common --> supprimer de l'appel
+! itype_aux tablo integer: deja defini - common --> supprimer de l'appel
+! itag tablo integer: deja defini - common --> supprimer de l'appel
+! xp_aux tablo real integer: deja defini - common --> supprimer de l'appel
+! up_aux tablo real integer: deja defini - common --> supprimer de l'appel
+! vx_aux tablo real integer: deja defini - common --> supprimer de l'appel
+
+! ---- > trouver les dimension max des tablo
+! ---- > trouver si les tablo doivent etre declarer en local ou en common
+
+ if (ntag.ne.0) then
+ do n=1,ntag
+ ini=itag(n)
+ xp(ini)=xp0(ini)+dt*vx(ini)
+ if (xp(ini).gt.xmax) xp(ini)=xp(ini)-xmax+xmin
+ if (xp(ini).lt.xmin) xp(ini)=xp(ini)+xmax-xmin !! TO CHECK WITH GH
+ end do
+ end if
+
+ do n=1,np_aux
+ xp_aux(n)=xp_aux(n)+dt*vx_aux(n)
+ if (xp_aux(n).gt.xmax) xp_aux(n)=xp_aux(n)-xmax+xmin
+ if (xp_aux(n).lt.xmin) xp_aux(n)=xp_aux(n)+xmax-xmin
+ end do
+ ir=i
+ jr=j
+ call remeshz(np_aux,ir,jr)
+ if (ntag.ne.0) call remeshz_tag(ntag,ir,jr)
+
+!! fin modif
+
+ endif
+ enddo
+ enddo
+
+ print*,'ntag z_advec :', ntag_total
+
+end subroutine z_advec
+
+
+
+subroutine velox_z(np,j,k)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: np, j, k, i, nx3, ny3, nz3, jp1, jp2, kp1, kp2, ip1, ip2
+ real(WP) :: dx3, dy3, dz3, dxinv, yy, zz, x0, y0, z0, dx2
+ real(WP) :: yy1, yy2, zz1, zz2, b1, b2, c1, c2, a1, a2
+ real(WP) :: x, xx1, xx2
+
+ do i=1,np
+ vx(i)=0.
+ end do
+
+ nx3=nxsc
+ ny3=nx3
+ nz3=nx3
+ dx3=xmax/float(nx3)
+ dy3=dx3
+ dz3=dx3
+
+ dxinv=1./dx3
+
+ x0=xmin
+ y0=ymin
+ z0=zmin
+
+ dx2 = L/nxsc
+
+ yy=+float(j-1)*dx2
+ zz=+float(k-1)*dx2
+
+ jp1 = int((yy)/dx3)
+ jp2 = jp1 + 1
+ kp1 = int((zz)/dx3)
+ kp2 = kp1 + 1
+ yy1 = (yy-float(jp1)*dx3)/dx3
+ yy2 = 1.0 - yy1
+ zz1 = (zz-float(kp1)*dx3)/dx3
+ zz2 = 1.0 - zz1
+ b1=yy2
+ b2=yy1
+ c1=zz2
+ c2=zz1
+
+ jp1=mod(jp1+nx3,nx3) +1
+ jp2=mod(jp2+nx3,nx3) +1
+ kp1=mod(kp1+nx3,nx3) +1
+ kp2=mod(kp2+nx3,nx3) +1
+
+ do i = 1,np
+
+ x = xp(i)
+
+ ip1 = int((x-x0)*dxinv)
+ ip2 = ip1 + 1
+
+ xx1 = (x - float(ip1)*dx3-x0)*dxinv
+ xx2=1-xx1
+
+ ip1=mod(ip1+nx3,nx3) +1
+ ip2=mod(ip2+nx3,nx3) +1
+
+ a1 = xx2
+ a2 = xx1
+
+ vx(i)= vx(i) + Wsc(jp1,kp1,ip1)*a1*b1*c1
+ vx(i)= vx(i) + Wsc(jp1,kp1,ip2)*a2*b1*c1
+ vx(i)= vx(i) + Wsc(jp2,kp1,ip1)*a1*b2*c1
+ vx(i)= vx(i) + Wsc(jp2,kp1,ip2)*a2*b2*c1
+ vx(i)= vx(i) + Wsc(jp1,kp2,ip1)*a1*b1*c2
+ vx(i)= vx(i) + Wsc(jp1,kp2,ip2)*a2*b1*c2
+ vx(i)= vx(i) + Wsc(jp2,kp2,ip1)*a1*b2*c2
+ vx(i)= vx(i) + Wsc(jp2,kp2,ip2)*a2*b2*c2
+
+ end do
+
+end subroutine velox_z
+
+subroutine remeshz(np1,jj,kk)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: np1, jj, kk
+ integer :: i, nx2, n, ip0, ip1, ip2, k, j, ip3, ip4
+ real(WP) :: dxinv, dx2, x0, g1, x, xx0, xx1, xx2, a0, a1, a2, xx3, xx4, a3, a4
+
+ nx2 = nxsc
+ dx2 = L/nx2
+
+ dxinv=1./dx2
+
+ do i=1,nx2
+ SC(jj,kk,i)=0.
+ enddo
+
+ x0=xmin
+
+ do n = 1,np1
+ g1 = up_aux(n)
+ x = xp_aux(n)
+
+ if (itype_aux(n).eq.0) then
+ ip1 = int((x-x0)*dxinv)
+ ip0 = ip1 - 1
+ ip2 = ip1 + 1
+
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=1-xx1
+
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+
+ a0=-0.5*xx1*xx2
+ a1=1.-xx1**2
+ a2=0.5*xx0*xx1
+
+ SC(jj,kk,ip0) = SC(jj,kk,ip0) + g1*a0
+ SC(jj,kk,ip1) = SC(jj,kk,ip1) + g1*a1
+ SC(jj,kk,ip2) = SC(jj,kk,ip2) + g1*a2
+
+ else
+
+ ip1 = nint((x-x0)*dxinv)
+ ip0 = ip1 - 1
+ ip2 = ip1 + 1
+ ip3 =ip1-2
+ ip4=ip1+2
+
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=xx1-1
+ xx3=xx1+2
+ xx4=xx1-2
+
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ ip3=mod(ip3+nx2,nx2) +1
+ ip4=mod(ip4+nx2,nx2) +1
+
+ a0=-xx1*xx2*xx3*xx4/6.
+ a1=xx0*xx2*xx3*xx4/4.
+ a2=-xx0*xx1*xx3*xx4/6.
+ a3=xx0*xx1*xx2*xx4/24.
+ a4=xx0*xx1*xx2*xx3/24.
+
+ SC(jj,kk,ip0) = SC(jj,kk,ip0) + g1*a0
+ SC(jj,kk,ip1) = SC(jj,kk,ip1) + g1*a1
+ SC(jj,kk,ip2) = SC(jj,kk,ip2) + g1*a2
+ SC(jj,kk,ip3) = SC(jj,kk,ip3) + g1*a3
+ SC(jj,kk,ip4) = SC(jj,kk,ip4) + g1*a4
+
+ end if
+ end do
+
+!!! go to 222 !! -> ca veut dire quoi goto ?
+!!! do k=3,nx2-2
+!!! do j=1,nx2
+!!! do i=1,nx2
+!!! if (SC(i,j,k).gt.1.2) then
+!!! print*,'BOUM', i,j,k,SC(i,j,k)
+!!! goto 222
+!!! end if
+!!! end do
+!!! end do
+!!! end do
+
+end subroutine remeshz
+
+subroutine remeshz_tag(ntag,kk,jj)
+
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: ntag, jj, kk
+ integer :: i, nx2, n, ip0, ip1, ip2, k, j, ini, jp1, ip3, ip4
+ real(WP) :: dxinv, dx2, x0, g1, x, xx0, xx1, xx2, a0, a1, a2, b1, b2, y, u1, u2, yy0, yy1, yy2
+
+ nx2 = nxsc
+ dx2 = L/nx2
+
+ dxinv=1./dx2
+ x0=xmin
+
+ do n=1,ntag,2
+ i=itag(n)
+ ini=itag(n+1)
+ x=xp(i)
+ y=xp(ini)
+ u1=up(i)
+ u2=up(ini)
+ if (itype(i).eq.0) then
+ ip1 = int((x-x0)*dxinv)
+ jp1 = nint((y-x0)*dxinv)
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ yy1 = (y - float(jp1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=1-xx1
+ yy0=yy1+1
+ yy2=1-yy1
+ ip0=ip1-1
+ ip2=ip1+1
+ ip3=ip1+2
+ ip4=ip1+3
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ ip3=mod(ip3+nx2,nx2) +1
+ ip4=mod(ip4+nx2,nx2) +1
+ a0=-xx1*xx2/2.
+ a1=xx0*xx2
+ b1=yy0*yy2
+ b2=yy0*yy1/2.
+ SC(kk,jj,ip0)=SC(kk,jj,ip0)+a0*u1
+ SC(kk,jj,ip1)=SC(kk,jj,ip1)+a1*u1+(1.+yy1-b1-b2)*u2
+ SC(kk,jj,ip2)=SC(kk,jj,ip2)+xx1*u1-yy1*u2
+ SC(kk,jj,ip3)=SC(kk,jj,ip3)+(1.-a0-a1-xx1)*u1+b1*u2
+ SC(kk,jj,ip4)=SC(kk,jj,ip4)+b2*u2
+ else
+ ip1 = nint((x-x0)*dxinv)
+ jp1 = int((y-x0)*dxinv)
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ yy1 = (y - float(jp1)*dx2-x0)*dxinv
+ xx2=1-xx1
+ yy0=yy1+1
+ ip0=ip1-1
+ ip2=ip1+1
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ a0=-0.5*xx1*xx2
+ b2=0.5*yy0*yy1
+ SC(kk,jj,ip0)=SC(kk,jj,ip0)+a0*u1
+ SC(kk,jj,ip1)=SC(kk,jj,ip1)+(1.-a0)*u1+(1.-b2)*u2
+ SC(kk,jj,ip2)=SC(kk,jj,ip2)+b2*u2
+ endif
+ enddo
+
+end subroutine remeshz_tag
diff --git a/LEGI/AutresRemaillage/advec_O4v2_x.f90 b/LEGI/AutresRemaillage/advec_O4v2_x.f90
new file mode 100644
index 000000000..d6b472f7c
--- /dev/null
+++ b/LEGI/AutresRemaillage/advec_O4v2_x.f90
@@ -0,0 +1,640 @@
+module x_advec_m
+
+ use solver
+ use data
+
+ implicit none
+
+ integer :: npart, ntag_total, npg
+ real(WP) :: xmin, ymin, zmin, xmax, ymax, zmax, cfl
+ integer, dimension(:), allocatable :: itype
+ integer, dimension(:), allocatable :: itype_aux
+ integer, dimension(:), allocatable :: itag
+ integer, dimension(:), allocatable :: icfl
+ real(WP), dimension(:), allocatable :: up
+ real(WP), dimension(:), allocatable :: up_aux
+ real(WP), dimension(:), allocatable :: xp
+ real(WP), dimension(:), allocatable :: xp0
+ real(WP), dimension(:), allocatable :: xp_aux
+ real(WP), dimension(:), allocatable :: vx
+ real(WP), dimension(:), allocatable :: vx_aux
+
+ real(WP) :: circlim
+
+end module x_advec_m
+
+
+
+subroutine x_advec_init
+
+ use x_advec_m
+
+ implicit none
+
+ npg = nxsc !100000
+
+ allocate(up(npg))
+ allocate(up_aux(npg))
+ allocate(xp(npg)) !!-> code original npg = 256 mais ca marche pas
+ allocate(xp0(npg))
+ allocate(xp_aux(npg))
+ allocate(vx(npg))
+ allocate(vx_aux(npg))
+ allocate(itype(npg))
+ allocate(itype_aux(npg))
+ allocate(itag(npg))
+ allocate(icfl(npg))
+
+ up = 0.0
+ up_aux = 0.0
+ xp = 0.0
+ xp0 = 0.0
+ xp_aux = 0.0
+ vx = 0.0
+ vx_aux = 0.0
+ cfl = 0.0
+
+ itype = 0
+ itype_aux = 0
+ itag = 0
+ icfl = 0
+
+ npart = 0
+
+ xmin = 0.
+ ymin = 0.
+ zmin = 0.
+ xmax = L
+ ymax = L
+ zmax = L
+
+ circlim = -1 !10**(-5)
+
+end subroutine x_advec_init
+
+
+
+subroutine x_advec
+
+ use x_advec_m
+
+ implicit none
+
+ real(WP) :: dx, yy, zz,tm
+ integer :: i,j,k,np,n, jr, kr, np_aux, ini, ntag
+ integer :: jj, j1, j2, j3, j4, ic, m
+
+ dx = L/nxsc
+ cfl = dt/dx
+
+ npart = 0
+ ntag_total = 0
+ ntag = 0
+ tm = 0.
+ m = 4 !np_bl+1 (et np_bl=1)
+
+ do k=1,nxsc
+ zz=xmin+float(k-1)*dx
+ do j=1,nxsc
+ np=0
+ yy=xmin+float(j-1)*dx
+ jj=3
+ ic=0
+ do i=3,nxsc+2
+ jj=jj+ic
+ if (jj.le.nxsc+2) then
+ j1=mod(jj-1,nxsc)+1
+ j2=mod(jj,nxsc)+1
+ j3=mod(jj+1,nxsc)+1
+ j4=mod(jj+2,nxsc)+1
+ if ((mod(j1,m).eq.m-1).and. &
+ & (abs(SC(j1,j,k))+abs(SC(j2,j,k))+abs(SC(j3,j,k))+abs(SC(j4,j,k)).ge.circlim)) then !!! QUESTION : UG = SC ???
+ np=np+1
+ up(np)=SC(j1,j,k)
+ tm=tm+SC(j1,j,k)
+ xp(np)=xmin+float(j1-1)*dx
+ np=np+1
+ up(np)=SC(j2,j,k)
+ tm=tm+SC(j2,j,k)
+ xp(np)=xmin+float(j2-1)*dx
+ np=np+1
+ up(np)=SC(j3,j,k)
+ tm=tm+SC(j3,j,k)
+ xp(np)=xmin+float(j3-1)*dx
+ np=np+1
+ up(np)=SC(j4,j,k)
+ tm=tm+SC(j4,j,k)
+ xp(np)=xmin+float(j4-1)*dx
+ ic=4
+ elseif ((abs(SC(j1,j,k)).ge.circlim)) then
+ np=np+1
+ up(np)=SC(j1,j,k)
+ tm=tm+SC(j1,j,k)
+ xp(np)=xmin+float(j1-1)*dx
+ ic=1
+ else
+ ic=1
+ endif
+ endif
+ enddo
+! print*,i,j,k,np
+ if (np.ne.0) then
+ call velox_x(np,j,k)
+ do n=1,np
+ xp0(n)=xp(n)
+ xp(n)=xp(n)+0.5*dt*vx(n)
+ if (xp(n).gt.xmax) xp(n)=xp(n)-xmax+xmin
+ if (xp(n).lt.xmin) xp(n)=xp(n)+xmax-xmin
+ end do
+ call velox_x(np,j,k)
+
+!! debut modif
+
+
+ call tag_particles(np,np_aux,ntag)
+ ntag_total = ntag_total + ntag
+
+! print*,i,j,k,ntag_total,ntag
+
+ ! test
+ !itype(1:np_aux) = 1
+
+! np integer: deja defini (input) - local
+! np_aux integer: deja defini (output) - local
+! ntag integer: deja defini (output) - local
+! np_bl=1 --> supprimer de l'appel
+! icfl tablo integer: deja defini - common --> supprimer de l'appel
+! itype tablo integer: deja defini - common --> supprimer de l'appel
+! itype_aux tablo integer: deja defini - common --> supprimer de l'appel
+! itag tablo integer: deja defini - common --> supprimer de l'appel
+! xp_aux tablo real integer: deja defini - common --> supprimer de l'appel
+! up_aux tablo real integer: deja defini - common --> supprimer de l'appel
+! vx_aux tablo real integer: deja defini - common --> supprimer de l'appel
+
+! ---- > trouver les dimension max des tablo
+! ---- > trouver si les tablo doivent etre declarer en local ou en common
+
+ if (ntag.ne.0) then
+ do n=1,ntag
+ ini=itag(n)
+ xp(ini)=xp0(ini)+dt*vx(ini)
+ if (xp(ini).gt.xmax) xp(ini)=xp(ini)-xmax+xmin
+ if (xp(ini).lt.xmin) xp(ini)=xp(ini)+xmax-xmin !! TO CHECK WITH GH
+ end do
+ end if
+
+ do n=1,np_aux
+ xp_aux(n)=xp_aux(n)+dt*vx_aux(n)
+ if (xp_aux(n).gt.xmax) xp_aux(n)=xp_aux(n)-xmax+xmin
+ if (xp_aux(n).lt.xmin) xp_aux(n)=xp_aux(n)+xmax-xmin
+ end do
+ jr=j
+ kr=k
+ call remeshx(np_aux,jr,kr)
+ if (ntag.ne.0) call remeshx_tag(ntag,jr,kr)
+
+!! fin modif
+
+ npart=npart+np
+ end if
+ end do
+ end do
+
+ print*, 'NPART, NTAG ', npart,ntag_total
+ print*,'ntag x_advec :', ntag_total
+
+
+end subroutine x_advec
+
+
+
+subroutine velox_x(np,j,k)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: np, j, k, i, nx3, ny3, nz3, jp1, jp2, kp1, kp2, ip1, ip2
+ real(WP) :: dx3, dy3, dz3, dxinv, yy, zz, x0, y0, z0, dx2
+ real(WP) :: yy1, yy2, zz1, zz2, b1, b2, c1, c2, a1, a2
+ real(WP) :: x, xx1, xx2
+
+ do i=1,np
+ vx(i)=0.
+ end do
+
+ nx3=nxsc
+ ny3=nx3
+ nz3=nx3
+ dx3=xmax/float(nx3)
+ dy3=dx3
+ dz3=dx3
+
+ dxinv=1./dx3
+
+ x0=xmin
+ y0=ymin
+ z0=zmin
+
+ dx2 = L/nxsc
+
+ yy=+float(j-1)*dx2
+ zz=+float(k-1)*dx2
+
+ jp1 = int((yy)/dx3)
+ jp2 = jp1 + 1
+ kp1 = int((zz)/dx3)
+ kp2 = kp1 + 1
+ yy1 = (yy-float(jp1)*dx3)/dx3
+ yy2 = 1.0 - yy1
+ zz1 = (zz-float(kp1)*dx3)/dx3
+ zz2 = 1.0 - zz1
+ b1=yy2
+ b2=yy1
+ c1=zz2
+ c2=zz1
+
+ jp1=mod(jp1+nx3,nx3) +1
+ jp2=mod(jp2+nx3,nx3) +1
+ kp1=mod(kp1+nx3,nx3) +1
+ kp2=mod(kp2+nx3,nx3) +1
+
+ do i = 1,np
+
+ x = xp(i)
+
+ ip1 = int((x-x0)*dxinv)
+ ip2 = ip1 + 1
+
+ xx1 = (x - float(ip1)*dx3-x0)*dxinv
+ xx2=1-xx1
+
+ ip1=mod(ip1+nx3,nx3) +1
+ ip2=mod(ip2+nx3,nx3) +1
+
+ a1 = xx2
+ a2 = xx1
+
+ vx(i)= vx(i) + Usc(ip1,jp1,kp1)*a1*b1*c1
+ vx(i)= vx(i) + Usc(ip2,jp1,kp1)*a2*b1*c1
+ vx(i)= vx(i) + Usc(ip1,jp2,kp1)*a1*b2*c1
+ vx(i)= vx(i) + Usc(ip2,jp2,kp1)*a2*b2*c1
+ vx(i)= vx(i) + Usc(ip1,jp1,kp2)*a1*b1*c2
+ vx(i)= vx(i) + Usc(ip2,jp1,kp2)*a2*b1*c2
+ vx(i)= vx(i) + Usc(ip1,jp2,kp2)*a1*b2*c2
+ vx(i)= vx(i) + Usc(ip2,jp2,kp2)*a2*b2*c2
+
+ end do
+
+end subroutine velox_x
+
+subroutine remeshx(np1,jj,kk)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: np1, jj, kk
+ integer :: i, nx2, n, ip0, ip1, ip2, k, j, ip3, ip4
+ real(WP) :: dxinv, dx2, x0, g1, x, xx0, xx1, xx2, a0, a1, a2, xx3, xx4, a3, a4
+
+ nx2 = nxsc
+ dx2 = L/nx2
+
+ dxinv=1./dx2
+
+ do i=1,nx2
+ SC(i,jj,kk)=0.
+ enddo
+
+ x0=xmin
+
+ do n = 1,np1
+ g1 = up_aux(n)
+ x = xp_aux(n)
+
+ if (itype_aux(n).eq.0) then !!! QUESTION : TOUJOURS D'ACTUALIT2 ? N'APPARAIT PLUS DANS REMESHX DE SRC_COTTET
+ ip1 = int((x-x0)*dxinv)
+ else
+ ip1 = nint((x-x0)*dxinv)
+ end if
+ ip0 = ip1 - 1
+ ip2 = ip1 + 1
+ ip3 =ip1-2
+ ip4=ip1+2
+
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=xx1-1
+ xx3=xx1+2
+ xx4=xx1-2
+
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ ip3=mod(ip3+nx2,nx2) +1
+ ip4=mod(ip4+nx2,nx2) +1
+
+ a0=-xx1*xx2*xx3*xx4/6.
+ a1=xx0*xx2*xx3*xx4/4.
+ a2=-xx0*xx1*xx3*xx4/6.
+ a3=xx0*xx1*xx2*xx4/24.
+ a4=xx0*xx1*xx2*xx3/24.
+
+ SC(ip0,jj,kk) = SC(ip0,jj,kk) + g1*a0
+ SC(ip1,jj,kk) = SC(ip1,jj,kk) + g1*a1
+ SC(ip2,jj,kk) = SC(ip2,jj,kk) + g1*a2
+ SC(ip3,jj,kk) = SC(ip3,jj,kk) + g1*a3
+ SC(ip4,jj,kk) = SC(ip4,jj,kk) + g1*a4
+
+ end do
+
+!!! go to 222 !! -> ca veut dire quoi goto ?
+!!! do k=3,nx2-2
+!!! do j=1,nx2
+!!! do i=1,nx2
+!!! if (SC(i,j,k).gt.1.2) then
+!!! print*,'BOUM', i,j,k,SC(i,j,k)
+!!! goto 222
+!!! end if
+!!! end do
+!!! end do
+!!! end do
+
+end subroutine remeshx
+
+subroutine tag_particles(npt,npart_aux,ntag)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: npt, npart_aux, ntag
+
+ integer :: ntype(npg),ncfl(npg),npart_bl(npg),i_nbl(npg)
+ real(WP) :: amin_lambda(npg), cflm
+
+ real(WP) :: x0, dx, dx_bl, dx_bl_inv
+ integer :: nblock, m, nbl, i, j, ini, jj, jc, j1, jjj
+
+ x0=xmin
+
+ cflm = cfl
+
+ m = 4 !np_bl+1
+ nblock=nxsc/(m)
+ dx = L/nxsc
+ dx_bl=float(m)*dx
+ dx_bl_inv=1./dx_bl
+
+ do nbl=1,nblock
+ amin_lambda(nbl)=111.
+ npart_bl(nbl)=0
+ i_nbl(nbl)=0
+ enddo
+
+ do i=1,npt
+ nbl=1+int((xp0(i)-x0+0.0001)*dx_bl_inv)
+ amin_lambda(nbl)=amin1(amin_lambda(nbl),vx(i)*cflm)
+ npart_bl(nbl)=npart_bl(nbl)+1
+ i_nbl(nbl)=i
+ enddo
+
+ do nbl=1,nblock-1
+ if (i_nbl(nbl).ne.0) then
+ ini=i_nbl(nbl)
+ if ((ini.lt.npt).and.(xp0(ini+1).lt.xp0(ini)+1.5*dx)) then
+ amin_lambda(nbl)=amin1(amin_lambda(nbl),vx(ini+1)*cflm)
+ end if
+ end if
+ end do
+
+ nbl=nblock
+ if (i_nbl(nbl).ne.0) then
+ ini = i_nbl(nbl)
+ if ((xp0(ini).ge.xmax-1.5*dx).and.(xp0(ini+1).lt.xmin+0.5*dx)) then
+ amin_lambda(nbl)=amin1(amin_lambda(nbl),vx(1)*cflm)
+ end if
+ end if
+
+ do nbl=1,nblock
+ if (amin_lambda(nbl).lt.nint(amin_lambda(nbl))) then
+ ntype(nbl)=1
+ ncfl(nbl)=nint(amin_lambda(nbl))
+ else
+ ntype(nbl)=0
+ ncfl(nbl)=int(amin_lambda(nbl))
+ if (amin_lambda(nbl).lt.0) ncfl(nbl)=int(amin_lambda(nbl))-1
+ endif
+ end do
+
+ do i=1,npt
+ nbl=1+int((xp0(i)-x0+0.0001)*dx_bl_inv)
+ itype(i)=ntype(nbl)
+ icfl(i)=ncfl(nbl)
+ end do
+
+ ntag=0
+ npart_aux=0
+ j=1
+ jc=0
+ do i=1,npt !-1
+ j=j+jc
+ if (j.le.npt) then
+ j1=j
+ jj=j+1
+ jjj=j+2
+ if ((jjj.le.npt).and.(icfl(jj).ne.icfl(jjj)) &
+! & .and.1.eq.2 & !! TEST
+ & .and.(itype(jj).ne.itype(jjj)) &
+ & .and.(xp0(jjj).le.xmin+mod(xp0(jj)+1.5*dx-xmin,xmax-xmin))) then
+ ntag=ntag+1
+ itag(ntag)=j1
+ ntag=ntag+1
+ itag(ntag)=jj
+ ntag=ntag+1
+ itag(ntag)=jjj
+ ntag=ntag+1
+ itag(ntag)=j+3
+ jc=4
+! print*, j1,jj,jjj
+! print*,xp0(j1),xp0(jj),xp0(jjj)
+! print*,icfl(jj),icfl(jjj)
+! print*,itype(jj),itype(jjj)
+! print*,xp0(jjj),xmin+mod(xp0(jj)+1.5*dx-xmin,xmax-xmin)
+! stop
+ else
+ npart_aux=npart_aux+1
+ xp_aux(npart_aux)=xp(j1)
+ up_aux(npart_aux)=up(j1)
+ vx_aux(npart_aux)=vx(j1)
+ itype_aux(npart_aux)=itype(j1)
+ jc=1
+ endif
+ endif
+ enddo
+
+end subroutine tag_particles
+
+subroutine remeshx_tag(ntag,jj,kk)
+
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: ntag, jj, kk
+ integer :: i, nx2, n, ip1, ip2, k, j, ini, inni, innni,jp1, ip3, ip4
+ integer :: ip1p, ip2p, ip0, ip1m
+ integer :: ip3m, ip2m, ip5p, ip4p, ip3p
+ real(WP) :: dxinv, dx2, x0, g1, x, xx0, xx1m, xx1p, xx2p, xx2, a0, a1, a2, b1, b2, y, u1, u2, yy0, yy1, yy2
+ real(WP) :: x1p, x2p, u1p, u2p, x1m, u1m, u0
+ real(WP) :: coef1m, coef0, coef1p, coef2p
+
+ nx2 = nxsc
+ dx2 = L/nx2
+
+ dxinv=1./dx2
+
+ do n=1,ntag-1,4
+ i=itag(n)
+ ini=itag(n+1)
+ inni=itag(n+2)
+ innni=itag(n+3)
+ x1m=xp(i)
+ x0=xp(ini)
+ x1p=xp(inni)
+ x2p=xp(innni)
+ u1m=up(i)
+ u0=up(ini)
+ u1p=up(inni)
+ u2p=up(innni)
+ if (itype(i).eq.0) then
+ ip0 = int((x0-xmin)*dxinv)
+ ip1m = int((x1m-xmin)*dxinv)
+ ip1p = nint((x1p-xmin)*dxinv)
+ ip2p = nint((x2p-xmin)*dxinv)
+ else
+ ip0 = nint((x0-xmin)*dxinv)
+ ip1m = nint((x1m-xmin)*dxinv)
+ ip1p = int((x1p-xmin)*dxinv)
+ ip2p = int((x2p-xmin)*dxinv)
+ end if
+ xx0 = (x0 - float(ip0)*dx2-xmin)*dxinv
+ xx1m = (x1m - float(ip1m)*dx2-xmin)*dxinv
+ xx1p = (x1p - float(ip1p)*dx2-xmin)*dxinv
+ xx2p = (x2p - float(ip2p)*dx2-xmin)*dxinv
+ ip3m=ip0-3
+ ip2m=ip0-2
+ ip1m=ip0-1
+ ip5p=ip0+5
+ ip4p=ip0+4
+ ip3p=ip0+3
+ ip2p=ip0+2
+ ip1p=ip0+1
+ ip0=mod(ip0+nxsc,nxsc) +1
+ ip1p=mod(ip1p+nxsc,nxsc) +1
+ ip2p=mod(ip2p+nxsc,nxsc) +1
+ ip3p=mod(ip3p+nxsc,nxsc) +1
+ ip4p=mod(ip4p+nxsc,nxsc) +1
+ ip5p=mod(ip5p+nxsc,nxsc) +1
+ ip1m=mod(ip1m+nxsc,nxsc) +1
+ ip2m=mod(ip2m+nxsc,nxsc) +1
+ ip3m=mod(ip3m+nxsc,nxsc) +1
+
+ if (itype(i).eq.0) then
+
+ coef1m=(xx1m-2.)*(xx1m-1.)*xx1m*(xx1m+1)/24.
+ SC(ip3m,jj,kk)=SC(ip3m,jj,kk)+u1m*coef1m
+
+ coef1m=-(xx1m-2.)*(xx1m-1.)*xx1m*(xx1m+2)/6.
+ coef0=(xx0-2.)*(xx0-1.)*xx0*(xx0+1)/24.
+ SC(ip2m,jj,kk)=SC(ip2m,jj,kk)+u1m*coef1m+u0*coef0
+
+ coef1m=(xx1m-2.)*(xx1m-1.)*(xx1m+1.)*(xx1m+2)/4.
+ coef0=-(xx0-2.)*(xx0-1.)*xx0*(xx0+2)/6.
+ coef1p=(xx1p-1.)*(xx1p-0.)*(xx1p+1.)*(xx1p+2)/24.
+ SC(ip1m,jj,kk)=SC(ip1m,jj,kk)+u1m*coef1m+u0*coef0+u1p*coef1p
+
+ coef1m=-(xx1m-2.)*(xx1m-0.)*(xx1m+1.)*(xx1m+2)/6.
+ coef0=(xx0-2.)*(xx0-1.)*(xx0+1.)*(xx0+2)/4.
+ coef1p=-(xx1p-1.)*(xx1p-0.)*(xx1p+1.)*(xx1p+3)/6.
+ coef2p=(xx2p-1.)*(xx2p-0.)*(xx2p+1.)*(xx2p+2)/24.
+ SC(ip0,jj,kk)=SC(ip0,jj,kk)+u1m*coef1m+u0*coef0+u1p*coef1p+u2p*coef2p
+
+ coef1m=(xx1m-1.)*(xx1m-0.)*(xx1m+1.)*(xx1m+2)/24.
+ coef1m=coef1m-(xx1m-2.)*(xx1m-1.)*(xx1m+0.)*(xx1m+1)/24.
+ coef0=-(xx0-2.)*(xx0-0.)*(xx0+1.)*(xx0+2)/6.
+ coef0=coef0+(xx0-3.)*(xx0-1.)*(xx0+0.)*(xx0+1)/6.
+ coef0=coef0+(xx0-1.)*(xx0-0.)*(xx0+1.)*(xx0+2)/24.
+ coef0=coef0-(xx0-2.)*(xx0-1.)*(xx0+0.)*(xx0+1)/24.
+ coef1p=(xx1p-2.)*(xx1p-1.)*(xx1p+0.)*(xx1p+1)/24.
+ coef1p=coef1p-(xx1p-1.)*(xx1p-0.)*(xx1p+1.)*(xx1p+2)/24.
+ coef1p=coef1p-(xx1p-2.)*(xx1p-1.)*(xx1p+0.)*(xx1p+2)/6.
+ coef1p=coef1p+(xx1p-1.)*(xx1p-0.)*(xx1p+1.)*(xx1p+3)/6.
+ coef2p=(xx2p-2.)*(xx2p-1.)*(xx2p+0.)*(xx2p+1)/24.
+ coef2p=coef2p-(xx2p-1.)*(xx2p-0.)*(xx2p+1.)*(xx2p+2)/24.
+ SC(ip1p,jj,kk)=SC(ip1p,jj,kk)+u1m*coef1m+u0*coef0+u1p*coef1p+u2p*coef2p
+
+ coef1m=(xx1m-2.)*(xx1m-1.)*(xx1m+0.)*(xx1m+1)/24.
+ coef0=-(xx0-3.)*(xx0-1.)*(xx0+0.)*(xx0+1)/6.
+ coef1p=(xx1p-2.)*(xx1p-1.)*(xx1p+1.)*(xx1p+2)/4.
+ coef2p=-(xx2p-2.)*(xx2p-1.)*(xx2p+0.)*(xx2p+2)/6.
+ SC(ip2p,jj,kk)=SC(ip2p,jj,kk)+u1m*coef1m+u0*coef0+u1p*coef1p+u2p*coef2p
+
+ coef0=(xx0-2.)*(xx0-1.)*(xx0+0.)*(xx0+1)/24.
+ coef1p=-(xx1p-2.)*(xx1p-0.)*(xx1p+1.)*(xx1p+2)/6.
+ coef2p=(xx2p-2.)*(xx2p-1.)*(xx2p+1.)*(xx2p+2)/4.
+ SC(ip3p,jj,kk)=SC(ip3p,jj,kk)+u0*coef0+u1p*coef1p+u2p*coef2p
+
+ coef1p=(xx1p-1.)*(xx1p-0.)*(xx1p+1.)*(xx1p+2)/24.
+ coef2p=-(xx2p-2.)*(xx2p-0.)*(xx2p+1.)*(xx2p+2)/6.
+ SC(ip4p,jj,kk)=SC(ip4p,jj,kk)+u1p*coef1p+u2p*coef2p
+
+ coef2p=(xx2p-1.)*(xx2p-0.)*(xx2p+1.)*(xx2p+2)/24.
+ SC(ip5p,jj,kk)=SC(ip5p,jj,kk)+u2p*coef2p
+
+ else
+
+ coef1m=(xx1m-2.)*(xx1m-1.)*(xx1m+0.)*(xx1m+1)/24.
+ SC(ip3m,jj,kk)=SC(ip3m,jj,kk)+u1m*coef1m
+
+ coef1m=-(xx1m-2.)*(xx1m-1.)*(xx1m+0.)*(xx1m+2)/6.
+ coef0=(xx0-2.)*(xx0-1.)*xx0*(xx0+1)/24.
+ SC(ip2m,jj,kk)=SC(ip2m,jj,kk)+u1m*coef1m+u0*coef0
+
+ coef1m=(xx1m-2.)*(xx1m-1.)*(xx1m+1.)*(xx1m+2)/4.
+ coef0=-(xx0-2.)*(xx0-1.)*(xx0+0.)*(xx0+2)/6.
+ coef1p=(xx1p-3.)*(xx1p-2.)*(xx1p-1.)*(xx1p+0)/24.
+ SC(ip1m,jj,kk)=SC(ip1m,jj,kk)+u1m*coef1m+u0*coef0+u1p*coef1p
+
+ coef1m=-(xx1m-2.)*(xx1m-0.)*(xx1m+1.)*(xx1m+2)/6.
+ coef1m=coef1m+(xx1m-1.)*(xx1m-0.)*(xx1m+1.)*(xx1m+2)/24.
+ coef0=(xx0-2.)*(xx0-1.)*(xx0+1.)*(xx0+2)/4.
+ coef0=coef0-(xx0-2.)*(xx0-0.)*(xx0+1.)*(xx0+2)/6.
+ coef0=coef0+(xx0-1.)*(xx0-0.)*(xx0+1.)*(xx0+2)/24.
+ coef0=coef0-(xx0-0.)*(xx0+1.)*(xx0+2.)*(xx0+3)/24.
+ coef1p=(xx1p-2.)*(xx1p-1.)*(xx1p+0.)*(xx1p+1)/24.
+ coef1p=coef1p-(xx1p-3.)*(xx1p-2.)*(xx1p-1.)*(xx1p+0)/24.
+ coef1p=coef1p-(xx1p-2.)*(xx1p-1.)*(xx1p+0.)*(xx1p+2)/6.
+ coef1p=coef1p+(xx1p-2.)*(xx1p-1.)*(xx1p+1.)*(xx1p+2)/4.
+ coef2p=(xx2p-2.)*(xx2p-1.)*(xx2p+0.)*(xx2p+1)/24.
+ coef2p=coef2p-(xx2p-2.)*(xx2p-1.)*(xx2p+0.)*(xx2p+2)/6.
+ SC(ip0,jj,kk)=SC(ip0,jj,kk)+u1m*coef1m+u0*coef0+u1p*coef1p+u2p*coef2p
+
+ coef0=(xx0-0.)*(xx0+1.)*(xx0+2.)*(xx0+3)/24.
+ coef1p=-(xx1p-2.)*(xx1p-0.)*(xx1p+1.)*(xx1p+2)/6.
+ coef2p=(xx2p-2.)*(xx2p-1.)*(xx2p+1.)*(xx2p+2)/4.
+ SC(ip1p,jj,kk)=SC(ip1p,jj,kk)+u0*coef0+u1p*coef1p+u2p*coef2p
+
+ coef1p=(xx1p-1.)*(xx1p-0.)*(xx1p+1.)*(xx1p+2)/24.
+ coef2p=-(xx2p-2.)*(xx2p-0.)*(xx2p+1.)*(xx2p+2)/6.
+ SC(ip2p,jj,kk)=SC(ip2p,jj,kk)+u1p*coef1p+u2p*coef2p
+
+ coef2p=(xx2p-1.)*(xx2p-0.)*(xx2p+1.)*(xx2p+2)/24.
+ SC(ip3p,jj,kk)=SC(ip3p,jj,kk)+u2p*coef2p
+
+ endif
+ enddo
+
+end subroutine remeshx_tag
diff --git a/LEGI/AutresRemaillage/advec_O4v2_y.f90 b/LEGI/AutresRemaillage/advec_O4v2_y.f90
new file mode 100644
index 000000000..76a152ccf
--- /dev/null
+++ b/LEGI/AutresRemaillage/advec_O4v2_y.f90
@@ -0,0 +1,445 @@
+subroutine y_advec
+
+ use x_advec_m
+
+ implicit none
+
+ real(WP) :: dx, yy, zz, tm
+ integer :: i,j,k,np,n, ir, kr, ini, np_aux, ntag
+ integer :: jj, j1, j2, j3, j4, ic, m
+
+ ntag = 0
+ ntag_total = 0
+ np_aux = 0
+ tm = 0.
+ m = 4 !np_bl+1 (et np_bl=1)
+
+ dx = L/nxsc
+ cfl = dt/dx
+
+ do k=1,nxsc
+ zz=xmin+float(k-1)*dx
+ do i=1,nxsc
+ np=0
+ yy=xmin+float(i-1)*dx
+ np=0
+ jj=3
+ ic=0
+ do j=3,nxsc+2
+ jj=jj+ic
+ if (jj.le.nxsc+2) then
+ j1=mod(jj-1,nxsc)+1
+ j2=mod(jj,nxsc)+1
+ j3=mod(jj+1,nxsc)+1
+ j4=mod(jj+2,nxsc)+1
+ if ((mod(j1,m).eq.m-1).and. &
+ & (abs(SC(i,j1,k))+abs(SC(i,j2,k))+abs(SC(i,j3,k))+abs(SC(i,j4,k)).ge.circlim)) then
+ np=np+1
+ up(np)=SC(i,j1,k)
+ tm=tm+SC(i,j1,k)
+ xp(np)=xmin+float(j1-1)*dx
+ np=np+1
+ up(np)=SC(i,j2,k)
+ tm=tm+SC(i,j2,k)
+ xp(np)=xmin+float(j2-1)*dx
+ np=np+1
+ up(np)=SC(i,j3,k)
+ tm=tm+SC(i,j3,k)
+ xp(np)=xmin+float(j3-1)*dx
+ np=np+1
+ up(np)=SC(i,j4,k)
+ tm=tm+SC(i,j4,k)
+ xp(np)=xmin+float(j4-1)*dx
+ ic=4
+ elseif ((abs(SC(i,j1,k)).ge.circlim)) then
+ np=np+1
+ up(np)=SC(i,j1,k)
+ tm=tm+SC(i,j1,k)
+ xp(np)=xmin+float(j1-1)*dx
+ ic=1
+ else
+ ic=1
+ endif
+ endif
+ enddo
+ if (np.ne.0) then
+ call velox_y(np,i,k)
+ do n=1,np
+ xp0(n)=xp(n)
+ xp(n)=xp(n)+0.5*dt*vx(n)
+ if (xp(n).gt.xmax) xp(n)=xp(n)-xmax+xmin
+ if (xp(n).lt.xmin) xp(n)=xp(n)+xmax-xmin
+ enddo
+ call velox_y(np,i,k)
+
+!! debut modif
+
+ call tag_particles(np,np_aux,ntag)
+ ntag_total = ntag_total + ntag
+ ! test
+ !itype(1:np_aux) = 1
+
+! np integer: deja defini (input) - local
+! np_aux integer: deja defini (output) - local
+! ntag integer: deja defini (output) - local
+! np_bl=1 --> supprimer de l'appel
+! icfl tablo integer: deja defini - common --> supprimer de l'appel
+! itype tablo integer: deja defini - common --> supprimer de l'appel
+! itype_aux tablo integer: deja defini - common --> supprimer de l'appel
+! itag tablo integer: deja defini - common --> supprimer de l'appel
+! xp_aux tablo real integer: deja defini - common --> supprimer de l'appel
+! up_aux tablo real integer: deja defini - common --> supprimer de l'appel
+! vx_aux tablo real integer: deja defini - common --> supprimer de l'appel
+
+! ---- > trouver les dimension max des tablo
+! ---- > trouver si les tablo doivent etre declarer en local ou en common
+
+ if (ntag.ne.0) then
+ do n=1,ntag
+ ini=itag(n)
+ xp(ini)=xp0(ini)+dt*vx(ini)
+ if (xp(ini).gt.xmax) xp(ini)=xp(ini)-xmax+xmin
+ if (xp(ini).lt.xmin) xp(ini)=xp(ini)+xmax-xmin !! TO CHECK WITH GH
+ end do
+ end if
+
+ do n=1,np_aux
+ xp_aux(n)=xp_aux(n)+dt*vx_aux(n)
+ if (xp_aux(n).gt.xmax) xp_aux(n)=xp_aux(n)-xmax+xmin
+ if (xp_aux(n).lt.xmin) xp_aux(n)=xp_aux(n)+xmax-xmin
+ end do
+ ir=i
+ kr=k
+ call remeshy(np_aux,ir,kr)
+ if (ntag.ne.0) call remeshy_tag(ntag,ir,kr)
+
+!! fin modif
+
+ endif
+ enddo
+ enddo
+
+ print*,'ntag y_advec :', ntag_total
+
+
+end subroutine y_advec
+
+
+
+subroutine velox_y(np,j,k)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: np, j, k, i, nx3, ny3, nz3, jp1, jp2, kp1, kp2, ip1, ip2
+ real(WP) :: dx3, dy3, dz3, dxinv, yy, zz, x0, y0, z0, dx2
+ real(WP) :: yy1, yy2, zz1, zz2, b1, b2, c1, c2, a1, a2
+ real(WP) :: x, xx1, xx2
+
+ do i=1,np
+ vx(i)=0.
+ end do
+
+ nx3=nxsc
+ ny3=nx3
+ nz3=nx3
+ dx3=xmax/float(nx3)
+ dy3=dx3
+ dz3=dx3
+
+ dxinv=1./dx3
+
+ x0=xmin
+ y0=ymin
+ z0=zmin
+
+ dx2 = L/nxsc
+
+ yy=+float(j-1)*dx2
+ zz=+float(k-1)*dx2
+
+ jp1 = int((yy)/dx3)
+ jp2 = jp1 + 1
+ kp1 = int((zz)/dx3)
+ kp2 = kp1 + 1
+ yy1 = (yy-float(jp1)*dx3)/dx3
+ yy2 = 1.0 - yy1
+ zz1 = (zz-float(kp1)*dx3)/dx3
+ zz2 = 1.0 - zz1
+ b1=yy2
+ b2=yy1
+ c1=zz2
+ c2=zz1
+
+ jp1=mod(jp1+nx3,nx3) +1
+ jp2=mod(jp2+nx3,nx3) +1
+ kp1=mod(kp1+nx3,nx3) +1
+ kp2=mod(kp2+nx3,nx3) +1
+
+ do i = 1,np
+
+ x = xp(i)
+
+ ip1 = int((x-x0)*dxinv)
+ ip2 = ip1 + 1
+
+ xx1 = (x - float(ip1)*dx3-x0)*dxinv
+ xx2=1-xx1
+
+ ip1=mod(ip1+nx3,nx3) +1
+ ip2=mod(ip2+nx3,nx3) +1
+
+ a1 = xx2
+ a2 = xx1
+
+ vx(i)= vx(i) + Vsc(jp1,ip1,kp1)*a1*b1*c1
+ vx(i)= vx(i) + Vsc(jp1,ip2,kp1)*a2*b1*c1
+ vx(i)= vx(i) + Vsc(jp2,ip1,kp1)*a1*b2*c1
+ vx(i)= vx(i) + Vsc(jp2,ip2,kp1)*a2*b2*c1
+ vx(i)= vx(i) + Vsc(jp1,ip1,kp2)*a1*b1*c2
+ vx(i)= vx(i) + Vsc(jp1,ip2,kp2)*a2*b1*c2
+ vx(i)= vx(i) + Vsc(jp2,ip1,kp2)*a1*b2*c2
+ vx(i)= vx(i) + Vsc(jp2,ip2,kp2)*a2*b2*c2
+
+ end do
+
+end subroutine velox_y
+
+subroutine remeshy(np1,jj,kk)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: np1, jj, kk
+ integer :: i, nx2, n, ip0, ip1, ip2, k, j, ip3, ip4
+ real(WP) :: dxinv, dx2, x0, g1, x, xx0, xx1, xx2, a0, a1, a2, xx3, xx4, a3, a4
+
+ nx2 = nxsc
+ dx2 = L/nx2
+
+ dxinv=1./dx2
+
+ do i=1,nx2
+ SC(jj,i,kk)=0.
+ enddo
+
+ x0=xmin
+
+ do n = 1,np1
+ g1 = up_aux(n)
+ x = xp_aux(n)
+
+ if (itype_aux(n).eq.0) then
+ ip1 = int((x-x0)*dxinv)
+ else
+ ip1 = nint((x-x0)*dxinv)
+ endif
+
+ ip0 = ip1 - 1
+ ip2 = ip1 + 1
+ ip3 =ip1-2
+ ip4=ip1+2
+
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=xx1-1
+ xx3=xx1+2
+ xx4=xx1-2
+
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ ip3=mod(ip3+nx2,nx2) +1
+ ip4=mod(ip4+nx2,nx2) +1
+
+ a0=-xx1*xx2*xx3*xx4/6.
+ a1=xx0*xx2*xx3*xx4/4.
+ a2=-xx0*xx1*xx3*xx4/6.
+ a3=xx0*xx1*xx2*xx4/24.
+ a4=xx0*xx1*xx2*xx3/24.
+
+ SC(jj,ip0,kk) = SC(jj,ip0,kk) + g1*a0
+ SC(jj,ip1,kk) = SC(jj,ip1,kk) + g1*a1
+ SC(jj,ip2,kk) = SC(jj,ip2,kk) + g1*a2
+ SC(jj,ip3,kk) = SC(jj,ip3,kk) + g1*a3
+ SC(jj,ip4,kk) = SC(jj,ip4,kk) + g1*a4
+
+ end do
+
+!!! go to 222 !! -> ca veut dire quoi goto ?
+!!! do k=3,nx2-2
+!!! do j=1,nx2
+!!! do i=1,nx2
+!!! if (SC(i,j,k).gt.1.2) then
+!!! print*,'BOUM', i,j,k,SC(i,j,k)
+!!! goto 222
+!!! end if
+!!! end do
+!!! end do
+!!! end do
+
+end subroutine remeshy
+
+subroutine remeshy_tag(ntag,jj,kk)
+
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: jj, kk, ntag
+ integer :: i, nx2, n, ip1, ip2, k, j, ini, inni, innni, jp1, ip3, ip4
+ integer :: ip1p, ip2p, ip0, ip1m
+ integer :: ip3m, ip2m, ip5p, ip4p, ip3p
+ real(WP) :: dxinv, dx2, x0, g1, x, xx0, xx1m, xx1p, xx2p, xx2, a0, a1, a2, b1, b2, y, u1, u2, yy0, yy1, yy2
+ real(WP) :: x1p, x2p, u1p, u2p, x1m, u1m, u0
+ real(WP) :: coef1m, coef0, coef1p, coef2p
+
+ nx2 = nxsc
+ dx2 = L/nx2
+
+ dxinv=1./dx2
+
+ do n=1,ntag-1,4
+ i=itag(n)
+ ini=itag(n+1)
+ inni=itag(n+2)
+ innni=itag(n+3)
+ x1m=xp(i)
+ x0=xp(ini)
+ x1p=xp(inni)
+ x2p=xp(innni)
+ u1m=up(i)
+ u0=up(ini)
+ u1p=up(inni)
+ u2p=up(innni)
+ if (itype(i).eq.0) then
+ ip0 = int((x0-xmin)*dxinv)
+ ip1m = int((x1m-xmin)*dxinv)
+ ip1p = nint((x1p-xmin)*dxinv)
+ ip2p = nint((x2p-xmin)*dxinv)
+ else
+ ip0 = nint((x0-xmin)*dxinv)
+ ip1m = nint((x1m-xmin)*dxinv)
+ ip1p = int((x1p-xmin)*dxinv)
+ ip2p = int((x2p-xmin)*dxinv)
+ end if
+ xx0 = (x0 - float(ip0)*dx2-xmin)*dxinv
+ xx1m = (x1m - float(ip1m)*dx2-xmin)*dxinv
+ xx1p = (x1p - float(ip1p)*dx2-xmin)*dxinv
+ xx2p = (x2p - float(ip2p)*dx2-xmin)*dxinv
+ ip3m=ip0-3
+ ip2m=ip0-2
+ ip1m=ip0-1
+ ip5p=ip0+5
+ ip4p=ip0+4
+ ip3p=ip0+3
+ ip2p=ip0+2
+ ip1p=ip0+1
+ ip0=mod(ip0+nxsc,nxsc) +1
+ ip1p=mod(ip1p+nxsc,nxsc) +1
+ ip2p=mod(ip2p+nxsc,nxsc) +1
+ ip3p=mod(ip3p+nxsc,nxsc) +1
+ ip4p=mod(ip4p+nxsc,nxsc) +1
+ ip5p=mod(ip5p+nxsc,nxsc) +1
+ ip1m=mod(ip1m+nxsc,nxsc) +1
+ ip2m=mod(ip2m+nxsc,nxsc) +1
+ ip3m=mod(ip3m+nxsc,nxsc) +1
+
+ if (itype(i).eq.0) then
+ coef1m=(xx1m-2.)*(xx1m-1.)*xx1m*(xx1m+1)/24.
+ SC(jj,ip3m,kk)=SC(jj,ip3m,kk)+u1m*coef1m
+
+ coef1m=-(xx1m-2.)*(xx1m-1.)*xx1m*(xx1m+2)/6.
+ coef0=(xx0-2.)*(xx0-1.)*xx0*(xx0+1)/24.
+ SC(jj,ip2m,kk)=SC(jj,ip2m,kk)+u1m*coef1m+u0*coef0
+
+ coef1m=(xx1m-2.)*(xx1m-1.)*(xx1m+1.)*(xx1m+2)/4.
+ coef0=-(xx0-2.)*(xx0-1.)*xx0*(xx0+2)/6.
+ coef1p=(xx1p-1.)*(xx1p-0.)*(xx1p+1.)*(xx1p+2)/24.
+ SC(jj,ip1m,kk)=SC(jj,ip1m,kk)+u1m*coef1m+u0*coef0+u1p*coef1p
+
+ coef1m=-(xx1m-2.)*(xx1m-0.)*(xx1m+1.)*(xx1m+2)/6.
+ coef0=(xx0-2.)*(xx0-1.)*(xx0+1.)*(xx0+2)/4.
+ coef1p=-(xx1p-1.)*(xx1p-0.)*(xx1p+1.)*(xx1p+3)/6.
+ coef2p=(xx2p-1.)*(xx2p-0.)*(xx2p+1.)*(xx2p+2)/24.
+ SC(jj,ip0,kk)=SC(jj,ip0,kk)+u1m*coef1m+u0*coef0+u1p*coef1p+u2p*coef2p
+
+ coef1m=(xx1m-1.)*(xx1m-0.)*(xx1m+1.)*(xx1m+2)/24.
+ coef1m=coef1m-(xx1m-2.)*(xx1m-1.)*(xx1m+0.)*(xx1m+1)/24.
+ coef0=-(xx0-2.)*(xx0-0.)*(xx0+1.)*(xx0+2)/6.
+ coef0=coef0+(xx0-3.)*(xx0-1.)*(xx0+0.)*(xx0+1)/6.
+ coef0=coef0+(xx0-1.)*(xx0-0.)*(xx0+1.)*(xx0+2)/24.
+ coef0=coef0-(xx0-2.)*(xx0-1.)*(xx0+0.)*(xx0+1)/24.
+ coef1p=(xx1p-2.)*(xx1p-1.)*(xx1p+0.)*(xx1p+1)/24.
+ coef1p=coef1p-(xx1p-1.)*(xx1p-0.)*(xx1p+1.)*(xx1p+2)/24.
+ coef1p=coef1p-(xx1p-2.)*(xx1p-1.)*(xx1p+0.)*(xx1p+2)/6.
+ coef1p=coef1p+(xx1p-1.)*(xx1p-0.)*(xx1p+1.)*(xx1p+3)/6.
+ coef2p=(xx2p-2.)*(xx2p-1.)*(xx2p+0.)*(xx2p+1)/24.
+ coef2p=coef2p-(xx2p-1.)*(xx2p-0.)*(xx2p+1.)*(xx2p+2)/24.
+ SC(jj,ip1p,kk)=SC(jj,ip1p,kk)+u1m*coef1m+u0*coef0+u1p*coef1p+u2p*coef2p
+
+ coef1m=(xx1m-2.)*(xx1m-1.)*(xx1m+0.)*(xx1m+1)/24.
+ coef0=-(xx0-3.)*(xx0-1.)*(xx0+0.)*(xx0+1)/6.
+ coef1p=(xx1p-2.)*(xx1p-1.)*(xx1p+1.)*(xx1p+2)/4.
+ coef2p=-(xx2p-2.)*(xx2p-1.)*(xx2p+0.)*(xx2p+2)/6.
+ SC(jj,ip2p,kk)=SC(jj,ip2p,kk)+u1m*coef1m+u0*coef0+u1p*coef1p+u2p*coef2p
+
+ coef0=(xx0-2.)*(xx0-1.)*(xx0+0.)*(xx0+1)/24.
+ coef1p=-(xx1p-2.)*(xx1p-0.)*(xx1p+1.)*(xx1p+2)/6.
+ coef2p=(xx2p-2.)*(xx2p-1.)*(xx2p+1.)*(xx2p+2)/4.
+ SC(jj,ip3p,kk)=SC(jj,ip3p,kk)+u0*coef0+u1p*coef1p+u2p*coef2p
+
+ coef1p=(xx1p-1.)*(xx1p-0.)*(xx1p+1.)*(xx1p+2)/24.
+ coef2p=-(xx2p-2.)*(xx2p-0.)*(xx2p+1.)*(xx2p+2)/6.
+ SC(jj,ip4p,kk)=SC(jj,ip4p,kk)+u1p*coef1p+u2p*coef2p
+
+ coef2p=(xx2p-1.)*(xx2p-0.)*(xx2p+1.)*(xx2p+2)/24.
+ SC(jj,ip5p,kk)=SC(jj,ip5p,kk)+u2p*coef2p
+
+ else
+
+ coef1m=(xx1m-2.)*(xx1m-1.)*(xx1m+0.)*(xx1m+1)/24.
+ SC(jj,ip3m,kk)=SC(jj,ip3m,kk)+u1m*coef1m
+
+ coef1m=-(xx1m-2.)*(xx1m-1.)*(xx1m+0.)*(xx1m+2)/6.
+ coef0=(xx0-2.)*(xx0-1.)*xx0*(xx0+1)/24.
+ SC(jj,ip2m,kk)=SC(jj,ip2m,kk)+u1m*coef1m+u0*coef0
+
+ coef1m=(xx1m-2.)*(xx1m-1.)*(xx1m+1.)*(xx1m+2)/4.
+ coef0=-(xx0-2.)*(xx0-1.)*(xx0+0.)*(xx0+2)/6.
+ coef1p=(xx1p-3.)*(xx1p-2.)*(xx1p-1.)*(xx1p+0)/24.
+ SC(jj,ip1m,kk)=SC(jj,ip1m,kk)+u1m*coef1m+u0*coef0+u1p*coef1p
+
+ coef1m=-(xx1m-2.)*(xx1m-0.)*(xx1m+1.)*(xx1m+2)/6.
+ coef1m=coef1m+(xx1m-1.)*(xx1m-0.)*(xx1m+1.)*(xx1m+2)/24.
+ coef0=(xx0-2.)*(xx0-1.)*(xx0+1.)*(xx0+2)/4.
+ coef0=coef0-(xx0-2.)*(xx0-0.)*(xx0+1.)*(xx0+2)/6.
+ coef0=coef0+(xx0-1.)*(xx0-0.)*(xx0+1.)*(xx0+2)/24.
+ coef0=coef0-(xx0-0.)*(xx0+1.)*(xx0+2.)*(xx0+3)/24.
+ coef1p=(xx1p-2.)*(xx1p-1.)*(xx1p+0.)*(xx1p+1)/24.
+ coef1p=coef1p-(xx1p-3.)*(xx1p-2.)*(xx1p-1.)*(xx1p+0)/24.
+ coef1p=coef1p-(xx1p-2.)*(xx1p-1.)*(xx1p+0.)*(xx1p+2)/6.
+ coef1p=coef1p+(xx1p-2.)*(xx1p-1.)*(xx1p+1.)*(xx1p+2)/4.
+ coef2p=(xx2p-2.)*(xx2p-1.)*(xx2p+0.)*(xx2p+1)/24.
+ coef2p=coef2p-(xx2p-2.)*(xx2p-1.)*(xx2p+0.)*(xx2p+2)/6.
+ SC(jj,ip0,kk)=SC(jj,ip0,kk)+u1m*coef1m+u0*coef0+u1p*coef1p+u2p*coef2p
+
+ coef0=(xx0-0.)*(xx0+1.)*(xx0+2.)*(xx0+3)/24.
+ coef1p=-(xx1p-2.)*(xx1p-0.)*(xx1p+1.)*(xx1p+2)/6.
+ coef2p=(xx2p-2.)*(xx2p-1.)*(xx2p+1.)*(xx2p+2)/4.
+ SC(jj,ip1p,kk)=SC(jj,ip1p,kk)+u0*coef0+u1p*coef1p+u2p*coef2p
+
+ coef1p=(xx1p-1.)*(xx1p-0.)*(xx1p+1.)*(xx1p+2)/24.
+ coef2p=-(xx2p-2.)*(xx2p-0.)*(xx2p+1.)*(xx2p+2)/6.
+ SC(jj,ip2p,kk)=SC(jj,ip2p,kk)+u1p*coef1p+u2p*coef2p
+
+ coef2p=(xx2p-1.)*(xx2p-0.)*(xx2p+1.)*(xx2p+2)/24.
+ SC(jj,ip3p,kk)=SC(jj,ip3p,kk)+u2p*coef2p
+
+ endif
+ enddo
+
+end subroutine remeshy_tag
diff --git a/LEGI/AutresRemaillage/advec_O4v2_z.f90 b/LEGI/AutresRemaillage/advec_O4v2_z.f90
new file mode 100644
index 000000000..0165c1cb4
--- /dev/null
+++ b/LEGI/AutresRemaillage/advec_O4v2_z.f90
@@ -0,0 +1,447 @@
+subroutine z_advec
+
+ use x_advec_m
+
+ implicit none
+
+ real(WP) :: dx, yy, zz, tm
+ integer :: i,j,k,np,n, ir, jr, ini, ntag, np_aux
+ integer :: jj, j1, j2, j3, j4, ic, m
+
+
+
+ dx = L/nxsc
+ cfl = dt/dx
+
+ ntag = 0
+ ntag_total = 0
+ np_aux = 0
+ tm = 0.
+ m = 4 !np_bl+1 (et np_bl=1)
+
+
+ do j=1,nxsc
+ zz=xmin+float(k-1)*dx
+ do i=1,nxsc
+ np=0
+ yy=xmin+float(i-1)*dx
+ np=0
+ jj=3
+ ic=0
+ do k=3,nxsc+2
+ jj=jj+ic
+ if (jj.le.nxsc+2) then
+ j1=mod(jj-1,nxsc)+1
+ j2=mod(jj,nxsc)+1
+ j3=mod(jj+1,nxsc)+1
+ j4=mod(jj+2,nxsc)+1
+ if ((mod(j1,m).eq.m-1).and. &
+ & (abs(SC(i,j,j1))+abs(SC(i,j,j2))+abs(SC(i,j,j3))+abs(SC(i,j,j4)).ge.circlim)) then
+ np=np+1
+ up(np)=SC(i,j,j1)
+ tm=tm+SC(i,j,j1)
+ xp(np)=xmin+float(j1-1)*dx
+ np=np+1
+ up(np)=SC(i,j,j2)
+ tm=tm+SC(i,j,j2)
+ xp(np)=xmin+float(j2-1)*dx
+ np=np+1
+ up(np)=SC(i,j,j3)
+ tm=tm+SC(i,j,j3)
+ xp(np)=xmin+float(j3-1)*dx
+ np=np+1
+ up(np)=SC(i,j,j4)
+ tm=tm+SC(i,j,j4)
+ xp(np)=xmin+float(j4-1)*dx
+ ic=4
+ elseif ((abs(SC(i,j,j1)).ge.circlim)) then
+ np=np+1
+ up(np)=SC(i,j,j1)
+ tm=tm+SC(i,j,j1)
+ xp(np)=xmin+float(j1-1)*dx
+ ic=1
+ else
+ ic=1
+ endif
+ endif
+ enddo
+ if (np.ne.0) then
+ call velox_z(np,i,j)
+ do n=1,np
+ xp0(n)=xp(n)
+ xp(n)=xp(n)+0.5*dt*vx(n)
+ if (xp(n).gt.xmax) xp(n)=xp(n)-xmax+xmin
+ if (xp(n).lt.xmin) xp(n)=xp(n)+xmax-xmin
+ enddo
+ call velox_z(np,i,j)
+
+!! debut modif
+
+ call tag_particles(np,np_aux,ntag)
+ ntag_total = ntag_total + ntag
+ ! test
+ !itype(1:np_aux) = 1
+
+! np integer: deja defini (input) - local
+! np_aux integer: deja defini (output) - local
+! ntag integer: deja defini (output) - local
+! np_bl=1 --> supprimer de l'appel
+! icfl tablo integer: deja defini - common --> supprimer de l'appel
+! itype tablo integer: deja defini - common --> supprimer de l'appel
+! itype_aux tablo integer: deja defini - common --> supprimer de l'appel
+! itag tablo integer: deja defini - common --> supprimer de l'appel
+! xp_aux tablo real integer: deja defini - common --> supprimer de l'appel
+! up_aux tablo real integer: deja defini - common --> supprimer de l'appel
+! vx_aux tablo real integer: deja defini - common --> supprimer de l'appel
+
+! ---- > trouver les dimension max des tablo
+! ---- > trouver si les tablo doivent etre declarer en local ou en common
+
+ if (ntag.ne.0) then
+ do n=1,ntag
+ ini=itag(n)
+ xp(ini)=xp0(ini)+dt*vx(ini)
+ if (xp(ini).gt.xmax) xp(ini)=xp(ini)-xmax+xmin
+ if (xp(ini).lt.xmin) xp(ini)=xp(ini)+xmax-xmin !! TO CHECK WITH GH
+ end do
+ end if
+
+ do n=1,np_aux
+ xp_aux(n)=xp_aux(n)+dt*vx_aux(n)
+ if (xp_aux(n).gt.xmax) xp_aux(n)=xp_aux(n)-xmax+xmin
+ if (xp_aux(n).lt.xmin) xp_aux(n)=xp_aux(n)+xmax-xmin
+ end do
+ ir=i
+ jr=j
+ call remeshz(np_aux,ir,jr)
+ if (ntag.ne.0) call remeshz_tag(ntag,ir,jr)
+
+!! fin modif
+
+ endif
+ enddo
+ enddo
+
+ print*,'ntag z_advec :', ntag_total
+
+end subroutine z_advec
+
+
+
+subroutine velox_z(np,j,k)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: np, j, k, i, nx3, ny3, nz3, jp1, jp2, kp1, kp2, ip1, ip2
+ real(WP) :: dx3, dy3, dz3, dxinv, yy, zz, x0, y0, z0, dx2
+ real(WP) :: yy1, yy2, zz1, zz2, b1, b2, c1, c2, a1, a2
+ real(WP) :: x, xx1, xx2
+
+ do i=1,np
+ vx(i)=0.
+ end do
+
+ nx3=nxsc
+ ny3=nx3
+ nz3=nx3
+ dx3=xmax/float(nx3)
+ dy3=dx3
+ dz3=dx3
+
+ dxinv=1./dx3
+
+ x0=xmin
+ y0=ymin
+ z0=zmin
+
+ dx2 = L/nxsc
+
+ yy=+float(j-1)*dx2
+ zz=+float(k-1)*dx2
+
+ jp1 = int((yy)/dx3)
+ jp2 = jp1 + 1
+ kp1 = int((zz)/dx3)
+ kp2 = kp1 + 1
+ yy1 = (yy-float(jp1)*dx3)/dx3
+ yy2 = 1.0 - yy1
+ zz1 = (zz-float(kp1)*dx3)/dx3
+ zz2 = 1.0 - zz1
+ b1=yy2
+ b2=yy1
+ c1=zz2
+ c2=zz1
+
+ jp1=mod(jp1+nx3,nx3) +1
+ jp2=mod(jp2+nx3,nx3) +1
+ kp1=mod(kp1+nx3,nx3) +1
+ kp2=mod(kp2+nx3,nx3) +1
+
+ do i = 1,np
+
+ x = xp(i)
+
+ ip1 = int((x-x0)*dxinv)
+ ip2 = ip1 + 1
+
+ xx1 = (x - float(ip1)*dx3-x0)*dxinv
+ xx2=1-xx1
+
+ ip1=mod(ip1+nx3,nx3) +1
+ ip2=mod(ip2+nx3,nx3) +1
+
+ a1 = xx2
+ a2 = xx1
+
+ vx(i)= vx(i) + Wsc(jp1,kp1,ip1)*a1*b1*c1
+ vx(i)= vx(i) + Wsc(jp1,kp1,ip2)*a2*b1*c1
+ vx(i)= vx(i) + Wsc(jp2,kp1,ip1)*a1*b2*c1
+ vx(i)= vx(i) + Wsc(jp2,kp1,ip2)*a2*b2*c1
+ vx(i)= vx(i) + Wsc(jp1,kp2,ip1)*a1*b1*c2
+ vx(i)= vx(i) + Wsc(jp1,kp2,ip2)*a2*b1*c2
+ vx(i)= vx(i) + Wsc(jp2,kp2,ip1)*a1*b2*c2
+ vx(i)= vx(i) + Wsc(jp2,kp2,ip2)*a2*b2*c2
+
+ end do
+
+end subroutine velox_z
+
+subroutine remeshz(np1,jj,kk)
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: np1, jj, kk
+ integer :: i, nx2, n, ip0, ip1, ip2, k, j, ip3, ip4
+ real(WP) :: dxinv, dx2, x0, g1, x, xx0, xx1, xx2, a0, a1, a2, xx3, xx4, a3, a4
+
+ nx2 = nxsc
+ dx2 = L/nx2
+
+ dxinv=1./dx2
+
+ do i=1,nx2
+ SC(jj,kk,i)=0.
+ enddo
+
+ x0=xmin
+
+ do n = 1,np1
+ g1 = up_aux(n)
+ x = xp_aux(n)
+
+ if (itype_aux(n).eq.0) then
+ ip1 = int((x-x0)*dxinv)
+ else
+ ip1 = nint((x-x0)*dxinv)
+ endif
+
+ ip0 = ip1 - 1
+ ip2 = ip1 + 1
+ ip3 =ip1-2
+ ip4=ip1+2
+
+ xx1 = (x - float(ip1)*dx2-x0)*dxinv
+ xx0=xx1+1
+ xx2=xx1-1
+ xx3=xx1+2
+ xx4=xx1-2
+
+ ip1=mod(ip1+nx2,nx2) +1
+ ip0=mod(ip0+nx2,nx2) +1
+ ip2=mod(ip2+nx2,nx2) +1
+ ip3=mod(ip3+nx2,nx2) +1
+ ip4=mod(ip4+nx2,nx2) +1
+
+ a0=-xx1*xx2*xx3*xx4/6.
+ a1=xx0*xx2*xx3*xx4/4.
+ a2=-xx0*xx1*xx3*xx4/6.
+ a3=xx0*xx1*xx2*xx4/24.
+ a4=xx0*xx1*xx2*xx3/24.
+
+ SC(jj,kk,ip0) = SC(jj,kk,ip0) + g1*a0
+ SC(jj,kk,ip1) = SC(jj,kk,ip1) + g1*a1
+ SC(jj,kk,ip2) = SC(jj,kk,ip2) + g1*a2
+ SC(jj,kk,ip3) = SC(jj,kk,ip3) + g1*a3
+ SC(jj,kk,ip4) = SC(jj,kk,ip4) + g1*a4
+
+ end do
+
+!!! go to 222 !! -> ca veut dire quoi goto ?
+!!! do k=3,nx2-2
+!!! do j=1,nx2
+!!! do i=1,nx2
+!!! if (SC(i,j,k).gt.1.2) then
+!!! print*,'BOUM', i,j,k,SC(i,j,k)
+!!! goto 222
+!!! end if
+!!! end do
+!!! end do
+!!! end do
+
+end subroutine remeshz
+
+subroutine remeshz_tag(ntag,kk,jj)
+
+
+ use x_advec_m
+
+ implicit none
+
+ integer :: jj, kk, ntag
+ integer :: i, nx2, n, ip1, ip2, k, j, ini, inni, innni, jp1, ip3, ip4
+ integer :: ip1p, ip2p, ip0, ip1m
+ integer :: ip3m, ip2m, ip5p, ip4p, ip3p
+ real(WP) :: dxinv, dx2, x0, g1, x, xx0, xx1m, xx1p, xx2p, xx2, a0, a1, a2, b1, b2, y, u1, u2, yy0, yy1, yy2
+ real(WP) :: x1p, x2p, u1p, u2p, x1m, u1m, u0
+ real(WP) :: coef1m, coef0, coef1p, coef2p
+
+ nx2 = nxsc
+ dx2 = L/nx2
+
+ dxinv=1./dx2
+
+ do n=1,ntag-1,4
+ i=itag(n)
+ ini=itag(n+1)
+ inni=itag(n+2)
+ innni=itag(n+3)
+ x1m=xp(i)
+ x0=xp(ini)
+ x1p=xp(inni)
+ x2p=xp(innni)
+ u1m=up(i)
+ u0=up(ini)
+ u1p=up(inni)
+ u2p=up(innni)
+ if (itype(i).eq.0) then
+ ip0 = int((x0-xmin)*dxinv)
+ ip1m = int((x1m-xmin)*dxinv)
+ ip1p = nint((x1p-xmin)*dxinv)
+ ip2p = nint((x2p-xmin)*dxinv)
+ else
+ ip0 = nint((x0-xmin)*dxinv)
+ ip1m = nint((x1m-xmin)*dxinv)
+ ip1p = int((x1p-xmin)*dxinv)
+ ip2p = int((x2p-xmin)*dxinv)
+ end if
+ xx0 = (x0 - float(ip0)*dx2-xmin)*dxinv
+ xx1m = (x1m - float(ip1m)*dx2-xmin)*dxinv
+ xx1p = (x1p - float(ip1p)*dx2-xmin)*dxinv
+ xx2p = (x2p - float(ip2p)*dx2-xmin)*dxinv
+ ip3m=ip0-3
+ ip2m=ip0-2
+ ip1m=ip0-1
+ ip5p=ip0+5
+ ip4p=ip0+4
+ ip3p=ip0+3
+ ip2p=ip0+2
+ ip1p=ip0+1
+ ip0=mod(ip0+nxsc,nxsc) +1
+ ip1p=mod(ip1p+nxsc,nxsc) +1
+ ip2p=mod(ip2p+nxsc,nxsc) +1
+ ip3p=mod(ip3p+nxsc,nxsc) +1
+ ip4p=mod(ip4p+nxsc,nxsc) +1
+ ip5p=mod(ip5p+nxsc,nxsc) +1
+ ip1m=mod(ip1m+nxsc,nxsc) +1
+ ip2m=mod(ip2m+nxsc,nxsc) +1
+ ip3m=mod(ip3m+nxsc,nxsc) +1
+
+ if (itype(i).eq.0) then
+ coef1m=(xx1m-2.)*(xx1m-1.)*xx1m*(xx1m+1)/24.
+ SC(kk,jj,ip3m)=SC(kk,jj,ip3m)+u1m*coef1m
+
+ coef1m=-(xx1m-2.)*(xx1m-1.)*xx1m*(xx1m+2)/6.
+ coef0=(xx0-2.)*(xx0-1.)*xx0*(xx0+1)/24.
+ SC(kk,jj,ip2m)=SC(kk,jj,ip2m)+u1m*coef1m+u0*coef0
+
+ coef1m=(xx1m-2.)*(xx1m-1.)*(xx1m+1.)*(xx1m+2)/4.
+ coef0=-(xx0-2.)*(xx0-1.)*xx0*(xx0+2)/6.
+ coef1p=(xx1p-1.)*(xx1p-0.)*(xx1p+1.)*(xx1p+2)/24.
+ SC(kk,jj,ip1m)=SC(kk,jj,ip1m)+u1m*coef1m+u0*coef0+u1p*coef1p
+
+ coef1m=-(xx1m-2.)*(xx1m-0.)*(xx1m+1.)*(xx1m+2)/6.
+ coef0=(xx0-2.)*(xx0-1.)*(xx0+1.)*(xx0+2)/4.
+ coef1p=-(xx1p-1.)*(xx1p-0.)*(xx1p+1.)*(xx1p+3)/6.
+ coef2p=(xx2p-1.)*(xx2p-0.)*(xx2p+1.)*(xx2p+2)/24.
+ SC(kk,jj,ip0)=SC(kk,jj,ip0)+u1m*coef1m+u0*coef0+u1p*coef1p+u2p*coef2p
+
+ coef1m=(xx1m-1.)*(xx1m-0.)*(xx1m+1.)*(xx1m+2)/24.
+ coef1m=coef1m-(xx1m-2.)*(xx1m-1.)*(xx1m+0.)*(xx1m+1)/24.
+ coef0=-(xx0-2.)*(xx0-0.)*(xx0+1.)*(xx0+2)/6.
+ coef0=coef0+(xx0-3.)*(xx0-1.)*(xx0+0.)*(xx0+1)/6.
+ coef0=coef0+(xx0-1.)*(xx0-0.)*(xx0+1.)*(xx0+2)/24.
+ coef0=coef0-(xx0-2.)*(xx0-1.)*(xx0+0.)*(xx0+1)/24.
+ coef1p=(xx1p-2.)*(xx1p-1.)*(xx1p+0.)*(xx1p+1)/24.
+ coef1p=coef1p-(xx1p-1.)*(xx1p-0.)*(xx1p+1.)*(xx1p+2)/24.
+ coef1p=coef1p-(xx1p-2.)*(xx1p-1.)*(xx1p+0.)*(xx1p+2)/6.
+ coef1p=coef1p+(xx1p-1.)*(xx1p-0.)*(xx1p+1.)*(xx1p+3)/6.
+ coef2p=(xx2p-2.)*(xx2p-1.)*(xx2p+0.)*(xx2p+1)/24.
+ coef2p=coef2p-(xx2p-1.)*(xx2p-0.)*(xx2p+1.)*(xx2p+2)/24.
+ SC(kk,jj,ip1p)=SC(kk,jj,ip1p)+u1m*coef1m+u0*coef0+u1p*coef1p+u2p*coef2p
+
+ coef1m=(xx1m-2.)*(xx1m-1.)*(xx1m+0.)*(xx1m+1)/24.
+ coef0=-(xx0-3.)*(xx0-1.)*(xx0+0.)*(xx0+1)/6.
+ coef1p=(xx1p-2.)*(xx1p-1.)*(xx1p+1.)*(xx1p+2)/4.
+ coef2p=-(xx2p-2.)*(xx2p-1.)*(xx2p+0.)*(xx2p+2)/6.
+ SC(kk,jj,ip2p)=SC(kk,jj,ip2p)+u1m*coef1m+u0*coef0+u1p*coef1p+u2p*coef2p
+
+ coef0=(xx0-2.)*(xx0-1.)*(xx0+0.)*(xx0+1)/24.
+ coef1p=-(xx1p-2.)*(xx1p-0.)*(xx1p+1.)*(xx1p+2)/6.
+ coef2p=(xx2p-2.)*(xx2p-1.)*(xx2p+1.)*(xx2p+2)/4.
+ SC(kk,jj,ip3p)=SC(kk,jj,ip3p)+u0*coef0+u1p*coef1p+u2p*coef2p
+
+ coef1p=(xx1p-1.)*(xx1p-0.)*(xx1p+1.)*(xx1p+2)/24.
+ coef2p=-(xx2p-2.)*(xx2p-0.)*(xx2p+1.)*(xx2p+2)/6.
+ SC(kk,jj,ip4p)=SC(kk,jj,ip4p)+u1p*coef1p+u2p*coef2p
+
+ coef2p=(xx2p-1.)*(xx2p-0.)*(xx2p+1.)*(xx2p+2)/24.
+ SC(kk,jj,ip5p)=SC(kk,jj,ip5p)+u2p*coef2p
+
+ else
+
+ coef1m=(xx1m-2.)*(xx1m-1.)*(xx1m+0.)*(xx1m+1)/24.
+ SC(kk,jj,ip3m)=SC(kk,jj,ip3m)+u1m*coef1m
+
+ coef1m=-(xx1m-2.)*(xx1m-1.)*(xx1m+0.)*(xx1m+2)/6.
+ coef0=(xx0-2.)*(xx0-1.)*xx0*(xx0+1)/24.
+ SC(kk,jj,ip2m)=SC(kk,jj,ip2m)+u1m*coef1m+u0*coef0
+
+ coef1m=(xx1m-2.)*(xx1m-1.)*(xx1m+1.)*(xx1m+2)/4.
+ coef0=-(xx0-2.)*(xx0-1.)*(xx0+0.)*(xx0+2)/6.
+ coef1p=(xx1p-3.)*(xx1p-2.)*(xx1p-1.)*(xx1p+0)/24.
+ SC(kk,jj,ip1m)=SC(kk,jj,ip1m)+u1m*coef1m+u0*coef0+u1p*coef1p
+
+ coef1m=-(xx1m-2.)*(xx1m-0.)*(xx1m+1.)*(xx1m+2)/6.
+ coef1m=coef1m+(xx1m-1.)*(xx1m-0.)*(xx1m+1.)*(xx1m+2)/24.
+ coef0=(xx0-2.)*(xx0-1.)*(xx0+1.)*(xx0+2)/4.
+ coef0=coef0-(xx0-2.)*(xx0-0.)*(xx0+1.)*(xx0+2)/6.
+ coef0=coef0+(xx0-1.)*(xx0-0.)*(xx0+1.)*(xx0+2)/24.
+ coef0=coef0-(xx0-0.)*(xx0+1.)*(xx0+2.)*(xx0+3)/24.
+ coef1p=(xx1p-2.)*(xx1p-1.)*(xx1p+0.)*(xx1p+1)/24.
+ coef1p=coef1p-(xx1p-3.)*(xx1p-2.)*(xx1p-1.)*(xx1p+0)/24.
+ coef1p=coef1p-(xx1p-2.)*(xx1p-1.)*(xx1p+0.)*(xx1p+2)/6.
+ coef1p=coef1p+(xx1p-2.)*(xx1p-1.)*(xx1p+1.)*(xx1p+2)/4.
+ coef2p=(xx2p-2.)*(xx2p-1.)*(xx2p+0.)*(xx2p+1)/24.
+ coef2p=coef2p-(xx2p-2.)*(xx2p-1.)*(xx2p+0.)*(xx2p+2)/6.
+ SC(kk,jj,ip0)=SC(kk,jj,ip0)+u1m*coef1m+u0*coef0+u1p*coef1p+u2p*coef2p
+
+ coef0=(xx0-0.)*(xx0+1.)*(xx0+2.)*(xx0+3)/24.
+ coef1p=-(xx1p-2.)*(xx1p-0.)*(xx1p+1.)*(xx1p+2)/6.
+ coef2p=(xx2p-2.)*(xx2p-1.)*(xx2p+1.)*(xx2p+2)/4.
+ SC(kk,jj,ip1p)=SC(kk,jj,ip1p)+u0*coef0+u1p*coef1p+u2p*coef2p
+
+ coef1p=(xx1p-1.)*(xx1p-0.)*(xx1p+1.)*(xx1p+2)/24.
+ coef2p=-(xx2p-2.)*(xx2p-0.)*(xx2p+1.)*(xx2p+2)/6.
+ SC(kk,jj,ip2p)=SC(kk,jj,ip2p)+u1p*coef1p+u2p*coef2p
+
+ coef2p=(xx2p-1.)*(xx2p-0.)*(xx2p+1.)*(xx2p+2)/24.
+ SC(kk,jj,ip3p)=SC(kk,jj,ip3p)+u2p*coef2p
+
+ endif
+ enddo
+
+end subroutine remeshz_tag
diff --git a/LEGI/advec.F90 b/LEGI/advec.F90
new file mode 100644
index 000000000..863c4020c
--- /dev/null
+++ b/LEGI/advec.F90
@@ -0,0 +1,286 @@
+module advec
+
+ use string
+
+
+ character(len=str_short), private :: type_solveur ! numerical method use to advect the scalar
+
+! public ! XXX En fortran 2003, les public seraient remplacé par des
+! ! protected (sauf pour la procedure advec_init)
+!
+! integer :: nb_part_bl ! number of particles in a block
+! integer :: nb_bound_part_bl ! number of particles in a block boundary wich used a specific remeshing
+! ! formula
+!
+! !integer :: nb_part ! number of particles
+! !integer :: nb_tag_total, nb_tag ! number of tagged particles (the totale one and the one on a given line)
+!
+! ! Variables to switch to "PROTECTED" when we will use fortran2003 norm
+! real(WP) :: xmin, ymin, zmin ! minimal coordinate in the num. domain
+! real(WP) :: xmax, ymax, zmax ! maximal coordinate in the num. domain
+!
+! integer :: Nx_sc,Ny_sc,Nz_sc ! number of mesh for the scalar
+! integer :: Nx_V,Ny_V,Nz_V ! number of mesh for the velocity
+!
+! real(WP) :: dx_sc, dy_sc, dz_sc ! space step for the scalar meshing
+! real(WP) :: dx_V, dy_V, dz_V ! space step for the velocity
+!
+!
+! real(WP) :: sc_min, sc_moy ! particles are create
+! ! only if the scalar reache
+! ! a given threshold
+! real(WP) :: eps_bl ! a small number in comparaison of the block dimensionned lengh
+
+
+ ! Procedure publiques
+ public :: advecX_init
+ character(len=str_short), public :: type_solv
+
+ ! Procedure privées
+ public :: tag_particles_O2
+
+
+ contains
+
+
+! ===== Méthodes publiques et génériques =====
+
+! Initialisation du solveur d'avection (générique, indépendant du solveur utilisé)
+subroutine advec_init(order)
+
+ use data
+ use parallel
+ use advecVAR
+
+ ! @param order = to choose the remeshing method (and thus the order)
+ character(len=str_short), optional, intent(in) :: order
+
+ ! Remplissage de la valeur par defaut en l'absence de choix explicite de la méthode numérique
+ if(present(order)) then
+ type_solveur = order
+ else
+ type_solveur = 'p_O2COR'
+ end if
+
+ ! Common initialisation
+ Nx_sc=nxsc ! nxsc is defined in the "parallel" module (mpi_int.F90)
+ Ny_sc=nysc
+ Nz_sc=nysc
+ Nx_V=nxsc
+ Ny_V=nysc
+ Nz_V=nysc
+
+ xmin = 0.
+ ymin = 0.
+ zmin = 0.
+ xmax = L
+ ymax = L
+ zmax = L
+
+ dx_sc = (xmax-xmin)/Nx_sc ! xmax, xmin = maximal and minimal x-coordinate of the domain
+ dy_sc = (ymax-ymin)/Ny_sc ! ymax, ymin = maximal and minimal y-coordinate on the domain
+ dz_sc = (zmax-zmin)/Nz_sc ! zmax, zmin = maximal and minimal z-coordinate on the domain
+
+ dx_V = (xmax-xmin)/Nx_V ! xmax, xmin = maximal and minimal x-coordinate of the domain
+ dy_V = (ymax-ymin)/Ny_V ! ymax, ymin = maximal and minimal y-coordinate on the domain
+ dz_V = (zmax-zmin)/Nz_V ! zmax, zmin = maximal and minimal z-coordinate on the domain
+
+
+ nb_part_max_bl = 0
+ sc_min = -1 !10**(-5)
+ sc_moy = sc_min/4.
+
+
+
+ ! Initialisation part adapted to each method
+ select case(type_solveur)
+ case('p_O2COR')
+ nb_part_max_bl = 2
+ nb_part_bound_bl = 1
+ case('p_O2LIM')
+ nb_part_max_bl = 2
+ nb_part_bound_bl = 1
+ case('p_O4')
+ nb_part_max_bl = 4
+ nb_part_bound_bl = 2
+ case default
+ nb_part_max_bl = 2
+ nb_part_bound_bl = 1
+ end select
+
+end subroutine advec_init
+
+! Return the name of the particular method used for the advection
+function type_solv()
+ character(len=str_short) :: type_solv
+
+ type_solv = type_solveur
+end function
+
+!> To tag particles which need specific remeshing formula
+subroutine tag_particles_O2(nb_part, N_sc, d_sc, range_min, range_max, p_pos, &
+ & p_V, p_type, nb_tag, nb_no_tag, p_ind_tag, p_ind_no_tag)
+
+ use advecVAR
+
+ ! @param[in] nb_part = number of particles on the current line
+ ! @param[in] N_sc = number of mesh point in the current direction (Nx_sc, Ny_sc or Nz_sc)
+ ! @param[in] d_sc = space step (for the scalar discretisation) in the current direction (dx_sc, dy_sc, dz_sc)
+ ! @param[in] range_min = minimal value of the coordinate in the current direction (xmin, ymin or zmin)
+ ! @param[in] range_max = maximal value of the coordinate in the current direction (xmax, ymax or zmax)
+ ! @param[in] p_pos = particles position(before any advection)
+ ! @param[in] p_V = particles velocity
+ ! @param[out] p_type = particles type for remeshing (left or center)
+ ! @param[out] nb_tag = number of tagged particles on the current line
+ ! @param[out] nb_no_tag = number of untagged particles on the current line
+ ! @param[out] p_ind_tag = indice of tagged particles on the current line
+ ! @param[out] p_ind_no_tag = indice of untagged particles on the current line
+
+ integer, intent(in) :: nb_part, N_sc
+ real(WP), intent(in) :: d_sc, range_min, range_max
+ real(WP), dimension(N_sc), intent(in) :: p_pos
+ real(WP), dimension(nb_part), intent(in):: p_V
+ integer, intent(out) :: nb_tag, nb_no_tag
+ integer,dimension(nb_part),intent(out) :: p_ind_tag, p_ind_no_tag, p_type
+
+ ! information about block
+ integer, dimension(:),allocatable :: bl_nb_part ! nb of particles created on each block
+ integer, dimension(:),allocatable :: bl_final_ind ! indice of the last mesh point of the
+ ! block wich contains a particle
+ real(WP), dimension(:),allocatable :: bl_lambdaMin ! lamda = V*dt/dx,
+ ! bl_lambdaMin =min of lambda in the current block
+ ! the block type depend of its value
+ integer, dimension(:),allocatable :: bl_ind ! integer as lambda in (bl_N,bl_N+1) for a left block
+ ! and lambda in (bl_N-1/2, bl_N+1/2) for a right block
+ ! (called "index" in the article)
+ integer, dimension(:),allocatable :: bl_type ! bloc type
+ integer :: nb_bl ! number of block on the current line
+ real(WP) :: d_bl ! dimensionned length of a block
+ ! Information about particles
+ integer,dimension(nb_part) :: p_ind ! algebric number of mesh crossed over by the particle during its advection
+
+ real(WP) :: cfl ! = d_sc/dt
+ integer :: ind_bl ! number of the current block
+ integer :: ind ! indice used during computation
+ integer :: i, increment ! used in the "while"
+ ! and "do" structure
+
+
+ nb_bl=N_sc/nb_part_max_bl
+ d_bl=float(nb_part_max_bl)*d_sc
+ eps_bl=d_sc/10.0
+ cfl = dt/d_sc
+
+ allocate(bl_nb_part(nb_bl))
+ allocate(bl_final_ind(nb_bl))
+ allocate(bl_lambdaMin(nb_bl))
+ allocate(bl_ind(nb_bl))
+ allocate(bl_type(nb_bl))
+
+ bl_nb_part = 0
+ bl_final_ind = 0
+ bl_lambdaMin = maxval(p_V)*cfl ! or any big value. We can perhaps used
+ ! the cfl condition from the NS solver
+
+ ! Update information on each block
+ do i=1,nb_part
+ ind_bl=1+int((p_pos(i)-range_min+eps_bl)/d_bl) ! XXX le eps_bl
+ !pour éviter les erreurs d'arrondi
+ bl_lambdaMin(ind_bl)=min(bl_lambdaMin(ind_bl),p_V(i)*cfl)
+ bl_nb_part(ind_bl)=bl_nb_part(ind_bl)+1
+ bl_final_ind(ind_bl)=i
+ enddo
+
+ ! bl_lambda_min depends on the velocity of the first particle of the next block too
+ do ind_bl=1,nb_bl-1
+ if (bl_final_ind(ind_bl)/=0) then
+ ind=bl_final_ind(ind_bl)
+ ! the block type depends on the deplacement of the first particle of
+ ! a following block too
+ if ((ind<nb_part).and.(p_pos(ind+1)<p_pos(ind)+1.5*d_sc)) then
+ bl_lambdaMin(ind_bl)=min(bl_lambdaMin(ind_bl),p_V(ind+1)*cfl)
+ end if
+ end if
+ end do
+ ! And at the end of the line
+ ind_bl=nb_bl
+ if (bl_final_ind(ind_bl)/=0) then
+ if ((p_pos(nb_part) >= range_max-1.5*d_sc).and.(p_pos(1) < range_min+0.5*d_sc)) then
+ bl_lambdaMin(ind_bl)=min(bl_lambdaMin(ind_bl),p_V(1)*cfl)
+ end if
+ end if
+
+ ! To determine the block type
+ do ind_bl=1,nb_block
+ ind=nint(bl_lambdaMin(ind_bl))
+ if (bl_lambdaMin(ind_bl)<ind) then
+ ! => center type
+ bl_type(ind_bl)=1
+ bl_ind(ind_bl)=ind
+ else
+ ! => left type
+ bl_type(ind_bl)=0
+ bl_ind(ind_bl)=ind
+ endif
+ end do
+
+ ! Send the information to each particles
+ do i=1,nb_part
+ ind_bl=1+int((p_pos(i)-range_min+eps_bl)/d_bl) ! XXX toujours la même chose.
+ p_type(i)=bl_type(ind_bl)
+ p_ind(i)=bl_ind(ind_bl)
+ end do
+
+ ! Particles tag
+ nb_tag=0
+ nb_no_tag=0
+ ! Far from the boundary
+ i=2
+ do while (i<nb_part)
+ ! NB : if a particle j and j+1 are in the same block they always have
+ ! the same indice p_ind_bl
+ if ((p_ind(i)/=p_ind(i+1)).and.(p_type(i)/=p_type(i+1)).and.(p_pos(i+1)<=p_pos(i)+1.5*d_sc)) then
+ ! The current particle and the following one are tagged
+ nb_tag=nb_tag+1
+ p_ind_tag(nb_tag)=i
+ nb_tag=nb_tag+1
+ p_ind_tag(nb_tag)=i+1
+ increment = 2
+ else
+ nb_no_tag=nb_no_tag+1
+ p_ind_no_tag(nb_no_tag)=i
+ increment=1
+ end if
+ i = i + increment
+ end do
+
+ ! On the boundary
+ if (nb_part>1) then
+ if ((p_ind(1)/=p_ind(nb_part)).and.(p_type(1)/=p_type(nb_part))&
+ & .and.(p_pos(1)+range_max-range_min<=p_pos(nb_part)+1.5*d_sc)) then
+ ! avant on vérifiait que la dernière maille n'était pas déjà taggée :
+ ! "& .and.(itag(ntag).ne.npt)) then" avec itag = p_ind_tag et npt=nb_part
+ nb_tag=nb_tag+1
+ p_ind_tag(nb_tag)=nb_part
+ nb_tag=nb_tag+1
+ p_ind_tag(nb_tag)=1
+ else
+ nb_no_tag = nb_no_tag +1
+ p_ind_no_tag(nb_no_tag)=nb_part
+ nb_no_tag = nb_no_tag +1
+ p_ind_no_tag(nb_no_tag)=1
+ end if
+ end if
+
+ deallocate(bl_nb_part)
+ deallocate(bl_final_ind)
+ deallocate(bl_lambdaMin)
+ deallocate(bl_ind)
+ deallocate(bl_type)
+
+end subroutine tag_particles_O2
+
+!> ===== Private procedure =====
+
+
+end module advec
diff --git a/LEGI/advecVAR.F90 b/LEGI/advecVAR.F90
new file mode 100644
index 000000000..e3f94e3ea
--- /dev/null
+++ b/LEGI/advecVAR.F90
@@ -0,0 +1,36 @@
+module advecVAR
+
+ use precision
+
+! XXX Si passage au fortran 2003 : basculer toutes ces variables dans le module
+! advec (fichier advec.F90) et mettre toutes les variables en protected.
+! Seul la procédure "advec_init" doit pouvoir les modifier, mais de nombreuses
+! procédures doivent pouvoir y accéder.
+
+ public ! XXX En fortran 2003, les public seraient remplacé par des
+ ! protected (sauf pour la procedure advec_init)
+
+ integer :: nb_part_bl ! number of particles in a block
+ integer :: nb_bound_part_bl ! number of particles in a block boundary wich used a specific remeshing
+ ! formula
+
+ !integer :: nb_part ! number of particles
+ !integer :: nb_tag_total, nb_tag ! number of tagged particles (the totale one and the one on a given line)
+
+ ! Variables to switch to "PROTECTED" when we will use fortran2003 norm
+ real(WP) :: xmin, ymin, zmin ! minimal coordinate in the num. domain
+ real(WP) :: xmax, ymax, zmax ! maximal coordinate in the num. domain
+
+ integer :: Nx_sc,Ny_sc,Nz_sc ! number of mesh for the scalar
+ integer :: Nx_V,Ny_V,Nz_V ! number of mesh for the velocity
+
+ real(WP) :: dx_sc, dy_sc, dz_sc ! space step for the scalar meshing
+ real(WP) :: dx_V, dy_V, dz_V ! space step for the velocity
+
+
+ real(WP) :: sc_min, sc_moy ! particles are create
+ ! only if the scalar reache
+ ! a given threshold
+ real(WP) :: eps_bl ! a small number in comparaison of the block dimensionned lengh
+
+end module advecVAR
diff --git a/LEGI/advecX.F90 b/LEGI/advecX.F90
new file mode 100644
index 000000000..f7027b78e
--- /dev/null
+++ b/LEGI/advecX.F90
@@ -0,0 +1,527 @@
+module advecX
+
+ implicit none
+
+ !private
+
+
+ ! Public procedure : generique procedure to advect the scalar with a particles solver
+ public :: advecX_calc
+
+ ! Private porcedures
+ ! particles solver with remeshing method - all of them are corrected to allow large CFL number
+ private :: advecX_calc_O2COR ! order 2
+ !private :: advecX_calc_O2SM ! order 2 with limitaters
+ ! generic procedure used by each method
+ private :: veloX ! calcul de la vitesse
+ private :: countX ! to count the number of particles on a given line
+ ! remeshing method
+ private :: remeshX_O2 ! to remesh untagged particles
+ private :: remeshX_O2_tag ! to remesh tagged particles
+
+
+ contains
+
+! ===== Public procedure =====
+
+
+
+!> Scalar advection (this procedure call the right solver, depending on the
+ !! simulation setup)
+subroutine advecX_calc(dt, Vx, scal)
+
+ use advecVAR
+ use advec
+
+ ! @param[in] dt = time step
+ ! @param[in] Vx = velocity along x (could be discretised on a bigger mesh then the scalar)
+ ! @param[in,out] SC = scalar field to advect
+
+ real(WP), intent(in) :: dt
+ real(WP), dimension(Nx_V, Ny_V, Nz_V), intent(in) :: Vx
+ real(WP), dimension(Nx_sc,Ny_sc,Nz_sc), intent(inout) :: scal
+
+ character(len=str_short) :: type_solveur
+ type_solveur = type_solv()
+
+ select case(type_solveur)
+ case('p_O2COR')
+ call advecX_calc_O2COR(dt, Vx, scal)
+ case default
+ call advecX_calc_O2COR(dt, Vx, scal)
+ end select
+
+end subroutine advecX_calc
+
+
+
+! ===== Private procedure =====
+
+! ---- Order 2 solver, with correction for large time step ----
+
+
+!> Advection during a time step dt - order 2
+subroutine advecX_calc_O2COR(dt, Vx, scal)
+
+ use advecVAR
+ use advec
+
+ ! @param[in] dt = time step
+ ! @param[in] Vx = velocity along x (could be discretised on a bigger mesh then the scalar)
+ ! @param[in,out] scal= scalar field to advect
+
+ real(WP), intent(in) :: dt
+ real(WP), dimension(Nx_V, Ny_V, Nz_V), intent(in) :: Vx
+ real(WP), dimension(Nx_sc,Ny_sc,Nz_sc), intent(inout) :: scal
+
+
+ integer :: i,j,k ! indice of the currend mesh point
+ integer :: n ! indice of the currend particle
+ integer :: nb_tag ! number of tagged particles on the current line
+ integer :: nb_no_tag ! number of untagged particles on the current line
+ integer :: nb_tag_total ! total number of tagged particles
+ integer :: nb_part ! number of particles on the current line
+ integer :: nb_part_total! total number of particles
+
+ real(WP), dimension(Nx_sc) :: p_pos ! particles position
+ real(WP), dimension(Nx_sc) :: p_SC ! scalar field advect by the particles
+ real(WP), dimension(:), allocatable :: p_V ! particles velocity
+ real(WP), dimension(:), allocatable :: p_pos_0 ! initial position of particles before any advection
+ integer, dimension(:), allocatable :: p_type ! does the particle belong to a center block or a left one ?
+ integer, dimension(:), allocatable :: p_ind_tag ! indice of tagged particles
+ integer, dimension(:), allocatable :: p_ind_no_tag ! indice of untagged particles
+
+ nb_part_total = 0
+
+ ! Some particles on the end of bloc need specific remeshing formula. To identify them, we tag them.
+ nb_tag = 0 ! number of tagged particles on the current line !
+ nb_tag_total = 0 ! total number of tagged particles
+
+ do k=1,Nz_sc
+ do j=1,Ny_sc
+ nb_part=0
+ ! Creation of particles
+ call countX(nb_part, j, k, p_pos, p_SC, scal)
+ ! If no particles, no more computations
+ if (nb_part.ne.0) then
+ ! Dynamical allocation of tables
+ allocate(p_V(nb_part))
+ allocate(p_type(nb_part))
+ allocate(p_ind_tag(nb_part))
+ allocate(p_ind_no_tag(nb_part))
+ allocate(p_pos_0(nb_part))
+ ! Advection
+ ! A RK2 scheme is used :
+ ! x(t_n+1)= x(t_n) + dt*v[x(t_n)+dt/2*v(x(t_n))].
+ ! Computation of the intermediar point used in the RK scheme
+ call veloX(nb_part,j,k,p_pos, p_V, Vx)
+ ! XXX la boucle peut être remplacée par des manips directement
+ ! sur les tableaux.
+ do n=1,nb_part
+ ! Save the initial position
+ p_pos_0(n)=p_pos(n)
+ ! Update the position to compute the intermediar point used in
+ ! the RK2 scheme
+ p_pos(n)=p_pos(n)+0.5*dt*p_V(n)
+ ! If the velocity is too high, the particles go out of the box.
+ ! Thus we use the periodical condition.
+ ! We assume that dt*abs(V) < xmax-xmin.
+ if (p_pos(n).gt.xmax) p_pos(n)=p_pos(n)-xmax+xmin
+ if (p_pos(n).lt.xmin) p_pos(n)=p_pos(n)+xmax-xmin
+ end do
+ ! Computation of the "final" velocity
+ call veloX(nb_part,j,k, p_pos, p_V, Vx)
+
+ ! As the velocity is known, we could now determinate :
+ ! - the bloc type (right or left)
+ ! - if the particles on a block boundary need a specific remeshing
+ ! formula (then the particles will be tagged)
+ ! All this is done by the procedure tag_particles.
+ call tag_particles_O2(nb_part, Nx_sc, dx_sc, xmin, xmax, p_pos_0, p_V, &
+ & p_type, nb_tag, nb_no_tag, p_ind_tag, p_ind_no_tag)
+ nb_tag_total = nb_tag_total + nb_tag
+
+ ! End of the advection
+ do n=1, nb_part
+ p_pos(n)=p_pos_0(n)*dt*p_V(n)
+ ! If the velocity is too high, the particles go out of the box.
+ if (p_pos(n).gt.xmax) p_pos(n)=p_pos(n)-xmax+xmin
+ if (p_pos(n).lt.xmin) p_pos(n)=p_pos(n)+xmax-xmin
+ end do
+ ! it is possible to write : p_pos(1:nb_part)=p_pos+dt*p_V
+ ! do not forget that size of p_pos_0, p_V = nb_part but
+ ! size(p_pos)=Nx_sc
+ deallocate(p_pos_0)
+ deallocate(p_V)
+ ! Remeshing the particles
+ ! For not tagged particles :
+ call remeshX_O2(p_pos, p_SC, nb_no_tag, p_ind_no_tag, p_type, j,k,scal)
+ ! For tagged particles :
+ if (nb_tag.ne.0) call remeshX_O2_tag(p_pos, p_SC, nb_tag, p_ind_tag, p_type,j,k,scal)
+
+ nb_part_total = nb_part_total + nb_part
+ end if
+ end do
+ end do
+
+ !print*, 'NPART, NTAG ', npart,ntag_total
+ !print*,'ntag advecX :', ntag_total
+
+ deallocate(p_type)
+ deallocate(p_ind_tag)
+ deallocate(p_ind_no_tag)
+
+end subroutine advecX_calc_O2COR
+
+
+!> To count the number of particles on a line
+subroutine countX(nb_part, j, k, p_pos, p_SC, scal)
+
+ use advecVAR
+
+ ! This method does only depend on the block size and not of the order of the
+ ! method used
+
+ ! @param[out] nb_part = number of particles
+ ! @param[in] j = y-indice of the line
+ ! @param[in] k = z-indice of the line
+ ! @param[out] p_pos = particles position
+ ! @param[out] p_SC = value of scalar field on the particles
+ ! @param[in] scal = scalar field
+
+ integer, intent(out) :: nb_part
+ integer, intent(in) :: j,k
+ real(WP), dimension(Nx_sc), intent(out) :: p_pos, p_SC
+ real(WP), dimension(Nx_sc,Ny_sc,Nz_sc), intent(in) :: scal
+
+ integer :: increment ! increment ensure reading each mesh
+ ! point one time (as we can add more than one particles)
+ integer :: stencil ! number of particles to add when we
+ ! when we are at a block end
+ integer :: i ! x-indice
+ integer :: ind ! temp. int used in some "do"
+ real(WP) :: temp
+
+ ! The following variables are defined in another module
+ ! min_SC = minimal value of the scalar needed to create a particles (if the scalar is null, there's nothing to advect)
+
+ nb_part = 0
+ increment = 1
+ stencil = 2*nb_bound_part_bl
+
+ ! There two cases
+ ! "Interior" particles (ie inside a block) : could be isolated
+ ! Particles on a block boundary : as we use specific remshing formula, they
+ ! could not be isolated to ensure consistance
+ i = nb_bound_part_bl ! table indice = i+1
+ do while(i < Nx_sc+nb_bound_part_bl-1)
+ ! On a block boundary : either no particles either particle on all the
+ ! end of the block and all the beginning of the following one
+ ! end of block <-> mod(i,nb_part_bl) +1 = nb_part_bl - nb_bound_part_bl +1
+ if (mod(i,nb_part_bl)== nb_part_bl - nb_bound_part_bl) then
+ ! as anounced, if the scalar is big enough we create particles at the block end and at
+ ! the next block begining. Therefor, particles are created if the average scalar value
+ ! on these mesh point is big enough. The "sum" is used in the purpose.
+ if (sum((/(abs(scal(mod(i+ind,Nx_sc)+1,j,k)),ind=0,stencil-1)/))>sc_moy*stencil) then
+ do ind = 0, stencil-1
+ nb_part = nb_part +1
+ p_pos(nb_part) = xmin+float(mod(i+ind,Nx_sc))*dx_sc
+ p_SC(nb_part) = scal(mod(i+ind,Nx_sc)+1,j,k)
+ end do
+ end if
+ increment = stencil
+ ! Inside a block
+ ! NB : if we use the minimal block size, there is no "inside" particles
+ else
+ if (abs(scal(i,j,k)).gt.sc_min) then
+ nb_part = nb_part +1
+ p_pos(nb_part) = xmin+float(mod(i+ind,Nx_sc))*dx_sc
+ p_SC(nb_part) = scal(mod(i+ind,Nx_sc)+1,j,k)
+ end if
+ increment = 1
+ end if
+ i = i + increment
+ end do
+
+end subroutine countX
+
+
+!> Compute the velocity on each particles of a line of indice (j_sc, k_sc) on the particles grid
+subroutine veloX(nb_part,j_sc,k_sc, p_pos, p_V, Vx)
+
+ use advecVAR
+
+ ! @param[in] nb_part = number of particles on each line
+ ! @param[in] j = indice corresponding to the line y-coordinate on the scalar grid
+ ! @param[in] k = indice corresponding to the line z-coordinate on the scalar grid
+ ! @param[in] p_pos = particles position
+ ! @param[out] p_V = particles velocity
+ ! @param[in] Vx = velocity along x-axis
+
+
+ integer, intent(in) :: nb_part, j_sc, k_sc ! see above
+ real(WP), dimension(Nx_sc), intent(in) :: p_pos
+ real(WP), dimension(nb_part), intent(out) :: p_V
+ real(WP), dimension(Nx_V,Ny_V,Nz_V), intent(in) :: Vx
+
+
+ integer :: i_Vm, i_Vp ! mesh point on the velocity grid such as:
+ ! x(i_Vm) <= x(i_sc) < x(i_Vm+1=i_Vp)
+ integer :: j_Vm, j_Vp ! as i_Vm, i_Vp but along the y-axis
+ integer :: k_Vm, k_Vp ! as j_Vm, j_Vp but along the z-axis
+ integer :: i
+
+ real(WP) :: y, z ! coordinate of the line
+ real(WP) :: yym, zzm ! adimentioned distance from the previous
+ ! mesh point : yym = [y(j_sc)-y(j_Vm)]/dy_V
+ real(WP) :: yyp, zzp ! adimentioned distance from the next mesh point
+ real(WP) :: am, ap, bm, bp, cm, cp ! weight for the interpolation (a along the x-axis,
+ ! b along the y-axis and c along the z-axis)
+
+ real(WP) :: x, xxm, xxp ! as above but for one particles of the line
+
+
+
+ ! Computation of the (y- and z-) coordinate of the line
+ y=+float(j_sc-1)*dy_sc ! = distance from ymin (could be different from y-coordinate if ymin != 0)
+ z=+float(k_sc-1)*dz_sc ! = distance from zmin
+
+ j_Vm = int(y/dy_V)
+ j_Vp = j_Vm + 1
+ k_Vm = int(z/dz_V)
+ k_Vp = k_Vm + 1
+
+ yym = (y-(j_Vm-1.)*dy_V)/dy_V
+ yyp = 1.0 - yym
+ zzm = (z-(k_Vm-1.)*dz_V)/dz_V
+ zzp = 1.0 - zzm
+
+ ! In fortran table indice start from one. If a particles is out the domain, we use periodic boundaries conditions
+ j_Vm = mod(j_Vm+Ny_V, Ny_V)+1
+ j_Vp = mod(j_Vp+Ny_V, Ny_V)+1
+ k_Vm = mod(k_Vm+Nz_V, Nz_V)+1
+ k_Vp = mod(k_Vp+Nz_V, Nz_V)+1
+
+ bm=yyp
+ bp=yym
+ cm=zzp
+ cp=zzm
+
+ do i = 1,nb_part
+
+ x = p_pos(i)
+
+ i_Vm = int((x-xmin)/dx_V)
+ i_Vp = i_Vm + 1
+
+ xxm = (x-(i_Vm-1.)*dx_V)/dx_V
+ xxp = 1.0 - xxm
+
+ ! In fortran table indice start from one. If a particles is out the domain, we use periodic boundaries conditions
+ i_Vm = mod(i_Vm+Nx_V, Nx_V)+1
+ i_Vp = mod(i_Vp+Nx_V, Nx_V)+1
+
+
+ am = xxp
+ ap = xxm
+
+ ! velocity interpolation
+ p_V(i)= Vx(i_Vm,j_Vm,k_Vm)*am*bm*cm
+ p_V(i)= p_V(i) + Vx(i_Vp,j_Vm,k_Vm)*ap*bm*cm
+ p_V(i)= p_V(i) + Vx(i_Vm,j_Vp,k_Vm)*am*bp*cm
+ p_V(i)= p_V(i) + Vx(i_Vp,j_Vp,k_Vm)*ap*bp*cm
+ p_V(i)= p_V(i) + Vx(i_Vm,j_Vm,k_Vp)*am*bm*cp
+ p_V(i)= p_V(i) + Vx(i_Vp,j_Vm,k_Vp)*ap*bm*cp
+ p_V(i)= p_V(i) + Vx(i_Vm,j_Vp,k_Vp)*am*bp*cp
+ p_V(i)= p_V(i) + Vx(i_Vp,j_Vp,k_Vp)*ap*bp*cp
+
+ end do
+
+end subroutine veloX
+
+!> remeshing with an order 2 method, corrected to allow large CFL number - untagged particles
+subroutine remeshX_O2(p_pos, p_sc, nb_no_tag, p_ind_no_tag, p_type,j,k, scal)
+
+ use advecVAR
+
+ ! @param[in] p_pos = particles position
+ ! @param[in] p_sc = scalar advected by the particles
+ ! @param[in] nb_no_tag = number of untagged particles
+ ! @param[in] p_ind_tag = indice of untagged particles
+ ! @param[in] j,k = indice of of the current line (y-coordinate and z-coordinate)
+ ! @param[in] p_type = equal 0 (resp 1) if the particle belong to a centered (resp left) block
+ ! @param[in,out] scal = scalar field
+
+ integer, intent(in) :: nb_no_tag, j, k
+ integer, dimension(:), intent(in) :: p_ind_no_tag, p_type
+ real(WP), dimension(:), intent(in) :: p_pos, p_sc
+ real(WP), dimension(Nx_sc,Ny_sc,Nz_sc), intent(inout) :: scal
+
+ integer :: i, ind ! indice of the current untagged particles and
+ ! number of the current particles in all
+ ! the (untagged and tagged) particles
+ integer :: iM, i0, iP ! indice of the the mesh point use in the remeshing
+ ! (they depend on the block type)
+ real(WP) :: x0, xP, xM, a0, aP, aM ! interpolation weight
+ real(WP) :: pos, sca ! position of the particle and scalar value which is advected
+
+
+
+ scal(:,j,k)=0.
+
+
+ do i = 1,nb_no_tag
+ ind = p_ind_no_tag(i)
+ sca = p_sc(ind)
+ pos = p_pos(ind)
+
+ ! Mesh point used in remeshing formula
+ if (p_type(ind)==0) then
+ ! the particle belong to a left bloc
+ i0 = int((pos-xmin)/dx_sc)
+ else
+ ! the particle belong to a centered bloc
+ i0 = nint((pos-xmin)/dx_sc)
+ end if
+ iM = i0 - 1
+ iP = i0 + 1
+
+ ! Distance to mesh points
+ x0 = (pos - float(i0)*dx_sc-xmin)/dx_sc
+ xP=x0+1
+ xM=x0-1
+
+ ! Table indice strats from 1 + periodic boundary conditions
+ iM=mod(iM+Nx_sc,Nx_sc) +1
+ i0=mod(i0+Nx_sc,Nx_sc) +1
+ iP=mod(iP+Nx_sc,Nx_sc) +1
+
+ ! Interpolation weights
+ aM=0.5*x0*xM
+ a0=1.-x0**2
+ aP=0.5*x0*xP
+
+ ! Remeshing
+ scal(iM,j,k) = scal(iM,j,k) + sca*aM
+ scal(i0,j,k) = scal(i0,j,k) + sca*a0
+ scal(iP,j,k) = scal(iP,j,k) + sca*aP
+
+ end do
+
+end subroutine remeshX_O2
+
+
+!> remeshing with an order 2 method, corrected to allow large CFL number - tagged particles
+subroutine remeshX_O2_tag(p_pos, p_sc, nb_tag, p_ind_tag, p_type,j,k,scal)
+
+ use advecVAR
+
+ ! @param[in] p_pos = particles position
+ ! @param[in] p_sc = scalar advected by the particles
+ ! @param[in] nb_no_tag = number of tagged particles
+ ! @param[in] p_ind_tag = indice of tagged particles
+ ! @param[in] j,k = indice of of the current line (y-coordinate and z-coordinate)
+ ! @param[in] p_type = equal 0 (resp 1) if the particle belong to a centered (resp left) block
+ ! @param[in,out] scal = scalar field
+
+ integer, intent(in) :: nb_tag, j, k
+ integer, dimension(:), intent(in) :: p_ind_tag, p_type
+ real(WP), dimension(:), intent(in) :: p_pos, p_sc
+ real(WP), dimension(Nx_sc,Ny_sc,Nz_sc), intent(inout) :: scal
+
+ integer :: i, ind, indP ! indice of the current particles in the group of tagged
+ ! particles and in the group of all particles.
+ ! ind1 is the indice of the next (tagged) particles
+ real(WP) :: pos, sca ! position of the particle and scalar value which is advected
+ real(WP) :: posP, scaP ! position of the next particle and scalar value which is advected
+ integer :: iM, i0, iP, iP2, iP3 ! indice of the the nearest mesh points
+ ! (they depend on the block type)
+ integer :: i0_bis ! indice of the the nearest mesh point for the indP=ind+1 particle
+ real(WP) :: x0, xM, xP, aM, a0, aP, aP2 ! interpolation weight for the particles
+ real(WP) :: xM_bis, x0_bis, xP_bis ! distance to mesh points for the next particle (ind+1=indP)
+ real(WP) :: b0, bP, bP2, bP3 ! interpolation weight for the next particle (ind+1=indP)
+
+ ! The particles are tagged by group of 2 : the one at a block end and the
+ ! following on the beginning of the next block
+ do i=1,nb_tag-1,2
+ ind=p_ind_tag(i)
+ pos = p_pos(ind)
+ sca = p_sc(ind)
+ indP=p_ind_tag(i+1)
+ posP= p_pos(indP)
+ scaP = p_sc(indP)
+ if (p_type(ind)==0) then
+ ! the current block is a left one. As it is tagged, the following
+ ! one is a centered one
+ ! Indice of mesh point used to remesh
+ i0 = int((pos-xmin)/dx_sc)
+ i0_bis = nint((posP-xmin)/dx_sc)
+ iM=i0-1
+ iP=i0+1
+ iP2=i0+2
+ iP3=i0+3
+ iM=mod(iM+Nx_sc,Nx_sc) +1
+ i0=mod(i0+Nx_sc,Nx_sc) +1
+ iP=mod(iP+Nx_sc,Nx_sc) +1
+ iP2=mod(iP2+Nx_sc,Nx_sc) +1
+ iP3=mod(iP3+Nx_sc,Nx_sc) +1
+
+ ! Distance to mesh point
+ x0 = (pos - float(i0)*dx_sc-xmin)/dx_sc
+ x0_bis = (posP - float(i0_bis)*dx_sc-xmin)/dx_sc
+ xM=x0-1
+ xP=x0+1
+ xP_bis=x0_bis+1
+ xM_bis=x0_bis-1
+
+ ! Interpolation weights
+ aM=0.5*x0*xM
+ a0=1.-x0**2
+ aP=x0
+ aP2=0.5*xM*x0
+ b0=xP_bis*x0_bis/2.
+ bP=-x0_bis
+ bP2=1.-x0_bis**2
+ bP3=x0_bis*xP_bis/2.
+
+ ! Remeshing
+ scal(iM,j,k)= scal(iM,j,k)+aM*sca
+ scal(i0,j,k)= scal(i0,j,k)+a0*sca+b0*scaP
+ scal(iP,j,k)= scal(iP,j,k)+aP*sca+bP*scaP
+ scal(iP2,j,k)=scal(iP2,j,k)+aP2*sca+bP2*scaP
+ scal(iP3,j,k)=scal(iP3,j,k)+bP3*scaP
+
+ ! if the particules belong to a left block (and therefore the following
+ ! particles belong to a centered one - verwise it wouldn't be tagged)
+ else
+ ! Compute the weight
+ i0 = nint((pos-xmin)/dx_sc)
+ i0_bis = int((posP-xmin)/dx_sc)
+ iM=i0-1
+ iP=i0+1
+ iM=mod(iM+Nx_sc,Nx_sc) +1
+ i0=mod(i0+Nx_sc,Nx_sc) +1
+ iP=mod(iP+Nx_sc,Nx_sc) +1
+
+ x0 = (pos - float(i0)*dx_sc-xmin)/dx_sc
+ x0_bis = (posP - float(i0_bis)*dx_sc-xmin)/dx_sc
+ xM = x0-1
+ xP_bis = x0_bis+1
+
+
+ aM=0.5*x0*xP
+ a0=1-aM
+ bP=0.5*x0_bis*xP_bis
+ b0=1-bP
+
+ ! Update the scalar field
+ scal(iM,j,k)=scal(iM,j,k)+aM*sca
+ scal(i0,j,k)=scal(i0,j,k)+a0*sca+b0*scaP
+ scal(iP,j,k)=scal(iP,j,k)+bP*scaP
+ endif
+ enddo
+
+end subroutine remeshX_O2_tag
+
+
+end module advecX
diff --git a/LEGI/advecY.F90 b/LEGI/advecY.F90
new file mode 100644
index 000000000..65a9f0cb6
--- /dev/null
+++ b/LEGI/advecY.F90
@@ -0,0 +1,529 @@
+module advecY
+
+ ! XXX Si on déplace ce use dans les subroutines qui en ont besoin, on est
+ ! pas gagnant ? (les variables publiques de advec ne seraient plus visibles
+ ! pas tous !?!)
+ use advec
+
+ implicit none
+
+
+ ! Public procedure : generique procedure to advect the scalar with a particles solver
+ public :: advecY_calc
+
+ ! Private porcedures
+ ! particles solver with remeshing method - all of them are corrected to allow large CFL number
+ private :: advecY_calc_O2COR ! order 2
+ !private :: advecX_calc_O2SM ! order 2 with limitaters
+ ! generic procedure used by each method
+ private :: veloY ! calcul de la vitesse
+ private :: countY ! to count the number of particles on a given line
+ ! remeshing method
+ private :: remeshY_O2
+ private :: remeshY_O2_tag
+
+
+ contains
+
+! ===== Public procedure =====
+
+
+
+!> Scalar advection (this procedure call the right solver, depending on the
+ !! simulaiton setup)
+subroutine advecY_calc(dt, Vy, SC)
+
+ use advecVAR
+
+ ! @param dt = time step
+ ! @param Vy = velocity along y (could be discretised on a bigger mesh then the scalar)
+ ! @param[in,out] SC = scalar field to advect
+
+ real(WP), intent(in) :: dt
+ real(WP), dimension(Nx_V, Ny_V, Nz_V), intent(in) :: Vy
+ real(WP), dimension(Nx_sc,Ny_sc,Nz_sc), intent(inout) :: SC
+
+ character(len=str_short) :: type_solveur
+ type_solveur = type_solv()
+
+ select case(type_solveur)
+ case('p_O2COR')
+ call advecY_calc_O2COR(dt, Vy, SC)
+ case default
+ call advecY_calc_O2COR(dt, Vy, SC)
+ end select
+
+end subroutine advecY_calc
+
+! ===== Private procedure =====
+
+! ---- Order 2 solver, with correction for large time step ----
+
+
+!> Advection during a time step dt - order 2
+subroutine advecY_calc_O2COR(dt,Vy,SC)
+
+ use advecVAR
+
+ ! @param dt = time step
+ ! @param Vy = velocity along y (could be discretised on a bigger mesh then the scalar)
+ ! @param[in,out] SC = scalar field to advect
+
+ real(WP), intent(in) :: dt
+ real(WP), dimension(Nx_V, Ny_V, Nz_V), intent(in) :: Vy
+ real(WP), dimension(Nx_sc,Ny_sc,Nz_sc), intent(inout) :: SC
+
+ integer :: i,j,k ! indice of the currend mesh point
+ integer :: n ! indice of the currend particle
+ integer :: jr, kr, np_aux, ini
+ integer :: nb_tag ! number of tagged particles on the current line
+ integer :: nb_no_tag ! number of untagged particles on the current line
+ integer :: nb_tag_total ! total number of tagged
+ integer :: nb_part ! number of particles on the current line
+ integer :: nb_part_total! total number of particles
+
+ real(WP), dimension(Ny_sc) :: p_pos ! particles position
+ real(WP), dimension(Ny_sc) :: p_SC ! scalar field advect by the particles
+ real(WP), dimension(:), allocatable :: p_V ! particles velocity
+ real(WP), dimension(:), allocatable :: p_pos_0 ! initial position of particles before any advection
+ integer, dimension(:), allocatable :: p_type ! does the particle belong to a center block or a left one ?
+ integer, dimension(:), allocatable :: p_ind_tag ! indice of tagged particles
+ integer, dimension(:), allocatable :: p_ind_no_tag ! indice of untagged particles
+
+ nb_part_total = 0
+
+ ! Some particles on block boundaries need specific remeshing formula.
+ ! Therefor they are tagged
+ nb_tag = 0
+ nb_tag_total = 0
+
+ do k=1,Nz_sc
+ do i=1,Nx_sc
+ nb_part=0
+ ! Creation of particles
+ call countY(nb_part, i, k, p_pos, p_SC, SC)
+ ! If no particles, no more computations
+ if (nb_part.ne.0) then
+ ! Dynamical allocation of tables
+ allocate(p_V(nb_part))
+ allocate(p_type(nb_part))
+ allocate(p_ind_tag(nb_part))
+ allocate(p_ind_no_tag(nb_part))
+ allocate(p_pos_0(nb_part))
+ ! Advection
+ ! A RK2 scheme is used :
+ ! y(t_n+1)= y(t_n) + dt*v[y(t_n)+dt/2*v(y(t_n))].
+ ! Computation of the intermediar point used in the RK scheme
+ call veloY(nb_part,i,k,p_pos, p_V, Vy)
+ ! XXX la boucle peut être remplacée par des manips directement
+ ! sur les tableaux.
+ do n=1,nb_part
+ ! Save the initial position
+ p_pos_0(n)=p_pos(n)
+ ! Update the position to compute the intermediar point used in
+ ! the RK2 scheme
+ p_pos(n)=p_pos(n)+0.5*dt*p_V(n)
+ ! If the velocity is too high, the particles go out of the box.
+ ! Thus we use the periodical condition.
+ ! We assume that dt*abs(V) < ymax-ymin.
+ if (p_pos(n).gt.ymax) p_pos(n)=p_pos(n)-ymax+ymin
+ if (p_pos(n).lt.ymin) p_pos(n)=p_pos(n)+ymax-ymin
+ end do
+ ! Computation of the "final" velocity
+ call veloY(nb_part,i,k, p_pos, p_SC, Vy)
+
+ ! As the velocity is known, we could now determinate :
+ ! - the bloc type (right or left)
+ ! - if the particles on a block boundary need a specific remeshing
+ ! formula (then the particles will be tagged)
+ ! All this is done by the procedure tag_particles.
+ call tag_particles_O2(nb_part, Ny_sc, dy_sc, ymin, ymax, p_pos_0, p_V,&
+ & p_type, nb_tag, nb_no_tag, p_ind_tag, p_ind_no_tag)
+ nb_tag_total = nb_tag_total + nb_tag
+
+ ! End of the advection
+ do n=1, nb_part
+ p_pos(n)=p_pos_0(n)*dt*p_V(n)
+ ! If the velocity is too high, the particles go out of the box.
+ if (p_pos(n).gt.ymax) p_pos(n)=p_pos(n)-ymax+ymin
+ if (p_pos(n).lt.ymin) p_pos(n)=p_pos(n)+ymax-ymin
+ end do
+ ! it is possible to write : p_pos(1:nb_part)=p_pos+dt*p_V
+ ! do not forget that size of p_pos_0, p_V = nb_part but
+ ! size(p_pos)=Ny_sc
+ deallocate(p_pos_0)
+ deallocate(p_V)
+ ! Remeshing the particles
+ ! For not tagged particles :
+ call remeshY_O2(p_pos, p_sc, nb_no_tag, p_ind_no_tag, p_type, i,k, SC)
+ ! For tagged particles :
+ if (nb_tag.ne.0) call remeshY_O2_tag(p_pos, p_SC, nb_tag, p_ind_tag, p_type,j,k, SC)
+
+ nb_part_total = nb_part_total + nb_part
+ end if
+ end do
+ end do
+
+ !print*, 'NPART, NTAG ', npart,ntag_total
+ !print*,'ntag advecY :', ntag_total
+
+
+ deallocate(p_type)
+ deallocate(p_ind_tag)
+ deallocate(p_ind_no_tag)
+
+
+end subroutine advecY_calc_O2COR
+
+!> To count the number of particles on a line (along the y-axis)
+subroutine countY(nb_part, i, k, p_pos, p_SC, SC)
+
+ use advecVAR
+
+ ! This method does only depend on the block size and not of the order of the
+ ! method used
+
+ ! @param[out] nb_part = number of particles
+ ! @param[in] i = x-indice of the line
+ ! @param[in] k = z-indice of the line
+ ! @param[out] p_pos = particles position
+ ! @param[out] p_SC = value of scalar field on the particles
+ ! @param[in] SC = scalar field to advect
+
+ integer, intent(out) :: nb_part
+ integer, intent(in) :: i,k
+ real(WP), dimension(Ny_sc), intent(out) :: p_pos, p_SC
+ real(WP), dimension(Nx_sc,Ny_sc,Nz_sc), intent(in) :: SC
+
+ integer :: increment ! increment ensure reading each mesh
+ ! point one time (as we can add more than one particles)
+ integer :: stencil ! number of particles to add when we
+ ! when we are at a block end
+ integer :: j ! y-indice
+ integer :: ind ! temp. int used in some "do"
+
+ ! The following variables are defined in another module
+ ! SC = the scalar field
+ ! min_SC = minimal value of the scalar needed to create a particles (if the scalar is null, there's nothing to advect)
+
+ nb_part = 0
+ p_pos = 0
+ p_SC = 0
+
+ increment = 1
+ stencil = 2*nb_bound_part_bl
+
+ ! There two cases
+ ! "Interior" particles (ie inside a block) : could be isolated
+ ! Particles on a block boundary : as we use specific remshing formula, they could
+ ! not be isolated to ensure consistance. If a block end contains some particle
+ ! the begining of the following block will contain particles. It explains why
+ ! the counter do not start a 0. The first particles are initialized with the one
+ ! at the line end
+ j = nb_bound_part_bl
+ do while(j < Ny_sc+nb_bound_part_bl-1)
+ ! On a block boundary : either no particles either particle on all the
+ ! end of the block and all the beginning of the following one
+ ! end of block <-> mod(i,nb_part_bl) +1 = nb_part_bl - nb_bound_part_bl +1
+ if (mod(j,nb_part_bl)== nb_part_bl - nb_bound_part_bl) then
+ ! as anounced, if the scalar is big enough we create particles at the block end and at
+ ! the next block begining. Therefor, particles are created if the average scalar value
+ ! on these mesh point is big enough. The "sum" is used in the purpose.
+ if (sum((/(abs(SC(i,mod(j+ind,Ny_sc)+1,k)),ind=0,stencil-1)/))>sc_moy*stencil) then
+ do ind = 0, stencil-1
+ nb_part = nb_part +1
+ p_pos(nb_part) = ymin+float(mod(j+ind,Ny_sc))*dy_sc
+ p_SC(nb_part) = SC(i,mod(j+ind,Ny_sc)+1,k)
+ end do
+ end if
+ increment = stencil
+ ! Inside a block
+ ! NB : if we use the minimal block size, there is no "inside" particles
+ else
+ if (abs(SC(i,j,k)).gt.sc_min) then
+ nb_part = nb_part +1
+ p_pos(nb_part) = ymin+float(mod(j+ind,Ny_sc))*dy_sc
+ p_SC(nb_part) = SC(i,mod(j+ind,Ny_sc)+1,k)
+ end if
+ increment = 1
+ end if
+ j = j + increment
+ end do
+end subroutine countY
+
+
+!> Compute the velocity on each particles of a line of indice (i_sc, k_sc) on the particles grid
+subroutine veloY(nb_part,i_sc,k_sc, p_pos, p_V, Vy)
+
+ use advecVAR
+
+ ! @param[in] nb_part = number of particles on each line
+ ! @param[in] i = indice corresponding to the line y-coordinate on the scalar grid
+ ! @param[in] k = indice corresponding to the line y-coordinate on the scalar grid
+ ! @param[in] p_pos = particles position
+ ! @param[out] p_V = particles velocity
+ ! @param[in] Vy = velocity field (component along y-axis)
+
+ integer, intent(in) :: nb_part, i_sc, k_sc ! see above
+ real(WP), dimension(Ny_sc), intent(in) :: p_pos
+ real(WP), dimension(nb_part), intent(out) :: p_V
+ real(WP), dimension(Nx_V, Ny_V, Nz_V), intent(in) :: Vy
+
+
+ integer :: i_Vm, i_Vp ! mesh point on the velocity grid such as:
+ ! x(i_Vm) <= x(i_sc) < x(i_Vm+1=i_Vp)
+ integer :: j_Vm, j_Vp ! as i_Vm, i_Vp but along the y-axis
+ integer :: k_Vm, k_Vp ! as j_Vm, j_Vp but along the z-axis
+ integer :: j ! indice of the current particle of the line
+
+ real(WP) :: x, z ! coordinate of the line
+ real(WP) :: xxm, zzm ! adimentioned distance from the previous
+ ! mesh point : yym = [y(j_sc)-y(j_Vm)]/dx_V
+ real(WP) :: xxp, zzp ! adimentioned distance from the next mesh point
+ real(WP) :: am, ap, bm, bp, cm, cp ! weight for the interpolation (a along the x-axis,
+ ! b along the y-axis and c along the z-axis)
+
+ real(WP) :: y, yym, yyp ! as above but for one particles of the line
+
+
+
+ ! Computation of the (y- and z-) coordinate of the line
+ x=+float(i_sc-1)*dx_sc ! = distance from ymin (could be different from y-coordinate)
+ z=+float(k_sc-1)*dz_sc ! = distance from zmin
+
+ i_Vm = int(x/dx_V)
+ i_Vp = i_Vm + 1
+ k_Vm = int(z/dz_V)
+ k_Vp = k_Vm + 1
+
+ xxm = (x-(i_Vm-1.)*dx_V)/dx_V
+ xxp = 1.0 - xxm
+ zzm = (z-(k_Vm-1.)*dz_V)/dz_V
+ zzp = 1.0 - zzm
+
+ ! In fortran table indice start from one. If a particles is out the domain, we use periodic boundaries conditions
+ i_Vm = mod(i_Vm+Nx_V, Nx_V)+1
+ i_Vp = mod(i_Vp+Nx_V, Nx_V)+1
+ k_Vm = mod(k_Vm+Nz_V, Nz_V)+1
+ k_Vp = mod(k_Vp+Nz_V, Nz_V)+1
+
+ am=xxp
+ ap=xxm
+ cm=zzp
+ cp=zzm
+
+ do j = 1,nb_part
+
+ y = p_pos(j)
+
+ j_Vm = int((y-ymin)/dy_V)
+ j_Vp = j_Vm + 1
+
+ yym = (y-(j_Vm-1.)*dy_V)/dy_V
+ yyp = 1.0 - yym
+
+ ! In fortran table indice start from one. If a particles is out the domain, we use periodic boundaries conditions
+ j_Vm = mod(j_Vm+Ny_V, Ny_V)+1
+ j_Vp = mod(j_Vp+Ny_V, Ny_V)+1
+
+
+ am = xxp
+ ap = xxm
+
+ ! velocity interpolation
+ p_V(j)= Vy(i_Vm,j_Vm,k_Vm)*am*bm*cm
+ p_V(j)= p_V(j) + Vy(i_Vp,j_Vm,k_Vm)*ap*bm*cm
+ p_V(j)= p_V(j) + Vy(i_Vm,j_Vp,k_Vm)*am*bp*cm
+ p_V(j)= p_V(j) + Vy(i_Vp,j_Vp,k_Vm)*ap*bp*cm
+ p_V(j)= p_V(j) + Vy(i_Vm,j_Vm,k_Vp)*am*bm*cp
+ p_V(j)= p_V(j) + Vy(i_Vp,j_Vm,k_Vp)*ap*bm*cp
+ p_V(j)= p_V(j) + Vy(i_Vm,j_Vp,k_Vp)*am*bp*cp
+ p_V(j)= p_V(j) + Vy(i_Vp,j_Vp,k_Vp)*ap*bp*cp
+
+ end do
+
+end subroutine veloY
+
+!> remeshing with an order 2 method, corrected to allow large CFL number - untagged particles
+subroutine remeshY_O2(p_pos, p_sc, nb_no_tag, p_ind_no_tag, p_type,i,k,SC)
+
+ use advecVAR
+
+ ! @param[in] p_pos = particles position
+ ! @param[in] p_sc = scalar advected by the particles
+ ! @param[in] nb_no_tag = number of untagged particles
+ ! @param[in] p_ind_tag = indice of untagged particles
+ ! @param[in] i,k = indice of of the current line (x-coordinate and z-coordinate)
+ ! @param[in] p_type = equal 0 (resp 1) if the particle belong to a centered (resp left) block
+ ! @param[in,out] SC = scalar field to advect
+
+ integer, intent(in) :: nb_no_tag, i, k
+ integer, dimension(:), intent(in) :: p_ind_no_tag, p_type
+ real(WP), dimension(:), intent(in) :: p_pos, p_sc
+ real(WP), dimension(Nx_sc,Ny_sc,Nz_sc), intent(inout) :: SC
+
+ integer :: j, ind ! indice of the current untagged particles and
+ ! number of the current particles in all
+ ! the (untagged and tagged) particles
+ integer :: jM, j0, jP ! indice of the the mesh point use in the remeshing
+ ! (they depend on the block type)
+ real(WP) :: yM, y0, yP ! adimensionned distance to mesh points
+ real(WP) :: bM, b0, bP ! interpolation weight
+ real(WP) :: pos, sca ! position of the particle and scalar it advects
+
+
+ SC(i,:,k)=0.
+
+ do j = 1, nb_no_tag
+ ind = p_ind_no_tag(j)
+ sca = p_sc(ind)
+ pos = p_pos(ind)
+
+ ! Mesh point used in remeshing formula
+ if (p_type(ind)==0) then
+ ! the particle belong to a centered bloc
+ j0 = int((pos-ymin)/dy_sc)
+ else
+ ! the particle belong to a left block
+ j0 = nint((pos-ymin)/dy_sc)
+ end if
+ jM = j0 - 1
+ jP = j0 + 1
+
+ ! Distance to mesh points
+ y0 = (pos - float(j0)*dy_sc-ymin)/dy_sc
+ yM = y0-1
+ yP = y0+1
+
+ jM=mod(jM+Ny_sc,Ny_sc) +1
+ j0=mod(j0+Ny_sc,Ny_sc) +1
+ jP=mod(jP+Ny_sc,Ny_sc) +1
+
+ ! Interpolation weights
+ bM=0.5*y0*yM
+ b0=1.-y0**2
+ bP=0.5*y0*yP
+
+ ! remeshing
+ SC(i,jM,k) = SC(i,jM,k) + sca*bM
+ SC(i,j0,k) = SC(i,j0,k) + sca*b0
+ SC(i,jP,k) = SC(i,jP,k) + sca*bP
+
+ end do
+
+end subroutine remeshY_O2
+
+!> remeshing with an order 2 method, corrected to allow large CFL number - tagged particles
+subroutine remeshY_O2_tag(p_pos, p_sc, nb_tag, p_ind_tag, p_type,i,k,SC)
+
+ use advecVAR
+
+ ! @param[in] p_pos = particles position
+ ! @param[in] p_sc = scalar advected by the particles
+ ! @param[in] nb_no_tag = number of tagged particles
+ ! @param[in] p_ind_tag = indice of tagged particles
+ ! @param[in] i,k = indice of of the current line (x-coordinate and y-coordinate)
+ ! @param[in] ptype = equal 0 (resp 1) if the particle belong to a centered (resp left) block
+ ! @param[in,out] SC = scalar field to advect
+
+ integer, intent(in) :: nb_tag, i, k
+ integer, dimension(:), intent(in) :: p_ind_tag, p_type
+ real(WP), dimension(:), intent(in) :: p_pos, p_sc
+ real(WP), dimension(Nx_sc,Ny_sc,Nz_sc), intent(inout) :: SC
+
+ integer :: j, ind, indP ! indice of the current particles in the group of tagged
+ ! particles and in the group of all particles.
+ real(WP) :: pos, sca ! position of the particle and scalar it advects
+ real(WP) :: posP, scaP ! position of the next particle and scalar it advects
+ ! indP is the indice of the next (tagged) particles
+ integer :: jM, j0, jP, jP2, jP3 ! indice of the the nearest mesh points
+ ! (they depend on the block type)
+ integer :: j0_bis ! indice of the the nearest mesh point for the indP=ind+1 particle
+ real(WP) :: yM, y0, yP, aM, a0, aP,aP2 ! interpolation weight for the particles
+ real(WP) :: yM_bis, y0_bis, yP_bis ! adimensionned distance to mesh points
+ real(wp) :: b0, bP, bP2, bP3 ! interpolation weight for the next particle (ind+1=indP)
+
+ ! The particles are tagged by group of 2 : the one at a block end and the
+ ! following on the beginning of the next block
+ do j=1,nb_tag-1,2
+ ind=p_ind_tag(j)
+ pos = p_pos(ind)
+ sca = p_sc(ind)
+ indP=p_ind_tag(j+1)
+ posP= p_pos(indP)
+ scaP = p_sc(indP)
+
+ if (p_type(ind)==0) then
+ ! the current block is a left one. As it is tagged, the following
+ ! one is a centered one
+
+ ! Indice of mesh point used in order to remesh
+ j0 = int((pos-ymin)/dy_sc)
+ j0_bis = nint((posP-ymin)/dy_sc)
+ jM=j0-1
+ jP=j0+1
+ jP2=j0+2
+ jP3=j0+3
+ jM=mod(jM+Ny_sc,Ny_sc) +1
+ j0=mod(j0+Ny_sc,Ny_sc) +1
+ jP=mod(jP+Ny_sc,Ny_sc) +1
+ jP2=mod(jP2+Ny_sc,Ny_sc) +1
+ jP3=mod(jP3+Ny_sc,Ny_sc) +1
+
+ ! Distance to mesh point
+ y0 = (pos - float(j0)*dy_sc-ymin)/dy_sc
+ y0_bis = (posP - float(j0_bis)*dy_sc-ymin)/dy_sc
+ yP=y0+1
+ yM=y0-1
+ yP_bis=y0_bis+1
+ yM_bis=y0_bis-1
+
+ ! Interpolation weight
+ aM=y0*yM/2.
+ a0=1-y0**2
+ aP=y0
+ aP2=y0*yM/2.
+ b0=y0_bis*yP_bis/2.
+ bP=-y0_bis
+ bP2=1-y0_bis**2
+ bP3=y0_bis*yP_bis/2.
+
+ ! Remeshing
+ SC(i,jM,k)=SC(i,jM,k)+aM*sca
+ SC(i,j0,k)=SC(i,j0,k)+a0*sca+b0*scaP
+ SC(i,jP,k)=SC(i,jP,k)+aP*sca+bP*scaP
+ SC(i,jP2,k)=SC(i,jP2,k)+aP2*sca+bP2*scaP
+ SC(i,jP3,k)=SC(i,jP3,k)+bP3*scaP
+
+ else
+ ! the particules belong to a centered block (and therefore the following
+ ! particles belong to a left one - verwise it wouldn't be tagged)
+
+ j0 = nint((pos-ymin)/dy_sc)
+ j0_bis = int((posP-ymin)/dy_sc)
+ jM=j0-1
+ jP=j0+1
+ jM=mod(jM+Ny_sc,Ny_sc) +1
+ j0=mod(j0+Ny_sc,Ny_sc) +1
+ jP=mod(jP+Ny_sc,Ny_sc) +1
+
+ y0 = (pos - float(j0)*dy_sc-ymin)/dy_sc
+ y0_bis = (posP - float(j0_bis)*dy_sc-ymin)/dy_sc
+ yM=y0-1.
+ yP_bis=y0_bis+1
+
+ aM=0.5*y0*yP
+ a0=1-aM
+ bP=0.5*y0_bis*yP_bis
+ b0=1-bP
+
+ ! Remeshing
+ SC(i,jM,k)=SC(i,jM,k)+aM*sca
+ SC(i,j0,k)=SC(i,j0,k)+a0*sca+b0*scaP
+ SC(i,jP,k)=SC(i,jP,k)+bP*scaP
+ endif
+ enddo
+
+end subroutine remeshY_O2_tag
+
+end module advecY
diff --git a/LEGI/advecZ.F90 b/LEGI/advecZ.F90
new file mode 100644
index 000000000..c2ad4f90b
--- /dev/null
+++ b/LEGI/advecZ.F90
@@ -0,0 +1,514 @@
+module advecZ
+
+ use advec
+
+ implicit none
+
+ !private
+
+
+ ! Public procedure : generique procedure to advect the scalar with a particles solver
+ public :: advecZ_calc
+
+ ! Private porcedures
+ ! particles solver with remeshing method - all of them are corrected to allow large CFL number
+ private :: advecZ_calc_O2COR ! order 2
+ !private :: advecZ_calc_O2SM ! order 2 with limitaters
+ ! generic procedure used by each method
+ private :: veloZ ! calcul de la vitesse
+ private :: countZ ! to count the number of particles on a given line
+ ! remeshing method
+ private :: remeshZ_O2
+ private :: remeshZ_O2_tag
+
+
+ contains
+
+! ===== Public procedure =====
+
+
+
+!> Scalar advection (this procedure call the right solver, depending on the
+ ! simulation setup)
+subroutine advecZ_calc(dt, Vz, SC)
+
+ use advecVAR
+
+ ! @param[in] dt = time step
+ ! @param[in] Vx = velocity along x (could be discretised on a bigger mesh then the scalar)
+ ! @param[in,out] SC = scalar field to advect
+
+ real(WP), intent(in) :: dt
+ real(WP), dimension(Nx_V, Ny_V, Nz_V), intent(in) :: Vz
+ real(WP), dimension(Nx_sc,Ny_sc,Nz_sc), intent(inout) :: SC
+
+
+ character(len=str_short) :: type_solveur
+ type_solveur = type_solv()
+
+ select case(type_solveur)
+ case('p_O2COR')
+ call advecZ_calc_O2COR(dt, Vz, SC)
+ case default
+ call advecZ_calc_O2COR(dt, Vz, SC)
+ end select
+
+end subroutine advecZ_calc
+
+
+
+! ===== Private procedure =====
+
+! ---- Order 2 solver, with correction for large time step ----
+!> Advection during a time step dt - order 2
+subroutine advecZ_calc_O2COR(dt, Vz, SC)
+
+ use advecVAR
+
+ ! @param[in] dt = time step
+ ! @param[in] Vx = velocity along x (could be discretised on a bigger mesh then the scalar)
+ ! @param[in,out] SC = scalar field to advect
+
+ real(WP), intent(in) :: dt
+ real(WP), dimension(Nx_V, Ny_V, Nz_V), intent(in) :: Vz
+ real(WP), dimension(Nx_sc,Ny_sc,Nz_sc), intent(inout) :: SC
+
+
+ integer :: i,j,k ! indice of the currend mesh point
+ integer :: n ! indice of the currend particle
+ integer :: nb_tag ! number of tagged particles on the current line
+ integer :: nb_no_tag ! number of untagged particles on the current line
+ integer :: nb_tag_total ! total number of tagged
+ integer :: nb_part ! number of particles on the current line
+ integer :: nb_part_total! total number of particles
+
+ real(WP), dimension(Nz_sc) :: p_pos ! particles position
+ real(WP), dimension(Nz_Sc) :: p_SC ! scalar field advect by the particles
+ real(WP), dimension(:), allocatable :: p_V ! particles velocity
+ real(WP), dimension(:), allocatable :: p_pos_0 ! initial position of particles before any advection
+ integer, dimension(:), allocatable :: p_type ! does the particle belong to a center block or a left one ?
+ integer, dimension(:), allocatable :: p_ind_tag ! indice of tagged particles
+ integer, dimension(:), allocatable :: p_ind_no_tag ! indice of untagged particles
+
+ nb_part_total = 0
+
+ ! Some particles on the end of bloc need specific remeshing formula. To identify them, we tag them.
+ nb_tag = 0 ! number of tagged particles on the current line !
+ nb_tag_total = 0 ! total number of tagged particles
+
+ do j=1,Ny_sc
+ do i=1,Nx_sc
+ nb_part=0
+ ! Creation of particles
+ call countZ(nb_part, i, j, p_pos, p_SC, SC)
+ ! If no particles, no more computations
+ if (nb_part.ne.0) then
+ ! Dynamical allocation of tables
+ allocate(p_V(nb_part))
+ allocate(p_type(nb_part))
+ allocate(p_ind_tag(nb_part))
+ allocate(p_ind_no_tag(nb_part))
+ allocate(p_pos_0(nb_part))
+ ! Advection
+ ! A RK2 scheme is used :
+ ! z(t_n+1)= z(t_n) + dt*v[z(t_n)+dt/2*v(z(t_n))].
+ ! Computation of the intermediar point used in the RK scheme
+ call veloZ(nb_part,i,j,p_pos, p_V, Vz)
+ ! XXX la boucle peut être remplacée par des manips directement
+ ! sur les tableaux.
+ do n=1,nb_part
+ ! Save the initial position
+ p_pos_0(n)=p_pos(n)
+ ! Update the position to compute the intermediar point used in
+ ! the RK2 scheme
+ p_pos(n)=p_pos(n)+0.5*dt*p_V(n)
+ ! If the velocity is too high, the particles go out of the box.
+ ! Thus we use the periodical condition.
+ ! We assume that dt*abs(V) < xmax-xmin.
+ if (p_pos(n).gt.zmax) p_pos(n)=p_pos(n)-zmax+zmin
+ if (p_pos(n).lt.zmin) p_pos(n)=p_pos(n)+zmax-zmin
+ end do
+ ! Computation of the "final" velocity
+ call veloZ(nb_part,i,j, p_pos, p_SC, Vz)
+
+ ! As the velocity is known, we could now determinate :
+ ! - the bloc type (right or left)
+ ! - if the particles on a block boundary need a specific remeshing
+ ! formula (then the particles will be tagged)
+ ! All this is done by the procedure tag_particles.
+ call tag_particles_O2(nb_part, Nz_sc, dz_sc, zmin, zmax, p_pos_0, &
+ & p_V, p_type, nb_tag, nb_no_tag, p_ind_tag, p_ind_no_tag)
+ nb_tag_total = nb_tag_total + nb_tag
+
+ ! End of the advection
+ do n=1, nb_part
+ p_pos(n)=p_pos_0(n)*dt*p_V(n)
+ ! If the velocity is too high, the particles go out of the box.
+ if (p_pos(n).gt.zmax) p_pos(n)=p_pos(n)-zmax+zmin
+ if (p_pos(n).lt.zmin) p_pos(n)=p_pos(n)+zmax-zmin
+ end do
+ ! it is possible to write : p_pos(1:nb_part)=p_pos+dt*p_V
+ ! do not forget that size of p_pos_0, p_V = nb_part but
+ ! size(p_pos)=Nz_sc
+ deallocate(p_pos_0)
+ deallocate(p_V)
+ ! Remeshing the particles
+ ! For not tagged particles :
+ call remeshZ_O2(p_pos, p_SC, nb_no_tag, p_ind_no_tag, p_type, j,k, SC)
+ ! For tagged particles :
+ if (nb_tag.ne.0) call remeshZ_O2_tag(p_pos, p_SC, nb_tag, p_ind_tag, p_type,j,k, SC)
+
+ nb_part_total = nb_part_total + nb_part
+ end if
+ end do
+ end do
+
+ !print*, 'NPART, NTAG ', npart,ntag_total
+ !print*,'ntag advecZ :', ntag_total
+
+ deallocate(p_type)
+ deallocate(p_ind_tag)
+ deallocate(p_ind_no_tag)
+
+end subroutine advecZ_calc_O2COR
+
+
+!> To count the number of particles on a line
+subroutine countZ(nb_part, i, j, p_pos, p_SC, SC)
+
+ use advecVAR
+
+ ! This method does only depend on the block size and not of the order of the
+ ! method used
+
+ ! @param[out] nb_part = number of particles
+ ! @param[in] i = x-indice of the line
+ ! @param[in] j = y-indice of the line
+ ! @param[out] p_pos = particles position
+ ! @param[out] p_SC = value of scalar field on the particles
+ ! @param[in] SC = scalar field
+
+ integer, intent(out) :: nb_part
+ integer, intent(in) :: i,j
+ real(WP), dimension(Nz_sc), intent(out) :: p_pos, p_SC
+ real(WP), dimension(Nx_sc,Ny_sc,Nz_sc), intent(in) :: SC
+
+ integer :: increment ! increment ensure reading each mesh
+ ! point one time (as we can add more than one particles)
+ integer :: stencil ! number of particles to add when we
+ ! when we are at a block end
+ integer :: k ! z-indice
+ integer :: ind ! temp. int used in some "do"
+
+ ! The following variables are defined in another module
+ ! min_SC = minimal value of the scalar needed to create a particles (if the scalar is null, there's nothing to advect)
+
+ nb_part = 0
+ increment = 1
+ stencil = 2*nb_bound_part_bl
+
+ ! There two cases
+ ! "Interior" particles (ie inside a block) : could be isolated
+ ! Particles on a block boundary : as we use specific remshing formula, they
+ ! could not be isolated to ensure consistance
+ k = nb_bound_part_bl ! indice in the table = k+1
+ do while(k < Nz_sc+nb_bound_part_bl-1)
+ ! On a block boundary : either no particles either particle on all the
+ ! end of the block and all the beginning of the following one
+ ! end of block <-> mod(i,nb_part_bl) +1 = nb_part_bl - nb_bound_part_bl +1
+ if (mod(k,nb_part_bl)== nb_part_bl - nb_bound_part_bl) then
+ ! as anounced, if the scalar is big enough we create particles at the block end and at
+ ! the next block begining. Therefor, particles are created if the average scalar value
+ ! on these mesh point is big enough. The "sum" is used in the purpose.
+ if (sum((/(abs(SC(i,j,mod(k+ind,Nz_sc)+1)),ind=0,stencil-1)/))>sc_moy*stencil) then
+ do ind = 0, stencil-1
+ nb_part = nb_part +1
+ p_pos(nb_part) = zmin+float(mod(k+ind,Nz_sc))*dz_sc
+ p_SC(nb_part) = SC(i,j,mod(k+ind,Nz_sc)+1)
+ end do
+ end if
+ increment = stencil
+ ! Inside a block
+ ! NB : if we use the minimal block size, there is no "inside" particles
+ else
+ if (abs(SC(i,j,k)).gt.sc_min) then
+ nb_part = nb_part +1
+ p_pos(nb_part) = zmin+float(mod(k+ind,Nz_sc))*dz_sc
+ p_SC(nb_part) = SC(i,j,mod(k+ind,Nz_sc)+1)
+ end if
+ increment = 1
+ end if
+ k = k + increment
+ end do
+
+end subroutine countZ
+
+!> Compute the velocity on each particles of a line of indice (i_sc, j_sc) on the particles grid
+subroutine veloZ(nb_part,j_sc,k_sc, p_pos, p_V, Vz)
+
+ use advecVAR
+
+ ! @param[in] nb_part = number of particles on each line
+ ! @param[in] i = indice corresponding to the line x-coordinate on the scalar grid
+ ! @param[in] j = indice corresponding to the line y-coordinate on the scalar grid
+ ! @param[in] p_pos = particles position
+ ! @param[out] p_V = particles velocity
+ ! @param[in] Vz = velocity along z-axis
+
+ integer, intent(in) :: nb_part, j_sc, k_sc ! see above
+ real(WP), dimension(Ny_sc), intent(in) :: p_pos
+ real(WP), dimension(nb_part), intent(out) :: p_V
+ real(WP), dimension(Nx_V,Ny_V,Nz_V), intent(in) :: Vz
+
+
+ integer :: i_Vm, i_Vp ! mesh point on the velocity grid such as:
+ ! x(i_Vm) <= x(i_sc) < x(i_Vm+1=i_Vp)
+ integer :: j_Vm, j_Vp ! as i_Vm, i_Vp but along the y-axis
+ integer :: k_Vm, k_Vp ! as j_Vm, j_Vp but along the z-axis
+ integer :: k
+
+ real(WP) :: x, y ! coordinate of the line
+ real(WP) :: xxm, yym ! adimentioned distance from the previous
+ ! mesh point : xxm = [x(j_sc)-x(j_Vm)]/dx_V
+ real(WP) :: xxp, yyp ! adimentioned distance from the next mesh point
+ real(WP) :: am, ap, bm, bp, cm, cp ! weight for the interpolation (a along the x-axis,
+ ! b along the y-axis and c along the z-axis)
+
+ real(WP) :: z, zzm, zzp ! as above but for one particles of the line
+
+
+
+ ! Computation of the (x- and y-) coordinate of the line
+ x=+float(j_sc-1)*dx_sc ! = distance from xmin (could be different from x-coordinate if xmin != 0)
+ y=+float(k_sc-1)*dy_sc ! = distance from ymin
+
+ i_Vm = int(x/dx_V)
+ i_Vp = i_Vm + 1
+ j_Vm = int(y/dy_V)
+ j_Vp = j_Vm + 1
+
+ xxm = (x-(i_Vm-1.)*dx_V)/dx_V
+ xxp = 1.0 - xxm
+ yym = (y-(j_Vm-1.)*dy_V)/dy_V
+ yyp = 1.0 - yym
+
+ am=xxp
+ ap=xxm
+ bm=yyp
+ bp=yym
+
+ do k = 1,nb_part
+
+ z = p_pos(k)
+
+ k_Vm = int((z-zmin)/dz_V)
+ k_Vp = k_Vm + 1
+
+ zzm = (z-(k_Vm-1.)*dz_V)/dz_V
+ zzp = 1.0 - zzm
+
+ ! In fortran table indice start from one. If a particles is out the domain, we use periodic boundaries conditions
+ k_Vm = mod(k_Vm+Nz_V, Nz_V)+1
+ k_Vp = mod(k_Vp+Nz_V, Nz_V)+1
+
+
+ cm = zzp
+ cp = zzm
+
+ ! velocity interpolation
+ p_V(k)= Vz(i_Vm,j_Vm,k_Vm)*am*bm*cm
+ p_V(k)= p_V(k) + Vz(i_Vp,j_Vm,k_Vm)*ap*bm*cm
+ p_V(k)= p_V(k) + Vz(i_Vm,j_Vp,k_Vm)*am*bp*cm
+ p_V(k)= p_V(k) + Vz(i_Vp,j_Vp,k_Vm)*ap*bp*cm
+ p_V(k)= p_V(k) + Vz(i_Vm,j_Vm,k_Vp)*am*bm*cp
+ p_V(k)= p_V(k) + Vz(i_Vp,j_Vm,k_Vp)*ap*bm*cp
+ p_V(k)= p_V(k) + Vz(i_Vm,j_Vp,k_Vp)*am*bp*cp
+ p_V(k)= p_V(k) + Vz(i_Vp,j_Vp,k_Vp)*ap*bp*cp
+
+ end do
+
+end subroutine veloZ
+
+!> remeshing with an order 2 method, corrected to allow large CFL number - untagged particles
+subroutine remeshZ_O2(p_pos, p_sc, nb_no_tag, p_ind_no_tag, p_type,i,j,SC)
+
+ use advecVAR
+
+ ! @param[in] p_pos = particles position
+ ! @param[in] p_sc = scalar advected by the particles
+ ! @param[in] nb_no_tag = number of untagged particles
+ ! @param[in] p_ind_tag = indice of untagged particles
+ ! @param[in] i,j = indice of of the current line (x-coordinate and y-coordinate)
+ ! @param[in] p_type = equal 0 (resp 1) if the particle belong to a centered (resp left) block
+ ! @param[in,out] SC = scalar field
+
+ integer, intent(in) :: nb_no_tag, i, j
+ integer, dimension(:), intent(in) :: p_ind_no_tag, p_type
+ real(WP), dimension(:), intent(in) :: p_pos, p_sc
+ real(WP), dimension(Nx_sc,Ny_sc,Nz_sc), intent(inout) :: SC
+
+ integer :: k, ind ! indice of the current untagged particles and
+ ! number of the current particles in all
+ ! the (untagged and tagged) particles
+ real(WP) :: pos, sca ! particle position and scalar value associated
+ integer :: kM, k0, kP ! indice of the the mesh point use in the remeshing
+ ! (they depend on the block type)
+ real(WP) :: zM, z0, zP, c0, cM, cP ! interpolction weight
+
+
+
+ SC(i,j,:)=0.
+
+
+ do k = 1,nb_no_tag
+ ind = p_ind_no_tag(k)
+ sca = p_sc(ind)
+ pos = p_pos(ind)
+
+ ! Mesh point used in remeshing formula
+ if (p_type(ind)==0) then
+ ! the particle belong to a left bloc
+ k0 = int((pos-zmin)/dz_sc)
+ else
+ ! the particle belong to a left bloc
+ k0 = nint((pos-zmin)/dz_sc)
+ end if
+ kM = k0 - 1
+ kP = k0 + 1
+
+ ! Distance to mesh points
+ z0 = (pos - float(k0)*dz_sc-zmin)/dz_sc
+ zM = z0-1
+ zP = z0+1
+
+ kM=mod(kM+Nz_sc,Nz_sc) +1
+ k0=mod(k0+Nz_sc,Nz_sc) +1
+ kP=mod(kP+Nz_sc,Nz_sc) +1
+
+ ! Interpolation weights
+ cM=0.5*z0*zM
+ c0=1.-z0**2
+ cP=0.5*z0*zP
+
+ ! Remeshing
+ SC(kM,j,k) = SC(kM,j,k) + sca*cM
+ SC(k0,j,k) = SC(k0,j,k) + sca*c0
+ SC(kP,j,k) = SC(kP,j,k) + sca*cP
+
+ end do
+
+end subroutine remeshZ_O2
+
+!> remeshing with an order 2 method, corrected to allow large CFL number - tagged particles
+subroutine remeshZ_O2_tag(p_pos, p_sc, nb_tag, p_ind_tag, p_type,i,j,SC)
+
+ use advecVAR
+
+ ! @param[in] p_pos = particles position
+ ! @param[in] p_sc = scalar advected by the particles
+ ! @param[in] nb_no_tag = number of tagged particles
+ ! @param[in] p_ind_tag = indice of tagged particles
+ ! @param[in] i,j = indice of of the current line (y-coordinate and z-coordinate)
+ ! @param[in] ptype = equal 0 (resp 1) if the particle belong to a centered (resp left) block
+ ! @param[in,out] SC = scalar field
+
+ integer, intent(in) :: nb_tag, i, j
+ integer, dimension(:), intent(in) :: p_ind_tag, p_type
+ real(WP), dimension(:), intent(in) :: p_pos, p_sc
+ real(WP), dimension(Nx_sc,Ny_sc,Nz_sc), intent(inout) :: SC
+
+ integer :: k, ind, indP ! indice of the current particles in the group of tagged
+ ! particles and in the group of all particles.
+ ! ind1 is the indice of the next (tagged) particles
+ real(WP) :: pos, sca ! particle position and scalar value associated
+ real(WP) :: posP, scaP ! position of the next particle and scalar value associated
+ integer :: k0, kM, kP, kP2, kP3 ! indice of the the nearest mesh points
+ ! (they depend on the block type)
+ integer :: k0_bis ! indice of the the nearest mesh point for the indP=ind+1 particle
+ real(WP) :: zM, z0, zP, aM, a0, aP, aP2 ! interpolation weight for the particles
+ real(WP) :: zM_bis, z0_bis, zP_bis ! distance to mesh points
+ real(WP) :: bM, b0, bP, bP2, bP3 ! interpolation weight for the next particle (ind+1=indP)
+
+ ! The particles are tagged by group of 2 : the one at a block end and the
+ ! following on the beginning of the next block
+ do k=1,nb_tag-1,2
+ ind = p_ind_tag(k)
+ pos = p_pos(ind)
+ sca = p_sc(ind)
+ indP=p_ind_tag(k+1)
+ posP= p_pos(indP)
+ scaP = p_sc(indP)
+ if (p_type(ind)==0) then
+ ! the current block is a left one. As it is tagged, the following
+ ! one is a centered one
+ ! Indice of mesh point used to remesh
+ k0 = int((pos-zmin)/dz_sc)
+ k0_bis = nint((posP-zmin)/dz_sc)
+ kM=k0-1
+ kP=k0+1
+ kP2=k0+2
+ kP3=k0+3
+ kM=mod(kM+Nz_sc,Nz_sc) +1
+ k0=mod(k0+Nz_sc,Nz_sc) +1
+ kP=mod(kP+Nz_sc,Nz_sc) +1
+ kP2=mod(kP2+Nz_sc,Nz_sc) +1
+ kP3=mod(kP3+Nz_sc,Nz_sc) +1
+
+ ! Distance to mesh point
+ z0 = (pos - float(k0)*dz_sc-zmin)/dz_sc
+ z0_bis = (posP - float(k0_bis)*dz_sc-zmin)/dz_sc
+ zM=z0-1
+ zP=z0+1
+ zP_bis=z0_bis+1
+ zM_bis=z0_bis-1
+
+ ! Interpolation weights
+ aM=0.5*z0*zM
+ a0=1.-z0**2
+ aP=z0
+ aP2=0.5*zM*z0
+ b0=zP_bis*z0_bis/2.
+ bP=-z0_bis
+ bP2=1.-z0_bis**2
+ bP3=z0_bis*zP_bis/2.
+
+ ! Remeshing
+ SC(i,j,kM)= SC(i,j,kM) +aM*sca
+ SC(i,j,k0)= SC(i,j,k0) +a0*sca+b0*scaP
+ SC(i,j,kP)= SC(i,j,kP) +aP*sca+bP*scaP
+ SC(i,j,kP2)=SC(i,j,kP2)+aP2*sca+bP2*scaP
+ SC(i,j,kP3)=SC(i,j,kP3)+bP3*scaP
+
+ ! if the particules belong to a left block (and therefore the following
+ ! particles belong to a centered one - verwise it wouldn't be tagged)
+ else
+ ! Compute the weight
+ k0 = nint((pos-zmin)/dz_sc)
+ k0_bis = int((posP-zmin)/dz_sc)
+ kM=k0-1
+ kP=k0+1
+ kM=mod(kM+Nz_sc,Nz_sc) +1
+ k0=mod(k0+Nz_sc,Nz_sc) +1
+ kP=mod(kP+Nz_sc,Nz_sc) +1
+
+ z0 = (pos - float(k0)*dz_sc-zmin)/dz_sc
+ z0_bis = (posP - float(k0_bis)*dz_sc-xmin)/dz_sc
+ zM = z0-1
+ zP_bis = z0_bis+1
+
+
+ aM=0.5*z0*zP
+ a0=1-aM
+ bP=0.5*z0_bis*zP_bis
+ b0=1-bP
+
+ ! Update the scalar field
+ SC(i,j,kM)=SC(i,j,kM)+aM*sca
+ SC(i,j,k0)=SC(i,j,k0)+a0*sca+b0*scaP
+ SC(i,j,kP)=SC(i,j,kP)+bP*scaP
+ endif
+ enddo
+
+end subroutine remeshZ_O2_tag
+
+end module advecZ
--
GitLab