Index: /trunk/SOURCES/Ant40_files/bmelt-ant-regions-oce_mod.f90
===================================================================
--- /trunk/SOURCES/Ant40_files/bmelt-ant-regions-oce_mod.f90	(revision 34)
+++ /trunk/SOURCES/Ant40_files/bmelt-ant-regions-oce_mod.f90	(revision 34)
@@ -0,0 +1,389 @@
+!> \file bmelt-ant-regions-oce_mod.f90
+!! Module qui calcule la fusion basale en fonction de la temp oceanique (grounded ou ice shelves).
+!<
+
+!> \namespace bmelt_ant_region_oce
+!! Module qui calcule la fusion basale (grounded ou ice shelves)
+!! \author C. Dumas
+!! \date juillet
+!! @note pour les ice shelves Antarctique, tient compte des différentes régions
+!! A choisir dans le module_choix
+!! @note Used module
+!! @note   - use module3D_phy
+!<
+MODULE  bmelt_ant_regions_oce         ! cat juillet 2005
+
+       use module3D_phy
+       use interface_input
+       use io_netcdf_grisli
+
+implicit none
+
+REAL,dimension(nx,ny) ::  bmgrz       !< fusion basale a la grounding zone
+real,dimension(nx,ny) ::  bmshelf     !< fusion basale sous shelf
+
+REAL,dimension(nx,ny) :: dist_talu    !< distance du point au talu continental
+REAL,dimension(nx,ny) :: typeshelf    !< Type de shelf Ronne->1 Ross ->2 ....
+                                      !< utilises pour moduler la fusion sous le shelf
+integer, dimension(10) :: region      !< pour écrire dans le fichier param
+character(len=30),dimension(10) :: regname !< nom des régions 
+real :: bsupshelf
+integer :: grd                        !< pour une sortie 
+
+real :: bmelt_Ross,bmgrz_Ross         !< fusion basale Ross
+real :: bmelt_FRis,bmgrz_FRis         !< fusion basale filchner-Ronne ice shelf
+real :: bmelt_Amery,bmgrz_Amery       !< Amery ice shelves
+real :: bmelt_PIG,bmgrz_PIG           !< Pine Island glacier
+real :: bmelt_Pen,bmgrz_Pen           !< Petits ice shelves Peninsule
+real :: bmelt_other,bmgrz_other       !< Autre ice shelves
+real :: bmelt_talus,bmgrz_talus       !< Au dela du talus continental
+
+real,dimension(nx,ny) :: temp_oce1, temp_oce2 ! temerature oceanique fichier 1 et 2
+real,dimension(nx,ny) :: diff_temp_oce ! variation temperature ocean => pour calcul variation fusion basale
+real :: beta_coefbmelt                ! coef variation fusion basale fonction de variation de temp oceanique (m.yr-1.K-1)
+
+
+CONTAINS
+!-------------------------------------------------------------------------------
+subroutine init_bmelt
+
+  character(len=100) :: file_temp_oce1      ! fichier de ref (ctrl) contenant la carte des temperatures de l'ocean de subsurface
+  character(len=100) :: file_temp_oce2      ! 2eme fichier : c'est la difference qui permet de calculer la variation de fusion basale
+  character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat
+  integer :: nbr_pts ! nbr de points utilises pour calcul temperature oceanique sur points sans valeurs
+  logical,dimension(nx,ny) :: mask_oce1,mask_oce2 ! masque des zones avec temp ocean (true si existe valeur)
+  real :: flpath ! distance de recherche (rayon pour les temp oceanique non presentes dans le fichier)
+  integer :: n,imin,imax,jmin,jmax
+  logical :: test ! pour stoper boucle
+
+! Cette routine fait l'initialisation pour la fusion basale.
+! Elle est appelée par inputfile-vec-0.5.f90
+
+
+namelist/bmelt_ant_reg/bmelt_Ross,bmgrz_Ross,bmelt_FRis,bmgrz_FRis,  & 
+                       bmelt_Amery,bmgrz_Amery,bmelt_PIG,bmgrz_PIG,  &
+                       bmelt_Pen,bmgrz_Pen,bmelt_other,bmgrz_other,  &       
+                       bmelt_talus,bmgrz_talus      
+
+namelist/bmelt_ocean_temp/file_temp_oce1,file_temp_oce2,beta_coefbmelt
+
+
+rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
+read(num_param,bmelt_ant_reg)
+rewind(num_param)
+read(num_param,bmelt_ocean_temp)
+
+
+!    ecriture dans le fichier parametres
+     write(num_rep_42,*)'fusion basale sous les ice shelves : bmelt-ant-regions_mod'
+     write(num_rep_42,*)'----------------------------------------------------------'
+
+! lecture du fichier contenant les distances au talu continental (m)
+      open(88,file=TRIM(DIRNAMEINP)//'distance_talu-40km.xy')
+
+      do j=1,ny
+         do i=1,nx
+            read(88,'(i3,1x,i3,1x,f10.2)') k,k,dist_talu(i,j)
+         enddo
+      enddo
+      close(88)
+
+! les bords de la grille sont mis a dist_talus=0 de force
+
+
+      dist_talu(1,:)=0.   !a gauche
+      dist_talu(2,:)=0.
+      
+      dist_talu(nx-1,:)=0. !  a droite 
+      dist_talu(nx,:)=0.  
+     
+      dist_talu(:,1)=0.    ! en bas   
+      dist_talu(:,2)=0.
+
+      dist_talu(:,ny-1)=0.    ! en haut
+      dist_talu(:,ny)=0.
+
+
+      debug_3D(:,:,32)=dist_talu(:,:)
+
+! lecture du fichier contenant les types de shelves
+!  Ronne-Flichner ->1, Ross -> 2 , Amery -> 3, 
+!  PIG-> 4, les petits shelves au dessus de PIG -> 5
+
+      open(88,file=TRIM(DIRNAMEINP)//'numer-ice-shelves-juil07.dat',status='OLD')
+      typeshelf(:,:)=100
+      do k=1,nx*ny
+        read(88,*,end=36) i,j,typeshelf(i,j)
+      end do
+36      close(88)
+
+
+
+  region(:)=0
+  regname(1)='Ronne-Fichner'
+  regname(2)='Ross'
+  regname(3)='Amery'
+  regname(4)='Pig'
+  regname(5)='Petits ice shelves peninsule'
+  regname(6)='autres ice shelves'
+  regname(7)='Au dela du talus continental'
+ 
+bms_init:    do j=1,ny
+         do i=1,nx
+
+talus:    if (dist_talu(i,j).GT.5000.) then
+            if (nint(typeshelf(i,j)).eq.1) then ! Ronne-Filchner FRis
+                    bmshelf(i,j)=bmelt_FRis         
+                    bmgrz(i,j)=bmgrz_FRis         
+
+
+               if (region(1).eq.0) then          ! pour l'impression dans 42
+                  region(1)=1
+                  write(num_rep_42,80)regname(1),bmshelf(i,j),bmgrz(i,j)
+               endif
+
+80    format(a32,' bmshelf=',f8.4,'  bmgrz=',f8.4)
+
+            else  if (nint(typeshelf(i,j)).eq.2) then ! Ross
+               bmshelf(i,j)=bmelt_Ross  
+               bmgrz(i,j)=bmgrz_Ross      
+
+               if (region(2).eq.0) then
+                  region(2)=1
+                  write(num_rep_42,80)regname(2),bmshelf(i,j),bmgrz(i,j)
+               endif
+
+
+            else  if (nint(typeshelf(i,j)).eq.3) then ! Amery
+                    bmshelf(i,j)=bmelt_Amery   
+                    bmgrz(i,j)=bmgrz_Amery          
+               if (region(3).eq.0) then
+                  region(3)=1
+                  write(num_rep_42,80)regname(3),bmshelf(i,j),bmgrz(i,j)
+               endif
+
+
+            else  if (nint(typeshelf(i,j)).eq.4) then ! Pig
+                    bmshelf(i,j)=bmelt_PIG  
+                    bmgrz(i,j)=bmgrz_PIG       
+
+               if (region(4).eq.0) then
+                  region(4)=1
+                  write(num_rep_42,80)regname(4),bmshelf(i,j),bmgrz(i,j)
+               endif
+
+            else  if (nint(typeshelf(i,j)).eq.5) then ! petits shelves peninsule
+                    bmshelf(i,j)=bmelt_Pen     
+                    bmgrz(i,j)=bmgrz_Pen       
+
+               if (region(5).eq.0) then
+                  region(5)=1
+                  write(num_rep_42,80)regname(5),bmshelf(i,j),bmgrz(i,j)
+               endif
+
+
+            else !
+                    typeshelf(i,j)=6
+                    bmshelf(i,j)=bmelt_other     
+                    bmgrz(i,j)=bmgrz_other       
+
+               if (region(6).eq.0) then
+                  region(6)=1
+                  write(num_rep_42,80)regname(6),bmshelf(i,j),bmgrz(i,j)
+               endif
+
+
+            endif
+
+
+       else   ! au dela du talus continental
+
+                    bmshelf(i,j)=bmelt_talus  
+                    bmgrz(i,j)=bmgrz_talus   
+ 
+               if (region(7).eq.0) then
+                  region(7)=1
+                  write(num_rep_42,80)regname(7),bmshelf(i,j),bmgrz(i,j)
+               endif
+
+
+            endif talus   ! fin du test sur dist_talus
+
+         enddo
+      enddo bms_init
+
+
+    call lect_input(1,'VOTEMP',1,temp_oce1,TRIM(DIRNAMEINP)//trim(file_temp_oce1),file_ncdf)    ! temp ocean ctrl
+    call lect_input(1,'VOTEMP',1,temp_oce2,TRIM(DIRNAMEINP)//trim(file_temp_oce2),file_ncdf)    ! temp ocean 
+
+
+
+
+    ! recherche de la valeur du plus proche voisin
+! initialisation
+    mask_oce1(:,:)=.true.
+    where (temp_oce1(:,:).lt.-9999.)
+       mask_oce1(:,:)=.false.
+    endwhere
+
+!    print*, 'temp_oce1(141,141)',temp_oce1(141,141)
+!    print*, 'mask_oce1(141,141)',mask_oce1(141,141)
+
+!    print*, 'temp_oce1(5,5)',temp_oce1(5,5)
+!    print*, 'mask_oce1(5,5)',mask_oce1(5,5)
+
+
+    do j=1,ny
+       do i=1,nx
+          if (.not.mask_oce1(i,j)) then
+!             print*,'boucle mask i j ',i,j,temp_oce1(i,j)
+             test=.true.
+             n=1
+             do while (test)
+                imin=max(1,i-n)
+                imax=min(NX,i+n)
+                jmin=max(1,j-n)
+                jmax=min(NY,j+n)
+                nbr_pts = count(mask_oce1(imin:imax,jmin:jmax))
+!                print*,'nbr_pts', nbr_pts
+                if (nbr_pts.GT.1) then
+                   temp_oce1(i,j)=sum(temp_oce1(imin:imax,jmin:jmax),mask=mask_oce1(imin:imax,jmin:jmax))/nbr_pts ! calcul la moyenne de temp_oce1 sur les pts non masques
+!                   print*,'fin boucle mask i j n ',i,j,n,temp_oce1(i,j)
+                   test=.false.
+                endif
+                n=n+1
+             enddo
+          endif
+       enddo
+    enddo
+
+    
+!    print*, 'temp_oce1 :',temp_oce1(5,5),TRIM(DIRNAMEINP)//trim(file_temp_oce1)
+
+!    stop
+! initialisation
+    mask_oce2(:,:)=.true.
+    where (temp_oce2(:,:).lt.-9999.)
+       mask_oce2(:,:)=.false.
+    endwhere
+
+    do j=1,ny
+       do i=1,nx
+          if (.not.mask_oce2(i,j)) then
+             test=.true.
+             n=1
+             do while (test)
+                imin=max(1,i-n)
+                imax=min(NX,i+n)
+                jmin=max(1,j-n)
+                jmax=min(NY,j+n)
+                nbr_pts = count(mask_oce2(imin:imax,jmin:jmax))
+                if (nbr_pts.GT.1) then
+                   temp_oce2(i,j)=sum(temp_oce2(imin:imax,jmin:jmax),mask=mask_oce2(imin:imax,jmin:jmax))/nbr_pts ! calcul la moyenne de temp_oce1 sur les pts non masques
+                   test=.false.
+                endif
+                n=n+1
+             enddo
+          endif
+       enddo
+    enddo
+
+
+    diff_temp_oce(:,:)=temp_oce2(:,:)-temp_oce1(:,:)
+
+    return
+  end subroutine init_bmelt
+
+
+!________________________________________________________________________________
+subroutine bmeltshelf
+
+
+! cette routine calcule la fusion basale proprement dite
+
+integer :: ngr           ! nombre de voisins flottants
+
+    do j=2,ny-1
+      do i=2,nx-1
+ 
+
+         shelf:   if (flot(i,j)) then    ! partie flottante
+            
+            bmelt(i,j)=diff_temp_oce(i,j)*beta_coefbmelt + bmshelf(i,j)
+
+            if (fbm(i,j)) then 
+               bmelt(i,j)=diff_temp_oce(i,j)*beta_coefbmelt + bmgrz(i,j)
+            endif
+
+
+! ATTENTION le bloc suivant sert a determiner la fusion basale d'equilibre
+! pour les shelves stationnaires
+
+            if (igrdline.eq.1) then
+               corrbmelt(i,j)=hdot(i,j)+bmelt(i,j)   ! le bmelt d'equilibre
+               debug_3D(i,j,28)=corrbmelt(i,j)
+            endif
+
+
+         else                   ! point posé, on compte le nombre de voisins flottants
+            ngr=0
+            if (flot(i+1,j)) ngr=ngr+1
+            if (flot(i-1,j)) ngr=ngr+1
+            if (flot(i,j+1)) ngr=ngr+1
+            if (flot(i,j-1)) ngr=ngr+1
+
+!   la fusion des points limites est une combinaison entre valeur posée et valeur flottante
+!   en fonction du nombre de points flottants
+
+            bmelt(i,j)= ngr/4.* (bmgrz(i,j)+diff_temp_oce(i,j)*beta_coefbmelt) + (1.-ngr/4.)*bmelt(i,j)
+         endif shelf
+      end do
+   end do
+
+! les bords de la grille sont mis a dist_talus=0 de force
+! attention, pourrait poser des problemes avec AGRIF
+
+   bmelt(1,:)=bmelt_talus   !a gauche
+   bmelt(2,:)=bmelt_talus
+      
+   bmelt(nx-1,:)=bmelt_talus !  a droite 
+   bmelt(nx,:)=bmelt_talus  
+     
+   bmelt(:,1)=bmelt_talus    ! en bas   
+   bmelt(:,2)=bmelt_talus
+
+   bmelt(:,ny-1)=bmelt_talus    ! en haut
+   bmelt(:,ny)=bmelt_talus
+
+!cdc test 
+!      bmelt(:,:)=diff_temp_oce(:,:)
+!      bmelt(:,:)=bmelt(:,:) + beta_coefbmelt*diff_temp_oce(:,:)
+
+
+if (igrdline.eq.1) then
+    bmelt(:,:)=0.      ! hdot donne alors le -bmelt
+endif
+
+! ecriture des bmelt pour en faire des statistiques
+!--------------------------------------------------
+!open(144,file='bmelt-test')
+!do j=1,ny
+!   do i=1,nx
+!      if (fbm(i,j)) then
+!         grd=1
+!      else
+!         grd=0
+!      endif
+!      if ((flot(i,j)).and.(H(i,j).gt.1.)) then
+!         write(144,999) i,j,grd,nint(typeshelf(i,j)), & 
+!              bmelt(i,j),bmgrz(i,j),bmshelf(i,j),hdot(i,j),dist_talu(i,j)
+!      endif
+!   end do
+!end do
+!999 format(4(i0,2x),5(e12.3,2x))
+!close(144)
+
+    return
+    end subroutine  bmeltshelf
+
+END MODULE bmelt_ant_regions_oce
Index: /trunk/SOURCES/Ant40_files/lect-anteis_mod.f90
===================================================================
--- /trunk/SOURCES/Ant40_files/lect-anteis_mod.f90	(revision 33)
+++ /trunk/SOURCES/Ant40_files/lect-anteis_mod.f90	(revision 34)
@@ -115,4 +115,9 @@
 !    endwhere
 
+
+    where (BSOC(:,:).LT.-9999.)
+       BSOC(:,:)=-9999.
+    endwhere
+
     sealevel0=0.  ! voir a passer dans le fichier parametre
     S0(:,:)=max(S0(:,:),sealevel0)   ! pour etre au niveau des mers : ATTENTION si SEALEV <0
@@ -231,10 +236,13 @@
     open(88,file=ghf_fich)
 
-    write(num_rep_42,*) 'flux geothermique Shapiro : ', TRIM(DIRNAMEINP)//'ijphi-40km-ant.dat'
-
-    do k=1,nx*ny
-       read(88,*) i,j,ghf(i,j)
-    end do
-    close(88)
+!    write(num_rep_42,*) 'flux geothermique Shapiro : ', TRIM(DIRNAMEINP)//'ijphi-40km-ant.dat'
+!    do k=1,nx*ny
+!       read(88,*) i,j,ghf(i,j)
+!    end do
+!    close(88)
+
+    call lect_input(1,'ghf',1,ghf,ghf_fich,file_ncdf)    ! socle
+
+
 
     ! pour passer les flux des mW/m2 au J/m2/an      
Index: /trunk/SOURCES/Ant40_files/module_choix-antar40-0.4.f90
===================================================================
--- /trunk/SOURCES/Ant40_files/module_choix-antar40-0.4.f90	(revision 33)
+++ /trunk/SOURCES/Ant40_files/module_choix-antar40-0.4.f90	(revision 34)
@@ -45,6 +45,7 @@
 
 !--------------Lecture climat forcage---------------
-!use climat_forcage_mois_mod ! forcage mensuel GCM 1 Snapshot Fev 2015
-use climat_perturb_mod   ! climat perturbe a reverifier Dec 2015
+use climat_forcage_mois_mod ! forcage mensuel GCM 1 Snapshot Fev 2015
+!use climat_perturb_mod   ! climat perturbe a reverifier Dec 2015
+!use climat_Grice2sea_years_mod  ! climat force par fichier SMB directement (grice2sea)
 
 !use climat_forcage_mod
@@ -105,5 +106,6 @@
 
 !--------------Fusion basale------------------------
-use bmelt_ant_regions ! pour l'Antarctique avec régions
+!use bmelt_ant_regions ! pour l'Antarctique avec régions
+use bmelt_ant_regions_oce ! pour Ant region et temp ocean
 !use bmelt_nor_regions ! pour le nord avec régions
 !use bmelt_nor_depth   ! pour le nord avec profondeur d'eau
Index: /trunk/SOURCES/Fichiers-parametres/anteis1_LISTE-VAR-NETCDF.dat
===================================================================
--- /trunk/SOURCES/Fichiers-parametres/anteis1_LISTE-VAR-NETCDF.dat	(revision 33)
+++ /trunk/SOURCES/Fichiers-parametres/anteis1_LISTE-VAR-NETCDF.dat	(revision 34)
@@ -144,4 +144,8 @@
 1	1	1			
 -----------------------------------------------------------
+ghf
+23
+1       1      -1
+-----------------------------------------------------------
 phid		
 24	
Index: /trunk/SOURCES/Fichiers-parametres/greeneem15_LISTE-VAR-NETCDF.dat
===================================================================
--- /trunk/SOURCES/Fichiers-parametres/greeneem15_LISTE-VAR-NETCDF.dat	(revision 33)
+++ /trunk/SOURCES/Fichiers-parametres/greeneem15_LISTE-VAR-NETCDF.dat	(revision 34)
@@ -144,4 +144,8 @@
 1	1	1			
 -----------------------------------------------------------
+ghf
+23
+1       1      -1
+-----------------------------------------------------------
 phid		
 24	
Index: /trunk/SOURCES/Fichiers-parametres/hemin15_LISTE-VAR-NETCDF.dat
===================================================================
--- /trunk/SOURCES/Fichiers-parametres/hemin15_LISTE-VAR-NETCDF.dat	(revision 33)
+++ /trunk/SOURCES/Fichiers-parametres/hemin15_LISTE-VAR-NETCDF.dat	(revision 34)
@@ -144,4 +144,8 @@
 1	1	1			
 -----------------------------------------------------------
+ghf
+23
+1       1      -1
+-----------------------------------------------------------
 phid		
 24	
Index: /trunk/SOURCES/Fichiers-parametres/hemin40_LISTE-VAR-NETCDF.dat
===================================================================
--- /trunk/SOURCES/Fichiers-parametres/hemin40_LISTE-VAR-NETCDF.dat	(revision 33)
+++ /trunk/SOURCES/Fichiers-parametres/hemin40_LISTE-VAR-NETCDF.dat	(revision 34)
@@ -144,4 +144,8 @@
 1	1	1			
 -----------------------------------------------------------
+ghf
+23
+1       1      -1
+-----------------------------------------------------------
 phid		
 24	
Index: /trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_mod.f90
===================================================================
--- /trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_mod.f90	(revision 33)
+++ /trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_mod.f90	(revision 34)
@@ -16,5 +16,5 @@
 
 
-use module3d_phy,only: nx,ny,S,S0,H,Bsoc,acc,abl,BM,Tann,Tjuly,Ts,time,dt,num_param,num_rep_42,dirnameinp
+use module3d_phy,only: nx,ny,S,S0,H,Bsoc,acc,abl,BM,Tann,Tjuly,Ts,time,dt,num_param,num_rep_42,dirnameinp,coefbmshelf
 !use lect_climref_Ice2sea
 use netcdf
@@ -104,5 +104,5 @@
   smb_file  = trim(dirnameinp)//trim(smb_file)
 
-  call Read_Ncdf_var('z',smb_file,tab)
+  call Read_Ncdf_var('smb',smb_file,tab)
 
 !  call lect_input(3,'smb',1,bm,smb_file,trim(dirnameinp)//trim(runname)//'.nc')
@@ -132,5 +132,5 @@
   temp_annual_file = trim(dirnameinp)//trim(temp_annual_file)
 
-  call Read_Ncdf_var('z',temp_annual_file,tab)
+  call Read_Ncdf_var('Tann',temp_annual_file,tab)
   Tann(:,:)  = tab(:,:)
 !  call lect_input(3,'Tann',1,Tann,temp_annual_file,trim(dirnameinp)//trim(runname)//'.nc')
@@ -284,4 +284,7 @@
      call init_grad_smb
   endif
+! coefbmshelf coefficient pour la fusion basale sous les ice shelves
+! ici pas d'evolution temporelle donc coefbmshelf est fixe
+  coefbmshelf=1.
 
   return
Index: /trunk/SOURCES/Makefile
===================================================================
--- /trunk/SOURCES/Makefile	(revision 33)
+++ /trunk/SOURCES/Makefile	(revision 34)
@@ -10,5 +10,5 @@
 
 # Choice for compiler (ifort = 1, gfortran = 0)
-ifort ?= 0
+ifort ?= 1
 
 # Choice for librairies (mkl_c = 1 = MKL,  = 0 = BLAS)
Index: /trunk/SOURCES/Makefile.grisli.inc
===================================================================
--- /trunk/SOURCES/Makefile.grisli.inc	(revision 33)
+++ /trunk/SOURCES/Makefile.grisli.inc	(revision 34)
@@ -214,4 +214,5 @@
         lect-anteis_mod.o \
         bmelt-ant-regions_mod.o \
+	bmelt-ant-regions-oce_mod.o \
         fake-routines-ant_mod.o module_choix-antar40-0.4.o \
         track_ant40_mod.o
@@ -515,5 +516,8 @@
 
 bmelt-ant-regions_mod.o : Ant40_files/bmelt-ant-regions_mod.f90
-	$(FT)  Ant40_files/bmelt-ant-regions_mod.f90
+	$(FT) Ant40_files/bmelt-ant-regions_mod.f90
+
+bmelt-ant-regions-oce_mod.o : Ant40_files/bmelt-ant-regions-oce_mod.f90
+	$(FT) Ant40_files/bmelt-ant-regions-oce_mod.f90
 
 bmelt-nor-regions_mod.o : Heminord_files/bmelt-nor-regions_mod.f90
Index: /trunk/SOURCES/Netcdf-routines/Description_Variables.dat
===================================================================
--- /trunk/SOURCES/Netcdf-routines/Description_Variables.dat	(revision 33)
+++ /trunk/SOURCES/Netcdf-routines/Description_Variables.dat	(revision 34)
@@ -225,4 +225,12 @@
 "celsius"	
 "Basal temperature "
+------------------------------------------------------
+"xxx"
+23      "ghf"  1
+"o"
+"geothermal_heat_flux"
+"geothermal_heat_flux"
+"m.w/mÂ²"
+"Geothermal heat flux"
 ------------------------------------------------------
 "xxx"
Index: /trunk/SOURCES/dragging_jorge_mod.f90
===================================================================
--- /trunk/SOURCES/dragging_jorge_mod.f90	(revision 34)
+++ /trunk/SOURCES/dragging_jorge_mod.f90	(revision 34)
@@ -0,0 +1,411 @@
+module dragging_jorge
+
+!jalv: module inspire par le dragging_hwat_contmaj mais simplifie
+! Definie les zones de stream avec uniquement:
+! * un critere sur la hauteur d'eau
+! * presence des sediments 
+ 
+
+  use module3d_phy
+!  use sedim_declar
+
+implicit none
+logical,dimension(nx,ny) :: fleuvemx
+logical,dimension(nx,ny) :: fleuvemy
+logical,dimension(nx,ny) :: fleuve
+logical,dimension(nx,ny) :: cote
+
+real :: seuil_sedim              ! seuil sur hwater pour avoir du glissement
+real :: valmax
+integer :: imax,jmax
+integer :: i_moins1,i_plus1,j_moins1,j_plus1
+integer :: lmax=20
+integer :: idep,jdep,iloc,jloc
+
+integer :: itestd
+
+real,dimension(nx,ny) :: betamx_temp      ! pour lissage du betamx
+real,dimension(nx,ny) :: betamy_temp      ! pour lissage du betamy
+
+real,dimension(nx,ny) :: mksedim      ! jalv, hÃ©ritÃ© de HUdson
+
+real,dimension(nx,ny) :: h_sedim      ! epaisseur de sediment
+
+
+real :: toto
+real :: neffstream  !jalv valeur en dessous de la quelle onpasse en stream
+real :: tostick   ! pour la glace posee
+real :: tob_ile    ! pour les iles
+real :: cry_lim=50.  ! courbure limite pour le suivi des fleuves
+contains
+!-------------------------------------------------------------------------------
+subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
+
+implicit none
+ 
+namelist/drag_jorge/hwatstream,cf,tobmax,toblim,seuil_sedim    
+
+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 dra_hudson_jorge
+!--------------------------------------------------------------------
+
+rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
+read(num_param,drag_jorge)
+
+write(num_rep_42,428)'!___________________________________________________________' 
+write(num_rep_42,428) '&drag_jorge              ! nom du bloc drag_jorge '
+write(num_rep_42,*)
+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_sedim         = ', seuil_sedim
+write(num_rep_42,*)'/'                            
+write(num_rep_42,428) '! hwatstream (m) :  critere de passage en stream en partant de la cote'
+write(num_rep_42,428) '!  si hwater > hwatstream '
+write(num_rep_42,428) '! cf coefficient de la loi de frottement fonction Neff'
+write(num_rep_42,428) '! seulement pour les points cotiers'
+write(num_rep_42,428) '! tobmax : (Pa) frottement maxi sous les streams '
+write(num_rep_42,428) '! toblim : (Pa)  pour les iles '
+write(num_rep_42,*)
+
+tostick=1.e5   ! valeurs pour les points non flgzmx
+tob_ile=tobmax/2.
+moteurmax=toblim
+
+
+neffstream=1.e6  !jalv
+
+!-------------------------------------------------------------------
+! masque stream
+
+      mstream_mx(:,:)=1
+      mstream_my(:,:)=1
+
+! coefficient permettant de modifier le basal drag.
+      drag_mx(:,:)=1.
+      drag_my(:,:)=1.
+
+
+
+
+
+! lecture du fichier sediment :
+
+open(20,file=TRIM(DIRNAMEINP)//'Sediment_map_Gabi_Laske/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.1)  ! pas d'epaisseurs de sediments negatives               
+
+
+where (h_sedim(:,:).gt.250) 
+      mksedim(:,:)=2
+endwhere
+
+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 gzmx initial a ete calculÃ© dans flottab pour les points a cheval sur la grounding line
+
+
+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
+
+
+! critere de passage en stream
+do  j=2,ny-1  
+   do i=2,nx-1  
+      if (((hwater(i,j).gt.hwatstream).or.((mksedim(i,j).eq.2).and.(hwater(i,j).ge.seuil_sedim))).and.(.not.flot(i,j))) then
+
+!IF ((neffmx(i,j).LT.neffstream).AND.(neffmy(i,j).LT.neffstream).and.(.not.flot(i,j))) THEN
+!        IF ((((neffmx(i,j).LT.neffstream).AND.(neffmy(i,j).LT.neffstream)).or.((mksedim(i,j).eq.2.).and.(hwater(i,j).gt.seuil_sedim))).AND.(.NOT.flot(i,j))) THEN
+
+         fleuve(i,j)=.true.
+      else
+         fleuve(i,j)=.false.
+      endif
+   enddo
+enddo
+
+do j=1,ny-1
+   do i=1,nx-1
+      if (fleuve(i,j)) then
+         fleuvemx(i,j)=.true.
+         fleuvemx(i+1,j)=.true.
+         fleuvemy(i,j)=.true.
+         fleuvemy(i,j+1)=.true.
+      end if
+   end do
+end do
+
+! pas de fleuve dans les endroits trop en pente
+
+!fleuvemx(:,:)=fleuvemx(:,:).and.(abs(rog*Hmx(:,:)*sdx(:,:)).lt.toblim)
+!fleuvemy(:,:)=fleuvemy(:,:).and.(abs(rog*Hmy(:,:)*sdy(:,:)).lt.toblim)
+
+!.....................
+
+!jalv
+!call detect_assym(nx,ny,0,76,1,0,1,0,betamx,itestd)
+!if (itestd.gt.0) then
+!   write(6,*) 'avant calcul betax asymetrie sur betamx  pour time=',time
+!   stop
+!else
+!   write(6,*) 'avant calcul betax pas d asymetrie sur betamx  pour time=',time
+!end if
+
+!call detect_assym(nx,ny,0,76,1,0,1,0,betamy,itestd)
+!if (itestd.gt.0) then
+!    write(6,*) 'avant calcul betay asymetrie sur betamy  pour time=',time
+!    stop
+!else
+!    write(6,*) 'avant calcul betay pas d asymetrie sur betamy  pour time=',time
+!end if
+
+
+
+!...........................................!   
+! frottement selon x (avant dans neffect)   !    ! test jalv: tous les cas ont la meme loi pour beta: beta=cf*neff
+!...........................................!    ! mais le sueil maxi est different pour chaque cas
+ 
+
+
+DO j=1,ny
+  DO i=1,nx
+
+    IF ((.NOT.ilemx(i,j)).AND.(fleuvemx(i,j))) gzmx(i,j)=.TRUE. 
+    
+    IF (cotemx(i,j)) THEN                        ! point cotier: frott maxi=80, mini=8
+
+!        betamx(i,j)=cf*neffmx(i,j)*1.      ! clim3_39
+!        betamx(i,j)=cf*neffmx(i,j)*1.      ! clim3_36
+!        betamx(i,j)=cf*neffmx(i,j)*1.      ! clim3_35
+!         betamx(i,j)=cf*neffmx(i,j)*1.    ! clim3_22         
+!         betamx(i,j)=cf*neffmx(i,j)*0.1  !clim3_21 
+        betamx(i,j)=cf*neffmx(i,j)*0.1 ! standar
+!         betamx(i,j)=cf*neffmx(i,j)*0.1/h_sedim(i,j) !clim3_20
+!         betamx(i,j)=min(betamx(i,j),80.)
+!         betamx(i,j)=max(betamx(i,j),8.)
+    ELSE IF ((gzmx(i,j)).AND.(.NOT.cotemx(i,j))) THEN  ! stream normal: frott maxi=160, mini=16                  
+
+!        betamx(i,j)=cf*neffmx(i,j)*50.      ! clim3_39
+!         betamx(i,j)=cf*neffmx(i,j)*500. ! clim3_36
+!         betamx(i,j)=cf*neffmx(i,j)*100. ! clim3_22
+!         betamx(i,j)=cf*neffmx(i,j)*10.   ! clim3_21
+        betamx(i,j)=cf*neffmx(i,j)*10. ! standar
+!         betamx(i,j)=cf*neffmx(i,j)*10./h_sedim(i,j)
+!         betamx(i,j)=min(betamx(i,j),160.)           
+!         betamx(i,j)=max(betamx(i,j),16.)
+         
+        IF ((mksedim(i,j).EQ.2).AND.(hwater(i,j).GE.seuil_sedim)) THEN  ! jalv: stream  en zone sediment
+
+!            betamx(i,j)=cf*neffmx(i,j)*5.      ! clim3_39
+!            betamx(i,j)=cf*neffmx(i,j)*50. ! clim3_36
+!            betamx(i,j)=cf*neffmx(i,j)*10. ! clim3_22
+!            betamx(i,j)=cf*neffmx(i,j)*1. !clim3_21
+            betamx(i,j)=cf*neffmx(i,j)*1. ! standar
+!            betamx(i,j)=cf*neffmx(i,j)*1./h_sedim(i,j)
+!            betamx(i,j)=min(betamx(i,j),120.)          !stream sediment: frott maxi=120, mini=12 
+!           betamx(i,j)=max(betamx(i,j),12.)   
+        ENDIF
+
+    ELSE IF (ilemx(i,j)) THEN                    ! ile : frott maxi=40, mini=4
+
+!        betamx(i,j)=cf*neffmx(i,j)*0.1      ! clim3_39
+!         betamx(i,j)=cf*neffmx(i,j)*0.1   ! clim3_22
+        betamx(i,j)=cf*neffmx(i,j)*0.1 ! standar
+!         betamx(i,j)=cf*neffmx(i,j)*0.1  !clim3_21
+!         betamx(i,j)=cf*neffmx(i,j)*0.1/h_sedim(i,j)
+!         betamx(i,j)=min(betamx(i,j),40.)
+!         betamx(i,j)=max(betamx(i,j),4.)
+    ELSE IF (flotmx(i,j)) THEN                   ! flottant : frott 0.
+        betamx(i,j)=0.
+    ELSE                                         ! grounded, SIA
+!         betamx(i,j)=tostick                      ! frottement glace posee (1 bar) divise par mil !jalv
+!         betamx(i,j)=200.                          ! jalv: le cas restant glace posee a priori: frott=150
+        betamx(i,j)=cf*neffmx(i,j)*1000.
+!        betamx(i,j)=cf*neffmx(i,j)*1000./h_sedim(i,j)
+    ENDIF
+
+!         if (cf*neffmx(i,j).gt.1500.) then
+!            gzmx(i,j)=.false.
+!            betamx(i,j)=tostick
+!         endif
+
+   end do
+end do
+
+
+!................................................!
+!...test..jalv..lissage..des..beta...............!
+!...on..calcule..la..moyenne..ponderee..des......!
+!...points..aux..alentours.......................!
+
+
+!do j=1,ny
+!   do i=3,nx-2
+
+!      betamx_temp(i,j)=(betamx(i,j)+(0.3*betamx(i-1,j))+(0.2*betamx(i-2,j))+(0.3*betamx(i+1,j))+(0.2*betamx(i+2,j)))/2
+     
+!   end do
+!end do
+
+
+
+!betamx(:,:)=betamx_temp(:,:)
+
+ 
+
+!...........................................!
+! frottement selon y (avant dans neffect)   !   ! test jalv: tous les cas ont la meme loi pour beta: beta=cf*neff
+!...........................................!   ! mais le seuil maxi est different pour chaque cas
+
+
+do j=1,ny                 
+   do i=1,nx   
+
+      if ((.not.ilemy(i,j)).and.(fleuvemy(i,j))) gzmy(i,j)=.true. 
+
+      
+      IF (cotemy(i,j)) THEN                        ! point cotier: froot maxi=40, mini=4
+
+!          betamy(i,j)=cf*neffmy(i,j)*1.      ! clim3_39
+!          betamy(i,j)=cf*neffmy(i,j)*1.      ! clim3_36        
+!          betamy(i,j)=cf*neffmy(i,j)*1.      ! clim3_35
+!          betamy(i,j)=cf*neffmy(i,j)*1.    ! clim3_22
+!          betamy(i,j)=cf*neffmy(i,j)*0.1  !clim3_21
+          betamy(i,j)=cf*neffmy(i,j)*0.1  ! standar
+!         betamy(i,j)=cf*neffmy(i,j)*0.1/h_sedim(i,j) 
+!         betamy(i,j)=min(betamy(i,j),80.)
+!         betamy(i,j)=max(betamy(i,j),8.)
+      else if ((gzmy(i,j)).and.(.not.cotemy(i,j))) then  ! stream normal: frott maxi=80, mini=8
+
+!          betamy(i,j)=cf*neffmy(i,j)*50.      ! clim3_39
+!          betamy(i,j)=cf*neffmy(i,j)*500. ! clim3_36
+!          betamy(i,j)=cf*neffmy(i,j)*100. ! clim3_22
+!          betamy(i,j)=cf*neffmy(i,j)*10. !clim3_21
+          betamy(i,j)=cf*neffmy(i,j)*10. ! standar
+!          betamy(i,j)=cf*neffmy(i,j)*10./h_sedim(i,j)
+!          betamy(i,j)=min(betamy(i,j),160.)            
+!          betamy(i,j)=max(betamy(i,j),16.)
+
+         if ((mksedim(i,j).eq.2).and.(hwater(i,j).ge.seuil_sedim)) then              ! jalv: stream en zone sediment
+
+!             betamy(i,j)=cf*neffmy(i,j)*5.      ! clim3_39
+!            betamy(i,j)=cf*neffmy(i,j)*50. ! clim3_36
+!            betamy(i,j)=cf*neffmy(i,j)*10. ! clim3_22
+!             betamy(i,j)=cf*neffmy(i,j)*1. !clim3_21  
+             betamy(i,j)=cf*neffmy(i,j)*1. ! standar 
+!            betamy(i,j)=cf*neffmy(i,j)*1./h_sedim(i,j)
+!            betamy(i,j)=min(betamy(i,j),120.)      !stream sediment: frott maxi=60, mini=0.
+!            betamy(i,j)=max(betamy(i,j),12.)  
+         endif
+      
+      else if (ilemy(i,j)) then                   ! ile : frott maxi=20, mini=2
+
+!          betamy(i,j)=cf*neffmy(i,j)*0.1      ! clim3_39
+!         betamy(i,j)=cf*neffmy(i,j)*0.1   ! clim3_22  
+          betamy(i,j)=cf*neffmy(i,j)*0.1   ! standar
+!         betamy(i,j)=cf*neffmy(i,j)*0.1  !clim3_21
+!         betamy(i,j)=cf*neffmy(i,j)*0.1/h_sedim(i,j)
+!         betamy(i,j)=min(betamy(i,j),40.)
+!         betamy(i,j)=max(betamy(i,j),4.)
+      else if (flotmy(i,j)) then                  ! flottant : frott 0
+         betamy(i,j)=0.
+      else                                        ! grounded, SIA
+!         betamy(i,j)=tostick/1000                ! frottement glace posee divise par mil !jalv
+!         betamy(i,j)=200.                         ! jalv: le cas restant glace posee a priori: frott maxi 200
+         betamy(i,j)=cf*neffmy(i,j)*1000.
+!         betamy(i,j)=cf*neffmy(i,j)*1000./h_sedim(i,j)
+      endif
+
+!         if (cf*neffmy(i,j).gt.1500.) then
+!            gzmy(i,j)=.false.
+!            betamy(i,j)=tostick
+!         endif
+
+   end do
+end do
+
+
+
+!................................................!
+!...test..jalv..lissage..des..beta...............!
+!...on..calcule..la..moyenne..ponderee..des......!
+!...points..aux..alentours.......................!
+
+
+!do j=3,ny-2
+!   do i=1,nx
+      
+!      betamy_temp(i,j)=(betamy(i,j)+(0.3*betamy(i,j-1))+(0.2*betamy(i,j-2))+(0.3*betamy(i,j+1))+(0.2*betamy(i,j+2)))/2
+   
+!   end do
+!end do
+
+!betamy(:,:)=betamy_temp(:,:)
+
+
+
+!jalv
+!call detect_assym(nx,ny,0,76,1,0,1,0,betamx,itestd)
+!if (itestd.gt.0) then
+!   write(6,*) 'apres dragging asymetrie sur betamx  pour time=',time
+!   stop
+!else
+!   write(6,*) 'apres dragging pas d asymetrie sur betamx  pour time=',time
+!end if
+
+!call detect_assym(nx,ny,0,76,1,0,1,0,betamy,itestd)
+!if (itestd.gt.0) then
+!    write(6,*) 'apres dragging asymetrie sur betamy  pour time=',time
+!    stop
+!else
+!    write(6,*) 'apres dragging pas d asymetrie sur betamy  pour time=',time
+!end if
+
+
+
+
+
+return
+end subroutine dragging
+
+end module dragging_jorge
Index: /trunk/SOURCES/initial-0.3.f90
===================================================================
--- /trunk/SOURCES/initial-0.3.f90	(revision 33)
+++ /trunk/SOURCES/initial-0.3.f90	(revision 34)
@@ -71,5 +71,5 @@
   endif
 
-!cdc call sortie_ncdf_cat
+  call sortie_ncdf_cat
 
   ! calcul de Hmx et Hmy -> shift=-1, dim=1 -> H(i-1,j)
Index: /trunk/SOURCES/lect-eis.f90
===================================================================
--- /trunk/SOURCES/lect-eis.f90	(revision 33)
+++ /trunk/SOURCES/lect-eis.f90	(revision 34)
@@ -18,5 +18,5 @@
 
   implicit none
-  character(len=17)  :: DIRNAMEINP
+  character(len=20)  :: DIRNAMEINP
   integer nx,ny,i,j,kdum
   real data(nx,ny)
