Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 1544)
+++ trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 1549)
@@ -27,4 +27,5 @@
       use mod_interface_dyn_phys, only: init_interface_dyn_phys
       USE phys_state_var_init_mod, ONLY: phys_state_var_init
+      USE physiq_mod, ONLY: physiq
       IMPLICIT NONE
 
@@ -95,5 +96,5 @@
       REAL du(nlayer),dv(nlayer),dtemp(nlayer)
       REAL dudyn(nlayer),dvdyn(nlayer),dtempdyn(nlayer)
-      REAL dpsurf    
+      REAL dpsurf(1)    
       REAL,ALLOCATABLE :: dq(:,:)
       REAL,ALLOCATABLE :: dqdyn(:,:)
@@ -792,5 +793,5 @@
 c       ----------------------------------------------------------
 
-           psurf=psurf+dtphys*dpsurf   ! surface pressure change
+           psurf=psurf+dtphys*dpsurf(1)   ! surface pressure change
            DO ilevel=1,nlevel
              plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
Index: trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 1544)
+++ 	(revision )
@@ -1,2642 +1,0 @@
-      SUBROUTINE physiq(
-     $            ngrid,nlayer,nq
-     $            ,firstcall,lastcall
-     $            ,pday,ptime,ptimestep
-     $            ,pplev,pplay,pphi
-     $            ,pu,pv,pt,pq
-     $            ,flxw
-     $            ,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 geometry_mod, only: longitude, latitude, cell_area
-      use comgeomfi_h, only: 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 time_phylmdz_mod, only: iphysiq, day_step, ecritstart, daysec
-      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
-      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 planetwide_mod
-      use phyredem, only: physdem0, physdem1
-      use eofdump_mod, only: eofdump
-      USE comvert_mod, ONLY: ap,bp,aps,bps
-#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    flxw(ngrid,nlayer)      vertical mass flux (kg/s) at layer lower boundary
-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 "chimiedata.h"
-#include "netcdf.inc"
-
-c Arguments :
-c -----------
-
-c   inputs:
-c   -------
-      INTEGER,INTENT(in) :: ngrid ! number of atmospheric columns
-      INTEGER,INTENT(in) :: nlayer ! number of atmospheric layers
-      INTEGER,INTENT(in) :: nq ! number of tracers
-      LOGICAL,INTENT(in) :: firstcall ! signals first call to physics
-      LOGICAL,INTENT(in) :: lastcall ! signals last call to physics
-      REAL,INTENT(in) :: pday ! number of elapsed sols since reference Ls=0
-      REAL,INTENT(in) :: ptime ! "universal time", given as fraction of sol (e.g.: 0.5 for noon)
-      REAL,INTENT(in) :: ptimestep ! physics timestep (s)
-      REAL,INTENT(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
-      REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
-      REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential at mid-layer (m2s-2)
-      REAL,INTENT(in) :: pu(ngrid,nlayer) ! zonal wind component (m/s)
-      REAL,INTENT(in) :: pv(ngrid,nlayer) ! meridional wind component (m/s)
-      REAL,INTENT(in) :: pt(ngrid,nlayer) ! temperature (K)
-      REAL,INTENT(in) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
-      REAL,INTENT(in) :: flxw(ngrid,nlayer) ! vertical mass flux (ks/s)
-                                            ! at lower boundary of layer
-
-c   outputs:
-c   --------
-c     physical tendencies
-      REAL,INTENT(out) :: pdu(ngrid,nlayer) ! zonal wind tendency (m/s/s)
-      REAL,INTENT(out) :: pdv(ngrid,nlayer) ! meridional wind tendency (m/s/s)
-      REAL,INTENT(out) :: pdt(ngrid,nlayer) ! temperature tendency (K/s)
-      REAL,INTENT(out) :: pdq(ngrid,nlayer,nq) ! tracer tendencies (../kg/s)
-      REAL,INTENT(out) :: pdpsrf(ngrid) ! surface pressure tendency (Pa/s)
-      LOGICAL,INTENT(out) :: tracerdyn ! signal to the dynamics to advect tracers or not
-
-
-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.
-#ifdef DUSTSTORM
-      REAL pq_tmp(ngrid, nlayer, 2) ! To compute tendencies due the dust bomb
-#endif
-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 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 photochemistry
-      logical :: asis             ! true  : adaptative semi-implicit symmetric (asis) chemical solver
-                                  ! false : euler backward chemical solver
-      REAL surfdust(ngrid,nlayer) ! dust surface area (m2/m3, if photochemistry)
-      REAL surfice(ngrid,nlayer)  !  ice surface area (m2/m3, if photochemistry)
-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 dsodust(ngrid,nlayer)   ! density-scaled opacity (in infrared) 
-      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)
-      REAL zh(ngrid,nlayer)      ! potential temperature (K)
-      REAL pw(ngrid,nlayer) ! vertical velocity (m/s) (>0 when downwards)
-      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 mass(nq)             ! global mass of tracers (g)
-      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",longitude,latitude,
-     &                 nsoilmx,ngrid,nlayer,nq,
-     &                 ptimestep,pday,time_phys,cell_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
-#ifdef DUSTSTORM
-      pq_tmp(:,:,:)=0
-#endif
-      igout=ngrid/2+1 
-
-
-      zday=pday+ptime ! compute time, in sols (and fraction thereof)
-
-c     Compute Solar Longitude (Ls) :
-c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      if (season) then
-         call solarlong(zday,zls)
-      else
-         call solarlong(float(day_ini),zls)
-      end if
-
-c     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
-
-      ! Compute vertical velocity (m/s) from vertical mass flux
-      ! w = F / (rho*area) and rho = P/(r*T)
-      ! but first linearly interpolate mass flux to mid-layers
-      do l=1,nlayer-1
-       pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
-      enddo
-      pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
-      do l=1,nlayer
-       pw(1:ngrid,l)=(pw(1:ngrid,l)*r*pt(1:ngrid,l)) / 
-     &               (pplay(1:ngrid,l)*cell_area(1:ngrid))
-       ! NB: here we use r and not rnew since this diagnostic comes
-       ! from the dynamics
-      enddo
-
-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,latitude,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
-
-#ifdef DUSTSTORM
-!! specific case: save the quantity of dust before adding perturbation
-       if (firstcall) then
-        pq_tmp(1:ngrid,1:nlayer,1)=pq(1:ngrid,1:nlayer,igcm_dust_mass)
-        pq_tmp(1:ngrid,1:nlayer,2)=pq(1:ngrid,1:nlayer,igcm_dust_number)
-       endif
-#endif
-         
-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,dsodust,tauscaling,taucloudtes,rdust,rice,
-     $     nuice,co2ice)
-
-
-#ifdef DUSTSTORM
-!! specific case: compute the added quantity of dust for perturbation
-       if (firstcall) then
-         pdq(1:ngrid,1:nlayer,igcm_dust_mass)= 
-     &      pdq(1:ngrid,1:nlayer,igcm_dust_mass)
-     &      - pq_tmp(1:ngrid,1:nlayer,1)
-     &      + pq(1:ngrid,1:nlayer,igcm_dust_mass)
-         pdq(1:ngrid,1:nlayer,igcm_dust_number)=
-     &      pdq(1:ngrid,1:nlayer,igcm_dust_number)
-     &      - pq_tmp(1:ngrid,1:nlayer,2)
-     &      + pq(1:ngrid,1:nlayer,igcm_dust_number)
-       endif
-#endif
-
-c          Outputs for basic check (middle of domain)
-c          ------------------------------------------
-           write(*,'("Ls =",f11.6," check lat =",f10.6,
-     &               " lon =",f11.6)')
-     &           zls*180./pi,latitude(igout)*180/pi,
-     &                       longitude(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.*longitude(ig)/pi/15.
-                 sl_ra  = pi*(1.0-sl_lct/12.)
-                 sl_lat = 180.*latitude(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
-
-            asis = .false.
-
-            if (asis) then
-            call calchim_asis(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)
-            else
-            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)
-            end if
-
-           ! 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.*longitude/pi, 180.*latitude/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     &        ,latitude(ig)*180./pi, longitude(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)
-              ! also store vmr_ice*rice for better diagnostics of rice
-               vmr(1:ngrid,1:nlayer)=vmr(1:ngrid,1:nlayer)*
-     &                               rice(1:ngrid,1:nlayer)
-               call wstats(ngrid,"vmr_h2oice_rice",
-     &                "H2O ice mixing ratio times ice particule size",
-     &                    "(mol/mol)*m",
-     &                    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") then
-
-!                   volume mixing ratio
-
-                    vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
-     &                            *mmean(1:ngrid,1:nlayer)/mmol(iq)
-
-                    call wstats(ngrid,"vmr_"//trim(noms(iq)),
-     $                        "Volume mixing ratio","mol/mol",3,vmr)
-                    if ((noms(iq).eq."o")
-     $             .or. (noms(iq).eq."co2")
-     $             .or. (noms(iq).eq."o3")
-     $             .or. (noms(iq).eq."ar")
-     $             .or. (noms(iq).eq."o2")
-     $             .or. (noms(iq).eq."h2o_vap") ) then
-                      call writediagfi(ngrid,"vmr_"//trim(noms(iq)),
-     $                         "Volume mixing ratio","mol/mol",3,vmr)
-                    end if
-
-!                   number density (molecule.cm-3)
-
-                    rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
-     &                          *rho(1:ngrid,1:nlayer)*n_avog/
-     &                           (1000*mmol(iq))
-
-!                   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)
-
-!                   vertical column (molecule.cm-2)
-
-                    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))  
-                    call writediagfi(ngrid,"c_"//trim(noms(iq)),
-     $                          "column","mol cm-2",2,colden(1,iq))
-
-!                   global mass (g)
-               
-                    call planetwide_sumval(colden(:,iq)/6.022e23
-     $                            *mmol(iq)*1.e4*cell_area(:),mass(iq))
-
-                    call writediagfi(ngrid,"mass_"//trim(noms(iq)),
-     $                              "global mass","g",0,mass(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"        ""
-      if (activice) then
-        comm_TAU_ICE(1:ngrid) = taucloudtes(1:ngrid)
-      else
-        comm_TAU_ICE(1:ngrid) = tauTES(1:ngrid)
-      endif
-      !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)
-            call WRITEDIAGFI(ngrid,'saturation',
-     &           'h2o vap saturation ratio','dimless',3,satu)
-            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)
-             call WRITEDIAGFI(ngrid,'dsodust',
-     &                        'dust density scaled opacity',
-     &                        'm2.kg-1',3,dsodust)
-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)
-
-         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)
-        call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
-     &                   'w.m-2',1,zdtsw)
-        call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
-     &                   'w.m-2',1,zdtlw)
-         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
Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 1549)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 1549)
@@ -0,0 +1,2650 @@
+      MODULE physiq_mod
+
+      IMPLICIT NONE
+
+      CONTAINS
+
+      SUBROUTINE physiq(
+     $            ngrid,nlayer,nq
+     $            ,firstcall,lastcall
+     $            ,pday,ptime,ptimestep
+     $            ,pplev,pplay,pphi
+     $            ,pu,pv,pt,pq
+     $            ,flxw
+     $            ,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 geometry_mod, only: longitude, latitude, cell_area
+      use comgeomfi_h, only: 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 time_phylmdz_mod, only: iphysiq, day_step, ecritstart, daysec
+      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
+      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 planetwide_mod
+      use phyredem, only: physdem0, physdem1
+      use eofdump_mod, only: eofdump
+      USE comvert_mod, ONLY: ap,bp,aps,bps
+#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    flxw(ngrid,nlayer)      vertical mass flux (kg/s) at layer lower boundary
+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 "chimiedata.h"
+#include "netcdf.inc"
+
+c Arguments :
+c -----------
+
+c   inputs:
+c   -------
+      INTEGER,INTENT(in) :: ngrid ! number of atmospheric columns
+      INTEGER,INTENT(in) :: nlayer ! number of atmospheric layers
+      INTEGER,INTENT(in) :: nq ! number of tracers
+      LOGICAL,INTENT(in) :: firstcall ! signals first call to physics
+      LOGICAL,INTENT(in) :: lastcall ! signals last call to physics
+      REAL,INTENT(in) :: pday ! number of elapsed sols since reference Ls=0
+      REAL,INTENT(in) :: ptime ! "universal time", given as fraction of sol (e.g.: 0.5 for noon)
+      REAL,INTENT(in) :: ptimestep ! physics timestep (s)
+      REAL,INTENT(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
+      REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
+      REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential at mid-layer (m2s-2)
+      REAL,INTENT(in) :: pu(ngrid,nlayer) ! zonal wind component (m/s)
+      REAL,INTENT(in) :: pv(ngrid,nlayer) ! meridional wind component (m/s)
+      REAL,INTENT(in) :: pt(ngrid,nlayer) ! temperature (K)
+      REAL,INTENT(in) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
+      REAL,INTENT(in) :: flxw(ngrid,nlayer) ! vertical mass flux (ks/s)
+                                            ! at lower boundary of layer
+
+c   outputs:
+c   --------
+c     physical tendencies
+      REAL,INTENT(out) :: pdu(ngrid,nlayer) ! zonal wind tendency (m/s/s)
+      REAL,INTENT(out) :: pdv(ngrid,nlayer) ! meridional wind tendency (m/s/s)
+      REAL,INTENT(out) :: pdt(ngrid,nlayer) ! temperature tendency (K/s)
+      REAL,INTENT(out) :: pdq(ngrid,nlayer,nq) ! tracer tendencies (../kg/s)
+      REAL,INTENT(out) :: pdpsrf(ngrid) ! surface pressure tendency (Pa/s)
+      LOGICAL,INTENT(out) :: tracerdyn ! signal to the dynamics to advect tracers or not
+
+
+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.
+#ifdef DUSTSTORM
+      REAL pq_tmp(ngrid, nlayer, 2) ! To compute tendencies due the dust bomb
+#endif
+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 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 photochemistry
+      logical :: asis             ! true  : adaptative semi-implicit symmetric (asis) chemical solver
+                                  ! false : euler backward chemical solver
+      REAL surfdust(ngrid,nlayer) ! dust surface area (m2/m3, if photochemistry)
+      REAL surfice(ngrid,nlayer)  !  ice surface area (m2/m3, if photochemistry)
+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 dsodust(ngrid,nlayer)   ! density-scaled opacity (in infrared) 
+      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)
+      REAL zh(ngrid,nlayer)      ! potential temperature (K)
+      REAL pw(ngrid,nlayer) ! vertical velocity (m/s) (>0 when downwards)
+      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 mass(nq)             ! global mass of tracers (g)
+      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",longitude,latitude,
+     &                 nsoilmx,ngrid,nlayer,nq,
+     &                 ptimestep,pday,time_phys,cell_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
+#ifdef DUSTSTORM
+      pq_tmp(:,:,:)=0
+#endif
+      igout=ngrid/2+1 
+
+
+      zday=pday+ptime ! compute time, in sols (and fraction thereof)
+
+c     Compute Solar Longitude (Ls) :
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      if (season) then
+         call solarlong(zday,zls)
+      else
+         call solarlong(float(day_ini),zls)
+      end if
+
+c     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
+
+      ! Compute vertical velocity (m/s) from vertical mass flux
+      ! w = F / (rho*area) and rho = P/(r*T)
+      ! but first linearly interpolate mass flux to mid-layers
+      do l=1,nlayer-1
+       pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
+      enddo
+      pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
+      do l=1,nlayer
+       pw(1:ngrid,l)=(pw(1:ngrid,l)*r*pt(1:ngrid,l)) / 
+     &               (pplay(1:ngrid,l)*cell_area(1:ngrid))
+       ! NB: here we use r and not rnew since this diagnostic comes
+       ! from the dynamics
+      enddo
+
+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,latitude,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
+
+#ifdef DUSTSTORM
+!! specific case: save the quantity of dust before adding perturbation
+       if (firstcall) then
+        pq_tmp(1:ngrid,1:nlayer,1)=pq(1:ngrid,1:nlayer,igcm_dust_mass)
+        pq_tmp(1:ngrid,1:nlayer,2)=pq(1:ngrid,1:nlayer,igcm_dust_number)
+       endif
+#endif
+         
+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,dsodust,tauscaling,taucloudtes,rdust,rice,
+     $     nuice,co2ice)
+
+
+#ifdef DUSTSTORM
+!! specific case: compute the added quantity of dust for perturbation
+       if (firstcall) then
+         pdq(1:ngrid,1:nlayer,igcm_dust_mass)= 
+     &      pdq(1:ngrid,1:nlayer,igcm_dust_mass)
+     &      - pq_tmp(1:ngrid,1:nlayer,1)
+     &      + pq(1:ngrid,1:nlayer,igcm_dust_mass)
+         pdq(1:ngrid,1:nlayer,igcm_dust_number)=
+     &      pdq(1:ngrid,1:nlayer,igcm_dust_number)
+     &      - pq_tmp(1:ngrid,1:nlayer,2)
+     &      + pq(1:ngrid,1:nlayer,igcm_dust_number)
+       endif
+#endif
+
+c          Outputs for basic check (middle of domain)
+c          ------------------------------------------
+           write(*,'("Ls =",f11.6," check lat =",f10.6,
+     &               " lon =",f11.6)')
+     &           zls*180./pi,latitude(igout)*180/pi,
+     &                       longitude(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.*longitude(ig)/pi/15.
+                 sl_ra  = pi*(1.0-sl_lct/12.)
+                 sl_lat = 180.*latitude(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
+
+            asis = .false.
+
+            if (asis) then
+            call calchim_asis(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)
+            else
+            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)
+            end if
+
+           ! 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.*longitude/pi, 180.*latitude/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     &        ,latitude(ig)*180./pi, longitude(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)
+              ! also store vmr_ice*rice for better diagnostics of rice
+               vmr(1:ngrid,1:nlayer)=vmr(1:ngrid,1:nlayer)*
+     &                               rice(1:ngrid,1:nlayer)
+               call wstats(ngrid,"vmr_h2oice_rice",
+     &                "H2O ice mixing ratio times ice particule size",
+     &                    "(mol/mol)*m",
+     &                    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") then
+
+!                   volume mixing ratio
+
+                    vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
+     &                            *mmean(1:ngrid,1:nlayer)/mmol(iq)
+
+                    call wstats(ngrid,"vmr_"//trim(noms(iq)),
+     $                        "Volume mixing ratio","mol/mol",3,vmr)
+                    if ((noms(iq).eq."o")
+     $             .or. (noms(iq).eq."co2")
+     $             .or. (noms(iq).eq."o3")
+     $             .or. (noms(iq).eq."ar")
+     $             .or. (noms(iq).eq."o2")
+     $             .or. (noms(iq).eq."h2o_vap") ) then
+                      call writediagfi(ngrid,"vmr_"//trim(noms(iq)),
+     $                         "Volume mixing ratio","mol/mol",3,vmr)
+                    end if
+
+!                   number density (molecule.cm-3)
+
+                    rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
+     &                          *rho(1:ngrid,1:nlayer)*n_avog/
+     &                           (1000*mmol(iq))
+
+!                   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)
+
+!                   vertical column (molecule.cm-2)
+
+                    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))  
+                    call writediagfi(ngrid,"c_"//trim(noms(iq)),
+     $                          "column","mol cm-2",2,colden(1,iq))
+
+!                   global mass (g)
+               
+                    call planetwide_sumval(colden(:,iq)/6.022e23
+     $                            *mmol(iq)*1.e4*cell_area(:),mass(iq))
+
+                    call writediagfi(ngrid,"mass_"//trim(noms(iq)),
+     $                              "global mass","g",0,mass(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"        ""
+      if (activice) then
+        comm_TAU_ICE(1:ngrid) = taucloudtes(1:ngrid)
+      else
+        comm_TAU_ICE(1:ngrid) = tauTES(1:ngrid)
+      endif
+      !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)
+            call WRITEDIAGFI(ngrid,'saturation',
+     &           'h2o vap saturation ratio','dimless',3,satu)
+            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)
+             call WRITEDIAGFI(ngrid,'dsodust',
+     &                        'dust density scaled opacity',
+     &                        'm2.kg-1',3,dsodust)
+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)
+
+         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)
+        call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
+     &                   'w.m-2',1,zdtsw)
+        call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
+     &                   'w.m-2',1,zdtlw)
+         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
+
+      END SUBROUTINE physiq
+
+      END MODULE physiq_mod
