      subroutine facAPE
c
c  Calculate new solutions for 
c   Animal Permanent Environment, random
c  Too many animals to store all diagonal blocks
c
      include 'TThead.f'
c
      real*8 ry(no),zrz(nog),dif,c,w,y,z,d,s,x
      real*8 wv(nog+10),ay(no),zpg,ssa(nog),chi
      real*8 acc(no),rhs(no),dg,dx,p(no,no),u
      real*8 vnois(no),hh(nog),qqq(no),q(no)
      integer iv(nog+10),jsr(nam),levs(nrec),kpg,ndf

      mlev=0
      levs = 0
      iprev=0
      ssa=0.d0
      noz=0
      qqq=0.d0
      spg=0.d0
      kpg=0
      zpg = 0.d0
      itdr = 0
      iprog = 0
      jsr=0
      mfac = 0
c      print *,' TTanm '
c  ////////////////////////////////////////////////
      do 10 krec=1,mrec
         iam = an(krec)
         if(iam.gt.mfac)mfac=iam
         levs(iam)=levs(iam)+1
 10   continue
      kend = 0
      mlev=maxm
      do 11 Liam = 2001,mfac
         jrec = levs(Liam)
         rhs=0.d0
         if(jrec.lt.1)go to 50
          kstr = kend + 1
          kend = kend + jrec
        do 12 lrec = kstr,kend
             krec = panm(lrec)
             icod = kcod(krec)
             ihys = hys(krec)
             iam = an(krec)
             iage = bage(krec)
             iymj = kymj(krec)
             itdr(iam) = itdr(iam)+1
c             if(Liam.ne.iam)print 1948,Liam,iam
c 1948  format(1x,'ANM ',4i10)
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             ay=0.d0
          do 20 ktr=1,ntr
             y = obs(krec,ktr)
             if(y.gt.0.d0)then
            y = y-sage(iage,ktr)-shys(ihys,ktr)
            y = y-sanm(iam,ktr) -symj(iymj,ktr)
            ay(ktr)=y
         endif
 20     continue
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
        ry=0.d0
        do 30 itr=1,ntr
         z = 0.d0
         do 35 itc=1,ntr
            m=ihmssf(itr,itc,ntr)
            z = z + ric(icod,m)*ay(itc)
 35      continue
         ry(itr) = z
 30      continue
c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
c
c  form  rhs
c
      do 40 itr=1,ntr
        rhs(itr)=rhs(itr)+ry(itr)
        do 741 itc=itr,ntr
           m=ihmssf(itr,itc,ntr)
           zrz(m)=zrz(m)+ric(icod,m)
 741     continue
 40    continue
 12    continue

c/////////////////////////////////////////////////////////////

c  ******************************************************
c  Add matrix adjustments, if any
c *******************************************************

c       print 1226,Liam,dg
c 1226  format(5x,'A-inv diag',i8,f15.6)
 50    do 830 kog=1,nog
          zrz(kog) = zrz(kog) + Pi(kog)
 830   continue
       call dkmvhf(zrz,ntr,wv,iv)
c
c  solve
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=zrz(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 = zrz(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 kr=1,no
         z=0.d0
         do 140 kc=1,no
            m=ihmssf(kr,kc,no)
            z=z + zrz(m)*rhs(kc)
 140     continue
         z = z + vnois(kr)
         dif = z - sape(Liam,kr)
         ccn = ccn + (dif*dif)
         ccd = ccd + (z*z)
         sape(Liam,kr) = z 
 102  continue
c
      if(itest.gt.0.and.Liam.eq.100000)then
        print 2586,Liam
       print 2587,(sape(Liam,L),L=1,no)
 2586 format(1x,'APE ',i6)
 2587 format(1x,10f8.3)
      endif
 11   continue
c
      if(bayes.gt.0)then
      ssa=0.d0
      noz=0
      do 350 Liam=2001,mfac
      if(itdr(Liam).gt.0)then
         m=0
         do 231 ir=1,no
            c=sape(Liam,ir)
            do 231 ic=ir,no
               m=m+1
               ssa(m)=ssa(m)+c*sape(Liam,ic)
 231    continue
           noz=noz+1
        endif
 350    continue

         Pi=0.d0
         ndf=noz+2
         call fgchi1(ndf,chi)
         dif=1.d0/chi

         do 238 m=1,nog
            Pi(m)=ssa(m)*dif
 238     continue
         do 239 i=1,no
            m=ihmssf(i,i,no)
            hh(i)=Pi(m)
 239     continue
         if(itest.gt.0)then
         print *,'APE ',ndf,chi
         print 3088,(hh(j),j=1,no)
 3088    format(1x,5f12.4)
         endif
         nm=2
         write(19)iter,nm,Pi
         call dkmvhf(Pi,no,wv,iv)
      endif
c
      return
      end
