Index: trunk/SOURCES/Fichiers-parametres/Makefile.tof-lsce3130.inc
===================================================================
--- trunk/SOURCES/Fichiers-parametres/Makefile.tof-lsce3130.inc	(revision 13)
+++ trunk/SOURCES/Fichiers-parametres/Makefile.tof-lsce3130.inc	(revision 14)
@@ -36,5 +36,5 @@
 IFORT= ifort
 
-ARITHM    = -O2 -fp-model precise  -heap-arrays -traceback -diag-disable warn # -traceback   -CB  -g  # options pour une meilleure arithmetique
+ARITHM    = -O2 -fp-model precise  -heap-arrays -traceback -diag-disable warn -mcmodel=medium # -traceback   -CB  -g  # options pour une meilleure arithmetique
                                                                       # (normalement reproductible)
 #ARITHM    =  -fp-model precise -warn all -CU -CA # options pour une meilleure arithmetique (normalement reproductible)
Index: trunk/SOURCES/Fichiers-parametres/hemin15_LISTE-VAR-NETCDF.dat
===================================================================
--- trunk/SOURCES/Fichiers-parametres/hemin15_LISTE-VAR-NETCDF.dat	(revision 14)
+++ trunk/SOURCES/Fichiers-parametres/hemin15_LISTE-VAR-NETCDF.dat	(revision 14)
@@ -0,0 +1,702 @@
+LISTE-VAR-NETCDF.dat   nom du fichier
+
+Liste des variables a sortir sur les fichier netcdf.
+Toutes les variables pour lesquelle une sortie est prevue sont dans cette liste.
+Les variables dans ce fichier sont regroupees par classe
+
+Une variable est sortie si :
+le module correspondant est active (exemple eau basale)
+et isortie=1
+
+La periodicite avec  laquelle ces variables sont sorties 
+ainsi que les temps predefinis eventuels sont donnes dans
+TEMPS-NETCDF.dat 
+
+classnumber : classe de variable
+
+typenod : 'o' noeuds majeurs, '>' noeuds mineurs x, '^' mineurs y
+	      'x' noeuds milieux
+
+snapmax : nombre d'enregistrements max dans un fichier netcdf, 
+si plus, il reouvre un nouveau fichier.
+
+2 lignes par classe sous la forme :
+------------------------------------------
+classnumber snapmax  ntypenode
+typenod1 typenod2 typenod3 ...
+
+
+ex :
+============================================ classe 1 
+1	100    4                   ! class 1  snapmax=10   4 type de noeud 
+'o' 	'>'    	  '^'  	    'x'	   ! les types de noeuds
+------------------------------------------
+2 lignes par variable  sous la forme :
+varname1 !attention varname1 5 caract maxi et pas de TAB apres	
+numero    		
+isortie  dtncdf interv
+!!! Possibilite de mettre des commentaires
+ligne vide
+-------------------------
+varname2 
+numero
+isortie  dtncdf interv
+-------------------------
+
+Entre les classes, le separateur est la ligne   ==========
+Entre les variables, le separateur est la ligne ----------
+Il ne teste que les 10 premiers caracteres de ========== et
+de ----------
+
+
+classnumber :   identificateur de classe de variables
+	    	0 = coordonnees (x,y,lat,long...)
+		1 = Var 2D a sortir regulierement (variables etats, geometrie ...)
+		2 = Var 2D  a sortir+ Debug
+		3 = Var 3D 
+		4 = Var de debug
+		5 = sortie juste une fois au debut et une fois a la fin
+
+snapmax :       nombre de snapshot max aÂ  sortir   
+
+varname : 5 caracteres maxi (un blanc sera ajoutÃÂ© au dÃÂ©but par le code)
+	  a aligner a gauche, pas de TAB apres
+
+isortie : 1-> oui, 0->non (pas de sortie du tout)
+
+dtncdf  : donne la periodicite avec laquelle ces variables sont sorties
+          indice pointe dans dtsortie_ncdf definie dans TEMPS-NETCDF.dat
+          si 0 ne sort que sur les pas de temps intermediaires 
+
+interv : pour cette variable  
+  		     interv=-1  premier+dernier
+                     interv= 0  ne sort que le premier pas de temps
+		     interv= 1  premier + dernier + predefinis
+
+ex. 
+1      1     0    sort avec periodicite 1 + premier pas de temps
+1      0     -1   sort seulement premier et dernier
+0      1     1    ne sort pas du tout
+1      0     0    ne sort que le premier
+
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ fin du header commentaire  							       	   
+============================================ classe 1 
+1	220     4           ! class 1 nombre de snapmax 100 4 type de noeud 
+'o' 	'>'    	  '^'  	    'x'					       
+-----------------------------------------------------------
+S	                   ! nom de la variable dans GRISLI (celui du tableau)
+1                          ! numero de la variable      
+1	1	1          ! isortie  ,  dtncdf (periodicite) , interv
+!!!! commentaire
+!!! comm
+Ligne vide 
+-----------------------------------------------------------
+H	
+2                
+1	1	1
+-----------------------------------------------------------
+Bsoc	
+3
+1	1	1
+-----------------------------------------------------------
+Mk	
+4
+1	1	1			
+-----------------------------------------------------------
+Hdot	
+5
+1	1	1			
+-----------------------------------------------------------
+B		
+7		
+1	1	1		
+-----------------------------------------------------------
+Bm		
+10
+1	1	1				
+-----------------------------------------------------------
+Acc		
+11	
+1	4	1				
+-----------------------------------------------------------
+Abl		
+12
+1	1	1		
+-----------------------------------------------------------
+Calv		
+13
+0	1	1				
+-----------------------------------------------------------
+dHdt		
+14
+0	1	1			
+-----------------------------------------------------------
+Tann		
+18
+1	1	1			
+-----------------------------------------------------------
+Tjja		
+19
+1	1	1			
+-----------------------------------------------------------
+Tb		
+20
+1	1	1			
+-----------------------------------------------------------
+phid		
+24	
+1	1       1				
+-----------------------------------------------------------
+Bmelt		
+25
+1	1	1				
+-----------------------------------------------------------
+Uxbar_maj		
+30
+1	1	1				
+-----------------------------------------------------------
+Uybar_maj		
+31
+1	1	1			
+-----------------------------------------------------------
+Umag		
+36
+1	1	1				
+-----------------------------------------------------------
+Uslmag		
+37
+1	1	1				
+-----------------------------------------------------------
+bm_ref_t
+197
+1	1	1
+-----------------------------------------------------------
+grad_bm
+198
+1	1	1
+-----------------------------------------------------------
+bm_time
+199
+1	1	1
+-----------------------------------------------------------
+============================================ classe 2 
+2       220     3  ! class 2 nombre de snapmax 100  3 type de noeud a sortir
+'>'    	'^'   	'o'   	
+-----------------------------------------------------------
+dS	         ! ecart par rapport a la surface observee	
+6
+1	1	1						
+-----------------------------------------------------------
+Uxbar	        ! sur noeud mineur	
+32	
+1	1	1				
+-----------------------------------------------------------
+Uybar	        ! sur noeud mineur	
+33
+1	1	1			
+-----------------------------------------------------------
+Uxsli	        ! sur noeud mineur	
+34
+1	1	1			
+-----------------------------------------------------------
+Uysli		 ! sur noeud mineur	
+35
+1	1	1				
+-----------------------------------------------------------
+Uxdef		 ! sur noeud mineur	
+38
+1	1	1	
+-----------------------------------------------------------
+Uydef		 ! sur noeud mineur	
+39
+1	1	1				
+-----------------------------------------------------------
+frtmx		
+40
+1	1	1				
+-----------------------------------------------------------
+frtmy	         ! sur noeud mineur	
+41
+1	1	1				
+-----------------------------------------------------------
+tobmx		 ! sur noeud mineur	
+42
+1	1	1					
+-----------------------------------------------------------
+tobmy		
+43
+1	1	1					
+-----------------------------------------------------------
+toshe		
+44
+0	1	1					
+-----------------------------------------------------------
+epsxx	        !'o'	
+45
+1	3	-1					
+-----------------------------------------------------------
+epsyy	        !'o'	
+46
+1	3	-1					
+-----------------------------------------------------------
+epsxy	        !'x'	
+47
+1	3	1					
+-----------------------------------------------------------
+eps	        'o'	
+48
+1	3	1					
+-----------------------------------------------------------
+pvi	         ! 'o'	
+50
+1	1	1					
+-----------------------------------------------------------
+pvm		!  'x'
+51
+1	1	1					
+-----------------------------------------------------------
+betamx		
+52
+1	1	1					
+-----------------------------------------------------------
+betamy		
+53
+1	1	1					
+-----------------------------------------------------------
+beta_c		
+54
+1	1	1					
+-----------------------------------------------------------
+Hwat		
+60
+1	1	1					
+-----------------------------------------------------------
+Hdw		
+61
+0	1	1					
+-----------------------------------------------------------
+Pgx		
+62
+0	1	1				
+-----------------------------------------------------------
+Pgy		
+63
+0	1	1					
+-----------------------------------------------------------
+Kond		
+64
+0	1	1					
+-----------------------------------------------------------
+phx		
+65
+0	1	1					
+-----------------------------------------------------------
+phy		
+66
+0	1	1					
+-----------------------------------------------------------
+Nefmx		
+68
+1	1	1					
+-----------------------------------------------------------
+Nefmy		
+69
+1	2	1					
+-----------------------------------------------------------
+posx		
+70
+1	1	1					
+-----------------------------------------------------------
+posy		
+71
+1	1	1				
+-----------------------------------------------------------
+fronx		
+72
+0	1	1					
+-----------------------------------------------------------
+frony		
+73
+0	1	1	
+-----------------------------------------------------------
+H_uw_x
+192
+1	1	1
+-----------------------------------------------------------
+H_uw_y
+193
+1	1	1
+-----------------------------------------------------------		
+============================================ classe 3 
+3       100       4   !class 3 nombre de snapmax 100 4 type de noeud
+'o' 	'>'    	  '^'  	    'x'
+-----------------------------------------------------------
+CP		
+74
+0	0	-1	
+-----------------------------------------------------------	
+CT		
+75
+0	0	-1	
+-----------------------------------------------------------
+SUX		
+76
+0	0	-1	
+-----------------------------------------------------------
+SUY		
+77
+0	0	-1	
+-----------------------------------------------------------
+TPMP		
+78
+1	0	-1
+-----------------------------------------------------------
+UX		
+79
+1	3	1	
+-----------------------------------------------------------
+UY		
+80
+1	3	1	
+-----------------------------------------------------------
+UZR		
+81
+1	3	1	
+-----------------------------------------------------------
+Chaldef_maj		
+82
+1	0	-1	
+-----------------------------------------------------------
+T		
+83
+1	3	1 
+-----------------------------------------------------------
+xdep_out		
+84
+0	0	0 	
+-----------------------------------------------------------
+ydep_out		
+85
+0	0	0 	
+-----------------------------------------------------------
+tdep_out		
+86
+0	0	0 	
+-----------------------------------------------------------
+!! Variables 3 D a ajouter
+!! SA, S2A,SA_mx,S2A_mx,SA_my,S2A_my
+
+============================================ classe 4 
+4       20     4   !class 4 nombre de snapmax 100 4 type de noeud
+'o' 	'>'    	  '^'  	    'x'
+-----------------------------------------------------------
+btot		
+15
+0	2	-1			
+-----------------------------------------------------------
+bshelftot	
+16
+0	1	-1			
+-----------------------------------------------------------
+fleuve		
+101	
+0	1	-1				
+-----------------------------------------------------------
+cote		
+102
+0	1	-1					
+-----------------------------------------------------------
+fleuvy		
+103
+1	1	-1				
+-----------------------------------------------------------
+fleuvy		
+104
+1	1	-1				
+-----------------------------------------------------------
+pot_w		
+105	
+0	1	-1					
+-----------------------------------------------------------
+pot_f		
+106	
+0	1	-1					
+-----------------------------------------------------------
+pot_tot		
+107
+0	1	-1				
+-----------------------------------------------------------
+pot_hw		
+108
+0	1	-1					
+-----------------------------------------------------------
+uxdef		
+109
+1	1	1			
+-----------------------------------------------------------
+uydef		
+110
+1	1	1				
+-----------------------------------------------------------
+uxflgz		
+111
+1	1	-1			
+-----------------------------------------------------------
+uyflgz		
+112
+1	1	1				
+-----------------------------------------------------------
+uxsl		
+113
+1	1	1					
+-----------------------------------------------------------
+uysl		
+114
+1	1	1					
+-----------------------------------------------------------
+B_sm		
+115
+0	1	1					
+-----------------------------------------------------------
+Bslop		
+116
+0	1	-1					
+-----------------------------------------------------------
+Bcrx		
+117
+0	1	-1			
+-----------------------------------------------------------
+Bcry		
+118
+0	1	-1					
+-----------------------------------------------------------
+Bcrxy		
+119
+0	1	-1				
+-----------------------------------------------------------
+Bteta		
+120
+0	1	-1				
+-----------------------------------------------------------
+sdx		
+121	
+1	1	-1					
+-----------------------------------------------------------
+sdy		
+122
+1	1	-1					
+-----------------------------------------------------------
+motmx		
+123
+0	1	1					
+-----------------------------------------------------------
+motmy		
+124
+0	1	1				
+-----------------------------------------------------------
+str		
+125
+1	1	1					
+-----------------------------------------------------------
+TUIJ		
+126
+0	1	-1					
+-----------------------------------------------------------
+ppvi		
+127
+0	1	-1				
+-----------------------------------------------------------
+bms		
+128
+0	1	-1				
+-----------------------------------------------------------
+dTann		
+129
+0	1	-1					
+-----------------------------------------------------------
+rapac		
+130
+0	1	-1		
+-----------------------------------------------------------
+bmtot		
+131
+1	1	-1					
+-----------------------------------------------------------
+dist_talus		
+132
+0	1	-1			
+-----------------------------------------------------------
+typeshelf	
+133	
+0	1	-1				
+-----------------------------------------------------------
+bmshelf		
+134
+0	1	-1			
+-----------------------------------------------------------
+opposx		
+135
+1	1	1		
+-----------------------------------------------------------
+opposy		
+136
+1	1	1					
+-----------------------------------------------------------
+uxnew		
+137
+1	1	-1				
+-----------------------------------------------------------
+uynew		
+138
+0	1	-1					
+-----------------------------------------------------------
+okumat		
+139
+0	1	-1		
+-----------------------------------------------------------
+okvmat		
+140
+0	1	-1					
+-----------------------------------------------------------
+ghost_x	
+141
+0	1	-1					
+-----------------------------------------------------------
+ghost_y			
+142
+0	1	-1					
+-----------------------------------------------------------
+tbout		
+143
+0	1	-1				
+-----------------------------------------------------------
+bm-bmelt		
+144
+1	1	-1		
+-----------------------------------------------------------		
+difx			
+145
+0	1	-1				
+-----------------------------------------------------------
+dify		
+146	
+0	1	-1					
+-----------------------------------------------------------
+advmx		
+147
+0	1	-1					
+-----------------------------------------------------------
+advmy		
+148
+0	1	-1					
+-----------------------------------------------------------
+deltaH		
+149
+0	1	-1					
+-----------------------------------------------------------
+arelax		
+150
+0	1	-1					
+-----------------------------------------------------------
+brelax		
+151
+0	1	-1				
+-----------------------------------------------------------
+crelax		
+152
+0	1	-1				
+-----------------------------------------------------------
+drelax		
+153
+0	1	-1				
+-----------------------------------------------------------
+erelax		
+154
+0	1	-1			
+-----------------------------------------------------------		
+frelax  	
+155
+0	1	-1		
+-----------------------------------------------------------		
+nb_ile		
+156								
+0	1	-1			
+-----------------------------------------------------------		
+berg    	
+157 											
+0	1	-1			
+-----------------------------------------------------------
+i_hp	   	
+185 											
+1	1	1			
+-----------------------------------------------------------
+hp	   	
+186 											
+1	1	1			
+-----------------------------------------------------------
+coefdef_mx	
+161
+0	1	1
+-----------------------------------------------------------
+coefdef_my
+162
+0	1	1
+-----------------------------------------------------------		
+Vslmag	
+168
+1	1	1
+-----------------------------------------------------------
+Vsl_x		
+169
+1	1	1					
+-----------------------------------------------------------
+Vsl_y	
+170
+1	1	1
+-----------------------------------------------------------
+stress_x		
+182
+1	1	1					
+-----------------------------------------------------------
+stress_y
+183
+1	1	1
+-----------------------------------------------------------
+S2a_law1
+190
+1	1	1
+-----------------------------------------------------------
+S2a_law2
+191
+1	1	1
+-----------------------------------------------------------
+============================================ classe 5
+5       100     4       ! class 5 nombre de snapmax 100 4 type de noeud a sortir
+'o' 	'>'    	  '^'  	    'x'
+-----------------------------------------------------------
+B_cry		
+8
+0	1	-1				
+-----------------------------------------------------------
+phi		
+23	
+1	0	-1				
+-----------------------------------------------------------
+Mk_init		
+9	
+1	0	-1		
+-----------------------------------------------------------
+Vcolmag	
+158
+1	0	-1
+-----------------------------------------------------------
+Vcol_x		
+159
+1	0	-1					
+-----------------------------------------------------------
+Vcol_y	
+160
+1	0	-1
+-----------------------------------------------------------
+					
+==========================================================
Index: trunk/SOURCES/Fichiers-parametres/hemin15_TEMPS-NETCDF.dat
===================================================================
--- trunk/SOURCES/Fichiers-parametres/hemin15_TEMPS-NETCDF.dat	(revision 14)
+++ trunk/SOURCES/Fichiers-parametres/hemin15_TEMPS-NETCDF.dat	(revision 14)
@@ -0,0 +1,43 @@
+TEMPS-HZ.dat
+
+Permet de définir les temps auxquels ont lieu les sorties hz :
+Deux méthodes complémentaires: 
+
+1) sortie régulière. dtsortie_hz(ndtsortie) est un tableau qui contient 
+les pas de temps de sortie. Pour chaque variable on choisi celui qu'on
+veut par indice dans dt_out_hz
+
+2) predef_tsort(npredeft)  est un tableaux de temps prédéfinis pour les sorties
+
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ fin des commentaires~~~~~~~~~~~~~
+2			nbr de classes a sortir	dans ce run	    	
+1  2	class   (voir les classes et variables dans LISTE_VAR_NCDF.dat)
+-------------------------------------------------------------------------------------------
+4	ndtsortie      lecture de dtsortie_hz
+10000.
+50000.
+100000.
+1.e10
+-------------------------------------------------------------------------------------------
+15	npredeft        lecture de predef_tsort
+10
+100
+200
+500
+100010
+100100
+200010
+200100
+-20990
+-20900
+-20500
+-19500
+-121990
+-125990
+-125900
+-419990
+-419500
+-419000
+
+
+-------------------------------------------------------------------------------------------
Index: trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_mod.f90
===================================================================
--- trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_mod.f90	(revision 13)
+++ trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_mod.f90	(revision 14)
@@ -112,4 +112,8 @@
 !  end where
 
+!cdc test debug Hemin15
+!  where (bm(:,:).lt.-20) bm(:,:)=-10.
+
+
   acc(:,:) = 0.
   abl(:,:) = 0.
@@ -129,4 +133,8 @@
   Tann(:,:)  = tab(:,:)
 !  call lect_input(3,'Tann',1,Tann,temp_annual_file,trim(dirnameinp)//trim(runname)//'.nc')
+
+!cdc test debug Hemin15
+!  where (Tann(:,:).lt.-100.) Tann(:,:)=10.
+
 
   ta0(:,:)   = Tann(:,:)
Index: trunk/SOURCES/Greeneem_files/module_choix-greeneem.f90
===================================================================
--- trunk/SOURCES/Greeneem_files/module_choix-greeneem.f90	(revision 13)
+++ trunk/SOURCES/Greeneem_files/module_choix-greeneem.f90	(revision 14)
@@ -67,8 +67,10 @@
 
 !use dragging_vitbil
-! use dragging_hwatstream
+!use dragging_hwatstream
 !use dragging_hwat_cont
-use dragging_hwat_contmaj          !version neem sans streams
+!use dragging_hwat_contmaj          !version neem sans streams validee greeneem15
 !use dragging_neem               ! neem avec streamsn,masque offline
+use dragging_hwat_sedim          ! version C. Dumas
+
 
 !------------ spinup -----------------------------------
Index: trunk/SOURCES/Hemin15_files/fake-routines-hemin15_mod.f90
===================================================================
--- trunk/SOURCES/Hemin15_files/fake-routines-hemin15_mod.f90	(revision 14)
+++ trunk/SOURCES/Hemin15_files/fake-routines-hemin15_mod.f90	(revision 14)
@@ -0,0 +1,48 @@
+!> \file fake-routines-hemin15_mod.f90
+!! Module qui contient des routines vides
+!<
+
+!> \namespace fake_nor
+!! Module qui contient des routines vides
+!! \author ...
+!! \date ...
+!! @note Used module
+!! @note   - use module3D_phy
+!! @note   - use deform_declar
+!<
+
+
+module fake_hemin15
+
+use module3d_phy
+implicit none
+
+integer :: fake
+
+contains
+!____________________________________________________
+subroutine initial_heino
+
+if (itracebug.eq.1)  call tracebug(' Entree dans routine fake initial_heino')
+fake=1
+
+return
+end subroutine initial_heino
+
+!____________________________________________________
+subroutine track_change_T
+
+if (itracebug.eq.1)  call tracebug(' Entree dans routine fake track_change_T')
+fake=1
+return
+end subroutine track_change_T
+
+!____________________________________________________
+subroutine time_step_recul
+
+if (itracebug.eq.1)  call tracebug(' Entree dans routine fake time_step_recul')
+! quand on n'utilise pas le recul
+end subroutine time_step_recul
+
+end module fake_hemin15
+!
Index: trunk/SOURCES/Hemin15_files/geography-hemin15_mod.f90
===================================================================
--- trunk/SOURCES/Hemin15_files/geography-hemin15_mod.f90	(revision 14)
+++ trunk/SOURCES/Hemin15_files/geography-hemin15_mod.f90	(revision 14)
@@ -0,0 +1,30 @@
+!> \file geography-hemin40_mod.f90
+!! Ce module permet de choisir le domaine geographique qu'on veut utiliser
+!<
+
+!> \namespace geography
+!! Module pour choisir le domaine geographique qu'on veut utiliser
+!! \author ...
+!! \date ...
+!! @note Used module
+!! @note   - use module_nord
+!<
+
+module geography
+
+! quand on veut choisir l eurasie
+! avec plusieurs resolutions
+! use module_euras40 
+! use module_euras20 
+! use module_euras10 
+
+! quand on veut choisir le nord
+  use module_hemin15  
+
+! L'antarctique en 40 km
+! use  module_ant40
+
+! l'expérience Heino d'ISMIP
+!  use module_Heino
+
+end module geography
Index: trunk/SOURCES/Hemin15_files/lect-hemin15_mod.f90
===================================================================
--- trunk/SOURCES/Hemin15_files/lect-hemin15_mod.f90	(revision 14)
+++ trunk/SOURCES/Hemin15_files/lect-hemin15_mod.f90	(revision 14)
@@ -0,0 +1,128 @@
+module lect_topo_hemin15
+
+  use module3D_phy
+  
+    character(len=50) :: FILE1 ! CHARACTER(LEN=30) :: FILE1, FILE2,  FILE3
+    character(len=50) :: FILE2
+    character(len=80) :: filin
+    real, dimension(nx,ny,5) :: bidon          ! pour l'appel a courbure
+
+contains
+ 
+subroutine input_topo
+
+namelist/topo_file/file1,file2
+rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
+read(num_param,topo_file)
+! formats pour les ecritures dans 42
+428 format(A)
+write(num_rep_42,428)'!___________________________________________________________' 
+write(num_rep_42,428) '&topo_file                                  ! nom du bloc '
+write(num_rep_42,*)
+write(num_rep_42,*) 'file1 = ', file1
+write(num_rep_42,*) 'file2 = ', file2
+write(num_rep_42,*)'/' 
+write(num_rep_42,428) '! file1 : topo de depart'
+write(num_rep_42,428) '! file2 : topo de reference'             
+write(num_rep_42,*)
+
+!====================================== La reponse est 42 ===========
+! write(42,*)
+! write(42,*)' Fichiers en entree'
+! write(42,*)'----------------------'
+!==================================================================== 
+
+! dans param :
+!      file1=TRIM(DIRNAMEINP)//'topo-21k.g40'     ! topo LGM ICE_5G (1=topo de depart)
+!      file1=TRIM(DIRNAMEINP)//'hemin2.g40'
+!      file2=TRIM(DIRNAMEINP)//'hemin2.g40'       ! topo actuelle
+!      write(42,*) 'topo de depart', file1
+!      write(42,*) 'topo reference', file2
+
+      
+! lecture adaptee aux fichiers intercomparaison EISMINT
+       nxx=nx
+       nyy=ny
+
+!     lecture de la topo actuelle
+!     ---------------------------
+     open (20,file=TRIM(DIRNAMEINP)//file2,status='old')
+       
+     read(20,'(A80)') TITRE
+     read(20,*) NI,NJ,NXX,NYY,STEP
+     read(20,*)
+         do J=1,ny 
+          do I=1,nx
+             read (20,*)  S0(I,J),H0(I,J),BSOC0(I,J)
+             S0(i,j)=max(S0(i,j),0.)
+          end do
+        end do
+     close(20)
+
+     
+!     lecture de la topo de depart
+!     ---------------------------
+     open (20,file=TRIM(DIRNAMEINP)//file1,status='old')
+!         open (20,file='../INPUT-DATA/hemin.g50')
+     read(20,'(A80)') TITRE
+     read(20,*) NI,NJ,NXX,NYY,STEP
+     read(20,*)
+         do J=1,ny 
+          do I=1,nx
+             read (20,*) S(I,J),H(I,J),BSOC(I,J)
+          end do
+        end do
+     close(20)
+
+! calcul des courbures du socle
+
+     call courbure(nx,ny,dx,Bsoc,bidon(:,:,1),bidon(:,:,2),bidon(:,:,3), &
+          bidon(:,:,4),socle_cry,bidon(:,:,5))
+     socle_cry(:,:)=socle_cry(:,:)*dx*dx
+
+! lecture des coordonnées geographiques
+
+    filin=TRIM(DIRNAMEINP)//'coord_grisli_HN_15km.dat'
+
+! les coordonnees sont calculees en °dec avec GMT,
+! les longitudes sont comprises entre -180 et +180 (negative a l'Ouest de
+! Greenwich et positive a l'Est)
+    open(unit=2004,file=filin,iostat=ios)
+    do k=1,nx*ny
+       read(2004,*) i,j,XCC(i,j),YCC(i,j),XLONG(i,j),YLAT(i,j)
+    enddo
+    close(2004)
+    write(42,*) 'fichier grille: ', filin
+              
+    xmin=xcc(1,1)/1000.
+    ymin=ycc(1,1)/1000.
+    xmax=xcc(nx,ny)/1000.
+    ymax=ycc(nx,ny)/1000.
+                         
+! lecture du flux geothermique de Shapiro
+    open(88,file=TRIM(DIRNAMEINP)//'ijphi_hemin15.dat')
+
+    write(42,*) 'flux geothermique Shapiro : ',TRIM(DIRNAMEINP)//'ijphi_hemin15.dat'
+
+    do k=1,nx*ny
+       read(88,*) i,j,ghf(i,j)
+!        print*, i,j,ghf(i,j)
+    end do
+    close(88)
+! pour passer les flux des mW/m2 au J/m2/an      
+    ghf(:,:)=-SECYEAR/1000.*ghf(:,:)
+!     write(42,*) 'flux geothermique fixe : 55 mW/m2' 
+!     ghf(:,:)=-SECYEAR/1000.*55. !B6norcg2
+
+! print*,'lect topo'
+! print*,'shb',S(101,91),H(101,91),B(101,91)
+! print*,'shb0',S0(101,91),H0(101,91),BSOC0(101,91)
+!    Initialisation du Masque
+!------------------------------------------------
+! pour l'Hemisphere Nord mko vrai partout (version 2006)
+    MK0(:,:)=1
+
+!------------------------------------------------     
+end subroutine input_topo
+
+end module lect_topo_hemin15
Index: trunk/SOURCES/Hemin15_files/module_choix-hemin15.f90
===================================================================
--- trunk/SOURCES/Hemin15_files/module_choix-hemin15.f90	(revision 14)
+++ trunk/SOURCES/Hemin15_files/module_choix-hemin15.f90	(revision 14)
@@ -0,0 +1,140 @@
+!> \file module_choix-hemin40-0.4.f90
+!! Ce module permet de choisir quels package on veut utiliser
+!! 
+!<
+
+!> \namespace MODULE_CHOIX
+!! Ce module permet de choisir quels package on veut utiliser
+!! \author CatRitz
+!! \date ...
+!! @note Used modules:
+!! @note   - use lect_topo_hemin40    
+!! @note   - use lect_clim_act_hemin40
+!! @note   - use climat_forcage_mod
+!! @note   - use no_lakes
+!! @note   - use isostasie_mod  
+!! @note   - use deformation_mod
+!! @note   - use  prop_thermiques_ice
+!! @note   - use eau_basale
+!! @note   - use sliding_Bindschadler
+!! @note   - use dragging_hwat_contmaj
+!! @note   - use eq_ellip_sgbsv_mod
+!! @note   - use calving_frange
+!! @note   - use equat_adv_diff_2D        
+!! @note   - use bmelt_nor_regions 
+!! @note   - use fake_nor
+!! @note   - use out_profile  
+!! @note   - use out_cptr    
+!! @note   - use out_hz       
+!! @note   - use output_hemin40_mod
+!! @note   - use printtable
+!<
+
+MODULE MODULE_CHOIX
+
+!--------------Lecture topo ------------------------
+! use lect_topo_heino      ! pour Heino
+! use lect_topo_anteis    ! pour l'antarctique 40 km
+use lect_topo_hemin15    ! pour l'hemisphere nord 40 km
+!use lect_topo_eurasie   ! pour l'eurasie  (toute resolutions)
+
+!--------------Lecture climat ref ------------------
+! Le climat de ref si necessaire est maintenant lu par le module de climat directement
+!use lect_clim_act_anteis
+!use lect_clim_act_hemin40 ! pour l'hemisphere nord et l'eurasie
+! use  climat_heino
+
+!--------------Module climat ---------------
+!use climat_forcage_mois_mod ! forcage mensuel GCM 1 Snapshot Fev 2015
+use climat_Grice2sea_years_mod ! forcage avec SMB type MAR Fev2015
+
+! Anciens modules pas encore valides
+!!!!!!use climat_perturb_mod ! pour simule climat equilibre type Loveclim
+!!!!!!use climat_forcage_stat_mois_mod  ! climat constant mensuel GCMs
+!!!!!!use climat_forcage_mod ! pour cycle force GCM avec index
+!!!!!!use climat_forcage_insolation_mod ! methode JB multi-snapshots mensuelle mais pour un etat stationnaire avec correction topo GCM => GRISLI
+!!!!!!use climat_synthes_mod
+!!!!!!use climat_profil_mod
+!!!!!!use climat_regions_delta
+
+!use ablation_mod ! calcul de l'ablation (PDD ou autre methode)
+use no_ablation  ! pas de calcul de l'ablation => lecture fichier SMB (necessaire avec climat_Grice2sea_years_mod)
+
+! pas de lacs proglaciaires
+use no_lakes
+
+! suivi des traceurs ou pas ?
+!use tracer_mod  ! probablement pas compatible avec Hemin-40 actuellement
+use notracer_mod
+
+!--------------Choix isostasie----------------------
+use isostasie_mod    ! module permettant de calculer la deflexion isostasique
+!use noisostasie_mod ! module pour ne pas avoir d'isostasie
+
+
+!--------------Module Physique---------------------
+!use deformation_mod         ! module concernant les lois de deformation
+use deformation_mod_2lois   ! module concernant les lois de deformation
+
+!--------------Module propritete thermique de la glace
+use prop_thermiques_ice
+! use prop_therm_ice_heino
+
+! eau basale
+use eau_basale
+
+! sliding- dragging
+!use module sliding_vitbal
+use sliding_Bindschadler
+!use sliding_dragging_heino    ! loi de glissement heino : mis dans diffusiv : 
+
+
+!use dragging_vitbil
+!use dragging_hwatstream
+!use dragging_hwat_contmaj ! frottement basal catherine grisliv8 ice stream peu actifs
+use dragging_hwat_sedim
+!use dragging_hwat_cont
+
+!------------ spinup -----------------------------------
+use no_spinup           ! spinup=0
+!use spinup_vitbil         ! plusieurs variantes de spinup : inclue le dragging
+                          ! enlever le use dragging dans ce cas là
+
+! use eq_elliptique_mod ! ancienne version (remplimat-5)
+use eq_ellip_sgbsv_mod ! nouvelle version juillet 2008
+
+! calving
+!use calving_vincent ! ancienne version
+use calving_frange
+
+! choix resolution de la conservation de la masse
+!-----------------------------------------------
+!use equat_adv_diff_2D          ! conservation masse avec advection-diffusion
+use equat_adv_diff_2D_vect      ! le vecteur est maintenant dans l'appel
+
+
+!--------------Fusion basale------------------------
+!use bmelt_ant_regions ! pour l'Antarctique avec régions
+!use bmelt_nor_regions ! pour le nord avec régions pas encore valide pour hemin15
+use bmelt_nor_depth   ! pour le nord avec profondeur d'eau
+
+
+! fake -------fausses routines pour compatibilités
+! use fake_heino
+! use fake_ant
+use fake_hemin15
+!-------------- Outputs-----------------------------
+!
+use out_profile  !profils
+!use out_cptr     ! fichier de reprise 
+use out_hz       ! module de sorties snapshot dans un plan horizontal
+
+!-----Suivi temporel
+use output_hemin40_mod
+! use output_antarcti_mod
+! use output_heino
+!-----debugage...
+use printtable
+
+
+end module MODULE_CHOIX
Index: trunk/SOURCES/Hemin15_files/paradim-hemin15_mod.f90
===================================================================
--- trunk/SOURCES/Hemin15_files/paradim-hemin15_mod.f90	(revision 14)
+++ trunk/SOURCES/Hemin15_files/paradim-hemin15_mod.f90	(revision 14)
@@ -0,0 +1,37 @@
+!> \file paradim-hemin15_mod.f90
+!! Def des dimensions du grille hemin 15 km  
+!<
+
+!> \namespace module_nord
+!! Def des dimensions du grille hemin a 15 km  
+!! \author C. Dumas
+!! \date 15/05/2015
+!<
+
+
+module module_hemin15
+  
+  implicit none
+
+  character(len=7), parameter :: geoplace='hemin15'    ! geoplace
+  character(len=20), parameter  :: dirnameinp='../../INPUT/HEMIN15/' ! input directory
+  character(len=20), parameter  :: dirforcage='../../INPUT/Forcage/' ! input directory
+
+! dimensionnement grilles
+!-------------------------
+  integer, parameter :: nx=642           !< dimension selon x (indice i)
+  integer, parameter :: ny=642           !< dimension selon y (indice j)
+  integer, parameter :: nlmin=141        !< min(nx,ny), pour matrice elliptique  
+  integer, parameter :: nz=21            !< dimension selon z (indice k) 
+  integer, parameter :: nzm=4            !< nb de points verticaux dans le socle
+  integer, parameter :: nn=241           !< pour des sorties
+!cdc integer, parameter :: nx1=nx-1,ny1=ny-1,nz1=nz-1
+
+
+! pas horizontal
+!-----------------
+  real, parameter :: dx=15000.0             !< pas en X en m
+  real, parameter :: dy=15000.0             !< pas en Y en m
+
+
+end module module_hemin15
Index: trunk/SOURCES/Hemin15_files/track_hemin15_mod.f90
===================================================================
--- trunk/SOURCES/Hemin15_files/track_hemin15_mod.f90	(revision 14)
+++ trunk/SOURCES/Hemin15_files/track_hemin15_mod.f90	(revision 14)
@@ -0,0 +1,28 @@
+!> \file track_hemin15_mod.f90
+!! Ce module permet d'utiliser des routines de traces et de debug 
+!<
+
+!> \namespace track_debug
+!! Permet d'utiliser des routines de traces et de debug
+!! \author ...
+!! \date ...
+!! @note Ce module est different dans chaque geographie     
+!! @note faire use dans le main         
+!! @note Used module
+!! @note   - use module3d_phy
+!! @note   - use module_choix
+
+module track_debug                    ! ce module est different dans chaque geographie et 
+                                      ! permet d'utiliser des routines de traces et de debug   
+                                      ! faire use dans le main                                 
+                                     
+
+
+use module3d_phy
+use module_choix
+
+contains
+
+subroutine bidon_debug            
+end subroutine bidon_debug
+end module track_debug
Index: trunk/SOURCES/Hemin40_files/lect-clim-act-hemin40_mod.f90
===================================================================
--- trunk/SOURCES/Hemin40_files/lect-clim-act-hemin40_mod.f90	(revision 13)
+++ 	(revision )
@@ -1,51 +1,0 @@
-!> \file lect-clim-act-hemin40_mod.f90
-!! Module pour la lecture du climat
-!<
-
-!> \namespace lect_clim_act_hemin40
-!! Module pour la lecture du climat
-!! \author ...
-!! \date ...
-!! @note Used module
-!! @note   - use module3D_phy
-!<
-
-module lect_clim_act_hemin40
-
-  use module3D_phy
-  
-    CHARACTER(LEN=35) :: FILE3
-
-contains
- 
-subroutine input_climat_ref
-
-  if (geoplace.eq.'hemin40') then
-!     file3=TRIM(DIRNAMEINP)//'clim_hemin00.g40' ! climat actuel (reference)
-      file3=TRIM(DIRNAMEINP)//'climat2_hn00.g40' ! climat actuel (mod Udson...)
-  elseif (geoplace(1:5).eq.'euras') then  
-      file3=TRIM(DIRNAMEINP//'climat_'//geoplace//'.dat')
-  else
-      print*,'dans lect_clim_act_nord40 wrong geoplace :',geoplace
-  endif
-  
-        write(42,*) 'climat de reference ',file3
-
-!  lecture du fichier climat actuel
-!  --------------------------------
-    open (11,file=file3)
-       read(11,*) 
-       read(11,*)
-       read(11,*)      
-          do J=1,NY 
-             do I=1,NX
-             read(11,*) MK(I,J),PRECIP0(I,J),TA0(I,J),TJ0(I,J)
-             enddo
-          enddo
-    close(11)  
-
-!----------------------------------------------
-end subroutine input_climat_ref
-
-
-end module  lect_clim_act_hemin40
Index: trunk/SOURCES/Hemin40_files/module_choix-hemin40-0.4.f90
===================================================================
--- trunk/SOURCES/Hemin40_files/module_choix-hemin40-0.4.f90	(revision 13)
+++ trunk/SOURCES/Hemin40_files/module_choix-hemin40-0.4.f90	(revision 14)
@@ -92,5 +92,6 @@
 !use dragging_vitbil
 !use dragging_hwatstream
-use dragging_hwat_contmaj
+!use dragging_hwat_contmaj ! frottement basal catherine grisliv8 ice stream peu actifs
+use dragging_hwat_sedim
 !use dragging_hwat_cont
 
Index: trunk/SOURCES/Makefile.grisli.inc
===================================================================
--- trunk/SOURCES/Makefile.grisli.inc	(revision 13)
+++ trunk/SOURCES/Makefile.grisli.inc	(revision 14)
@@ -67,9 +67,9 @@
 	interface_input.o deformation_type.o autre_pr_temp_type.o  \
 	iso_declar_mod-0.3.o isostasie_mod-0.3.o noisostasie_mod-0.3.o \
-        deformation_mod_2lois.o deformation_mod_2lois_isotherme.o      \
+        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 \
-	prop-thermiques_mod.o  \
+	dragging_hwat_sedim_mod.o  prop-thermiques_mod.o  \
 	calving_frange.o no_calving.o no_lakes.o \
 	out_profile_mod.o printtable_mod.o mix-SIA-L1_mod.o \
@@ -101,5 +101,5 @@
 
 
-# Modules communs Post geographiques
+# Modules communs Post geographiques 
 
 mod_post_geo = step_temp_type.o
@@ -203,4 +203,10 @@
         track_hemin40_mod.o
 
+Liste_hemin15 = output_hemin40_mod.o \
+        lect-hemin15_mod.o  \
+        bmelt-hemin40-depth_mod.o \
+        fake-routines-hemin15_mod.o  module_choix-hemin15.o \
+        track_hemin15_mod.o
+
 Liste_greeneem15 = output_greeneem_mod-0.4.o \
 	lect-greeneem_mod.o \
@@ -268,4 +274,6 @@
 Dim_hemin40 = paradim-hemin40_mod.o geography-hemin40_mod.o
 
+Dim_hemin15 = paradim-hemin15_mod.o geography-hemin15_mod.o
+
 Dim_eura40 = paradim-euras40_mod.o geography-euras40_mod.o
 
@@ -690,4 +698,7 @@
 	$(FT) dragging_hwat_contmaj_mod.f90
 
+dragging_hwat_sedim_mod.o : dragging_hwat_sedim_mod.f90
+	$(FT) dragging_hwat_sedim_mod.f90
+
 dragging_calc_beta_mod.o : Antarctique_general_files/dragging_calc_beta_mod.f90
 	$(FT) Antarctique_general_files/dragging_calc_beta_mod.f90
@@ -1095,5 +1106,5 @@
 
 moy_mxmy.o : moy_mxmy.f90
-	$(FT)  moy_mxmy.f90
+	$(FT) moy_mxmy.f90
 
 neffect-0.3.o : neffect-0.3.f90 #module3d.mod
@@ -1220,7 +1231,4 @@
 	$(FT) Hemin40_files/output_hemin40_mod.f90
 
-lect-clim-act-hemin40_mod.o : Hemin40_files/lect-clim-act-hemin40_mod.f90
-	$(FT) Hemin40_files/lect-clim-act-hemin40_mod.f90
-
 lect-hemin40_mod.o : Hemin40_files/lect-hemin40_mod.f90
 	$(FT) Hemin40_files/lect-hemin40_mod.f90
@@ -1241,4 +1249,24 @@
 	$(FT) Hemin40_files/track_hemin40_mod.f90
 
+# Routines Hemin-15 :
+
+paradim-hemin15_mod.o : Hemin15_files/paradim-hemin15_mod.f90
+	$(FT) Hemin15_files/paradim-hemin15_mod.f90
+
+geography-hemin15_mod.o : Hemin15_files/geography-hemin15_mod.f90
+	$(FT) Hemin15_files/geography-hemin15_mod.f90
+
+lect-hemin15_mod.o : Hemin15_files/lect-hemin15_mod.f90
+	$(FT) Hemin15_files/lect-hemin15_mod.f90
+
+fake-routines-hemin15_mod.o : Hemin15_files/fake-routines-hemin15_mod.f90
+	$(FT) Hemin15_files/fake-routines-hemin15_mod.f90
+
+module_choix-hemin15.o : Hemin15_files/module_choix-hemin15.f90
+	$(FT) Hemin15_files/module_choix-hemin15.f90
+
+track_hemin15_mod.o : Hemin15_files/track_hemin15_mod.f90
+	$(FT) Hemin15_files/track_hemin15_mod.f90
+
 # routines Eurasie-40
 
@@ -1269,20 +1297,20 @@
 
 icetemp_declar-0.5_mod.o : Temperature-routines/icetemp_declar-0.5_mod.f90
-	$(FT)  Temperature-routines/icetemp_declar-0.5_mod.f90
+	$(FT) Temperature-routines/icetemp_declar-0.5_mod.f90
 
 interface_advec_icetemp.o :  Temperature-routines/interface_advec_icetemp.f90
-	$(FT)  Temperature-routines/interface_advec_icetemp.f90
+	$(FT) Temperature-routines/interface_advec_icetemp.f90
 
 interface_icetempmod.o :  Temperature-routines/interface_icetempmod.f90
-	$(FT)  Temperature-routines/interface_icetempmod.f90
+	$(FT) Temperature-routines/interface_icetempmod.f90
 
 interface_prop_th_icetemp.o:  Temperature-routines/interface_prop_th_icetemp.f90
-	$(FT)  Temperature-routines/interface_prop_th_icetemp.f90
+	$(FT) Temperature-routines/interface_prop_th_icetemp.f90
 
 interface_Qprod_icetemp.o:  Temperature-routines/interface_Qprod_icetemp.f90
-	$(FT)  Temperature-routines/interface_Qprod_icetemp.f90
+	$(FT) Temperature-routines/interface_Qprod_icetemp.f90
 
 interface_temp_col.o :  Temperature-routines/interface_temp_col.f90
-	$(FT)  Temperature-routines/interface_temp_col.f90
+	$(FT) Temperature-routines/interface_temp_col.f90
 
 # Routines Agrif
@@ -1318,5 +1346,5 @@
 
 graphique_L2.o : New-remplimat/graphique_L2.f90
-	$(FT)  New-remplimat/graphique_L2.f90
+	$(FT) New-remplimat/graphique_L2.f90
 
 graph-sgbsv.o : New-remplimat/graph-sgbsv.f90
@@ -1375,11 +1403,11 @@
 
 toy_declar_3D_grisli_mod.o : Recul_force_grounding_line/toy_declar_3D_grisli_mod.f90
-	$(FT)  Recul_force_grounding_line/toy_declar_3D_grisli_mod.f90
-
-declar_io_ncdf_toy_mod_fake.o :  Recul_force_grounding_line/declar_io_ncdf_toy_mod_fake.f90
-	$(FT)  Recul_force_grounding_line/declar_io_ncdf_toy_mod_fake.f90
+	$(FT) Recul_force_grounding_line/toy_declar_3D_grisli_mod.f90
+
+declar_io_ncdf_toy_mod_fake.o : Recul_force_grounding_line/declar_io_ncdf_toy_mod_fake.f90
+	$(FT) Recul_force_grounding_line/declar_io_ncdf_toy_mod_fake.f90
 
 declar_toy_retreat.o : Recul_force_grounding_line/declar_toy_retreat.f90
-	$(FT)  Recul_force_grounding_line/declar_toy_retreat.f90
+	$(FT) Recul_force_grounding_line/declar_toy_retreat.f90
 
 toy_retreat_mod.o  : Recul_force_grounding_line/toy_retreat_mod.f90
@@ -1526,4 +1554,29 @@
 	$(routine_elliptiques) $(NCDF_LIB)  $(MKL_LIB)
 
+Hemin-15 : $(Dim_hemin15) $(mod_dim_communs) \
+	$(toy_recul)  \
+	$(mod_communs) \
+	$(mod_clim_tof) \
+	$(mod_no_tracers) \
+	$(mod_ell) $(Liste_hemin15) \
+	$(mod_post_geo) \
+	$(diagnoshelf) \
+	$(Liste_Netcdf) \
+	$(routines_communes) steps_time_loop.o \
+	$(routine_elliptiques)
+
+	$(LK) -o  ../bin/Hemin-15 \
+	$(Dim_hemin15) $(mod_dim_communs) \
+	$(toy_recul)  \
+	$(mod_communs) \
+	$(mod_clim_tof) \
+	$(mod_no_tracers) \
+	$(mod_ell) $(Liste_hemin15) \
+	$(mod_post_geo) \
+	$(diagnoshelf) \
+	$(Liste_Netcdf) \
+	$(routines_communes) steps_time_loop.o \
+	$(routine_elliptiques) $(NCDF_LIB)  $(MKL_LIB)
+
 Greeneem15 : $(Dim_greeneem15) $(mod_dim_communs) \
 	$(toy_recul)  \
Index: trunk/SOURCES/Netcdf-routines/sortie_netcdf_GRISLI_mod.0.2-hassine.f90
===================================================================
--- trunk/SOURCES/Netcdf-routines/sortie_netcdf_GRISLI_mod.0.2-hassine.f90	(revision 13)
+++ trunk/SOURCES/Netcdf-routines/sortie_netcdf_GRISLI_mod.0.2-hassine.f90	(revision 14)
@@ -178,5 +178,5 @@
     !------------------------------------
     ! open(num_file,file='../'//trim(dirsource)//'/Fichiers-parametres/'//geoplace//'_TEMPS-NETCDF.dat')
-    open(num_file,file=trim(dirsource)//'/Fichiers-parametres/'//geoplace//'_TEMPS-NETCDF.dat')
+    open(num_file,file=trim(dirsource)//'/Fichiers-parametres/'//geoplace//'_TEMPS-NETCDF.dat',status='old')
 
 
Index: trunk/SOURCES/dragging_hwat_sedim_mod.f90
===================================================================
--- trunk/SOURCES/dragging_hwat_sedim_mod.f90	(revision 14)
+++ trunk/SOURCES/dragging_hwat_sedim_mod.f90	(revision 14)
@@ -0,0 +1,475 @@
+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
+
+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 :: tobmax    ! tobmax : (Pa) frottement maxi sous les streams
+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,tobmax,toblim,seuil_hwater,seuil_sedim,seuil_neff,coef_sedim,coef_gz,coef_ile    
+
+if (itracebug.eq.1) write(num_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_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,*) 'tobmax         = ', tobmax
+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) '! tobmax : (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=tobmax/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
+   
+      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)*0.1 
+! 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)*0.1 
+! 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
+
