c
c  This program provides an interactive example of use of the OPAL opacity
c  calculation subroutines.  Compile by, for example,
c
c    f77 -o z_xcotrin_example z_xcotrin_example.f z14xcotrin21.f
c
c  and then run it via
c
c    ./z_xcotrin_example
c
c  It will ask you for the input values needed to read and calculate opacities.
c
c  (Note that internal real variables are single precision.)
c
      program z_xcotrin_example
c     -------------------------
c
      character*255 cf_hz, cf_ofe, cbeg_ferg, cf_d
      character*255 cdirin, cdir_mol, cdir_cond, cf_used
      character*1 cyn
c
      dimension xiz(19), fninz(19), xmet(5)
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( / 'Example of use of OPAL opacity routines in',
     $     ' z14xcotrin21.f' / /
     $     '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,'("------ 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,'("------ Input error ... try again:",$)')
 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,'("=== Bad Z_READ value ... try again:",$)')
               goto 190
            endif
         endif
c
         goto 220
 210     write(6,'("------ Input error ... try again:",$)')
 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,'("=== No filename was input ... try again:",$)')
            goto 220
         endif
c
         goto 250
 240     write(6,'("------ Input error ... try again:",$)')
 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,'("=== Bad OFE_BRACKET value ... try again:",$)')
            goto 250
         else if ( ofe_bracket .ne. 0.0 ) then
            goto 280
 270        write(6,'("------ Input error ... try again:",$)')
 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,'("------ Input error ... try again:",$)')
 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,'("=== Bad I_MOL value ... try again:",$)')
            goto 90
         else if ( i_mol .eq. 991 ) then
            goto 120
 110        write(6,'("------ 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,'("=== No name was input ... try again:",$)')
               goto 120
            endif
         endif
c
         goto 60
 50      write(6,'("------ Input error ... try again:",$)')
 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,'("=== Bad I_COND value ... try again:",$)')
            goto 60
         endif
c
         goto 150
 140     write(6,'("------ Input error ... try again:",$)')
 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,'("=== Bad I_CNO value ... try again:",$)')
            goto 150
         endif
c
         goto 310
 300     write(6,'("------ Input error ... try again:",$)')
 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,'("------ Input error ... try again:",$)')
 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,'("------ Input error ... try again:")')
 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,'("------ Input error ... try again:",$)')
 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,'("------ Input error ... try again:")')
 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,'("------ Input error ... try again:",$)')
 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
            write(6,1010) list_level
 1010       format( / '      LIST_LEVEL =', i12 /
     $           '      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
         write(6,1020) cdirin(:max(1,lnblnk(cdirin)))
 1020    format( / '      CDIRIN = ''', a, '''' /
     $        '      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
            write(6,1030) cdir_mol(:max(1,lnblnk(cdir_mol)))
 1030       format( / '      CDIR_MOL = ''', a, '''' /
     $           '      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
            write(6,1050) cbeg_ferg(:max(1,lnblnk(cbeg_ferg)))
 1050       format( / '      CBEG_FERG = ''', a, '''' /
     $           '     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
            write(6,1060) iu, z_read, cf_hz(:max(1,lnblnk(cf_hz))),
     $           ofe_bracket, cf_ofe(:max(1,lnblnk(cf_ofe)))
 1060       format( / '      IU =', i3 /
     $           '      Z_READ =', f10.7 /
     $           '      CF_HZ = ''', a, '''' /
     $           '      OFE_BRACKET =', f10.6 /
     $           '      CF_OFE = ''', a, '''' /
     $           '      CALL READ_BASIC_OPAL_OPAC( IU,',
     $           ' 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
            write(6,1070) iu, z_read, cf_hz(:max(1,lnblnk(cf_hz))),
     $           ofe_bracket, cf_ofe(:max(1,lnblnk(cf_ofe))),
     $           i_mol, i_cond, i_cno
 1070       format( / '      IU =', i3 /
     $           '      Z_READ =', f10.7 /
     $           '      CF_HZ = ''', a, '''' /
     $           '      OFE_BRACKET =', f10.6 /
     $           '      CF_OFE = ''', a, '''' /
     $           '      I_MOL =',i4 /
     $           '      I_COND =',i3 /
     $           '      I_CNO =',i4 /
     $           '      CALL READ_EXTENDED_OPAC( IU,',
     $           ' Z_READ, CF_HZ, OFE_BRACKET, CF_OFE,' /
     $           '     $                         I_MOL,',
     $           ' 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,'("------ Input error ... try again:",$)')
 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,'("------ Input error ... try again:",$)')
 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,
     $              '("=== No filename was input ... try again:",$)')
               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,'("------ Input error ... try again:",$)')
 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,'("=== No filename was input ... try again:",$)')
            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) ( xiz(i), fninz(i), i = 1, 19 ),
     $     xiz_heavy, fninz_heavy
 2010 format( / 'The OPAL 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,'("GN93 mix).")')
         else if ( itype .eq. 3 ) then
            write(6,'("GS98 mix).")')
         else if ( itype .eq. 4 ) then
            write(6,'("L03 mix).")')
         else if ( itype .eq. 5 ) then
            write(6,'("AGS04 mix).")')
         else if ( itype .eq. 6 ) then
            write(6,'("S92 mix).")')
         else if ( itype .eq. 7 ) then
            write(6,'("S92AE mix).")')
         else if ( itype .ne. 99 .or. cbeg_ferg .eq. ' ' ) then
            write(6,'("unknown mix).")')
         else
            write(6,'(a," mix).")')
     $           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,'(" ")')
      goto 2810
 2800 write(6,'("------ Input error ... try again:",$)')
 2810 write(6,2820)
 2820 format(/ 'Do you wish to compute opacity values? (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 temperature, density, and',
     $        ' composition for which to compute opacity:' )
c
         if ( i_cno .le. 0 ) then
            cyn = 'n'
         else
            goto 3020
 3010       write(6,'("------ Input error ... try again:",$)')
 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
            goto 3050
 3040       write(6,'("------ Input error ... try again:",$)')
 3050       write(6,3060)
 3060       format( / '  Give  Tlog   (log10{T}, the log',
     $           ' of the temperature in K),' /
     $           '        RHOlog (log10{RHO}, the log',
     $           ' of the density in g/cm^3),' /
     $           '        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) tlog, rholog, x, z, exc, exo
c
            y = 1.0 - x - z - exc - exo
c
            if ( tlog .lt. 2.0 .or. tlog .gt. 11.0 ) then
               write(6,'("=== Bad value of logT ... try again:",$)')
               goto 3050
            else if ( rholog .lt. -20.0 .or. rholog .gt. 15.0 ) then
               write(6,'("=== Bad value of logRHO ... try again:",$)')
               goto 3050
            else if ( x .lt. 0.0 .or. x .gt. 1.0 ) then
               write(6,'("=== Bad value of X ... try again:",$)')
               goto 3050
            else if ( z .lt. 0.0 .or. z .gt. 1.0 ) then
               write(6,'("=== Bad value of Z ... try again:",$)')
               goto 3050
            else if ( exc .lt. -z .or. exc .gt. 1.0 ) then
               write(6,'("=== Bad value of exC ... try again:",$)')
               goto 3050
            else if ( exo .lt. -z .or. exc .gt. 1.0 ) then
               write(6,'("=== Bad value of exO ... try again:",$)')
               goto 3050
            else if ( exc + exo .lt. -z ) then
               write(6,3067)
 3067          format( '=== Bad values: exC + exO < -Z',
     $              ' ... try again:', $ )
               goto 3050
            else if ( y .lt. -1.e-6 ) then
               write(6,3068)
 3068          format( '=== Bad values: X + Z + exC + exO > 1',
     $              ' ... try again:', $ )
               goto 3050
            endif
c
         else
c
            goto 3080
 3070       write(6,'("------ Input error ... try again:",$)')
 3080       write(6,3090)
 3090       format( / '  Give  Tlog   (log10{T}, the log',
     $           ' of the temperature in K),' /
     $           '        RHOlog (log10{RHO}, the log',
     $           ' of the density in g/cm^3),' /
     $           '        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) tlog, rholog, x,
     $           ( xmet(i), i = 1, 5 )
c
            y = 1.0 - x - xmet(1) - xmet(2) - xmet(3)
     $           - xmet(4) - xmet(5)
c
            if ( tlog .lt. 2.0 .or. tlog .gt. 11.0 ) then
               write(6,'("=== Bad value of logT ... try again:",$)')
               goto 3080
            else if ( rholog .lt. -20.0 .or. rholog .gt. 15.0 ) then
               write(6,'("=== Bad value of logRHO ... try again:",$)')
               goto 3080
            else if ( x .lt. 0.0 .or. x .gt. 1.0 ) then
               write(6,'("=== Bad value of X ... try again:",$)')
               goto 3080
            else if ( xmet(1) .lt. 0.0 .or. xmet(1) .gt. 1.0 ) then
               write(6,'("=== Bad value of xC ... try again:",$)')
               goto 3080
            else if ( xmet(2) .lt. 0.0 .or. xmet(2) .gt. 1.0 ) then
               write(6,'("=== Bad value of xN ... try again:",$)')
               goto 3080
            else if ( xmet(3) .lt. 0.0 .or. xmet(3) .gt. 1.0 ) then
               write(6,'("=== Bad value of xO ... try again:",$)')
               goto 3080
            else if ( xmet(4) .lt. 0.0 .or. xmet(4) .gt. 1.0 ) then
               write(6,'("=== Bad value of xNe ... try again:",$)')
               goto 3080
            else if ( xmet(5) .lt. 0.0 .or. xmet(5) .gt. 1.0 ) then
               write(6,'("=== Bad value of xHeavy ... try again:",$)')
            else if ( y .lt. -1.e-6 ) then
               write(6,
     $              '("=== Bad values: sum{Xi} > 1 ... try again:",$)')
               goto 3080
            endif
c
         endif
c
         slt = tlog - 6.0
         slr = rholog - 3.0 * slt
c
         write(6,3200) y
 3200    format( / '[this implies a helium mass fraction of  Y =',
     $        f10.6, '; also, recall that' / ' SLT = Tlog - 6.0  and'
     $        '  SLR = RHOlog - 3.0 * SLT].  Now: compute opacity:' )
c
         if ( cyn .ne. 'y' .and. cyn .ne. 'Y' ) then
c
            write(6,3300) z, x, exc, exo, slt, slr
 3300       format( / '      Z   =', f11.7 /
     $           '      X   =', f11.7 /
     $           '      exC =', f11.7 /
     $           '      exO =', f11.7 /
     $           '      SLT =', f11.6 /
     $           '      SLR =', f11.6 /
     $           '      CALL OPAL( Z, X, exC, exO, SLT, SLR )' )
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
            write(6,3400) x, slt, slr, ( xmet(i), i = 1, 5 )
 3400       format( /
     $           '      X       =', f11.7 /
     $           '      SLT     =', f11.6 /
     $           '      SLR     =', f11.6 /
     $           '      XMET(1) =', f11.7 /
     $           '      XMET(2) =', f11.7 /
     $           '      XMET(3) =', f11.7 /
     $           '      XMET(4) =', f11.7 /
     $           '      XMET(5) =', f11.7 /
     $           '      NMET    = 5' /
     $           '      FU      = 0.0' /
     $           '      CALL OPAL_X_CNO_FU( X,',
     $           ' SLT, SLR, XMET, NMET, FU )' )
c
            call opal_x_cno_fu( x, slt, slr, xmet, 5, 0.0 )
c
         endif
c
         write(6,3450)
 3450    format( '      CALL ASK_LAST_OPAC( OPAC,',
     $        ' DKT, DKR, DKTRHO, FEDGE, FTR, FZ )' )
c
         call ask_last_opac( opac, dkt, dkr, dktrho,
     $                       fedge, ftr, fz )
c
         write(6,3500)
 3500    format( / 'The above returns the following values',
     $        ' for the edge-tests and the opacity:' )
c
         write(6,3510) '    FZ', fz
 3510    format( '> ', a6, ' =', f11.7, '  ', $ )
c
         if ( fz .le. 0.0 ) then
c
            write(6,3520)
 3520       format( '[far beyond allowed Z-range:',
     $           ' no extrapolation possible]' )
            write(6,3510) '   FTR', ftr
            write(6,3530)
 3530       format( '[this may not have been computed',
     $           ' (because FZ was 0.0)]' )
c
         else
c
            if ( fz .le. 0.999999 ) then
               write(6,3540)
 3540          format( '[just outside allowed Z-range:',
     $              ' was able to extrapolate]' )
            else
               write(6,3550)
 3550          format( '[inside the allowed Z-range:',
     $              ' was able to interpolate]' )
            endif
c
            write(6,3510) '   FTR', ftr
c
            if ( ftr .le. 0.0 ) then
               write(6,3560)
 3560          format( '[far outside allowed (T6,R):',
     $              ' no extrapolation possible]' )
            else if ( ftr .le. 0.999999 ) then
               write(6,3570)
 3570          format( '[just outside allowed (T6,R):',
     $              ' was able to extrapolate]' )
            else
               write(6,3580)
 3580          format( '[inside the allowed (T6,R):',
     $              ' was able to interpolate]' )
            endif
c
         endif
c
         write(6,3510) ' FEDGE', fedge
         if ( fedge .le. 0.0 ) then
            write(6,'("[far outside",$)')
         else if ( fedge .le. 0.999999 ) then
            write(6,'("[just outside",$)')
         else
            write(6,'("[inside",$)')
         endif
         write(6,3590)
 3590    format( ' the (Z,T6,R) boundaries of the matrices]' )
c
         if ( abs( opac ) .gt. 99.9999 ) then
            write(6,3600) '  OPAC', opac
 3600       format( '> ', a6, ' =', 1p, e11.3, '  ', $ )
         else
            write(6,3610) '  OPAC', opac
 3610       format( '> ', a6, ' =', f11.6, '  ', $ )
         endif
c
         if ( opac .gt. 9.99e+34 ) then
c
            write(6,3620)
 3620       format( '[the value of log10(Kappa)',
     $           ' was not computed at all]' )
c
         else
c
            if ( abs( opac ) .gt. 29.9999 ) then
               write(6,3630)
 3630          format( '[a bad value of log10(Kappa):',
     $              ' this should not happen]' )
            else
               write(6,3640)
 3640          format( '[log10(Kappa): the opacity',
     $              ' value that was computed]' )
            endif
c
            if ( abs( dkt ) .gt. 99.9999 ) then
               write(6,3600) '   DKT', dkt
            else
               write(6,3610) '   DKT', dkt
            endif
            write(6,3650)
 3650       format( '[dLog(Kappa)/dLog(T) at constant R',
     $           ' (not constant RHO)]' )
c
            if ( abs( dkr ) .gt. 99.9999 ) then
               write(6,3600) '   DKR', dkr
            else
               write(6,3610) '   DKR', dkr
            endif
            write(6,3660)
 3660       format( '[dLog(Kappa)/dLog(RHO) at constant T]' )
c
            if ( abs( dktrho ) .gt. 99.9999 ) then
               write(6,3600) 'DKTRHO', dktrho
            else
               write(6,3610) 'DKTRHO', dktrho
            endif
            write(6,3670)
 3670       format( '[dLog(Kappa)/dLog(T) at constant RHO]' )
c
         endif
c
c  Get the composition actually used when obtaining the opacity (which may
c  not always be exactly the same as the input composition), and any
c  CNO-interpolation factors that may have been used.
c
         call ask_last_xcnou( z_use, x_use, exc_use, exo_use,
     $                        slt_use, slr_use, fcn, fcon, fcnone, fu )
c
         write(6,3800) z_use, x_use, exc_use, exo_use
 3800    format( 'The above was obtained using the composition:' /
     $        '  Z =', f11.7, ',  X =', f11.7,
     $        ',  exC =', f11.7, ',  exO =', f11.7 )
c
         if ( ( cyn .eq. 'y' .or. cyn .eq. 'Y' ) .and.
     $        opac .lt. 9.99e+34 ) then
            write(6,3810) fcn, fcon, fcnone
 3810       format( 'and the following CNO-interpolation factors:' /
     $           '  f(C->N) =', f11.6, ',  f(CO->N) =', f11.6,
     $           ',  f(CNO->Ne) =', f11.6 )
         endif
c
         goto 3910
 3900    write(6,'("------ Input error ... try again:",$)')
 3910    write(6,3920)
 3920    format( / 'Do you wish to compute another',
     $        ' opacity value? (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
