Index: /branches/iLoveclim/SOURCES/3D-physique-gen_mod.f90
===================================================================
--- /branches/iLoveclim/SOURCES/3D-physique-gen_mod.f90	(revision 90)
+++ /branches/iLoveclim/SOURCES/3D-physique-gen_mod.f90	(revision 91)
@@ -253,4 +253,5 @@
   real (kind = 8),dimension(nx,ny) :: calvingGRIS
   real (kind = 8),dimension(nx,ny),target :: calvin_GRIS ! calving flux[m3/s] given to CLIO
+  real (kind = 8),dimension(nx,ny,20),target :: bmshelfCLIO  ! basal melting given by CLIO  
   real, dimension(nx,ny) :: coef_drag   !< coefficient de la loi de friction non lineaire : depend de la valeur de alpha_drag
                                         !< si alpha_drag = 1, coef_drag = drag_centre  
Index: /branches/iLoveclim/SOURCES/Hemin40_files/module_choix-hemin40-0.4.f90
===================================================================
--- /branches/iLoveclim/SOURCES/Hemin40_files/module_choix-hemin40-0.4.f90	(revision 90)
+++ /branches/iLoveclim/SOURCES/Hemin40_files/module_choix-hemin40-0.4.f90	(revision 91)
@@ -117,6 +117,6 @@
 !use bmelt_ant_regions ! pour l'Antarctique avec régions
 !use bmelt_nor_regions ! pour le nord avec régions
-use bmelt_nor_depth   ! pour le nord avec profondeur d'eau
-
+!use bmelt_nor_depth   ! pour le nord avec profondeur d'eau
+use bmelt_clio_coupl_mod ! coupled version: CLIO bm rates (needs SHELFMELT==1 iLOVECLIM flag)
 
 ! fake -------fausses routines pour compatibilités
Index: /branches/iLoveclim/SOURCES/Makefile.grisli-gfortran.inc
===================================================================
--- /branches/iLoveclim/SOURCES/Makefile.grisli-gfortran.inc	(revision 90)
+++ /branches/iLoveclim/SOURCES/Makefile.grisli-gfortran.inc	(revision 91)
@@ -82,5 +82,6 @@
 	conserv-mass-adv-diff_sept2009_mod.o  \
 	eaubasale-0.5_mod.o spinup_mod.o no_spinup_mod.o \
-	tridiag_mod-0.3.o out_horiz_mod.o bmelt-seuil-profondeur_mod.o 
+	tridiag_mod-0.3.o out_horiz_mod.o \
+	bmelt_clio_coupl_mod.o bmelt-seuil-profondeur_mod.o 
 
 mod_communs_RGL =   util_recovery.o \
Index: /branches/iLoveclim/SOURCES/Makefile.grisli.inc
===================================================================
--- /branches/iLoveclim/SOURCES/Makefile.grisli.inc	(revision 90)
+++ /branches/iLoveclim/SOURCES/Makefile.grisli.inc	(revision 91)
@@ -81,5 +81,6 @@
 	conserv-mass-adv-diff_sept2009_mod.o  \
 	eaubasale-0.5_mod.o spinup_mod.o no_spinup_mod.o \
-	tridiag_mod-0.3.o out_horiz_mod.o bmelt-seuil-profondeur_mod.o 
+	tridiag_mod-0.3.o out_horiz_mod.o \
+	bmelt_clio_coupl_mod.o bmelt-seuil-profondeur_mod.o 
 
 mod_communs_RGL =   util_recovery.o \
Index: /branches/iLoveclim/SOURCES/bmelt_clio_coupl_mod.f90
===================================================================
--- /branches/iLoveclim/SOURCES/bmelt_clio_coupl_mod.f90	(revision 91)
+++ /branches/iLoveclim/SOURCES/bmelt_clio_coupl_mod.f90	(revision 91)
@@ -0,0 +1,236 @@
+!> \file bmelt_clio_coupl_mod.f90
+!! Read the bmelt computed in CLIO (thersf.f) for the 20 vertical layers
+!<
+
+!> \namespace  BMELT_CLIO_COUPL_MOD
+!! Read the bmelt computed in CLIO and compute grounded basal melt
+!! @note Only for coupled version
+!! \author Aurelien Quiquet
+!! \date October 2016
+!! @note Used module
+!! @note   - use module3D_phy
+!<
+
+
+MODULE  BMELT_CLIO_COUPL_MOD
+
+
+  USE module3D_phy,only:nx,ny,i,j,ro,row,h,hdot,bmelt,corrbmelt,igrdline,flot,fbm,bmshelfCLIO
+
+implicit none
+
+real, parameter       :: bmgrz_fact = 2. !< bmshelf/bmgrz ratio (fixed)
+
+REAL,dimension(nx,ny) ::  bmgrz     !< fusion basale a la grounding zone
+real,dimension(nx,ny) ::  bmshelf   !< fusion basale sous shelf
+
+real,dimension(20)    ::  z_CLIO    !< depth of the CLIO layers (center of layer)
+real,dimension(20)    ::  dz_CLIO   !< thickness of the CLIO layers
+
+
+CONTAINS
+!-------------------------------------------------------------------------------
+
+!> SUBROUTINE: init_bmelt
+!! This routine does the initialisation of basal melting rates of ice shelves
+!<
+subroutine init_bmelt
+
+! local variables:
+integer :: locdepth ! local depth of the shelf
+integer :: noc      ! loop integer on the vertical oceanic layers
+
+! This routine does the initialisation of basal melting rates of ice shelves
+!  -  z_CLIO
+!  - dz_CLIO
+!  - bmshelf and bmgrz
+! Called by initial-0.3.f90
+
+!    ecriture dans le fichier parametres
+     write(42,*)'fusion basale sous les ice shelves : bmelt from CLIO   '
+     write(42,*)'-------------------------------------------------------'
+
+
+!the depths of the centers of CLIO vertical layers:
+     z_CLIO(20)  = 5.00
+     z_CLIO(19)  = 15.98
+     z_CLIO(18)  = 29.17
+     z_CLIO(17)  = 45.20
+     z_CLIO(16)  = 64.96
+     z_CLIO(15)  = 89.75
+     z_CLIO(14)  = 121.52
+     z_CLIO(13)  = 163.28
+     z_CLIO(12)  = 219.86
+     z_CLIO(11)  = 299.26
+     z_CLIO(10)  = 415.07
+     z_CLIO(9)   = 588.88
+     z_CLIO(8)   = 850.19
+     z_CLIO(7)   = 1225.11
+     z_CLIO(6)   = 1717.90
+     z_CLIO(5)   = 2307.36
+     z_CLIO(4)   = 2963.25
+     z_CLIO(3)   = 3661.11
+     z_CLIO(2)   = 4385.22
+     z_CLIO(1)   = 5126.18
+     
+!the thicknesses of the CLIO vertical layers:
+     dz_CLIO(20) = 10.00
+     dz_CLIO(19) = 11.96
+     dz_CLIO(18) = 14.42
+     dz_CLIO(17) = 17.64
+     dz_CLIO(16) = 21.88
+     dz_CLIO(15) = 27.70
+     dz_CLIO(14) = 35.84
+     dz_CLIO(13) = 47.68
+     dz_CLIO(12) = 65.48
+     dz_CLIO(11) = 93.38
+     dz_CLIO(10) = 138.18
+     dz_CLIO(9)  = 209.44
+     dz_CLIO(8)  = 313.18
+     dz_CLIO(7)  = 436.66
+     dz_CLIO(6)  = 548.92
+     dz_CLIO(5)  = 630.00
+     dz_CLIO(4)  = 681.78
+     dz_CLIO(3)  = 713.94
+     dz_CLIO(2)  = 734.28
+     dz_CLIO(1)  = 747.64
+
+     
+     do j=1,ny
+        do i=1,nx
+
+! Init: everywhere to surface CLIO melt.
+! In CLIO, when we don't have ocean, I assumed a very high bm
+!          In this case, in GRISLI, we put constant values (0.2)
+           if (bmshelfCLIO(i,j,20).lt.999d0) then! 
+              bmshelf(i,j) = real(bmshelfCLIO(i,j,20))
+              bmgrz(i,j) = bmshelf(i,j) * bmgrz_fact
+           else
+              bmshelf(i,j) = 0.2
+              bmgrz(i,j) = bmshelf(i,j) * bmgrz_fact
+           endif
+
+! Now, we look at the depth of the ice shelves to use the right bm
+           !if (flot(i,j).and.(H(i,j).gt.1.)) then
+           if (flot(i,j)) then
+              
+              bmelt(i,j) = bmshelf(i,j) ! init
+              
+              locdepth=ro/row*H(i,j)
+              do noc=20,2,-1
+                 if ( locdepth .gt. z_CLIO(noc)+dz_CLIO(noc)/2. ) then
+                    if (bmshelfCLIO(i,j,noc+1).lt.999d0) then
+                       bmelt(i,j) = real(bmshelfCLIO(i,j,noc+1))
+                    endif
+                 else
+                    exit
+                 endif
+              enddo
+              
+           endif
+            
+        enddo
+     enddo
+
+
+      return
+      end subroutine init_bmelt
+
+!________________________________________________________________________________
+
+!> SUBROUTINE: bmeltshelf
+!! This routine computes the actual basal melting rates
+!<
+subroutine bmeltshelf
+
+! local variables:
+integer :: locdepth ! local depth of the shelf
+integer :: noc      ! loop integer on the vertical oceanic layers
+integer :: ngr      ! number of floating points, neighbours of a grounded points
+real    :: bmsum    ! temporary bm 
+
+do j=1,ny
+   do i=1,nx
+
+! Init: everywhere to surface CLIO melt.
+! In CLIO, when we don't have ocean, I assumed a very high bm
+!          In this case, in GRISLI, we put constant values (0.2)
+      if (bmshelfCLIO(i,j,20).lt.999d0) then! 
+         bmshelf(i,j) = real(bmshelfCLIO(i,j,20))
+         bmgrz(i,j) = bmshelf(i,j) * bmgrz_fact
+      else
+         bmshelf(i,j) = 0.2
+         bmgrz(i,j) = bmshelf(i,j) * bmgrz_fact
+      endif
+
+! Now, we look at the depth of the ice shelves to use the right bm
+      if (flot(i,j)) then    ! partie flottante
+
+         bmelt(i,j) = bmshelf(i,j) ! init
+         
+         locdepth=ro/row*H(i,j)
+         do noc=20,2,-1
+            if ( locdepth.gt. z_CLIO(noc)+dz_CLIO(noc)/2. ) then
+               if (bmshelfCLIO(i,j,noc+1).lt.999d0) then
+                  bmelt(i,j) = real(bmshelfCLIO(i,j,noc+1))
+               endif
+            else
+               exit
+            endif
+         enddo
+         
+         if (fbm(i,j))  bmelt(i,j)=bmgrz(i,j)
+
+
+! ATTENTION LE BLOC SUIVANT SERT A FAIRE DES ICE SHELVES STATIONNAIRES
+! igrdline est défini dans itnitial1
+! afq -- not tested for coupled applications...
+         if (igrdline.eq.1) then
+            corrbmelt(i,j)=corrbmelt(i,j)+hdot(i,j)*0.8
+            bmelt(i,j)=bmelt(i,j)+corrbmelt(i,j)
+         endif
+
+      endif ! on flot
+   enddo
+enddo
+
+do j=2,ny-1
+   do i=2,nx-1
+      
+      if (.not.flot(i,j)) then   ! grounded point, we account for the floating neighbours
+
+         bmsum=0.
+         ngr=0
+         if (flot(i+1,j)) then
+            ngr=ngr+1
+            bmsum= bmsum+bmelt(i+1,j)
+         endif
+         if (flot(i-1,j)) then
+            ngr=ngr+1
+            bmsum= bmsum+bmelt(i-1,j)
+         endif
+         if (flot(i,j+1)) then
+            ngr=ngr+1
+            bmsum= bmsum+bmelt(i,j+1)
+         endif
+         if (flot(i,j-1)) then
+            ngr=ngr+1
+            bmsum= bmsum+bmelt(i,j-1)
+         endif
+
+! Grounding point basal melting rate is a combined effect of grounded and floating values:
+
+           bmelt(i,j)= ngr/4.*bmgrz(i,j)+(1.-ngr/4.)*bmsum
+
+        endif
+
+     end do
+  end do
+
+
+    return
+    end subroutine  bmeltshelf
+
+
+
+END MODULE  BMELT_CLIO_COUPL_MOD
Index: /branches/iLoveclim/SOURCES/main3D-0.4-40km.f90
===================================================================
--- /branches/iLoveclim/SOURCES/main3D-0.4-40km.f90	(revision 90)
+++ /branches/iLoveclim/SOURCES/main3D-0.4-40km.f90	(revision 91)
@@ -137,4 +137,6 @@
   end do time_loop
 
+  bmshelfCLIO(:,:,:) = 0d0
+
   if (itracebug.eq.1)  call tracebug('dans main avant call out_recovery ')
   call out_recovery(iout)
Index: /branches/iLoveclim/SOURCES/steps_time_loop.f90
===================================================================
--- /branches/iLoveclim/SOURCES/steps_time_loop.f90	(revision 90)
+++ /branches/iLoveclim/SOURCES/steps_time_loop.f90	(revision 91)
@@ -246,5 +246,7 @@
      call ablation
 
-
+     ! oceanic forcing
+     !=====================
+     call bmeltshelf  ! afq -- icetemp is not called before the 4th dtt, so a first call to bmelt here
 
      ! SIA dynamic velocities  (deformation + loi de glissement )
@@ -289,4 +291,7 @@
         call icetemp
         ! subroutines pour le calcul de la fusion basale
+        ! afq -- 28-10-2016: there is a reset of bmelt in temp_col in the
+        ! current version of the code, so shelf and grounded melt have to 
+        ! be called here. TODO: remove reset in icetemp routines??
         call bmeltshelf
         call bmelt_grounded
