329 lines
9.3 KiB
Fortran
329 lines
9.3 KiB
Fortran
SUBROUTINE CONREF
|
|
C =================
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ARRAY1.FOR'
|
|
COMMON/CUBCON/ACNV,BCNV,DEL,GRDADB,DELMDE,RHO,FLXTOT,GRAVD
|
|
common/imucnn/imucon
|
|
dimension idcon(mdepth),
|
|
* flxtt(mdepth),delta0(mdepth)
|
|
C
|
|
IF(ICONV.LE.0.AND.INDL.EQ.0) RETURN
|
|
FLXTO0=SIG4P*TEFF**4
|
|
DLTND=DELTA(ND)
|
|
IFNDM1=0
|
|
C
|
|
ICBEG=0
|
|
ICEND=0
|
|
DO ID=2,ND
|
|
T=TEMP(ID)
|
|
P=PTOTAL(ID)
|
|
PG=PGS(ID)
|
|
PRAD=P-PG-HALF*DENS(ID)*VTURB(ID)**2
|
|
TM=TEMP(ID-1)
|
|
PM=PTOTAL(ID-1)
|
|
PGM=PGS(ID-1)
|
|
PRADM=PM-PGM-HALF*DENS(ID-1)*VTURB(ID-1)**2
|
|
IF(ID.EQ.ND.AND.IFNDM1.EQ.1) THEN
|
|
FAC=DLTND*(P-PM)/(P+PM)
|
|
T=TM*(UN+FAC)/(UN-FAC)
|
|
END IF
|
|
if(ilgder.eq.0) then
|
|
T0=HALF*(T+TM)
|
|
P0=HALF*(P+PM)
|
|
PG0=HALF*(PG+PGM)
|
|
PR0=HALF*(PRAD+PRADM)
|
|
AB0=HALF*(ABROSD(ID)+ABROSD(ID-1))
|
|
DLT=(T-TM)/(P-PM)*P0/T0
|
|
else
|
|
T0=SQRT(T*TM)
|
|
P0=SQRT(P*PM)
|
|
PG0=SQRT(PG*PGM)
|
|
PR0=SQRT(PRAD*PRADM)
|
|
AB0=SQRT(ABROSD(ID)*ABROSD(ID-1))
|
|
DLT=LOG(T/TM)/LOG(P/PM)
|
|
end if
|
|
DELTA(ID)=DLT
|
|
if(idisk.eq.0) then
|
|
flxtot=flxto0
|
|
gravd=grav
|
|
else
|
|
flxtot=flxto0*(1.d0-thetav(id))
|
|
gravd=qgrav*zd(id)
|
|
end if
|
|
flxtt(id)=flxtot
|
|
C
|
|
C convective flux
|
|
C
|
|
CALL CONVEC(ID,T0,P0,PG0,PR0,AB0,DLT,FLXCNV,VCON)
|
|
FLXC(ID)=FLXCNV
|
|
idcon(id)=0
|
|
if(flxc(id).gt.0.) idcon(id)=1
|
|
IF(ICBEG.EQ.0.AND.FLXC(ID).GT.0..AND.FLXC(ID-1).EQ.0..
|
|
* AND.ID.GT.IDCONZ) ICBEG=ID
|
|
if(icbeg.gt.0.and.flxc(id).gt.0.) icend=id
|
|
END DO
|
|
C
|
|
c correction algorithm - if in the convetion zone one has
|
|
c some depth point have gradient slightly smaller that the
|
|
c adaibatic gradient, recompute gradients to have convection there.
|
|
C in this case, the radiation flux is held fixed
|
|
c
|
|
c
|
|
icbeg0=icbeg
|
|
if(ideepc.gt.0) then
|
|
if(icend.eq.nd-1.and.ideepc.eq.2) icend=nd
|
|
icbegd=icend
|
|
do id=icend,icbeg,-1
|
|
if(idcon(id).gt.0) then
|
|
icbegd=id
|
|
else
|
|
igap=0
|
|
do idd=id-1,id-ndcgap,-1
|
|
if(idcon(idd).gt.0) igap=1
|
|
end do
|
|
if(igap.gt.0) then
|
|
icbegd=id
|
|
else
|
|
go to 10
|
|
end if
|
|
end if
|
|
end do
|
|
10 continue
|
|
icbeg0=icbegd
|
|
end if
|
|
c
|
|
if(ideepc.eq.3) icend=nd
|
|
if(ideepc.ge.4) then
|
|
if(icend.le.icbegp) then
|
|
icbeg0=icbegp
|
|
icend=icendp
|
|
end if
|
|
end if
|
|
icbegp=icbeg0
|
|
icendp=icend
|
|
c
|
|
if(icbeg0.gt.0.and.icend.gt.0) then
|
|
write(6,601) icbeg0,icend
|
|
601 format(/' convective refinement between depths ',2i4/)
|
|
c
|
|
c check the temperature at the depth just above the base of
|
|
c convection zone; if there is an oscillation there, re-adjust
|
|
c the temperature
|
|
c
|
|
if(temp(icbeg0-1).lt.temp(icbeg0-2).and.
|
|
* temp(icbeg0-1).lt.temp(icbeg0)) then
|
|
temp(icbeg0-1)=half*(temp(icbeg0)+temp(icbeg0-2))
|
|
end if
|
|
c
|
|
do id=icbeg0,icend
|
|
T=TEMP(ID)
|
|
P=PTOTAL(ID)
|
|
PG=PGS(ID)
|
|
PRAD=P-PG-HALF*DENS(ID)*VTURB(ID)**2
|
|
TM=TEMP(ID-1)
|
|
PM=PTOTAL(ID-1)
|
|
PGM=PGS(ID-1)
|
|
PRADM=PM-PGM-HALF*DENS(ID-1)*VTURB(ID-1)**2
|
|
if(ilgder.eq.0) then
|
|
P0=HALF*(P+PM)
|
|
T0=HALF*(T+TM)
|
|
PG0=HALF*(PG+PGM)
|
|
PR0=HALF*(PRAD+PRADM)
|
|
AB0=HALF*(ABROSD(ID)+ABROSD(ID-1))
|
|
DLT=(T-TM)/(P-PM)*P0/T0
|
|
ppd=(p+pm)/(p-pm)
|
|
else
|
|
T0=SQRT(T*TM)
|
|
P0=SQRT(P*PM)
|
|
PG0=SQRT(PG*PGM)
|
|
PR0=SQRT(PRAD*PRADM)
|
|
AB0=SQRT(ABROSD(ID)*ABROSD(ID-1))
|
|
DLT=LOG(T/TM)/LOG(P/PM)
|
|
end if
|
|
if(idisk.eq.0) then
|
|
flxtot=flxto0
|
|
gravd=grav
|
|
else
|
|
flxtot=flxto0*(1.d0-thetav(id))
|
|
gravd=qgrav*zd(id)
|
|
end if
|
|
fcnv=flxtot-flrd(id)
|
|
tor=t
|
|
if(fcnv.le.0.) fcnv=flxtot*(un-flrd(id)/flrd(2))
|
|
if(fcnv.lt.0.) fcnv=0.001*flxtot
|
|
c
|
|
c iteration loop to correct temperature
|
|
c
|
|
iic=0
|
|
20 iic=iic+1
|
|
if(ilgder.eq.0) then
|
|
T0=HALF*(T+TM)
|
|
AB0=HALF*(ABROSD(ID)+ABROSD(ID-1))
|
|
DLT=(T-TM)/(P-PM)*P0/T0
|
|
else
|
|
T0=SQRT(T*TM)
|
|
AB0=SQRT(ABROSD(ID)*ABROSD(ID-1))
|
|
DLT=LOG(T/TM)/LOG(P/PM)
|
|
end if
|
|
DELTA(ID)=DLT
|
|
if(flxc(id)/flxtot.gt.crflim) then
|
|
CALL CONVC1(ID,T0,P0,PG0,PR0,AB0,DLT,FLXCNV,FC0)
|
|
deltae=(fcnv/fc0)**0.666666666666667
|
|
deltaa=deltae+bcnv*sqrt(deltae)
|
|
dlt=deltaa+grdadb
|
|
end if
|
|
told=t
|
|
if(ilgder.eq.0) then
|
|
dlp=dlt/ppd
|
|
t=tm*(un+dlp)/(un-dlp)
|
|
else
|
|
t=tm*(P/PM)**DLT
|
|
end if
|
|
flxcnv=fc0*deltae**1.5
|
|
ff=flxcnv/flxtot
|
|
dtt=(t-told)/told
|
|
erfl=(flrd(id)+flxcnv)/flxtot
|
|
if(ilgder.eq.0) then
|
|
T0=HALF*(T+TM)
|
|
DLT=(T-TM)/(P-PM)*P0/T0
|
|
else
|
|
T0=SQRT(T*TM)
|
|
DLT=LOG(T/TM)/LOG(P/PM)
|
|
end if
|
|
CALL CONVEC(ID,T0,P0,PG0,PR0,AB0,DLT,FLXCN0,VCON)
|
|
erfl=(flrd(id)+flxcn0)/flxtot
|
|
if(abs(dtt).gt.1.e-9.and.iic.lt.10) go to 20
|
|
delta(id)=dlt
|
|
temp(id)=t
|
|
c write(6,667) id,iic,tor,t,flxcn0/flxtot,dlt,grdadb
|
|
end do
|
|
c
|
|
c new refinement procedure
|
|
c
|
|
if(iter.ge.imucon) then
|
|
icbeg0=icbeg
|
|
icend=nd
|
|
icendp=nd
|
|
write(6,674) imucon,icbeg0,icend
|
|
674 format(/' modification with imucon: icbeg0,icend',3i4)
|
|
write(6,677) icbeg0,icend
|
|
677 format(/' new refinement procedure: icbeg0, icend ',2i4/
|
|
*' entries are: id,itrc,t,dlt,grdadb,flrd/ft,flr/ft,fcn0/ft,
|
|
*(fcn0+flr)/ft'/)
|
|
do id=icbeg0,icend
|
|
t=temp(id)
|
|
p=ptotal(id)
|
|
tm=temp(id-1)
|
|
pm=ptotal(id-1)
|
|
pg=pgs(id)
|
|
PG=PGS(ID)
|
|
PRAD=P-PG-HALF*DENS(ID)*VTURB(ID)**2
|
|
PGM=PGS(ID)
|
|
PRADM=PM-PGM-HALF*DENS(ID-1)*VTURB(ID-1)**2
|
|
pg0=sqrt(pg*pgm)
|
|
pr0=sqrt(prad*pradm)
|
|
told=t
|
|
fcnv=flxtt(id)-flrd(id)
|
|
t0=sqrt(t*tm)
|
|
p0=sqrt(p*pm)
|
|
ab0=sqrt(abrosd(id)*abrosd(id-1))
|
|
dlt=log(t/tm)/log(p/pm)
|
|
call convc1(id,t0,p0,pg0,prad0,ab0,dlt,flxcn0,fc0)
|
|
alp=min(flrd(id),flxtt(id))/t0**4/dlt
|
|
c
|
|
if(fcnv.le.0.) go to 200
|
|
if(flxcn0.le.0.) go to 200
|
|
c
|
|
bet=bcnv/t0**3
|
|
deltae=(fcnv/fc0)**twothr
|
|
deltaa=deltae+bcnv*sqrt(deltae)
|
|
dlt=deltaa+grdadb
|
|
t=tm*(P/PM)**DLT
|
|
t0=sqrt(t*tm)
|
|
c
|
|
itrnrc=0
|
|
100 continue
|
|
itrnrc=itrnrc+1
|
|
t1=un/t
|
|
dltp=t1/log(p/pm)
|
|
t0p=half*t0*t1
|
|
dele=(flxtt(id)-alp*t0**4*dlt)/fc0
|
|
dele3= dele**third
|
|
delep=-alp*t0**4*dlt/fc0*(two*t1+d1tp/dlt)
|
|
vl=dlt-grdadb-dele3*(dele3+bet*t0**3)
|
|
bb=dltp-twothr*delep/dele3-
|
|
* bet*t0**3*dele3*(1.5d0*t1+third*delep/dele)
|
|
dt=-vl/bb*t1
|
|
t=t*(un+dt)
|
|
t0=sqrt(t*tm)
|
|
dlt=log(t/tm)/log(p/pm)
|
|
call convc1(id,t0,p0,pg0,prad0,ab0,dlt,flxcn0,fc0)
|
|
645 format(2i4,1pe11.3,0pf8.2,1p3e13.5)
|
|
if(abs(dt).lt.1.e-9.or.itrnrc.gt.20) go to 110
|
|
go to 100
|
|
110 continue
|
|
go to 230
|
|
c
|
|
200 continue
|
|
if(flxtt(id).lt.flrd(id)) then
|
|
alp=flxtt(id)/flrd(id)*t0**4*delta0(id)
|
|
itrnrc=0
|
|
210 continue
|
|
itrnrc=itrnrc+1
|
|
t1=un/t
|
|
dltp=t1/log(p/pm)
|
|
t0p=half*t0*t1
|
|
dele=alp-t0**4*dlt
|
|
delep=-t0**4*dlt*(two*t1+dltp/dlt)
|
|
dt=-dele/delep*t1
|
|
t=t*(un+dt)
|
|
write(6,645) id,itrnrc,dt,t
|
|
if(abs(dt).lt.1.e-6.or.itrnrc.gt.20) go to 220
|
|
t0=sqrt(t*tm)
|
|
dlt=log(t/tm)/log(p/pm)
|
|
go to 210
|
|
220 continue
|
|
alp=flrd(id)/t0**4/dlt
|
|
end if
|
|
c
|
|
230 continue
|
|
delta(id)=dlt
|
|
temp(id)=t
|
|
flr=alp*t0**4*dlt
|
|
if(dlt.ge.grdadb) then
|
|
flc=fc0*dele
|
|
write(6,646) id,itrnrc,t,dlt,grdadb,flrd(id)/flxtt(id),
|
|
* flr/flxtt(id),flc/flxtt(id),flxcn0/flxtt(id),
|
|
* (flr+flxcn0)/flxtt(id)
|
|
else
|
|
itrnrc=0
|
|
write(6,646) id,itrnrc,t,dlt,grdadb,flrd(id)/flxtt(id),
|
|
* flr/flxtt(id)
|
|
end if
|
|
646 format(2i4,f8.1,1p7e12.4)
|
|
end do
|
|
end if
|
|
c
|
|
if(ioptab.ge.-1.and.ifryb.gt.0) then
|
|
do id=1,nd
|
|
t=temp(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)
|
|
END DO
|
|
END IF
|
|
c
|
|
call tdpini
|
|
call conout(1,ipconf)
|
|
end if
|
|
c
|
|
return
|
|
end
|