subroutine cshell
*
*  cshell steps from the center shell to the next shell.
*  gshell steps from shell n to shell n+1
*  See Schwarzschild & Harm 1965, ApJ, 142, 855
*  for a partial explanation of variable meanings
*
      implicit double precision(a-h,o-z)

      common/coef/hrp(400),hrt(400),hra(400),hbp(400),hbt(400),hba(400),
     1            hpp(400),hpt(400),hpa(400),htp(400),htt(400),hta(400)
      common/cal/ut1,ot1,et1,op1,o1,ep1,qnp1,qnt1,p1,q1,f1,t1,w1,b1,e1,
     1 qe1,ea1,qn1,up1,qe2,qn2,qnp2,qnt2
      common/temp/s1,r1,s2,r2,b2,p2,t2,ea2,xc2,xo2,fca2,f2,q2,w2,c
      common/thermo/u2,up2,ut2,e2,ep2,et2,psi2,pg,o2,op2,ot2,fp,ft,
     1 en2(5),fn,fcc,ce,ci,cif,w,nu
      common/contrl/ds,g,sm,wc,it,nite,ja,jb,j,k,l
 
      data fnat/0.4342945/

      hrp(1)=-(up2)*.333333
      hrt(1)=-(ut2)*.333333
      hra(1)=-.333333*(u2+3.*r2-s2+.622089-sm)
      hpp(1)=0
      hpt(1)=0
      hpa(1)=0
      hbp(1)=(-qe2*ep2+qnp2)*fnat
      hbt(1)=-qe2*(e2-ea2+et2*fnat)+qnt2*fnat
      hba(1)=-b2+(-qe2*(e2-ea2)+qn2)*fnat
      htp(1)=0
      htt(1)=0
      hta(1)=0
      it=-1
      return
*......................................................................*

      entry gshell

      gs=2.*fnat/(s2-s1)
      ar=gs*(r1-r2)+fnat*(q1+q2)
      ap=gs*(p1-p2)-fnat*(f1+f2)
      gp=2.*fnat/(p1-p2)
      gt=gp*(t1-t2)/(p1-p2)
      at=gp*(t1-t2)-fnat*(w1+w2)
      ab=gs*(b2-b1)+fnat*(qe1*(e1-ea1)+qe2*(e2-ea2)-qn1-qn2)
      c1=q1*up1+hrp(k-1)*(-gs+3.*q1)
      c2=q1*ut1+hrt(k-1)*(-gs+3.*q1)
      c3=gs+3.*q2
      c4=q2*up2
      c5=q2*ut2
      car=ar+hra(k-1)*(gs-3.*q1)
      c6=-gs-f1-4.*f1*hrp(k-1)
      c7=-4.*f1*hrt(k-1)
      c8=-4.*f2
      c9=gs-f2
      cap=ap+4.*f1*hra(k-1)
      c10=gt+w1+fnat*w1*(op1/o1+hbp(k-1)/b1)
      c11=-gp-4.*w1+fnat*w1*(ot1/o1+hbt(k-1)/b1)
      c12=fnat*w2/b2
      c13=-gt+w2+w2*op2*fnat/o2
      c14=gp-4.*w2+w2*ot2*fnat/o2
      cat=at-fnat*w1*hba(k-1)/b1
      c15=gs*hbp(k-1)+fnat*(qnp1-qe1*ep1)
      c16=-qe1*(e1-ea1)+fnat*(qnt1-qe1*et1)+gs*hbt(k-1)
      c17=-gs
      c18=fnat*(qnp2-qe2*ep2) 
      c19=-qe2*(e2-ea2)+fnat*(qnt2-qe2*et2)
      cab=ab-gs*hba(k-1)
      c20=c1*c8-c3*c6
      c21=c2*c8-c3*c7
      c22=c4*c8-c3*c9
      c23=c5*c8
      carp=c8*car-c3*cap
      c24=c12*c15-c10*c17
      c25=c12*c16-c11*c17
      c26=c12*c18-c13*c17
      c27=c12*c19-c14*c17
      cabt=c12*cab-c17*cat
      c28=c20*c25-c21*c24
      c29=c22*c25-c21*c26
      c30=c23*c25-c21*c27
      cp=c25*carp-c21*cabt
      c32=c20*c26-c22*c24
      c33=c20*c27-c23*c24
      ct=c20*cabt-c24*carp
      cd1=-1./c28
      cd2=-1./c17
      cd3=-1./c3
      hpp(k)=c29*cd1
      hpt(k)=c30*cd1
      hpa(k)=-cp*cd1
      htp(k)=c32*cd1
      htt(k)=c33*cd1
      hta(k)=-ct*cd1
      hbp(k)=( c15*hpp(k)+c16*htp(k)+c18)*cd2
      hbt(k)=( c15*hpt(k)+c16*htt(k)+c19)*cd2
      hba(k)=( c15*hpa(k)+c16*hta(k)-cab)*cd2
      hrp(k)=( c1*hpp(k)+c2*htp(k)+c4)*cd3
      hrt(k)=( c1*hpt(k)+c2*htt(k)+c5)*cd3
      hra(k)=( c1*hpa(k)+c2*hta(k)-car)*cd3

      return
      end 

************************************************************************