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 ************************************************************************