143 lines
3.8 KiB
Fortran
143 lines
3.8 KiB
Fortran
SUBROUTINE RYBHEQ
|
|
C =================
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
common/grdpra/GRD(MDEPTH),pra(mdepth),pgs0(mdepth),ANTP(MDEPTH)
|
|
common/rybpgs/CS(MDEPTH),PRAD2D(MDEPTH),F1HE
|
|
dimension pradfc(mdepth)
|
|
c
|
|
DO ID=1,ND
|
|
GRD(ID)=0.
|
|
pra(id)=0.
|
|
pradt(id)=0.
|
|
pradfc(id)=0.
|
|
END DO
|
|
PRD0=0.
|
|
c
|
|
if(ifprad.gt.0) then
|
|
CALL OPAINI(1)
|
|
DO IJ=1,NFREQ
|
|
CALL OPACF1(IJ)
|
|
CALL RTEFR1(IJ)
|
|
if(.not.lskip(1,ij))
|
|
* GRD(1)=GRD(1)+W(IJ)*ABSO1(1)*FH(IJ)*RAD1(1)
|
|
DO ID=2,ND
|
|
if(.not.lskip(id,ij)) then
|
|
GRD(ID)=GRD(ID)+(RAD1(ID)*FAK1(ID)-
|
|
* RAD1(ID-1)*FAK1(ID-1))*W(IJ)
|
|
pra(id)=pra(id)+RAD1(ID)*FAK1(ID)*W(IJ)
|
|
end if
|
|
END DO
|
|
pra(1)=pra(1)+RAD1(1)*FAK1(1)*W(IJ)-ABSO1(1)*W(IJ)*
|
|
* (RAD1(1)*FH(IJ)-HEXTRD(IJ))
|
|
END DO
|
|
if(idisk.eq.0) then
|
|
GRD(1)=PCK*GRD(1)/DENS(1)
|
|
else
|
|
grd(1)=pck*grd(2)
|
|
end if
|
|
PRD0=PRD0*DENS(1)*DM(1)*PCK
|
|
c
|
|
do id=1,nd
|
|
pra(id)=pra(id)*pck
|
|
pradt(id)=pra(id)
|
|
pradfc(id)=pra(id)/(2.5213e-15*temp(id)**4)
|
|
end do
|
|
C
|
|
do id=2,nd-1
|
|
dmm=un/(dm(id)-dm(id-1))
|
|
dmp=un/(dm(id+1)-dm(id))
|
|
dm0=two/(dm(id+1)-dm(id-1))
|
|
qq=((pra(id+1)-pra(id))*dmp-(pra(id)-pra(id-1))*dmm)*dm0
|
|
prad2d(id)=qq
|
|
end do
|
|
prad2d(1)=prad2d(2)
|
|
prad2d(nd)=0.
|
|
end if
|
|
c
|
|
if(idisk.eq.0) then
|
|
PGS0(1)=DM(1)*(GRAV-GRD(1))
|
|
DO ID=2,ND
|
|
PGS0(ID)=PGS0(ID-1)-PCK*GRD(ID)+GRAV*(DM(ID)-DM(ID-1))
|
|
END DO
|
|
else
|
|
pgs0(1)=pgs(1)
|
|
do id=2,nd
|
|
grv=(dm(id)-dm(id-1))*qgrav*(zd(id)+zd(id-1))*half
|
|
pgs0(id)=pgs0(id-1)-pck*grd(id)+grv
|
|
end do
|
|
end if
|
|
if(iprybh.gt.0) write(6,603) iter
|
|
603 format(/' rybheq in iter',i4,
|
|
* ' dm,dens,pgs0,pra,pradfc,grd,grad,ggrav'/)
|
|
c
|
|
if(iprybh.gt.0) then
|
|
do id=1,nd
|
|
if(id.eq.1) then
|
|
gra=grd(id)
|
|
ggr=qgrav*zd(id)
|
|
else
|
|
gra=pck*grd(id)/(dm(id)-dm(id-1))
|
|
ggr=qgrav*(zd(id)+zd(id-1))*half
|
|
end if
|
|
if(idisk.eq.0) then
|
|
write(6,604) id,dm(id),dens(id),pgs0(id),pra(id),
|
|
* pradfc(id),grd(id),gra
|
|
else
|
|
write(6,604) id,dm(id),dens(id),pgs0(id),pra(id),
|
|
* pradfc(id),grd(id),gra,ggr,zd(id)
|
|
end if
|
|
end do
|
|
end if
|
|
604 format(i4,1p9e13.4)
|
|
C
|
|
if(idisk.eq.0) then
|
|
DO ID=1,ND
|
|
T=TEMP(ID)
|
|
pgs(id)=pgs0(id)
|
|
AN=PGS(ID)/BOLK/T
|
|
CALL ELDENS(ID,T,AN,ANE,ENRG,ENTT,WM,1)
|
|
RHO=WMM(ID)*(AN-ANE)
|
|
DENS(ID)=RHO
|
|
ELEC(ID)=ANE
|
|
CALL WNSTOR(ID)
|
|
CALL STEQEQ(ID,POP,1)
|
|
c IF(.NOT.LCHC.and.iter.lt.ielcor) CALL ELCOR(ID)
|
|
END DO
|
|
else
|
|
c
|
|
c --------------------------
|
|
c the rest is only for disks
|
|
c --------------------------
|
|
c
|
|
cprad=2.5213e-15
|
|
do id=1,nd
|
|
cs(id)=pgs0(id)/dens(id)/temp(id)
|
|
c pi(id)=cprad*pradfc(id)
|
|
end do
|
|
c
|
|
hr1=grd(1)/qgrav
|
|
hg1=sqrt(two*cs(1)*temp(1)/qgrav)
|
|
x=(zd(1)-hr1)/hg1
|
|
IF(X.LT.3.) THEN
|
|
IF(X.LT.0.) X=0.
|
|
F1HE=8.86226925D-1*EXP(X*X)*ERFCX(X)
|
|
ELSE
|
|
F1HE=HALF*(UN-HALF/X/X)/X
|
|
END IF
|
|
c
|
|
ntemp=2
|
|
call pgset(ntemp)
|
|
c
|
|
do id=1,nd
|
|
pgs(id)=pgs0(id)
|
|
ptotal(id)=pgs0(id)+pra(id)
|
|
end do
|
|
c
|
|
END IF
|
|
C
|
|
RETURN
|
|
END
|