subroutine ccontr (icall,c1,c2,c3,ier1,ier2)
      parameter (lth=40,lvr=30,no=50)
      common /rocm0/ theta (lth)
      common /rocm4/ y
      double precision c1,c2,c3,theta,y,w,one,two
      dimension c2(*),c3(lvr,*)
      data one,two/1.,2./
      if (icall.le.1) return
      w=y
      
         y=(y**theta(3)-one)/theta(3)
         
      
      
      call cels (c1,c2,c3,ier1,ier2)
      y=w
      c1=c1-two*(theta(3)-one)*log(y)
      
      return
      end