      subroutine facAGE
c
c  Calculate new solutions for 
c   Age at Calving
c
      include 'TThead.f'
c
      real*8 wv(nog+10),rhs(no),diags(nog),dif,c,y,z,w,x,dx
      real*8 ry(ntr),ay(ntr),vnois(no),hh(no)
      real*8 p(no,no),q(no),u
      integer iv(nog+10),levs(nrec)

      mfac = 0
      levs=0
c      print *,' TTAge '
c////////////////////////////////////////////////
      do 10 krec=1,mrec
        iage = bage(krec)
        if(iage.gt.mfac)mfac=iage
        levs(iage)=levs(iage)+1
 10   continue  
c////////////////////////////////////////////////
      mlev = 0
      kend=0
      do 11 lage=1,mfac
         jrec = levs(lage)
         if(jrec.lt.1)go to 11
c         print *,' Age ',lage,jrec
         rhs=0.d0
         diags=0.d0
         kstr = kend + 1
         kend = kend + jrec
      do 12 lrec=kstr,kend
        krec = page(lrec)
c
        icod = kcod(krec)
         iam  = an(krec)
         ihys = hys(krec)
         iage = bage(krec)
         iymj = kymj(krec)
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         ay=0.d0
         do 20 ktr=1,ntr
            y=obs(krec,ktr)
            if(y.gt.0.d0)then
            y = y - shys(ihys,ktr)-sanm(iam,ktr)-sape(iam,ktr)
            y = y - symj(iymj,ktr)
            ay(ktr) = y
         endif
 20     continue
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
        ry=0.d0
        do 30 ir=1,ntr
         z = 0.d0
         do 35 ic=1,ntr
            m=ihmssf(ir,ic,ntr)
            z = z + ric(icod,m)*ay(ic)
 35      continue
         ry(ir) = z
 30      continue
c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
c
c  form  X'R^-1X  and X'R^-1y
c
      do 40 ir=1,ntr
         w = ry(ir)
        rhs(ir)=rhs(ir)+w
        do 41 ic=ir,ntr
           m=ihmssf(ir,ic,ntr)
        diags(m)=diags(m)+ric(icod,m)
 41   continue
 40   continue    
 12    continue
c
c  solve, invert diagonal block
c
 100   mlev = mlev + 1
c       do 119 ir=1,no
c          mc=ihmssf(ir,ir,no)
c          diags(mc)=diags(mc)+0.4d0
c 119  continue    
      call dkmvhf(diags,no,wv,iv)
c
      vnois = 0.d0
      if(bayes.gt.0)then
        P=0.d0
        hh=0.d0
        do 89 i=1,no
           call fgnor3(u)
           hh(i)=u
           x=diags(ihmssf(i,i,no))
           if(x.gt.0.d0)then
              im=i-1
              if(i.gt.1)then
                 do 90 j=1,im
                    x=x - P(j,i)*P(j,i)
 90              continue
              endif
           else
              x=0.d0
           endif
           z=0.d0
           if(x.gt.0.d0)then
              P(i,i)=dsqrt(x)
              z=1.d0/P(i,i)
           endif
           ip=i+1
           if(ip.gt.no)go to 89
           do 81 j=ip,no
              x=diags(ihmssf(j,i,no))
              if(i.gt.1)then
                 do 92 k=1,im
                    x=x - P(j,k)*P(i,k)
 92              continue
              endif
              P(j,i)=x*z
 81        continue
           
 89     continue
c
         do 93 i=1,no
            x=0.d0
            do 94 j=1,i
               x=x +P(i,j)*hh(j)
 94         continue
            vnois(i)=x
 93      continue
         
      endif
c
      do 102 ir=1,no
         z=0.d0
         do 140 ic=1,no
            m=ihmssf(ir,ic,no)
            z=z+diags(m)*rhs(ic)
 140     continue
         z = z + vnois(ir)
         dif = z - sage(lage,ir)
         ccn = ccn + dif*dif
         ccd = ccd + z*z
         sage(lage,ir) = z 
 102  continue
c
      if(itest.gt.0)then
      if(lage.eq.3)then
       print 2873,lage
       print 2874,(sage(lage,k),k=1,no)
 2873  format(1x,'AGE ',i6)
 2874  format(1x,10f8.3)
      endif
      endif
 11   continue
      return
      end
