subroutine kastprof_fn(tsurf,temp,popsk,play,plev,qvap)
      implicit none

!     creates a kasting temperature profile

#include "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "callkeys.h"

      REAL play(nlayermx)   ! Pressure at the middle of the layers (Pa)
      REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
      REAL psurf,tsurf      
      REAL tcond(nlayermx)
      REAL qvap(nlayermx)
      REAL temp(nlayermx)   ! temperature at the middle of the layers

      real, parameter :: Rstar   = 8.314
      real, parameter :: ptriple = 518000.0
      real, parameter :: latcond = 5.9e5     ! as in condense_co2
      real, parameter :: mn      = 0.02801   ! molar mass of background gas (N2)
      real, parameter :: mv      = 0.04401   ! molar mass of background gas (CO2)

      real popsk(nlayermx)  ! this is dodgy as in rest of code, zpopsk(ngridmx,nlayermx)
                            ! but kastprof requires 1D, so its ok I guess

      real ppco2
      real vmr_local,q_local,gamma,dqdlnT,dlnqdlnT,dlnpdlnT,dpdT,dln_p

      integer l,iter

      logical co2mixedcond

      co2mixedcond=.true.

      if(ngridmx.gt.1)then
         print*,'DO NOT RUN KASTPROF IN 3D!!!'
         call abort
      end if

      ! 1st, we make pot. temperature constant everywhere
      do l=1,nlayermx
            temp(l)=TsurfK*popsk(l)
      enddo

!     then do all the other stuff
      vmr_local = 1.-n2mixratio
      q_local   = vmr_local*mv/mn

      qvap(1)   = q_local

      ppco2=vmr_local*play(1)

      if(play(1).lt.ptriple)then !
         tcond(1) = (-3167.8)/(log(.01*ppco2)-23.23)         ! Fanale's formula
      else
         tcond(1) = 684.2-92.3*log(ppco2)+4.32*log(ppco2)**2 ! liquid-vapour (CRC handbook 2003 data)
      endif

      do l=1,nlayermx-1

         ppco2=vmr_local*play(l) ! this shld be implicit i.e. vmr_local(l+1) ??

         if(play(l).le.ptriple)then ! ERROR
            tcond(l) = (-3167.8)/(log(.01*ppco2)-23.23)         ! Fanale's formula
         else
            tcond(l) = 684.2-92.3*log(ppco2)+4.32*log(ppco2)**2 ! liquid-vapour (CRC handbook 2003 data)
         endif

         if(co2mixedcond)then
            if (temp(l) .le. Tstrat)then
               temp(l)  = Tstrat            ! isothermal stratosphere
               temp(l+1) = Tstrat
               qvap(l+1) = qvap(l)
               !print*,'strato'

            elseif(temp(l) .le. tcond(l))then

               dln_p     = log(play(l+1))-log(play(l))
               gamma     = cpp - q_local*latcond/temp(l)
               dqdlnT    = (mv*latcond/(mn*temp(l)) - gamma) & 
                            / (latcond/temp(l) + Rstar/(mn*q_local))
               dlnqdlnT  = (1./q_local)*dqdlnT
               dlnpdlnT  = mv*latcond/(Rstar*temp(l)) - dlnqdlnT/(1. + q_local*mn/mv)
               dpdT      = (play(l)/temp(l))*dlnpdlnT
               q_local   = q_local + dqdlnT*dln_p/dlnpdlnT

               vmr_local = q_local*mn/mv


               if(1.eq.1)then
               print*,'l         =',l
               print*,'gamma     =',gamma
               print*,'gamma bit =',q_local*latcond/temp(l)
               print*,'cpp       =',cpp
               print*,'latcond   =',latcond
               print*,'q_local   =',q_local
               print*,'T         =',temp(l)
               print*,'p         =',play(l)
               print*,'dqdlnT    =',dqdlnT
               print*,'dlnqdlnT  =',dlnqdlnT
               print*,'dlnpdlnT  =',dlnpdlnT
               print*,'dpdT      =',dpdT
               call abort
               endif

               qvap(l+1) = q_local
               temp(l+1) = temp(l) + temp(l)*dln_p/(dlnpdlnT)
               !print*,'CO2co'
            else
               !print*,'tropo'


               ! CO2 is not condensing; therefore dry convection
               qvap(l+1) = qvap(l)
            endif

         else

            if (temp(l) .le. Tstrat)then
               temp(l) = Tstrat            
               ! isothermal stratosphere
            elseif(temp(l) .le. tcond(l))then
               temp(l)=tcond(l)
               ! pure CO2 condensation
            else
               ! CO2 is not condensing; therefore dry convection 
            endif

         endif

      enddo ! vertical loop

      temp(nlayermx) = Tstrat
      qvap(nlayermx) = qvap(nlayermx-1)

      print*,'---fin---'

      call abort

      return
    end subroutine kastprof_fn