**==PRINTB.FOR
      SUBROUTINE PRINTB(DTY,PER,NWRT1,NWRT2,NWRT3,NWRT4,NMONT,NMOD,IEND)
      IMPLICIT REAL*8(A-H, L, M, O-Z)
      SAVE
      INTEGER MAXMSH
      PARAMETER (MAXMSH = 2000)
      COMMON H(60,MAXMSH),DH(60,MAXMSH),EPS,DEL,DH0,NMESH,JIN,IW(200)
      COMMON /NUCMAT/ HNUC(100,MAXMSH), DHNUC(100,MAXMSH)
      COMMON /TRANS / HT(26,MAXMSH,2)
      COMMON /SODDS / ALPHA, RML, CMG, CSI, CFE, CT(10), AGE, DT, M1,
     :  EC, BM, ANG, CM, MTA, MTB, TM(2), T0, M0, TC(2), OS, AC, RCD,
     :  RMG, RHL, XF, DR, AK1 ,RMT, AK2, ITH, IX, IY, IZ, IB, ISX(45),
     :  TRB
* extra common for mesh-spacing
      COMMON /PMESH / PMH(2), PME(2), IAGB
      COMMON /INF   / Q(60)
      COMMON /DINF  / QD(60)
      COMMON /OP    / ZS, LEDD, VM, GR, GRAD, ETH, RLF, EGR, R, QQ
      COMMON /ABUND / WA1(10), YA(10), WA(3), AVM
c GR is ov value, XX is Schwarzschild value, of gradr-grada
      COMMON /STAT2 / PL, RL, U, P, RHO, FK, T, SF, ST, ZT, GRADA, CP,
c     :                CH, S, PR, PG, PF, PT, EN, RR(20), EX, WX, EXX(3),
c     :                RW(16)
     :                CH, S, PR, PG, PF, PT, EN, RPP, R33, R34, RBE, RBP,
     :                RPC, RPN, RPO, R3A, RAC, RAN, RAO, RANE, RCC, RCO,
     :                ROO, RGNE, RGMG, RCCG, RPNG, EX, ENX, EXX(3), AAWT,
     :                WT, RW(14)
      COMMON /YUK1  / PX(34), WMH, WMHE, VMH, VME, VMC, VMG, BE, VLH,
     :                VLE, VLC, VLN, VLT, MCB(12),WWW(100)
      COMMON /CNSTS / CPI, PI4, CLN10, CA, CB, CL, CD, CG, CR, CC, CEVB,
     &                CEN, CPL, CMEVMU, CSECYR, LSUN, MSUN, RSUN, TSUNYR
      COMMON /ATDATA/ DH2(4), CHI(26,9), OMG(27), AM(10), BN(10), JZ(10)
      COMMON /NCDATA/ QPP, Q33, Q34, QBE, QBP, QPC, QPNA, QPO, Q3A, QAC,
     :                QAN, QAO, QANE, QCCA, QCO, QOO, QGNE, QGMG, QCCG,
     :                QPNG, CNUC(390)
      COMMON /WT    / PWT,SG,MK,MT2,LQP,LKP
*
* Extra COMMON for main-sequence evolution.
*
      COMMON /ZAMS  / TKH(2)
      COMMON /MESH  / TRC1,TRC2,DD,DT1,DT2,MWT,MWTS,IVMC,IVMS
      COMMON /MIX   / KBICZ,KICZ
      COMMON /DHBLOC/ IDREDGE
      COMMON /WIND  / WIND,SWIND,PERIOD
      COMMON /EVMODE/ IMODE
      COMMON /INERTI/ VI(2)
      COMMON /MONTAG/ WCV, WL
C Extra common for extra timestep control stuff - P for previous, C for current
      COMMON /DTCONT/ VLHP(2), VLEP(2), VLCP(2), RLFP(2), TOTMP(2), VLHC(2),
     :     VLEC(2), VLCC(2), RLFC(2), TOTMC(2)
      COMMON /WINDS / WINDML(2), FAKEWIND(2), BE2
      COMMON /TIDES / MENV(2), RENV(2)
      PS(VX) = 0.5D0*(VX+DABS(VX))
      DIMENSION XA(10), TCB(12), RCB(12)
      DIMENSION SX(45,MAXMSH+1)
      DIMENSION XINIT(44),YIELD(44)
      CHARACTER*5 AX(34), BX(16), CX(16), DX(17)
      INTEGER KCB(12), KMX(3), KEX(12), KENV(2)
      REAL*8 EMAX(3), VMX(3), MEX(12), THB(12) ! EXX(3),
      INTEGER IPX(25)
      DATA IPX/9, 17, 2, 3, 4, 5, 6, 8, 10, 11, 12, 13, 14, 15, 16,
     &     18, 19, 20, 21, 28, 27, 7, 24, 25, 26/
C Data required to compute yields - Anders & Grevasse '89 abunds
C Used A&G for H,He,C,N,O,Ne rather than input for convenience
      DATA XINIT/4.801d-5,2.929d-5,9.353d-9,0d0,4.725d-9,3.650d-5,0d0,
     :     4.363d-6,3.887d-6,2.167d-5,4.051d-7,4.127d-6,1.302d-4,0d0,
     :     3.339d-5,5.148d-4,6.766d-5,7.760d-5,0d0,0d0,5.798d-5,6.530d-4,
     :     3.426d-5,2.352d-5,8.155d-5,3.958d-4,3.222d-6,1.866d-5,1.169d-3,
     :     2.855d-5,3.697d-6,0d0,0d0,3.358d-6,4.944d-5,0d0,1.958d-5,8.594d-7,
     :     7.057d-1,2.752d-1,3.032d-3,1.105d-3,9.592d-3,1.619d-3/
      DATA YIELD/0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,
     :     0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,
     :     0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0/
      DATA IDREDGE/0/
      DATA AX/' psi ','  P  ',' rho ','  T  ','kappa','grada',' grad',
     :'gr-ga','  m  ','  H1 ',' He4 ',' C12 ',' N14 ',' O16 ',' Ne20',
     :' He3 ','  r  ','  L  ',' Eth ',' Enuc',' Eneu','  dm ',' k**2',
     :'n/n+1','lr/lp','lm/lp',' U   ',' S   ','L/Edd',' mu  ','  mu ',
     :' LK ','Ebind','beta '/
      DATA BX/'  g  ','  n  ','  D  ',' He3 ',' Li7 ',' Be7 ',' B11 ',
     :     ' C13 ',' C14 ',' N15 ',' O17 ',' O18 ',' F19 ',' Ne21',
     :     ' Ne22','  m  '/
      DATA CX/' Na22',' Na23',' Mg24',' Mg25',' Mg26','Al26M',
     :     'Al26G',' Al27',' Si28',' Si29',' Si30',' P31 ',' S32 ',
     :     ' S33 ',' S34 ','  m  '/
      DATA DX/' Fe56',' Fe57',' Fe58',' Fe59',' Fe60',' Co59',
     :     ' Ni58',' Ni59',' Ni60',' Ni61', '  p  ',' He4 ',' C12 ',
     :     ' N14 ',' O16 ',' Ne20','  m  '/
99001 FORMAT (/, '  K', 17(4X, A5, 2X),/)
99002 FORMAT (I4, 1P, 15(1X,E10.3))
99003 FORMAT (/, 1X, I7, '  dt/age/MH/MHe  tn/tKH/Mb  P/rlf/dM LH/LHe/',
     :'LC Lth/Lnu/m  H1      He4     C12     N14     O16     Ne20    ',
     :'He3       psi    dens    temp', 2(/, F8.4, 1P, E16.9, 4E10.3,
     :0P, 7F8.5, F8.3, 2F8.4, 2X, A4), /, 3F8.4, F10.4, 1P, 2E10.3, 0P,
     :F10.4, 7F8.5, F8.3, 2F8.4, 2X, A4, /, F8.4, 12F8.3, 5F8.4, I6,
     :/,2I4,12I8,/,15F8.4,1P,3E10.3,0P,/,15I8)
99004 FORMAT (I13, 1P, 9E13.5, 3(/, 10E13.5))
99005 FORMAT (/, '  K', 15(6X, A5, 3X),/)
99006 FORMAT (I4, 1P, 11E14.7)
C Store previous values of extra timestep control variables
      VMHP = VMH
      DO I=1,2
         VLHP(I) = VLHC(I)
         VLEP(I) = VLEC(I)
         VLCP(I) = VLCC(I)
         RLFP(I) = RLFC(I)
         TOTMP(I) = TOTMC(I)
      END DO
            DO ISTAR = 1,IMODE
               JC = 1
               JTHB = 1
               JE = 1
               KTM = 0
               KCE = 1
               KMH = 0
               KME = 0
               TMAX = 0.0D0
               VI(ISTAR) = 0d0
               DO K=1,NMESH - 1
                  VI(ISTAR) = VI(ISTAR) + 2.0/3.0*(DEXP(H(4+15*(ISTAR-1),K)+DH(4+15*(ISTAR-1),K))
     : - DEXP(H(4+15*(ISTAR-1),K+1)+DH(4+15*(ISTAR-1),K)))*DEXP(H(7+15*(ISTAR-1),K)
     :                 +DH(7+15*(ISTAR-1),K))**2.0
               END DO
               VI(ISTAR) = VI(ISTAR) + 2.0/5.0*DEXP(H(4+15*(ISTAR-1),NMESH)+DH(4+15*(ISTAR-1),NMESH))*
     :              DEXP(H(7+15*(ISTAR-1),NMESH)+DH(7+15*(ISTAR-1),NMESH))**2.0
               DO J = 1, 3
                  KMX(J) = 1
                  VMX(J) = 0.0D0
C     EXX(J) = 0.0D0
                  EMAX(J) = 0.0D0
               END DO
               EXLIM = 10.0*H(8,1)/EXP(H(4,1))
               DO J = 1, 34
                  PX(J) = 0.0D0
                  SX(J,1) = 1.0D-8
                  IF (J .LE. 12) KCB(J) = 0
                  IF (J .LE. 12) MCB(J) = 0d0
                  IF (J .LE. 12) KEX(J) = 0
                  IF (J .LE. 12) MEX(J) = 0.0D0
                  IF (J .LE. 12) THB(J) = 0.0D0
               END DO
               WMH = 0d0
               WMHE = 0d0
               VMH = 0d0
               VME = 0d0
               VMC = 0d0
               VMG = 0d0
               BE = 0d0
               VLH = 0d0
               VLE = 0d0
               VLC = 0d0
               VLN = 0d0
               VLT = 0d0
               NPRINT = 0
C     Decide whether to print the interior or not
               IF ( MOD(NMOD,NWRT1).EQ.0 ) NPRINT = 1               
               IF ( NPRINT.EQ.1 ) WRITE(32+20*(ISTAR-1), 99001) (AX(ISX(J)), J=1,15)
               IF (NPRINT.EQ.1.AND.IW(102).NE.0) WRITE (36+20*(ISTAR-1),*) NMOD
               IF (NPRINT.EQ.1.AND.IW(102).NE.0) WRITE (37+20*(ISTAR-1),*) NMOD
               IF (NPRINT.EQ.1.AND.IW(102).NE.0) WRITE (38+20*(ISTAR-1),*) NMOD
      IF ( NPRINT.EQ.1.AND.IW(102).NE.0 ) WRITE(36+20*(ISTAR-1), 99001) (BX(J), J=1,16)
      IF ( NPRINT.EQ.1.AND.IW(102).NE.0 ) WRITE(37+20*(ISTAR-1), 99001) (CX(J), J=1,16)
      IF ( NPRINT.EQ.1.AND.IW(102).NE.0 ) WRITE(38+20*(ISTAR-1), 99001) (DX(J), J=1,17)
         DO KK = 1, NMESH
            K = NMESH + 1 - KK
            DO J = 1, 20
               IF ( IEND.LT.0 ) QD(J) = 0.0D0
               IF ( IEND.LT.0 ) Q(J) = H(J+15*(ISTAR-1), K)
               IF ( IEND.GE.0 ) QD(J) = DH(J+15*(ISTAR-1), K)
               IF ( IEND.GE.0 ) Q(J) = H(J+15*(ISTAR-1), K) + QD(J)
            END DO
            IF (ISTAR.EQ.2) Q(28) = H(13, K)
            CALL FUNCS1 (-1, K, 1, NMESH, ISTAR)
C call mass-loss to get surface conditions
            IF (ISTAR.EQ.1.AND.IMODE.EQ.2.AND.K.EQ.1) THEN
               DO J=1,30
                  IF ( IEND.LT.0 ) QD(J) = 0.0D0
                  IF ( IEND.LT.0 ) Q(J) = H(J, K)
                  IF ( IEND.GE.0 ) QD(J) = DH(J, K)
                  IF ( IEND.GE.0 ) Q(J) = H(J, K) + QD(J)
               END DO
               CALL MASSLOSS(-1,A1,A2,A3,A4,A5)
            END IF
*
* Correct abundances for actual atomic masses
*
            DO I = 1, 10
               XA(I) = AM(I)*YA(I)/AVM
            END DO
            XH = XA(1)
            XHE = XA(2)
            XC = XA(3)
            XN = XA(4)
            XO = XA(5)
            XNE = XA(6)
            XMG = XA(7)
            XSI = XA(8)
            XFE = XA(9)
            XHE3 = XA(10)
C Evaluate the functions to be printed
            WF=DSQRT(1.0D0+DEXP(Q(1)))
            PX(1) = 2.0D0*(WF-DLOG(WF+1.0D0))+Q(1)
            PX(2) = P
            PX(3) = RHO
            PX(4) = T
            TMAX = MAX(TMAX, SX(4,KK))
            IF ( PX(4).GE.TMAX ) KTM = KK + 1
            PX(5) = FK
            PX(6) = GRADA
            PX(7) = GRAD
            PX(8) = DMIN1(9.99D0, GR)
            PX(9) = VM/MSUN
C locate convective/radiative boundaries (GR=0)
            IF (KK.LT.2.OR.JC.GT.12.OR.SX(8,KK)*PX(8).GT.0.0D0
     &           .OR. AAWT.LT.1.0D-1) GO TO 1
C            IF (K.GT.KBICZ+1) GOTO 1
               MCB(JC) = (SX(9,KK)*PX(8)-PX(9)*SX(8,KK))/(PX(8)-SX(8,KK))
               KCB(JC) = K
               TCB(JC) = (SX(4,KK)*PX(8)-PX(4)*SX(8,KK))/(PX(8)-SX(8,KK)) !T
               RCB(JC) = (SX(3,KK)*PX(8)-PX(3)*SX(8,KK))/(PX(8)-SX(8,KK)) !RHO
C     Find T at base of conv. envelope
            IF (K.GT.50.AND.K.LT.300) THEN
               TBCE = (SX(4,KK)*PX(8)-PX(4)*SX(8,KK))/(PX(8)-SX(8,KK))
            END IF
C Find T and rho of convective pulse extrema
C         IF (K.LT.750) THEN
C            TCB(JC) = T
C            RCB(JC) = RHO
C         END IF
            IF ( PX(4).GT.1.0D5 ) KCE = KK
            JC = JC + 1
C locate convective/semiconvective boundaries (loosely, GR=0.01)
 1          IF ((PX(8)-DR)*(SX(8,KK)-DR).GT.0.0D0.OR.KK.LT.2.OR.JC.GT.12
     &           .OR. AAWT.LT.1.0D-1) GO TO 2
C            IF (K.GT.KBICZ+1) GOTO 2
               MCB(JC)=-(SX(9,KK)*(PX(8)-DR)-PX(9)*(SX(8,KK)-DR))
     :              /(PX(8)-SX(8,KK))
               KCB(JC) = -K
               JC=JC + 1
C locate burning shell boundaries
 2          IF (XH.GT.XF.AND.SX(10,KK).LT.XF) THEN
               VMH = (PX(9)*(XF-SX(10,KK))+SX(9,KK)*(XH-XF))/(XH-SX(10,KK))
               KMH = K
               PMH(ISTAR) = (PX(2)*(XF-SX(10,KK))+SX(2,KK)*(XH-XF))/(XH-SX(10,KK))
            END IF
            IF (XHE.GT.XF.AND.SX(11,KK).LT.XF.AND.K.GT.NMESH/2) THEN
               VME = (PX(9)*(XF-SX(11,KK))+SX(9,KK)*(XHE-XF))/
     &              (XHE-SX(11,KK))
               KME = K
               PME(ISTAR) = (PX(2)*(XF-SX(11,KK))+SX(2,KK)*(XHE-XF))/
     &              (XHE-SX(11,KK))
            END IF
            IF (XC.GT.XF.AND.SX(12,KK).LT.XF) VMC = (PX(9)*
     :           (XF-SX(12,KK))+SX(9,KK)*(XC-XF))/(XC-SX(12,KK))
C N abundance in H-burning ash
            IF (XH.GT.1d-10.AND.SX(10,KK).LT.1d-10) XASH = XN
C     Find T_hmid (XH = 0.345)
            IF (XH.GT.0.345.AND.SX(10,KK).LT.0.345) THEN
               THMID = (PX(4)*(0.345-SX(10,KK))+SX(4,KK)*(XH-0.345))/
     :              (XH-SX(10,KK))
            END IF
C Find T_hbase (XH = 0.05)
            IF (XH.GT.0.05.AND.SX(10,KK).LT.0.05) THEN
               THBASE = (PX(4)*(0.05-SX(10,KK))+SX(4,KK)*(XH-0.05))/
     :              (XH-SX(10,KK))
            END IF
C Find T_intmid (Xhe = 0.45)
            IF (XHE.GT.0.45.AND.SX(11,KK).LT.0.45) THEN
               TINTMID = (PX(4)*(0.45-SX(11,KK))+SX(4,KK)*(XHE-0.45))/
     :              (XHE-SX(11,KK))
            END IF
C Find T_intbase (Xhe = 0.05)
            IF (XHE.GT.0.05.AND.SX(11,KK).LT.0.05) THEN
               TINTBASE = (PX(4)*(0.05-SX(11,KK))+SX(4,KK)*(XHE-0.05))/
     :              (XHE-SX(11,KK))
            END IF
C find points of max energy generation
            EXX(1) = CMEVMU*((QPP + 0.5*Q33)*RPP + QPC*RPC + QPO*RPO
     :           + QPNA*RPN + QPNG*RPNG)
            EXX(2) = CMEVMU*(Q3A*R3A + QAC*RAC + QAN*RAN + QAO*RAO
     :           + QANE*RANE)
            EXX(3) = CMEVMU*(QCCA*RCC + QCCG*RCCG + QCO*RCO +QOO*ROO
     :           + QGNE*RGNE + QGMG*RGMG)
            DO J = 1, 3
               IF (EXX(J).GE.EMAX(J)) THEN
                  IF (EXX(J).GE.EXLIM) KMX(J) = K
                  IF (EXX(J).GE.EXLIM) VMX(J) = PX(9)
                  EMAX(J) = EXX(J)
               END IF
            END DO
            PX(10) = XH
            PX(11) = XHE
            PX(12) = XC
            PX(13) = XN
            PX(14) = XO
            PX(15) = XNE
            PX(16) = XHE3 !XMG
            PX(17) = R/RSUN
            PX(18) = Q(8)/LSUN
            PX(19) = ETH
            PX(20) = EX
C locate boundaries of burning regions (EX > EXLIM)
            IF((PX(20)-EXLIM)*(SX(20,KK)-EXLIM).LE.0.0D0.AND.JE.LT.12.AND.
     &           KK.GE.2) THEN
               MEX(JE)= (SX(9,KK)*(PX(20)-EXLIM)-PX(9)*(SX(20,KK)-EXLIM))
     :              /(PX(20)-SX(20,KK))
               KEX(JE) = K
               JE = JE  + 1
            END IF
            PX(21) = -EN
            PX(22) = PX(9)-SX(9,KK)
            WF = SX(17,KK)/PX(17)
            PX(23) = (SX(9,KK)*SX(23,KK)*WF**2+PX(22)*0.4D0*(1.0D0-WF**5)
     :           /(1.0D0-WF**3))/PX(9)
            WF = 1.0D0/DLOG10(PX(2)/SX(2,KK))
            IF ( KK.EQ.1 ) WF = 0.0D0
            PX(24) = WF*DLOG10(PX(3)/SX(3,KK))
            PX(25) = -WF*DLOG10(PX(17)/SX(17,KK))
            PX(26) = -WF*DLOG10(PX(9)/SX(9,KK))
            PX(27) = U
            PX(28) = S
            PX(29) = Q(8)/LEDD
            TADJ = FK*CP*T*(1.0D11*PX(22)/(4*PI4*R*R))**2/(3*CL*PR)/CSECYR
            PX(30) = LKP !TADJ
            PX(31) = CR*RHO*T/PG
C         PX(32) = HT(2,K)*HT(4,K)*1.0D66
            PX(32) = MT2*PS(PX(30) - SX(30,KK))! Thermohaline mixing coeff ! LKP
C RJS 4/4/07 Locate TH boundaries
            IF (KK.GT.2.AND.JTHB.LT.12.AND.((PX(32).LE.-1d-20.AND
     :   .SX(32,KK).GT.-1d-20).OR.(PX(32).GT.-1d-20.AND.SX(32,KK).LE.-1d-20))) THEN
               THB(JTHB) = (SX(9,KK)*PX(32)-PX(9)*SX(32,KK))/(PX(32)-SX(32,KK))
               JTHB = JTHB + 1
            END IF
C Print the interior details on first `page', if required
            IF ( NPRINT.EQ.1.AND.MOD(K-1,NWRT2).EQ.0 ) WRITE(32+20*(ISTAR-1),99002) K,
     :           (PX(ISX(J)), J=1,15)
99007       FORMAT (I4, 1P, 17(1X,E10.3))
            IF ( NPRINT.EQ.1.AND.MOD(K,10*NWRT2).EQ.0 ) WRITE(32+20*(ISTAR-1),99002)
            IF ( NPRINT.EQ.1.AND.MOD(K-1,NWRT2).EQ.0.AND.IW(102).NE.0) 
     :           WRITE(36+20*(ISTAR-1),99007) K,
     :           HNUC(1+50*(ISTAR-1),K), HNUC(2+50*(ISTAR-1),K), HNUC(3+50*(ISTAR-1),K),
     :           HNUC(4+50*(ISTAR-1),K), HNUC(5+50*(ISTAR-1),K), HNUC(6+50*(ISTAR-1),K),
     :           HNUC(7+50*(ISTAR-1),K), HNUC(8+50*(ISTAR-1),K), HNUC(9+50*(ISTAR-1),K),
     :           HNUC(10+50*(ISTAR-1),K),HNUC(11+50*(ISTAR-1),K),HNUC(12+50*(ISTAR-1),K),
     :           HNUC(13+50*(ISTAR-1),K),HNUC(14+50*(ISTAR-1),K),HNUC(15+50*(ISTAR-1),K),
     :           PX(9)
            IF ( NPRINT.EQ.1.AND.MOD(K,10*NWRT2).EQ.0.AND.IW(102).NE.0)
     :           WRITE(36+20*(ISTAR-1),99002)
            IF ( NPRINT.EQ.1.AND.MOD(K-1,NWRT2).EQ.0.AND.IW(102).NE.0) 
     :           WRITE(37+20*(ISTAR-1),99007) K,
     :           HNUC(16+50*(ISTAR-1),K),HNUC(17+50*(ISTAR-1),K),HNUC(18+50*(ISTAR-1),K),
     :           HNUC(19+50*(ISTAR-1),K),HNUC(20+50*(ISTAR-1),K),HNUC(21+50*(ISTAR-1),K),
     :           HNUC(22+50*(ISTAR-1),K),HNUC(23+50*(ISTAR-1),K),HNUC(24+50*(ISTAR-1),K),
     :           HNUC(25+50*(ISTAR-1),K),HNUC(26+50*(ISTAR-1),K),HNUC(27+50*(ISTAR-1),K),
     :           HNUC(28+50*(ISTAR-1),K),HNUC(29+50*(ISTAR-1),K),HNUC(30+50*(ISTAR-1),K),
     :           PX(9)
            IF ( NPRINT.EQ.1.AND.MOD(K,10*NWRT2).EQ.0.AND.IW(102).NE.0) 
     :           WRITE(37+20*(ISTAR-1),99002)
            IF ( NPRINT.EQ.1.AND.MOD(K-1,NWRT2).EQ.0.AND.IW(102).NE.0 )
     :           WRITE(38+20*(ISTAR-1),99007) K,
     :           HNUC(31+50*(ISTAR-1),K),HNUC(32+50*(ISTAR-1),K),HNUC(33+50*(ISTAR-1),K),
     :           HNUC(34+50*(ISTAR-1),K),HNUC(35+50*(ISTAR-1),K),HNUC(36+50*(ISTAR-1),K),
     :           HNUC(37+50*(ISTAR-1),K),HNUC(38+50*(ISTAR-1),K),HNUC(39+50*(ISTAR-1),K),
     :           HNUC(40+50*(ISTAR-1),K),HNUC(41+50*(ISTAR-1),K),HNUC(42+50*(ISTAR-1),K),
     :           HNUC(43+50*(ISTAR-1),K),HNUC(44+50*(ISTAR-1),K),HNUC(45+50*(ISTAR-1),K),
     :           HNUC(46+50*(ISTAR-1),K), PX(9)         
            IF ( NPRINT.EQ.1.AND.MOD(K,10*NWRT2).EQ.0.AND.IW(102).NE.0)
     :           WRITE(38+20*(ISTAR-1),99002)
            DO J=1,34
               SX(J,KK+1)=PX(J)
            END DO
C Some integrated quantities, to be printed in the short summary
            WMH = WMH + PX(22)*XH
            WMHE = WMHE + PX(22)*XHE
            IF (PX(9).GT.VMHP) THEN
               BE = BE - PX(22)*EGR
            ELSE
               BE = 0d0
            END IF
            BE2 = BE+0.01
         PX(33) = BE
         PX(34) = PR/(PG+PR)
            WF = PX(22)*MSUN/LSUN
            VLH = VLH + WF*EXX(1)
            VLE = VLE + WF*EXX(2)
            VLC = VLC + WF*EXX(3)
            VLN = VLN - WF*EN
            VLT = VLT + WF*ETH
            IF (NMONT.NE.0.AND.MOD(NMOD,NMONT).EQ.0) THEN
C Output for use in MONTAGE -- may need smart output control RJS 29/5/08
99008          FORMAT (1P,12(1X,E14.7))
99018          FORMAT (I4, 1X, I6, 1X, 1P, E20.13, 2(1X,E12.5))
               IF (K.EQ.NMESH) THEN
                  WRITE (42+20*(ISTAR-1),99018) NMESH, NMOD, AGE*CSECYR, VLH, VLE
               END IF
               WRITE (42+20*(ISTAR-1),99008) DLOG10(PX(4)), DLOG10(PX(3)), DLOG10(WL),
     :              DLOG10(DMAX1(WCV,1d-30)), DLOG10(PX(17)*RSUN*1d11), PX(9), 
     :              PX(10), PX(16), PX(11), PX(12), PX(13), PX(14) 
            END IF
         END DO
         IF (NWRT3.EQ.1) GO TO 3
C Write further `pages' for each detailed model, if required
         DO I=2, NWRT3
            IJ = 15*I-15
            IF ( NPRINT.EQ.1 ) WRITE(32+20*(ISTAR-1), 99001) (AX(ISX(J+IJ)), J=1,15)
            DO KK = 1, NMESH
               K = NMESH + 1 - KK
C Have to choose the SX's that are wanted
               IF ( NPRINT.EQ.1.AND.MOD(K-1,NWRT2).EQ.0 ) WRITE(32+20*(ISTAR-1),99002)
     :              K, (SX(ISX(J+IJ),KK+1), J=1,15)
               IF ( NPRINT.EQ.1.AND.MOD(K,10*NWRT2).EQ.0 ) WRITE(32+20*(ISTAR-1),99002)
            END DO
         END DO
 3       TM(ISTAR) = VM
C         PER = BM/3.55223D0*(ANG/(TM(1)*(BM-TM(1))))**3
         IF (ISTAR.EQ.1) THEN
            PER = BM/3.55223D0*(Q(13)/(TM(1)*(BM-TM(1))))**3
         ELSE
C            PER = BM/3.55223D0*(H(13,1)/(TM(1)*(BM-TM(1))))**3
            PER = (TM(1)+TM(2))/3.55223D0*(H(13,1)/(TM(1)*TM(2)))**3
         END IF
         DMT = PX(9)*DH(4+15*(ISTAR-1), 1)/DTY
         TKH(ISTAR) = DABS(BE)/(Q(8)*CSECYR)
         TN = 4.0D10*PX(9)/Q(8)
C convective mixing timescale needed in FUNCS1
         TC(ISTAR) = RCD/(CSECYR*TN*DSQRT(VM))
*
* Write to the numerical data storage unit for plotting purposes.
*
C Separation
         IF (IMODE.EQ.2) THEN
            SEP = (TM(1)+TM(2))*(H(13,1)/(TM(1)*TM(2)))**2.0
         ELSE
            SEP = BM*(H(13,1)/(TM(1)*(BM-TM(1))))**2.0
         END IF
C Total angular momentum
         HTOT = H(14,1) + H(29,1) + H(13,1)
C Total spin angular momentum
         HSPINTOT = H(14,1) + H(29,1)
C Orbital angular velocity
         OORB = H(13,1)*(TM(1)+TM(2))/(TM(1)*TM(2)*SEP**2.0)
         VIORB = (TM(1)*TM(2)/(TM(1)+TM(2)))*SEP**2.0
         OSPIN = Q(14)/VI(ISTAR)
         OCRIT = DSQRT(6.67d-11*PX(9)*2d30/(6.96d8*PX(17))**3.0) !/DSQRT(CG)
         SEP = SEP/RSUN 
         WRITE (33+20*(ISTAR-1),115) NMOD,AGE,LOG10(PX(17)),LOG10(PX(4)),
     &        LOG10(PX(18)),PX(9),VMH,VME,LOG10(MAX(VLH,1.01D-10)),
     &        LOG10(MAX(VLE,1.01D-10)), LOG10(MAX(VLC,1.01D-10)),
     &        MCB,VMX(1),VMX(2),LOG10(FK),DTY,(PX(JJ),JJ=10,14), PX(16),
     &        RLF,Q(14),PER,SEP,BM/MSUN,H(13,1),HSPINTOT,HTOT,
     &        OORB*DSQRT(CG), OSPIN*DSQRT(CG), VI(ISTAR),OCRIT, DMT,MEX,THB,MENV(ISTAR)/MSUN,
     &        RENV(ISTAR)/RSUN, DLOG10(SX(3,2)), DLOG10(SX(4,2))
C There are 99 things in the format statement and 74 have been used.
 115     FORMAT (I6,1P,E16.9,0P,24F10.5,1P,3E13.6,18(1X,E12.5),0P,52F9.5)
         CALL FLUSH(33+20*(ISTAR-1))
C Write output for nucleosynthesis stuff
         write (41+20*(ISTAR - 1),116) NMOD,AGE,VMH,VME,VMH - VME,PX(9),XASH,(TCB(I)/1d8,
     :        I=1,12),RCB,
     :        TBCE,THMID,THBASE,TINTMID,TINTBASE,PERIOD,DABS(DMT)
 116     format (I6,1P,E16.9,0P,17(1X,E13.6),1P,19(1X,E13.6),15I5)
         CALL FLUSH(41+20*(ISTAR - 1))
         IF (IW(102).NE.0) THEN
C Write surface/centre abundances, yields -- but only if using NS code!
            write (39+20*(ISTAR - 1),117) NMOD,AGE,(HNUC(I+50*(ISTAR-1),1),I=1,50),
     :           WINDML(ISTAR)/MSUN*CSECYR, FAKEWIND(ISTAR)/MSUN*CSECYR, DTY
 117        format (I6,1P,E16.9,54(1X,E13.6))
            CALL FLUSH(39+20*(ISTAR - 1))
            write (40+20*(ISTAR - 1),117) NMOD,AGE,(HNUC(I+50*(ISTAR-1),NMESH),I=1,50)
     :           , H(5+15*(ISTAR - 1),NMESH)
            call flush(40+20*(ISTAR - 1))
         END IF
*
* Write detailed model data for plotting purposes. -- Note new unit number needed! RJS
*
C      NWRT5 = 21
C      IF (IEND.LT.0) WRITE (11,'(2I6)') NMESH, NWRT5
C      IF (NPRINT.EQ.1) THEN
C         WRITE (11,10001) NMOD, AGE
C         DO KK = 1, NMESH
C            WRITE (11,10002) (SX(IPX(J),KK+1), J=1, NWRT5)
C         END DO
C      END IF
10001    FORMAT (I6,1P,E17.9,0P)
10002    FORMAT (1P,E13.6,4E11.4,25E11.3,0P)
*
C Print the short summary, for every NWRT4'th model
         IF (MOD(NMOD,NWRT4).NE.0) RETURN
         BMS = BM/MSUN
         IF (IMODE.EQ.1) THEN
            BM = BM - (WINDML(1)+FAKEWIND(1))*DT
            BMS = BM/MSUN
         END IF
         FR = DLOG10(PX(17))
         FL = DLOG10(PX(18))
* print at max He energy generation, rather than Tmax
* Restore to Tmax
         KHE = NMESH + 2 - KMX(2)
         SDC = DLOG10(SX(3,2))
         SDM = DLOG10(SX(3,KTM))
         SDS = DLOG10(PX(3))
         STC = DLOG10(SX(4,2))
         STM = DLOG10(SX(4,KTM))
         STS = DLOG10(PX(4))
         PER0 =  1.5D0*FR - 0.5D0*DLOG10(8.157D0*PX(9))
         WMI = DLOG10(SX(9,KCE)*SX(17,KCE)**2*SX(23,KCE))
C RJS - Viscous mesh
         IF (IVMC.EQ.1) THEN
C            MWT = DMAX1(0d0, 0.75*(DTY/1d-4 - 1))
            MWT = DMAX1(0d0, 0.25*(DTY/1d-4 - 1))
            MWT = DEXP(-DSQRT(MWT))
            MWT = DMAX1(0.1d0,DMIN1(1d0,MWT))
         ELSE
            MWT = 0d0
         END IF
C Timestep control for TP-AGB
C Switch on mechanism for dealing with dredge up when necessary
         VLEC(ISTAR) = VLE
         LHEF = (VLEC(ISTAR) - VLEP(ISTAR))/VLEP(ISTAR)
         IF (IAGB.EQ.1) THEN
            IF (IDREDGE.EQ.0.AND.VLEC(ISTAR).GT.5d5.AND.LHEF.LT.0d0) IDREDGE = 1
C Find out where next to outer conv boundary is
            IF (IDREDGE.EQ.1) THEN
               MAXB = 0d0
               DO II = 1,12
                  IF (MCB(II).GT.MAXB.AND.MCB(II).LT.1d0) THEN
                     MAXB = MCB(II)
                  END IF
               END DO
               IDREDGE = 2
            END IF
C If boundary moves in, reduce DD
            IF (IDREDGE.EQ.2.AND.VLEC(ISTAR).LT.1d5) THEN
               MAXB2 = 0d0
               DO II = 1,12
                  IF (MCB(II).GT.MAXB2.AND.MCB(II).LT.1d0) MAXB2 = MCB(II)
               END DO
               IF (MAXB2.LT.MAXB) IDREDGE = 3
            END IF
            IF (IDREDGE.EQ.3) THEN
               IF (VLEC(ISTAR).GT.2d4) THEN
                  DD = 5.0
               ELSE
                  DD = 1.0
               END IF
            END IF
            IF (IDREDGE.EQ.3.AND.VLEC(ISTAR).LT.3d3) THEN
               DD = 5.0
               IDREDGE = 0
            END IF
            KICZ = KCB(3)
            KBICZ = KMX(2)
         END IF
         WRITE(32+20*(ISTAR-1),99003) NMOD,PX(9),DTY,TN,PER,VLH,VLT,SX(10,2),SX(11,2),
     :        SX(12,2),SX(13,2),SX(14,2),SX(15,2),SX(16,2),SX(1,2),SDC,STC,
     :        'cntr',VMH,AGE,TKH(ISTAR),RLF,VLE,VLN,XH,XHE,XC,XN,XO,XNE,XHE3,PX(1),
     :        SDS,STS,'srfc',VME,WMH,WMHE,BMS,DMT,VLC,SX(9,KTM),SX(10,KTM),
     :        SX(11,KTM),SX(12,KTM),SX(13,KTM),SX(14,KTM),SX(15,KTM),SX(16,KTM),
     :        SX(1,KTM),SDM,STM,'Tmax',VMC,MCB,WMI,PX(23),PER0,FR,FL,NMOD,
     :        KMH,KME,KCB,VMX,MEX,EMAX,KMX,KEX
         CALL FLUSH (32+20*(ISTAR-1))
         CALL FLUSH (36+20*(ISTAR-1))
         CALL FLUSH (37+20*(ISTAR-1))
         CALL FLUSH (38+20*(ISTAR-1))
C Find convective envelope base
         MENV(ISTAR) = 0d0
         RENV(ISTAR) = 0d0
        DO K=NMESH-1,1,-1 !surf to centre
            IF (SX(8,K+1).GT.0d0.AND.SX(8,K).LT.0d0.AND.SX(9,K).GT.MENV(ISTAR)/MSUN
     :          .AND.K.LT.NMESH-20) THEN
               MENV(ISTAR) = MSUN*(SX(9,K)*SX(8,K+1)-SX(9,K+1)*SX(8,K))/(SX(8,K+1)-SX(8,K))
               RENV(ISTAR) = RSUN*(SX(17,K)*SX(8,K+1)-SX(17,K+1)*SX(8,K))/(SX(8,K+1)-SX(8,K))
               KENV(ISTAR) = NMESH - K
            END IF
         END DO
C Update current values for extra timestep control variables
         VLHC(ISTAR) = VLH
         VLEC(ISTAR) = VLE
         VLCC(ISTAR) = VLC
         RLFC(ISTAR) = RLF
         TOTMC(ISTAR) = TM(ISTAR)
      END DO
      RETURN
      END
