      subroutine facHYS
c
c  Calculate new solutions for 
c   Herd-year-seasons,  Random factor
c
      include 'TThead.f'
c
      real*8 wv(nog+10),rhs(no),diags(nog),dif,c,y,z,w
      real*8 ry(ntr),ay(ntr),ssh(nog),p(no,no)
      real*8 vnois(no),hh(nog),dx,chi,x,q(no),u
      integer iv(nog+10),levs(nrec),ndf,noz

      levs = 0
      mfac = 0
      ssh=0.d0
      noz=0
c      print *,' TTHYS '
c////////////////////////////////////////////////
      do 10 krec=1,mrec
         ihys = hys(krec)
         if(ihys.gt.mfac)mfac=ihys
         levs(ihys)=levs(ihys) + 1
 10   continue
c////////////////////////////////////////////////
      kend = 0
      mlev=0
      do 11 Lhys=1,mfac
         jrec = levs(Lhys)
         rhs = 0.d0
         diags = 0.d0
         if(jrec.lt.1)go to 11
         kstr = kend + 1
         kend = kend + jrec
      do 12 lrec=kstr,kend
         krec = phys(lrec)
c
         icod = kcod(krec)
         iam = an(krec)
         iage = bage(krec)
         ihys=hys(krec)
         iymj=kymj(krec)
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         ay=0.d0
         do 20 ktr=1,no
            y = obs(krec,ktr)
            if(y.gt.0.d0)then
            y=y-sage(iage,ktr)-sanm(iam,ktr)-symj(iymj,ktr)
            y=y-sape(iam,ktr)
            ay(ktr) = y
            endif
 20     continue
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
        ry=0.d0
        do 30 ir=1,no
         z = 0.d0
         do 35 ic=1,no
            m=ihmssf(ir,ic,no)
            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,no
         rhs(ir)=rhs(ir)+ry(ir)
c
         do 41 ic=ir,no
            m=ihmssf(ir,ic,no)
            diags(m)=diags(m)+ric(icod,m)
 41   continue
 40   continue    
 12   continue
c
c  solve, invert each diagonal block
c
 100  mlev=mlev + 1
       do 201 ir=1,no
          do 201 ic=ir,no
             mom=ihmssf(ir,ic,no)
          diags(mom)=diags(mom) + Hi(mom)
 201   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
          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
       
       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 - shys(Lhys,ir)
         ccn = ccn + (dif*dif)
         ccd = ccd + (z*z)
         shys(Lhys,ir) = z 
 102  continue
c
      if(itest.gt.0.and.Lhys.eq.3)then
      print 2221,Lhys
      print 2223,(shys(Lhys,j),j=1,ntr)
 2221 format(1x,'HYS',i6)
 2223 format(1x,10f8.4)
      endif
c      if(levs(Lhys).lt.5)go to 11
      if(bayes.gt.0)then
      do 210 ir=1,no
         c=shys(Lhys,ir)
         do 210 ic=ir,no
            m= ihmssf(ir,ic,no)
            ssh(m)=ssh(m)+c*shys(Lhys,ic)
 210     continue
      noz = noz + 1   
      endif
 11   continue
      if(bayes.gt.0)then
         Hi=0.d0
         ndf=noz+2
         call fgchi1(ndf,chi)
         dif = 1.d0/chi
         m=0
         do 218 kr=1,no
          do 218 kc=kr,no
             m=m+1
            Hi(m)= ssh(m)*dif
 218     continue
         nm=3
         do 219 ir=1,no
            hh(ir)=Hi(ihmssf(ir,ir,no))
 219     continue
         if(itest.gt.0)then
         print 2221,ndf
         print 2227,(hh(k),k=1,no)
 2227    format(1x,5f12.4)
         endif
         write(19)iter,nm,Hi
         call dkmvhf(Hi,no,wv,iv)
      endif
      return
      end
