Index: branches/GRISLIv3/SOURCES/furst_schoof_mod.f90
===================================================================
--- branches/GRISLIv3/SOURCES/furst_schoof_mod.f90	(revision 364)
+++ branches/GRISLIv3/SOURCES/furst_schoof_mod.f90	(revision 397)
@@ -3,489 +3,500 @@
 module furst_schoof_mod
 
-use module3d_phy
-use deformation_mod_2lois
-
-implicit none
-
-real :: frot_coef ! afq: fixed in Tsai but = beta in Schoof (for m_weert = 1)
-real :: alpha_flot
-integer :: gr_select
-real,dimension(nx,ny) :: back_force_x
-real,dimension(nx,ny) :: back_force_y
-real, parameter :: m_weert = 1.0 ! afq: is this defined elsewhere in the model???
-real, parameter :: inv_mweert = 1./m_weert
+  use module3d_phy, only:nx,ny
+  use deformation_mod_2lois, only:
+
+  implicit none
+
+  real :: frot_coef ! afq: fixed in Tsai but = beta in Schoof (for m_weert = 1)
+  real :: alpha_flot
+  integer :: gr_select
+  real,dimension(nx,ny) :: back_force_x
+  real,dimension(nx,ny) :: back_force_y
+  real, parameter :: m_weert = 1.0 ! afq: is this defined elsewhere in the model???
+  real, parameter :: inv_mweert = 1./m_weert
 
 contains
 
-subroutine init_furst_schoof
-
-namelist/furst_schoof/frot_coef,gr_select
-
-rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
-read(num_param,furst_schoof)
-write(num_rep_42,'(A)')'!___________________________________________________________' 
-write(num_rep_42,'(A)') '&furst_schoof                                ! nom du bloc '
-write(num_rep_42,*)
-write(num_rep_42,*)'frot_coef  = ',frot_coef
-write(num_rep_42,*)'gr_select  = ',gr_select
-write(num_rep_42,*)'/'
-write(num_rep_42,*)'! frot_coef : solid friction law coef 0.6 in GMD 2018'
-write(num_rep_42,*)'! gr_select = 1 : Tsai , 2 : Schoof'
-
-!  frot = 0.035 ! f in the solid friction law, 0.035 in Ritz et al. 2001	
-
-  ! choix Tsai (loi de Coulomb) ou Schoof (loi de Weertman)
-  alpha_flot= ro/row
-  back_force_x(:,:)=0.1
-  back_force_y(:,:)=0.1
-  
-end subroutine init_furst_schoof  
-  
-subroutine interpol_glflux
-
-  implicit none
-  real :: xpos,ypos
-  real :: Hgl
-  !  integer :: igl
-
-  real :: phi_im2, phi_im1, phi_i, phi_ip1, phi_ip2, phi_ip3
-  real :: phi_jm2, phi_jm1, phi_j, phi_jp1, phi_jp2, phi_jp3
-  real,dimension(nx,ny) :: phi_xgl, phi_ygl ! flux gr line sur maille stag
-  integer,dimension(nx,ny) :: countx, county   ! how often do we modify ux/uy
-  real :: phi_prescr
-  real :: toutpetit = 1e-6
-  real :: denom, prodscal
-
-  real :: bfx, bfy
-
-  !debug
-  real,dimension(nx,ny) :: xpos_tab=9999.
-  real,dimension(nx,ny) :: ypos_tab=9999.
-  real,dimension(nx,ny) :: Hglx_tab=9999.
-  real,dimension(nx,ny) :: Hgly_tab=9999.
-  real,dimension(nx,ny) :: Hmxloc, Hmyloc ! pour eviter division par 0
-  real,dimension(nx,ny) :: phi_prescr_tabx, phi_prescr_taby
-
-  real,dimension(nx,ny) :: uxbartemp, uybartemp
-  
-  real,dimension(nx,ny) :: archimtab
-  
-  real,dimension(nx,ny) :: imx_diag_in
-
-  phi_prescr_tabx=0.
-  phi_prescr_taby=0.
-
-  Hmxloc(:,:)= max(Hmx(:,:),1.)
-  Hmyloc(:,:)=max(Hmy(:,:),1.)
-
-  uxbartemp(:,:)=0.
-  uybartemp(:,:)=0.
-  countx(:,:)=0.
-  county(:,:)=0.
-  gr_line_schoof(:,:)=0
-
-  archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot -sealevel_2d(:,:)
-
-  imx_diag_in(:,:) = imx_diag(:,:)
-  
-  ! detection des noeuds grounding line
-  do j= 3,ny-3
-     do i=3,nx-3
-        !archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel
-        if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel_2d(i,j))) then    ! grounded and with ice
-           !           print*,'schoof gr line ij',i,j
-           ! selon x (indice i)
-           
-           if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0..and..not.flot_marais(i-1,j)) then
-           !if (flot(i-1,j).and.Uxbar(i,j).LT.0.) then    ! grounding line between i-1,i    CAS WEST (ecoulement vers grille West)
-
-              call  calc_xgl4schoof(-dx,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i-1,j)),real(Bsoc(i-1,j)),sealevel_2d(i-1,j),xpos,Hgl)
-              xpos_tab(i,j)=xpos
-              Hglx_tab(i,j)=Hgl
-              
-              ! afq: the back force is on the staggered grid, the GL can be either West or East to this point.
-              if (xpos .lt. -dx/2.) then
-                 bfx = back_force_x(i,j)
-              else
-                 bfx = back_force_x(i+1,j)
-              endif
-              if (gr_select.eq.1) then ! flux de Tsai
-                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfx,phi_prescr)
-              else if (gr_select.eq.2) then ! flux de Schoof
-                 ! afq: the dragging coef. is on the staggered grid, the GL can be either West or East to this point.
-                 if (xpos .lt. -dx/2.) then
-                    frot_coef = betamx(i,j)
-                 else
-                    frot_coef = betamx(i+1,j)
-                 endif
-                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfx,phi_prescr)
-              else
-                 print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
-                 stop
-              endif
-              phi_prescr=-phi_prescr ! Doit Ãªtre negatif dans le cas West
-              phi_prescr_tabx(i,j)=phi_prescr
-
-              if (xpos .lt. -dx/2.) then  !  GL a west du i staggered,                          o centre, x stag
-
-                 !  grille    ! .....x......o......x......o...G..x......O......x......o......x......o
-                 !            !            i-2           i-1            i            i+1           i+2
-                 !flux        !     in            out         G out            in
-
-
-!afq --                 if (.not.(ilemx(i,j)).and..not.(ilemx(i-1,j))) then
-                    
-                 ! extrapolation pour avoir uxbar(i-1,j) qui sera ensuite prescrit dans call rempli_L2
-
-                 phi_im2 = Uxbar(i-2,j) * Hmxloc(i-2,j)
-                 phi_xgl(i-1,j)  = phi_im2 + dx * (phi_prescr-phi_im2)/(2.5 * dx + xpos)
-                 uxbartemp(i-1,j)   = uxbartemp(i-1,j)+ phi_xgl(i-1,j) / Hmxloc(i-1,j) ! attention division par 0 possible ?
-                 countx   (i-1,j)   = countx(i-1,j)+1
-                 imx_diag (i-1,j)   = 1
-                 gr_line_schoof(i-1,j) = 1
-
-                 ! extrapolation pour avoir uxbar(i,j)
-
-                 phi_ip1 = Uxbar(i+1,j) * Hmxloc(i+1,j)
-                 phi_xgl(i,j) = phi_ip1 + dx * (phi_prescr - phi_ip1)/(0.5 * dx - xpos)
-                 !                 print*,'uxbar(i,j)=',uxbar(i,j)
-                 uxbartemp(i,j)   = uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j)
-                 countx   (i,j)   = countx(i,j)+1
-                 imx_diag (i,j)   = 1
-                 gr_line_schoof(i,j) = 1
-                 !                 print*,'schoof uxbar(i,j)=',uxbar(i,j)
-                 !                 print*,'Hgl',Hgl
-                 !                 print*,'phi_prescr',phi_prescr, phi_prescr/Hgl
-                 !                 print*
-
-!afq --                 end if
-              
-              else            !  GL a Est du i staggered,                          o centre, x stag
-
-                 !  grille    ! .....x......o......x......o......x...G..O......x......o......x......o
-                 !            i-2           i-1            i            i+1           i+2
-                 !flux        !     in             in           out  G         out           in
-
-!afq --                 if (.not.(ilemx(i,j)).and..not.(ilemx(i+1,j))) then
-
-                 ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2
-
-                 phi_im1 = Uxbar(i-1,j) * Hmxloc(i-1,j)
-                 phi_xgl(i,j)  = phi_im1 + dx * (phi_prescr-phi_im1)/(1.5 * dx + xpos)
-                 uxbartemp(i,j)   = uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j) ! attention division par 0 possible ?
-                 countx   (i,j)   = countx(i,j)+1
-                 imx_diag (i,j) = 1
-                 gr_line_schoof(i,j) = 1
-
-                 ! extrapolation pour avoir uxbar(i+1,j)
-
-                 phi_ip2 = Uxbar(i+2,j) * Hmxloc(i+2,j)
-                 phi_xgl(i+1,j) = phi_ip2 + dx * (phi_prescr - phi_ip2)/(1.5 * dx  - xpos)
-                 uxbartemp(i+1,j)   = uxbartemp(i+1,j) +  phi_xgl(i+1,j) / Hmxloc(i+1,j)
-                 countx   (i+1,j)   = countx(i+1,j)+1
-                 imx_diag (i+1,j)   = 1
-                 gr_line_schoof(i+1,j) = 1
-
-!afq --                 end if
-
-              end if
-
-           else if (archimtab(i+1,j).LT.0..and.Uxbar(i+1,j).GT.0..and..not.flot_marais(i+1,j)) then
-           !else if (flot(i+1,j).and.Uxbar(i+1,j).GT.0.) then    ! grounding line between i,i+1    CAS EST (ecoulement vers grille Est)
-
-              call  calc_xgl4schoof(dx,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i+1,j)),real(Bsoc(i+1,j)),sealevel_2d(i+1,j),xpos,Hgl)
-              xpos_tab(i,j)=xpos
-              Hglx_tab(i,j)=Hgl
-              
-              ! afq: the back force is on the staggered grid, the GL can be either West or East to this point.
-              if (xpos .lt. dx/2.) then
-                 bfx = back_force_x(i,j)
-              else
-                 bfx = back_force_x(i+1,j)
-              endif
-              if (gr_select.eq.1) then ! flux de Tsai
-                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfx,phi_prescr)
-              else if (gr_select.eq.2) then ! flux de Schoof
-                 ! afq: the dragging coef. is on the staggered grid, the GL can be either West or East to this point.
-                 if (xpos .lt. dx/2.) then
-                    frot_coef = betamx(i,j)
-                 else
-                    frot_coef = betamx(i+1,j)
-                 endif
-                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfx,phi_prescr)
-              else
-                 print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
-                 stop
-              endif
-              phi_prescr_tabx(i,j)=phi_prescr
-              if (xpos .lt. dx/2.) then  !  GL a west du i staggered,                          o centre, x stag
-
-                 !  grille    ! .....x......o......x......o......x......O..G...x......o......x......o
-                 !            !            i-2           i-1            i            i+1           i+2
-                 !flux        !                    in           out        G  out            in
-
-!afq --                 if (.not.(ilemx(i,j)).and..not.(ilemx(i+1,j))) then
-                    
-                 ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2
-
-                 phi_im1 = Uxbar(i-1,j) * Hmxloc(i-1,j)
-                 phi_xgl(i,j)  = phi_im1 + dx * (phi_prescr-phi_im1)/(1.5 * dx + xpos)
-                 uxbartemp(i,j)   =  uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j) ! attention division par 0 possible ?
-                 countx   (i,j)   = countx(i,j)+1
-                 imx_diag (i,j) = 1
-                 gr_line_schoof(i,j) = 1
-
-                 ! extrapolation pour avoir uxbar(i+1,j)
-
-                 phi_ip2 = Uxbar(i+2,j) * Hmxloc(i+2,j)
-                 phi_xgl(i+1,j) = phi_ip2 + dx * (phi_prescr - phi_ip2)/(1.5* dx - xpos)
-                 uxbartemp(i+1,j)   = uxbartemp(i+1,j) + phi_xgl(i+1,j) / Hmxloc(i+1,j)
-                 countx   (i+1,j)   = countx(i+1,j)+1
-                 imx_diag (i+1,j) = 1
-                 gr_line_schoof(i+1,j) = 1
-
-!afq --                 end if
-              
-              else   !  GL a west du i staggered,                          o centre, x stag
-
-                 !  grille    ! ......x......o......x......O......x...G..o......x......o......x.....o
-                 !            !             i-1            i            i+1           i+2          i+3
-                 !flux        !                     in           out  G        out            in
-
-!afq --                 if (.not.(ilemx(i+1,j)).and..not.(ilemx(i+2,j))) then
-                    
-                 ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2
-
-                 phi_i = Uxbar(i,j) * Hmxloc(i,j)
-                 phi_xgl(i+1,j)  = phi_i + dx * (phi_prescr-phi_i)/(0.5 * dx + xpos)
-                 uxbartemp(i+1,j)   = uxbartemp(i+1,j) + phi_xgl(i+1,j) / Hmxloc(i+1,j) ! attention division par 0 possible ?
-                 countx   (i+1,j)   = countx(i+1,j)+1
-                 imx_diag (i+1,j) = 1
-                 gr_line_schoof(i+1,j) = 1
-
-                 ! extrapolation pour avoir uxbar(i+1,j)
-
-                 phi_ip3 = Uxbar(i+3,j) * Hmxloc(i+3,j)
-                 phi_xgl(i+2,j) = phi_ip3 + dx * (phi_prescr - phi_ip3)/(2.5 * dx  - xpos)
-                 uxbartemp(i+2,j)   = uxbartemp(i+2,j) + phi_xgl(i+2,j) / Hmxloc(i+2,j)
-                 countx   (i+2,j)   = countx(i+2,j)+1
-                 imx_diag (i+2,j) = 1
-                 gr_line_schoof(i+2,j) = 1
-
-!afq --                 end if
-              
-              end if
-
-           end if
-
-           !*******************************************************************************
-           ! selon y (indice j)
-           if (archimtab(i,j-1).LT.0..and.Uybar(i,j).LT.0..and..not.flot_marais(i,j-1)) then 
-           !if (flot(i,j-1).and.Uybar(i,j).LT.0.) then    ! grounding line between j-1,j    CAS SUD (ecoulement vers grille SUD)
-
-
-              call  calc_xgl4schoof(-dy,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i,j-1)),real(Bsoc(i,j-1)),sealevel_2d(i,j-1),ypos,Hgl)
-              ypos_tab(i,j)=ypos
-              Hgly_tab(i,j)=Hgl
-              
-              ! afq: the back force is on the staggered grid, the GL can be either South or North to this point.
-              if (ypos .lt. -dy/2.) then
-                 bfy = back_force_y(i,j)
-              else
-                 bfy = back_force_y(i,j+1)
-              endif
-              if (gr_select.eq.1) then ! flux de Tsai
-                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfy,phi_prescr)
-              else if (gr_select.eq.2) then ! flux de Schoof
-                 ! afq: the dragging coef. is on the staggered grid, the GL can be either South or North to this point.
-                 if (ypos .lt. -dy/2.) then
-                    frot_coef = betamy(i,j)
-                 else
-                    frot_coef = betamy(i,j+1)
-                 endif
-                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfy,phi_prescr)
-              else
-                 print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
-                 stop
-              endif
-              phi_prescr=-phi_prescr ! Doit Ãªtre negatif dans le cas Sud
-              phi_prescr_taby(i,j)=phi_prescr
-              if (ypos .lt. -dy/2.) then  !  GL au sud du j staggered,                          o centre, x stag
-
-                 !  grille    ! .....x......o......x......o...G..x......O......x......o......x......o
-                 !            !            j-2           j-1            j            j+1           j+2
-                 !flux        !     in            out         G out            in
-
-!afq --                 if (.not.(ilemy(i,j)).and..not.(ilemy(i,j-1))) then
-
-                 ! extrapolation pour avoir uybar(i,j-1) qui sera ensuite prescrit dans call rempli_L2
-
-                 phi_jm2 = Uybar(i,j-2) * Hmyloc(i,j-2)
-                 phi_ygl(i,j-1)  = phi_jm2 + dy * (phi_prescr-phi_jm2)/(2.5 * dy + ypos)
-                 uybartemp(i,j-1)   = uybartemp(i,j-1) + phi_ygl(i,j-1) / Hmyloc(i,j-1) ! attention division par 0 possible ?
-                 county   (i,j-1)   = county(i,j-1)+1
-                 imy_diag (i,j-1) = 1
-                 gr_line_schoof(i,j-1) = 1
-
-                 ! extrapolation pour avoir uybar(i,j)
-
-                 phi_jp1 = Uybar(i,j+1) * Hmyloc(i,j+1)
-                 phi_ygl(i,j) = phi_jp1 + dy * (phi_prescr - phi_jp1)/(0.5 * dy - ypos)
-                 uybartemp(i,j)   = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j)
-                 county   (i,j)   = county(i,j)+1
-                 imy_diag (i,j) = 1
-                 gr_line_schoof(i,j) = 1
-
-
-!afq --                 end if
-              
-              else            !  GL au nord du j staggered,                          o centre, x stag
-
-                 !  grille    ! .....x......o......x......o......x...G..O......x......o......x......o
-                 !            j-2           j-1            j            j+1           j+2
-                 !flux        !     in             in           out  G         out           in
-
-!afq --                 if (.not.(ilemy(i,j)).and..not.(ilemy(i,j+1))) then
-                    
-                 ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2
-
-                 phi_jm1 = Uybar(i,j-1) * Hmyloc(i,j-1)
-                 phi_ygl(i,j)  = phi_jm1 + dy * (phi_prescr-phi_jm1)/(1.5 * dy + ypos)
-                 uybartemp(i,j)   = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j) ! attention division par 0 possible ?
-                 county   (i,j)   = county(i,j)+1
-                 imy_diag(i,j) = 1
-                 gr_line_schoof(i,j) = 1
-
-                 ! extrapolation pour avoir uybar(i,j+1)
-
-                 phi_jp2 = Uybar(i,j+2) * Hmyloc(i,j+2)
-                 phi_ygl(i,j+1) = phi_jp2 + dy * (phi_prescr - phi_jp2)/(1.5 * dy  - ypos)
-                 uybartemp(i,j+1)   = uybartemp(i,j+1) + phi_ygl(i,j+1) / Hmyloc(i,j+1)
-                 county   (i,j+1)   = county(i,j+1)+1
-                 imy_diag (i,j+1) = 1
-                 gr_line_schoof(i,j+1) = 1
-
-!afq --                 end if
-              
-              end if
-              
-           else if (archimtab(i,j+1).LT.0..and.Uybar(i,j+1).GT.0..and..not.flot_marais(i,j+1)) then
-           !else if (flot(i,j+1).and.Uybar(i,j+1).GT.0.) then    ! grounding line between j,j+1    CAS NORD (ecoulement vers grille Nord)
-
-              call  calc_xgl4schoof(dy,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i,j+1)),real(Bsoc(i,j+1)),sealevel_2d(i,j+1),ypos,Hgl)
-              ypos_tab(i,j)=ypos
-              Hgly_tab(i,j)=Hgl
-              
-              ! afq: the back force is on the staggered grid, the GL can be either South or North to this point.
-              if (ypos .lt. dy/2.) then
-                 bfy = back_force_y(i,j)
-              else
-                 bfy = back_force_y(i,j+1)
-              endif
-              if (gr_select.eq.1) then ! flux de Tsai
-                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfy,phi_prescr)
-              else if (gr_select.eq.2) then ! flux de Schoof
-                 ! afq: the dragging coef. is on the staggered grid, the GL can be either South or North to this point.
-                 if (ypos .lt. dy/2.) then
-                    frot_coef = betamy(i,j)
-                 else
-                    frot_coef = betamy(i,j+1)
-                 endif
-                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfy,phi_prescr)
-              else
-                 print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
-                 stop
-              endif
-              phi_prescr_taby(i,j)=phi_prescr
-              if (ypos .lt. dy/2.) then  !  GL au sud du j staggered,                          o centre, x stag
-
-                 !  grille    ! .....x......o......x......o......x......O..G...x......o......x......o
-                 !            !            j-2           j-1            j            j+1           j+2
-                 !flux        !                    in           out        G  out            in
-
-!afq --                 if (.not.(ilemy(i,j)).and..not.(ilemy(i,j+1))) then
-
-                 ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2
-
-                 phi_jm1 = Uybar(i,j-1) * Hmyloc(i,j-1)
-                 phi_ygl(i,j)  = phi_jm1 + dy * (phi_prescr-phi_jm1)/(1.5 * dy + ypos)
-                 uybartemp(i,j)   = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j) ! attention division par 0 possible ?
-                 county   (i,j)   = county(i,j)+1
-                 imy_diag (i,j) = 1
-                 gr_line_schoof(i,j) = 1
-
-                 ! extrapolation pour avoir uybar(i,j+1)
-
-                 phi_jp2 = Uybar(i,j+2) * Hmyloc(i,j+2)
-                 phi_ygl(i,j+1) = phi_jp2 + dy * (phi_prescr - phi_jp2)/(1.5* dy - ypos)
-                 uybartemp(i,j+1)   = uybartemp(i,j+1) + phi_ygl(i,j+1) / Hmyloc(i,j+1)
-                 county   (i,j+1)   = county(i,j+1)+1
-                 imy_diag (i,j+1) = 1
-                 gr_line_schoof(i,j+1) = 1
-
-!afq --                 end if
-              
-              else   !  GL au nord du j staggered,                          o centre, x stag
-
-                 !  grille    ! ......x......o......x......O......x...G..o......x......o......x.....o
-                 !            !             j-1            j            j+1           j+2          j+3
-                 !flux        !                     in           out  G        out            in
-
-!afq --                 if (.not.(ilemy(i,j+1)).and..not.(ilemy(i,j+2))) then
-
-                 ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2
-
-                 phi_j = Uybar(i,j) * Hmyloc(i,j)
-                 phi_ygl(i,j+1)  = phi_j + dy * (phi_prescr-phi_j)/(0.5 * dy + ypos)
-                 uybartemp(i,j+1)   = uybartemp(i,j+1) +  phi_ygl(i,j+1) / Hmyloc(i,j+1) ! attention division par 0 possible ?
-                 county   (i,j+1)   = county(i,j+1)+1
-                 imy_diag (i,j+1) = 1
-                 gr_line_schoof(i,j+1) = 1
-
-                 ! extrapolation pour avoir uybar(i,j+1)
-
-                 phi_jp3 = Uybar(i,j+3) * Hmyloc(i,j+3)
-                 phi_ygl(i,j+2) = phi_jp3 + dy * (phi_prescr - phi_jp3)/(2.5 * dy  - ypos)
-                 uybartemp(i,j+2)   = uybartemp(i,j+2) +  phi_ygl(i,j+2) / Hmyloc(i,j+2)
-                 county   (i,j+2)   = county(i,j+2)+1
-                 imy_diag (i,j+2) = 1
-                 gr_line_schoof(i,j+2) = 1
-
-!afq --                 end if
-
-              end if
-
-           end if
-        end if
-     end do
-  end do
-
-  ! afq -- if we discard the points with multiple fluxes coming, uncomment:
-  where (countx(:,:).ge.2)
-     uxbartemp(:,:)=uxbar(:,:)*countx(:,:)
-     imx_diag(:,:) = imx_diag_in(:,:)
-  end where
-  where (county(:,:).ge.2)
-     uybartemp(:,:)=uybar(:,:)*county(:,:)
-     imy_diag(:,:) = imx_diag_in(:,:)
-  end where
-  ! afq --
-  
-  where (countx(:,:).ne.0) uxbar(:,:)= max(min(uxbartemp(:,:)/countx(:,:),V_limit),-V_limit)
-  where (county(:,:).ne.0) uybar(:,:)= max(min(uybartemp(:,:)/county(:,:),V_limit),-V_limit)
-           
-  do j= 3,ny-3
-     do i=3,nx-3
-        denom = sqrt((uxbartemp(i,j)**2 + uybartemp(i,j)**2))*sqrt((uxbar(i,j)**2 + uybar(i,j)**2))
-        if (denom.gt.toutpetit) then
-           ! prodscal is cos theta in : U . V = u v cos theta
-           ! this number is between -1 and 1, 1 = same direction & -1 = opposite direction
-           ! we could use this as an artificial back force... 
-           prodscal = (uxbartemp(i,j)*uxbar(i,j) + uybartemp(i,j)*uybar(i,j))/  &
-                (sqrt((uxbartemp(i,j)**2 + uybartemp(i,j)**2))*sqrt((uxbar(i,j)**2 + uybar(i,j)**2)))
-        endif
-     end do
-  end do
-  
+  subroutine init_furst_schoof
+    
+    use module3d_phy, only:num_param,num_rep_42
+    use param_phy_mod, only:row,ro
+
+    implicit none
+
+    namelist/furst_schoof/frot_coef,gr_select
+
+    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
+    read(num_param,furst_schoof)
+    write(num_rep_42,'(A)')'!___________________________________________________________' 
+    write(num_rep_42,'(A)') '&furst_schoof                                ! nom du bloc '
+    write(num_rep_42,*)
+    write(num_rep_42,*)'frot_coef  = ',frot_coef
+    write(num_rep_42,*)'gr_select  = ',gr_select
+    write(num_rep_42,*)'/'
+    write(num_rep_42,*)'! frot_coef : solid friction law coef 0.6 in GMD 2018'
+    write(num_rep_42,*)'! gr_select = 1 : Tsai , 2 : Schoof'
+
+    !  frot = 0.035 ! f in the solid friction law, 0.035 in Ritz et al. 2001	
+
+    ! choix Tsai (loi de Coulomb) ou Schoof (loi de Weertman)
+    alpha_flot= ro/row
+    back_force_x(:,:)=0.1
+    back_force_y(:,:)=0.1
+
+  end subroutine init_furst_schoof
+
+  subroutine interpol_glflux
+    
+    use module3d_phy, only:hmx,hmy,gr_line_schoof,Bsoc,H,sealevel_2d,imx_diag,imy_diag,&
+         uxbar,uybar,flot_marais,dx,dy,Abar,betamx,betamy,V_limit,debug_3D
+    use deformation_mod_2lois, only:glen
+    
+    implicit none
+
+    integer :: i,j
+    real :: xpos,ypos
+    real :: Hgl
+    !  integer :: igl
+
+    real :: phi_im2, phi_im1, phi_i, phi_ip1, phi_ip2, phi_ip3
+    real :: phi_jm2, phi_jm1, phi_j, phi_jp1, phi_jp2, phi_jp3
+    real,dimension(nx,ny) :: phi_xgl, phi_ygl ! flux gr line sur maille stag
+    integer,dimension(nx,ny) :: countx, county   ! how often do we modify ux/uy
+    real :: phi_prescr
+    real :: toutpetit = 1e-6
+    real :: denom, prodscal
+
+    real :: bfx, bfy
+
+    !debug
+    real,dimension(nx,ny) :: xpos_tab=9999.
+    real,dimension(nx,ny) :: ypos_tab=9999.
+    real,dimension(nx,ny) :: Hglx_tab=9999.
+    real,dimension(nx,ny) :: Hgly_tab=9999.
+    real,dimension(nx,ny) :: Hmxloc, Hmyloc ! pour eviter division par 0
+    real,dimension(nx,ny) :: phi_prescr_tabx, phi_prescr_taby
+
+    real,dimension(nx,ny) :: uxbartemp, uybartemp
+
+    real,dimension(nx,ny) :: archimtab
+
+    real,dimension(nx,ny) :: imx_diag_in
+
+    phi_prescr_tabx=0.
+    phi_prescr_taby=0.
+
+    Hmxloc(:,:)= max(Hmx(:,:),1.)
+    Hmyloc(:,:)=max(Hmy(:,:),1.)
+
+    uxbartemp(:,:)=0.
+    uybartemp(:,:)=0.
+    countx(:,:)=0.
+    county(:,:)=0.
+    gr_line_schoof(:,:)=0
+
+    archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot -sealevel_2d(:,:)
+
+    imx_diag_in(:,:) = imx_diag(:,:)
+
+    ! detection des noeuds grounding line
+    do j= 3,ny-3
+       do i=3,nx-3
+          !archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel
+          if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel_2d(i,j))) then    ! grounded and with ice
+             !           print*,'schoof gr line ij',i,j
+             ! selon x (indice i)
+
+             if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0..and..not.flot_marais(i-1,j)) then
+                !if (flot(i-1,j).and.Uxbar(i,j).LT.0.) then    ! grounding line between i-1,i    CAS WEST (ecoulement vers grille West)
+
+                call  calc_xgl4schoof(-dx,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i-1,j)),real(Bsoc(i-1,j)),sealevel_2d(i-1,j),xpos,Hgl)
+                xpos_tab(i,j)=xpos
+                Hglx_tab(i,j)=Hgl
+
+                ! afq: the back force is on the staggered grid, the GL can be either West or East to this point.
+                if (xpos .lt. -dx/2.) then
+                   bfx = back_force_x(i,j)
+                else
+                   bfx = back_force_x(i+1,j)
+                endif
+                if (gr_select.eq.1) then ! flux de Tsai
+                   call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfx,phi_prescr)
+                else if (gr_select.eq.2) then ! flux de Schoof
+                   ! afq: the dragging coef. is on the staggered grid, the GL can be either West or East to this point.
+                   if (xpos .lt. -dx/2.) then
+                      frot_coef = betamx(i,j)
+                   else
+                      frot_coef = betamx(i+1,j)
+                   endif
+                   call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfx,phi_prescr)
+                else
+                   print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
+                   stop
+                endif
+                phi_prescr=-phi_prescr ! Doit Ãªtre negatif dans le cas West
+                phi_prescr_tabx(i,j)=phi_prescr
+
+                if (xpos .lt. -dx/2.) then  !  GL a west du i staggered,                          o centre, x stag
+
+                   !  grille    ! .....x......o......x......o...G..x......O......x......o......x......o
+                   !            !            i-2           i-1            i            i+1           i+2
+                   !flux        !     in            out         G out            in
+
+
+                   !afq --                 if (.not.(ilemx(i,j)).and..not.(ilemx(i-1,j))) then
+
+                   ! extrapolation pour avoir uxbar(i-1,j) qui sera ensuite prescrit dans call rempli_L2
+
+                   phi_im2 = Uxbar(i-2,j) * Hmxloc(i-2,j)
+                   phi_xgl(i-1,j)  = phi_im2 + dx * (phi_prescr-phi_im2)/(2.5 * dx + xpos)
+                   uxbartemp(i-1,j)   = uxbartemp(i-1,j)+ phi_xgl(i-1,j) / Hmxloc(i-1,j) ! attention division par 0 possible ?
+                   countx   (i-1,j)   = countx(i-1,j)+1
+                   imx_diag (i-1,j)   = 1
+                   gr_line_schoof(i-1,j) = 1
+
+                   ! extrapolation pour avoir uxbar(i,j)
+
+                   phi_ip1 = Uxbar(i+1,j) * Hmxloc(i+1,j)
+                   phi_xgl(i,j) = phi_ip1 + dx * (phi_prescr - phi_ip1)/(0.5 * dx - xpos)
+                   !                 print*,'uxbar(i,j)=',uxbar(i,j)
+                   uxbartemp(i,j)   = uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j)
+                   countx   (i,j)   = countx(i,j)+1
+                   imx_diag (i,j)   = 1
+                   gr_line_schoof(i,j) = 1
+                   !                 print*,'schoof uxbar(i,j)=',uxbar(i,j)
+                   !                 print*,'Hgl',Hgl
+                   !                 print*,'phi_prescr',phi_prescr, phi_prescr/Hgl
+                   !                 print*
+
+                   !afq --                 end if
+
+                else            !  GL a Est du i staggered,                          o centre, x stag
+
+                   !  grille    ! .....x......o......x......o......x...G..O......x......o......x......o
+                   !            i-2           i-1            i            i+1           i+2
+                   !flux        !     in             in           out  G         out           in
+
+                   !afq --                 if (.not.(ilemx(i,j)).and..not.(ilemx(i+1,j))) then
+
+                   ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2
+
+                   phi_im1 = Uxbar(i-1,j) * Hmxloc(i-1,j)
+                   phi_xgl(i,j)  = phi_im1 + dx * (phi_prescr-phi_im1)/(1.5 * dx + xpos)
+                   uxbartemp(i,j)   = uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j) ! attention division par 0 possible ?
+                   countx   (i,j)   = countx(i,j)+1
+                   imx_diag (i,j) = 1
+                   gr_line_schoof(i,j) = 1
+
+                   ! extrapolation pour avoir uxbar(i+1,j)
+
+                   phi_ip2 = Uxbar(i+2,j) * Hmxloc(i+2,j)
+                   phi_xgl(i+1,j) = phi_ip2 + dx * (phi_prescr - phi_ip2)/(1.5 * dx  - xpos)
+                   uxbartemp(i+1,j)   = uxbartemp(i+1,j) +  phi_xgl(i+1,j) / Hmxloc(i+1,j)
+                   countx   (i+1,j)   = countx(i+1,j)+1
+                   imx_diag (i+1,j)   = 1
+                   gr_line_schoof(i+1,j) = 1
+
+                   !afq --                 end if
+
+                end if
+
+             else if (archimtab(i+1,j).LT.0..and.Uxbar(i+1,j).GT.0..and..not.flot_marais(i+1,j)) then
+                !else if (flot(i+1,j).and.Uxbar(i+1,j).GT.0.) then    ! grounding line between i,i+1    CAS EST (ecoulement vers grille Est)
+
+                call  calc_xgl4schoof(dx,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i+1,j)),real(Bsoc(i+1,j)),sealevel_2d(i+1,j),xpos,Hgl)
+                xpos_tab(i,j)=xpos
+                Hglx_tab(i,j)=Hgl
+
+                ! afq: the back force is on the staggered grid, the GL can be either West or East to this point.
+                if (xpos .lt. dx/2.) then
+                   bfx = back_force_x(i,j)
+                else
+                   bfx = back_force_x(i+1,j)
+                endif
+                if (gr_select.eq.1) then ! flux de Tsai
+                   call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfx,phi_prescr)
+                else if (gr_select.eq.2) then ! flux de Schoof
+                   ! afq: the dragging coef. is on the staggered grid, the GL can be either West or East to this point.
+                   if (xpos .lt. dx/2.) then
+                      frot_coef = betamx(i,j)
+                   else
+                      frot_coef = betamx(i+1,j)
+                   endif
+                   call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfx,phi_prescr)
+                else
+                   print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
+                   stop
+                endif
+                phi_prescr_tabx(i,j)=phi_prescr
+                if (xpos .lt. dx/2.) then  !  GL a west du i staggered,                          o centre, x stag
+
+                   !  grille    ! .....x......o......x......o......x......O..G...x......o......x......o
+                   !            !            i-2           i-1            i            i+1           i+2
+                   !flux        !                    in           out        G  out            in
+
+                   !afq --                 if (.not.(ilemx(i,j)).and..not.(ilemx(i+1,j))) then
+
+                   ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2
+
+                   phi_im1 = Uxbar(i-1,j) * Hmxloc(i-1,j)
+                   phi_xgl(i,j)  = phi_im1 + dx * (phi_prescr-phi_im1)/(1.5 * dx + xpos)
+                   uxbartemp(i,j)   =  uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j) ! attention division par 0 possible ?
+                   countx   (i,j)   = countx(i,j)+1
+                   imx_diag (i,j) = 1
+                   gr_line_schoof(i,j) = 1
+
+                   ! extrapolation pour avoir uxbar(i+1,j)
+
+                   phi_ip2 = Uxbar(i+2,j) * Hmxloc(i+2,j)
+                   phi_xgl(i+1,j) = phi_ip2 + dx * (phi_prescr - phi_ip2)/(1.5* dx - xpos)
+                   uxbartemp(i+1,j)   = uxbartemp(i+1,j) + phi_xgl(i+1,j) / Hmxloc(i+1,j)
+                   countx   (i+1,j)   = countx(i+1,j)+1
+                   imx_diag (i+1,j) = 1
+                   gr_line_schoof(i+1,j) = 1
+
+                   !afq --                 end if
+
+                else   !  GL a west du i staggered,                          o centre, x stag
+
+                   !  grille    ! ......x......o......x......O......x...G..o......x......o......x.....o
+                   !            !             i-1            i            i+1           i+2          i+3
+                   !flux        !                     in           out  G        out            in
+
+                   !afq --                 if (.not.(ilemx(i+1,j)).and..not.(ilemx(i+2,j))) then
+
+                   ! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2
+
+                   phi_i = Uxbar(i,j) * Hmxloc(i,j)
+                   phi_xgl(i+1,j)  = phi_i + dx * (phi_prescr-phi_i)/(0.5 * dx + xpos)
+                   uxbartemp(i+1,j)   = uxbartemp(i+1,j) + phi_xgl(i+1,j) / Hmxloc(i+1,j) ! attention division par 0 possible ?
+                   countx   (i+1,j)   = countx(i+1,j)+1
+                   imx_diag (i+1,j) = 1
+                   gr_line_schoof(i+1,j) = 1
+
+                   ! extrapolation pour avoir uxbar(i+1,j)
+
+                   phi_ip3 = Uxbar(i+3,j) * Hmxloc(i+3,j)
+                   phi_xgl(i+2,j) = phi_ip3 + dx * (phi_prescr - phi_ip3)/(2.5 * dx  - xpos)
+                   uxbartemp(i+2,j)   = uxbartemp(i+2,j) + phi_xgl(i+2,j) / Hmxloc(i+2,j)
+                   countx   (i+2,j)   = countx(i+2,j)+1
+                   imx_diag (i+2,j) = 1
+                   gr_line_schoof(i+2,j) = 1
+
+                   !afq --                 end if
+
+                end if
+
+             end if
+
+             !*******************************************************************************
+             ! selon y (indice j)
+             if (archimtab(i,j-1).LT.0..and.Uybar(i,j).LT.0..and..not.flot_marais(i,j-1)) then 
+                !if (flot(i,j-1).and.Uybar(i,j).LT.0.) then    ! grounding line between j-1,j    CAS SUD (ecoulement vers grille SUD)
+
+
+                call  calc_xgl4schoof(-dy,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i,j-1)),real(Bsoc(i,j-1)),sealevel_2d(i,j-1),ypos,Hgl)
+                ypos_tab(i,j)=ypos
+                Hgly_tab(i,j)=Hgl
+
+                ! afq: the back force is on the staggered grid, the GL can be either South or North to this point.
+                if (ypos .lt. -dy/2.) then
+                   bfy = back_force_y(i,j)
+                else
+                   bfy = back_force_y(i,j+1)
+                endif
+                if (gr_select.eq.1) then ! flux de Tsai
+                   call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfy,phi_prescr)
+                else if (gr_select.eq.2) then ! flux de Schoof
+                   ! afq: the dragging coef. is on the staggered grid, the GL can be either South or North to this point.
+                   if (ypos .lt. -dy/2.) then
+                      frot_coef = betamy(i,j)
+                   else
+                      frot_coef = betamy(i,j+1)
+                   endif
+                   call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfy,phi_prescr)
+                else
+                   print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
+                   stop
+                endif
+                phi_prescr=-phi_prescr ! Doit Ãªtre negatif dans le cas Sud
+                phi_prescr_taby(i,j)=phi_prescr
+                if (ypos .lt. -dy/2.) then  !  GL au sud du j staggered,                          o centre, x stag
+
+                   !  grille    ! .....x......o......x......o...G..x......O......x......o......x......o
+                   !            !            j-2           j-1            j            j+1           j+2
+                   !flux        !     in            out         G out            in
+
+                   !afq --                 if (.not.(ilemy(i,j)).and..not.(ilemy(i,j-1))) then
+
+                   ! extrapolation pour avoir uybar(i,j-1) qui sera ensuite prescrit dans call rempli_L2
+
+                   phi_jm2 = Uybar(i,j-2) * Hmyloc(i,j-2)
+                   phi_ygl(i,j-1)  = phi_jm2 + dy * (phi_prescr-phi_jm2)/(2.5 * dy + ypos)
+                   uybartemp(i,j-1)   = uybartemp(i,j-1) + phi_ygl(i,j-1) / Hmyloc(i,j-1) ! attention division par 0 possible ?
+                   county   (i,j-1)   = county(i,j-1)+1
+                   imy_diag (i,j-1) = 1
+                   gr_line_schoof(i,j-1) = 1
+
+                   ! extrapolation pour avoir uybar(i,j)
+
+                   phi_jp1 = Uybar(i,j+1) * Hmyloc(i,j+1)
+                   phi_ygl(i,j) = phi_jp1 + dy * (phi_prescr - phi_jp1)/(0.5 * dy - ypos)
+                   uybartemp(i,j)   = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j)
+                   county   (i,j)   = county(i,j)+1
+                   imy_diag (i,j) = 1
+                   gr_line_schoof(i,j) = 1
+
+
+                   !afq --                 end if
+
+                else            !  GL au nord du j staggered,                          o centre, x stag
+
+                   !  grille    ! .....x......o......x......o......x...G..O......x......o......x......o
+                   !            j-2           j-1            j            j+1           j+2
+                   !flux        !     in             in           out  G         out           in
+
+                   !afq --                 if (.not.(ilemy(i,j)).and..not.(ilemy(i,j+1))) then
+
+                   ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2
+
+                   phi_jm1 = Uybar(i,j-1) * Hmyloc(i,j-1)
+                   phi_ygl(i,j)  = phi_jm1 + dy * (phi_prescr-phi_jm1)/(1.5 * dy + ypos)
+                   uybartemp(i,j)   = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j) ! attention division par 0 possible ?
+                   county   (i,j)   = county(i,j)+1
+                   imy_diag(i,j) = 1
+                   gr_line_schoof(i,j) = 1
+
+                   ! extrapolation pour avoir uybar(i,j+1)
+
+                   phi_jp2 = Uybar(i,j+2) * Hmyloc(i,j+2)
+                   phi_ygl(i,j+1) = phi_jp2 + dy * (phi_prescr - phi_jp2)/(1.5 * dy  - ypos)
+                   uybartemp(i,j+1)   = uybartemp(i,j+1) + phi_ygl(i,j+1) / Hmyloc(i,j+1)
+                   county   (i,j+1)   = county(i,j+1)+1
+                   imy_diag (i,j+1) = 1
+                   gr_line_schoof(i,j+1) = 1
+
+                   !afq --                 end if
+
+                end if
+
+             else if (archimtab(i,j+1).LT.0..and.Uybar(i,j+1).GT.0..and..not.flot_marais(i,j+1)) then
+                !else if (flot(i,j+1).and.Uybar(i,j+1).GT.0.) then    ! grounding line between j,j+1    CAS NORD (ecoulement vers grille Nord)
+
+                call  calc_xgl4schoof(dy,alpha_flot,real(H(i,j)),real(Bsoc(i,j)),sealevel_2d(i,j),real(H(i,j+1)),real(Bsoc(i,j+1)),sealevel_2d(i,j+1),ypos,Hgl)
+                ypos_tab(i,j)=ypos
+                Hgly_tab(i,j)=Hgl
+
+                ! afq: the back force is on the staggered grid, the GL can be either South or North to this point.
+                if (ypos .lt. dy/2.) then
+                   bfy = back_force_y(i,j)
+                else
+                   bfy = back_force_y(i,j+1)
+                endif
+                if (gr_select.eq.1) then ! flux de Tsai
+                   call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfy,phi_prescr)
+                else if (gr_select.eq.2) then ! flux de Schoof
+                   ! afq: the dragging coef. is on the staggered grid, the GL can be either South or North to this point.
+                   if (ypos .lt. dy/2.) then
+                      frot_coef = betamy(i,j)
+                   else
+                      frot_coef = betamy(i,j+1)
+                   endif
+                   call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfy,phi_prescr)
+                else
+                   print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
+                   stop
+                endif
+                phi_prescr_taby(i,j)=phi_prescr
+                if (ypos .lt. dy/2.) then  !  GL au sud du j staggered,                          o centre, x stag
+
+                   !  grille    ! .....x......o......x......o......x......O..G...x......o......x......o
+                   !            !            j-2           j-1            j            j+1           j+2
+                   !flux        !                    in           out        G  out            in
+
+                   !afq --                 if (.not.(ilemy(i,j)).and..not.(ilemy(i,j+1))) then
+
+                   ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2
+
+                   phi_jm1 = Uybar(i,j-1) * Hmyloc(i,j-1)
+                   phi_ygl(i,j)  = phi_jm1 + dy * (phi_prescr-phi_jm1)/(1.5 * dy + ypos)
+                   uybartemp(i,j)   = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j) ! attention division par 0 possible ?
+                   county   (i,j)   = county(i,j)+1
+                   imy_diag (i,j) = 1
+                   gr_line_schoof(i,j) = 1
+
+                   ! extrapolation pour avoir uybar(i,j+1)
+
+                   phi_jp2 = Uybar(i,j+2) * Hmyloc(i,j+2)
+                   phi_ygl(i,j+1) = phi_jp2 + dy * (phi_prescr - phi_jp2)/(1.5* dy - ypos)
+                   uybartemp(i,j+1)   = uybartemp(i,j+1) + phi_ygl(i,j+1) / Hmyloc(i,j+1)
+                   county   (i,j+1)   = county(i,j+1)+1
+                   imy_diag (i,j+1) = 1
+                   gr_line_schoof(i,j+1) = 1
+
+                   !afq --                 end if
+
+                else   !  GL au nord du j staggered,                          o centre, x stag
+
+                   !  grille    ! ......x......o......x......O......x...G..o......x......o......x.....o
+                   !            !             j-1            j            j+1           j+2          j+3
+                   !flux        !                     in           out  G        out            in
+
+                   !afq --                 if (.not.(ilemy(i,j+1)).and..not.(ilemy(i,j+2))) then
+
+                   ! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2
+
+                   phi_j = Uybar(i,j) * Hmyloc(i,j)
+                   phi_ygl(i,j+1)  = phi_j + dy * (phi_prescr-phi_j)/(0.5 * dy + ypos)
+                   uybartemp(i,j+1)   = uybartemp(i,j+1) +  phi_ygl(i,j+1) / Hmyloc(i,j+1) ! attention division par 0 possible ?
+                   county   (i,j+1)   = county(i,j+1)+1
+                   imy_diag (i,j+1) = 1
+                   gr_line_schoof(i,j+1) = 1
+
+                   ! extrapolation pour avoir uybar(i,j+1)
+
+                   phi_jp3 = Uybar(i,j+3) * Hmyloc(i,j+3)
+                   phi_ygl(i,j+2) = phi_jp3 + dy * (phi_prescr - phi_jp3)/(2.5 * dy  - ypos)
+                   uybartemp(i,j+2)   = uybartemp(i,j+2) +  phi_ygl(i,j+2) / Hmyloc(i,j+2)
+                   county   (i,j+2)   = county(i,j+2)+1
+                   imy_diag (i,j+2) = 1
+                   gr_line_schoof(i,j+2) = 1
+
+                   !afq --                 end if
+
+                end if
+
+             end if
+          end if
+       end do
+    end do
+
+    ! afq -- if we discard the points with multiple fluxes coming, uncomment:
+    where (countx(:,:).ge.2)
+       uxbartemp(:,:)=uxbar(:,:)*countx(:,:)
+       imx_diag(:,:) = imx_diag_in(:,:)
+    end where
+    where (county(:,:).ge.2)
+       uybartemp(:,:)=uybar(:,:)*county(:,:)
+       imy_diag(:,:) = imx_diag_in(:,:)
+    end where
+    ! afq --
+
+    where (countx(:,:).ne.0) uxbar(:,:)= max(min(uxbartemp(:,:)/countx(:,:),V_limit),-V_limit)
+    where (county(:,:).ne.0) uybar(:,:)= max(min(uybartemp(:,:)/county(:,:),V_limit),-V_limit)
+
+    do j= 3,ny-3
+       do i=3,nx-3
+          denom = sqrt((uxbartemp(i,j)**2 + uybartemp(i,j)**2))*sqrt((uxbar(i,j)**2 + uybar(i,j)**2))
+          if (denom.gt.toutpetit) then
+             ! prodscal is cos theta in : U . V = u v cos theta
+             ! this number is between -1 and 1, 1 = same direction & -1 = opposite direction
+             ! we could use this as an artificial back force... 
+             prodscal = (uxbartemp(i,j)*uxbar(i,j) + uybartemp(i,j)*uybar(i,j))/  &
+                  (sqrt((uxbartemp(i,j)**2 + uybartemp(i,j)**2))*sqrt((uxbar(i,j)**2 + uybar(i,j)**2)))
+          endif
+       end do
+    end do
+
 !!$  do j=1,ny
 !!$     do i=1,nx
@@ -500,132 +511,135 @@
 !!$  enddo
 
-  debug_3D(:,:,66) = phi_prescr_tabx(:,:)
-  debug_3D(:,:,67) = phi_prescr_taby(:,:)
-  
-end subroutine interpol_glflux
-  
-!-------------------------------------------------------------------------
-! calcule la position sous maille de la ligne dechouage
- subroutine calc_xgl4schoof(dy,alpha,H_0,B_0,SL_0,H_1,B_1,SL_1,xpos,Hgl)
-                
-
-  implicit none
-
-  real,intent(in)        ::    dy        !< longueur orientee de la maille 
-  real,intent(in)        ::    alpha     !< coefficient de flottaison = rho/rhow
-  real,intent(in)        ::    H_0       !< epaisseur au point pose
-  real,intent(in)        ::    B_0       !< altitude socle au point pose
-  real,intent(in)        ::    SL_0      !< sea level au point pose
-  real,intent(in)        ::    H_1       !< epaisseur au point flottant
-  real,intent(in)        ::    B_1       !< altitude socle au point flottant
-  real,intent(in)        ::    SL_1      !< sea level au point flottant
-  real,intent(out)       ::    xpos      !< position de la ligne (en distance depuis le point pose)
-  real,intent(out)       ::    Hgl       !< epaisseur de glace a la grounding line
-   
-  
-
-  real                   ::    Cgl       ! variable de travail  
-
-!~   Cgl   = (-alpha * (H_1-H_0) - (B_1-B_0)) ! -(rho/rhow (H1-H0) + (B1-B0))
-
-!~   if (abs(Cgl).gt.1.e-5) then   
-!~      xpos    = (B_0 + alpha * H_0 - SL) / Cgl 
-!~   else
-!~      xpos    = 1. - 1./abs(dy)   ! verifier
-!~   end if
-  
-!~   if (xpos.LT.0..OR.xpos.GT.1.) then
-!~ 		print*,'calc_xgl4schoof : xpos value error i,j=',i,j
-!~ 		print*, 'xpos,Cgl', xpos,Cgl
-!~ 		print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL
-!~ 		stop
-!~ 	endif
-		
-  Cgl   = (B_1-B_0) + alpha * (H_1-H_0) - (SL_1-SL_0)
-
-  if (abs(Cgl).gt.1.e-5) then   
-     xpos    = (SL_0 - (B_0 + alpha * H_0)) / Cgl 
-  else
-     xpos    = 0.5  ! verifier
-  end if
-  
-!  print*
-!  print*, 'xpos', xpos
-!	print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL
-!	print*,'B_1, H_1', B_1, H_1
-
-	
-  if (xpos.LT.0..OR.xpos.GT.1.) then
-     print*,'calc_xgl4schoof : xpos value error i,j=',i,j
-     print*, 'xpos,Cgl', xpos,Cgl
-     print*,'B_0, alpha, H_0, SL_0', B_0, alpha, H_0, SL_0
-     print*,'archim=',B_1+H_1*alpha - SL_1
-     !stop
-     xpos = min(max(0.,xpos),1.)
-  endif
-	
-  xpos = xpos * dy
-  
-  Hgl= H_0 + (H_1 - H_0) * xpos/dy
-
-end subroutine calc_xgl4schoof
-
-
-!--------------------------------------------------------------------
-! Schoof_flux : calcule le coefficient et le flux de Schoof 
-!-------------------------------------------------------------------- 
-
-! Dans un premier temps, pas de terme horizontal (regarder dans toy_recul une tentative de confinement)
-
-subroutine flux_Schoof4Schoof (Hgl,A_mean,C_frot,alpha,n_glen,m_weert,back_force_coef,phi_schoof)
-
-  
-implicit none
-
-real,intent(in) :: Hgl
-real,intent(in) :: A_mean
-real,intent(in) :: C_frot
-real,intent(in) :: alpha
-real,intent(in) :: n_glen
-real,intent(in) :: m_weert
-real,intent(in) :: back_force_coef
-real,intent(out) :: phi_schoof
-
-! afq: in standard GRISLI (19/05/17), tob= -beta ub => m_weert = 1 and C_frot = beta
-
-!phi_schoof = (((A_mean * rog**(n_glen+1) * (1 - alpha)**n_glen) / (4**n_glen *C_frot)) **(1./(1.+1./m_weert)))*back_force_coef**(n_glen/(1+1./m_weert))
-!phi_schoof = phi_schoof * Hgl**((n_glen+3.+1/m_weert)/(1.+1./m_weert))
-phi_schoof = (((A_mean * rog**(n_glen+1.) * (1. - alpha)**n_glen) / (4.**n_glen*C_frot)) **(1./(1.+inv_mweert)))*back_force_coef**(n_glen/(1.+inv_mweert))
-phi_schoof = phi_schoof * Hgl**((n_glen+3.+inv_mweert)/(1.+inv_mweert))
-
-
-
-end subroutine flux_Schoof4Schoof
-
-!--------------------------------------------------------------------
-! Tsai_flux : calcule le coefficient et le flux de Schoof 
-!-------------------------------------------------------------------- 
-
-! Dans un premier temps, pas de terme horizontal (regarder dans toy_recul une tentative de confinement)
-
-subroutine flux_Tsai4Schoof (Hgl,A_mean,f_frot,alpha,n_glen,back_force_coef,phi_Tsai)
-
-implicit none
-
-real,intent(in) :: Hgl
-real,intent(in) :: A_mean
-real,intent(in) :: f_frot
-real,intent(in) :: alpha
-real,intent(in) :: n_glen
-real,intent(in) :: back_force_coef
-real,intent(out) :: phi_Tsai
-
-real, parameter :: Q0 = 0.61
-
-
-phi_Tsai = Q0 * ((8. * A_mean * rog**(n_glen) * (1. - alpha)**(n_glen-1.)) / (4.**n_glen *f_frot))*back_force_coef**(n_glen-1.)
-phi_Tsai = phi_Tsai * Hgl**(n_glen+2.)
-
-end subroutine flux_Tsai4Schoof
+    debug_3D(:,:,66) = phi_prescr_tabx(:,:)
+    debug_3D(:,:,67) = phi_prescr_taby(:,:)
+
+  end subroutine interpol_glflux
+
+  !-------------------------------------------------------------------------
+  ! calcule la position sous maille de la ligne dechouage
+  subroutine calc_xgl4schoof(dy,alpha,H_0,B_0,SL_0,H_1,B_1,SL_1,xpos,Hgl)
+    
+    implicit none
+
+    integer :: i,j
+    real,intent(in)        ::    dy        !< longueur orientee de la maille 
+    real,intent(in)        ::    alpha     !< coefficient de flottaison = rho/rhow
+    real,intent(in)        ::    H_0       !< epaisseur au point pose
+    real,intent(in)        ::    B_0       !< altitude socle au point pose
+    real,intent(in)        ::    SL_0      !< sea level au point pose
+    real,intent(in)        ::    H_1       !< epaisseur au point flottant
+    real,intent(in)        ::    B_1       !< altitude socle au point flottant
+    real,intent(in)        ::    SL_1      !< sea level au point flottant
+    real,intent(out)       ::    xpos      !< position de la ligne (en distance depuis le point pose)
+    real,intent(out)       ::    Hgl       !< epaisseur de glace a la grounding line
+
+
+
+    real                   ::    Cgl       ! variable de travail  
+
+    !~   Cgl   = (-alpha * (H_1-H_0) - (B_1-B_0)) ! -(rho/rhow (H1-H0) + (B1-B0))
+
+    !~   if (abs(Cgl).gt.1.e-5) then   
+    !~      xpos    = (B_0 + alpha * H_0 - SL) / Cgl 
+    !~   else
+    !~      xpos    = 1. - 1./abs(dy)   ! verifier
+    !~   end if
+
+    !~   if (xpos.LT.0..OR.xpos.GT.1.) then
+    !~ 		print*,'calc_xgl4schoof : xpos value error i,j=',i,j
+    !~ 		print*, 'xpos,Cgl', xpos,Cgl
+    !~ 		print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL
+    !~ 		stop
+    !~ 	endif
+
+    Cgl   = (B_1-B_0) + alpha * (H_1-H_0) - (SL_1-SL_0)
+
+    if (abs(Cgl).gt.1.e-5) then   
+       xpos    = (SL_0 - (B_0 + alpha * H_0)) / Cgl 
+    else
+       xpos    = 0.5  ! verifier
+    end if
+
+    !  print*
+    !  print*, 'xpos', xpos
+    !	print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL
+    !	print*,'B_1, H_1', B_1, H_1
+
+
+    if (xpos.LT.0..OR.xpos.GT.1.) then
+       print*,'calc_xgl4schoof : xpos value error i,j=',i,j
+       print*, 'xpos,Cgl', xpos,Cgl
+       print*,'B_0, alpha, H_0, SL_0', B_0, alpha, H_0, SL_0
+       print*,'archim=',B_1+H_1*alpha - SL_1
+       !stop
+       xpos = min(max(0.,xpos),1.)
+    endif
+
+    xpos = xpos * dy
+
+    Hgl= H_0 + (H_1 - H_0) * xpos/dy
+
+  end subroutine calc_xgl4schoof
+
+
+  !--------------------------------------------------------------------
+  ! Schoof_flux : calcule le coefficient et le flux de Schoof 
+  !-------------------------------------------------------------------- 
+
+  ! Dans un premier temps, pas de terme horizontal (regarder dans toy_recul une tentative de confinement)
+
+  subroutine flux_Schoof4Schoof (Hgl,A_mean,C_frot,alpha,n_glen,m_weert,back_force_coef,phi_schoof)
+
+    use param_phy_mod, only:rog
+
+    implicit none
+
+    real,intent(in) :: Hgl
+    real,intent(in) :: A_mean
+    real,intent(in) :: C_frot
+    real,intent(in) :: alpha
+    real,intent(in) :: n_glen
+    real,intent(in) :: m_weert
+    real,intent(in) :: back_force_coef
+    real,intent(out) :: phi_schoof
+
+    ! afq: in standard GRISLI (19/05/17), tob= -beta ub => m_weert = 1 and C_frot = beta
+
+    !phi_schoof = (((A_mean * rog**(n_glen+1) * (1 - alpha)**n_glen) / (4**n_glen *C_frot)) **(1./(1.+1./m_weert)))*back_force_coef**(n_glen/(1+1./m_weert))
+    !phi_schoof = phi_schoof * Hgl**((n_glen+3.+1/m_weert)/(1.+1./m_weert))
+    phi_schoof = (((A_mean * rog**(n_glen+1.) * (1. - alpha)**n_glen) / (4.**n_glen*C_frot)) **(1./(1.+inv_mweert)))*back_force_coef**(n_glen/(1.+inv_mweert))
+    phi_schoof = phi_schoof * Hgl**((n_glen+3.+inv_mweert)/(1.+inv_mweert))
+
+
+
+  end subroutine flux_Schoof4Schoof
+
+  !--------------------------------------------------------------------
+  ! Tsai_flux : calcule le coefficient et le flux de Schoof 
+  !-------------------------------------------------------------------- 
+
+  ! Dans un premier temps, pas de terme horizontal (regarder dans toy_recul une tentative de confinement)
+
+  subroutine flux_Tsai4Schoof (Hgl,A_mean,f_frot,alpha,n_glen,back_force_coef,phi_Tsai)
+
+    use param_phy_mod, only:rog
+
+    implicit none
+
+    real,intent(in) :: Hgl
+    real,intent(in) :: A_mean
+    real,intent(in) :: f_frot
+    real,intent(in) :: alpha
+    real,intent(in) :: n_glen
+    real,intent(in) :: back_force_coef
+    real,intent(out) :: phi_Tsai
+
+    real, parameter :: Q0 = 0.61
+
+
+    phi_Tsai = Q0 * ((8. * A_mean * rog**(n_glen) * (1. - alpha)**(n_glen-1.)) / (4.**n_glen *f_frot))*back_force_coef**(n_glen-1.)
+    phi_Tsai = phi_Tsai * Hgl**(n_glen+2.)
+
+  end subroutine flux_Tsai4Schoof
 
 end module furst_schoof_mod
