Index: /branches/GRISLIv3/SOURCES/bmelt-beckmann-gcm_mod.f90
===================================================================
--- /branches/GRISLIv3/SOURCES/bmelt-beckmann-gcm_mod.f90	(revision 460)
+++ /branches/GRISLIv3/SOURCES/bmelt-beckmann-gcm_mod.f90	(revision 461)
@@ -35,5 +35,5 @@
   real*8 :: K_t
   real*8 :: n_tour
-  real*8, dimension (:,:,:), allocatable :: temp_ocean , mask_oce2, temp_ocean_2 !< thermal forcing            , input
+  real*8, dimension (:,:,:), allocatable :: temp_ocean     !< thermal forcing            , input
   real*8, dimension (:,:,:), allocatable :: salinity_ocean !< thermal forcing            , input
   
@@ -41,4 +41,7 @@
       
 !  integer, dimension (nx,ny) :: profondeur
+  real*8, dimension (:), allocatable :: mean_TF     !< mean thermal forcing for a given basin (non-local formula)
+  real*8, dimension (:), allocatable :: IS_area     !< ice shelf area for a given basin (non-local formula) 
+  real*8, dimension(nx,ny) :: deltaT_basin
   integer, dimension(nx,ny) :: bassin
   real*8, dimension (:), allocatable :: zoc             !< depth of oceanic levels    , input
@@ -47,4 +50,6 @@
   character(len=200) :: TO_file, SO_file, Bassin_file   ! fichiers de forcage
   logical, dimension(:), allocatable :: mask_bassin_vide  ! mask true if bassin without any GCM data
+  integer :: nonlocal_bool              !< non local ISMIP parametrisation? (1=yes)
+  integer :: grlmelt_bool               !< do we apply ocean melt at the grounded grounding line points (1=yes)
   
 contains
@@ -63,4 +68,6 @@
     character(len=100) :: SO_file                   !< files read
     character(len=100) :: Bassin_file               !< files read
+    character(len=100) :: deltaT_file               !< files read, spatial correction for Toc
+    integer :: deltaT_bool                          !< Do we apply a spatial correction for Toc (1=yes)
     integer            :: i,j,k               
     integer :: n,imin,imax,jmin,jmax
@@ -74,5 +81,6 @@
     integer :: nmax ! to avoid inifite loop 
     integer :: err
-    namelist/bmelt_beckmann_gcm_mod/K_t,TO_file,SO_file,Bassin_file,nbassins,bmelt_empty
+
+    namelist/bmelt_beckmann_gcm_mod/K_t,TO_file,SO_file,Bassin_file,nbassins,deltaT_bool,deltaT_file,bmelt_empty,nonlocal_bool,grlmelt_bool
     
     rewind(num_param)        ! loop back at the beginning of the param_list.dat file
@@ -112,18 +120,4 @@
        end if
     end if
-    if (.not.allocated(mask_oce2)) then
-      allocate(mask_oce2(nx,ny,nzoc),stat=err)
-      if (err/=0) then
-          print *,"erreur a l'allocation du tableau mask_oce2",err
-          stop 4
-       end if
-    end if
-    if (.not.allocated(temp_ocean_2)) then
-      allocate(temp_ocean_2(nx,ny,nzoc),stat=err)
-      if (err/=0) then
-          print *,"erreur a l'allocation du tableau temp_ocean_2",err
-          stop 4
-       end if
-    end if
     if (.not.allocated(salinity_ocean)) then
       allocate(salinity_ocean(nx,ny,nzoc),stat=err)
@@ -151,4 +145,18 @@
       if (err/=0) then
           print *,"erreur a l'allocation du tableau max_depth_bassin",err
+          stop 4
+       end if
+    end if
+    if (.not.allocated(mean_TF)) then
+      allocate(mean_TF(nbassins),stat=err)
+      if (err/=0) then
+          print *,"erreur a l'allocation du tableau mean_TF",err
+          stop 4
+       end if
+    end if
+    if (.not.allocated(IS_area)) then
+      allocate(IS_area(nbassins),stat=err)
+      if (err/=0) then
+          print *,"erreur a l'allocation du tableau IS_area",err
           stop 4
        end if
@@ -163,5 +171,13 @@
     call Read_Ncdf_var('so',file_inputs,tab3d)
     salinity_ocean(:,:,:)  = tab3d(:,:,:)
-    
+
+    if (deltaT_bool.eq.1) then
+       file_inputs=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(deltaT_file)
+       call Read_Ncdf_var('deltaT_basin',file_inputs,tab2d)
+       deltaT_basin(:,:) = tab2d(:,:)
+    else
+       deltaT_basin(:,:) = 0.
+    endif
+
     file_inputs=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(Bassin_file)
     call Read_Ncdf_var('mask',file_inputs,tab2d)
@@ -226,5 +242,5 @@
    where (temp_ocean(:,:,:).gt.1.e10)
       mask_gcm_init(:,:,:) = 0
-	endwhere
+   endwhere
    
 ! mask_gcm = 1 si valeur dans gcm, 0 si pas de valeur  => tableau update lorsqu'on rempli la grille   
@@ -232,5 +248,5 @@
    where (temp_ocean(:,:,:).gt.1.e10)
       mask_gcm(:,:,:) = 0
-	endwhere
+   endwhere
    !$OMP END WORKSHARE  
    
@@ -352,4 +368,6 @@
     !$OMP WORKSHARE
     ice_draft(:,:) = S(:,:)-H(:,:)-sealevel_2d(:,:)
+    mean_TF(:) = 0.d0
+    IS_area(:) = 0.d0
     !$OMP END WORKSHARE
 
@@ -370,11 +388,11 @@
                 kinf = ksup - 1
                 if ( ice_draft(i,j) .gt. zoc(1) ) then 
-                  TO_draft(i,j) = temp_ocean(i,j,1)
-                  SO_draft(i,j) = salinity_ocean(i,j,1)
-                  T_freez(i,j) = 0.0939 - 0.057 * salinity_ocean(i,j,1) + 0.000764 * zoc(1) 
+                   TO_draft(i,j) = temp_ocean(i,j,1)
+                   SO_draft(i,j) = salinity_ocean(i,j,1)
+                   T_freez(i,j) = 0.0939 - 0.057 * salinity_ocean(i,j,1) + 0.000764 * zoc(1) 
                 elseif ( ice_draft(i,j) .lt. zoc(nzoc) ) then
-                  TO_draft(i,j) = temp_ocean(i,j,nzoc)
-                  SO_draft(i,j) = salinity_ocean(i,j,nzoc)
-                  T_freez(i,j) = 0.0939 - 0.057 * salinity_ocean(i,j,nzoc) + 0.000764 * zoc(nzoc) 
+                   TO_draft(i,j) = temp_ocean(i,j,nzoc)
+                   SO_draft(i,j) = salinity_ocean(i,j,nzoc)
+                   T_freez(i,j) = 0.0939 - 0.057 * salinity_ocean(i,j,nzoc) + 0.000764 * zoc(nzoc) 
                 else
                   !TO_draft(i,j) = (   (zoc(ksup)-ice_draft(i,j)) * temp_ocean(i,j,kinf)   &
@@ -382,7 +400,12 @@
                   !SO_draft(i,j) = (   (zoc(ksup)-ice_draft(i,j)) * salinity_ocean(i,j,kinf)   &
                   !    &                   + (ice_draft(i,j)-zoc(kinf)) * salinity_ocean(i,j,ksup) ) / (zoc(ksup)-zoc(kinf))
-                  TO_draft(i,j) = temp_ocean(i,j,ksup) + ((ice_draft(i,j)-zoc(ksup))*(temp_ocean(i,j,kinf)-temp_ocean(i,j,ksup))/(zoc(kinf)-zoc(ksup))) 
-                  SO_draft(i,j) = salinity_ocean(i,j,ksup) + ((ice_draft(i,j)-zoc(ksup))*(salinity_ocean(i,j,kinf)-salinity_ocean(i,j,ksup))/(zoc(kinf)-zoc(ksup)))
-                  T_freez(i,j) = 0.0939 - 0.057 * SO_draft(i,j) + 0.000764 * ice_draft(i,j)
+                   TO_draft(i,j) = temp_ocean(i,j,ksup) + ((ice_draft(i,j)-zoc(ksup))*(temp_ocean(i,j,kinf)-temp_ocean(i,j,ksup))/(zoc(kinf)-zoc(ksup))) 
+                   SO_draft(i,j) = salinity_ocean(i,j,ksup) + ((ice_draft(i,j)-zoc(ksup))*(salinity_ocean(i,j,kinf)-salinity_ocean(i,j,ksup))/(zoc(kinf)-zoc(ksup)))
+                   T_freez(i,j) = 0.0939 - 0.057 * SO_draft(i,j) + 0.000764 * ice_draft(i,j)
+                endif
+
+                if (nonlocal_bool.eq.1) then
+                   mean_TF( bassin(i,j) ) = mean_TF( bassin(i,j) ) + mesh_area(i,j) * (TO_draft(i,j)-T_freez(i,j))
+                   IS_area( bassin(i,j) ) = IS_area( bassin(i,j) ) + mesh_area(i,j)
                 endif
                 
@@ -401,4 +424,10 @@
     enddo
     !$OMP END DO
+
+    where ( IS_area(:).gt.0.d0)     ! for all basins that have ice shelves
+       mean_TF(:) = mean_TF(:) / IS_area(:)
+    elsewhere                       ! we have no floating points over the whole basin
+       mean_TF(:) = -9999.d0
+    endwhere
     
     ! 3  - Calculation of melting rate :
@@ -409,5 +438,6 @@
     do j=1,ny
        do i=1,nx
-          if ( TO_draft(i,j) .gt. -9000.d0 ) then ! floating points
+          !if ( TO_draft(i,j) .gt. -9000.d0 ) then ! floating points
+          if ( flot(i,j) ) then ! floating points
             !if (mask_bassin_vide(bassin(i,j)) .or. (bassin(i,j).eq.1)) then ! bassins without any GCM data or bassin 1
             !if (mask_bassin_vide(bassin(i,j)) .or. (bassin(i,j).eq.1)) then !bassins without any GCM data
@@ -415,5 +445,9 @@
                bmelt(i,j) = bmelt_empty
             else  ! Bassins with GCM data
-               bmelt(i,j) = coef_OM * ( TO_draft(i,j) - T_freez(i,j) )* abs( TO_draft(i,j) - T_freez(i,j) )
+               if (nonlocal_bool.eq.1) then
+                  bmelt(i,j) = coef_OM * ( TO_draft(i,j) + deltaT_basin(i,j) - T_freez(i,j) )* abs( mean_TF(bassin(i,j)) + deltaT_basin(i,j) )
+               else
+                  bmelt(i,j) = coef_OM * ( TO_draft(i,j) + deltaT_basin(i,j) - T_freez(i,j) )* abs( TO_draft(i,j) + deltaT_basin(i,j) - T_freez(i,j) )
+               endif
             endif
           else
@@ -427,28 +461,32 @@
     do j=2,ny-1
        do i=2,nx-1
-          if ( TO_draft(i,j) .le. -9000.d0 ) then !grounded points
-             ngr=0
-             bmloc=0.d0
-             !if (TO_draft(i+1,j) .le. -9000.d0 ) then
-             if (flot(i+1,j)) then
-                ngr=ngr+1
-                bmloc=bmloc+bmelt(i+1,j)
+          !if ( TO_draft(i,j) .le. -9000.d0 ) then !grounded points
+          if ( .not. flot(i,j) ) then
+             if (grlmelt_bool.eq.1) then
+                ngr=0
+                bmloc=0.d0
+                !if (TO_draft(i+1,j) .le. -9000.d0 ) then
+                if (flot(i+1,j)) then
+                   ngr=ngr+1
+                   bmloc=bmloc+bmelt(i+1,j)
+                endif
+                !if (TO_draft(i-1,j) .le. -9000.d0 ) then
+                if (flot(i-1,j)) then
+                   ngr=ngr+1
+                   bmloc=bmloc+bmelt(i-1,j)
+                endif
+                !if (TO_draft(i,j+1) .le. -9000.d0 ) then
+                if (flot(i,j+1)) then
+                   ngr=ngr+1
+                   bmloc=bmloc+bmelt(i,j+1)
+                endif
+                !if (TO_draft(i,j-1) .le. -9000.d0 ) then
+                if (flot(i,j-1)) then
+                   ngr=ngr+1
+                   bmloc=bmloc+bmelt(i,j-1)
+                endif
+             else
+                ngr=0   ! afq -- I do not allow subgrid ocean melt
              endif
-             !if (TO_draft(i-1,j) .le. -9000.d0 ) then
-             if (flot(i-1,j)) then
-                ngr=ngr+1
-                bmloc=bmloc+bmelt(i-1,j)
-             endif
-             !if (TO_draft(i,j+1) .le. -9000.d0 ) then
-             if (flot(i,j+1)) then
-                ngr=ngr+1
-                bmloc=bmloc+bmelt(i,j+1)
-             endif
-             !if (TO_draft(i,j-1) .le. -9000.d0 ) then
-             if (flot(i,j-1)) then
-                ngr=ngr+1
-                bmloc=bmloc+bmelt(i,j-1)
-             endif
-             ngr=0   ! afq -- I do not allow subgrid ocean melt
              bmelt(i,j)=(ngr/4.)*bmloc+(1-ngr/4.)*bmelt(i,j) 
           endif
