Index: trunk/LMDZ.GENERIC/README
===================================================================
--- trunk/LMDZ.GENERIC/README	(revision 747)
+++ trunk/LMDZ.GENERIC/README	(revision 749)
@@ -750,2 +750,5 @@
 - Bug fix in physiq : the size of OLR_nu is L_NSPECTI and not L_NSPECTV
 - A more robust aerosol_mod + iniaerosol : problems with ifort+parallel solved, while still OK with other compilers and seq runs.
+
+== 01/08/2012 == AS
+- Changed name of physiq to physiqg to allow for easy interfacing with latest LMDz dynamical core.
Index: trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 747)
+++ 	(revision )
@@ -1,1760 +1,0 @@
-      subroutine physiq(ngrid,nlayer,nq,   &
-                  firstcall,lastcall,      &
-                  pday,ptime,ptimestep,    &
-                  pplev,pplay,pphi,        &
-                  pu,pv,pt,pq,             &
-                  pw,                      &
-                  pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
- 
-      use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind
-      use watercommon_h, only : RLVTT, Psat_water
-      use gases_h
-      use radcommon_h, only: sigma
-      use radii_mod, only: h2o_reffrad, co2_reffrad
-      use aerosol_mod
-      implicit none
-
-
-!==================================================================
-!     
-!     Purpose
-!     -------
-!     Central subroutine for all the physics parameterisations in the
-!     universal model. Originally adapted from the Mars LMDZ model.
-!
-!     The model can be run without or with tracer transport
-!     depending on the value of "tracer" in file "callphys.def".
-!
-!
-!   It includes:
-!
-!      1.  Initialization:
-!      1.1 Firstcall initializations
-!      1.2 Initialization for every call to physiq
-!      1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
-!      2. Compute radiative transfer tendencies
-!         (longwave and shortwave).
-!      4. Vertical diffusion (turbulent mixing):
-!      5. Convective adjustment
-!      6. Condensation and sublimation of gases (currently just CO2).
-!      7. TRACERS
-!       7a. water and water ice
-!       7c. other schemes for tracer transport (lifting, sedimentation)
-!       7d. updates (pressure variations, surface budget)
-!      9. Surface and sub-surface temperature calculations
-!     10. Write outputs :
-!           - "startfi", "histfi" if it's time
-!           - Saving statistics if "callstats = .true."
-!           - Output any needed variables in "diagfi" 
-!     10. Diagnostics: mass conservation of tracers, radiative energy balance etc.
-!
-!   arguments
-!   ---------
-!
-!   input
-!   -----
-!    ecri                  period (in dynamical timestep) to write output
-!    ngrid                 Size of the horizontal grid.
-!                          All internal loops are performed on that grid.
-!    nlayer                Number of vertical layers.
-!    nq                    Number of advected fields
-!    firstcall             True at the first call
-!    lastcall              True at the last call
-!    pday                  Number of days counted from the North. Spring
-!                          equinoxe.
-!    ptime                 Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
-!    ptimestep             timestep (s)
-!    pplay(ngrid,nlayer)   Pressure at the middle of the layers (Pa)
-!    pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
-!    pphi(ngrid,nlayer)    Geopotential at the middle of the layers (m2s-2)
-!    pu(ngrid,nlayer)      u component of the wind (ms-1)
-!    pv(ngrid,nlayer)      v component of the wind (ms-1)
-!    pt(ngrid,nlayer)      Temperature (K)
-!    pq(ngrid,nlayer,nq)   Advected fields
-!    pudyn(ngrid,nlayer)    \ 
-!    pvdyn(ngrid,nlayer)     \ Dynamical temporal derivative for the
-!    ptdyn(ngrid,nlayer)     / corresponding variables
-!    pqdyn(ngrid,nlayer,nq) /
-!    pw(ngrid,?)           vertical velocity
-!
-!   output
-!   ------
-!
-!    pdu(ngrid,nlayermx)        \
-!    pdv(ngrid,nlayermx)         \  Temporal derivative of the corresponding
-!    pdt(ngrid,nlayermx)         /  variables due to physical processes.
-!    pdq(ngrid,nlayermx)        /
-!    pdpsrf(ngrid)             /
-!    tracerdyn                 call tracer in dynamical part of GCM ?
-!
-!
-!     Authors
-!     -------
-!           Frederic Hourdin	15/10/93
-!           Francois Forget	1994
-!           Christophe Hourdin	02/1997 
-!           Subroutine completely rewritten by F. Forget (01/2000)
-!           Water ice clouds: Franck Montmessin (update 06/2003)
-!           Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
-!           New correlated-k radiative scheme: R. Wordsworth (2009)
-!           Many specifically Martian subroutines removed: R. Wordsworth (2009)       
-!           Improved water cycle: R. Wordsworth / B. Charnay (2010)
-!           To F90: R. Wordsworth (2010)
-!           New turbulent diffusion scheme: J. Leconte (2012)
-!           Loops converted to F90 matrix format: J. Leconte (2012)
-!     
-!==================================================================
-
-
-!    0.  Declarations :
-!    ------------------
-
-#include "dimensions.h"
-#include "dimphys.h"
-#include "comgeomfi.h"
-#include "surfdat.h"
-#include "comsoil.h"
-#include "comdiurn.h"
-#include "callkeys.h"
-#include "comcstfi.h"
-#include "planete.h"
-#include "comsaison.h"
-#include "control.h"
-#include "tracer.h"
-#include "watercap.h"
-#include "netcdf.inc"
-
-! Arguments :
-! -----------
-
-!   inputs:
-!   -------
-      integer ngrid,nlayer,nq
-      real ptimestep
-      real pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
-      real pphi(ngridmx,nlayer)
-      real pu(ngridmx,nlayer),pv(ngridmx,nlayer)
-      real pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
-      real pw(ngridmx,nlayer)        ! pvervel transmitted by dyn3d
-      real zh(ngridmx,nlayermx)      ! potential temperature (K)
-      logical firstcall,lastcall
-
-      real pday
-      real ptime 
-      logical tracerdyn
-
-!   outputs:
-!   --------
-!     physical tendencies
-      real pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
-      real pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
-      real pdpsrf(ngridmx) ! surface pressure tendency
-
-
-! Local saved variables:
-! ----------------------
-!     aerosol (dust or ice) extinction optical depth  at reference wavelength 
-!     "longrefvis" set in dimradmars.h , for one of the "naerkind"  kind of
-!      aerosol optical properties:
-!      real aerosol(ngridmx,nlayermx,naerkind) 
-!     this is now internal to callcorrk and hence no longer needed here
-
-      integer day_ini                ! Initial date of the run (sol since Ls=0) 
-      integer icount                 ! counter of calls to physiq during the run.
-      real tsurf(ngridmx)            ! Surface temperature (K)
-      real tsoil(ngridmx,nsoilmx)    ! sub-surface temperatures (K)
-      real albedo(ngridmx)           ! Surface albedo
-
-      real albedo0(ngridmx)          ! Surface albedo
-      integer rnat(ngridmx)          ! added by BC
-      save rnat
-
-      real emis(ngridmx)             ! Thermal IR surface emissivity
-      real dtrad(ngridmx,nlayermx)   ! Net atm. radiative heating rate (K.s-1)
-      real fluxrad_sky(ngridmx)      ! rad. flux from sky absorbed by surface (W.m-2)
-      real fluxrad(ngridmx)          ! Net radiative surface flux (W.m-2)
-      real capcal(ngridmx)           ! surface heat capacity (J m-2 K-1)
-      real fluxgrd(ngridmx)          ! surface conduction flux (W.m-2)
-      real qsurf(ngridmx,nqmx)       ! tracer on surface (e.g. kg.m-2)
-      real q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy 
-
-      save day_ini, icount
-      save tsurf,tsoil
-      save albedo0,albedo,emis,q2
-      save capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
-
-! Local variables :
-! -----------------
-
-!     aerosol (dust or ice) extinction optical depth  at reference wavelength 
-!     for the "naerkind" optically active aerosols:
-      real aerosol(ngridmx,nlayermx,naerkind) 
-
-      character*80 fichier 
-      integer l,ig,ierr,iq,iaer
-
-      real fluxsurf_lw(ngridmx)      ! incident LW (IR) surface flux (W.m-2)
-      real fluxsurf_sw(ngridmx)      ! incident SW (stellar) surface flux (W.m-2)
-      real fluxtop_lw(ngridmx)       ! Outgoing LW (IR) flux to space (W.m-2)
-      real fluxabs_sw(ngridmx)       ! Absorbed SW (stellar) flux (W.m-2)
-
-      real fluxtop_dn(ngridmx)
-      real fluxdyn(ngridmx)          ! horizontal heat transport by dynamics      
-      real OLR_nu(ngridmx,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
-      real OSR_nu(ngridmx,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
- 
-      real,save :: sensibFlux(ngridmx)   ! turbulent flux given by the atm to the surface
- 
-      save fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu
-
-
-      real zls                       ! solar longitude (rad)
-      real zday                      ! date (time since Ls=0, in martian days)
-      real zzlay(ngridmx,nlayermx)   ! altitude at the middle of the layers
-      real zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
-      real latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
-
-!     Tendencies due to various processes:
-      real dqsurf(ngridmx,nqmx)
-      real,save :: zdtlw(ngridmx,nlayermx)                    ! (K/s)
-      real,save :: zdtsw(ngridmx,nlayermx)                    ! (K/s)
-
-      real cldtlw(ngridmx,nlayermx)                           ! (K/s) LW heating rate for clear areas
-      real cldtsw(ngridmx,nlayermx)                           ! (K/s) SW heating rate for clear areas
-      real zdtsurf(ngridmx)                                   ! (K/s)
-      real dtlscale(ngridmx,nlayermx)
-      real zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx)  ! (m.s-2)
-      real zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx)         ! (K/s)
-      real zdtdif(ngridmx,nlayermx)			      ! (K/s)
-      real zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx)  ! (m.s-2)
-      real zdhadj(ngridmx,nlayermx)                           ! (K/s)
-      real zdtgw(ngridmx,nlayermx)                            ! (K/s)
-      real zdtmr(ngridmx,nlayermx)
-      real zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx)    ! (m.s-2)
-      real zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
-      real zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
-      real zdumr(ngridmx,nlayermx),zdvmr(ngridmx,nlayermx)
-      real zdtsurfmr(ngridmx)
-      
-      real zdmassmr(ngridmx,nlayermx),zdpsrfmr(ngridmx)
-
-      real zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
-      real zdqevap(ngridmx,nlayermx)
-      real zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
-      real zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
-      real zdqadj(ngridmx,nlayermx,nqmx)
-      real zdqc(ngridmx,nlayermx,nqmx)
-      real zdqmr(ngridmx,nlayermx,nqmx),zdqsurfmr(ngridmx,nqmx)
-      real zdqlscale(ngridmx,nlayermx,nqmx)
-      real zdqslscale(ngridmx,nqmx)
-      real zdqchim(ngridmx,nlayermx,nqmx)
-      real zdqschim(ngridmx,nqmx)
-
-      real zdteuv(ngridmx,nlayermx)    ! (K/s)
-      real zdtconduc(ngridmx,nlayermx) ! (K/s)
-      real zdumolvis(ngridmx,nlayermx)
-      real zdvmolvis(ngridmx,nlayermx)
-      real zdqmoldiff(ngridmx,nlayermx,nqmx)
-
-!     Local variables for local calculations:
-      real zflubid(ngridmx)
-      real zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
-      real zdum1(ngridmx,nlayermx)
-      real zdum2(ngridmx,nlayermx)
-      real ztim1,ztim2,ztim3, z1,z2
-      real ztime_fin
-      real zdh(ngridmx,nlayermx)
-      integer length
-      parameter (length=100)
-
-! local variables only used for diagnostics (output in file "diagfi" or "stats")
-! ------------------------------------------------------------------------------
-      real ps(ngridmx), zt(ngridmx,nlayermx)
-      real zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
-      real zq(ngridmx,nlayermx,nqmx)
-      character*2 str2
-      character*5 str5
-      real zdtadj(ngridmx,nlayermx)
-      real zdtdyn(ngridmx,nlayermx),ztprevious(ngridmx,nlayermx)
-      save ztprevious 
-      real reff(ngridmx,nlayermx) ! effective dust radius (used if doubleq=T)
-      real qtot1,qtot2            ! total aerosol mass
-      integer igmin, lmin
-      logical tdiag
-
-      real zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
-      real zstress(ngrid), cd
-      real hco2(nqmx), tmean, zlocal(nlayermx)
-      real vmr(ngridmx,nlayermx) ! volume mixing ratio
-
-      real time_phys
-
-!     reinstated by RW for diagnostic
-      real tau_col(ngridmx)
-      save tau_col
-      
-!     included by RW to reduce insanity of code
-      real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS
-
-!     included by RW to compute tracer column densities
-      real qcol(ngridmx,nqmx)
-
-!     included by RW for H2O precipitation
-      real zdtrain(ngridmx,nlayermx)
-      real zdqrain(ngridmx,nlayermx,nqmx)
-      real zdqsrain(ngridmx)
-      real zdqssnow(ngridmx)
-      real rainout(ngridmx)
-
-!     included by RW for H2O Manabe scheme
-      real dtmoist(ngridmx,nlayermx)
-      real dqmoist(ngridmx,nlayermx,nqmx)
-
-      real qvap(ngridmx,nlayermx)
-      real dqvaplscale(ngridmx,nlayermx)
-      real dqcldlscale(ngridmx,nlayermx)
-      real rneb_man(ngridmx,nlayermx)
-      real rneb_lsc(ngridmx,nlayermx)
-
-!     included by RW to account for surface cooling by evaporation
-      real dtsurfh2olat(ngridmx)
-
-
-!     to test energy conservation (RW+JL)
-      real mass(ngridmx,nlayermx),massarea(ngridmx,nlayermx)
-      real dEtot, dEtots, AtmToSurf_TurbFlux
-      real dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
-      real dEzRadsw(ngridmx,nlayermx),dEzRadlw(ngridmx,nlayermx),dEzdiff(ngridmx,nlayermx)
-      real dEdiffs(ngridmx),dEdiff(ngridmx)
-      real madjdE(ngridmx), lscaledE(ngridmx)
-!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
-      
-      real dItot, dVtot
-
-!     included by BC for evaporation     
-      real qevap(ngridmx,nlayermx,nqmx)
-      real tevap(ngridmx,nlayermx)
-      real dqevap1(ngridmx,nlayermx)
-      real dtevap1(ngridmx,nlayermx)
-
-!     included by BC for hydrology
-      real hice(ngridmx)
-
-!     included by RW to test water conservation (by routine)
-      real h2otot
-      real dWtot, dWtots
-      real h2o_surf_all
-      logical watertest
-      save watertest
-
-!     included by RW for RH diagnostic
-      real qsat(ngridmx,nlayermx), RH(ngridmx,nlayermx), H2Omaxcol(ngridmx),psat_tmp
-
-!     included by RW for hydrology
-      real dqs_hyd(ngridmx,nqmx)
-      real zdtsurf_hyd(ngridmx)
-
-!     included by RW for water cycle conservation tests
-      real icesrf,liqsrf,icecol,vapcol
-
-!     included by BC for double radiative transfer call
-      logical clearsky
-      real zdtsw1(ngridmx,nlayermx)
-      real zdtlw1(ngridmx,nlayermx)
-      real fluxsurf_lw1(ngridmx)     
-      real fluxsurf_sw1(ngridmx)  
-      real fluxtop_lw1(ngridmx)    
-      real fluxabs_sw1(ngridmx)  
-      real tau_col1(ngrid)
-      real OLR_nu1(ngrid,L_NSPECTI)
-      real OSR_nu1(ngrid,L_NSPECTV)
-      real tf, ntf
-
-!     included by BC for cloud fraction computations
-      real,save :: cloudfrac(ngridmx,nlayermx)
-      real,save :: totcloudfrac(ngridmx)
-
-!     included by RW for vdifc water conservation test
-      real nconsMAX
-      real vdifcncons(ngridmx)
-      real cadjncons(ngridmx)
-
-!      double precision qsurf_hist(ngridmx,nqmx)
-      real qsurf_hist(ngridmx,nqmx)
-      save qsurf_hist
-
-!     included by RW for temp convadj conservation test
-      real playtest(nlayermx)
-      real plevtest(nlayermx)
-      real ttest(nlayermx)
-      real qtest(nlayermx)
-      integer igtest
-
-!     included by RW for runway greenhouse 1D study
-      real muvar(ngridmx,nlayermx+1)
-
-!     included by RW for variable H2O particle sizes
-      real,save :: reffrad(ngridmx,nlayermx,naerkind) ! aerosol effective radius (m)
-      real :: reffrad1(ngridmx,nlayermx,naerkind) ! for clearsky call to callcorrk
-      real,save :: nueffrad(ngridmx,nlayermx,naerkind) ! aerosol effective radius variance
-      real :: reffh2oliq(ngridmx,nlayermx) ! liquid water particles effective radius (m)
-      real :: reffh2oice(ngridmx,nlayermx) ! water ice particles effective radius (m)
-      real reffH2O(ngridmx,nlayermx)
-      real reffcol(ngridmx,naerkind)
-
-!     included by RW for sourceevol
-      real, save :: ice_initial(ngridmx)!, isoil
-      real delta_ice,ice_tot
-      real, save :: ice_min(ngridmx)
-      
-      integer num_run
-      logical,save :: ice_update
-      
-!=======================================================================
-      
-
-! 1. Initialisation
-! -----------------
-
-!  1.1   Initialisation only at first call
-!  ---------------------------------------
-      if (firstcall) then
-
-
-!        variables set to 0
-!        ~~~~~~~~~~~~~~~~~~
-         dtrad(:,:) = 0.0
-         fluxrad(:) = 0.0
-         tau_col(:) = 0.0
-         zdtsw(:,:) = 0.0
-         zdtlw(:,:) = 0.0
-
-!        initialize aerosol indexes 
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-            call iniaerosol()
-
-      
-!        initialize tracer names, indexes and properties
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         tracerdyn=tracer
-         if (tracer) then
-            call initracer()
-         endif ! end tracer
-
-! 
-
-!        read startfi (initial parameters)
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         call phyetat0("startfi.nc",0,0,nsoilmx,nq,           &
-               day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,   &
-               cloudfrac,totcloudfrac,hice)
-
-         if (pday.ne.day_ini) then
-           write(*,*) "ERROR in physiq.F90:"
-           write(*,*) "bad synchronization between physics and dynamics"
-           write(*,*) "dynamics day: ",pday
-           write(*,*) "physics day:  ",day_ini
-           stop
-         endif
-
-         write (*,*) 'In physiq day_ini =', day_ini
-
-!        Initialize albedo and orbital calculation
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         call surfini(ngrid,qsurf,albedo0)
-         call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
-
-         albedo(:)=albedo0(:)
-
-         if(tlocked)then
-            print*,'Planet is tidally locked at resonance n=',nres
-            print*,'Make sure you have the right rotation rate!!!'
-         endif
-
-!        initialize soil 
-!        ~~~~~~~~~~~~~~~
-         if (callsoil) then
-            call soil(ngrid,nsoilmx,firstcall,inertiedat,     &
-                ptimestep,tsurf,tsoil,capcal,fluxgrd)
-         else
-            print*,'WARNING! Thermal conduction in the soil turned off'
-            capcal(:)=1.e6
-            fluxgrd(:)=0.
-            if(noradsurf)then
-                  fluxgrd(:)=10.0
-            endif
-            print*,'Flux from ground = ',fluxgrd,' W m^-2'
-         endif
-         icount=1
-
-!        decide whether to update ice at end of run
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         ice_update=.false.
-         if(sourceevol)then
-            open(128,file='num_run',form='formatted')
-            read(128,*) num_run 
-            close(128)
-        
-            if(num_run.ne.0.and.mod(num_run,2).eq.0)then
-            !if(num_run.ne.0.and.mod(num_run,3).eq.0)then
-               print*,'Updating ice at end of this year!'
-               ice_update=.true.
-               ice_min(:)=1.e4
-            endif  
-         endif
-
-!        define surface as continent or ocean
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         rnat(:)=1
-         do ig=1,ngridmx
-!            if(iceball.or.oceanball.or.(inertiedat(ig,1).gt.1.E4))then
-            if(inertiedat(ig,1).gt.1.E4)then
-            rnat(ig)=0
-            endif
-         enddo
-
-         print*,'WARNING! Surface type currently decided by surface inertia'
-         print*,'This should be improved e.g. in newstart.F'
-
-
-!        initialise surface history variable 
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         qsurf_hist(:,:)=qsurf(:,:)
-
-!        initialise variable for dynamical heating diagnostic
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         ztprevious(:,:)=pt(:,:)
-
-!        Set temperature just above condensation temperature (for Early Mars)
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         if(nearco2cond) then
-            write(*,*)' WARNING! Starting at Tcond+1K'
-            do l=1, nlayer
-               do ig=1,ngrid
-                  pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4     &
-                      -pt(ig,l)) / ptimestep
-               enddo
-            enddo
-         endif
-
-         if(meanOLR)then
-            ! to record global radiative balance
-            call system('rm -f rad_bal.out')
-            ! to record global mean/max/min temperatures
-            call system('rm -f tem_bal.out')
-            ! to record global hydrological balance
-            call system('rm -f h2o_bal.out')
-         endif
-
-         watertest=.false.
-         if(water)then
-            ! initialise variables for water cycle
-
-            if(enertest)then
-               watertest = .true.
-            endif
-
-            if(ice_update)then
-               ice_initial(:)=qsurf(:,igcm_h2o_ice)
-            endif
-
-         endif
-         call su_watercycle ! even if we don't have a water cycle, we might
-                            ! need epsi for the wvp definitions in callcorrk.F
-
-      endif        !  (end of "if firstcall")
-
-! ---------------------------------------------------
-! 1.2   Initializations done at every physical timestep:
-! ---------------------------------------------------
-
-      if (ngrid.NE.ngridmx) then
-         print*,'STOP in PHYSIQ'
-         print*,'Probleme de dimensions :'
-         print*,'ngrid     = ',ngrid
-         print*,'ngridmx   = ',ngridmx
-         stop
-      endif
-
-!     Initialize various variables
-!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-      pdu(1:ngridmx,1:nlayermx) = 0.0
-      pdv(1:ngridmx,1:nlayermx) = 0.0
-      if ( .not.nearco2cond ) then
-         pdt(1:ngridmx,1:nlayermx) = 0.0
-      endif
-      pdq(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
-      pdpsrf(1:ngridmx)       = 0.0
-      zflubid(1:ngridmx)      = 0.0
-      zdtsurf(1:ngridmx)      = 0.0
-      dqsurf(1:ngridmx,1:nqmx)= 0.0
-
-      zday=pday+ptime ! compute time, in sols (and fraction thereof)
-
-!     Compute Stellar Longitude (Ls)
-!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      if (season) then
-         call stellarlong(zday,zls)
-      else
-         call stellarlong(float(day_ini),zls)
-      end if
-
-!     Compute geopotential between layers
-!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-      zzlay(1:ngridmx,1:nlayermx)=pphi(1:ngridmx,1:nlayermx)/g
-
-      zzlev(1:ngridmx,1)=0.
-      zzlev(1:ngridmx,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
-
-      do l=2,nlayer
-         do ig=1,ngrid
-            z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
-            z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
-            zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
-         enddo
-      enddo
-!     Potential temperature calculation may not be the same in physiq and dynamic...
-
-!     Compute potential temperature
-!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-      do l=1,nlayer         
-         do ig=1,ngrid 
-            zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
-            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
-            mass(ig,l)  = (pplev(ig,l) - pplev(ig,l+1))/g
-	    massarea(ig,l)=mass(ig,l)*area(ig)
-         enddo
-      enddo
-
-!-----------------------------------------------------------------------
-!    2. Compute radiative tendencies
-!-----------------------------------------------------------------------
-
-      if (callrad) then
-         if( mod(icount-1,iradia).eq.0.or.lastcall) then
-
-!          Compute local stellar zenith angles
-!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-           call orbite(zls,dist_star,declin)
-
-           if (tlocked) then
-              ztim1=SIN(declin)
-              ztim2=COS(declin)*COS(2.*pi*(zday/year_day) - zls*nres)
-              ztim3=-COS(declin)*SIN(2.*pi*(zday/year_day) - zls*nres)
-
-              call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
-               ztim1,ztim2,ztim3,mu0,fract)
-
-           elseif (diurnal) then
-               ztim1=SIN(declin)
-               ztim2=COS(declin)*COS(2.*pi*(zday-.5))
-               ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
-
-               call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
-               ztim1,ztim2,ztim3,mu0,fract)
-
-           else
-
-               call mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
-               ! WARNING: this function appears not to work in 1D
-
-           endif
-
-           if (corrk) then
-
-!          a) Call correlated-k radiative transfer scheme
-!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-              if(kastprof)then
-                 print*,'kastprof should not = true here'
-                 call abort
-              endif
-              muvar(:,:)=0.0 ! only used for climate evolution studies (kcm1d) for now
-      
-!             standard callcorrk
-              clearsky=.false.
-              call callcorrk(ngrid,nlayer,pq,nq,qsurf,                   &
-                      albedo,emis,mu0,pplev,pplay,pt,                    &
-                      tsurf,fract,dist_star,aerosol,muvar,               &
-                      zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
-                      fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,               &
-                      reffrad,tau_col,cloudfrac,totcloudfrac,            &
-                      clearsky,firstcall,lastcall)      
-
-!             Option to call scheme once more for clear regions
-              if(CLFvarying)then
-
-                 ! ---> PROBLEMS WITH ALLOCATED ARRAYS 
-                 ! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
-                 clearsky=.true.
-                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,                  &
-                      albedo,emis,mu0,pplev,pplay,pt,                      &
-                      tsurf,fract,dist_star,aerosol,muvar,                 &
-                      zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
-                      fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,              &
-                      reffrad1,tau_col1,cloudfrac,totcloudfrac,             &
-                      clearsky,firstcall,lastcall)
-                 clearsky = .false.  ! just in case.     
-
-                 ! Sum the fluxes and heating rates from cloudy/clear cases
-                 do ig=1,ngrid
-                    tf=totcloudfrac(ig)
-                    ntf=1.-tf
-                    
-                    fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw(ig)
-                    fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw(ig)
-                    fluxtop_lw(ig)  = ntf*fluxtop_lw1(ig)  + tf*fluxtop_lw(ig)
-                    fluxabs_sw(ig)  = ntf*fluxabs_sw1(ig)  + tf*fluxabs_sw(ig)
-                    tau_col(ig)     = ntf*tau_col1(ig)     + tf*tau_col(ig)
-                    
-                    zdtlw(ig,1:nlayermx) = ntf*zdtlw1(ig,1:nlayermx) + tf*zdtlw(ig,1:nlayermx)
-                    zdtsw(ig,1:nlayermx) = ntf*zdtsw1(ig,1:nlayermx) + tf*zdtsw(ig,1:nlayermx)
-
-                    OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV)		       
-                    OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI)		       
-
-                 enddo
-
-              endif !CLFvarying
-
-!             Radiative flux from the sky absorbed by the surface (W.m-2)
-              GSR=0.0
-              fluxrad_sky(1:ngridmx)=emis(1:ngridmx)*fluxsurf_lw(1:ngridmx)+fluxsurf_sw(1:ngridmx)*(1.-albedo(1:ngridmx))
-
-              if(noradsurf)then ! no lower surface; SW flux just disappears
-                 GSR = SUM(fluxsurf_sw(1:ngridmx)*area(1:ngridmx))/totarea
-                 fluxrad_sky(1:ngridmx)=emis(1:ngridmx)*fluxsurf_lw(1:ngridmx)
-                 print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
-              endif
-
-!             Net atmospheric radiative heating rate (K.s-1)
-              dtrad(1:ngridmx,1:nlayermx)=zdtsw(1:ngridmx,1:nlayermx)+zdtlw(1:ngridmx,1:nlayermx)
-
-           elseif(newtonian)then
-
-!          b) Call Newtonian cooling scheme
-!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-              call newtrelax(mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
-
-              zdtsurf(1:ngridmx) = +(pt(1:ngridmx,1)-tsurf(1:ngridmx))/ptimestep
-              ! e.g. surface becomes proxy for 1st atmospheric layer ?
-
-           else
-
-!          c) Atmosphere has no radiative effect
-!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-              fluxtop_dn(1:ngridmx)  = fract(1:ngridmx)*mu0(1:ngridmx)*Fat1AU/dist_star**2
-              if(ngrid.eq.1)then ! / by 4 globally in 1D case...
-                 fluxtop_dn(1)  = fract(1)*Fat1AU/dist_star**2/2.0
-              endif
-              fluxsurf_sw(1:ngridmx) = fluxtop_dn(1:ngridmx)
-              fluxrad_sky(1:ngridmx) = fluxtop_dn(1:ngridmx)*(1.-albedo(1:ngridmx))
-              fluxtop_lw(1:ngridmx)  = emis(1:ngridmx)*sigma*tsurf(1:ngridmx)**4
-              ! radiation skips the atmosphere entirely
-
-
-              dtrad(1:ngridmx,1:nlayermx)=0.0
-              ! hence no atmospheric radiative heating
-
-           endif                ! if corrk
-
-        endif ! of if(mod(icount-1,iradia).eq.0)
-
-
-!    Transformation of the radiative tendencies
-!    ------------------------------------------
-
-        zplanck(1:ngridmx)=tsurf(1:ngridmx)*tsurf(1:ngridmx)
-        zplanck(1:ngridmx)=emis(1:ngridmx)*sigma*zplanck(1:ngridmx)*zplanck(1:ngridmx)
-        fluxrad(1:ngridmx)=fluxrad_sky(1:ngridmx)-zplanck(1:ngridmx)
-        pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+dtrad(1:ngridmx,1:nlayermx)
-
-!-------------------------
-! test energy conservation
-         if(enertest)then
-            dEtotSW  = cpp*SUM(massarea(:,:)*zdtsw(:,:))/totarea
-            dEtotLW  = cpp*SUM(massarea(:,:)*zdtlw(:,:))/totarea
-            dEtotsSW = SUM(fluxsurf_sw(:)*(1.-albedo(:))*area(:))/totarea
-            dEtotsLW = SUM((fluxsurf_lw(:)*emis(:)-zplanck(:))*area(:))/totarea
-	    dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:)
-	    dEzRadlw(:,:)=cpp*mass(:,:)*zdtlw(:,:)
-            print*,'---------------------------------------------------------------'
-            print*,'In corrk SW atmospheric heating       =',dEtotSW,' W m-2'
-            print*,'In corrk LW atmospheric heating       =',dEtotLW,' W m-2'
-            print*,'atmospheric net rad heating (SW+LW)   =',dEtotLW+dEtotSW,' W m-2'
-            print*,'In corrk SW surface heating           =',dEtotsSW,' W m-2'
-            print*,'In corrk LW surface heating           =',dEtotsLW,' W m-2'
-            print*,'surface net rad heating (SW+LW)       =',dEtotsLW+dEtotsSW,' W m-2'
-         endif
-!-------------------------
-
-      endif ! of if (callrad)
-
-!-----------------------------------------------------------------------
-!    4. Vertical diffusion (turbulent mixing):
-!    -----------------------------------------
-
-      if (calldifv) then
-      
-         zflubid(1:ngridmx)=fluxrad(1:ngridmx)+fluxgrd(1:ngridmx)
-
-         zdum1(1:ngridmx,1:nlayermx)=0.0
-         zdum2(1:ngridmx,1:nlayermx)=0.0
-
-
-!JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff	 
-	 if (UseTurbDiff) then
-	 
-           call turbdiff(ngrid,nlayer,nq,rnat,       &
-              ptimestep,capcal,lwrite,               &
-              pplay,pplev,zzlay,zzlev,z0,            &
-              pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf,   &
-              zdum1,zdum2,pdt,pdq,zflubid,           &
-              zdudif,zdvdif,zdtdif,zdtsdif,          &
-	      sensibFlux,q2,zdqdif,zdqevap,zdqsdif,lastcall)
-
-	 else
-	 
-           zdh(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)/zpopsk(1:ngridmx,1:nlayermx)
- 
-           call vdifc(ngrid,nlayer,nq,rnat,zpopsk,   &
-              ptimestep,capcal,lwrite,               &
-              pplay,pplev,zzlay,zzlev,z0,            &
-              pu,pv,zh,pq,tsurf,emis,qsurf,          &
-              zdum1,zdum2,zdh,pdq,zflubid,           &
-              zdudif,zdvdif,zdhdif,zdtsdif,          &
-	      sensibFlux,q2,zdqdif,zdqsdif,lastcall)
-
-           zdtdif(1:ngridmx,1:nlayermx)=zdhdif(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx) ! for diagnostic only
-	   zdqevap(1:ngridmx,1:nlayermx)=0.
-
-         end if !turbdiff
-
-         pdv(1:ngridmx,1:nlayermx)=pdv(1:ngridmx,1:nlayermx)+zdvdif(1:ngridmx,1:nlayermx)
-         pdu(1:ngridmx,1:nlayermx)=pdu(1:ngridmx,1:nlayermx)+zdudif(1:ngridmx,1:nlayermx)
-         pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+zdtdif(1:ngridmx,1:nlayermx)
-         zdtsurf(1:ngridmx)=zdtsurf(1:ngridmx)+zdtsdif(1:ngridmx)
-         if (tracer) then 
-           pdq(1:ngridmx,1:nlayermx,1:nqmx)=pdq(1:ngridmx,1:nlayermx,1:nqmx)+ zdqdif(1:ngridmx,1:nlayermx,1:nqmx) 
-           dqsurf(1:ngridmx,1:nqmx)=dqsurf(1:ngridmx,1:nqmx) + zdqsdif(1:ngridmx,1:nqmx)
-         end if ! of if (tracer)
-
-         !-------------------------
-         ! test energy conservation
-         if(enertest)then
-	    dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:)
-            do ig = 1, ngrid
-	       dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
-	       dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
-            enddo
-	    dEtot = SUM(dEdiff(:)*area(:))/totarea
-	    dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
-            dEtots = SUM(dEdiffs(:)*area(:))/totarea
-	    AtmToSurf_TurbFlux=SUM(sensibFlux(:)*area(:))/totarea
-            if (UseTurbDiff) then
-	       print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
-               print*,'In TurbDiff non-cons atm nrj change   =',dEtot,' W m-2'
-               print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
-	    else
-	       print*,'In vdifc sensible flux (atm=>surf)    =',AtmToSurf_TurbFlux,' W m-2'
-               print*,'In vdifc non-cons atm nrj change      =',dEtot,' W m-2'
-               print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
-	    end if
-! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme
-!    but not given back elsewhere
-         endif
-         !-------------------------
-
-         !-------------------------
-         ! test water conservation
-         if(watertest.and.water)then
-            dWtot = SUM(massarea(:,:)*zdqdif(:,:,igcm_h2o_vap))*ptimestep/totarea
-            dWtots =  SUM(zdqsdif(:,igcm_h2o_vap)*area(:))*ptimestep/totarea
-            do ig = 1, ngrid
-	       vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap))
-	    Enddo
-            dWtot = dWtot + SUM(massarea(:,:)*zdqdif(:,:,igcm_h2o_ice))*ptimestep/totarea
-            dWtots = dWtots + SUM(zdqsdif(:,igcm_h2o_ice)*area(:))*ptimestep/totarea
-            do ig = 1, ngrid
-	       vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice))
-	    Enddo	    
-            nconsMAX=MAXVAL(vdifcncons(:))
-
-            print*,'---------------------------------------------------------------'
-            print*,'In difv atmospheric water change        =',dWtot,' kg m-2'
-            print*,'In difv surface water change            =',dWtots,' kg m-2'
-            print*,'In difv non-cons factor                 =',dWtot+dWtots,' kg m-2'
-            print*,'In difv MAX non-cons factor             =',nconsMAX,' kg m-2 s-1'
-
-         endif
-         !-------------------------
-
-      else    
-
-         if(.not.newtonian)then
-
-            zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + (fluxrad(1:ngridmx) + fluxgrd(1:ngridmx))/capcal(1:ngridmx)
-
-         endif
-
-      endif ! of if (calldifv)
-
-
-!-----------------------------------------------------------------------
-!   5. Dry convective adjustment:
-!   -----------------------------
-
-      if(calladj) then
-
-         zdh(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)/zpopsk(1:ngridmx,1:nlayermx)
-         zduadj(1:ngridmx,1:nlayermx)=0.0
-         zdvadj(1:ngridmx,1:nlayermx)=0.0
-         zdhadj(1:ngridmx,1:nlayermx)=0.0
-
-
-         call convadj(ngrid,nlayer,nq,ptimestep,    &
-              pplay,pplev,zpopsk,                   &
-              pu,pv,zh,pq,                          &
-              pdu,pdv,zdh,pdq,                      &
-              zduadj,zdvadj,zdhadj,                 &
-              zdqadj)
-
-         pdu(1:ngridmx,1:nlayermx) = pdu(1:ngridmx,1:nlayermx) + zduadj(1:ngridmx,1:nlayermx)
-         pdv(1:ngridmx,1:nlayermx) = pdv(1:ngridmx,1:nlayermx) + zdvadj(1:ngridmx,1:nlayermx)
-         pdt(1:ngridmx,1:nlayermx)    = pdt(1:ngridmx,1:nlayermx) + zdhadj(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx)
-         zdtadj(1:ngridmx,1:nlayermx) = zdhadj(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx) ! for diagnostic only
- 
-         if(tracer) then 
-            pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqadj(1:ngridmx,1:nlayermx,1:nqmx) 
-         end if
-
-         !-------------------------
-         ! test energy conservation
-         if(enertest)then
-            dEtot=cpp*SUM(massarea(:,:)*zdtadj(:,:))/totarea
-          print*,'In convadj atmospheric energy change  =',dEtot,' W m-2'
-         endif
-         !-------------------------
-
-         !-------------------------
-         ! test water conservation
-         if(watertest)then
-            dWtot = SUM(massarea(:,:)*zdqadj(:,:,igcm_h2o_vap))*ptimestep/totarea
-            do ig = 1, ngrid
-	       cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap))
-	    Enddo
-            dWtot = dWtot + SUM(massarea(:,:)*zdqadj(:,:,igcm_h2o_ice))*ptimestep/totarea
-            do ig = 1, ngrid
-	       cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice))
-	    Enddo	    
-            nconsMAX=MAXVAL(cadjncons(:))
-
-            print*,'In convadj atmospheric water change     =',dWtot,' kg m-2'
-            print*,'In convadj MAX non-cons factor          =',nconsMAX,' kg m-2 s-1'
-         endif
-         !-------------------------
-               
-      endif ! of if(calladj)
-
-!-----------------------------------------------------------------------
-!   6. Carbon dioxide condensation-sublimation:
-!   -------------------------------------------
-
-      if (co2cond) then
-         if (.not.tracer) then
-            print*,'We need a CO2 ice tracer to condense CO2'
-            call abort
-         endif
-     print*, 'we are in co2cond!!!!!'
-         call condense_cloud(ngrid,nlayer,nq,ptimestep,   &
-              capcal,pplay,pplev,tsurf,pt,                  &
-              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,        &
-              qsurf(1,igcm_co2_ice),albedo,emis,            &
-              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,               &
-              zdqc,reffrad)
-
-
-         pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+zdtc(1:ngridmx,1:nlayermx)
-         pdv(1:ngridmx,1:nlayermx)=pdv(1:ngridmx,1:nlayermx)+zdvc(1:ngridmx,1:nlayermx)
-         pdu(1:ngridmx,1:nlayermx)=pdu(1:ngridmx,1:nlayermx)+zduc(1:ngridmx,1:nlayermx)
-         zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + zdtsurfc(1:ngridmx)
-
-         pdq(1:ngridmx,1:nlayermx,1:nqmx)=pdq(1:ngridmx,1:nlayermx,1:nqmx)+ zdqc(1:ngridmx,1:nlayermx,1:nqmx) 
-         ! Note: we do not add surface co2ice tendency
-         ! because qsurf(:,igcm_co2_ice) is updated in condens_co2cloud
-
-         !-------------------------
-         ! test energy conservation
-         if(enertest)then
-            dEtot = cpp*SUM(massarea(:,:)*zdtc(:,:))/totarea
-            dEtots = SUM(capcal(:)*zdtsurfc(:)*area(:))/totarea
-            print*,'In co2cloud atmospheric energy change   =',dEtot,' W m-2'
-            print*,'In co2cloud surface energy change       =',dEtots,' W m-2'
-         endif
-         !-------------------------
-
-      endif  ! of if (co2cond)
-
-
-!-----------------------------------------------------------------------
-!   7. Specific parameterizations for tracers 
-!   -----------------------------------------
-
-      if (tracer) then 
-
-!   7a. Water and ice
-!     ---------------
-         if (water) then
-
-!        ----------------------------------------
-!        Water ice condensation in the atmosphere
-!        ----------------------------------------
-            if(watercond.and.(RLVTT.gt.1.e-8))then
-
-!             ----------------
-!             Moist convection
-!             ----------------
-               ! Re-evaporate cloud water/ice
-!	       dqevap1(1:ngridmx,1:nlayermx)=0.
-!	       dtevap1(1:ngridmx,1:nlayermx)=0.
-!               call evap(ptimestep,pt,pq,pdq,pdt,dqevap1,dtevap1,qevap,tevap)
-               
-	       dqmoist(1:ngridmx,1:nlayermx,1:nqmx)=0.
-	       dtmoist(1:ngridmx,1:nlayermx)=0.
-
-               call moistadj(pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
-
-               pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap)   &
-	                  +dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_vap)
-               pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) =pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice)     &
-                          +dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_ice)
-               pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+dtmoist(1:ngridmx,1:nlayermx)
-
-!               do l=1,nlayermx
-!	         do ig=1,ngridmx
-!	          if(isnan(dqmoist(ig,l,igcm_h2o_ice))) print*,'Nan in dqmoist ice, abort,ig,l',ig,l
-!	          if(isnan(dqmoist(ig,l,igcm_h2o_vap))) print*,'Nan in dqmoist vap, abort,ig,l',ig,l
-!	          if(isnan(dtmoist(ig,l))) print*,'Nan in dtmoist, abort,ig,l',ig,l		  
-!                  if(isnan(dqmoist(ig,l,igcm_h2o_ice)).or.isnan(dqmoist(ig,l,igcm_h2o_vap)).or.isnan(dtmoist(ig,l))) STOP
-!	         enddo
-!	       enddo
-!	       print*,'after moist',dqmoist(185,1,igcm_h2o_ice),dqmoist(185,1,igcm_h2o_vap),dtmoist(185,1)*86400.
-
-               !-------------------------
-               ! test energy conservation
-               if(enertest)then
-                  dEtot=cpp*SUM(massarea(:,:)*dtmoist(:,:))/totarea
-                  do ig = 1, ngrid
-                     madjdE(ig) = cpp*SUM(mass(:,:)*dtmoist(:,:))
-                  enddo
-                  print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
-                  print*,'In moistadj MAX atmospheric energy change   =',MAXVAL(pdt(:,:))*86400.,'K/day'
-		  
-		! test energy conservation
-	          dWtot = SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap))*ptimestep/totarea
-                  dWtot = dWtot + SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_ice))*ptimestep/totarea
-                  print*,'In moistadj atmospheric water change    =',dWtot,' kg m-2'
-               endif
-               !-------------------------
-	       
-
-!        --------------------------------
-!        Large scale condensation/evaporation
-!        --------------------------------
-
-               call largescale(ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc)
-
-               pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+dtlscale(1:ngridmx,1:nlayermx)
-               pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap)+dqvaplscale(1:ngridmx,1:nlayermx)
-               pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice)+dqcldlscale(1:ngridmx,1:nlayermx)
-
-               !-------------------------
-               ! test energy conservation
-               if(enertest)then
-                  do ig = 1, ngrid
-                     lscaledE(ig) = cpp*SUM(mass(:,:)*dtlscale(:,:))
-                  enddo
-                  dEtot=cpp*SUM(massarea(:,:)*(dtlscale(:,:)))/totarea
-		  if(isnan(dEtot)) then
-		     print*,'Nan in largescale, abort'
-                     STOP
-		  endif
-		  print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
-
-               ! test water conservation
-	          dWtot = SUM(massarea(:,:)*dqvaplscale(:,:))*ptimestep/totarea
-                  dWtot = dWtot + SUM(massarea(:,:)*dqcldlscale(:,:))*ptimestep/totarea
-                  print*,'In largescale atmospheric water change  =',dWtot,' kg m-2'
-               endif
-               !-------------------------
-
-               ! compute cloud fraction
-               do l = 1, nlayer
-                  do ig = 1,ngrid
-                     cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l)) 
-                  enddo
-               enddo
-
-            endif                ! of if (watercondense)
-           
-
-!        --------------------------------
-!        Water ice / liquid precipitation
-!        --------------------------------
-            if(waterrain)then
-
-               zdqrain(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
-               zdqsrain(1:ngridmx)    = 0.0
-               zdqssnow(1:ngridmx)    = 0.0
-
-               call rain(ptimestep,pplev,pplay,pt,pdt,pq,pdq,            &
-                   zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
-
-               pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) &
-	             +zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_vap)
-               pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) &
-	             +zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_ice)
-               pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+zdtrain(1:ngridmx,1:nlayermx)
-               dqsurf(1:ngridmx,igcm_h2o_vap) = dqsurf(1:ngridmx,igcm_h2o_vap)+zdqsrain(1:ngridmx) ! a bug was here
-               dqsurf(1:ngridmx,igcm_h2o_ice) = dqsurf(1:ngridmx,igcm_h2o_ice)+zdqssnow(1:ngridmx)
-               rainout(1:ngridmx)             = zdqsrain(1:ngridmx)+zdqssnow(1:ngridmx) ! diagnostic   
-
-
-               !-------------------------
-               ! test energy conservation
-               if(enertest)then
-                  dEtot=cpp*SUM(massarea(:,:)*zdtrain(:,:))/totarea
-                 print*,'In rain atmospheric T energy change       =',dEtot,' W m-2'
-                  dItot = SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_ice))/totarea*RLVTT/cpp
-		  dItot = dItot + SUM(area(:)*zdqssnow(:))/totarea*RLVTT/cpp
-                  dVtot = SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap))*ptimestep/totarea
-		  dVtot = dVtot + SUM(area(:)*zdqsrain(:))/totarea*RLVTT/cpp
-		  dEtot = dItot + dVtot
-                 print*,'In rain dItot =',dItot,' W m-2'
-                 print*,'In rain dVtot =',dVtot,' W m-2'
-                 print*,'In rain atmospheric L energy change       =',dEtot,' W m-2'
-
-               ! test water conservation
-	          dWtot = SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap))*ptimestep/totarea
-                  dWtot = dWtot + SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_ice))*ptimestep/totarea
-		  dWtots =  SUM((zdqsrain(:)+zdqssnow(:))*area(:))*ptimestep/totarea
-                 print*,'In rain atmospheric water change        =',dWtot,' kg m-2'
-                 print*,'In rain surface water change            =',dWtots,' kg m-2'
-                 print*,'In rain non-cons factor                 =',dWtot+dWtots,' kg m-2'
-              endif
-              !-------------------------
-
-            end if                 ! of if (waterrain)
-         end if                    ! of if (water)
-
-
-!   7c. Aerosol particles
-!     -------------------
-!        ------------- 
-!        Sedimentation
-!        ------------- 
-        if (sedimentation) then 
-           zdqsed(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
-           zdqssed(1:ngridmx,1:nqmx)  = 0.0
-
-
-           !-------------------------
-           ! find qtot
-           if(watertest)then
-              iq=3
-	      dWtot = SUM(massarea(:,:)*pq(:,:,iq))*ptimestep/totarea
-	      dWtots = SUM(massarea(:,:)*pdq(:,:,iq))*ptimestep/totarea
-              print*,'Before sedim pq  =',dWtot,' kg m-2'
-              print*,'Before sedim pdq =',dWtots,' kg m-2'
-           endif
-           !-------------------------
-
-           call callsedim(ngrid,nlayer,ptimestep,           &
-                pplev,zzlev,pt,pq,pdq,zdqsed,zdqssed,nq,reffH2O)
-
-           !-------------------------
-           ! find qtot
-           if(watertest)then
-              iq=3
-	      dWtot = SUM(massarea(:,:)*pq(:,:,iq))*ptimestep/totarea
-	      dWtots = SUM(massarea(:,:)*pdq(:,:,iq))*ptimestep/totarea
-              print*,'After sedim pq  =',dWtot,' kg m-2'
-              print*,'After sedim pdq =',dWtots,' kg m-2'
-           endif
-           !-------------------------
-
-           ! for now, we only allow H2O ice to sediment
-              ! and as in rain.F90, whether it falls as rain or snow depends
-              ! only on the surface temperature
-           pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqsed(1:ngridmx,1:nlayermx,1:nqmx)
-           dqsurf(1:ngridmx,1:nqmx) = dqsurf(1:ngridmx,1:nqmx) + zdqssed(1:ngridmx,1:nqmx)
-
-           !-------------------------
-           ! test water conservation
-           if(watertest)then
-              dWtot = SUM(massarea(:,:)*(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice)))*ptimestep/totarea
-	      dWtots = SUM((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))*area(:))*ptimestep/totarea
-              print*,'In sedim atmospheric ice change         =',dWtot,' kg m-2'
-              print*,'In sedim surface ice change             =',dWtots,' kg m-2'
-              print*,'In sedim non-cons factor                =',dWtot+dWtots,' kg m-2'
-           endif
-           !-------------------------
-
-        end if   ! of if (sedimentation)
-
-
-!   7d. Updates
-!     ---------
-
-!       -----------------------------------
-!       Updating atm mass and tracer budget
-!       -----------------------------------
-
-         if(mass_redistrib) then
-
-            zdmassmr(1:ngridmx,1:nlayermx) = mass(1:ngridmx,1:nlayermx) * &
-	       (   zdqevap(1:ngridmx,1:nlayermx)                          &
-	         + zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_vap)             &
-		 + dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_vap)             &
-		 + dqvaplscale(1:ngridmx,1:nlayermx) )
-		 
-            
-	    call writediagfi(ngridmx,"mass_evap","mass gain"," ",3,zdmassmr)
-	    call writediagfi(ngridmx,"mass","mass"," ",3,mass)
-
-	    call mass_redistribution(ngrid,nlayer,nq,ptimestep,   &
-                       rnat,capcal,pplay,pplev,pt,tsurf,pq,qsurf,     &
-		       pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr,  &
-		       zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
-	 
-
-            pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqmr(1:ngridmx,1:nlayermx,1:nqmx)
-            dqsurf(1:ngridmx,1:nqmx)         = dqsurf(1:ngridmx,1:nqmx) + zdqsurfmr(1:ngridmx,1:nqmx)
-            pdt(1:ngridmx,1:nlayermx)        = pdt(1:ngridmx,1:nlayermx) + zdtmr(1:ngridmx,1:nlayermx)
-            pdu(1:ngridmx,1:nlayermx)        = pdu(1:ngridmx,1:nlayermx) + zdumr(1:ngridmx,1:nlayermx)
-            pdv(1:ngridmx,1:nlayermx)        = pdv(1:ngridmx,1:nlayermx) + zdvmr(1:ngridmx,1:nlayermx)
-	    pdpsrf(1:ngridmx)                = pdpsrf(1:ngridmx) + zdpsrfmr(1:ngridmx)
-            zdtsurf(1:ngridmx)               = zdtsurf(1:ngridmx) + zdtsurfmr(1:ngridmx)
-	    
-!	    print*,'after mass redistrib, q=',pq(211,1:nlayermx,igcm_h2o_vap)+ptimestep*pdq(211,1:nlayermx,igcm_h2o_vap) 
-	 endif
-
-
-
-!       ---------------------------------
-!       Updating tracer budget on surface
-!       ---------------------------------
-
-         if(hydrology)then
-
-            call hydrol(ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd,    &
-               capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice)
-            ! note: for now, also changes albedo in the subroutine
-
-            zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + zdtsurf_hyd(1:ngridmx)
-            qsurf(1:ngridmx,1:nqmx) = qsurf(1:ngridmx,1:nqmx) + ptimestep*dqs_hyd(1:ngridmx,1:nqmx)
-            ! when hydrology is used, other dqsurf tendencies are all added to dqs_hyd inside
-
-            !-------------------------
-            ! test energy conservation
-            if(enertest)then
-               dEtots = SUM(area(:)*capcal(:)*zdtsurf_hyd(:))/totarea
-               print*,'In hydrol surface energy change     =',dEtots,' W m-2'
-            endif
-            !-------------------------
-
-            !-------------------------
-            ! test water conservation
-            if(watertest)then
-	       dWtots =  SUM(dqs_hyd(:,igcm_h2o_ice)*area(:))*ptimestep/totarea
-               print*,'In hydrol surface ice change            =',dWtots,' kg m-2'
-	       dWtots =  SUM(dqs_hyd(:,igcm_h2o_vap)*area(:))*ptimestep/totarea
-               print*,'In hydrol surface water change          =',dWtots,' kg m-2'
-               print*,'---------------------------------------------------------------'
-            endif
-            !-------------------------
-
-         ELSE     ! of if (hydrology)
-
-            qsurf(1:ngridmx,1:nqmx)=qsurf(1:ngridmx,1:nqmx)+ptimestep*dqsurf(1:ngridmx,1:nqmx)
-
-         END IF   ! of if (hydrology)
-
-         ! Add qsurf to qsurf_hist, which is what we save in
-         ! diagfi.nc etc. At the same time, we set the water 
-         ! content of ocean gridpoints back to zero, in order
-         ! to avoid rounding errors in vdifc, rain
-         qsurf_hist(:,:) = qsurf(:,:)
-
-         if(ice_update)then
-            ice_min(1:ngridmx)=min(ice_min(1:ngridmx),qsurf(1:ngridmx,igcm_h2o_ice))
-         endif
-
-      endif   !  of if (tracer) 
-
-!-----------------------------------------------------------------------
-!   9. Surface and sub-surface soil temperature
-!-----------------------------------------------------------------------
-
-
-!     Increment surface temperature
-      tsurf(1:ngridmx)=tsurf(1:ngridmx)+ptimestep*zdtsurf(1:ngridmx) 
-
-!     Compute soil temperatures and subsurface heat flux
-      if (callsoil) then
-         call soil(ngrid,nsoilmx,.false.,inertiedat,        &
-                ptimestep,tsurf,tsoil,capcal,fluxgrd)
-      endif
-
-!-------------------------
-! test energy conservation
-      if(enertest)then
-         dEtots = SUM(area(:)*capcal(:)*zdtsurf(:))/totarea
-         print*,'Surface energy change                 =',dEtots,' W m-2'
-      endif
-!-------------------------
-
-!-----------------------------------------------------------------------
-!  10. Perform diagnostics and write output files
-!-----------------------------------------------------------------------
-
-!    -------------------------------
-!    Dynamical fields incrementation
-!    -------------------------------
-!    For output only: the actual model integration is performed in the dynamics
- 
-      ! temperature, zonal and meridional wind
-      zt(1:ngridmx,1:nlayermx) = pt(1:ngridmx,1:nlayermx) + pdt(1:ngridmx,1:nlayermx)*ptimestep
-      zu(1:ngridmx,1:nlayermx) = pu(1:ngridmx,1:nlayermx) + pdu(1:ngridmx,1:nlayermx)*ptimestep
-      zv(1:ngridmx,1:nlayermx) = pv(1:ngridmx,1:nlayermx) + pdv(1:ngridmx,1:nlayermx)*ptimestep
-
-      ! diagnostic
-      zdtdyn(1:ngridmx,1:nlayermx)     = pt(1:ngridmx,1:nlayermx)-ztprevious(1:ngridmx,1:nlayermx)
-      ztprevious(1:ngridmx,1:nlayermx) = zt(1:ngridmx,1:nlayermx)
-
-      if(firstcall)then
-         zdtdyn(1:ngridmx,1:nlayermx)=0.0
-      endif
-
-      ! dynamical heating diagnostic
-      do ig=1,ngrid
-         fluxdyn(ig)= SUM(zdtdyn(ig,:) *mass(ig,:))*cpp/ptimestep
-      enddo
-
-      ! tracers
-      zq(1:ngridmx,1:nlayermx,1:nqmx) = pq(1:ngridmx,1:nlayermx,1:nqmx) + pdq(1:ngridmx,1:nlayermx,1:nqmx)*ptimestep
-
-      ! surface pressure
-      ps(1:ngridmx) = pplev(1:ngridmx,1) + pdpsrf(1:ngridmx)*ptimestep
-
-      ! pressure
-      do l=1,nlayer
-          zplev(1:ngridmx,l) = pplev(1:ngridmx,l)/pplev(1:ngridmx,1)*ps(:)
-          zplay(1:ngridmx,l) = pplay(1:ngridmx,l)/pplev(1:ngridmx,1)*ps(1:ngridmx)
-      enddo
-
-!     ---------------------------------------------------------
-!     Surface and soil temperature information
-!     ---------------------------------------------------------
-
-      Ts1 = SUM(area(:)*tsurf(:))/totarea
-      Ts2 = MINVAL(tsurf(:))
-      Ts3 = MAXVAL(tsurf(:))
-      if(callsoil)then
-         TsS = SUM(area(:)*tsoil(:,nsoilmx))/totarea        ! mean temperature at bottom soil layer
-         print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]     ave[Tdeep]'
-         print*,Ts1,Ts2,Ts3,TsS
-      else
-         print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]'
-         print*,Ts1,Ts2,Ts3
-      endif
-
-!     ---------------------------------------------------------
-!     Check the energy balance of the simulation during the run
-!     ---------------------------------------------------------
-
-      if(corrk)then
-
-         ISR = SUM(area(:)*fluxtop_dn(:))/totarea
-         ASR = SUM(area(:)*fluxabs_sw(:))/totarea
-         OLR = SUM(area(:)*fluxtop_lw(:))/totarea
-         GND = SUM(area(:)*fluxgrd(:))/totarea
-         DYN = SUM(area(:)*fluxdyn(:))/totarea
-         do ig=1,ngrid             
-            if(fluxtop_dn(ig).lt.0.0)then
-               print*,'fluxtop_dn has gone crazy'
-               print*,'fluxtop_dn=',fluxtop_dn(ig)
-               print*,'tau_col=',tau_col(ig)
-               print*,'aerosol=',aerosol(ig,:,:)
-               print*,'temp=   ',pt(ig,:)
-               print*,'pplay=  ',pplay(ig,:)
-               call abort
-            endif
-         end do
-                     
-         if(ngridmx.eq.1)then
-            DYN=0.0
-         endif
-
-         print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
-         print*, ISR,ASR,OLR,GND,DYN
-         
-         if(enertest)then
-            print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
-            print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
-            print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
-         endif
-
-         if(meanOLR)then
-            if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
-               ! to record global radiative balance
-               open(92,file="rad_bal.out",form='formatted',position='append')
-               write(92,*) zday,ISR,ASR,OLR
-               close(92)
-               open(93,file="tem_bal.out",form='formatted',position='append')
-               write(93,*) zday,Ts1,Ts2,Ts3,TsS
-               close(93)
-            endif
-         endif
-
-      endif
-
-
-!     ------------------------------------------------------------------
-!     Diagnostic to test radiative-convective timescales in code
-!     ------------------------------------------------------------------
-      if(testradtimes)then
-         open(38,file="tau_phys.out",form='formatted',position='append')
-         ig=1
-         do l=1,nlayer
-            write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
-         enddo
-         close(38)
-         print*,'As testradtimes enabled,'
-         print*,'exiting physics on first call'
-         call abort
-      endif
-
-!     ---------------------------------------------------------
-!     Compute column amounts (kg m-2) if tracers are enabled
-!     ---------------------------------------------------------
-      if(tracer)then
-         qcol(1:ngridmx,1:nqmx)=0.0
-         do iq=1,nq
-           do ig=1,ngridmx
-              qcol(ig,iq) = SUM( zq(ig,1:nlayermx,iq) * mass(ig,1:nlayermx))
-           enddo
-         enddo
-
-         ! Generalised for arbitrary aerosols now. (LK)
-         reffcol(1:ngridmx,1:naerkind)=0.0
-         if(co2cond.and.(iaero_co2.ne.0))then
-            call co2_reffrad(zq,reffrad)
-            do ig=1,ngridmx
-               reffcol(ig,iaero_co2) = SUM(zq(ig,1:nlayermx,igcm_co2_ice)*reffrad(ig,1:nlayermx,iaero_co2)*mass(ig,1:nlayermx))
-            enddo
-         endif
-         if(water.and.(iaero_h2o.ne.0))then
-            call h2o_reffrad(zq,zt,reffrad,nueffrad)
-            do ig=1,ngridmx
-               reffcol(ig,iaero_h2o) = SUM(zq(ig,1:nlayermx,igcm_h2o_ice)*reffrad(ig,1:nlayermx,iaero_h2o)*mass(ig,1:nlayermx))
-            enddo
-         endif
-
-      endif
-
-!     ---------------------------------------------------------
-!     Test for water conservation if water is enabled
-!     ---------------------------------------------------------
-
-      if(water)then
-
-         icesrf = SUM(area(:)*qsurf_hist(:,igcm_h2o_ice))/totarea
-         liqsrf = SUM(area(:)*qsurf_hist(:,igcm_h2o_vap))/totarea
-         icecol = SUM(area(:)*qcol(:,igcm_h2o_ice))/totarea
-         vapcol = SUM(area(:)*qcol(:,igcm_h2o_vap))/totarea
-
-         h2otot = icesrf + liqsrf + icecol + vapcol
-
-         print*,' Total water amount [kg m^-2]: ',h2otot
-         print*,' Surface ice    Surface liq.   Atmos. con.     Atmos. vap. [kg m^-2] '
-         print*, icesrf,liqsrf,icecol,vapcol
-
-         if(meanOLR)then
-            if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
-               ! to record global water balance
-               open(98,file="h2o_bal.out",form='formatted',position='append')
-               write(98,*) zday,icesrf,liqsrf,icecol,vapcol
-               close(98)
-            endif
-         endif
-
-      endif
-
-!     ---------------------------------------------------------
-!     Calculate RH for diagnostic if water is enabled
-!     ---------------------------------------------------------
-
-      if(water)then
-         do l = 1, nlayer
-            do ig = 1, ngrid
-!               call watersat(pt(ig,l),pplay(ig,l),qsat(ig,l))
-               call Psat_water(zt(ig,l),pplay(ig,l),psat_tmp,qsat(ig,l))
-
-               RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l)
-            enddo
-         enddo
-
-         ! compute maximum possible H2O column amount (100% saturation)
-         do ig=1,ngrid
-               H2Omaxcol(ig) = SUM( qsat(ig,:) * mass(ig,:))
-         enddo
-
-      endif
-
-
-         print*,''
-         print*,'--> Ls =',zls*180./pi
-!        -------------------------------------------------------------------
-!        Writing NetCDF file  "RESTARTFI" at the end of the run
-!        -------------------------------------------------------------------
-!        Note: 'restartfi' is stored just before dynamics are stored
-!              in 'restart'. Between now and the writting of 'restart',
-!              there will have been the itau=itau+1 instruction and
-!              a reset of 'time' (lastacll = .true. when itau+1= itaufin)
-!              thus we store for time=time+dtvr
-
-         if(lastcall) then
-            ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
-
-
-!           Update surface ice distribution to iterate to steady state if requested
-            if(ice_update)then
-
-               do ig = 1, ngrid
-
-                  delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
-
-                  ! add multiple years of evolution
-                  qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep 
-
-                  ! if ice has gone -ve, set to zero
-                  if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
-                     qsurf_hist(ig,igcm_h2o_ice) = 0.0 
-                  endif
-
-                  ! if ice is seasonal, set to zero (NEW)
-                  if(ice_min(ig).lt.0.01)then
-                     qsurf_hist(ig,igcm_h2o_ice) = 0.0 
-                  endif
-
-               enddo
-
-               ! enforce ice conservation
-               ice_tot= SUM(qsurf_hist(:,igcm_h2o_ice)*area(:) )
-               qsurf_hist(:,igcm_h2o_ice) = qsurf_hist(:,igcm_h2o_ice)*(icesrf/ice_tot)
-
-            endif
-
-            write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
-            call physdem1("restartfi.nc",long,lati,nsoilmx,nq,            &
-                    ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
-                    area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,      &
-                    cloudfrac,totcloudfrac,hice)
-         endif
-
-!        -----------------------------------------------------------------
-!        Saving statistics :
-!        -----------------------------------------------------------------
-!        ("stats" stores and accumulates 8 key variables in file "stats.nc"
-!        which can later be used to make the statistic files of the run:
-!        "stats")          only possible in 3D runs !
-
-         
-         if (callstats) then
-
-            call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
-            call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
-            call wstats(ngrid,"fluxsurf_lw",                               &
-                        "Thermal IR radiative flux to surface","W.m-2",2,  &
-                         fluxsurf_lw)
-!            call wstats(ngrid,"fluxsurf_sw",                               &
-!                        "Solar radiative flux to surface","W.m-2",2,       &
-!                         fluxsurf_sw_tot)
-            call wstats(ngrid,"fluxtop_lw",                                &
-                        "Thermal IR radiative flux to space","W.m-2",2,    &
-                        fluxtop_lw)
-!            call wstats(ngrid,"fluxtop_sw",                                &
-!                        "Solar radiative flux to space","W.m-2",2,         &
-!                        fluxtop_sw_tot)
-
-            call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
-            call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
-            call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
-
-            call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
-            call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
-            call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
-            call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
-            call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
-
-           if (tracer) then
-             do iq=1,nq
-                call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
-                call wstats(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
-                          'kg m^-2',2,qsurf(1,iq) )
-
-                call wstats(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    & 
-                          'kg m^-2',2,qcol(1,iq) )
-!                call wstats(ngridmx,trim(noms(iq))//'_reff',                          & 
-!                          trim(noms(iq))//'_reff',                                   & 
-!                          'm',3,reffrad(1,1,iq))
-              end do
-            if (water) then
-               vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
-               call wstats(ngrid,"vmr_h2ovapor",                           &
-                          "H2O vapour volume mixing ratio","mol/mol",       & 
-                          3,vmr)
-            endif ! of if (water)
-
-           endif !tracer 
-
-            if(lastcall) then
-              write (*,*) "Writing stats..."
-              call mkstats(ierr)
-            endif
-          endif !if callstats
-
-
-!        ----------------------------------------------------------------------
-!        output in netcdf file "DIAGFI", containing any variable for diagnostic
-!        (output with  period "ecritphy", set in "run.def")
-!        ----------------------------------------------------------------------
-!        writediagfi can also be called from any other subroutine for any variable.
-!        but its preferable to keep all the calls in one place...
-
-        call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
-        call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
-        call writediagfi(ngrid,"temp","temperature","K",3,zt)
-        call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
-        call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
-        call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
-        call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
-        call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
-
-!     Total energy balance diagnostics
-        if(callrad.and.(.not.newtonian))then
-           call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo)
-           call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
-           call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
-           call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
-           call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
-           call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
-        endif
-	
-	if(enertest) then
-	  if (calldifv) then
-	     call writediagfi(ngridmx,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
-!	     call writediagfi(ngridmx,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
-!	     call writediagfi(ngridmx,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
-!	     call writediagfi(ngridmx,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
-	     call writediagfi(ngridmx,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
-	  endif
-	  if (corrk) then
-!	     call writediagfi(ngridmx,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
-!	     call writediagfi(ngridmx,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
-	  endif 
-          if(watercond) then
-	     call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE) 
-	     call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
-	     call writediagfi(ngridmx,"RH","relative humidity"," ",3,RH)
-!	     call writediagfi(ngridmx,"qsatatm","atm qsat"," ",3,qsat)
-!	     call writediagfi(ngridmx,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
-          endif 
-        endif
-
-!     Temporary inclusions for heating diagnostics
-!        call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
-!        call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
-!        call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
-!        call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
-
-        ! debugging
-        !call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
-        !call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
-
-!     Output aerosols
-        if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
-	  call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
-        if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
-	  call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
-        if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
-	  call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
-        if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
-	  call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
-
-!     Output tracers
-       if (tracer) then
-        do iq=1,nq
-          call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
-!          call writediagfi(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
-!                          'kg m^-2',2,qsurf(1,iq) )
-          call writediagfi(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
-                          'kg m^-2',2,qsurf_hist(1,iq) )
-          call writediagfi(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    & 
-                          'kg m^-2',2,qcol(1,iq) )
-
-          if(watercond.or.CLFvarying)then
-!             call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
-!             call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
-!             call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
-             call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
-          endif
-
-          if(waterrain)then
-             call writediagfi(ngridmx,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
-             call writediagfi(ngridmx,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
-          endif
-
-          if(hydrology)then
-             call writediagfi(ngridmx,"hice","oceanic ice height","m",2,hice)
-          endif
-
-          if(ice_update)then
-             call writediagfi(ngridmx,"ice_min","min annual ice","m",2,ice_min)
-             call writediagfi(ngridmx,"ice_ini","initial annual ice","m",2,ice_initial)
-          endif
-
-          do ig=1,ngrid
-             if(tau_col(ig).gt.1.e3)then
-                print*,'WARNING: tau_col=',tau_col(ig)
-                print*,'at ig=',ig,'in PHYSIQ'
-             endif
-          end do
-
-          call writediagfi(ngridmx,"tau_col","Total aerosol optical depth","[]",2,tau_col)
-
-        enddo
-       endif
-
-!      output spectrum 
-      if(specOLR.and.corrk)then      
-         call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
-         call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
-      endif
-
-
-      icount=icount+1
-
-      ! deallocate gas variables
-      if (lastcall) then
-        IF ( ALLOCATED( gnom ) ) DEALLOCATE( gnom )
-        IF ( ALLOCATED( gfrac ) ) DEALLOCATE( gfrac )  ! both allocated in su_gases.F90
-      endif
-
-
-      return
-    end subroutine physiq
Index: trunk/LMDZ.GENERIC/libf/phystd/physiqg.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiqg.F90	(revision 749)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiqg.F90	(revision 749)
@@ -0,0 +1,1759 @@
+      subroutine physiqg(ngrid,nlayer,nq,   &
+                  firstcall,lastcall,      &
+                  pday,ptime,ptimestep,    &
+                  pplev,pplay,pphi,        &
+                  pu,pv,pt,pq,             &
+                  pw,                      &
+                  pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
+ 
+      use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind
+      use watercommon_h, only : RLVTT, Psat_water
+      use gases_h
+      use radcommon_h, only: sigma
+      use radii_mod, only: h2o_reffrad, co2_reffrad
+      use aerosol_mod
+      implicit none
+
+
+!==================================================================
+!     
+!     Purpose
+!     -------
+!     Central subroutine for all the physics parameterisations in the
+!     universal model. Originally adapted from the Mars LMDZ model.
+!
+!     The model can be run without or with tracer transport
+!     depending on the value of "tracer" in file "callphys.def".
+!
+!
+!   It includes:
+!
+!      1.  Initialization:
+!      1.1 Firstcall initializations
+!      1.2 Initialization for every call to physiq
+!      1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
+!      2. Compute radiative transfer tendencies
+!         (longwave and shortwave).
+!      4. Vertical diffusion (turbulent mixing):
+!      5. Convective adjustment
+!      6. Condensation and sublimation of gases (currently just CO2).
+!      7. TRACERS
+!       7a. water and water ice
+!       7c. other schemes for tracer transport (lifting, sedimentation)
+!       7d. updates (pressure variations, surface budget)
+!      9. Surface and sub-surface temperature calculations
+!     10. Write outputs :
+!           - "startfi", "histfi" if it's time
+!           - Saving statistics if "callstats = .true."
+!           - Output any needed variables in "diagfi" 
+!     10. Diagnostics: mass conservation of tracers, radiative energy balance etc.
+!
+!   arguments
+!   ---------
+!
+!   input
+!   -----
+!    ecri                  period (in dynamical timestep) to write output
+!    ngrid                 Size of the horizontal grid.
+!                          All internal loops are performed on that grid.
+!    nlayer                Number of vertical layers.
+!    nq                    Number of advected fields
+!    firstcall             True at the first call
+!    lastcall              True at the last call
+!    pday                  Number of days counted from the North. Spring
+!                          equinoxe.
+!    ptime                 Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
+!    ptimestep             timestep (s)
+!    pplay(ngrid,nlayer)   Pressure at the middle of the layers (Pa)
+!    pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
+!    pphi(ngrid,nlayer)    Geopotential at the middle of the layers (m2s-2)
+!    pu(ngrid,nlayer)      u component of the wind (ms-1)
+!    pv(ngrid,nlayer)      v component of the wind (ms-1)
+!    pt(ngrid,nlayer)      Temperature (K)
+!    pq(ngrid,nlayer,nq)   Advected fields
+!    pudyn(ngrid,nlayer)    \ 
+!    pvdyn(ngrid,nlayer)     \ Dynamical temporal derivative for the
+!    ptdyn(ngrid,nlayer)     / corresponding variables
+!    pqdyn(ngrid,nlayer,nq) /
+!    pw(ngrid,?)           vertical velocity
+!
+!   output
+!   ------
+!
+!    pdu(ngrid,nlayermx)        \
+!    pdv(ngrid,nlayermx)         \  Temporal derivative of the corresponding
+!    pdt(ngrid,nlayermx)         /  variables due to physical processes.
+!    pdq(ngrid,nlayermx)        /
+!    pdpsrf(ngrid)             /
+!    tracerdyn                 call tracer in dynamical part of GCM ?
+!
+!
+!     Authors
+!     -------
+!           Frederic Hourdin	15/10/93
+!           Francois Forget	1994
+!           Christophe Hourdin	02/1997 
+!           Subroutine completely rewritten by F. Forget (01/2000)
+!           Water ice clouds: Franck Montmessin (update 06/2003)
+!           Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
+!           New correlated-k radiative scheme: R. Wordsworth (2009)
+!           Many specifically Martian subroutines removed: R. Wordsworth (2009)       
+!           Improved water cycle: R. Wordsworth / B. Charnay (2010)
+!           To F90: R. Wordsworth (2010)
+!           New turbulent diffusion scheme: J. Leconte (2012)
+!           Loops converted to F90 matrix format: J. Leconte (2012)
+!     
+!==================================================================
+
+
+!    0.  Declarations :
+!    ------------------
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comgeomfi.h"
+#include "surfdat.h"
+#include "comsoil.h"
+#include "comdiurn.h"
+#include "callkeys.h"
+#include "comcstfi.h"
+#include "planete.h"
+#include "comsaison.h"
+#include "control.h"
+#include "tracer.h"
+#include "watercap.h"
+#include "netcdf.inc"
+
+! Arguments :
+! -----------
+
+!   inputs:
+!   -------
+      integer ngrid,nlayer,nq
+      real ptimestep
+      real pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
+      real pphi(ngridmx,nlayer)
+      real pu(ngridmx,nlayer),pv(ngridmx,nlayer)
+      real pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
+      real pw(ngridmx,nlayer)        ! pvervel transmitted by dyn3d
+      real zh(ngridmx,nlayermx)      ! potential temperature (K)
+      logical firstcall,lastcall
+
+      real pday
+      real ptime 
+      logical tracerdyn
+
+!   outputs:
+!   --------
+!     physical tendencies
+      real pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
+      real pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
+      real pdpsrf(ngridmx) ! surface pressure tendency
+
+
+! Local saved variables:
+! ----------------------
+!     aerosol (dust or ice) extinction optical depth  at reference wavelength 
+!     "longrefvis" set in dimradmars.h , for one of the "naerkind"  kind of
+!      aerosol optical properties:
+!      real aerosol(ngridmx,nlayermx,naerkind) 
+!     this is now internal to callcorrk and hence no longer needed here
+
+      integer day_ini                ! Initial date of the run (sol since Ls=0) 
+      integer icount                 ! counter of calls to physiq during the run.
+      real tsurf(ngridmx)            ! Surface temperature (K)
+      real tsoil(ngridmx,nsoilmx)    ! sub-surface temperatures (K)
+      real albedo(ngridmx)           ! Surface albedo
+
+      real albedo0(ngridmx)          ! Surface albedo
+      integer rnat(ngridmx)          ! added by BC
+      save rnat
+
+      real emis(ngridmx)             ! Thermal IR surface emissivity
+      real dtrad(ngridmx,nlayermx)   ! Net atm. radiative heating rate (K.s-1)
+      real fluxrad_sky(ngridmx)      ! rad. flux from sky absorbed by surface (W.m-2)
+      real fluxrad(ngridmx)          ! Net radiative surface flux (W.m-2)
+      real capcal(ngridmx)           ! surface heat capacity (J m-2 K-1)
+      real fluxgrd(ngridmx)          ! surface conduction flux (W.m-2)
+      real qsurf(ngridmx,nqmx)       ! tracer on surface (e.g. kg.m-2)
+      real q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy 
+
+      save day_ini, icount
+      save tsurf,tsoil
+      save albedo0,albedo,emis,q2
+      save capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
+
+! Local variables :
+! -----------------
+
+!     aerosol (dust or ice) extinction optical depth  at reference wavelength 
+!     for the "naerkind" optically active aerosols:
+      real aerosol(ngridmx,nlayermx,naerkind) 
+
+      character*80 fichier 
+      integer l,ig,ierr,iq,iaer
+
+      real fluxsurf_lw(ngridmx)      ! incident LW (IR) surface flux (W.m-2)
+      real fluxsurf_sw(ngridmx)      ! incident SW (stellar) surface flux (W.m-2)
+      real fluxtop_lw(ngridmx)       ! Outgoing LW (IR) flux to space (W.m-2)
+      real fluxabs_sw(ngridmx)       ! Absorbed SW (stellar) flux (W.m-2)
+
+      real fluxtop_dn(ngridmx)
+      real fluxdyn(ngridmx)          ! horizontal heat transport by dynamics      
+      real OLR_nu(ngridmx,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
+      real OSR_nu(ngridmx,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
+ 
+      real,save :: sensibFlux(ngridmx)   ! turbulent flux given by the atm to the surface
+ 
+      save fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu
+
+
+      real zls                       ! solar longitude (rad)
+      real zday                      ! date (time since Ls=0, in martian days)
+      real zzlay(ngridmx,nlayermx)   ! altitude at the middle of the layers
+      real zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
+      real latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
+
+!     Tendencies due to various processes:
+      real dqsurf(ngridmx,nqmx)
+      real,save :: zdtlw(ngridmx,nlayermx)                    ! (K/s)
+      real,save :: zdtsw(ngridmx,nlayermx)                    ! (K/s)
+
+      real cldtlw(ngridmx,nlayermx)                           ! (K/s) LW heating rate for clear areas
+      real cldtsw(ngridmx,nlayermx)                           ! (K/s) SW heating rate for clear areas
+      real zdtsurf(ngridmx)                                   ! (K/s)
+      real dtlscale(ngridmx,nlayermx)
+      real zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx)  ! (m.s-2)
+      real zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx)         ! (K/s)
+      real zdtdif(ngridmx,nlayermx)			      ! (K/s)
+      real zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx)  ! (m.s-2)
+      real zdhadj(ngridmx,nlayermx)                           ! (K/s)
+      real zdtgw(ngridmx,nlayermx)                            ! (K/s)
+      real zdtmr(ngridmx,nlayermx)
+      real zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx)    ! (m.s-2)
+      real zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
+      real zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
+      real zdumr(ngridmx,nlayermx),zdvmr(ngridmx,nlayermx)
+      real zdtsurfmr(ngridmx)
+      
+      real zdmassmr(ngridmx,nlayermx),zdpsrfmr(ngridmx)
+
+      real zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
+      real zdqevap(ngridmx,nlayermx)
+      real zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
+      real zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
+      real zdqadj(ngridmx,nlayermx,nqmx)
+      real zdqc(ngridmx,nlayermx,nqmx)
+      real zdqmr(ngridmx,nlayermx,nqmx),zdqsurfmr(ngridmx,nqmx)
+      real zdqlscale(ngridmx,nlayermx,nqmx)
+      real zdqslscale(ngridmx,nqmx)
+      real zdqchim(ngridmx,nlayermx,nqmx)
+      real zdqschim(ngridmx,nqmx)
+
+      real zdteuv(ngridmx,nlayermx)    ! (K/s)
+      real zdtconduc(ngridmx,nlayermx) ! (K/s)
+      real zdumolvis(ngridmx,nlayermx)
+      real zdvmolvis(ngridmx,nlayermx)
+      real zdqmoldiff(ngridmx,nlayermx,nqmx)
+
+!     Local variables for local calculations:
+      real zflubid(ngridmx)
+      real zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
+      real zdum1(ngridmx,nlayermx)
+      real zdum2(ngridmx,nlayermx)
+      real ztim1,ztim2,ztim3, z1,z2
+      real ztime_fin
+      real zdh(ngridmx,nlayermx)
+      integer length
+      parameter (length=100)
+
+! local variables only used for diagnostics (output in file "diagfi" or "stats")
+! ------------------------------------------------------------------------------
+      real ps(ngridmx), zt(ngridmx,nlayermx)
+      real zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
+      real zq(ngridmx,nlayermx,nqmx)
+      character*2 str2
+      character*5 str5
+      real zdtadj(ngridmx,nlayermx)
+      real zdtdyn(ngridmx,nlayermx),ztprevious(ngridmx,nlayermx)
+      save ztprevious 
+      real reff(ngridmx,nlayermx) ! effective dust radius (used if doubleq=T)
+      real qtot1,qtot2            ! total aerosol mass
+      integer igmin, lmin
+      logical tdiag
+
+      real zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
+      real zstress(ngrid), cd
+      real hco2(nqmx), tmean, zlocal(nlayermx)
+      real vmr(ngridmx,nlayermx) ! volume mixing ratio
+
+      real time_phys
+
+!     reinstated by RW for diagnostic
+      real tau_col(ngridmx)
+      save tau_col
+      
+!     included by RW to reduce insanity of code
+      real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS
+
+!     included by RW to compute tracer column densities
+      real qcol(ngridmx,nqmx)
+
+!     included by RW for H2O precipitation
+      real zdtrain(ngridmx,nlayermx)
+      real zdqrain(ngridmx,nlayermx,nqmx)
+      real zdqsrain(ngridmx)
+      real zdqssnow(ngridmx)
+      real rainout(ngridmx)
+
+!     included by RW for H2O Manabe scheme
+      real dtmoist(ngridmx,nlayermx)
+      real dqmoist(ngridmx,nlayermx,nqmx)
+
+      real qvap(ngridmx,nlayermx)
+      real dqvaplscale(ngridmx,nlayermx)
+      real dqcldlscale(ngridmx,nlayermx)
+      real rneb_man(ngridmx,nlayermx)
+      real rneb_lsc(ngridmx,nlayermx)
+
+!     included by RW to account for surface cooling by evaporation
+      real dtsurfh2olat(ngridmx)
+
+
+!     to test energy conservation (RW+JL)
+      real mass(ngridmx,nlayermx),massarea(ngridmx,nlayermx)
+      real dEtot, dEtots, AtmToSurf_TurbFlux
+      real dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
+      real dEzRadsw(ngridmx,nlayermx),dEzRadlw(ngridmx,nlayermx),dEzdiff(ngridmx,nlayermx)
+      real dEdiffs(ngridmx),dEdiff(ngridmx)
+      real madjdE(ngridmx), lscaledE(ngridmx)
+!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
+      
+      real dItot, dVtot
+
+!     included by BC for evaporation     
+      real qevap(ngridmx,nlayermx,nqmx)
+      real tevap(ngridmx,nlayermx)
+      real dqevap1(ngridmx,nlayermx)
+      real dtevap1(ngridmx,nlayermx)
+
+!     included by BC for hydrology
+      real hice(ngridmx)
+
+!     included by RW to test water conservation (by routine)
+      real h2otot
+      real dWtot, dWtots
+      real h2o_surf_all
+      logical watertest
+      save watertest
+
+!     included by RW for RH diagnostic
+      real qsat(ngridmx,nlayermx), RH(ngridmx,nlayermx), H2Omaxcol(ngridmx),psat_tmp
+
+!     included by RW for hydrology
+      real dqs_hyd(ngridmx,nqmx)
+      real zdtsurf_hyd(ngridmx)
+
+!     included by RW for water cycle conservation tests
+      real icesrf,liqsrf,icecol,vapcol
+
+!     included by BC for double radiative transfer call
+      logical clearsky
+      real zdtsw1(ngridmx,nlayermx)
+      real zdtlw1(ngridmx,nlayermx)
+      real fluxsurf_lw1(ngridmx)     
+      real fluxsurf_sw1(ngridmx)  
+      real fluxtop_lw1(ngridmx)    
+      real fluxabs_sw1(ngridmx)  
+      real tau_col1(ngrid)
+      real OLR_nu1(ngrid,L_NSPECTI)
+      real OSR_nu1(ngrid,L_NSPECTV)
+      real tf, ntf
+
+!     included by BC for cloud fraction computations
+      real,save :: cloudfrac(ngridmx,nlayermx)
+      real,save :: totcloudfrac(ngridmx)
+
+!     included by RW for vdifc water conservation test
+      real nconsMAX
+      real vdifcncons(ngridmx)
+      real cadjncons(ngridmx)
+
+!      double precision qsurf_hist(ngridmx,nqmx)
+      real qsurf_hist(ngridmx,nqmx)
+      save qsurf_hist
+
+!     included by RW for temp convadj conservation test
+      real playtest(nlayermx)
+      real plevtest(nlayermx)
+      real ttest(nlayermx)
+      real qtest(nlayermx)
+      integer igtest
+
+!     included by RW for runway greenhouse 1D study
+      real muvar(ngridmx,nlayermx+1)
+
+!     included by RW for variable H2O particle sizes
+      real,save :: reffrad(ngridmx,nlayermx,naerkind) ! aerosol effective radius (m)
+      real,save :: nueffrad(ngridmx,nlayermx,naerkind) ! aerosol effective radius variance
+      real :: reffh2oliq(ngridmx,nlayermx) ! liquid water particles effective radius (m)
+      real :: reffh2oice(ngridmx,nlayermx) ! water ice particles effective radius (m)
+      real reffH2O(ngridmx,nlayermx)
+      real reffcol(ngridmx,naerkind)
+
+!     included by RW for sourceevol
+      real, save :: ice_initial(ngridmx)!, isoil
+      real delta_ice,ice_tot
+      real, save :: ice_min(ngridmx)
+      
+      integer num_run
+      logical,save :: ice_update
+      
+!=======================================================================
+      
+
+! 1. Initialisation
+! -----------------
+
+!  1.1   Initialisation only at first call
+!  ---------------------------------------
+      if (firstcall) then
+
+
+!        variables set to 0
+!        ~~~~~~~~~~~~~~~~~~
+         dtrad(:,:) = 0.0
+         fluxrad(:) = 0.0
+         tau_col(:) = 0.0
+         zdtsw(:,:) = 0.0
+         zdtlw(:,:) = 0.0
+
+!        initialize aerosol indexes 
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            call iniaerosol()
+
+      
+!        initialize tracer names, indexes and properties
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         tracerdyn=tracer
+         if (tracer) then
+            call initracer()
+         endif ! end tracer
+
+! 
+
+!        read startfi (initial parameters)
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         call phyetat0("startfi.nc",0,0,nsoilmx,nq,           &
+               day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,   &
+               cloudfrac,totcloudfrac,hice)
+
+         if (pday.ne.day_ini) then
+           write(*,*) "ERROR in physiq.F90:"
+           write(*,*) "bad synchronization between physics and dynamics"
+           write(*,*) "dynamics day: ",pday
+           write(*,*) "physics day:  ",day_ini
+           stop
+         endif
+
+         write (*,*) 'In physiq day_ini =', day_ini
+
+!        Initialize albedo and orbital calculation
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         call surfini(ngrid,qsurf,albedo0)
+         call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
+
+         albedo(:)=albedo0(:)
+
+         if(tlocked)then
+            print*,'Planet is tidally locked at resonance n=',nres
+            print*,'Make sure you have the right rotation rate!!!'
+         endif
+
+!        initialize soil 
+!        ~~~~~~~~~~~~~~~
+         if (callsoil) then
+            call soil(ngrid,nsoilmx,firstcall,inertiedat,     &
+                ptimestep,tsurf,tsoil,capcal,fluxgrd)
+         else
+            print*,'WARNING! Thermal conduction in the soil turned off'
+            capcal(:)=1.e6
+            fluxgrd(:)=0.
+            if(noradsurf)then
+                  fluxgrd(:)=10.0
+            endif
+            print*,'Flux from ground = ',fluxgrd,' W m^-2'
+         endif
+         icount=1
+
+!        decide whether to update ice at end of run
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         ice_update=.false.
+         if(sourceevol)then
+            open(128,file='num_run',form='formatted')
+            read(128,*) num_run 
+            close(128)
+        
+            if(num_run.ne.0.and.mod(num_run,2).eq.0)then
+            !if(num_run.ne.0.and.mod(num_run,3).eq.0)then
+               print*,'Updating ice at end of this year!'
+               ice_update=.true.
+               ice_min(:)=1.e4
+            endif  
+         endif
+
+!        define surface as continent or ocean
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         rnat(:)=1
+         do ig=1,ngridmx
+!            if(iceball.or.oceanball.or.(inertiedat(ig,1).gt.1.E4))then
+            if(inertiedat(ig,1).gt.1.E4)then
+            rnat(ig)=0
+            endif
+         enddo
+
+         print*,'WARNING! Surface type currently decided by surface inertia'
+         print*,'This should be improved e.g. in newstart.F'
+
+
+!        initialise surface history variable 
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         qsurf_hist(:,:)=qsurf(:,:)
+
+!        initialise variable for dynamical heating diagnostic
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         ztprevious(:,:)=pt(:,:)
+
+!        Set temperature just above condensation temperature (for Early Mars)
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         if(nearco2cond) then
+            write(*,*)' WARNING! Starting at Tcond+1K'
+            do l=1, nlayer
+               do ig=1,ngrid
+                  pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4     &
+                      -pt(ig,l)) / ptimestep
+               enddo
+            enddo
+         endif
+
+         if(meanOLR)then
+            ! to record global radiative balance
+            call system('rm -f rad_bal.out')
+            ! to record global mean/max/min temperatures
+            call system('rm -f tem_bal.out')
+            ! to record global hydrological balance
+            call system('rm -f h2o_bal.out')
+         endif
+
+         watertest=.false.
+         if(water)then
+            ! initialise variables for water cycle
+
+            if(enertest)then
+               watertest = .true.
+            endif
+
+            if(ice_update)then
+               ice_initial(:)=qsurf(:,igcm_h2o_ice)
+            endif
+
+         endif
+         call su_watercycle ! even if we don't have a water cycle, we might
+                            ! need epsi for the wvp definitions in callcorrk.F
+
+      endif        !  (end of "if firstcall")
+
+! ---------------------------------------------------
+! 1.2   Initializations done at every physical timestep:
+! ---------------------------------------------------
+
+      if (ngrid.NE.ngridmx) then
+         print*,'STOP in PHYSIQ'
+         print*,'Probleme de dimensions :'
+         print*,'ngrid     = ',ngrid
+         print*,'ngridmx   = ',ngridmx
+         stop
+      endif
+
+!     Initialize various variables
+!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      pdu(1:ngridmx,1:nlayermx) = 0.0
+      pdv(1:ngridmx,1:nlayermx) = 0.0
+      if ( .not.nearco2cond ) then
+         pdt(1:ngridmx,1:nlayermx) = 0.0
+      endif
+      pdq(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
+      pdpsrf(1:ngridmx)       = 0.0
+      zflubid(1:ngridmx)      = 0.0
+      zdtsurf(1:ngridmx)      = 0.0
+      dqsurf(1:ngridmx,1:nqmx)= 0.0
+
+      zday=pday+ptime ! compute time, in sols (and fraction thereof)
+
+!     Compute Stellar Longitude (Ls)
+!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      if (season) then
+         call stellarlong(zday,zls)
+      else
+         call stellarlong(float(day_ini),zls)
+      end if
+
+!     Compute geopotential between layers
+!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      zzlay(1:ngridmx,1:nlayermx)=pphi(1:ngridmx,1:nlayermx)/g
+
+      zzlev(1:ngridmx,1)=0.
+      zzlev(1:ngridmx,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
+
+      do l=2,nlayer
+         do ig=1,ngrid
+            z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
+            z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
+            zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
+         enddo
+      enddo
+!     Potential temperature calculation may not be the same in physiq and dynamic...
+
+!     Compute potential temperature
+!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      do l=1,nlayer         
+         do ig=1,ngrid 
+            zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
+            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
+            mass(ig,l)  = (pplev(ig,l) - pplev(ig,l+1))/g
+	    massarea(ig,l)=mass(ig,l)*area(ig)
+         enddo
+      enddo
+
+!-----------------------------------------------------------------------
+!    2. Compute radiative tendencies
+!-----------------------------------------------------------------------
+
+      if (callrad) then
+         if( mod(icount-1,iradia).eq.0.or.lastcall) then
+
+!          Compute local stellar zenith angles
+!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           call orbite(zls,dist_star,declin)
+
+           if (tlocked) then
+              ztim1=SIN(declin)
+              ztim2=COS(declin)*COS(2.*pi*(zday/year_day) - zls*nres)
+              ztim3=-COS(declin)*SIN(2.*pi*(zday/year_day) - zls*nres)
+
+              call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
+               ztim1,ztim2,ztim3,mu0,fract)
+
+           elseif (diurnal) then
+               ztim1=SIN(declin)
+               ztim2=COS(declin)*COS(2.*pi*(zday-.5))
+               ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
+
+               call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
+               ztim1,ztim2,ztim3,mu0,fract)
+
+           else
+
+               call mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
+               ! WARNING: this function appears not to work in 1D
+
+           endif
+
+           if (corrk) then
+
+!          a) Call correlated-k radiative transfer scheme
+!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+              if(kastprof)then
+                 print*,'kastprof should not = true here'
+                 call abort
+              endif
+              muvar(:,:)=0.0 ! only used for climate evolution studies (kcm1d) for now
+      
+!             standard callcorrk
+              clearsky=.false.
+              call callcorrk(ngrid,nlayer,pq,nq,qsurf,                   &
+                      albedo,emis,mu0,pplev,pplay,pt,                    &
+                      tsurf,fract,dist_star,aerosol,muvar,               &
+                      zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
+                      fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,               &
+                      reffrad,tau_col,cloudfrac,totcloudfrac,            &
+                      clearsky,firstcall,lastcall)      
+
+!             Option to call scheme once more for clear regions
+              if(CLFvarying)then
+
+                 ! ---> PROBLEMS WITH ALLOCATED ARRAYS 
+                 ! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
+                 clearsky=.true.
+                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,                  &
+                      albedo,emis,mu0,pplev,pplay,pt,                      &
+                      tsurf,fract,dist_star,aerosol,muvar,                 &
+                      zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
+                      fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,              &
+                      reffrad,tau_col1,cloudfrac,totcloudfrac,             &
+                      clearsky,firstcall,lastcall)
+                 clearsky = .false.  ! just in case.     
+
+                 ! Sum the fluxes and heating rates from cloudy/clear cases
+                 do ig=1,ngrid
+                    tf=totcloudfrac(ig)
+                    ntf=1.-tf
+                    
+                    fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw(ig)
+                    fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw(ig)
+                    fluxtop_lw(ig)  = ntf*fluxtop_lw1(ig)  + tf*fluxtop_lw(ig)
+                    fluxabs_sw(ig)  = ntf*fluxabs_sw1(ig)  + tf*fluxabs_sw(ig)
+                    tau_col(ig)     = ntf*tau_col1(ig)     + tf*tau_col(ig)
+                    
+                    zdtlw(ig,1:nlayermx) = ntf*zdtlw1(ig,1:nlayermx) + tf*zdtlw(ig,1:nlayermx)
+                    zdtsw(ig,1:nlayermx) = ntf*zdtsw1(ig,1:nlayermx) + tf*zdtsw(ig,1:nlayermx)
+
+                    OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV)		       
+                    OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI)		       
+
+                 enddo
+
+              endif !CLFvarying
+
+!             Radiative flux from the sky absorbed by the surface (W.m-2)
+              GSR=0.0
+              fluxrad_sky(1:ngridmx)=emis(1:ngridmx)*fluxsurf_lw(1:ngridmx)+fluxsurf_sw(1:ngridmx)*(1.-albedo(1:ngridmx))
+
+              if(noradsurf)then ! no lower surface; SW flux just disappears
+                 GSR = SUM(fluxsurf_sw(1:ngridmx)*area(1:ngridmx))/totarea
+                 fluxrad_sky(1:ngridmx)=emis(1:ngridmx)*fluxsurf_lw(1:ngridmx)
+                 print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
+              endif
+
+!             Net atmospheric radiative heating rate (K.s-1)
+              dtrad(1:ngridmx,1:nlayermx)=zdtsw(1:ngridmx,1:nlayermx)+zdtlw(1:ngridmx,1:nlayermx)
+
+           elseif(newtonian)then
+
+!          b) Call Newtonian cooling scheme
+!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+              call newtrelax(mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
+
+              zdtsurf(1:ngridmx) = +(pt(1:ngridmx,1)-tsurf(1:ngridmx))/ptimestep
+              ! e.g. surface becomes proxy for 1st atmospheric layer ?
+
+           else
+
+!          c) Atmosphere has no radiative effect
+!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+              fluxtop_dn(1:ngridmx)  = fract(1:ngridmx)*mu0(1:ngridmx)*Fat1AU/dist_star**2
+              if(ngrid.eq.1)then ! / by 4 globally in 1D case...
+                 fluxtop_dn(1)  = fract(1)*Fat1AU/dist_star**2/2.0
+              endif
+              fluxsurf_sw(1:ngridmx) = fluxtop_dn(1:ngridmx)
+              fluxrad_sky(1:ngridmx) = fluxtop_dn(1:ngridmx)*(1.-albedo(1:ngridmx))
+              fluxtop_lw(1:ngridmx)  = emis(1:ngridmx)*sigma*tsurf(1:ngridmx)**4
+              ! radiation skips the atmosphere entirely
+
+
+              dtrad(1:ngridmx,1:nlayermx)=0.0
+              ! hence no atmospheric radiative heating
+
+           endif                ! if corrk
+
+        endif ! of if(mod(icount-1,iradia).eq.0)
+
+
+!    Transformation of the radiative tendencies
+!    ------------------------------------------
+
+        zplanck(1:ngridmx)=tsurf(1:ngridmx)*tsurf(1:ngridmx)
+        zplanck(1:ngridmx)=emis(1:ngridmx)*sigma*zplanck(1:ngridmx)*zplanck(1:ngridmx)
+        fluxrad(1:ngridmx)=fluxrad_sky(1:ngridmx)-zplanck(1:ngridmx)
+        pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+dtrad(1:ngridmx,1:nlayermx)
+
+!-------------------------
+! test energy conservation
+         if(enertest)then
+            dEtotSW  = cpp*SUM(massarea(:,:)*zdtsw(:,:))/totarea
+            dEtotLW  = cpp*SUM(massarea(:,:)*zdtlw(:,:))/totarea
+            dEtotsSW = SUM(fluxsurf_sw(:)*(1.-albedo(:))*area(:))/totarea
+            dEtotsLW = SUM((fluxsurf_lw(:)*emis(:)-zplanck(:))*area(:))/totarea
+	    dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:)
+	    dEzRadlw(:,:)=cpp*mass(:,:)*zdtlw(:,:)
+            print*,'---------------------------------------------------------------'
+            print*,'In corrk SW atmospheric heating       =',dEtotSW,' W m-2'
+            print*,'In corrk LW atmospheric heating       =',dEtotLW,' W m-2'
+            print*,'atmospheric net rad heating (SW+LW)   =',dEtotLW+dEtotSW,' W m-2'
+            print*,'In corrk SW surface heating           =',dEtotsSW,' W m-2'
+            print*,'In corrk LW surface heating           =',dEtotsLW,' W m-2'
+            print*,'surface net rad heating (SW+LW)       =',dEtotsLW+dEtotsSW,' W m-2'
+         endif
+!-------------------------
+
+      endif ! of if (callrad)
+
+!-----------------------------------------------------------------------
+!    4. Vertical diffusion (turbulent mixing):
+!    -----------------------------------------
+
+      if (calldifv) then
+      
+         zflubid(1:ngridmx)=fluxrad(1:ngridmx)+fluxgrd(1:ngridmx)
+
+         zdum1(1:ngridmx,1:nlayermx)=0.0
+         zdum2(1:ngridmx,1:nlayermx)=0.0
+
+
+!JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff	 
+	 if (UseTurbDiff) then
+	 
+           call turbdiff(ngrid,nlayer,nq,rnat,       &
+              ptimestep,capcal,lwrite,               &
+              pplay,pplev,zzlay,zzlev,z0,            &
+              pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf,   &
+              zdum1,zdum2,pdt,pdq,zflubid,           &
+              zdudif,zdvdif,zdtdif,zdtsdif,          &
+	      sensibFlux,q2,zdqdif,zdqevap,zdqsdif,lastcall)
+
+	 else
+	 
+           zdh(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)/zpopsk(1:ngridmx,1:nlayermx)
+ 
+           call vdifc(ngrid,nlayer,nq,rnat,zpopsk,   &
+              ptimestep,capcal,lwrite,               &
+              pplay,pplev,zzlay,zzlev,z0,            &
+              pu,pv,zh,pq,tsurf,emis,qsurf,          &
+              zdum1,zdum2,zdh,pdq,zflubid,           &
+              zdudif,zdvdif,zdhdif,zdtsdif,          &
+	      sensibFlux,q2,zdqdif,zdqsdif,lastcall)
+
+           zdtdif(1:ngridmx,1:nlayermx)=zdhdif(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx) ! for diagnostic only
+	   zdqevap(1:ngridmx,1:nlayermx)=0.
+
+         end if !turbdiff
+
+         pdv(1:ngridmx,1:nlayermx)=pdv(1:ngridmx,1:nlayermx)+zdvdif(1:ngridmx,1:nlayermx)
+         pdu(1:ngridmx,1:nlayermx)=pdu(1:ngridmx,1:nlayermx)+zdudif(1:ngridmx,1:nlayermx)
+         pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+zdtdif(1:ngridmx,1:nlayermx)
+         zdtsurf(1:ngridmx)=zdtsurf(1:ngridmx)+zdtsdif(1:ngridmx)
+         if (tracer) then 
+           pdq(1:ngridmx,1:nlayermx,1:nqmx)=pdq(1:ngridmx,1:nlayermx,1:nqmx)+ zdqdif(1:ngridmx,1:nlayermx,1:nqmx) 
+           dqsurf(1:ngridmx,1:nqmx)=dqsurf(1:ngridmx,1:nqmx) + zdqsdif(1:ngridmx,1:nqmx)
+         end if ! of if (tracer)
+
+         !-------------------------
+         ! test energy conservation
+         if(enertest)then
+	    dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:)
+            do ig = 1, ngrid
+	       dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
+	       dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
+            enddo
+	    dEtot = SUM(dEdiff(:)*area(:))/totarea
+	    dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
+            dEtots = SUM(dEdiffs(:)*area(:))/totarea
+	    AtmToSurf_TurbFlux=SUM(sensibFlux(:)*area(:))/totarea
+            if (UseTurbDiff) then
+	       print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
+               print*,'In TurbDiff non-cons atm nrj change   =',dEtot,' W m-2'
+               print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
+	    else
+	       print*,'In vdifc sensible flux (atm=>surf)    =',AtmToSurf_TurbFlux,' W m-2'
+               print*,'In vdifc non-cons atm nrj change      =',dEtot,' W m-2'
+               print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
+	    end if
+! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme
+!    but not given back elsewhere
+         endif
+         !-------------------------
+
+         !-------------------------
+         ! test water conservation
+         if(watertest.and.water)then
+            dWtot = SUM(massarea(:,:)*zdqdif(:,:,igcm_h2o_vap))*ptimestep/totarea
+            dWtots =  SUM(zdqsdif(:,igcm_h2o_vap)*area(:))*ptimestep/totarea
+            do ig = 1, ngrid
+	       vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap))
+	    Enddo
+            dWtot = dWtot + SUM(massarea(:,:)*zdqdif(:,:,igcm_h2o_ice))*ptimestep/totarea
+            dWtots = dWtots + SUM(zdqsdif(:,igcm_h2o_ice)*area(:))*ptimestep/totarea
+            do ig = 1, ngrid
+	       vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice))
+	    Enddo	    
+            nconsMAX=MAXVAL(vdifcncons(:))
+
+            print*,'---------------------------------------------------------------'
+            print*,'In difv atmospheric water change        =',dWtot,' kg m-2'
+            print*,'In difv surface water change            =',dWtots,' kg m-2'
+            print*,'In difv non-cons factor                 =',dWtot+dWtots,' kg m-2'
+            print*,'In difv MAX non-cons factor             =',nconsMAX,' kg m-2 s-1'
+
+         endif
+         !-------------------------
+
+      else    
+
+         if(.not.newtonian)then
+
+            zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + (fluxrad(1:ngridmx) + fluxgrd(1:ngridmx))/capcal(1:ngridmx)
+
+         endif
+
+      endif ! of if (calldifv)
+
+
+!-----------------------------------------------------------------------
+!   5. Dry convective adjustment:
+!   -----------------------------
+
+      if(calladj) then
+
+         zdh(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)/zpopsk(1:ngridmx,1:nlayermx)
+         zduadj(1:ngridmx,1:nlayermx)=0.0
+         zdvadj(1:ngridmx,1:nlayermx)=0.0
+         zdhadj(1:ngridmx,1:nlayermx)=0.0
+
+
+         call convadj(ngrid,nlayer,nq,ptimestep,    &
+              pplay,pplev,zpopsk,                   &
+              pu,pv,zh,pq,                          &
+              pdu,pdv,zdh,pdq,                      &
+              zduadj,zdvadj,zdhadj,                 &
+              zdqadj)
+
+         pdu(1:ngridmx,1:nlayermx) = pdu(1:ngridmx,1:nlayermx) + zduadj(1:ngridmx,1:nlayermx)
+         pdv(1:ngridmx,1:nlayermx) = pdv(1:ngridmx,1:nlayermx) + zdvadj(1:ngridmx,1:nlayermx)
+         pdt(1:ngridmx,1:nlayermx)    = pdt(1:ngridmx,1:nlayermx) + zdhadj(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx)
+         zdtadj(1:ngridmx,1:nlayermx) = zdhadj(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx) ! for diagnostic only
+ 
+         if(tracer) then 
+            pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqadj(1:ngridmx,1:nlayermx,1:nqmx) 
+         end if
+
+         !-------------------------
+         ! test energy conservation
+         if(enertest)then
+            dEtot=cpp*SUM(massarea(:,:)*zdtadj(:,:))/totarea
+          print*,'In convadj atmospheric energy change  =',dEtot,' W m-2'
+         endif
+         !-------------------------
+
+         !-------------------------
+         ! test water conservation
+         if(watertest)then
+            dWtot = SUM(massarea(:,:)*zdqadj(:,:,igcm_h2o_vap))*ptimestep/totarea
+            do ig = 1, ngrid
+	       cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap))
+	    Enddo
+            dWtot = dWtot + SUM(massarea(:,:)*zdqadj(:,:,igcm_h2o_ice))*ptimestep/totarea
+            do ig = 1, ngrid
+	       cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice))
+	    Enddo	    
+            nconsMAX=MAXVAL(cadjncons(:))
+
+            print*,'In convadj atmospheric water change     =',dWtot,' kg m-2'
+            print*,'In convadj MAX non-cons factor          =',nconsMAX,' kg m-2 s-1'
+         endif
+         !-------------------------
+               
+      endif ! of if(calladj)
+
+!-----------------------------------------------------------------------
+!   6. Carbon dioxide condensation-sublimation:
+!   -------------------------------------------
+
+      if (co2cond) then
+         if (.not.tracer) then
+            print*,'We need a CO2 ice tracer to condense CO2'
+            call abort
+         endif
+     print*, 'we are in co2cond!!!!!'
+         call condense_cloud(ngrid,nlayer,nq,ptimestep,   &
+              capcal,pplay,pplev,tsurf,pt,                  &
+              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,        &
+              qsurf(1,igcm_co2_ice),albedo,emis,            &
+              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,               &
+              zdqc,reffrad)
+
+
+         pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+zdtc(1:ngridmx,1:nlayermx)
+         pdv(1:ngridmx,1:nlayermx)=pdv(1:ngridmx,1:nlayermx)+zdvc(1:ngridmx,1:nlayermx)
+         pdu(1:ngridmx,1:nlayermx)=pdu(1:ngridmx,1:nlayermx)+zduc(1:ngridmx,1:nlayermx)
+         zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + zdtsurfc(1:ngridmx)
+
+         pdq(1:ngridmx,1:nlayermx,1:nqmx)=pdq(1:ngridmx,1:nlayermx,1:nqmx)+ zdqc(1:ngridmx,1:nlayermx,1:nqmx) 
+         ! Note: we do not add surface co2ice tendency
+         ! because qsurf(:,igcm_co2_ice) is updated in condens_co2cloud
+
+         !-------------------------
+         ! test energy conservation
+         if(enertest)then
+            dEtot = cpp*SUM(massarea(:,:)*zdtc(:,:))/totarea
+            dEtots = SUM(capcal(:)*zdtsurfc(:)*area(:))/totarea
+            print*,'In co2cloud atmospheric energy change   =',dEtot,' W m-2'
+            print*,'In co2cloud surface energy change       =',dEtots,' W m-2'
+         endif
+         !-------------------------
+
+      endif  ! of if (co2cond)
+
+
+!-----------------------------------------------------------------------
+!   7. Specific parameterizations for tracers 
+!   -----------------------------------------
+
+      if (tracer) then 
+
+!   7a. Water and ice
+!     ---------------
+         if (water) then
+
+!        ----------------------------------------
+!        Water ice condensation in the atmosphere
+!        ----------------------------------------
+            if(watercond.and.(RLVTT.gt.1.e-8))then
+
+!             ----------------
+!             Moist convection
+!             ----------------
+               ! Re-evaporate cloud water/ice
+!	       dqevap1(1:ngridmx,1:nlayermx)=0.
+!	       dtevap1(1:ngridmx,1:nlayermx)=0.
+!               call evap(ptimestep,pt,pq,pdq,pdt,dqevap1,dtevap1,qevap,tevap)
+               
+	       dqmoist(1:ngridmx,1:nlayermx,1:nqmx)=0.
+	       dtmoist(1:ngridmx,1:nlayermx)=0.
+
+               call moistadj(pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
+
+               pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap)   &
+	                  +dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_vap)
+               pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) =pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice)     &
+                          +dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_ice)
+               pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+dtmoist(1:ngridmx,1:nlayermx)
+
+!               do l=1,nlayermx
+!	         do ig=1,ngridmx
+!	          if(isnan(dqmoist(ig,l,igcm_h2o_ice))) print*,'Nan in dqmoist ice, abort,ig,l',ig,l
+!	          if(isnan(dqmoist(ig,l,igcm_h2o_vap))) print*,'Nan in dqmoist vap, abort,ig,l',ig,l
+!	          if(isnan(dtmoist(ig,l))) print*,'Nan in dtmoist, abort,ig,l',ig,l		  
+!                  if(isnan(dqmoist(ig,l,igcm_h2o_ice)).or.isnan(dqmoist(ig,l,igcm_h2o_vap)).or.isnan(dtmoist(ig,l))) STOP
+!	         enddo
+!	       enddo
+!	       print*,'after moist',dqmoist(185,1,igcm_h2o_ice),dqmoist(185,1,igcm_h2o_vap),dtmoist(185,1)*86400.
+
+               !-------------------------
+               ! test energy conservation
+               if(enertest)then
+                  dEtot=cpp*SUM(massarea(:,:)*dtmoist(:,:))/totarea
+                  do ig = 1, ngrid
+                     madjdE(ig) = cpp*SUM(mass(:,:)*dtmoist(:,:))
+                  enddo
+                  print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
+                  print*,'In moistadj MAX atmospheric energy change   =',MAXVAL(pdt(:,:))*86400.,'K/day'
+		  
+		! test energy conservation
+	          dWtot = SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap))*ptimestep/totarea
+                  dWtot = dWtot + SUM(massarea(:,:)*dqmoist(:,:,igcm_h2o_ice))*ptimestep/totarea
+                  print*,'In moistadj atmospheric water change    =',dWtot,' kg m-2'
+               endif
+               !-------------------------
+	       
+
+!        --------------------------------
+!        Large scale condensation/evaporation
+!        --------------------------------
+
+               call largescale(ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc)
+
+               pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+dtlscale(1:ngridmx,1:nlayermx)
+               pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap)+dqvaplscale(1:ngridmx,1:nlayermx)
+               pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice)+dqcldlscale(1:ngridmx,1:nlayermx)
+
+               !-------------------------
+               ! test energy conservation
+               if(enertest)then
+                  do ig = 1, ngrid
+                     lscaledE(ig) = cpp*SUM(mass(:,:)*dtlscale(:,:))
+                  enddo
+                  dEtot=cpp*SUM(massarea(:,:)*(dtlscale(:,:)))/totarea
+		  if(isnan(dEtot)) then
+		     print*,'Nan in largescale, abort'
+                     STOP
+		  endif
+		  print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
+
+               ! test water conservation
+	          dWtot = SUM(massarea(:,:)*dqvaplscale(:,:))*ptimestep/totarea
+                  dWtot = dWtot + SUM(massarea(:,:)*dqcldlscale(:,:))*ptimestep/totarea
+                  print*,'In largescale atmospheric water change  =',dWtot,' kg m-2'
+               endif
+               !-------------------------
+
+               ! compute cloud fraction
+               do l = 1, nlayer
+                  do ig = 1,ngrid
+                     cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l)) 
+                  enddo
+               enddo
+
+            endif                ! of if (watercondense)
+           
+
+!        --------------------------------
+!        Water ice / liquid precipitation
+!        --------------------------------
+            if(waterrain)then
+
+               zdqrain(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
+               zdqsrain(1:ngridmx)    = 0.0
+               zdqssnow(1:ngridmx)    = 0.0
+
+               call rain(ptimestep,pplev,pplay,pt,pdt,pq,pdq,            &
+                   zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
+
+               pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) &
+	             +zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_vap)
+               pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) &
+	             +zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_ice)
+               pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+zdtrain(1:ngridmx,1:nlayermx)
+               dqsurf(1:ngridmx,igcm_h2o_vap) = dqsurf(1:ngridmx,igcm_h2o_vap)+zdqsrain(1:ngridmx) ! a bug was here
+               dqsurf(1:ngridmx,igcm_h2o_ice) = dqsurf(1:ngridmx,igcm_h2o_ice)+zdqssnow(1:ngridmx)
+               rainout(1:ngridmx)             = zdqsrain(1:ngridmx)+zdqssnow(1:ngridmx) ! diagnostic   
+
+
+               !-------------------------
+               ! test energy conservation
+               if(enertest)then
+                  dEtot=cpp*SUM(massarea(:,:)*zdtrain(:,:))/totarea
+                 print*,'In rain atmospheric T energy change       =',dEtot,' W m-2'
+                  dItot = SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_ice))/totarea*RLVTT/cpp
+		  dItot = dItot + SUM(area(:)*zdqssnow(:))/totarea*RLVTT/cpp
+                  dVtot = SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap))*ptimestep/totarea
+		  dVtot = dVtot + SUM(area(:)*zdqsrain(:))/totarea*RLVTT/cpp
+		  dEtot = dItot + dVtot
+                 print*,'In rain dItot =',dItot,' W m-2'
+                 print*,'In rain dVtot =',dVtot,' W m-2'
+                 print*,'In rain atmospheric L energy change       =',dEtot,' W m-2'
+
+               ! test water conservation
+	          dWtot = SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap))*ptimestep/totarea
+                  dWtot = dWtot + SUM(massarea(:,:)*zdqrain(:,:,igcm_h2o_ice))*ptimestep/totarea
+		  dWtots =  SUM((zdqsrain(:)+zdqssnow(:))*area(:))*ptimestep/totarea
+                 print*,'In rain atmospheric water change        =',dWtot,' kg m-2'
+                 print*,'In rain surface water change            =',dWtots,' kg m-2'
+                 print*,'In rain non-cons factor                 =',dWtot+dWtots,' kg m-2'
+              endif
+              !-------------------------
+
+            end if                 ! of if (waterrain)
+         end if                    ! of if (water)
+
+
+!   7c. Aerosol particles
+!     -------------------
+!        ------------- 
+!        Sedimentation
+!        ------------- 
+        if (sedimentation) then 
+           zdqsed(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
+           zdqssed(1:ngridmx,1:nqmx)  = 0.0
+
+
+           !-------------------------
+           ! find qtot
+           if(watertest)then
+              iq=3
+	      dWtot = SUM(massarea(:,:)*pq(:,:,iq))*ptimestep/totarea
+	      dWtots = SUM(massarea(:,:)*pdq(:,:,iq))*ptimestep/totarea
+              print*,'Before sedim pq  =',dWtot,' kg m-2'
+              print*,'Before sedim pdq =',dWtots,' kg m-2'
+           endif
+           !-------------------------
+
+           call callsedim(ngrid,nlayer,ptimestep,           &
+                pplev,zzlev,pt,pq,pdq,zdqsed,zdqssed,nq,reffH2O)
+
+           !-------------------------
+           ! find qtot
+           if(watertest)then
+              iq=3
+	      dWtot = SUM(massarea(:,:)*pq(:,:,iq))*ptimestep/totarea
+	      dWtots = SUM(massarea(:,:)*pdq(:,:,iq))*ptimestep/totarea
+              print*,'After sedim pq  =',dWtot,' kg m-2'
+              print*,'After sedim pdq =',dWtots,' kg m-2'
+           endif
+           !-------------------------
+
+           ! for now, we only allow H2O ice to sediment
+              ! and as in rain.F90, whether it falls as rain or snow depends
+              ! only on the surface temperature
+           pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqsed(1:ngridmx,1:nlayermx,1:nqmx)
+           dqsurf(1:ngridmx,1:nqmx) = dqsurf(1:ngridmx,1:nqmx) + zdqssed(1:ngridmx,1:nqmx)
+
+           !-------------------------
+           ! test water conservation
+           if(watertest)then
+              dWtot = SUM(massarea(:,:)*(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice)))*ptimestep/totarea
+	      dWtots = SUM((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))*area(:))*ptimestep/totarea
+              print*,'In sedim atmospheric ice change         =',dWtot,' kg m-2'
+              print*,'In sedim surface ice change             =',dWtots,' kg m-2'
+              print*,'In sedim non-cons factor                =',dWtot+dWtots,' kg m-2'
+           endif
+           !-------------------------
+
+        end if   ! of if (sedimentation)
+
+
+!   7d. Updates
+!     ---------
+
+!       -----------------------------------
+!       Updating atm mass and tracer budget
+!       -----------------------------------
+
+         if(mass_redistrib) then
+
+            zdmassmr(1:ngridmx,1:nlayermx) = mass(1:ngridmx,1:nlayermx) * &
+	       (   zdqevap(1:ngridmx,1:nlayermx)                          &
+	         + zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_vap)             &
+		 + dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_vap)             &
+		 + dqvaplscale(1:ngridmx,1:nlayermx) )
+		 
+            
+	    call writediagfi(ngridmx,"mass_evap","mass gain"," ",3,zdmassmr)
+	    call writediagfi(ngridmx,"mass","mass"," ",3,mass)
+
+	    call mass_redistribution(ngrid,nlayer,nq,ptimestep,   &
+                       rnat,capcal,pplay,pplev,pt,tsurf,pq,qsurf,     &
+		       pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr,  &
+		       zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
+	 
+
+            pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqmr(1:ngridmx,1:nlayermx,1:nqmx)
+            dqsurf(1:ngridmx,1:nqmx)         = dqsurf(1:ngridmx,1:nqmx) + zdqsurfmr(1:ngridmx,1:nqmx)
+            pdt(1:ngridmx,1:nlayermx)        = pdt(1:ngridmx,1:nlayermx) + zdtmr(1:ngridmx,1:nlayermx)
+            pdu(1:ngridmx,1:nlayermx)        = pdu(1:ngridmx,1:nlayermx) + zdumr(1:ngridmx,1:nlayermx)
+            pdv(1:ngridmx,1:nlayermx)        = pdv(1:ngridmx,1:nlayermx) + zdvmr(1:ngridmx,1:nlayermx)
+	    pdpsrf(1:ngridmx)                = pdpsrf(1:ngridmx) + zdpsrfmr(1:ngridmx)
+            zdtsurf(1:ngridmx)               = zdtsurf(1:ngridmx) + zdtsurfmr(1:ngridmx)
+	    
+!	    print*,'after mass redistrib, q=',pq(211,1:nlayermx,igcm_h2o_vap)+ptimestep*pdq(211,1:nlayermx,igcm_h2o_vap) 
+	 endif
+
+
+
+!       ---------------------------------
+!       Updating tracer budget on surface
+!       ---------------------------------
+
+         if(hydrology)then
+
+            call hydrol(ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd,    &
+               capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice)
+            ! note: for now, also changes albedo in the subroutine
+
+            zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + zdtsurf_hyd(1:ngridmx)
+            qsurf(1:ngridmx,1:nqmx) = qsurf(1:ngridmx,1:nqmx) + ptimestep*dqs_hyd(1:ngridmx,1:nqmx)
+            ! when hydrology is used, other dqsurf tendencies are all added to dqs_hyd inside
+
+            !-------------------------
+            ! test energy conservation
+            if(enertest)then
+               dEtots = SUM(area(:)*capcal(:)*zdtsurf_hyd(:))/totarea
+               print*,'In hydrol surface energy change     =',dEtots,' W m-2'
+            endif
+            !-------------------------
+
+            !-------------------------
+            ! test water conservation
+            if(watertest)then
+	       dWtots =  SUM(dqs_hyd(:,igcm_h2o_ice)*area(:))*ptimestep/totarea
+               print*,'In hydrol surface ice change            =',dWtots,' kg m-2'
+	       dWtots =  SUM(dqs_hyd(:,igcm_h2o_vap)*area(:))*ptimestep/totarea
+               print*,'In hydrol surface water change          =',dWtots,' kg m-2'
+               print*,'---------------------------------------------------------------'
+            endif
+            !-------------------------
+
+         ELSE     ! of if (hydrology)
+
+            qsurf(1:ngridmx,1:nqmx)=qsurf(1:ngridmx,1:nqmx)+ptimestep*dqsurf(1:ngridmx,1:nqmx)
+
+         END IF   ! of if (hydrology)
+
+         ! Add qsurf to qsurf_hist, which is what we save in
+         ! diagfi.nc etc. At the same time, we set the water 
+         ! content of ocean gridpoints back to zero, in order
+         ! to avoid rounding errors in vdifc, rain
+         qsurf_hist(:,:) = qsurf(:,:)
+
+         if(ice_update)then
+            ice_min(1:ngridmx)=min(ice_min(1:ngridmx),qsurf(1:ngridmx,igcm_h2o_ice))
+         endif
+
+      endif   !  of if (tracer) 
+
+!-----------------------------------------------------------------------
+!   9. Surface and sub-surface soil temperature
+!-----------------------------------------------------------------------
+
+
+!     Increment surface temperature
+      tsurf(1:ngridmx)=tsurf(1:ngridmx)+ptimestep*zdtsurf(1:ngridmx) 
+
+!     Compute soil temperatures and subsurface heat flux
+      if (callsoil) then
+         call soil(ngrid,nsoilmx,.false.,inertiedat,        &
+                ptimestep,tsurf,tsoil,capcal,fluxgrd)
+      endif
+
+!-------------------------
+! test energy conservation
+      if(enertest)then
+         dEtots = SUM(area(:)*capcal(:)*zdtsurf(:))/totarea
+         print*,'Surface energy change                 =',dEtots,' W m-2'
+      endif
+!-------------------------
+
+!-----------------------------------------------------------------------
+!  10. Perform diagnostics and write output files
+!-----------------------------------------------------------------------
+
+!    -------------------------------
+!    Dynamical fields incrementation
+!    -------------------------------
+!    For output only: the actual model integration is performed in the dynamics
+ 
+      ! temperature, zonal and meridional wind
+      zt(1:ngridmx,1:nlayermx) = pt(1:ngridmx,1:nlayermx) + pdt(1:ngridmx,1:nlayermx)*ptimestep
+      zu(1:ngridmx,1:nlayermx) = pu(1:ngridmx,1:nlayermx) + pdu(1:ngridmx,1:nlayermx)*ptimestep
+      zv(1:ngridmx,1:nlayermx) = pv(1:ngridmx,1:nlayermx) + pdv(1:ngridmx,1:nlayermx)*ptimestep
+
+      ! diagnostic
+      zdtdyn(1:ngridmx,1:nlayermx)     = pt(1:ngridmx,1:nlayermx)-ztprevious(1:ngridmx,1:nlayermx)
+      ztprevious(1:ngridmx,1:nlayermx) = zt(1:ngridmx,1:nlayermx)
+
+      if(firstcall)then
+         zdtdyn(1:ngridmx,1:nlayermx)=0.0
+      endif
+
+      ! dynamical heating diagnostic
+      do ig=1,ngrid
+         fluxdyn(ig)= SUM(zdtdyn(ig,:) *mass(ig,:))*cpp/ptimestep
+      enddo
+
+      ! tracers
+      zq(1:ngridmx,1:nlayermx,1:nqmx) = pq(1:ngridmx,1:nlayermx,1:nqmx) + pdq(1:ngridmx,1:nlayermx,1:nqmx)*ptimestep
+
+      ! surface pressure
+      ps(1:ngridmx) = pplev(1:ngridmx,1) + pdpsrf(1:ngridmx)*ptimestep
+
+      ! pressure
+      do l=1,nlayer
+          zplev(1:ngridmx,l) = pplev(1:ngridmx,l)/pplev(1:ngridmx,1)*ps(:)
+          zplay(1:ngridmx,l) = pplay(1:ngridmx,l)/pplev(1:ngridmx,1)*ps(1:ngridmx)
+      enddo
+
+!     ---------------------------------------------------------
+!     Surface and soil temperature information
+!     ---------------------------------------------------------
+
+      Ts1 = SUM(area(:)*tsurf(:))/totarea
+      Ts2 = MINVAL(tsurf(:))
+      Ts3 = MAXVAL(tsurf(:))
+      if(callsoil)then
+         TsS = SUM(area(:)*tsoil(:,nsoilmx))/totarea        ! mean temperature at bottom soil layer
+         print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]     ave[Tdeep]'
+         print*,Ts1,Ts2,Ts3,TsS
+      else
+         print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]'
+         print*,Ts1,Ts2,Ts3
+      endif
+
+!     ---------------------------------------------------------
+!     Check the energy balance of the simulation during the run
+!     ---------------------------------------------------------
+
+      if(corrk)then
+
+         ISR = SUM(area(:)*fluxtop_dn(:))/totarea
+         ASR = SUM(area(:)*fluxabs_sw(:))/totarea
+         OLR = SUM(area(:)*fluxtop_lw(:))/totarea
+         GND = SUM(area(:)*fluxgrd(:))/totarea
+         DYN = SUM(area(:)*fluxdyn(:))/totarea
+         do ig=1,ngrid             
+            if(fluxtop_dn(ig).lt.0.0)then
+               print*,'fluxtop_dn has gone crazy'
+               print*,'fluxtop_dn=',fluxtop_dn(ig)
+               print*,'tau_col=',tau_col(ig)
+               print*,'aerosol=',aerosol(ig,:,:)
+               print*,'temp=   ',pt(ig,:)
+               print*,'pplay=  ',pplay(ig,:)
+               call abort
+            endif
+         end do
+                     
+         if(ngridmx.eq.1)then
+            DYN=0.0
+         endif
+
+         print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
+         print*, ISR,ASR,OLR,GND,DYN
+         
+         if(enertest)then
+            print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2'
+            print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2'
+            print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2'
+         endif
+
+         if(meanOLR)then
+            if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
+               ! to record global radiative balance
+               open(92,file="rad_bal.out",form='formatted',position='append')
+               write(92,*) zday,ISR,ASR,OLR
+               close(92)
+               open(93,file="tem_bal.out",form='formatted',position='append')
+               write(93,*) zday,Ts1,Ts2,Ts3,TsS
+               close(93)
+            endif
+         endif
+
+      endif
+
+
+!     ------------------------------------------------------------------
+!     Diagnostic to test radiative-convective timescales in code
+!     ------------------------------------------------------------------
+      if(testradtimes)then
+         open(38,file="tau_phys.out",form='formatted',position='append')
+         ig=1
+         do l=1,nlayer
+            write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l)
+         enddo
+         close(38)
+         print*,'As testradtimes enabled,'
+         print*,'exiting physics on first call'
+         call abort
+      endif
+
+!     ---------------------------------------------------------
+!     Compute column amounts (kg m-2) if tracers are enabled
+!     ---------------------------------------------------------
+      if(tracer)then
+         qcol(1:ngridmx,1:nqmx)=0.0
+         do iq=1,nq
+           do ig=1,ngridmx
+              qcol(ig,iq) = SUM( zq(ig,1:nlayermx,iq) * mass(ig,1:nlayermx))
+           enddo
+         enddo
+
+         ! Generalised for arbitrary aerosols now. (LK)
+         reffcol(1:ngridmx,1:naerkind)=0.0
+         if(co2cond.and.(iaero_co2.ne.0))then
+            call co2_reffrad(zq,reffrad)
+            do ig=1,ngridmx
+               reffcol(ig,iaero_co2) = SUM(zq(ig,1:nlayermx,igcm_co2_ice)*reffrad(ig,1:nlayermx,iaero_co2)*mass(ig,1:nlayermx))
+            enddo
+         endif
+         if(water.and.(iaero_h2o.ne.0))then
+            call h2o_reffrad(zq,zt,reffrad,nueffrad)
+            do ig=1,ngridmx
+               reffcol(ig,iaero_h2o) = SUM(zq(ig,1:nlayermx,igcm_h2o_ice)*reffrad(ig,1:nlayermx,iaero_h2o)*mass(ig,1:nlayermx))
+            enddo
+         endif
+
+      endif
+
+!     ---------------------------------------------------------
+!     Test for water conservation if water is enabled
+!     ---------------------------------------------------------
+
+      if(water)then
+
+         icesrf = SUM(area(:)*qsurf_hist(:,igcm_h2o_ice))/totarea
+         liqsrf = SUM(area(:)*qsurf_hist(:,igcm_h2o_vap))/totarea
+         icecol = SUM(area(:)*qcol(:,igcm_h2o_ice))/totarea
+         vapcol = SUM(area(:)*qcol(:,igcm_h2o_vap))/totarea
+
+         h2otot = icesrf + liqsrf + icecol + vapcol
+
+         print*,' Total water amount [kg m^-2]: ',h2otot
+         print*,' Surface ice    Surface liq.   Atmos. con.     Atmos. vap. [kg m^-2] '
+         print*, icesrf,liqsrf,icecol,vapcol
+
+         if(meanOLR)then
+            if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
+               ! to record global water balance
+               open(98,file="h2o_bal.out",form='formatted',position='append')
+               write(98,*) zday,icesrf,liqsrf,icecol,vapcol
+               close(98)
+            endif
+         endif
+
+      endif
+
+!     ---------------------------------------------------------
+!     Calculate RH for diagnostic if water is enabled
+!     ---------------------------------------------------------
+
+      if(water)then
+         do l = 1, nlayer
+            do ig = 1, ngrid
+!               call watersat(pt(ig,l),pplay(ig,l),qsat(ig,l))
+               call Psat_water(zt(ig,l),pplay(ig,l),psat_tmp,qsat(ig,l))
+
+               RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l)
+            enddo
+         enddo
+
+         ! compute maximum possible H2O column amount (100% saturation)
+         do ig=1,ngrid
+               H2Omaxcol(ig) = SUM( qsat(ig,:) * mass(ig,:))
+         enddo
+
+      endif
+
+
+         print*,''
+         print*,'--> Ls =',zls*180./pi
+!        -------------------------------------------------------------------
+!        Writing NetCDF file  "RESTARTFI" at the end of the run
+!        -------------------------------------------------------------------
+!        Note: 'restartfi' is stored just before dynamics are stored
+!              in 'restart'. Between now and the writting of 'restart',
+!              there will have been the itau=itau+1 instruction and
+!              a reset of 'time' (lastacll = .true. when itau+1= itaufin)
+!              thus we store for time=time+dtvr
+
+         if(lastcall) then
+            ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
+
+
+!           Update surface ice distribution to iterate to steady state if requested
+            if(ice_update)then
+
+               do ig = 1, ngrid
+
+                  delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
+
+                  ! add multiple years of evolution
+                  qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep 
+
+                  ! if ice has gone -ve, set to zero
+                  if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
+                     qsurf_hist(ig,igcm_h2o_ice) = 0.0 
+                  endif
+
+                  ! if ice is seasonal, set to zero (NEW)
+                  if(ice_min(ig).lt.0.01)then
+                     qsurf_hist(ig,igcm_h2o_ice) = 0.0 
+                  endif
+
+               enddo
+
+               ! enforce ice conservation
+               ice_tot= SUM(qsurf_hist(:,igcm_h2o_ice)*area(:) )
+               qsurf_hist(:,igcm_h2o_ice) = qsurf_hist(:,igcm_h2o_ice)*(icesrf/ice_tot)
+
+            endif
+
+            write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
+            call physdem1("restartfi.nc",long,lati,nsoilmx,nq,            &
+                    ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
+                    area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,      &
+                    cloudfrac,totcloudfrac,hice)
+         endif
+
+!        -----------------------------------------------------------------
+!        Saving statistics :
+!        -----------------------------------------------------------------
+!        ("stats" stores and accumulates 8 key variables in file "stats.nc"
+!        which can later be used to make the statistic files of the run:
+!        "stats")          only possible in 3D runs !
+
+         
+         if (callstats) then
+
+            call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
+            call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
+            call wstats(ngrid,"fluxsurf_lw",                               &
+                        "Thermal IR radiative flux to surface","W.m-2",2,  &
+                         fluxsurf_lw)
+!            call wstats(ngrid,"fluxsurf_sw",                               &
+!                        "Solar radiative flux to surface","W.m-2",2,       &
+!                         fluxsurf_sw_tot)
+            call wstats(ngrid,"fluxtop_lw",                                &
+                        "Thermal IR radiative flux to space","W.m-2",2,    &
+                        fluxtop_lw)
+!            call wstats(ngrid,"fluxtop_sw",                                &
+!                        "Solar radiative flux to space","W.m-2",2,         &
+!                        fluxtop_sw_tot)
+
+            call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
+            call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
+            call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
+
+            call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
+            call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
+            call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
+            call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
+            call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
+
+           if (tracer) then
+             do iq=1,nq
+                call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
+                call wstats(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
+                          'kg m^-2',2,qsurf(1,iq) )
+
+                call wstats(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    & 
+                          'kg m^-2',2,qcol(1,iq) )
+!                call wstats(ngridmx,trim(noms(iq))//'_reff',                          & 
+!                          trim(noms(iq))//'_reff',                                   & 
+!                          'm',3,reffrad(1,1,iq))
+              end do
+            if (water) then
+               vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
+               call wstats(ngrid,"vmr_h2ovapor",                           &
+                          "H2O vapour volume mixing ratio","mol/mol",       & 
+                          3,vmr)
+            endif ! of if (water)
+
+           endif !tracer 
+
+            if(lastcall) then
+              write (*,*) "Writing stats..."
+              call mkstats(ierr)
+            endif
+          endif !if callstats
+
+
+!        ----------------------------------------------------------------------
+!        output in netcdf file "DIAGFI", containing any variable for diagnostic
+!        (output with  period "ecritphy", set in "run.def")
+!        ----------------------------------------------------------------------
+!        writediagfi can also be called from any other subroutine for any variable.
+!        but its preferable to keep all the calls in one place...
+
+        call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
+        call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
+        call writediagfi(ngrid,"temp","temperature","K",3,zt)
+        call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
+        call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
+        call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
+        call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
+        call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
+
+!     Total energy balance diagnostics
+        if(callrad.and.(.not.newtonian))then
+           call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo)
+           call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
+           call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
+           call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
+           call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
+           call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
+        endif
+	
+	if(enertest) then
+	  if (calldifv) then
+	     call writediagfi(ngridmx,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
+!	     call writediagfi(ngridmx,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
+!	     call writediagfi(ngridmx,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
+!	     call writediagfi(ngridmx,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
+	     call writediagfi(ngridmx,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
+	  endif
+	  if (corrk) then
+!	     call writediagfi(ngridmx,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
+!	     call writediagfi(ngridmx,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
+	  endif 
+          if(watercond) then
+	     call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE) 
+	     call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
+	     call writediagfi(ngridmx,"RH","relative humidity"," ",3,RH)
+!	     call writediagfi(ngridmx,"qsatatm","atm qsat"," ",3,qsat)
+!	     call writediagfi(ngridmx,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
+          endif 
+        endif
+
+!     Temporary inclusions for heating diagnostics
+!        call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
+!        call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
+!        call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
+!        call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
+
+        ! debugging
+        !call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
+        !call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
+
+!     Output aerosols
+        if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
+	  call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
+        if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
+	  call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
+        if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
+	  call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
+        if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
+	  call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
+
+!     Output tracers
+       if (tracer) then
+        do iq=1,nq
+          call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
+!          call writediagfi(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
+!                          'kg m^-2',2,qsurf(1,iq) )
+          call writediagfi(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
+                          'kg m^-2',2,qsurf_hist(1,iq) )
+          call writediagfi(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    & 
+                          'kg m^-2',2,qcol(1,iq) )
+
+          if(watercond.or.CLFvarying)then
+!             call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
+!             call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
+!             call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
+             call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
+          endif
+
+          if(waterrain)then
+             call writediagfi(ngridmx,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
+             call writediagfi(ngridmx,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
+          endif
+
+          if(hydrology)then
+             call writediagfi(ngridmx,"hice","oceanic ice height","m",2,hice)
+          endif
+
+          if(ice_update)then
+             call writediagfi(ngridmx,"ice_min","min annual ice","m",2,ice_min)
+             call writediagfi(ngridmx,"ice_ini","initial annual ice","m",2,ice_initial)
+          endif
+
+          do ig=1,ngrid
+             if(tau_col(ig).gt.1.e3)then
+                print*,'WARNING: tau_col=',tau_col(ig)
+                print*,'at ig=',ig,'in PHYSIQ'
+             endif
+          end do
+
+          call writediagfi(ngridmx,"tau_col","Total aerosol optical depth","[]",2,tau_col)
+
+        enddo
+       endif
+
+!      output spectrum 
+      if(specOLR.and.corrk)then      
+         call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
+         call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
+      endif
+
+
+      icount=icount+1
+
+      ! deallocate gas variables
+      if (lastcall) then
+        IF ( ALLOCATED( gnom ) ) DEALLOCATE( gnom )
+        IF ( ALLOCATED( gfrac ) ) DEALLOCATE( gfrac )  ! both allocated in su_gases.F90
+      endif
+
+
+      return
+    end subroutine physiqg