131 lines
4.0 KiB
Fortran
131 lines
4.0 KiB
Fortran
SUBROUTINE CONOUT(IMOD,IPRIN)
|
|
C =============================
|
|
C
|
|
C Diagnostic outprint of temperature gradients, convective flux,
|
|
C and their derivatives
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
COMMON/CUBCON/A,B,DEL,GRDADB,DELMDE,RHO,FLXTOT,GRAVD
|
|
C
|
|
IF(IPRIN.GT.0) WRITE(6,600)
|
|
ANEREL=ELEC(1)/(DENS(1)/WMM(1)+ELEC(1))
|
|
ICBEG=0
|
|
FLXTO0=SIG4P*TEFF**4
|
|
DO ID=1,ND
|
|
T=TEMP(ID)
|
|
PTOT=PTOTAL(ID)
|
|
PG=PGS(ID)
|
|
PRAD=PTOT-PG-HALF*DENS(ID)*VTURB(ID)**2
|
|
if(prad.lt.0.) prad=0.
|
|
IF(IMOD.EQ.2) THEN
|
|
if(ioptab.ge.0) then
|
|
CALL OPACF0(ID,NFREQ)
|
|
CALL MEANOP(T,ABSO,SCAT,OPROS,OPPLA)
|
|
ABROSD(ID)=OPROS/DENS(ID)
|
|
else
|
|
call meanopt(t,id,dens(id),opros,oppla)
|
|
if(hmix0.gt.0.) abrosd(id)=opros
|
|
end if
|
|
END IF
|
|
FLXTOT=flxto0
|
|
if(idisk.eq.1) then
|
|
flxtot=flxto0*(un-thetav(id))
|
|
gravd=zd(id)*qgrav
|
|
prad=pradt(id)
|
|
end if
|
|
IF(ID.EQ.1) THEN
|
|
TAU=DM(ID)*ABROSD(ID)
|
|
FLXCR=0.
|
|
GRDADB=0.
|
|
DELTA(ID)=0.
|
|
FLXC(ID)=0.
|
|
FLXCNV=0.
|
|
ELSE
|
|
TM=TEMP(ID-1)
|
|
TAU=TAUM+HALF*(DM(ID)-DM(ID-1))*(ABROSD(ID)+ABROSD(ID-1))
|
|
PTOTM=PTOTAL(ID-1)
|
|
PGM=PGS(ID-1)
|
|
PRADM=PTOTM-PGM-HALF*DENS(ID-1)*VTURB(ID-1)**2
|
|
if(idisk.eq.1) pradm=pradt(id-1)
|
|
if(pradm.lt.0.) pradm=0.
|
|
if(ilgder.eq.0) then
|
|
T0=HALF*(T+TM)
|
|
PT0=HALF*(PTOT+PTOTM)
|
|
PG0=HALF*(PG+PGM)
|
|
PR0=HALF*(PRAD+PRADM)
|
|
AB0=HALF*(ABROSD(ID)+ABROSD(ID-1))
|
|
DLT=(T-TM)/(PTOT-PTOTM)*PT0/T0
|
|
else
|
|
T0=SQRT(T*TM)
|
|
PT0=SQRT(PTOT*PTOTM)
|
|
PG0=SQRT(PG*PGM)
|
|
PR0=SQRT(PRAD*PRADM)
|
|
AB0=SQRT(ABROSD(ID)*ABROSD(ID-1))
|
|
DLT=LOG(T/TM)/LOG(PTOT/PTOTM)
|
|
end if
|
|
DELTA(ID)=DLT
|
|
flxcnv=0.
|
|
if(idisk.ne.1.or.id.lt.nd)
|
|
* CALL CONVEC(ID,T0,PT0,PG0,PR0,AB0,DLT,FLXCNV,VCON)
|
|
if(hmix0.gt.0.) FLXC(ID)=FLXCNV
|
|
IF(ICBEG.EQ.0.AND.FLXC(ID).GT.0..AND.FLXC(ID-1).EQ.0..
|
|
* AND.ID.GT.25) ICBEG=ID
|
|
if(icbeg.gt.0.and.flxc(id).gt.0.) icend=id
|
|
END IF
|
|
PRADR=PRAD/PTOT
|
|
conrel=0.
|
|
radrel=1.
|
|
if(flxtot.gt.0.) then
|
|
conrel=FLXCNV/FLXTOT
|
|
radrel=flrd(id)/flxtot
|
|
end if
|
|
IF(IPRIN.GT.0) WRITE(6,601) ID,TAU,T,DELTA(ID),
|
|
* GRDADB,conrel,radrel,conrel+radrel
|
|
|
|
TAUM=TAU
|
|
END DO
|
|
if(iprin.gt.0) write(6,603) icbeg,icend
|
|
603 format(/' convective zone between depths (inclusive) ',2i4/)
|
|
C
|
|
c if ICONV=3, then:
|
|
c in the convective zone the radiative (+ convective) equilibrium
|
|
c is taken obligatorily in the differential form, i.e.
|
|
c NDRE is modified to have the value just below the beginning of
|
|
c convection zone, and
|
|
c rediff(id) has to be set to unity, and reint(id) to 0 for id => ndre
|
|
c
|
|
IF(ICBEG.GT.3.AND.ICONV.EQ.3) THEN
|
|
NDRE=ICBEG-1
|
|
DO ID=1,ND
|
|
IF(ID.GE.NDRE) THEN
|
|
REINT(ID)=0.
|
|
REDIF(ID)=1.
|
|
else
|
|
reint(id)=1.
|
|
redif(id)=0.
|
|
END IF
|
|
END DO
|
|
WRITE(6,602) NDRE
|
|
END IF
|
|
c
|
|
IF(ICBEG.GT.3.AND.ICONV.EQ.2) THEN
|
|
NDRE=ICBEG-1
|
|
DO ID=1,ND
|
|
IF(ID.GE.NDRE) THEN
|
|
REDIF(ID)=1.
|
|
END IF
|
|
END DO
|
|
WRITE(6,602) NDRE
|
|
END IF
|
|
c
|
|
600 FORMAT(//' ID',3X,'TAUR',7X,'TEMP',6X,
|
|
* 'DELTA',2X,'DELTA(AD)',2X,'CON/TOT RAD/TOT (C+R)/TOT'//)
|
|
601 FORMAT(1H ,I4,1PD9.2,0PF9.1,1P5D10.2)
|
|
602 FORMAT(//' NDRE IS RESET IN CONOUT DUE TO THE EXISTENCE OF'
|
|
* ,' CONVECTIVE ZONE'/' NDRE= ',I3/)
|
|
RETURN
|
|
END
|