source: branches/GRISLIv3/SOURCES/diffusiv-polyn-0.6.f90 @ 446

!> \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, only: V_limit,iglen,diffmx,diffmy,flgzmx,flgzmy,ubx,uby,&
       uxflgz,uyflgz,ddbx,ddby,sdx,sdy,flotmx,flotmy,slope2mx,slope2my,hmx,hmy,uxdef,uydef,&
       uxbar,uybar,H,flot
  use geography, only: nx,ny,dx,dy
  use runparam, only: itracebug
  use param_phy_mod, only:rog
  use deformation_mod_2lois, only:n1poly,n2poly,glen,ddx,ddy,s2a_mx,s2a_my
  use module_choix, only:sliding

  implicit none

  integer :: i,j
  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