PRO dspec, pause=pause

	; Runs one version of rdev39 to make a spectrum with a preturbed parameter
	; Then does a retrieval using the default version

   ; uflag
   ;  1 - change line strength
   ;  2 - change air broadended halfwidth
   ; purturbation value is a fractional scaling

	; Requires a valid t15 file...

	FORWARD_FUNCTION readstat, readctl, readprfs, setfilenames, mksyn, readt15

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

	uflag = 0
	ugas = 0
	uval = 0.0

;	READ, uflag, ugas, uval, PROMPT=' Enter parameter flag (1=linestrength, 2=air broadened HW, gas id and perturbation value (fractional) : '

;Rather than reading this in, set to do both and get the gas ID number from the files - assuming that we're looking at the gas of interest

rc = readstat( stdef, fn.stfile )
case stdef.gas(0) of
 'H2O': ugas=1
 'CO2': ugas=2
 'O3' : ugas=3
 'N2O': ugas=4
 'CO' : ugas=5
 'CH4': ugas=6
 'O2' : ugas=7
 'NO' : ugas=8
 'SO2' : ugas=9
 'NO2' : ugas=10
 'NH3' : ugas=11
 'HNO3' : ugas=12
 'OH' : ugas=13
 'HF' : ugas=14
 'HCL' : ugas=15
 'HBR' : ugas=16
 'HI' : ugas=17
 'CLO' : ugas=18
 'OCS' : ugas=19
 'H2CO' : ugas=20
 'HOCL' : ugas=21
 'HO2' : ugas=22
 'H2O2' : ugas=23
 'HONO' : ugas=24
 'HO2NO2' : ugas=25
 'N2O5' : ugas=26
 'CLONO2' : ugas=27
 'HCN' : ugas=28
 'CH3F' : ugas=29
 'CH3CL' : ugas=30
 'CF4' : ugas=31
 'CCL2F2' : ugas=32
 'CCL3F' : ugas=33
 'CH3CCL3' : ugas=34
 'CCL4' : ugas=35
 'COF2' : ugas=36
 'COCLF' : ugas=37
 'C2H6' : ugas=38
 'C2H4' : ugas=39
 'C2H2' : ugas=40
 'N2' : ugas=41
 'CHF2CL' : ugas=42
 'COCL2' : ugas=43
 'CH3BR' : ugas=44
 'CH3I' : ugas=45
 else:  begin
  print, 'Please enter the gas id number into the table in this code'
  stop
 endelse
 endcase

print, 'gas =', stdef.gas(0)
print, 'ugas =', ugas


	; 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 = '/project/ya4/bin/fastc.k fasc'
  sfit = '/project/ya4/bin/sfit2.k -i -v 394xo'
  rdrv = '/project/ya4/sfitsrc/394lp/rdrv39'

	; Step 1
	; Do a retrieval on a default synthetic spectrum and save the profile & column

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

	; Make the forward model rdrv.ctl
	mksyn, file = fn.t15one

	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

	; Compute the default airmass & spectrum
	SPAWN, fast
	SPAWN, rdrv

	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, rdrv

;	; 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

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

	;read in the unperturbed pbp file

  unpertpbp=readspec(fn.pbfile)

	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...'

	; Step 2
	; Make a perturbed spectrum & do a retrieval - 1st the linestrength
	PRINT, 'Doing line strength calculation'
	uflag=1
	uval=0.05 ;a 5% perturbation


	; Make the forward model rdrv.ctl - overwrite t15syn
	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.'		; retrieval flag
		PRINTF, lun, fn.t15syn
		PRINTF, lun, '.FALSE.'		; write k-mat flag
		PRINTF, lun, '.FALSE.'		; write Sa flag
		PRINTF, lun, '.FALSE.'		; write gas output files flag
		PRINTF, lun, uflag, ugas, uval
	FREE_LUN, lun

	; Make perturbed spectrum
	SPAWN, rdrv

	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...'

	;; 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
;
; 	; Do default retrieval on perturbed spectrum
;	SPAWN, rdrv

	; Read in Statevector file & get perturbed vmr/column
;	rc = 0
;	rc = readstat( stat, fn.stfile )

	;READ, fpert, PROMPT=' Enter fractional value of perturbation : '
	fpert = uval


;read in pbp
  newpbp=readspec(fn.pbfile)

 ;calculate Kb : change in spectrum/change in lineparam/perturbation

dif = unpertpbp - newpbp
;fpert=20 ;20% perturbation in the line shape

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

Kb=dif/fpert

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


; Step 2
  ; Make a perturbed spectrum & do a retrieval - now the air broadened half width
  PRINT, 'Doing air broadened half width calculations'
  uflag=2
  uval=0.05 ;a 5% perturbation


  ; Make the forward model rdrv.ctl - overwrite t15syn
  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.'    ; retrieval flag
    PRINTF, lun, fn.t15syn
    PRINTF, lun, '.FALSE.'    ; write k-mat flag
    PRINTF, lun, '.FALSE.'    ; write Sa flag
    PRINTF, lun, '.FALSE.'    ; write gas output files flag
    PRINTF, lun, uflag, ugas, uval
  FREE_LUN, lun

  ; Make perturbed spectrum
  SPAWN, rdrv

  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...'

;  ; 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

;  ; Do default retrieval on perturbed spectrum
;  SPAWN, rdrv

;  ; Read in Statevector file & get perturbed vmr/column
;  rc = 0
;  rc = readstat( stat, fn.stfile )

  fpert = uval

 ;read in pbp
  newpbp2=readspec(fn.pbfile)

 ;calculate Kb : change in spectrum/change in lineparam/perturbation

dif2 = unpertpbp - newpbp2

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

Kb2=dif2/fpert

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


END