204 lines
6.7 KiB
Fortran
204 lines
6.7 KiB
Fortran
SUBROUTINE RYBENE
|
|
C =================
|
|
C
|
|
C complementing partial matrices of complete linearization
|
|
c corresponding to the contribution of ALI frequencies and the
|
|
c convection to the radiative/convective equilibrium equation in
|
|
c the Rybicki formalism
|
|
c
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
INCLUDE 'ARRAY1.FOR'
|
|
COMMON/CUBCON/ACNV,BCNV,DEL,GRDADB,DELMDE,RHO,FLXTOT,GRAVD
|
|
COMMON/RYBMTX/RA(MDEPTH),RB(MDEPTH),RC(MDEPTH),VR(MDEPTH),
|
|
* UA(MDEPTH),UB(MDEPTH),UC(MDEPTH),
|
|
* VA(MDEPTH),VB(MDEPTH),VC(MDEPTH),WR(MDEPTH),
|
|
* WM(MDEPTH,MDEPTH)
|
|
common/deridt/dert
|
|
C
|
|
C contribution from the radiative equilibrium part of the energy eq.
|
|
c
|
|
flxto0=sig4p*teff**4
|
|
DO ID=1,ND
|
|
WR(ID)=WR(ID)+FCOOL(ID)
|
|
IF(IDISK.EQ.1) WR(ID)=WR(ID)-reint(id)*TVISC(ID)
|
|
if(reint(id).gt.0.) then
|
|
IF(ID.GT.1) WM(ID,ID-1)
|
|
* =WM(ID,ID-1)+AREIT(ID)*dens(id)*reint(id)
|
|
WM(ID,ID) =WM(ID,ID) + REIT(ID)*dens(id)*reint(id)
|
|
IF(IDISK.EQ.1) WM(ID,ID)=WM(ID,ID)+
|
|
* DTVIST(ID)*reint(id)
|
|
IF(ID.LT.ND) WM(ID,ID+1)
|
|
* =WM(ID,ID+1)+CREIT(ID)*dens(id)*reint(id)
|
|
end if
|
|
c
|
|
if(idisk.eq.0) then
|
|
flxtot=flxto0
|
|
gravd=grav
|
|
else
|
|
flxtot=flxto0*(1.d0-thetav(id))
|
|
gravd=qgrav*zd(id)
|
|
end if
|
|
c
|
|
if(redif(id).gt.0) then
|
|
WR(ID)=WR(ID)+flxtot*redif(id)
|
|
IF(ID.GT.1) WM(ID,ID-1)=WM(ID,ID-1)+REDTM(ID)*REDIF(ID)
|
|
WM(ID,ID)= WM(ID,ID)+ REDT (ID)*REDIF(ID)
|
|
IF(ID.LT.ND) WM(ID,ID+1)=WM(ID,ID+1)+REDTP(ID)*REDIF(ID)
|
|
end if
|
|
END DO
|
|
c
|
|
C contribution from convection
|
|
C
|
|
IF(HMIX0.LE.0.OR.ICONV.LE.0) RETURN
|
|
if(dert.eq.0.) dert=0.01
|
|
DO 10 ID=idconz,ND
|
|
T=TEMP(ID)
|
|
P=PTOTAL(ID)
|
|
TM=TEMP(ID-1)
|
|
PM=PTOTAL(ID-1)
|
|
pr0=0.
|
|
if(icentr.eq.0.or.id.eq.nd) then
|
|
T0=SQRT(T*TM)
|
|
P0=SQRT(P*PM)
|
|
AB0=SQRT(ABROSD(ID)*ABROSD(ID-1))
|
|
pr0=sqrt(pradt(id)*pradt(id-1))*pck
|
|
DLP=UN/LOG(P/PM)
|
|
DLT=LOG(T/TM)*DLP
|
|
DELTA(ID)=DLT
|
|
DDT0= DLP/T
|
|
DDTM=-DLP/TM
|
|
CALL CONVEC(ID,T0,P0,P0-pr0,pr0,AB0,DLT,FLXCNV,VCON)
|
|
FLXC(ID)=FLXCNV
|
|
DHCDD=0.
|
|
IF(DELMDE.GT.0.) THEN
|
|
bb=bcnv*half
|
|
cord=un-bb/sqrt(bb*bb+(dlt-grdadb))
|
|
DHCDD=1.5D0/DELMDE*FLXCNV*cord
|
|
END IF
|
|
if(id.lt.icbegp-2) then
|
|
flxc(id)=0.
|
|
go to 10
|
|
end if
|
|
T1=(un+dert)*T0
|
|
CALL CONVEC(ID,T1,P0,p0-pr0,pr0,AB0,DLT,FLXC1,VCON)
|
|
DHCDT0=(FLXC1-FLXCNV)*HALF/dert
|
|
DHCDT =DHCDT0/T+DHCDD*DDT0
|
|
DHCDTM=DHCDT0/TM+DHCDD*DDTM
|
|
else
|
|
dlm=log(p/pm)
|
|
dlp=log(ptotal(id+1)/p)
|
|
dl0=dlm+dlp
|
|
palf=dlp/dlm/dl0
|
|
pgam=dlm/dlp/dl0
|
|
dlt=palf*log(t/tm)+pgam*log(temp(id+1)/t)
|
|
delta(id)=dlt
|
|
ddtm=-palf/tm
|
|
ddtp=pgam/temp(id+1)
|
|
ddt0=(palf-pgam)/t
|
|
pr0=sqrt(pradt(id)*pradt(id-1))*pck
|
|
call convec(id,t,p,p-pr0,pr0,abrosd(id),dlt,flxcnv,vcon)
|
|
flxc(id)=flxcnv
|
|
dhcdd=0.
|
|
IF(DELMDE.GT.0.) then
|
|
bb=bcnv*half
|
|
cord=un-bb/sqrt(bb*bb+(dlt-grdadb))
|
|
DHCDD=1.5D0/DELMDE*FLXCNV*cord
|
|
end if
|
|
if(id.lt.icbegp-2) then
|
|
flxc(id)=0.
|
|
go to 10
|
|
end if
|
|
T1=(un+dert)*T
|
|
CALL CONVEC(ID,T1,P,p-pr0,pr0,abrosd(id),DLT,FLXC1,VCON)
|
|
DHCDT0=(FLXC1-FLXCNV)*HALF/dert
|
|
DHCDT =DHCDT0/T+DHCDD*DDT0
|
|
DHCDTM=DHCDD*DDTM
|
|
DHCDTP=DHCDD*DDTP
|
|
end if
|
|
|
|
C
|
|
C additional terms in matrix WM
|
|
C
|
|
C ** 1. differential equation form
|
|
C
|
|
if(redif(id).gt.0) then
|
|
WM(ID,ID-1)=WM(ID,ID-1)+DHCDTM*redif(id)
|
|
WM(ID,ID)= WM(ID,ID)+DHCDT*redif(id)
|
|
WR(ID)=WR(ID)-FLXC(ID)*redif(id)
|
|
if(icentr.ne.0.and.id.lt.nd) then
|
|
WM(ID,ID+1)=WM(ID,ID+1)+DHCDTP*redif(id)
|
|
end if
|
|
END IF
|
|
C
|
|
C ** 2. integral equation form
|
|
C
|
|
if(reint(id).gt.0.) then
|
|
if(id.lt.nd) then
|
|
TP=TEMP(ID+1)
|
|
PTP=PTOTAL(ID+1)
|
|
T0=SQRT(T*TP)
|
|
P0=SQRT(P*PTP)
|
|
pr0=sqrt(pradt(id)*pradt(id-1))*pck
|
|
AB0=SQRT(ABROSD(ID)*ABROSD(ID+1))
|
|
DLP=UN/LOG(PTP/P)
|
|
DLT=LOG(TP/T)*DLP
|
|
DDTP0= DLP/TP
|
|
DDTPM=-DLP/T
|
|
CALL CONVEC(ID,T0,P0,p0-pr0,pr0,AB0,DLT,FLXCNV,VCON)
|
|
DHCDDP=0.
|
|
IF(DELMDE.GT.0.) DHCDDP=1.5D0/DELMDE*FLXCNV
|
|
T1=1.001D0*T0
|
|
CALL CONVEC(ID,T1,P0,p0-pr0,pr0,AB0,DLT,FLXC1,VCON)
|
|
DHCDT0=(FLXC1-FLXCNV)*500.
|
|
DHCDTP=DHCDT0/TP+DHCDDP*DDTP0
|
|
DHCDTU=DHCDT0/T+DHCDDP*DDTPM
|
|
C
|
|
C additional terms in matrices A and B (the row corresponding to
|
|
C energy balance, i.e. T) due to convection
|
|
C
|
|
DELM=(DM(ID+1)-DM(ID-1))*HALF
|
|
RDELM=DENS(ID)/DELM*reint(id)
|
|
IF(ICONV.GT.0) THEN
|
|
WM(ID,ID-1)=WM(ID,ID-1)-RDELM*DHCDTM
|
|
WM(ID,ID)=WM(ID,ID)+RDELM*(DHCDTP-DHCDTU)
|
|
WM(ID,ID+1)=WM(ID,ID+1)+RDELM*DHCDTP
|
|
END IF
|
|
WR(ID)=WR(ID)-RDELM*(FLXCNV-FLXC(ID))
|
|
ELSE
|
|
TP=TEMP(ID-2)
|
|
PTP=PTOTAL(ID-2)
|
|
T0=SQRT(T*TP)
|
|
P0=SQRT(P*PTP)
|
|
pr0=sqrt(pradt(id)*pradt(id-1))*pck
|
|
AB0=SQRT(ABROSD(ID)*ABROSD(ID-2))
|
|
DLP=UN/LOG(PTP/P)
|
|
DLT=LOG(TP/T)*DLP
|
|
DDTP0= DLP/TP
|
|
DDTPM=-DLP/T
|
|
CALL CONVEC(ID,T0,P0,p0-pr0,pr0,AB0,DLT,FLXCNV,VCON)
|
|
if(flxcnv.le.0.) go to 10
|
|
DHCDDP=0.
|
|
IF(DELMDE.GT.0.) DHCDDP=1.5D0/DELMDE*FLXCNV
|
|
T1=1.001D0*T0
|
|
CALL CONVEC(ID,T1,P0,p0-pr0,pr0,AB0,DLT,FLXC1,VCON)
|
|
DHCDTP=(FLXC1-FLXCNV)*1.D3/T0*HALF+DHCDDP*DDTP0
|
|
C
|
|
C additional terms in matrices A and B (the row corresponding to
|
|
C energy balance, i.e. T) due to convection
|
|
C
|
|
DELM=(DM(ID)-DM(ID-2))*HALF
|
|
RDELM=DENS(ID)/DELM*reint(id)
|
|
DELHC=WMM(ID)/DELM*(FLXCNV-FLXC(ID))
|
|
WM(ID,ID-1)=WM(ID,ID-1)+RDELM*(DHCDT-DHCDTP)
|
|
WM(ID,ID)=WM(ID,ID)+RDELM*DHCDT
|
|
WR(ID)=WR(ID)-RDELM*(FLXC(ID)-FLXCNV)
|
|
END IF
|
|
END IF
|
|
10 CONTINUE
|
|
c
|
|
RETURN
|
|
END
|