Index: unk/SOURCES/Antarctique_general_files/dragging_prescr_beta_mod_Ant.f90
===================================================================
--- /trunk/SOURCES/Antarctique_general_files/dragging_prescr_beta_mod_Ant.f90	(revision 330)
+++ 	(revision )
@@ -1,530 +1,0 @@
-!> \file dragging_prescr_beta_mod.f90
-!! Module pour calculer le beta a partir de vitesses de bilan
-!<
-
-!> \namespace dragging_prescr_beta
-!! Calcule le beta a partir de vitesses de bilan
-!! @note Il faut partir d'un cptr 
-!! \author ...
-!! \date ...
-!! @note Used module
-!! @note   - use module3D_phy
-!<
-
-module dragging_prescr_beta
-
-  ! Calcule le beta a partir de vitesses de bilan
-  ! il faut partir d'un cptr.
-
-  use module3d_phy
-  use interface_input
-  implicit none
-  real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite
-  real, dimension(nx,ny) :: Vcol_y           !<
-  real, dimension(nx,ny) :: Vsl_x            !<
-  real, dimension(nx,ny) :: Vsl_y            !<
-
-
-  real, dimension(nx,ny) :: drag_centre      !< beta on major node (average)
-  !hassine !real, dimension(nx,ny) :: beta_centre      ! beta on major node (average)
-  real, dimension(nx,ny) :: betamoy_x        !< pour faire les moyennes glissantes
-  real, dimension(nx,ny) :: betamoy_y        !<
-  integer, dimension(nx,ny) :: nbvoisx       !<
-  integer, dimension(nx,ny) :: nbvoisy       !<
-
-  real, dimension(nx,ny) ::  corrdS_centre   !< correction dS 
-  real, dimension(nx,ny) ::  corrdS_mx       !<
-  real, dimension(nx,ny) ::  corrdS_my       !<
-  real, dimension(nx,ny) ::  fleuves_nommes  !< each ice stream is numbered
-
-  real, dimension(nx,ny) ::  corrstream_centre    !< pour correction stream
-  real, dimension(nx,ny) ::  corrstream_mx        !<
-  real, dimension(nx,ny) ::  corrstream_my        !<
-  integer :: numstream
-
-  integer,parameter             :: n_g = 2        !< taille de cette matrice (-n_g:n_g,-n_g:n_g)
-  real, dimension(-2:2,-2:2)    :: M_gauss        !< matrice qui contient la gaussienne pour iterations beta
-  real                          :: sigma_beta=1   !< sigma de la Gaussienne
-
-
-  real, dimension(nx,ny) :: drag_instab_topo !< regions detected as topo instability
-  real :: beta_limgz                         !< when beta gt beta_limgz -> not gzmx
-  real :: beta_min                           !< minimum value of beta
-  real :: beta_mult                          !< coefficient of beta field
-  integer  :: ill,jll,nmoy
-
-  real :: maxi                               !< calcul correction dS a la grounding line
-  real :: mini
-
-  logical :: corr_def = .False.               !< for deformation correction (compatibility)
-
-contains
-
-  !-----------------------------------------------------------------------------------------
-  !> SUBROUTINE: init_dragging 
-  !! Cette routine fait l'initialisation du dragging.
-  !<
-  subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
-                                ! nouvelle version : lit les fichiers nc
-    implicit none
-    character(len=100) :: beta_c_file               !  beta on centered grid
-    character(len=100) :: betamx_file               !  beta on staggered grid mx
-    character(len=100) :: betamy_file               !  beta on staggered grid mx
-
-    namelist/beta_prescr/beta_c_file,beta_limgz,beta_min,beta_mult
-
-    if (itracebug.eq.1)  call tracebug(' Prescr_beta       subroutine init_dragging')
-
-    iter_beta=0
-
-    ! lecture des parametres du run                    
-    !--------------------------------------------------------------------
-
-    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
-428 format(A)
-    read(num_param,beta_prescr)
-
-    write(num_rep_42,428) '!___________________________________________________________' 
-    write(num_rep_42,428) '!  read beta on centered grid'
-    write(num_rep_42,beta_prescr)                    
-    write(num_rep_42,428) '!  beta_file : nom des fichier qui contiennent les beta_c'
-    write(num_rep_42,428) '!  above beta_limgz, gzmx is false'
-    write(num_rep_42,428) '!  if grounded, beta_min minimum value '
-    write(num_rep_42,428) '!  beta_mult : coefficient just after reading' 
-    write(num_rep_42,*)
-    write(num_rep_42,428) '!___________________________________________________________' 
-   
-
-    beta_c_file = trim(dirnameinp)//trim(beta_c_file)
-    betamx_file = trim(dirnameinp)//trim(betamx_file)
-    betamy_file = trim(dirnameinp)//trim(betamy_file)
-    betamax = beta_limgz
-
-!  read the beta value on centered and staggered grids
-!    call lect_input(1,'beta_c',1,drag_centre,beta_c_file,trim(dirnameinp)//trim(runname)//'.nc')
-!    call lect_input(1,'betamx',1,drag_mx,betamx_file,trim(dirnameinp)//trim(runname)//'.nc')
-!    call lect_input(1,'betamy',1,drag_my,betamy_file,trim(dirnameinp)//trim(runname)//'.nc')
-
-
-    call lect_input(1,'beta_c',1,drag_centre,beta_c_file,'')
-
-! multiplication par le coefficient beta_mult
-    if (itracebug.eq.1)  call tracebug(' Prescr_beta       apres lecture')
-    drag_centre(:,:) = drag_centre(:,:)*beta_mult
-
-
-!    call lect_input(1,'betamx',1,drag_mx,betamx_file,'')
-!    call lect_input(1,'betamy',1,drag_my,betamy_file,'')
-
-
-     call  beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre)         ! redistribute on staggered grid
-
-     if  (beta_limgz.gt.0.) then
-        beta_centre(:,:) =  drag_centre(:,:)
-        betamx(:,:)      =  drag_mx(:,:)
-        betamy(:,:)      =  drag_my(:,:)
-
-        if (itracebug.eq.1)  call tracebug(' Prescr_beta       apres call beta_c2stag')
-
-     else if (beta_limgz.lt.0.) then                              ! attention sans doute pour plastic
-
-        beta_centre(:,:) = - beta_limgz
-        betamx(:,:)      = - beta_limgz
-        betamy(:,:)      = - beta_limgz
-        drag_centre(:,:) = - beta_limgz 
-        drag_mx(:,:)     = - beta_limgz
-        drag_my(:,:)     = - beta_limgz 
-        beta_limgz = 1.
-     end if
-
-
-! divers traitements qui peuvent etre appliques a beta : mis sous forme de routine
-
-    ! call beta_plateau_slope   !  give value beta_limgz to places where slope is >1.e-2
-    ! call beta_stag2c          !  average the staggered values on central ones
-    ! call beta_plateau_S       !  give value beta_limgz to places where S > 2500 m
-
-   call beta_min_value(beta_min)      !  valeur mini que peut avoir beta (en Pa an m-1)
-                                ! on peut envisager une valeur ~ 10 
-                                ! rappel : beta doit être positif.
-
-    ! call  beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre)         ! redistribute on staggered grid
-
-
-    !  mstream attribution
-    !________________________
-    ! mstream -> regions where streaming is allowed
-    ! drag_mx -> value of beta if an ice stream is allowed
-    ! betamx  -> value of dragging (may depend on other conditions)
-
-    ! to ensure having some "Dirichlet" points, there must be at least 
-    ! a mstream=0 in every island
-
-
-! un tour de lissage
-!    call beta_stag2c(nx,ny,drag_mx,drag_my,drag_centre)                  ! average  staggered 
-!                                                                         ! on central 
-
-!    call beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre)                  ! centered distributed
-!                                                                         ! on staggered
-!
-    beta_centre(:,:) =  drag_centre(:,:)                                 ! surtout pour sorties
-
-
-
-    where (drag_mx(:,:).ge.beta_limgz)
-       mstream_mx(:,:) = 0
-       betamx(:,:)     = beta_limgz
-       drag_mx(:,:)    = beta_limgz
-    elsewhere
-       mstream_mx(:,:) = 1
-       betamx(:,:)     = drag_mx(:,:)
-    endwhere
-
-    where (drag_my(:,:).ge.beta_limgz)
-       mstream_my(:,:) = 0
-       betamy(:,:)     = beta_limgz
-       drag_my(:,:)    = beta_limgz
-    elsewhere
-       mstream_my(:,:) = 1
-       betamy(:,:)     = drag_my(:,:)
-    endwhere
-
-    if (itracebug.eq.1)  call tracebug(' Prescr_beta  apres lecture')
-    mstream_mx(:,:) = 0
-    mstream_my(:,:) = 0
-
-    if (itracebug.eq.1)  call tracebug(' Prescr_beta       avant call gzm')
-    call gzm_beta_prescr
-!    call write_datfile2(nx,ny,betamx,betamy,'beta-LBq15-08.dat')
-
-    call make_gaussienne(n_g,sigma_beta,M_gauss)
-STOP
-    return
-
-  end subroutine init_dragging
-  !________________________________________________________________________________
-
-  !> SUBROUTINE: dragging 
-  !! Defini la localisation des streams et le frottement basal
-  !<
-  !-------------------------------------------------------------------------
-  subroutine dragging   ! defini la localisation des streams et le frottement basal
-
-    implicit none
-
-
-
-    if (itracebug.eq.1)  call tracebug(' beta_prescr         subroutine dragging')
-
-
-!   tentative d'iteration
-
-!    call beta_stag2c(nx,ny,drag_mx,drag_my,drag_centre)                  ! average  staggered 
-                                                                          ! on central 
-!   call iter_gauss(n_g,nx,ny,M_gauss,-1.,drag_centre,corrstream_centre,hdot)  ! correct central
-!   call beta_c2stag(nx,ny,drag_mx,drag_my,corrstream_centre)             ! centered distributed
-                                                                          ! on staggered
-!    drag_centre(:,:) = corrstream_centre(:,:)
-!    beta_centre(:,:) =  drag_centre(:,:)                                 ! surtout pour sorties
-
-    ! pour l'instant pas d'autre condition que mstream
-    where (mstream_mx(:,:).eq.1)
-       betamx(:,:) = drag_mx(:,:)
-    elsewhere 
-       betamx(:,:)= beta_limgz
-    end where
-
-    where (mstream_my(:,:).eq.1)
-       betamy(:,:) = drag_my(:,:)
-    elsewhere 
-       betamy(:,:)= beta_limgz
-    end where
-
-    betamx(:,:)=drag_mx(:,:)
-    betamy(:,:)=drag_my(:,:)
-
-    call  gzm_beta_prescr        ! determine gz et flgz et met a 0 le beta des points flottants
-
-
-    return
-  end subroutine dragging
-
-  !----------------------------------------------------------------------------------
-  !>SUBROUTINE: gzm_beta_prescr
-  !! Calcul de gzmx
-  !<
-
-  subroutine gzm_beta_prescr     ! attribue flgzmx et gzmx (et my)
-    
-    logical ::  test_Tx          ! test if there is a basal melting point in the surrounding
-    logical ::  test_Ty          ! test if there is a basal melting point in the surrounding
-    logical ::  test_beta_x      ! test if there is a low drag point in the surrounding
-    logical ::  test_beta_y      ! test if there is a low drag point in the surrounding
-
-    real    ::  beta_forc_fleuve ! below this value -> ice stream
-
-    beta_forc_fleuve = 5.e3
-
-    ! calcul de gzmx 
-
-
-    gzmx(:,:)=.true.
-    gzmy(:,:)=.true.
-    flgzmx(:,:)=.false.
-    flgzmy(:,:)=.false.
-
-
-    ! points flottants : flgzmx mais pas gzmx
-    do j=2,ny
-       do i=2,nx
-          if (flot(i,j).and.(flot(i-1,j))) then
-             flotmx(i,j)=.true.
-             flgzmx(i,j)=.true.
-             gzmx(i,j)=.false.
-             betamx(i,j)=0.
-          end if
-          if (flot(i,j).and.(flot(i,j-1))) then
-             flotmy(i,j)=.true.
-             flgzmy(i,j)=.true.
-             gzmy(i,j)=.false.
-             betamy(i,j)=0.
-          end if
-       end do
-    end do
-
-    !--------- autres criteres
-
-    ! points poses
-    gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.beta_limgz)   !  Pas de calcul pour les points
-    gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.beta_limgz)   !  au fort beta
-
-
-
-    ! critere (lache) sur la temperature
-    do j = 2, ny-1
-       do i =2, nx-1 
-
-!   test s'il y a un point tempere dans les environs
-          test_Tx = any( ibase( i-1:i   , j-1:j+1 ).gt.1)
-          test_Ty = any( ibase( i-1:i+1 , j-1:j   ).gt.1)
-
-!   test s'il y a un point low drag dans les environs
-          test_beta_x = any( drag_centre( i-1:i   , j-1:j+1 ) .lt. beta_forc_fleuve )
-          test_beta_y = any( drag_centre( i-1:i+1 , j-1:j   ) .lt. beta_forc_fleuve )
-
-!   critere pour gzmx
-          gzmx(i,j) =  gzmx(i,j) .and. (test_Tx.or.test_beta_x)
-          gzmy(i,j) =  gzmy(i,j) .and. (test_Ty.or.test_beta_y)
-
-       end do
-    end do
-
-    flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
-    flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:)
-
-    where ((.not.flotmx(:,:)).and.(.not.gzmx(:,:)))
-       betamx(:,:) = beta_limgz
-    end where
-
-    where ((.not.flotmy(:,:)).and.(.not.gzmy(:,:)))
-       betamy(:,:) = beta_limgz
-    end where
-
-
-    
-  end subroutine gzm_beta_prescr
-
-
-  !----------------------------------------------------------------------------------
-  ! Routines qui permettent des manip simples sur beta.
-
-
-!______________________________________________________________________________________
-!>SUBROUTINE:  beta_plateau_slope 
-!! give value beta_limgz to places where slope is >1.e-2 
-!<
-
-  subroutine beta_plateau_slope      ! give value beta_limgz to places where slope is >1.e-2
-
-    call slope_surf
-
-    where (abs(sdx(:,:)).gt.1.e-2)
-       drag_mx(:,:) = beta_limgz
-    end where
-
-    where (abs(sdy(:,:)).gt.1.e-2)
-       drag_my(:,:) = beta_limgz
-    end where
-  end subroutine beta_plateau_slope
-!______________________________________________________________________________________
-!> SUBROUTINE beta_plateau_S 
-!! give value beta_limgz to places where S > 2500 m
-!<
-  subroutine beta_plateau_S     !  give value beta_limgz to places where S > 2500 m
-
-    where (S(:,:).gt.2500.)     ! moyen simple de nettoyer l'interieur. WAIS pas touche
-       drag_centre(:,:)=beta_limgz
-       drag_mx(:,:)=beta_limgz
-       drag_my(:,:)=beta_limgz
-    end where
-  end subroutine beta_plateau_S
-!______________________________________________________________________________________
-!>SUBROUTINE:   beta_min_value 
-!! valeur mini que peut avoir beta (en Pa an m-1) 
-!<
-
-  subroutine beta_min_value(val_betamin)
-
-  real :: val_betamin     ! valeur mini que peut avoir beta (en Pa an m-1)
-
-    drag_centre(:,:) = max(drag_centre(:,:),val_betamin)
-    drag_mx(:,:) = max(drag_mx(:,:),val_betamin)
-    drag_my(:,:) = max(drag_my(:,:),val_betamin)
-    where(flot(:,:))
-       drag_centre(:,:) = 0.
-    end where
-
-  end subroutine beta_min_value
-!______________________________________________________________________________________
-!>SUBROUTINE:   beta_stag2c 
-!! average the staggered values on central ones
-!<
-
-  subroutine beta_stag2c(nxx,nyy,R_mx,R_my,R_centre)      ! average the staggered values on central ones
-
-    implicit none
-    integer                                ::  nxx,nyy                ! dimension des tableaux
-    real, dimension(nxx,nyy), intent(in)   ::  R_mx                   ! valeur du beta sur maille mx
-    real, dimension(nxx,nyy), intent(in)   ::  R_my                   ! valeur du beta sur maille my
-    real, dimension(nxx,nyy), intent(out)  ::  R_centre               ! valeur du beta sur maille centree   
-    
-    R_centre(:,:)=0.
-    do j=2,ny-1
-       do i=2,nx-1
-          R_centre(i,j)= ((R_mx(i,j)+R_mx(i+1,j)) &
-               + (R_my(i,j)+R_my(i,j+1)))*0.25
-       end do
-    end do
-  end subroutine beta_stag2c
-!______________________________________________________________________________________
-!>SUBROUTINE:   beta_c2stag_min  
-!! redistribute central values on staggered grid
-!<
-
-  subroutine beta_c2stag(nxx,nyy,R_mx,R_my,R_centre)              ! redistribute central values on staggered grid
-
-    implicit none
-    integer                                ::  nxx,nyy                ! dimension des tableaux
-    real, dimension(nxx,nyy), intent(out)  ::  R_mx                   ! valeur du beta sur maille mx
-    real, dimension(nxx,nyy), intent(out)  ::  R_my                   ! valeur du beta sur maille my
-    real, dimension(nxx,nyy), intent(in)   ::  R_centre               ! valeur du beta sur maille centree
-
-    do j=2,nyy 
-       do i=2,nxx
-          R_mx(i,j)= (R_centre(i-1,j)+R_centre(i,j))*0.5
-          R_my(i,j)= (R_centre(i,j-1)+R_centre(i,j))*0.5
-       end do
-    end do
-
-  end subroutine beta_c2stag
-!______________________________________________________________________________________
-!>SUBROUTINE:   beta_c2stag_min  
-!! redistribute central values on staggered grid
-!<
-
-  subroutine beta_c2stag_min(nxx,nyy,R_mx,R_my,R_centre)              ! redistribute central values on staggered grid
-
-    implicit none
-    integer                                ::  nxx,nyy                ! dimension des tableaux
-    real, dimension(nxx,nyy), intent(inout)::  R_mx                   ! valeur du beta sur maille mx
-    real, dimension(nxx,nyy), intent(inout)::  R_my                   ! valeur du beta sur maille my
-    real, dimension(nxx,nyy), intent(in)   ::  R_centre               ! valeur du beta sur maille centree
-
-    do j=2,nyy 
-       do i=2,nxx
-          R_mx(i,j)= min(R_mx(i,j),(R_centre(i-1,j)+R_centre(i,j))*0.5)
-          R_my(i,j)= min(R_my(i,j),(R_centre(i,j-1)+R_centre(i,j))*0.5)
-       end do
-    end do
-
-  end subroutine beta_c2stag_min
-!______________________________________________________________________________________
-!>SUBROUTINE:   make_gaussienne 
-!! cree une matrice gaussienne normalisee a 1
-!<
-
-  subroutine  make_gaussienne(n,sigma,M)
-
-    implicit none
-    integer,intent(in)                       :: n          ! demi taille du tableau gaussienne
-    real, intent(in)                         :: sigma      ! sigma gaussienne exprime en dx
-    real,dimension(-n:n,-n:n),intent(inout)  :: M          ! tableau gaussienne
-
-
-!   variables locales
-    real                                     :: dist2      ! distance au centre (au carre)
-    real                                     :: som        ! pour calibration
-    real                                     :: sigma2     ! pour gaussienne
-
-    som = 0.
-    sigma2 = sigma * sigma * 2
-    do j = -n, n
-       do i = -n, n
-          dist2 = (i*i+j*j)
-          M(i,j) = exp(-dist2/sigma2)
-          som=som+M(i,j)
-       end do
-    end do
-
-    M(:,:) = M(:,:) / som
-
-!    do j = -n, n
-!       write(6,'(5(f0.4,1x))') (M(i,j),i=-n,n)
-!    end do
-
-  end subroutine make_gaussienne
-
-!______________________________________________________________________________________
-!>SUBROUTINE:  iter_gauss
-!! utilise hdot pour faire les iterations sur beta_centre
-!<
-  subroutine iter_gauss(n,nxx,nyy,M,coef,beta_in,beta_out,Ecart)
-
-    integer,                  intent(in)    :: n          ! taille matrice gaussienne 
-    integer,                  intent(in)    :: nxx,nyy    ! taille du domaine
-    real,dimension(-n:n,-n:n),intent(in)    :: M          ! matrice gaussienne 
-    real,                     intent(in)    :: coef       ! coefficient multiplicateur
-    real,dimension(nx,ny),    intent(in)    :: beta_in    ! la valeur initiale
-    real,dimension(nx,ny),    intent(out)   :: beta_out   ! la valeur apres iteration
-    real,dimension(nx,ny),    intent(in)    :: Ecart      ! l'ecart a l'obs
-
-! variables locales
-    integer                                 :: ii,jj
-
-    do j= 1 + n, nyy - n
-       do i = 1 + n, nxx - n
-          beta_out(i,j) = beta_in(i,j)
-
-          do jj = -n, n
-             do ii = -n, n 
-                beta_out(i,j) = beta_out(i,j) + coef * M(ii,jj) * ecart(i+ii,j+jj)
-             end do
-          end do
-
-       end do
-    end do
-
-    beta_out(:,:) = max (beta_out(:,:),0.)
-  end subroutine iter_gauss
-
-end module dragging_prescr_beta
-
-!!$! add the nodes of topographical instability (on major nodes)
-! a garder eventuellement pour runs instabilite
-!!$beta_file = trim(dirnameinp)//'sensib_H_Gael-15km.dat'
-!!$call lect_datfile(nx,ny,drag_instab_topo,1,beta_file)   
-!!$
-!!$where (drag_instab_topo(:,:).gt.0.)
-!!$   drag_centre(:,:)=min(drag_centre(:,:),400.)
-!!$endwhere
Index: /trunk/SOURCES/Makefile.grisli.inc
===================================================================
--- /trunk/SOURCES/Makefile.grisli.inc	(revision 330)
+++ /trunk/SOURCES/Makefile.grisli.inc	(revision 331)
@@ -71,9 +71,6 @@
 	sliding_Bindshadler_mod.o no_sliding.o			    \
 	fake_beta_iter_vitbil_mod.o \
-	dragging_hwatermax_0.2_mod.o \
-	dragging_hwat-contigu_mod.o dragging_hwat_contmaj_mod.o \
-	dragging_hwat_sedim_mod.o dragging_neff_slope_mod.o \
 	dragging_param_beta_mod.o dragging_param_beta_sedim_mod.o \
-	dragging_param_beta_nolin_mod.o bilan_eau_mod.o \
+	bilan_eau_mod.o \
 	calving_frange.o calving_frange_abuk.o \
 	calving_frange_ISMIP_glaciers.o calving_frange_ISMIP_fracture.o \
@@ -94,8 +91,5 @@
   deformation_mod_2lois.o deformation_mod_2lois_isotherme.o      \
   sliding_Bindshadler_mod.o no_sliding.o                      \
-  dragging_hwatermax_0.2_mod.o dragging_calc_beta_mod.o  \
-  dragging_hwat-contigu_mod.o dragging_hwat_contmaj_mod.o \
-  dragging_hwat_sedim_mod.o dragging_neff_slope_mod.o \
-	dragging_param_beta_mod.o \
+  dragging_param_beta_mod.o dragging_param_beta_sedim_mod.o \
   calving_frange.o calving_frange_abuk.o calving_frange_ISMIP_fracture.o \
 	no_calving.o no_rsl.o imposed_regions_rsl.o lake_rsl_mod.o \
Index: /trunk/SOURCES/Old-sources/dragging_hwat-contigu_mod.f90
===================================================================
--- /trunk/SOURCES/Old-sources/dragging_hwat-contigu_mod.f90	(revision 331)
+++ /trunk/SOURCES/Old-sources/dragging_hwat-contigu_mod.f90	(revision 331)
@@ -0,0 +1,279 @@
+!> \file dragging_hwat-contigu_mod.f90
+!! Defini les zones de stream avec 2 criteres
+!<
+
+!> \namespace  dragging_hwat_cont
+!! Defini les zones de stream avec 2 criteres
+!! \author ...
+!! \date ... 
+!! @note Les 2 criteres sont: 
+!! @note   * un critere sur la hauteur d'eau
+!! @note   * un critere de continuite
+!! depuis la cote
+!! @note Used module
+!! @note   - use module3D_phy
+!<
+module dragging_hwat_cont
+
+! Defini les zones de stream avec un critere sur la hauteur d'eau + un critere de continuite
+! depuis la cote
+
+use module3d_phy
+use fake_beta_iter_vitbil_mod
+
+implicit none
+!logical,dimension(nx,ny) :: fleuvemx
+!logical,dimension(nx,ny) :: fleuvemy
+real,dimension(nx,ny) :: hwatmx
+real,dimension(nx,ny) :: hwatmY
+
+real :: valmax
+integer :: imax,jmax
+integer :: i_moins1,i_plus1,j_moins1,j_plus1
+integer :: lmax=20
+integer :: idep,jdep,iloc,jloc
+
+real :: tostick
+
+contains
+!-------------------------------------------------------------------------------
+!> SUBROUTINE: init_dragging
+!! Cette routine fait l'initialisation du dragging.
+!>
+subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
+
+implicit none
+ 
+namelist/drag_hwat_cont/hwatstream,cf,betamax,toblim,tostick    
+
+if (itracebug.eq.1)  call tracebug(' Dragging avec hwatermax')
+
+
+! formats pour les ecritures dans 42
+428 format(A)
+
+! lecture des parametres du run                      block draghwat
+!--------------------------------------------------------------------
+
+rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
+read(num_param,drag_hwat_cont)
+
+write(num_rep_42,428)'!___________________________________________________________' 
+write(num_rep_42,428) '&drag_hwat_cont              ! nom du bloc hwatermax '
+write(num_rep_42,*)
+write(num_rep_42,*) 'hwatsream      = ',hwatstream 
+write(num_rep_42,*) 'cf             = ',cf
+write(num_rep_42,*) 'betamax         = ', betamax
+write(num_rep_42,*) 'toblim         = ', toblim
+write(num_rep_42,*)'/'                            
+write(num_rep_42,428) '! hwatstream (m) :  critere de passage en stream en partant de la cote'
+write(num_rep_42,428) '!  si hwater > hwatstream '
+write(num_rep_42,428) '! cf coefficient de la loi de frottement fonction Neff'
+write(num_rep_42,428) '! betamax : (Pa) frottement maxi '
+write(num_rep_42,428) '! toblim : (Pa) tes pour les iles '
+write(num_rep_42,*)
+
+inv_beta=0
+tostick=1.e5   ! valeurs pour les points non flgzmx
+
+!-------------------------------------------------------------------
+! masque stream
+
+      mstream_mx(:,:)=1
+      mstream_my(:,:)=1
+
+
+! coefficient permettant de modifier le basal drag.
+      drag_mx(:,:)=1.
+      drag_my(:,:)=1.
+
+   
+      return
+      end subroutine init_dragging
+!________________________________________________________________________________
+
+!> SUBROUTINE: dragging
+!! Defini la localisation des streams et le frottement basal
+!>
+!-------------------------------------------------------------------------
+subroutine dragging   ! defini la localisation des streams et le frottement basal
+
+
+!         les iles n'ont pas de condition neff mais ont la condition toblim
+!         (ce bloc etait avant dans flottab)
+
+
+do j=2,ny
+   do i=2,nx
+      ilemx(i,j)=ilemx(i,j).and.(abs(tobmx(i,j)).lt.toblim)
+      ilemy(i,j)=ilemy(i,j).and.(abs(tobmy(i,j)).lt.toblim)
+   end do
+end do
+
+! calcul de la hauteur d'eau sur les mailles mineures
+do j=1,ny
+   do i=2,nx
+      hwatmx(i,j)=0.5*(hwater(i,j)+hwater(i-1,j))
+   end do
+end do
+
+do j=2,ny
+   do i=1,nx
+      hwatmy(i,j)=0.5*(hwater(i,j)+hwater(i,j-1))
+   end do
+end do
+
+
+debug_3D(:,:,1)=hwatmx(:,:)
+debug_3D(:,:,2)=hwatmy(:,:)
+
+
+! le gzmx initial a ete calculé dans flottab pour les points a cheval sur la grounding line
+
+! Partant de la cote, le fleuve est determine en remontant vers le point de hauteur d'eau max
+! tant que ce point est superieur a hwatermax.
+
+! Attention : ce systeme ne permet pas d'avoir des fleuves convergents
+!             et les fleuves n'ont une épaisseur que de 1 noeud.
+
+fleuvemx(:,:)=.false.
+fleuvemy(:,:)=.false.
+
+fleuve_mx: do j=1,ny  
+   do i=2,nx                      ! le noeud 1 n'est pas attribue
+
+cote_mx :  if (gzmx(i,j)) then    ! a ce niveau les gzmx sont sur la grounding line
+         idep=i
+         jdep=j
+
+suit_x : do l=1,lmax         ! debut de la boucle de suivi, lmax longueur maxi des fleuves
+           i_moins1=max(idep-1,2)
+           j_moins1=max(jdep-1,1)
+           i_plus1=min(idep+1,nx)
+           j_plus1=min(jdep+1,ny)
+           
+           valmax=0.
+
+! recherche du max en excluant les points flottants
+         do jloc=j_moins1,j_plus1
+            do iloc=i_moins1,i_plus1
+ 
+               if ((hwatmx(iloc,jloc).gt.valmax).and.(.not.flotmx(iloc,jloc))) then
+                  imax=iloc
+                  jmax=jloc
+                  valmax=hwatmx(iloc,jloc)
+               endif
+            end do
+         end do
+
+         if (hwatmx(imax,jmax).gt.hwatstream) then
+            fleuvemx(imax,jmax)=.true.
+            idep=imax
+            jdep=jmax
+         end if
+      end do suit_x
+
+   end if cote_mx
+
+if ((.not.ilemx(i,j)).and.(fleuvemx(i,j))) gzmx(i,j)=.true. 
+
+
+
+! calcul du frottement basal (ce bloc etait avant dans neffect)
+         if (gzmx(i,j)) then                          ! stream
+            betamx(i,j)=cf*neffmx(i,j) 
+         else if (flotmx(i,j).or.ilemx(i,j)) then     ! flottant ou ile
+            betamx(i,j)=0.
+         else                                         ! grounded, SIA
+            betamx(i,j)=tostick                       ! frottement glace posee (1 bar)
+         endif      
+
+ end do
+end do fleuve_mx
+
+
+
+
+fleuve_my: do j=2,ny                  ! le noeud 1 n'est pas attribue
+   do i=1,nx   
+
+cote_my :  if (gzmy(i,j)) then    ! a ce niveau les gzmx sont sur la grounding line
+         idep=i
+         jdep=j
+
+suit_y : do l=1,lmax         ! debut de la boucle de suivi, lmax longueur maxi des fleuves
+           i_moins1=max(idep-1,1)
+           j_moins1=max(jdep-1,2)
+           i_plus1=min(idep+1,nx)
+           j_plus1=min(jdep+1,ny)
+           
+           valmax=0.
+
+! recherche du max en excluant les points flottants
+         do jloc=j_moins1,j_plus1
+            do iloc=i_moins1,i_plus1
+ 
+               if ((hwatmy(iloc,jloc).gt.valmax).and.(.not.flotmy(iloc,jloc))) then
+                  imax=iloc
+                  jmax=jloc
+                  valmax=hwatmy(iloc,jloc)
+               endif
+            end do
+         end do
+
+         if (hwatmy(imax,jmax).gt.hwatstream) then
+            fleuvemy(imax,jmax)=.true.
+            idep=imax
+            jdep=jmax
+         end if
+      end do suit_y
+
+   end if cote_my
+
+if ((.not.ilemy(i,j)).and.(fleuvemy(i,j))) gzmy(i,j)=.true. 
+
+
+ ! calcul du frottement basal (ce bloc etait avant dans neffect)
+   if (gzmy(i,j)) then                ! stream
+      betamy(i,j)=cf*neffmy(i,j) 
+   else if (flotmy(i,j).or.ilemy(i,j)) then     ! flottant ou ile
+      betamy(i,j)=0.
+   else                               ! grounded, SIA
+      betamy(i,j)=tostick             ! frottement glace posee
+   endif
+
+
+ end do
+end do fleuve_my
+
+
+
+where (fleuvemx(:,:))
+   debug_3D(:,:,3)=1
+elsewhere
+   debug_3D(:,:,3)=0
+endwhere
+
+where (flotmx(:,:))
+   debug_3D(:,:,3)=-1
+endwhere
+
+!_____________________
+
+
+where (fleuvemy(:,:))
+   debug_3D(:,:,4)=1
+elsewhere
+   debug_3D(:,:,4)=0
+endwhere
+
+where (flotmy(:,:))
+   debug_3D(:,:,4)=-1
+end where
+
+
+return
+end subroutine dragging
+
+end module dragging_hwat_cont
+
Index: /trunk/SOURCES/Old-sources/dragging_hwat_contmaj_mod.f90
===================================================================
--- /trunk/SOURCES/Old-sources/dragging_hwat_contmaj_mod.f90	(revision 331)
+++ /trunk/SOURCES/Old-sources/dragging_hwat_contmaj_mod.f90	(revision 331)
@@ -0,0 +1,391 @@
+!> \file dragging_hwat_contmaj_mod.f90
+!! Defini les zones de stream avec 3 criteres 
+!<
+
+!> \namespace  dragging_hwat_contmaj
+!! Defini les zones de stream avec 3 criteres
+!! \author ...
+!! \date ...
+!! @note Les trois criteres sont:
+!! @note   * un critere sur la hauteur d'eau 
+!! @note   * un critere de continuite depuis la cote
+!! @note   * un critere sur la courbure du socle (si negatif, vallees) 
+!! @note Used module
+!! @note   - use module3D_phy
+!! @note   - use param_phy_mod
+!<
+
+module dragging_hwat_contmaj
+
+! Defini les zones de stream avec :
+! * un critere sur la hauteur d'eau 
+! * un critere de continuite depuis la cote
+! * un critere sur la courbure du socle (si negatif, vallees) 
+
+use module3d_phy
+use param_phy_mod
+use fake_beta_iter_vitbil_mod
+
+implicit none
+!logical,dimension(nx,ny) :: fleuvemx
+!logical,dimension(nx,ny) :: fleuvemy
+logical,dimension(nx,ny) :: fleuve
+logical,dimension(nx,ny) :: cote
+
+
+real :: valmax
+integer :: imax,jmax
+integer :: i_moins1,i_plus1,j_moins1,j_plus1
+integer :: lmax=20
+integer :: idep,jdep,iloc,jloc
+
+
+real :: tostick   ! pour la glace posee
+real :: tob_ile    ! pour les iles
+real :: cry_lim=50.  ! courbure limite pour le suivi des fleuves
+
+
+real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite avec spinup cat
+real, dimension(nx,ny) :: Vcol_y           !< uniquement pour compatibilite avec spinup cat
+real, dimension(nx,ny) :: Vsl_x            !< uniquement pour compatibilite avec spinup cat
+real, dimension(nx,ny) :: Vsl_y            !< uniquement pour compatibilite avec spinup cat
+logical :: corr_def = .false.              !< for deformation correction, pour compatibilite beta
+
+contains
+!-------------------------------------------------------------------------------
+
+!> SUBROUTINE: init_dragging
+!! Cette routine fait l'initialisation du dragging.
+!>
+subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
+
+implicit none
+ 
+namelist/drag_hwat_cont/hwatstream,cf,betamax,toblim,tostick    
+
+if (itracebug.eq.1)  call tracebug(' dragging avec hwatermax')
+
+
+! formats pour les ecritures dans 42
+428 format(A)
+
+! lecture des parametres du run                      block draghwat
+!--------------------------------------------------------------------
+
+rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
+read(num_param,drag_hwat_cont)
+
+write(num_rep_42,428)'!___________________________________________________________' 
+write(num_rep_42,428) '&drag_hwat_cont              ! nom du bloc hwatermax '
+write(num_rep_42,*)
+write(num_rep_42,*) 'hwatstream      = ',hwatstream 
+write(num_rep_42,*) 'cf             = ',cf
+write(num_rep_42,*) 'betamax         = ', betamax
+write(num_rep_42,*) 'toblim         = ', toblim
+write(num_rep_42,*)'/'                            
+write(num_rep_42,428) '! hwatstream (m) :  critere de passage en stream en partant de la cote'
+write(num_rep_42,428) '!  si hwater > hwatstream '
+write(num_rep_42,428) '! cf coefficient de la loi de frottement fonction Neff'
+write(num_rep_42,428) '! seulement pour les points cotiers'
+write(num_rep_42,428) '! betamax : (Pa) frottement maxi sous les streams '
+write(num_rep_42,428) '! toblim : (Pa)  pour les iles '
+write(num_rep_42,*)
+
+inv_beta=0
+tostick=1.e5   ! valeurs pour les points non flgzmx
+tob_ile=betamax/2.
+moteurmax=toblim
+betamax_2d(:,:)=betamax
+
+!-------------------------------------------------------------------
+! masque stream
+
+      mstream_mx(:,:)=1
+      mstream_my(:,:)=1
+
+
+! coefficient permettant de modifier le basal drag.
+      drag_mx(:,:)=1.
+      drag_my(:,:)=1.
+
+where (socle_cry(:,:).lt.cry_lim)
+   debug_3D(:,:,25)=1
+elsewhere
+   debug_3D(:,:,25)=0
+endwhere
+
+return
+end subroutine init_dragging
+!________________________________________________________________________________
+
+!> SUBROUTINE: dragging
+!! Defini la localisation des streams et le frottement basal
+!>
+!-------------------------------------------------------------------------
+subroutine dragging   ! defini la localisation des streams et le frottement basal
+
+
+!         les iles n'ont pas de condition neff mais ont la condition toblim
+!         (ce bloc etait avant dans flottab)
+
+
+do j=2,ny
+   do i=2,nx
+      ilemx(i,j)=ilemx(i,j).and.(abs(rog*Hmx(i,j)*sdx(i,j)).lt.toblim)
+      ilemy(i,j)=ilemy(i,j).and.(abs(rog*Hmy(i,j)*sdy(i,j)).lt.toblim)
+   end do
+end do
+
+! le calcul des fleuves se fait sur les mailles majeures en partant de la cote
+
+
+
+
+! le gzmx initial a ete calculé dans flottab pour les points a cheval sur la grounding line
+
+! Partant de la cote, le fleuve est determine en remontant vers le point de hauteur d'eau max
+! tant que ce point est superieur a hwatermax.
+
+! Attention : ce systeme ne permet pas d'avoir des fleuves convergents
+!             et les fleuves n'ont une épaisseur que de 1 noeud.
+
+fleuvemx(:,:)=.false.
+fleuvemy(:,:)=.false.
+fleuve(:,:)=.false.
+cote(:,:)=.false.
+
+! detection des cotes
+do  j=2,ny-1  
+   do i=2,nx-1  
+      if ((.not.flot(i,j)).and. & 
+      ((flot(i+1,j)).or.(flot(i,j+1)).or.(flot(i-1,j)).or.(flot(i,j-1)))) then
+         cote(i,j)=.true.
+      endif
+   end do
+end do
+
+fleuve_maj: do j=2,ny-1  
+ifleuve:   do i=2,nx-1                     
+
+cote_detect :  if (cote(i,j)) then
+         idep=i
+         jdep=j
+
+!         if (hwater(i,j).gt.hwatstream) then
+         if (socle_cry(i,j).lt.0.) then    ! dans une vallee
+            fleuve(i,j)=.true. 
+!           write(166,*)
+!           write(166,*)'cote',i,j,hwater(i,j),hwatstream
+         else
+            cote(i,j)=.false.
+            cycle ifleuve
+         endif
+
+suit : do l=1,lmax         ! debut de la boucle de suivi, lmax longueur maxi des fleuves
+           i_moins1=max(idep-1,2)
+           j_moins1=max(jdep-1,1)
+           i_plus1=min(idep+1,nx)
+           j_plus1=min(jdep+1,ny)
+           
+
+
+!!$! recherche du max en suivant la hauteur d'eau maxi chez les voisins
+!!$! * en excluant les points flottants
+!!$! * et ceux qui sont deja tagges fleuves   !version 18 juillet 2007
+
+!!$           valmax=0.
+!!$
+!!$         do jloc=j_moins1,j_plus1
+!!$            do iloc=i_moins1,i_plus1
+!!$ 
+!!$               if ((hwater(iloc,jloc).gt.valmax) & 
+!!$                    .and.(.not.flot(iloc,jloc))  &
+!!$                    .and.(.not.fleuve(iloc,jloc))) then
+!!$                  imax=iloc
+!!$                  jmax=jloc
+!!$                  valmax=hwater(iloc,jloc)
+!!$               endif
+!!$            end do
+!!$         end do
+
+! recherche du max en suivant le socle le plus profond
+! * en excluant les points flottants
+! * et ceux qui sont deja tagges fleuves
+
+           valmax=1000.
+
+        do jloc=j_moins1,j_plus1
+            do iloc=i_moins1,i_plus1
+ 
+               if ((B(iloc,jloc).lt.valmax)      & 
+                    .and.(.not.flot(iloc,jloc))  &
+                    .and.(.not.fleuve(iloc,jloc)).and.(socle_cry(iloc,jloc).lt.cry_lim)) then
+                  imax=iloc
+                  jmax=jloc
+                  valmax=B(iloc,jloc)
+               endif
+            end do
+         end do
+
+
+
+!         if (hwater(imax,jmax).gt.max(hwatstream-L*40.,200.)) then  ! version 18 juil
+         if ((hwater(imax,jmax).gt.hwatstream).and.(socle_cry(i,j).lt.cry_lim)) then
+            fleuve(imax,jmax)=.true.
+            idep=imax
+            jdep=jmax
+!            write(166,*) imax,jmax,hwater(imax,jmax)
+         else
+            fleuve(imax,jmax)=.false.
+            exit suit
+         end if
+ 
+      end do suit
+
+   end if cote_detect
+
+end do ifleuve
+end do fleuve_maj
+
+
+! detection de problemes
+!write(166,*)
+!write(166,*) 'detection de pb dans fleuve   time=',time
+
+do j=2,ny
+   do i=2,nx
+      if (fleuve(i,j).and.(socle_cry(i,j).gt.cry_lim)) then
+         write(166,*) 'pb fleuve',time,i,j,hwater(i,j),socle_cry(i,j)
+      endif
+   end do
+end do
+
+do j=1,ny-1
+   do i=1,nx-1
+      if (fleuve(i,j)) then
+         fleuvemx(i,j)=.true.
+         fleuvemx(i+1,j)=.true.
+         fleuvemy(i,j)=.true.
+         fleuvemy(i,j+1)=.true.
+      end if
+   end do
+end do
+
+! pas de fleuve dans les endroits trop en pente
+
+fleuvemx(:,:)=fleuvemx(:,:).and.(abs(rog*Hmx(:,:)*sdx(:,:)).lt.toblim)
+fleuvemy(:,:)=fleuvemy(:,:).and.(abs(rog*Hmy(:,:)*sdy(:,:)).lt.toblim)
+
+do j=1,ny
+   do i=1,nx
+      if ((.not.ilemx(i,j)).and.(fleuvemx(i,j))) gzmx(i,j)=.true. 
+
+! calcul du frottement basal (ce bloc etait avant dans neffect)
+
+      if (cotemx(i,j)) then                        ! point cotier
+         betamx(i,j)=cf*neffmx(i,j) 
+         betamx(i,j)=min(betamx(i,j),betamax)
+      else if ((gzmx(i,j)).and.(.not.cotemx(i,j))) then  ! stream
+         betamx(i,j)=betamax
+         betamx(i,j)=cf*neffmx(i,j)*20.                  ! attention a cette modif juillet 2007
+         betamx(i,j)=min(betamx(i,j),betamax)
+         betamx(i,j)=max(betamx(i,j),10.)
+
+         if (cf*neffmx(i,j).gt.1500.) then
+            gzmx(i,j)=.false.
+            betamx(i,j)=tostick
+         endif
+
+
+      else if (ilemx(i,j)) then
+         betamx(i,j)=cf*neffmx(i,j) 
+         betamx(i,j)=min(betamx(i,j),tob_ile)       
+      else if (flotmx(i,j)) then     ! flottant ou ile
+         betamx(i,j)=0.
+      else                                         ! grounded, SIA
+         betamx(i,j)=tostick                       ! frottement glace posee (1 bar)
+      endif
+
+   end do
+end do
+
+
+do j=1,ny                  ! le noeud 1 n'est pas attribue
+   do i=1,nx   
+
+      if ((.not.ilemy(i,j)).and.(fleuvemy(i,j))) gzmy(i,j)=.true. 
+
+
+ ! calcul du frottement basal (ce bloc etait avant dans neffect)
+
+      if (cotemy(i,j)) then                        ! point cotier
+         betamy(i,j)=cf*neffmy(i,j) 
+         betamy(i,j)=min(betamy(i,j),betamax)
+      else if ((gzmy(i,j)).and.(.not.cotemy(i,j))) then  ! stream
+         betamy(i,j)=betamax
+         betamy(i,j)=cf*neffmy(i,j)*20. 
+         betamy(i,j)=min(betamy(i,j),betamax)
+         betamy(i,j)=max(betamy(i,j),10.)
+
+         if (cf*neffmy(i,j).gt.1500.) then
+            gzmy(i,j)=.false.
+            betamy(i,j)=tostick
+         endif
+
+      else if (ilemy(i,j)) then
+         betamy(i,j)=cf*neffmy(i,j) 
+         betamy(i,j)=min(betamy(i,j),tob_ile)     
+      else if (flotmy(i,j)) then     ! flottant ou ile
+         betamy(i,j)=0.
+      else                               ! grounded, SIA
+         betamy(i,j)=tostick             ! frottement glace posee
+      endif
+   end do
+end do
+
+
+where (fleuve(:,:))
+   debug_3D(:,:,1)=1
+elsewhere (flot(:,:))
+   debug_3D(:,:,1)=2
+elsewhere
+   debug_3D(:,:,1)=0
+endwhere
+
+where (cote(:,:))
+   debug_3D(:,:,2)=1
+elsewhere
+   debug_3D(:,:,2)=0
+endwhere
+
+where (fleuvemx(:,:))
+   debug_3D(:,:,3)=1
+elsewhere
+   debug_3D(:,:,3)=0
+endwhere
+
+where (flotmx(:,:))
+   debug_3D(:,:,3)=-1
+endwhere
+
+!_____________________
+
+
+where (fleuvemy(:,:))
+   debug_3D(:,:,4)=1
+elsewhere
+   debug_3D(:,:,4)=0
+endwhere
+
+where (flotmy(:,:))
+   debug_3D(:,:,4)=-1
+end where
+
+debug_3D(:,:,23)= abs(RO*G*HMX(:,:)*sdx(:,:)*1.e-5) 
+debug_3D(:,:,24)= abs(RO*G*HMY(:,:)*sdy(:,:)*1.e-5) 
+
+return
+end subroutine dragging
+
+end module dragging_hwat_contmaj
+
Index: /trunk/SOURCES/Old-sources/dragging_hwat_sedim_mod.f90
===================================================================
--- /trunk/SOURCES/Old-sources/dragging_hwat_sedim_mod.f90	(revision 331)
+++ /trunk/SOURCES/Old-sources/dragging_hwat_sedim_mod.f90	(revision 331)
@@ -0,0 +1,477 @@
+module dragging_hwat_sedim
+
+! Defini les zones de stream avec :
+! * un critere sur la hauteur d'eau 
+! * un critere sur la courbure du socle (si negatif, vallees)
+! * un critere sur l'epaisseur de sediments (carte de Gabi & Laske) 
+
+use module3d_phy
+use fake_beta_iter_vitbil_mod
+
+implicit none
+!logical,dimension(nx,ny) :: fleuvemx
+!logical,dimension(nx,ny) :: fleuvemy
+logical,dimension(nx,ny) :: fleuve
+logical,dimension(nx,ny) :: fleuve_sedim
+logical,dimension(nx,ny) :: fleuve_sedimx
+logical,dimension(nx,ny) :: fleuve_sedimy
+logical,dimension(nx,ny) :: cote
+
+
+real :: valmax
+integer :: imax,jmax
+integer :: i_moins1,i_plus1,j_moins1,j_plus1
+integer :: lmax=20
+integer :: idep,jdep,iloc,jloc
+
+
+real :: tostick   ! pour la glace posee
+real :: tob_ile    ! pour les iles
+real :: cry_lim=50.  ! courbure limite pour le suivi des fleuves
+
+real,dimension(nx,ny) :: h_sedim    ! epaisseur de sediment
+real,dimension(nx,ny) :: neff ! pression effective noeuds majeurs
+real :: seuil_hwater  ! seuil sur hwater pour avoir du glissement en zone sediment
+real :: seuil_sedim   ! seuil sur h_sedim pour avoir du glissement en zone sediment
+real :: seuil_neff    ! seuil sur la pression effective pour avoir glissement en zone sediment
+real :: coef_sedim    ! coef frottement zones sediments  (1)
+real :: coef_gz       ! coef frottement zones stream std  (10)
+real :: coef_ile      ! coef frottement zones iles   (0.1)
+real :: toto
+
+real :: betacotemax=100  ! valeur max pour betamx et betamy sur les cote
+real :: betacotemin=4
+real :: betasedmax=400    ! 400 dans c-HadCM3-evol-basal-melt-sed : cHadCM3s  et 250 dans c-HadCM3-evol-basal-melt-sed-betasedimin6
+real :: betasedmin=10   ! 10 dans c-HadCM3-evol-basal-melt-sed : cHadCM3s et 6 dans c-HadCM3-evol-basal-melt-sed-betasedimin6
+real :: betastreamax=500
+real :: betastreamin=10
+real :: betailemax=250
+real :: betailemin=6
+
+real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite avec spinup cat
+real, dimension(nx,ny) :: Vcol_y           !< uniquement pour compatibilite avec spinup cat
+real, dimension(nx,ny) :: Vsl_x            !< uniquement pour compatibilite avec spinup cat
+real, dimension(nx,ny) :: Vsl_y            !< uniquement pour compatibilite avec spinup cat
+logical :: corr_def = .false.              !< for deformation correction, pour compatibilite beta
+
+character(len=150) :: file_sedim   ! fichier epaisseur sediment
+
+contains
+!-------------------------------------------------------------------------------
+subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
+
+implicit none
+ 
+namelist/drag_hwat_sedim/file_sedim,hwatstream,cf,betamax,toblim,seuil_hwater,seuil_sedim,seuil_neff,coef_sedim,coef_gz,coef_ile    
+
+if (itracebug.eq.1) write(num_tracebug,*)' dragging avec hwatermax'
+
+inv_beta=0
+
+! formats pour les ecritures dans 42
+428 format(A)
+
+! lecture des parametres du run                      block draghwat
+!--------------------------------------------------------------------
+
+rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
+read(num_param,drag_hwat_sedim)
+
+write(num_rep_42,428)'!___________________________________________________________' 
+write(num_rep_42,428) '&drag_hwat_sedim              ! nom du bloc hwat_sedim '
+write(num_rep_42,*)
+write(num_rep_42,*) 'file_sedim     = ', trim(file_sedim)
+write(num_rep_42,*) 'hwatstream     = ', hwatstream 
+write(num_rep_42,*) 'cf             = ', cf
+write(num_rep_42,*) 'betamax        = ', betamax
+write(num_rep_42,*) 'toblim         = ', toblim
+write(num_rep_42,*) 'seuil_hwater   = ', seuil_hwater
+write(num_rep_42,*) 'seuil_sedim    = ', seuil_sedim
+write(num_rep_42,*) 'seuil_neff     = ', seuil_neff
+write(num_rep_42,*) 'coef_sedim     = ', coef_sedim
+write(num_rep_42,*) 'coef_gz        = ', coef_gz
+write(num_rep_42,*) 'coef_ile       = ', coef_ile
+
+write(num_rep_42,*) '/'
+write(num_rep_42,428) '! file_sedim : fichier carte epaisseur sediment'                            
+write(num_rep_42,428) '! hwatstream (m) :  critere de passage en stream en partant de la cote si hwater > hwatstream '
+write(num_rep_42,428) '! cf coefficient de la loi de frottement fonction Neff'
+write(num_rep_42,428) '! seulement pour les points cotiers'
+write(num_rep_42,428) '! betamax : (Pa) frottement maxi sous les streams '
+write(num_rep_42,428) '! toblim : (Pa)  pour les iles '
+write(num_rep_42,428) '! seuil_hwater (m) : seuil hwater pour avoir glissement sur zone sediment'
+write(num_rep_42,428) '! seuil_sedim (m) : seuil epaisseur sediment pour avoir glissement sur zone sediment'
+write(num_rep_42,428) '! seuil_neff (Pa) seuil sur la pression effective pour avoir glissement en zone sediment'
+write(num_rep_42,428) '! coef_sedim : coef frottement zones stream sediments'
+write(num_rep_42,428) '! coef_gz : coef frottement zones stream std'
+write(num_rep_42,428) '! coef_ile : coef frottement zones iles'
+write(num_rep_42,*)
+
+tostick=1.e5   ! valeurs pour les points non flgzmx
+tob_ile=betamax/2.
+moteurmax=toblim
+
+!-------------------------------------------------------------------
+! masque stream
+
+      mstream_mx(:,:)=1
+      mstream_my(:,:)=1
+
+
+! coefficient permettant de modifier le basal drag.
+      drag_mx(:,:)=1.
+      drag_my(:,:)=1.
+
+where (socle_cry(:,:).lt.cry_lim)
+   debug_3D(:,:,25)=1
+elsewhere
+   debug_3D(:,:,25)=0
+endwhere
+
+! lecture du fichier sediment :
+!open(20,file=TRIM(DIRNAMEINP)//trim(file_sedim)) !   'sediment_ij_hemin40.dat')
+!do j=1,ny
+!   do i=1,nx
+!      read(20,*),toto,toto,h_sedim(i,j)
+!   enddo
+!enddo
+!h_sedim(:,:)=max(h_sedim(:,:),0.)  ! pas d'epaisseurs de sediments negatives
+h_sedim(:,:)=10. 
+  
+      return
+      end subroutine init_dragging
+!________________________________________________________________________________
+
+
+!-------------------------------------------------------------------------
+subroutine dragging   ! defini la localisation des streams et le frottement basal
+
+
+!         les iles n'ont pas de condition neff mais ont la condition toblim
+!         (ce bloc etait avant dans flottab)
+
+
+do j=2,ny
+   do i=2,nx
+      ilemx(i,j)=ilemx(i,j).and.(abs(rog*Hmx(i,j)*sdx(i,j)).lt.toblim)
+      ilemy(i,j)=ilemy(i,j).and.(abs(rog*Hmy(i,j)*sdy(i,j)).lt.toblim)
+   end do
+end do
+
+! le calcul des fleuves se fait sur les mailles majeures en partant de la cote
+
+
+
+
+! le gzmx initial a ete calculé dans flottab pour les points a cheval sur la grounding line
+
+! Partant de la cote, le fleuve est determine en remontant vers le point de hauteur d'eau max
+! tant que ce point est superieur a hwatermax.
+
+! Attention : ce systeme ne permet pas d'avoir des fleuves convergents
+!             et les fleuves n'ont une épaisseur que de 1 noeud.
+
+fleuvemx(:,:)=.false.
+fleuvemy(:,:)=.false.
+fleuve(:,:)=.false.
+cote(:,:)=.false.
+
+! detection des cotes
+do  j=2,ny-1  
+   do i=2,nx-1  
+      if ((.not.flot(i,j)).and. & 
+      ((flot(i+1,j)).or.(flot(i,j+1)).or.(flot(i-1,j)).or.(flot(i,j-1)))) then
+         cote(i,j)=.true.
+      endif
+   end do
+end do
+
+! calcul de neff (pression effective sur noeuds majeurs)
+do j=1,ny-1
+   do i=1,nx-1
+        neff(i,j)=(neffmx(i,j)+neffmx(i,j+1)+neffmy(i,j)+neffmy(i+1,j))/4
+   enddo
+enddo
+
+! test sur hauteur d'eau et epaisseur sediment :
+!where (h_sedim(:,:).GT.seuil_sedim.and.hwater(:,:).GT.seuil_hwater.and.(.not.flot(:,:))) fleuve(:,:)=.true.
+where (h_sedim(:,:).GT.seuil_sedim.and.neff(:,:).LE.seuil_neff.and.hwater(:,:).GT.seuil_hwater.and.(.not.flot(:,:))) fleuve(:,:)=.true.
+
+fleuve_sedim(:,:)=fleuve(:,:)  ! points uniquement fleuves sediment
+
+
+fleuve_maj: do j=2,ny-1  
+ifleuve:   do i=2,nx-1                     
+
+cote_detect :  if (cote(i,j)) then
+   idep=i
+   jdep=j
+         
+      if (socle_cry(i,j).lt.0.) then    ! dans une vallee
+         fleuve(i,j)=.true. 
+!         write(166,*)
+!         write(166,*)'cote',i,j,hwater(i,j),hwatstream
+      else
+         cote(i,j)=.false.
+         cycle ifleuve
+      endif
+
+suit : do l=1,lmax         ! debut de la boucle de suivi, lmax longueur maxi des fleuves
+           i_moins1=max(idep-1,2)
+           j_moins1=max(jdep-1,1)
+           i_plus1=min(idep+1,nx)
+           j_plus1=min(jdep+1,ny)
+           
+
+
+!!$! recherche du max en suivant la hauteur d'eau maxi chez les voisins
+!!$! * en excluant les points flottants
+!!$! * et ceux qui sont deja tagges fleuves   !version 18 juillet 2007
+
+!!$           valmax=0.
+!!$
+!!$         do jloc=j_moins1,j_plus1
+!!$            do iloc=i_moins1,i_plus1
+!!$ 
+!!$               if ((hwater(iloc,jloc).gt.valmax) & 
+!!$                    .and.(.not.flot(iloc,jloc))  &
+!!$                    .and.(.not.fleuve(iloc,jloc))) then
+!!$                  imax=iloc
+!!$                  jmax=jloc
+!!$                  valmax=hwater(iloc,jloc)
+!!$               endif
+!!$            end do
+!!$         end do
+
+! recherche du max en suivant le socle le plus profond
+! * en excluant les points flottants
+! * et ceux qui sont deja tagges fleuves
+
+           valmax=1000.
+
+        do jloc=j_moins1,j_plus1
+            do iloc=i_moins1,i_plus1
+ 
+               if ((B(iloc,jloc).lt.valmax)      & 
+                    .and.(.not.flot(iloc,jloc))  &
+                    .and.(.not.fleuve(iloc,jloc)).and.(socle_cry(iloc,jloc).lt.cry_lim)) then
+                  imax=iloc
+                  jmax=jloc
+                 valmax=B(iloc,jloc)
+               endif
+            end do
+         end do
+
+
+
+
+         if ((hwater(imax,jmax).gt.hwatstream).and.(socle_cry(i,j).lt.cry_lim)) then
+            fleuve(imax,jmax)=.true.
+            idep=imax
+            jdep=jmax
+!            write(166,*) imax,jmax,hwater(imax,jmax)
+         else
+!            fleuve(imax,jmax)=.false.
+            exit suit
+         end if
+ 
+      end do suit
+
+   end if cote_detect
+
+end do ifleuve
+end do fleuve_maj
+
+
+! detection de problemes
+!write(166,*)
+!write(166,*) 'detection de pb dans fleuve   time=',time
+!do j=2,ny
+!   do i=2,ny
+!      if (fleuve(i,j).and.(socle_cry(i,j).gt.cry_lim)) then
+!         write(166,*) 'pb fleuve',time,i,j,hwater(i,j),socle_cry(i,j)
+!      endif
+!   end do
+!end do
+
+do j=1,ny-1
+   do i=1,nx-1
+      if (fleuve(i,j)) then
+         fleuvemx(i,j)=.true.
+         fleuvemx(i+1,j)=.true.
+         fleuvemy(i,j)=.true.
+         fleuvemy(i,j+1)=.true.
+      end if
+      if (fleuve_sedim(i,j)) then
+         fleuve_sedimx(i,j)=.true.
+         fleuve_sedimx(i+1,j)=.true.
+         fleuve_sedimy(i,j)=.true.
+         fleuve_sedimy(i,j+1)=.true.
+      endif
+   end do
+end do
+
+! pas de fleuve dans les endroits trop en pente
+!fleuvemx(:,:)=fleuvemx(:,:).and.(abs(rog*Hmx(:,:)*sdx(:,:)).lt.toblim)
+!fleuvemy(:,:)=fleuvemy(:,:).and.(abs(rog*Hmy(:,:)*sdy(:,:)).lt.toblim)
+!fleuvemx(:,:)=fleuvemx(:,:).and.(abs(rog*Hmx(:,:)*sdx(:,:)).lt.toblim)
+!fleuvemy(:,:)=fleuvemy(:,:).and.(abs(rog*Hmy(:,:)*sdy(:,:)).lt.toblim)
+
+
+do j=1,ny
+   do i=1,nx
+      if ((.not.ilemx(i,j)).and.(fleuvemx(i,j))) gzmx(i,j)=.true. 
+
+! calcul du frottement basal (ce bloc etait avant dans neffect)
+
+      if (cotemx(i,j)) then                        ! point cotier
+         betamx(i,j)=cf*neffmx(i,j) 
+! cat         betamx(i,j)=min(betamx(i,j),tobmax)
+! version Jorge
+! point cotier: frott maxi=40, mini=4
+!         betamx(i,j)=min(betamx(i,j),betacotemax)
+!         betamx(i,j)=max(betamx(i,j),betacotemin)
+! version Jorge stream  en zone sediment :
+      else if (fleuve_sedimx(i,j).and.gzmx(i,j).and.(.not.cotemx(i,j))) then  ! stream  en zone sediment
+! stream sediment: frott maxi=60, mini=6 
+!cdc            betamx(i,j)=cf*neffmx(i,j)*3
+             betamx(i,j)=cf*neffmx(i,j)*coef_sedim
+!            betamx(i,j)=min(betamx(i,j),betasedmax)          
+!            betamx(i,j)=max(betamx(i,j),betasedmin)
+      else if ((gzmx(i,j)).and.(.not.cotemx(i,j))) then  ! stream normal
+! cat         betamx(i,j)=tobmax
+! cat         betamx(i,j)=cf*neffmx(i,j)*20.                  ! attention a cette modif juillet 2007
+! cat         betamx(i,j)=min(betamx(i,j),tobmax)
+! cat         betamx(i,j)=max(betamx(i,j),10.)
+! version Jorge
+! stream normal: frott maxi=80, mini=8
+         betamx(i,j)=cf*neffmx(i,j)*coef_gz         
+!         betamx(i,j)=min(betamx(i,j),betastreamax)           
+!         betamx(i,j)=max(betamx(i,j),betastreamin)
+
+! cat         if (cf*neffmx(i,j).gt.1500.) then
+! cat          gzmx(i,j)=.false.
+! cat           betamx(i,j)=tostick
+! cat        endif
+
+      else if (ilemx(i,j)) then
+         betamx(i,j)=cf*neffmx(i,j)*coef_ile
+! cat         betamx(i,j)=min(betamx(i,j),tob_ile)    
+! version Jorge
+! ile : frott maxi=20, mini=2   
+!         betamx(i,j)=min(betamx(i,j),betailemax)
+!         betamx(i,j)=max(betamx(i,j),betailemin)
+      else if (flotmx(i,j)) then     ! flottant ou ile
+         betamx(i,j)=0.
+      else                                         ! grounded, SIA
+         betamx(i,j)=tostick                       ! frottement glace posee (1 bar)
+      endif
+   end do
+end do
+
+
+do j=1,ny                  ! le noeud 1 n'est pas attribue
+   do i=1,nx   
+
+      if ((.not.ilemy(i,j)).and.(fleuvemy(i,j))) gzmy(i,j)=.true. 
+
+
+ ! calcul du frottement basal (ce bloc etait avant dans neffect)
+
+      if (cotemy(i,j)) then                        ! point cotier
+         betamy(i,j)=cf*neffmy(i,j) 
+! cat         betamy(i,j)=min(betamy(i,j),tobmax)
+! version Jorge
+! point cotier: frott maxi=40, mini=4
+!         betamy(i,j)=min(betamy(i,j),betacotemax)
+!         betamy(i,j)=max(betamy(i,j),betacotemin)
+! version Jorge stream  en zone sediment :
+      else if ((fleuve_sedimy(i,j).and.gzmx(i,j).and.(.not.cotemx(i,j)))) then  ! stream  en zone sediment
+! stream sediment: frott maxi=60, mini=6 
+!cdc            betamy(i,j)=cf*neffmy(i,j)*3
+             betamy(i,j)=cf*neffmy(i,j)*coef_sedim
+!            betamy(i,j)=min(betamy(i,j),betasedmax)          
+!            betamy(i,j)=max(betamy(i,j),betasedmin)
+      else if ((gzmy(i,j)).and.(.not.cotemy(i,j))) then  ! stream normal
+! cat         betamy(i,j)=tobmax
+! cat         betamy(i,j)=cf*neffmy(i,j)*20. 
+! cat         betamy(i,j)=min(betamy(i,j),tobmax)
+! cat         betamy(i,j)=max(betamy(i,j),10.)
+! version Jorge
+! stream normal: frott maxi=80, mini=8
+         betamy(i,j)=cf*neffmy(i,j)*coef_gz                     
+!         betamy(i,j)=min(betamy(i,j),betastreamax)           
+!         betamy(i,j)=max(betamy(i,j),betastreamin)
+
+! cat         if (cf*neffmy(i,j).gt.1500.) then
+! cat            gzmy(i,j)=.false.
+! cat            betamy(i,j)=tostick
+! cat         endif
+
+      else if (ilemy(i,j)) then
+         betamy(i,j)=cf*neffmy(i,j)*coef_ile
+! cat         betamy(i,j)=min(betamy(i,j),tob_ile)     
+! version Jorge
+! ile : frott maxi=20, mini=2   
+!         betamy(i,j)=min(betamy(i,j),betailemax)
+!         betamy(i,j)=max(betamy(i,j),betailemin)
+      else if (flotmy(i,j)) then     ! flottant ou ile
+         betamy(i,j)=0.
+      else                               ! grounded, SIA
+         betamy(i,j)=tostick             ! frottement glace posee
+      endif
+   end do
+end do
+
+!print*,'dragging debug h_sedim(78,144),hwater,seuil_sedim,seuil_hwater',h_sedim(78,144),seuil_sedim,hwater(78,144),seuil_hwater,fleuve(78,144),gzmx(78,144)
+!print*,'dragging debug h_sedim(90,157),hwater,seuil_sedim,seuil_hwater',h_sedim(90,157),seuil_sedim,hwater(90,157),seuil_hwater,fleuve(90,157),gzmx(90,157),betamx(90,157),betamx(78,144)
+
+
+where (fleuve(:,:))
+   debug_3D(:,:,1)=1
+elsewhere (flot(:,:))
+   debug_3D(:,:,1)=2
+elsewhere
+   debug_3D(:,:,1)=0
+endwhere
+
+where (cote(:,:))
+   debug_3D(:,:,2)=1
+elsewhere
+   debug_3D(:,:,2)=0
+endwhere
+
+where (fleuvemx(:,:))
+   debug_3D(:,:,3)=1
+elsewhere
+   debug_3D(:,:,3)=0
+endwhere
+
+where (flotmx(:,:))
+   debug_3D(:,:,3)=-1
+endwhere
+
+!_____________________
+
+
+where (fleuvemy(:,:))
+   debug_3D(:,:,4)=1
+elsewhere
+   debug_3D(:,:,4)=0
+endwhere
+
+where (flotmy(:,:))
+   debug_3D(:,:,4)=-1
+end where
+
+debug_3D(:,:,23)= abs(RO*G*HMX(:,:)*sdx(:,:)*1.e-5) 
+debug_3D(:,:,24)= abs(RO*G*HMY(:,:)*sdy(:,:)*1.e-5) 
+
+
+!print*,'point (79,145)',hwater(79,145),h_sedim(79,145),flot(79,145),fleuve_sedim(79,145),fleuve(79,145),betamx(79,145),fleuve_sedimx(79,145),gzmx(79,145)
+
+return
+end subroutine dragging
+
+end module dragging_hwat_sedim
+
Index: /trunk/SOURCES/Old-sources/dragging_hwatermax_0.2_mod.f90
===================================================================
--- /trunk/SOURCES/Old-sources/dragging_hwatermax_0.2_mod.f90	(revision 331)
+++ /trunk/SOURCES/Old-sources/dragging_hwatermax_0.2_mod.f90	(revision 331)
@@ -0,0 +1,151 @@
+!> \file dragging_hwatermax_0.2_mod.f90
+!! Defini les zones de stream avec un critere
+!<
+
+!> \namespace  dragging_hwatstream
+!! Defini les zones de stream avec 1 critere
+!! \author ...
+!! \date ...
+!! @note Le critere est:
+!! @note   * un critere sur la hauteur d'eau 
+!! @note Used module
+!! @note   - use module3D_phy
+!<
+module dragging_hwatstream
+  ! Defini les zones de stream avec un critere sur la hauteur d'eau
+
+  use module3D_phy
+  use fake_beta_iter_vitbil_mod
+
+
+  implicit none
+  real :: tostick
+
+contains
+  !-------------------------------------------------------------------------------
+  !> SUBROUTINE: init_dragging
+  !!Cette routine fait l'initialisation du dragging.
+  !>
+  subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
+
+    implicit none
+
+    namelist/drag_hwat/hwatstream,cf,betamax,toblim,tostick    
+
+    if (itracebug.eq.1)  call tracebug(' Dragging avec hwatermax')
+
+
+    ! formats pour les ecritures dans 42
+428 format(A)
+
+    ! lecture des parametres du run                      block draghwat
+    !--------------------------------------------------------------------
+
+    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
+    read(num_param,drag_hwat)
+
+    write(num_rep_42,428)'!___________________________________________________________' 
+    write(num_rep_42,428) '&drag_hwat              ! nom du bloc hwatermax '
+    write(num_rep_42,*)
+    write(num_rep_42,*) 'hwatsream      = ',hwatstream 
+    write(num_rep_42,*) 'cf             = ',cf
+    write(num_rep_42,*) 'betamax         = ', betamax
+    write(num_rep_42,*) 'toblim         = ', toblim
+    write(num_rep_42,*)'/'                            
+    write(num_rep_42,428) '! hwatstream (m) :  critere de passage en stream si hwater > hwatstream '
+    write(num_rep_42,428) '! cf coefficient de la loi de frottement fonction Neff'
+    write(num_rep_42,428) '! betamax : (Pa) frottement maxi '
+    write(num_rep_42,428) '! toblim : (Pa) tes pour les iles '
+    write(num_rep_42,*)
+
+    inv_beta=0
+    tostick=1.e5   ! valeurs pour les points non flgzmx
+
+    !-------------------------------------------------------------------
+    ! masque stream
+
+    mstream_mx(:,:)=1
+    mstream_my(:,:)=1
+
+
+    ! coefficient permettant de modifier le basal drag.
+    drag_mx(:,:)=1.
+    drag_my(:,:)=1.
+
+
+    return
+  end subroutine init_dragging
+  !________________________________________________________________________________
+
+  !> SUBROUTINE: dragging
+  !!Defini la localisation des streams et le frottement basal
+  !>
+  !-------------------------------------------------------------------------
+  subroutine dragging   ! defini la localisation des streams et le frottement basal
+
+
+    !         les iles n'ont pas de condition neff mais ont la condition toblim
+    !         (ce bloc etait avant dans flottab)
+
+    do j=2,ny
+       do i=2,nx
+          ilemx(i,j)=ilemx(i,j).and.(abs(tobmx(i,j)).lt.toblim)
+          ilemy(i,j)=ilemy(i,j).and.(abs(tobmy(i,j)).lt.toblim)
+       end do
+    end do
+
+
+    gzm_x: do j=1,ny  ! 2,ny
+       do i=2,nx ! -1
+
+          ! gzmx pour les hauteurs d'eau superieures a un seuil
+          ! attention, il semble qu'il y avait une erreur dans la precedente
+          ! formulation car le test etait seulement sur hwater(i,j).
+          ! le gzmx initial a ete calculé dans flottab pour les points a cheval sur la grounding line
+
+          gzmx(i,j)=gzmx(i,j).or. & 
+               (((hwater(i,j)+hwater(i-1,j))*0.5.gt.hwatstream) &
+               .and..not.flotmx(i,j))
+		  fleuvemx(i,j)=gzmx(i,j)
+
+          ! calcul du frottement basal (ce bloc etait avant dans neffect)
+          if (gzmx(i,j)) then                          ! stream
+             betamx(i,j)=cf*neffmx(i,j) 
+          else if (flotmx(i,j).or.ilemx(i,j)) then     ! flottant ou ile
+             betamx(i,j)=0.
+          else                                         ! grounded, SIA
+             betamx(i,j)=tostick                       ! frottement glace posee (1 bar)
+          endif
+
+       end do
+    end do gzm_x
+
+    gzm_y: do j=2,ny  !-1
+       do i=1,nx  ! 2,nx
+
+          ! gzmy pour les hauteurs d'eau superieures a un seuil
+          ! attention, il semble qu'il y avait une erreur dans la precedente
+          ! formulation car le test etait seulement sur hwater(i,j).
+
+          gzmy(i,j)=gzmy(i,j).or. &
+               (((hwater(i,j)+hwater(i,j-1))*0.5.gt.hwatstream) &
+               .and..not.flotmy(i,j))
+          fleuvemy(i,j)=gzmx(i,j)
+
+          ! calcul du frottement basal (ce bloc etait avant dans neffect)
+          if (gzmy(i,j)) then                ! stream
+             betamy(i,j)=cf*neffmy(i,j) 
+          else if (flotmy(i,j).or.ilemy(i,j)) then     ! flottant ou ile
+             betamy(i,j)=0.
+          else                               ! grounded, SIA
+             betamy(i,j)=tostick             ! frottement glace posee
+          endif
+
+       end do
+    end do gzm_y
+
+    return
+  end subroutine dragging
+
+end module dragging_hwatstream
+
Index: /trunk/SOURCES/Old-sources/dragging_neff_slope_mod.f90
===================================================================
--- /trunk/SOURCES/Old-sources/dragging_neff_slope_mod.f90	(revision 331)
+++ /trunk/SOURCES/Old-sources/dragging_neff_slope_mod.f90	(revision 331)
@@ -0,0 +1,422 @@
+!> \file dragging_neff_contmaj_mod.f90
+!! Defini les zones de stream avec 3 criteres 
+!<
+
+!> \namespace  dragging_neff_contmaj
+!! Defini les zones de stream avec 3 criteres
+!! \author ...
+!! \date ...
+!! @note Les trois criteres sont:
+!! @note   * un critere sur la pression effective 
+!! @note   * un critere de continuite depuis la cote
+!! @note   * un critere sur la courbure du socle (si negatif, vallees) 
+!! @note Used module
+!! @note   - use module3D_phy
+!! @note   - use param_phy_mod
+!<
+
+module dragging_neff_slope
+
+! Defini les zones de stream avec :
+! * un critere sur la pression effective 
+! * un critere de continuite depuis la cote
+! * un critere sur la courbure du socle (si negatif, vallees) 
+
+use module3d_phy
+use param_phy_mod
+use interface_input
+use io_netcdf_grisli
+use fake_beta_iter_vitbil_mod
+
+implicit none
+
+logical,dimension(nx,ny) :: fleuve
+logical,dimension(nx,ny) :: cote
+logical,dimension(nx,ny) :: slowssamx ! not actual stream, but ssa used as Ub
+logical,dimension(nx,ny) :: slowssamy ! not actual stream, but ssa used as Ub
+logical,dimension(nx,ny) :: slowssa   ! not actual stream, but ssa used as Ub
+
+real,dimension(nx,ny) :: hires_slope ! slope comupted on a high resolution topography
+character(len=100) :: slope_fich     ! fichier grille
+character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat
+
+real :: valmax
+integer :: imax,jmax
+integer :: i_moins1,i_plus1,j_moins1,j_plus1
+integer :: lmax=20
+integer :: idep,jdep,iloc,jloc
+
+
+real :: tostick   ! pour la glace posee
+real :: betamin   ! betamin : (Pa) frottement mini sous les streams
+real :: tob_ile    ! pour les iles
+real :: cry_lim=50.  ! courbure limite pour le suivi des fleuves
+
+real,dimension(nx,ny) :: neff ! pression effective noeuds majeurs
+real :: seuil_neff    ! seuil sur la pression effective pour avoir glissement en zone sediment
+real :: coef_gz       ! coef frottement zones stream std  (10)
+real :: coef_ile      ! coef frottement zones iles   (0.1)
+
+real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite avec spinup cat
+real, dimension(nx,ny) :: Vcol_y           !< uniquement pour compatibilite avec spinup cat
+real, dimension(nx,ny) :: Vsl_x            !< uniquement pour compatibilite avec spinup cat
+real, dimension(nx,ny) :: Vsl_y            !< uniquement pour compatibilite avec spinup cat
+logical :: corr_def = .false.              !< for deformation correction, pour compatibilite beta
+
+contains
+!-------------------------------------------------------------------------------
+
+!> SUBROUTINE: init_dragging
+!! Cette routine fait l'initialisation du dragging.
+!>
+subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
+
+implicit none
+
+! local variables, defining the rugosity-enhanced flow
+real :: expo_slope
+real :: pente_min, pente_max
+
+namelist/drag_neff_slope/cf,betamax,betamin,toblim,tostick,seuil_neff,coef_gz,coef_ile,slope_fich,expo_slope,pente_min,pente_max
+
+if (itracebug.eq.1)  call tracebug(' dragging avec neff et slope')
+
+
+! formats pour les ecritures dans 42
+428 format(A)
+
+! lecture des parametres du run                      block drag neff
+!--------------------------------------------------------------------
+
+rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
+read(num_param,drag_neff_slope)
+
+write(num_rep_42,428)'!___________________________________________________________' 
+write(num_rep_42,428) '&drag_neff_slope          ! nom du bloc dragging neff slope'
+write(num_rep_42,*)
+write(num_rep_42,*) 'cf               = ',cf
+write(num_rep_42,*) 'betamax          = ', betamax
+write(num_rep_42,*) 'betamin          = ', betamin
+write(num_rep_42,*) 'toblim           = ', toblim
+write(num_rep_42,*) 'tostick          = ', tostick
+write(num_rep_42,*) 'seuil_neff       = ', seuil_neff
+write(num_rep_42,*) 'coef_gz          = ', coef_gz
+write(num_rep_42,*) 'coef_ile         = ', coef_ile
+write(num_rep_42,'(A,A)') 'slope_fich = ', slope_fich
+write(num_rep_42,*) 'expo_slope       = ', expo_slope
+write(num_rep_42,*) 'pente_min        = ', pente_min
+write(num_rep_42,*) 'pente_max        = ', pente_max
+write(num_rep_42,*)'/'                            
+write(num_rep_42,428) '! cf coefficient de la loi de frottement fonction Neff'
+write(num_rep_42,428) '! betamax : (Pa) frottement maxi sous les streams '
+write(num_rep_42,428) '! betamin : (Pa) frottement mini sous les streams '
+write(num_rep_42,428) '! toblim : (Pa)  pour les iles '
+write(num_rep_42,428) '! tostick : (Pa)  pour les points non flgzmx '
+write(num_rep_42,428) '! seuil_neff (Pa) seuil sur la pression effective pour avoir glissement'
+write(num_rep_42,428) '! coef_gz : coef frottement zones stream std'
+write(num_rep_42,428) '! coef_ile : coef frottement zones iles'
+write(num_rep_42,428) '! slope_fich : fichier pente bedrock'
+write(num_rep_42,428) '! expo_slope : exposant fonction pente'
+write(num_rep_42,428) '! pente_min : no flow enhancement below this'
+write(num_rep_42,428) '! pente_max : maximum flow enhancement above this'
+write(num_rep_42,*)
+
+inv_beta=0
+tob_ile=betamax/2.
+moteurmax=toblim
+
+! betamax_2d depends on sub-grid slopes.
+slope_fich=trim(dirnameinp)//trim(slope_fich)
+call lect_input(1,'slope',1,hires_slope(:,:),slope_fich,file_ncdf)
+
+! from slopes, we create an index between 0 & 1
+! 1 is very mountainous, 0 is flat
+hires_slope(:,:)=(max(min(hires_slope(:,:),pente_max),pente_min)-pente_min)/(pente_max-pente_min)
+! now we compute the actual betamax used by the remplimat routines
+! /!\ the slope is used to modify the beta where we have a temperate base,
+! but NO ice stream... -> Slow SSA zone (SSA used to compute Ub)
+betamax_2d(:,:)=max ( tostick * (1. - (1 - betamax / tostick ) * hires_slope(:,:)**(1./expo_slope)) , betamax )
+!do j=1,ny
+!   do i=1,nx
+!      write(18745,*) betamax_2d (i,j)
+!   enddo
+!enddo
+
+!-------------------------------------------------------------------
+! masque stream
+      mstream_mx(:,:)=1
+      mstream_my(:,:)=1
+
+
+! coefficient permettant de modifier le basal drag.
+      drag_mx(:,:)=1.
+      drag_my(:,:)=1.
+
+where (socle_cry(:,:).lt.cry_lim)
+   debug_3D(:,:,25)=1
+elsewhere
+   debug_3D(:,:,25)=0
+endwhere
+
+   
+      return
+      end subroutine init_dragging
+!________________________________________________________________________________
+
+!> SUBROUTINE: dragging
+!! Defini la localisation des streams et le frottement basal
+!>
+!-------------------------------------------------------------------------
+subroutine dragging   ! defini la localisation des streams et le frottement basal
+!$ USE OMP_LIB
+
+!         les iles n'ont pas de condition neff mais ont la condition toblim
+!         (ce bloc etait avant dans flottab)
+
+!$OMP PARALLEL
+!$OMP DO
+do j=2,ny
+   do i=2,nx
+      ilemx(i,j)=ilemx(i,j).and.(abs(rog*Hmx(i,j)*sdx(i,j)).lt.toblim)
+      ilemy(i,j)=ilemy(i,j).and.(abs(rog*Hmy(i,j)*sdy(i,j)).lt.toblim)
+   end do
+end do
+!$OMP END DO
+
+! le calcul des fleuves se fait sur les mailles majeures en partant de la cote
+
+
+
+
+! le gzmx initial a ete calculé dans flottab pour les points a cheval sur la grounding line
+
+! Partant de la cote, le fleuve est determine en remontant vers le point de hauteur d'eau max
+! tant que ce point est superieur a hwatermax.
+
+! Attention : ce systeme ne permet pas d'avoir des fleuves convergents
+!             et les fleuves n'ont une épaisseur que de 1 noeud.
+!$OMP WORKSHARE
+fleuvemx(:,:)=.false.
+fleuvemy(:,:)=.false.
+fleuve(:,:)=.false.
+cote(:,:)=.false.
+!$OMP END WORKSHARE
+
+! detection des cotes
+!$OMP DO
+do  j=2,ny-1  
+   do i=2,nx-1  
+      if ((.not.flot(i,j)).and. & 
+      ((flot(i+1,j)).or.(flot(i,j+1)).or.(flot(i-1,j)).or.(flot(i,j-1)))) then
+         cote(i,j)=.true.
+      endif
+   end do
+end do
+!$OMP END DO
+
+! calcul de neff (pression effective sur noeuds majeurs)
+!$OMP DO
+do j=1,ny-1
+   do i=1,nx-1
+        neff(i,j)=(neffmx(i,j)+neffmx(i+1,j)+neffmy(i,j)+neffmy(i,j+1))/4
+   enddo
+enddo
+!$OMP END DO
+!aurel, for the borders:
+!$OMP WORKSHARE
+neff(:,ny)=1.e5
+neff(nx,:)=1.e5
+
+
+! aurel, we add the neff threshold:
+where ((neff(:,:).le.seuil_neff).and.(.not.flot(:,:)).and.(H(:,:).gt.1.)) fleuve(:,:)=.true.
+!$OMP END WORKSHARE
+
+!$OMP DO
+do j=1,ny-1
+   do i=1,nx-1
+      if (fleuve(i,j)) then
+         fleuvemx(i,j)=.true.
+         fleuvemx(i+1,j)=.true.
+         fleuvemy(i,j)=.true.
+         fleuvemy(i,j+1)=.true.
+      end if
+   end do
+end do
+!$OMP END DO
+
+! we look for the warm based points that will not be treated as stream (ub from SSA):
+!$OMP WORKSHARE
+slowssa(:,:)=.false.
+slowssamx(:,:)=.false.
+slowssamy(:,:)=.false.
+!$OMP END WORKSHARE
+!$OMP DO
+do  j=1,ny
+   do i=1,nx
+      !if ((not(flot(i,j))).and.(hwater(i,j).gt.0.1)) slowssa(i,j)=.true.
+      if ((.not.(flot(i,j))).and.(ibase(i,j).ne.1).and.(H(i,j).gt.1.)) slowssa(i,j)=.true.
+   end do
+end do
+!$OMP END DO
+!$OMP DO
+do j=1,ny-1
+   do i=1,nx-1
+      if (slowssa(i,j)) then
+         slowssamx(i,j)=.true.
+         slowssamx(i+1,j)=.true.
+         slowssamy(i,j)=.true.
+         slowssamy(i,j+1)=.true.
+      end if
+   end do
+end do
+!$OMP END DO
+!$OMP DO
+do j=1,ny
+   do i=1,nx
+
+      ! the actual streams and the warm based points are gzmx now:
+      if ( ((.not.ilemx(i,j)).and.(fleuvemx(i,j))) .or. ((.not.ilemx(i,j)).and.(slowssamx(i,j))) ) gzmx(i,j)=.true. 
+
+      
+! calcul du frottement basal (ce bloc etait avant dans neffect)
+
+      if (cotemx(i,j)) then                        ! point cotier
+         betamx(i,j)=cf*neffmx(i,j) 
+         betamx(i,j)=min(betamx(i,j),betamax)
+
+      else if ((gzmx(i,j)).and.(.not.cotemx(i,j))) then ! tous les points SSA
+
+         if (fleuvemx(i,j)) then                 ! the actual streams
+            betamx(i,j)=cf*neffmx(i,j)*coef_gz
+            betamx(i,j)=min(betamx(i,j),betamax)
+            betamx(i,j)=max(betamx(i,j),betamin)
+
+            if (cf*neffmx(i,j).gt.betamax*2.) then  ! a stream that's becoming grounded...
+               if (slowssamx(i,j)) then
+                  betamx(i,j)=betamax_2d(i,j)
+               else
+                  gzmx(i,j)=.false.
+                  betamx(i,j)=tostick
+               endif
+            endif
+
+         else                    ! not an actual stream, SSA is used to compute Ub
+            betamx(i,j)=betamax_2d(i,j)
+         endif
+         
+      else if (ilemx(i,j)) then
+         betamx(i,j)=cf*neffmx(i,j)*coef_ile 
+         betamx(i,j)=min(betamx(i,j),tob_ile)       
+      else if (flotmx(i,j)) then     ! flottant ou ile
+         betamx(i,j)=0.
+      else                                         ! grounded, SIA
+         betamx(i,j)=tostick                       ! frottement glace posee (1 bar)
+      endif
+
+   end do
+end do
+!$OMP END DO
+
+!$OMP DO
+do j=1,ny
+   do i=1,nx
+
+      ! the actual streams and the warm based points are gzmx now:
+      if ( ((.not.ilemy(i,j)).and.(fleuvemy(i,j))) .or. ((.not.ilemy(i,j)).and.(slowssamy(i,j))) ) gzmy(i,j)=.true. 
+
+      
+! calcul du frottement basal (ce bloc etait avant dans neffect)
+
+      if (cotemy(i,j)) then                        ! point cotier
+         betamy(i,j)=cf*neffmy(i,j) 
+         betamy(i,j)=min(betamy(i,j),betamax)
+
+      else if ((gzmy(i,j)).and.(.not.cotemy(i,j))) then ! tous les points SSA
+
+         if (fleuvemy(i,j)) then                 ! the actual streams
+            betamy(i,j)=cf*neffmy(i,j)*coef_gz
+            betamy(i,j)=min(betamy(i,j),betamax)
+            betamy(i,j)=max(betamy(i,j),betamin)
+
+            if (cf*neffmy(i,j).gt.betamax*2.) then  ! a stream that's becoming grounded...
+               if (slowssamy(i,j)) then
+                  betamy(i,j)=betamax_2d(i,j)
+               else
+                  gzmy(i,j)=.false.
+                  betamy(i,j)=tostick
+               endif
+            endif
+
+         else                    ! not an actual stream, SSA is used to compute Ub
+            betamy(i,j)=betamax_2d(i,j)
+         endif
+         
+      else if (ilemy(i,j)) then
+         betamy(i,j)=cf*neffmy(i,j)*coef_ile 
+         betamy(i,j)=min(betamy(i,j),tob_ile)       
+      else if (flotmy(i,j)) then     ! flottant ou ile
+         betamy(i,j)=0.
+      else                                         ! grounded, SIA
+         betamy(i,j)=tostick                       ! frottement glace posee (1 bar)
+      endif
+
+   end do
+end do
+!$OMP END DO
+
+!$OMP DO
+do j=2,ny-1
+   do i=2,nx-1
+      beta_centre(i,j) = ((betamx(i,j)+betamx(i+1,j)) &
+          + (betamy(i,j)+betamy(i,j+1)))*0.25
+   end do
+end do
+!$OMP END DO
+!$OMP END PARALLEL
+
+!~ where (fleuve(:,:))
+!~    debug_3D(:,:,1)=1
+!~ elsewhere (flot(:,:))
+!~    debug_3D(:,:,1)=2
+!~ elsewhere
+!~    debug_3D(:,:,1)=0
+!~ endwhere
+!~ 
+!~ where (cote(:,:))
+!~    debug_3D(:,:,2)=1
+!~ elsewhere
+!~    debug_3D(:,:,2)=0
+!~ endwhere
+!~ 
+!~ where (fleuvemx(:,:))
+!~    debug_3D(:,:,3)=1
+!~ elsewhere
+!~    debug_3D(:,:,3)=0
+!~ endwhere
+!~ 
+!~ where (flotmx(:,:))
+!~    debug_3D(:,:,3)=-1
+!~ endwhere
+
+!_____________________
+
+
+!~ where (fleuvemy(:,:))
+!~    debug_3D(:,:,4)=1
+!~ elsewhere
+!~    debug_3D(:,:,4)=0
+!~ endwhere
+!~ 
+!~ where (flotmy(:,:))
+!~    debug_3D(:,:,4)=-1
+!~ end where
+!~ 
+!~ debug_3D(:,:,23)= abs(RO*G*HMX(:,:)*sdx(:,:)*1.e-5) 
+!~ debug_3D(:,:,24)= abs(RO*G*HMY(:,:)*sdy(:,:)*1.e-5) 
+
+return
+end subroutine dragging
+
+end module dragging_neff_slope
+
Index: /trunk/SOURCES/Old-sources/dragging_prescr_beta_mod_Ant.f90
===================================================================
--- /trunk/SOURCES/Old-sources/dragging_prescr_beta_mod_Ant.f90	(revision 331)
+++ /trunk/SOURCES/Old-sources/dragging_prescr_beta_mod_Ant.f90	(revision 331)
@@ -0,0 +1,530 @@
+!> \file dragging_prescr_beta_mod.f90
+!! Module pour calculer le beta a partir de vitesses de bilan
+!<
+
+!> \namespace dragging_prescr_beta
+!! Calcule le beta a partir de vitesses de bilan
+!! @note Il faut partir d'un cptr 
+!! \author ...
+!! \date ...
+!! @note Used module
+!! @note   - use module3D_phy
+!<
+
+module dragging_prescr_beta
+
+  ! Calcule le beta a partir de vitesses de bilan
+  ! il faut partir d'un cptr.
+
+  use module3d_phy
+  use interface_input
+  implicit none
+  real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite
+  real, dimension(nx,ny) :: Vcol_y           !<
+  real, dimension(nx,ny) :: Vsl_x            !<
+  real, dimension(nx,ny) :: Vsl_y            !<
+
+
+  real, dimension(nx,ny) :: drag_centre      !< beta on major node (average)
+  !hassine !real, dimension(nx,ny) :: beta_centre      ! beta on major node (average)
+  real, dimension(nx,ny) :: betamoy_x        !< pour faire les moyennes glissantes
+  real, dimension(nx,ny) :: betamoy_y        !<
+  integer, dimension(nx,ny) :: nbvoisx       !<
+  integer, dimension(nx,ny) :: nbvoisy       !<
+
+  real, dimension(nx,ny) ::  corrdS_centre   !< correction dS 
+  real, dimension(nx,ny) ::  corrdS_mx       !<
+  real, dimension(nx,ny) ::  corrdS_my       !<
+  real, dimension(nx,ny) ::  fleuves_nommes  !< each ice stream is numbered
+
+  real, dimension(nx,ny) ::  corrstream_centre    !< pour correction stream
+  real, dimension(nx,ny) ::  corrstream_mx        !<
+  real, dimension(nx,ny) ::  corrstream_my        !<
+  integer :: numstream
+
+  integer,parameter             :: n_g = 2        !< taille de cette matrice (-n_g:n_g,-n_g:n_g)
+  real, dimension(-2:2,-2:2)    :: M_gauss        !< matrice qui contient la gaussienne pour iterations beta
+  real                          :: sigma_beta=1   !< sigma de la Gaussienne
+
+
+  real, dimension(nx,ny) :: drag_instab_topo !< regions detected as topo instability
+  real :: beta_limgz                         !< when beta gt beta_limgz -> not gzmx
+  real :: beta_min                           !< minimum value of beta
+  real :: beta_mult                          !< coefficient of beta field
+  integer  :: ill,jll,nmoy
+
+  real :: maxi                               !< calcul correction dS a la grounding line
+  real :: mini
+
+  logical :: corr_def = .False.               !< for deformation correction (compatibility)
+
+contains
+
+  !-----------------------------------------------------------------------------------------
+  !> SUBROUTINE: init_dragging 
+  !! Cette routine fait l'initialisation du dragging.
+  !<
+  subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
+                                ! nouvelle version : lit les fichiers nc
+    implicit none
+    character(len=100) :: beta_c_file               !  beta on centered grid
+    character(len=100) :: betamx_file               !  beta on staggered grid mx
+    character(len=100) :: betamy_file               !  beta on staggered grid mx
+
+    namelist/beta_prescr/beta_c_file,beta_limgz,beta_min,beta_mult
+
+    if (itracebug.eq.1)  call tracebug(' Prescr_beta       subroutine init_dragging')
+
+    iter_beta=0
+
+    ! lecture des parametres du run                    
+    !--------------------------------------------------------------------
+
+    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
+428 format(A)
+    read(num_param,beta_prescr)
+
+    write(num_rep_42,428) '!___________________________________________________________' 
+    write(num_rep_42,428) '!  read beta on centered grid'
+    write(num_rep_42,beta_prescr)                    
+    write(num_rep_42,428) '!  beta_file : nom des fichier qui contiennent les beta_c'
+    write(num_rep_42,428) '!  above beta_limgz, gzmx is false'
+    write(num_rep_42,428) '!  if grounded, beta_min minimum value '
+    write(num_rep_42,428) '!  beta_mult : coefficient just after reading' 
+    write(num_rep_42,*)
+    write(num_rep_42,428) '!___________________________________________________________' 
+   
+
+    beta_c_file = trim(dirnameinp)//trim(beta_c_file)
+    betamx_file = trim(dirnameinp)//trim(betamx_file)
+    betamy_file = trim(dirnameinp)//trim(betamy_file)
+    betamax = beta_limgz
+
+!  read the beta value on centered and staggered grids
+!    call lect_input(1,'beta_c',1,drag_centre,beta_c_file,trim(dirnameinp)//trim(runname)//'.nc')
+!    call lect_input(1,'betamx',1,drag_mx,betamx_file,trim(dirnameinp)//trim(runname)//'.nc')
+!    call lect_input(1,'betamy',1,drag_my,betamy_file,trim(dirnameinp)//trim(runname)//'.nc')
+
+
+    call lect_input(1,'beta_c',1,drag_centre,beta_c_file,'')
+
+! multiplication par le coefficient beta_mult
+    if (itracebug.eq.1)  call tracebug(' Prescr_beta       apres lecture')
+    drag_centre(:,:) = drag_centre(:,:)*beta_mult
+
+
+!    call lect_input(1,'betamx',1,drag_mx,betamx_file,'')
+!    call lect_input(1,'betamy',1,drag_my,betamy_file,'')
+
+
+     call  beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre)         ! redistribute on staggered grid
+
+     if  (beta_limgz.gt.0.) then
+        beta_centre(:,:) =  drag_centre(:,:)
+        betamx(:,:)      =  drag_mx(:,:)
+        betamy(:,:)      =  drag_my(:,:)
+
+        if (itracebug.eq.1)  call tracebug(' Prescr_beta       apres call beta_c2stag')
+
+     else if (beta_limgz.lt.0.) then                              ! attention sans doute pour plastic
+
+        beta_centre(:,:) = - beta_limgz
+        betamx(:,:)      = - beta_limgz
+        betamy(:,:)      = - beta_limgz
+        drag_centre(:,:) = - beta_limgz 
+        drag_mx(:,:)     = - beta_limgz
+        drag_my(:,:)     = - beta_limgz 
+        beta_limgz = 1.
+     end if
+
+
+! divers traitements qui peuvent etre appliques a beta : mis sous forme de routine
+
+    ! call beta_plateau_slope   !  give value beta_limgz to places where slope is >1.e-2
+    ! call beta_stag2c          !  average the staggered values on central ones
+    ! call beta_plateau_S       !  give value beta_limgz to places where S > 2500 m
+
+   call beta_min_value(beta_min)      !  valeur mini que peut avoir beta (en Pa an m-1)
+                                ! on peut envisager une valeur ~ 10 
+                                ! rappel : beta doit être positif.
+
+    ! call  beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre)         ! redistribute on staggered grid
+
+
+    !  mstream attribution
+    !________________________
+    ! mstream -> regions where streaming is allowed
+    ! drag_mx -> value of beta if an ice stream is allowed
+    ! betamx  -> value of dragging (may depend on other conditions)
+
+    ! to ensure having some "Dirichlet" points, there must be at least 
+    ! a mstream=0 in every island
+
+
+! un tour de lissage
+!    call beta_stag2c(nx,ny,drag_mx,drag_my,drag_centre)                  ! average  staggered 
+!                                                                         ! on central 
+
+!    call beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre)                  ! centered distributed
+!                                                                         ! on staggered
+!
+    beta_centre(:,:) =  drag_centre(:,:)                                 ! surtout pour sorties
+
+
+
+    where (drag_mx(:,:).ge.beta_limgz)
+       mstream_mx(:,:) = 0
+       betamx(:,:)     = beta_limgz
+       drag_mx(:,:)    = beta_limgz
+    elsewhere
+       mstream_mx(:,:) = 1
+       betamx(:,:)     = drag_mx(:,:)
+    endwhere
+
+    where (drag_my(:,:).ge.beta_limgz)
+       mstream_my(:,:) = 0
+       betamy(:,:)     = beta_limgz
+       drag_my(:,:)    = beta_limgz
+    elsewhere
+       mstream_my(:,:) = 1
+       betamy(:,:)     = drag_my(:,:)
+    endwhere
+
+    if (itracebug.eq.1)  call tracebug(' Prescr_beta  apres lecture')
+    mstream_mx(:,:) = 0
+    mstream_my(:,:) = 0
+
+    if (itracebug.eq.1)  call tracebug(' Prescr_beta       avant call gzm')
+    call gzm_beta_prescr
+!    call write_datfile2(nx,ny,betamx,betamy,'beta-LBq15-08.dat')
+
+    call make_gaussienne(n_g,sigma_beta,M_gauss)
+STOP
+    return
+
+  end subroutine init_dragging
+  !________________________________________________________________________________
+
+  !> SUBROUTINE: dragging 
+  !! Defini la localisation des streams et le frottement basal
+  !<
+  !-------------------------------------------------------------------------
+  subroutine dragging   ! defini la localisation des streams et le frottement basal
+
+    implicit none
+
+
+
+    if (itracebug.eq.1)  call tracebug(' beta_prescr         subroutine dragging')
+
+
+!   tentative d'iteration
+
+!    call beta_stag2c(nx,ny,drag_mx,drag_my,drag_centre)                  ! average  staggered 
+                                                                          ! on central 
+!   call iter_gauss(n_g,nx,ny,M_gauss,-1.,drag_centre,corrstream_centre,hdot)  ! correct central
+!   call beta_c2stag(nx,ny,drag_mx,drag_my,corrstream_centre)             ! centered distributed
+                                                                          ! on staggered
+!    drag_centre(:,:) = corrstream_centre(:,:)
+!    beta_centre(:,:) =  drag_centre(:,:)                                 ! surtout pour sorties
+
+    ! pour l'instant pas d'autre condition que mstream
+    where (mstream_mx(:,:).eq.1)
+       betamx(:,:) = drag_mx(:,:)
+    elsewhere 
+       betamx(:,:)= beta_limgz
+    end where
+
+    where (mstream_my(:,:).eq.1)
+       betamy(:,:) = drag_my(:,:)
+    elsewhere 
+       betamy(:,:)= beta_limgz
+    end where
+
+    betamx(:,:)=drag_mx(:,:)
+    betamy(:,:)=drag_my(:,:)
+
+    call  gzm_beta_prescr        ! determine gz et flgz et met a 0 le beta des points flottants
+
+
+    return
+  end subroutine dragging
+
+  !----------------------------------------------------------------------------------
+  !>SUBROUTINE: gzm_beta_prescr
+  !! Calcul de gzmx
+  !<
+
+  subroutine gzm_beta_prescr     ! attribue flgzmx et gzmx (et my)
+    
+    logical ::  test_Tx          ! test if there is a basal melting point in the surrounding
+    logical ::  test_Ty          ! test if there is a basal melting point in the surrounding
+    logical ::  test_beta_x      ! test if there is a low drag point in the surrounding
+    logical ::  test_beta_y      ! test if there is a low drag point in the surrounding
+
+    real    ::  beta_forc_fleuve ! below this value -> ice stream
+
+    beta_forc_fleuve = 5.e3
+
+    ! calcul de gzmx 
+
+
+    gzmx(:,:)=.true.
+    gzmy(:,:)=.true.
+    flgzmx(:,:)=.false.
+    flgzmy(:,:)=.false.
+
+
+    ! points flottants : flgzmx mais pas gzmx
+    do j=2,ny
+       do i=2,nx
+          if (flot(i,j).and.(flot(i-1,j))) then
+             flotmx(i,j)=.true.
+             flgzmx(i,j)=.true.
+             gzmx(i,j)=.false.
+             betamx(i,j)=0.
+          end if
+          if (flot(i,j).and.(flot(i,j-1))) then
+             flotmy(i,j)=.true.
+             flgzmy(i,j)=.true.
+             gzmy(i,j)=.false.
+             betamy(i,j)=0.
+          end if
+       end do
+    end do
+
+    !--------- autres criteres
+
+    ! points poses
+    gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.beta_limgz)   !  Pas de calcul pour les points
+    gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.beta_limgz)   !  au fort beta
+
+
+
+    ! critere (lache) sur la temperature
+    do j = 2, ny-1
+       do i =2, nx-1 
+
+!   test s'il y a un point tempere dans les environs
+          test_Tx = any( ibase( i-1:i   , j-1:j+1 ).gt.1)
+          test_Ty = any( ibase( i-1:i+1 , j-1:j   ).gt.1)
+
+!   test s'il y a un point low drag dans les environs
+          test_beta_x = any( drag_centre( i-1:i   , j-1:j+1 ) .lt. beta_forc_fleuve )
+          test_beta_y = any( drag_centre( i-1:i+1 , j-1:j   ) .lt. beta_forc_fleuve )
+
+!   critere pour gzmx
+          gzmx(i,j) =  gzmx(i,j) .and. (test_Tx.or.test_beta_x)
+          gzmy(i,j) =  gzmy(i,j) .and. (test_Ty.or.test_beta_y)
+
+       end do
+    end do
+
+    flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
+    flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:)
+
+    where ((.not.flotmx(:,:)).and.(.not.gzmx(:,:)))
+       betamx(:,:) = beta_limgz
+    end where
+
+    where ((.not.flotmy(:,:)).and.(.not.gzmy(:,:)))
+       betamy(:,:) = beta_limgz
+    end where
+
+
+    
+  end subroutine gzm_beta_prescr
+
+
+  !----------------------------------------------------------------------------------
+  ! Routines qui permettent des manip simples sur beta.
+
+
+!______________________________________________________________________________________
+!>SUBROUTINE:  beta_plateau_slope 
+!! give value beta_limgz to places where slope is >1.e-2 
+!<
+
+  subroutine beta_plateau_slope      ! give value beta_limgz to places where slope is >1.e-2
+
+    call slope_surf
+
+    where (abs(sdx(:,:)).gt.1.e-2)
+       drag_mx(:,:) = beta_limgz
+    end where
+
+    where (abs(sdy(:,:)).gt.1.e-2)
+       drag_my(:,:) = beta_limgz
+    end where
+  end subroutine beta_plateau_slope
+!______________________________________________________________________________________
+!> SUBROUTINE beta_plateau_S 
+!! give value beta_limgz to places where S > 2500 m
+!<
+  subroutine beta_plateau_S     !  give value beta_limgz to places where S > 2500 m
+
+    where (S(:,:).gt.2500.)     ! moyen simple de nettoyer l'interieur. WAIS pas touche
+       drag_centre(:,:)=beta_limgz
+       drag_mx(:,:)=beta_limgz
+       drag_my(:,:)=beta_limgz
+    end where
+  end subroutine beta_plateau_S
+!______________________________________________________________________________________
+!>SUBROUTINE:   beta_min_value 
+!! valeur mini que peut avoir beta (en Pa an m-1) 
+!<
+
+  subroutine beta_min_value(val_betamin)
+
+  real :: val_betamin     ! valeur mini que peut avoir beta (en Pa an m-1)
+
+    drag_centre(:,:) = max(drag_centre(:,:),val_betamin)
+    drag_mx(:,:) = max(drag_mx(:,:),val_betamin)
+    drag_my(:,:) = max(drag_my(:,:),val_betamin)
+    where(flot(:,:))
+       drag_centre(:,:) = 0.
+    end where
+
+  end subroutine beta_min_value
+!______________________________________________________________________________________
+!>SUBROUTINE:   beta_stag2c 
+!! average the staggered values on central ones
+!<
+
+  subroutine beta_stag2c(nxx,nyy,R_mx,R_my,R_centre)      ! average the staggered values on central ones
+
+    implicit none
+    integer                                ::  nxx,nyy                ! dimension des tableaux
+    real, dimension(nxx,nyy), intent(in)   ::  R_mx                   ! valeur du beta sur maille mx
+    real, dimension(nxx,nyy), intent(in)   ::  R_my                   ! valeur du beta sur maille my
+    real, dimension(nxx,nyy), intent(out)  ::  R_centre               ! valeur du beta sur maille centree   
+    
+    R_centre(:,:)=0.
+    do j=2,ny-1
+       do i=2,nx-1
+          R_centre(i,j)= ((R_mx(i,j)+R_mx(i+1,j)) &
+               + (R_my(i,j)+R_my(i,j+1)))*0.25
+       end do
+    end do
+  end subroutine beta_stag2c
+!______________________________________________________________________________________
+!>SUBROUTINE:   beta_c2stag_min  
+!! redistribute central values on staggered grid
+!<
+
+  subroutine beta_c2stag(nxx,nyy,R_mx,R_my,R_centre)              ! redistribute central values on staggered grid
+
+    implicit none
+    integer                                ::  nxx,nyy                ! dimension des tableaux
+    real, dimension(nxx,nyy), intent(out)  ::  R_mx                   ! valeur du beta sur maille mx
+    real, dimension(nxx,nyy), intent(out)  ::  R_my                   ! valeur du beta sur maille my
+    real, dimension(nxx,nyy), intent(in)   ::  R_centre               ! valeur du beta sur maille centree
+
+    do j=2,nyy 
+       do i=2,nxx
+          R_mx(i,j)= (R_centre(i-1,j)+R_centre(i,j))*0.5
+          R_my(i,j)= (R_centre(i,j-1)+R_centre(i,j))*0.5
+       end do
+    end do
+
+  end subroutine beta_c2stag
+!______________________________________________________________________________________
+!>SUBROUTINE:   beta_c2stag_min  
+!! redistribute central values on staggered grid
+!<
+
+  subroutine beta_c2stag_min(nxx,nyy,R_mx,R_my,R_centre)              ! redistribute central values on staggered grid
+
+    implicit none
+    integer                                ::  nxx,nyy                ! dimension des tableaux
+    real, dimension(nxx,nyy), intent(inout)::  R_mx                   ! valeur du beta sur maille mx
+    real, dimension(nxx,nyy), intent(inout)::  R_my                   ! valeur du beta sur maille my
+    real, dimension(nxx,nyy), intent(in)   ::  R_centre               ! valeur du beta sur maille centree
+
+    do j=2,nyy 
+       do i=2,nxx
+          R_mx(i,j)= min(R_mx(i,j),(R_centre(i-1,j)+R_centre(i,j))*0.5)
+          R_my(i,j)= min(R_my(i,j),(R_centre(i,j-1)+R_centre(i,j))*0.5)
+       end do
+    end do
+
+  end subroutine beta_c2stag_min
+!______________________________________________________________________________________
+!>SUBROUTINE:   make_gaussienne 
+!! cree une matrice gaussienne normalisee a 1
+!<
+
+  subroutine  make_gaussienne(n,sigma,M)
+
+    implicit none
+    integer,intent(in)                       :: n          ! demi taille du tableau gaussienne
+    real, intent(in)                         :: sigma      ! sigma gaussienne exprime en dx
+    real,dimension(-n:n,-n:n),intent(inout)  :: M          ! tableau gaussienne
+
+
+!   variables locales
+    real                                     :: dist2      ! distance au centre (au carre)
+    real                                     :: som        ! pour calibration
+    real                                     :: sigma2     ! pour gaussienne
+
+    som = 0.
+    sigma2 = sigma * sigma * 2
+    do j = -n, n
+       do i = -n, n
+          dist2 = (i*i+j*j)
+          M(i,j) = exp(-dist2/sigma2)
+          som=som+M(i,j)
+       end do
+    end do
+
+    M(:,:) = M(:,:) / som
+
+!    do j = -n, n
+!       write(6,'(5(f0.4,1x))') (M(i,j),i=-n,n)
+!    end do
+
+  end subroutine make_gaussienne
+
+!______________________________________________________________________________________
+!>SUBROUTINE:  iter_gauss
+!! utilise hdot pour faire les iterations sur beta_centre
+!<
+  subroutine iter_gauss(n,nxx,nyy,M,coef,beta_in,beta_out,Ecart)
+
+    integer,                  intent(in)    :: n          ! taille matrice gaussienne 
+    integer,                  intent(in)    :: nxx,nyy    ! taille du domaine
+    real,dimension(-n:n,-n:n),intent(in)    :: M          ! matrice gaussienne 
+    real,                     intent(in)    :: coef       ! coefficient multiplicateur
+    real,dimension(nx,ny),    intent(in)    :: beta_in    ! la valeur initiale
+    real,dimension(nx,ny),    intent(out)   :: beta_out   ! la valeur apres iteration
+    real,dimension(nx,ny),    intent(in)    :: Ecart      ! l'ecart a l'obs
+
+! variables locales
+    integer                                 :: ii,jj
+
+    do j= 1 + n, nyy - n
+       do i = 1 + n, nxx - n
+          beta_out(i,j) = beta_in(i,j)
+
+          do jj = -n, n
+             do ii = -n, n 
+                beta_out(i,j) = beta_out(i,j) + coef * M(ii,jj) * ecart(i+ii,j+jj)
+             end do
+          end do
+
+       end do
+    end do
+
+    beta_out(:,:) = max (beta_out(:,:),0.)
+  end subroutine iter_gauss
+
+end module dragging_prescr_beta
+
+!!$! add the nodes of topographical instability (on major nodes)
+! a garder eventuellement pour runs instabilite
+!!$beta_file = trim(dirnameinp)//'sensib_H_Gael-15km.dat'
+!!$call lect_datfile(nx,ny,drag_instab_topo,1,beta_file)   
+!!$
+!!$where (drag_instab_topo(:,:).gt.0.)
+!!$   drag_centre(:,:)=min(drag_centre(:,:),400.)
+!!$endwhere
Index: unk/SOURCES/dragging_hwat-contigu_mod.f90
===================================================================
--- /trunk/SOURCES/dragging_hwat-contigu_mod.f90	(revision 330)
+++ 	(revision )
@@ -1,279 +1,0 @@
-!> \file dragging_hwat-contigu_mod.f90
-!! Defini les zones de stream avec 2 criteres
-!<
-
-!> \namespace  dragging_hwat_cont
-!! Defini les zones de stream avec 2 criteres
-!! \author ...
-!! \date ... 
-!! @note Les 2 criteres sont: 
-!! @note   * un critere sur la hauteur d'eau
-!! @note   * un critere de continuite
-!! depuis la cote
-!! @note Used module
-!! @note   - use module3D_phy
-!<
-module dragging_hwat_cont
-
-! Defini les zones de stream avec un critere sur la hauteur d'eau + un critere de continuite
-! depuis la cote
-
-use module3d_phy
-use fake_beta_iter_vitbil_mod
-
-implicit none
-!logical,dimension(nx,ny) :: fleuvemx
-!logical,dimension(nx,ny) :: fleuvemy
-real,dimension(nx,ny) :: hwatmx
-real,dimension(nx,ny) :: hwatmY
-
-real :: valmax
-integer :: imax,jmax
-integer :: i_moins1,i_plus1,j_moins1,j_plus1
-integer :: lmax=20
-integer :: idep,jdep,iloc,jloc
-
-real :: tostick
-
-contains
-!-------------------------------------------------------------------------------
-!> SUBROUTINE: init_dragging
-!! Cette routine fait l'initialisation du dragging.
-!>
-subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
-
-implicit none
- 
-namelist/drag_hwat_cont/hwatstream,cf,betamax,toblim,tostick    
-
-if (itracebug.eq.1)  call tracebug(' Dragging avec hwatermax')
-
-
-! formats pour les ecritures dans 42
-428 format(A)
-
-! lecture des parametres du run                      block draghwat
-!--------------------------------------------------------------------
-
-rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
-read(num_param,drag_hwat_cont)
-
-write(num_rep_42,428)'!___________________________________________________________' 
-write(num_rep_42,428) '&drag_hwat_cont              ! nom du bloc hwatermax '
-write(num_rep_42,*)
-write(num_rep_42,*) 'hwatsream      = ',hwatstream 
-write(num_rep_42,*) 'cf             = ',cf
-write(num_rep_42,*) 'betamax         = ', betamax
-write(num_rep_42,*) 'toblim         = ', toblim
-write(num_rep_42,*)'/'                            
-write(num_rep_42,428) '! hwatstream (m) :  critere de passage en stream en partant de la cote'
-write(num_rep_42,428) '!  si hwater > hwatstream '
-write(num_rep_42,428) '! cf coefficient de la loi de frottement fonction Neff'
-write(num_rep_42,428) '! betamax : (Pa) frottement maxi '
-write(num_rep_42,428) '! toblim : (Pa) tes pour les iles '
-write(num_rep_42,*)
-
-inv_beta=0
-tostick=1.e5   ! valeurs pour les points non flgzmx
-
-!-------------------------------------------------------------------
-! masque stream
-
-      mstream_mx(:,:)=1
-      mstream_my(:,:)=1
-
-
-! coefficient permettant de modifier le basal drag.
-      drag_mx(:,:)=1.
-      drag_my(:,:)=1.
-
-   
-      return
-      end subroutine init_dragging
-!________________________________________________________________________________
-
-!> SUBROUTINE: dragging
-!! Defini la localisation des streams et le frottement basal
-!>
-!-------------------------------------------------------------------------
-subroutine dragging   ! defini la localisation des streams et le frottement basal
-
-
-!         les iles n'ont pas de condition neff mais ont la condition toblim
-!         (ce bloc etait avant dans flottab)
-
-
-do j=2,ny
-   do i=2,nx
-      ilemx(i,j)=ilemx(i,j).and.(abs(tobmx(i,j)).lt.toblim)
-      ilemy(i,j)=ilemy(i,j).and.(abs(tobmy(i,j)).lt.toblim)
-   end do
-end do
-
-! calcul de la hauteur d'eau sur les mailles mineures
-do j=1,ny
-   do i=2,nx
-      hwatmx(i,j)=0.5*(hwater(i,j)+hwater(i-1,j))
-   end do
-end do
-
-do j=2,ny
-   do i=1,nx
-      hwatmy(i,j)=0.5*(hwater(i,j)+hwater(i,j-1))
-   end do
-end do
-
-
-debug_3D(:,:,1)=hwatmx(:,:)
-debug_3D(:,:,2)=hwatmy(:,:)
-
-
-! le gzmx initial a ete calculé dans flottab pour les points a cheval sur la grounding line
-
-! Partant de la cote, le fleuve est determine en remontant vers le point de hauteur d'eau max
-! tant que ce point est superieur a hwatermax.
-
-! Attention : ce systeme ne permet pas d'avoir des fleuves convergents
-!             et les fleuves n'ont une épaisseur que de 1 noeud.
-
-fleuvemx(:,:)=.false.
-fleuvemy(:,:)=.false.
-
-fleuve_mx: do j=1,ny  
-   do i=2,nx                      ! le noeud 1 n'est pas attribue
-
-cote_mx :  if (gzmx(i,j)) then    ! a ce niveau les gzmx sont sur la grounding line
-         idep=i
-         jdep=j
-
-suit_x : do l=1,lmax         ! debut de la boucle de suivi, lmax longueur maxi des fleuves
-           i_moins1=max(idep-1,2)
-           j_moins1=max(jdep-1,1)
-           i_plus1=min(idep+1,nx)
-           j_plus1=min(jdep+1,ny)
-           
-           valmax=0.
-
-! recherche du max en excluant les points flottants
-         do jloc=j_moins1,j_plus1
-            do iloc=i_moins1,i_plus1
- 
-               if ((hwatmx(iloc,jloc).gt.valmax).and.(.not.flotmx(iloc,jloc))) then
-                  imax=iloc
-                  jmax=jloc
-                  valmax=hwatmx(iloc,jloc)
-               endif
-            end do
-         end do
-
-         if (hwatmx(imax,jmax).gt.hwatstream) then
-            fleuvemx(imax,jmax)=.true.
-            idep=imax
-            jdep=jmax
-         end if
-      end do suit_x
-
-   end if cote_mx
-
-if ((.not.ilemx(i,j)).and.(fleuvemx(i,j))) gzmx(i,j)=.true. 
-
-
-
-! calcul du frottement basal (ce bloc etait avant dans neffect)
-         if (gzmx(i,j)) then                          ! stream
-            betamx(i,j)=cf*neffmx(i,j) 
-         else if (flotmx(i,j).or.ilemx(i,j)) then     ! flottant ou ile
-            betamx(i,j)=0.
-         else                                         ! grounded, SIA
-            betamx(i,j)=tostick                       ! frottement glace posee (1 bar)
-         endif      
-
- end do
-end do fleuve_mx
-
-
-
-
-fleuve_my: do j=2,ny                  ! le noeud 1 n'est pas attribue
-   do i=1,nx   
-
-cote_my :  if (gzmy(i,j)) then    ! a ce niveau les gzmx sont sur la grounding line
-         idep=i
-         jdep=j
-
-suit_y : do l=1,lmax         ! debut de la boucle de suivi, lmax longueur maxi des fleuves
-           i_moins1=max(idep-1,1)
-           j_moins1=max(jdep-1,2)
-           i_plus1=min(idep+1,nx)
-           j_plus1=min(jdep+1,ny)
-           
-           valmax=0.
-
-! recherche du max en excluant les points flottants
-         do jloc=j_moins1,j_plus1
-            do iloc=i_moins1,i_plus1
- 
-               if ((hwatmy(iloc,jloc).gt.valmax).and.(.not.flotmy(iloc,jloc))) then
-                  imax=iloc
-                  jmax=jloc
-                  valmax=hwatmy(iloc,jloc)
-               endif
-            end do
-         end do
-
-         if (hwatmy(imax,jmax).gt.hwatstream) then
-            fleuvemy(imax,jmax)=.true.
-            idep=imax
-            jdep=jmax
-         end if
-      end do suit_y
-
-   end if cote_my
-
-if ((.not.ilemy(i,j)).and.(fleuvemy(i,j))) gzmy(i,j)=.true. 
-
-
- ! calcul du frottement basal (ce bloc etait avant dans neffect)
-   if (gzmy(i,j)) then                ! stream
-      betamy(i,j)=cf*neffmy(i,j) 
-   else if (flotmy(i,j).or.ilemy(i,j)) then     ! flottant ou ile
-      betamy(i,j)=0.
-   else                               ! grounded, SIA
-      betamy(i,j)=tostick             ! frottement glace posee
-   endif
-
-
- end do
-end do fleuve_my
-
-
-
-where (fleuvemx(:,:))
-   debug_3D(:,:,3)=1
-elsewhere
-   debug_3D(:,:,3)=0
-endwhere
-
-where (flotmx(:,:))
-   debug_3D(:,:,3)=-1
-endwhere
-
-!_____________________
-
-
-where (fleuvemy(:,:))
-   debug_3D(:,:,4)=1
-elsewhere
-   debug_3D(:,:,4)=0
-endwhere
-
-where (flotmy(:,:))
-   debug_3D(:,:,4)=-1
-end where
-
-
-return
-end subroutine dragging
-
-end module dragging_hwat_cont
-
Index: unk/SOURCES/dragging_hwat_contmaj_mod.f90
===================================================================
--- /trunk/SOURCES/dragging_hwat_contmaj_mod.f90	(revision 330)
+++ 	(revision )
@@ -1,391 +1,0 @@
-!> \file dragging_hwat_contmaj_mod.f90
-!! Defini les zones de stream avec 3 criteres 
-!<
-
-!> \namespace  dragging_hwat_contmaj
-!! Defini les zones de stream avec 3 criteres
-!! \author ...
-!! \date ...
-!! @note Les trois criteres sont:
-!! @note   * un critere sur la hauteur d'eau 
-!! @note   * un critere de continuite depuis la cote
-!! @note   * un critere sur la courbure du socle (si negatif, vallees) 
-!! @note Used module
-!! @note   - use module3D_phy
-!! @note   - use param_phy_mod
-!<
-
-module dragging_hwat_contmaj
-
-! Defini les zones de stream avec :
-! * un critere sur la hauteur d'eau 
-! * un critere de continuite depuis la cote
-! * un critere sur la courbure du socle (si negatif, vallees) 
-
-use module3d_phy
-use param_phy_mod
-use fake_beta_iter_vitbil_mod
-
-implicit none
-!logical,dimension(nx,ny) :: fleuvemx
-!logical,dimension(nx,ny) :: fleuvemy
-logical,dimension(nx,ny) :: fleuve
-logical,dimension(nx,ny) :: cote
-
-
-real :: valmax
-integer :: imax,jmax
-integer :: i_moins1,i_plus1,j_moins1,j_plus1
-integer :: lmax=20
-integer :: idep,jdep,iloc,jloc
-
-
-real :: tostick   ! pour la glace posee
-real :: tob_ile    ! pour les iles
-real :: cry_lim=50.  ! courbure limite pour le suivi des fleuves
-
-
-real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite avec spinup cat
-real, dimension(nx,ny) :: Vcol_y           !< uniquement pour compatibilite avec spinup cat
-real, dimension(nx,ny) :: Vsl_x            !< uniquement pour compatibilite avec spinup cat
-real, dimension(nx,ny) :: Vsl_y            !< uniquement pour compatibilite avec spinup cat
-logical :: corr_def = .false.              !< for deformation correction, pour compatibilite beta
-
-contains
-!-------------------------------------------------------------------------------
-
-!> SUBROUTINE: init_dragging
-!! Cette routine fait l'initialisation du dragging.
-!>
-subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
-
-implicit none
- 
-namelist/drag_hwat_cont/hwatstream,cf,betamax,toblim,tostick    
-
-if (itracebug.eq.1)  call tracebug(' dragging avec hwatermax')
-
-
-! formats pour les ecritures dans 42
-428 format(A)
-
-! lecture des parametres du run                      block draghwat
-!--------------------------------------------------------------------
-
-rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
-read(num_param,drag_hwat_cont)
-
-write(num_rep_42,428)'!___________________________________________________________' 
-write(num_rep_42,428) '&drag_hwat_cont              ! nom du bloc hwatermax '
-write(num_rep_42,*)
-write(num_rep_42,*) 'hwatstream      = ',hwatstream 
-write(num_rep_42,*) 'cf             = ',cf
-write(num_rep_42,*) 'betamax         = ', betamax
-write(num_rep_42,*) 'toblim         = ', toblim
-write(num_rep_42,*)'/'                            
-write(num_rep_42,428) '! hwatstream (m) :  critere de passage en stream en partant de la cote'
-write(num_rep_42,428) '!  si hwater > hwatstream '
-write(num_rep_42,428) '! cf coefficient de la loi de frottement fonction Neff'
-write(num_rep_42,428) '! seulement pour les points cotiers'
-write(num_rep_42,428) '! betamax : (Pa) frottement maxi sous les streams '
-write(num_rep_42,428) '! toblim : (Pa)  pour les iles '
-write(num_rep_42,*)
-
-inv_beta=0
-tostick=1.e5   ! valeurs pour les points non flgzmx
-tob_ile=betamax/2.
-moteurmax=toblim
-betamax_2d(:,:)=betamax
-
-!-------------------------------------------------------------------
-! masque stream
-
-      mstream_mx(:,:)=1
-      mstream_my(:,:)=1
-
-
-! coefficient permettant de modifier le basal drag.
-      drag_mx(:,:)=1.
-      drag_my(:,:)=1.
-
-where (socle_cry(:,:).lt.cry_lim)
-   debug_3D(:,:,25)=1
-elsewhere
-   debug_3D(:,:,25)=0
-endwhere
-
-return
-end subroutine init_dragging
-!________________________________________________________________________________
-
-!> SUBROUTINE: dragging
-!! Defini la localisation des streams et le frottement basal
-!>
-!-------------------------------------------------------------------------
-subroutine dragging   ! defini la localisation des streams et le frottement basal
-
-
-!         les iles n'ont pas de condition neff mais ont la condition toblim
-!         (ce bloc etait avant dans flottab)
-
-
-do j=2,ny
-   do i=2,nx
-      ilemx(i,j)=ilemx(i,j).and.(abs(rog*Hmx(i,j)*sdx(i,j)).lt.toblim)
-      ilemy(i,j)=ilemy(i,j).and.(abs(rog*Hmy(i,j)*sdy(i,j)).lt.toblim)
-   end do
-end do
-
-! le calcul des fleuves se fait sur les mailles majeures en partant de la cote
-
-
-
-
-! le gzmx initial a ete calculé dans flottab pour les points a cheval sur la grounding line
-
-! Partant de la cote, le fleuve est determine en remontant vers le point de hauteur d'eau max
-! tant que ce point est superieur a hwatermax.
-
-! Attention : ce systeme ne permet pas d'avoir des fleuves convergents
-!             et les fleuves n'ont une épaisseur que de 1 noeud.
-
-fleuvemx(:,:)=.false.
-fleuvemy(:,:)=.false.
-fleuve(:,:)=.false.
-cote(:,:)=.false.
-
-! detection des cotes
-do  j=2,ny-1  
-   do i=2,nx-1  
-      if ((.not.flot(i,j)).and. & 
-      ((flot(i+1,j)).or.(flot(i,j+1)).or.(flot(i-1,j)).or.(flot(i,j-1)))) then
-         cote(i,j)=.true.
-      endif
-   end do
-end do
-
-fleuve_maj: do j=2,ny-1  
-ifleuve:   do i=2,nx-1                     
-
-cote_detect :  if (cote(i,j)) then
-         idep=i
-         jdep=j
-
-!         if (hwater(i,j).gt.hwatstream) then
-         if (socle_cry(i,j).lt.0.) then    ! dans une vallee
-            fleuve(i,j)=.true. 
-!           write(166,*)
-!           write(166,*)'cote',i,j,hwater(i,j),hwatstream
-         else
-            cote(i,j)=.false.
-            cycle ifleuve
-         endif
-
-suit : do l=1,lmax         ! debut de la boucle de suivi, lmax longueur maxi des fleuves
-           i_moins1=max(idep-1,2)
-           j_moins1=max(jdep-1,1)
-           i_plus1=min(idep+1,nx)
-           j_plus1=min(jdep+1,ny)
-           
-
-
-!!$! recherche du max en suivant la hauteur d'eau maxi chez les voisins
-!!$! * en excluant les points flottants
-!!$! * et ceux qui sont deja tagges fleuves   !version 18 juillet 2007
-
-!!$           valmax=0.
-!!$
-!!$         do jloc=j_moins1,j_plus1
-!!$            do iloc=i_moins1,i_plus1
-!!$ 
-!!$               if ((hwater(iloc,jloc).gt.valmax) & 
-!!$                    .and.(.not.flot(iloc,jloc))  &
-!!$                    .and.(.not.fleuve(iloc,jloc))) then
-!!$                  imax=iloc
-!!$                  jmax=jloc
-!!$                  valmax=hwater(iloc,jloc)
-!!$               endif
-!!$            end do
-!!$         end do
-
-! recherche du max en suivant le socle le plus profond
-! * en excluant les points flottants
-! * et ceux qui sont deja tagges fleuves
-
-           valmax=1000.
-
-        do jloc=j_moins1,j_plus1
-            do iloc=i_moins1,i_plus1
- 
-               if ((B(iloc,jloc).lt.valmax)      & 
-                    .and.(.not.flot(iloc,jloc))  &
-                    .and.(.not.fleuve(iloc,jloc)).and.(socle_cry(iloc,jloc).lt.cry_lim)) then
-                  imax=iloc
-                  jmax=jloc
-                  valmax=B(iloc,jloc)
-               endif
-            end do
-         end do
-
-
-
-!         if (hwater(imax,jmax).gt.max(hwatstream-L*40.,200.)) then  ! version 18 juil
-         if ((hwater(imax,jmax).gt.hwatstream).and.(socle_cry(i,j).lt.cry_lim)) then
-            fleuve(imax,jmax)=.true.
-            idep=imax
-            jdep=jmax
-!            write(166,*) imax,jmax,hwater(imax,jmax)
-         else
-            fleuve(imax,jmax)=.false.
-            exit suit
-         end if
- 
-      end do suit
-
-   end if cote_detect
-
-end do ifleuve
-end do fleuve_maj
-
-
-! detection de problemes
-!write(166,*)
-!write(166,*) 'detection de pb dans fleuve   time=',time
-
-do j=2,ny
-   do i=2,nx
-      if (fleuve(i,j).and.(socle_cry(i,j).gt.cry_lim)) then
-         write(166,*) 'pb fleuve',time,i,j,hwater(i,j),socle_cry(i,j)
-      endif
-   end do
-end do
-
-do j=1,ny-1
-   do i=1,nx-1
-      if (fleuve(i,j)) then
-         fleuvemx(i,j)=.true.
-         fleuvemx(i+1,j)=.true.
-         fleuvemy(i,j)=.true.
-         fleuvemy(i,j+1)=.true.
-      end if
-   end do
-end do
-
-! pas de fleuve dans les endroits trop en pente
-
-fleuvemx(:,:)=fleuvemx(:,:).and.(abs(rog*Hmx(:,:)*sdx(:,:)).lt.toblim)
-fleuvemy(:,:)=fleuvemy(:,:).and.(abs(rog*Hmy(:,:)*sdy(:,:)).lt.toblim)
-
-do j=1,ny
-   do i=1,nx
-      if ((.not.ilemx(i,j)).and.(fleuvemx(i,j))) gzmx(i,j)=.true. 
-
-! calcul du frottement basal (ce bloc etait avant dans neffect)
-
-      if (cotemx(i,j)) then                        ! point cotier
-         betamx(i,j)=cf*neffmx(i,j) 
-         betamx(i,j)=min(betamx(i,j),betamax)
-      else if ((gzmx(i,j)).and.(.not.cotemx(i,j))) then  ! stream
-         betamx(i,j)=betamax
-         betamx(i,j)=cf*neffmx(i,j)*20.                  ! attention a cette modif juillet 2007
-         betamx(i,j)=min(betamx(i,j),betamax)
-         betamx(i,j)=max(betamx(i,j),10.)
-
-         if (cf*neffmx(i,j).gt.1500.) then
-            gzmx(i,j)=.false.
-            betamx(i,j)=tostick
-         endif
-
-
-      else if (ilemx(i,j)) then
-         betamx(i,j)=cf*neffmx(i,j) 
-         betamx(i,j)=min(betamx(i,j),tob_ile)       
-      else if (flotmx(i,j)) then     ! flottant ou ile
-         betamx(i,j)=0.
-      else                                         ! grounded, SIA
-         betamx(i,j)=tostick                       ! frottement glace posee (1 bar)
-      endif
-
-   end do
-end do
-
-
-do j=1,ny                  ! le noeud 1 n'est pas attribue
-   do i=1,nx   
-
-      if ((.not.ilemy(i,j)).and.(fleuvemy(i,j))) gzmy(i,j)=.true. 
-
-
- ! calcul du frottement basal (ce bloc etait avant dans neffect)
-
-      if (cotemy(i,j)) then                        ! point cotier
-         betamy(i,j)=cf*neffmy(i,j) 
-         betamy(i,j)=min(betamy(i,j),betamax)
-      else if ((gzmy(i,j)).and.(.not.cotemy(i,j))) then  ! stream
-         betamy(i,j)=betamax
-         betamy(i,j)=cf*neffmy(i,j)*20. 
-         betamy(i,j)=min(betamy(i,j),betamax)
-         betamy(i,j)=max(betamy(i,j),10.)
-
-         if (cf*neffmy(i,j).gt.1500.) then
-            gzmy(i,j)=.false.
-            betamy(i,j)=tostick
-         endif
-
-      else if (ilemy(i,j)) then
-         betamy(i,j)=cf*neffmy(i,j) 
-         betamy(i,j)=min(betamy(i,j),tob_ile)     
-      else if (flotmy(i,j)) then     ! flottant ou ile
-         betamy(i,j)=0.
-      else                               ! grounded, SIA
-         betamy(i,j)=tostick             ! frottement glace posee
-      endif
-   end do
-end do
-
-
-where (fleuve(:,:))
-   debug_3D(:,:,1)=1
-elsewhere (flot(:,:))
-   debug_3D(:,:,1)=2
-elsewhere
-   debug_3D(:,:,1)=0
-endwhere
-
-where (cote(:,:))
-   debug_3D(:,:,2)=1
-elsewhere
-   debug_3D(:,:,2)=0
-endwhere
-
-where (fleuvemx(:,:))
-   debug_3D(:,:,3)=1
-elsewhere
-   debug_3D(:,:,3)=0
-endwhere
-
-where (flotmx(:,:))
-   debug_3D(:,:,3)=-1
-endwhere
-
-!_____________________
-
-
-where (fleuvemy(:,:))
-   debug_3D(:,:,4)=1
-elsewhere
-   debug_3D(:,:,4)=0
-endwhere
-
-where (flotmy(:,:))
-   debug_3D(:,:,4)=-1
-end where
-
-debug_3D(:,:,23)= abs(RO*G*HMX(:,:)*sdx(:,:)*1.e-5) 
-debug_3D(:,:,24)= abs(RO*G*HMY(:,:)*sdy(:,:)*1.e-5) 
-
-return
-end subroutine dragging
-
-end module dragging_hwat_contmaj
-
Index: unk/SOURCES/dragging_hwat_sedim_mod.f90
===================================================================
--- /trunk/SOURCES/dragging_hwat_sedim_mod.f90	(revision 330)
+++ 	(revision )
@@ -1,477 +1,0 @@
-module dragging_hwat_sedim
-
-! Defini les zones de stream avec :
-! * un critere sur la hauteur d'eau 
-! * un critere sur la courbure du socle (si negatif, vallees)
-! * un critere sur l'epaisseur de sediments (carte de Gabi & Laske) 
-
-use module3d_phy
-use fake_beta_iter_vitbil_mod
-
-implicit none
-!logical,dimension(nx,ny) :: fleuvemx
-!logical,dimension(nx,ny) :: fleuvemy
-logical,dimension(nx,ny) :: fleuve
-logical,dimension(nx,ny) :: fleuve_sedim
-logical,dimension(nx,ny) :: fleuve_sedimx
-logical,dimension(nx,ny) :: fleuve_sedimy
-logical,dimension(nx,ny) :: cote
-
-
-real :: valmax
-integer :: imax,jmax
-integer :: i_moins1,i_plus1,j_moins1,j_plus1
-integer :: lmax=20
-integer :: idep,jdep,iloc,jloc
-
-
-real :: tostick   ! pour la glace posee
-real :: tob_ile    ! pour les iles
-real :: cry_lim=50.  ! courbure limite pour le suivi des fleuves
-
-real,dimension(nx,ny) :: h_sedim    ! epaisseur de sediment
-real,dimension(nx,ny) :: neff ! pression effective noeuds majeurs
-real :: seuil_hwater  ! seuil sur hwater pour avoir du glissement en zone sediment
-real :: seuil_sedim   ! seuil sur h_sedim pour avoir du glissement en zone sediment
-real :: seuil_neff    ! seuil sur la pression effective pour avoir glissement en zone sediment
-real :: coef_sedim    ! coef frottement zones sediments  (1)
-real :: coef_gz       ! coef frottement zones stream std  (10)
-real :: coef_ile      ! coef frottement zones iles   (0.1)
-real :: toto
-
-real :: betacotemax=100  ! valeur max pour betamx et betamy sur les cote
-real :: betacotemin=4
-real :: betasedmax=400    ! 400 dans c-HadCM3-evol-basal-melt-sed : cHadCM3s  et 250 dans c-HadCM3-evol-basal-melt-sed-betasedimin6
-real :: betasedmin=10   ! 10 dans c-HadCM3-evol-basal-melt-sed : cHadCM3s et 6 dans c-HadCM3-evol-basal-melt-sed-betasedimin6
-real :: betastreamax=500
-real :: betastreamin=10
-real :: betailemax=250
-real :: betailemin=6
-
-real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite avec spinup cat
-real, dimension(nx,ny) :: Vcol_y           !< uniquement pour compatibilite avec spinup cat
-real, dimension(nx,ny) :: Vsl_x            !< uniquement pour compatibilite avec spinup cat
-real, dimension(nx,ny) :: Vsl_y            !< uniquement pour compatibilite avec spinup cat
-logical :: corr_def = .false.              !< for deformation correction, pour compatibilite beta
-
-character(len=150) :: file_sedim   ! fichier epaisseur sediment
-
-contains
-!-------------------------------------------------------------------------------
-subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
-
-implicit none
- 
-namelist/drag_hwat_sedim/file_sedim,hwatstream,cf,betamax,toblim,seuil_hwater,seuil_sedim,seuil_neff,coef_sedim,coef_gz,coef_ile    
-
-if (itracebug.eq.1) write(num_tracebug,*)' dragging avec hwatermax'
-
-inv_beta=0
-
-! formats pour les ecritures dans 42
-428 format(A)
-
-! lecture des parametres du run                      block draghwat
-!--------------------------------------------------------------------
-
-rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
-read(num_param,drag_hwat_sedim)
-
-write(num_rep_42,428)'!___________________________________________________________' 
-write(num_rep_42,428) '&drag_hwat_sedim              ! nom du bloc hwat_sedim '
-write(num_rep_42,*)
-write(num_rep_42,*) 'file_sedim     = ', trim(file_sedim)
-write(num_rep_42,*) 'hwatstream     = ', hwatstream 
-write(num_rep_42,*) 'cf             = ', cf
-write(num_rep_42,*) 'betamax        = ', betamax
-write(num_rep_42,*) 'toblim         = ', toblim
-write(num_rep_42,*) 'seuil_hwater   = ', seuil_hwater
-write(num_rep_42,*) 'seuil_sedim    = ', seuil_sedim
-write(num_rep_42,*) 'seuil_neff     = ', seuil_neff
-write(num_rep_42,*) 'coef_sedim     = ', coef_sedim
-write(num_rep_42,*) 'coef_gz        = ', coef_gz
-write(num_rep_42,*) 'coef_ile       = ', coef_ile
-
-write(num_rep_42,*) '/'
-write(num_rep_42,428) '! file_sedim : fichier carte epaisseur sediment'                            
-write(num_rep_42,428) '! hwatstream (m) :  critere de passage en stream en partant de la cote si hwater > hwatstream '
-write(num_rep_42,428) '! cf coefficient de la loi de frottement fonction Neff'
-write(num_rep_42,428) '! seulement pour les points cotiers'
-write(num_rep_42,428) '! betamax : (Pa) frottement maxi sous les streams '
-write(num_rep_42,428) '! toblim : (Pa)  pour les iles '
-write(num_rep_42,428) '! seuil_hwater (m) : seuil hwater pour avoir glissement sur zone sediment'
-write(num_rep_42,428) '! seuil_sedim (m) : seuil epaisseur sediment pour avoir glissement sur zone sediment'
-write(num_rep_42,428) '! seuil_neff (Pa) seuil sur la pression effective pour avoir glissement en zone sediment'
-write(num_rep_42,428) '! coef_sedim : coef frottement zones stream sediments'
-write(num_rep_42,428) '! coef_gz : coef frottement zones stream std'
-write(num_rep_42,428) '! coef_ile : coef frottement zones iles'
-write(num_rep_42,*)
-
-tostick=1.e5   ! valeurs pour les points non flgzmx
-tob_ile=betamax/2.
-moteurmax=toblim
-
-!-------------------------------------------------------------------
-! masque stream
-
-      mstream_mx(:,:)=1
-      mstream_my(:,:)=1
-
-
-! coefficient permettant de modifier le basal drag.
-      drag_mx(:,:)=1.
-      drag_my(:,:)=1.
-
-where (socle_cry(:,:).lt.cry_lim)
-   debug_3D(:,:,25)=1
-elsewhere
-   debug_3D(:,:,25)=0
-endwhere
-
-! lecture du fichier sediment :
-!open(20,file=TRIM(DIRNAMEINP)//trim(file_sedim)) !   'sediment_ij_hemin40.dat')
-!do j=1,ny
-!   do i=1,nx
-!      read(20,*),toto,toto,h_sedim(i,j)
-!   enddo
-!enddo
-!h_sedim(:,:)=max(h_sedim(:,:),0.)  ! pas d'epaisseurs de sediments negatives
-h_sedim(:,:)=10. 
-  
-      return
-      end subroutine init_dragging
-!________________________________________________________________________________
-
-
-!-------------------------------------------------------------------------
-subroutine dragging   ! defini la localisation des streams et le frottement basal
-
-
-!         les iles n'ont pas de condition neff mais ont la condition toblim
-!         (ce bloc etait avant dans flottab)
-
-
-do j=2,ny
-   do i=2,nx
-      ilemx(i,j)=ilemx(i,j).and.(abs(rog*Hmx(i,j)*sdx(i,j)).lt.toblim)
-      ilemy(i,j)=ilemy(i,j).and.(abs(rog*Hmy(i,j)*sdy(i,j)).lt.toblim)
-   end do
-end do
-
-! le calcul des fleuves se fait sur les mailles majeures en partant de la cote
-
-
-
-
-! le gzmx initial a ete calculé dans flottab pour les points a cheval sur la grounding line
-
-! Partant de la cote, le fleuve est determine en remontant vers le point de hauteur d'eau max
-! tant que ce point est superieur a hwatermax.
-
-! Attention : ce systeme ne permet pas d'avoir des fleuves convergents
-!             et les fleuves n'ont une épaisseur que de 1 noeud.
-
-fleuvemx(:,:)=.false.
-fleuvemy(:,:)=.false.
-fleuve(:,:)=.false.
-cote(:,:)=.false.
-
-! detection des cotes
-do  j=2,ny-1  
-   do i=2,nx-1  
-      if ((.not.flot(i,j)).and. & 
-      ((flot(i+1,j)).or.(flot(i,j+1)).or.(flot(i-1,j)).or.(flot(i,j-1)))) then
-         cote(i,j)=.true.
-      endif
-   end do
-end do
-
-! calcul de neff (pression effective sur noeuds majeurs)
-do j=1,ny-1
-   do i=1,nx-1
-        neff(i,j)=(neffmx(i,j)+neffmx(i,j+1)+neffmy(i,j)+neffmy(i+1,j))/4
-   enddo
-enddo
-
-! test sur hauteur d'eau et epaisseur sediment :
-!where (h_sedim(:,:).GT.seuil_sedim.and.hwater(:,:).GT.seuil_hwater.and.(.not.flot(:,:))) fleuve(:,:)=.true.
-where (h_sedim(:,:).GT.seuil_sedim.and.neff(:,:).LE.seuil_neff.and.hwater(:,:).GT.seuil_hwater.and.(.not.flot(:,:))) fleuve(:,:)=.true.
-
-fleuve_sedim(:,:)=fleuve(:,:)  ! points uniquement fleuves sediment
-
-
-fleuve_maj: do j=2,ny-1  
-ifleuve:   do i=2,nx-1                     
-
-cote_detect :  if (cote(i,j)) then
-   idep=i
-   jdep=j
-         
-      if (socle_cry(i,j).lt.0.) then    ! dans une vallee
-         fleuve(i,j)=.true. 
-!         write(166,*)
-!         write(166,*)'cote',i,j,hwater(i,j),hwatstream
-      else
-         cote(i,j)=.false.
-         cycle ifleuve
-      endif
-
-suit : do l=1,lmax         ! debut de la boucle de suivi, lmax longueur maxi des fleuves
-           i_moins1=max(idep-1,2)
-           j_moins1=max(jdep-1,1)
-           i_plus1=min(idep+1,nx)
-           j_plus1=min(jdep+1,ny)
-           
-
-
-!!$! recherche du max en suivant la hauteur d'eau maxi chez les voisins
-!!$! * en excluant les points flottants
-!!$! * et ceux qui sont deja tagges fleuves   !version 18 juillet 2007
-
-!!$           valmax=0.
-!!$
-!!$         do jloc=j_moins1,j_plus1
-!!$            do iloc=i_moins1,i_plus1
-!!$ 
-!!$               if ((hwater(iloc,jloc).gt.valmax) & 
-!!$                    .and.(.not.flot(iloc,jloc))  &
-!!$                    .and.(.not.fleuve(iloc,jloc))) then
-!!$                  imax=iloc
-!!$                  jmax=jloc
-!!$                  valmax=hwater(iloc,jloc)
-!!$               endif
-!!$            end do
-!!$         end do
-
-! recherche du max en suivant le socle le plus profond
-! * en excluant les points flottants
-! * et ceux qui sont deja tagges fleuves
-
-           valmax=1000.
-
-        do jloc=j_moins1,j_plus1
-            do iloc=i_moins1,i_plus1
- 
-               if ((B(iloc,jloc).lt.valmax)      & 
-                    .and.(.not.flot(iloc,jloc))  &
-                    .and.(.not.fleuve(iloc,jloc)).and.(socle_cry(iloc,jloc).lt.cry_lim)) then
-                  imax=iloc
-                  jmax=jloc
-                 valmax=B(iloc,jloc)
-               endif
-            end do
-         end do
-
-
-
-
-         if ((hwater(imax,jmax).gt.hwatstream).and.(socle_cry(i,j).lt.cry_lim)) then
-            fleuve(imax,jmax)=.true.
-            idep=imax
-            jdep=jmax
-!            write(166,*) imax,jmax,hwater(imax,jmax)
-         else
-!            fleuve(imax,jmax)=.false.
-            exit suit
-         end if
- 
-      end do suit
-
-   end if cote_detect
-
-end do ifleuve
-end do fleuve_maj
-
-
-! detection de problemes
-!write(166,*)
-!write(166,*) 'detection de pb dans fleuve   time=',time
-!do j=2,ny
-!   do i=2,ny
-!      if (fleuve(i,j).and.(socle_cry(i,j).gt.cry_lim)) then
-!         write(166,*) 'pb fleuve',time,i,j,hwater(i,j),socle_cry(i,j)
-!      endif
-!   end do
-!end do
-
-do j=1,ny-1
-   do i=1,nx-1
-      if (fleuve(i,j)) then
-         fleuvemx(i,j)=.true.
-         fleuvemx(i+1,j)=.true.
-         fleuvemy(i,j)=.true.
-         fleuvemy(i,j+1)=.true.
-      end if
-      if (fleuve_sedim(i,j)) then
-         fleuve_sedimx(i,j)=.true.
-         fleuve_sedimx(i+1,j)=.true.
-         fleuve_sedimy(i,j)=.true.
-         fleuve_sedimy(i,j+1)=.true.
-      endif
-   end do
-end do
-
-! pas de fleuve dans les endroits trop en pente
-!fleuvemx(:,:)=fleuvemx(:,:).and.(abs(rog*Hmx(:,:)*sdx(:,:)).lt.toblim)
-!fleuvemy(:,:)=fleuvemy(:,:).and.(abs(rog*Hmy(:,:)*sdy(:,:)).lt.toblim)
-!fleuvemx(:,:)=fleuvemx(:,:).and.(abs(rog*Hmx(:,:)*sdx(:,:)).lt.toblim)
-!fleuvemy(:,:)=fleuvemy(:,:).and.(abs(rog*Hmy(:,:)*sdy(:,:)).lt.toblim)
-
-
-do j=1,ny
-   do i=1,nx
-      if ((.not.ilemx(i,j)).and.(fleuvemx(i,j))) gzmx(i,j)=.true. 
-
-! calcul du frottement basal (ce bloc etait avant dans neffect)
-
-      if (cotemx(i,j)) then                        ! point cotier
-         betamx(i,j)=cf*neffmx(i,j) 
-! cat         betamx(i,j)=min(betamx(i,j),tobmax)
-! version Jorge
-! point cotier: frott maxi=40, mini=4
-!         betamx(i,j)=min(betamx(i,j),betacotemax)
-!         betamx(i,j)=max(betamx(i,j),betacotemin)
-! version Jorge stream  en zone sediment :
-      else if (fleuve_sedimx(i,j).and.gzmx(i,j).and.(.not.cotemx(i,j))) then  ! stream  en zone sediment
-! stream sediment: frott maxi=60, mini=6 
-!cdc            betamx(i,j)=cf*neffmx(i,j)*3
-             betamx(i,j)=cf*neffmx(i,j)*coef_sedim
-!            betamx(i,j)=min(betamx(i,j),betasedmax)          
-!            betamx(i,j)=max(betamx(i,j),betasedmin)
-      else if ((gzmx(i,j)).and.(.not.cotemx(i,j))) then  ! stream normal
-! cat         betamx(i,j)=tobmax
-! cat         betamx(i,j)=cf*neffmx(i,j)*20.                  ! attention a cette modif juillet 2007
-! cat         betamx(i,j)=min(betamx(i,j),tobmax)
-! cat         betamx(i,j)=max(betamx(i,j),10.)
-! version Jorge
-! stream normal: frott maxi=80, mini=8
-         betamx(i,j)=cf*neffmx(i,j)*coef_gz         
-!         betamx(i,j)=min(betamx(i,j),betastreamax)           
-!         betamx(i,j)=max(betamx(i,j),betastreamin)
-
-! cat         if (cf*neffmx(i,j).gt.1500.) then
-! cat          gzmx(i,j)=.false.
-! cat           betamx(i,j)=tostick
-! cat        endif
-
-      else if (ilemx(i,j)) then
-         betamx(i,j)=cf*neffmx(i,j)*coef_ile
-! cat         betamx(i,j)=min(betamx(i,j),tob_ile)    
-! version Jorge
-! ile : frott maxi=20, mini=2   
-!         betamx(i,j)=min(betamx(i,j),betailemax)
-!         betamx(i,j)=max(betamx(i,j),betailemin)
-      else if (flotmx(i,j)) then     ! flottant ou ile
-         betamx(i,j)=0.
-      else                                         ! grounded, SIA
-         betamx(i,j)=tostick                       ! frottement glace posee (1 bar)
-      endif
-   end do
-end do
-
-
-do j=1,ny                  ! le noeud 1 n'est pas attribue
-   do i=1,nx   
-
-      if ((.not.ilemy(i,j)).and.(fleuvemy(i,j))) gzmy(i,j)=.true. 
-
-
- ! calcul du frottement basal (ce bloc etait avant dans neffect)
-
-      if (cotemy(i,j)) then                        ! point cotier
-         betamy(i,j)=cf*neffmy(i,j) 
-! cat         betamy(i,j)=min(betamy(i,j),tobmax)
-! version Jorge
-! point cotier: frott maxi=40, mini=4
-!         betamy(i,j)=min(betamy(i,j),betacotemax)
-!         betamy(i,j)=max(betamy(i,j),betacotemin)
-! version Jorge stream  en zone sediment :
-      else if ((fleuve_sedimy(i,j).and.gzmx(i,j).and.(.not.cotemx(i,j)))) then  ! stream  en zone sediment
-! stream sediment: frott maxi=60, mini=6 
-!cdc            betamy(i,j)=cf*neffmy(i,j)*3
-             betamy(i,j)=cf*neffmy(i,j)*coef_sedim
-!            betamy(i,j)=min(betamy(i,j),betasedmax)          
-!            betamy(i,j)=max(betamy(i,j),betasedmin)
-      else if ((gzmy(i,j)).and.(.not.cotemy(i,j))) then  ! stream normal
-! cat         betamy(i,j)=tobmax
-! cat         betamy(i,j)=cf*neffmy(i,j)*20. 
-! cat         betamy(i,j)=min(betamy(i,j),tobmax)
-! cat         betamy(i,j)=max(betamy(i,j),10.)
-! version Jorge
-! stream normal: frott maxi=80, mini=8
-         betamy(i,j)=cf*neffmy(i,j)*coef_gz                     
-!         betamy(i,j)=min(betamy(i,j),betastreamax)           
-!         betamy(i,j)=max(betamy(i,j),betastreamin)
-
-! cat         if (cf*neffmy(i,j).gt.1500.) then
-! cat            gzmy(i,j)=.false.
-! cat            betamy(i,j)=tostick
-! cat         endif
-
-      else if (ilemy(i,j)) then
-         betamy(i,j)=cf*neffmy(i,j)*coef_ile
-! cat         betamy(i,j)=min(betamy(i,j),tob_ile)     
-! version Jorge
-! ile : frott maxi=20, mini=2   
-!         betamy(i,j)=min(betamy(i,j),betailemax)
-!         betamy(i,j)=max(betamy(i,j),betailemin)
-      else if (flotmy(i,j)) then     ! flottant ou ile
-         betamy(i,j)=0.
-      else                               ! grounded, SIA
-         betamy(i,j)=tostick             ! frottement glace posee
-      endif
-   end do
-end do
-
-!print*,'dragging debug h_sedim(78,144),hwater,seuil_sedim,seuil_hwater',h_sedim(78,144),seuil_sedim,hwater(78,144),seuil_hwater,fleuve(78,144),gzmx(78,144)
-!print*,'dragging debug h_sedim(90,157),hwater,seuil_sedim,seuil_hwater',h_sedim(90,157),seuil_sedim,hwater(90,157),seuil_hwater,fleuve(90,157),gzmx(90,157),betamx(90,157),betamx(78,144)
-
-
-where (fleuve(:,:))
-   debug_3D(:,:,1)=1
-elsewhere (flot(:,:))
-   debug_3D(:,:,1)=2
-elsewhere
-   debug_3D(:,:,1)=0
-endwhere
-
-where (cote(:,:))
-   debug_3D(:,:,2)=1
-elsewhere
-   debug_3D(:,:,2)=0
-endwhere
-
-where (fleuvemx(:,:))
-   debug_3D(:,:,3)=1
-elsewhere
-   debug_3D(:,:,3)=0
-endwhere
-
-where (flotmx(:,:))
-   debug_3D(:,:,3)=-1
-endwhere
-
-!_____________________
-
-
-where (fleuvemy(:,:))
-   debug_3D(:,:,4)=1
-elsewhere
-   debug_3D(:,:,4)=0
-endwhere
-
-where (flotmy(:,:))
-   debug_3D(:,:,4)=-1
-end where
-
-debug_3D(:,:,23)= abs(RO*G*HMX(:,:)*sdx(:,:)*1.e-5) 
-debug_3D(:,:,24)= abs(RO*G*HMY(:,:)*sdy(:,:)*1.e-5) 
-
-
-!print*,'point (79,145)',hwater(79,145),h_sedim(79,145),flot(79,145),fleuve_sedim(79,145),fleuve(79,145),betamx(79,145),fleuve_sedimx(79,145),gzmx(79,145)
-
-return
-end subroutine dragging
-
-end module dragging_hwat_sedim
-
Index: unk/SOURCES/dragging_hwatermax_0.2_mod.f90
===================================================================
--- /trunk/SOURCES/dragging_hwatermax_0.2_mod.f90	(revision 330)
+++ 	(revision )
@@ -1,151 +1,0 @@
-!> \file dragging_hwatermax_0.2_mod.f90
-!! Defini les zones de stream avec un critere
-!<
-
-!> \namespace  dragging_hwatstream
-!! Defini les zones de stream avec 1 critere
-!! \author ...
-!! \date ...
-!! @note Le critere est:
-!! @note   * un critere sur la hauteur d'eau 
-!! @note Used module
-!! @note   - use module3D_phy
-!<
-module dragging_hwatstream
-  ! Defini les zones de stream avec un critere sur la hauteur d'eau
-
-  use module3D_phy
-  use fake_beta_iter_vitbil_mod
-
-
-  implicit none
-  real :: tostick
-
-contains
-  !-------------------------------------------------------------------------------
-  !> SUBROUTINE: init_dragging
-  !!Cette routine fait l'initialisation du dragging.
-  !>
-  subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
-
-    implicit none
-
-    namelist/drag_hwat/hwatstream,cf,betamax,toblim,tostick    
-
-    if (itracebug.eq.1)  call tracebug(' Dragging avec hwatermax')
-
-
-    ! formats pour les ecritures dans 42
-428 format(A)
-
-    ! lecture des parametres du run                      block draghwat
-    !--------------------------------------------------------------------
-
-    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
-    read(num_param,drag_hwat)
-
-    write(num_rep_42,428)'!___________________________________________________________' 
-    write(num_rep_42,428) '&drag_hwat              ! nom du bloc hwatermax '
-    write(num_rep_42,*)
-    write(num_rep_42,*) 'hwatsream      = ',hwatstream 
-    write(num_rep_42,*) 'cf             = ',cf
-    write(num_rep_42,*) 'betamax         = ', betamax
-    write(num_rep_42,*) 'toblim         = ', toblim
-    write(num_rep_42,*)'/'                            
-    write(num_rep_42,428) '! hwatstream (m) :  critere de passage en stream si hwater > hwatstream '
-    write(num_rep_42,428) '! cf coefficient de la loi de frottement fonction Neff'
-    write(num_rep_42,428) '! betamax : (Pa) frottement maxi '
-    write(num_rep_42,428) '! toblim : (Pa) tes pour les iles '
-    write(num_rep_42,*)
-
-    inv_beta=0
-    tostick=1.e5   ! valeurs pour les points non flgzmx
-
-    !-------------------------------------------------------------------
-    ! masque stream
-
-    mstream_mx(:,:)=1
-    mstream_my(:,:)=1
-
-
-    ! coefficient permettant de modifier le basal drag.
-    drag_mx(:,:)=1.
-    drag_my(:,:)=1.
-
-
-    return
-  end subroutine init_dragging
-  !________________________________________________________________________________
-
-  !> SUBROUTINE: dragging
-  !!Defini la localisation des streams et le frottement basal
-  !>
-  !-------------------------------------------------------------------------
-  subroutine dragging   ! defini la localisation des streams et le frottement basal
-
-
-    !         les iles n'ont pas de condition neff mais ont la condition toblim
-    !         (ce bloc etait avant dans flottab)
-
-    do j=2,ny
-       do i=2,nx
-          ilemx(i,j)=ilemx(i,j).and.(abs(tobmx(i,j)).lt.toblim)
-          ilemy(i,j)=ilemy(i,j).and.(abs(tobmy(i,j)).lt.toblim)
-       end do
-    end do
-
-
-    gzm_x: do j=1,ny  ! 2,ny
-       do i=2,nx ! -1
-
-          ! gzmx pour les hauteurs d'eau superieures a un seuil
-          ! attention, il semble qu'il y avait une erreur dans la precedente
-          ! formulation car le test etait seulement sur hwater(i,j).
-          ! le gzmx initial a ete calculé dans flottab pour les points a cheval sur la grounding line
-
-          gzmx(i,j)=gzmx(i,j).or. & 
-               (((hwater(i,j)+hwater(i-1,j))*0.5.gt.hwatstream) &
-               .and..not.flotmx(i,j))
-		  fleuvemx(i,j)=gzmx(i,j)
-
-          ! calcul du frottement basal (ce bloc etait avant dans neffect)
-          if (gzmx(i,j)) then                          ! stream
-             betamx(i,j)=cf*neffmx(i,j) 
-          else if (flotmx(i,j).or.ilemx(i,j)) then     ! flottant ou ile
-             betamx(i,j)=0.
-          else                                         ! grounded, SIA
-             betamx(i,j)=tostick                       ! frottement glace posee (1 bar)
-          endif
-
-       end do
-    end do gzm_x
-
-    gzm_y: do j=2,ny  !-1
-       do i=1,nx  ! 2,nx
-
-          ! gzmy pour les hauteurs d'eau superieures a un seuil
-          ! attention, il semble qu'il y avait une erreur dans la precedente
-          ! formulation car le test etait seulement sur hwater(i,j).
-
-          gzmy(i,j)=gzmy(i,j).or. &
-               (((hwater(i,j)+hwater(i,j-1))*0.5.gt.hwatstream) &
-               .and..not.flotmy(i,j))
-          fleuvemy(i,j)=gzmx(i,j)
-
-          ! calcul du frottement basal (ce bloc etait avant dans neffect)
-          if (gzmy(i,j)) then                ! stream
-             betamy(i,j)=cf*neffmy(i,j) 
-          else if (flotmy(i,j).or.ilemy(i,j)) then     ! flottant ou ile
-             betamy(i,j)=0.
-          else                               ! grounded, SIA
-             betamy(i,j)=tostick             ! frottement glace posee
-          endif
-
-       end do
-    end do gzm_y
-
-    return
-  end subroutine dragging
-
-end module dragging_hwatstream
-
Index: unk/SOURCES/dragging_neff_slope_mod.f90
===================================================================
--- /trunk/SOURCES/dragging_neff_slope_mod.f90	(revision 330)
+++ 	(revision )
@@ -1,422 +1,0 @@
-!> \file dragging_neff_contmaj_mod.f90
-!! Defini les zones de stream avec 3 criteres 
-!<
-
-!> \namespace  dragging_neff_contmaj
-!! Defini les zones de stream avec 3 criteres
-!! \author ...
-!! \date ...
-!! @note Les trois criteres sont:
-!! @note   * un critere sur la pression effective 
-!! @note   * un critere de continuite depuis la cote
-!! @note   * un critere sur la courbure du socle (si negatif, vallees) 
-!! @note Used module
-!! @note   - use module3D_phy
-!! @note   - use param_phy_mod
-!<
-
-module dragging_neff_slope
-
-! Defini les zones de stream avec :
-! * un critere sur la pression effective 
-! * un critere de continuite depuis la cote
-! * un critere sur la courbure du socle (si negatif, vallees) 
-
-use module3d_phy
-use param_phy_mod
-use interface_input
-use io_netcdf_grisli
-use fake_beta_iter_vitbil_mod
-
-implicit none
-
-logical,dimension(nx,ny) :: fleuve
-logical,dimension(nx,ny) :: cote
-logical,dimension(nx,ny) :: slowssamx ! not actual stream, but ssa used as Ub
-logical,dimension(nx,ny) :: slowssamy ! not actual stream, but ssa used as Ub
-logical,dimension(nx,ny) :: slowssa   ! not actual stream, but ssa used as Ub
-
-real,dimension(nx,ny) :: hires_slope ! slope comupted on a high resolution topography
-character(len=100) :: slope_fich     ! fichier grille
-character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat
-
-real :: valmax
-integer :: imax,jmax
-integer :: i_moins1,i_plus1,j_moins1,j_plus1
-integer :: lmax=20
-integer :: idep,jdep,iloc,jloc
-
-
-real :: tostick   ! pour la glace posee
-real :: betamin   ! betamin : (Pa) frottement mini sous les streams
-real :: tob_ile    ! pour les iles
-real :: cry_lim=50.  ! courbure limite pour le suivi des fleuves
-
-real,dimension(nx,ny) :: neff ! pression effective noeuds majeurs
-real :: seuil_neff    ! seuil sur la pression effective pour avoir glissement en zone sediment
-real :: coef_gz       ! coef frottement zones stream std  (10)
-real :: coef_ile      ! coef frottement zones iles   (0.1)
-
-real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite avec spinup cat
-real, dimension(nx,ny) :: Vcol_y           !< uniquement pour compatibilite avec spinup cat
-real, dimension(nx,ny) :: Vsl_x            !< uniquement pour compatibilite avec spinup cat
-real, dimension(nx,ny) :: Vsl_y            !< uniquement pour compatibilite avec spinup cat
-logical :: corr_def = .false.              !< for deformation correction, pour compatibilite beta
-
-contains
-!-------------------------------------------------------------------------------
-
-!> SUBROUTINE: init_dragging
-!! Cette routine fait l'initialisation du dragging.
-!>
-subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
-
-implicit none
-
-! local variables, defining the rugosity-enhanced flow
-real :: expo_slope
-real :: pente_min, pente_max
-
-namelist/drag_neff_slope/cf,betamax,betamin,toblim,tostick,seuil_neff,coef_gz,coef_ile,slope_fich,expo_slope,pente_min,pente_max
-
-if (itracebug.eq.1)  call tracebug(' dragging avec neff et slope')
-
-
-! formats pour les ecritures dans 42
-428 format(A)
-
-! lecture des parametres du run                      block drag neff
-!--------------------------------------------------------------------
-
-rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
-read(num_param,drag_neff_slope)
-
-write(num_rep_42,428)'!___________________________________________________________' 
-write(num_rep_42,428) '&drag_neff_slope          ! nom du bloc dragging neff slope'
-write(num_rep_42,*)
-write(num_rep_42,*) 'cf               = ',cf
-write(num_rep_42,*) 'betamax          = ', betamax
-write(num_rep_42,*) 'betamin          = ', betamin
-write(num_rep_42,*) 'toblim           = ', toblim
-write(num_rep_42,*) 'tostick          = ', tostick
-write(num_rep_42,*) 'seuil_neff       = ', seuil_neff
-write(num_rep_42,*) 'coef_gz          = ', coef_gz
-write(num_rep_42,*) 'coef_ile         = ', coef_ile
-write(num_rep_42,'(A,A)') 'slope_fich = ', slope_fich
-write(num_rep_42,*) 'expo_slope       = ', expo_slope
-write(num_rep_42,*) 'pente_min        = ', pente_min
-write(num_rep_42,*) 'pente_max        = ', pente_max
-write(num_rep_42,*)'/'                            
-write(num_rep_42,428) '! cf coefficient de la loi de frottement fonction Neff'
-write(num_rep_42,428) '! betamax : (Pa) frottement maxi sous les streams '
-write(num_rep_42,428) '! betamin : (Pa) frottement mini sous les streams '
-write(num_rep_42,428) '! toblim : (Pa)  pour les iles '
-write(num_rep_42,428) '! tostick : (Pa)  pour les points non flgzmx '
-write(num_rep_42,428) '! seuil_neff (Pa) seuil sur la pression effective pour avoir glissement'
-write(num_rep_42,428) '! coef_gz : coef frottement zones stream std'
-write(num_rep_42,428) '! coef_ile : coef frottement zones iles'
-write(num_rep_42,428) '! slope_fich : fichier pente bedrock'
-write(num_rep_42,428) '! expo_slope : exposant fonction pente'
-write(num_rep_42,428) '! pente_min : no flow enhancement below this'
-write(num_rep_42,428) '! pente_max : maximum flow enhancement above this'
-write(num_rep_42,*)
-
-inv_beta=0
-tob_ile=betamax/2.
-moteurmax=toblim
-
-! betamax_2d depends on sub-grid slopes.
-slope_fich=trim(dirnameinp)//trim(slope_fich)
-call lect_input(1,'slope',1,hires_slope(:,:),slope_fich,file_ncdf)
-
-! from slopes, we create an index between 0 & 1
-! 1 is very mountainous, 0 is flat
-hires_slope(:,:)=(max(min(hires_slope(:,:),pente_max),pente_min)-pente_min)/(pente_max-pente_min)
-! now we compute the actual betamax used by the remplimat routines
-! /!\ the slope is used to modify the beta where we have a temperate base,
-! but NO ice stream... -> Slow SSA zone (SSA used to compute Ub)
-betamax_2d(:,:)=max ( tostick * (1. - (1 - betamax / tostick ) * hires_slope(:,:)**(1./expo_slope)) , betamax )
-!do j=1,ny
-!   do i=1,nx
-!      write(18745,*) betamax_2d (i,j)
-!   enddo
-!enddo
-
-!-------------------------------------------------------------------
-! masque stream
-      mstream_mx(:,:)=1
-      mstream_my(:,:)=1
-
-
-! coefficient permettant de modifier le basal drag.
-      drag_mx(:,:)=1.
-      drag_my(:,:)=1.
-
-where (socle_cry(:,:).lt.cry_lim)
-   debug_3D(:,:,25)=1
-elsewhere
-   debug_3D(:,:,25)=0
-endwhere
-
-   
-      return
-      end subroutine init_dragging
-!________________________________________________________________________________
-
-!> SUBROUTINE: dragging
-!! Defini la localisation des streams et le frottement basal
-!>
-!-------------------------------------------------------------------------
-subroutine dragging   ! defini la localisation des streams et le frottement basal
-!$ USE OMP_LIB
-
-!         les iles n'ont pas de condition neff mais ont la condition toblim
-!         (ce bloc etait avant dans flottab)
-
-!$OMP PARALLEL
-!$OMP DO
-do j=2,ny
-   do i=2,nx
-      ilemx(i,j)=ilemx(i,j).and.(abs(rog*Hmx(i,j)*sdx(i,j)).lt.toblim)
-      ilemy(i,j)=ilemy(i,j).and.(abs(rog*Hmy(i,j)*sdy(i,j)).lt.toblim)
-   end do
-end do
-!$OMP END DO
-
-! le calcul des fleuves se fait sur les mailles majeures en partant de la cote
-
-
-
-
-! le gzmx initial a ete calculé dans flottab pour les points a cheval sur la grounding line
-
-! Partant de la cote, le fleuve est determine en remontant vers le point de hauteur d'eau max
-! tant que ce point est superieur a hwatermax.
-
-! Attention : ce systeme ne permet pas d'avoir des fleuves convergents
-!             et les fleuves n'ont une épaisseur que de 1 noeud.
-!$OMP WORKSHARE
-fleuvemx(:,:)=.false.
-fleuvemy(:,:)=.false.
-fleuve(:,:)=.false.
-cote(:,:)=.false.
-!$OMP END WORKSHARE
-
-! detection des cotes
-!$OMP DO
-do  j=2,ny-1  
-   do i=2,nx-1  
-      if ((.not.flot(i,j)).and. & 
-      ((flot(i+1,j)).or.(flot(i,j+1)).or.(flot(i-1,j)).or.(flot(i,j-1)))) then
-         cote(i,j)=.true.
-      endif
-   end do
-end do
-!$OMP END DO
-
-! calcul de neff (pression effective sur noeuds majeurs)
-!$OMP DO
-do j=1,ny-1
-   do i=1,nx-1
-        neff(i,j)=(neffmx(i,j)+neffmx(i+1,j)+neffmy(i,j)+neffmy(i,j+1))/4
-   enddo
-enddo
-!$OMP END DO
-!aurel, for the borders:
-!$OMP WORKSHARE
-neff(:,ny)=1.e5
-neff(nx,:)=1.e5
-
-
-! aurel, we add the neff threshold:
-where ((neff(:,:).le.seuil_neff).and.(.not.flot(:,:)).and.(H(:,:).gt.1.)) fleuve(:,:)=.true.
-!$OMP END WORKSHARE
-
-!$OMP DO
-do j=1,ny-1
-   do i=1,nx-1
-      if (fleuve(i,j)) then
-         fleuvemx(i,j)=.true.
-         fleuvemx(i+1,j)=.true.
-         fleuvemy(i,j)=.true.
-         fleuvemy(i,j+1)=.true.
-      end if
-   end do
-end do
-!$OMP END DO
-
-! we look for the warm based points that will not be treated as stream (ub from SSA):
-!$OMP WORKSHARE
-slowssa(:,:)=.false.
-slowssamx(:,:)=.false.
-slowssamy(:,:)=.false.
-!$OMP END WORKSHARE
-!$OMP DO
-do  j=1,ny
-   do i=1,nx
-      !if ((not(flot(i,j))).and.(hwater(i,j).gt.0.1)) slowssa(i,j)=.true.
-      if ((.not.(flot(i,j))).and.(ibase(i,j).ne.1).and.(H(i,j).gt.1.)) slowssa(i,j)=.true.
-   end do
-end do
-!$OMP END DO
-!$OMP DO
-do j=1,ny-1
-   do i=1,nx-1
-      if (slowssa(i,j)) then
-         slowssamx(i,j)=.true.
-         slowssamx(i+1,j)=.true.
-         slowssamy(i,j)=.true.
-         slowssamy(i,j+1)=.true.
-      end if
-   end do
-end do
-!$OMP END DO
-!$OMP DO
-do j=1,ny
-   do i=1,nx
-
-      ! the actual streams and the warm based points are gzmx now:
-      if ( ((.not.ilemx(i,j)).and.(fleuvemx(i,j))) .or. ((.not.ilemx(i,j)).and.(slowssamx(i,j))) ) gzmx(i,j)=.true. 
-
-      
-! calcul du frottement basal (ce bloc etait avant dans neffect)
-
-      if (cotemx(i,j)) then                        ! point cotier
-         betamx(i,j)=cf*neffmx(i,j) 
-         betamx(i,j)=min(betamx(i,j),betamax)
-
-      else if ((gzmx(i,j)).and.(.not.cotemx(i,j))) then ! tous les points SSA
-
-         if (fleuvemx(i,j)) then                 ! the actual streams
-            betamx(i,j)=cf*neffmx(i,j)*coef_gz
-            betamx(i,j)=min(betamx(i,j),betamax)
-            betamx(i,j)=max(betamx(i,j),betamin)
-
-            if (cf*neffmx(i,j).gt.betamax*2.) then  ! a stream that's becoming grounded...
-               if (slowssamx(i,j)) then
-                  betamx(i,j)=betamax_2d(i,j)
-               else
-                  gzmx(i,j)=.false.
-                  betamx(i,j)=tostick
-               endif
-            endif
-
-         else                    ! not an actual stream, SSA is used to compute Ub
-            betamx(i,j)=betamax_2d(i,j)
-         endif
-         
-      else if (ilemx(i,j)) then
-         betamx(i,j)=cf*neffmx(i,j)*coef_ile 
-         betamx(i,j)=min(betamx(i,j),tob_ile)       
-      else if (flotmx(i,j)) then     ! flottant ou ile
-         betamx(i,j)=0.
-      else                                         ! grounded, SIA
-         betamx(i,j)=tostick                       ! frottement glace posee (1 bar)
-      endif
-
-   end do
-end do
-!$OMP END DO
-
-!$OMP DO
-do j=1,ny
-   do i=1,nx
-
-      ! the actual streams and the warm based points are gzmx now:
-      if ( ((.not.ilemy(i,j)).and.(fleuvemy(i,j))) .or. ((.not.ilemy(i,j)).and.(slowssamy(i,j))) ) gzmy(i,j)=.true. 
-
-      
-! calcul du frottement basal (ce bloc etait avant dans neffect)
-
-      if (cotemy(i,j)) then                        ! point cotier
-         betamy(i,j)=cf*neffmy(i,j) 
-         betamy(i,j)=min(betamy(i,j),betamax)
-
-      else if ((gzmy(i,j)).and.(.not.cotemy(i,j))) then ! tous les points SSA
-
-         if (fleuvemy(i,j)) then                 ! the actual streams
-            betamy(i,j)=cf*neffmy(i,j)*coef_gz
-            betamy(i,j)=min(betamy(i,j),betamax)
-            betamy(i,j)=max(betamy(i,j),betamin)
-
-            if (cf*neffmy(i,j).gt.betamax*2.) then  ! a stream that's becoming grounded...
-               if (slowssamy(i,j)) then
-                  betamy(i,j)=betamax_2d(i,j)
-               else
-                  gzmy(i,j)=.false.
-                  betamy(i,j)=tostick
-               endif
-            endif
-
-         else                    ! not an actual stream, SSA is used to compute Ub
-            betamy(i,j)=betamax_2d(i,j)
-         endif
-         
-      else if (ilemy(i,j)) then
-         betamy(i,j)=cf*neffmy(i,j)*coef_ile 
-         betamy(i,j)=min(betamy(i,j),tob_ile)       
-      else if (flotmy(i,j)) then     ! flottant ou ile
-         betamy(i,j)=0.
-      else                                         ! grounded, SIA
-         betamy(i,j)=tostick                       ! frottement glace posee (1 bar)
-      endif
-
-   end do
-end do
-!$OMP END DO
-
-!$OMP DO
-do j=2,ny-1
-   do i=2,nx-1
-      beta_centre(i,j) = ((betamx(i,j)+betamx(i+1,j)) &
-          + (betamy(i,j)+betamy(i,j+1)))*0.25
-   end do
-end do
-!$OMP END DO
-!$OMP END PARALLEL
-
-!~ where (fleuve(:,:))
-!~    debug_3D(:,:,1)=1
-!~ elsewhere (flot(:,:))
-!~    debug_3D(:,:,1)=2
-!~ elsewhere
-!~    debug_3D(:,:,1)=0
-!~ endwhere
-!~ 
-!~ where (cote(:,:))
-!~    debug_3D(:,:,2)=1
-!~ elsewhere
-!~    debug_3D(:,:,2)=0
-!~ endwhere
-!~ 
-!~ where (fleuvemx(:,:))
-!~    debug_3D(:,:,3)=1
-!~ elsewhere
-!~    debug_3D(:,:,3)=0
-!~ endwhere
-!~ 
-!~ where (flotmx(:,:))
-!~    debug_3D(:,:,3)=-1
-!~ endwhere
-
-!_____________________
-
-
-!~ where (fleuvemy(:,:))
-!~    debug_3D(:,:,4)=1
-!~ elsewhere
-!~    debug_3D(:,:,4)=0
-!~ endwhere
-!~ 
-!~ where (flotmy(:,:))
-!~    debug_3D(:,:,4)=-1
-!~ end where
-!~ 
-!~ debug_3D(:,:,23)= abs(RO*G*HMX(:,:)*sdx(:,:)*1.e-5) 
-!~ debug_3D(:,:,24)= abs(RO*G*HMY(:,:)*sdy(:,:)*1.e-5) 
-
-return
-end subroutine dragging
-
-end module dragging_neff_slope
-
Index: unk/SOURCES/dragging_param_beta_nolin_mod.f90
===================================================================
--- /trunk/SOURCES/dragging_param_beta_nolin_mod.f90	(revision 330)
+++ 	(revision )
@@ -1,269 +1,0 @@
-!> \file dragging_neff_contmaj_mod.f90
-!! Defini les zones de stream avec 3 criteres 
-!<
-
-!> \namespace  dragging_neff_contmaj
-!! Defini les zones de stream avec 3 criteres
-!! \author ...
-!! \date ...
-!! @note Les trois criteres sont:
-!! @note   * un critere sur la pression effective 
-!! @note   * un critere de continuite depuis la cote
-!! @note   * un critere sur la courbure du socle (si negatif, vallees) 
-!! @note Used module
-!! @note   - use module3D_phy
-!! @note   - use param_phy_mod
-!<
-
-module dragging_param_beta_nolin
-
-! Defini les zones de stream avec :
-! * un critere sur la pression effective 
-! * un critere de continuite depuis la cote
-! * un critere sur la courbure du socle (si negatif, vallees) 
-
-use module3d_phy
-use param_phy_mod
-use interface_input
-use io_netcdf_grisli
-use fake_beta_iter_vitbil_mod
-
-implicit none
-
-logical,dimension(nx,ny) :: cote
-
-real,dimension(nx,ny) :: beta_param ! local variable
-
-real :: betamin   ! betamin : (Pa) frottement mini sous les streams
-
-real :: beta_slope        ! = A in:     B = A x Neff ** K
-real :: beta_expo         ! = K in:     B = A x Neff ** K
-
-real,dimension(nx,ny) :: neff   ! pression effective noeuds majeurs
-real,dimension(nx,ny) :: hwatmx ! hauteur d'eau staggered grille - afq
-real,dimension(nx,ny) :: hwatmy ! hauteur d'eau staggered grille - afq
-
-real,dimension(nx,ny) :: uslmag ! vitesse basale - afq
-
-real :: coef_ile      ! coef frottement zones iles   (0.1)
-
-real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite avec spinup cat
-real, dimension(nx,ny) :: Vcol_y           !< uniquement pour compatibilite avec spinup cat
-real, dimension(nx,ny) :: Vsl_x            !< uniquement pour compatibilite avec spinup cat
-real, dimension(nx,ny) :: Vsl_y            !< uniquement pour compatibilite avec spinup cat
-logical :: corr_def = .false.              !< for deformation correction, pour compatibilite beta
-
-contains
-!-------------------------------------------------------------------------------
-
-!> SUBROUTINE: init_dragging
-!! Cette routine fait l'initialisation du dragging.
-!>
-subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
-
-implicit none
-
-namelist/drag_param_beta/beta_slope,beta_expo,betamax,betamin,coef_ile
-
-if (itracebug.eq.1)  call tracebug(' dragging avec neff et slope')
-
-
-! formats pour les ecritures dans 42
-428 format(A)
-
-! lecture des parametres du run                      block drag neff
-!--------------------------------------------------------------------
-
-rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
-read(num_param,drag_param_beta)
-
-write(num_rep_42,428)'!___________________________________________________________' 
-write(num_rep_42,428) '&drag_param_beta_nolin    ! nom du bloc dragging param beta'
-write(num_rep_42,*)
-write(num_rep_42,*) 'beta_slope            = ', beta_slope
-write(num_rep_42,*) 'beta_expo             = ', beta_expo
-write(num_rep_42,*) 'betamax               = ', betamax
-write(num_rep_42,*) 'betamin               = ', betamin
-write(num_rep_42,*)'/'                            
-write(num_rep_42,428) '! Slope & expo of beta = (- slope x Neff ** expo ) / Uslmag'
-
-inv_beta=0
-!-------------------------------------------------------------------
-! masque stream
-mstream_mx(:,:)=1
-mstream_my(:,:)=1
-
-
-! coefficient permettant de modifier le basal drag.
-drag_mx(:,:)=1.
-drag_my(:,:)=1.
-      
-betamax_2d(:,:) = betamax
-
-! afq -- toblim = 0. !0.7e5 ! afq -- pour les iles, mais est-ce vraiment utile?
-  
-return
-end subroutine init_dragging
-!________________________________________________________________________________
-
-!> SUBROUTINE: dragging
-!! Defini la localisation des streams et le frottement basal
-!>
-!-------------------------------------------------------------------------
-subroutine dragging   ! defini la localisation des streams et le frottement basal
-!$ USE OMP_LIB
-
-!         les iles n'ont pas de condition neff mais ont la condition toblim
-!         (ce bloc etait avant dans flottab)
-
-!$OMP PARALLEL
-
-!    -- !$OMP DO
-!    -- do j=2,ny
-!    --    do i=2,nx
-!    --       ilemx(i,j)=ilemx(i,j).and.(abs(rog*Hmx(i,j)*sdx(i,j)).lt.toblim)
-!    --       ilemy(i,j)=ilemy(i,j).and.(abs(rog*Hmy(i,j)*sdy(i,j)).lt.toblim)
-!    --    end do
-!    -- end do
-!    -- !$OMP END DO
-
-
-! le gzmx initial a ete calculé dans flottab pour les points a cheval sur la grounding line
-
-
-!$OMP WORKSHARE
-fleuvemx(:,:)=.false.
-fleuvemy(:,:)=.false.
-cote(:,:)=.false.
-gzmx(:,:)=.true.
-gzmy(:,:)=.true.
-flgzmx(:,:)=.false.
-flgzmy(:,:)=.false.
-!$OMP END WORKSHARE
-
-! detection des cotes
-!$OMP DO
-do  j=2,ny-1  
-   do i=2,nx-1  
-      if ((.not.flot(i,j)).and. & 
-      ((flot(i+1,j)).or.(flot(i,j+1)).or.(flot(i-1,j)).or.(flot(i,j-1)))) then
-         cote(i,j)=.true.
-      endif
-   end do
-end do
-!$OMP END DO
-
-! calcul de neff (pression effective sur noeuds majeurs)
-if (sum(neffmx(:,:)).le.0.) neffmx(:,:) =1.e8
-if (sum(neffmy(:,:)).le.0.) neffmy(:,:) =1.e8
-
-!$OMP DO
-do j=1,ny-1
-   do i=1,nx-1
-        neff(i,j)=(neffmx(i,j)+neffmx(i+1,j)+neffmy(i,j)+neffmy(i,j+1))/4
-   enddo
-enddo
-!$OMP END DO
-!aurel, for the borders:
-!$OMP WORKSHARE
-neff(:,ny)=1.e5
-neff(nx,:)=1.e5
-! calcul de hwat (staggered)
-hwatmx(:,:)=0.
-hwatmy(:,:)=0.
-!$OMP END WORKSHARE
-!$OMP DO
-do i=2,nx
-   hwatmx(i,:)=(hwater(i-1,:)+hwater(i,:))/2.
-enddo
-!$OMP END DO
-!$OMP DO
-do j=2,ny
-   hwatmy(:,j)=(hwater(:,j-1)+hwater(:,j))/2.
-enddo
-!$OMP END DO
-
-!$OMP WORKSHARE
-
-! new parametrisation of beta on Neff:
-betamx(:,:)= beta_slope*(neffmx(:,:)**beta_expo)
-betamy(:,:)= beta_slope*(neffmy(:,:)**beta_expo)
-
-where (ilemx(:,:)) betamx(:,:) = betamx(:,:) * coef_ile
-where (ilemy(:,:)) betamy(:,:) = betamy(:,:) * coef_ile
-
-!uslmag(:,:) = (ux(:,:,nz)**2+uy(:,:,nz)**2)**(0.5)
-!where (uslmag(:,:).gt.0)
-!   betamx(:,:) = betamx(:,:) / uslmag(:,:)
-!   betamy(:,:) = betamy(:,:) / uslmag(:,:)
-!endwhere
-
-where (abs(ux(:,:,nz)).gt.0)
-   betamx(:,:) = betamx(:,:) / min(abs(ux(:,:,nz)),100.)
-endwhere
-where (abs(uy(:,:,nz)).gt.0)
-   betamy(:,:) = betamy(:,:) / min(abs(uy(:,:,nz)),100.)
-endwhere
-
-betamx(:,:)=max(betamx(:,:),betamin)
-betamy(:,:)=max(betamy(:,:),betamin)
-betamx(:,:)=min(betamx(:,:),betamax)
-betamy(:,:)=min(betamy(:,:),betamax)
-
-where ( hwatmx(:,:) .lt. 0.5 ) betamx(:,:) = betamax
-where ( hwatmy(:,:) .lt. 0.5 ) betamy(:,:) = betamax
-
-!$OMP END WORKSHARE
-
-
-! calcul de gzmx 
-
-! points flottants : flgzmx mais pas gzmx
-!$OMP DO
-do j=2,ny
-   do i=2,nx
-      if (flot(i,j).and.(flot(i-1,j))) then
-         flotmx(i,j)=.true.
-         flgzmx(i,j)=.true.
-         gzmx(i,j)=.false.
-         betamx(i,j)=0.
-      end if
-      if (flot(i,j).and.(flot(i,j-1))) then
-         flotmy(i,j)=.true.
-         flgzmy(i,j)=.true.
-         gzmy(i,j)=.false.
-         betamy(i,j)=0.
-      end if
-   end do
-end do
-!$OMP END DO
-
-!--------- autres criteres
-
-! points poses
-!$OMP WORKSHARE
-gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.betamax)   !  Pas de calcul pour les points
-gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.betamax)   !  au fort beta
-
-flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
-flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:)
-
-fleuvemx(:,:)=gzmx(:,:)
-fleuvemy(:,:)=gzmy(:,:)
-!$OMP END WORKSHARE
-
-!$OMP DO
-do j=2,ny-1
-   do i=2,nx-1
-      beta_centre(i,j) = ((betamx(i,j)+betamx(i+1,j)) &
-          + (betamy(i,j)+betamy(i,j+1)))*0.25
-   end do
-end do
-!$OMP END DO
-!$OMP END PARALLEL
-
-return
-end subroutine dragging
-
-end module dragging_param_beta_nolin
-
