      subroutine write4
*
*     print the structure of envelopes and convective transfer properties
*     both scaled down,  also interpolate between shells to find exact
*     physical properties when eta=0,20, at the bottom of the convection
*     zone and when the mass of the envelope is=1e-6* stellar mass
*
      implicit double precision(a-h,o-z)

      double precision m
      logical writeit

      common/shells/r(800),m(800),p(800),t(800),el1(800),ol1(800),
     1 rl1(800),atg(800),rtg(800),phs(800)
      common/tenv/rhok(3,60,35),eta(3,60,35),pp(3,60,35),datg(3,60,35),
     1 od(3,18,35),tk(3,60),uu(3,60,35),xtt(3,60,35),xrt(3,60,35)
      common/je/je
      common/comp/amhyhe,amheca,x(4)
      common/xcomp/xcomp(400,3)
      common/slr/sl(11)
      common/c/az(3),z(3),taux,ts,ps,sm,rs,err,te,gs,rho,ol1x,rl1x,el1x
     1 ,atgx,rtgx,tl,pl,xmass,acu,am,bk,an0,pi
      common/d/ ml(3,60),nmax,kind,lmax,nsuff,nlim,nlim1,j,iter,lmam, 
     1 nlum,ko,nhp
      common/dirac/ psi(141),f12(141),f32(141)
      common/grad/ drdtp(400),cp(400),ttg(400),drdtpx,cpx,pecx,ttgx,
     1 deptx,dpdr,dtdr,dmdr,dr,xr,xt,cv,drdptx
      common/surf/u(2,3),v(2,3),w(2,3),rm,bm,is,ks,ls,ms,kstart,npass 
      common/terp/stpms,rat1
      common/wprop/ bn2(400),ac2(400)
      common/writeit/writeit
      common/oldval/bold,teold 
      common/mixl/aml,mls
      dimension s(400),d(400),nbn2(400),amfrac(400)
*
* Format statements
*
12    format(/13h log(m/msun)=,f6.3,2x,13h log(r/rsun)=,f6.3,2x,11h log(
     1grav)=,f6.3,2x,11h log(teff)=,f6.3,2x,13h log(l/lsun)=,f6.3,2x,12h
     2 log(l*/m*)=,f6.3//                          12h iterations=,i2,
     359h  nc=-1 if convv + nondeg, 0 if convv + deg, +1 if nonconvv/)
33    format(' too few shells')
37    format(4x,6hradius,5x,11hlog(1-r/rs),2x,11hlog(1-m/ms),2x,4hlogp,
     1 4x,4hlogt,4x,6hlogrho,4x,7hopacity,5x,3heta,5x,3httg,8x,3hatg,8x,
     2 3hrtg,6x,2hnc,3x,4hne/z)
900   format(/,4x,6hradius,5x,11hlog(1-r/rs),2x,11hlog(1-m/ms),10x,3hbn2
     1,11x,3hac2,8x,2hje,4x,4hnbn2)
38    format(2x,1pe11.5,3x,0pf8.4,5x,f8.4,1x,2(1x,f7.3),2x,f7.3,3x,
     1 1pe9.3,2x,0pf6.2,3(2x,e9.3),2x,i2,3x,f5.3) 
902   format(2x,1pe11.5,3x,0pf8.4,5x,f8.4,2(2x,1pe16.9),0p,3f12.9,i5)
101   format(/1x,8henvelope,2x,9hstructure,16x,6hscale=,i2,5x,
     1 ' total shells: ',i5)
103   format('log(1-r/rs)=',f8.4,2x,'log(1-m/ms)=',f8.4,2x,'logp=',
     1 f7.3,2x,'logt=',f6.3,2x,'logrho=',f6.3,2x,'logkappa=',f6.3,2x,
     2 'eta=',f6.2)
104   format(//57h physical conditions at the bottom of the convection z
     1one/)
106   format(///////31h convective transfer properties,16x,6hscale=i2)
107   format(/4x,6hz=rs-r,5x,11hlog(1-r/rs),2x,11hlog(1-m/ms),3x,
     1    7hlocal g,6x,        6hpeclet,6x,3httg,6x,5hsuper,5x,10hmix le
     2ngth,2x,5hfc/ft,3x,9hconv velo,3x,5hv/vth)
108   format(2x,1pe11.5,3x,0pf8.4,5x,f8.4,4x,1pe9.3,4x,e9.3,3x,0pf5.3,
     1 3x,1pe9.3,3x,e9.3,3x,0pf5.3,3x,1pe9.3,3x,0pf5.3)
109   format(//' physical conditions when m(env)/m(star)=',1pd14.6/)
*
*     envelope structure
*
      rads=r(1)
      nmod=1
      alc4pig = -6.07666 - 2. * dlog10(gs)
      amass = 10.**sm * 1.989e33
      do 4 n=1,nlim 
987     format('m(n): ',f20.17)
        if( 1.-r(n)/rads .ge. 1.d-18 ) then
           d(n) = dlog10( 1.-r(n)/rads )
        else
           d(n) = -10000.
        endif
        amfrac(n)=10.**m(n)/sm
        ame1=1.d0-amfrac(n)
        if (ame1 .ge. 1.d-18) then
           s(n)=dlog10(ame1)
        else
            s(n) = -10000.
        endif
4     continue

      xs=dlog10(sl(j))
      xr=dlog10(rs) 
      xg=dlog10(gs) 
      xt=dlog10(te) 
      xm=dlog10(sm) 
      ratio=xs-xm

      ion=-1
      do 11 n=1,nlim,nmod
         if ( atg(n) .gt. rtg(n) ) then 
            nc=1
         else
            if ( el1(n) .lt. 0. ) then
               nc=-1
            else
               nc=0 
            endif
         endif
         opy=10.**ol1(n)
*
*     compute degree of ionization xion 
*
         if ( ion .gt. 0. ) then
            xion=1.0
         else
            zzz = xcomp(n,1)*z(1) + xcomp(n,2)*z(2) + xcomp(n,3)*z(3) 
            xion=10.**(1.5*t(n)-rl1(n))*am/bk/zzz/1.329153
     1             *f12(el1(n))
            if ( xion .gt. 1.0) then
               xion=1.0
               ion=1
            endif
         endif
11    continue
*
* end of loop
*
*     physical properties at places of interest
*
      it=1
      stop=dlog10(1.-10.**stpms)
      do 90 n=1,nlim
         l=n
         if ( s(n) .gt. stop) goto 91
90    continue
91    slope=(s(l-1)-stop)/(s(l-1)-s(l)) 
*
* returns to 22 from below
*
22    dfin=d(l-1)+(d(l)-d(l-1))*slope
      sfin=s(l-1)+(s(l)-s(l-1))*slope
      pfin=p(l-1)+(p(l)-p(l-1))*slope
      tfin=t(l-1)+(t(l)-t(l-1))*slope
      rfin=rl1(l-1)+(rl1(l)-rl1(l-1))*slope
      ofin=ol1(l-1)+(ol1(l)-ol1(l-1))*slope
      efin=el1(l-1)+(el1(l)-el1(l-1))*slope
      if ( it .le. 1 ) then
         stp = 10.**stop
         slum=sl(1) 
         u(1,is)=pfin
         u(2,is)=tfin
         v(1,is)=dlog10(slum) 
         v(2,is)=dlog10(1-10.**dfin)+w(2,is)
         it=it+1
*
*  skip over the next section with the following statement. 
*  the next section computes the convective properties.
*  it's caused some problems in the hot models
*
         if ( nlim1 .le. 0 ) go to 30
         do 25 n=kind,nlim
            l=n
            if ( atg(n) .gt. rtg(n) ) goto 26
25       continue
26       slope=(atg(l-1)-rtg(l-1))/(atg(l-1)-rtg(l-1)-atg(l)+rtg(l))
         go to 22
      endif
*
* bottom of loop
*
* convection zone
*
      nmod1=((nlim1-kind)/100)+1
*
* compute geometrical depth zc, local gravity glocl, superadiabaticity
* super, ratio of convective to total flux fcf, convective velocity
* conv, ratio of convective to thermal velocities vcvt and peclet
* number
*
      do 52 n=kind,nlim1,nmod1
         zc=rads-r(n)
         glocl=27398.*10.**m(n)*(6.96e+10/r(n))**2
         super=ttg(n)-atg(n)
         fcf=1.-ttg(n)/rtg(n) 
         if ( fcf .le. 0.) then
            fcf=0.
            conv=0. 
         else
            conv=27398.*(te/5800.)**4*10.**(m(n)-t(n)-rl1(n))
     1           *3.90e+33*phs(n)
     1           *drdtp(n)*fcf/(16.*pi*cp(n)*r(n)**2)
            conv=conv*aml
*
* choose mixing-length version here
* (mls) =0 Bohm-Vitense or (mls=1) Bohm & Cassinelli
* default is Bohm-Vitense version
*
            if(mls.eq.1)then
*
* ml2 or ml3 version
*
              conv=2.*dsqrt(2.d0)*conv
            endif
            if ( conv .gt. 0.0 ) conv = conv**0.333333
         endif
         azave = xcomp(n,1)*az(1)+xcomp(n,2)*az(2)+xcomp(n,3)*az(3) 
         vcvt=conv/(3.*bk*an0/azave*10.**t(n))**0.5
*
* also need to adjust the peclet #. default is Bohm-Vitense.
*
         peclet=cp(n)*phs(n)*conv/4.535688e-04*10.**(2.*rl1(n)+ol1(n)-
     1         3.*t(n))
         peclet=peclet*aml*aml
*
* Bohm & Cassinelli version
*
         if(mls.eq.1)then
           peclet=peclet/dsqrt(2.d0)
         endif
52    continue

30    continue

      do 901 n=1,nlim,nmod
         nbn2(n)=bn2(n)/dabs(bn2(n))
         if (bn2(n) .ge. 0.) then
            bn2(n)=dlog10(bn2(n))
         else
            bn2(n) = -7.
         endif
         if (ac2(n) .ge. 0.) then
            ac2(n)=dlog10(ac2(n))
         else
            ac2(n) = -7.
         endif
901   continue
1901  format(5(2x,f15.10))

      return
      end 

************************************************************************

