Index: trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90	(revision 3411)
+++ trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90	(revision 3412)
@@ -54,8 +54,9 @@
                               fast,fasthaze,haze,metcloud,monoxcloud,&
                               n2cond,nearn2cond,noseason_day,conservn2, &
-                              convergeps,kbo,triton,paleo,paleoyears, &
+                              convergeps,kbo,triton,paleo,paleoyears,glaflow, &
                               carbox, methane,&
                               oldplutovdifc,oldplutocorrk,oldplutosedim, &
-                              aerohaze,haze_proffix,source_haze,&
+                              aerohaze,haze_proffix,source_haze, tsurfmax, &
+                              albmin_ch4, &
                               season, sedimentation,generic_condensation, &
                               specOLR, &
@@ -291,4 +292,7 @@
       REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
 
+      REAL,SAVE :: ptime0    ! store the first time
+      REAL dstep
+      REAL,SAVE :: glastep=20   ! step in pluto day to spread glacier
 
 
@@ -617,4 +621,8 @@
          albedo_snow_SPECTV(:)=0.0
          albedo_n2_ice_SPECTV(:)=0.0
+
+         ptime0=ptime
+         write (*,*) 'In physiq ptime0 =', ptime
+
          call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_n2_ice_SPECTV)
 
@@ -759,4 +767,9 @@
       taux(1:ngrid) = 0.0
       tauy(1:ngrid) = 0.0
+
+      if (conservn2) then
+         write(*,*) 'conservn2 iniloop'
+         call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
+      endif
 
       zday=pday+ptime ! Compute time, in sols (and fraction thereof).
@@ -848,4 +861,10 @@
        enddo
       endif
+
+      if (conservn2) then
+         write(*,*) 'conservn2 thermo'
+         call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
+      endif
+
 !---------------------------------
 ! II. Compute radiative tendencies
@@ -1094,4 +1113,8 @@
 !!      call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
 
+   if (conservn2) then
+      write(*,*) 'conservn2 radiat'
+      call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
+   endif
 
 !  --------------------------------------------
@@ -1229,4 +1252,9 @@
       endif ! end of 'calldifv'
 
+      if (conservn2) then
+        write(*,*) 'conservn2 diff'
+        call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)+ &
+                                       dqsurf(:,1)*ptimestep)
+      endif
 
 !-------------------
@@ -1322,4 +1350,9 @@
       endif  ! end of 'n2cond'
 
+      if (conservn2) then
+       write(*,*) 'conservn2 n2cond'
+       call testconservmass(ngrid,nlayer,pplev(:,1)+ &
+           pdpsrf(:)*ptimestep,qsurf(:,1)+dqsurf(:,1)*ptimestep)
+      endif
 
 !---------------------------------------------
@@ -1592,6 +1625,6 @@
       if (conservn2) then
          write(*,*) 'conservn2 tracer'
-         ! call testconservmass(ngrid,nlayer,pplev(:,1)+   &
-         !    pdpsrf(:)*ptimestep,qsurf(:,1))
+         call testconservmass(ngrid,nlayer,pplev(:,1)+   &
+            pdpsrf(:)*ptimestep,qsurf(:,1))
       endif
 
@@ -1646,9 +1679,21 @@
 
 
-      ! ! Increment surface temperature
-      ! if(ok_slab_ocean)then  !AF24: removed
-
+! VII.1 Increment surface temperature
+!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
-      ! Compute soil temperatures and subsurface heat flux.
+
+      ! Prevent surface (.e.g. non volatile ch4) to exceed max temperature
+      ! Lellouch et al., 2000,2011
+      IF (tsurfmax) THEN
+        DO ig=1,ngrid
+         if (albedo_equivalent(ig).gt.albmin_ch4.and. &
+                           qsurf(ig,igcm_n2).eq.0.) then
+              tsurf(ig)=min(tsurf(ig),54.)
+         endif
+        ENDDO
+      ENDIF
+
+! VII.2 Compute soil temperatures and subsurface heat flux.
+!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       if (callsoil) then
          call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat,   &
@@ -1656,6 +1701,9 @@
       endif
 
-
-!       if (ok_slab_ocean) then !AF24: removed
+      ! ! For output :
+      ! tidat_out(:,:)=0.
+      ! DO l=1,min(nlayermx,nsoilmx)
+      !    tidat_out(:,l)=tidat(:,l)
+      ! ENDDO
 
       ! Test energy conservation
@@ -1665,4 +1713,49 @@
       endif
 
+
+
+!   VII.3 multiply tendencies of cond/subli for paleo loop only in the
+!       last Pluto year of the simulation
+!       Year day must be adapted in the startfi for each object
+!       Paleo uses year_day to calculate the annual mean tendancies
+!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      IF (paleo) then
+         if (zday.gt.day_ini+ptime0+nday-year_day) then
+            DO iq=1,nq
+             DO ig=1,ngrid
+               qsurfyear(ig,iq)=qsurfyear(ig,iq)+ &
+                             (qsurf(ig,iq)-qsurf1(ig,iq))  !kg m-2 !ptimestep
+             ENDDO
+            ENDDO
+         endif
+      endif
+
+!   VII.4 Glacial flow at each timestep glastep or at lastcall
+!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      IF (fast.and.glaflow) THEN
+         if ((mod(zday-day_ini-ptime0,glastep)).lt.1. &
+                                                   .or.lastcall) then
+           IF (lastcall) then
+            dstep=mod(zday-day_ini-ptime0,glastep)*daysec
+           else
+            dstep=glastep*daysec
+           endif
+           zdqflow(:,:)=qsurf(:,:)
+           IF (paleo) then
+             call spreadglacier_paleo(ngrid,nq,qsurf, &
+                                    phisfinew,dstep,tsurf)
+           else
+             call spreadglacier_simple(ngrid,nq,qsurf,dstep)
+           endif
+           zdqflow(:,:)=(zdqflow(:,:)-qsurf(:,:))/dstep
+
+           if (conservn2) then
+            write(*,*) 'conservn2 glaflow'
+            call testconservmass(ngrid,nlayer,pplev(:,1)+ &
+            pdpsrf(:)*ptimestep,qsurf(:,1))
+           endif
+
+         endif
+      ENDIF
 
 !---------------------------------------------------
@@ -1816,5 +1909,4 @@
 
 
-      ! ! Test for water conservation. !AF24: removed
 
       ! Calculate RH_generic (Generic Relative Humidity) for diagnostic.
@@ -2271,9 +2363,9 @@
 
          do iq=1,nq
-          ! call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
+          call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
           !  call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
           !                   'kg m^-2',2,qsurf_hist(1,iq) )
-          !  call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
-          !                  'kg m^-2',2,qcol(1,iq) )
+           call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
+                           'kg m^-2',2,qcol(1,iq) )
 !          call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
 !                          'kg m^-2',2,qsurf(1,iq) )
