c
c  This program computes an opacity table in the same format as the OPAL
c  opacity tables but at an arbitrary compomposition.
c
c  Compile by, for example,
c
c    f77 -o z_xcotrin_ZXtable z_xcotrin_ZXtable.f z14xcotrin21.f
c
c  and then run it via
c
c    ./z_xcotrin_ZXtable
c
c  It will ask you for the input values needed to read and calculate opacities.
c
c  Note that, due to the smoothing recommended by OPAL when reading in the
c  opacities, any tables created will not be exactly identical to those in
c  the input OPAL file (e.g., 'GN93hz') even at the same composition.
c
c  (Note that internal real variables are single precision.)
c
      program z_xcotrin_ZXtable
c     -------------------------
c
      parameter ( ntopal = 70, nropal = 19, ntferg = 69 )
c
      character*255 cf_hz, cf_ofe, cbeg_ferg, cf_d, fnalt
      character*255 cdirin, cdir_mol, cdir_cond, cf_used
      character*1 cyn
c
      logical lxst
c
      dimension xiz(19), fninz(19), xmet(5)
c
      dimension fltopal(ntopal), fltferg(ntferg)
c
      data fltopal / 3.75, 3.80, 3.85, 3.90, 3.95,
     $     4.00, 4.05, 4.10, 4.15, 4.20,  4.25, 4.30, 4.35, 4.40, 4.45,
     $     4.50, 4.55, 4.60, 4.65, 4.70,  4.75, 4.80, 4.85, 4.90, 4.95,
     $     5.00, 5.05, 5.10, 5.15, 5.20,  5.25, 5.30, 5.35, 5.40, 5.45,
     $     5.50, 5.55, 5.60, 5.65, 5.70,  5.75, 5.80, 5.85, 5.90, 5.95,
     $     6.00, 6.10, 6.20, 6.30, 6.40,  6.50, 6.60, 6.70, 6.80, 6.90,
     $     7.00, 7.10, 7.20, 7.30, 7.40,  7.50, 7.60, 7.70, 7.80, 7.90,
     $     8.00, 8.10, 8.30, 8.50, 8.70 /
c
      data fltferg / 3.70, 3.65, 3.60, 3.55, 3.50,
     $     3.49, 3.48, 3.47, 3.46, 3.45,  3.44, 3.43, 3.42, 3.41, 3.40,
     $     3.39, 3.38, 3.37, 3.36, 3.35,  3.34, 3.33, 3.32, 3.31, 3.30,
     $     3.29, 3.28, 3.27, 3.26, 3.25,  3.24, 3.23, 3.22, 3.21, 3.20,
     $     3.19, 3.18, 3.17, 3.16, 3.15,  3.14, 3.13, 3.12, 3.11, 3.10,
     $     3.09, 3.08, 3.07, 3.06, 3.05,  3.04, 3.03, 3.02, 3.01, 3.00,
     $     2.99, 2.98, 2.97, 2.96, 2.95,  2.94, 2.93, 2.92, 2.91, 2.90,
     $     2.85, 2.80, 2.75, 2.70 /
c
      call ask_z_limits( nzmax, zmin, zmax )
c
      iu = 7
      i_mol = 1
      i_cond = 1
      i_cno = 1
      cf_hz = 'AGS04'
      cf_ofe = ' '
      cbeg_ferg = ' '
      cdirin = ' '
      cdir_mol = ' '
      cdir_cond = ' '
c
c  The following just asks the user to input the relevant variables
c
      write(6,10)
 10   format( / 'Create an OPAL-format opacity table using',
     $     ' opacity routines in z14xcotrin21.f' / )
c
      call ask_smooth( initsmooth, lowcosmooth, interpcosmooth )
c
      if ( initsmooth .ne. 0 ) write(6,11)
 11   format( / 'NOTE that, due to the smoothing',
     $     ' recommended by OPAL when reading in the' /
     $     'opacities, any tables created will not',
     $     ' be exactly identical to those in' /
     $     'the input OPAL file (e.g., ''GN93hz'') even',
     $     ' at the same composition.' / )
c
      write(6,12)
 12   format( / 'Note that this interactive example',
     $     ' program allows some of the more complex' /
     $     'opacity interpolation options, but does',
     $     ' not require that you use them.' / /
     $     'For example, CNO-interpolation (to account',
     $     ' for the effects of the nuclear' /
     $     'reactions that convert C and O to N,',
     $     ' and N to Ne) usually have little effect' /
     $     'on the opacity, at least at the temperatures',
     $     ' where these effects can be' /
     $     'computed: you can avoid this complication',
     $     ' by setting ''I_CNO = 0'' when asked.' /
     $     'If electron conduction at high density',
     $     ' is irrelevant to you, set ''I_COND = 0''.' )
      write(6,15)
 15   format( 'If molecular opacities at low temperature',
     $     ' are unnecessary, set ''I_MOL = 0'';' /
     $     'and even if you need them, the program''s',
     $     ' ''best choice'' (''I_MOL = 1'') is likely' /
     $     'to suffice (specific choices should',
     $     ' not be needed, let alone ''non-standard''' /
     $     'cversions).  Finally, when actually',
     $     ' computing the opacities, note that large' /
     $     'amounts of carbon and oxygen only result',
     $     ' from helium-burning reactions; for' /
     $     'the Sun (and for most stellar envelopes),',
     $     ' ''exC = 0.0'' and ''exO = 0.0'' are' /
     $     'appropriate (i.e., no C or O has been',
     $     ' added to the initial metallicity Z).' / /
     $     'You will be prompted for the',
     $     ' required variable values...' / )
c
      goto 30
 20   write(6,9999)
 9999 format('------ Input error ... try again:',$)
 30   write(6,40)
 40   format( / 'Do you wish to read in a previously-computed',
     $     ' opacity dump-file? (y/n): ', $ )
      read(5,'(a1)',err=20,end=20) cyn
      if ( cyn .ne. 'y' .and. cyn .ne. 'Y' .and.
     $     cyn .ne. 'n' .and. cyn .ne. 'N' ) goto 20
c
c  If we are not reading a binary opacity dump-file: must read the ASCCI files:
c
      if ( cyn .ne. 'y' .and. cyn .ne. 'Y' ) then
c
         if ( nzmax .ge. 14 ) then
            z_read = 0.02
         else
            write(6,170) nzmax
 170        format( / '[Note that the OPAL Z-range will be only',
     $           ' partially available, since you have' /
     $           'compiled a version of these routines with only',
     $           i3, ' allowed Z-storage values.]', $ )
            goto 190
 180        write(6,9999)
 190        write(6,200)
 200        format( / '  Give  Z_READ (a Z-value from 0.0 to 0.1:',
     $           ' opacities will be read' /
     $           '                in for as wide a range',
     $           'around Z_READ as possible): ', $ )
            read(5,*,err=180,end=180) z_read
            if ( z_read .lt. 0.0 .or. z_read .gt. 0.1 ) then
               write(6,9997) 'Bad Z_READ value'
 9997          format('=== ',a,' ... try again:',$)
               goto 190
            endif
         endif
c
         goto 220
 210     write(6,9999)
 220     write(6,230)
 230     format( / 'Which set of OPAL opacities (i.e.,',
     $        ' Z-composition file) do you wish to use?' /
     $        '  Give  CF_HZ (typically  GN93hz  ,',
     $        '  GS98hz  , or  AGS04hz): ', $ )
         read(5,'(a255)',err=210,end=210) cf_hz
         if ( cf_hz .eq. ' ' ) then
            write(6,9997) 'No filename was input'
            goto 220
         endif
c
         goto 250
 240     write(6,9999)
 250     write(6,260)
 260     format( / 'Do you wish to use a non-solar O/Fe',
     $        ' ratio (i.e., [O/Fe] non-zero)?' /
     $        '  Give  OFE_BRACKET ([O/Fe]: log',
     $        ' of relative ratio; 0.0 = solar): ' $ )
         read(5,*,err=240,end=240) ofe_bracket
         if ( ofe_bracket .lt. -0.5 .or. ofe_bracket .gt. 1.0 ) then
            write(6,9997) 'Bad OFE_BRACKET value'
            goto 250
         else if ( ofe_bracket .ne. 0.0 ) then
            goto 280
 270        write(6,9999)
 280        write(6,290)
 290        format( / 'Do you wish to use a non-default OPAL',
     $           ' opacity file for [O/Fe]-interpolation?' /
     $            ' Give  CF_OFE (this filename;',
     $           ' blank means use default): ', $ )
            read(5,'(a255)',err=270,end=270) cf_ofe
         endif
c
         goto 90
 80      write(6,9999)
 90      write(6,100)
 100     format( / 'Do you wish to use molecular opacities at low',
     $        ' temperatures?' / '  Give  I_MOL (0 = no, 1 = yes,',
     $        ' 21 = require GN93, 51 = require AGS04,' /
     $        '               991 = require a non-standard',
     $        ' set of molecular opacities): ', $ )
         read(5,*,err=80,end=80) i_mol
         if ( i_mol .ne. 0 .and. i_mol .ne. 1 .and. i_mol .ne. 21 .and.
     $        i_mol .ne. 31 .and. i_mol .ne. 41 .and.
     $        i_mol .ne. 51 .and. i_mol .ne. 61 .and.
     $        i_mol .ne. 71 .and. i_mol .ne. 991 ) then
            write(6,9997) 'Bad I_MOL value'
            goto 90
         else if ( i_mol .eq. 991 ) then
            goto 120
 110        write(6,9998)
 9998       format('------ Input error ... try again:')
 120        write(6,130)
 130        format( '  Give the non-standard molecular-opacity',
     $           ' file-name-beginning [NOTE THAT' /
     $           '     g  g98.  l03.  ags04.  s92.  s92ae.',
     $           '  ARE THE 6 STANDARD CASES]: ', $ )
            read(5,'(a255)',err=110,end=110) cbeg_ferg
            if ( cbeg_ferg .eq. ' ' ) then
               write(6,9997) 'No name was input'
               goto 120
            endif
         endif
c
         goto 60
 50      write(6,9999)
 60      write(6,70)
 70      format( / 'Do you wish to use conductive opacities at high',
     $        ' densities?' / '  Give  I_COND (0 = no, 1 = yes,',
     $        ' 2 = require Potekhin et al. 2006 ones): ', $ )
         read(5,*,err=50,end=50) i_cond
         if ( i_cond .lt. 0 .or. i_cond .gt. 2 ) then
            write(6,9997) 'I_COND'
            goto 60
         endif
c
         goto 150
 140     write(6,9999)
 150     write(6,160)
 160     format( / 'Will you wish to vary the CNO-makeup of the',
     $        ' metallicity Z (above 10^4 K)?' /
     $        '  Give  I_CNO (0 = no [usual case],',
     $        ' 1 = yes [more test case than useful]): ', $ )
         read(5,*,err=140,end=140) i_cno
         if ( i_cno .lt. 0 .or. i_cno .gt. 1 ) then
            write(6,9997) 'Bad I_CNO value'
            goto 150
         endif
c
         goto 310
 300     write(6,9999)
 310     write(6,320)
 320     format( / 'What directory contains the OPAL opacity',
     $        ' files (e.g, /home/user/opaldir/)?' /
     $        '  Give  CDIRIN (this directory name;',
     $        ' blank means the local directory):' )
         read(5,'(a255)',err=300,end=300) cdirin
c
         cdir_mol = cdirin
         cdir_cond = cdirin
c
         if ( i_mol .gt. 0 ) then
            goto 340
 330        write(6,9999)
 340        write(6,350)
 350        format( / 'Are the molecular opacity files in the',
     $           ' same directory as OPAL files? (y/n): ', $ )
            read(5,'(a1)',err=330,end=330) cyn
            if ( cyn .ne. 'y' .and. cyn .ne. 'Y' .and.
     $           cyn .ne. 'n' .and. cyn .ne. 'N' ) goto 330
            if ( cyn .ne. 'y' .and. cyn .ne. 'Y' ) then
               goto 370
 360           write(6,9998)
 370           write(6,380)
 380           format('  Give  CDIR_MOL (the directory name;',
     $              ' blank means the local directory):' )
               read(5,'(a255)',err=360,end=360) cdir_mol
            endif
         endif
c
         if ( i_cond .gt. 0 ) then
            goto 400
 390        write(6,9999)
 400        write(6,410)
 410        format( / 'Are the conductive opacity files in the',
     $           ' same directory as OPAL files? (y/n): ', $ )
            read(5,'(a1)',err=390,end=390) cyn
            if ( cyn .ne. 'y' .and. cyn .ne. 'Y' .and.
     $           cyn .ne. 'n' .and. cyn .ne. 'N' ) goto 390
            if ( cyn .ne. 'y' .and. cyn .ne. 'Y' ) then
               goto 430
 420           write(6,9998)
 430           write(6,440)
 440           format('  Give  CDIR_COND (the directory name;',
     $              ' blank means the local directory):' )
               read(5,'(a255)',err=420,end=420) cdir_cond
            endif
         endif
c
         goto 460
 450     write(6,9999)
 460     write(6,470)
 470     format( / 'How much logging-output do you wish',
     $        ' describing the opacities that are read in?' /
     $        '  Give  LIST_LEVEL (0 = none, 1 = print',
     $        ' first of each type of opacity file,' /
     $        '                    9999 = print EVERY',
     $        ' opacity file that is read in): ', $ )
         read(5,*,err=450,end=450) list_level
c
c  Now we call the subroutines that read in OPAL opacities from ASCII files ...
c
         write(6,1000)
 1000    format( / / 'Now we call the subroutines that',
     $        ' read in OPAL opacities from ASCII files:' /
     $        '  This may take a while ...' / )
c
         if ( list_level .gt. 0 ) then
c
c            write(6,1010) list_level
c 1010       format( / '      LIST_LEVEL =', i12 /
c     $           '      CALL SET_OPAL_LIST_LEVEL( LIST_LEVEL )' )
c
c  Set the amount of logging output, if there should be any
c
            call set_opal_list_level( list_level )
c
         endif
c
c         write(6,1020) cdirin(:max(1,lnblnk(cdirin)))
c 1020    format( / '      CDIRIN = ''', a, '''' /
c     $        '      CALL SET_OPAL_DIR( CDIRIN )' )
c
c  Set the directory where the OPAL opacity files are to be found
c
         call set_opal_dir( cdirin )
c
         if ( i_mol .gt. 0 .and. cdir_mol .ne. cdirin ) then
c
c            write(6,1030) cdir_mol(:max(1,lnblnk(cdir_mol)))
c 1030       format( / '      CDIR_MOL = ''', a, '''' /
c     $           '      CALL SET_MOL_DIR( CDIR_MOL )' )
c
c  If necessary, set the directory where molecular opacity files are found
c
            call set_mol_dir( cdir_mol )
c
         endif
c
         if ( i_cond .gt. 0 .and. cdir_cond .ne. cdirin ) then
c
            write(6,1040) cdir_cond(:max(1,lnblnk(cdir_cond)))
 1040       format( / '      CDIR_COND = ''', a, '''' /
     $           '     CALL SET_COND_DIR( CDIR_COND )' )
c
c  If necessary, set the directory where conductive opacity files are found
c
            call set_cond_dir( cdir_cond )
c
         endif
c
         if ( i_mol .eq. 991 ) then
c
c            write(6,1050) cbeg_ferg(:max(1,lnblnk(cbeg_ferg)))
c 1050       format( / '      CBEG_FERG = ''', a, '''' /
c     $           '     CALL SET_FERG_USER( CBEG_FERG )' )
c
c  If necessary, specify a non-standard molecular opacity case
c    (in general, you should not do this):
c
            call set_ferg_user( cbeg_ferg )
c
         endif
c
         if ( i_mol .eq. 0 .and.
     $        i_cond .eq. 0 .and. i_cno .eq. 0 ) then
c
c            write(6,1060) iu, z_read, cf_hz(:max(1,lnblnk(cf_hz))),
c     $           ofe_bracket, cf_ofe(:max(1,lnblnk(cf_ofe)))
c 1060       format( / '      IU =', i3 /
c     $           '      Z_READ =', f10.7 /
c     $           '      CF_HZ = ''', a, '''' /
c     $           '      OFE_BRACKET =', f10.6 /
c     $           '      CF_OFE = ''', a, '''' /
c     $           '      CALL READ_BASIC_OPAL_OPAC( IU,',
c     $           ' Z_READ, CF_HZ, OFE_BRACKET, CF_OFE )' / )
c
c  Read in the opal opacities using the basic interface, if that will work
c  (note that this will produce log output, if you have set LIST_LEVEL > 0):
c
            call read_basic_opal_opac( iu, z_read, cf_hz,
     $                                 ofe_bracket, cf_ofe )
c
         else
c
c            write(6,1070) iu, z_read, cf_hz(:max(1,lnblnk(cf_hz))),
c     $           ofe_bracket, cf_ofe(:max(1,lnblnk(cf_ofe))),
c     $           i_mol, i_cond, i_cno
c 1070       format( / '      IU =', i3 /
c     $           '      Z_READ =', f10.7 /
c     $           '      CF_HZ = ''', a, '''' /
c     $           '      OFE_BRACKET =', f10.6 /
c     $           '      CF_OFE = ''', a, '''' /
c     $           '      I_MOL =',i4 /
c     $           '      I_COND =',i3 /
c     $           '      I_CNO =',i4 /
c     $           '      CALL READ_EXTENDED_OPAC( IU,',
c     $           ' Z_READ, CF_HZ, OFE_BRACKET, CF_OFE,' /
c     $           '     $                         I_MOL,',
c     $           ' I_COND, I_CNO, '' '' )' / )
c
c  Otherwise, read in the opal opacities using the extended interface
c  (note that this will produce log output, if you have set LIST_LEVEL > 0):
c
            call read_extended_opac( iu, z_read, cf_hz,
     $                               ofe_bracket, cf_ofe,
     $                               i_mol, i_cond, i_cno, ' ' )
c
         endif
c
c  We are done reading in the OPAL opacities; ask if a dumpfile is desired.
c
         write(6,1200)
 1200    format( / 'We are now done reading in the OPAL opacities.' / )
c
         goto 1220
 1210    write(6,9999)
 1220    write(6,1230)
 1230    format( / 'Do you wish to create a binary',
     $        ' opacity-dumpfile for future use? (y/n): ', $ )
         read(5,'(a1)',err=1210,end=1210) cyn
         if ( cyn .ne. 'y' .and. cyn .ne. 'Y' .and.
     $        cyn .ne. 'n' .and. cyn .ne. 'N' ) goto 1210
c
         if ( cyn .eq. 'y' .or. cyn .eq. 'Y' ) then
c
            goto 1250
 1240       write(6,9999)
 1250       write(6,1260)
 1260       format( / '  Give  CF_D (the name of',
     $           ' the binary opacity dump-file):' )
            read(5,'(a255)',err=1240,end=1240) cf_d
            if ( cf_d .eq. ' ' ) then
               write(6,9997) 'No filename was input'
               goto 1250
            endif
c
            write(6,1270)
 1270       format( / 'We will now write the binary opacity',
     $           ' dump-file.  This should be fairly fast...' )
c
            write(6,1280) iu, cf_d(:max(1,lnblnk(cf_d)))
 1280       format( / '      IU =', i3 /
     $           '      CF_D = ''', a, '''' /
     $           '      CALL DUMP_OPAL_OPAC( IU, CF_D )' )
c
c  If so specified, write the binary opacity dump-file:
c
            call dump_opal_opac( iu, cf_d ) ;
c
            write(6,1290)
 1290       format( / 'We are done writing the',
     $           ' binary opacity dump-file.' / )
c
         endif
c
c  Below is the other branch of the IF statement, for the case when we should
c  read in opacities from a binary opacity dump-file (rather than ASCII files).
c
      else
c
c  Ask for the name of the opacity dump-file to read from.
c
         goto 1550
 1540    write(6,9999)
 1550    write(6,1560)
 1560    format( / '  Give  CF_D (the name of',
     $        ' the binary opacity dump-file):' )
         read(5,'(a255)',err=1540,end=1540) cf_d
         if ( cf_d .eq. ' ' ) then
            write(6,9997) 'No filename was input'
            goto 1550
         endif
c
         write(6,1700)
 1700    format( / / 'Now we call the subroutine that',
     $        ' reads in opacities from a binary dump-file:' /
     $        '  This should be fairly fast ...' / )
c
         write(6,1780) iu, cf_d(:max(1,lnblnk(cf_d)))
 1780    format( / '      IU =', i3 /
     $           '      CF_D = ''', a, '''' /
     $           '      CALL READ_OPAL_DUMP( IU, CF_D )' )
c
c  Read in the opacities from the binary opacity dump-file:
c
         call read_opal_dump( iu, cf_d )
c
         write(6,1790)
 1790    format( / 'We are done reading the',
     $        ' binary opacity dump-file.' / )
c
      endif
c
c  We are now done reading in the opacities.
c
c
c  The following fairly-obscure subroutines just allow this program to report
c  back to you on some of the characteristics of the opacities now available.
c
      call ask_opal_file_used( 1, cf_used )
c
      write(6,2000) cf_used(:max(1,lnblnk(cf_used)))
 2000 format( / '[The stored OPAL opacities were based on the file ''',
     $     a, ''',' / '  supplemented by other files',
     $     ' as needed/specified by the user input.]' )
c
      call ask_opal_z_mix( 9, xiz, 19, fninz, 19 )
c
      xiz_heavy = 1.0 - xiz(1) - xiz(2) - xiz(3) - xiz(4)
      fninz_heavy = 1.0 - fninz(1) - fninz(2) - fninz(3) - fninz(4)
c
      write(6,2010) 'OPAL', ( xiz(i), fninz(i), i = 1, 19 ),
     $     xiz_heavy, fninz_heavy
 2010 format( / 'The ', a, ' mixture has the following RELATIVE',
     $     ' composition for the metallicity Z:' / /
     $     '         mass      number          mass      number ',
     $     '         mass      number' /
     $     '  el   fraction   fraction  el   fraction   fraction',
     $     '  el   fraction   fraction' /
     $     '  --  ---------  ---------  --  ---------  ---------',
     $     '  --  ---------  ---------' /
     $     '   C', 2f11.7,'   N', 2f11.7,'   O', 2f11.7 /
     $     '  Ne', 2f11.7,'  Na', 2f11.7,'  Mg', 2f11.7 /
     $     '  Al', 2f11.7,'  Si', 2f11.7,'   P', 2f11.7 /
     $     '   S', 2f11.7,'  Cl', 2f11.7,'  Ar', 2f11.7 /
     $     '   K', 2f11.7,'  Ca', 2f11.7,'  Ti', 2f11.7 /
     $     '  Cr', 2f11.7,'  Mn', 2f11.7,'  Fe', 2f11.7 /
     $     '  Ni', 2f11.7,'     total_heavies(above_Ne)=', 2f11.7 )
c
      call ask_z_use( nzuse, zlo, zmid, zhi, zloex, zhiex )
c
      write(6,2020) zlo, zhi
 2020 format( / 'OPAL opacities are available for',
     $     ' metallicities Z from', f10.7, ' to', f10.7 )
c
      call ask_alex_use( kalex, kalex_avail, itype )
c
      if ( kalex_avail .le. 0 ) then
         i_mol = 0
         write(6,2100)
 2100    format( / 'Molecular opacities are NOT available.' )
      else if ( kalex_avail .eq. 1 ) then
         write(6,2110)
 2110    format( / 'The OLD molecular opacities are available,',
     $        ' from Alexander & Ferguson (1994).' )
      else
         write(6,2120)
 2120    format( / 'Molecular opacities of Ferguson et al. (2005)',
     $        ' are available (for ', $ )
         if ( itype .eq. 2 ) then
            write(6,2125) 'GN93'
 2125       format(a,' mix).')
         else if ( itype .eq. 3 ) then
            write(6,2125) 'GS98'
         else if ( itype .eq. 4 ) then
            write(6,2125) 'L03'
         else if ( itype .eq. 5 ) then
            write(6,2125) 'AGS04'
         else if ( itype .eq. 6 ) then
            write(6,2125) 'S92'
         else if ( itype .eq. 7 ) then
            write(6,2125) 'S92AE'
         else if ( itype .ne. 99 .or. cbeg_ferg .eq. ' ' ) then
            write(6,2125) 'unknown'
         else
            write(6,2125) cbeg_ferg(:max(1,lnblnk(cbeg_ferg)))
         endif
      endif
      if ( kalex_avail .gt. 0 ) then
         if ( kalex .gt. 0 ) then
            write(6,2130)
 2130       format( '    These molecular opacities will be',
     $           ' used automatically at low temperatures.' )
         else
            write(6,2140)
 2140       format( '    However, by default these molecular',
     $           ' opacities will be IGNORED (never used).' )
         endif
      endif
c
      call ask_cond_use( kcond, kcond_avail, kreplace_itoh )
c
      if ( kcond_avail .le. 0 ) then
         i_cond = 0
         write(6,2200)
 2200    format( / 'Conductive opacities are NOT available.' )
      else if ( kcond_avail .eq. 1 ) then
         write(6,2210)
 2210    format( / 'The OLD conductive opacities are available,',
     $        ' from Hubbard & Lampe (1969)' )
         if ( kreplace_itoh .ge. 0 ) then
            write(6,2220)
 2220       format( '  [supplemented by the formulae',
     $           ' of Itoh et al. (1983,1984)].' )
         endif
      else
         write(6,2230)
 2230    format( / 'Conductive opacities of Potekhin et al. (2006)',
     $        ' are available [using', $ )
         if ( kreplace_itoh .gt. 1 ) then
            write(6,2240)
 2240       format( / '  biquadratic interpolation in',
     $           ' logT, logRHO, and RMS nuclear charge].' )
         else if ( kreplace_itoh .eq. 1 ) then
            write(6,2250)
 2250       format( ' cubic' / '  interpolation in logT and logRHO,',
     $           ' and linear in the RMS nuclear charge].' )
         else if ( kreplace_itoh .eq. 0 ) then
            write(6,2260)
 2260       format( / '  biquadratic interpolation in',
     $           ' logT, logRHO, and mean nuclear charge].' )
         else if ( kreplace_itoh .lt. 0 ) then
            write(6,2270)
 2270       format( ' cubic' / '  interpolation in logT and logRHO,',
     $           ' and linear in the mean nuclear charge].' )
         endif
      endif
      if ( kcond_avail .gt. 0 ) then
         if ( kcond .gt. 0 ) then
            write(6,2280)
 2280       format( '    These conductive opacities will be',
     $           ' used automatically at high densities.' )
         else
            write(6,2290)
 2290       format( '    However, by default these conductive',
     $           ' opacities will be IGNORED (never used).' )
         endif
      endif
c
      call ask_cno_interp( kcno, kuser, kcno_avail, kuser_avail )
c
      if ( kcno_avail .le. 0 ) then
         i_cno = 0
         write(6,2400)
 2400    format( / 'CNO-interpolation (in the make-up of Z)',
     $        ' is NOT available.' )
      else
         write(6,2410)
 2410    format( / 'CNO-interpolation (in the make-up of Z)',
     $        ' is available if desired.' )
         if ( kcno .gt. 0 ) then
            write(6,2420)
 2420       format( '    The opacity effects from CNO-interpolation',
     $           ' will be used above about 10^4 K.' )
         else
            write(6,2430)
 2430       format( '    However, any opacity effects from',
     $           ' CNO-interpolation will be IGNORED.' )
         endif
      endif
c
c  Now we get to the point where actual computation of opacities takes place.
c
      write(6,5100)
      goto 2810
 2800 write(6,9999)
 2810 write(6,2820)
 2820 format(/ 'Do you wish to compute and output a',
     $     ' (text) opacity table? (y/n): ', $ )
      read(5,'(a1)',err=2800,end=2800) cyn
      if ( cyn .ne. 'y' .and. cyn .ne. 'Y' .and.
     $     cyn .ne. 'n' .and. cyn .ne. 'N' ) goto 2800
c
c  Loop as long as opacity calculations are desired
c
      do while ( cyn .eq. 'y' .or. cyn .eq. 'Y' )
c
         write(6,3000)
 3000    format( / / 'Specify filename and composition for which',
     $        ' to compute the opacity table...' )
c
         goto 4020
 4010    write(6,9999)
 4020    write(6,4030)
 4030    format( / '  Give the name of the output file in which',
     $        ' to store the opacity table:' )
         read(5,'(a255)',err=4010,end=4010) cf_d
c
         goto 4110
 4100    write(6,9999)
 4110    write(6,4120)
 4120    format( / '  How many header lines do you wish',
     $        ' before table values (e.g., 0 or 7): ', $ )
         read(5,*,err=4100,end=4100) n_headers
         if ( n_headers .lt. 0 .or. n_headers .gt. 300 ) then
            write(6,9997) 'Bad value of N_headers'
            goto 4110
         endif
c
         k_mol = 0
         kt_cond = 0
         kr_cond = 0
c
         if ( i_mol .gt. 0 ) then
c
            goto 4050
 4040       write(6,9999)
 4050       write(6,4060) fltopal(1)
 4060       format( / '  Do you wish to extend the table below'
     $           ' lowest OPAL logT =',f5.2,' ? (y/n): ', $ )
            read(5,'(a1)',err=4040,end=4040) cyn
            if ( cyn .ne. 'y' .and. cyn .ne. 'Y' .and.
     $           cyn .ne. 'n' .and. cyn .ne. 'N' ) goto 4040
c
            if ( cyn .ne. 'n' .and. cyn .ne. 'N' ) then
               goto 4080
 4070          write(6,9999)
 4080          write(6,4090)
 4090          format( / '    Give the lowest temperature logT',
     $              ' you wish to have in the table: ', $ )
               read(5,*,err=4070,end=4070) tlog
               if ( tlog .lt. 0.0 .or. tlog .gt. 9.999 ) then
                  write(6,9997) 'Bad value of logT'
                  goto 4080
               endif
               if ( tlog .lt. 3.74999 ) k_mol = 1
               do while ( k_mol .lt. ntferg .and.
     $              tlog .lt. fltferg(max(1,k_mol)) - 1.e-5 )
                  k_mol = k_mol + 1
               enddo
               if ( k_mol .gt. 0 ) then
                  write(6,4095) 'reset to', fltferg(k_mol)
 4095             format( '  ...Table lower limit is ', a,
     $                 ' logT =', f5.2 )
               else
                  write(6,4095) 'still', fltopal(1)
               endif
            endif
c
         endif
c
         if ( i_cond .gt. 0 ) then
c
            goto 4310
 4300       write(6,9999)
 4310       write(6,4320) fltopal(ntopal), 0.5 * ( nropal - 17 )
 4320       format( / '  Do you wish to extend the table',
     $           ' (using conductive opacities)' /
     $           '    above the highest OPAL logT =',f5.2,
     $           ' and/or logR =', f4.1, ' ? (y/n): ', $ )
            read(5,'(a1)',err=4300,end=4300) cyn
            if ( cyn .ne. 'y' .and. cyn .ne. 'Y' .and.
     $           cyn .ne. 'n' .and. cyn .ne. 'N' ) goto 4300
c
            if ( cyn .ne. 'n' .and. cyn .ne. 'N' ) then
               goto 4340
 4330          write(6,9999)
 4340          write(6,4350)
 4350          format( / '    Give the maximum logT and logR',
     $              ' values you want: ', $ )
               read(5,*,err=4330,end=4330) tlog, rlog
               if ( tlog .lt. 0.0 .or. tlog .gt. 99.999 ) then
                  write(6,9997) 'Bad value of logT'
                  goto 4330
               else if ( rlog .lt. -9.0 .or. rlog .gt. 99.999 ) then
                  write(6,9997) 'Bad value of logR'
                  goto 4330
               endif
               tlog = min( tlog , 9.99 ) - 1.e-5
               do while ( fltopal(ntopal) + 0.1 * kt_cond .lt. tlog )
                  kt_cond = kt_cond + 1
               enddo
               rlog = min( rlog , 16.0 ) - 1.e-5
               do while ( 0.5 * ( nropal + kr_cond - 17 ) .lt. rlog )
                  kr_cond = kr_cond + 1
               enddo
               if ( kt_cond .gt. 0 .or. kr_cond .gt. 0 ) then
                  write(6,4360) 'reset to',
     $                 min( 9.99 , fltopal(ntopal) + 0.1 * kt_cond ),
     $                 0.5 * ( nropal + kr_cond - 17 )
 4360             format( '  ...Table upper limts are ', a,
     $                 ' logT =', f5.2, ' , logR =', f5.1 )
               else
                  write(6,4360) 'still', fltopal(ntopal),
     $                 0.5 * ( nropal - 17 )
               endif
            endif
c
         endif
c
         if ( i_cno .le. 0 ) then
            cyn = 'n'
         else
            goto 3020
 3010       write(6,9999)
 3020       write(6,3030)
 3030       format( / '  Do you wish to explicitly specify',
     $           ' the component abundances in Z? (y/n): ', $ )
            read(5,'(a1)',err=3010,end=3010) cyn
            if ( cyn .ne. 'y' .and. cyn .ne. 'Y' .and.
     $           cyn .ne. 'n' .and. cyn .ne. 'N' ) goto 3010
         endif
c
         if ( cyn .ne. 'y' .and. cyn .ne. 'Y' ) then
c
            k_cno = 0
c
            goto 3050
 3040       write(6,9999)
 3050       write(6,3060)
 3060       format( / '  Give  X      (the hydrogen mass fraction),' /
     $           '        Z      (the mass fraction of metals),' /
     $           '        exC    (any excess carbon abundance',
     $           ' [beyond that contained in Z]),' /
     $           '        exO    (any excess oxygen abundance',
     $           ' [beyond that contained in Z]):' )
            read(5,*,err=3040,end=3040) x, z, exc, exo
c
            y = 1.0 - x - z - exc - exo
c
            if ( x .lt. 0.0 .or. x .gt. 1.0 ) then
               write(6,9997) 'Bad value of X'
               goto 3050
            else if ( z .lt. 0.0 .or. z .gt. 1.0 ) then
               write(6,9997) 'Bad value of Z'
               goto 3050
            else if ( exc .lt. -z .or. exc .gt. 1.0 ) then
               write(6,9997) 'Bad value of exC'
               goto 3050
            else if ( exo .lt. -z .or. exc .gt. 1.0 ) then
               write(6,9997) 'Bad value of exO'
               goto 3050
            else if ( exc + exo .lt. -z ) then
               write(6,9997) 'Bad values: exC + exO < -Z'
               goto 3050
            else if ( y .lt. -1.e-6 ) then
               write(6,9997) 'Bad values: X + Z + exC + exO > 1'
               goto 3050
            endif
c
            call opal( z, x, exc, exo, 0.0, -4.0 )
c
            call ask_last_opac( opac, dkt, dkr, dktrho,
     $                          fedge, ftr, fz )
c
            call ask_last_xcnou( z_use, x_use, exc_use, exo_use,
     $                           slt_use, slr_use,
     $                           fcn, fcon, fcnone, fu )
c
            y_use = 1.0 - x_use - z_use - exc_use - exo_use
c
            write(6,3097) y_use, z_use, x_use, exc_use, exo_use
c
            if ( fz .le. 1.e-6 ) then
               write(6,3069) z_use
 3069          format( '...Error: Z =', f10.7, ' gives FZ = 0.0',
     $              ' (out of Z-range): table NOT computed.' )
               goto 3910
            endif
c
         else
c
            k_cno = 1
c
            goto 3080
 3070       write(6,9999)
 3080       write(6,3090)
 3090       format( / '  Give  X      (the hydrogen mass fraction),' /
     $           '        xC     (the [total] carbon mass fraction),' /
     $           '        xN     (the nitrogen mass fraction),' /
     $           '        xO     (the [total] oxygen mass fraction),' /
     $           '        xNe    (the neon mass fraction),' /
     $           '        xHeavy (the total mass fraction',
     $           ' of elements heavier than neon):' )
            read(5,*,err=3070,end=3070) x, ( xmet(i), i = 1, 5 )
c
            y = 1.0 - x - xmet(1) - xmet(2) - xmet(3)
     $           - xmet(4) - xmet(5)
c
            if ( x .lt. 0.0 .or. x .gt. 1.0 ) then
               write(6,9997) 'Bad value of X'
               goto 3080
            else if ( xmet(1) .lt. 0.0 .or. xmet(1) .gt. 1.0 ) then
               write(6,9997) 'Bad value of xC'
               goto 3080
            else if ( xmet(2) .lt. 0.0 .or. xmet(2) .gt. 1.0 ) then
               write(6,9997) 'Bad value of xN'
               goto 3080
            else if ( xmet(3) .lt. 0.0 .or. xmet(3) .gt. 1.0 ) then
               write(6,9997) 'Bad value of xO'
               goto 3080
            else if ( xmet(4) .lt. 0.0 .or. xmet(4) .gt. 1.0 ) then
               write(6,9997) 'Bad value of xNe'
               goto 3080
            else if ( xmet(5) .lt. 0.0 .or. xmet(5) .gt. 1.0 ) then
               write(6,9997) 'Bad value of xHeavy'
            else if ( y .lt. -1.e-6 ) then
               write(6,9997) 'Bad values: sum{Xi} > 1'
               goto 3080
            endif
c
            call opal_x_cno_fu( x, 0.0, -4.0, xmet, 5, 0.0 )
c
            call ask_last_opac( opac, dkt, dkr, dktrho,
     $                          fedge, ftr, fz )
c
            call ask_last_xcnou( z_use, x_use, exc_use, exo_use,
     $                           slt_use, slr_use,
     $                           fcn, fcon, fcnone, fu )
c
            y_use = 1.0 - x_use - z_use - exc_use - exo_use
c
            write(6,3097) y_use, z_use, x_use, exc_use, exo_use
 3097       format( / 'The above yields Y =', f11.7,
     $           ' from the composition:' /
     $           '  Z =', f11.7, ',  X =', f11.7,
     $           ',  exC =', f11.7, ',  exO =', f11.7 )
c
            if ( fz .le. 1.e-6 ) then
               write(6,3098)
 3098          format( '...Error: this gives FZ = 0.0',
     $              ' (out of allowed Z-range): table NOT computed.' )
               goto 3910
            endif
c
            write(6,3099) fcn, fcon, fcnone
 3099       format( 'and the following CNO-interpolation factors:' /
     $           '  f(C->N) =', f11.6, ',  f(CO->N) =', f11.6,
     $           ',  f(CNO->Ne) =', f11.6 )
c
         endif
c
         goto 4210
 4200    write(6,9999)
 4210    write(6,4220) cf_d(:max(1,lnblnk(cf_d))), n_headers
 4220    format( / 'Create new opacity output file called' /
     $        '  ', a / 'with', i4, ' header lines', $ )
         if ( k_mol .gt. 0 ) write(6,4230) fltferg(k_mol)
 4230    format( ' (extend table down to logT =', f5.2, ')', $ )
         if ( kt_cond .gt. 0 .or. kr_cond .gt. 0 ) write(6,4235)
     $        min( 9.99 , fltopal(ntopal) + 0.1 * kt_cond ),
     $        0.5 * ( nropal + kr_cond - 17 )
 4235    format( / '(extend table up to logT =', f5.2,
     $        ' logR =', f5.1, ')', $ )
         write(6,4240)
 4240    format( / '  Confirm? (y/n): ', $ )
         read(5,'(a1)',err=4200,end=4200) cyn
         if ( cyn .ne. 'y' .and. cyn .ne. 'Y' .and.
     $        cyn .ne. 'n' .and. cyn .ne. 'N' ) goto 4200
c
         if ( cyn .ne. 'y' .and. cyn .ne. 'Y' ) goto 3910
c
c  Open output file and store any header lines
c
         iu = 33
c
         call inqfil( cf_d, lxst )
c
         if ( lxst ) then
            goto 4980
 4970       write(6,9999)
 4980       write(6,4990) cf_d(:lnblnk(cf_d))
 4990       format( / 'This file already exists:' / '  ', a /
     $           'Do you wish to overwrite it? (y/n): ', $ )
            read(5,'(a1)',err=4970,end=4970) cyn
            if ( cyn .ne. 'y' .and. cyn .ne. 'Y' .and.
     $           cyn .ne. 'n' .and. cyn .ne. 'N' ) goto 4970
            if ( cyn .ne. 'y' .and. cyn .ne. 'Y' ) then
 4995          format('Table NOT created.')
               goto 3910
            endif
         endif
c
         call linux_get_home_dir(cf_d,fnalt,ialt)
c
         if ( ialt .gt. 0 ) then
            open(iu,file=fnalt,form='FORMATTED',status='UNKNOWN',
     $           iostat=ioperr,err=4900)
         else
            open(iu,file=cf_d,form='FORMATTED',status='UNKNOWN',
     $           iostat=ioperr,err=4900)
         endif
c
         goto 5000
c
 4900    write(6,4910) ioperr, iu, cf_d(:lnblnk(cf_d))
 4910    format(/'Error',i12,' opening unit',i3,' for new file:'/
     $        '  ', a / 'Table NOT created.' )
         goto 3910
c
 5000    continue
c
         if ( n_headers .ge. 20 ) then
            write(iu,2010) cf_used(:max(1,lnblnk(cf_used))),
     $           ( xiz(i), fninz(i), i = 1, 19 ),
     $           xiz_heavy, fninz_heavy
            n_headers = n_headers - 13
         endif
c
         if ( n_headers .ge. 9 ) then
            write(iu,5010) z_use, x_use, exc_use, exo_use, y_use
 5010       format( / '  Z =', f11.7, ',  X =', f11.7,
     $           ',  exC =', f11.7, ',  exO =', f11.7,
     $           '  [Y =', f11.7, ']' )
            n_headers = n_headers - 2
            if ( k_cno .gt. 0 .and. n_headers .ge. 8 ) then
               write(iu,5020) fcn, fcon, fcnone
 5020          format( 'with  f(C->N) =', f11.6, ',  f(CO->N) =',
     $              f11.6, ',  f(CNO->Ne) =', f11.6 )
               n_headers = n_headers - 1
            endif
         endif
c
         if ( k_mol .gt. 0 .and. n_headers .ge. 9 ) then
            write(iu,5022) fltferg(k_mol)
 5022       format( / '  ...Table extended via molecular opacities,',
     $           ' down to logT =', f5.2 )
            n_headers = n_headers - 2
         endif
c
         if ( ( kt_cond .gt. 0 .or. kr_cond .gt. 0 ) .and.
     $        n_headers .ge. 9 ) then
            write(iu,5024)
     $           min( 9.99 , fltopal(ntopal) + 0.1 * kt_cond ),
     $           0.5 * ( nropal + kr_cond - 17 )
 5024       format( / '  ...Table extended via conductive opacities,',
     $           ' up to logT =', f5.2, ' , logR =', f5.1 )
            n_headers = n_headers - 2
         endif
c
         if ( n_headers .ge. 7 ) then
            do i = 7, n_headers
c               write(iu,'("")')
               write(iu,5100)
 5100          format( ' ' )
            enddo
         endif
c
         if ( n_headers .ge. 6 )
     $        write(iu,5030) cf_used(:max(1,lnblnk(cf_used))),
     $        x_use, y_use, z_use, exc_use, exo_use
 5030    format( 'TABLE #  1     ', a16, '   X=', f6.4, ' Y=', f6.4,
     $        ' Z=', f6.4, ' dXc=', f6.4, ' dXo=', f6.4 )
c
         if ( n_headers .ge. 5 ) write(iu,5100)
c
         if ( n_headers .ge. 4 ) write(iu,5040)
 5040    format( '                                                  ',
     $        'log R' )
c
         if ( n_headers .ge. 3 ) write(iu,5100)
c
         if ( n_headers .ge. 2 ) write(iu,5050)
     $        ( 0.5 * ( i - 17 ), i = 1, nropal + kr_cond )
 5050    format( 'logT', f6.1, 99f7.1 )
c
         if ( n_headers .ge. 1 ) write(iu,5100)
c
c  Loop over T to create lines in table...
c
         do k = 1 - k_mol, ntopal + kt_cond
c
            if ( k .le. 0 ) then
               tlog = fltferg(1-k)
            else if ( k .le. ntopal ) then
               tlog = fltopal(k)
            else
               tlog = min( 9.99 ,
     $              ( k - ntopal ) * 0.1 + fltopal(ntopal) )
            endif
c
            write(iu,'(f4.2,$)') tlog
c
            slt = tlog - 6.0
c
c  Loop over RHO to fill this line of the table...
c
            j = 0
c
            do while ( j .lt. nropal + kr_cond .and.
     $           ( j .lt. nropal .or. k .gt. 0 ) )
c
               j = j + 1
c
               slr = 0.5 * ( j - 17 )
c
               if ( k_cno .eq. 0 ) then
c
c  Call the simpler OPAL-opacity interface, if no CNO-interpolation is needed:
c
                  call opal( z, x, exc, exo, slt, slr )
c
               else
c
c  otherwise, call the more complicated OPAL-opacity interface:
c
                  call opal_x_cno_fu( x, slt, slr, xmet, 5, 0.0 )
c
               endif
c
               call ask_last_opac( opac, dkt, dkr, dktrho,
     $                             fedge, ftr, fz )
c
               if ( ftr .ge. 0.9 .and. opac .gt. -29.999 .and.
     $              ( ( k .le. ntopal  .and. opac .lt. 9.999 ) .or.
     $              ( k .gt. ntopal  .and. opac .lt. 19.999 ) ) ) then
c
                  if ( j .lt. nropal + kr_cond .and.
     $                 ( j .lt. nropal .or. k .gt. 0 ) ) then
                     write(iu,5110) opac
 5110                format(f7.3,$)
                  else
                     write(iu,5120) opac
 5120                format(f7.3)
                  endif
c
               else if ( ( j .le. nropal .or. k .le. 0 ) .and.
     $                 ( tlog .lt. 7.0 .or.
     $                 ( j .le. 14 .and. k .le. ntopal ) ) ) then
c
                  if ( j .lt. nropal + kr_cond .and.
     $                 ( j .lt. nropal .or. k .gt. 0 ) ) then
                     write(iu,5130)
 5130                format('  9.999',$)
                  else
                     write(iu,5140)
 5140                format('  9.999')
                  endif
c
               else if ( k .gt. ntopal .or. j .gt. nropal ) then
c
                  if ( j .lt. nropal + kr_cond ) then
                     write(iu,5150)
 5150                format(' 99.999',$)
                  else
                     write(iu,5160)
 5160                format(' 99.999')
                  endif
c
               else
c
                  write(iu,5100)
                  j = nropal + kr_cond
c
               endif
c
            enddo
c
c  End of table line
c
         enddo
c
         close(iu)
c
         goto 3910
 3900    write(6,9999)
 3910    write(6,3920)
 3920    format( / 'Do you wish to compute another',
     $        ' (text) opacity table? (y/n): ', $ )
         read(5,'(a1)',err=3900,end=3900) cyn
         if ( cyn .ne. 'y' .and. cyn .ne. 'Y' .and.
     $        cyn .ne. 'n' .and. cyn .ne. 'N' ) goto 3900
c
      enddo
c
      stop
      end
