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

!     creates a temperature profile for an N2-CO2 atmosphere
!     a fine resolution grid is used for accuracy

#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 gamma,dqdlnT,dlnqdlnT,dlnpdlnT,dpdT,dln_p

      integer l,iter

      real sigma_loc,q_loc,vmr_loc,ppvar_loc,T_loc,p_loc
      real Tcond
      real Told,qold,pold

      integer, parameter :: Niter = 2000

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

      l=1

      vmr_loc = 0.05 ! value at ground
      q_loc   = vmr_loc*mv/mn
      T_loc   = Tsurf

      dln_p = log(pceil/plev(1))/Niter ! please check
      do iter=1,Niter

         sigma_loc = (pceil/plev(1))**(real(iter)/real(Niter))
         p_loc     = plev(1)*sigma_loc

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

            if (T_loc.le.Tstrat)then
               T_loc     = Tstrat            ! isothermal stratosphere
            elseif(T_loc.le.Tcond)then

               gamma     = cpp - q_loc*latcond/T_loc
               dqdlnT    = (mv*latcond/(mn*T_loc) - gamma) & 
                            / (latcond/T_loc + Rstar/(mn*q_loc))
               dlnqdlnT  = (1./q_loc)*dqdlnT
               dlnpdlnT  = mv*latcond/(Rstar*T_loc) - dlnqdlnT/(1. + q_loc*mn/mv)
               dpdT      = (p_loc/T_loc)*dlnpdlnT
               q_loc     = q_loc  + dqdlnT*dln_p/dlnpdlnT
               vmr_loc   = q_loc*mn/mv ! needed for Tcond calc.
               T_loc     = T_loc*(1 + dln_p/dlnpdlnT)
               print*,T_loc
            else
               ! nothing doing; so dry convection
               T_loc     = tsurf*(p_loc/plev(1))**rcp
            endif

            ! update true gridpoints when we reach them
            if((p_loc.le.play(l)) .and. (l.le.nlayermx))then
               temp(l)=T_loc
               qvap(l)=q_loc
               ! do an avg you lazy bastard
!               temp(l)=T_loc*(p_loc-play(l))/(p_loc-pold)+Told*(play(l)-pold)/(p_loc-pold)
!               qvap(l)=q_loc*(p_loc-play(l))/(p_loc-pold)+qold*(play(l)-pold)/(p_loc-pold)

!               f=fA*(y2-y)/(y2-y1)+fB*(y-y1)/(y2-y1)

!               print*,'l=',l
!               print*,'temp(l)=',temp(l)
!               print*,'qvap(l)=',qvap(l)


               qvap(l)=0.05*mv/mn ! suface amount everywhere to test

               l=l+1
            endif

            Told=T_loc
            qold=q_loc
            pold=p_loc

      enddo ! vertical loop

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

      return
    end subroutine N2broadprof_fn