Index: branches/GRISLIv3/SOURCES/spinup_mod.f90
===================================================================
--- branches/GRISLIv3/SOURCES/spinup_mod.f90	(revision 364)
+++ branches/GRISLIv3/SOURCES/spinup_mod.f90	(revision 429)
@@ -19,10 +19,6 @@
 module spinup_vitbil
 
-  use module3D_phy
-  use param_phy_mod
-  use deformation_mod_2lois
-  use interface_input
-  use io_netcdf_grisli
-  
+  use geography, only: nx,ny,nz
+
   implicit none
 
@@ -36,10 +32,10 @@
   real, dimension(nx,ny) :: init_coefmy      !< rescalling coefficient before limitation
 
-! compatibilites avec remplimat
+  ! compatibilites avec remplimat
   logical :: corr_def = .False.               !< for deformation correction (compatibility)
 
   real, dimension(nx,ny) :: Vsl_x             !<
   real, dimension(nx,ny) :: Vsl_y             !<
-  
+
   integer :: type_vitbil ! defini le type de fichier a lire : lect_vitbil ou lect_vitbil_Lebrocq
 
@@ -51,5 +47,9 @@
   !<
   subroutine init_spinup
-  namelist/spinup/ispinup,type_vitbil
+
+    use module3D_phy, only: ispinup,num_param,num_rep_42
+    use runparam, only : itracebug
+    
+    namelist/spinup/ispinup,type_vitbil
 
     ! put here which type of spinup
@@ -63,47 +63,47 @@
     ! ispinup = 3                          ! fait le calcul des temperatures 
     !                                      ! avec les vitesses de bilan
-  
-
-
-
-! lecture des parametres
-
-rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
-read(num_param,spinup)
-
-write(num_rep_42,*)'!___________________________________________________________'
-write(num_rep_42,*)'!        spinup          module spinup_vitbil               '
-write(num_rep_42,spinup)                   ! pour ecrire les parametres lus
-write(num_rep_42,*)
-write(num_rep_42,*)
-write(num_rep_42,*)'! ispinup = 0     run standard voir no_spinup'
-write(num_rep_42,*)'! ispinup = 1     equilibre temperature avec vitesses grisli voir no_spinup'
-write(num_rep_42,*)
-write(num_rep_42,*)'! ispinup >1      use spinup_vitbil'
-write(num_rep_42,*)'! ispinup = 2     conservation de la masse avec vitesses bilan '
-write(num_rep_42,*)'! ispinup = 3     equilibre temperature avec vitesses bilan'
-write(num_rep_42,*)'! type_vitbil type de lecture des vitesses de bilan 1: lect_vitbil, 2 : lect_vitbil_Lebrocq'
-write(num_rep_42,*)
-
- 
-if (ispinup.eq.0) then
-   write(6,*)' ispinup = 0 et ispinup = 1 doivent etre appeles par no_spinup '
-   write(6,*) 'et il faut rajouter un dragging'
-   stop
-endif
-
-if (ispinup.eq.3) then
-  select case (type_vitbil)
-  case (1)
-!cdc Methode Cat lecture fichier selon x et y sur grille stag
-    call lect_vitbil
-  case(2)
-!cdc Methode Tof lecture fichier Lebrocq vitesse et direction
-    call lect_vitbil_Lebrocq
-  case default
-    print*,'type_vitbil valeur invalide dans spinup_vitbil'
-    stop
-  end select
-endif
+
+
+
+
+    ! lecture des parametres
+
+    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
+    read(num_param,spinup)
+
+    write(num_rep_42,*)'!___________________________________________________________'
+    write(num_rep_42,*)'!        spinup          module spinup_vitbil               '
+    write(num_rep_42,spinup)                   ! pour ecrire les parametres lus
+    write(num_rep_42,*)
+    write(num_rep_42,*)
+    write(num_rep_42,*)'! ispinup = 0     run standard voir no_spinup'
+    write(num_rep_42,*)'! ispinup = 1     equilibre temperature avec vitesses grisli voir no_spinup'
+    write(num_rep_42,*)
+    write(num_rep_42,*)'! ispinup >1      use spinup_vitbil'
+    write(num_rep_42,*)'! ispinup = 2     conservation de la masse avec vitesses bilan '
+    write(num_rep_42,*)'! ispinup = 3     equilibre temperature avec vitesses bilan'
+    write(num_rep_42,*)'! type_vitbil type de lecture des vitesses de bilan 1: lect_vitbil, 2 : lect_vitbil_Lebrocq'
+    write(num_rep_42,*)
+
+
+    if (ispinup.eq.0) then
+       write(6,*)' ispinup = 0 et ispinup = 1 doivent etre appeles par no_spinup '
+       write(6,*) 'et il faut rajouter un dragging'
+       stop
+    endif
+
+    if (ispinup.eq.3) then
+       select case (type_vitbil)
+       case (1)
+          !cdc Methode Cat lecture fichier selon x et y sur grille stag
+          call lect_vitbil
+       case(2)
+          !cdc Methode Tof lecture fichier Lebrocq vitesse et direction
+          call lect_vitbil_Lebrocq
+       case default
+          print*,'type_vitbil valeur invalide dans spinup_vitbil'
+          stop
+       end select
+    endif
 
     if (itracebug.eq.1)  call tracebug(' fin routine init_spinup de spinup_vitbil')
@@ -111,9 +111,14 @@
 
   subroutine lect_vitbil
-
+    
+    use module3D_phy, only: debug_3D,num_rep_42,num_param
+    use io_netcdf_grisli, only: Read_Ncdf_var
+    use geography, only: dirnameinp
+    use runparam, only : itracebug
+    
     character(len=100) :: balance_Ux_file        !  balance velocity file Ux
     character(len=100) :: balance_Uy_file        !  balance velocity file Uy
-	  real*8, dimension(:,:),   pointer    :: tab               !< tableau 2d real pointer
-	  
+    real*8, dimension(:,:),   pointer    :: tab               !< tableau 2d real pointer
+
     namelist/vitbil_upwind/balance_Ux_file, balance_Uy_file  
 
@@ -135,8 +140,8 @@
 
     ! read balance velocities on staggered nodes
-    
+
     balance_Ux_file = trim(dirnameinp)//trim(balance_Ux_file)
     balance_Uy_file = trim(dirnameinp)//trim(balance_Uy_file)
-    
+
     call Read_Ncdf_var('Vcol_x',balance_Ux_file,tab)
     Vcol_x(:,:)=tab(:,:)
@@ -144,6 +149,6 @@
     Vcol_y(:,:)=tab(:,:)
 
-!    call lect_input(2,'Vcol',1,Vcol_x,balance_vel_file,trim(dirnameinp)//trim(runname)//'.nc')  !read Vcol_x
-!    call lect_input(2,'Vcol',2,Vcol_y,balance_vel_file,trim(dirnameinp)//trim(runname)//'.nc')  !read Vcol_y
+    !    call lect_input(2,'Vcol',1,Vcol_x,balance_vel_file,trim(dirnameinp)//trim(runname)//'.nc')  !read Vcol_x
+    !    call lect_input(2,'Vcol',2,Vcol_y,balance_vel_file,trim(dirnameinp)//trim(runname)//'.nc')  !read Vcol_y
     !call lect_datfile(nx,ny,Vcol_x,1,balance_vel_file)   ! read Vcol_x   
     !call lect_datfile(nx,ny,Vcol_y,2,balance_vel_file)   ! read Vcol_y 
@@ -153,8 +158,12 @@
 
   end subroutine lect_vitbil
-  
+
   subroutine lect_vitbil_Lebrocq
 
-
+    use module3D_phy, only: num_rep_42,mstream_mx,itracebug,num_param,pi
+    use io_netcdf_grisli, only: Read_Ncdf_var
+    use geography, only: dirnameinp
+    use runparam, only : itracebug
+    
     character(len=100) :: vit_balance_file        !  balance velocity file
     character(len=100) :: flowdir_balance_file    !  balance flow direction file
@@ -164,5 +173,6 @@
     real, dimension(nx,ny) :: Vx_Lebrocq, Vy_Lebrocq ! vitesses selon x et y
     integer, dimension(nx,ny) :: v_good  !points avec vitesse Lebrocq OK
-    
+    integer :: i,j
+
     namelist/vitbil_Lebrocq/vit_balance_file, flowdir_balance_file  
 
@@ -184,8 +194,8 @@
 
     ! read balance velocities on staggered nodes
-    
+
     vit_balance_file = trim(dirnameinp)//trim(vit_balance_file)
     flowdir_balance_file = trim(dirnameinp)//trim(flowdir_balance_file)
-    
+
     call Read_Ncdf_var('z',vit_balance_file,tab)
     V_Lebrocq(:,:)=tab(:,:)
@@ -194,45 +204,45 @@
 
     where (V_Lebrocq(:,:).ge.0.)
-      v_good(:,:)=1
+       v_good(:,:)=1
     elsewhere
-      v_good(:,:)=0
+       v_good(:,:)=0
     endwhere
-    
+
     do j=2,ny-1
-      do i=2,nx-1
-        if (V_Lebrocq(i,j).lt.0.) then
-          
-          V_Lebrocq(i,j)= (V_Lebrocq(i-1,j)*v_good(i-1,j) + V_Lebrocq(i+1,j)*v_good(i+1,j) + &
-                         V_Lebrocq(i,j-1)*v_good(i,j-1) + V_Lebrocq(i,j+1)*v_good(i,j+1)) / &
-                         (v_good(i-1,j)+v_good(i+1,j)+v_good(i,j-1)+v_good(i,j-1))
-        endif
-      enddo
+       do i=2,nx-1
+          if (V_Lebrocq(i,j).lt.0.) then
+
+             V_Lebrocq(i,j)= (V_Lebrocq(i-1,j)*v_good(i-1,j) + V_Lebrocq(i+1,j)*v_good(i+1,j) + &
+                  V_Lebrocq(i,j-1)*v_good(i,j-1) + V_Lebrocq(i,j+1)*v_good(i,j+1)) / &
+                  (v_good(i-1,j)+v_good(i+1,j)+v_good(i,j-1)+v_good(i,j-1))
+          endif
+       enddo
     enddo
 
 
     where (V_Lebrocq(:,:).ge.100.)
-      V_Lebrocq(:,:)=100.
+       V_Lebrocq(:,:)=100.
     endwhere
 
-		! calcul des vitesses selon x et selon y :
-		where (V_Lebrocq(:,:).ge.0.)
-			Vx_Lebrocq(:,:) = V_Lebrocq(:,:)*sin(flowdir_Lebrocq(:,:))
-			Vy_Lebrocq(:,:) = V_Lebrocq(:,:)*cos(flowdir_Lebrocq(:,:))
-		elsewhere
-			Vx_Lebrocq(:,:) = 0.
-			Vy_Lebrocq(:,:) = 0.
-		endwhere
-
-		
-		! calcul des vitesses sur les points staggered
-		
-		do j=2,ny
-			do i=2,nx
-				Vcol_x(i,j) = (Vx_Lebrocq(i-1,j) + Vx_Lebrocq(i,j))/2.
-				Vcol_y(i,j) = (Vy_Lebrocq(i,j-1) + Vy_Lebrocq(i,j))/2.
-			enddo
-		enddo
-
-  end subroutine lect_vitbil_Lebrocq 
+    ! calcul des vitesses selon x et selon y :
+    where (V_Lebrocq(:,:).ge.0.)
+       Vx_Lebrocq(:,:) = V_Lebrocq(:,:)*sin(flowdir_Lebrocq(:,:))
+       Vy_Lebrocq(:,:) = V_Lebrocq(:,:)*cos(flowdir_Lebrocq(:,:))
+    elsewhere
+       Vx_Lebrocq(:,:) = 0.
+       Vy_Lebrocq(:,:) = 0.
+    endwhere
+
+
+    ! calcul des vitesses sur les points staggered
+
+    do j=2,ny
+       do i=2,nx
+          Vcol_x(i,j) = (Vx_Lebrocq(i-1,j) + Vx_Lebrocq(i,j))/2.
+          Vcol_y(i,j) = (Vy_Lebrocq(i,j-1) + Vy_Lebrocq(i,j))/2.
+       enddo
+    enddo
+
+  end subroutine lect_vitbil_Lebrocq
   !____________________________________________________________________________________
 
@@ -242,4 +252,6 @@
   subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
 
+    use module3D_phy, only: mstream_mx,mstream_my,drag_mx,drag_my,gzmx,gzmy,flgzmx,flgzmy,flotmx,flotmy
+    
     implicit none
 
@@ -265,5 +277,7 @@
 
   subroutine force_balance_vel
-
+    
+    use module3D_phy, only: uxbar, uybar,debug_3D
+    use runparam, only : itracebug
     if (itracebug.eq.1)  call tracebug(' Subroutine force_balance_vel')
 
@@ -281,5 +295,12 @@
   !<
   subroutine calc_coef_vitbil
-
+    
+    use module3D_phy, only: ddbx,ddby,gzmx,gzmy,flgzmx,flgzmy,fleuvemx,fleuvemy,flot,flotmx,flotmy,&
+         coefmxbmelt,coefmybmelt,sdx,sdy,uxbar,uybar,uxdef,uydef,ubx,uby,debug_3D,iglen
+    use deformation_mod_2lois, only: n1poly,n2poly,sa_mx,sa_my,s2a_mx,s2a_my
+    use runparam, only : itracebug
+    
+    integer :: i,j,k
+    
     ddbx(:,:)=0.
     ddby(:,:)=0.
@@ -290,25 +311,25 @@
     fleuvemx(:,:)=.false.
     fleuvemy(:,:)=.false.
-    
+
     do j=2,ny
-      do i=2,nx
+       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.
+             flotmx(i,j)=.true.
+             flgzmx(i,j)=.true.
+             gzmx(i,j)=.false.
           end if
           if (flot(i,j).and.(flot(i,j-1))) then
-            flotmy(i,j)=.true.
-            flgzmy(i,j)=.true.
-            gzmy(i,j)=.false.
+             flotmy(i,j)=.true.
+             flgzmy(i,j)=.true.
+             gzmy(i,j)=.false.
           end if
-      end do
-    end do   
-    
+       end do
+    end do
+
     where (coefmxbmelt(:,:).le.0.)
-      gzmx(:,:)=.false.
+       gzmx(:,:)=.false.
     endwhere
     where (coefmybmelt(:,:).le.0.)
-      gzmy(:,:)=.false.
+       gzmy(:,:)=.false.
     endwhere
     flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
@@ -316,5 +337,5 @@
     fleuvemx(:,:)=gzmx(:,:)
     fleuvemy(:,:)=gzmy(:,:)
-    
+
     if (itracebug.eq.1)  call tracebug(' Subroutine  calc_coef_vitbil')
     call slope_surf
@@ -341,7 +362,7 @@
              coef_defmx(i,j) = 1.
              Uxbar(i,j)      = Vcol_x(i,j)
-! dmr below is a cast from a real*4 to logical*4
-! dmr cannot be implicit in gfortran
-! dmr              flgzmx(i,j)     = Vcol_x(i,j)
+             ! dmr below is a cast from a real*4 to logical*4
+             ! dmr cannot be implicit in gfortran
+             ! dmr              flgzmx(i,j)     = Vcol_x(i,j)
              flgzmx(i,j)     = (nint(Vcol_x(i,j)).ne.0)
              uxdef(i,j)      = 0.
@@ -365,5 +386,5 @@
                    coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j)
                    ddbx(i,j) = 0
-                    Ubx(i,j)  = 0.
+                   Ubx(i,j)  = 0.
                 else
                    coef_defmx(i,j)=1.                                  ! no rescaling of deformation
@@ -385,7 +406,7 @@
              coef_defmy(i,j) = 1.
              Uybar(i,j)      = Vcol_y(i,j)
-! dmr below is a cast from a real*4 to logical*4
-! dmr cannot be implicit in gfortran
-! dmr             flgzmy(i,j)     = Vcol_y(i,j)
+             ! dmr below is a cast from a real*4 to logical*4
+             ! dmr cannot be implicit in gfortran
+             ! dmr             flgzmy(i,j)     = Vcol_y(i,j)
              flgzmy(i,j)     = (nint(Vcol_y(i,j)).ne.0)
 
@@ -393,5 +414,5 @@
              Uby(i,j)        = Vcol_y(i,j)
 
- if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test0')
+             if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test0')
           else
              coldy:      if ((coefmybmelt(i,j).le.0.).and.(.not.gzmy(i,j))) then   !base froide 
@@ -402,7 +423,7 @@
                 if (abs(Uy_deformation(i,j)).gt.0.01) then             ! calcule par calc_ubar_def
                    coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j)
- if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test1')
+                   if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test1')
                 else
- if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test2')
+                   if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test2')
                    coef_defmy(i,j)=1.
                 end if
@@ -412,10 +433,10 @@
                    ddby(i,j) = 0.
                    Uby(i,j)  = 0.
- if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test3')
+                   if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test3')
                 else
                    coef_defmy(i,j)=1.                                  ! pas de rescaling de la deformation
                    ddby(i,j)=-(Vcol_y(i,j)-Uy_deformation(i,j))/sdy(i,j)        ! pb si directions opposees
                    Uby(i,j) =  Vcol_y(i,j)-Uy_deformation(i,j)
-if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.144))  call tracebug(' test4')
+                   if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.144))  call tracebug(' test4')
 
                 end if
@@ -426,5 +447,5 @@
        end do
     end do
- if (itracebug.eq.1)  call tracebug(' apres test floty')
+    if (itracebug.eq.1)  call tracebug(' apres test floty')
 
     ! remarque coef peut être <1  a cause de la pente locale, mais il joue toujours son role 
@@ -453,5 +474,12 @@
 
   subroutine limit_coef_vitbil
-
+    
+    use module3D_phy, only: sdx,sdy,flotmx,flotmy,coefmxbmelt,coefmybmelt,gzmx,gzmy,ddbx,ddby,debug_3D,&
+         iglen
+    use deformation_mod_2lois, only: n1poly,n2poly,sa_mx,sa_my,s2a_mx,s2a_my
+    use runparam, only : itracebug
+    
+    integer :: i,j,k
+    
     call slope_surf
 
@@ -583,4 +611,10 @@
   subroutine calc_ubar_def   ! calculate velocity due to deformation before calibration
 
+    use module3D_phy, only: iglen,flotmx,flotmy,hmx,hmy,sdx,sdy,slope2mx,slope2my
+    use deformation_mod_2lois, only: n1poly,n2poly,glen,s2a_mx,s2a_my,ddx,ddy
+    use geography, only: dx,dy
+    use param_phy_mod, only: rog
+    use runparam, only: itracebug
+    
     implicit none              ! extrait de diffusiv pour calculer la partie deformation
 
@@ -588,4 +622,5 @@
     real :: inv_4dx ! inverse de dx pour eviter les divisions =1/(4*dx)
     real :: inv_4dy ! inverse de dy pour eviter les divisions =1/(4*dy)
+    integer :: i,j
 
     inv_4dx=1./(4*dx)
@@ -646,9 +681,14 @@
 
   subroutine ajust_ghf
-
+    
+    use module3D_phy, only: debug_3D,flot,ibase,ghf,ghf0,secyear
+    
     implicit none
+    
     real,dimension(nx,ny) :: coefdef_maj     !< coefficient deformation
     real :: increment_ghf
     real :: ghf_min
+    integer :: i,j
+    
     increment_ghf=0.00001                    !< exprime en  W/m2   
     ghf_min=0.025                                           
