Index: trunk/LMDZ.GENERIC/README
===================================================================
--- trunk/LMDZ.GENERIC/README	(revision 1470)
+++ trunk/LMDZ.GENERIC/README	(revision 1477)
@@ -1083,2 +1083,4 @@
   aerosol properties files and surface files in datadir) still works.
 
+== 21/09/2015 == MT
+- Cleanup of physiq.F90 and changed condense_cloud.F90 to condense_co2.F90
Index: trunk/LMDZ.GENERIC/libf/phystd/callsedim.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callsedim.F	(revision 1470)
+++ trunk/LMDZ.GENERIC/libf/phystd/callsedim.F	(revision 1477)
@@ -107,5 +107,5 @@
       do iq=1,nq
        if((radius(iq).gt.1.e-9).and.(iq.ne.igcm_co2_ice)) then   
-!         (no sedim for gases, and co2_ice sedim is done in condense_cloud)      
+!         (no sedim for gases, and co2_ice sedim is done in condense_co2)      
 
 ! store locally updated tracers
Index: trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90	(revision 1470)
+++ 	(revision )
@@ -1,522 +1,0 @@
-      subroutine condense_cloud(ngrid,nlayer,nq,ptimestep, &
-          pcapcal,pplay,pplev,ptsrf,pt, &
-          pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, &
-          piceco2,psolaralb,pemisurf, &
-          pdtc,pdtsrfc,pdpsrf,pduc,pdvc, &
-          pdqc)
-
-      use radinc_h, only : naerkind
-      use gases_h, only: gfrac, igas_co2
-      use radii_mod, only : co2_reffrad
-      use aerosol_mod, only : iaero_co2
-      USE surfdat_h, only: albedodat, albedice, emisice, emissiv
-      USE comgeomfi_h, only: lati
-      USE tracer_h, only: noms, rho_co2
-      use comcstfi_mod, only: g, r, cpp
-
-      implicit none
-
-!==================================================================
-!     Purpose
-!     -------
-!     Condense and/or sublime CO2 ice on the ground and in the
-!     atmosphere, and sediment the ice.
-!  
-!     Inputs
-!     ------
-!     ngrid                 Number of vertical columns
-!     nlayer                Number of layers
-!     pplay(ngrid,nlayer)   Pressure layers
-!     pplev(ngrid,nlayer+1) Pressure levels 
-!     pt(ngrid,nlayer)      Temperature (in K)
-!     ptsrf(ngrid)          Surface temperature
-!     
-!     pdt(ngrid,nlayer)   Time derivative before condensation/sublimation of pt
-!     pdtsrf(ngrid)         Time derivative before condensation/sublimation of ptsrf
-!     pqsurf(ngrid,nq)      Sedimentation flux at the surface (kg.m-2.s-1)
-!     
-!     Outputs
-!     -------
-!     pdpsrf(ngrid)       \  Contribution of condensation/sublimation
-!     pdtc(ngrid,nlayer)  /  to the time derivatives of Ps, pt, and ptsrf
-!     pdtsrfc(ngrid)     /
-!     
-!     Both
-!     ----
-!     piceco2(ngrid)        CO2 ice at the surface (kg/m2)
-!     psolaralb(ngrid)      Albedo at the surface
-!     pemisurf(ngrid)       Emissivity of the surface
-!
-!     Authors
-!     ------- 
-!     Francois Forget (1996)
-!     Converted to Fortran 90 and slightly modified by R. Wordsworth (2009)
-!     Includes simplifed nucleation by J. Leconte (2011)
-!     
-!==================================================================
-
-!-----------------------------------------------------------------------
-!     Arguments
-
-      INTEGER,INTENT(IN) :: ngrid
-      INTEGER,INTENT(IN) :: nlayer
-      INTEGER,INTENT(IN) :: nq
-      REAL,INTENT(IN) :: ptimestep 
-      REAL,INTENT(IN) :: pcapcal(ngrid)
-      REAL,INTENT(IN) :: pplay(ngrid,nlayer)
-      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1)
-      REAL,INTENT(IN) :: ptsrf(ngrid)
-      REAL,INTENT(IN) :: pt(ngrid,nlayer)
-      REAL,INTENT(IN) :: pphi(ngrid,nlayer)
-      REAL,INTENT(IN) :: pdt(ngrid,nlayer)
-      REAL,INTENT(IN) :: pdu(ngrid,nlayer)
-      REAL,INTENT(IN) :: pdv(ngrid,nlayer)
-      REAL,INTENT(IN) :: pdtsrf(ngrid)
-      REAL,INTENT(IN) :: pu(ngrid,nlayer)
-      REAL,INTENT(IN) :: pv(ngrid,nlayer)
-      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
-      REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq)
-      REAL,INTENT(INOUT) :: piceco2(ngrid)
-      REAL,INTENT(OUT) :: psolaralb(ngrid)
-      REAL,INTENT(OUT) :: pemisurf(ngrid)
-      REAL,INTENT(OUT) :: pdtc(ngrid,nlayer)
-      REAL,INTENT(OUT) :: pdtsrfc(ngrid)
-      REAL,INTENT(OUT) :: pdpsrf(ngrid)
-      REAL,INTENT(OUT) :: pduc(ngrid,nlayer)
-      REAL,INTENT(OUT) :: pdvc(ngrid,nlayer)
-      REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq)
-
-!-----------------------------------------------------------------------
-!     Local variables
-
-      INTEGER l,ig,icap,ilay,it,iq
-
-      REAL reffrad(ngrid,nlayer) ! radius (m) of the co2 ice particles
-      REAL*8 zt(ngrid,nlayer)
-      REAL zq(ngrid,nlayer,nq)
-      REAL zcpi
-      REAL ztcond (ngrid,nlayer)
-      REAL ztnuc (ngrid,nlayer)
-      REAL ztcondsol(ngrid) 
-      REAL zdiceco2(ngrid)
-      REAL zcondicea(ngrid,nlayer), zcondices(ngrid)
-      REAL zfallice(ngrid), Mfallice(ngrid) 
-      REAL zmflux(nlayer+1)
-      REAL zu(nlayer),zv(nlayer)
-      REAL ztsrf(ngrid)
-      REAL ztc(nlayer), ztm(nlayer+1) 
-      REAL zum(nlayer+1) , zvm(nlayer+1)
-      LOGICAL condsub(ngrid)
-      REAL subptimestep
-      Integer Ntime
-      real masse (ngrid,nlayer), w(ngrid,nlayer,nq)
-      real wq(ngrid,nlayer+1)
-      real vstokes,reff
-
-!     Special diagnostic variables
-      real tconda1(ngrid,nlayer)
-      real tconda2(ngrid,nlayer)
-      real zdtsig (ngrid,nlayer)
-      real zdt (ngrid,nlayer)
-
-!-----------------------------------------------------------------------
-!     Saved local variables
-
-      REAL,SAVE :: latcond=5.9e5
-      REAL,SAVE :: ccond
-      REAL,SAVE :: cpice=1000.
-      REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: emisref
-!$OMP THREADPRIVATE(latcond,ccond,cpice,emisref)
-
-      LOGICAL,SAVE :: firstcall=.true.
-!$OMP THREADPRIVATE(firstcall)
-      REAL,EXTERNAL :: SSUM
-
-      REAL,EXTERNAL :: CBRT
-
-      INTEGER,SAVE :: i_co2ice=0      ! co2 ice
-!$OMP THREADPRIVATE(i_co2ice)
-      CHARACTER(LEN=20) :: tracername ! to temporarily store text
-
-      integer igas
-
-      real ppco2
-
-!-----------------------------------------------------------------------
-!     Initializations
-
-      pdqc(1:ngrid,1:nlayer,1:nq)=0
-      pdtc(1:ngrid,1:nlayer)=0
-      zq(1:ngrid,1:nlayer,1:nq)=0
-      zt(1:ngrid,1:nlayer)=0
-
-!     Initialisation
-      IF (firstcall) THEN
-
-         ALLOCATE(emisref(ngrid)) !! this should be deallocated in lastcall...
-
-         ! find CO2 ice tracer
-         do iq=1,nq
-            tracername=noms(iq)
-            if (tracername.eq."co2_ice") then
-               i_co2ice=iq
-            endif
-         enddo
-         
-        write(*,*) "condense_cloud: i_co2ice=",i_co2ice       
-
-        if((i_co2ice.lt.1))then
-           print*,'In condens_cloud but no CO2 ice tracer, exiting.'
-           print*,'Still need generalisation to arbitrary species!'
-           stop
-        endif
-
-         ccond=cpp/(g*latcond)
-         print*,'In condens_cloud: ccond=',ccond,' latcond=',latcond
-
-!          Prepare special treatment if gas is not pure CO2
-             !if (addn2) then
-             !   m_co2   = 44.01E-3 ! CO2 molecular mass (kg/mol)   
-             !   m_noco2 = 28.02E-3 ! N2 molecular mass (kg/mol)  
-!               Compute A and B coefficient use to compute
-!               mean molecular mass Mair defined by
-!               1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
-!               1/Mair = A*q(ico2) + B
-             !   A = (1/m_co2 - 1/m_noco2)
-             !   B = 1/m_noco2
-             !endif
-
-!          Minimum CO2 mixing ratio below which mixing occurs with layer above: 
-           !qco2min =0.75  
-
-           firstcall=.false.
-      ENDIF
-      zcpi=1./cpp
-
-!-----------------------------------------------------------------------
-!     Calculation of CO2 condensation / sublimation 
-!
-!     Variables used:
-!     piceco2(ngrid)       amount of co2 ice on the ground  (kg/m2)
-!     zcondicea(ngrid,l)   condensation rate in layer l     (kg/m2/s)
-!     zcondices(ngrid)     condensation rate on the ground  (kg/m2/s)
-!     zfallice(ngrid)      flux of ice falling on surface   (kg/m2/s)
-!     pdtc(ngrid,nlayer) dT/dt due to phase changes       (K/s)
-     
-
-!     Tendencies initially set to 0 (except pdtc)
-      DO l=1,nlayer
-         DO ig=1,ngrid
-            zcondicea(ig,l) = 0.
-            pduc(ig,l)  = 0
-            pdvc(ig,l)  = 0
-            pdqc(ig,l,i_co2ice)  = 0
-         END DO
-      END DO
-      
-      DO ig=1,ngrid
-         Mfallice(ig) = 0.
-         zfallice(ig) = 0.
-         zcondices(ig) = 0.
-         pdtsrfc(ig) = 0.
-         pdpsrf(ig) = 0.
-         condsub(ig) = .false.
-         zdiceco2(ig) = 0.
-      ENDDO
-
-!-----------------------------------------------------------------------
-!     Atmospheric condensation
-
-
-!     Compute CO2 Volume mixing ratio
-!     -------------------------------
-!      if (addn2) then
-!         DO l=1,nlayer
-!            DO ig=1,ngrid
-!              qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep
-!!             Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
-!              mmean=1/(A*qco2 +B)
-!              vmr_co2(ig,l) = qco2*mmean/m_co2 
-!            ENDDO
-!         ENDDO
-!      else
-!         DO l=1,nlayer
-!            DO ig=1,ngrid
-!              vmr_co2(ig,l)=0.5
-!            ENDDO
-!         ENDDO
-!      end if
-
-!     Forecast the atmospheric frost temperature 'ztcond' and nucleation temperature 'ztnuc'
-      DO l=1,nlayer
-         DO ig=1,ngrid
-            ppco2=gfrac(igas_CO2)*pplay(ig,l)
-            call get_tcond_co2(ppco2,ztcond(ig,l))
-            call get_tnuc_co2(ppco2,ztnuc(ig,l))
-         ENDDO
-      ENDDO
-      
-!     Initialize zq and zt at the beginning of the sub-timestep loop
-      DO l=1,nlayer
-         DO ig=1,ngrid
-            zt(ig,l)=pt(ig,l)
-            zq(ig,l,i_co2ice)=pq(ig,l,i_co2ice)
-            IF( zq(ig,l,i_co2ice).lt.-1.e-6 ) THEN
-               print*,'Uh-oh, zq = ',zq(ig,l,i_co2ice),'at ig,l=',ig,l
-               if(l.eq.1)then
-                  print*,'Perhaps the atmosphere is collapsing on surface...?'
-               endif
-            END IF
-         ENDDO
-      ENDDO
-
-!     Calculate the mass of each atmospheric layer (kg.m-2)
-      do  ilay=1,nlayer
-         DO ig=1,ngrid
-            masse(ig,ilay)=(pplev(ig,ilay) - pplev(ig,ilay+1)) /g
-         end do
-      end do
-
-!     -----------------------------------------------
-!     START CONDENSATION/SEDIMENTATION SUB-TIME LOOP
-!     -----------------------------------------------
-      Ntime =  20               ! number of sub-timestep 
-      subptimestep = ptimestep/float(Ntime)           
-
-      DO it=1,Ntime
-
-!     Add the tendencies from other physical processes at each subtimstep
-         DO l=1,nlayer
-            DO ig=1,ngrid
-               zt(ig,l)   = zt(ig,l)   + pdt(ig,l)   * subptimestep
-               zq(ig,l,i_co2ice) = zq(ig,l,i_co2ice) + pdq(ig,l,i_co2ice) * subptimestep
-            END DO
-         END DO
-
-
-!     Gravitational sedimentation
-            
-!     sedimentation computed from radius computed from q in module radii_mod
-	 call co2_reffrad(ngrid,nlayer,nq,zq,reffrad)
-	 
-         do  ilay=1,nlayer
-            DO ig=1,ngrid
-
-               reff = reffrad(ig,ilay)
-
-               call stokes                      &
-                   (pplev(ig,ilay),pt(ig,ilay), &
-                    reff,vstokes,rho_co2)
-
-               !w(ig,ilay,i_co2ice) = 0.0
-               w(ig,ilay,i_co2ice) = vstokes *  subptimestep * &
-                   pplev(ig,ilay)/(r*pt(ig,ilay))
-
-            end do
-         end do
-
-!     Computing q after sedimentation
-
-         call vlz_fi(ngrid,nlayer,zq(1,1,i_co2ice),2.,masse,w(1,1,i_co2ice),wq)
-
-
-!     Progressively accumulating the flux to the ground
-!     Mfallice is the total amount of ice fallen to the ground
-         DO ig=1,ngrid
-            Mfallice(ig) =  Mfallice(ig) + wq(ig,i_co2ice)
-         end do
-
-
-!     Condensation / sublimation in the atmosphere
-!     --------------------------------------------
-!     (calculation of zcondicea, zfallice and pdtc)
-!     (MODIFICATIONS FOR EARLY MARS: falling heat neglected, condensation
-!     of CO2 into tracer i_co2ice)
-
-         DO l=nlayer , 1, -1
-            DO ig=1,ngrid
-               pdtc(ig,l)=0.
-
-
-   !     ztcond-> ztnuc in test beneath to nucleate only when super saturation occurs(JL 2011)
-               IF ((zt(ig,l).LT.ztnuc(ig,l)).or.(zq(ig,l,i_co2ice).gt.1.E-10)) THEN 
-                  condsub(ig)=.true.
-                  pdtc(ig,l)   = (ztcond(ig,l) - zt(ig,l))/subptimestep
-                  pdqc(ig,l,i_co2ice) = pdtc(ig,l)*ccond*g
-
-!     Case when the ice from above sublimes entirely
-                  IF ((zq(ig,l,i_co2ice).lt.-pdqc(ig,l,i_co2ice)*subptimestep) &
-                      .AND. (zq(ig,l,i_co2ice).gt.0)) THEN
-
-                     pdqc(ig,l,i_co2ice) = -zq(ig,l,i_co2ice)/subptimestep
-                     pdtc(ig,l)   =-zq(ig,l,i_co2ice)/(ccond*g*subptimestep)
-
-                  END IF
-
-!     Temperature and q after condensation
-                  zt(ig,l)   = zt(ig,l)   + pdtc(ig,l)   * subptimestep
-                  zq(ig,l,i_co2ice) = zq(ig,l,i_co2ice) + pdqc(ig,l,i_co2ice) * subptimestep
-               END IF
-
-            ENDDO
-         ENDDO
-      ENDDO                     ! end of subtimestep loop
-
-!     Computing global tendencies after the subtimestep
-      DO l=1,nlayer
-         DO ig=1,ngrid
-            pdtc(ig,l) = &
-              (zt(ig,l) - (pt(ig,l) + pdt(ig,l)*ptimestep))/ptimestep
-            pdqc(ig,l,i_co2ice) = &
-              (zq(ig,l,i_co2ice)-(pq(ig,l,i_co2ice)+pdq(ig,l,i_co2ice)*ptimestep))/ptimestep
-         END DO
-      END DO
-      DO ig=1,ngrid
-         zfallice(ig) = Mfallice(ig)/ptimestep
-      END DO
-
-
-!-----------------------------------------------------------------------
-!     Condensation/sublimation on the ground
-!     (calculation of zcondices and pdtsrfc)
-
-!     forecast of ground temperature ztsrf and frost temperature ztcondsol
-      DO ig=1,ngrid
-         ppco2=gfrac(igas_CO2)*pplay(ig,1)
-         call get_tcond_co2(ppco2,ztcondsol(ig))
-         
-         ztsrf(ig) = ptsrf(ig)
-
-         if((ztsrf(ig).le.ztcondsol(ig)+2.0).and.(ngrid.eq.1))then
-            print*,'CO2 is condensing on the surface in 1D. This atmosphere is doomed.'
-            print*,'T_surf = ',ztsrf,'K'
-            print*,'T_cond = ',ztcondsol,'K'
-            open(116,file='surf_vals.out')
-            write(116,*) 0.0, pplev(1,1), 0.0, 0.0
-            close(116)
-            call abort
-         endif
-
-         ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
-
-      ENDDO
-     
-      DO ig=1,ngrid
-         IF(ig.GT.ngrid/2+1) THEN
-            icap=2
-         ELSE
-            icap=1
-         ENDIF
-
-!     Loop over where we have condensation / sublimation
-         IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR.           &        ! ground condensation
-                    (zfallice(ig).NE.0.) .OR.              &        ! falling snow
-                    ((ztsrf(ig) .GT. ztcondsol(ig)) .AND.  &        ! ground sublimation
-                    ((piceco2(ig)+zfallice(ig)*ptimestep) .NE. 0.))) THEN
-            condsub(ig) = .true.
-
-!     Condensation or partial sublimation of CO2 ice
-            zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
-                /(latcond*ptimestep)         
-            pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
-
-!     If the entire CO_2 ice layer sublimes
-!     (including what has just condensed in the atmosphere)
-            IF((piceco2(ig)/ptimestep+zfallice(ig)).LE. &
-                -zcondices(ig))THEN
-               zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig)
-               pdtsrfc(ig)=(latcond/pcapcal(ig))*       &
-                   (zcondices(ig))
-            END IF
-
-!     Changing CO2 ice amount and pressure
-
-            zdiceco2(ig) = zcondices(ig) + zfallice(ig)
-            piceco2(ig)  = piceco2(ig) + zdiceco2(ig)*ptimestep
-            pdpsrf(ig)   = -zdiceco2(ig)*g
-
-            IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
-               PRINT*,'STOP in condens'
-               PRINT*,'condensing more than total mass'
-               PRINT*,'Grid point ',ig
-               PRINT*,'Ps = ',pplev(ig,1)
-               PRINT*,'d Ps = ',pdpsrf(ig)
-               STOP
-            ENDIF
-         END IF
-      ENDDO
-
-!     Surface albedo and emissivity of the ground below the snow (emisref)
-!     --------------------------------------------------------------------
-      DO ig=1,ngrid
-         IF(lati(ig).LT.0.) THEN
-            icap=2 ! Southern Hemisphere
-         ELSE
-            icap=1 ! Nortnern hemisphere
-         ENDIF
-
-         if(.not.piceco2(ig).ge.0.) THEN
-            if(piceco2(ig).le.-1.e-8) print*,   &
-                'WARNING in condense_co2cloud: piceco2(',ig,')=', piceco2(ig)
-            piceco2(ig)=0.
-         endif
-         if (piceco2(ig).gt.0) then
-            psolaralb(ig) = albedice(icap)
-            emisref(ig)   = emisice(icap)
-         else
-            psolaralb(ig) = albedodat(ig)
-            emisref(ig)   = emissiv
-            pemisurf(ig)  = emissiv
-         end if
-      end do
-
-      return
-    end subroutine condense_cloud
-
-!-------------------------------------------------------------------------
-    subroutine get_tcond_co2(p,tcond)
-!   Calculates the condensation temperature for CO2
-
-      implicit none
-
-      real p, peff, tcond
-      real, parameter :: ptriple=518000.0
-
-      peff=p
-
-      if(peff.lt.ptriple)then
-         tcond = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula
-      else
-         tcond = 684.2-92.3*log(peff)+4.32*log(peff)**2 
-         ! liquid-vapour transition (based on CRC handbook 2003 data)
-      endif
-      return
-
-    end subroutine get_tcond_co2
-    
-    
-    
- 
-!-------------------------------------------------------------------------
-    subroutine get_tnuc_co2(p,tnuc)
-!   Calculates the nucleation temperature for CO2, based on a simple super saturation criterion 
-!     (JL 2011)
-
-      use callkeys_mod, only: co2supsat
-
-      implicit none
-
-      real p, peff, tnuc
-      real, parameter :: ptriple=518000.0
-
-      peff=p/co2supsat
-
-      if(peff.lt.ptriple)then
-         tnuc = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula
-      else
-         tnuc = 684.2-92.3*log(peff)+4.32*log(peff)**2 
-         ! liquid-vapour transition (based on CRC handbook 2003 data)
-      endif
-      return
-
-    end subroutine get_tnuc_co2
Index: trunk/LMDZ.GENERIC/libf/phystd/condense_co2.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/condense_co2.F90	(revision 1477)
+++ trunk/LMDZ.GENERIC/libf/phystd/condense_co2.F90	(revision 1477)
@@ -0,0 +1,522 @@
+      subroutine condense_co2(ngrid,nlayer,nq,ptimestep, &
+          pcapcal,pplay,pplev,ptsrf,pt, &
+          pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, &
+          piceco2,psolaralb,pemisurf, &
+          pdtc,pdtsrfc,pdpsrf,pduc,pdvc, &
+          pdqc)
+
+      use radinc_h, only : naerkind
+      use gases_h, only: gfrac, igas_co2
+      use radii_mod, only : co2_reffrad
+      use aerosol_mod, only : iaero_co2
+      USE surfdat_h, only: albedodat, albedice, emisice, emissiv
+      USE comgeomfi_h, only: lati
+      USE tracer_h, only: noms, rho_co2
+      use comcstfi_mod, only: g, r, cpp
+
+      implicit none
+
+!==================================================================
+!     Purpose
+!     -------
+!     Condense and/or sublime CO2 ice on the ground and in the
+!     atmosphere, and sediment the ice.
+!  
+!     Inputs
+!     ------
+!     ngrid                 Number of vertical columns
+!     nlayer                Number of layers
+!     pplay(ngrid,nlayer)   Pressure layers
+!     pplev(ngrid,nlayer+1) Pressure levels 
+!     pt(ngrid,nlayer)      Temperature (in K)
+!     ptsrf(ngrid)          Surface temperature
+!     
+!     pdt(ngrid,nlayer)   Time derivative before condensation/sublimation of pt
+!     pdtsrf(ngrid)         Time derivative before condensation/sublimation of ptsrf
+!     pqsurf(ngrid,nq)      Sedimentation flux at the surface (kg.m-2.s-1)
+!     
+!     Outputs
+!     -------
+!     pdpsrf(ngrid)       \  Contribution of condensation/sublimation
+!     pdtc(ngrid,nlayer)  /  to the time derivatives of Ps, pt, and ptsrf
+!     pdtsrfc(ngrid)     /
+!     
+!     Both
+!     ----
+!     piceco2(ngrid)        CO2 ice at the surface (kg/m2)
+!     psolaralb(ngrid)      Albedo at the surface
+!     pemisurf(ngrid)       Emissivity of the surface
+!
+!     Authors
+!     ------- 
+!     Francois Forget (1996)
+!     Converted to Fortran 90 and slightly modified by R. Wordsworth (2009)
+!     Includes simplifed nucleation by J. Leconte (2011)
+!     
+!==================================================================
+
+!-----------------------------------------------------------------------
+!     Arguments
+
+      INTEGER,INTENT(IN) :: ngrid
+      INTEGER,INTENT(IN) :: nlayer
+      INTEGER,INTENT(IN) :: nq
+      REAL,INTENT(IN) :: ptimestep 
+      REAL,INTENT(IN) :: pcapcal(ngrid)
+      REAL,INTENT(IN) :: pplay(ngrid,nlayer)
+      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1)
+      REAL,INTENT(IN) :: ptsrf(ngrid)
+      REAL,INTENT(IN) :: pt(ngrid,nlayer)
+      REAL,INTENT(IN) :: pphi(ngrid,nlayer)
+      REAL,INTENT(IN) :: pdt(ngrid,nlayer)
+      REAL,INTENT(IN) :: pdu(ngrid,nlayer)
+      REAL,INTENT(IN) :: pdv(ngrid,nlayer)
+      REAL,INTENT(IN) :: pdtsrf(ngrid)
+      REAL,INTENT(IN) :: pu(ngrid,nlayer)
+      REAL,INTENT(IN) :: pv(ngrid,nlayer)
+      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
+      REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq)
+      REAL,INTENT(INOUT) :: piceco2(ngrid)
+      REAL,INTENT(OUT) :: psolaralb(ngrid)
+      REAL,INTENT(OUT) :: pemisurf(ngrid)
+      REAL,INTENT(OUT) :: pdtc(ngrid,nlayer)
+      REAL,INTENT(OUT) :: pdtsrfc(ngrid)
+      REAL,INTENT(OUT) :: pdpsrf(ngrid)
+      REAL,INTENT(OUT) :: pduc(ngrid,nlayer)
+      REAL,INTENT(OUT) :: pdvc(ngrid,nlayer)
+      REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq)
+
+!-----------------------------------------------------------------------
+!     Local variables
+
+      INTEGER l,ig,icap,ilay,it,iq
+
+      REAL reffrad(ngrid,nlayer) ! radius (m) of the co2 ice particles
+      REAL*8 zt(ngrid,nlayer)
+      REAL zq(ngrid,nlayer,nq)
+      REAL zcpi
+      REAL ztcond (ngrid,nlayer)
+      REAL ztnuc (ngrid,nlayer)
+      REAL ztcondsol(ngrid) 
+      REAL zdiceco2(ngrid)
+      REAL zcondicea(ngrid,nlayer), zcondices(ngrid)
+      REAL zfallice(ngrid), Mfallice(ngrid) 
+      REAL zmflux(nlayer+1)
+      REAL zu(nlayer),zv(nlayer)
+      REAL ztsrf(ngrid)
+      REAL ztc(nlayer), ztm(nlayer+1) 
+      REAL zum(nlayer+1) , zvm(nlayer+1)
+      LOGICAL condsub(ngrid)
+      REAL subptimestep
+      Integer Ntime
+      real masse (ngrid,nlayer), w(ngrid,nlayer,nq)
+      real wq(ngrid,nlayer+1)
+      real vstokes,reff
+
+!     Special diagnostic variables
+      real tconda1(ngrid,nlayer)
+      real tconda2(ngrid,nlayer)
+      real zdtsig (ngrid,nlayer)
+      real zdt (ngrid,nlayer)
+
+!-----------------------------------------------------------------------
+!     Saved local variables
+
+      REAL,SAVE :: latcond=5.9e5
+      REAL,SAVE :: ccond
+      REAL,SAVE :: cpice=1000.
+      REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: emisref
+!$OMP THREADPRIVATE(latcond,ccond,cpice,emisref)
+
+      LOGICAL,SAVE :: firstcall=.true.
+!$OMP THREADPRIVATE(firstcall)
+      REAL,EXTERNAL :: SSUM
+
+      REAL,EXTERNAL :: CBRT
+
+      INTEGER,SAVE :: i_co2ice=0      ! co2 ice
+!$OMP THREADPRIVATE(i_co2ice)
+      CHARACTER(LEN=20) :: tracername ! to temporarily store text
+
+      integer igas
+
+      real ppco2
+
+!-----------------------------------------------------------------------
+!     Initializations
+
+      pdqc(1:ngrid,1:nlayer,1:nq)=0
+      pdtc(1:ngrid,1:nlayer)=0
+      zq(1:ngrid,1:nlayer,1:nq)=0
+      zt(1:ngrid,1:nlayer)=0
+
+!     Initialisation
+      IF (firstcall) THEN
+
+         ALLOCATE(emisref(ngrid)) !! this should be deallocated in lastcall...
+
+         ! find CO2 ice tracer
+         do iq=1,nq
+            tracername=noms(iq)
+            if (tracername.eq."co2_ice") then
+               i_co2ice=iq
+            endif
+         enddo
+         
+        write(*,*) "condense_co2: i_co2ice=",i_co2ice       
+
+        if((i_co2ice.lt.1))then
+           print*,'In condens_cloud but no CO2 ice tracer, exiting.'
+           print*,'Still need generalisation to arbitrary species!'
+           stop
+        endif
+
+         ccond=cpp/(g*latcond)
+         print*,'In condens_cloud: ccond=',ccond,' latcond=',latcond
+
+!          Prepare special treatment if gas is not pure CO2
+             !if (addn2) then
+             !   m_co2   = 44.01E-3 ! CO2 molecular mass (kg/mol)   
+             !   m_noco2 = 28.02E-3 ! N2 molecular mass (kg/mol)  
+!               Compute A and B coefficient use to compute
+!               mean molecular mass Mair defined by
+!               1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
+!               1/Mair = A*q(ico2) + B
+             !   A = (1/m_co2 - 1/m_noco2)
+             !   B = 1/m_noco2
+             !endif
+
+!          Minimum CO2 mixing ratio below which mixing occurs with layer above: 
+           !qco2min =0.75  
+
+           firstcall=.false.
+      ENDIF
+      zcpi=1./cpp
+
+!-----------------------------------------------------------------------
+!     Calculation of CO2 condensation / sublimation 
+!
+!     Variables used:
+!     piceco2(ngrid)       amount of co2 ice on the ground  (kg/m2)
+!     zcondicea(ngrid,l)   condensation rate in layer l     (kg/m2/s)
+!     zcondices(ngrid)     condensation rate on the ground  (kg/m2/s)
+!     zfallice(ngrid)      flux of ice falling on surface   (kg/m2/s)
+!     pdtc(ngrid,nlayer) dT/dt due to phase changes       (K/s)
+     
+
+!     Tendencies initially set to 0 (except pdtc)
+      DO l=1,nlayer
+         DO ig=1,ngrid
+            zcondicea(ig,l) = 0.
+            pduc(ig,l)  = 0
+            pdvc(ig,l)  = 0
+            pdqc(ig,l,i_co2ice)  = 0
+         END DO
+      END DO
+      
+      DO ig=1,ngrid
+         Mfallice(ig) = 0.
+         zfallice(ig) = 0.
+         zcondices(ig) = 0.
+         pdtsrfc(ig) = 0.
+         pdpsrf(ig) = 0.
+         condsub(ig) = .false.
+         zdiceco2(ig) = 0.
+      ENDDO
+
+!-----------------------------------------------------------------------
+!     Atmospheric condensation
+
+
+!     Compute CO2 Volume mixing ratio
+!     -------------------------------
+!      if (addn2) then
+!         DO l=1,nlayer
+!            DO ig=1,ngrid
+!              qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep
+!!             Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
+!              mmean=1/(A*qco2 +B)
+!              vmr_co2(ig,l) = qco2*mmean/m_co2 
+!            ENDDO
+!         ENDDO
+!      else
+!         DO l=1,nlayer
+!            DO ig=1,ngrid
+!              vmr_co2(ig,l)=0.5
+!            ENDDO
+!         ENDDO
+!      end if
+
+!     Forecast the atmospheric frost temperature 'ztcond' and nucleation temperature 'ztnuc'
+      DO l=1,nlayer
+         DO ig=1,ngrid
+            ppco2=gfrac(igas_CO2)*pplay(ig,l)
+            call get_tcond_co2(ppco2,ztcond(ig,l))
+            call get_tnuc_co2(ppco2,ztnuc(ig,l))
+         ENDDO
+      ENDDO
+      
+!     Initialize zq and zt at the beginning of the sub-timestep loop
+      DO l=1,nlayer
+         DO ig=1,ngrid
+            zt(ig,l)=pt(ig,l)
+            zq(ig,l,i_co2ice)=pq(ig,l,i_co2ice)
+            IF( zq(ig,l,i_co2ice).lt.-1.e-6 ) THEN
+               print*,'Uh-oh, zq = ',zq(ig,l,i_co2ice),'at ig,l=',ig,l
+               if(l.eq.1)then
+                  print*,'Perhaps the atmosphere is collapsing on surface...?'
+               endif
+            END IF
+         ENDDO
+      ENDDO
+
+!     Calculate the mass of each atmospheric layer (kg.m-2)
+      do  ilay=1,nlayer
+         DO ig=1,ngrid
+            masse(ig,ilay)=(pplev(ig,ilay) - pplev(ig,ilay+1)) /g
+         end do
+      end do
+
+!     -----------------------------------------------
+!     START CONDENSATION/SEDIMENTATION SUB-TIME LOOP
+!     -----------------------------------------------
+      Ntime =  20               ! number of sub-timestep 
+      subptimestep = ptimestep/float(Ntime)           
+
+      DO it=1,Ntime
+
+!     Add the tendencies from other physical processes at each subtimstep
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               zt(ig,l)   = zt(ig,l)   + pdt(ig,l)   * subptimestep
+               zq(ig,l,i_co2ice) = zq(ig,l,i_co2ice) + pdq(ig,l,i_co2ice) * subptimestep
+            END DO
+         END DO
+
+
+!     Gravitational sedimentation
+            
+!     sedimentation computed from radius computed from q in module radii_mod
+	 call co2_reffrad(ngrid,nlayer,nq,zq,reffrad)
+	 
+         do  ilay=1,nlayer
+            DO ig=1,ngrid
+
+               reff = reffrad(ig,ilay)
+
+               call stokes                      &
+                   (pplev(ig,ilay),pt(ig,ilay), &
+                    reff,vstokes,rho_co2)
+
+               !w(ig,ilay,i_co2ice) = 0.0
+               w(ig,ilay,i_co2ice) = vstokes *  subptimestep * &
+                   pplev(ig,ilay)/(r*pt(ig,ilay))
+
+            end do
+         end do
+
+!     Computing q after sedimentation
+
+         call vlz_fi(ngrid,nlayer,zq(1,1,i_co2ice),2.,masse,w(1,1,i_co2ice),wq)
+
+
+!     Progressively accumulating the flux to the ground
+!     Mfallice is the total amount of ice fallen to the ground
+         DO ig=1,ngrid
+            Mfallice(ig) =  Mfallice(ig) + wq(ig,i_co2ice)
+         end do
+
+
+!     Condensation / sublimation in the atmosphere
+!     --------------------------------------------
+!     (calculation of zcondicea, zfallice and pdtc)
+!     (MODIFICATIONS FOR EARLY MARS: falling heat neglected, condensation
+!     of CO2 into tracer i_co2ice)
+
+         DO l=nlayer , 1, -1
+            DO ig=1,ngrid
+               pdtc(ig,l)=0.
+
+
+   !     ztcond-> ztnuc in test beneath to nucleate only when super saturation occurs(JL 2011)
+               IF ((zt(ig,l).LT.ztnuc(ig,l)).or.(zq(ig,l,i_co2ice).gt.1.E-10)) THEN 
+                  condsub(ig)=.true.
+                  pdtc(ig,l)   = (ztcond(ig,l) - zt(ig,l))/subptimestep
+                  pdqc(ig,l,i_co2ice) = pdtc(ig,l)*ccond*g
+
+!     Case when the ice from above sublimes entirely
+                  IF ((zq(ig,l,i_co2ice).lt.-pdqc(ig,l,i_co2ice)*subptimestep) &
+                      .AND. (zq(ig,l,i_co2ice).gt.0)) THEN
+
+                     pdqc(ig,l,i_co2ice) = -zq(ig,l,i_co2ice)/subptimestep
+                     pdtc(ig,l)   =-zq(ig,l,i_co2ice)/(ccond*g*subptimestep)
+
+                  END IF
+
+!     Temperature and q after condensation
+                  zt(ig,l)   = zt(ig,l)   + pdtc(ig,l)   * subptimestep
+                  zq(ig,l,i_co2ice) = zq(ig,l,i_co2ice) + pdqc(ig,l,i_co2ice) * subptimestep
+               END IF
+
+            ENDDO
+         ENDDO
+      ENDDO                     ! end of subtimestep loop
+
+!     Computing global tendencies after the subtimestep
+      DO l=1,nlayer
+         DO ig=1,ngrid
+            pdtc(ig,l) = &
+              (zt(ig,l) - (pt(ig,l) + pdt(ig,l)*ptimestep))/ptimestep
+            pdqc(ig,l,i_co2ice) = &
+              (zq(ig,l,i_co2ice)-(pq(ig,l,i_co2ice)+pdq(ig,l,i_co2ice)*ptimestep))/ptimestep
+         END DO
+      END DO
+      DO ig=1,ngrid
+         zfallice(ig) = Mfallice(ig)/ptimestep
+      END DO
+
+
+!-----------------------------------------------------------------------
+!     Condensation/sublimation on the ground
+!     (calculation of zcondices and pdtsrfc)
+
+!     forecast of ground temperature ztsrf and frost temperature ztcondsol
+      DO ig=1,ngrid
+         ppco2=gfrac(igas_CO2)*pplay(ig,1)
+         call get_tcond_co2(ppco2,ztcondsol(ig))
+         
+         ztsrf(ig) = ptsrf(ig)
+
+         if((ztsrf(ig).le.ztcondsol(ig)+2.0).and.(ngrid.eq.1))then
+            print*,'CO2 is condensing on the surface in 1D. This atmosphere is doomed.'
+            print*,'T_surf = ',ztsrf,'K'
+            print*,'T_cond = ',ztcondsol,'K'
+            open(116,file='surf_vals.out')
+            write(116,*) 0.0, pplev(1,1), 0.0, 0.0
+            close(116)
+            call abort
+         endif
+
+         ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
+
+      ENDDO
+     
+      DO ig=1,ngrid
+         IF(ig.GT.ngrid/2+1) THEN
+            icap=2
+         ELSE
+            icap=1
+         ENDIF
+
+!     Loop over where we have condensation / sublimation
+         IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR.           &        ! ground condensation
+                    (zfallice(ig).NE.0.) .OR.              &        ! falling snow
+                    ((ztsrf(ig) .GT. ztcondsol(ig)) .AND.  &        ! ground sublimation
+                    ((piceco2(ig)+zfallice(ig)*ptimestep) .NE. 0.))) THEN
+            condsub(ig) = .true.
+
+!     Condensation or partial sublimation of CO2 ice
+            zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
+                /(latcond*ptimestep)         
+            pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
+
+!     If the entire CO_2 ice layer sublimes
+!     (including what has just condensed in the atmosphere)
+            IF((piceco2(ig)/ptimestep+zfallice(ig)).LE. &
+                -zcondices(ig))THEN
+               zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig)
+               pdtsrfc(ig)=(latcond/pcapcal(ig))*       &
+                   (zcondices(ig))
+            END IF
+
+!     Changing CO2 ice amount and pressure
+
+            zdiceco2(ig) = zcondices(ig) + zfallice(ig)
+            piceco2(ig)  = piceco2(ig) + zdiceco2(ig)*ptimestep
+            pdpsrf(ig)   = -zdiceco2(ig)*g
+
+            IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
+               PRINT*,'STOP in condens'
+               PRINT*,'condensing more than total mass'
+               PRINT*,'Grid point ',ig
+               PRINT*,'Ps = ',pplev(ig,1)
+               PRINT*,'d Ps = ',pdpsrf(ig)
+               STOP
+            ENDIF
+         END IF
+      ENDDO
+
+!     Surface albedo and emissivity of the ground below the snow (emisref)
+!     --------------------------------------------------------------------
+      DO ig=1,ngrid
+         IF(lati(ig).LT.0.) THEN
+            icap=2 ! Southern Hemisphere
+         ELSE
+            icap=1 ! Nortnern hemisphere
+         ENDIF
+
+         if(.not.piceco2(ig).ge.0.) THEN
+            if(piceco2(ig).le.-1.e-8) print*,   &
+                'WARNING in condense_co2cloud: piceco2(',ig,')=', piceco2(ig)
+            piceco2(ig)=0.
+         endif
+         if (piceco2(ig).gt.0) then
+            psolaralb(ig) = albedice(icap)
+            emisref(ig)   = emisice(icap)
+         else
+            psolaralb(ig) = albedodat(ig)
+            emisref(ig)   = emissiv
+            pemisurf(ig)  = emissiv
+         end if
+      end do
+
+      return
+    end subroutine condense_co2
+
+!-------------------------------------------------------------------------
+    subroutine get_tcond_co2(p,tcond)
+!   Calculates the condensation temperature for CO2
+
+      implicit none
+
+      real p, peff, tcond
+      real, parameter :: ptriple=518000.0
+
+      peff=p
+
+      if(peff.lt.ptriple)then
+         tcond = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula
+      else
+         tcond = 684.2-92.3*log(peff)+4.32*log(peff)**2 
+         ! liquid-vapour transition (based on CRC handbook 2003 data)
+      endif
+      return
+
+    end subroutine get_tcond_co2
+    
+    
+    
+ 
+!-------------------------------------------------------------------------
+    subroutine get_tnuc_co2(p,tnuc)
+!   Calculates the nucleation temperature for CO2, based on a simple super saturation criterion 
+!     (JL 2011)
+
+      use callkeys_mod, only: co2supsat
+
+      implicit none
+
+      real p, peff, tnuc
+      real, parameter :: ptriple=518000.0
+
+      peff=p/co2supsat
+
+      if(peff.lt.ptriple)then
+         tnuc = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula
+      else
+         tnuc = 684.2-92.3*log(peff)+4.32*log(peff)**2 
+         ! liquid-vapour transition (based on CRC handbook 2003 data)
+      endif
+      return
+
+    end subroutine get_tnuc_co2
Index: trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 1470)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 1477)
@@ -52,54 +52,67 @@
 !   It includes:
 !
-!      1.  Initialization:
-!      1.1 Firstcall initializations
-!      1.2 Initialization for every call to physiq
-!      1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
-!      2. Compute radiative transfer tendencies
-!         (longwave and shortwave).
-!      4. Vertical diffusion (turbulent mixing):
-!      5. Convective adjustment
-!      6. Condensation and sublimation of gases (currently just CO2).
-!      7. TRACERS
-!       7a. water and water ice
-!       7c. other schemes for tracer transport (lifting, sedimentation)
-!       7d. updates (pressure variations, surface budget)
-!      9. Surface and sub-surface temperature calculations
-!     10. Write outputs :
-!           - "startfi", "histfi" if it's time
-!           - Saving statistics if "callstats = .true."
-!           - Output any needed variables in "diagfi" 
-!     10. Diagnostics: mass conservation of tracers, radiative energy balance etc.
+!      I. Initialization :
+!         I.1 Firstcall initializations.
+!         I.2 Initialization for every call to physiq.
+!
+!      II. Compute radiative transfer tendencies (longwave and shortwave) :
+!         II.a Option 1 : Call correlated-k radiative transfer scheme.
+!         II.b Option 2 : Call Newtonian cooling scheme.
+!         II.c Option 3 : Atmosphere has no radiative effect.
+!
+!      III. Vertical diffusion (turbulent mixing) :
+!
+!      IV. Dry Convective adjustment :
+!
+!      V. Condensation and sublimation of gases (currently just CO2).
+!
+!      VI. Tracers
+!         VI.1. Water and water ice.
+!         VI.2. Aerosols and particles.
+!         VI.3. Updates (pressure variations, surface budget).
+!         VI.4. Slab Ocean.
+!         VI.5. Surface Tracer Update.
+!
+!      VII. Surface and sub-surface soil temperature.
+!
+!      VIII. Perform diagnostics and write output files.
+!
 !
 !   arguments
 !   ---------
 !
-!   input
+!   INPUT
 !   -----
-!    ecri                  period (in dynamical timestep) to write output
+!
 !    ngrid                 Size of the horizontal grid.
-!                          All internal loops are performed on that grid.
 !    nlayer                Number of vertical layers.
-!    nq                    Number of advected fields
-!    firstcall             True at the first call
-!    lastcall              True at the last call
-!    pday                  Number of days counted from the North. Spring
-!                          equinoxe.
-!    ptime                 Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
-!    ptimestep             timestep (s)
-!    pplay(ngrid,nlayer)   Pressure at the middle of the layers (Pa)
-!    pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
-!    pphi(ngrid,nlayer)    Geopotential at the middle of the layers (m2s-2)
-!    pu(ngrid,nlayer)      u component of the wind (ms-1)
-!    pv(ngrid,nlayer)      v component of the wind (ms-1)
-!    pt(ngrid,nlayer)      Temperature (K)
-!    pq(ngrid,nlayer,nq)   Advected fields
+!    nq                    Number of advected fields.
+!    nametrac              Name of corresponding advected fields.
+!
+!    firstcall             True at the first call.
+!    lastcall              True at the last call.
+!
+!    pday                  Number of days counted from the North. Spring equinoxe.
+!    ptime                 Universal time (0<ptime<1): ptime=0.5 at 12:00 UT.
+!    ptimestep             timestep (s).
+!
+!    pplay(ngrid,nlayer)   Pressure at the middle of the layers (Pa).
+!    pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa).
+!    pphi(ngrid,nlayer)    Geopotential at the middle of the layers (m2.s-2).
+!
+!    pu(ngrid,nlayer)      u, zonal component of the wind (ms-1).
+!    pv(ngrid,nlayer)      v, meridional component of the wind (ms-1).
+!
+!    pt(ngrid,nlayer)      Temperature (K).
+!
+!    pq(ngrid,nlayer,nq)   Advected fields.
+!
 !    pudyn(ngrid,nlayer)    \
 !    pvdyn(ngrid,nlayer)     \ Dynamical temporal derivative for the
-!    ptdyn(ngrid,nlayer)     / corresponding variables
+!    ptdyn(ngrid,nlayer)     / corresponding variables.
 !    pqdyn(ngrid,nlayer,nq) /
 !    flxw(ngrid,nlayer)      vertical mass flux (kg/s) at layer lower boundary
 !
-!   output
+!   OUTPUT
 !   ------
 !
@@ -127,4 +140,5 @@
 !           Loops converted to F90 matrix format: J. Leconte (2012)
 !           No more ngridmx/nqmx, F90 commons and adaptation to parallel: A. Spiga (2012)
+!           Purge of the code : M. Turbet (2015)
 !==================================================================
 
@@ -138,66 +152,68 @@
 ! -----------
 
-!   inputs:
+!   INPUTS:
 !   -------
-      integer,intent(in) :: ngrid ! number of atmospheric columns
-      integer,intent(in) :: nlayer ! number of atmospheric layers
-      integer,intent(in) :: nq ! number of tracers
-      character*20,intent(in) :: nametrac(nq) ! name of the tracer from dynamics
-      logical,intent(in) :: firstcall ! signals first call to physics
-      logical,intent(in) :: lastcall ! signals last call to physics
-      real,intent(in) :: pday ! number of elapsed sols since reference Ls=0
-      real,intent(in) :: ptime ! "universal time", given as fraction of sol (e.g.: 0.5 for noon)
-      real,intent(in) :: ptimestep ! physics timestep (s)
-      real,intent(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
-      real,intent(in) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
-      real,intent(in) :: pphi(ngrid,nlayer) ! geopotential at mid-layer (m2s-2)
-      real,intent(in) :: pu(ngrid,nlayer) ! zonal wind component (m/s)
-      real,intent(in) :: pv(ngrid,nlayer) ! meridional wind component (m/s)
-      real,intent(in) :: pt(ngrid,nlayer) ! temperature (K)
-      real,intent(in) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
-      real,intent(in) :: flxw(ngrid,nlayer) ! vertical mass flux (ks/s)
-                                            ! at lower boundary of layer
-
-
-!   outputs:
+
+      integer,intent(in) :: ngrid             ! Number of atmospheric columns.
+      integer,intent(in) :: nlayer            ! Number of atmospheric layers.
+      integer,intent(in) :: nq                ! Number of tracers.
+      character*20,intent(in) :: nametrac(nq) ! Names of the tracers taken from dynamics.
+      
+      logical,intent(in) :: firstcall ! Signals first call to physics.
+      logical,intent(in) :: lastcall  ! Signals last call to physics.
+      
+      real,intent(in) :: pday                  ! Number of elapsed sols since reference Ls=0.
+      real,intent(in) :: ptime                 ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
+      real,intent(in) :: ptimestep             ! Physics timestep (s).
+      real,intent(in) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
+      real,intent(in) :: pplay(ngrid,nlayer)   ! Mid-layer pressure (Pa).
+      real,intent(in) :: pphi(ngrid,nlayer)    ! Geopotential at mid-layer (m2s-2).
+      real,intent(in) :: pu(ngrid,nlayer)      ! Zonal wind component (m/s).
+      real,intent(in) :: pv(ngrid,nlayer)      ! Meridional wind component (m/s).
+      real,intent(in) :: pt(ngrid,nlayer)      ! Temperature (K).
+      real,intent(in) :: pq(ngrid,nlayer,nq)   ! Tracers (kg/kg_of_air).
+      real,intent(in) :: flxw(ngrid,nlayer)    ! Vertical mass flux (ks/s) at lower boundary of layer
+
+!   OUTPUTS:
 !   --------
-!     physical tendencies
-      real,intent(out) :: pdu(ngrid,nlayer) ! zonal wind tendency (m/s/s)
-      real,intent(out) :: pdv(ngrid,nlayer) ! meridional wind tendency (m/s/s)
-      real,intent(out) :: pdt(ngrid,nlayer) ! temperature tendency (K/s)
-      real,intent(out) :: pdq(ngrid,nlayer,nq) ! tracer tendencies (../kg/s)
-      real,intent(out) :: pdpsrf(ngrid) ! surface pressure tendency (Pa/s)
-      logical,intent(out) :: tracerdyn ! signal to the dynamics to advect tracers or not
+
+!     Physical tendencies :
+
+      real,intent(out) :: pdu(ngrid,nlayer)    ! Zonal wind tendencies (m/s/s).
+      real,intent(out) :: pdv(ngrid,nlayer)    ! Meridional wind tendencies (m/s/s).
+      real,intent(out) :: pdt(ngrid,nlayer)    ! Temperature tendencies (K/s).
+      real,intent(out) :: pdq(ngrid,nlayer,nq) ! Tracer tendencies (kg/kg_of_air/s).
+      real,intent(out) :: pdpsrf(ngrid)        ! Surface pressure tendency (Pa/s).
+      logical,intent(out) :: tracerdyn         ! Signal to the dynamics to advect tracers or not.     
 
 ! Local saved variables:
 ! ----------------------
-!     aerosol (dust or ice) extinction optical depth  at reference wavelength 
-!     "longrefvis" set in dimradmars.h , for one of the "naerkind"  kind of
-!      aerosol optical properties:
-!      real aerosol(ngrid,nlayer,naerkind) 
-!     this is now internal to callcorrk and hence no longer needed here
-
-      integer day_ini                ! Initial date of the run (sol since Ls=0) 
-      integer icount                 ! counter of calls to physiq during the run.
-      real, dimension(:),allocatable,save ::  tsurf  ! Surface temperature (K)
-      real, dimension(:,:),allocatable,save ::  tsoil  ! sub-surface temperatures (K)
-      real, dimension(:),allocatable,save :: albedo ! Surface albedo
+
+      integer day_ini                                  ! Initial date of the run (sol since Ls=0).
+      integer icount                                   ! Counter of calls to physiq during the run.
+      real, dimension(:),allocatable,save ::  tsurf    ! Surface temperature (K).
+      real, dimension(:,:),allocatable,save ::  tsoil  ! Sub-surface temperatures (K).
+      real, dimension(:),allocatable,save :: albedo    ! Surface albedo.
+      
 !$OMP THREADPRIVATE(tsurf,tsoil,albedo)
 
-      real,dimension(:),allocatable,save :: albedo0 ! Surface albedo
-      real,dimension(:),allocatable,save :: rnat ! added by BC
+      real,dimension(:),allocatable,save :: albedo0 ! Initial surface albedo.
+      real,dimension(:),allocatable,save :: rnat    ! Defines the type of the grid (ocean,continent,...). By BC.
+      
 !$OMP THREADPRIVATE(albedo0,rnat)
 
-      real,dimension(:),allocatable,save :: emis ! Thermal IR surface emissivity
-      real,dimension(:,:),allocatable,save :: dtrad ! Net atm. radiative heating rate (K.s-1)
-      real,dimension(:),allocatable,save :: fluxrad_sky ! rad. flux from sky absorbed by surface (W.m-2)
-      real,dimension(:),allocatable,save :: fluxrad ! Net radiative surface flux (W.m-2)
-      real,dimension(:),allocatable,save :: capcal ! surface heat capacity (J m-2 K-1)
-      real,dimension(:),allocatable,save :: fluxgrd ! surface conduction flux (W.m-2)
-      real,dimension(:,:),allocatable,save :: qsurf ! tracer on surface (e.g. kg.m-2)
-      real,dimension(:,:),allocatable,save :: q2 ! Turbulent Kinetic Energy 
+      real,dimension(:),allocatable,save :: emis        ! Thermal IR surface emissivity.
+      real,dimension(:,:),allocatable,save :: dtrad     ! Net atmospheric radiative heating rate (K.s-1).
+      real,dimension(:),allocatable,save :: fluxrad_sky ! Radiative flux from sky absorbed by surface (W.m-2).
+      real,dimension(:),allocatable,save :: fluxrad     ! Net radiative surface flux (W.m-2).
+      real,dimension(:),allocatable,save :: capcal      ! Surface heat capacity (J m-2 K-1).
+      real,dimension(:),allocatable,save :: fluxgrd     ! Surface conduction flux (W.m-2).
+      real,dimension(:,:),allocatable,save :: qsurf     ! Tracer on surface (e.g. kg.m-2).
+      real,dimension(:,:),allocatable,save :: q2        ! Turbulent Kinetic Energy.
+      
 !$OMP THREADPRIVATE(emis,dtrad,fluxrad_sky,fluxrad,capcal,fluxgrd,qsurf,q2)
 
       save day_ini, icount
+      
 !$OMP THREADPRIVATE(day_ini,icount)
 
@@ -205,80 +221,96 @@
 ! -----------------
 
-!     aerosol (dust or ice) extinction optical depth  at reference wavelength 
+!     Aerosol (dust or ice) extinction optical depth  at reference wavelength 
 !     for the "naerkind" optically active aerosols:
-      real aerosol(ngrid,nlayer,naerkind) 
-      real zh(ngrid,nlayer)      ! potential temperature (K)
-      real pw(ngrid,nlayer) ! vertical velocity (m/s) (>0 when downwards)
-
-      character*80 fichier 
-      integer l,ig,ierr,iq,iaer
-      
-      !!! this is saved for diagnostic
-      real,dimension(:),allocatable,save :: fluxsurf_lw ! incident LW (IR) surface flux (W.m-2)
-      real,dimension(:),allocatable,save :: fluxsurf_sw ! incident SW (stellar) surface flux (W.m-2)
-      real,dimension(:),allocatable,save :: fluxtop_lw ! Outgoing LW (IR) flux to space (W.m-2)
-      real,dimension(:),allocatable,save :: fluxabs_sw ! Absorbed SW (stellar) flux (W.m-2)
-      real,dimension(:),allocatable,save :: fluxtop_dn
-      real,dimension(:),allocatable,save :: fluxdyn ! horizontal heat transport by dynamics  
-      real,dimension(:,:),allocatable,save :: OLR_nu ! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
-      real,dimension(:,:),allocatable,save :: OSR_nu ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
-      real,dimension(:,:),allocatable,save :: zdtlw ! (K/s)
-      real,dimension(:,:),allocatable,save :: zdtsw ! (K/s)
-      real,dimension(:),allocatable,save :: sensibFlux ! turbulent flux given by the atm to the surface
+
+      real aerosol(ngrid,nlayer,naerkind) ! Aerosols.
+      real zh(ngrid,nlayer)               ! Potential temperature (K).
+      real pw(ngrid,nlayer)               ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
+
+      integer l,ig,ierr,iq
+      
+      ! FOR DIAGNOSTIC :
+      
+      real,dimension(:),allocatable,save :: fluxsurf_lw  ! Incident Long Wave (IR) surface flux (W.m-2).
+      real,dimension(:),allocatable,save :: fluxsurf_sw  ! incident Short Wave (stellar) surface flux (W.m-2).
+      real,dimension(:),allocatable,save :: fluxtop_lw   ! Outgoing LW (IR) flux to space (W.m-2).
+      real,dimension(:),allocatable,save :: fluxabs_sw   ! Absorbed SW (stellar) flux (W.m-2).
+      real,dimension(:),allocatable,save :: fluxtop_dn   ! Incoming SW (stellar) radiation at the top of the atmosphere (W.m-2).
+      real,dimension(:),allocatable,save :: fluxdyn      ! Horizontal heat transport by dynamics (W.m-2).
+      real,dimension(:,:),allocatable,save :: OLR_nu     ! Outgoing LW radiation in each band (Normalized to the band width (W/m2/cm-1)).
+      real,dimension(:,:),allocatable,save :: OSR_nu     ! Outgoing SW radiation in each band (Normalized to the band width (W/m2/cm-1)).
+      real,dimension(:,:),allocatable,save :: zdtlw      ! LW heating tendencies (K/s).
+      real,dimension(:,:),allocatable,save :: zdtsw      ! SW heating tendencies (K/s).
+      real,dimension(:),allocatable,save :: sensibFlux   ! Turbulent flux given by the atmosphere to the surface (W.m-2).
+      
 !$OMP THREADPRIVATE(fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu,&
-	!$OMP zdtlw,zdtsw,sensibFlux)
-
-      real zls                       ! solar longitude (rad)
-      real zlss                      ! sub solar point longitude (rad)
-      real zday                      ! date (time since Ls=0, in martian days)
-      real zzlay(ngrid,nlayer)   ! altitude at the middle of the layers
-      real zzlev(ngrid,nlayer+1) ! altitude at layer boundaries
-      real latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
-
-!     Tendencies due to various processes:
-      real dqsurf(ngrid,nq)
-      real cldtlw(ngrid,nlayer)                           ! (K/s) LW heating rate for clear areas
-      real cldtsw(ngrid,nlayer)                           ! (K/s) SW heating rate for clear areas
-      real zdtsurf(ngrid)                                   ! (K/s)
-      real dtlscale(ngrid,nlayer)
-      real zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer)  ! (m.s-2)
-      real zdhdif(ngrid,nlayer), zdtsdif(ngrid)         ! (K/s)
-      real zdtdif(ngrid,nlayer)			      ! (K/s)
-      real zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer)  ! (m.s-2)
-      real zdhadj(ngrid,nlayer)                           ! (K/s)
-      real zdtgw(ngrid,nlayer)                            ! (K/s)
-      real zdtmr(ngrid,nlayer)
-      real zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer)    ! (m.s-2)
-      real zdtc(ngrid,nlayer),zdtsurfc(ngrid)
-      real zdvc(ngrid,nlayer),zduc(ngrid,nlayer)
-      real zdumr(ngrid,nlayer),zdvmr(ngrid,nlayer)
-      real zdtsurfmr(ngrid)
-      
-      real zdmassmr(ngrid,nlayer),zdpsrfmr(ngrid)
-      real zdmassmr_col(ngrid)
-
-      real zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq)
-      real zdqevap(ngrid,nlayer)
-      real zdqsed(ngrid,nlayer,nq), zdqssed(ngrid,nq)
-      real zdqdev(ngrid,nlayer,nq), zdqsdev(ngrid,nq)
-      real zdqadj(ngrid,nlayer,nq)
-      real zdqc(ngrid,nlayer,nq)
-      real zdqmr(ngrid,nlayer,nq),zdqsurfmr(ngrid,nq)
-      real zdqlscale(ngrid,nlayer,nq)
-      real zdqslscale(ngrid,nq)
-      real zdqchim(ngrid,nlayer,nq)
-      real zdqschim(ngrid,nq)
-
-      real zdteuv(ngrid,nlayer)    ! (K/s)
-      real zdtconduc(ngrid,nlayer) ! (K/s)
-      real zdumolvis(ngrid,nlayer)
-      real zdvmolvis(ngrid,nlayer)
-      real zdqmoldiff(ngrid,nlayer,nq)
-
-!     Local variables for local calculations:
+	
+        !$OMP zdtlw,zdtsw,sensibFlux)
+
+      real zls                       ! Solar longitude (radians).
+      real zlss                      ! Sub solar point longitude (radians).
+      real zday                      ! Date (time since Ls=0, calculated in sols).
+      real zzlay(ngrid,nlayer)       ! Altitude at the middle of the atmospheric layers.
+      real zzlev(ngrid,nlayer+1)     ! Altitude at the atmospheric layer boundaries.
+
+! TENDENCIES due to various processes :
+
+      ! For Surface Temperature : (K/s)
+      real zdtsurf(ngrid)     ! Cumulated tendencies.
+      real zdtsurfmr(ngrid)   ! Mass_redistribution routine.
+      real zdtsurfc(ngrid)    ! Condense_co2 routine.
+      real zdtsdif(ngrid)     ! Turbdiff/vdifc routines.
+      real zdtsurf_hyd(ngrid) ! Hydrol routine.
+            
+      ! For Atmospheric Temperatures : (K/s)    
+      real dtlscale(ngrid,nlayer)                             ! Largescale routine.
+      real zdtc(ngrid,nlayer)                                 ! Condense_co2 routine.
+      real zdtdif(ngrid,nlayer)	                              ! Turbdiff/vdifc routines.
+      real zdtmr(ngrid,nlayer)                                ! Mass_redistribution routine.
+      real zdtrain(ngrid,nlayer)                              ! Rain routine.
+      real dtmoist(ngrid,nlayer)                              ! Moistadj routine.
+      real dt_ekman(ngrid,noceanmx), dt_hdiff(ngrid,noceanmx) ! Slab_ocean routine.
+      real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer)         ! Callcorrk routine.
+                              
+      ! For Surface Tracers : (kg/m2/s)
+      real dqsurf(ngrid,nq)                 ! Cumulated tendencies.
+      real zdqslscale(ngrid,nq)             ! Largescale routine.
+      real zdqsdif(ngrid,nq)                ! Turbdiff/vdifc routines.
+      real zdqssed(ngrid,nq)                ! Callsedim routine.
+      real zdqsurfmr(ngrid,nq)              ! Mass_redistribution routine.
+      real zdqsrain(ngrid), zdqssnow(ngrid) ! Rain routine.
+      real dqs_hyd(ngrid,nq)                ! Hydrol routine.
+                  
+      ! For Tracers : (kg/kg_of_air/s)
+      real zdqc(ngrid,nlayer,nq)      ! Condense_co2 routine.
+      real zdqadj(ngrid,nlayer,nq)    ! Convadj routine.
+      real zdqdif(ngrid,nlayer,nq)    ! Turbdiff/vdifc routines.
+      real zdqevap(ngrid,nlayer)      ! Turbdiff routine.
+      real zdqsed(ngrid,nlayer,nq)    ! Callsedim routine.
+      real zdqmr(ngrid,nlayer,nq)     ! Mass_redistribution routine.
+      real zdqrain(ngrid,nlayer,nq)   ! Rain routine.
+      real dqmoist(ngrid,nlayer,nq)   ! Moistadj routine.
+      real dqvaplscale(ngrid,nlayer)  ! Largescale routine.
+      real dqcldlscale(ngrid,nlayer)  ! Largescale routine.
+                        
+      ! For Winds : (m/s/s)
+      real zdvc(ngrid,nlayer),zduc(ngrid,nlayer)     ! Condense_co2 routine.
+      real zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer) ! Convadj routine.
+      real zdumr(ngrid,nlayer),zdvmr(ngrid,nlayer)   ! Mass_redistribution routine.
+      real zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer) ! Turbdiff/vdifc routines.
+      real zdhdif(ngrid,nlayer)                      ! Turbdiff/vdifc routines.
+      real zdhadj(ngrid,nlayer)                      ! Convadj routine.
+     
+      ! For Pressure and Mass :
+      real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
+      real zdmassmr_col(ngrid)    ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s).
+      real zdpsrfmr(ngrid)        ! Pressure tendency for mass_redistribution routine (Pa/s).
+
+      
+    
+! Local variables for LOCAL CALCULATIONS:
+! ---------------------------------------
       real zflubid(ngrid)
       real zplanck(ngrid),zpopsk(ngrid,nlayer)
-      real zdum1(ngrid,nlayer)
-      real zdum2(ngrid,nlayer)
       real ztim1,ztim2,ztim3, z1,z2
       real ztime_fin
@@ -287,59 +319,29 @@
       real taux(ngrid),tauy(ngrid)
 
-      integer length
-      parameter (length=100)
-
-! local variables only used for diagnostics (output in file "diagfi" or "stats")
-! ------------------------------------------------------------------------------
-      real ps(ngrid), zt(ngrid,nlayer)
-      real zu(ngrid,nlayer),zv(ngrid,nlayer)
-      real zq(ngrid,nlayer,nq)
-      character*2 str2
-      character*5 str5
-      real zdtadj(ngrid,nlayer)
-      real zdtdyn(ngrid,nlayer)
-      real,allocatable,dimension(:,:),save :: ztprevious
+
+! local variables for DIAGNOSTICS : (diagfi & stat)
+! -------------------------------------------------
+      real ps(ngrid)                                     ! Surface Pressure.
+      real zt(ngrid,nlayer)                              ! Atmospheric Temperature.
+      real zu(ngrid,nlayer),zv(ngrid,nlayer)             ! Zonal and Meridional Winds.
+      real zq(ngrid,nlayer,nq)                           ! Atmospheric Tracers.
+      real zdtadj(ngrid,nlayer)                          ! Convadj Diagnostic.
+      real zdtdyn(ngrid,nlayer)                          ! Dynamical Heating (K/s).
+      real,allocatable,dimension(:,:),save :: ztprevious ! Previous loop Atmospheric Temperature (K)                                                         ! Useful for Dynamical Heating calculation.
 !$OMP THREADPRIVATE(ztprevious)
-      real reff(ngrid,nlayer) ! effective dust radius (used if doubleq=T)
-      real qtot1,qtot2            ! total aerosol mass
-      integer igmin, lmin
-      logical tdiag
-
-      real zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
-      real zstress(ngrid), cd
-      real hco2(nq), tmean, zlocal(nlayer)
-      real vmr(ngrid,nlayer) ! volume mixing ratio
-
+
+      real reff(ngrid,nlayer)                       ! Effective dust radius (used if doubleq=T).
+      real vmr(ngrid,nlayer)                        ! volume mixing ratio
       real time_phys
-
-!     reinstated by RW for diagnostic
-      real,allocatable,dimension(:),save :: tau_col
+      real,allocatable,dimension(:),save :: tau_col ! Total Aerosol Optical Depth.
 !$OMP THREADPRIVATE(tau_col)
       
-!     included by RW to reduce insanity of code
-      real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS
-
-!     included by RW to compute tracer column densities
-      real qcol(ngrid,nq)
-
-!     included by RW for H2O precipitation
-      real zdtrain(ngrid,nlayer)
-      real zdqrain(ngrid,nlayer,nq)
-      real zdqsrain(ngrid)
-      real zdqssnow(ngrid)
-      real rainout(ngrid)
+      real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
+      
+      real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2).
 
 !     included by RW for H2O Manabe scheme
-      real dtmoist(ngrid,nlayer)
-      real dqmoist(ngrid,nlayer,nq)
-
-      real qvap(ngrid,nlayer)
-      real dqvaplscale(ngrid,nlayer)
-      real dqcldlscale(ngrid,nlayer)
-      real rneb_man(ngrid,nlayer)
-      real rneb_lsc(ngrid,nlayer)
-
-!     included by RW to account for surface cooling by evaporation
-      real dtsurfh2olat(ngrid)
+      real rneb_man(ngrid,nlayer) ! H2O cloud fraction (moistadj).
+      real rneb_lsc(ngrid,nlayer) ! H2O cloud fraction (large scale).
 
 
@@ -352,94 +354,57 @@
       real dEdiffs(ngrid),dEdiff(ngrid)
       real madjdE(ngrid), lscaledE(ngrid),madjdEz(ngrid,nlayer), lscaledEz(ngrid,nlayer)
+      
 !JL12 conservation test for mean flow kinetic energy has been disabled temporarily
-      real dtmoist_max,dtmoist_min
-      
+
+      real dtmoist_max,dtmoist_min     
       real dItot, dItot_tmp, dVtot, dVtot_tmp
 
-!     included by BC for evaporation     
-      real qevap(ngrid,nlayer,nq)
-      real tevap(ngrid,nlayer)
-      real dqevap1(ngrid,nlayer)
-      real dtevap1(ngrid,nlayer)
-
-!     included by BC for hydrology
-      real,allocatable,save :: hice(:)
+      real,allocatable,save :: hice(:) ! Oceanic Ice height. by BC
  !$OMP THREADPRIVATE(hice)     
 
-!     included by RW to test water conservation (by routine)
-      real h2otot
+      real h2otot                      ! Total Amount of water. For diagnostic.
+      real icesrf,liqsrf,icecol,vapcol ! Total Amounts of water (ice,liq,vap). For diagnostic.
       real dWtot, dWtot_tmp, dWtots, dWtots_tmp
-      real h2o_surf_all
-      logical watertest
-      save watertest
+      logical,save :: watertest
 !$OMP THREADPRIVATE(watertest)
 
-!     included by RW for RH diagnostic
-      real qsat(ngrid,nlayer), RH(ngrid,nlayer), H2Omaxcol(ngrid),psat_tmp
-
-!     included by RW for hydrology
-      real dqs_hyd(ngrid,nq)
-      real zdtsurf_hyd(ngrid)
-
-!     included by RW for water cycle conservation tests
-      real icesrf,liqsrf,icecol,vapcol
-
-!     included by BC for double radiative transfer call
-      logical clearsky
-      real zdtsw1(ngrid,nlayer)
-      real zdtlw1(ngrid,nlayer)
-      real fluxsurf_lw1(ngrid)     
-      real fluxsurf_sw1(ngrid)  
-      real fluxtop_lw1(ngrid)    
-      real fluxabs_sw1(ngrid)  
-      real tau_col1(ngrid)
-      real OLR_nu1(ngrid,L_NSPECTI)
-      real OSR_nu1(ngrid,L_NSPECTV)
+      real qsat(ngrid,nlayer) ! Water Vapor Volume Mixing Ratio at saturation (kg/kg_of_air).
+      real RH(ngrid,nlayer)   ! Relative humidity.
+      real H2Omaxcol(ngrid)   ! Maximum possible H2O column amount (at 100% saturation) (kg/m2).
+      real psat_tmp
+      
+      logical clearsky ! For double radiative transfer call. By BC
+      
+      real fluxsurf_lw1(ngrid), fluxsurf_sw1(ngrid), fluxtop_lw1(ngrid), fluxabs_sw1(ngrid) ! For SW/LW flux diagnostics.
+      real tau_col1(ngrid)                                                                  ! For aerosol optical depth diagnostic.
+      real OLR_nu1(ngrid,L_NSPECTI), OSR_nu1(ngrid,L_NSPECTV)                               ! For Outgoing Radiation diagnostics.
       real tf, ntf
 
-!     included by BC for cloud fraction computations
-      real,allocatable,dimension(:,:),save :: cloudfrac
-      real,allocatable,dimension(:),save :: totcloudfrac
+      real,allocatable,dimension(:,:),save :: cloudfrac  ! Fraction of clouds (%).
+      real,allocatable,dimension(:),save :: totcloudfrac ! Column fraction of clouds (%).
 !$OMP THREADPRIVATE(cloudfrac,totcloudfrac)
 
-!     included by RW for vdifc water conservation test
-      real nconsMAX
-      real vdifcncons(ngrid)
-      real cadjncons(ngrid)
-
-!      double precision qsurf_hist(ngrid,nq)
+      real nconsMAX, vdifcncons(ngrid), cadjncons(ngrid) ! Vdfic water conservation test. By RW
+
       real,allocatable,dimension(:,:),save :: qsurf_hist
 !$OMP THREADPRIVATE(qsurf_hist)
 
-!     included by RW for temp convadj conservation test
-      real playtest(nlayer)
-      real plevtest(nlayer)
-      real ttest(nlayer)
-      real qtest(nlayer)
-      integer igtest
-
-!     included by RW for runway greenhouse 1D study
-      real muvar(ngrid,nlayer+1)
-
-!     included by RW for variable H2O particle sizes
-      real,allocatable,dimension(:,:,:),save :: reffrad ! aerosol effective radius (m)
-      real,allocatable,dimension(:,:,:),save :: nueffrad ! aerosol effective radius variance
+      real muvar(ngrid,nlayer+1) ! For Runaway Greenhouse 1D study. By RW
+
+      real,allocatable,dimension(:,:,:),save :: reffrad  ! Aerosol effective radius (m). By RW
+      real,allocatable,dimension(:,:,:),save :: nueffrad ! Aerosol effective radius variance. By RW
 !$OMP THREADPRIVATE(reffrad,nueffrad)
-!      real :: nueffrad_dummy(ngrid,nlayer,naerkind) !! AS. This is temporary. Check below why.
-      real :: reffh2oliq(ngrid,nlayer) ! liquid water particles effective radius (m)
-      real :: reffh2oice(ngrid,nlayer) ! water ice particles effective radius (m)
-   !   real reffH2O(ngrid,nlayer)
+
       real reffcol(ngrid,naerkind)
 
-!     included by RW for sourceevol
+!     Sourceevol for 'accelerated ice evolution'. By RW
       real,allocatable,dimension(:),save :: ice_initial
       real delta_ice,ice_tot
-      real,allocatable,dimension(:),save :: ice_min
-      
+      real,allocatable,dimension(:),save :: ice_min      
       integer num_run
       logical,save :: ice_update
 !$OMP THREADPRIVATE(ice_initial,ice_min,ice_update)
 
-!     included by BC for slab ocean
+!     For slab ocean. By BC
       real, dimension(:),allocatable,save ::  pctsrf_sic
       real, dimension(:,:),allocatable,save :: tslab
@@ -452,19 +417,20 @@
       real :: tsurf2(ngrid)
       real :: flux_o(ngrid),flux_g(ngrid),fluxgrdocean(ngrid)
-      real :: dt_ekman(ngrid,noceanmx),dt_hdiff(ngrid,noceanmx)
       real :: flux_sens_lat(ngrid)
       real :: qsurfint(ngrid,nq)
 
 
-!=======================================================================
-
-! 1. Initialisation
+!==================================================================================================
+
 ! -----------------
-
-!  1.1   Initialisation only at first call
-!  ---------------------------------------
+! I. INITIALISATION
+! -----------------
+
+! --------------------------------
+! I.1   First Call Initialisation.
+! --------------------------------
       if (firstcall) then
 
-        !!! ALLOCATE SAVED ARRAYS
+        ! Allocate saved arrays.
         ALLOCATE(tsurf(ngrid))
         ALLOCATE(tsoil(ngrid,nsoilmx))    
@@ -507,16 +473,13 @@
         ALLOCATE(zmasq(ngrid))
         ALLOCATE(knindex(ngrid))
-!        ALLOCATE(qsurfint(ngrid,nqmx))
-
-
-        !! this is defined in comsaison_h
+
+        ! This is defined in comsaison_h
         ALLOCATE(mu0(ngrid))
-        ALLOCATE(fract(ngrid))
-
-            
-         !! this is defined in radcommon_h
-         ALLOCATE(glat(ngrid))
-
-!        variables set to 0
+        ALLOCATE(fract(ngrid))            
+         ! This is defined in radcommon_h
+        ALLOCATE(glat(ngrid))
+        
+
+!        Variables set to 0
 !        ~~~~~~~~~~~~~~~~~~
          dtrad(:,:) = 0.0
@@ -526,26 +489,25 @@
          zdtlw(:,:) = 0.0
 
-!        initialize aerosol indexes 
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-            call iniaerosol()
-
-      
-!        initialize tracer names, indexes and properties
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+!        Initialize aerosol indexes.
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         call iniaerosol()
+
+      
+!        Initialize tracer names, indexes and properties.
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          tracerdyn=tracer
-         IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) !! because noms is an argument of physdem1
-                                                      !! whether or not tracer is on
+         IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
          if (tracer) then
             call initracer(ngrid,nq,nametrac)
-         endif ! end tracer
-
-! 
-
-!        read startfi (initial parameters)
+         endif
+
+
+!        Read 'startfi.nc' file.
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         call phyetat0(ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,   &
-               day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,   &
-               cloudfrac,totcloudfrac,hice,  &
-               rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
+         call phyetat0(ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
+                       day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,   &
+                       cloudfrac,totcloudfrac,hice,                   &
+                       rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
 
          if (pday.ne.day_ini) then
@@ -559,5 +521,5 @@
          write (*,*) 'In physiq day_ini =', day_ini
 
-!        Initialize albedo and orbital calculation
+!        Initialize albedo and orbital calculation.
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          call surfini(ngrid,nq,qsurf,albedo0)
@@ -571,15 +533,15 @@
          endif
 
-!        Initialize oceanic variables
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!        Initialize oceanic variables.
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
          if (ok_slab_ocean)then
 
            call ocean_slab_init(ngrid,ptimestep, tslab, &
-                sea_ice, pctsrf_sic)
-
-            knindex(:) = 0
-
-            do ig=1,ngrid
+                                sea_ice, pctsrf_sic)
+
+           knindex(:) = 0
+
+           do ig=1,ngrid
               zmasq(ig)=1
               knindex(ig) = ig
@@ -587,46 +549,45 @@
                  zmasq(ig)=0
               endif
-            enddo
-
+           enddo
 
            CALL init_masquv(ngrid,zmasq)
 
-
-         endif
-
-
-!        initialize soil 
-!        ~~~~~~~~~~~~~~~
+         endif ! end of 'ok_slab_ocean'.
+
+
+!        Initialize soil.
+!        ~~~~~~~~~~~~~~~~
          if (callsoil) then
+         
             call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
-                ptimestep,tsurf,tsoil,capcal,fluxgrd)
+                      ptimestep,tsurf,tsoil,capcal,fluxgrd)
 
             if (ok_slab_ocean) then
                do ig=1,ngrid
                   if (nint(rnat(ig)).eq.2) then
-                    capcal(ig)=capcalocean
-                    if (pctsrf_sic(ig).gt.0.5) then
-                      capcal(ig)=capcalseaice
-                      if (qsurf(ig,igcm_h2o_ice).gt.0.) then
-                        capcal(ig)=capcalsno
-                      endif
-                    endif
+                     capcal(ig)=capcalocean
+                     if (pctsrf_sic(ig).gt.0.5) then
+                        capcal(ig)=capcalseaice
+                        if (qsurf(ig,igcm_h2o_ice).gt.0.) then
+                           capcal(ig)=capcalsno
+                        endif
+                     endif
                   endif
                enddo
-            endif
-
-
-
-
-         else
+            endif ! end of 'ok_slab_ocean'.
+
+         else ! else of 'callsoil'.
+         
             print*,'WARNING! Thermal conduction in the soil turned off'
             capcal(:)=1.e6
             fluxgrd(:)=intheat
             print*,'Flux from ground = ',intheat,' W m^-2'
-         endif
+            
+         endif ! end of 'callsoil'.
+         
          icount=1
 
-!        decide whether to update ice at end of run
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!        Decide whether to update ice at end of run.
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          ice_update=.false.
          if(sourceevol)then
@@ -635,7 +596,7 @@
                      status="old",iostat=ierr)
             if (ierr.ne.0) then
-              write(*,*) "physiq: Error! No num_run file!"
-              write(*,*) " (which is needed for sourceevol option)"
-              stop
+               write(*,*) "physiq: Error! No num_run file!"
+               write(*,*) " (which is needed for sourceevol option)"
+               stop
             endif
             read(128,*) num_run 
@@ -645,32 +606,35 @@
         
             if(num_run.ne.0.and.mod(num_run,2).eq.0)then
-            !if(num_run.ne.0.and.mod(num_run,3).eq.0)then
                print*,'Updating ice at end of this year!'
                ice_update=.true.
                ice_min(:)=1.e4
-            endif  
-         endif
-
-!  In newstart now, will have to be remove (BC 2014)
-!        define surface as continent or ocean
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            endif
+            
+         endif ! end of 'sourceevol'.
+
+
+         ! Here is defined the type of the surface : Continent or Ocean.
+         ! BC2014 : This is now already done in newstart.
+         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          if (.not.ok_slab_ocean) then
+         
            rnat(:)=1.
            do ig=1,ngrid
-!             if(iceball.or.oceanball.or.(inertiedat(ig,1).gt.1.E4))then
-             if(inertiedat(ig,1).gt.1.E4)then
-               rnat(ig)=0
-             endif
+              if(inertiedat(ig,1).gt.1.E4)then
+                 rnat(ig)=0
+              endif
            enddo
 
            print*,'WARNING! Surface type currently decided by surface inertia'
            print*,'This should be improved e.g. in newstart.F'
-         endif!(.not.ok_slab_ocean)
-
-!        initialise surface history variable 
+           
+         endif ! end of 'ok_slab_ocean'.
+
+
+!        Initialize surface history variable.
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          qsurf_hist(:,:)=qsurf(:,:)
 
-!        initialise variable for dynamical heating diagnostic
+!        Initialize variable for dynamical heating diagnostic.
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          ztprevious(:,:)=pt(:,:)
@@ -688,17 +652,14 @@
          endif
 
-         if(meanOLR)then
-            ! to record global radiative balance
-            call system('rm -f rad_bal.out')
-            ! to record global mean/max/min temperatures
-            call system('rm -f tem_bal.out')
-            ! to record global hydrological balance
-            call system('rm -f h2o_bal.out')
-         endif
-
-         watertest=.false.
-         if(water)then
-            ! initialise variables for water cycle
-
+         if(meanOLR)then          
+            call system('rm -f rad_bal.out') ! to record global radiative balance.           
+            call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.            
+            call system('rm -f h2o_bal.out') ! to record global hydrological balance.
+         endif
+
+
+         watertest=.false.        
+         if(water)then ! Initialize variables for water cycle
+            
             if(enertest)then
                watertest = .true.
@@ -710,42 +671,41 @@
 
          endif
+         
          call su_watercycle ! even if we don't have a water cycle, we might
                             ! need epsi for the wvp definitions in callcorrk.F
 
-         if (ngrid.ne.1) then ! no need to create a restart file in 1d
-           call physdem0("restartfi.nc",long,lati,nsoilmx,ngrid,nlayer,nq, &
-                         ptimestep,pday+nday,time_phys,area, &
-                         albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
+         if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
+            call physdem0("restartfi.nc",long,lati,nsoilmx,ngrid,nlayer,nq, &
+                          ptimestep,pday+nday,time_phys,area, &
+                          albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
          endif
          
-      endif        !  (end of "if firstcall")
-
-! ---------------------------------------------------
-! 1.2   Initializations done at every physical timestep:
-! ---------------------------------------------------
-!     Initialize various variables
-!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      
+      endif ! end of 'firstcall'
+
+! ------------------------------------------------------
+! I.2   Initializations done at every physical timestep:
+! ------------------------------------------------------
+
+      ! Initialize various variables
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      
+      if ( .not.nearco2cond ) then
+         pdt(1:ngrid,1:nlayer) = 0.0
+      endif      
+      zdtsurf(1:ngrid)      = 0.0
+      pdq(1:ngrid,1:nlayer,1:nq) = 0.0
+      dqsurf(1:ngrid,1:nq)= 0.0
       pdu(1:ngrid,1:nlayer) = 0.0
       pdv(1:ngrid,1:nlayer) = 0.0
-      if ( .not.nearco2cond ) then
-         pdt(1:ngrid,1:nlayer) = 0.0
-      endif
-      pdq(1:ngrid,1:nlayer,1:nq) = 0.0
       pdpsrf(1:ngrid)       = 0.0
-      zflubid(1:ngrid)      = 0.0
-      zdtsurf(1:ngrid)      = 0.0
-      dqsurf(1:ngrid,1:nq)= 0.0
+      zflubid(1:ngrid)      = 0.0      
       flux_sens_lat(1:ngrid) = 0.0
       taux(1:ngrid) = 0.0
       tauy(1:ngrid) = 0.0
 
-
-
-
-      zday=pday+ptime ! compute time, in sols (and fraction thereof)
-
-!     Compute Stellar Longitude (Ls), and orbital parameters
-!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      zday=pday+ptime ! Compute time, in sols (and fraction thereof).
+
+      ! Compute Stellar Longitude (Ls), and orbital parameters.
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       if (season) then
          call stellarlong(zday,zls)
@@ -755,4 +715,5 @@
 
       call orbite(zls,dist_star,declin,right_ascen)
+      
       if (tlocked) then
               zlss=Mod(-(2.*pi*(zday/year_day)*nres - right_ascen),2.*pi)
@@ -764,34 +725,27 @@
 
 
-!    Compute variations of g with latitude (oblate case)
-!
-        if (oblate .eqv. .false.) then
-           glat(:) = g
-        else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then
-        print*,'I need values for flatten, J2, Rmean and MassPlanet to compute glat (else set oblate=.false.)'
-
-        call abort
-        else
-
-        gmplanet = MassPlanet*grav*1e24
-        do ig=1,ngrid
-           glat(ig)= gmplanet/(Rmean**2) * (1.D0 + 0.75 *J2 - 2.0*flatten/3. + (2.*flatten - 15./4.* J2) * cos(2. * (pi/2. - lati(ig)))) 
-        end do
-        endif
-
-!!      write(*,*) 'lati, glat =',lati(1)*180./pi,glat(1), lati(ngrid/2)*180./pi, glat(ngrid/2)
-
-
-
-!     Compute geopotential between layers
-!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
+      ! Compute variations of g with latitude (oblate case).
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      if (oblate .eqv. .false.) then      
+         glat(:) = g         
+      else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then      
+         print*,'I need values for flatten, J2, Rmean and MassPlanet to compute glat (else set oblate=.false.)'
+         call abort
+      else
+         gmplanet = MassPlanet*grav*1e24
+         do ig=1,ngrid
+            glat(ig)= gmplanet/(Rmean**2) * (1.D0 + 0.75 *J2 - 2.0*flatten/3. + (2.*flatten - 15./4.* J2) * cos(2. * (pi/2. - lati(ig)))) 
+         end do
+      endif
+
+      ! Compute geopotential between layers.
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
       do l=1,nlayer         
-      zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
+         zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
       enddo
 
       zzlev(1:ngrid,1)=0.
-      zzlev(1:ngrid,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
+      zzlev(1:ngrid,nlayer+1)=1.e7 ! Dummy top of last layer above 10000 km...
 
       do l=2,nlayer
@@ -801,10 +755,9 @@
             zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
          enddo
-      enddo
-!     Potential temperature calculation may not be the same in physiq and dynamic...
-
-!     Compute potential temperature
-!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
+      enddo     
+
+      ! Compute potential temperature
+      ! Note : Potential temperature calculation may not be the same in physiq and dynamic...
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       do l=1,nlayer         
          do ig=1,ngrid
@@ -818,49 +771,47 @@
      ! Compute vertical velocity (m/s) from vertical mass flux
      ! w = F / (rho*area) and rho = P/(r*T)
-     ! but first linearly interpolate mass flux to mid-layers
-     do l=1,nlayer-1
-       pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
-     enddo
-     pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
-     do l=1,nlayer
-       pw(1:ngrid,l)=(pw(1:ngrid,l)*r*pt(1:ngrid,l)) /  &
-                     (pplay(1:ngrid,l)*area(1:ngrid))
-     enddo
-
-!-----------------------------------------------------------------------
-!    2. Compute radiative tendencies
-!-----------------------------------------------------------------------
+     ! But first linearly interpolate mass flux to mid-layers
+      do l=1,nlayer-1
+         pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
+      enddo
+      pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
+      do l=1,nlayer
+         pw(1:ngrid,l)=(pw(1:ngrid,l)*r*pt(1:ngrid,l)) /  &
+                       (pplay(1:ngrid,l)*area(1:ngrid))
+      enddo
+
+!---------------------------------
+! II. Compute radiative tendencies
+!---------------------------------
 
       if (callrad) then
          if( mod(icount-1,iradia).eq.0.or.lastcall) then
 
-!          Compute local stellar zenith angles
-!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-           if (tlocked) then
-! JL14 corrects tidally resonant (and inclined) cases. nres=omega_rot/omega_orb
-              ztim1=SIN(declin)
-              ztim2=COS(declin)*COS(zlss)
-              ztim3=COS(declin)*SIN(zlss)
-
-              call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
-               ztim1,ztim2,ztim3,mu0,fract, flatten)
-
-           elseif (diurnal) then
-              ztim1=SIN(declin)
-              ztim2=COS(declin)*COS(2.*pi*(zday-.5))
-              ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
-
-              call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
-               ztim1,ztim2,ztim3,mu0,fract, flatten)
-           else if(diurnal .eqv. .false.) then
+          ! Compute local stellar zenith angles
+          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            if (tlocked) then
+            ! JL14 corrects tidally resonant (and inclined) cases. nres=omega_rot/omega_orb
+               ztim1=SIN(declin)
+               ztim2=COS(declin)*COS(zlss)
+               ztim3=COS(declin)*SIN(zlss)
+
+               call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
+                            ztim1,ztim2,ztim3,mu0,fract, flatten)
+
+            elseif (diurnal) then
+               ztim1=SIN(declin)
+               ztim2=COS(declin)*COS(2.*pi*(zday-.5))
+               ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
+
+               call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
+                            ztim1,ztim2,ztim3,mu0,fract, flatten)
+            else if(diurnal .eqv. .false.) then
 
                call mucorr(ngrid,declin,lati,mu0,fract,10000.,rad,flatten)
                ! WARNING: this function appears not to work in 1D
 
-           endif 
+            endif 
            
-           !! Eclipse incoming sunlight (e.g. Saturn ring shadowing)
-        
+            ! Eclipse incoming sunlight (e.g. Saturn ring shadowing).       
             if(rings_shadow) then
                 call call_rings(ngrid, ptime, pday, diurnal)
@@ -868,138 +819,132 @@
 
 
-           if (corrk) then
-
-!          a) Call correlated-k radiative transfer scheme
-!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-              if(kastprof)then
-                 print*,'kastprof should not = true here'
-                 call abort
-              endif
-	      if(water) then
+            if (corrk) then
+            
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! II.a Call correlated-k radiative transfer scheme
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+               if(kastprof)then
+                  print*,'kastprof should not = true here'
+                  call abort
+               endif
+	       if(water) then
                   muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayer,igcm_h2o_vap)) 
 		  muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayer,igcm_h2o_vap)) 
 		  ! take into account water effect on mean molecular weight
-	      else
+	       else
                   muvar(1:ngrid,1:nlayer+1)=mugaz
-	      endif 
-      
-!             standard callcorrk
-
-              if(ok_slab_ocean) then
-                tsurf2(:)=tsurf(:) 
-                do ig=1,ngrid
-                  if (nint(rnat(ig))==0) then
-                    tsurf(ig)=((1.-pctsrf_sic(ig))*tslab(ig,1)**4+pctsrf_sic(ig)*tsea_ice(ig)**4)**0.25
-                  endif
-                enddo
-              endif!(ok_slab_ocean)
+	       endif 
+      
+
+               if(ok_slab_ocean) then
+                  tsurf2(:)=tsurf(:) 
+                  do ig=1,ngrid
+                     if (nint(rnat(ig))==0) then
+                        tsurf(ig)=((1.-pctsrf_sic(ig))*tslab(ig,1)**4+pctsrf_sic(ig)*tsea_ice(ig)**4)**0.25
+                     endif
+                  enddo
+               endif !(ok_slab_ocean)
                 
-
-                clearsky=.false.
-                call callcorrk(ngrid,nlayer,pq,nq,qsurf,                   &
-                      albedo,emis,mu0,pplev,pplay,pt,                    &
-                      tsurf,fract,dist_star,aerosol,muvar,               &
-                      zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
-                      fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,               &
-                      tau_col,cloudfrac,totcloudfrac,                    &
-                      clearsky,firstcall,lastcall)      
-
-!             Option to call scheme once more for clear regions
-                if(CLFvarying)then
-
-                 ! ---> PROBLEMS WITH ALLOCATED ARRAYS 
-                 ! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
-                   clearsky=.true.
-                   call callcorrk(ngrid,nlayer,pq,nq,qsurf,                  &
-                      albedo,emis,mu0,pplev,pplay,pt,                      &
-                      tsurf,fract,dist_star,aerosol,muvar,                 &
-                      zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
-                      fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,              &
-                      tau_col1,cloudfrac,totcloudfrac,                     &
-                      clearsky,firstcall,lastcall)
-                   clearsky = .false.  ! just in case.     
-
-
-                 ! Sum the fluxes and heating rates from cloudy/clear cases
-                   do ig=1,ngrid
-                      tf=totcloudfrac(ig)
-                      ntf=1.-tf
+               ! standard callcorrk
+               clearsky=.false.
+               call callcorrk(ngrid,nlayer,pq,nq,qsurf,                          &
+                              albedo,emis,mu0,pplev,pplay,pt,                    &
+                              tsurf,fract,dist_star,aerosol,muvar,               &
+                              zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
+                              fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,               &
+                              tau_col,cloudfrac,totcloudfrac,                    &
+                              clearsky,firstcall,lastcall)      
+
+                ! Option to call scheme once more for clear regions
+               if(CLFvarying)then
+
+                  ! ---> PROBLEMS WITH ALLOCATED ARRAYS : temporary solution in callcorrk: do not deallocate if CLFvarying ...
+                  clearsky=.true.
+                  call callcorrk(ngrid,nlayer,pq,nq,qsurf,                            &
+                                 albedo,emis,mu0,pplev,pplay,pt,                      &
+                                 tsurf,fract,dist_star,aerosol,muvar,                 &
+                                 zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
+                                 fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,              &
+                                 tau_col1,cloudfrac,totcloudfrac,                     &
+                                 clearsky,firstcall,lastcall)
+                  clearsky = .false.  ! just in case.     
+
+
+                  ! Sum the fluxes and heating rates from cloudy/clear cases
+                  do ig=1,ngrid
+                     tf=totcloudfrac(ig)
+                     ntf=1.-tf
                     
-                    fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw(ig)
-                    fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw(ig)
-                    fluxtop_lw(ig)  = ntf*fluxtop_lw1(ig)  + tf*fluxtop_lw(ig)
-                    fluxabs_sw(ig)  = ntf*fluxabs_sw1(ig)  + tf*fluxabs_sw(ig)
-                    tau_col(ig)     = ntf*tau_col1(ig)     + tf*tau_col(ig)
+                     fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw(ig)
+                     fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw(ig)
+                     fluxtop_lw(ig)  = ntf*fluxtop_lw1(ig)  + tf*fluxtop_lw(ig)
+                     fluxabs_sw(ig)  = ntf*fluxabs_sw1(ig)  + tf*fluxabs_sw(ig)
+                     tau_col(ig)     = ntf*tau_col1(ig)     + tf*tau_col(ig)
                     
-                    zdtlw(ig,1:nlayer) = ntf*zdtlw1(ig,1:nlayer) + tf*zdtlw(ig,1:nlayer)
-                    zdtsw(ig,1:nlayer) = ntf*zdtsw1(ig,1:nlayer) + tf*zdtsw(ig,1:nlayer)
-
-                    OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV)		       
-                    OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI)		       
-
+                     zdtlw(ig,1:nlayer) = ntf*zdtlw1(ig,1:nlayer) + tf*zdtlw(ig,1:nlayer)
+                     zdtsw(ig,1:nlayer) = ntf*zdtsw1(ig,1:nlayer) + tf*zdtsw(ig,1:nlayer)
+
+                     OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV)		       
+                     OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI)		       
                  enddo
 
-              endif !CLFvarying
-
-              if(ok_slab_ocean) then
-                tsurf(:)=tsurf2(:)
-              endif!(ok_slab_ocean)
-
-
-!             Radiative flux from the sky absorbed by the surface (W.m-2)
-              GSR=0.0
-              fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+fluxsurf_sw(1:ngrid)*(1.-albedo(1:ngrid))
-
-              !if(noradsurf)then ! no lower surface; SW flux just disappears
-              !   GSR = SUM(fluxsurf_sw(1:ngrid)*area(1:ngrid))/totarea
-              !   fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
-              !   print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
-              !endif
-
-!             Net atmospheric radiative heating rate (K.s-1)
-              dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
-
-           elseif(newtonian)then
-
-!          b) Call Newtonian cooling scheme
-!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-              call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
-
-              zdtsurf(1:ngrid) = +(pt(1:ngrid,1)-tsurf(1:ngrid))/ptimestep
-              ! e.g. surface becomes proxy for 1st atmospheric layer ?
-
-           else
-
-!          c) Atmosphere has no radiative effect
-!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-              fluxtop_dn(1:ngrid)  = fract(1:ngrid)*mu0(1:ngrid)*Fat1AU/dist_star**2
-              if(ngrid.eq.1)then ! / by 4 globally in 1D case...
-                 fluxtop_dn(1)  = fract(1)*Fat1AU/dist_star**2/2.0
-              endif
-              fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
-              fluxrad_sky(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid))
-              fluxtop_lw(1:ngrid)  = emis(1:ngrid)*sigma*tsurf(1:ngrid)**4
-              ! radiation skips the atmosphere entirely
-
-
-              dtrad(1:ngrid,1:nlayer)=0.0
-              ! hence no atmospheric radiative heating
-
-           endif                ! if corrk
-
-        endif ! of if(mod(icount-1,iradia).eq.0)
+               endif ! end of CLFvarying.
+
+               if(ok_slab_ocean) then
+                  tsurf(:)=tsurf2(:)
+               endif
+
+
+               ! Radiative flux from the sky absorbed by the surface (W.m-2)
+               GSR=0.0
+               fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+fluxsurf_sw(1:ngrid)*(1.-albedo(1:ngrid))
+
+                            !if(noradsurf)then ! no lower surface; SW flux just disappears
+                            !   GSR = SUM(fluxsurf_sw(1:ngrid)*area(1:ngrid))/totarea
+                            !   fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
+                            !   print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
+                            !endif
+
+               ! Net atmospheric radiative heating rate (K.s-1)
+               dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
+
+            elseif(newtonian)then
+
+               ! II.b Call Newtonian cooling scheme
+               ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+               call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
+
+               zdtsurf(1:ngrid) = +(pt(1:ngrid,1)-tsurf(1:ngrid))/ptimestep
+               ! e.g. surface becomes proxy for 1st atmospheric layer ?
+
+            else
+            
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! II.c Atmosphere has no radiative effect 
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+               fluxtop_dn(1:ngrid)  = fract(1:ngrid)*mu0(1:ngrid)*Fat1AU/dist_star**2
+               if(ngrid.eq.1)then ! / by 4 globally in 1D case...
+                  fluxtop_dn(1)  = fract(1)*Fat1AU/dist_star**2/2.0
+               endif
+               fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
+               fluxrad_sky(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid))
+               fluxtop_lw(1:ngrid)  = emis(1:ngrid)*sigma*tsurf(1:ngrid)**4
+
+               dtrad(1:ngrid,1:nlayer)=0.0 ! no atmospheric radiative heating
+
+            endif ! end of corrk
+
+         endif ! of if(mod(icount-1,iradia).eq.0)
        
 
-!    Transformation of the radiative tendencies
-!    ------------------------------------------
-
-        zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
-        zplanck(1:ngrid)=emis(1:ngrid)*sigma*zplanck(1:ngrid)*zplanck(1:ngrid)
-        fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
-        pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
-
-!-------------------------
-! test energy conservation
+         ! Transformation of the radiative tendencies
+         ! ------------------------------------------
+         zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
+         zplanck(1:ngrid)=emis(1:ngrid)*sigma*zplanck(1:ngrid)*zplanck(1:ngrid)
+         fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
+         pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
+         
+         ! Test of energy conservation
+         !----------------------------
          if(enertest)then
 	    call planetwide_sumval(cpp*massarea(:,:)*zdtsw(:,:)/totarea_planet,dEtotSW)
@@ -1010,20 +955,21 @@
 	    dEzRadlw(:,:)=cpp*mass(:,:)*zdtlw(:,:)
 	    if (is_master) then
-            	print*,'---------------------------------------------------------------'
-            	print*,'In corrk SW atmospheric heating       =',dEtotSW,' W m-2'
-            	print*,'In corrk LW atmospheric heating       =',dEtotLW,' W m-2'
-            	print*,'atmospheric net rad heating (SW+LW)   =',dEtotLW+dEtotSW,' W m-2'
-            	print*,'In corrk SW surface heating           =',dEtotsSW,' W m-2'
-            	print*,'In corrk LW surface heating           =',dEtotsLW,' W m-2'
-            	print*,'surface net rad heating (SW+LW)       =',dEtotsLW+dEtotsSW,' W m-2'
+               print*,'---------------------------------------------------------------'
+               print*,'In corrk SW atmospheric heating       =',dEtotSW,' W m-2'
+               print*,'In corrk LW atmospheric heating       =',dEtotLW,' W m-2'
+               print*,'atmospheric net rad heating (SW+LW)   =',dEtotLW+dEtotSW,' W m-2'
+               print*,'In corrk SW surface heating           =',dEtotsSW,' W m-2'
+               print*,'In corrk LW surface heating           =',dEtotsLW,' W m-2'
+               print*,'surface net rad heating (SW+LW)       =',dEtotsLW+dEtotsSW,' W m-2'
 	    endif
-         endif
-!-------------------------
+         endif ! end of 'enertest'
 
       endif ! of if (callrad)
 
-!-----------------------------------------------------------------------
-!    4. Vertical diffusion (turbulent mixing):
-!    -----------------------------------------
+
+
+!  --------------------------------------------
+!  III. Vertical diffusion (turbulent mixing) :
+!  --------------------------------------------
 
       if (calldifv) then
@@ -1031,37 +977,34 @@
          zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
 
-         zdum1(1:ngrid,1:nlayer)=0.0
-         zdum2(1:ngrid,1:nlayer)=0.0
-
-
-!JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff	 
+         ! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
 	 if (UseTurbDiff) then
 	 
-           call turbdiff(ngrid,nlayer,nq,rnat,       &
-              ptimestep,capcal,lwrite,               &
-              pplay,pplev,zzlay,zzlev,z0,            &
-              pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf,   &
-              zdum1,zdum2,pdt,pdq,zflubid,           &
-              zdudif,zdvdif,zdtdif,zdtsdif,          &
-	      sensibFlux,q2,zdqdif,zdqevap,zdqsdif,  &
-              taux,tauy,lastcall)
+            call turbdiff(ngrid,nlayer,nq,rnat,                  &
+                          ptimestep,capcal,lwrite,               &
+                          pplay,pplev,zzlay,zzlev,z0,            &
+                          pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf,   &
+                          pdt,pdq,zflubid,                       &
+                          zdudif,zdvdif,zdtdif,zdtsdif,          &
+	                  sensibFlux,q2,zdqdif,zdqevap,zdqsdif,  &
+                          taux,tauy,lastcall)
 
 	 else
 	 
-           zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
+            zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
  
-           call vdifc(ngrid,nlayer,nq,rnat,zpopsk,   &
-              ptimestep,capcal,lwrite,               &
-              pplay,pplev,zzlay,zzlev,z0,            &
-              pu,pv,zh,pq,tsurf,emis,qsurf,          &
-              zdum1,zdum2,zdh,pdq,zflubid,           &
-              zdudif,zdvdif,zdhdif,zdtsdif,          &
-	      sensibFlux,q2,zdqdif,zdqsdif,          &
-              taux,tauy,lastcall)
-
-           zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
-	   zdqevap(1:ngrid,1:nlayer)=0.
-
-         end if !turbdiff
+            call vdifc(ngrid,nlayer,nq,rnat,zpopsk,           &
+                       ptimestep,capcal,lwrite,               &
+                       pplay,pplev,zzlay,zzlev,z0,            &
+                       pu,pv,zh,pq,tsurf,emis,qsurf,          &
+                       zdh,pdq,zflubid,                       &
+                       zdudif,zdvdif,zdhdif,zdtsdif,          &
+	               sensibFlux,q2,zdqdif,zdqsdif,          &
+                       taux,tauy,lastcall)
+
+            zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
+	    zdqevap(1:ngrid,1:nlayer)=0.
+
+         end if !end of 'UseTurbDiff'
+
 
          pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)
@@ -1070,8 +1013,8 @@
          zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
 
+
          if(ok_slab_ocean)then
             flux_sens_lat(1:ngrid)=(zdtsdif(1:ngrid)*capcal(1:ngrid)-fluxrad(1:ngrid))
          endif
-
 
 
@@ -1081,7 +1024,9 @@
          end if ! of if (tracer)
 
+
+         ! test energy conservation
          !-------------------------
-         ! test energy conservation
          if(enertest)then
+         
 	    dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:)
             do ig = 1, ngrid
@@ -1089,32 +1034,35 @@
 	       dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
             enddo
+            
 	    call planetwide_sumval(dEdiff(:)*area(:)/totarea_planet,dEtot)
 	    dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
 	    call planetwide_sumval(dEdiffs(:)*area(:)/totarea_planet,dEtots)
 	    call planetwide_sumval(sensibFlux(:)*area(:)/totarea_planet,AtmToSurf_TurbFlux)
+            
 	    if (is_master) then
-             if (UseTurbDiff) then
-	       	print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
-               	print*,'In TurbDiff non-cons atm nrj change   =',dEtot,' W m-2'
-               	print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
-	     else
-	       	print*,'In vdifc sensible flux (atm=>surf)    =',AtmToSurf_TurbFlux,' W m-2'
-               	print*,'In vdifc non-cons atm nrj change      =',dEtot,' W m-2'
-               	print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
-	     end if
-	    endif ! of if (is_master)
-! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme
-!    but not given back elsewhere
-         endif
-         !-------------------------
-
-         !-------------------------
-         ! test water conservation
+            
+               if (UseTurbDiff) then
+	       	  print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
+               	  print*,'In TurbDiff non-cons atm nrj change   =',dEtot,' W m-2'
+                  print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
+	       else
+	       	  print*,'In vdifc sensible flux (atm=>surf)    =',AtmToSurf_TurbFlux,' W m-2'
+               	  print*,'In vdifc non-cons atm nrj change      =',dEtot,' W m-2'
+               	  print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
+	       end if
+	    endif ! end of 'is_master'
+            
+         ! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
+         endif ! end of 'enertest'
+
+
+         ! Test water conservation.
          if(watertest.and.water)then
+         
 	    call planetwide_sumval(massarea(:,:)*zdqdif(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dWtot_tmp)
 	    call planetwide_sumval(zdqsdif(:,igcm_h2o_vap)*area(:)*ptimestep/totarea_planet,dWtots_tmp)
             do ig = 1, ngrid
 	       vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap))
-	    Enddo
+	    enddo
 	    call planetwide_sumval(massarea(:,:)*zdqdif(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
 	    call planetwide_sumval(zdqsdif(:,igcm_h2o_ice)*area(:)*ptimestep/totarea_planet,dWtots)
@@ -1123,19 +1071,19 @@
             do ig = 1, ngrid
 	       vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice))
-	    Enddo	    
+	    enddo	    
 	    call planetwide_maxval(vdifcncons(:),nconsMAX)
 
 	    if (is_master) then
-            	print*,'---------------------------------------------------------------'
-            	print*,'In difv atmospheric water change        =',dWtot,' kg m-2'
-            	print*,'In difv surface water change            =',dWtots,' kg m-2'
-            	print*,'In difv non-cons factor                 =',dWtot+dWtots,' kg m-2'
-            	print*,'In difv MAX non-cons factor             =',nconsMAX,' kg m-2 s-1'
+               print*,'---------------------------------------------------------------'
+               print*,'In difv atmospheric water change        =',dWtot,' kg m-2'
+               print*,'In difv surface water change            =',dWtots,' kg m-2'
+               print*,'In difv non-cons factor                 =',dWtot+dWtots,' kg m-2'
+               print*,'In difv MAX non-cons factor             =',nconsMAX,' kg m-2 s-1'
 	    endif
 
-         endif
+         endif ! end of 'watertest'
          !-------------------------
 
-      else    
+      else ! calldifv
 
          if(.not.newtonian)then
@@ -1145,10 +1093,10 @@
          endif
 
-      endif ! of if (calldifv)
-
-
-!-----------------------------------------------------------------------
-!   5. Dry convective adjustment:
-!   -----------------------------
+      endif ! end of 'calldifv'
+
+
+!----------------------------------
+!   IV. Dry convective adjustment :
+!----------------------------------
 
       if(calladj) then
@@ -1160,10 +1108,10 @@
 
 
-         call convadj(ngrid,nlayer,nq,ptimestep,    &
-              pplay,pplev,zpopsk,                   &
-              pu,pv,zh,pq,                          &
-              pdu,pdv,zdh,pdq,                      &
-              zduadj,zdvadj,zdhadj,                 &
-              zdqadj)
+         call convadj(ngrid,nlayer,nq,ptimestep,            &
+                      pplay,pplev,zpopsk,                   &
+                      pu,pv,zh,pq,                          &
+                      pdu,pdv,zdh,pdq,                      &
+                      zduadj,zdvadj,zdhadj,                 &
+                      zdqadj)
 
          pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer)
@@ -1176,50 +1124,48 @@
          end if
 
-         !-------------------------
-         ! test energy conservation
+         ! Test energy conservation
          if(enertest)then
 	    call planetwide_sumval(cpp*massarea(:,:)*zdtadj(:,:)/totarea_planet,dEtot)
             if (is_master) print*,'In convadj atmospheric energy change  =',dEtot,' W m-2'
          endif
-         !-------------------------
-
-         !-------------------------
-         ! test water conservation
+
+         ! Test water conservation
          if(watertest)then
 	    call planetwide_sumval(massarea(:,:)*zdqadj(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dWtot_tmp)
             do ig = 1, ngrid
 	       cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap))
-	    Enddo
+	    enddo
 	    call planetwide_sumval(massarea(:,:)*zdqadj(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
 	    dWtot = dWtot + dWtot_tmp
             do ig = 1, ngrid
 	       cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice))
-	    Enddo	    
+	    enddo	    
 	    call planetwide_maxval(cadjncons(:),nconsMAX)
 
             if (is_master) then
-	    	print*,'In convadj atmospheric water change     =',dWtot,' kg m-2'
-            	print*,'In convadj MAX non-cons factor          =',nconsMAX,' kg m-2 s-1'
+	       print*,'In convadj atmospheric water change     =',dWtot,' kg m-2'
+               print*,'In convadj MAX non-cons factor          =',nconsMAX,' kg m-2 s-1'
 	    endif
-         endif
-         !-------------------------
+            
+         endif ! end of 'watertest'
          
-      endif ! of if(calladj)
-
-!-----------------------------------------------------------------------
-!   6. Carbon dioxide condensation-sublimation:
-!   -------------------------------------------
+      endif ! end of 'calladj'
+
+!-----------------------------------------------
+!   V. Carbon dioxide condensation-sublimation :
+!-----------------------------------------------
 
       if (co2cond) then
+      
          if (.not.tracer) then
             print*,'We need a CO2 ice tracer to condense CO2'
             call abort
          endif
-         call condense_cloud(ngrid,nlayer,nq,ptimestep,   &
-              capcal,pplay,pplev,tsurf,pt,                  &
-              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,        &
-              qsurf(1,igcm_co2_ice),albedo,emis,            &
-              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,               &
-              zdqc)
+         call condense_co2(ngrid,nlayer,nq,ptimestep,                    &
+                           capcal,pplay,pplev,tsurf,pt,                  &
+                           pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,        &
+                           qsurf(1,igcm_co2_ice),albedo,emis,            &
+                           zdtc,zdtsurfc,pdpsrf,zduc,zdvc,               &
+                           zdqc)
 
          pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer)
@@ -1229,8 +1175,5 @@
 
          pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq) 
-         ! Note: we do not add surface co2ice tendency
-         ! because qsurf(:,igcm_co2_ice) is updated in condens_co2cloud
-
-         !-------------------------
+
          ! test energy conservation
          if(enertest)then
@@ -1238,46 +1181,42 @@
 	    call planetwide_sumval(capcal(:)*zdtsurfc(:)*area(:)/totarea_planet,dEtots)
 	    if (is_master) then
-	    	print*,'In co2cloud atmospheric energy change   =',dEtot,' W m-2'
-            	print*,'In co2cloud surface energy change       =',dEtots,' W m-2'
+	       print*,'In co2cloud atmospheric energy change   =',dEtot,' W m-2'
+               print*,'In co2cloud surface energy change       =',dEtots,' W m-2'
 	    endif
          endif
-         !-------------------------
-
-      endif  ! of if (co2cond)
-
-
-!-----------------------------------------------------------------------
-!   7. Specific parameterizations for tracers 
-!   -----------------------------------------
-
-      if (tracer) then 
-
-!   7a. Water and ice
-!     ---------------
+
+      endif  ! end of 'co2cond'
+
+
+!---------------------------------------------
+!   VI. Specific parameterizations for tracers 
+!---------------------------------------------
+
+      if (tracer) then
+      
+  ! ---------------------
+  !   VI.1. Water and ice
+  ! ---------------------
          if (water) then
 
-!        ----------------------------------------
-!        Water ice condensation in the atmosphere
-!        ----------------------------------------
+            ! Water ice condensation in the atmosphere
             if(watercond.and.(RLVTT.gt.1.e-8))then
-
-!             ----------------
-!             Moist convection
-!             ----------------
 
 	       dqmoist(1:ngrid,1:nlayer,1:nq)=0.
 	       dtmoist(1:ngrid,1:nlayer)=0.
-
+               
+               ! Moist Convective Adjustment.
+               ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 	       call moistadj(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
 
-               pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap)   &
-	                  +dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)
-               pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =pdq(1:ngrid,1:nlayer,igcm_h2o_ice)     &
-                          +dqmoist(1:ngrid,1:nlayer,igcm_h2o_ice)
+               pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap)     &
+	                                          + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)
+               pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice)     &
+                                                  + dqmoist(1:ngrid,1:nlayer,igcm_h2o_ice)
                pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtmoist(1:ngrid,1:nlayer)
 
-               !-------------------------
-               ! test energy conservation
+               ! Test energy conservation.
                if(enertest)then
+               
 		  call planetwide_sumval(cpp*massarea(:,:)*dtmoist(:,:)/totarea_planet,dEtot)
 		  call planetwide_maxval(dtmoist(:,:),dtmoist_max)
@@ -1287,21 +1226,20 @@
                      madjdE(ig) = cpp*SUM(mass(:,:)*dtmoist(:,:))
                   enddo
+                  
 		  if (is_master) then
-                  	print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
-                  	print*,'In moistadj MAX atmospheric energy change   =',dtmoist_max*ptimestep,'K/step'
-                  	print*,'In moistadj MIN atmospheric energy change   =',dtmoist_min*ptimestep,'K/step'
+                     print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
+                     print*,'In moistadj MAX atmospheric energy change   =',dtmoist_max*ptimestep,'K/step'
+                     print*,'In moistadj MIN atmospheric energy change   =',dtmoist_min*ptimestep,'K/step'
 		  endif
 		  
-		! test energy conservation
 		  call planetwide_sumval(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+	&
-		  	massarea(:,:)*dqmoist(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
+		  	                 massarea(:,:)*dqmoist(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
 		  if (is_master) print*,'In moistadj atmospheric water change    =',dWtot,' kg m-2'
-               endif
-               !-------------------------
+                  
+               endif ! end of 'enertest'
 	       
 
-!        --------------------------------
-!        Large scale condensation/evaporation
-!        --------------------------------
+               ! Large scale condensation/evaporation.
+               ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                call largescale(ngrid,nlayer,nq,ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc)
 
@@ -1310,6 +1248,5 @@
                pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice)+dqcldlscale(1:ngrid,1:nlayer)
 
-               !-------------------------
-               ! test energy conservation
+               ! Test energy conservation.
                if(enertest)then
                   lscaledEz(:,:) = cpp*mass(:,:)*dtlscale(:,:)
@@ -1318,18 +1255,15 @@
                   enddo
 		  call planetwide_sumval(cpp*massarea(:,:)*dtlscale(:,:)/totarea_planet,dEtot)
-!		  if(isnan(dEtot)) then ! NB: isnan() is not a standard function...
-!		     print*,'Nan in largescale, abort'
-!                     STOP
-!		  endif
+
 		  if (is_master) print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
 
-               ! test water conservation
+                  ! Test water conservation.
 		  call planetwide_sumval(massarea(:,:)*dqvaplscale(:,:)*ptimestep/totarea_planet+	&
-		  	SUM(massarea(:,:)*dqcldlscale(:,:))*ptimestep/totarea_planet,dWtot)
+		  	                 SUM(massarea(:,:)*dqcldlscale(:,:))*ptimestep/totarea_planet,dWtot)
+                        
 		  if (is_master) print*,'In largescale atmospheric water change  =',dWtot,' kg m-2'
-               endif
-               !-------------------------
-
-               ! compute cloud fraction
+               endif ! end of 'enertest'
+
+               ! Compute cloud fraction.
                do l = 1, nlayer
                   do ig=1,ngrid
@@ -1338,10 +1272,8 @@
                enddo
 
-            endif                ! of if (watercondense)
+            endif ! end of 'watercond'
            
-
-!        --------------------------------
-!        Water ice / liquid precipitation
-!        --------------------------------
+            ! Water ice / liquid precipitation.
+            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
             if(waterrain)then
 
@@ -1351,18 +1283,18 @@
 
                call rain(ngrid,nlayer,nq,ptimestep,pplev,pplay,pt,pdt,pq,pdq,            &
-                   zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
+                         zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
 
                pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) &
-	             +zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)
+	                                          + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)
                pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice) &
-	             +zdqrain(1:ngrid,1:nlayer,igcm_h2o_ice)
+	                                          + zdqrain(1:ngrid,1:nlayer,igcm_h2o_ice)
                pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtrain(1:ngrid,1:nlayer)
-               dqsurf(1:ngrid,igcm_h2o_vap) = dqsurf(1:ngrid,igcm_h2o_vap)+zdqsrain(1:ngrid) ! a bug was here
-               dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid)
-               rainout(1:ngrid)             = zdqsrain(1:ngrid)+zdqssnow(1:ngrid) ! diagnostic   
-
-               !-------------------------
-               ! test energy conservation
+               
+               dqsurf(1:ngrid,igcm_h2o_vap) = dqsurf(1:ngrid,igcm_h2o_vap)+zdqsrain(1:ngrid)
+               dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid) 
+
+               ! Test energy conservation.
                if(enertest)then
+               
 		  call planetwide_sumval(cpp*massarea(:,:)*zdtrain(:,:)/totarea_planet,dEtot)
 		  if (is_master) print*,'In rain atmospheric T energy change       =',dEtot,' W m-2'
@@ -1374,14 +1306,16 @@
 		  dVtot = dVtot + dVtot_tmp
 		  dEtot = dItot + dVtot
+                  
 		  if (is_master) then
-                  	print*,'In rain dItot =',dItot,' W m-2'
-                  	print*,'In rain dVtot =',dVtot,' W m-2'
-                  	print*,'In rain atmospheric L energy change       =',dEtot,' W m-2'
+                     print*,'In rain dItot =',dItot,' W m-2'
+                     print*,'In rain dVtot =',dVtot,' W m-2'
+                     print*,'In rain atmospheric L energy change       =',dEtot,' W m-2'
 		  endif
 
-               ! test water conservation
+                  ! Test water conservation
 		  call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+	&
 		  	massarea(:,:)*zdqrain(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
 		  call planetwide_sumval((zdqsrain(:)+zdqssnow(:))*area(:)*ptimestep/totarea_planet,dWtots)
+                  
 		  if (is_master) then
                   	print*,'In rain atmospheric water change        =',dWtot,' kg m-2'
@@ -1389,82 +1323,72 @@
                   	print*,'In rain non-cons factor                 =',dWtot+dWtots,' kg m-2'
 		  endif
-              endif
-              !-------------------------
-
-            end if                 ! of if (waterrain)
-         end if                    ! of if (water)
-
-
-!   7c. Aerosol particles
-!     -------------------
-!        ------------- 
-!        Sedimentation
-!        ------------- 
-        if (sedimentation) then 
-           zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0
-           zdqssed(1:ngrid,1:nq)  = 0.0
-
-
-           !-------------------------
-           ! find qtot
-           if(watertest)then
-              iq=igcm_h2o_ice
-	      call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot)
-	      call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots)
-	      if (is_master) then
-	      	print*,'Before sedim pq  =',dWtot,' kg m-2'
-              	print*,'Before sedim pdq =',dWtots,' kg m-2'
-	      endif
-           endif
-           !-------------------------
-
-           call callsedim(ngrid,nlayer,ptimestep,           &
-                pplev,zzlev,pt,pdt,pq,pdq,zdqsed,zdqssed,nq)
-
-           !-------------------------
-           ! find qtot
-           if(watertest)then
-              iq=igcm_h2o_ice
-	      call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot)
-	      call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots)
-	      if (is_master) then
-              	print*,'After sedim pq  =',dWtot,' kg m-2'
-              	print*,'After sedim pdq =',dWtots,' kg m-2'
-  	      endif
-           endif
-           !-------------------------
-
-           ! for now, we only allow H2O ice to sediment
-              ! and as in rain.F90, whether it falls as rain or snow depends
-              ! only on the surface temperature
-           pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
-           dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
-
-           !-------------------------
-           ! test water conservation
-           if(watertest)then
-	      call planetwide_sumval(massarea(:,:)*(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice))*ptimestep/totarea_planet,dWtot)
-	      call planetwide_sumval((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))*area(:)*ptimestep/totarea_planet,dWtots)
-	      if (is_master) then
-              	print*,'In sedim atmospheric ice change         =',dWtot,' kg m-2'
-              	print*,'In sedim surface ice change             =',dWtots,' kg m-2'
-              	print*,'In sedim non-cons factor                =',dWtot+dWtots,' kg m-2'
-	      endif
-           endif
-           !-------------------------
-
-        end if   ! of if (sedimentation)
-
-
-!   7d. Updates
-!     ---------
-
-!       -----------------------------------
-!       Updating atm mass and tracer budget
-!       -----------------------------------
-
+                  
+               endif ! end of 'enertest'
+
+            end if ! enf of 'waterrain'
+            
+         end if ! end of 'water'
+
+  ! -------------------------
+  !   VI.2. Aerosol particles
+  ! -------------------------
+
+         ! Sedimentation.
+         if (sedimentation) then
+        
+            zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0
+            zdqssed(1:ngrid,1:nq)  = 0.0
+
+            if(watertest)then
+           
+               iq=igcm_h2o_ice
+	       call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot)
+	       call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots)
+	       if (is_master) then
+	      	  print*,'Before sedim pq  =',dWtot,' kg m-2'
+                  print*,'Before sedim pdq =',dWtots,' kg m-2'
+	       endif
+            endif
+           
+            call callsedim(ngrid,nlayer,ptimestep,           &
+                          pplev,zzlev,pt,pdt,pq,pdq,        &
+                          zdqsed,zdqssed,nq)
+
+            if(watertest)then
+               iq=igcm_h2o_ice
+	       call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot)
+	       call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots)
+	       if (is_master) then
+              	  print*,'After sedim pq  =',dWtot,' kg m-2'
+              	  print*,'After sedim pdq =',dWtots,' kg m-2'
+  	       endif
+            endif
+
+            ! Whether it falls as rain or snow depends only on the surface temperature
+            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
+            dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
+
+            ! Test water conservation
+            if(watertest)then
+	       call planetwide_sumval(massarea(:,:)*(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice))*ptimestep/totarea_planet,dWtot)
+	       call planetwide_sumval((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))*area(:)*ptimestep/totarea_planet,dWtots)
+	       if (is_master) then
+              	  print*,'In sedim atmospheric ice change         =',dWtot,' kg m-2'
+              	  print*,'In sedim surface ice change             =',dWtots,' kg m-2'
+              	  print*,'In sedim non-cons factor                =',dWtot+dWtots,' kg m-2'
+	       endif
+            endif
+
+         end if ! end of 'sedimentation'
+
+
+  ! ---------------
+  !   VI.3. Updates
+  ! ---------------
+
+         ! Updating Atmospheric Mass and Tracers budgets.
          if(mass_redistrib) then
 
-            zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * &
+            zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) *     &
 	       (   zdqevap(1:ngrid,1:nlayer)                          &
 	         + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)             &
@@ -1478,109 +1402,95 @@
 	    call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
 	    call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
-	    call writediagfi(ngrid,"mass","mass"," ",3,mass)
-
-	    call mass_redistribution(ngrid,nlayer,nq,ptimestep,   &
-                       rnat,capcal,pplay,pplev,pt,tsurf,pq,qsurf,     &
-		       pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr,  &
-		       zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
+	    call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
+
+	    call mass_redistribution(ngrid,nlayer,nq,ptimestep,                     &
+                                     rnat,capcal,pplay,pplev,pt,tsurf,pq,qsurf,     &
+		                     pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr,        &
+		                     zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
 	 
-
             pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
-            dqsurf(1:ngrid,1:nq)         = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
-            pdt(1:ngrid,1:nlayer)        = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer)
-            pdu(1:ngrid,1:nlayer)        = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
-            pdv(1:ngrid,1:nlayer)        = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
-	    pdpsrf(1:ngrid)                = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
-            zdtsurf(1:ngrid)               = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
+            dqsurf(1:ngrid,1:nq)       = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
+            pdt(1:ngrid,1:nlayer)      = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer)
+            pdu(1:ngrid,1:nlayer)      = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
+            pdv(1:ngrid,1:nlayer)      = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
+	    pdpsrf(1:ngrid)            = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
+            zdtsurf(1:ngrid)           = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
 	    
-!	    print*,'after mass redistrib, q=',pq(211,1:nlayer,igcm_h2o_vap)+ptimestep*pdq(211,1:nlayer,igcm_h2o_vap) 
 	 endif
 
-
-!   7e. Ocean
-!     ---------
-
-!       ---------------------------------
-!       Slab_ocean
-!       ---------------------------------
-      if (ok_slab_ocean)then
-
-         do ig=1,ngrid
-            qsurfint(:,igcm_h2o_ice)=qsurf(:,igcm_h2o_ice)
-         enddo
-
-         call ocean_slab_ice(ptimestep, &
-               ngrid, knindex, tsea_ice, fluxrad, &
-               zdqssnow, qsurf(:,igcm_h2o_ice), &
-               -zdqsdif(:,igcm_h2o_vap), &
-               flux_sens_lat,tsea_ice, pctsrf_sic, &
-               taux,tauy,icount)
-
-
-         call ocean_slab_get_vars(ngrid,tslab, &
-               sea_ice, flux_o, &
-               flux_g, dt_hdiff, &
-               dt_ekman)
-
+  ! ------------------
+  !   VI.4. Slab Ocean
+  ! ------------------
+
+         if (ok_slab_ocean)then
+
+            do ig=1,ngrid
+               qsurfint(:,igcm_h2o_ice)=qsurf(:,igcm_h2o_ice)
+            enddo
+
+            call ocean_slab_ice(ptimestep,                          &
+                                ngrid, knindex, tsea_ice, fluxrad,  &
+                                zdqssnow, qsurf(:,igcm_h2o_ice),    &
+                              - zdqsdif(:,igcm_h2o_vap),            &
+                                flux_sens_lat,tsea_ice, pctsrf_sic, &
+                                taux,tauy,icount)
+
+
+            call ocean_slab_get_vars(ngrid,tslab,      &
+                                     sea_ice, flux_o,  &
+                                     flux_g, dt_hdiff, &
+                                     dt_ekman)
+   
             do ig=1,ngrid
                if (nint(rnat(ig)).eq.1)then
-               tslab(ig,1) = 0.
-               tslab(ig,2) = 0.
-               tsea_ice(ig) = 0.
-               sea_ice(ig) = 0.
-               pctsrf_sic(ig) = 0.
-               qsurf(ig,igcm_h2o_ice)=qsurfint(ig,igcm_h2o_ice)
+                  tslab(ig,1) = 0.
+                  tslab(ig,2) = 0.
+                  tsea_ice(ig) = 0.
+                  sea_ice(ig) = 0.
+                  pctsrf_sic(ig) = 0.
+                  qsurf(ig,igcm_h2o_ice) = qsurfint(ig,igcm_h2o_ice)
                endif
             enddo
 
-
-      endif! (ok_slab_ocean)
-
-!       ---------------------------------
-!       Updating tracer budget on surface
-!       ---------------------------------
+         endif ! end of 'ok_slab_ocean'
+
+  ! -----------------------------
+  !   VI.5. Surface Tracer Update
+  ! -----------------------------
 
          if(hydrology)then
             
-            call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd,    &
-               capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic,  &
-               sea_ice)
-            ! note: for now, also changes albedo in the subroutine
-
-            zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurf_hyd(1:ngrid)
+            call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd,             &
+                        capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic,  &
+                        sea_ice)
+
+            zdtsurf(1:ngrid)    = zdtsurf(1:ngrid) + zdtsurf_hyd(1:ngrid)
             qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqs_hyd(1:ngrid,1:nq)
-            ! when hydrology is used, other dqsurf tendencies are all added to dqs_hyd inside
-
-            !-------------------------
-            ! test energy conservation
+
+            ! Test energy conservation
             if(enertest)then
 	       call planetwide_sumval(area(:)*capcal(:)*zdtsurf_hyd(:)/totarea_planet,dEtots)
 	       if (is_master) print*,'In hydrol surface energy change     =',dEtots,' W m-2'
             endif
-            !-------------------------
-
-            !-------------------------
+
             ! test water conservation
             if(watertest)then
 	       call planetwide_sumval(dqs_hyd(:,igcm_h2o_ice)*area(:)*ptimestep/totarea_planet,dWtots)
 	       if (is_master) print*,'In hydrol surface ice change            =',dWtots,' kg m-2'
-	       call planetwide_sumval(dqs_hyd(:,igcm_h2o_vap)*area(:)*ptimestep/totarea_planet,dWtots)
+	          call planetwide_sumval(dqs_hyd(:,igcm_h2o_vap)*area(:)*ptimestep/totarea_planet,dWtots)
 	       if (is_master) then
-               	print*,'In hydrol surface water change          =',dWtots,' kg m-2'
-               	print*,'---------------------------------------------------------------'
+                  print*,'In hydrol surface water change          =',dWtots,' kg m-2'
+                  print*,'---------------------------------------------------------------'
 	       endif
             endif
-            !-------------------------
-
-         ELSE     ! of if (hydrology)
+
+         else ! of if (hydrology)
 
             qsurf(1:ngrid,1:nq)=qsurf(1:ngrid,1:nq)+ptimestep*dqsurf(1:ngrid,1:nq)
 
-         END IF   ! of if (hydrology)
-
-         ! Add qsurf to qsurf_hist, which is what we save in
-         ! diagfi.nc etc. At the same time, we set the water 
-         ! content of ocean gridpoints back to zero, in order
-         ! to avoid rounding errors in vdifc, rain
+         end if ! of if (hydrology)
+
+         ! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water 
+         ! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain.
          qsurf_hist(:,:) = qsurf(:,:)
 
@@ -1589,12 +1499,13 @@
          endif
 
-      endif   !  of if (tracer) 
-
-!-----------------------------------------------------------------------
-!   9. Surface and sub-surface soil temperature
-!-----------------------------------------------------------------------
-
-!     Increment surface temperature
-
+      endif! end of if 'tracer'
+
+
+!------------------------------------------------
+!   VII. Surface and sub-surface soil temperature 
+!------------------------------------------------
+
+
+      ! Increment surface temperature
       if(ok_slab_ocean)then
          do ig=1,ngrid
@@ -1608,56 +1519,60 @@
       else
         tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid) 
-      endif!(ok_slab_ocean)
-
-!     Compute soil temperatures and subsurface heat flux
+      endif ! end of 'ok_slab_ocean'
+
+
+      ! Compute soil temperatures and subsurface heat flux.
       if (callsoil) then
          call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat,   &
-                ptimestep,tsurf,tsoil,capcal,fluxgrd)     
+                   ptimestep,tsurf,tsoil,capcal,fluxgrd)     
       endif
 
 
       if (ok_slab_ocean) then
-            do ig=1,ngrid
-               fluxgrdocean(ig)=fluxgrd(ig)
-               if (nint(rnat(ig)).eq.0) then
+      
+         do ig=1,ngrid
+         
+            fluxgrdocean(ig)=fluxgrd(ig)
+            if (nint(rnat(ig)).eq.0) then
                capcal(ig)=capcalocean
                fluxgrd(ig)=0.
                fluxgrdocean(ig)=pctsrf_sic(ig)*flux_g(ig)+(1-pctsrf_sic(ig))*(dt_hdiff(ig,1)+dt_ekman(ig,1))
-                 do iq=1,nsoilmx
-                    tsoil(ig,iq)=tsurf(ig)
-                 enddo
-                 if (pctsrf_sic(ig).gt.0.5) then
-                   capcal(ig)=capcalseaice
-                   if (qsurf(ig,igcm_h2o_ice).gt.0.) then
-                   capcal(ig)=capcalsno
-                   endif
-                 endif
-               endif
-            enddo
-      endif !(ok_slab_ocean)
-
-!-------------------------
-! test energy conservation
+               do iq=1,nsoilmx
+                  tsoil(ig,iq)=tsurf(ig)
+               enddo
+               if (pctsrf_sic(ig).gt.0.5) then
+                  capcal(ig)=capcalseaice
+                  if (qsurf(ig,igcm_h2o_ice).gt.0.) then
+                     capcal(ig)=capcalsno
+                  endif
+               endif               
+            endif
+            
+         enddo
+         
+      endif !end of 'ok_slab_ocean'
+
+
+      ! Test energy conservation
       if(enertest)then
 	 call planetwide_sumval(area(:)*capcal(:)*zdtsurf(:)/totarea_planet,dEtots)	 
 	 if (is_master) print*,'Surface energy change                 =',dEtots,' W m-2'
       endif
-!-------------------------
-
-!-----------------------------------------------------------------------
-!  10. Perform diagnostics and write output files
-!-----------------------------------------------------------------------
-
-!    -------------------------------
-!    Dynamical fields incrementation
-!    -------------------------------
-!    For output only: the actual model integration is performed in the dynamics
+
+
+!---------------------------------------------------
+!   VIII. Perform diagnostics and write output files
+!---------------------------------------------------
+
+      ! Note : For output only: the actual model integration is performed in the dynamics.
+
+
  
-      ! temperature, zonal and meridional wind
+      ! Temperature, zonal and meridional winds.
       zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
       zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
       zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
 
-      ! diagnostic
+      ! Diagnostic.
       zdtdyn(1:ngrid,1:nlayer)     = pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)
       ztprevious(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
@@ -1667,25 +1582,18 @@
       endif
 
-      ! dynamical heating diagnostic
+      ! Dynamical heating diagnostic.
       do ig=1,ngrid
          fluxdyn(ig)= SUM(zdtdyn(ig,:) *mass(ig,:))*cpp/ptimestep
       enddo
 
-      ! tracers
+      ! Tracers.
       zq(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq) + pdq(1:ngrid,1:nlayer,1:nq)*ptimestep
 
-      ! surface pressure
+      ! Surface pressure.
       ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
 
-      ! pressure
-      do l=1,nlayer
-          zplev(1:ngrid,l) = pplev(1:ngrid,l)/pplev(1:ngrid,1)*ps(:)
-          zplay(1:ngrid,l) = pplay(1:ngrid,l)/pplev(1:ngrid,1)*ps(1:ngrid)
-      enddo
-
-!     ---------------------------------------------------------
-!     Surface and soil temperature information
-!     ---------------------------------------------------------
-
+
+
+      ! Surface and soil temperature information
       call planetwide_sumval(area(:)*tsurf(:)/totarea_planet,Ts1)
       call planetwide_minval(tsurf(:),Ts2)
@@ -1693,17 +1601,15 @@
       if(callsoil)then
          TsS = SUM(area(:)*tsoil(:,nsoilmx))/totarea        ! mean temperature at bottom soil layer
-           print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]     ave[Tdeep]'
-           print*,Ts1,Ts2,Ts3,TsS
+         print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]     ave[Tdeep]'
+         print*,Ts1,Ts2,Ts3,TsS
       else
-      	if (is_master) then
-           print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]'
-           print*,Ts1,Ts2,Ts3
-	endif
+      	 if (is_master) then
+            print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]'
+            print*,Ts1,Ts2,Ts3
+	 endif
       end if
 
-!     ---------------------------------------------------------
-!     Check the energy balance of the simulation during the run
-!     ---------------------------------------------------------
-
+
+      ! Check the energy balance of the simulation during the run
       if(corrk)then
 
@@ -1730,6 +1636,6 @@
 	 
 	 if (is_master) then
-         	print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
-         	print*, ISR,ASR,OLR,GND,DYN
+            print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
+            print*, ISR,ASR,OLR,GND,DYN
 	 endif
 
@@ -1748,7 +1654,7 @@
                open(93,file="tem_bal.out",form='formatted',position='append')
                if(callsoil)then
-               	write(93,*) zday,Ts1,Ts2,Ts3,TsS
+               	  write(93,*) zday,Ts1,Ts2,Ts3,TsS
 	       else
-	        write(93,*) zday,Ts1,Ts2,Ts3
+	          write(93,*) zday,Ts1,Ts2,Ts3
 	       endif
                close(93)
@@ -1756,10 +1662,8 @@
          endif
 
-      endif
-
-
-!     ------------------------------------------------------------------
-!     Diagnostic to test radiative-convective timescales in code
-!     ------------------------------------------------------------------
+      endif ! end of 'corrk'
+
+
+      ! Diagnostic to test radiative-convective timescales in code.
       if(testradtimes)then
          open(38,file="tau_phys.out",form='formatted',position='append')
@@ -1774,16 +1678,15 @@
       endif
 
-!     ---------------------------------------------------------
-!     Compute column amounts (kg m-2) if tracers are enabled
-!     ---------------------------------------------------------
+
+      ! Compute column amounts (kg m-2) if tracers are enabled.
       if(tracer)then
          qcol(1:ngrid,1:nq)=0.0
          do iq=1,nq
-           do ig=1,ngrid
-              qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer))
-           enddo
+            do ig=1,ngrid
+               qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer))
+            enddo
          enddo
 
-         ! Generalised for arbitrary aerosols now. (LK)
+         ! Generalised for arbitrary aerosols now. By LK
          reffcol(1:ngrid,1:naerkind)=0.0
          if(co2cond.and.(iaero_co2.ne.0))then
@@ -1801,10 +1704,8 @@
          endif
 
-      endif
-
-!     ---------------------------------------------------------
-!     Test for water conservation if water is enabled
-!     ---------------------------------------------------------
-
+      endif ! end of 'tracer'
+
+
+      ! Test for water conservation.
       if(water)then
 
@@ -1817,7 +1718,7 @@
 	 
 	 if (is_master) then
-         	print*,' Total water amount [kg m^-2]: ',h2otot
-         	print*,' Surface ice    Surface liq.   Atmos. con.     Atmos. vap. [kg m^-2] '
-         	print*, icesrf,liqsrf,icecol,vapcol
+            print*,' Total water amount [kg m^-2]: ',h2otot
+            print*,' Surface ice    Surface liq.   Atmos. con.     Atmos. vap. [kg m^-2] '
+            print*, icesrf,liqsrf,icecol,vapcol
 	 endif
 
@@ -1833,9 +1734,8 @@
       endif
 
-!     ---------------------------------------------------------
-!     Calculate RH for diagnostic if water is enabled
-!     ---------------------------------------------------------
-
+
+      ! Calculate RH (Relative Humidity) for diagnostic.
       if(water)then
+      
          do l = 1, nlayer
             do ig=1,ngrid
@@ -1845,16 +1745,19 @@
          enddo
 
-         ! compute maximum possible H2O column amount (100% saturation)
+         ! Compute maximum possible H2O column amount (100% saturation).
          do ig=1,ngrid
-               H2Omaxcol(ig) = SUM( qsat(ig,:) * mass(ig,:))
+            H2Omaxcol(ig) = SUM( qsat(ig,:) * mass(ig,:))
          enddo
 
-      endif
-
-
-           print*,'--> Ls =',zls*180./pi
-!        -------------------------------------------------------------------
+      endif ! end of 'water'
+
+
+      print*,'--> Ls =',zls*180./pi
+      
+      
+!----------------------------------------------------------------------
 !        Writing NetCDF file  "RESTARTFI" at the end of the run
-!        -------------------------------------------------------------------
+!----------------------------------------------------------------------
+
 !        Note: 'restartfi' is stored just before dynamics are stored
 !              in 'restart'. Between now and the writting of 'restart',
@@ -1863,132 +1766,126 @@
 !              thus we store for time=time+dtvr
 
-         if(lastcall) then
-            ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
-
-
-!           Update surface ice distribution to iterate to steady state if requested
-            if(ice_update)then
-
-               do ig=1,ngrid
-
-                  delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
-
-                  ! add multiple years of evolution
-                  qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep 
-
-                  ! if ice has gone -ve, set to zero
-                  if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
-                     qsurf_hist(ig,igcm_h2o_ice) = 0.0 
-                  endif
-
-                  ! if ice is seasonal, set to zero (NEW)
-                  if(ice_min(ig).lt.0.01)then
-                     qsurf_hist(ig,igcm_h2o_ice) = 0.0 
-                  endif
-
-               enddo
-
-               ! enforce ice conservation
-               ice_tot= SUM(qsurf_hist(:,igcm_h2o_ice)*area(:) )
-               qsurf_hist(:,igcm_h2o_ice) = qsurf_hist(:,igcm_h2o_ice)*(icesrf/ice_tot)
-
-            endif
-
-            if (ngrid.ne.1) then
-              write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
-!#ifdef CPP_PARA
-!! for now in parallel we use a different routine to write restartfi.nc
-!            call phyredem(ngrid,"restartfi.nc",           &
-!                    ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
-!                    cloudfrac,totcloudfrac,hice)
-!#else
-!            call physdem1(ngrid,"restartfi.nc",long,lati,nsoilmx,nq,            &
-!                    ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
-!                    area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,      &
-!                    cloudfrac,totcloudfrac,hice,noms)
-!#endif
-              call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
-                      ptimestep,ztime_fin, &
-                      tsurf,tsoil,emis,q2,qsurf_hist, &
-                      cloudfrac,totcloudfrac,hice, &
-                      rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
-            endif
-
-            if(ok_slab_ocean) then
-              call ocean_slab_final!(tslab, seaice)
-            end if
+
+
+      if(lastcall) then
+         ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
+
+         ! Update surface ice distribution to iterate to steady state if requested
+         if(ice_update)then
+
+            do ig=1,ngrid
+
+               delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
+               
+               ! add multiple years of evolution
+               qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep 
+
+               ! if ice has gone -ve, set to zero
+               if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
+                  qsurf_hist(ig,igcm_h2o_ice) = 0.0 
+               endif
+
+               ! if ice is seasonal, set to zero (NEW)
+               if(ice_min(ig).lt.0.01)then
+                  qsurf_hist(ig,igcm_h2o_ice) = 0.0 
+               endif
+
+            enddo
+
+            ! enforce ice conservation
+            ice_tot= SUM(qsurf_hist(:,igcm_h2o_ice)*area(:) )
+            qsurf_hist(:,igcm_h2o_ice) = qsurf_hist(:,igcm_h2o_ice)*(icesrf/ice_tot)
+
+         endif
+
+         if (ngrid.ne.1) then
+            write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
+            
+            call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
+                          ptimestep,ztime_fin,                    &
+                          tsurf,tsoil,emis,q2,qsurf_hist,         &
+                          cloudfrac,totcloudfrac,hice,            &
+                          rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
+         endif
+
+         if(ok_slab_ocean) then
+            call ocean_slab_final!(tslab, seaice)
+         end if
 
          
-         endif ! of if (lastcall)
-
-
-!        -----------------------------------------------------------------
+      endif ! end of 'lastcall'
+
+
+!-----------------------------------
 !        Saving statistics :
-!        -----------------------------------------------------------------
-!        ("stats" stores and accumulates 8 key variables in file "stats.nc"
-!        which can later be used to make the statistic files of the run:
-!        "stats")          only possible in 3D runs !
+!-----------------------------------
+
+!    Note :("stats" stores and accumulates 8 key variables in file "stats.nc"
+!           which can later be used to make the statistic files of the run:
+!           "stats")          only possible in 3D runs !!!
 
          
-         if (callstats) then
-
-            call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
-            call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
-            call wstats(ngrid,"fluxsurf_lw",                               &
-                        "Thermal IR radiative flux to surface","W.m-2",2,  &
-                         fluxsurf_lw)
+      if (callstats) then
+
+         call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
+         call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
+         call wstats(ngrid,"fluxsurf_lw",                               &
+                     "Thermal IR radiative flux to surface","W.m-2",2,  &
+                     fluxsurf_lw)
+         call wstats(ngrid,"fluxtop_lw",                                &
+                     "Thermal IR radiative flux to space","W.m-2",2,    &
+                     fluxtop_lw)
+                     
 !            call wstats(ngrid,"fluxsurf_sw",                               &
 !                        "Solar radiative flux to surface","W.m-2",2,       &
-!                         fluxsurf_sw_tot)
-            call wstats(ngrid,"fluxtop_lw",                                &
-                        "Thermal IR radiative flux to space","W.m-2",2,    &
-                        fluxtop_lw)
+!                         fluxsurf_sw_tot)                     
 !            call wstats(ngrid,"fluxtop_sw",                                &
 !                        "Solar radiative flux to space","W.m-2",2,         &
 !                        fluxtop_sw_tot)
 
-            call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
-            call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
-            call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
-            call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo)
-            call wstats(ngrid,"p","Pressure","Pa",3,pplay)
-            call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
-            call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
-            call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
-            call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
-            call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
-
-           if (tracer) then
-             do iq=1,nq
-                call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
-                call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
-                          'kg m^-2',2,qsurf(1,iq) )
-
-                call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    & 
+
+         call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
+         call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
+         call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
+         call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo)
+         call wstats(ngrid,"p","Pressure","Pa",3,pplay)
+         call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
+         call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
+         call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
+         call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
+         call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
+
+         if (tracer) then
+            do iq=1,nq
+               call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
+               call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
+                           'kg m^-2',2,qsurf(1,iq) )
+               call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    & 
                           'kg m^-2',2,qcol(1,iq) )
-!                call wstats(ngrid,trim(noms(iq))//'_reff',                          & 
+                          
+!              call wstats(ngrid,trim(noms(iq))//'_reff',                          & 
 !                          trim(noms(iq))//'_reff',                                   & 
 !                          'm',3,reffrad(1,1,iq))
-              end do
+
+            end do
+            
             if (water) then
                vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
-               call wstats(ngrid,"vmr_h2ovapor",                           &
-                          "H2O vapour volume mixing ratio","mol/mol",       & 
-                          3,vmr)
-            endif ! of if (water)
-
-           endif !tracer 
-
-          if(watercond.and.CLFvarying)then
-             call wstats(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
-             call wstats(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
-             call wstats(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
-             call wstats(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
-             call wstats(ngrid,"RH","relative humidity"," ",3,RH)
-          endif
-
-
-
-           if (ok_slab_ocean) then
+               call wstats(ngrid,"vmr_h2ovapor",                             &
+                           "H2O vapour volume mixing ratio","mol/mol",       & 
+                           3,vmr)
+            endif
+
+         endif ! end of 'tracer' 
+
+         if(watercond.and.CLFvarying)then
+            call wstats(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
+            call wstats(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
+            call wstats(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
+            call wstats(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
+            call wstats(ngrid,"RH","relative humidity"," ",3,RH)
+         endif
+
+         if (ok_slab_ocean) then
             call wstats(ngrid,"dt_hdiff1","dt_hdiff1","K/s",2,dt_hdiff(:,1))
             call wstats(ngrid,"dt_hdiff2","dt_hdiff2","K/s",2,dt_hdiff(:,2))
@@ -2001,32 +1898,36 @@
             call wstats(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
             call wstats(ngrid,"rnat","nature of the surface","",2,rnat)
-           endif! (ok_slab_ocean)
-
-            if(lastcall) then
-              write (*,*) "Writing stats..."
-              call mkstats(ierr)
-            endif
-          endif !if callstats
-
-
-!        ----------------------------------------------------------------------
-!        output in netcdf file "DIAGFI", containing any variable for diagnostic
-!        (output with  period "ecritphy", set in "run.def")
-!        ----------------------------------------------------------------------
-!        writediagfi can also be called from any other subroutine for any variable.
-!        but its preferable to keep all the calls in one place...
-
-        call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
-	call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi)
-	call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
-	call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
-        call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
-        call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
-        call writediagfi(ngrid,"temp","temperature","K",3,zt)
-        call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
-        call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
-        call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
-        call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
-        call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
+         endif
+
+         if(lastcall) then
+            write (*,*) "Writing stats..."
+            call mkstats(ierr)
+         endif
+         
+      endif ! end of 'callstats'
+
+
+!-----------------------------------------------------------------------------------------------------
+!           OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic
+!
+!             Note 1 : output with  period "ecritphy", set in "run.def"
+!
+!             Note 2 : writediagfi can also be called from any other subroutine for any variable,
+!                      but its preferable to keep all the calls in one place ...
+!-----------------------------------------------------------------------------------------------------
+
+
+      call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
+      call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi)
+      call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
+      call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
+      call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
+      call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
+      call writediagfi(ngrid,"temp","temperature","K",3,zt)
+      call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
+      call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
+      call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
+      call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
+      call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
 
 !     Subsurface temperatures
@@ -2034,27 +1935,30 @@
 !        call writediagsoil(ngrid,"temp","temperature","K",3,tsoil)
 
-!     Oceanic layers
-        if(ok_slab_ocean) then
-          call writediagfi(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic)
-          call writediagfi(ngrid,"tsea_ice","tsea_ice","K",2,tsea_ice)
-          call writediagfi(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
-          call writediagfi(ngrid,"tslab1","tslab1","K",2,tslab(:,1))
-          call writediagfi(ngrid,"tslab2","tslab2","K",2,tslab(:,2))
-          call writediagfi(ngrid,"dt_hdiff1","dt_hdiff1","K/s",2,dt_hdiff(:,1))
-          call writediagfi(ngrid,"dt_hdiff2","dt_hdiff2","K/s",2,dt_hdiff(:,2))
-          call writediagfi(ngrid,"dt_ekman1","dt_ekman1","K/s",2,dt_ekman(:,1))
-          call writediagfi(ngrid,"dt_ekman2","dt_ekman2","K/s",2,dt_ekman(:,2))
-          call writediagfi(ngrid,"rnat","nature of the surface","",2,rnat)
-          call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
-          call writediagfi(ngrid,"latentFlux","latent heat flux","w.m^-2",2,zdqsdif(:,igcm_h2o_vap)*RLVTT)
-       endif! (ok_slab_ocean)
-
-!     Total energy balance diagnostics
-        if(callrad.and.(.not.newtonian))then
-           call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo)
-           call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
-           call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
-           call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
-           call writediagfi(ngrid,"shad","rings"," ", 2, fract)
+      ! Oceanic layers
+      if(ok_slab_ocean) then
+         call writediagfi(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic)
+         call writediagfi(ngrid,"tsea_ice","tsea_ice","K",2,tsea_ice)
+         call writediagfi(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
+         call writediagfi(ngrid,"tslab1","tslab1","K",2,tslab(:,1))
+         call writediagfi(ngrid,"tslab2","tslab2","K",2,tslab(:,2))
+         call writediagfi(ngrid,"dt_hdiff1","dt_hdiff1","K/s",2,dt_hdiff(:,1))
+         call writediagfi(ngrid,"dt_hdiff2","dt_hdiff2","K/s",2,dt_hdiff(:,2))
+         call writediagfi(ngrid,"dt_ekman1","dt_ekman1","K/s",2,dt_ekman(:,1))
+         call writediagfi(ngrid,"dt_ekman2","dt_ekman2","K/s",2,dt_ekman(:,2))
+         call writediagfi(ngrid,"rnat","nature of the surface","",2,rnat)
+         call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
+         call writediagfi(ngrid,"latentFlux","latent heat flux","w.m^-2",2,zdqsdif(:,igcm_h2o_vap)*RLVTT)
+      endif
+      
+
+      ! Total energy balance diagnostics
+      if(callrad.and.(.not.newtonian))then
+      
+         call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo)
+         call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
+         call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
+         call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
+         call writediagfi(ngrid,"shad","rings"," ", 2, fract)
+         
 !           call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
 !           call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
@@ -2063,100 +1967,119 @@
 !           call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1)
 !           call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
-           if(ok_slab_ocean) then
-             call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrdocean)
-           else
-             call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
-           endif!(ok_slab_ocean)
-           call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
-	 endif
+
+         if(ok_slab_ocean) then
+            call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrdocean)
+         else
+            call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
+         endif
+         
+         call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
+         
+      endif ! end of 'callrad'
 	
-	if(enertest) then
-	  if (calldifv) then
-	     call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
+      if(enertest) then
+      
+	 if (calldifv) then
+         
+	    call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
+            call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
+            
 !	     call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
 !	     call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
 !	     call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
-	     call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
-	  endif
-	  if (corrk) then
-	     call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
-	     call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
-	  endif 
-          if(watercond) then
+	     
+	 endif
+          
+	 if (corrk) then
+	    call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
+	    call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
+	 endif
+          
+         if(watercond) then
+
+	    call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE) 
+	    call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
+	    call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat)
+             
 !	     call writediagfi(ngrid,"lscaledEz","heat from largescale","W m-2",3,lscaledEz) 
-!	     call writediagfi(ngrid,"madjdEz","heat from moistadj","W m-2",3,madjdEz)
-	     call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE) 
-	     call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
-	     call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat)
+!	     call writediagfi(ngrid,"madjdEz","heat from moistadj","W m-2",3,madjdEz)             
 !	     call writediagfi(ngrid,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
-          endif 
-        endif
-
-!     Temporary inclusions for heating diagnostics
-!        call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
-        call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
-        call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
-        call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
-
-        ! debugging
+
+         endif
+          
+      endif ! end of 'enertest'
+
+
+        ! Temporary inclusions for heating diagnostics.
+        !call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
+        !call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
+        !call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
+        ! call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
+        
+        ! For Debugging.
         !call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
         !call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
-
-!     Output aerosols
-        if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
-	  call writediagfi(ngrid,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
-        if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
-	  call writediagfi(ngrid,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
-        if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
-	  call writediagfi(ngrid,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
-        if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
-	  call writediagfi(ngrid,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
-
-!     Output tracers
-       if (tracer) then
-        do iq=1,nq
-          call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
+        
+
+      ! Output aerosols.
+      if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
+	 call writediagfi(ngrid,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
+      if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
+	 call writediagfi(ngrid,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
+      if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
+	 call writediagfi(ngrid,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
+      if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
+	 call writediagfi(ngrid,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
+
+      ! Output tracers.
+      if (tracer) then
+      
+         do iq=1,nq
+            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))//'_surf',trim(noms(iq))//'_surf',  & 
-!                          'kg m^-2',2,qsurf(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) )
-
-          if(watercond.or.CLFvarying)then
-             call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
-             call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
-             call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
-             call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
-	     call writediagfi(ngrid,"RH","relative humidity"," ",3,RH)
-          endif
-
-          if(waterrain)then
-             call writediagfi(ngrid,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
-             call writediagfi(ngrid,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
-          endif
-
-          if((hydrology).and.(.not.ok_slab_ocean))then
-             call writediagfi(ngrid,"hice","oceanic ice height","m",2,hice)
-          endif
-
-          if(ice_update)then
-             call writediagfi(ngrid,"ice_min","min annual ice","m",2,ice_min)
-             call writediagfi(ngrid,"ice_ini","initial annual ice","m",2,ice_initial)
-          endif
-
-          do ig=1,ngrid
-             if(tau_col(ig).gt.1.e3)then
-!                print*,'WARNING: tau_col=',tau_col(ig)
-!                print*,'at ig=',ig,'in PHYSIQ'
-             endif
-          end do
-
-          call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col)
-
-        enddo
-       endif
-
-!      output spectrum 
+!                          'kg m^-2',2,qsurf(1,iq) )                          
+
+            if(watercond.or.CLFvarying)then
+               call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
+               call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
+               call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
+               call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
+	       call writediagfi(ngrid,"RH","relative humidity"," ",3,RH)
+            endif
+
+            if(waterrain)then
+               call writediagfi(ngrid,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
+               call writediagfi(ngrid,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
+            endif
+
+            if((hydrology).and.(.not.ok_slab_ocean))then
+               call writediagfi(ngrid,"hice","oceanic ice height","m",2,hice)
+            endif
+
+            if(ice_update)then
+               call writediagfi(ngrid,"ice_min","min annual ice","m",2,ice_min)
+               call writediagfi(ngrid,"ice_ini","initial annual ice","m",2,ice_initial)
+            endif
+
+            do ig=1,ngrid
+               if(tau_col(ig).gt.1.e3)then
+                  print*,'WARNING: tau_col=',tau_col(ig)
+                  print*,'at ig=',ig,'in PHYSIQ'
+               endif         
+            end do
+
+            call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col)
+
+         enddo ! end of 'nq' loop
+         
+       endif ! end of 'tracer'
+
+
+      ! Output spectrum.
       if(specOLR.and.corrk)then      
          call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
@@ -2166,100 +2089,3 @@
       icount=icount+1
 
-      if (lastcall) then
-
-        ! deallocate gas variables
-!$OMP BARRIER
-!$OMP MASTER
-        IF ( ALLOCATED( gnom ) ) DEALLOCATE( gnom )
-        IF ( ALLOCATED( gfrac ) ) DEALLOCATE( gfrac )  ! both allocated in su_gases.F90
-!$OMP END MASTER
-!$OMP BARRIER
-
-        ! deallocate saved arrays
-        ! this is probably not that necessary
-        ! ... but probably a good idea to clean the place before we leave
-        IF ( ALLOCATED(tsurf)) DEALLOCATE(tsurf)
-        IF ( ALLOCATED(tsoil)) DEALLOCATE(tsoil)
-        IF ( ALLOCATED(albedo)) DEALLOCATE(albedo)
-        IF ( ALLOCATED(albedo0)) DEALLOCATE(albedo0)
-        IF ( ALLOCATED(rnat)) DEALLOCATE(rnat)
-        IF ( ALLOCATED(emis)) DEALLOCATE(emis)
-        IF ( ALLOCATED(dtrad)) DEALLOCATE(dtrad)
-        IF ( ALLOCATED(fluxrad_sky)) DEALLOCATE(fluxrad_sky)
-        IF ( ALLOCATED(fluxrad)) DEALLOCATE(fluxrad)
-        IF ( ALLOCATED(capcal)) DEALLOCATE(capcal)
-        IF ( ALLOCATED(fluxgrd)) DEALLOCATE(fluxgrd)
-        IF ( ALLOCATED(qsurf)) DEALLOCATE(qsurf)
-        IF ( ALLOCATED(q2)) DEALLOCATE(q2)
-        IF ( ALLOCATED(ztprevious)) DEALLOCATE(ztprevious)
-	IF ( ALLOCATED(hice)) DEALLOCATE(hice)
-        IF ( ALLOCATED(cloudfrac)) DEALLOCATE(cloudfrac)
-        IF ( ALLOCATED(totcloudfrac)) DEALLOCATE(totcloudfrac)
-        IF ( ALLOCATED(qsurf_hist)) DEALLOCATE(qsurf_hist)
-        IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad)
-        IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad)
-        IF ( ALLOCATED(ice_initial)) DEALLOCATE(ice_initial)
-        IF ( ALLOCATED(ice_min)) DEALLOCATE(ice_min)
-
-        IF ( ALLOCATED(fluxsurf_lw)) DEALLOCATE(fluxsurf_lw)
-        IF ( ALLOCATED(fluxsurf_sw)) DEALLOCATE(fluxsurf_sw)
-        IF ( ALLOCATED(fluxtop_lw)) DEALLOCATE(fluxtop_lw)
-        IF ( ALLOCATED(fluxabs_sw)) DEALLOCATE(fluxabs_sw)
-        IF ( ALLOCATED(fluxtop_dn)) DEALLOCATE(fluxtop_dn)
-        IF ( ALLOCATED(fluxdyn)) DEALLOCATE(fluxdyn)
-        IF ( ALLOCATED(OLR_nu)) DEALLOCATE(OLR_nu)
-        IF ( ALLOCATED(OSR_nu)) DEALLOCATE(OSR_nu)
-        IF ( ALLOCATED(sensibFlux)) DEALLOCATE(sensibFlux)
-        IF ( ALLOCATED(zdtlw)) DEALLOCATE(zdtlw)
-        IF ( ALLOCATED(zdtsw)) DEALLOCATE(zdtsw)
-        IF ( ALLOCATED(tau_col)) DEALLOCATE(tau_col)
-	IF ( ALLOCATED(pctsrf_sic)) DEALLOCATE(pctsrf_sic)
-        IF ( ALLOCATED(tslab)) DEALLOCATE(tslab)
-        IF ( ALLOCATED(tsea_ice)) DEALLOCATE(tsea_ice)
-        IF ( ALLOCATED(sea_ice)) DEALLOCATE(sea_ice)
-        IF ( ALLOCATED(zmasq)) DEALLOCATE(zmasq)
-        IF ( ALLOCATED(knindex)) DEALLOCATE(knindex)
-
-        !! this is defined in comsaison_h
-        IF ( ALLOCATED(mu0)) DEALLOCATE(mu0)
-        IF ( ALLOCATED(fract)) DEALLOCATE(fract)
-
-        !! this is defined in comsoil_h
-        IF ( ALLOCATED(layer)) DEALLOCATE(layer)
-        IF ( ALLOCATED(mlayer)) DEALLOCATE(mlayer)
-        IF ( ALLOCATED(inertiedat)) DEALLOCATE(inertiedat)
-
-        !! this is defined in surfdat_h
-        IF ( ALLOCATED(albedodat)) DEALLOCATE(albedodat)
-        IF ( ALLOCATED(phisfi)) DEALLOCATE(phisfi)
-        IF ( ALLOCATED(zmea)) DEALLOCATE(zmea)
-        IF ( ALLOCATED(zstd)) DEALLOCATE(zstd)
-        IF ( ALLOCATED(zsig)) DEALLOCATE(zsig)
-        IF ( ALLOCATED(zgam)) DEALLOCATE(zgam)
-        IF ( ALLOCATED(zthe)) DEALLOCATE(zthe)
-        IF ( ALLOCATED(dryness)) DEALLOCATE(dryness)
-        IF ( ALLOCATED(watercaptag)) DEALLOCATE(watercaptag)
- 
-        !! this is defined in tracer_h
-        IF ( ALLOCATED(noms)) DEALLOCATE(noms)
-        IF ( ALLOCATED(mmol)) DEALLOCATE(mmol)
-        IF ( ALLOCATED(radius)) DEALLOCATE(radius)
-        IF ( ALLOCATED(rho_q)) DEALLOCATE(rho_q)
-        IF ( ALLOCATED(qext)) DEALLOCATE(qext)
-        IF ( ALLOCATED(alpha_lift)) DEALLOCATE(alpha_lift)
-        IF ( ALLOCATED(alpha_devil)) DEALLOCATE(alpha_devil)
-        IF ( ALLOCATED(qextrhor)) DEALLOCATE(qextrhor)
-        IF ( ALLOCATED(igcm_dustbin)) DEALLOCATE(igcm_dustbin)
-
-        !! this is defined in comgeomfi_h
-        IF ( ALLOCATED(lati)) DEALLOCATE(lati)
-        IF ( ALLOCATED(long)) DEALLOCATE(long)
-        IF ( ALLOCATED(area)) DEALLOCATE(area)
-        IF ( ALLOCATED(sinlat)) DEALLOCATE(sinlat)
-        IF ( ALLOCATED(coslat)) DEALLOCATE(coslat)
-        IF ( ALLOCATED(sinlon)) DEALLOCATE(sinlon)
-        IF ( ALLOCATED(coslon)) DEALLOCATE(coslon)
-      endif
-
-      return
     end subroutine physiq
Index: trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90	(revision 1470)
+++ trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90	(revision 1477)
@@ -3,5 +3,5 @@
           pplay,pplev,pzlay,pzlev,pz0,                 &
           pu,pv,pt,ppopsk,pq,ptsrf,pemis,pqsurf,       &
-          pdufi,pdvfi,pdtfi,pdqfi,pfluxsrf,            &
+          pdtfi,pdqfi,pfluxsrf,            &
           Pdudif,pdvdif,pdtdif,pdtsrf,sensibFlux,pq2,  &
           pdqdif,pdqevap,pdqsdif,flux_u,flux_v,lastcall)
@@ -54,5 +54,4 @@
       REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
       REAL,INTENT(IN) :: pemis(ngrid)
-      REAL,INTENT(IN) :: pdufi(ngrid,nlay),pdvfi(ngrid,nlay)
       REAL,INTENT(IN) :: pdtfi(ngrid,nlay)
       REAL,INTENT(IN) :: pfluxsrf(ngrid)
@@ -189,6 +188,6 @@
       DO ilev=1,nlay
          DO ig=1,ngrid
-            zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
-            zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
+            zu(ig,ilev)=pu(ig,ilev)
+            zv(ig,ilev)=pv(ig,ilev)
             zt(ig,ilev)=pt(ig,ilev)+pdtfi(ig,ilev)*ptimestep
             zh(ig,ilev)=pt(ig,ilev)*zovExner(ig,ilev) !for call vdif_kc, but could be moved and computed there
@@ -689,6 +688,6 @@
       do ilev = 1, nlay
          do ig=1,ngrid
-            pdudif(ig,ilev)=(zu(ig,ilev)-(pu(ig,ilev)+pdufi(ig,ilev)*ptimestep))/ptimestep
-            pdvdif(ig,ilev)=(zv(ig,ilev)-(pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep))/ptimestep
+            pdudif(ig,ilev)=(zu(ig,ilev)-(pu(ig,ilev)))/ptimestep
+            pdvdif(ig,ilev)=(zv(ig,ilev)-(pv(ig,ilev)))/ptimestep
             pdtdif(ig,ilev)=( zt(ig,ilev)- pt(ig,ilev))/ptimestep-pdtfi(ig,ilev)
          enddo
Index: trunk/LMDZ.GENERIC/libf/phystd/vdifc.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/vdifc.F	(revision 1470)
+++ trunk/LMDZ.GENERIC/libf/phystd/vdifc.F	(revision 1477)
@@ -3,5 +3,5 @@
      &     pplay,pplev,pzlay,pzlev,pz0,
      &     pu,pv,ph,pq,ptsrf,pemis,pqsurf,
-     &     pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
+     &     pdhfi,pdqfi,pfluxsrf,
      &     pdudif,pdvdif,pdhdif,pdtsrf,sensibFlux,pq2,
      &     pdqdif,pdqsdif,lastcall)
@@ -49,5 +49,5 @@
       REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
       REAL ptsrf(ngrid),pemis(ngrid)
-      REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
+      REAL pdhfi(ngrid,nlay)
       REAL pfluxsrf(ngrid)
       REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
@@ -205,6 +205,6 @@
       DO ilev=1,nlay
          DO ig=1,ngrid
-            zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
-            zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
+            zu(ig,ilev)=pu(ig,ilev)
+            zv(ig,ilev)=pv(ig,ilev)
             zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
          ENDDO
@@ -662,7 +662,7 @@
          do ig=1,ngrid
             pdudif(ig,ilev)=(zu(ig,ilev)-
-     &           (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep))/ptimestep
+     &           (pu(ig,ilev)))/ptimestep
             pdvdif(ig,ilev)=(zv(ig,ilev)-
-     &           (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep))/ptimestep
+     &           (pv(ig,ilev)))/ptimestep
             hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep 
 
