
	program corasc

c   ===================================================================
c		Stephane ERARD, octobre 1998


c	Corrects saturation of some channels in the Ascraeus session.
c	The new values are computed from the correlation observed with 
c	the neighbouring channels at lower flux.
c	Translated from IDL to Fortran, to preserve consistency of 
c	files format.

c   ===================================================================


	real*4 	temp, temp2,limi(12)
	integer*2 	matrice(120,25,0:63), atm(120,25,0:63),
     &			i,j,k,spec(72),spec1(72),
     &			sto(120,25,8),nbf
	character*50 ch


c	Reads calibrated data
	ch='$1$DUA1:[PCCOMMON.BASEISM.CAL_DATA.MARS]'
	open(1,file=ch//'asceven.atm',form='unformatted',recl=36,
     &access='direct',status='old')
	open(2,file=ch//'asceven.ca0',form='unformatted',recl=36,
     &access='direct',status='old')

	nbf=0
	coef=32768.*2.
	do i=1,120
	 do j=1,25
	  nbf=nbf+1
	  read(2'nbf)spec		!even channels
	  read(1'nbf)spec1		! same with atm correction

	  do k=1,8		!Saves instrumental parameters
	   sto(i,j,k)=spec(k)
	  enddo
	  do k=1,64		! On stocke les pairs seuls en ordre
	   matrice(i,j,k-1)=spec(k+8)	!/coef
	   atm(i,j,k-1)=spec1(k+8)	!/coef
	  enddo

	 enddo
	enddo
	type*,'nb of raw spectra read',nbf

	close(2)
	close(1)


c     ---------- Interpolate saturated values ------------

	limi(4)=.31*coef		! channel 10
	limi(6)=.30*coef		! channel 14
	limi(8)=.305*coef		! channel 18
	limi(10)=.34*coef		! channel 22
	limi(11)=.35*coef		! channel 23

	do i=1,120
	 do j=1,25

c 	Channel 10
	  if (matrice(i,j,4).gt.limi(4)) then 
		temp=matrice(i,j,3)/2.+matrice(i,j,5)/2.
		temp2=temp * .992950 - 0.000143766*coef
		matrice(i,j,4) = int(temp2)
		atm(i,j,4) = int(temp2)
	  endif

c	Channel 14
	  if (matrice(i,j,6).gt.limi(6)) then 
		temp=matrice(i,j,5)/2.+matrice(i,j,7)/2.
		temp2=temp * 1.00783 - 0.001288*coef
		matrice(i,j,6) = int(temp2)
		atm(i,j,6) = int(temp2)
	  endif

c	Channel 18
	  if (matrice(i,j,8).gt.limi(8)) then 
		temp=matrice(i,j,7)/2.+matrice(i,j,9)/2.
		temp2=temp * 1.00614 - 0.00073248*coef
		matrice(i,j,8) = int(temp2)
		atm(i,j,8) = int(temp2)
	  endif

c	Channel 22
	  if (matrice(i,j,10).gt.limi(10)) then 
		temp=matrice(i,j,9)
		temp2=temp * 0.998776 + 0.00164536*coef
		matrice(i,j,10) = int(temp2)
		atm(i,j,10) = int(temp2)
	  endif

c	Channel 23
	  if (matrice(i,j,11).gt.limi(11)) then 
		temp=matrice(i,j,12)
		temp2=temp * 1.0025 - 0.00157386*coef
		matrice(i,j,11) = int(temp2)
		atm(i,j,11) = int(temp2)
	  endif

	 end do
	end do


c		----- Writes the corrected data -----

	nbf=0
	open(1,file='asceven.atm',form='unformatted',recl=36,
     &access='direct',status='new')  ! do not keep previous version
	open(2,file='asceven.cal',form='unformatted',recl=36,
     &access='direct',status='new')

	do i=1,120
	 do j=1,25
	  nbf=nbf+1

	  do k=1,8		!Save instrumental parameters
	   spec(k)=sto(i,j,k)
	   spec1(k)=sto(i,j,k)
	  enddo
	  do k=1,64
	   spec(k+8)=matrice(i,j,k-1)
	   spec1(k+8)=atm(i,j,k-1)
	  enddo

	  write(2'nbf)spec	!even channels
	  write(1'nbf)spec1	! atm correction

	 enddo
	enddo
	type*,'nb of registered spectra',nbf

	close(2)
	close(1)

5000	stop
      end
