!
!     quadratic interpolation subroutine , output is alpgt-->interpolated lnP
!
       SUBROUTINE quadInterpLnPr(h1,h2,h3,alp1,alp2,alp3,h,alpgt)
 
       implicit none
       REAL::h1,h2,h3,h,alpgt,alp1,alp2,alp3
!!!DRAGAN        REAL*8::h3m2,h3m1,a,b,c,d,dlt
      REAL::h3m2,h3m1,a,b,c,d,dlt
!       logical debug

         dlt=(alp1-alp3)*(alp2-alp1)*(alp3-alp2)
         h3m2=h3-h2
         h3m1=h3-h1
         b=(h3m2*(alp3*alp3-alp1*alp1)   &
            -h3m1*(alp3*alp3-alp2*alp2))/dlt
         c=(h3m1*(alp3-alp2)-h3m2*(alp3-alp1))/dlt
        
	
	
	 d=h1-c*alp1*alp1-b*alp1-h


         if (abs(c).le.0.01) then
            print *,"h3m1= ",h3m1," alp3=", alp3
            print *,"alp2= ",alp2," h3m2=", h3m2
            print *,"alp1= ",alp1," dlt =", dlt
            print *,'(h3m1*(alp3-alp2)-h3m2*(alp3-alp1)), dlt = ', &
                     (h3m1*(alp3-alp2)-h3m2*(alp3-alp1)), dlt 
            print *,'c = ', (h3m1*(alp3-alp2)-h3m2*(alp3-alp1))/dlt
            print *,'h3m1,alp3-alp2=',h3m1,alp3-alp2,h3m1*(alp3-alp2)
            print *,'h3m2,alp3-alp1=',h3m2,alp3-alp1,h3m2*(alp3-alp1)
         c=0.02   !!!!to avoid FPE - DRAGAN
	 endif
        
	
	 d=h1-c*alp1*alp1-b*alp1-h
	
	
	
	 if((b*b-4*c*d).ge.0) then
           alpgt=(-b-sqrt(b*b-4*c*d))/2/c
         else
           a=(h1-h2)/(alp1-alp2)
           b=h1-a*alp1
           alpgt=(h-b)/a
         endif

         alpgt = alp1*(h-h2)*(h-h3)/((h1-h2)*(h1-h3)) &
          + alp2*(h-h1)*(h-h3)/((h2-h1)*(h2-h3)) &
               + alp3*(h-h1)*(h-h2)/((h3-h1)*(h3-h2))

      END SUBROUTINE quadInterpLnPr