       subroutine csvals(ri,zi,ci,cc,si,cri,czi,sri,szi,
     & crri,czzi,crzi)

c***********************************************************************
c      Sound speed and slowness interpolation subroutine
c      Written by Orlando Camargo Rodriguez
c      orodrig@ualg.pt
c      Universidade do Algarve
c      Physics Department
c      Signal Processing Laboratory
c      Faro, 19/08/2009 at 18:30
c****&******************************************************************

c      Include global (common) variables:

       include 'global.for'

c      Define local variables:

       character*4 cclss

       real*8 zi,ri,ci,si,cc
       real*8 cri,czi,sri,szi
       real*8 crri,czzi,crzi
       real*8 drdn,dzdn
       real*8 k,root,root32,root52,a
       real*8 bmunk,bmunk2,epsiln,eta

c-----------------------------------------------------------------------

       cclss = cclass(1:4) 

       epsiln = 7.4e-3;
       bmunk  = 1300.0
       bmunk2 = bmunk*bmunk

c***********************************************************************
c     Sound speed profiles:
c***********************************************************************
       if (cdist.eq.'c(z,z)') then

     	  cri = 0.0
	 crri = 0.0 
	 crzi = 0.0

c=======================================================================
c       Analytical sound speed profiles:
c=======================================================================
	    if (cclss.eq.'ISOV') then

	         ci = c0(1)
	        czi = 0.0 
	       czzi = 0.0 

	    elseif (cclss.eq.'LINP') then

	        k = (c0(2)-c0(1))/(z0(2)-z0(1))
	       ci = c0(1) + k*( zi - z0(1) )
	      czi = k
	     czzi = 0.0

	    elseif (cclss.eq.'PARP') then

	        k = (c0(2)-c0(1))/( (z0(2)-z0(1))**2.0 )
	       ci = c0(1) + k*( zi - z0(1) )**2.0
	      czi = 2*k*( zi - z0(1) )
	     czzi = 2*k	

	    elseif (cclss.eq.'EXPP') then

	        k = log( c0(1)/c0(2) )/( z0(2) - z0(1) )
	       ci = c0(1)*exp( -k*( zi - z0(1) ) )
	      czi = -k*ci
	     czzi = k*k*ci

	    elseif (cclss.eq.'N2LP') then

	        k = ( ( c0(1)/c0(2) )**2.0 - 1 )/( z0(2) - z0(1) )
	     root = sqrt( 1 + k*( zi - z0(1) ) )
	   root32 = root**( 3.0/2.0 )
	   root52 = root**( 5.0/2.0 )  
	       ci = c0(1)/sqrt( 1 + k*( zi - z0(1) ) )
	      czi =    -k*c0(1)/( 2*root32 )
	     czzi = 3*k*k*c0(1)/( 4*root52 )

	    elseif (cclss.eq.'ISQP') then

	        a = ( c0(2)/c0(1) -1 )**2.0
	     root = sqrt( a/(1-a) )
	        k = root/( z0(2) - z0(1) )
	     root = sqrt( 1 + ( k*( zi - z0(1) ) )**(2.0) )
	   root32 = root**( 3.0/2.0 )
	   root52 = root**( 5.0/2.0 )
	       ci = c0(1)*( 1 + k*( zi - z0(1) )/root )
	      czi = c0(1)*k/root32
	     czzi = -3*c0(1)*k**3.0/root52

	    elseif (cclss.eq.'MUNK') then 

	    eta = 2*( zi - z0(1) )/bmunk
	     ci = c0(1)*( 1 + epsiln*( eta + exp(-eta) - 1 ) )
	    czi = 2*epsiln*c0(1)*( 1 - exp(-eta) )/bmunk
	   czzi = 4*epsiln*c0(1)*exp( -eta )/bmunk2

c=======================================================================
c       Tabulated sound speed profile:
c=======================================================================
	    elseif (cclss.eq.'TABL') then

	    call cvals1(nz0,z0,c0,zi,ci,czi,czzi)

	    else

            write(6,*) 'Unknown sound speed profile'
	    write(6,*) 'aborting calculations...'
            stop

	    end if

c***********************************************************************
c     Sound speed field:
c***********************************************************************

       elseif (cdist.eq.'c(r,z)') then

       call cvals2(nr0,nz0,r0,z0,c02d,ri,zi,ci,cri,czi,crri,czzi,crzi)

c-----------------------------------------------------------------------
       else

           write(6,*) 'Unknown sound speed distribution'
	   write(6,*) 'aborting calculations...'
           stop

       end if

        cc = ci*ci

        si =  1.0/ci 
       sri = -cri/cc
       szi = -czi/cc


c***********************************************************************
c     Back to main:
c*********************************************************************** 	

       return

       end 
