PRO dtemplbl, fastlayers=fastlayers, nlayers=nlayers, pause = pause ;, tof = toff

	; Purturbs the temperature profile of a synthetic spectrum layer by layer,
	; & builds up the Kb matrix. Starts from Jim's deltaTemp program, but adapted
	; to look at the spectrum instead of the vmr.
	; RB July 2010

	; assumes we're in a working rdrv retrieval directory

	IF( NOT KEYWORD_SET( fastlayers )) THEN BEGIN
    READ, fastlayers, PROMPT=' Enter number of fastcode layers (TAB=43): '
  ENDIF
  IF( NOT KEYWORD_SET( nlayers )) THEN BEGIN
    READ, nlayers, PROMPT=' Enter number of retrieval layers (TAB=47): '
  ENDIF

	FORWARD_FUNCTION readstat, readctl, setfilenames, readspec, mksyn, readt15


	; Set up file names from LAPACK version sfit2 run
	rc = setfilenames( fn )

  fpert=1 ;force the temperature perturbation to be 1 degree


	; Use default ILS for all calcs
	cmd = 'cp ' + fn.ilsfit + ' ' + fn.ilsnam
	SPAWN, cmd
	; Use default cinput for all calcs
	cmd = 'cp ' + fn.cindef + ' ' + fn.cinnam
	SPAWN, cmd
	; Use default fasc.inp for all calcs
	cmd = 'cp ' + fn.fasdef + ' ' + fn.fasnam
	SPAWN, cmd
	; Use default fasc.prf for all calcs
	cmd = 'cp ' + fn.prfdef + ' ' + fn.prfnam
	SPAWN, cmd
	; Set default rdrv.ctl at start
	cmd = 'cp ' + fn.rdfdef + ' ' + fn.rdfnam
	SPAWN, cmd

	; Define codes to do the work
	fast = '/p/ya4/bin/fast2.k -ffasc -v394'
	sfit = '/p/ya4/bin/sfit2.k -s0 -v394lp'
	;rdrv = '/p/ya4/sfitsrc/394lp/rdrv39'

	; STEP 1

	; Make the default spectrum & do a retrieval
	; Make the dummy spectrum
	mksyn, file = fn.t15one

	; Create a forward model rdrv.ctl file from the existing one
	rc = 0
	rc = readctl( ctl, fn.rdfnam )

	PRINT, 'Overwriting : ', fn.rdfnam
	OPENW, lun, fn.rdfnam, /GET_LUN
		PRINTF, lun, ctl.binf
		PRINTF, lun, ctl.ptfl
		PRINTF, lun, ctl.msfl
		PRINTF, lun, ctl.mxfl
		PRINTF, lun, fn.t15one
		PRINTF, lun, '.FALSE.'
		PRINTF, lun, fn.t15syn
		PRINTF, lun, '.FALSE.'
		PRINTF, lun, '.FALSE.'
	FREE_LUN, lun

	SPAWN, fast
	SPAWN, sfit

	junk=''
	IF( KEYWORD_SET( pause )) THEN BEGIN
		READ, junk, PROMPT=' Finished default spectrum - Hit enter to continue or q(uit): '
		IF( STRUPCASE( junk ) EQ 'Q' )THEN STOP, '   ...Quit by command.'
	ENDIF ELSE PRINT,' Finished default spectrum..'

	; Compute a retrieval from the default spectrum
	; Make a retrieval rdrv.ctl
	PRINT, 'Overwriting : ', fn.rdfnam
	OPENW, lun, fn.rdfnam, /GET_LUN
		PRINTF, lun, ctl.binf
		PRINTF, lun, ctl.ptfl
		PRINTF, lun, ctl.msfl
		PRINTF, lun, ctl.mxfl
		PRINTF, lun, fn.t15syn
		PRINTF, lun, '.TRUE.'
		PRINTF, lun, '.TRUE.'
		PRINTF, lun, 'K.out'
		PRINTF, lun, '.TRUE.'
		PRINTF, lun, 'Sa.complete'
	FREE_LUN, lun

	; Compute the default retrieval, vmr, column ...
	SPAWN, sfit

	; Save the necessary files
	;cmd = 'cp PRFS.out PRFS.out.def'
	;SPAWN, cmd
	;cmd = 'cp statevec statevec.def'
	;SPAWN, cmd
	;cmd = 'cp cdetail cdetail.def'
	;SPAWN, cmd
  ;cmd = 'cp pbpfile pbpfile.def'
  ;SPAWN, cmd

  unpertpbp=readspec(fn.pbfile)

  ; Read in Statevector file & get neutral vmr
  rc = 0
  rc = readstat( stdef, fn.stfile )

 ;set up the Kb matrix
 temp=size(unpertpbp)
 npoints=temp(1)

 Kb=fltarr(fastlayers, npoints); we're going to fill this column by column


	IF( KEYWORD_SET( pause )) THEN BEGIN
		READ, junk, PROMPT=' Finished retrieval of default - Hit enter to continue or q(uit): '
		IF( STRUPCASE( junk ) EQ 'Q' )THEN STOP, '   ...Quit by command.'
	ENDIF ELSE PRINT, ' Finished retrieval of default... going into loop through layers'


	; STEP 2: Do a temperature perturbation of each layer, get the new spectrum and make the column of Kb


;put a loop around this whole thing for the number of layers that we have
for mmm=0, fastlayers-1 do begin
  print, 'doing layer '+string(mmm)+' of '+string(fastlayers)

  ; Use default ILS for all calcs
  cmd = 'cp ' + fn.ilsfit + ' ' + fn.ilsnam
  SPAWN, cmd
  ; Use default cinput for all calcs
  cmd = 'cp ' + fn.cindef + ' ' + fn.cinnam
  SPAWN, cmd
  ; Use default fasc.inp for all calcs
  cmd = 'cp ' + fn.fasdef + ' ' + fn.fasnam
  SPAWN, cmd
  ; Use default fasc.prf for all calcs
  cmd = 'cp ' + fn.prfdef + ' ' + fn.prfnam
  SPAWN, cmd
  ; Set default rdrv.ctl at start
  cmd = 'cp ' + fn.rdfdef + ' ' + fn.rdfnam
  SPAWN, cmd

  ; Perturb the temperature profile and make a synthetic spectrum

	; Make a forward model rdrv.ctl
	PRINT, 'Overwriting : ', fn.rdfnam
	OPENW, lun, fn.rdfnam, /GET_LUN
		PRINTF, lun, ctl.binf
		PRINTF, lun, ctl.ptfl
		PRINTF, lun, ctl.msfl
		PRINTF, lun, ctl.mxfl
		PRINTF, lun, fn.t15one
		PRINTF, lun, '.FALSE.'
		PRINTF, lun, fn.t15syn
		PRINTF, lun, '.FALSE.'
		PRINTF, lun, '.FALSE.'
	FREE_LUN, lun

	; Read in the neutral temperature profile
	profs = 70
	lys   = 0

	refmodrd, ref, lys, fn.prfdef, profs
	neutT = ref[2].prf[*]

	; Perturb temperature profile
	ref[2].prf[mmm] = neutT(mmm) + fpert

	; Write out the purturbed fasc.prf file
	refmodwr, ref, lys, fn.prfnam, profs
	; Compute spectrum from the atmosphere
	SPAWN, fast

   ; sometimes the sza changes after going through this process. this corrects for that
	openr,1,'fasc.ms'
	readf,1,junk
	close,1

	tempvalue=strsplit(junk,/extract)
	newsza=fix(tempvalue(0))
	sedcommand='sed s/NNNNN/'+tempvalue(0)+'/g cinput.def.NNNNN > cinput'
	spawn, sedcommand

	SPAWN, sfit

  newpbp=readspec(fn.pbfile)

	IF( KEYWORD_SET( pause )) THEN BEGIN
		READ, junk, PROMPT='Finished perturbed spectrum -  Hit enter to continue or q(uit): '
		IF( STRUPCASE( junk ) EQ 'Q' )THEN STOP, '   ...Quit by command.'
	ENDIF ELSE PRINT, 'Finished perturbed spectrum...'+string(mmm)


	; Read in Statevector file
	rc = 0
	rc = readstat( stat, fn.stfile )

	IF( rc NE 0 ) THEN BEGIN
		PRINTF, -2,'Could not read st file: ', fn.stfile
		STOP
	ENDIF

;calculate Kb : change in spectrum/change in Temperature (old-new)/perturbation

dif = unpertpbp - newpbp

Kb(mmm,*)=dif/fpert

endfor

;print out the file. This will be used in errorcalcncar.
openw,1,'Kb_temp.cov'
for i=0, npoints-1 do printf,1,Kb(*,i)
close,1

end