PRO COSMO_PIPE

  COMMON psm, param
  COMMON sky, model

  cosmodir = GETPSMDIR('cosmo')
  dir_exists = FILE_TEST(cosmodir)
  IF ((dir_exists EQ 0)) THEN FILE_MKDIR, cosmodir

  cmbdir = GETPSMDIR('cmb')
  dir_exists = FILE_TEST(cmbdir)
  IF ((dir_exists EQ 0)) THEN FILE_MKDIR, cmbdir

  ;; Clear the directories if necessary   
  ;;------------------------------------
  docosmo = ((STRLOWCASE(CMB_MODEL()) NE 'no_cmb') OR (STRLOWCASE(SZ_MODEL()) NE 'no_sz')) AND SKY_TASK() EQ 'new'
  IF docosmo THEN CLEAR_PSM_DIRECTORY, 'cosmo'
  docmb = (STRLOWCASE(CMB_MODEL()) NE 'no_cmb') AND (SKY_TASK() EQ 'new')
  IF docmb THEN CLEAR_PSM_DIRECTORY, 'cmb'

  ;; Filenames to write the theory CMB cl
  ;;--------------------------------------
  cmb_lensedcl_file = cmbdir + 'cmb' + '_lensed_cl.fits' 
  cmb_unlensedcl_file = cmbdir + 'cmb' + '_unlensed_cl.fits' 

  ;; Set sky resolution
  skyres = SKY_RESOLUTION()
  ;; Set sky nside
  skynside = SKY_NSIDE()
  ;; Set sky lmax
  skylmax = SKY_LMAX() 
  ;; Set polar
  IF STRUPCASE(param.global.fields) EQ 'TP' THEN polar = 1 ELSE polar = 0  
  ;; Set sky pixelisation
  skypix = SKY_PIX()

  ;;--------------------------------------------------------------------------------------------------------------------;;
  ;; Get input cl 
  ;;--------------------------------------------------------------------------------------------------------------------;;

  CASE model.cmb.model OF

     'gaussian': BEGIN
        CASE model.cmb.cmb_cl_source OF
           ;; Input cl read from file
           'standard-LCDM': BEGIN
              unlensedcl = CMB_BESTFIT_CL(/UNLENSED) 
              lensedcl = CMB_BESTFIT_CL()
           END
           ;; Input cl obtained by running CAMB
           'CAMB': BEGIN

              ;; We have to change the current working directory to the
              ;; camb directory to run post July 2011 versions of
              ;; camb. This is done below.

              ;; Find camb path, and cd to that directory
              cambfullpath = FILE_DIRNAME(FILE_WHICH(GETENV('PATH'), 'camb'))
              CD, cambfullpath, current=currentdir

              ;; Get an absolute filename for the cosmodir directory
              isrelative = FIRSTCHAR(cosmodir,2) EQ './'
              IF isrelative THEN cdir = FIX_SEPARATOR(currentdir)+STRMID(cosmodir,2,STRLEN(cosmodir)-2) ELSE cdir=cosmodir

              ;; write the CAMB parameter file, run CAMB, store CAMB ouputs in cmb directory
              WRITE_CAMB_PARAM_FILE, model.cosmo, cdir+'camb.ini', lmax=skylmax, polar=polar, output_root=FIX_SEPARATOR(cdir)+'camb'
              IF PSM_VERBOSITY() GE 1 THEN PSM_INFO_SUBTITLE, 'RUNNING CAMB'

              SPAWN, 'camb ' +cdir+'camb.ini'

              ;; Get back to original working directory
              CD, currentdir
              
              ;; unlensed cl
              ;;-------------
              READ_CAMB_CL, cosmodir+'camb_scalCls.txt', cl, polar=polar
              IF model.cosmo.r GT 0 THEN BEGIN
                 READ_CAMB_CL, cosmodir+'camb_tensCls.txt', cl_tens, polar=polar
                 cl.C_1_1[cl.L] = cl.C_1_1[cl.L] + cl_tens.C_1_1[cl.L]
                 IF polar EQ 1 THEN BEGIN
                    cl.C_2_2[cl.L] = cl.C_2_2[cl.L] + cl_tens.C_2_2[cl.L]
                    cl.C_3_3[cl.L] = cl.C_3_3[cl.L] + cl_tens.C_3_3[cl.L]
                    cl.C_1_2[cl.L] = cl.C_1_2[cl.L] + cl_tens.C_1_2[cl.L]
                 ENDIF
              ENDIF
              unlensedcl = cl
              
              ;; lensed cl
              ;;-----------
              IF model.cosmo.r GT 0 THEN BEGIN
                 READ_CAMB_CL, cosmodir+'camb_lensedtotCls.txt', lensedcl, polar=polar
              ENDIF ELSE BEGIN
                 READ_CAMB_CL, cosmodir+'camb_lensedCls.txt', lensedcl, polar=polar
              ENDELSE
           END
           
        ENDCASE
        PSM_ID_KEYGEN, cmb_id
        hbase = PSM_BASE_HDR(datatype='comp',dataform='cl', /setval)
        hcomp = PSM_COMP_HDR(compname='cmb',compid=cmb_id)
        MK_PSMHDR, xhdr, hbase, hcomp 
        WRITE_CL, lensedcl, cmb_lensedcl_file, pxwin=0, beam=0, xhdr=xhdr
        WRITE_CL, unlensedcl, cmb_unlensedcl_file, pxwin=0, beam=0, xhdr=xhdr
        ;;stop
     END 

     'nongaussian_fnl': BEGIN
        ;;stop
        WRITE_CAMB_PARAM_FILE, model.cosmo, cosmodir+'camb.ini', lmax=skylmax, polar=polar
        IF PSM_VERBOSITY() GE 1 THEN PSM_INFO_SUBTITLE, 'RUNNING CAMB'
        SPAWN, 'camb ' +cosmodir+'camb.ini'
        
        ;; unlensed cl
        ;;-------------
        READ_CAMB_CL, cosmodir+'camb_scalCls.txt', cl, polar=polar
        IF model.cosmo.r GT 0 THEN BEGIN
           READ_CAMB_CL, cosmodir+'camb_tensCls.txt', cl_tens, polar=polar
           cl.C_1_1[cl.L] = cl.C_1_1[cl.L] + cl_tens.C_1_1[cl.L]
           IF polar EQ 1 THEN BEGIN
              cl.C_2_2[cl.L] = cl.C_2_2[cl.L] + cl_tens.C_2_2[cl.L]
              cl.C_3_3[cl.L] = cl.C_3_3[cl.L] + cl_tens.C_3_3[cl.L]
              cl.C_1_2[cl.L] = cl.C_1_2[cl.L] + cl_tens.C_1_2[cl.L]
           ENDIF
        ENDIF
        unlensedcl = cl
        
        ;; lensed cl
        ;;-----------
        IF model.cosmo.r GT 0 THEN BEGIN
           READ_CAMB_CL, cosmodir+'camb_lensedtotCls.txt', lensedcl, polar=polar
        ENDIF ELSE BEGIN
           READ_CAMB_CL, cosmodir+'camb_lensedCls.txt', lensedcl, polar=polar
        ENDELSE

        ;; write cl
        ;;----------
        PSM_ID_KEYGEN, cmb_id
        hbase = PSM_BASE_HDR(datatype='comp',dataform='cl', /setval)
        hcomp = PSM_COMP_HDR(compname='cmb',compid=cmb_id)
        MK_PSMHDR, xhdr, hbase, hcomp 
        WRITE_CL, lensedcl, cmb_lensedcl_file, pxwin=0, beam=0, xhdr=xhdr
        WRITE_CL, unlensedcl, cmb_unlensedcl_file, pxwin=0, beam=0, xhdr=xhdr
     END
     
     ELSE: BEGIN
     END

  ENDCASE

END