subroutine physiq(ngrid,nlayer,nq,
     *            firstcall,lastcall,
     *            pday,ptime,ptimestep,
     *            pplev,pplay,pphi,
     *            pu,pv,pt,pq,
     *            pw,
     *            pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)

      use radinc_h, only : naerkind
      use planet_h
      USE ioipsl_getincom
      use comgeomfi_h
      use comsoil_h
      use aerosol_mod
      implicit none

!==================================================================
!     Purpose
!     -------
!     Central subroutine for all the physics parameterisations in the
!     universal model. Originally adapted from the Mars LMDZ model.
!
!     The model can be run without or with tracer transport
!     depending on the value of "tracer" in file "callphys.def".
!
!   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 :
!       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.
!
!   arguments
!   ---------
!
!   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
!    pudyn(ngrid,nlayer)    \
!    pvdyn(ngrid,nlayer)     \ Dynamical temporal derivative for the
!    ptdyn(ngrid,nlayer)     / corresponding variables
!    pqdyn(ngrid,nlayer,nq) /
!    pw(ngrid,?)           vertical velocity
!
!   output
!   ------
!
!    pdu(ngrid,nlayermx)        \
!    pdv(ngrid,nlayermx)         \  Temporal derivative of the corresponding
!    pdt(ngrid,nlayermx)         /  variables due to physical processes.
!    pdq(ngrid,nlayermx)        /
!    pdpsrf(ngrid)             /
!    tracerdyn                 call tracer in dynamical part of GCM ?
!
!     Authors
!     -------
!           Frederic Hourdin	15/10/93
!           Francois Forget	1994
!           Christophe Hourdin	02/1997
!           Subroutine completly rewritten by F.Forget (01/2000)
!           Water ice clouds: Franck Montmessin (update 06/2003)
!           Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
!           New correlated-k radiative scheme: R. Wordsworth (2009)
!           Many specifically Martian subroutines removed: R. Wordsworth (2009)
!           Pluto version improved :T. Bertrand (2015,2020)
!==================================================================

c    0.  Declarations :
c    ------------------
#include "dimensions.h"
#include "dimphys.h"
#include "surfdat.h"
#include "comdiurn.h"
#include "callkeys.h"
#include "comcstfi.h"
#include "comsaison.h"
#include "control.h"
#include "comg1d.h"
#include "tracer.h"
#include "netcdf.inc"

c Arguments :
c -----------

c   inputs:
c   -------
      INTEGER ngrid,nlayer,nq
      REAL ptimestep
      REAL pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
      REAL pphi(ngridmx,nlayer)
      REAL pu(ngridmx,nlayer),pv(ngridmx,nlayer)
      REAL pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
      REAL pw(ngridmx,nlayer)        ! pvervel transmitted by dyn3d
      REAL zh(ngridmx,nlayermx)      ! potential temperature (K)
      LOGICAL firstcall,lastcall
      REAL pday
      REAL ptime
      logical tracerdyn

c   outputs:
c   --------
c     physical tendencies
      REAL pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
      REAL pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
      REAL pdpsrf(ngridmx) ! surface pressure tendency
      REAL picen2(ngrid)

c Local saved variables:
c ----------------------
      INTEGER day_ini  ! Initial date of the run (sol since Ls=0)
      INTEGER icount     ! counter of calls to physiq during the run.
      REAL tsurf(ngridmx)            ! Surface temperature (K)
      REAL tsoil(ngridmx,nsoilmx)    ! sub-surface temperatures (K)
      REAL,SAVE :: tidat(ngridmx,nsoilmx)    ! thermal inertia soil
      REAL tidat_out(ngridmx,nlayermx)  ! thermal inertia soil but output on vertical grid
      REAL tsub(ngridmx)             ! temperature of the deepest layer
      REAL albedo(ngridmx)           ! Surface albedo

      REAL emis(ngridmx)             ! Thermal IR surface emissivity
      REAL,SAVE :: dtrad(ngridmx,nlayermx)   ! Net atm. radiative heating rate (K.s-1)
      REAL,SAVE :: fluxrad_sky(ngridmx)      ! rad. flux from sky absorbed by surface (W.m-2)
      REAL,SAVE :: fluxrad(ngridmx)          ! Net radiative surface flux (W.m-2)
      REAL,SAVE :: capcal(ngridmx)           ! surface heat capacity (J m-2 K-1)
      REAL dplanck(ngridmx)          ! for output (W.s/m2/K)
      REAL,SAVE :: fluxgrd(ngridmx)          ! surface conduction flux (W.m-2)
      REAL,SAVE :: qsurf(ngridmx,nqmx)       ! tracer on surface (e.g. kg.m-2)
      REAL q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy
      REAL qsurf1(ngridmx,nqmx)      ! saving qsurf to calculate flux over long timescales kg.m-2
      REAL flusurf(ngridmx,nqmx)     ! flux cond/sub kg.m-2.s-1
      REAL,SAVE :: flusurfold(ngridmx,nqmx)  ! old flux cond/sub kg.m-2.s-1
      REAL totmass                   ! total mass of n2 (atm + surf)
      REAL globave                   ! globalaverage 2D ps
      REAL globaveice(nqmx)          ! globalaverage 2D ice
      REAL globavenewice(nqmx)          ! globalaverage 2D ice

      REAL,SAVE :: ptime0    ! store the first time

      REAL dstep
      REAL,SAVE :: glastep=20   ! step en pluto day pour etaler le glacier

      SAVE day_ini,icount
      SAVE tsurf,tsoil,tsub
      SAVE albedo,emis,q2


      REAL stephan
      DATA stephan/5.67e-08/  ! Stephan Boltzman constant
      SAVE stephan

      REAL acond,bcond
      SAVE acond,bcond
      REAL tcond1p4Pa
      DATA tcond1p4Pa/38/

c Local variables :
c -----------------
!     Tendencies for the paleoclimate mode
      REAL,SAVE :: qsurfyear(ngridmx,nqmx)  ! kg.m-2 averaged mass of ice lost/gained in the last Pluto year of the run
      REAL,SAVE :: phisfinew(ngridmx)       ! geopotential of the bedrock (= phisfi-qsurf/1000*g)
      REAL qsurfpal(ngridmx,nqmx)           ! qsurf after a paleoclimate step : for physdem1 and restartfi
      REAL phisfipal(ngridmx)               ! geopotential after a paleoclimate step : for physdem1 and restartfi
      REAL oblipal                   ! change of obliquity
      REAL peri_daypal               ! new periday
      REAL eccpal                    ! change of eccentricity
      REAL tpalnew                   ! change of time
      REAL adjustnew                 ! change in N2 ice albedo
      REAL pdaypal                   ! new pday = day_ini + step
      REAL zdt_tot                   ! time range corresponding to the flux of qsurfyear
      REAL massacc(nqmx)             ! accumulated mass
      REAL masslost(nqmx)            ! accumulated mass

!     aerosol (ice or haze) extinction optical depth  at reference wavelength
!     for the "naerkind" optically active aerosols:
      REAL aerosol(ngridmx,nlayermx,naerkind)

      CHARACTER*80 fichier
      INTEGER l,ig,ierr,igout,iq,i, tapphys
      INTEGER lecttsoil     ! lecture of tsoil from proftsoil
!      INTEGER iqmin                     ! Used if iceparty engaged

      REAL fluxsurf_lw(ngridmx)      ! incident LW (IR) surface flux (W.m-2)
      REAL fluxsurf_sw(ngridmx)      ! incident SW (solar) surface flux (W.m-2)
      REAL,SAVE :: fluxtop_lw(ngridmx)       ! Outgoing LW (IR) flux to space (W.m-2)
      REAL,SAVE :: fluxtop_sw(ngridmx)       ! Outgoing SW (solar) flux to space (W.m-2)

!     included by RW for equilibration test
      real,save :: fluxtop_dn(ngridmx)
      real fluxabs_sw(ngridmx) ! absorbed shortwave radiation
      real fluxdyn(ngridmx)    ! horizontal heat transport by dynamics

      REAL zls                       ! solar longitude (rad)
      REAL zfluxuv                     ! Lyman alpha flux at 1AU
                                     !  ph/cm2/s
      REAL zday                      ! date (time since Ls=0, in martian days)
      REAL saveday                      ! saved date
      SAVE saveday
      REAL savedeclin                   ! saved declin
      SAVE savedeclin
      !REAL adjust                      ! adjustment for surface albedo
      REAL zzlay(ngridmx,nlayermx)   ! altitude at the middle of the layers
      REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries

c     Aerosol effective radius used for radiative transfer (units=meters)
      REAL reffrad(ngridmx,nlayermx,naerkind)

c     Tendencies due to various processes:
      REAL dqsurf(ngridmx,nqmx)
      REAL zdtlw(ngridmx,nlayermx)      ! (K/s)
      REAL zdtsw(ngridmx,nlayermx)      ! (K/s)
      REAL cldtlw(ngridmx,nlayermx)     ! (K/s) LW heating rate for clear areas
      REAL cldtsw(ngridmx,nlayermx)     ! (K/s) SW heating rate for clear areas
      REAL zdtsurf(ngridmx)             ! (K/s)
      REAL zdtlscale(ngridmx,nlayermx)
      REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx)  ! (m.s-2)
      REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx)         ! (K/s)
      REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx)  ! (m.s-2)
      REAL zdhadj(ngridmx,nlayermx)                           ! (K/s)
      REAL zdtgw(ngridmx,nlayermx)                            ! (K/s)
      REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx)    ! (m.s-2)
      REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
      REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
      REAL tconds(ngridmx,nlayermx)

      REAL zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
      REAL zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
      REAL zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
      REAL zdqadj(ngridmx,nlayermx,nqmx)
      REAL zdqc(ngridmx,nlayermx,nqmx)
      REAL zdqlscale(ngridmx,nlayermx,nqmx)
      REAL zdqslscale(ngridmx,nqmx), zdqsc(ngridmx,nqmx)
      REAL zdqchim(ngridmx,nlayermx,nqmx)
      REAL zdqschim(ngridmx,nqmx)
      REAL zdqch4cloud(ngridmx,nlayermx,nqmx)
      REAL zdqsch4cloud(ngridmx,nqmx)
      REAL zdtch4cloud(ngridmx,nlayermx)
      REAL zdqcocloud(ngridmx,nlayermx,nqmx)
      REAL zdqscocloud(ngridmx,nqmx)
      REAL zdtcocloud(ngridmx,nlayermx)
      REAL rice_ch4(ngridmx,nlayermx)    ! Methane ice geometric mean radius (m)
      REAL rice_co(ngridmx,nlayermx)     ! CO ice geometric mean radius (m)

      REAL zdqsch4fast(ngridmx)    ! used only for fast model nogcm
      REAL zdqch4fast(ngridmx)    ! used only for fast model nogcm
      REAL zdqscofast(ngridmx)    ! used only for fast model nogcm
      REAL zdqcofast(ngridmx)    ! used only for fast model nogcm
      REAL zdqflow(ngridmx,nqmx)

      REAL zdteuv(ngridmx,nlayermx)    ! (K/s)
      REAL zdtconduc(ngridmx,nlayermx) ! (K/s)
      REAL zdumolvis(ngridmx,nlayermx)
      REAL zdvmolvis(ngridmx,nlayermx)
      real zdqmoldiff(ngridmx,nlayermx,nqmx)

      ! Haze relatated tendancies
      REAL zdqhaze(ngridmx,nlayermx,nqmx)
      REAL zdqprodhaze(ngridmx,nqmx)
      REAL zdqsprodhaze(ngridmx)
      REAL zdqhaze_col(ngridmx)
      REAL zdqphot_prec(ngrid,nlayer)
      REAL zdqphot_ch4(ngrid,nlayer)
      REAL zdqconv_prec(ngrid,nlayer)
      REAL zdq_source(ngridmx,nlayermx,nqmx)
      ! Fast Haze relatated tendancies
      REAL fluxbot(ngridmx)
      REAL gradflux(ngridmx)
      REAL fluxlym_sol_bot(ngridmx)      ! Solar flux Lyman alpha ph.m-2.s-1 reaching the surface
      REAL fluxlym_ipm_bot(ngridmx)      ! IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 reaching the surface
      REAL flym_sol(ngridmx)      ! Incident Solar flux Lyman alpha ph.m-2.s-1
      REAL flym_ipm(ngridmx)      ! Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1

      real nconsMAX
      real vdifcncons(ngridmx)
      real cadjncons(ngridmx)
      real ch4csMAX
      real ch4cncons(ngridmx)
      real ch4sedncons(ngridmx)
      real ch4sedsMAX
      real condncons(ngridmx)
      real dWtot, dWtots
      real dWtotn2, dWtotsn2
      real nconsMAXn2
      real vdifcnconsn2(ngridmx)
      real cadjnconsn2(ngridmx)
      real condnconsn2(ngridmx)

c     Local variable for local intermediate calcul:
      REAL zflubid(ngridmx)
      REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
      REAL zdum1(ngridmx,nlayermx)
      REAL zdum2(ngridmx,nlayermx)
      REAL ztim1,ztim2,ztim3, z1,z2
      REAL ztime_fin
      REAL zdh(ngridmx,nlayermx)
      INTEGER length
      PARAMETER (length=100)

c local variables only used for diagnostic (output in file "diagfi" or "stats")
c -----------------------------------------------------------------------------
      REAL ps(ngridmx), zt(ngridmx,nlayermx)
      REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
      REAL zq(ngridmx,nlayermx,nqmx)
      REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx)
      character*2 str2
      character*5 str5
      real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
      real zdtdyn(ngridmx,nlayermx),ztprevious(ngridmx,nlayermx)
      save ztprevious
      real qtot1,qtot2            ! total aerosol mass
      integer igmin, lmin
      logical tdiag

      REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
      REAL cd
      real vmr(ngridmx,nlayermx) ! volume mixing ratio

      REAL time_phys

      REAL tau_col(ngrid)
      real flux1,flux2,flux3,ts1,ts2,ts3

      real qcol(ngridmx,nqmx)
!     Pluto adding variables
      real vmr_ch4(ngridmx)  ! vmr ch4
      real vmr_co(ngridmx)  ! vmr co
      real rho(ngridmx,nlayermx) ! density
      real zrho_ch4(ngridmx,nlayermx) ! density methane kg.m-3
      real zrho_co(ngridmx,nlayermx) ! density CO kg.m-3
      real zrho_haze(ngridmx,nlayermx) ! density haze kg.m-3
      real zdqrho_photprec(ngridmx,nlayermx) !photolysis rate kg.m-3.s-1
      real zq1temp_ch4(ngridmx) !
      real qsat_ch4(ngridmx) !
      real qsat_ch4_l1(ngridmx) !

!      CHARACTER(LEN=20) :: txt ! to temporarly store text for eddy tracers
      real profmmr(ngridmx,nlayermx) ! fixed profile of haze if haze_proffix
      real sensiblehf1(ngridmx) ! sensible heat flux
      real sensiblehf2(ngridmx) ! sensible heat flux

!     included by RW to test energy conservation
      REAL dEtot, dEtots, masse, vabs, dvabs, ang0

!     included by RW to allow variations in cp with location
      REAL cpp3D(ngridmx,nlayermx)   ! specific heat capacity at const. pressure
      REAL rcp3D(ngridmx,nlayermx)   ! R / specific heat capacity
      real cppNI, rcpNI, nnu ! last one just for Seb version
      REAL zpopskNI(ngridmx,nlayermx)

!     included by RW to make 1D saves not every timestep
      integer countG1D
!     integer countG3D,saveG3D
      save countG1D
!     save countG3D,saveG3D
      CHARACTER(len=10) :: tname

c=======================================================================

c 1. Initialisation:
c -----------------

c  1.1   Initialisation only at first call
c  ---------------------------------------
      IF (firstcall) THEN
         if(ngrid.eq.1)then
            saveG1D=1
            countG1D=1
         endif

c        variables set to 0
c        ~~~~~~~~~~~~~~~~~~
         call zerophys(ngrid*nlayer,dtrad)
         call zerophys(ngrid,fluxrad)
         call zerophys(ngrid*nlayer*nq,zdqsed)
         call zerophys(ngrid*nq, zdqssed)
         call zerophys(ngrid*nlayer*nq,zdqadj)

c        initialize tracer names, indexes and properties
c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         tracerdyn=tracer   ! variable tracer in dynamics
         IF (tracer) THEN
            CALL initracer()
         ENDIF  ! end tracer

!        Initialize aerosol indexes: done in initracer
!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        call iniaerosol()

c        read startfi (initial parameters)
c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         CALL phyetat0 ("startfi.nc",0,0,
     &         nsoilmx,nq,
     &         day_ini,time_phys,
     &         tsurf,tsoil,emis,q2,qsurf)

         if (pday.ne.day_ini) then
           write(*,*) "PHYSIQ: ERROR: bad synchronization between ",
     &                "physics and dynamics"
           write(*,*) "dynamics day: ",pday
           write(*,*) "physics day:  ",day_ini
           stop
         endif

         write (*,*) 'In physiq day_ini =', day_ini
         ptime0=ptime
         write (*,*) 'In physiq ptime0 =', ptime

c        Initialize albedo and orbital calculation
c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         CALL surfini(ngrid,albedo)
         CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)

         saveday=pday
         savedeclin=0.
         !adjust=0. ! albedo adjustment for convergeps

c        Initialize surface inventory and geopotential for paleo mode
c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         if (paleo) then
           write(*,*) 'Paleo time tpal = ',tpal
           DO ig=1,ngrid
             phisfinew(ig)=phisfi(ig)-qsurf(ig,1)*g/1000. ! topo of bedrock below the ice
           ENDDO
         endif

         if (conservn2) then
            write(*,*) 'conservn2 initial'
            call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
         endif

c        initialize soil
c        ~~~~~~~~~~~~~~~
         IF (callsoil) THEN
            CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
     s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
         ELSE
            PRINT*,'WARNING! Thermal conduction in the soil turned off'
            DO ig=1,ngrid
               capcal(ig)=1.e5
               fluxgrd(ig)=0.
            ENDDO
         ENDIF
         icount=1

      ENDIF        !  (end of "if firstcall")

c ---------------------------------------------------
c 1.2   Initializations done at every physical timestep:
c ---------------------------------------------------
c
      IF (ngrid.NE.ngridmx) THEN
         PRINT*,'STOP in PHYSIQ'
         PRINT*,'Probleme de dimensions :'
         PRINT*,'ngrid     = ',ngrid
         PRINT*,'ngridmx   = ',ngridmx
         STOP
      ENDIF

c     Initialize various variables
c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      call zerophys(ngrid*nlayer, pdv)
      call zerophys(ngrid*nlayer, pdu)
      call zerophys(ngrid*nlayer,pdt)
      call zerophys(ngrid*nlayer*nq, pdq)
      call zerophys(ngrid, pdpsrf)
      call zerophys(ngrid, zflubid)
      call zerophys(ngrid, zdtsurf)
      zdtsdif(1:ngrid)=0
      call zerophys(ngrid*nq, dqsurf)

      if (conservn2) then
         write(*,*) 'conservn2 iniloop'
         call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
      endif

      igout=ngrid/2+1

      zday=pday+ptime ! compute time, in sols (and fraction thereof)

c     Compute Solar Longitude (Ls) :
c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      if (season) then
         call solarlong(zday,zls)
      else
         call solarlong(float(day_ini),zls)
      end if

c     Get Lyman alpha flux at specific Ls
      if (haze) then
         call lymalpha(zls,zfluxuv)
         print*, 'Haze lyman-alpha zls,zfluxuv=',zls,zfluxuv
      end if

c     Compute geopotential between layers
c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DO l=1,nlayer
         DO ig=1,ngrid
            zzlay(ig,l)=pphi(ig,l)/g
         ENDDO
      ENDDO
      DO ig=1,ngrid
         zzlev(ig,1)=0.
         zzlev(ig,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
      ENDDO
      DO l=2,nlayer
         DO ig=1,ngrid
            z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
            z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
            zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
         ENDDO
      ENDDO

!     Potential temperature calculation not the same in physiq and dynamic...
c     Compute potential temperature zh
      DO l=1,nlayer
         DO ig=1,ngrid
            zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
         ENDDO
      ENDDO

c --------------------------------------------------------
c    1.3 thermosphere
c --------------------------------------------------------

c ajout de la conduction depuis la thermosphere
      IF (callconduct) THEN

          call conduction (ngrid,nlayer,ptimestep,
     $                  pplay,pt,zzlay,zzlev,zdtconduc,tsurf)
          DO l=1,nlayer
             DO ig=1,ngrid
                pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
             ENDDO
          ENDDO

      ENDIF

! ajout de la viscosite moleculaire
      IF (callmolvis) THEN
          call molvis(ngrid,nlayer,ptimestep,
     $                   pplay,pt,zzlay,zzlev,
     $                   zdtconduc,pu,tsurf,zdumolvis)
          call molvis(ngrid,nlayer,ptimestep,
     $                   pplay,pt,zzlay,zzlev,
     $                   zdtconduc,pv,tsurf,zdvmolvis)

          DO l=1,nlayer
             DO ig=1,ngrid
             ! pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l)
                pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
                pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
             ENDDO
          ENDDO
      ENDIF

      IF (callmoldiff) THEN
           call moldiff_red(ngrid,nlayer,nq,
     &                   pplay,pplev,pt,pdt,pq,pdq,ptimestep,
     &                   zzlay,zdtconduc,zdqmoldiff)

           DO l=1,nlayer
              DO ig=1,ngrid
                 DO iq=1, nq
                  pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
                 ENDDO
              ENDDO
           ENDDO
      ENDIF

      if (conservn2) then
         write(*,*) 'conservn2 thermo'
         call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
      endif

c-----------------------------------------------------------------------
c    2. Compute radiative tendencies :
c---------  --------------------------------------------------------------
c     Saving qsurf to compute paleo flux condensation/sublimation
      DO iq=1, nq
         DO ig=1,ngrid
             IF (qsurf(ig,iq).lt.0.) then
                qsurf(ig,iq)=0.
             ENDIF
             qsurf1(ig,iq)=qsurf(ig,iq)
         ENDDO
      ENDDO

      IF (callrad) THEN
         IF( MOD(icount-1,iradia).EQ.0) THEN

c          Local stellar zenith angle
c          ~~~~~~~~~~~~~~~~~~~~~~~~~~
           IF (triton) then
              CALL orbitetriton(zls,zday,dist_sol,declin)
           ELSE
              CALL orbite(zls,dist_sol,declin)
           ENDIF

           IF (diurnal) THEN
               ztim1=SIN(declin)
               ztim2=COS(declin)*COS(2.*pi*(zday-.5))
               ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))

               CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
     s         ztim1,ztim2,ztim3,mu0,fract)

           ELSE
               CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
               ! WARNING: this function appears not to work in 1D
           ENDIF

c         Albedo /IT changes depending of surface ices (only in 3D)
c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
           if (ngrid.ne.1) then

            !! Specific to change albedo of N2 so that Psurf
            !! converges toward 1.4 Pa in "1989" seasons for Triton
            !! converges toward 1.1 Pa in "2015" seasons for Pluto
            if (convergeps) then
              if (triton) then
                ! 1989 declination
                if (declin*180./pi.gt.-46..and.declin*180./pi.lt.-45.
     &          .and.zday.gt.saveday+1000.
     &          .and.declin.lt.savedeclin) then
                  call globalaverage2d(ngrid,pplev(:,1),globave)
                  if (globave.gt.1.5) then
                        adjust=adjust+0.005
                  else if (globave.lt.1.3) then
                        adjust=adjust-0.005
                  endif
                  saveday=zday
                endif
                call surfprop(ngrid,nq,fract,qsurf,tsurf,tidat,capcal,
     &           adjust,dist_sol,albedo,emis,flusurfold,ptimestep,zls)
                savedeclin=declin
              else
                ! Pluto : 2015 declination current epoch
                if (declin*180./pi.gt.50.and.declin*180./pi.lt.51.
     &          .and.zday.gt.saveday+10000.
     &          .and.declin.gt.savedeclin) then
                  call globalaverage2d(ngrid,pplev(:,1),globave)
                  if (globave.gt.1.2) then
                        adjust=adjust+0.005
                  else if (globave.lt.1.) then
                        adjust=adjust-0.005
                  endif
                  saveday=zday
                endif
                call surfprop(ngrid,nq,fract,qsurf,tsurf,tidat,capcal,
     &           adjust,dist_sol,albedo,emis,flusurfold,ptimestep,zls)
                savedeclin=declin
              endif
            else
              ! Surface properties
              call surfprop(ngrid,nq,fract,qsurf,tsurf,tidat,
     &        capcal,adjust,dist_sol,albedo,emis,flusurfold,ptimestep,
     &         zls)
            end if
!TB22
           else  ! here ngrid=1
              call surfprop(ngrid,nq,fract,qsurf,tsurf,tidat,
     &        capcal,adjust,dist_sol,albedo,emis,flusurfold,ptimestep,
     &         zls)

           end if ! if ngrid ne 1

c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c          Fixed zenith angle in 1D
c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
           if(ngrid.eq.1.and.globmean1d) then
             !global mean 1D radiative equiilium
             ig=1
             mu0(ig)     = 0.5
             fract(ig)= 1/(4*mu0(ig))
             !print*,'WARNING GLOBAL MEAN CALCULATION'
           endif

c          Call main radiative transfer scheme
c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

c          Radiative transfer
c          ------------------
           if (corrk) then
c             print*,'TB22 start corrk',tsurf   ! Bug testphys1d nodebug
             CALL callcorrk(icount,ngrid,nlayer,pq,nq,qsurf,
     &             albedo,emis,mu0,pplev,pplay,pt,
     &             tsurf,fract,dist_sol,igout,aerosol,cpp3D,
     &             zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,
     &             fluxtop_sw,fluxtop_dn,reffrad,tau_col,ptime,pday,
     &             firstcall,lastcall,zzlay)

c            Radiative flux from the sky absorbed by the surface (W.m-2)
c            ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
             DO ig=1,ngrid
                fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
     &               +fluxsurf_sw(ig)*(1.-albedo(ig))
             ENDDO

c            Net atmospheric radiative heating/cooling rate (K.s-1):
c            ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
             DO l=1,nlayer
                DO ig=1,ngrid
                 dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
                ENDDO
             ENDDO

           else ! corrk = F

c            Atmosphere has no radiative effect
c            ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
             DO ig=1,ngrid
                 fluxtop_dn(ig)  = fract(ig)*mu0(ig)*Fat1AU/dist_sol**2
                 fluxrad_sky(ig) = fluxtop_dn(ig)*(1.-albedo(ig))
                 fluxtop_sw(ig)  = fluxtop_dn(ig)*albedo(ig)
                 fluxtop_lw(ig)  = emis(ig)*stephan*tsurf(ig)**4
             ENDDO ! radiation skips the atmosphere entirely
             DO l=1,nlayer
                 DO ig=1,ngrid
                    dtrad(ig,l)=0.0
                 ENDDO
             ENDDO ! hence no atmospheric radiative heating

           endif                ! if corrk

         ENDIF ! of if(mod(icount-1,iradia).eq.0) ie radiative timestep

!        Transformation of the radiative tendencies:
!        -------------------------------------------

!        Net radiative surface flux (W.m-2)
!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

         DO ig=1,ngrid
           zplanck(ig)=tsurf(ig)*tsurf(ig)
           zplanck(ig)=emis(ig)*
     &          stephan*zplanck(ig)*zplanck(ig)
           fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
         ENDDO

         DO l=1,nlayer
           DO ig=1,ngrid
              pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
           ENDDO
         ENDDO

      ENDIF ! of IF (callrad)

      if (conservn2) then
         write(*,*) 'conservn2 radiat'
         call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
      endif

!-----------------------------------------------------------------------
!    4. Vertical diffusion (turbulent mixing):
!-----------------------------------------------------------------------

      IF (calldifv) THEN

         DO ig=1,ngrid
            zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
         ENDDO

         CALL zerophys(ngrid*nlayer,zdum1)
         CALL zerophys(ngrid*nlayer,zdum2)
         do l=1,nlayer
            do ig=1,ngrid
               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
            enddo
         enddo

c        Calling vdif (Martian version WITH N2 condensation)
         CALL vdifc(ngrid,nlayer,nq,zpopsk,
     &        ptimestep,capcal,lwrite,
     &        pplay,pplev,zzlay,zzlev,z0,
     &        pu,pv,zh,pq,pt,tsurf,emis,qsurf,
     &        zdum1,zdum2,zdh,pdq,pdt,zflubid,
     &        zdudif,zdvdif,zdhdif,zdtsdif,q2,
     &        zdqdif,zdqsdif,qsat_ch4,qsat_ch4_l1) !,zq1temp_ch4,qsat_ch4)

         DO l=1,nlayer
            DO ig=1,ngrid
               pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
               pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
               pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)

               zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
            ENDDO
         ENDDO

         DO ig=1,ngrid
            zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
         ENDDO

         bcond=1./tcond1p4Pa
         acond=r/lw_n2

         if (tracer) then
          DO iq=1, nq
           DO l=1,nlayer
            DO ig=1,ngrid
             pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
            ENDDO
           ENDDO
          ENDDO

          DO iq=1, nq
              DO ig=1,ngrid
                 dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
              ENDDO
          ENDDO

         end if ! of if (tracer)

!------------------------------------------------------------------
! test methane conservation
c         if(methane)then
c           call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
c     &          ptimestep,pplev,zdqdif,zdqsdif,'CH4',' vdifc ')
c         endif  ! methane
!------------------------------------------------------------------
! test CO conservation
c         if(carbox)then
c           call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
c     &          ptimestep,pplev,zdqdif,zdqsdif,'CO ',' vdifc ')
c         endif  ! carbox
!------------------------------------------------------------------

      ELSE   ! case with calldifv=F
         ztim1=4.*stephan*ptimestep
         DO ig=1,ngrid
           ztim2=ztim1*emis(ig)*tsurf(ig)**3
           z1=capcal(ig)*tsurf(ig)+
     s     ztim2*tsurf(ig)+ (fluxrad(ig)+fluxgrd(ig))*ptimestep
           z2= capcal(ig)+ztim2
           zdtsurf(ig)=(z1/z2 - tsurf(ig))/ptimestep

           ! for output:
           !dplanck(ig)=4.*stephan*ptimestep*emis(ig)*tsurf(ig)**3
         ENDDO

c        ------------------------------------------------------------------
c        Methane surface sublimation and condensation in fast model (nogcm)
c        ------------------------------------------------------------------
         if ((methane).and.(fast).and.condmetsurf) THEN

           call ch4surf(ngrid,nlayer,nq,ptimestep,
     &     tsurf,zdtsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf,
     &     zdqch4fast,zdqsch4fast)

           DO ig=1,ngrid
              dqsurf(ig,igcm_ch4_ice)= dqsurf(ig,igcm_ch4_ice) +
     &                    zdqsch4fast(ig)
              pdq(ig,1,igcm_ch4_gas)= pdq(ig,1,igcm_ch4_gas) +
     &                    zdqch4fast(ig)
              zdtsurf(ig)=zdtsurf(ig)+lw_ch4*zdqsch4fast(ig)/capcal(ig)
           ENDDO
         end if
c        ------------------------------------------------------------------
c        CO surface sublimation and condensation in fast model (nogcm)
c        ------------------------------------------------------------------
         if ((carbox).and.(fast).and.condcosurf) THEN

           call cosurf(ngrid,nlayer,nq,ptimestep,
     &     tsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf,
     &     zdqcofast,zdqscofast)

           DO ig=1,ngrid
              dqsurf(ig,igcm_co_ice)= dqsurf(ig,igcm_co_ice) +
     &                    zdqscofast(ig)
              pdq(ig,1,igcm_co_gas)= pdq(ig,1,igcm_co_gas) +
     &                    zdqcofast(ig)
              zdtsurf(ig)=zdtsurf(ig)+lw_co*zdqscofast(ig)/capcal(ig)
           ENDDO
         end if

      ENDIF ! of IF (calldifv)

      if (conservn2) then
        write(*,*) 'conservn2 diff'
        call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)+
     &                                  dqsurf(:,1)*ptimestep)
      endif

!------------------------------------------------------------------
! test CO conservation
!         if(carbox)then
!           call testconservfast(ngrid,nlayer,nq,
!     &          ptimestep,pplev,zdqcofast,zdqscofast,'CO ',' vdifc ')
!         endif  ! carbox
!------------------------------------------------------------------

!-----------------------------------------------------------------------
!   5. Dry convective adjustment:
!   -----------------------------

      IF(calladj) THEN

         DO l=1,nlayer
            DO ig=1,ngrid
                  zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
            ENDDO
         ENDDO
         CALL zerophys(ngrid*nlayer,zduadj)
         CALL zerophys(ngrid*nlayer,zdvadj)
         CALL zerophys(ngrid*nlayer,zdhadj)

            CALL convadj(ngrid,nlayer,nq,ptimestep,
     &           pplay,pplev,zpopsk,
     &           pu,pv,zh,pq,
     &           pdu,pdv,zdh,pdq,
     &           zduadj,zdvadj,zdhadj,
     &           zdqadj)


         DO l=1,nlayer
            DO ig=1,ngrid
               pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
               pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
               pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
               zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
            ENDDO
         ENDDO

         if(tracer) then
           DO iq=1, nq
            DO l=1,nlayer
              DO ig=1,ngrid
                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
              ENDDO
            ENDDO
           ENDDO
         end if

      ENDIF ! of IF(calladj)

!------------------------------------------------------------------
! test methane conservation
c         if(methane)then
c           call testchange(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
c     &          ptimestep,pplev,zdqadj,'CH4','calladj')
c         endif  ! methane
!------------------------------------------------------------------
! test CO conservation
c         if(carbox)then
c           call testchange(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
c     &          ptimestep,pplev,zdqadj,'CO ','calladj')
c         endif  ! carbox
!------------------------------------------------------------------

!-----------------------------------------------------------------------
!   6. Nitrogen condensation-sublimation:
!   -------------------------------------------

      IF (N2cond) THEN
         if (.not.tracer.or.igcm_n2.eq.0) then
             print*,'We need a n2 ice tracer to condense n2'
             call abort
         endif

c          write(*,*) 'before N2 condens :'
c          call massn2(ngrid,nlayer,pplev,qsurf(:,1),dqsurf(:,igcm_n2),
c     &        pdpsrf,ptimestep)

         if (tracer) then
           zdqc(:,:,:)=0.
           zdqsc(:,:)=0.
           call condens_n2(ngrid,nlayer,nq,ptimestep,
     &        capcal,pplay,pplev,tsurf,pt,
     &        pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
     &        qsurf(1,igcm_n2),albedo,emis,
     &        zdtc,zdtsurfc,pdpsrf,zduc,zdvc,
     &        zdqc,zdqsc(1,igcm_n2))
         endif

         !!! TB temporaire
         !zdtc=0.
         !zdvc=0.
         !zduc=0.
         !zdqc=0.
         !zdqsc=0.
         !zdtsurfc=0.

         DO l=1,nlayer
           DO ig=1,ngrid
             pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
             pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
             pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
           ENDDO
         ENDDO
         DO ig=1,ngrid
            zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
         ENDDO

         DO iq=1, nq
c           TB: option eddy not condensed with N2
c           txt=noms(iq)
c           IF (txt(1:4).ne."eddy") THEN
            DO l=1,nlayer
              DO ig=1,ngrid
               pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
              ENDDO
            ENDDO
c           ENDIF
         ENDDO
         DO ig=1,ngrid
            dqsurf(ig,igcm_n2)= dqsurf(ig,igcm_n2) + zdqsc(ig,igcm_n2)
         ENDDO

      ENDIF  ! of IF (N2cond)

      if (conservn2) then
       write(*,*) 'conservn2 n2cond'
       call testconservmass(ngrid,nlayer,pplev(:,1)+
     &      pdpsrf(:)*ptimestep,qsurf(:,1)+dqsurf(:,1)*ptimestep)
      endif

!------------------------------------------------------------------
! test n2 conservation for one tendency / pplevis not updated here w pdpsrf
!         if(tracer)then
!          call testconserv(ngrid,nlayer,nq,igcm_n2,igcm_n2,
!     &          ptimestep,pplev,zdqc,zdqsc,'N2 ',' n2cond')
!          call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
!         endif  ! n2
!------------------------------------------------------------------
! test methane conservation
!         if(methane)then
!           call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
!     &          ptimestep,pplev,zdqc,zdqsc,'CH4',' n2cond')
!         endif  ! methane
!------------------------------------------------------------------
! test CO conservation
c         if(carbox)then
c           call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
c     &          ptimestep,pplev,zdqc,zdqsc,'CO ',' n2cond')
c         endif  ! carbox
!------------------------------------------------------------------

c-----------------------------------------------------------------------
c   7. Specific parameterizations for tracers
c   -----------------------------------------

      if (tracer) then

c   7a. Methane, CO, and ice
c      ---------------------------------------
c      Methane ice condensation in the atmosphere
c      ----------------------------------------
       zdqch4cloud(:,:,:)=0.
       if ((methane).and.(metcloud).and.(.not.fast)) THEN
           call ch4cloud(ngrid,nlayer,naerkind,ptimestep,
     &                pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
     &                pq,pdq,zdqch4cloud,zdqsch4cloud,zdtch4cloud,
     &                nq,rice_ch4)

           DO l=1,nlayer
              DO ig=1,ngrid
                 pdq(ig,l,igcm_ch4_gas)=pdq(ig,l,igcm_ch4_gas)+
     &                                   zdqch4cloud(ig,l,igcm_ch4_gas)
                 pdq(ig,l,igcm_ch4_ice)=pdq(ig,l,igcm_ch4_ice)+
     &                                   zdqch4cloud(ig,l,igcm_ch4_ice)
              ENDDO
           ENDDO

           ! Increment methane ice surface tracer tendency
           DO ig=1,ngrid
               dqsurf(ig,igcm_ch4_ice)=dqsurf(ig,igcm_ch4_ice)+
     &                               zdqsch4cloud(ig,igcm_ch4_ice)
           ENDDO

           ! update temperature tendancy
           DO ig=1,ngrid
             DO l=1,nlayer
               pdt(ig,l)=pdt(ig,l)+zdtch4cloud(ig,l)
             ENDDO
           ENDDO
       else
         rice_ch4(:,:)=0 ! initialization needed for callsedim
       end if

c      ---------------------------------------
c      CO ice condensation in the atmosphere
c      ----------------------------------------
       zdqcocloud(:,:,:)=0.
       IF ((carbox).and.(monoxcloud).and.(.not.fast)) THEN
           call cocloud(ngrid,nlayer,naerkind,ptimestep,
     &                pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
     &                pq,pdq,zdqcocloud,zdqscocloud,zdtcocloud,
     &                nq,rice_co,qsurf(1,igcm_n2),dqsurf(1,igcm_n2))

           DO l=1,nlayer
               DO ig=1,ngrid
                 pdq(ig,l,igcm_co_gas)=pdq(ig,l,igcm_co_gas)+
     &                                   zdqcocloud(ig,l,igcm_co_gas)
                 pdq(ig,l,igcm_co_ice)=pdq(ig,l,igcm_co_ice)+
     &                                   zdqcocloud(ig,l,igcm_co_ice)
               ENDDO
           ENDDO

           ! Increment CO ice surface tracer tendency
           DO ig=1,ngrid
            dqsurf(ig,igcm_co_ice)=dqsurf(ig,igcm_co_ice)+
     &                               zdqscocloud(ig,igcm_co_ice)
           ENDDO

           ! update temperature tendancy
           DO ig=1,ngrid
             DO l=1,nlayer
               pdt(ig,l)=pdt(ig,l)+zdtcocloud(ig,l)
             ENDDO
           ENDDO
       ELSE
         rice_co(:,:)=0 ! initialization needed for callsedim
       END IF  ! of IF (carbox)

!------------------------------------------------------------------
! test methane conservation
c         if(methane)then
c           call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
c     &         ptimestep,pplev,zdqch4cloud,zdqsch4cloud,'CH4','ch4clou')
c         endif  ! methane
!------------------------------------------------------------------
! test CO conservation
c         if(carbox)then
c           call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
c     &          ptimestep,pplev,zdqcocloud,zdqscocloud,'CO ','cocloud')
c         endif  ! carbox
!------------------------------------------------------------------

c   7b. Haze particle production
c     -------------------
       IF (haze) THEN

          zdqphot_prec(:,:)=0.
          zdqphot_ch4(:,:)=0.
          zdqhaze(:,:,:)=0.
          ! Forcing to a fixed haze profile if haze_proffix
          if (haze_proffix.and.i_haze.gt.0.) then
           call haze_prof(ngrid,nlayer,zzlay,pplay,pt,
     &                               reffrad,profmmr)
           zdqhaze(:,:,i_haze)=(profmmr(:,:)-pq(:,:,igcm_haze))
     &                                               /ptimestep
          else
           call hazecloud(ngrid,nlayer,nq,ptimestep,
     &         pplay,pplev,pq,pdq,dist_sol,mu0,zfluxuv,zdqhaze,
     &         zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
          endif

          DO iq=1, nq ! should be updated
             DO l=1,nlayer
               DO ig=1,ngrid
                    pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqhaze(ig,l,iq)
               ENDDO
             ENDDO
          ENDDO

       ENDIF

       IF (fast.and.fasthaze) THEN
         call prodhaze(ngrid,nlayer,nq,ptimestep,pplev,pq,pdq,dist_sol,
     &             mu0,declin,zdqprodhaze,zdqsprodhaze,gradflux,fluxbot,
     &             fluxlym_sol_bot,fluxlym_ipm_bot,flym_sol,flym_ipm)

         DO ig=1,ngrid
           pdq(ig,1,igcm_ch4_gas)=pdq(ig,1,igcm_ch4_gas)+
     &                                    zdqprodhaze(ig,igcm_ch4_gas)
           pdq(ig,1,igcm_prec_haze)=pdq(ig,1,igcm_prec_haze)+
     &                                  zdqprodhaze(ig,igcm_prec_haze)
           pdq(ig,1,igcm_haze)=abs(pdq(ig,1,igcm_haze)+
     &                                  zdqprodhaze(ig,igcm_haze))
           qsurf(ig,igcm_haze)= qsurf(ig,igcm_haze)+
     &                                      zdqsprodhaze(ig)*ptimestep
         ENDDO

       ENDIF

c   7c. Aerosol particles
c     -------------------

c      -------------
c      Sedimentation
c      -------------
       IF (sedimentation) THEN
         call zerophys(ngrid*nlayer*nq, zdqsed)
         call zerophys(ngrid*nq, zdqssed)

         call callsedim(ngrid,nlayer, ptimestep,
     &        pplev,zzlev, pt,rice_ch4,rice_co,
     &        pq, pdq, zdqsed, zdqssed,nq,pphi)

         DO iq=1, nq
           DO l=1,nlayer
             DO ig=1,ngrid
               pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
             ENDDO
           ENDDO
         ENDDO
         DO iq=2, nq
           DO ig=1,ngrid
              dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
           ENDDO
         ENDDO

       END IF   ! of IF (sedimentation)

!------------------------------------------------------------------
! test methane conservation
c         if(methane)then
c           call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice,
c     &          ptimestep,pplev,zdqsed,zdqssed,'CH4',' sedim ')
c         endif  ! methane
!------------------------------------------------------------------
! test CO conservation
c         if(carbox)then
c           call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice,
c     &          ptimestep,pplev,zdqsed,zdqssed,'CO ',' sedim ')
c         endif  ! carbox
!------------------------------------------------------------------

c   7d. Updates
c     ---------
!        write(*,*) 'before update qsurf:'
!        call massn2(ngrid,nlayer,pplev,qsurf(:,1),dqsurf(:,igcm_n2),
!     &        pdpsrf,ptimestep)

c       ---------------------------------
c       Updating tracer budget on surface (before spread of glacier)
c       ---------------------------------
        DO iq=1, nq
          DO ig=1,ngrid
            qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
          ENDDO
        ENDDO

      endif !  of if (tracer)

      if (conservn2) then
         write(*,*) 'conservn2 tracer'
         call testconservmass(ngrid,nlayer,pplev(:,1)+
     &       pdpsrf(:)*ptimestep,qsurf(:,1))
      endif

      DO ig=1,ngrid
        flusurf(ig,igcm_n2)=(qsurf(ig,igcm_n2)-
     &                              qsurf1(ig,igcm_n2))/ptimestep
        flusurfold(ig,igcm_n2)=flusurf(ig,igcm_n2)
        if (methane) then
          flusurf(ig,igcm_ch4_ice)=(qsurf(ig,igcm_ch4_ice)-
     &                              qsurf1(ig,igcm_ch4_ice))/ptimestep
          flusurfold(ig,igcm_ch4_ice)=flusurf(ig,igcm_ch4_ice)
        endif
        if (carbox) then
          flusurf(ig,igcm_co_ice)=(qsurf(ig,igcm_co_ice)-
     &                              qsurf1(ig,igcm_co_ice))/ptimestep
          !flusurfold(ig,igcm_co_ice)=flusurf(ig,igcm_co_ice)
        endif
      ENDDO

      !! Special source of haze particle !
      ! todo: should be placed in haze and use tendency of n2 instead of flusurf
      IF (source_haze) THEN
             call hazesource(ngrid,nlayer,nq,ptimestep,
     &                       pplev,flusurf,mu0,zdq_source)

             DO iq=1, nq
               DO l=1,nlayer
                 DO ig=1,ngrid
                    pdq(ig,l,iq)=pdq(ig,l,iq)+zdq_source(ig,l,iq)
                 ENDDO
               ENDDO
             ENDDO
      ENDIF

!-----------------------------------------------------------------------
!   9. Surface and sub-surface soil temperature
!-----------------------------------------------------------------------
!   For diagnostic
!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DO ig=1,ngrid
         rho(ig,1) = pplay(ig,1)/(r*pt(ig,1))
         sensiblehf1(ig)=rho(ig,1)*cpp*(0.4/log(zzlay(ig,1)/z0))**2*
     &                   (pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1))**0.5*
     &                   (tsurf(ig)-pt(ig,1))
         sensiblehf2(ig)=zflubid(ig)-capcal(ig)*zdtsdif(ig)

      ENDDO

!   9.1 Increment Surface temperature:
!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DO ig=1,ngrid
         tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
      ENDDO

      ! Prevent surface (.e.g. non volatile ch4) to exceed max temperature
      ! Lellouch et al., 2000,2011
      IF (tsurfmax) THEN
        DO ig=1,ngrid
         if (albedo(ig).gt.albmin_ch4.and.
     &                      qsurf(ig,igcm_n2).eq.0.) then
              tsurf(ig)=min(tsurf(ig),54.)
         endif
        ENDDO
      ENDIF
!
!   9.2 Compute soil temperatures and subsurface heat flux:
!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      IF (callsoil) THEN
         CALL soil(ngrid,nsoilmx,.false.,tidat,
     &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
      ENDIF

      ! For output :
      tidat_out(:,:)=0.
      DO l=1,min(nlayermx,nsoilmx)
         tidat_out(:,l)=tidat(:,l)
      ENDDO

!   9.3 multiply tendencies of cond/subli for paleo loop only in the
!       last Pluto year of the simulation
!       Year day must be adapted in the startfi for each object
!       Paleo uses year_day to calculate the annual mean tendancies
!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      IF (paleo) then
         if (zday.gt.day_ini+ptime0+nday-year_day) then
            DO iq=1,nq
             DO ig=1,ngrid
               qsurfyear(ig,iq)=qsurfyear(ig,iq)+
     &                        (qsurf(ig,iq)-qsurf1(ig,iq))  !kg m-2 !ptimestep
             ENDDO
            ENDDO
         endif
      endif

!   9.4 Glacial flow at each timestep glastep or at lastcall
!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      IF (fast.and.glaflow) THEN
         if ((mod(zday-day_ini-ptime0,glastep)).lt.1.
     &                                              .or.lastcall) then
           IF (lastcall) then
            dstep=mod(zday-day_ini-ptime0,glastep)*daysec
           else
            dstep=glastep*daysec
           endif
           zdqflow(:,:)=qsurf(:,:)
           IF (paleo) then
             call spreadglacier_paleo(ngrid,nq,qsurf,
     &                               phisfinew,dstep,tsurf)
           else
             call spreadglacier_simple(ngrid,nq,qsurf,dstep)
           endif
           zdqflow(:,:)=(zdqflow(:,:)-qsurf(:,:))/dstep

           if (conservn2) then
            write(*,*) 'conservn2 glaflow'
            call testconservmass(ngrid,nlayer,pplev(:,1)+
     &       pdpsrf(:)*ptimestep,qsurf(:,1))
           endif

         endif
      ENDIF
!-----------------------------------------------------------------------
!  10. Perform diagnostics and write output files
!-----------------------------------------------------------------------
!     ---------------------------------------------------------
!     Check the energy balance of the simulation during the run
!     ---------------------------------------------------------
      flux1 = 0
      flux2 = 0
      flux3 = 0
      do ig=1,ngrid
	    flux1 = flux1 +
     &        area(ig)*(fluxtop_dn(ig) - fluxtop_sw(ig))
	    flux2 = flux2 + area(ig)*fluxtop_lw(ig)
	    flux3 = flux3 + area(ig)*fluxtop_dn(ig)
            fluxabs_sw(ig)=fluxtop_dn(ig) - fluxtop_sw(ig)

      end do
      print*,'Incident solar flux, absorbed solar flux, OLR (W/m2)'
      print*, flux3/totarea, '      ',  flux1/totarea ,
     & '           = ', flux2/totarea

      if(meanOLR)then
         ! to record global radiative balance
         open(92,file="rad_bal.out",form='formatted',access='append')
         write(92,*) zday,flux3/totarea,flux1/totarea,flux2/totarea
         close(92)
      endif

!     Surface temperature information
      ts1 = 0
      ts2 = 999
      ts3 = 0
      do ig=1,ngrid
         ts1 = ts1 + area(ig)*tsurf(ig)
         ts2 = min(ts2,tsurf(ig))
      	 ts3 = max(ts3,tsurf(ig))
      end do
!      print*,'Mean Tsurf =',ts1/totarea , ' Min Tsurf=',ts2,
!     &          ' Max Tsurf =',ts3
!      print*,'Max Tsurf=',ts3,'Min Tsurf=',ts2

!    -------------------------------
!    Dynamical fields incrementation
!    -------------------------------
! (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
      ! temperature, zonal and meridional wind
      if(firstcall) ztprevious(:,:)=pt(:,:)
      DO l=1,nlayer
        DO ig=1,ngrid
          zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
          zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
          zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
      ! diagnostic
            zdtdyn(ig,l)=ztprevious(ig,l)-pt(ig,l)
            ztprevious(ig,l)=zt(ig,l)
        ENDDO
      ENDDO
!      if(firstcall)call zerophys(zdtdyn)
!      if(firstcall) zdtdyn(:,:)=0.
      ! dynamical heating diagnostic
      DO ig=1,ngrid
         fluxdyn(ig)=0.
         DO l=1,nlayer
               fluxdyn(ig)=fluxdyn(ig) - (zdtdyn(ig,l)/ptimestep)
     &              *(pplev(ig,l)-pplev(ig,l+1))*cpp/g
         ENDDO
      ENDDO

!    -------------------------------
!    Other diagnostics
!    -------------------------------
     !! diagnostic qsat_lev1 (ps, T1)
!      CALL methanesat(ngridmx,zt(:,1),pplev(1,1),qsat_lev1(:),
!     &                                          qsurf(:,igcm_n2))

     !! eddy tracers : decay time
!      DO l=1,nlayer
!       DO ig=1,ngrid
!        pdq(ig,l,igcm_eddy1e6)=pdq(ig,l,igcm_eddy1e6)-
!     &          pq(ig,l,igcm_eddy1e6)*(1-exp(-ptimestep/1.e6))/ptimestep
!        pdq(ig,l,igcm_eddy1e7)=pdq(ig,l,igcm_eddy1e7)-
!     &          pq(ig,l,igcm_eddy1e7)*(1-exp(-ptimestep/1.e7))/ptimestep
!        pdq(ig,l,igcm_eddy5e7)=pdq(ig,l,igcm_eddy5e7)-
!     &          pq(ig,l,igcm_eddy5e7)*(1-exp(-ptimestep/5.e7))/ptimestep
!        pdq(ig,l,igcm_eddy1e8)=pdq(ig,l,igcm_eddy1e8)-
!     &          pq(ig,l,igcm_eddy1e8)*(1-exp(-ptimestep/1.e8))/ptimestep
!        pdq(ig,l,igcm_eddy5e8)=pdq(ig,l,igcm_eddy5e8)-
!     &          pq(ig,l,igcm_eddy5e8)*(1-exp(-ptimestep/5.e8))/ptimestep
!       ENDDO
!      ENDDO

!    -------------------------------
!    Update tracers
!    -------------------------------
      DO iq=1, nq
        DO l=1,nlayer
          DO ig=1,ngrid
            zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
          ENDDO
        ENDDO
      ENDDO

      DO ig=1,ngrid
          ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
      ENDDO
      call globalaverage2d(ngrid,ps,globave)

      ! pressure density
      IF (.not.fast) then !
       DO ig=1,ngrid
        DO l=1,nlayer
             zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
             zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
             rho(ig,l) = zplay(ig,l)/(r*zt(ig,l))
        ENDDO
        zplev(ig,nlayer+1)=pplev(ig,nlayer+1)/pplev(ig,1)*ps(ig)
       ENDDO
      ENDIF

!     ---------------------------------------------------------
!     Compute column amounts (kg m-2) if tracers are enabled
!     ---------------------------------------------------------
      call zerophys(ngrid*nq,qcol)
      if(tracer.and..not.fast)then
         do iq=1,nq
            DO ig=1,ngrid
               DO l=1,nlayer
                  qcol(ig,iq) = qcol(ig,iq) + zq(ig,l,iq) *
     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
               enddo
            enddo
         enddo
      endif

      if (methane) then
        IF (fast) then ! zq is the mixing ratio supposingly mixed in all atmosphere
          DO ig=1,ngrid
            vmr_ch4(ig)=zq(ig,1,igcm_ch4_gas)*
     &                        mmol(igcm_n2)/mmol(igcm_ch4_gas)*100.
          ENDDO
        ELSE
          DO ig=1,ngrid
!           compute vmr methane
            vmr_ch4(ig)=qcol(ig,igcm_ch4_gas)*
     &             g/ps(ig)*mmol(igcm_n2)/mmol(igcm_ch4_gas)*100.
!           compute density methane
            DO l=1,nlayer
               zrho_ch4(ig,l)=zq(ig,l,igcm_ch4_gas)*rho(ig,l)
            ENDDO
          ENDDO
        ENDIF
      endif

      if (carbox) then
        IF (fast) then
          DO ig=1,ngrid
            vmr_co(ig)=zq(ig,1,igcm_co_gas)*
     &                        mmol(igcm_n2)/mmol(igcm_co_gas)*100.
          ENDDO
        ELSE
         DO ig=1,ngrid
!           compute vmr CO
            vmr_co(ig)=qcol(ig,igcm_co_gas)*
     &             g/ps(ig)*mmol(igcm_n2)/mmol(igcm_co_gas)*100.
!           compute density CO
            DO l=1,nlayer
               zrho_co(ig,l)=zq(ig,l,igcm_co_gas)*rho(ig,l)
            ENDDO
         ENDDO
        ENDIF
      endif

      zrho_haze(:,:)=0.
      zdqrho_photprec(:,:)=0.
      IF (haze.and.aerohaze) then
       DO ig=1,ngrid
        DO l=1,nlayer
               zrho_haze(ig,l)=zq(ig,l,igcm_haze)*rho(ig,l)
               zdqrho_photprec(ig,l)=zdqphot_prec(ig,l)*rho(ig,l)
            ENDDO
        ENDDO
      ENDIF

      IF (fasthaze) then
       DO ig=1,ngrid
          qcol(ig,igcm_haze)=zq(ig,1,igcm_haze)*pplev(ig,1)/g
          qcol(ig,igcm_prec_haze)=zq(ig,1,igcm_prec_haze)*pplev(ig,1)/g
       ENDDO
      ENDIF

c     Info about Ls, declin...
      IF (fast) THEN
        write(*,*),'Ls=',zls*180./pi,' dec=',declin*180./pi
        write(*,*),'zday=',zday,' ps=',globave
       IF (lastcall) then
         write(*,*),'lastcall'
       ENDIF
      ELSE
       write(*,*),'Ls=',zls*180./pi,'decli=',declin*180./pi,'zday=',zday
      ENDIF

      lecttsoil=0 ! default value for lecttsoil
      call getin("lecttsoil",lecttsoil)
      IF (lastcall.and.(ngrid.EQ.1).and.(lecttsoil.eq.1)) THEN
      ! save tsoil temperature profile for 1D profile
         OPEN(13,file='proftsoil.out',form='formatted')
         DO i=1,nsoilmx
            write(13,*) tsoil(1,i)
         ENDDO
         CLOSE(13)
      ENDIF


      IF (ngrid.NE.1) THEN
!        -------------------------------------------------------------------
!        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',
!              there will have been the itau=itau+1 instruction and
!              a reset of 'time' (lastacll = .true. when itau+1= itaufin)
!              thus we store for time=time+dtvr

         IF(lastcall) THEN

            ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
            write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin

            if (paleo) then
             ! time range for tendencies of ice flux qsurfyear
             zdt_tot=year_day   ! Last year of simulation

             masslost(:)=0.
             massacc(:)=0.

             DO ig=1,ngrid
                ! update new reservoir of ice on the surface
                DO iq=1,nq
                 ! kg/m2 to be sublimed or condensed during paleoyears
                 qsurfyear(ig,iq)=qsurfyear(ig,iq)*
     &                       paleoyears*365.25/(zdt_tot*daysec/86400.)

                ! special case if we sublime the entire reservoir
                 IF (-qsurfyear(ig,iq).gt.qsurf(ig,iq)) THEN
                   masslost(iq)=masslost(iq)+realarea(ig)*
     &                    (-qsurfyear(ig,iq)-qsurf(ig,iq))
                   qsurfyear(ig,iq)=-qsurf(ig,iq)
                 ENDIF

                 IF (qsurfyear(ig,iq).gt.0.) THEN
                   massacc(iq)=massacc(iq)+realarea(ig)*qsurfyear(ig,iq)
                 ENDIF

                ENDDO
             ENDDO

             DO ig=1,ngrid
                DO iq=1,nq
                  qsurfpal(ig,iq)=qsurf(ig,iq)+qsurfyear(ig,iq)
                  IF (qsurfyear(ig,iq).gt.0.) THEN
                   qsurfpal(ig,iq)=qsurfpal(ig,iq)-
     &                   qsurfyear(ig,iq)*masslost(iq)/massacc(iq)
                  ENDIF
                ENDDO
             ENDDO
             ! Finally ensure conservation of qsurf
             DO iq=1,nq
               call globalaverage2d(ngrid,qsurf(:,iq),globaveice(iq))
               call globalaverage2d(ngrid,qsurfpal(:,iq),
     &                                         globavenewice(iq))
               IF (globavenewice(iq).gt.0.) THEN
                  qsurfpal(:,iq)=qsurfpal(:,iq)*
     &                               globaveice(iq)/globavenewice(iq)
               ENDIF
             ENDDO

             ! update new geopotential depending on the ice reservoir
             phisfipal(:)=phisfinew(:)+qsurfpal(:,igcm_n2)*g/1000.
             !phisfipal(ig)=phisfi(ig)


             if (kbo.or.triton) then !  case of Triton : we do not change the orbital parameters

               pdaypal=pday ! no increment of pdaypal to keep same evolution of the subsolar point
               eccpal=1.-periheli/((periheli+aphelie)/2.)    !no change of ecc
               peri_daypal=peri_day ! no change
               oblipal=obliquit     ! no change
               tpalnew=tpal
               adjustnew=adjust

             else  ! Pluto
               ! update new pday and tpal (Myr) to be set in startfi controle
               pdaypal=int(day_ini+paleoyears*365.25/6.3872)
               tpalnew=tpal+paleoyears*1.e-6  ! Myrs

               ! update new N2 ice adjustment (not tested yet on Pluto)
               adjustnew=adjust

               ! update milankovitch parameters : obliquity,Lsp,ecc
               call calcmilank(tpalnew,oblipal,peri_daypal,eccpal)
               !peri_daypal=peri_day
               !eccpal=0.009

             endif

             write(*,*) "Paleo peri=",peri_daypal,"  tpal=",tpalnew
             write(*,*) "Paleo eccpal=",eccpal,"  tpal=",tpalnew

             ! create restartfi
             call physdem1pal("restartfi.nc",long,lati,nsoilmx,nq,
     .            ptimestep,pdaypal,
     .            ztime_fin,tsurf,tsoil,emis,q2,qsurfpal,
     .            area,albedodat,tidat,zmea,zstd,zsig,
     .            zgam,zthe,oblipal,eccpal,tpalnew,adjustnew,phisfipal,
     .            peri_daypal)
            else

               call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
     .              ptimestep,pday,
     .              ztime_fin,tsurf,tsoil,emis,q2,qsurf,
     .              area,albedodat,tidat,zmea,zstd,zsig,
     .              zgam,zthe)

            endif

         ENDIF

!        -----------------------------------------------------------------
!        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 !


         IF (callstats) THEN
             write (*,*) "stats have been turned off in the code.
     &                You can manage your own output in physiq.F "
             stop

!            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,"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)
!             call wstats(ngrid,"fluxtop_sw",
!     .                  "Solar radiative flux to space","W.m-2",2,
!     .                  fluxtop_sw_tot)

!            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)

            IF(lastcall) THEN
              write (*,*) "Writing stats..."
              call mkstats(ierr)
            ENDIF
          ENDIF !if callstats


!        ---------------------------------------------------------------------
!        3D OUTPUTS
!        ----------------------------------------------------------------------
!        output in netcdf file "DIAGFI", containing any variable for diagnostic
!        (output with  period "ecritphy", set in "run.def")
!        ----------------------------------------------------------------------

!          if(MOD(countG3D,saveG3D).eq.0)then
!          print*,'countG3D',countG3D

!!       BASIC

        call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
     &                  tsurf)
        call WRITEDIAGFI(ngrid,"rice_ch4","ch4 ice mass mean radius",
     &         "m",3,rice_ch4)
        call WRITEDIAGFI(ngrid,"area","area","m",2,area)
        call WRITEDIAGFI(ngrid,"mu0","cos sza","m",2,mu0)
        call WRITEDIAGFI(ngrid,"fract","fractional day","",2,fract)
        call WRITEDIAGFI(ngrid,"declin","declin","deg",0,declin*180./pi)
        call WRITEDIAGFI(ngrid,"ls","ls","deg",0,zls*180./pi)
        call WRITEDIAGFI(ngrid,"dist_sol","dist_sol","AU",0,dist_sol)
        call WRITEDIAGFI(ngrid,"phisfi","phisfi"," ",2,phisfi)
!        call WRITEDIAGFI(ngrid,"tsub","subsurface temperature","K",
!     &       1,tsub)
!        call WRITEDIAGFI(ngrid,"phitop","phitop"," ",2,phitop)

!        ENERGY - Total energy balance diagnostics

        call WRITEDIAGFI(ngrid,"albedo","albedo","sans dim",2,
     &                  albedo)
        call WRITEDIAGFI(ngrid,"emissivite","emissivite","sans dim",2,
     &                  emis)
        call WRITEDIAGFI(ngrid,"fluxtop_dn","fluxtop_dn","sans dim",2,
     &                  fluxtop_dn)
        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)

        IF (fast) then
          ! pressure reprocessed in nogcm.F
          call WRITEDIAGFI(ngrid,"globave","surf press","Pa",0,globave)
          call WRITEDIAGFI(ngrid,"ps","surf press","Pa",2,ps)
          call WRITEDIAGFI(ngrid,"fluxrad","fluxrad",
     &                                        "W m-2",2,fluxrad)
          call WRITEDIAGFI(ngrid,"fluxgrd","fluxgrd",
     &                                        "W m-2",2,fluxgrd)
          call WRITEDIAGFI(ngrid,"capcal","capcal",
     &                                        "W.s m-2 K-1",2,capcal)
!          call WRITEDIAGFI(ngrid,"dplanck","dplanck",
!     &                                        "W.s m-2 K-1",2,dplanck)
          call WRITEDIAGFI(ngrid,"tsoil","tsoil","K",3,tsoil)

          ! If we want to see tidat evolution with time:
          call WRITEDIAGFI(ngrid,"tidat","tidat","SI",3,tidat_out)


          IF (fasthaze) then
             call WRITEDIAGFI(ngrid,"gradflux","gradflux",
     &                                        "Ph m-2 s-1",2,gradflux)
             call WRITEDIAGFI(ngrid,"fluxbot","flux bottom",
     &                                        "Ph m-2 s-1",2,fluxbot)
             call WRITEDIAGFI(ngrid,"fluxbotsol","flux bottom sol",
     &                                  "Ph m-2 s-1",2,fluxlym_sol_bot)
             call WRITEDIAGFI(ngrid,"fluxbotipm","flux bottom ipm",
     &                                  "Ph m-2 s-1",2,fluxlym_ipm_bot)
             call WRITEDIAGFI(ngrid,"fluxlymipm","flux incid ipm",
     &                                  "Ph m-2 s-1",2,flym_ipm)
             call WRITEDIAGFI(ngrid,"fluxlymsol","flux incid sol",
     &                                  "Ph m-2 s-1",2,flym_sol)
             call WRITEDIAGFI(ngrid,"tend_prodhaze","prod haze",
     &                                    "kg m-2 s-1",2,zdqsprodhaze)
             call WRITEDIAGFI(ngrid,"tend_lostch4","tend ch4",
     &                     "kg kg-1 s-1",2,zdqprodhaze(1,igcm_ch4_gas))
             call WRITEDIAGFI(ngrid,"tend_prodhaze","tend prod haze",
     &                     "kg kg-1 s-1",2,zdqprodhaze(1,igcm_haze))
          ENDIF
          IF (paleo) then
            call WRITEDIAGFI(ngrid,"qsurfn2_year","qsurfn2_year",
     &                "kg m-2.plutoyear-1",2,qsurfyear(:,1))
            call WRITEDIAGFI(ngrid,"phisfipal","phisfipal",
     &                                        " ",2,phisfipal)
            call WRITEDIAGFI(ngrid,"zdqflow","zdqflow",
     &                                        " ",2,zdqflow(1,igcm_n2))
          ENDIF
        ELSE
          call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
          call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
          call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
          call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
          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,"p","Pression","Pa",3,pplay)
          call WRITEDIAGFI(ngrid,"fluxrad","fluxrad",
     &                                        "W m-2",2,fluxrad)
          call WRITEDIAGFI(ngrid,"fluxgrd","fluxgrd",
     &                                        "W m-2",2,fluxgrd)
        ENDIF

!       Usually not used :

!        call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
        call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
!        call WRITEDIAGFI(ngrid,"DYN","dynamical heat input","W m-2",
!     &       2,fluxdyn)

!        TENDANCIES

        call WRITEDIAGFI(ngrid,"dps","dps","K",2,pdpsrf)
        call WRITEDIAGFI(ngrid,"zdtsw","SW heating","K s-1",3,zdtsw)
        call WRITEDIAGFI(ngrid,"zdtlw","LW heating","K s-1",3,zdtlw)
        call WRITEDIAGFI(ngrid,"zdtc","tendancy T cond N2","K",3,zdtc)
        call WRITEDIAGFI(ngrid,"zdtch4cloud","tendancy T ch4cloud",
     &                                        "K",3,zdtch4cloud)
        call WRITEDIAGFI(ngrid,"zdtcocloud","tendancy T cocloud",
     &                                        "K",3,zdtcocloud)
        call WRITEDIAGFI(ngrid,"zdtsurfc","zdtsurfc","K",2,zdtsurfc)
        call WRITEDIAGFI(ngrid,"zdtdif","zdtdif","K",3,zdtdif)
        call WRITEDIAGFI(ngrid,"zdtconduc","tendancy T conduc",
     &                                                "K",3,zdtconduc)
        call WRITEDIAGFI(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
        call WRITEDIAGFI(ngrid,"zdtrad","rad heating","T s-1",3,dtrad)
        call WRITEDIAGFI(ngrid,"zdtadj","conv adj","T s-1",3,zdtadj)
!        call WRITEDIAGFI(ngrid,"zdqcloud","ch4 cloud","T s-1",
!     &                  3,zdqch4cloud(1,1,igcm_ch4_ice))
!        call WRITEDIAGFI(ngrid,"zdqcloudgas","ch4 cloud","T s-1",
!     &                  3,zdqch4cloud(1,1,igcm_ch4_gas))
!        call WRITEDIAGFI(ngrid,"zdqcn2_c","zdq condn2","",3,zdqc(:,:,3))
!        call WRITEDIAGFI(ngrid,"zdqdif_c","zdqdif","",3,zdqdif(:,:,3))
!        call WRITEDIAGFI(ngrid,"zdqadj_c","zdqadj","",3,zdqadj(:,:,3))
        !call WRITEDIAGFI(ngrid,"zdqsed_h","zdqsed","",3,zdqsed(:,:,7))
        !call WRITEDIAGFI(ngrid,"zdqssed","zdqssed","",2,zdqssed)
        call WRITEDIAGFI(ngrid,"zq1temp_ch4"," "," ",2,zq1temp_ch4)
        call WRITEDIAGFI(ngrid,"qsat_ch4"," "," ",2,qsat_ch4)
        call WRITEDIAGFI(ngrid,"qsat_ch4_l1"," "," ",2,qsat_ch4_l1)
        call WRITEDIAGFI(ngrid,"senshf1","senshf1"," ",2,sensiblehf1)
        call WRITEDIAGFI(ngrid,"senshf2","senshf2"," ",2,sensiblehf2)


!        OTHERS

!      call WRITEDIAGFI(ngrid,"dWtotdifv","dWtotdifv","kg m-2",1,dWtot)
!      call WRITEDIAGFI(ngrid,"dWtotsdifv","dWtotsdifv","kgm-2",1,dWtots)
!      call WRITEDIAGFI(ngrid,"nconsMX","nconsMX","kgm-2s-1",1,nconsMAX)


!       TRACERS

       if (tracer) then

!!!         afficher la tendance sur le flux de la glace
          call WRITEDIAGFI(ngridmx,'n2_iceflux',
     &                    'n2_iceflux',
     &                      "kg m^-2 s^-1",2,flusurf(1,igcm_n2) )

        do iq=1,nq
          call WRITEDIAGFI(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
          call WRITEDIAGFI(ngrid,trim(noms(iq))//'_Tcondn2',
     &                           trim(noms(iq))//'_Tcondn2',
     &                      'kg/kg',3,zdqc(1,1,iq))

          call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_surf',
     &                    trim(noms(iq))//'_surf',
     &                    'kg m^-2',2,qsurf(1,iq) )

          call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_col',
     &                    trim(noms(iq))//'_col',
     &                    'kg m^-2',2,qcol(1,iq) )

!         call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_aero',
!     &                    trim(noms(iq))//'_aero',
!     &                    'kg/kg',3,aerosol(1,1,iq))


        enddo

         call WRITEDIAGFI(ngridmx,'haze_reff',
     &                    'haze_reff',
     &                    'm',3,reffrad(1,1,1))

        if (methane) then

         call WRITEDIAGFI(ngridmx,'ch4_iceflux',
     &                    'ch4_iceflux',
     &                      "kg m^-2 s^-1",2,flusurf(1,igcm_ch4_ice) )

         call WRITEDIAGFI(ngrid,"vmr_ch4","vmr_ch4","%",
     &                                            2,vmr_ch4)
         if (.not.fast) THEN
          call WRITEDIAGFI(ngrid,"zrho_ch4","zrho_ch4","kg.m-3",
     &                                            3,zrho_ch4(:,:))
         endif

!        if (fast) THEN
!             call WRITEDIAGFI(ngrid,"zq1_ch4","zq ch4","kg m-2",
!    &                                      2,zq(1,1,igcm_ch4_gas))
!        endif

         ! Tendancies
         call WRITEDIAGFI(ngrid,"zdqch4cloud","ch4 cloud","T s-1",
     &                  3,zdqch4cloud(1,1,igcm_ch4_gas))
         call WRITEDIAGFI(ngrid,"zdqcn2_ch4","zdq condn2 ch4","",
     &                   3,zdqc(:,:,igcm_ch4_gas))
         call WRITEDIAGFI(ngrid,"zdqdif_ch4","zdqdif ch4","",
     &                   3,zdqdif(:,:,igcm_ch4_gas))
         call WRITEDIAGFI(ngrid,"zdqsdif_ch4_ice","zdqsdif ch4","",
     &                   2,zdqsdif(:,igcm_ch4_ice))
         call WRITEDIAGFI(ngrid,"zdqadj_ch4","zdqadj ch4","",
     &                   3,zdqadj(:,:,igcm_ch4_gas))
         call WRITEDIAGFI(ngrid,"zdqsed_ch4","zdqsed ch4","",
     &                   3,zdqsed(:,:,igcm_ch4_gas))
         call WRITEDIAGFI(ngrid,"zdqssed_ch4","zdqssed ch4","",
     &                   2,zdqssed(:,igcm_ch4_gas))
        call WRITEDIAGFI(ngrid,"zdtch4cloud","ch4 cloud","T s-1",
     &                  3,zdtch4cloud)

        endif

        if (carbox) then

         call WRITEDIAGFI(ngridmx,'co_iceflux',
     &                    'co_iceflux',
     &                      "kg m^-2 s^-1",2,flusurf(1,igcm_co_ice) )

         call WRITEDIAGFI(ngrid,"vmr_co","vmr_co","%",
     &                                            2,vmr_co)

         if (.not.fast) THEN
          call WRITEDIAGFI(ngrid,"zrho_co","zrho_co","kg.m-3",
     &                                            3,zrho_co(:,:))
         endif

!        if (fast) THEN
!        call WRITEDIAGFI(ngrid,"zq1_co","zq co","kg m-2",
!    &                                      2,zq(1,1,igcm_co_gas))
!        endif

         ! Tendancies
!        call WRITEDIAGFI(ngrid,"zdqcocloud","co cloud","T s-1",
!    &                  3,zdqcocloud(1,1,igcm_co_gas))
!        call WRITEDIAGFI(ngrid,"zdqcn2_co","zdq condn2 co","",
!    &                   3,zdqc(:,:,igcm_co_gas))
!        call WRITEDIAGFI(ngrid,"zdqdif_co","zdqdif co","",
!    &                   3,zdqdif(:,:,igcm_co_gas))
!        call WRITEDIAGFI(ngrid,"zdqsdif_co_ice","zdqsdif co","",
!    &                   2,zdqsdif(:,igcm_co_ice))
!        call WRITEDIAGFI(ngrid,"zdqadj_co","zdqadj co","",
!    &                   3,zdqadj(:,:,igcm_co_gas))
!        call WRITEDIAGFI(ngrid,"zdqsed_co","zdqsed co","",
!    &                   3,zdqsed(:,:,igcm_co_gas))
!        call WRITEDIAGFI(ngrid,"zdqssed_co","zdqssed co","",
!    &                   2,zdqssed(:,igcm_co_gas))
!       call WRITEDIAGFI(ngrid,"zdtcocloud","co cloud","T s-1",
!    &                  3,zdtcocloud)

       endif

        if (haze) then

!         call WRITEDIAGFI(ngrid,"zrho_haze","zrho_haze","kg.m-3",
!     &                                            3,zrho_haze(:,:))
         call WRITEDIAGFI(ngrid,"zdqrho_photprec","zdqrho_photprec",
     &                  "kg.m-3.s-1",3,zdqrho_photprec(:,:))
         call WRITEDIAGFI(ngrid,"zdqphot_prec","zdqphot_prec","",
     &                                         3,zdqphot_prec(:,:))
         call WRITEDIAGFI(ngrid,"zdqhaze_ch4","zdqhaze_ch4","",
     &             3,zdqhaze(:,:,igcm_ch4_gas))
         call WRITEDIAGFI(ngrid,"zdqhaze_prec","zdqhaze_prec","",
     &             3,zdqhaze(:,:,igcm_prec_haze))
         if (igcm_haze.ne.0) then
           call WRITEDIAGFI(ngrid,"zdqhaze_haze","zdqhaze_haze","",
     &               3,zdqhaze(:,:,igcm_haze))
           call WRITEDIAGFI(ngrid,"zdqssed_haze","zdqssed haze",
     &          "kg/m2/s",2,zdqssed(:,igcm_haze))
         endif
         call WRITEDIAGFI(ngrid,"zdqphot_ch4","zdqphot_ch4","",
     &                                         3,zdqphot_ch4(:,:))
         call WRITEDIAGFI(ngrid,"zdqconv_prec","zdqconv_prec","",
     &                                        3,zdqconv_prec(:,:))
!         call WRITEDIAGFI(ngrid,"zdqhaze_col","zdqhaze col","kg/m2/s",
!     &                   2,zdqhaze_col(:))
        endif

        if (aerohaze) then
         call WRITEDIAGFI(ngridmx,"tau_col",
     &         "Total aerosol optical depth","opacity",2,tau_col)
        endif

!         call WRITEDIAGFI(ngridmx,"tau_col",
!     &         "Total aerosol optical depth","[]",2,tau_col)
       endif

!       countG3D=1
!       else
!            countG3D=countG3D+1
!       endif ! if time to save

!      ----------------------------------------------------------------------
!      1D OUTPUTS
!      ----------------------------------------------------------------------
      ELSE         ! if(ngrid.eq.1)

       if(countG1D.eq.saveG1D)then

!      BASIC 1D ONLY

          call WRITEDIAGFI(ngrid,"ISR","incoming stellar rad.","W m-2",
     &       0,fluxtop_dn)
          call WRITEDIAGFI(ngrid,"ASR","absorbed stellar rad.","W m-2",
     &       0,fluxabs_sw)
          call WRITEDIAGFI(ngrid,"OLR","outgoing longwave rad.","W m-2",
     &       0,fluxtop_lw)
          call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
     &                  tsurf)
          call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",0,ps)
          call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
          call WRITEDIAGFI(ngrid,"p","Pression","Pa",3,pplay)

          call WRITEDIAGFI(ngrid,"fluxsurf_sw","sw surface flux",
     &                                         "K",0,fluxsurf_sw)
          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,"zdtdif","zdtdif","K",3,zdtdif)
          call WRITEDIAGFI(ngrid,"zdtconduc","tendancy T conduc",
     &                                                "K",3,zdtconduc)

          call WRITEDIAGFI(ngrid,"fluxsurf_sw","sw surface flux",
     &                                         "K",0,fluxsurf_sw)
          call WRITEDIAGFI(ngrid,"declin","lat subsolaire",
     &                                         "deg",0,declin*180./pi)
          call WRITEDIAGFI(ngrid,"ls","ls","deg",0,zls*180./pi)
          call WRITEDIAGFI(ngrid,"dist_sol","dist_sol","AU",0,dist_sol)
          call WRITEDIAGFI(ngrid,"mu0","solar zenith angle",
     &                                         "deg",0,mu0)
          call WRITEDIAGFI(ngrid,"albedo","albedo",
     &                                         "",0,albedo)
          call WRITEDIAGFI(ngrid,"emis","emis",
     &                                         "",0,emis)
          call WRITEDIAGFI(ngrid,"tsoil","tsoil","K",1,tsoil(1,:))
          call WRITEDIAGFI(ngrid,"tidat","tidat","SI",1,tidat_out(1,:))

!       OUTPUT OF TENDANCIES
!          call WRITEDIAGFI(ngrid,"zdqcloud","ch4 cloud","T s-1",
!     &                  3,zdqcloud(1,1,igcm_ch4_gas))
!          call WRITEDIAGFI(ngrid,"zdqcn2","zdq condn2","",1,zdqc(:,:,3))
!          call WRITEDIAGFI(ngrid,"zdqdif","zdqdif","",1,zdqdif(:,:,3))
!          call WRITEDIAGFI(ngrid,"zdqadj","zdqadj","",1,zdqadj(:,:,3))
!          call WRITEDIAGFI(ngrid,"zdqsed","zdqsed","",1,zdqsed(:,:,3))
!          call WRITEDIAGFI(ngrid,"zdqc","zdqc","",1,zdqc(:,:,3))
!          call WRITEDIAGFI(ngrid,"zdqhaze","zdqhaze","",1,
!     &               zdqhaze(:,:,3))

        ! 1D methane
        if (methane) then

         call WRITEDIAGFI(ngridmx,'ch4_iceflux',
     &                    'ch4_iceflux',
     &                      "kg m^-2 s^-1",0,flusurf(1,igcm_ch4_ice) )

         call WRITEDIAGFI(ngrid,"vmr_ch4","vmr_ch4","%",
     &                                            0,vmr_ch4(1))
         call WRITEDIAGFI(ngrid,"zrho_ch4","zrho_ch4","kg.m-3",
     &                                            1,zrho_ch4(1,:))

         ! Tendancies
         call WRITEDIAGFI(ngrid,"zdqch4cloud","ch4 cloud","T s-1",
     &                  1,zdqch4cloud(1,1,igcm_ch4_gas))
         call WRITEDIAGFI(ngrid,"zdqcn2_ch4","zdq condn2 ch4","",
     &                   1,zdqc(1,:,igcm_ch4_gas))
         call WRITEDIAGFI(ngrid,"zdqdif_ch4","zdqdif ch4","",
     &                   1,zdqdif(1,:,igcm_ch4_gas))
         call WRITEDIAGFI(ngrid,"zdqsdif_ch4_ice","zdqsdif ch4","",
     &                   1,zdqsdif(:,igcm_ch4_ice))
         call WRITEDIAGFI(ngrid,"zdqadj_ch4","zdqadj ch4","",
     &                   1,zdqadj(1,:,igcm_ch4_gas))
         call WRITEDIAGFI(ngrid,"zdqsed_ch4","zdqsed ch4","",
     &                   1,zdqsed(1,:,igcm_ch4_gas))
         call WRITEDIAGFI(ngrid,"zdqssed_ch4","zdqssed ch4","",
     &                   0,zdqssed(1,igcm_ch4_gas))

        endif

        ! 1D Haze
        if (haze) then

!         call WRITEDIAGFI(ngrid,"zrho_haze","zrho_haze","kg.m-3",
!     &                                            1,zrho_haze(:,:))
         call WRITEDIAGFI(ngrid,"zdqrho_photprec","zdqrho_photprec",
     &                  "kg.m-3.s-1",3,zdqrho_photprec)
         call WRITEDIAGFI(ngrid,"zdqphot_prec","zdqphot_prec","",
     &                                            3,zdqphot_prec)
         call WRITEDIAGFI(ngrid,"zdqhaze_ch4","zdqhaze_ch4","",
     &             1,zdqhaze(1,:,igcm_ch4_gas))
         call WRITEDIAGFI(ngrid,"zdqhaze_prec","zdqhaze_prec","",
     &             1,zdqhaze(1,:,igcm_prec_haze))
         call WRITEDIAGFI(ngrid,"zdqhaze_haze","zdqhaze_haze","",
     &             1,zdqhaze(1,:,igcm_haze))
         call WRITEDIAGFI(ngrid,"zdqphot_ch4","zdqphot_ch4","",
     &                                         3,zdqphot_ch4)
!         call WRITEDIAGFI(ngrid,"zdqconv_prec","zdqconv_prec","",
!     &                                         1,zdqconv_prec(1,:))
         call WRITEDIAGFI(ngrid,"zdqssed_haze","zdqssed haze","kg/m2/s",
     &                   0,zdqssed(1,igcm_haze))
!         call WRITEDIAGFI(ngrid,"zdqhaze_col","zdqhaze col","kg/m2/s",
!     &                   0,zdqhaze_col(:))
         if (haze_radproffix) then
           call WRITEDIAGFI(ngridmx,'haze_reffrad','haze_reffrad','m',
     &                   3,reffrad(1,1,1))
         endif
        endif

        if (aerohaze) then
         call WRITEDIAGFI(ngridmx,"tau_col",
     &         "Total aerosol optical depth","[]",0,tau_col(1))
         call WRITEDIAGFI(ngridmx,"tau_aero",
     &         "aerosol optical depth","[]",3,aerosol(1,1,1))
        endif

!      TRACERS 1D ONLY
        if (tracer) then
         do iq=1,nq
          call WRITEDIAGFI(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
          call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_surf',
     &                    trim(noms(iq))//'_surf',
     &                    'kg m^-2',2,qsurf(1,iq) )
          call WRITEDIAGFI(ngridmx,trim(noms(iq))//'_col',
     &                    trim(noms(iq))//'_col',
     &                    'kg m^-2',2,qcol(1,iq) )
         enddo
        endif

         countG1D=1
         else
            countG1D=countG1D+1
         endif ! if time to save

      END IF       ! if(ngrid.ne.1)

!      countG3D=countG3D+1
      icount=icount+1

      RETURN
      END