      subroutine imlc(lcparms,spotparms)
c
c  Modified for use in BINAFIT by C. Hummel.
c  Comment reads and get input from lcparms.
c  Write to outfile 7 instead of to screen 6.
c  Program written by R. E. Wilson, Florida
c
c  Main program for computing light and radial velocity curves                10
c                                                                             20
c  Version of May 1, 1992                                                     30
C                                                                             40
C     TO PRINT VELOCITIES IN KM/SEC, SET VUNIT=1.                             50
C     TO PRINT NORMALIZED VELOCITIES IN SAME COLUMNS, SET VUNIT EQUAL TO      60
C     DESIRED VELOCITY UNIT IN KM/SEC.                                        70
C                                                                             80
C     PARAMETER PSHIFT IS DEFINED AS THE PHASE AT WHICH PRIMARY               90
C     CONJUNCTION (STAR 1 AWAY FROM OBSERVER) WOULD OCCUR IF THE             100
C     ARGUMENT OF PERIASTRON WERE 90 DEGREES. SINCE THE NOMINAL VALUE        110
C     OF THIS QUANTITY IS ZERO, PSHIFT MAY BE USED TO INTRODUCE AN           120
C     ARBITRARY PHASE SHIFT.                                                 130
C                                                                            140
      dimension rad(4),drdo(4),xtha(4),xfia(4),po(2)                         150
      dimension rv(3011),grx(3011),gry(3011),grz(3011),rvq(3011),            160
     $grxq(3011),gryq(3011),grzq(3011),slump1(3011),slump2(3011),            170
     $fr1(3011),fr2(3011),glump1(3011),glump2(3011),xx1(3011),               180
     $xx2(3011),yy1(3011),yy2(3011),zz1(3011),zz2(3011),grv1(3011),          190
     $grv2(3011),rftemp(3011),rf1(3011),rf2(3011),csbt1(3011),               200
     $csbt2(3011),gmag1(3011),gmag2(3011),hld(3200),snfi(6400),              210
     $csfi(6400),tld(6400),snth(260),csth(260),theta(520),rho(520),          220
     $aa(07),bb(07),mmsave(124)                                              230
      DIMENSION XLAT(2,100),xlong(2,100)                                     240
      COMMON /FLVAR/ PER,PSHIFT,DP,DS,EF,EFC,ECOS,PERR,PHPER,PCONJ,          250
     $PHPERI,VSUM1,VSUM2,VRA1,VRA2,VKM1,VKM2,VUNIT,vfvu,trc                  260
      COMMON /DPDX/ DPDX1,DPDX2,PHSV,PCSV                                    270
      COMMON /ECCEN/ E,A,PERIOD,VGA,SINI,VF,VFAC,VGAM,VOL1,VOL2,IFC          280
      COMMON /KFAC/ KFF1,KFF2,kfo1,kfo2                                      290
      COMMON /INVAR/ KH,IPB,IRTE,NREF,IRVOL1,IRVOL2,mref,ifsmv1,ifsmv2,      300
     $icor1,icor2,ld                                                         310
      COMMON /SPOTS/SINLAT(2,100),COSLAT(2,100),SINLNG(2,100),COSLNG         320
     $(2,100),RADSP(2,100),temsp(2,100),xlng(2,100)                          330
      COMMON /NSPT/ NSP1,NSP2,IFAT1,IFAT2                                    340
      common /qrm/ rm

c     These three lines added by C. Hummel
      real*4 lcparms(54),spotparms(4,2,100)
      logical first      
      data first/.true./ 

      data xtha(1),xtha(2),xtha(3),xtha(4),xfia(1),xfia(2),xfia(3),          350
     $xfia(4)/0.,1.570796,1.570796,1.570796,0.,0.,1.570796,3.141593/         360
   96 FORMAT('   Phase   lite 1  lite 2  tot lt  nrm lt    dist   mag+K      370
     $  VR1    VR2    DV1    DV2  r1pol  r1pt r1sid r1bak r2pol  r2pt r2     380
     $sid r2bak')                                                            390
   45 FORMAT('   Phase    light 1    light 2    (1+2+3)  norm lite     d     400
     $ist    mag+K     V RAD 1   V RAD 2    del V1    del V2  V km/s 1       410
     $V km/s 2')                                                             420
   65 FORMAT('   Phase    light 1    light 2    (1+2+3)  norm lite     d     430
     $ist   mag+K    V RAD 1   V RAD 2    del V1    del V2   VNORM 1  VN     440
     $ORM 2')                                                                450
   97 FORMAT(1X,F8.4,6F8.4,4F7.3,8(1X,F5.4))                                 460
   47 FORMAT(' wv lth    L1       L2      x1     x2     y1     y2     el     470
     $3    m zero  factor')                                                  480
   48 FORMAT('  ecc  lng per   s-m axis     F1       F2      V gam  Pshi     490
     $ft    Incl      g1      g2    Nspot1    Nspot2')                       500
   54 FORMAT('    T1      T2     Alb 1   Alb 2    Pot 1    Pot 2    M2/M     510
     $1   x1(bolo)  x2(bolo)  y1(bolo)  y2(bolo)')                           520
   33 FORMAT(3X,I2,4X,I1,5X,I1,5X,I1,4X,I2,2X,I2,3X,E10.5,1X,F5.4,2X,        530
     $F8.3,2X,E14.4,2X,F4.3,3f10.4)                                          540
   74 FORMAT(' DIMENSIONLESS RADIAL VELOCITIES CONTAIN FACTOR P/(2PI*A)'     550
     $)                                                                      560
   46 FORMAT(4X,'GRID1/4    GRID2/4',7X,'POLAR SBR 1',8X,'POLAR SBR 2'       570
     $,7X,'SURF. AREA 1',7X,'SURF. AREA 2',6X,'PERI. PH.',6X,'CONJ. PH.      580
     $')                                                                     590
   50 FORMAT(40H PRIMARY COMPONENT EXCEEDS CRITICAL LOBE)                    600
   51 FORMAT(42H SECONDARY COMPONENT EXCEEDS CRITICAL LOBE)                  610
   43 FORMAT(1H1)                                                            620
   42 FORMAT(1H )                                                            630
   41 FORMAT(1X,9HCOMPONENT,3X,6HR POLE,8X,5HDERIV,8X,7HR POINT,8X,          640
     $5HDERIV,9X,6HR SIDE,8X,5HDERIV,9X,6HR BACK,8X,5HDERIV)                 650
   80 FORMAT(1X,F6.4,1X,F7.5)                                                660
    2 FORMAT(F5.4,1X,F6.2,1X,F7.3,1X,2(F6.3,1X),F7.4,1X,F5.4,1X,f7.3,        670
     $2(1X,F6.3))                                                            680
    6 FORMAT(2(F6.4,1X),2(F6.3,1X),2(F7.4,1X),F8.5,4(1X,F5.3))               690
    3 FORMAT(1X,F8.4,1X,4(F9.6,2X),2(F7.4,1x),4F10.6,2(1X,E10.4))            700
    1 FORMAT(I2,1X,I1,1X,I1,1X,I1,1X,I2,1X,I2,1X,E9.5,1X,F5.4,1X,F7.3,       710
     $1X,F4.3,3(1x,f9.4))                                                    720
    4 FORMAT(F6.4,2(1X,F7.4),4(1X,F5.3),1X,F6.4,1X,F6.3,1X,F6.4)             730
    5 FORMAT(1X,F5.4,1X,F7.2,2X,F9.4,2(2X,F7.4),2X,F7.4,2X,F6.4,2X,f7.3,     740
     $2(2X,F6.3),2(5X,I3))                                                   750
   34 FORMAT(1X,F6.4,2(2X,F7.4),4(2x,f5.3),2X,F6.4,2X,F6.3,2X,F6.4)          760
   49 FORMAT(' PROGRAM SHOULD NOT BE USED IN MODE 1 OR 3 WITH NON-ZERO E     770
     $CCENTRICITY')                                                          780
   10 FORMAT(' MODE   IPB  IFAT1 IFAT2  N1  N2    Period     THE    VUNI     790
     $T(km/s)   V FAC       PHN    PHstrt    PHstop     Phin')               800
  148 format('   ifrad  nref   mref   ifsmv1   ifsmv2   icor1   icor2        810
     $ld')                                                                   820
  149 format(5x,i1,2(6x,i1),7x,i1,8x,i1,i9,i8,i6)                            830
   40 FORMAT(4X,I1,8(2X,F12.5))                                              840
   94 FORMAT(1X,2(I8,2X),F20.9,F19.9,F19.9,F18.9,F14.5,F15.5)                850
    8 FORMAT(2(2X,F6.4),2(2X,F6.3),2(2X,F7.4),2X,F8.5,2f9.3,2f10.3)          860
   84 FORMAT(1X,I4,4F12.3)                                                   870
   85 FORMAT(3F7.2,F8.3)                                                     880
   83 FORMAT(1X,'STAR  CO-LATITUDE  LONGITUDE  SPOT RADIUS  TEMP. FACTOR     890
     $')                                                                     900
  150 format(' Star     M/Msun   (Mean Radius)/Rsun   M Bol    Log g (cg     910
     $s)')                                                                   920
  250 format(4x,I1,4x,f6.2,11x,f7.2,7x,f5.2,5x,f5.2)                         930
  350 format(' Primary star exceeds outer contact surface')                  940
  351 format(' Secondary star exceeds outer contact surface')                950
   22 format(8(i1,1x))                                                       960
      ot=1./3.                                                               970
      KH=17                                                                  980
      PI=3.141593                                                            990
      twopi=pi+pi                                                           1000
      DTR=1.745329E-2                                                       1010

ch    The following added by C. Hummel
      if(first)open(7,file='imlcout.dat')
      ifrad=lcparms(1)      
      nref=lcparms(2)
      mref=lcparms(3)   
      ifsmv1=lcparms(4) 
      ifsmv2=lcparms(5) 
      icor1=lcparms(6)  
      icor2=lcparms(7)  
      ld=lcparms(8)     
      MODE=lcparms(9)   
      IPB=lcparms(10)   
      IFAT1=lcparms(11) 
      IFAT2=lcparms(12) 
      N1=lcparms(13)    
      N2=lcparms(14)    
      PERIOD=lcparms(15)
      THE=lcparms(16)   
      VUNIT=lcparms(17) 
      PHN=lcparms(18)   
      PHSTRT=lcparms(19)
      PHSTOP=lcparms(20)
      PHIN=lcparms(21)  
      E=lcparms(22)     
      PER=lcparms(23)   
      A=lcparms(24)     
      F1=lcparms(25)    
      F2=lcparms(26)    
      VGA=lcparms(27)   
      PSHIFT=lcparms(28)
      XINCL=lcparms(29) 
      GR1=lcparms(30)   
      GR2=lcparms(31)   
      tavh=lcparms(32)  
      tavc=lcparms(33)  
      alb1=lcparms(34)  
      alb2=lcparms(35)  
      poth=lcparms(36)  
      potc=lcparms(37)  
      rm=lcparms(38)    
      xbol1=lcparms(39) 
      xbol2=lcparms(40) 
      ybol1=lcparms(41) 
      ybol2=lcparms(42) 
      WL=lcparms(43)    
      HLUM=lcparms(44)  
      CLUM=lcparms(45)  
      XH=lcparms(46)    
      xc=lcparms(47)    
      yh=lcparms(48)    
      yc=lcparms(49)    
      EL3=lcparms(50)   
      ZERO=lcparms(51)  
      FACTOR=lcparms(52)
      NSP1=lcparms(53)
      NSP2=lcparms(54)

ch
ch    DO 1000 IT=1,1000                                                     1020
ch    WRITE(7,43)                                                           1030
ch    read(5,22)ifrad,nref,mref,ifsmv1,ifsmv2,icor1,icor2,ld                1040
      if(ifrad.eq.9) stop                                                   1050
ch    READ(5,1)MODE,IPB,IFAT1,IFAT2,N1,N2,PERIOD,THE,VUNIT,PHN,PHSTRT,      1060
ch   $PHSTOP,PHIN                                                           1070
ch    READ(5,2)E,PER,A,F1,F2,VGA,PSHIFT,XINCL,GR1,GR2                       1080
ch    read(5,6) tavh,tavc,alb1,alb2,poth,potc,rm,xbol1,xbol2,ybol1,         1090
ch   $ybol2                                                                 1100
ch    READ(5,4)WL,HLUM,CLUM,XH,xc,yh,yc,EL3,ZERO,FACTOR                     1110
ch    NSP1=0                                                                1120
ch    NSP2=0                                                                1130
      DO 88 KP=1,2                                                          1140
ch    DO 87 I=1,100                                                         1150
      if(kp.eq.1)nsp=nsp1  ! added by C. Hummel
      if(kp.eq.2)nsp=nsp2  ! added by C. Hummel
      do 87 i=1,nsp        ! added by C. Hummel
ch    READ(5,85)XLAT(KP,I),XLONG(KP,I),RADSP(KP,I),TEMSP(KP,I)              1160
      xlat(kp,i)=spotparms(1,kp,i)    ! added by C. Hummel
      xlong(kp,i)=spotparms(2,kp,i)   ! added by C. Hummel
      radsp(kp,i)=spotparms(3,kp,i)   ! added by C. Hummel
      temsp(kp,i)=spotparms(4,kp,i)   ! added by C. Hummel
      xlng(kp,i)=dtr*xlong(kp,i)                                            1170
      IF(XLAT(KP,I).GE.200.) GOTO 88                                        1180
      SINLAT(KP,I)=SIN(DTR*XLAT(KP,I))                                      1190
      COSLAT(KP,I)=COS(DTR*XLAT(KP,I))                                      1200
      SINLNG(KP,I)=SIN(DTR*XLONG(KP,I))                                     1210
      COSLNG(KP,I)=COS(DTR*XLONG(KP,I))                                     1220
ch    IF(KP.EQ.1)NSP1=NSP1+1                                                1230
ch 87 IF(KP.EQ.2)NSP2=NSP2+1                                                1240
   87 continue  ! added by C. Hummel
   88 CONTINUE                                                              1250
      DINT1=PI*(1.-XBOL1/3.+2.*ybol1/9.)                                    1260
      DINT2=PI*(1.-XBOL2/3.+2.*ybol2/9.)                                    1270
      NSTOT=NSP1+NSP2                                                       1280
      NP1=N1+1                                                              1290
      NP2=N1+N2+2                                                           1300
      PERR=DTR*PER                                                          1310
      IRTE=0                                                                1320
      IRVOL1=0                                                              1330
      IRVOL2=0                                                              1340
      CALL SINCOS(1,N1,N1,SNTH,CSTH,SNFI,CSFI,MMSAVE)                       1350
      CALL SINCOS(2,N2,N1,SNTH,CSTH,SNFI,CSFI,MMSAVE)                       1360
      CALL modlog(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,     1370
     $rm,poth,potc,gr1,gr2,alb1,alb2,n1,n2,f1,f2,mod,xincl,the,mode,        1380
     $snth,csth,snfi,csfi,grv1,grv2,xx1,yy1,zz1,xx2,yy2,zz2,glump1,         1390
     $glump2,csbt1,csbt2,gmag1,gmag2)                                       1400
      CALL VOLUME(VOL1,RM,POTH,DP,F1,N1,N1,1,RV,GRX,GRY,GRZ,RVQ,            1410
     $GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMMD,SMD,      1420
     $GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,          1430
     $GMAG1,GMAG2,GR1,1)                                                    1440
      CALL VOLUME(VOL2,RM,POTC,DP,F2,N2,N1,2,RV,GRX,GRY,GRZ,RVQ,            1450
     $GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMMD,SMD,      1460
     $GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,          1470
     $GMAG1,GMAG2,GR2,1)                                                    1480
      if(e.eq.0.) goto 117                                                  1490
      DAP=1.+E                                                              1500
      P1AP=POTH-2.*E*RM/(1.-E*E)                                            1510
      VL1=VOL1                                                              1520
      CALL VOLUME(VL1,RM,P1AP,DAP,F1,N1,N1,1,RV,GRX,GRY,GRZ,RVQ,            1530
     $GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMMD,SMD,      1540
     $GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,          1550
     $GMAG1,GMAG2,GR1,2)                                                    1560
      DPDX1=(POTH-P1AP)*(1.-E*E)*.5/E                                       1570
      P2AP=POTC-2.*E/(1.-E*E)                                               1580
      VL2=VOL2                                                              1590
      CALL VOLUME(VL2,RM,P2AP,DAP,F2,N2,N1,2,RV,GRX,GRY,GRZ,RVQ,            1600
     $GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,SNTH,CSTH,SNFI,CSFI,SUMMD,SMD,      1610
     $GRV1,GRV2,XX1,YY1,ZZ1,XX2,YY2,ZZ2,CSBT1,CSBT2,GLUMP1,GLUMP2,          1620
     $GMAG1,GMAG2,GR2,2)                                                    1630
      DPDX2=(POTC-P2AP)*(1.-E*E)*.5/E                                       1640
  117 CONTINUE                                                              1650
      PHSV=POTH                                                             1660
      PCSV=POTC                                                             1670
      IF(E.EQ.0.) GOTO 61                                                   1680
      IF(MOD.EQ.1) WRITE(7,49)                                              1690
   61 CONTINUE                                                              1700
      CALL BBL(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,        1710
     $SLUMP1,SLUMP2,THETA,RHO,AA,BB,POTH,POTC,N1,N2,F1,F2,D,HLUM            1720
     $,clum,xh,xc,yh,yc,gr1,gr2,wl,sm1,sm2,tpolh,tpolc,sbrh,sbrc,ifat1,     1730
     $ifat2,tavh,tavc,alb1,alb2,xbol1,xbol2,ybol1,ybol2,phn,rm,xincl,       1740
     $hot,cool,snth,csth,snfi,csfi,tld,glump1,glump2,xx1,xx2,yy1,yy2,       1750
     $zz1,zz2,dint1,dint2,grv1,grv2,rftemp,rf1,rf2,csbt1,csbt2,gmag1,       1760
     $gmag2,mode)                                                           1770
      KH=0                                                                  1780
      if(kfo1.eq.0) goto 380                                                1790
      write(7,350)                                                          1800
      goto 381                                                              1810
  380 IF(KFF1.EQ.1) WRITE(7,50)                                             1820
  381 if(kfo2.eq.0) goto 382                                                1830
      write(7,351)                                                          1840
      goto 383                                                              1850
  382 IF(KFF2.EQ.1) WRITE(7,51)                                             1860

  383 if(first)then  ! added by C. Hummel
      IF((KFF1+KFF2+kfo1+kfo2).GT.0) WRITE(7,42)                            1870
      write(7,42)                                                           1880
      write(7,148)                                                          1890
      write(7,149) ifrad,nref,mref,ifsmv1,ifsmv2,icor1,icor2,ld             1900
      write(7,42)                                                           1910
      WRITE(7,10)                                                           1920
      WRITE(7,33)MODE,IPB,IFAT1,IFAT2,N1,N2,PERIOD,THE,VUNIT,VFAC,PHN,      1930
     $PHstrt,PHstop,Phin                                                    1940
      WRITE(7,42)                                                           1950
      WRITE(7,48)                                                           1960
      WRITE(7,5)E,PER,A,F1,F2,VGA,PSHIFT,XINCL,GR1,GR2,NSP1,NSP2            1970
      WRITE(7,42)                                                           1980
      WRITE(7,54)                                                           1990
      WRITE(7,8)TAVH,TAVC,ALB1,ALB2,PHSV,PCSV,RM,XBOL1,xbol2,ybol1,ybol2    2000
      WRITE(7,42)                                                           2010
      WRITE(7,47)                                                           2020
      WRITE(7,34)WL,HLUM,CLUM,XH,XC,yh,yc,el3,ZERO,FACTOR                   2030
      WRITE(7,42)                                                           2040
      endif ! added by C. Hummel

      IF(NSTOT.GT.0.and.first) WRITE(7,83) ! modified by C. Hummel          2050
      DO 188 KP=1,2                                                         2060
      IF((NSP1+KP-1).EQ.0) GOTO 188                                         2070
      IF((NSP2+(KP-2)**2).EQ.0) GOTO 188                                    2080
      NSPOT=NSP1                                                            2090
      IF(KP.EQ.2) NSPOT=NSP2                                                2100

      if(first)then    ! added by C. Hummel
      DO 187 I=1,NSPOT                                                      2110
  187 WRITE(7,84)KP,XLAT(KP,I),XLONG(KP,I),RADSP(KP,I),TEMSP(KP,I)          2120
      endif            ! added by C. Hummel

  188 if(first)WRITE(7,42) ! modified by C. Hummel                          2130
      if(first)write(7,150)! modified by C. Hummel                          2140

      rr1=.6203505*vol1**ot                                                 2150
      rr2=.6203505*vol2**ot                                                 2160
      tav1=10000.*tavh                                                      2170
      tav2=10000.*tavc                                                      2180
      call mlrg(a,period,rm,rr1,rr2,tav1,tav2,sms1,sms2,sr1,sr2,            2190
     $bolm1,bolm2,xlg1,xlg2)                                                2200
      ns1=1                                                                 2210
      ns2=2                                                                 2220

      if(first)then ! added by C. Hummel
      write(7,250) ns1,sms1,sr1,bolm1,xlg1                                  2230
      write(7,250) ns2,sms2,sr2,bolm2,xlg2                                  2240
      write(7,42)                                                           2250
      WRITE(7,46)                                                           2260
      WRITE(7,94) MMSAVE(NP1),MMSAVE(NP2),SBRH,SBRC,SM1,SM2,PHPERI,         2270
     $PCONJ                                                                 2280
      WRITE(7,42)                                                           2290
      endif ! added by C. Hummel

      ALL=HOT+COOL+EL3                                                      2300
      IF(MODE.EQ.-1) ALL=COOL+EL3                                           2310
      IF(IFRAD.EQ.0) GOTO 71                                                2320

      if(first)then ! added by C. Hummel
      WRITE(7,96)                                                           2330
      endif         ! added by C. Hummel

      GOTO 77                                                               2340
   71 IF(VUNIT.EQ.1..and.first) WRITE(7,45)  ! modified by C. Hummel        2350
      IF(VUNIT.NE.1..and.first) WRITE(7,65)  ! modified by C. Hummel        2360
   77 CONTINUE                                                              2370
      LL1=MMSAVE(N1)+1                                                      2380
      NPP2=NP2-1                                                            2390
      LL2=MMSAVE(NPP2)+1                                                    2400
      LLL1=MMSAVE(NP1)                                                      2410
      LLL2=MMSAVE(NP2)                                                      2420
      LLLL1=(LL1+LLL1)/2                                                    2430
      LLLL2=(LL2+LLL2)/2                                                    2440
      POTH=PHSV                                                             2450
      POTC=PCSV                                                             2460
      PO(1)=POTH                                                            2470
      PO(2)=POTC                                                            2480
      IF(E.EQ.0.) IRVOL1=1                                                  2490
      IF(E.EQ.0.) IRVOL2=1                                                  2500
      IF(E.EQ.0.) IRTE=1                                                    2510
      do 20,phas=phstrt,phstop,phin                                         2520
      CALL BBL(RV,GRX,GRY,GRZ,RVQ,GRXQ,GRYQ,GRZQ,MMSAVE,FR1,FR2,HLD,        2530
     $SLUMP1,SLUMP2,THETA,RHO,AA,BB,POTH,POTC,N1,N2,F1,F2,D,hlum,           2540
     $clum,xh,xc,yh,yc,gr1,gr2,wl,sm1,sm2,tpolh,tpolc,sbrh,sbrc,ifat1,      2550
     $ifat2,tavh,tavc,alb1,alb2,xbol1,xbol2,ybol1,ybol2,phas,rm,xincl,      2560
     $hot,cool,snth,csth,snfi,csfi,tld,glump1,glump2,xx1,xx2,yy1,yy2,       2570
     $zz1,zz2,dint1,dint2,grv1,grv2,rftemp,rf1,rf2,csbt1,csbt2,gmag1,       2580
     $gmag2,mode)                                                           2590
      HTT=HOT                                                               2600
      IF(MODE.EQ.-1) HTT=0.                                                 2610
      TOTAL=HTT+COOL+EL3                                                    2620
      TOTALL=TOTAL/ALL                                                      2630
      TOT=TOTALL*FACTOR                                                     2640
      SMAGG=-1.085736*ALOG(TOTALL)+ZERO                                     2650
      IF(IFRAD.EQ.0) GOTO 82                                                2660

      if(first)then ! added by C. Hummel
      WRITE(7,97)PHAS,HTT,COOL,TOTAL,TOT,D,SMAGG,VSUM1,VSUM2,VRA1,VRA2,     2670
     $RV(1),RV(LL1),RV(LLLL1),RV(LLL1),RVQ(1),RVQ(LL2),RVQ(LLLL2),RVQ(LL    2680
     $L2)                                                                   2690
      endif         ! added by C. Hummel

      GOTO 20                                                               2700

   82 if(first)then ! added by C. Hummel
      WRITE(7,3)PHAS,HTT,COOL,TOTAL,TOT,D,SMAGG,VSUM1,VSUM2,VRA1,VRA2,      2710
     $VKM1,VKM2                                                             2720
      endif         ! added by C. Hummel

   20 CONTINUE                                                              2730
      if(first)then ! added by C. Hummel
      WRITE(7,42)                                                           2740
      WRITE(7,41)                                                           2750
      WRITE(7,42)                                                           2760
      endif         ! added by C. Hummel

      do 119 ii=1,2                                                         2770
      gt1=2-ii                                                              2780
      gt2=ii-1                                                              2790
      f=f1*gt1+f2*gt2                                                       2800
      do 118 i=1,4                                                          2810
      call romqsp(po(ii),rm,f,dp,e,xtha(i),xfia(i),rad(i),drdo(i),          2820
     $drdq,dodq,ii,mode)                                                    2830
  118 continue                                                              2840

      if(first)then ! added by C. Hummel
      write(7,40) ii,rad(1),drdo(1),rad(2),drdo(2),rad(3),drdo(3),          2850
     $rad(4),drdo(4)                                                        2860
      endif         ! added by C. Hummel

  119 continue                                                              2870

      if(first)then ! added by C. Hummel
      WRITE(7,42)                                                           2880
      WRITE(7,74)                                                           2890
      endif         ! added by C. Hummel

c1000 CONTINUE      ! no loop, modified by C. Hummel                        2900
ch    STOP                                                                  2910
      if(first)close(7)
      first=.false.
      return 
      END                                                                   2920

