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