!> \file diffusiv-polyn-0.5.f90
!! Contient la subroutine de calcul diffusivites
!<

!> SUBROUTINE: diffusiv
!! Calcule de  diffusivites 
!! \author Catherine 
!! \date   18 mars 2006
!! @note Used modules:
!! @note    - use module3D_phy
!! @note    - use module_choix
!!
!<
      subroutine diffusiv()

!     ===============================================================
!     **     DIFFUSIVITIES                    18 mars 2006   ***
!     glissement avec reserve d'eau 29 Avril 97
!     choix sur le type de discretisation (locale iloc=1 )    
!
!     Modification Vince --> 2 fev 95
!     pour couplage avec Shelf
!
!     Modification C. Dumas Dec 99
!     DDX, SA,S2A,... ont une dimension en plus (loi de deformation)
!
!     Modifications  de Cat du 11 juin 2000
!     Modif Christophe 24 janv 2002 : passage au fortran 90
!     pour optimisation de la routine
!     Modif Vincent dec 2003 : on utilise plus frontmx/y > commentaires
!
!     Modif Cat mars 2003
!                  modif mars 2007 : moyennes de S2a avec S2a_mx et S2a_my

!     ===============================================================
!$ USE OMP_LIB
USE module3D_phy

use module_choix

implicit none

REAL :: ulgliss,ultot,uldef,glenexp
REAL :: INV_4DX ! inverse de dx pour eviter les divisions =1/(4*dx)
REAL :: INV_4DY ! inverse de dy pour eviter les divisions =1/(4*dy)

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

!     la valeur de ff est donnee dans initial
!     limitation de la vitesse de glissement de déformation et totale
ulgliss = V_limit
ultot   = V_limit
uldef   = 2000.


INV_4DX=1/(4*dx)
INV_4DY=1/(4*dy)

! initialisation de la diffusion
!$OMP PARALLEL
!$OMP WORKSHARE
diffmx(:,:)=0.
diffmy(:,:)=0.
!$OMP END WORKSHARE
!$OMP END PARALLEL

! calcul de SDX, SDY et de la pente au carre en mx et en my :

call slope_surf

call sliding     ! au sens vitesse de glissement

!  le glissement est maintenant dans un module a part choisi dans le module choix
!  pour l'instant seules les lois Heino (loigliss=4) et Bindshadler(loigliss=2)
!  sont programmees.

!      ddbx et ddby termes dus au glissement 
!      relation avec la vitesse de glissement ubx et uby 
!      ubx=-ddbx*sdx et uby=-ddby*sdy
!      -------------------------------------------------

! Calcul de la vitesse de glissement, qu'elle soit calculee par diagno ou par sliding
!------------------------------------------------------------------------------------
!$OMP PARALLEL
!$OMP WORKSHARE
where (flgzmx(:,:))
    ubx(:,:) = uxflgz(:,:)
elsewhere
    ubx(:,:) = ddbx(:,:)* (-sdx(:,:))     ! on pourrait mettre une limitation
endwhere

where (flgzmy(:,:))
    uby(:,:) = uyflgz(:,:)
elsewhere
    uby(:,:) = ddby(:,:)* (-sdy(:,:))
endwhere
!$OMP END WORKSHARE
!$OMP END PARALLEL
if (itracebug.eq.1) call tracebug(' diffusiv  apres calcul glissement')



! Deformation SIA : Calcul de ddy  et ddx pour chaque valeur de iglen
!------------------------------------------------------------------
do iglen=n1poly,n2poly
   glenexp=max(0.,(glen(iglen)-1.)/2.)
  !$OMP PARALLEL 
  !$OMP DO
   do j=1,ny
      do i=1,nx
         if (.not.flotmy(i,j)) then !  On calcule quand la deformation même pour les ice streams
            ddy(i,j,iglen)=((slope2my(i,j)**  &  ! slope2my calcule en debut de routine
                 glenexp)*(rog)**glen(iglen)) &
                 *hmy(i,j)**(glen(iglen)+1)
         endif
      end do
   end do
   !$OMP END DO
   !$OMP END PARALLEL
end do



do iglen=n1poly,n2poly
   glenexp=max(0.,(glen(iglen)-1.)/2.)
   !$OMP PARALLEL
   !$OMP DO
   do j=1,ny  
      do i=1,nx
         if (.not.flotmx(i,j)) then      ! bug y->x corrige le 5 avril 2008 (
            ddx(i,j,iglen)=((slope2mx(i,j)**  & ! slope2mx calcule en debut de routine
                 glenexp)*(rog)**glen(iglen)) &
                 *hmx(i,j)**(glen(iglen)+1)
         endif
      end do
   end do
   !$OMP END DO
   !$OMP END PARALLEL
end do


! vitesse due a la déformation-----------------------------------------------
!$OMP PARALLEL
!$OMP DO
ydef: do j=2,ny
       do i=1,nx 

        if (.not.flotmy(i,j)) then !  On calcule la deformation même pour les ice streams

         uydef(i,j)=0.
         diffmy(i,j)=0.

         do iglen=n1poly,n2poly
           diffmy(i,j)=diffmy(i,j)+ddy(i,j,iglen)*(s2a_my(i,j,1,iglen))
         end do

         uydef(i,j)=diffmy(i,j)*(-sdy(i,j))    ! partie deformation
         uydef(i,j)=min(uydef(i,j),uldef)
         uydef(i,j)=max(uydef(i,j),-uldef)

         diffmy(i,j)=diffmy(i,j)+ddby(i,j)     ! pour utilisation comme diffusion, a multiplier par H

         uybar(i,j)=uby(i,j)+uydef(i,j)


         else    ! points flottants
            uydef(i,j)=0.   ! normalement uxbed est attribue          
         endif
    end do
   end do ydef
!$OMP END DO

!$OMP DO
xdef: do j=1,ny
       do i=2,nx 

        if (.not.flotmx(i,j)) then !  On calcule la deformation même pour les ice streams

         uxdef(i,j)=0.
         diffmx(i,j)=0.

         do iglen=n1poly,n2poly
           diffmx(i,j)=diffmx(i,j)+ddx(i,j,iglen)*(s2a_mx(i,j,1,iglen))
         end do

         uxdef(i,j)=diffmx(i,j)*(-sdx(i,j))    ! partie deformation
         uxdef(i,j)=min(uxdef(i,j),uldef)
         uxdef(i,j)=max(uxdef(i,j),-uldef)


         diffmx(i,j)=diffmx(i,j)+ddbx(i,j)     ! pour utilisation comme diffusion,
                                               ! a multiplier par H
         uxbar(i,j)=ubx(i,j)+uxdef(i,j)


         else    ! points flottants
            uxdef(i,j)=0.   ! normalement uxbed est attribue          
         endif
    end do
   end do xdef
!$OMP END DO

if (itracebug.eq.1) call tracebug(' diffusiv  avant limit')


! limitation par ultot----------------------------------------------------------

! la limitation selon x
!$OMP WORKSHARE
where(.not.flgzmx(:,:))
   uxbar(:,:)=min(uxbar(:,:),ultot)
   uxbar(:,:)=max(uxbar(:,:),-ultot)
end where

! la limitation selon y
where(.not.flgzmy(:,:))
   uybar(:,:)=min(uybar(:,:),ultot)
   uybar(:,:)=max(uybar(:,:),-ultot)
end where
!$OMP END WORKSHARE

!    cas particulier des sommets ou il ne reste plus de neige.
!$OMP DO
do j=2,ny-1
   do i=2,nx-1
      if ((H(i,j).le.1.).and.(.not.flot(i,j))) then
         uxbar(i+1,j)=min(uxbar(i+1,j),ultot)  ! n'agit que si ux -> 
         uxbar(i,j)=max(uxbar(i,j),-ultot)     ! n'agit que si ux <-
         uybar(i,j+1)=min(uybar(i,j+1),ultot)
         uybar(i,j)=max(uybar(i,j),-ultot)
      endif
   end do
end do
!$OMP END DO
!$OMP END PARALLEL
if (itracebug.eq.1)  call tracebug(' fin diffusiv ')

return
end subroutine diffusiv