       subroutine calpav(ctitle,nthtas)
       
c***********************************************************************
c      CALPAV.FOR
c      CALculation of coherent acoustic Pressure And particle Velocity
c      Faro, Tue Jan  7 21:29:40 WET 2014
c      Written by Orlando Camargo Rodriguez
c****&******************************************************************
c      Note: 
c      readin.for will stop calculations if the array type is incorrect.
c      Therefore, after testing different array types the remaining one 
c      can only be the one available from the list. 
c****&******************************************************************

c      Include global (common) variables:

       include 'global.for'

c      Define local variables:

       character(len=128) header       
       character(len=60 ) ctitle
       character(len=8  ) v1

       integer*4 irefl(np),jbdry(np)
       integer*4 i,imax,j,k,nthtas
              
       real*8 tbdry(2,np)
       real*8 xp(3)
       real*8 es(2),e1(2),deltar(2)
       
       real*8 thetai,ctheta
       real*8 rh,zh,n,nxn,dzdr,qray
       real*8 sray,dr,dz,q0,w
       real*8 drarray,dzarray
       
       complex*16 dpdr2d(nhyd2,nhyd2),dpdz2d(nhyd2,nhyd2)
       complex*16   pl2d(nhyd2,nhyd2),  pc2d(nhyd2,nhyd2)
       complex*16   pr2d(nhyd2,nhyd2)
       complex*16   pu2d(nhyd2,nhyd2),  pd2d(nhyd2,nhyd2)
       
       complex*16  paux(npress)
       complex*16   dpdr(np),  dpdz(np)
       complex*16 pressl(np),pressc(np),pressr(np),pressu(np),pressd(np)
       complex*16  decay(np)
       complex*16 zp(3)
       
       complex*16 aray,pc,pl,pr,pu,pd
       complex*16 dpdri,dpdzi
       complex*16 dummi

c***********************************************************************

c***********************************************************************

       header = 'Signal Processing LABoratory: TRACEO
     & acoustic ray tracing program - Pressure And Velocity'
     
      header(117:124) = char(ichar(' '))
      header(125:126) = char(0) // char(1)
       
      write(header(127:128),'(a)') int(19785,kind=2)
      
      open(unit=outfil,file='pav.mat',form='unformatted',
     & access='stream')
     
c     Open the file and write the header:

      do i = 1, 16
         j = 8*(i-1) + 1
         write(outfil) header(j:j+7)
      end do

c     Write the angles:

      v1 = 'ray_elev'
	   
      call dtelrv(outfil,1,nthtas,v1,thetas)

c***********************************************************************
c      Get the right array dr and dz:
c***********************************************************************

       call csvals(rx,zx,cx,cc,sigmai,cri,czi,sri,szi,
     & crri,czzi,crzi)

       omega  = 2.0*pi*freqx
       lambda =     cx/freqx

       dzarray = lambda/10.0
       drarray = lambda/10.0

       if (artype.eq.'HRY') then

           do j = 2,nra
           drarray = min( abs( rarray(j) - rarray(j-1) ) ,
     &     drarray )
           end do
	   
       elseif (artype.eq.'VRY')  then

	   do k = 2,nza
	   dzarray = min( abs( zarray(k) - zarray(k-1) ) ,
     &	   dzarray )
	   end do
	   
       elseif (artype.eq.'LRY')  then

           do j = 2,nra
           drarray = min( abs( rarray(j) - rarray(j-1) ) ,
     &     drarray )
           dzarray = min( abs( zarray(j) - zarray(j-1) ) ,
     &     dzarray )
           end do

       else

           do j = 2,nra
           drarray = min( abs( rarray(j) - rarray(j-1) ) ,
     &     drarray )
           end do

           do k = 2,nza
           dzarray = min( abs( zarray(k) - zarray(k-1) ) ,
     &     dzarray )
           end do

       end if

c***********************************************************************
c      Solve the EIKonal and the DYNamic sets of EQuations:
c***********************************************************************

        do i = 1,nthtas

           thetai = -thetas(i)*pi/180.0

           ctheta = abs( cos( thetai ) ) 

c          Trace a ray as long as it is neither 90 or -90:
           
           if (ctheta.gt.1.0e-7) then
  
           call seikeq(thetai,imax,irefl,decay,jbdry,tbdry)
           call sdyneq(thetai,imax,irefl,decay,jbdry,tbdry)
	   
	   q0 = cx/( pi*dtheta/180.0 )
	   
c          Now that the ray has been calculated let's determine the
c          ray influence at each point of the array:

           if (artype.eq.'HRY') then

	       zh = zarray(1)

	       do j = 1,nra

                  rh = rarray(j)

                  if ((rh.ge.rmin).and.(rh.lt.rmax)) then

		     if (iretrn.eq.0) then
	
 		     call star(imax,rh,zh,q0,irefl,drarray,dzarray,
     &		     pl,pc,pr,pu,pd)

 		     pressl(j) = pressl(j) + pl
 		     pressc(j) = pressc(j) + pc
 		     pressr(j) = pressr(j) + pr
 		     pressu(j) = pressu(j) + pu
 		     pressd(j) = pressd(j) + pd
		     
		     else
		      
		     call mstar(imax,rh,zh,q0,irefl,drarray,dzarray,
     &		     pl,pc,pr,pu,pd)

 		     pressl(j) = pressl(j) + pl
 		     pressc(j) = pressc(j) + pc
 		     pressr(j) = pressr(j) + pr
 		     pressu(j) = pressu(j) + pu
 		     pressd(j) = pressd(j) + pd

		     end if

		  end if

               end do

	   elseif (artype.eq.'VRY')  then
               
	       rh = rarray(1)
	       
               if ((rh.ge.rmin).and.(rh.lt.rmax)) then

	         if (iretrn.eq.0) then

		 do j = 1,nza
	
 		    zh = zarray(j)

 		    call star(imax,rh,zh,q0,irefl,drarray,dzarray,
     &		    pl,pc,pr,pu,pd)

 		    pressl(j) = pressl(j) + pl
 		    pressc(j) = pressc(j) + pc
 		    pressr(j) = pressr(j) + pr
 		    pressu(j) = pressu(j) + pu
 		    pressd(j) = pressd(j) + pd

 		 end do

		 else

		 do j = 1,nza

 		    zh = zarray(j)

 		    call mstar(imax,rh,zh,q0,irefl,drarray,dzarray,
     &		    pl,pc,pr,pu,pd)

 		    pressl(j) = pressl(j) + pl
 		    pressc(j) = pressc(j) + pc
 		    pressr(j) = pressr(j) + pr
 		    pressu(j) = pressu(j) + pu
 		    pressd(j) = pressd(j) + pd

 		 end do

		 end if

	       end if

           elseif (artype.eq.'LRY')  then

               do j = 1,nra

	          rh = rarray(j)
                  zh = zarray(j)

		  if ((rh.ge.rmin).and.(rh.lt.rmax)) then

		     if (iretrn.eq.0) then

 		     call star(imax,rh,zh,q0,irefl,drarray,dzarray,
     &		     pl,pc,pr,pu,pd)

 		     pressl(j) = pressl(j) + pl
 		     pressc(j) = pressc(j) + pc
 		     pressr(j) = pressr(j) + pr
 		     pressu(j) = pressu(j) + pu
 		     pressd(j) = pressd(j) + pd

		     else

		     call mstar(imax,rh,zh,q0,irefl,drarray,dzarray,
     &		     pl,pc,pr,pu,pd)

 		     pressl(j) = pressl(j) + pl
 		     pressc(j) = pressc(j) + pc
 		     pressr(j) = pressr(j) + pr
 		     pressu(j) = pressu(j) + pu
 		     pressd(j) = pressd(j) + pd

		     end if

		  end if

	       end do

	   else

c              Yes: this is slow, but the code is short!

	       do j = 1,nra

                  rh = rarray(j)

                  if ((rh.ge.rmin).and.(rh.lt.rmax)) then
		  
                     if (iretrn.eq.0) then

		     do k = 1,nza

 		     	zh = zarray(k)

 		     	call star(imax,rh,zh,q0,irefl,drarray,dzarray,
     &		     	pl,pc,pr,pu,pd)

 		     	pl2d(k,j) = pl2d(k,j) + pl
 		     	pc2d(k,j) = pc2d(k,j) + pc
 		     	pr2d(k,j) = pr2d(k,j) + pr
 		     	pu2d(k,j) = pu2d(k,j) + pu
 		     	pd2d(k,j) = pd2d(k,j) + pd

 		     end do

		     else

		     do k = 1,nza

 		     	zh = zarray(k)

		     	call mstar(imax,rh,zh,q0,irefl,drarray,dzarray,
     &		     	pl,pc,pr,pu,pd)
 
 		     	pl2d(k,j) = pl2d(k,j) + pl
 		     	pc2d(k,j) = pc2d(k,j) + pc
 		     	pr2d(k,j) = pr2d(k,j) + pr
 		     	pu2d(k,j) = pu2d(k,j) + pu
 		     	pd2d(k,j) = pd2d(k,j) + pd

 		     end do
		     
		     end if

		  end if

               end do

           end if
	   
	   end if

       end do
       
c***********************************************************************
c      Now that we have the acoustic pressure at the star we can 
c      calculate the particle velocity components:
c***********************************************************************

       if (artype.eq.'HRY') then
       
       zh = zarray(1)
               
	  do j = 1,nra

 	     rh = rarray(j)
             
	     xp(1) = rh - drarray
	     xp(2) = rh
	     xp(3) = rh + drarray
	     
	     zp(1) = pressl(j)
	     zp(2) = pressc(j)
	     zp(3) = pressr(j)
	     
	     call cbpa1d(xp,zp,rh,dummi,dpdri,dummi)
	     
	     dpdr(j) = -imunit*dpdri
	     
	     xp(1) = zh - dzarray
	     xp(2) = zh
	     xp(3) = zh + dzarray
	     
	     zp(1) = pressd(j)
	     zp(3) = pressu(j)
	     
	     call cbpa1d(xp,zp,zh,dummi,dpdzi,dummi)
	     
	     dpdz(j) = -imunit*dpdzi
	     
 	  end do
	  
       elseif (artype.eq.'VRY') then
        
       rh = rarray(1)
	       
	  do j = 1,nza
	
 	  zh = zarray(j)
	
	  xp(1) = rh - drarray
	  xp(2) = rh
	  xp(3) = rh + drarray
	
	  zp(1) = pressl(j)
	  zp(2) = pressc(j)
	  zp(3) = pressr(j)
	
	  call cbpa1d(xp,zp,rh,dummi,dpdri,dummi)
	
	  dpdr(j) = -imunit*dpdri
	
	  xp(1) = zh - dzarray
	  xp(2) = zh
	  xp(3) = zh + dzarray
	
	  zp(1) = pressd(j)
	  zp(3) = pressu(j)
	
	  call cbpa1d(xp,zp,zh,dummi,dpdzi,dummi)
	
	  dpdz(j) = -imunit*dpdzi

 	  end do
	  
       elseif (artype.eq.'LRY')  then
	   
 	  do j = 1,nra
 
	     rh = rarray(j)
 	     zh = zarray(j)
	     
	     xp(1) = rh - drarray
	     xp(2) = rh
	     xp(3) = rh + drarray
	     
	     zp(1) = pressl(j)
	     zp(2) = pressc(j)
	     zp(3) = pressr(j)
	     
	     call cbpa1d(xp,zp,rh,dummi,dpdri,dummi)
	     
	     dpdr(j) = -imunit*dpdri
	     
	     xp(1) = zh - dzarray
	     xp(2) = zh
	     xp(3) = zh + dzarray
	     
	     zp(1) = pressd(j)
	     zp(3) = pressu(j)
	     
	     call cbpa1d(xp,zp,zh,dummi,dpdzi,dummi)
	     
	     dpdz(j) = -imunit*dpdzi

	  end do

       else
	  
	  do j = 1,nra
	
 	     rh = rarray(j)

 	     do k = 1,nza

 	  	zh = zarray(k)
		
		xp(1) = rh - drarray
		xp(2) = rh
		xp(3) = rh + drarray
	
		zp(1) = pl2d(k,j)
		zp(2) = pc2d(k,j)
		zp(3) = pr2d(k,j)
		
		call cbpa1d(xp,zp,rh,dummi,dpdri,dummi)
	
		dpdr2d(k,j) = -imunit*dpdri
	
		xp(1) = zh - dzarray
		xp(2) = zh
		xp(3) = zh + dzarray
		
		zp(1) = pd2d(k,j)
		zp(3) = pu2d(k,j)
			
		call cbpa1d(xp,zp,zh,dummi,dpdzi,dummi)
	
		dpdz2d(k,j) = -imunit*dpdzi
	
 	     end do

 	  end do

       end if
	
c***********************************************************************
c      Write the data to the output file:
c***********************************************************************

        if (artype.ne.'RRY') then
	
	    jj = max(nra,nza)
	    
            v1 = 'pressure'
	    
	    call dtelcm(outfil,jj,1,v1,pressc)
	    
	    v1 = 'gradpr_r'
	    
	    call dtelcm(outfil,jj,1,v1,dpdr)
	    
	    v1 = 'gradpr_z'
	    
            call dtelcm(outfil,jj,1,v1,dpdz)
		    
	else
	
c         Yes, this is not elegant, but it works!!!!!

            jj = 0

	    do j = 1,nra
	    
	       do i = 1,nza
	       
	       jj = jj + 1
	       
	       paux(jj) = pc2d(i,j)
	       
	       end do
	    
	    end do
	    
	    v1 = 'pressure'
	    
            call dtelcm(outfil,nza,nra,v1,paux)
	    
            jj = 0

	    do j = 1,nra
	    
	       do i = 1,nza
	       
	       jj = jj + 1
	       
	       paux(jj) = dpdr2d(i,j)
	       
	       end do
	    
	    end do
	    
	    v1 = 'gradpr_r'
	    
	    call dtelcm(outfil,nza,nra,v1,paux)
	    
            jj = 0

	    do j = 1,nra
	    
	       do i = 1,nza
	       
	       jj = jj + 1
	       
	       paux(jj) = dpdz2d(i,j)
	       
	       end do
	    
	    end do
	    
	    v1 = 'gradpr_z'
	    
	    call dtelcm(outfil,nza,nra,v1,paux)

	end if

c***********************************************************************
c      The Happy End:
c***********************************************************************

       v1 = 'rarray2d'
	   
       call dtelrv(outfil,1,nra,v1,rarray)

       v1 = 'zarray2d'
	   
       call dtelrv(outfil,1,nza,v1,zarray)

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

       return

       end
