subroutine physiq(ngrid,nlayer,nq,   &
                  firstcall,lastcall,      &
                  pday,ptime,ptimestep,    &
                  pplev,pplay,pphi,        &
                  pu,pv,pt,pq,             &
                  pw,                      &
                  pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
 
      use radinc_h, only : naerkind,L_NSPECTI,L_NSPECTV
      use watercommon_h, only : mx_eau_sol, To, RLVTT, mH2O
      use gases_h
      use radcommon_h, only: sigma
      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)
!     
!==================================================================


!    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,i, tapphys,nw

      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)

      real,save :: reffrad(ngridmx,nlayermx,naerkind) ! aerosol effective radius (m)

!     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 zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx)    ! (m.s-2)
      real zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
      real zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)

      real zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
      real zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
      real zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
      real zdqadj(ngridmx,nlayermx,nqmx)
      real zdqc(ngridmx,nlayermx,nqmx)
      real zdqlscale(ngridmx,nlayermx,nqmx)
      real zdqslscale(ngridmx,nqmx)
      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 dEtot, dEtots, masse, AtmToSurf_TurbFlux
      real dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
      real dEzRadsw(ngridmx,nlayermx),dEzRadlw(ngridmx,nlayermx),dEzdif(ngridmx,nlayermx)
      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 dqevap(ngridmx,nlayermx)
      real dtevap(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)

!     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 cloudfrac(ngridmx,nlayermx)
      real totcloudfrac(ngridmx)

!     included by RW for vdifc water conservation test
      real nconsMAX
      real vdifcncons(ngridmx)
      real cadjncons(ngridmx)
      real qzero1D
      save qzero1D

!      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 reffH2O(ngridmx,nlayermx)
      real reffcol(ngridmx,2)

!     included by RW for sourceevol
      real ice_initial(ngridmx)!, isoil
      real delta_ice,ice_tot
      real ice_min(ngridmx)
      save ice_min
      save ice_initial

      integer num_run
      logical ice_update
      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 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)

         do ig=1,ngrid
            albedo(ig)=albedo0(ig)
         enddo

         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'
            do ig=1,ngrid
               capcal(ig)=1.e6
               fluxgrd(ig)=0.
               if(noradsurf)then
                  fluxgrd(ig)=10.0
               endif
            enddo
            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
!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         do ig=1,ngridmx
            rnat(ig)=1

!            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 
!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         do ig=1,ngridmx
            do iq=1,nqmx
               qsurf_hist(ig,iq)=qsurf(ig,iq)
            enddo
         enddo

!        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

            qzero1D=0.0

            if(enertest)then
               watertest = .true.
            endif

            if(ngrid.eq.1)then
               qzero1D               = 0.0
               qsurf(1,igcm_h2o_vap) = qzero1D
               do l=1, nlayer
                  pq(1,l,igcm_h2o_vap)=0.0
                  pq(1,l,igcm_h2o_ice)=0.0
               enddo
            endif

            do ig=1,ngrid
               qsurf_hist(ig,igcm_h2o_vap) = qsurf(ig,igcm_h2o_vap)
               if(ice_update)then
                  ice_initial(ig)=qsurf(ig,igcm_h2o_ice)
               endif
            enddo

         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(:,:) = 0.0
      pdv(:,:) = 0.0
!      if ( (.not.nearco2cond).and.(.not.firstcall) ) then
      if ( .not.nearco2cond ) then
         pdt(:,:) = 0.0
      endif ! this was the source of an evil bug...
      pdq(:,:,:)  = 0.0
      pdpsrf(:)   = 0.0
      zflubid(:)  = 0.0
      zdtsurf(:)  = 0.0
      dqsurf(:,:) = 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
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      do l=1,nlayer
         do ig=1,ngrid
            zzlay(ig,l)=pphi(ig,l)/g
         enddo
      enddo
      do ig=1,ngrid
         zzlev(ig,1)=0.
         zzlev(ig,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
      enddo
      do l=2,nlayer
         do ig=1,ngrid
            z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
            z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
            zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
         enddo
      enddo
!     Potential temperature calculation 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)
         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)
                    
                    do l=1,nlayer
                       zdtlw(ig,l) = ntf*zdtlw1(ig,l) + tf*zdtlw(ig,l)
                       zdtsw(ig,l) = ntf*zdtsw1(ig,l) + tf*zdtsw(ig,l)
                    enddo

                    do nw=1,L_NSPECTV
                       OSR_nu(ig,nw) = ntf*OSR_nu1(ig,nw) + tf*OSR_nu(ig,nw)		       
                    enddo
                    do nw=1,L_NSPECTI
                       OLR_nu(ig,nw) = ntf*OLR_nu1(ig,nw) + tf*OLR_nu(ig,nw)		       
                    enddo

                 enddo

              endif !CLFvarying

!             Radiative flux from the sky absorbed by the surface (W.m-2)
              GSR=0.0
              do ig=1,ngrid
                 fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)             &
                      +fluxsurf_sw(ig)*(1.-albedo(ig))

                 if(noradsurf)then ! no lower surface; SW flux just disappears
                    GSR = GSR + fluxsurf_sw(ig)*area(ig)
                    fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
                 endif

              enddo
              if(noradsurf)then
                 print*,'SW lost in deep atmosphere = ',GSR/totarea,' W m^-2'
              endif

!             Net atmospheric radiative heating rate (K.s-1)
              do l=1,nlayer
                 do ig=1,ngrid
                    dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
                 enddo
              enddo


           elseif(newtonian)then

!          b) Call Newtonian cooling scheme
!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
              call newtrelax(mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)

              do ig=1,ngrid
                 zdtsurf(ig) = +(pt(ig,1)-tsurf(ig))/ptimestep
              enddo
              ! e.g. surface becomes proxy for 1st atmospheric layer ?

           else

!          c) Atmosphere has no radiative effect
!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
              do ig=1,ngrid
                 fluxtop_dn(ig)  = fract(ig)*mu0(ig)*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(ig) = fluxtop_dn(ig)
                 fluxrad_sky(ig) = fluxtop_dn(ig)*(1.-albedo(ig))
                 fluxtop_lw(ig)  = emis(ig)*sigma*tsurf(ig)**4
              enddo ! radiation skips the atmosphere entirely

              do l=1,nlayer
                 do ig=1,ngrid
                    dtrad(ig,l)=0.0
                 enddo
              enddo ! hence no atmospheric radiative heating

           endif                ! if corrk

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


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

        do ig=1,ngrid
           zplanck(ig)=tsurf(ig)*tsurf(ig)
           zplanck(ig)=emis(ig)*sigma*zplanck(ig)*zplanck(ig)
           fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
        enddo

        do l=1,nlayer
           do ig=1,ngrid
              pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
           enddo
        enddo

!-------------------------
! test energy conservation
         if(enertest)then
            dEtotSW  = 0.0
            dEtotLW  = 0.0
            dEtotsSW = 0.0
            dEtotsLW = 0.0
	    dEzRadsw(:,:)=0.
	    dEzRadlw(:,:)=0.
            do ig = 1, ngrid
               do l = 1, nlayer
                  masse  = (pplev(ig,l) - pplev(ig,l+1))/g
                  dEtotSW = dEtotSW + cpp*masse*zdtsw(ig,l)*area(ig)
                  dEtotLW = dEtotLW + cpp*masse*zdtlw(ig,l)*area(ig)
		  dEzRadsw(ig,l)=cpp*masse*zdtsw(ig,l)
		  dEzRadlw(ig,l)=cpp*masse*zdtlw(ig,l)
               enddo
               dEtotsSW = dEtotsSW + fluxsurf_sw(ig)*(1.-albedo(ig))*area(ig)
               dEtotsLW = dEtotsLW + emis(ig)*fluxsurf_lw(ig)*area(ig)-zplanck(ig)*area(ig)
            enddo
            dEtotSW  = dEtotSW/totarea
            dEtotLW  = dEtotLW/totarea
            dEtotsSW = dEtotsSW/totarea
            dEtotsLW = dEtotsLW/totarea
            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
      
         do ig=1,ngrid
            zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
         enddo

         zdum1(:,:)=0.0
         zdum2(:,:)=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,zdqsdif,lastcall)

	 else
	 
           do l=1,nlayer
             do ig=1,ngrid 
               zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
             enddo
           enddo
 
           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)

           do l=1,nlayer
             do ig=1,ngrid
                zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
             enddo
           enddo

         end if !turbdiff

         do l=1,nlayer
            do ig=1,ngrid
               pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
               pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
               pdt(ig,l)=pdt(ig,l)+zdtdif(ig,l)
            enddo
         enddo

         do ig=1,ngrid
            zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
         enddo

         if (tracer) then 
           do iq=1, nq
            do l=1,nlayer
              do ig=1,ngrid
                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq) 
              enddo
            enddo
           enddo
           do iq=1, nq
              do ig=1,ngrid
                 dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
              enddo
           enddo

         end if ! of if (tracer)

         !-------------------------
         ! test energy conservation
         if(enertest)then
            dEtot=0.0
            dEtots=0.0
	    dEzdif(:,:)=0.
            AtmToSurf_TurbFlux=0.
	    
            do ig = 1, ngrid
               do l = 1, nlayer
                  masse = (pplev(ig,l) - pplev(ig,l+1))/g
                  dEtot = dEtot + cpp*masse*zdtdif(ig,l)*area(ig)
		  dEzdif (ig,l)= cpp*masse*(zdtdif(ig,l))
               enddo
	       dEzdif(ig,1)= dEzdif(ig,1)+ sensibFlux(ig)! subtract flux to the ground
               dEtot = dEtot + sensibFlux(ig)*area(ig)! subtract flux to the ground
               dEtots = dEtots + capcal(ig)*zdtsdif(ig)*area(ig)-zflubid(ig)*area(ig)-sensibFlux(ig)*area(ig)
	       AtmToSurf_TurbFlux=AtmToSurf_TurbFlux+sensibFlux(ig)*area(ig)
            enddo
            dEtot  = dEtot/totarea
            dEtots = dEtots/totarea
	    AtmToSurf_TurbFlux=AtmToSurf_TurbFlux/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=0.0
            dWtots=0.0
            nconsMAX=0.0
            do ig = 1, ngrid

               vdifcncons(ig)=0.0
               do l = 1, nlayer
                  masse = (pplev(ig,l) - pplev(ig,l+1))/g

                  iq    = igcm_h2o_vap
                  dWtot = dWtot + masse*zdqdif(ig,l,iq)*ptimestep*area(ig)
                  vdifcncons(ig)=vdifcncons(ig) + masse*zdqdif(ig,l,iq)

                  iq    = igcm_h2o_ice
                  dWtot = dWtot + masse*zdqdif(ig,l,iq)*ptimestep*area(ig)
                  vdifcncons(ig)=vdifcncons(ig) + masse*zdqdif(ig,l,iq)
               
               enddo

               iq     = igcm_h2o_vap
               dWtots = dWtots + zdqsdif(ig,iq)*ptimestep*area(ig)
               vdifcncons(ig)=vdifcncons(ig)+zdqsdif(ig,iq)

               iq     = igcm_h2o_ice
               dWtots = dWtots + zdqsdif(ig,iq)*ptimestep*area(ig)
               vdifcncons(ig)=vdifcncons(ig)+zdqsdif(ig,iq)

               if(vdifcncons(ig).gt.nconsMAX)then
                  nconsMAX=vdifcncons(ig)
               endif

            enddo

            dWtot  = dWtot/totarea
            dWtots = dWtots/totarea
            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

            do ig=1,ngrid
               zdtsurf(ig) = zdtsurf(ig) + (fluxrad(ig) + fluxgrd(ig))/capcal(ig)
            enddo

         endif

      endif ! of if (calldifv)


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

      if(calladj) then

         do l=1,nlayer
            do ig=1,ngrid
               zdh(ig,l) = pdt(ig,l)/zpopsk(ig,l)
            enddo
         enddo
         zduadj(:,:)=0.0
         zdvadj(:,:)=0.0
         zdhadj(:,:)=0.0


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

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

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

         !-------------------------
         ! test energy conservation
         if(enertest)then
            dEtot=0.0
            do ig = 1, ngrid
               do l = 1, nlayer
                  masse = (pplev(ig,l) - pplev(ig,l+1))/g
                  dEtot = dEtot + cpp*masse*zdtadj(ig,l)*area(ig)
               enddo
            enddo
            dEtot=dEtot/totarea
          print*,'In convadj atmospheric energy change  =',dEtot,' W m-2'
         endif
         !-------------------------

         !-------------------------
         ! test water conservation
         if(watertest)then
            dWtot=0.0
            do ig = 1, ngrid
               cadjncons(ig)=0.0
               do l = 1, nlayer
                  masse = (pplev(ig,l) - pplev(ig,l+1))/g
                  
                  iq    = igcm_h2o_vap
                  dWtot = dWtot + masse*zdqadj(ig,l,iq)*ptimestep*area(ig)
                  cadjncons(ig)=cadjncons(ig) + masse*zdqadj(ig,l,iq)
                  
                  iq    = igcm_h2o_ice
                  dWtot = dWtot + masse*zdqadj(ig,l,iq)*ptimestep*area(ig)
                  cadjncons(ig)=cadjncons(ig) + masse*zdqadj(ig,l,iq)
                  
               enddo

            enddo
            dWtot=dWtot/totarea
            print*,'In convadj atmospheric water change     =',dWtot,' kg m-2'
         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

         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)

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

         do iq=1,nq ! should use new notation here !
            do l=1,nlayer
               do ig=1,ngrid
                  pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq) 
               enddo
            enddo
         enddo
         ! 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=0.0
            dEtots=0.0
            do ig = 1, ngrid
               do l = 1, nlayer
                  masse = (pplev(ig,l) - pplev(ig,l+1))/g
                  dEtot = dEtot + cpp*masse*zdtc(ig,l)*area(ig)
               enddo
               dEtots = dEtots + capcal(ig)*zdtsurfc(ig)*area(ig)
            enddo
            dEtot=dEtot/totarea
            dEtots=dEtots/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)then

	       if(RLVTT.gt.1.e-8)then

               ! Re-evaporate cloud water/ice
               call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
               DO l = 1, nlayer  
                  DO ig = 1, ngrid 
                     pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqevap(ig,l)
                     pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)-dqevap(ig,l)
                     pdt(ig,l)              = pdt(ig,l)+dtevap(ig,l)
                  enddo
               enddo

               call moistadj(pt,qevap,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
               do l=1,nlayer
                  do ig=1,ngrid
                     pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqmoist(ig,l,igcm_h2o_vap)
                     pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+dqmoist(ig,l,igcm_h2o_ice)
                     pdt(ig,l)              = pdt(ig,l)+dtmoist(ig,l)
                  enddo
               enddo

               !-------------------------
               ! test energy conservation
               if(enertest)then
                  dEtot=0.0
                  madjdE(:)=0.0
                  do ig = 1, ngrid
                     do l = 1, nlayer
                        masse = (pplev(ig,l) - pplev(ig,l+1))/g
                        dEtot = dEtot + cpp*masse*(dtmoist(ig,l)+dtevap(ig,l))*area(ig)
                        madjdE(ig) = madjdE(ig) + cpp*masse*(dtmoist(ig,l)+dtevap(ig,l))
                     enddo
                  enddo
                  dEtot=dEtot/totarea
                  print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
               endif
               !-------------------------

               !-------------------------
               ! test water conservation
               if(watertest)then
                  dWtot=0.0
                  do ig = 1, ngrid
                     do iq = 1 , nq
                        do l = 1, nlayer
                           masse = (pplev(ig,l) - pplev(ig,l+1))/g
                           dWtot = dWtot + masse*dqmoist(ig,l,igcm_h2o_vap)*area(ig)*ptimestep
                           dWtot = dWtot + masse*dqmoist(ig,l,igcm_h2o_ice)*area(ig)*ptimestep
                        enddo
                     enddo
                  enddo
                  dWtot=dWtot/totarea
                  print*,'In moistadj atmospheric water change    =',dWtot,' kg m-2'
               endif
               !-------------------------
               

	       endif
	       

               ! Re-evaporate cloud water/ice
               call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
               do l = 1, nlayer  
                  do ig = 1, ngrid
                     pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqevap(ig,l)
                     pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)-dqevap(ig,l)
                     pdt(ig,l)              = pdt(ig,l)+dtevap(ig,l)
                  enddo
               enddo ! note: we use qevap but not tevap in largescale/moistadj
                     ! otherwise is a big mess

               call largescale(ptimestep,pplev,pplay,pt,qevap,        & ! a bug was here!
                    pdt,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc,reffH2O)
               do l=1,nlayer
                  do ig=1,ngrid
                     pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+dqvaplscale(ig,l)
                     pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+dqcldlscale(ig,l)
                     pdt(ig,l)              = pdt(ig,l)+dtlscale(ig,l)

                     if(.not.aerofixed)then
                        reffrad(ig,l,2)=reffH2O(ig,l)
                     endif

                  enddo
               enddo

               !-------------------------
               ! test energy conservation
               if(enertest)then
                  dEtot=0.0
                  lscaledE(:)=0.0
                  do ig = 1, ngrid
                     do l = 1, nlayer
                        masse = (pplev(ig,l) - pplev(ig,l+1))/g
                        dEtot = dEtot + cpp*masse*(dtlscale(ig,l)+dtevap(ig,l))*area(ig)
                        lscaledE(ig) = lscaledE(ig) + cpp*masse*(dtlscale(ig,l)+dtevap(ig,l))
                     enddo
                  enddo
                  dEtot=dEtot/totarea
                  print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
               endif
               !-------------------------

               !-------------------------
               ! test water conservation
               if(watertest)then
                  dWtot=0.0
                  do ig = 1, ngrid
                     do iq = 1 , nq
                        do l = 1, nlayer
                           masse = (pplev(ig,l) - pplev(ig,l+1))/g
                           dWtot = dWtot + masse*dqvaplscale(ig,l)*area(ig)*ptimestep
                           dWtot = dWtot + masse*dqcldlscale(ig,l)*area(ig)*ptimestep
                        enddo
                     enddo
                  enddo
                  dWtot=dWtot/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

               ! compute total cloud fraction in column
               call totalcloudfrac(cloudfrac,totcloudfrac)

           endif                ! of if (watercondense)
           

!        --------------------------------
!        Water ice / liquid precipitation
!        --------------------------------
           if(waterrain)then

              zdqrain(:,:,:) = 0.0
              zdqsrain(:)    = 0.0
              zdqssnow(:)    = 0.0

              call rain(ptimestep,pplev,pplay,pt,pdt,pq,pdq,            &
                   zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)

              do l=1,nlayer
                 do ig=1,ngrid
                    pdq(ig,l,igcm_h2o_vap) = pdq(ig,l,igcm_h2o_vap)+zdqrain(ig,l,igcm_h2o_vap)
                    pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+zdqrain(ig,l,igcm_h2o_ice)
                    pdt(ig,l)              = pdt(ig,l)+zdtrain(ig,l)
                 enddo
              enddo

              do ig=1,ngrid
                 dqsurf(ig,igcm_h2o_vap) = dqsurf(ig,igcm_h2o_vap)+zdqsrain(ig) ! a bug was here
                 dqsurf(ig,igcm_h2o_ice) = dqsurf(ig,igcm_h2o_ice)+zdqssnow(ig)
                 rainout(ig)             = zdqsrain(ig)+zdqssnow(ig) ! diagnostic   
              enddo

              !-------------------------
              ! test energy conservation
              if(enertest)then
                 dEtot=0.0
                 do ig = 1, ngrid
                    do l = 1, nlayer
                       masse = (pplev(ig,l) - pplev(ig,l+1))/g
                       dEtot = dEtot + cpp*masse*zdtrain(ig,l)*area(ig)
                    enddo
                 enddo
                 dEtot=dEtot/totarea
                 print*,'In rain atmospheric energy change       =',dEtot,' W m-2'
              endif
              !-------------------------


              !-------------------------
              ! test energy conservation
              if(enertest)then
                 dEtot=0.0
                 do ig = 1, ngrid
                    do l = 1, nlayer
                       masse = (pplev(ig,l) - pplev(ig,l+1))/g
                       dEtot = dEtot + cpp*masse*zdtrain(ig,l)*area(ig)
                    enddo
                 enddo
                 dEtot=dEtot/totarea
                 print*,'In rain atmospheric T energy change       =',dEtot,' W m-2'

                 dEtot=0.0
                 do ig = 1, ngrid
                    do l = 1, nlayer
                       masse = (pplev(ig,l) - pplev(ig,l+1))/g
                       dItot = dItot + masse*zdqrain(ig,l,igcm_h2o_ice)*area(ig)
                       dVtot = dVtot + masse*zdqrain(ig,l,igcm_h2o_vap)*area(ig)
                    enddo
                    dItot = dItot + zdqssnow(ig)*area(ig)
                    dVtot = dVtot + zdqsrain(ig)*area(ig)
                 enddo
                 dEtot=(dItot*RLVTT/cpp + dVtot*RLVTT/cpp)/totarea
                 print*,'In rain dItot =',dItot*RLVTT/(cpp*totarea),' W m-2'
                 print*,'In rain dVtot =',dVtot*RLVTT/(cpp*totarea),' W m-2'
                 print*,'In rain atmospheric L energy change       =',dEtot,' W m-2'
              endif
              !-------------------------

              !-------------------------
              ! test water conservation
              if(watertest)then
                 dWtot=0.0
                 dWtots=0.0
                 do ig = 1, ngrid
                    !do iq = 1 , nq
                       do l = 1, nlayer
                          masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
                          dWtot = dWtot + masse*zdqrain(ig,l,igcm_h2o_vap)*area(ig)*ptimestep
                          dWtot = dWtot + masse*zdqrain(ig,l,igcm_h2o_ice)*area(ig)*ptimestep
                       enddo
                    !enddo
                    dWtots = dWtots + (zdqsrain(ig)+zdqssnow(ig))*area(ig)*ptimestep
                 enddo
                 dWtot=dWtot/totarea
                 dWtots=dWtots/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(:,:,:) = 0.0
           zdqssed(:,:)  = 0.0


           !-------------------------
           ! find qtot
           if(watertest)then
              dWtot=0.0
              dWtots=0.0
              iq=3
              do ig = 1, ngrid
                 do l = 1, nlayer
                    masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
                    dWtot  = dWtot  + masse*pq(ig,l,iq)*area(ig)*ptimestep
                    dWtots = dWtots + masse*pdq(ig,l,iq)*area(ig)*ptimestep
                 enddo
              enddo
              dWtot=dWtot/totarea
              dWtots=dWtots/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
              dWtot=0.0
              dWtots=0.0
              iq=3
              do ig = 1, ngrid
                 do l = 1, nlayer
                    masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
                    dWtot  = dWtot  + masse*pq(ig,l,iq)*area(ig)*ptimestep
                    dWtots = dWtots + masse*pdq(ig,l,iq)*area(ig)*ptimestep
                 enddo
              enddo
              dWtot=dWtot/totarea
              dWtots=dWtots/totarea
              print*,'After sedim pq  =',dWtot,' kg m-2'
              print*,'After sedim pdq =',dWtots,' kg m-2'
           endif
           !-------------------------

           do iq=1,nq
              ! 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
              do ig=1,ngrid
                 do l=1,nlayer
                    pdq(ig,l,iq) = pdq(ig,l,iq) + zdqsed(ig,l,iq)
                 enddo
                 dqsurf(ig,iq) = dqsurf(ig,iq) + zdqssed(ig,iq)
              enddo
           enddo

           !-------------------------
           ! test water conservation
           if(watertest)then
              dWtot=0.0
              dWtots=0.0
              do iq=1,nq
              do ig = 1, ngrid
                 do l = 1, nlayer
                    masse = (pplev(ig,l) - pplev(ig,l+1))/g ! equiv to l2c in rain
                    dWtot = dWtot + masse*zdqsed(ig,l,iq)*area(ig)*ptimestep
                 enddo
                 dWtots   = dWtots + zdqssed(ig,iq)*area(ig)*ptimestep
              enddo
              enddo
              dWtot=dWtot/totarea
              dWtots=dWtots/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 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

            do ig=1,ngrid
               zdtsurf(ig) = zdtsurf(ig) + zdtsurf_hyd(ig)
               do iq=1,nq
                  qsurf(ig,iq) = qsurf(ig,iq)+ptimestep*dqs_hyd(ig,iq)
               enddo
            enddo
            ! when hydrology is used, other dqsurf tendencies are all added to dqs_hyd inside

            !-------------------------
            ! test energy conservation
            if(enertest)then
               dEtot=0.0
               do ig = 1, ngrid
                  dEtots = dEtots + capcal(ig)*zdtsurf_hyd(ig)*area(ig)
               enddo
               dEtot=dEtot/totarea
               print*,'In hydrol atmospheric energy change     =',dEtot,' W m-2'
            endif
            !-------------------------

            !-------------------------
            ! test water conservation
            if(watertest)then
               dWtots=0.0
               do ig = 1, ngrid
                  dWtots = dWtots + dqs_hyd(ig,igcm_h2o_ice)*area(ig)*ptimestep
               enddo
               dWtots=dWtots/totarea
               print*,'In hydrol surface ice change            =',dWtots,' kg m-2'
               dWtots=0.0
               do ig = 1, ngrid
                  dWtots = dWtots + dqs_hyd(ig,igcm_h2o_vap)*area(ig)*ptimestep
               enddo
               dWtots=dWtots/totarea
               print*,'In hydrol surface water change          =',dWtots,' kg m-2'
               print*,'---------------------------------------------------------------'
            endif
            !-------------------------

         ELSE     ! of if (hydrology)

            do iq=1,nq
               do ig=1,ngrid
                  qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
               enddo
            enddo           

         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
         do ig = 1, ngrid
            do iq = 1, nq
               if(iq.eq.igcm_h2o_vap .and. rnat(ig).eq.0)then ! if liquid water and terrain = ocean 
                  qsurf_hist(ig,iq) = qsurf(ig,iq)
                  !qsurf(ig,iq)      = qcol(ig,iq) 
                  ! the value of qsurf we choose here makes NO DIFFERENCE TO ANYTHING AT ALL
               else
                  qsurf_hist(ig,iq) = qsurf(ig,iq)
               endif
            enddo
         enddo

         if(ice_update)then
            do ig = 1, ngrid
               ice_min(ig)=min(ice_min(ig),qsurf(ig,igcm_h2o_ice))
            enddo
         endif

      endif   !  of if (tracer) 

!-----------------------------------------------------------------------
!   9. Surface and sub-surface soil temperature
!-----------------------------------------------------------------------


!     Increment surface temperature
      do ig=1,ngrid
         tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) 
      enddo

!     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=0.0
         do ig = 1, ngrid
            dEtots = dEtots + capcal(ig)*zdtsurf(ig)*area(ig)
         enddo
         dEtots=dEtots/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
      do l=1,nlayer
        do ig=1,ngrid
          zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
          zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
          zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep

          ! diagnostic
          zdtdyn(ig,l)     = ztprevious(ig,l)-pt(ig,l)
          ztprevious(ig,l) = zt(ig,l)
        enddo
      enddo

      if(firstcall)then
         zdtdyn(:,:)=0.0
      endif

      ! dynamical heating diagnostic
      fluxdyn(:)=0.
      do ig=1,ngrid
         do l=1,nlayer
            fluxdyn(ig)=fluxdyn(ig) - (zdtdyn(ig,l)/ptimestep)   &
                 *(pplev(ig,l)-pplev(ig,l+1))*cpp/g
         enddo
      enddo

      ! tracers
      do iq=1, nq
        do l=1,nlayer
          do ig=1,ngrid
            zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
          enddo
        enddo
      enddo

      ! surface pressure
      do ig=1,ngrid
          ps(ig) = pplev(ig,1) + pdpsrf(ig)*ptimestep
      enddo

      ! pressure
      do l=1,nlayer
        do ig=1,ngrid
             zplev(ig,l) = pplev(ig,l)/pplev(ig,1)*ps(ig)
             zplay(ig,l) = pplay(ig,l)/pplev(ig,1)*ps(ig)
        enddo
      enddo

!     ---------------------------------------------------------
!     Surface and soil temperature information
!     ---------------------------------------------------------

      Ts1 = 0.0
      Ts2 = 99999.9
      Ts3 = 0.0
      TsS = 0.0 ! mean temperature at bottom soil layer
      do ig=1,ngrid
         Ts1 = Ts1 + area(ig)*tsurf(ig)
         Ts2 = min(Ts2,tsurf(ig))  
         Ts3 = max(Ts3,tsurf(ig))
         TsS = TsS + area(ig)*tsoil(ig,nsoilmx)
      end do
      Ts1=Ts1/totarea
      TsS=TsS/totarea
      if(callsoil)then
         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 = 0.0
         ASR = 0.0
         OLR = 0.0
         GND = 0.0
         DYN = 0.0
         do ig=1,ngrid 
            ISR = ISR + area(ig)*fluxtop_dn(ig) 
            ASR = ASR + area(ig)*fluxabs_sw(ig)
            OLR = OLR + area(ig)*fluxtop_lw(ig) 
            GND = GND + area(ig)*fluxgrd(ig) 
            DYN = DYN + area(ig)*fluxdyn(ig) 
            
            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/totarea,ASR/totarea,OLR/totarea,GND/totarea,DYN/totarea
         
         if(enertest)then
            print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR/totarea,' W m-2'
            print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR/totarea,' W m-2'
            print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR/totarea,' 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/totarea,ASR/totarea,OLR/totarea
               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, exiting physics on first call'
         call abort
      endif

!     ---------------------------------------------------------
!     Compute column amounts (kg m-2) if tracers are enabled
!     ---------------------------------------------------------
      if(tracer)then
         qcol(:,:)=0.0
         do iq=1,nq
            do ig=1,ngrid
               do l=1,nlayer
                  qcol(ig,iq) = qcol(ig,iq) + zq(ig,l,iq) *    &
                            (pplev(ig,l) - pplev(ig,l+1)) / g
               enddo
            enddo
         enddo

         ! not generalised for arbitrary aerosols yet!!!
         reffcol(:,:)=0.0
         do ig=1,ngrid
            do l=1,nlayer
               if(co2cond)then
                  reffcol(ig,1) = reffcol(ig,1) + zq(ig,l,igcm_co2_ice) * &
                                  reffrad(ig,l,1) *    &
                                  (pplev(ig,l) - pplev(ig,l+1)) / g
               endif
               if(water)then
                  reffcol(ig,2) = reffcol(ig,2) + zq(ig,l,igcm_h2o_ice) * &
                                  reffrad(ig,l,2) *    &
                                  (pplev(ig,l) - pplev(ig,l+1)) / g
               endif
            enddo
         enddo

      endif

!     ---------------------------------------------------------
!     Test for water conservation if water is enabled
!     ---------------------------------------------------------

      if(water)then

         icesrf = 0.0
         liqsrf = 0.0
         icecol = 0.0
         vapcol = 0.0

         h2otot = 0.0
         do ig=1,ngrid

            icesrf = icesrf + area(ig)*qsurf_hist(ig,igcm_h2o_ice)
            liqsrf = liqsrf + area(ig)*qsurf_hist(ig,igcm_h2o_vap)
            icecol = icecol + area(ig)*qcol(ig,igcm_h2o_ice)
            vapcol = vapcol + area(ig)*qcol(ig,igcm_h2o_vap)

            h2otot = h2otot + area(ig)*                       &
            (qcol(ig,igcm_h2o_ice)+qcol(ig,igcm_h2o_vap)      &
            +qsurf_hist(ig,igcm_h2o_ice)+qsurf_hist(ig,igcm_h2o_vap))
         end do

         print*,' Total water amount [kg m^-2]: ',h2otot/totarea
         print*,' Surface ice    Surface liq.   Atmos. con.     Atmos. vap. [kg m^-2] '
         print*, icesrf/totarea,liqsrf/totarea,icecol/totarea,vapcol/totarea

         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/totarea,liqsrf/totarea,icecol/totarea,vapcol/totarea
               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))
               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)=0.0
            do l=1,nlayer
               H2Omaxcol(ig) = H2Omaxcol(ig) + qsat(ig,l) *    &
                    (pplev(ig,l) - pplev(ig,l+1))/g
            enddo
         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*100.0 
                     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 
                     !qsurf_hist(ig,igcm_h2o_vap) = 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 
                     !qsurf_hist(ig,igcm_h2o_vap) = 0.0
                  endif

               enddo

               ! enforce ice conservation
               ice_tot=0.0
               do ig = 1, ngrid
                  ice_tot = ice_tot + qsurf_hist(ig,igcm_h2o_ice)*area(ig)
               enddo
               do ig = 1, ngrid
                  qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice)*(icesrf/ice_tot)
               enddo

            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) call writediagfi(ngridmx,"dEzdif","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdif)
	  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) 
          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) call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,1))
        if (igcm_h2o_ice.ne.0) call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(1,1,2))
        if (igcm_co2_ice.ne.0) call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,1))
        if (igcm_h2o_ice.ne.0) call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,2))

!     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(water)then
             call writediagfi(ngridmx,"H2Omaxcol","max. poss. H2O column","kg m^-2",2,H2Omaxcol)
          endif

          if(watercond.or.CLFvarying)then
             !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 STUFF
      if (lastcall) then
        IF ( ALLOCATED( gnom ) ) DEALLOCATE( gnom )    !! this was allocated in su_gases.F90
        IF ( ALLOCATED( gfrac ) ) DEALLOCATE( gfrac )  !! this was allocated in su_gases.F90
      endif


      return
    end subroutine physiq