!> \file resol_adv_diff_2D-sept2009.f90
!!Resoud l equation adv-diffusion par une methode de relaxation
!<

!> \namespace reso_adv_diff_2D_vect
!! Resoud l equation adv-diffusion par une methode de relaxation
!! @note   - passage du vecteur et de l epaisseur  en dummy 
!! @note   - epaisseur prescrite selon un masque 
!! @note use module3D_phy
!! @todo essayer dans le futur la methode de Richard dUX/dx = H dU/dx + U dH/dx
!<

module reso_adv_diff_2D_vect

use module3D_phy



implicit none
real :: omega   !< parametre schema temporel de la resolution partie diffusion
real :: mu_adv  !< parametre schema temporel de la resolution advection
real :: upwind  !< schema spatial pour l'advection

! tableaux de travail. resolution M H = Frelax
real,dimension(nx,ny) :: crelax      !< diagnonale de M
real,dimension(nx,ny) :: arelax      !< sous diagonale selon x
real,dimension(nx,ny) :: brelax      !< sur diagonale selon x
real,dimension(nx,ny) :: drelax      !< sous diagonale selon y
real,dimension(nx,ny) :: erelax      !< sur diagonale selon y
real,dimension(nx,ny) :: frelax      !< vecteur
real,dimension(nx,ny) :: c_west      !< sur demi mailles Ux
real,dimension(nx,ny) :: c_east      !< sur demi mailles Ux
real,dimension(nx,ny) :: c_north     !< sur demi mailles Uy
real,dimension(nx,ny) :: c_south     !< sur demi mailles Uy

real,dimension(nx,ny) :: bdx         !< pente socle
real,dimension(nx,ny) :: bdy         !< pente socle

real,dimension(nx,ny) :: hdx         !< pente epaisseur
real,dimension(nx,ny) :: hdy         !< pente epaisseur




integer, dimension(2) :: ijmax    ! position de maxval
integer :: iFAIL

contains 

!> initialise  init_reso_adv_diff_2D
!! definition des parametres qui gerent la resolution
!! @return omega, mu_adv, upwind

subroutine init_reso_adv_diff_2D


write(num_rep_42,*)'partie diffusive'
write(num_rep_42,*)'le type de schema temporel diffusif depend de omega'
write(num_rep_42,*)'0 -> explicite, 0.5 -> Crank-Nicolson'
write(num_rep_42,*)'1 -> semi implicite,  >1 -> over-implicite' 

! over implicite propose par Hindmasrsh

omega = 2.5

write(num_rep_42,*)'omega = ',omega

write(num_rep_42,*)
write(num_rep_42,*)'partie advective'
write(num_rep_42,*)' le schéma temporel advectif dépend de mu'
write(num_rep_42,*)'0 -> explicite, 0.5 -> Crank-Nicolson'
write(num_rep_42,*)'1 -> semi implicite,  >1 -> over-implicite'

mu_adv=1.  ! l'over implicite ne s'applique pas a l'equation advective
write(num_rep_42,*)'mu_adv = ',mu_adv

write(num_rep_42,*)' le schéma spatial dépend de upwind'
write(num_rep_42,*)'1 -> schema upwind,  0.5 -> schema centre'

upwind=1

write(num_rep_42,*)'upwind = ',upwind


write(num_rep_42,*)'------------------------------------------------------'

return
end subroutine init_reso_adv_diff_2D
!------------------------------------------------------------------

!> subroutine resol_adv_diff_2D_Hpresc
!! definition des parametres qui gerent la resolution
!! @param Dfx,Dfy             termes diffusifs
!! @param advx,advy           termes advectifs
!! @param vieuxH              ancienne valeur de H
!! @param H_presc             la valeur H prescrit pour certains noeuds
!! @param i_Hpresc            le masque ou la valeur de H est prescrite
!! @param bil                 le bilan pour la colonne 
!! @return newH               la valeur de H apres le pas de temps

subroutine resol_adv_diff_2D_vect(Dfx,Dfy,advx,advy,H_presc,i_Hpresc,bil,vieuxH,newH)

!$ USE OMP_LIB
implicit none

real,dimension(nx,ny), intent(in) :: Dfx          !< terme diffusif selon x
real,dimension(nx,ny), intent(in) :: Dfy          !< terme diffusif selon y
real,dimension(nx,ny), intent(in) :: Advx         !< terme advectif selon x
real,dimension(nx,ny), intent(in) :: Advy         !< terme advectif selon y
real,dimension(nx,ny), intent(in) :: vieuxH       !< ancienne valeur de H
real,dimension(nx,ny), intent(out):: newH         !< nouvelle valeur de H

real,dimension(nx,ny), intent(in) :: bil          !< bilan de masse pour la colonne

real,dimension(nx,ny), intent(in) :: H_presc      !< H value if prescribed
integer,dimension(nx,ny), intent(in) :: i_Hpresc  !< 1 if H is prescribed on this node, else 0


! tableaux de travail. resolution M H = Frelax
!!$real,dimension(nx,ny) :: crelax      !< diagnonale de M
!!$real,dimension(nx,ny) :: arelax      !< sous diagonale selon x
!!$real,dimension(nx,ny) :: brelax      !< sur diagonale selon x
!!$real,dimension(nx,ny) :: drelax      !< sous diagonale selon y
!!$real,dimension(nx,ny) :: erelax      !< sur diagonale selon y
!!$real,dimension(nx,ny) :: frelax      !< vecteur
!!$real,dimension(nx,ny) :: c_west      !< sur demi mailles Ux
!!$real,dimension(nx,ny) :: c_east      !< sur demi mailles Ux
!!$real,dimension(nx,ny) :: c_north     !< sur demi mailles Uy
!!$real,dimension(nx,ny) :: c_south     !<sur demi mailles Ux

!!$real,dimension(nx,ny) :: bdx         !< pente socle
!!$real,dimension(nx,ny) :: bdy         !< pente socle
!!$
!!$real,dimension(nx,ny) :: hdx         !< pente epaisseur
!!$real,dimension(nx,ny) :: hdy         !< pente epaisseur

real :: frdx,frdy                    !< pour calcul frelax : termes diffusion
real :: fraxw,fraxe,frays,frayn      !< termes advection

real,dimension(nx,ny) :: deltah      ! dans calcul relax
real :: delh                         ! dans calcul relax
real :: testh                        ! dans calcul relax

logical :: stopp
integer :: ntour
real :: reste


if (itracebug.eq.1)  call tracebug(' Entree dans routine resolution_diffusion')

!$OMP PARALLEL
!$OMP WORKSHARE
! calcul de bdx et hdx
hdx(:,:)=dx1*(vieuxH(:,:)-eoshift(vieuxH(:,:),shift=-1,boundary=0.,dim=1))
hdy(:,:)=dx1*(vieuxH(:,:)-eoshift(vieuxH(:,:),shift=-1,boundary=0.,dim=2))
bdx(:,:)=dx1*(B(:,:)-eoshift(B(:,:),shift=-1,boundary=0.,dim=1))
bdy(:,:)=dx1*(B(:,:)-eoshift(B(:,:),shift=-1,boundary=0.,dim=2))


! initialisations (qui feront aussi les conditions aux limites)
Arelax(:,:)=0.
Brelax(:,:)=0.
Drelax(:,:)=0.
Erelax(:,:)=0.
Crelax(:,:)=1.
Frelax(:,:)=0.
DeltaH(:,:)=0.
!$OMP END WORKSHARE

! schema spatial

if (upwind.eq.1.) then                 !schema amont
!$OMP WORKSHARE
   where (Advx(:,:).ge.0.)
      c_west(:,:)=1.
      c_east(:,:)=0.
   elsewhere
      c_west(:,:)=0.
      c_east(:,:)=1.
   end where

   where (Advy(:,:).ge.0.)
      c_south(:,:)=1.
      c_north(:,:)=0.
   elsewhere
      c_south(:,:)=0.
      c_north(:,:)=1.
   end where
!$OMP END WORKSHARE
else if (upwind.lt.1.) then             ! schema centre
!$OMP WORKSHARE
      c_west(:,:)=0.5
      c_east(:,:)=0.5
      c_south(:,:)=0.5
      c_north(:,:)=0.5
!$OMP END WORKSHARE      
end if


! attribution des elements des diagonales
!$OMP DO PRIVATE(frdx,frdy,fraxw,fraxe,frays,frayn)
do j=2,ny-1
  do i=2,nx-1

!  sous diagonale en x
     arelax(i,j)=-omega*Dtdx2*Dfx(i,j)   &               ! partie diffusive en x
          -mu_adv*dtdx*c_west(i,j)*Advx(i,j)             ! partie advective en x

!  sur diagonale en x
     brelax(i,j)=-omega*Dtdx2*Dfx(i+1,j)  &              ! partie diffusive
          +mu_adv*dtdx*c_east(i+1,j)*Advx(i+1,j)         ! partie advective

!  sous diagonale en y
     drelax(i,j)=-omega*Dtdx2*Dfy(i,j)   &               ! partie diffusive en y
          -mu_adv*dtdx*c_south(i,j)*Advy(i,j)            ! partie advective en y

!  sur diagonale en y
     erelax(i,j)=-omega*Dtdx2*Dfy(i,j+1)  &              ! partie diffusive
          +mu_adv*dtdx*c_north(i,j+1)*Advy(i,j+1)        ! partie advective



! diagonale
     crelax(i,j)=omega*Dtdx2*((Dfx(i+1,j)+Dfx(i,j))  &                         ! partie diffusive en x
                             +(Dfy(i,j+1)+Dfy(i,j)))                           ! partie diffusive en y
     crelax(i,j)=crelax(i,j)+mu_adv*dtdx* &
                  ( (c_west(i+1,j)*Advx(i+1,j)-c_east(i,j)*Advx(i,j)) &        !partie advective en x
                  +(c_south(i,j+1)*Advy(i,j+1)-c_north(i,j)*Advy(i,j)))        !partie advective en y
     crelax(i,j)=1.+crelax(i,j)                                                ! partie temporelle


! terme du vecteur

     frdx= -Dfx(i,j) * (Bdx(i,j)  +(1.-omega)*Hdx(i,j))         &  ! partie diffusive en x
           +Dfx(i+1,j)*(Bdx(i+1,j)+(1.-omega)*Hdx(i+1,j))

     frdy= -Dfy(i,j) * (Bdy(i,j)  +(1.-omega)*Hdy(i,j))         &  ! partie diffusive en y
           +Dfy(i,j+1)*(Bdy(i,j+1)+(1.-omega)*Hdy(i,j+1))

     fraxw= -c_west(i,j)*  Advx(i,j) * vieuxH(i-1,j)           &  ! partie advective en x 
            +c_west(i+1,j)*Advx(i+1,j)*vieuxH(i,j)                ! venant de l'west

     fraxe= -c_east(i,j) * Advx(i,j) * vieuxH(i,j)             &  ! partie advective en x 
            +c_east(i+1,j)*Advx(i+1,j)*vieuxH(i+1,j)              ! venant de l'est

     frays= -c_south(i,j) * Advy(i,j) * vieuxH(i,j-1)          &  ! partie advective en y
            +c_south(i,j+1)*Advy(i,j+1)*vieuxH(i,j)               ! venant du sud

     frayn= -c_north(i,j) * Advy(i,j) * vieuxH(i,j)            &  ! partie advective en y
            +c_north(i,j+1)*Advy(i,j+1)*vieuxH(i,j+1)             ! venant du nord



     frelax(i,j) = vieuxH(i,j) + Dt*bil(i,j)    &
                 + dtdx*(frdx+frdy)        &  ! termes lies a la diffusion (socle, omega)
                 + (1.-mu_adv)*dtdx*((fraxw+fraxe)+(frays+frayn)) 
                                              ! le dernier terme serait pour 
                                              ! de l'over implicite en diffusion
                                              ! (en pratique pas utilise)
     end do
  end do
!$OMP END DO

!$OMP WORKSHARE
where (i_hpresc(:,:) .eq.1)          ! thickness prescribed
   frelax(:,:) = H_presc(:,:)
   arelax(:,:) = 0.
   brelax(:,:) = 0.
   crelax(:,:) = 1.
   drelax(:,:) = 0.
   erelax(:,:) = 0.
end where
!$OMP END WORKSHARE
!$OMP END PARALLEL

stopp = .false.
ntour=0

relax_loop: do while(.not.stopp)
ntour=ntour+1
!$OMP PARALLEL
!$OMP DO PRIVATE(reste)
do j=2,ny-1
   do i=2,nx-1

      reste = (((arelax(i,j)*newH(i-1,j) +drelax(i,j)*newH(i,j-1)) &
           + (brelax(i,j)*newH(i+1,j) + erelax(i,j)*newH(i,j+1))) &
           + crelax(i,j)*newH(i,j))- frelax(i,j)

!if (ntour.eq.1)  debug_3D(i,j,49)=(((arelax(i,j)*newH(i-1,j) +drelax(i,j)*newH(i,j-1)) &
!           + (brelax(i,j)*newH(i+1,j) + erelax(i,j)*newH(i,j+1))) &
!           + crelax(i,j)*newH(i,j))


      deltaH(i,j) = reste/crelax(i,j)

   end do
end do
!$OMP END DO

!debug_3D(:,:,50)=arelax(:,:)
!debug_3D(:,:,51)=brelax(:,:)
!debug_3D(:,:,52)=crelax(:,:)
!debug_3D(:,:,53)=drelax(:,:)
!debug_3D(:,:,54)=erelax(:,:)
!debug_3D(:,:,55)=frelax(:,:)


!$OMP WORKSHARE
newH(:,:)=newH(:,:)-deltaH(:,:)
!$OMP END WORKSHARE
!$OMP END PARALLEL


! critere d'arret:
! ----------------         

delh=0

!$OMP PARALLEL
!$OMP DO REDUCTION(+:delh)
do j=2,ny-1
   do i=2,nx-1
      delh=delh+deltaH(i,j)**2
   end do
end do
!$OMP END DO
!$OMP END PARALLEL

if (delh.gt.0.) then
   testh=sqrt(delh)/((nx-2)*(ny-2))
else
   testh=0.
endif
stopp = (testh.lt.1.e-4).or.(ntour.gt.100)



end do relax_loop


! thickness at the upwind node
!do j = 1, ny
!   do i = 2, nx
!      debug_3D(i,j,92) = c_west(i,j) * newH(i-1,j) + c_east(i,j) * newH(i,j)
!   end do
!end do
!do j = 2, ny
!   do i = 1, nx
!      debug_3D(i,j,93) = c_south(i,j) * newH(i,j-1) + c_north(i,j) * newH(i,j)
!   end do
!end do

if (itracebug.eq.1)  call tracebug(' fin routine resolution_diffusion')

 
return

end subroutine resol_adv_diff_2D_vect

end module reso_adv_diff_2D_vect