Index: /trunk/SOURCES/Ant40_files/module_choix-antar40-0.4.f90
===================================================================
--- /trunk/SOURCES/Ant40_files/module_choix-antar40-0.4.f90	(revision 40)
+++ /trunk/SOURCES/Ant40_files/module_choix-antar40-0.4.f90	(revision 41)
@@ -87,4 +87,5 @@
 
 use dragging_hwat_contmaj ! frottement basal catherine grisliv8 ice stream peu actifs
+! use dragging_neff_contmaj ! dragging Aurelien et Christophe Fev 2016
 ! use dragging_hwat_sedim version sediment HN
 
Index: /trunk/SOURCES/Fichiers-parametres/anteis1_TEMPS-NETCDF.dat
===================================================================
--- /trunk/SOURCES/Fichiers-parametres/anteis1_TEMPS-NETCDF.dat	(revision 40)
+++ /trunk/SOURCES/Fichiers-parametres/anteis1_TEMPS-NETCDF.dat	(revision 41)
@@ -20,24 +20,12 @@
 1.e10
 -------------------------------------------------------------------------------------------
-15	npredeft        lecture de predef_tsort
+8	npredeft        lecture de predef_tsort
 10
 100
 200
 500
-100010
-100100
-200010
-200100
--20990
--20900
--20500
--19500
--121990
--125990
--125900
--419990
--419500
--419000
-
-
+1000
+2000
+5000
+7000
 -------------------------------------------------------------------------------------------
Index: /trunk/SOURCES/Makefile.grisli-gfortran.inc
===================================================================
--- /trunk/SOURCES/Makefile.grisli-gfortran.inc	(revision 40)
+++ /trunk/SOURCES/Makefile.grisli-gfortran.inc	(revision 41)
@@ -72,5 +72,5 @@
 	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_hwat_sedim_mod.o dragging_neff_contmaj_mod.o \
 	calving_frange.o no_calving.o no_lakes.o \
 	out_profile_mod.o printtable_mod.o mix-SIA-L1_mod.o \
@@ -700,4 +700,7 @@
 dragging_hwat_sedim_mod.o : dragging_hwat_sedim_mod.f90
 	$(FT) dragging_hwat_sedim_mod.f90
+
+dragging_neff_contmaj_mod.o : dragging_neff_contmaj_mod.f90
+	$(FT) dragging_neff_contmaj_mod.f90
 
 dragging_calc_beta_mod.o : Antarctique_general_files/dragging_calc_beta_mod.f90
Index: /trunk/SOURCES/Makefile.grisli.inc
===================================================================
--- /trunk/SOURCES/Makefile.grisli.inc	(revision 40)
+++ /trunk/SOURCES/Makefile.grisli.inc	(revision 41)
@@ -72,5 +72,5 @@
 	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_hwat_sedim_mod.o dragging_neff_contmaj_mod.o \
 	calving_frange.o no_calving.o no_lakes.o \
 	out_profile_mod.o printtable_mod.o mix-SIA-L1_mod.o \
@@ -704,4 +704,7 @@
 dragging_hwat_sedim_mod.o : dragging_hwat_sedim_mod.f90
 	$(FT) dragging_hwat_sedim_mod.f90
+
+dragging_neff_contmaj_mod.o : dragging_neff_contmaj_mod.f90
+	$(FT) dragging_neff_contmaj_mod.f90
 
 dragging_calc_beta_mod.o : Antarctique_general_files/dragging_calc_beta_mod.f90
Index: /trunk/SOURCES/dragging_hwat_sedim_mod.f90
===================================================================
--- /trunk/SOURCES/dragging_hwat_sedim_mod.f90	(revision 40)
+++ /trunk/SOURCES/dragging_hwat_sedim_mod.f90	(revision 41)
@@ -325,5 +325,5 @@
 
       if (cotemx(i,j)) then                        ! point cotier
-         betamx(i,j)=cf*neffmx(i,j)*0.1 
+         betamx(i,j)=cf*neffmx(i,j) 
 ! cat         betamx(i,j)=min(betamx(i,j),tobmax)
 ! version Jorge
@@ -379,5 +379,5 @@
 
       if (cotemy(i,j)) then                        ! point cotier
-         betamy(i,j)=cf*neffmy(i,j)*0.1 
+         betamy(i,j)=cf*neffmy(i,j) 
 ! cat         betamy(i,j)=min(betamy(i,j),tobmax)
 ! version Jorge
Index: /trunk/SOURCES/dragging_neff_contmaj_mod.f90
===================================================================
--- /trunk/SOURCES/dragging_neff_contmaj_mod.f90	(revision 41)
+++ /trunk/SOURCES/dragging_neff_contmaj_mod.f90	(revision 41)
@@ -0,0 +1,374 @@
+!> \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_contmaj
+
+! 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
+
+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 :: 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
+ 
+namelist/drag_neff_cont/hwatstream,cf,betamax,betamin,toblim,seuil_neff,coef_gz,coef_ile
+
+if (itracebug.eq.1)  call tracebug(' dragging avec hwatermax')
+
+
+! 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_cont)
+
+write(num_rep_42,428)'!___________________________________________________________' 
+write(num_rep_42,428) '&drag_neff_cont              ! nom du bloc dragging neff cont'
+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,*) 'betamin         = ', betamin
+write(num_rep_42,*) 'toblim          = ', toblim
+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,*)'/'                            
+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 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) '! 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,*)
+
+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
+
+   
+      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
+
+! 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+1,j)+neffmy(i,j)+neffmy(i,j+1))/4
+   enddo
+enddo
+!aurel, for the borders:
+neff(:,ny)=1.e5
+neff(nx,:)=1.e5
+
+
+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. 
+         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 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
+         else
+            fleuve(imax,jmax)=.false.
+            exit suit
+         end if
+ 
+      end do suit
+
+   end if cote_detect
+
+end do ifleuve
+end do fleuve_maj
+
+! aurel, we add the neff threshold:
+where ((neff(:,:).le.seuil_neff).and.(H(:,:).gt.50.).and.(.not.flot(:,:))) fleuve(:,:)=.true.
+
+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)=cf*neffmx(i,j)*coef_gz
+         betamx(i,j)=min(betamx(i,j),betamax)
+         betamx(i,j)=max(betamx(i,j),betamin)
+
+         !aurel, function of betamax and not 1500 as before:
+         if (cf*neffmx(i,j).gt.betamax*2.) then
+            gzmx(i,j)=.false.
+            betamx(i,j)=tostick
+         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
+
+
+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)=cf*neffmy(i,j)*coef_gz
+         betamy(i,j)=min(betamy(i,j),betamax)
+         betamy(i,j)=max(betamy(i,j),betamin)
+
+         !aurel, function of betamax and not 1500 as before:
+         if (cf*neffmy(i,j).gt.betamax*2.) then
+            gzmy(i,j)=.false.
+            betamy(i,j)=tostick
+         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
+      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_neff_contmaj
+
Index: /trunk/SOURCES/relaxation_water_diffusion.f90
===================================================================
--- /trunk/SOURCES/relaxation_water_diffusion.f90	(revision 40)
+++ /trunk/SOURCES/relaxation_water_diffusion.f90	(revision 41)
@@ -128,5 +128,5 @@
           erelax(i,j)=-kmy(i,j+1)*dtwdx2      ! drelax : diagonale i,j+1 
 
-          crelax(i,j)=1.+((kmx(i,j)+kmx(i+1,j))+(kmy(i,j+1)+kmy(i,j+1)))*dtwdx2 
+          crelax(i,j)=1.+((kmx(i,j)+kmx(i+1,j))+(kmy(i,j)+kmy(i,j+1)))*dtwdx2 
                                               !crelax : diagonale i,j 
 
