MODULE co2condens_mod

      IMPLICIT NONE

      CONTAINS

      SUBROUTINE co2condens(ngrid,nlayer,nq,ptimestep,
     $                  pcapcal,pplay,pplev,ptsrf,pt,
     $                  pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq,
     $                  piceco2,psolaralb,pemisurf,
     $                  pdtc,pdtsrfc,pdpsrf,pduc,pdvc,pdqc,
     $                  fluxsurf_sw,zls,
     $                  zdqssed_co2,pcondicea_co2microp,
     $                  zdtcloudco2)
                                                   
       use tracer_mod, only: noms
       use surfdat_h, only: emissiv, phisfi
       use geometry_mod, only: latitude ! grid point latitudes (rad)
       use planete_h, only: obliquit
       use comcstfi_h, only: cpp, g, r, pi 
#ifndef MESOSCALE
       USE vertical_layers_mod, ONLY: bp
#endif
       IMPLICIT NONE
c=======================================================================
c   subject:
c   --------
c   Condensation/sublimation of CO2 ice on the ground and in the
c   atmosphere 
c   (Scheme described in Forget et al., Icarus, 1998)
c
c   author:   Francois Forget     1994-1996 ; updated 1996 -- 2018
c   ------
c             adapted to external CO2 ice clouds scheme by Deborah Bardet (2018) '
c
c=======================================================================
c
c    0.  Declarations :
c    ------------------
c
      include "callkeys.h"

c-----------------------------------------------------------------------
c    Arguments :
c    ---------
      INTEGER,INTENT(IN) :: ngrid  ! number of atmospheric columns
      INTEGER,INTENT(IN) :: nlayer ! number of vertical layers
      INTEGER,INTENT(IN) :: nq  ! number of tracers

      REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
      REAL,INTENT(IN) :: pcapcal(ngrid)
      REAL,INTENT(IN) :: pplay(ngrid,nlayer) !mid-layer pressure (Pa)
      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
      REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
      REAL,INTENT(IN) :: pt(ngrid,nlayer) ! atmospheric temperature (K)
      REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential (m2.s-2)
      REAL,INTENT(IN) :: pdt(ngrid,nlayer) ! tendency on temperature from
                                           ! previous physical processes (K/s)
      REAL,INTENT(IN) :: pdu(ngrid,nlayer) ! tendency on zonal wind (m/s2)
                                           ! from previous physical processes
      REAL,INTENT(IN) :: pdv(ngrid,nlayer) ! tendency on meridional wind (m/s2)
                                           ! from previous physical processes
      REAL,INTENT(IN) :: pdtsrf(ngrid) ! tendency on surface temperature from
                                       ! previous physical processes (K/s)
      REAL,INTENT(IN) :: pu(ngrid,nlayer) ! zonal wind (m/s)
      REAL,INTENT(IN) :: pv(ngrid,nlayer) ! meridional wind (m/s)
      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (../kg_air)
      REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq) ! tendency on tracers from
                                              ! previous physical processes

      REAL,INTENT(IN) :: zdqssed_co2(ngrid)  ! CO2 flux at the surface  (kg.m-2.s-1)
      REAL,INTENT(IN) :: pcondicea_co2microp(ngrid,nlayer)! tendency due to CO2 condensation (kg/kg.s-1)
      REAL,INTENT(IN) :: zdtcloudco2(ngrid,nlayer) ! tendency on temperature due to latent heat                 

      REAL,INTENT(INOUT) :: piceco2(ngrid) ! CO2 ice on the surface (kg.m-2)
      REAL,INTENT(INOUT) :: psolaralb(ngrid,2) ! albedo of the surface
      REAL,INTENT(INOUT) :: pemisurf(ngrid) ! emissivity of the surface
      
      ! tendencies due to CO2 condensation/sublimation:
      REAL,INTENT(OUT) :: pdtc(ngrid,nlayer) ! tendency on temperature (K/s)
      REAL,INTENT(OUT) :: pdtsrfc(ngrid) ! tendency on surface temperature (K/s)
      REAL,INTENT(OUT) :: pdpsrf(ngrid) ! tendency on surface pressure (Pa/s)
      REAL,INTENT(OUT) :: pduc(ngrid,nlayer) ! tendency on zonal wind (m.s-2)
      REAL,INTENT(OUT) :: pdvc(ngrid,nlayer) ! tendency on meridional wind (m.s-2)
      REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq) ! tendency on tracers
      
      ! added to calculate flux dependent albedo:
      REAL,intent(in) :: fluxsurf_sw(ngrid,2)
      real,intent(in) :: zls ! solar longitude (rad)

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

      INTEGER i,j
      INTEGER l,ig,iq,icap
      REAL zt(ngrid,nlayer)
      REAL zcpi
      REAL ztcond (ngrid,nlayer+1) ! CO2 condensation temperature (atm)
      REAL ztcondsol(ngrid) ! CO2 condensation temperature (surface)
      REAL zdiceco2(ngrid)
      REAL zcondicea(ngrid,nlayer) ! condensation rate in layer  l  (kg/m2/s)
      REAL zcondices(ngrid) ! condensation rate on the ground (kg/m2/s)
      REAL zfallice(ngrid,nlayer+1) ! amount of ice falling from layer l (kg/m2/s)
      REAL zfallheat
      REAL zmflux(nlayer+1)
      REAL zu(nlayer),zv(nlayer)
      REAL zq(nlayer,nq),zq1(nlayer)
      REAL ztsrf(ngrid)
      REAL ztc(nlayer), ztm(nlayer+1) 
      REAL zum(nlayer+1) , zvm(nlayer+1)
      REAL zqm(nlayer+1,nq),zqm1(nlayer+1)
      REAL masse(nlayer),w(nlayer+1)
      REAL Sm(nlayer),Smq(nlayer,nq),mixmas,qmix
      LOGICAL condsub(ngrid)

      real :: emisref(ngrid)

c variable speciale diagnostique
      real tconda1(ngrid,nlayer)
      real tconda2(ngrid,nlayer)
c     REAL zdiceco2a(ngrid) ! for diagnostic only
      real zdtsig (ngrid,nlayer)
      real zdt (ngrid,nlayer)
      real vmr_co2(ngrid,nlayer) ! co2 volume mixing ratio
! improved_ztcond flag: If set to .true. (AND running with a 'co2' tracer)
! then condensation temperature is computed using partial pressure of CO2
      logical,parameter :: improved_ztcond=.true.

c   local saved variables
      integer,save :: ico2 ! index of CO2 tracer
      real,save :: qco2min,qco2,mmean
      real,parameter :: latcond=5.9e5 ! (J/kg) Latent heat of solid CO2 ice
      real,parameter :: tcond1mb=136.27 ! condensation temperature (K) at 1 mbar
      real,parameter :: cpice=1000. ! (J.kg-1.K-1) specific heat of CO2 ice
      REAL,SAVE :: acond,bcond,ccond
      real,save :: m_co2, m_noco2, A , B

      LOGICAL,SAVE :: firstcall = .true. !,firstcall2=.true.

c D.BARDET: to debug 
      real ztc3D(ngrid,nlayer)
      REAL ztm3D(ngrid,nlayer)
      REAL zmflux3D(ngrid,nlayer)
c----------------------------------------------------------------------

c   Initialisation
c   --------------
c
      ! AS: firstcall OK absolute
      IF (firstcall) THEN
         
         bcond=1./tcond1mb
         ccond=cpp/(g*latcond)
         acond=r/latcond

         firstcall=.false.
         write(*,*) 'Newcondens: improved_ztcond=',improved_ztcond
         PRINT*,'In newcondens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
         PRINT*,'acond,bcond,ccond',acond,bcond,ccond

         ico2=0

         if (tracer) then
c          Prepare Special treatment if one of the tracer is CO2 gas
           do iq=1,nq
             if (noms(iq).eq."co2") then
                ico2=iq
                m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
                m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
c               Compute A and B coefficient use to compute
c               mean molecular mass Mair defined by
c               1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
c               1/Mair = A*q(ico2) + B
                A =(1/m_co2 - 1/m_noco2)
                B=1/m_noco2
             endif
           enddo
c          minimum CO2 mix. ratio below which mixing occur with layer above: 
           qco2min =0.75  
         end if
      ENDIF ! of IF (firstcall)
      zcpi=1./cpp

c
c======================================================================
c    Calcul of CO2 condensation sublimation 
c    ============================================================
c  
c    Used variable :
c       piceco2(ngrid)   :  amount of co2 ice on the ground (kg/m2)
c       zcondicea(ngrid,l):  condensation rate in layer  l  (kg/m2/s)
c       zcondices(ngrid):  condensation rate on the ground (kg/m2/s)
c       zfallice(ngrid,l):amount of ice falling from layer l (kg/m2/s)
c                           
c       pdtc(ngrid,nlayer)  : dT/dt due to cond/sub
c
c
c     Tendencies set to 0 (except pdtc)
c     -------------------------------------
      DO l=1,nlayer
         DO ig=1,ngrid
           zcondicea(ig,l) = 0.
           zfallice(ig,l) = 0.
           pduc(ig,l)  = 0
           pdvc(ig,l)  = 0
         END DO
      END DO
         
      DO iq=1,nq
      DO l=1,nlayer
         DO ig=1,ngrid
           pdqc(ig,l,iq)  = 0
         END DO
      END DO
      END DO

      DO ig=1,ngrid
         zfallice(ig,nlayer+1) = 0.
         zcondices(ig) = 0.
         pdtsrfc(ig) = 0.
         pdpsrf(ig) = 0.
         condsub(ig) = .false.
         zdiceco2(ig) = 0.
      ENDDO
      zfallheat=0

c     *************************
c     ATMOSPHERIC CONDENSATION
c     *************************

c     Compute CO2 Volume mixing ratio
c     -------------------------------
      if (improved_ztcond.and.(ico2.ne.0)) then
         DO l=1,nlayer
            DO ig=1,ngrid
              qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep
c             Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
              mmean=1/(A*qco2 +B)
              vmr_co2(ig,l) = qco2*mmean/m_co2 
            ENDDO
         ENDDO
      else
         DO l=1,nlayer
            DO ig=1,ngrid
              vmr_co2(ig,l)=0.95
            ENDDO
         ENDDO
      end if

      IF (.NOT. co2clouds) then
c     forecast of atmospheric temperature zt and frost temperature ztcond
c     --------------------------------------------------------------------

      DO l=1,nlayer
         DO ig=1,ngrid
            zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
!            ztcond(ig,l)=1./(bcond-acond*log(.0095*pplay(ig,l)))
            if (pplay(ig,l).ge.1e-4) then
              ztcond(ig,l)=
     &         1./(bcond-acond*log(.01*vmr_co2(ig,l)*pplay(ig,l)))
            else
              ztcond(ig,l)=0.0 !mars Monica
            endif
         ENDDO
      ENDDO

      ztcond(:,nlayer+1)=ztcond(:,nlayer)
 
c      Condensation/sublimation in the atmosphere
c      ------------------------------------------
c      (calcul of zcondicea , zfallice and pdtc)
c
      DO l=nlayer , 1, -1
         DO ig=1,ngrid
           pdtc(ig,l)=0.
           IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
               condsub(ig)=.true.
               IF (zfallice(ig,l+1).gt.0) then  
                 zfallheat=zfallice(ig,l+1)*
     &           (pphi(ig,l+1)-pphi(ig,l) +
     &          cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/latcond
               ELSE
                   zfallheat=0.
               ENDIF
               pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
               zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
     &                        *ccond*pdtc(ig,l)- zfallheat
c              Case when the ice from above sublimes entirely
c              """""""""""""""""""""""""""""""""""""""""""""""
               IF (zfallice(ig,l+1).lt.- zcondicea(ig,l)) then
                  pdtc(ig,l)=(-zfallice(ig,l+1)+zfallheat)/
     &              (ccond*(pplev(ig,l)-pplev(ig,l+1)))
                  zcondicea(ig,l)= -zfallice(ig,l+1)
               END IF

               zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
            END IF
         ENDDO
      ENDDO
      
      ELSE
         
        DO ig=1,ngrid 
         zfallice(ig,1) = zdqssed_co2(ig)
        ENDDO
        DO l=nlayer , 1, -1
            DO ig=1,ngrid 
         zcondicea(ig,l) = pcondicea_co2microp(ig,l)* 
     &                        (pplev(ig,l) - pplev(ig,l+1))/g
            ENDDO
        ENDDO
        DO l=nlayer, 1, -1
         DO ig=1, ngrid
            zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
            pdtc(ig,l)=0.
         ENDDO
      ENDDO 
      ENDIF ! if not co2clouds

      call WRITEdiagfi(ngrid,"pdtc_atm",
     &         "temperature tendency due to CO2 condensation",
     &         " ",3,pdtc) 
     
       call WRITEdiagfi(ngrid,"zcondicea",
     &         "",
     &         " ",3,zcondicea)         
      
       call WRITEdiagfi(ngrid,"zfallice",
     &         "",
     &         " ",2,zfallice(ngrid,1))

c     *************************
c       SURFACE  CONDENSATION
c     *************************

c     forecast of ground temperature ztsrf and frost temperature ztcondsol
c     --------------------------------------------------------------------
      DO ig=1,ngrid
         ztcondsol(ig)=
     &          1./(bcond-acond*log(.01*vmr_co2(ig,1)*pplev(ig,1)))
         ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
      ENDDO

c
c      Condensation/sublimation on the ground 
c      --------------------------------------
c      (compute zcondices and pdtsrfc)
c
      DO ig=1,ngrid
         IF(latitude(ig).lt.0) THEN
           ! Southern hemisphere
            icap=2
         ELSE
           ! Northern hemisphere
            icap=1
         ENDIF
       
c
c        Loop on where we have  condensation/ sublimation 
         IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR.   ! ground cond 
     $       (zfallice(ig,1).NE.0.) .OR.           ! falling snow
     $      ((ztsrf(ig) .GT. ztcondsol(ig)) .AND.  ! ground sublim.
     $      ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN
            condsub(ig) = .true. 

            IF (zfallice(ig,1).gt.0) then  
                 zfallheat=zfallice(ig,1)*
     &           (pphi(ig,1)- phisfi(ig) +
     &           cpice*(ztcond(ig,1)-ztcondsol(ig)))/latcond
            ELSE
                 zfallheat=0.
            ENDIF

c           condensation or partial sublimation of CO2 ice
c           """""""""""""""""""""""""""""""""""""""""""""""
            zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) 
     &      /(latcond*ptimestep)         - zfallheat
            pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
       
c           If the entire CO_2 ice layer sublimes
c           """"""""""""""""""""""""""""""""""""""""""""""""""""
c           (including what has just condensed in the atmosphere)

            IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE.
     &          -zcondices(ig))THEN
               zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1) 
               pdtsrfc(ig)=(latcond/pcapcal(ig))*
     &         (zcondices(ig)+zfallheat)
            END IF

c           Changing CO2 ice amount and pressure :
c           """"""""""""""""""""""""""""""""""""

            zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
            piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
            pdpsrf(ig) = -zdiceco2(ig)*g

            IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
               PRINT*,'STOP in condens'
               PRINT*,'condensing more than total mass'
               PRINT*,'Grid point ',ig
               PRINT*,'Ps = ',pplev(ig,1)
               PRINT*,'d Ps = ',pdpsrf(ig)
               STOP
            ENDIF
         END IF ! if there is condensation/sublimmation
      ENDDO ! of DO ig=1,ngrid

c  ********************************************************************
c  Surface albedo and emissivity of the surface below the snow (emisref)
c  ********************************************************************

!     Check that amont of CO2 ice is not problematic
      DO ig=1,ngrid
           if(.not.piceco2(ig).ge.0.) THEN
              if(piceco2(ig).le.-5.e-8) print*,
     $         'WARNING newcondens piceco2(',ig,')=', piceco2(ig)
               piceco2(ig)=0.
           endif
      ENDDO
      
!     Set albedo and emissivity of the surface
!     ----------------------------------------
      CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref)

! set pemisurf() to emissiv when there is bare surface (needed for co2snow)
      DO ig=1,ngrid
        if (piceco2(ig).eq.0) then
          pemisurf(ig)=emissiv
        endif
      ENDDO

!         firstcall2=.false.
c ***************************************************************
c  Correction to account for redistribution between sigma or hybrid 
c  layers when changing surface pressure (and warming/cooling
c  of the CO2 currently changing phase).
c *************************************************************

      DO ig=1,ngrid
        if (condsub(ig)) then
           do l=1,nlayer
             ztc(l)  =zt(ig,l)   +pdtc(ig,l)  *ptimestep
             zu(l)   =pu(ig,l)   +pdu( ig,l)  *ptimestep
             zv(l)   =pv(ig,l)   +pdv( ig,l)  *ptimestep
            do iq=1,nq
             zq(l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
            enddo
           end do

c  Mass fluxes through the sigma levels (kg.m-2.s-1)  (>0 when up)
c  """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

            zmflux(1) = -zcondices(ig)
            DO l=1,nlayer
             zmflux(l+1) = zmflux(l) -zcondicea(ig,l)
#ifndef MESOSCALE
     &        + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld  
          if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
#else
          if (abs(zmflux(l+1)).lt.1E-13) zmflux(l+1)=0.
#endif
            END DO

c Mass of each layer
c ------------------ 
            DO l=1,nlayer
              masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
            END DO


c  Corresponding fluxes for T,U,V,Q
c  """"""""""""""""""""""""""""""""

c           averaging operator for TRANSPORT  
c           """"""""""""""""""""""""""""""""
c           Value transfert at the surface interface when condensation
c           sublimation:
            ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
            zum(1) = 0  
            zvm(1) = 0  
            do iq=1,nq
              zqm(1,iq)=0. ! most tracer do not condense !
            enddo
c           Special case if one of the tracer is CO2 gas
            if (ico2.ne.0) zqm(1,ico2)=1. ! flux is 100% CO2

c           Van Leer scheme:
            DO l=1,nlayer+1
                w(l)=-zmflux(l)*ptimestep
            END DO
            call vl1d(nlayer,ztc,2.,masse,w,ztm) 
            call vl1d(nlayer,zu ,2.,masse,w,zum) 
            call vl1d(nlayer,zv ,2.,masse,w,zvm) 
            do iq=1,nq
             do l=1,nlayer
              zq1(l)=zq(l,iq)
             enddo
             zqm1(1)=zqm(1,iq)
             call vl1d(nlayer,zq1,2.,masse,w,zqm1)
             do l=2,nlayer
              zq( l,iq)=zq1(l)
              zqm(l,iq)=zqm1(l)
             enddo
            enddo

c           Surface condensation affects low winds
            if (zmflux(1).lt.0) then 
                zum(1)= zu(1) *  (w(1)/masse(1))
                zvm(1)= zv(1) *  (w(1)/masse(1))
                if (w(1).gt.masse(1)) then ! ensure numerical stability
                  zum(1)= (zu(1)-zum(2))*masse(1)/w(1) + zum(2)
                  zvm(1)= (zv(1)-zvm(2))*masse(1)/w(1) + zvm(2)
                end if
            end if
                    
            ztm(nlayer+1)= ztc(nlayer) ! should not be used, but... 
            zum(nlayer+1)= zu(nlayer)  ! should not be used, but...
            zvm(nlayer+1)= zv(nlayer)  ! should not be used, but...
            do iq=1,nq
             zqm(nlayer+1,iq)= zq(nlayer,iq)
            enddo

#ifdef MESOSCALE
!!!! AS: This part must be commented in the mesoscale model
!!!! AS: ... to avoid instabilities.
!!!! AS: you have to compile with -DMESOSCALE to do so
#else  
c           Tendencies on T, U, V, Q 
c           """"""""""""""""""""""""
            DO l=1,nlayer
               IF(.not. co2clouds) THEN
c             Tendencies on T
                zdtsig(ig,l) = (1/masse(l)) *
     &        ( zmflux(l)*(ztm(l) - ztc(l)) 
     &        - zmflux(l+1)*(ztm(l+1) - ztc(l))
     &        + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l))  )
               ELSE
                zdtsig(ig,l) = (1/masse(l)) *
     &        ( zmflux(l)*(ztm(l) - ztc(l)) 
     &        - zmflux(l+1)*(ztm(l+1) - ztc(l)))               
               ENDIF
c D.BARDET: for diagnotics 
                zmflux3D(ig,l)=zmflux(l)
                ztm3D(ig,l)=ztm(l)
                ztc3D(ig,l)=ztc(l)
                
                pdtc(ig,l) =  pdtc(ig,l) + zdtsig(ig,l)

c             Tendencies on U
                pduc(ig,l)   = (1/masse(l)) *
     &        ( zmflux(l)*(zum(l) - zu(l))
     &        - zmflux(l+1)*(zum(l+1) - zu(l)) )


c             Tendencies on V
                pdvc(ig,l)   = (1/masse(l)) *
     &        ( zmflux(l)*(zvm(l) - zv(l))
     &        - zmflux(l+1)*(zvm(l+1) - zv(l)) )

            END DO

#endif

c           Tendencies on Q
            do iq=1,nq
!              if (noms(iq).eq.'co2') then 
              if (iq.eq.ico2) then 
c               SPECIAL Case when the tracer IS CO2 :
                DO l=1,nlayer
                  pdqc(ig,l,iq)= (1/masse(l)) *
     &           ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
     &           - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) 
     &           + zcondicea(ig,l)*(zq(l,iq)-1.) )
                END DO
              else
                DO l=1,nlayer
                  pdqc(ig,l,iq)= (1/masse(l)) *
     &           ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
     &           - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) 
     &           + zcondicea(ig,l)*zq(l,iq) )
                END DO
              end if
            enddo

       end if ! if (condsub)
      END DO  ! loop on ig 

c ***************************************************************
c   CO2 snow / clouds scheme
c ***************************************************************

      call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
     &        zcondicea,zcondices,zfallice,pemisurf)

c ***************************************************************
c Ecriture des diagnostiques
c ***************************************************************

c     DO l=1,nlayer
c        DO ig=1,ngrid
c         Taux de cond en kg.m-2.pa-1.s-1
c          tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
c          Taux de cond en kg.m-3.s-1
c          tconda2(ig,l)=tconda1(ig,l)*pplay(ig,l)*g/(r*pt(ig,l))
c        END DO
c     END DO
c     call WRITEDIAGFI(ngrid,'tconda1',
c    &'Taux de condensation CO2 atmospherique /Pa',
c    & 'kg.m-2.Pa-1.s-1',3,tconda1)
c     call WRITEDIAGFI(ngrid,'tconda2',
c    &'Taux de condensation CO2 atmospherique /m',
c    & 'kg.m-3.s-1',3,tconda2)

! output falling co2 ice in 1st layer:
!      call WRITEDIAGFI(ngrid,'fallice',
!     &'Precipitation of co2 ice',
!     & 'kg.m-2.s-1',2,zfallice(1,1))

#ifndef MESOSCALE 
! Extra special case for surface temperature tendency pdtsrfc:
! we want to fix the south pole temperature to CO2 condensation temperature
         if (caps.and.(obliquit.lt.27.)) then
           ! check if last grid point is the south pole
           if (abs(latitude(ngrid)-(-pi/2.)).lt.1.e-5) then
           ! NB: Updated surface pressure, at grid point 'ngrid', is
           !     ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
!             write(*,*) "newcondens: South pole: latitude(ngrid)=",
!     &                                           latitude(ngrid)
             ztcondsol(ngrid)=
     &          1./(bcond-acond*log(.01*vmr_co2(ngrid,1)*
     &                    (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
             pdtsrfc(ngrid)=(ztcondsol(ngrid)-ztsrf(ngrid))/ptimestep
           endif
         endif
#endif

      END SUBROUTINE co2condens



c *****************************************************************
      SUBROUTINE vl1d(nlayer,q,pente_max,masse,w,qm)
c
c    
c     Operateur de moyenne inter-couche pour calcul de transport type
c     Van-Leer " pseudo amont " dans la verticale
c    q,w sont des arguments d'entree  pour le s-pg ....
c    masse : masse de la couche Dp/g
c    w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
c    pente_max = 2 conseillee
c
c
c   --------------------------------------------------------------------
      IMPLICIT NONE

c
c
c
c   Arguments:
c   ----------
      integer nlayer
      real masse(nlayer),pente_max
      REAL q(nlayer),qm(nlayer+1)
      REAL w(nlayer+1)
c
c      Local 
c   ---------
c
      INTEGER l
c
      real dzq(nlayer),dzqw(nlayer),adzqw(nlayer),dzqmax
      real sigw, Mtot, MQtot
      integer m 
c     integer ismax,ismin 


c    On oriente tout dans le sens de la pression 
c     W > 0 WHEN DOWN !!!!!!!!!!!!!

      do l=2,nlayer
            dzqw(l)=q(l-1)-q(l)
            adzqw(l)=abs(dzqw(l))
      enddo

      do l=2,nlayer-1
            if(dzqw(l)*dzqw(l+1).gt.0.) then
                dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
            else
                dzq(l)=0.
            endif
            dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
            dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
      enddo

         dzq(1)=0.
         dzq(nlayer)=0.

       do l = 1,nlayer-1

c         Regular scheme (transfered mass < layer mass)
c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
             sigw=w(l+1)/masse(l+1)
             qm(l+1)=(q(l+1)+0.5*(1.-sigw)*dzq(l+1))
          else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
             sigw=w(l+1)/masse(l)
             qm(l+1)=(q(l)-0.5*(1.+sigw)*dzq(l))

c         Extended scheme (transfered mass > layer mass)
c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          else if(w(l+1).gt.0.) then
             m=l+1
             Mtot = masse(m)
             MQtot = masse(m)*q(m)
             do while ((m.lt.nlayer).and.(w(l+1).gt.(Mtot+masse(m+1))))
                m=m+1
                Mtot = Mtot + masse(m)
                MQtot = MQtot + masse(m)*q(m)
             end do
             if (m.lt.nlayer) then
                sigw=(w(l+1)-Mtot)/masse(m+1)
                qm(l+1)= (1/w(l+1))*(MQtot + (w(l+1)-Mtot)*
     &          (q(m+1)+0.5*(1.-sigw)*dzq(m+1)) )
             else
                w(l+1) = Mtot
                qm(l+1) = Mqtot / Mtot
                write(*,*) 'top layer is disapearing !'
                stop
             end if
          else      ! if(w(l+1).lt.0) 
             m = l-1 
             Mtot = masse(m+1)
             MQtot = masse(m+1)*q(m+1)
             if (m.gt.0) then ! because some compilers will have problems
                              ! evaluating masse(0)
              do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
                m=m-1
                Mtot = Mtot + masse(m+1)
                MQtot = MQtot + masse(m+1)*q(m+1)
                if (m.eq.0) exit
              end do
             endif
             if (m.gt.0) then
                sigw=(w(l+1)+Mtot)/masse(m)
                qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)* 
     &          (q(m)-0.5*(1.+sigw)*dzq(m))  )
             else
                qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*qm(1))
             end if   
          end if
       enddo

c     boundary conditions (not used in newcondens !!)
c         qm(nlayer+1)=0.
c         if(w(1).gt.0.) then
c            qm(1)=q(1)
c         else 
c           qm(1)=0.
c         end if

      END SUBROUTINE vl1d
      
      END MODULE co2condens_mod