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 tracer_mod, only: noms, mmol, igcm_co2, igcm_n2,
     &                      igcm_co, igcm_o, igcm_h2o_vap, igcm_h2o_ice,
     &                      igcm_ccn_mass, igcm_ccn_number,
     &                      igcm_dust_mass, igcm_dust_number, igcm_h2o2,
     &                      nuice_ref, rho_ice, rho_dust, ref_r0
      use comsoil_h, only: inertiedat, ! soil thermal inertia
     &                     tsoil, nsoilmx ! number of subsurface layers
      use comgeomfi_h, only: long, lati, area,
     &                       sinlon, coslon, sinlat, coslat
      use surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam,
     &                     zthe, z0, albedo_h2o_ice,
     &                     frost_albedo_threshold,
     &                     tsurf, co2ice, emis,
     &                     capcal, fluxgrd, qsurf
      use comsaison_h, only: dist_sol, declin, mu0, fract
      use slope_mod, only: theta_sl, psi_sl
      use conc_mod, only: rnew, cpnew, mmean
      use control_mod, only: iphysiq, day_step, ecritstart
      use dimradmars_mod, only: tauscaling, aerosol,
     &                          dtrad, fluxrad_sky, fluxrad, albedo,
     &                          naerkind
      use turb_mod, only: q2, wstar, ustar, sensibFlux, 
     &                    zmax_th, hfmax_th, turb_resolved
      use planete_h, only: aphelie, periheli, year_day, peri_day,
     &                     obliquit
      USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad, daysec
      use param_v4_h, only: nreact,n_avog, 
     &                      fill_data_thermos, allocate_param_thermos
      use iono_h, only: allocate_param_iono
#ifdef MESOSCALE
      use comsoil_h, only: mlayer,layer
      use surfdat_h, only: z0_default
      use comm_wrf
#else
      use phyredem, only: physdem0, physdem1
      use eofdump_mod, only: eofdump
#endif

      IMPLICIT NONE
c=======================================================================
c
c   subject:
c   --------
c
c   Organisation of the physical parametrisations of the LMD 
c   martian atmospheric general circulation model.
c
c   The GCM can be run without or with tracer transport
c   depending on the value of Logical "tracer" in file  "callphys.def"
c   Tracers may be water vapor, ice OR chemical species OR dust particles
c
c   SEE comments in initracer.F about numbering of tracer species...
c
c   It includes:
c
c      1. Initialization:
c      1.1 First call initializations
c      1.2 Initialization for every call to physiq
c      1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
c      2. Compute radiative transfer tendencies
c         (longwave and shortwave) for CO2 and aerosols.
c      3. Gravity wave and subgrid scale topography drag :
c      4. Vertical diffusion (turbulent mixing):
c      5. Convective adjustment
c      6. Condensation and sublimation of carbon dioxide.
c      7.  TRACERS :
c       7a. water and water ice
c       7b. call for photochemistry when tracers are chemical species
c       7c. other scheme for tracer (dust) transport (lifting, sedimentation)
c       7d. updates (CO2 pressure variations, surface budget)
c      8. Contribution to tendencies due to thermosphere
c      9. Surface and sub-surface temperature calculations
c     10. Write outputs :
c           - "startfi", "histfi" (if it's time)
c           - Saving statistics (if "callstats = .true.")
c           - Dumping eof (if "calleofdump = .true.")
c           - Output any needed variables in "diagfi" 
c     11. Diagnostic: mass conservation of tracers
c 
c   author: 
c   ------- 
c           Frederic Hourdin	15/10/93
c           Francois Forget		1994
c           Christophe Hourdin	02/1997 
c           Subroutine completly rewritten by F.Forget (01/2000)
c           Introduction of the photochemical module: S. Lebonnois (11/2002)
c           Introduction of the thermosphere module: M. Angelats i Coll (2002)
c           Water ice clouds: Franck Montmessin (update 06/2003)
c           Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
c             Nb: See callradite.F for more information.
c           Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags
c           jul 2011 malv+fgg: Modified calls to NIR heating routine and 15 um cooling parameterization
c           
c   arguments:
c   ----------
c
c   input:
c   ------
c    ecri                  period (in dynamical timestep) to write output
c    ngrid                 Size of the horizontal grid.
c                          All internal loops are performed on that grid.
c    nlayer                Number of vertical layers.
c    nq                    Number of advected fields
c    firstcall             True at the first call
c    lastcall              True at the last call
c    pday                  Number of days counted from the North. Spring
c                          equinoxe.
c    ptime                 Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
c    ptimestep             timestep (s)
c    pplay(ngrid,nlayer)   Pressure at the middle of the layers (Pa)
c    pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
c    pphi(ngrid,nlayer)    Geopotential at the middle of the layers (m2s-2)
c    pu(ngrid,nlayer)      u component of the wind (ms-1)
c    pv(ngrid,nlayer)      v component of the wind (ms-1)
c    pt(ngrid,nlayer)      Temperature (K)
c    pq(ngrid,nlayer,nq)   Advected fields
c    pudyn(ngrid,nlayer)    | 
c    pvdyn(ngrid,nlayer)    | Dynamical temporal derivative for the
c    ptdyn(ngrid,nlayer)    | corresponding variables
c    pqdyn(ngrid,nlayer,nq) |
c    pw(ngrid,?)           vertical velocity
c
c   output:
c   -------
c
c    pdu(ngrid,nlayer)       |
c    pdv(ngrid,nlayer)       |  Temporal derivative of the corresponding
c    pdt(ngrid,nlayer)       |  variables due to physical processes.
c    pdq(ngrid,nlayer,nq)    |
c    pdpsrf(ngrid)           |
c    tracerdyn                 call tracer in dynamical part of GCM ?

c
c=======================================================================
c
c    0.  Declarations :
c    ------------------

#include "callkeys.h"
#include "comg1d.h"
#include "nlteparams.h"
#include "comvert.h"
#include "chimiedata.h"
#include "netcdf.inc"

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

c   inputs:
c   -------
      INTEGER ngrid,nlayer,nq
      REAL ptimestep
      REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
      REAL pphi(ngrid,nlayer)
      REAL pu(ngrid,nlayer),pv(ngrid,nlayer)
      REAL pt(ngrid,nlayer),pq(ngrid,nlayer,nq)
      REAL pw(ngrid,nlayer) !Mars pvervel transmit par dyn3d
      REAL zh(ngrid,nlayer)      ! potential temperature (K)
      LOGICAL firstcall,lastcall

      REAL pday
      REAL ptime
      logical tracerdyn

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


c Local saved variables:
c ----------------------
      INTEGER,SAVE :: day_ini ! Initial date of the run (sol since Ls=0)
      INTEGER,SAVE :: icount     ! counter of calls to physiq during the run.
     
c     Variables used by the water ice microphysical scheme:
      REAL rice(ngrid,nlayer)    ! Water ice geometric mean radius (m)
      REAL nuice(ngrid,nlayer)   ! Estimated effective variance
                                     !   of the size distribution
      real rsedcloud(ngrid,nlayer) ! Cloud sedimentation radius
      real rhocloud(ngrid,nlayer)  ! Cloud density (kg.m-3)
      REAL surfdust(ngrid,nlayer) ! dust surface area (m2/m3, if photochemistry)
      REAL surfice(ngrid,nlayer)  !  ice surface area (m2/m3, if photochemistry)
      REAL inertiesoil(ngrid,nsoilmx) ! Time varying subsurface
                                      ! thermal inertia (J.s-1/2.m-2.K-1)
                                      ! (used only when tifeedback=.true.)

c     Variables used by the slope model
      REAL sl_ls, sl_lct, sl_lat
      REAL sl_tau, sl_alb, sl_the, sl_psi
      REAL sl_fl0, sl_flu
      REAL sl_ra, sl_di0
      REAL sky

      REAL,PARAMETER :: stephan = 5.67e-08 ! Stephan Boltzman constant

c Local variables :
c -----------------

      REAL CBRT
      EXTERNAL CBRT

!      CHARACTER*80 fichier 
      INTEGER l,ig,ierr,igout,iq,tapphys

      REAL fluxsurf_lw(ngrid)      !incident LW (IR) surface flux (W.m-2)
      REAL fluxsurf_sw(ngrid,2)    !incident SW (solar) surface flux (W.m-2)
      REAL fluxtop_lw(ngrid)       !Outgoing LW (IR) flux to space (W.m-2)
      REAL fluxtop_sw(ngrid,2)     !Outgoing SW (solar) flux to space (W.m-2)
      REAL tauref(ngrid)           ! Reference column optical depth at odpref
      real,parameter :: odpref=610. ! DOD reference pressure (Pa)
      REAL tau(ngrid,naerkind)     ! Column dust optical depth at each point
                                   ! AS: TBD: this one should be in a module !
      REAL zls                       !  solar 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)

c     Tendancies due to various processes:
      REAL dqsurf(ngrid,nq)
      REAL zdtlw(ngrid,nlayer)     ! (K/s)
      REAL zdtsw(ngrid,nlayer)     ! (K/s)
!      REAL cldtlw(ngrid,nlayer)     ! (K/s) LW heating rate for clear area
!      REAL cldtsw(ngrid,nlayer)     ! (K/s) SW heating rate for clear area
      REAL zdtnirco2(ngrid,nlayer) ! (K/s)
      REAL zdtnlte(ngrid,nlayer)   ! (K/s)
      REAL zdtsurf(ngrid)            ! (K/s)
      REAL zdtcloud(ngrid,nlayer)
      REAL zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer)  ! (m.s-2)
      REAL zdhdif(ngrid,nlayer), zdtsdif(ngrid)         ! (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 zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer)    ! (m.s-2)
      REAL zdtc(ngrid,nlayer),zdtsurfc(ngrid)
      REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer)

      REAL zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq)
      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 zdqcloud(ngrid,nlayer,nq)
      REAL zdqscloud(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)

c     Local variable for local intermediate calcul:
      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
      REAL zdh(ngrid,nlayer)
      INTEGER length
      PARAMETER (length=100)

c local variables only used for diagnostic (output in file "diagfi" or "stats")
c -----------------------------------------------------------------------------
      REAL ps(ngrid), zt(ngrid,nlayer)
      REAL zu(ngrid,nlayer),zv(ngrid,nlayer)
      REAL zq(ngrid,nlayer,nq)
      REAL fluxtop_sw_tot(ngrid), fluxsurf_sw_tot(ngrid)
      character*2 str2
!      character*5 str5
      real zdtdif(ngrid,nlayer), zdtadj(ngrid,nlayer)
      real rdust(ngrid,nlayer) ! dust geometric mean radius (m)
      integer igmin, lmin
      logical tdiag

      real co2col(ngrid)        ! CO2 column
      ! pplev and pplay are dynamical inputs and must not be modified in the physics.
      ! instead, use zplay and zplev :
      REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer) 
!      REAL zstress(ngrid),cd
      real tmean, zlocal(nlayer)
      real rho(ngrid,nlayer)  ! density
      real vmr(ngrid,nlayer)  ! volume mixing ratio
      real rhopart(ngrid,nlayer) ! number density of a given species
      real colden(ngrid,nq)     ! vertical column of tracers
      REAL mtot(ngrid)          ! Total mass of water vapor (kg/m2)
      REAL icetot(ngrid)        ! Total mass of water ice (kg/m2)
      REAL Nccntot(ngrid)       ! Total number of ccn (nbr/m2)
      REAL Mccntot(ngrid)       ! Total mass of ccn (kg/m2)
      REAL rave(ngrid)          ! Mean water ice effective radius (m)
      REAL opTES(ngrid,nlayer)  ! abs optical depth at 825 cm-1
      REAL tauTES(ngrid)        ! column optical depth at 825 cm-1
      REAL Qabsice                ! Water ice absorption coefficient
      REAL taucloudtes(ngrid)  ! Cloud opacity at infrared
                               !   reference wavelength using
                               !   Qabs instead of Qext
                               !   (direct comparison with TES)
                               
      REAL dqdustsurf(ngrid) ! surface q dust flux (kg/m2/s)
      REAL dndustsurf(ngrid) ! surface n dust flux (number/m2/s)
      REAL ndust(ngrid,nlayer) ! true n dust (kg/kg)
      REAL qdust(ngrid,nlayer) ! true q dust (kg/kg)
      REAL nccn(ngrid,nlayer)  ! true n ccn (kg/kg)
      REAL qccn(ngrid,nlayer)  ! true q ccn (kg/kg)

c Test 1d/3d scavenging
      real h2otot(ngrid)
      REAL satu(ngrid,nlayer)  ! satu ratio for output
      REAL zqsat(ngrid,nlayer) ! saturation

      REAL,SAVE :: time_phys

! Added for new NLTE scheme

      real co2vmr_gcm(ngrid,nlayer)
      real n2vmr_gcm(ngrid,nlayer)
      real ovmr_gcm(ngrid,nlayer)
      real covmr_gcm(ngrid,nlayer)
      integer ierr_nlte
      real*8  varerr

c Variables for PBL
      REAL zz1(ngrid)
      REAL lmax_th_out(ngrid)
      REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
      REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
      INTEGER lmax_th(ngrid),dimout,n_out,n
      CHARACTER(50) zstring
      REAL dtke_th(ngrid,nlayer+1)
      REAL zcdv(ngrid), zcdh(ngrid)
      REAL, ALLOCATABLE, DIMENSION(:,:) :: T_out
      REAL, ALLOCATABLE, DIMENSION(:,:) :: u_out ! Interpolated teta and u at z_out
      REAL u_out1(ngrid)
      REAL T_out1(ngrid)
      REAL, ALLOCATABLE, DIMENSION(:) :: z_out     ! height of interpolation between z0 and z1 [meters]
      REAL tstar(ngrid)  ! friction velocity and friction potential temp
      REAL L_mo(ngrid),vhf(ngrid),vvv(ngrid)
!      REAL zu2(ngrid)

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

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

c  1.1   Initialisation only at first call
c  ---------------------------------------
      IF (firstcall) THEN

c        variables set to 0
c        ~~~~~~~~~~~~~~~~~~
         aerosol(:,:,:)=0
         dtrad(:,:)=0

#ifndef MESOSCALE
         fluxrad(:)=0
         wstar(:)=0.
#endif

c        read startfi 
c        ~~~~~~~~~~~~
#ifndef MESOSCALE
! GCM. Read netcdf initial physical parameters.
         CALL phyetat0 ("startfi.nc",0,0,
     &         nsoilmx,ngrid,nlayer,nq,
     &         day_ini,time_phys,
     &         tsurf,tsoil,emis,q2,qsurf,co2ice,tauscaling)

         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

#else
! MESOSCALE. Supposedly everything is already set in modules.
! So we just check. And we calculate day_ini + two param in dyn3d/control_mod.
      print*,"check: rad,cpp,g,r,rcp,daysec"
      print*,rad,cpp,g,r,rcp,daysec
      PRINT*,'check: tsurf ',tsurf(1),tsurf(ngrid)
      PRINT*,'check: tsoil ',tsoil(1,1),tsoil(ngrid,nsoilmx)
      PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx)
      PRINT*,'check: midlayer,layer ', mlayer(:),layer(:)
      PRINT*,'check: tracernames ', noms
      PRINT*,'check: emis ',emis(1),emis(ngrid)
      PRINT*,'check: q2 ',q2(1,1),q2(ngrid,nlayer+1)
      PRINT*,'check: qsurf ',qsurf(1,1),qsurf(ngrid,nq)
      PRINT*,'check: co2 ',co2ice(1),co2ice(ngrid)
      day_step=daysec/ptimestep
      PRINT*,'Call to LMD physics:',day_step,' per Martian day'
      iphysiq=ptimestep
      day_ini = pday
#endif

c        initialize tracers
c        ~~~~~~~~~~~~~~~~~~
         tracerdyn=tracer
         IF (tracer) THEN
            CALL initracer(ngrid,nq,qsurf)
         ENDIF  ! end tracer

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

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

#ifndef MESOSCALE
c        Initialize thermospheric parameters
c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

         if (callthermos) then
            call fill_data_thermos
            call allocate_param_thermos(nlayer)
            call allocate_param_iono(nlayer,nreact)
            if(solvarmod.eq.0) call param_read
            if(solvarmod.eq.1) call param_read_e107
         endif
#endif
c        Initialize R and Cp as constant

         if (.not.callthermos .and. .not.photochem) then
                 do l=1,nlayer
                  do ig=1,ngrid
                   rnew(ig,l)=r
                   cpnew(ig,l)=cpp
                   mmean(ig,l)=mugaz
                   enddo
                  enddo  
         endif         

         if(callnlte.and.nltemodel.eq.2) call nlte_setup
         if(callnirco2.and.nircorr.eq.1) call NIR_leedat
         if(thermochem) call chemthermos_readini

        IF (tracer.AND.water.AND.(ngrid.NE.1)) THEN
          write(*,*)"physiq: water_param Surface water ice albedo:", 
     .                  albedo_h2o_ice
        ENDIF

#ifndef MESOSCALE
         if (callslope) call getslopes(ngrid,phisfi)

         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,time_phys,area,
     .                 albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
         endif
#endif
                   
      ENDIF        !  (end of "if firstcall")


c ---------------------------------------------------
c 1.2   Initializations done at every physical timestep:
c ---------------------------------------------------
c

c     Initialize various variables
c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      pdv(:,:)=0
      pdu(:,:)=0
      pdt(:,:)=0
      pdq(:,:,:)=0
      pdpsrf(:)=0
      zflubid(:)=0
      zdtsurf(:)=0
      dqsurf(:,:)=0
      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     Initialize pressure levels
c     ~~~~~~~~~~~~~~~~~~
      zplev(:,:) = pplev(:,:)
      zplay(:,:) = pplay(:,:)
      ps(:) = pplev(:,1)


c     Compute geopotential at interlayers
c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c     ponderation des altitudes au niveau des couches en dp/p

      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=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
            z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(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
c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DO l=1,nlayer
         DO ig=1,ngrid 
            zpopsk(ig,l)=(zplay(ig,l)/zplev(ig,1))**rcp
            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
         ENDDO
      ENDDO

#ifndef MESOSCALE
c-----------------------------------------------------------------------
c    1.2.5 Compute mean mass, cp, and R
c    --------------------------------

      if(photochem.or.callthermos) then
         call concentrations(ngrid,nlayer,nq,
     &                       zplay,pt,pdt,pq,pdq,ptimestep)
      endif
#endif
c-----------------------------------------------------------------------
c    2. Compute radiative tendencies :
c------------------------------------


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

c          Local Solar zenith angle
c          ~~~~~~~~~~~~~~~~~~~~~~~~
           CALL orbite(zls,dist_sol,declin)

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

c          NLTE cooling from CO2 emission
c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
           IF(callnlte) then
              if(nltemodel.eq.0.or.nltemodel.eq.1) then
                 CALL nltecool(ngrid,nlayer,nq,zplay,pt,pq,zdtnlte)
              else if(nltemodel.eq.2) then
                co2vmr_gcm(1:ngrid,1:nlayer)=
     &                      pq(1:ngrid,1:nlayer,igcm_co2)*
     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_co2)
                n2vmr_gcm(1:ngrid,1:nlayer)=
     &                      pq(1:ngrid,1:nlayer,igcm_n2)*
     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_n2)
                covmr_gcm(1:ngrid,1:nlayer)=
     &                      pq(1:ngrid,1:nlayer,igcm_co)*
     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_co)
                ovmr_gcm(1:ngrid,1:nlayer)=
     &                      pq(1:ngrid,1:nlayer,igcm_o)*
     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_o)
                 
                 CALL nlte_tcool(ngrid,nlayer,zplay*9.869e-6,
     $                pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm, 
     $                ovmr_gcm,  zdtnlte,ierr_nlte,varerr )
                 if(ierr_nlte.gt.0) then
                    write(*,*)
     $                'WARNING: nlte_tcool output with error message',
     $                'ierr_nlte=',ierr_nlte,'varerr=',varerr
                    write(*,*)'I will continue anyway'
                 endif

                 zdtnlte(1:ngrid,1:nlayer)=
     &                             zdtnlte(1:ngrid,1:nlayer)/86400.
              endif
           else
              zdtnlte(:,:)=0.
           endif

c          Find number of layers for LTE radiation calculations
           IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
     &          CALL nlthermeq(ngrid,nlayer,zplev,zplay)

c          Note: Dustopacity.F has been transferred to callradite.F
         
c          Call main radiative transfer scheme
c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c          Transfer through CO2 (except NIR CO2 absorption)
c            and aerosols (dust and water ice)

c          Radiative transfer
c          ------------------
           CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
     $     emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
     $     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
     $     tauref,tau,aerosol,tauscaling,taucloudtes,rdust,rice,
     $     nuice,co2ice)

c          Outputs for basic check (middle of domain)
c          ------------------------------------------
           write(*,'("Ls =",f11.6," check lat =",f10.6,
     &               " lon =",f11.6)')
     &           zls*180./pi,lati(igout)*180/pi,long(igout)*180/pi
           write(*,'(" tauref(",f4.0," Pa) =",f9.6,
     &             " tau(",f4.0," Pa) =",f9.6)')
     &            odpref,tauref(igout),
     &            odpref,tau(igout,1)*odpref/zplev(igout,1)
c          ---------------------------------------------------------
c          Call slope parameterization for direct and scattered flux
c          ---------------------------------------------------------
           IF(callslope) THEN
            print *, 'Slope scheme is on and computing...'
            DO ig=1,ngrid  
              sl_the = theta_sl(ig)
              IF (sl_the .ne. 0.) THEN
                ztim1=fluxsurf_sw(ig,1)+fluxsurf_sw(ig,2)
                DO l=1,2
                 sl_lct = ptime*24. + 180.*long(ig)/pi/15.
                 sl_ra  = pi*(1.0-sl_lct/12.)
                 sl_lat = 180.*lati(ig)/pi
                 sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
                 sl_alb = albedo(ig,l)
                 sl_psi = psi_sl(ig)
                 sl_fl0 = fluxsurf_sw(ig,l)
                 sl_di0 = 0.
                 if (mu0(ig) .gt. 0.) then
                  sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
                  sl_di0 = sl_di0*1370./dist_sol/dist_sol
                  sl_di0 = sl_di0/ztim1
                  sl_di0 = fluxsurf_sw(ig,l)*sl_di0
                 endif
                 ! you never know (roundup concern...)
                 if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
                 !!!!!!!!!!!!!!!!!!!!!!!!!!
                 CALL param_slope( mu0(ig), declin, sl_ra, sl_lat, 
     &                             sl_tau, sl_alb, sl_the, sl_psi,
     &                             sl_di0, sl_fl0, sl_flu )
                 !!!!!!!!!!!!!!!!!!!!!!!!!!
                 fluxsurf_sw(ig,l) = sl_flu
                ENDDO
              !!! compute correction on IR flux as well
              sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
              fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky
              ENDIF
            ENDDO
           ENDIF

c          CO2 near infrared absorption
c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
           zdtnirco2(:,:)=0
           if (callnirco2) then
              call nirco2abs (ngrid,nlayer,zplay,dist_sol,nq,pq,
     .                       mu0,fract,declin, zdtnirco2)
           endif

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)*(1.-albedo(ig,1))
     $         +fluxsurf_sw(ig,2)*(1.-albedo(ig,2))
           ENDDO


c          Net atmospheric radiative heating rate (K.s-1)
c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
           IF(callnlte) THEN
              CALL blendrad(ngrid, nlayer, zplay,
     &             zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
           ELSE
              DO l=1,nlayer
                 DO ig=1,ngrid
                    dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
     &                          +zdtnirco2(ig,l)
                  ENDDO
              ENDDO
           ENDIF

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

c    Transformation of the radiative tendencies:
c    -------------------------------------------

c          Net radiative surface flux (W.m-2)
c          ~~~~~~~~~~~~~~~~~~~~~~~~~~
c
           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)
               IF(callslope) THEN
                 sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
                 fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig)
               ENDIF
           ENDDO

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

      ENDIF ! of IF (callrad)

c-----------------------------------------------------------------------
c    3. Gravity wave and subgrid scale topography drag :
c    -------------------------------------------------


      IF(calllott)THEN

        CALL calldrag_noro(ngrid,nlayer,ptimestep,
     &                 zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
 
        DO l=1,nlayer
          DO ig=1,ngrid
            pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
            pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
            pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
          ENDDO
        ENDDO
      ENDIF

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

      IF (calldifv) THEN

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

         zdum1(:,:)=0
         zdum2(:,:)=0
         do l=1,nlayer
            do ig=1,ngrid
               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
            enddo
         enddo

c ----------------------
c Treatment of a special case : using new surface layer (Richardson based)
c without using the thermals in gcm and mesoscale can yield problems in
c weakly unstable situations when winds are near to 0. For those cases, we add
c a unit subgrid gustiness. Remember that thermals should be used we using the 
c Richardson based surface layer model.
        IF ( .not.calltherm 
     .       .and. callrichsl 
     .       .and. .not.turb_resolved) THEN
          DO ig=1, ngrid
             IF (zh(ig,1) .lt. tsurf(ig)) THEN
               wstar(ig)=1.
               hfmax_th(ig)=0.2
             ELSE
               wstar(ig)=0.
               hfmax_th(ig)=0.
             ENDIF      
          ENDDO
        ENDIF
c ----------------------

         IF (tke_heat_flux .ne. 0.) THEN
             zz1(:)=(pt(:,1)+pdt(:,1)*ptimestep)*(r/g)*
     &                 (-alog(zplay(:,1)/zplev(:,1)))
             pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1)
         ENDIF

c        Calling vdif (Martian version WITH CO2 condensation)
         CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
     $        ptimestep,capcal,lwrite,
     $        zplay,zplev,zzlay,zzlev,z0,
     $        pu,pv,zh,pq,tsurf,emis,qsurf,
     $        zdum1,zdum2,zdh,pdq,zflubid,
     $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
     &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,sensibFlux)


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

         IF (.not.turb_resolved) THEN
          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

          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)
         ELSE
           write (*,*) '******************************************'
           write (*,*) '** LES mode: the difv part is only used to'
           write (*,*) '**  - provide HFX and UST to the dynamics'
           write (*,*) '**  - update TSURF'
           write (*,*) '******************************************'
           !! Specific treatment for lifting in turbulent-resolving mode (AC)
           IF (lifting .and. doubleq) THEN
             !! lifted dust is injected in the first layer. 
             !! Sedimentation must be called after turbulent mixing, i.e. on next step, after WRF. 
             !! => lifted dust is not incremented before the sedimentation step.
             zdqdif(1:ngrid,1,1:nq)=0.
             zdqdif(1:ngrid,1,igcm_dust_number) = 
     .                  -zdqsdif(1:ngrid,igcm_dust_number)
             zdqdif(1:ngrid,1,igcm_dust_mass) = 
     .                  -zdqsdif(1:ngrid,igcm_dust_mass)
             zdqdif(1:ngrid,2:nlayer,1:nq) = 0.
             DO iq=1, nq
              IF ((iq .ne. igcm_dust_mass)
     &        .and. (iq .ne. igcm_dust_number)) THEN
                zdqsdif(:,iq)=0.
              ENDIF
             ENDDO
           ELSE
             zdqdif(1:ngrid,1:nlayer,1:nq) = 0.
             zdqsdif(1:ngrid,1:nq) = 0.
           ENDIF
         ENDIF
      ELSE    
         DO ig=1,ngrid
            zdtsurf(ig)=zdtsurf(ig)+
     s        (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
         ENDDO
         IF (turb_resolved) THEN
            write(*,*) 'Turbulent-resolving mode !' 
            write(*,*) 'Please set calldifv to T in callphys.def'
            STOP
         ENDIF
      ENDIF ! of IF (calldifv)

c-----------------------------------------------------------------------
c   5. Thermals :
c   -----------------------------

      if(calltherm .and. .not.turb_resolved) then
 
        call calltherm_interface(ngrid,nlayer,nq,
     $ tracer,igcm_co2,
     $ zzlev,zzlay,
     $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
     $ zplay,zplev,pphi,zpopsk,
     $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
     $ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
     
         DO l=1,nlayer
           DO ig=1,ngrid
              pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
              pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
              pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
              q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
           ENDDO
        ENDDO
 
        DO ig=1,ngrid
          q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
        ENDDO      
    
        if (tracer) then
        DO iq=1,nq
         DO l=1,nlayer
           DO ig=1,ngrid
             pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
           ENDDO
         ENDDO
        ENDDO
        endif

        lmax_th_out(:)=real(lmax_th(:))

      else   !of if calltherm
        lmax_th(:)=0
        wstar(:)=0.
        hfmax_th(:)=0.
        lmax_th_out(:)=0.
      end if

c-----------------------------------------------------------------------
c   5. Dry convective adjustment:
c   -----------------------------

      IF(calladj) THEN

         DO l=1,nlayer
            DO ig=1,ngrid
               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
            ENDDO
         ENDDO
         zduadj(:,:)=0
         zdvadj(:,:)=0
         zdhadj(:,:)=0

         CALL convadj(ngrid,nlayer,nq,ptimestep,
     $                zplay,zplev,zpopsk,lmax_th,
     $                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)



c-----------------------------------------------------------------------
c   6. Specific parameterizations for tracers 
c:   -----------------------------------------

      if (tracer) then 

c   6a. Water and ice
c     ---------------

c        ---------------------------------------
c        Water ice condensation in the atmosphere
c        ----------------------------------------
         IF (water) THEN

           call watercloud(ngrid,nlayer,ptimestep,
     &                zplev,zplay,pdpsrf,zzlay, pt,pdt,
     &                pq,pdq,zdqcloud,zdtcloud,
     &                nq,tau,tauscaling,rdust,rice,nuice,
     &                rsedcloud,rhocloud)
     
c Temperature variation due to latent heat release
           if (activice) then
               pdt(1:ngrid,1:nlayer) = 
     &          pdt(1:ngrid,1:nlayer) + 
     &          zdtcloud(1:ngrid,1:nlayer)
           endif
           
! increment water vapour and ice atmospheric tracers tendencies
           pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = 
     &       pdq(1:ngrid,1:nlayer,igcm_h2o_vap) + 
     &       zdqcloud(1:ngrid,1:nlayer,igcm_h2o_vap)
           pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = 
     &       pdq(1:ngrid,1:nlayer,igcm_h2o_ice) + 
     &       zdqcloud(1:ngrid,1:nlayer,igcm_h2o_ice)

! increment dust and ccn masses and numbers
! We need to check that we have Nccn & Ndust > 0
! This is due to single precision rounding problems
           if (microphys) then
             pdq(1:ngrid,1:nlayer,igcm_ccn_mass) = 
     &         pdq(1:ngrid,1:nlayer,igcm_ccn_mass) + 
     &         zdqcloud(1:ngrid,1:nlayer,igcm_ccn_mass)
             pdq(1:ngrid,1:nlayer,igcm_ccn_number) = 
     &         pdq(1:ngrid,1:nlayer,igcm_ccn_number) + 
     &         zdqcloud(1:ngrid,1:nlayer,igcm_ccn_number)
             where (pq(:,:,igcm_ccn_mass) + 
     &       ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
               pdq(:,:,igcm_ccn_mass) = 
     &         - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
               pdq(:,:,igcm_ccn_number) = 
     &         - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
             end where
             where (pq(:,:,igcm_ccn_number) + 
     &       ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
               pdq(:,:,igcm_ccn_mass) = 
     &         - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
               pdq(:,:,igcm_ccn_number) = 
     &         - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
             end where
           endif

           if (scavenging) then
             pdq(1:ngrid,1:nlayer,igcm_dust_mass) = 
     &         pdq(1:ngrid,1:nlayer,igcm_dust_mass) + 
     &         zdqcloud(1:ngrid,1:nlayer,igcm_dust_mass)
             pdq(1:ngrid,1:nlayer,igcm_dust_number) = 
     &         pdq(1:ngrid,1:nlayer,igcm_dust_number) + 
     &         zdqcloud(1:ngrid,1:nlayer,igcm_dust_number)
             where (pq(:,:,igcm_dust_mass) + 
     &       ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
               pdq(:,:,igcm_dust_mass) = 
     &         - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
               pdq(:,:,igcm_dust_number) = 
     &         - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
             end where
             where (pq(:,:,igcm_dust_number) + 
     &       ptimestep*pdq(:,:,igcm_dust_number) < 0.)
               pdq(:,:,igcm_dust_mass) = 
     &         - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
               pdq(:,:,igcm_dust_number) = 
     &         - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
             end where
           endif ! of if scavenging
                      

         END IF  ! of IF (water)

c   6b. Aerosol particles
c     -------------------

c        ----------
c        Dust devil :
c        ----------
         IF(callddevil) then 
           call dustdevil(ngrid,nlayer,nq, zplev,pu,pv,pt, tsurf,q2,
     &                zdqdev,zdqsdev)
 
           if (dustbin.ge.1) then
              do iq=1,nq
                 DO l=1,nlayer
                    DO ig=1,ngrid
                       pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
                    ENDDO
                 ENDDO
              enddo
              do iq=1,nq
                 DO ig=1,ngrid
                    dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
                 ENDDO
              enddo
           endif  ! of if (dustbin.ge.1)

         END IF ! of IF (callddevil)

c        ------------- 
c        Sedimentation :   acts also on water ice
c        ------------- 
         IF (sedimentation) THEN 
           zdqsed(1:ngrid,1:nlayer,1:nq)=0
           zdqssed(1:ngrid,1:nq)=0

           call callsedim(ngrid,nlayer, ptimestep,
     &                zplev,zzlev, zzlay, pt, pdt, rdust, rice,
     &                rsedcloud,rhocloud,
     &                pq, pdq, zdqsed, zdqssed,nq, 
     &                tau,tauscaling)

           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=1, nq
             DO ig=1,ngrid
                dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
             ENDDO
           ENDDO
         END IF   ! of IF (sedimentation)
        
c Add lifted dust to tendancies after sedimentation in the LES (AC)
      IF (turb_resolved) 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
      ENDIF

c
c   6c. Chemical species
c     ------------------

#ifndef MESOSCALE
c        --------------
c        photochemistry :
c        --------------
         IF (photochem .or. thermochem) then

!           dust and ice surface area
            call surfacearea(ngrid, nlayer, naerkind,
     $                       ptimestep, zplay, zzlay, 
     $                       pt, pq, pdq, nq, 
     $                       rdust, rice, tau, tauscaling, 
     $                       surfdust, surfice)
!           call photochemistry
            call calchim(ngrid,nlayer,nq,
     &                   ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0,
     $                   zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
     $                   zdqcloud,zdqscloud,tauref,co2ice,
     $                   pu,pdu,pv,pdv,surfdust,surfice)

           ! increment values of tracers:
           DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
                      ! tracers is zero anyways
             DO l=1,nlayer
               DO ig=1,ngrid
                 pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
               ENDDO
             ENDDO
           ENDDO ! of DO iq=1,nq
           
           ! add condensation tendency for H2O2
           if (igcm_h2o2.ne.0) then
             DO l=1,nlayer
               DO ig=1,ngrid
                 pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
     &                                +zdqcloud(ig,l,igcm_h2o2)
               ENDDO
             ENDDO
           endif

           ! increment surface values of tracers:
           DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
                      ! tracers is zero anyways
             DO ig=1,ngrid
               dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
             ENDDO
           ENDDO ! of DO iq=1,nq

           ! add condensation tendency for H2O2
           if (igcm_h2o2.ne.0) then
             DO ig=1,ngrid
               dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
     &                                +zdqscloud(ig,igcm_h2o2)
             ENDDO
           endif

         END IF  ! of IF (photochem.or.thermochem)
#endif

c   6d. Updates
c     ---------

        DO iq=1, nq
          DO ig=1,ngrid

c       ---------------------------------
c       Updating tracer budget on surface
c       ---------------------------------
            qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)

          ENDDO  ! (ig)
        ENDDO    ! (iq)

      endif !  of if (tracer) 

#ifndef MESOSCALE
c-----------------------------------------------------------------------
c   7. THERMOSPHERE CALCULATION
c-----------------------------------------------------------------------

      if (callthermos) then
        call thermosphere(ngrid,nlayer,nq,zplev,zplay,dist_sol,
     $     mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
     &     pt,pq,pu,pv,pdt,pdq,
     $     zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)

        DO l=1,nlayer
          DO ig=1,ngrid
            dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
            pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
     &                         +zdteuv(ig,l)
            pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
            pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
            DO iq=1, nq
              pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
            ENDDO
          ENDDO
        ENDDO

      endif ! of if (callthermos)
#endif

c-----------------------------------------------------------------------
c   8. Carbon dioxide condensation-sublimation:
c     (should be the last atmospherical physical process to be computed)
c   -------------------------------------------

      IF (tituscap) THEN
         !!! get the actual co2 seasonal cap from Titus observations
         CALL geticecover( ngrid, 180.*zls/pi,
     .                  180.*long/pi, 180.*lati/pi, co2ice )
         co2ice = co2ice * 10000.
      ENDIF
      
      
      pdpsrf(:) = 0

      IF (callcond) THEN
         CALL newcondens(ngrid,nlayer,nq,ptimestep,
     $              capcal,zplay,zplev,tsurf,pt,
     $              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
     $              co2ice,albedo,emis,
     $              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
     $              fluxsurf_sw,zls) 

         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

         IF (tracer) THEN
           DO iq=1, nq
            DO l=1,nlayer
              DO ig=1,ngrid
                pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq) 
              ENDDO
            ENDDO
           ENDDO
         ENDIF ! of IF (tracer)

#ifndef MESOSCALE
        ! update surface pressure
        DO ig=1,ngrid
          ps(ig) = zplev(ig,1) + pdpsrf(ig)*ptimestep
        ENDDO
        
        ! update pressure levels
        DO l=1,nlayer
         DO ig=1,ngrid
          zplay(ig,l) = aps(l) + bps(l)*ps(ig)
          zplev(ig,l) = ap(l)  + bp(l)*ps(ig)
         ENDDO
        ENDDO
        zplev(:,nlayer+1) = 0.
        
        ! update layers altitude
        DO l=2,nlayer
          DO ig=1,ngrid
           z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
           z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
           zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
          ENDDO
        ENDDO
#endif
      
      ENDIF  ! of IF (callcond)
      

c-----------------------------------------------------------------------
c   9. Surface  and sub-surface soil temperature
c-----------------------------------------------------------------------
c
c
c   9.1 Increment Surface temperature:
c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      DO ig=1,ngrid
         tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) 
      ENDDO

c  Prescribe a cold trap at south pole (except at high obliquity !!)
c  Temperature at the surface is set there to be the temperature
c  corresponding to equilibrium temperature between phases of CO2


      IF (tracer.AND.water.AND.(ngrid.NE.1)) THEN
!#ifndef MESOSCALE 
!         if (caps.and.(obliquit.lt.27.)) then => now done in newcondens
           ! NB: Updated surface pressure, at grid point 'ngrid', is
           !     ps(ngrid)=zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
!           tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095*
!     &                     (zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
!           tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095*ps(ngrid)))
!         endif
!#endif
c       -------------------------------------------------------------
c       Change of surface albedo in case of ground frost
c       everywhere except on the north permanent cap and in regions
c       covered by dry ice. 
c              ALWAYS PLACE these lines after newcondens !!!
c       -------------------------------------------------------------
         do ig=1,ngrid
           if ((co2ice(ig).eq.0).and.
     &        (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then
              albedo(ig,1) = albedo_h2o_ice
              albedo(ig,2) = albedo_h2o_ice
c              write(*,*) "frost thickness", qsurf(ig,igcm_h2o_ice)
c              write(*,*) "physiq.F frost :"
c     &        ,lati(ig)*180./pi, long(ig)*180./pi
           endif
         enddo  ! of do ig=1,ngrid
      ENDIF  ! of IF (tracer.AND.water.AND.(ngrid.NE.1))

c
c   9.2 Compute soil temperatures and subsurface heat flux:
c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      IF (callsoil) THEN
c       Thermal inertia feedback
        IF (tifeedback) THEN
         CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil)
         CALL soil(ngrid,nsoilmx,.false.,inertiesoil,
     s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
        ELSE
         CALL soil(ngrid,nsoilmx,.false.,inertiedat,
     s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
        ENDIF
      ENDIF
      

c-----------------------------------------------------------------------
c  10. Write output files
c  ----------------------

c    -------------------------------
c    Dynamical fields incrementation
c    -------------------------------
c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
      ! temperature, zonal and meridional wind
      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
        ENDDO
      ENDDO

      ! 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

      ! Density 
      DO l=1,nlayer
         DO ig=1,ngrid
            rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
         ENDDO
      ENDDO

      ! Potential Temperature

       DO ig=1,ngrid
          DO l=1,nlayer
              zh(ig,l) = zt(ig,l)*(zplev(ig,1)/zplay(ig,l))**rcp
          ENDDO
       ENDDO


c    Compute surface stress : (NB: z0 is a common in surfdat.h)
c     DO ig=1,ngrid
c        cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
c        zstress(ig) = rho(ig,1)*cd*(zu(ig,1)**2 + zv(ig,1)**2)
c     ENDDO

c     Sum of fluxes in solar spectral bands (for output only)
      DO ig=1,ngrid
             fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
             fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
      ENDDO
c ******* TEST ******************************************************
      ztim1 = 999
      DO l=1,nlayer
        DO ig=1,ngrid
           if (pt(ig,l).lt.ztim1) then
               ztim1 = pt(ig,l)
               igmin = ig
               lmin = l 
           end if
        ENDDO
      ENDDO
      if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
        write(*,*) 'PHYSIQ: stability WARNING :'
        write(*,*) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
     &              'ig l =', igmin, lmin
      end if
c *******************************************************************

c     ---------------------
c     Outputs to the screen 
c     ---------------------

      IF (lwrite) THEN
         PRINT*,'Global diagnostics for the physics'
         PRINT*,'Variables and their increments x and dx/dt * dt'
         WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
         WRITE(*,'(2f10.5,2f15.5)')
     s   tsurf(igout),zdtsurf(igout)*ptimestep,
     s   zplev(igout,1),pdpsrf(igout)*ptimestep
         WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
         WRITE(*,'(i4,6f10.5)') (l,
     s   pu(igout,l),pdu(igout,l)*ptimestep,
     s   pv(igout,l),pdv(igout,l)*ptimestep,
     s   pt(igout,l),pdt(igout,l)*ptimestep,
     s   l=1,nlayer)
      ENDIF ! of IF (lwrite)

c        ----------------------------------------------------------
c        ----------------------------------------------------------
c        INTERPOLATIONS IN THE SURFACE-LAYER
c        ----------------------------------------------------------
c        ----------------------------------------------------------

           n_out=0 ! number of elements in the z_out array.
                   ! for z_out=[3.,2.,1.,0.5,0.1], n_out must be set
                   ! to 5
           IF (n_out .ne. 0) THEN

           IF(.NOT. ALLOCATED(z_out)) ALLOCATE(z_out(n_out))
           IF(.NOT. ALLOCATED(T_out)) ALLOCATE(T_out(ngrid,n_out))
           IF(.NOT. ALLOCATED(u_out)) ALLOCATE(u_out(ngrid,n_out))

           z_out(:)=[3.,2.,1.,0.5,0.1]
           u_out(:,:)=0.
           T_out(:,:)=0.

           call pbl_parameters(ngrid,nlayer,ps,zplay,z0,
     & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,n_out,
     & T_out,u_out,ustar,tstar,L_mo,vhf,vvv)
                   ! pourquoi ustar recalcule ici? fait dans vdifc.

#ifndef MESOSCALE
            IF (ngrid .eq. 1) THEN
               dimout=0
            ELSE
               dimout=2
            ENDIF
            DO n=1,n_out
               write(zstring, '(F8.6)') z_out(n)
               call WRITEDIAGFI(ngrid,'T_out_'//trim(zstring),
     &   'potential temperature at z_out','K',dimout,T_out(:,n))
               call WRITEDIAGFI(ngrid,'u_out_'//trim(zstring),
     &   'horizontal velocity norm at z_out','m/s',dimout,u_out(:,n))
            ENDDO
            call WRITEDIAGFI(ngrid,'u_star',
     &   'friction velocity','m/s',dimout,ustar)
            call WRITEDIAGFI(ngrid,'teta_star',
     &   'friction potential temperature','K',dimout,tstar)
            call WRITEDIAGFI(ngrid,'L',
     &   'Monin Obukhov length','m',dimout,L_mo)
            call WRITEDIAGFI(ngrid,'vvv',
     &   'Vertical velocity variance at zout','m',dimout,vvv)
            call WRITEDIAGFI(ngrid,'vhf',
     &   'Vertical heat flux at zout','m',dimout,vhf)
#else
         T_out1(:)=T_out(:,1)
         u_out1(:)=u_out(:,1)
#endif

         ENDIF

c        ----------------------------------------------------------
c        ----------------------------------------------------------
c        END OF SURFACE LAYER INTERPOLATIONS
c        ----------------------------------------------------------
c        ----------------------------------------------------------

      IF (ngrid.NE.1) THEN

#ifndef MESOSCALE
c        -------------------------------------------------------------------
c        Writing NetCDF file  "RESTARTFI" at the end of the run
c        -------------------------------------------------------------------
c        Note: 'restartfi' is stored just before dynamics are stored
c              in 'restart'. Between now and the writting of 'restart',
c              there will have been the itau=itau+1 instruction and
c              a reset of 'time' (lastacll = .true. when itau+1= itaufin)
c              thus we store for time=time+dtvr

         IF(   ((ecritstart.GT.0) .and. 
     .          (MOD(icount*iphysiq,ecritstart).EQ.0)) 
     .    .or. lastcall  ) THEN
         
          ztime_fin = pday + ptime  + ptimestep/(float(iphysiq)*daysec)
     .               - day_ini - time_phys
          print*, pday,ptime,day_ini, time_phys
          write(*,'(A,I7,A,F12.5)') 
     .         'PHYSIQ: Ecriture du fichier restartfi ; icount=',
     .          icount,' date=',ztime_fin
            
            
          call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,
     .                ptimestep,ztime_fin,
     .                tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling)
          
         ENDIF
#endif

c        -------------------------------------------------------------------
c        Calculation of diagnostic variables written in both stats and
c          diagfi files
c        -------------------------------------------------------------------

         if (tracer) then
         
           if(doubleq) then
              do ig=1,ngrid 
                dqdustsurf(ig) = 
     &                zdqssed(ig,igcm_dust_mass)*tauscaling(ig)
                dndustsurf(ig) = 
     &                zdqssed(ig,igcm_dust_number)*tauscaling(ig)
                ndust(ig,:) =
     &                pq(ig,:,igcm_dust_number)*tauscaling(ig)
                qdust(ig,:) =
     &                pq(ig,:,igcm_dust_mass)*tauscaling(ig)
              enddo
              if (scavenging) then
                do ig=1,ngrid 
                  dqdustsurf(ig) = dqdustsurf(ig) + 
     &                     zdqssed(ig,igcm_ccn_mass)*tauscaling(ig)
                  dndustsurf(ig) = dndustsurf(ig) + 
     &                     zdqssed(ig,igcm_ccn_number)*tauscaling(ig)
                  nccn(ig,:) =
     &                     pq(ig,:,igcm_ccn_number)*tauscaling(ig)
                  qccn(ig,:) =
     &                     pq(ig,:,igcm_ccn_mass)*tauscaling(ig)
                enddo
              endif
           endif
                    
           if (water) then
             mtot(:)=0
             icetot(:)=0
             rave(:)=0
             tauTES(:)=0
             do ig=1,ngrid 
               do l=1,nlayer
                 mtot(ig) = mtot(ig) + 
     &                      zq(ig,l,igcm_h2o_vap) * 
     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
                 icetot(ig) = icetot(ig) + 
     &                        zq(ig,l,igcm_h2o_ice) * 
     &                        (zplev(ig,l) - zplev(ig,l+1)) / g
c                Computing abs optical depth at 825 cm-1 in each
c                  layer to simulate NEW TES retrieval
                 Qabsice = min(
     &             max(0.4e6*rice(ig,l)*(1.+nuice_ref)-0.05 ,0.),1.2
     &                        )
                 opTES(ig,l)= 0.75 * Qabsice * 
     &             zq(ig,l,igcm_h2o_ice) *
     &             (zplev(ig,l) - zplev(ig,l+1)) / g
     &             / (rho_ice * rice(ig,l) * (1.+nuice_ref))
                 tauTES(ig)=tauTES(ig)+ opTES(ig,l) 
               enddo
c              rave(ig)=rave(ig)/max(icetot(ig),1.e-30)       ! mass weight
c               if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
             enddo
             call watersat(ngrid*nlayer,zt,zplay,zqsat)
             satu(:,:) = zq(:,:,igcm_h2o_vap)/zqsat(:,:)

             if (scavenging) then
               Nccntot(:)= 0
               Mccntot(:)= 0
               rave(:)=0
               do ig=1,ngrid 
                 do l=1,nlayer
                    Nccntot(ig) = Nccntot(ig) + 
     &              zq(ig,l,igcm_ccn_number)*tauscaling(ig)
     &              *(zplev(ig,l) - zplev(ig,l+1)) / g
                    Mccntot(ig) = Mccntot(ig) + 
     &              zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
     &              *(zplev(ig,l) - zplev(ig,l+1)) / g
cccc Column integrated effective ice radius 
cccc is weighted by total ice surface area (BETTER than total ice mass) 
                    rave(ig) = rave(ig) + 
     &                      tauscaling(ig) *
     &                      zq(ig,l,igcm_ccn_number) *
     &                      (zplev(ig,l) - zplev(ig,l+1)) / g * 
     &                      rice(ig,l) * rice(ig,l)*  (1.+nuice_ref)
                 enddo
               rave(ig)=(icetot(ig)/rho_ice+Mccntot(ig)/rho_dust)*0.75
     &                  /max(pi*rave(ig),1.e-30) ! surface weight
               if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
               enddo
             else ! of if (scavenging)
               rave(:)=0
               do ig=1,ngrid 
                 do l=1,nlayer
                 rave(ig) = rave(ig) + 
     &                      zq(ig,l,igcm_h2o_ice) *
     &                      (zplev(ig,l) - zplev(ig,l+1)) / g * 
     &                      rice(ig,l) * (1.+nuice_ref)
                 enddo
                 rave(ig) = max(rave(ig) / 
     &             max(icetot(ig),1.e-30),1.e-30) ! mass weight
               enddo
             endif ! of if (scavenging)

           endif ! of if (water)
         endif ! of if (tracer)

#ifndef MESOSCALE
c        -----------------------------------------------------------------
c        WSTATS: Saving statistics
c        -----------------------------------------------------------------
c        ("stats" stores and accumulates 8 key variables in file "stats.nc"
c        which can later be used to make the statistic files of the run:
c        "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,"co2ice","CO2 ice cover",
     &                "kg.m-2",2,co2ice)
        call wstats(ngrid,"tauref","reference dod at 610 Pa","NU",
     &                2,tauref)
        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)
        call wstats(ngrid,"w","Vertical (down-up) wind",
     &                "m.s-1",3,pw)
        call wstats(ngrid,"rho","Atmospheric density","kg/m3",3,rho)
        call wstats(ngrid,"pressure","Pressure","Pa",3,zplay)
          call wstats(ngrid,"q2",
     &                "Boundary layer eddy kinetic energy",
     &                "m2.s-2",3,q2)
          call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
     &                emis)
c          call wstats(ngrid,"ssurf","Surface stress","N.m-2",
c    &                2,zstress)
c          call wstats(ngrid,"sw_htrt","sw heat.rate",
c    &                 "W.m-2",3,zdtsw)
c          call wstats(ngrid,"lw_htrt","lw heat.rate",
c    &                 "W.m-2",3,zdtlw)

          if (calltherm) then
            call wstats(ngrid,"zmax_th","Height of thermals",
     &                "m",2,zmax_th)
            call wstats(ngrid,"hfmax_th","Max thermals heat flux",
     &                "K.m/s",2,hfmax_th)
            call wstats(ngrid,"wstar",
     &                "Max vertical velocity in thermals",
     &                "m/s",2,wstar)
          endif

           if (tracer) then
             if (water) then
               vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
               call wstats(ngrid,"vmr_h2ovap",
     &                    "H2O vapor volume mixing ratio","mol/mol",
     &                    3,vmr)
               vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
               call wstats(ngrid,"vmr_h2oice",
     &                    "H2O ice volume mixing ratio","mol/mol",
     &                    3,vmr)
               vmr=zqsat(1:ngrid,1:nlayer)
     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
               call wstats(ngrid,"vmr_h2osat",
     &                    "saturation volume mixing ratio","mol/mol",
     &                    3,vmr)
               call wstats(ngrid,"h2o_ice_s",
     &                    "surface h2o_ice","kg/m2",
     &                    2,qsurf(1,igcm_h2o_ice))
               call wstats(ngrid,'albedo',
     &                         'albedo',
     &                         '',2,albedo(1,1))
               call wstats(ngrid,"mtot",
     &                    "total mass of water vapor","kg/m2",
     &                    2,mtot)
               call wstats(ngrid,"icetot",
     &                    "total mass of water ice","kg/m2",
     &                    2,icetot)
               call wstats(ngrid,"reffice",
     &                    "Mean reff","m",
     &                    2,rave)
              call wstats(ngrid,"Nccntot",
     &                    "condensation nuclei","Nbr/m2",
     &                    2,Nccntot)
              call wstats(ngrid,"Mccntot",
     &                    "condensation nuclei mass","kg/m2",
     &                    2,Mccntot)
              call wstats(ngrid,"rice",
     &                    "Ice particle size","m",
     &                    3,rice)
               if (.not.activice) then
                 call wstats(ngrid,"tauTESap",
     &                      "tau abs 825 cm-1","",
     &                      2,tauTES)
               else
                 call wstats(ngrid,'tauTES',
     &                   'tau abs 825 cm-1',
     &                  '',2,taucloudtes)
               endif

             endif ! of if (water)
             
             
           if (dustbin.ne.0) then
          
             call wstats(ngrid,'tau','taudust','SI',2,tau(1,1))
             
             if (doubleq) then
c            call wstats(ngrid,'qsurf','qsurf',
c     &                       'kg.m-2',2,qsurf(1,igcm_dust_mass))
c            call wstats(ngrid,'Nsurf','N particles',
c     &                       'N.m-2',2,qsurf(1,igcm_dust_number))
c            call wstats(ngrid,'dqsdev','ddevil lift',
c    &                       'kg.m-2.s-1',2,zdqsdev(1,1))
c            call wstats(ngrid,'dqssed','sedimentation',
c     &                       'kg.m-2.s-1',2,zdqssed(1,1))
c            call wstats(ngrid,'dqsdif','diffusion',
c     &                       'kg.m-2.s-1',2,zdqsdif(1,1))
               call wstats(ngrid,'dqsdust',
     &                        'deposited surface dust mass',
     &                        'kg.m-2.s-1',2,dqdustsurf)
               call wstats(ngrid,'dqndust',
     &                        'deposited surface dust number',
     &                        'number.m-2.s-1',2,dndustsurf)
               call wstats(ngrid,'reffdust','reffdust',
     &                        'm',3,rdust*ref_r0)
               call wstats(ngrid,'dustq','Dust mass mr',
     &                        'kg/kg',3,qdust)
               call wstats(ngrid,'dustN','Dust number',
     &                        'part/kg',3,ndust)
             else
               do iq=1,dustbin
                 write(str2(1:2),'(i2.2)') iq
                 call wstats(ngrid,'q'//str2,'mix. ratio',
     &                         'kg/kg',3,zq(1,1,iq))
                 call wstats(ngrid,'qsurf'//str2,'qsurf',
     &                         'kg.m-2',2,qsurf(1,iq))
               end do
             endif ! (doubleq)

             if (scavenging) then
               call wstats(ngrid,'ccnq','CCN mass mr',
     &                        'kg/kg',3,qccn)
               call wstats(ngrid,'ccnN','CCN number',
     &                        'part/kg',3,nccn)
             endif ! (scavenging)
          
          endif ! (dustbin.ne.0)


             if (thermochem.or.photochem) then
                do iq=1,nq
                  if (noms(iq) .ne. "dust_mass" .and.
     $                 noms(iq) .ne. "dust_number" .and.
     $                 noms(iq) .ne. "ccn_mass" .and.
     $                 noms(iq) .ne. "ccn_number" .and.
     $                 noms(iq) .ne. "h2o_vap" .and.
     $                 noms(iq) .ne. "h2o_ice") then
                    vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
     &                          *mmean(1:ngrid,1:nlayer)/mmol(iq)
                    rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
     &                          *rho(1:ngrid,1:nlayer)*n_avog/
     &                           (1000*mmol(iq))
                   call wstats(ngrid,"vmr_"//trim(noms(iq)),
     $                         "Volume mixing ratio","mol/mol",3,vmr)
!                   call wstats(ngrid,"rho_"//trim(noms(iq)),
!     $                     "Number density","cm-3",3,rhopart)
!                   call writediagfi(ngrid,"rho_"//trim(noms(iq)),
!     $                     "Number density","cm-3",3,rhopart)
                   do ig = 1,ngrid
                      colden(ig,iq) = 0.                           
                   end do
                   do l=1,nlayer                                    
                      do ig=1,ngrid                                  
                         colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq)
     $                                  *(zplev(ig,l)-zplev(ig,l+1)) 
     $                                  *6.022e22/(mmol(iq)*g)       
                      end do                                        
                   end do                                          
                   call wstats(ngrid,"c_"//trim(noms(iq)),           
     $                         "column","mol cm-2",2,colden(1,iq))  
                  end if ! of if (noms(iq) .ne. "dust_mass" ...)
                end do ! of do iq=1,nq
             end if ! of if (thermochem.or.photochem)

           end if ! of if (tracer)

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

         ENDIF !if callstats

c        (Store EOF for Mars Climate database software)
         IF (calleofdump) THEN
            CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
         ENDIF
#endif
!endif of ifndef MESOSCALE

#ifdef MESOSCALE
      
      !! see comm_wrf. 
      !! not needed when an array is already in a shared module.
      !! --> example : hfmax_th, zmax_th

      !state  real  HR_SW     ikj   misc  1  -  h  "HR_SW"     "HEATING RATE SW"                 "K/s"
      comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
      !state  real  HR_LW     ikj   misc  1  -  h  "HR_LW"     "HEATING RATE LW"                 "K/s"
      comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
      !state  real  SWDOWNZ    ij   misc  1  -  h  "SWDOWNZ"   "DOWNWARD SW FLUX AT SURFACE"     "W m-2"
      comm_SWDOWNZ(1:ngrid) = fluxsurf_sw_tot(1:ngrid)
      !state  real  TAU_DUST   ij   misc  1  -  h  "TAU_DUST"  "REFERENCE VISIBLE DUST OPACITY"  ""
      comm_TAU_DUST(1:ngrid) = tauref(1:ngrid)
      !state  real  RDUST     ikj   misc  1  -  h  "RDUST"     "DUST RADIUS"                     "m"
      comm_RDUST(1:ngrid,1:nlayer) = rdust(1:ngrid,1:nlayer)
      !state  real  QSURFDUST  ij   misc  1  -  h  "QSURFDUST" "DUST MASS AT SURFACE"            "kg m-2"
      IF (igcm_dust_mass .ne. 0) THEN
        comm_QSURFDUST(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass)
      ELSE
        comm_QSURFDUST(1:ngrid) = 0.
      ENDIF
      !state  real  MTOT       ij   misc  1  -  h  "MTOT"      "TOTAL MASS WATER VAPOR in pmic"  "pmic"
      comm_MTOT(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
      !state  real  ICETOT     ij   misc  1  -  h  "ICETOT"    "TOTAL MASS WATER ICE"            "kg m-2"
      comm_ICETOT(1:ngrid) = icetot(1:ngrid) * 1.e6 / rho_ice
      !state  real  VMR_ICE   ikj   misc  1  -  h  "VMR_ICE"   "VOL. MIXING RATIO ICE"           "ppm"
      IF (igcm_h2o_ice .ne. 0) THEN
        comm_VMR_ICE(1:ngrid,1:nlayer) = 1.e6
     .    * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
     .    * mmean(1:ngrid,1:nlayer) / mmol(igcm_h2o_ice)
      ELSE
        comm_VMR_ICE(1:ngrid,1:nlayer) = 0.
      ENDIF
      !state  real  TAU_ICE    ij   misc  1  -  h  "TAU_ICE"   "CLOUD OD at 825 cm-1 TES"        ""
      comm_TAU_ICE(1:ngrid) = taucloudtes(1:ngrid)
      !state  real  RICE      ikj   misc  1  -  h  "RICE"      "ICE RADIUS"                      "m"
      comm_RICE(1:ngrid,1:nlayer) = rice(1:ngrid,1:nlayer)

   
      !! calculate sensible heat flux in W/m2 for outputs
      !! -- the one computed in vdifc is not the real one
      !! -- vdifc must have been called
      if (.not.callrichsl) then
        sensibFlux(1:ngrid) = zflubid(1:ngrid)
     .         - capcal(1:ngrid)*zdtsdif(1:ngrid)
      else
        sensibFlux(1:ngrid) = 
     &   (pplay(1:ngrid,1)/(r*pt(1:ngrid,1)))*cpp
     &   *sqrt(pu(1:ngrid,1)*pu(1:ngrid,1)+pv(1:ngrid,1)*pv(1:ngrid,1)
     &         +(log(1.+0.7*wstar(1:ngrid) + 2.3*wstar(1:ngrid)**2))**2)
     &   *zcdh(1:ngrid)*(tsurf(1:ngrid)-zh(1:ngrid,1))
      endif
          
#else
#ifndef MESOINI

c        ==========================================================
c        WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
c          any variable for diagnostic (output with period
c          "ecritphy", set in "run.def")
c        ==========================================================
c        WRITEDIAGFI can ALSO be called from any other subroutines
c        for any variables !!
c        call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
c    &                  emis)
c        call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
c        call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
         call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
     &                  tsurf)
         call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
         call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
     &                                         ,"kg.m-2",2,co2ice)

         call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
     &                  "K",2,zt(1,7))
         call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
     &                  fluxsurf_lw)
         call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
     &                  fluxsurf_sw_tot)
         call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
     &                  fluxtop_lw)
         call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
     &                  fluxtop_sw_tot)
         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,"w","Vertical wind","m.s-1",3,pw)
         call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
c        call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
c        call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
         call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,zplay)
c        call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
c    &                  zstress)
c        call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
c    &                   'w.m-2',3,zdtsw)
c        call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
c    &                   'w.m-2',3,zdtlw)
            if (.not.activice) then
               CALL WRITEDIAGFI(ngrid,'tauTESap',
     &                         'tau abs 825 cm-1',
     &                         '',2,tauTES)
             else
               CALL WRITEDIAGFI(ngrid,'tauTES',
     &                         'tau abs 825 cm-1',
     &                         '',2,taucloudtes)
             endif

#else
     !!! this is to ensure correct initialisation of mesoscale model
        call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
     &                  tsurf)
        call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps) 
        call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
     &                  co2ice)
        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,"emis","Surface emissivity","w.m-1",2,
     &                  emis)
        call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
     &                       "K",3,tsoil)
        call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
     &                       "K",3,inertiedat)
#endif


c        ----------------------------------------------------------
c        Outputs of the CO2 cycle
c        ----------------------------------------------------------

         if (tracer.and.(igcm_co2.ne.0)) then
!          call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
!    &                     "kg/kg",2,zq(1,1,igcm_co2))
          call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
     &                     "kg/kg",3,zq(1,1,igcm_co2))
        
         ! Compute co2 column
         co2col(:)=0
         do l=1,nlayer
           do ig=1,ngrid
             co2col(ig)=co2col(ig)+
     &                  zq(ig,l,igcm_co2)*(zplev(ig,l)-zplev(ig,l+1))/g
           enddo
         enddo
         call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
     &                  co2col)
         endif ! of if (tracer.and.(igcm_co2.ne.0))

c        ----------------------------------------------------------
c        Outputs of the water cycle
c        ----------------------------------------------------------
         if (tracer) then
           if (water) then

#ifdef MESOINI
            !!!! waterice = q01, voir readmeteo.F90
            call WRITEDIAGFI(ngrid,'q01',noms(igcm_h2o_ice),
     &                      'kg/kg',3,
     &                       zq(1:ngrid,1:nlayer,igcm_h2o_ice))
            !!!! watervapor = q02, voir readmeteo.F90
            call WRITEDIAGFI(ngrid,'q02',noms(igcm_h2o_vap),
     &                      'kg/kg',3,
     &                       zq(1:ngrid,1:nlayer,igcm_h2o_vap))
            !!!! surface waterice qsurf02 (voir readmeteo)
            call WRITEDIAGFI(ngrid,'qsurf02','surface tracer',
     &                      'kg.m-2',2,
     &                       qsurf(1:ngrid,igcm_h2o_ice))
#endif

            CALL WRITEDIAGFI(ngrid,'mtot',
     &                       'total mass of water vapor',
     &                       'kg/m2',2,mtot)
            CALL WRITEDIAGFI(ngrid,'icetot',
     &                       'total mass of water ice',
     &                       'kg/m2',2,icetot)
            vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
            call WRITEDIAGFI(ngrid,'vmr_h2oice','h2o ice vmr',
     &                       'mol/mol',3,vmr)
            vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
            call WRITEDIAGFI(ngrid,'vmr_h2ovap','h2o vap vmr',
     &                       'mol/mol',3,vmr)
            CALL WRITEDIAGFI(ngrid,'reffice',
     &                       'Mean reff',
     &                       'm',2,rave)
            if (scavenging) then
              CALL WRITEDIAGFI(ngrid,"Nccntot",
     &                    "condensation nuclei","Nbr/m2",
     &                    2,Nccntot)
              CALL WRITEDIAGFI(ngrid,"Mccntot",
     &                    "mass condensation nuclei","kg/m2",
     &                    2,Mccntot)
            endif
            call WRITEDIAGFI(ngrid,'rice','Ice particle size',
     &                       'm',3,rice)
            call WRITEDIAGFI(ngrid,'h2o_ice_s',
     &                       'surface h2o_ice',
     &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice))
            CALL WRITEDIAGFI(ngrid,'albedo',
     &                         'albedo',
     &                         '',2,albedo(1,1))
              if (tifeedback) then
                 call WRITEDIAGSOIL(ngrid,"soiltemp",
     &                              "Soil temperature","K",
     &                              3,tsoil)
                 call WRITEDIAGSOIL(ngrid,'soilti',
     &                       'Soil Thermal Inertia',
     &                       'J.s-1/2.m-2.K-1',3,inertiesoil)
              endif
           endif !(water)


           if (water.and..not.photochem) then
             iq=nq
c            write(str2(1:2),'(i2.2)') iq
c            call WRITEDIAGFI(ngrid,'dqs'//str2,'dqscloud',
c    &                       'kg.m-2',2,zdqscloud(1,iq))
c            call WRITEDIAGFI(ngrid,'dqch'//str2,'var chim',
c    &                       'kg/kg',3,zdqchim(1,1,iq))
c            call WRITEDIAGFI(ngrid,'dqd'//str2,'var dif',
c    &                       'kg/kg',3,zdqdif(1,1,iq))
c            call WRITEDIAGFI(ngrid,'dqa'//str2,'var adj',
c    &                       'kg/kg',3,zdqadj(1,1,iq))
c            call WRITEDIAGFI(ngrid,'dqc'//str2,'var c',
c    &                       'kg/kg',3,zdqc(1,1,iq))
           endif  !(water.and..not.photochem)
         endif

c        ----------------------------------------------------------
c        Outputs of the dust cycle
c        ----------------------------------------------------------

        call WRITEDIAGFI(ngrid,'tauref',
     &                    'Dust ref opt depth','NU',2,tauref)

         if (tracer.and.(dustbin.ne.0)) then

          call WRITEDIAGFI(ngrid,'tau','taudust','SI',2,tau(1,1))

#ifndef MESOINI
           if (doubleq) then
c            call WRITEDIAGFI(ngrid,'qsurf','qsurf',
c     &                       'kg.m-2',2,qsurf(1,igcm_dust_mass))
c            call WRITEDIAGFI(ngrid,'Nsurf','N particles',
c     &                       'N.m-2',2,qsurf(1,igcm_dust_number))
c            call WRITEDIAGFI(ngrid,'dqsdev','ddevil lift',
c    &                       'kg.m-2.s-1',2,zdqsdev(1,1))
c            call WRITEDIAGFI(ngrid,'dqssed','sedimentation',
c     &                       'kg.m-2.s-1',2,zdqssed(1,1))
c            call WRITEDIAGFI(ngrid,'dqsdif','diffusion',
c     &                       'kg.m-2.s-1',2,zdqsdif(1,1))
c             call WRITEDIAGFI(ngrid,'sedice','sedimented ice',
c     &                       'kg.m-2.s-1',2,zdqssed(:,igcm_h2o_ice))
c             call WRITEDIAGFI(ngrid,'subice','sublimated ice',
c     &                       'kg.m-2.s-1',2,zdqsdif(:,igcm_h2o_ice))
             call WRITEDIAGFI(ngrid,'dqsdust',
     &                        'deposited surface dust mass',
     &                        'kg.m-2.s-1',2,dqdustsurf)
             call WRITEDIAGFI(ngrid,'dqndust',
     &                        'deposited surface dust number',
     &                        'number.m-2.s-1',2,dndustsurf)
             call WRITEDIAGFI(ngrid,'reffdust','reffdust',
     &                        'm',3,rdust*ref_r0)
             call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
     &                        'kg/kg',3,qdust)
             call WRITEDIAGFI(ngrid,'dustN','Dust number',
     &                        'part/kg',3,ndust)
c             call WRITEDIAGFI(ngrid,"tauscaling",
c     &                    "dust conversion factor"," ",2,tauscaling)
           else
             do iq=1,dustbin
               write(str2(1:2),'(i2.2)') iq
               call WRITEDIAGFI(ngrid,'q'//str2,'mix. ratio',
     &                         'kg/kg',3,zq(1,1,iq))
               call WRITEDIAGFI(ngrid,'qsurf'//str2,'qsurf',
     &                         'kg.m-2',2,qsurf(1,iq))
             end do
           endif ! (doubleq)

           if (scavenging) then
             call WRITEDIAGFI(ngrid,'ccnq','CCN mass mr',
     &                        'kg/kg',3,qccn)
             call WRITEDIAGFI(ngrid,'ccnN','CCN number',
     &                        'part/kg',3,nccn)
           endif ! (scavenging)

c          if (submicron) then
c            call WRITEDIAGFI(ngrid,'dustsubm','subm mass mr',
c    &                        'kg/kg',3,pq(1,1,igcm_dust_submicron))
c          endif ! (submicron)

#else
!     !!! to initialize mesoscale we need scaled variables
!     !!! because this must correspond to starting point for tracers
!             call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
!     &           'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_mass))
!             call WRITEDIAGFI(ngrid,'dustN','Dust number',
!     &           'part/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_number))
!             call WRITEDIAGFI(ngrid,'ccn','Nuclei mass mr',
!     &           'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_mass))
!             call WRITEDIAGFI(ngrid,'ccnN','Nuclei number',
!     &           'part/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_number))
      if (freedust) then
             call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
     &                        'kg/kg',3,qdust)
             call WRITEDIAGFI(ngrid,'dustN','Dust number',
     &                        'part/kg',3,ndust)
             call WRITEDIAGFI(ngrid,'ccn','CCN mass mr',
     &                        'kg/kg',3,qccn)
             call WRITEDIAGFI(ngrid,'ccnN','CCN number',
     &                        'part/kg',3,nccn)
      else
             call WRITEDIAGFI(ngrid,'dustq','Dust mass mr',
     &                        'kg/kg',3,pq(1,1,igcm_dust_mass))
             call WRITEDIAGFI(ngrid,'dustN','Dust number',
     &                        'part/kg',3,pq(1,1,igcm_dust_number))
             call WRITEDIAGFI(ngrid,'ccn','Nuclei mass mr',
     &                        'kg/kg',3,pq(1,1,igcm_ccn_mass))
             call WRITEDIAGFI(ngrid,'ccnN','Nuclei number',
     &                        'part/kg',3,pq(1,1,igcm_ccn_number))
      endif
#endif

         end if  ! (tracer.and.(dustbin.ne.0))


c        ----------------------------------------------------------
c        Thermospheric outputs
c        ----------------------------------------------------------

         if(callthermos) then 

            call WRITEDIAGFI(ngrid,"q15um","15 um cooling","K/s",
     $           3,zdtnlte)
            call WRITEDIAGFI(ngrid,"quv","UV heating","K/s",
     $           3,zdteuv)
            call WRITEDIAGFI(ngrid,"cond","Thermal conduction","K/s",
     $           3,zdtconduc)
            call WRITEDIAGFI(ngrid,"qnir","NIR heating","K/s",
     $           3,zdtnirco2)


             if (thermochem.or.photochem) then
                do iq=1,nq
                  if (noms(iq) .ne. "dust_mass" .and.
     $                 noms(iq) .ne. "dust_number" .and.
     $                 noms(iq) .ne. "ccn_mass" .and.
     $                 noms(iq) .ne. "ccn_number" .and.
     $                 noms(iq) .ne. "h2o_vap" .and.
     $                 noms(iq) .ne. "h2o_ice") then
                    vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
     &                          *mmean(1:ngrid,1:nlayer)/mmol(iq)
                    rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
     &                          *rho(1:ngrid,1:nlayer)*n_avog/
     &                           (1000*mmol(iq))
                   call writediagfi(ngrid,"rho_"//trim(noms(iq)),
     $                     "Number density","cm-3",3,rhopart)
                   if ((noms(iq).eq."o") .or. (noms(iq).eq."co2").or.
     $                 (noms(iq).eq."o3")) then
                      call writediagfi(ngrid,"vmr_"//trim(noms(iq)),
     $                         "Volume mixing ratio","mol/mol",3,vmr)
                   end if
                   do ig = 1,ngrid
                      colden(ig,iq) = 0.
                   end do
                   do l=1,nlayer
                      do ig=1,ngrid
                         colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq)
     $                                  *(zplev(ig,l)-zplev(ig,l+1))
     $                                  *6.022e22/(mmol(iq)*g)
                      end do
                   end do
                   call writediagfi(ngrid,"c_"//trim(noms(iq)),
     $                             "column","mol cm-2",2,colden(1,iq))
                  end if ! of if (noms(iq) .ne. "dust_mass" ...)
                end do ! of do iq=1,nq
             end if ! of if (thermochem.or.photochem)


         endif  !(callthermos)

c        ----------------------------------------------------------
c        ----------------------------------------------------------
c        PBL OUTPUS
c        ----------------------------------------------------------
c        ----------------------------------------------------------

c        ----------------------------------------------------------
c        Outputs of thermals
c        ----------------------------------------------------------
         if (calltherm) then

!        call WRITEDIAGFI(ngrid,'dtke',
!     &              'tendance tke thermiques','m**2/s**2',
!     &                         3,dtke_th)
!        call WRITEDIAGFI(ngrid,'d_u_ajs',
!     &              'tendance u thermiques','m/s',
!     &                         3,pdu_th*ptimestep)
!        call WRITEDIAGFI(ngrid,'d_v_ajs',
!     &              'tendance v thermiques','m/s',
!     &                         3,pdv_th*ptimestep)
!        if (tracer) then
!        if (nq .eq. 2) then
!        call WRITEDIAGFI(ngrid,'deltaq_th',
!     &              'delta q thermiques','kg/kg',
!     &                         3,ptimestep*pdq_th(:,:,2))
!        endif
!        endif

        call WRITEDIAGFI(ngrid,'zmax_th',
     &              'hauteur du thermique','m',
     &                         2,zmax_th)
        call WRITEDIAGFI(ngrid,'hfmax_th',
     &              'maximum TH heat flux','K.m/s',
     &                         2,hfmax_th)
        call WRITEDIAGFI(ngrid,'wstar',
     &              'maximum TH vertical velocity','m/s',
     &                         2,wstar)

         endif

c        ----------------------------------------------------------
c        ----------------------------------------------------------
c        END OF PBL OUTPUS
c        ----------------------------------------------------------
c        ----------------------------------------------------------


c        ----------------------------------------------------------
c        Output in netcdf file "diagsoil.nc" for subterranean
c          variables (output every "ecritphy", as for writediagfi)
c        ----------------------------------------------------------

         ! Write soil temperature
!        call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
!    &                     3,tsoil)
         ! Write surface temperature
!        call writediagsoil(ngrid,"tsurf","Surface temperature","K",
!    &                     2,tsurf)

c        ==========================================================
c        END OF WRITEDIAGFI
c        ==========================================================
#endif
! of ifdef MESOSCALE

      ELSE     ! if(ngrid.eq.1)

#ifndef MESOSCALE
         write(*,'("Ls =",f11.6," tauref(",f4.0," Pa) =",f9.6)')
     &    zls*180./pi,odpref,tauref
c      ----------------------------------------------------------------------
c      Output in grads file "g1d" (ONLY when using testphys1d)
c      (output at every X physical timestep)
c      ----------------------------------------------------------------------
c
c        CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
c         CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
c         CALL writeg1d(ngrid,1,ps,'ps','Pa')
         
c         CALL writeg1d(ngrid,nlayer,zt,'T','K')
c        CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
c        CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
c        CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')

! THERMALS STUFF 1D
         if(calltherm) then

        call WRITEDIAGFI(ngrid,'lmax_th',
     &              'hauteur du thermique','point',
     &                         0,lmax_th_out)
        call WRITEDIAGFI(ngrid,'zmax_th',
     &              'hauteur du thermique','m',
     &                         0,zmax_th)
        call WRITEDIAGFI(ngrid,'hfmax_th',
     &              'maximum TH heat flux','K.m/s',
     &                         0,hfmax_th)
        call WRITEDIAGFI(ngrid,'wstar',
     &              'maximum TH vertical velocity','m/s',
     &                         0,wstar)

         co2col(:)=0.
         if (tracer) then
         do l=1,nlayer
           do ig=1,ngrid
             co2col(ig)=co2col(ig)+
     &                  zq(ig,l,1)*(zplev(ig,l)-zplev(ig,l+1))/g
         enddo
         enddo

         end if
         call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass'          &
     &                                      ,'kg/m-2',0,co2col)
         endif ! of if (calltherm)

         call WRITEDIAGFI(ngrid,'w','vertical velocity'                 &
     &                              ,'m/s',1,pw)
         call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
         call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
     &                  tsurf)
         call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
         call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)

         call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
         call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
         call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho)
!         call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate",              &
!     &              "K.s-1",1,dtrad/zpopsk)
!        call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
!     &                   'w.m-2',1,zdtsw/zpopsk)
!        call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
!     &                   'w.m-2',1,zdtlw/zpopsk)
         call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
     &                                   ,"kg.m-2",0,co2ice)

! or output in diagfi.nc (for testphys1d)
         call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps)
         call WRITEDIAGFI(ngrid,'temp','Temperature',
     &                       'K',1,zt)
     
         if(tracer) then
c           CALL writeg1d(ngrid,1,tau,'tau','SI')
            do iq=1,nq
c              CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg') 
               call WRITEDIAGFI(ngrid,trim(noms(iq)),
     &              trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
            end do
           if (doubleq) then
             call WRITEDIAGFI(ngrid,'rdust','rdust',
     &                        'm',1,rdust)
           endif
           if (water.AND.tifeedback) then
             call WRITEDIAGFI(ngrid,"soiltemp",
     &                              "Soil temperature","K",
     &                              1,tsoil)
             call WRITEDIAGFI(ngrid,'soilti',
     &                       'Soil Thermal Inertia',
     &                       'J.s-1/2.m-2.K-1',1,inertiesoil)
           endif
         end if
         
cccccccccccccccccc scavenging & water outputs 1D TN ccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         IF (water) THEN

           if (.not.activice) then

             tauTES=0
             do l=1,nlayer
               Qabsice = min(
     &             max(0.4e6*rice(1,l)*(1.+nuice_ref)-0.05 ,0.),1.2
     &                      )
               opTES(1,l)= 0.75 * Qabsice *
     &             zq(1,l,igcm_h2o_ice) *
     &             (zplev(1,l) - zplev(1,l+1)) / g
     &             / (rho_ice * rice(1,l) * (1.+nuice_ref))
               tauTES=tauTES+ opTES(1,l)
             enddo
             CALL WRITEDIAGFI(ngrid,'tauTESap',
     &                         'tau abs 825 cm-1',
     &                         '',0,tauTES)
           else

             CALL WRITEDIAGFI(ngrid,'tauTES',
     &                         'tau abs 825 cm-1',
     &                         '',0,taucloudtes)
           endif
     
           mtot = 0
           icetot = 0
           h2otot = qsurf(1,igcm_h2o_ice)

           do l=1,nlayer
             mtot = mtot +  zq(1,l,igcm_h2o_vap) 
     &                 * (zplev(1,l) - zplev(1,l+1)) / g
             icetot = icetot +  zq(1,l,igcm_h2o_ice) 
     &                 * (zplev(1,l) - zplev(1,l+1)) / g
           end do
           h2otot = h2otot+mtot+icetot

             CALL WRITEDIAGFI(ngrid,'h2otot',
     &                         'h2otot',
     &                         'kg/m2',0,h2otot)
             CALL WRITEDIAGFI(ngrid,'mtot',
     &                         'mtot',
     &                         'kg/m2',0,mtot)
             CALL WRITEDIAGFI(ngrid,'icetot',
     &                         'icetot',
     &                         'kg/m2',0,icetot)

           if (scavenging) then

             rave = 0
             do l=1,nlayer
cccc Column integrated effective ice radius
cccc is weighted by total ice surface area (BETTER) 
             rave = rave + tauscaling(1) *
     &              zq(1,l,igcm_ccn_number) *
     &              (zplev(1,l) - zplev(1,l+1)) / g * 
     &              rice(1,l) * rice(1,l)*  (1.+nuice_ref)
             enddo
             rave=icetot*0.75/max(rave*pi*rho_ice,1.e-30) ! surface weight

              Nccntot= 0
              call watersat(ngrid*nlayer,zt,zplay,zqsat)
               do l=1,nlayer
                   Nccntot = Nccntot + 
     &              zq(1,l,igcm_ccn_number)*tauscaling(1)
     &              *(zplev(1,l) - zplev(1,l+1)) / g
                   satu(1,l) = zq(1,l,igcm_h2o_vap)/zqsat(1,l)
                   satu(1,l) = (max(satu(1,l),float(1))-1)
!     &                      * zq(1,l,igcm_h2o_vap) *
!     &                      (zplev(1,l) - zplev(1,l+1)) / g
               enddo
               call WRITEDIAGFI(ngrid,"satu","vap in satu","kg/kg",1,
     &                    satu)
               CALL WRITEDIAGFI(ngrid,'Nccntot',
     &                         'Nccntot',
     &                         'nbr/m2',0,Nccntot)

             call WRITEDIAGFI(ngrid,'zdqsed_dustq'
     & ,'sedimentation q','kg.m-2.s-1',1,zdqsed(1,:,igcm_dust_mass))
             call WRITEDIAGFI(ngrid,'zdqsed_dustN'
     &,'sedimentation N','Nbr.m-2.s-1',1,
     &                                 zdqsed(1,:,igcm_dust_number))

           else ! of if (scavenging)

cccc Column integrated effective ice radius 
cccc is weighted by total ice mass         (LESS GOOD)
             rave = 0
             do l=1,nlayer
               rave = rave + zq(1,l,igcm_h2o_ice)
     &              * (zplev(1,l) - zplev(1,l+1)) / g
     &              *  rice(1,l) * (1.+nuice_ref)
             enddo
             rave=max(rave/max(icetot,1.e-30),1.e-30) ! mass weight
           endif ! of if (scavenging)

           CALL WRITEDIAGFI(ngrid,'reffice',
     &                      'reffice',
     &                      'm',0,rave)

           do iq=1,nq
             call WRITEDIAGFI(ngrid,trim(noms(iq))//'_s',
     &            trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq))
           end do
     
        call WRITEDIAGFI(ngrid,'zdqcloud_ice','cloud ice',
     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice))
        call WRITEDIAGFI(ngrid,'zdqcloud_vap','cloud vap',
     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_vap))
        call WRITEDIAGFI(ngrid,'zdqcloud','cloud ice',
     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)
     &                          +zdqcloud(1,:,igcm_h2o_vap))
     
          call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
     &                    rice)
         ENDIF ! of IF (water)
         
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


         zlocal(1)=-log(zplay(1,1)/zplev(1,1))* Rnew(1,1)*zt(1,1)/g

         do l=2,nlayer-1
            tmean=zt(1,l)
            if(zt(1,l).ne.zt(1,l-1))
     &        tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
              zlocal(l)= zlocal(l-1)
     &        -log(zplay(1,l)/zplay(1,l-1))*rnew(1,l)*tmean/g
         enddo
         zlocal(nlayer)= zlocal(nlayer-1)-
     &                   log(zplay(1,nlayer)/zplay(1,nlayer-1))*
     &                   rnew(1,nlayer)*tmean/g
#endif

      END IF       ! if(ngrid.ne.1)

      icount=icount+1
      RETURN
      END