source: trunk/SOURCES/Old-sources/diffusiv-polyn-0.5.F90 @ 159

!> \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 module3D_phy

use module_choix

implicit none

REAL :: ulgliss,ultot,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)
INTEGER :: ii
integer :: ll

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

!     la valeur de ff est donnee dans initial
!     limitation de la vitesse de glissement et totale
ulgliss=10000.
ultot=10000.


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

! initialisation de la diffusion
diffmx(:,:)=0.
diffmy(:,:)=0.


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

call slope_surf


! Mise ne réserve des vitesses flgzmx et flgzmy   ! non ca a l'air de se cumuler
where (flgzmx(:,:))
!   uxflgz(:,:)=uxbar(:,:)
elsewhere
   uxflgz(:,:)=0.
endwhere

where (flgzmy(:,:))
!   uyflgz(:,:)=uybar(:,:)
elsewhere
   uyflgz(:,:)=0.
endwhere


! Calcul de ddy  et ddx pour chaque valeur de iglen
!------------------------------------------------------------------------------------


do iglen=n1poly,n2poly
   glenexp=max(0.,(glen(iglen)-1.)/2.)


   do J=1,NY
      do I=1,NX ! -1

         if (.not.flotmy(i,j)) then !  On calcule quand la deformation même pour les ice streams
!         if (.not.flgzmy(i,j)) then !  ancienne version
                                    !  si marine les points de la grounding zone
                                    ! sont traites avec le shelf
                                    !  si pas marine, ils sont traites ici



         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
end do

do iglen=n1poly,n2poly
   glenexp=max(0.,(glen(iglen)-1.)/2.)


   do J=1,NY  !-1
      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
end do




!      -------------------------------------------------
!      GLISSEMENT
!      DDBX et DDBY termes dus au glissement 
!      relation avec la vitesse de glissement UBx et UBy 
!      UBx=-DDBX*SDX et UBy=-DDBY*SDY
!      -------------------------------------------------

! 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.

call sliding

! vitesse de déformation-----------------------------------------------


ydef: do j=2,NY
       do I=2,NX  ! -1


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

           if (abs(sdy(i,j)).gt.10.e-10) then !limitation de la vitesse de glissement
              ddby(i,j)=min(ddby(i,j),abs(ulgliss/sdy(i,j)))
           endif

          uybar(i,j)=ddby(i,j)

         do iglen=n1poly,n2poly

!            uybar(i,j)=uybar(i,j)+ddy(i,j,iglen)*(s2a(i,j,1,iglen) &
!                   +s2a(i,j-1,1,iglen))/2.
            uybar(i,j)=uybar(i,j)+ddy(i,j,iglen)*(s2a_my(i,j,1,iglen))

         end do
         diffmy(i,j)=uybar(i,j)     ! pour utilisation comme diffusion, a multiplier par H
         uybar(i,j)=uybar(i,j)*(-sdy(i,j))

         uby(i,j)=-(ddby(i,j)*sdy(i,j))
         uydef(i,j)=uybar(i,j)-uby(i,j)


         else    ! points flottants
            uydef(i,j)=0.
            uby(i,j)=0.          
         endif
    end do
   end do ydef




xdef: do J=2,NY  !-1
        do I=2,NX

          if (.not.flotmx(i,j)) then

             if (abs(sdx(i,j)).gt.10.e-10) then !limitation de la vitesse de glissement
                ddbx(i,j)=min(ddbx(i,j),abs(ulgliss/sdx(i,j)))
             endif

             uxbar(i,j)=ddbx(i,j)


             do iglen=n1poly,n2poly

!                uxbar(i,j)=uxbar(i,j)+ddx(i,j,iglen)*(s2a(i,j,1,iglen) &
!                        +s2a(i-1,j,1,iglen))/2.
                uxbar(i,j)=uxbar(i,j)+ddx(i,j,iglen)*s2a_mx(i,j,1,iglen)

             end do
             diffmx(i,j)=uxbar(i,j)     ! pour utilisation comme diffusion, a multiplier par H
             uxbar(i,j)=uxbar(i,j)*(-sdx(i,j))
             ubx(i,j)=-(ddbx(i,j)*sdx(i,j))
             uxdef(i,j)=uxbar(i,j)-ubx(i,j)

          else 
             uxdef(i,j)=0.
             ubx(i,j)=0.   
          endif

       end do
    end do xdef

!    Recherche des vitesses >1000m/an 
ll=0
ii=0


!    mise a 0 si pas marine et flottant
if (marine) goto 88
do j=2,ny
   do I=2,nx
      if (.not.marine.and.flotmx(i,j).and..not.gzmx(i,j) ) then
         uxbar(i,j)=0.
         do k=1,nz
            ux(i,j,k)=uxbar(i,j)
         end do
      else if (.not.marine.and.gzmx(i,j)) then
         do k=1,nz
            ux(i,j,k)=uxbar(i,j)
         end do
      endif
      if (.not.marine.and.flotmy(i,j).and..not.gzmy(i,j)) then
         uybar(i,j)=0.
         do k=1,nz
            uy(i,j,k)=uybar(i,j)
         end do
      else if (.not.marine.and.gzmy(i,j)) then
         do k=1,nz
            uy(i,j,k)=uybar(i,j)
         end do
      endif
   end do
end do


88 continue

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

! la limitation selon x
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

!    cas particulier des sommets ou il ne reste plus de neige.
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

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

! On remet les vitesses flgz
where (flgzmx(:,:))
   uxbar(:,:)=uxflgz(:,:)
endwhere

where (flgzmy(:,:))
   uybar(:,:)=uyflgz(:,:)
endwhere
if (itracebug.eq.1)  call tracebug(' fin diffusiv ')

return
end subroutine diffusiv