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