Hot Fixes to C90.03
These are minor changes that will be incorporated in the next revision of the code. These only apply to the current version of the code stored on this site. In the following, the corrections that should be made to version 90.03 are indicated by line numbers within specific subroutines.
The code could crash when the length scale for a specified column density exceeded 10^38 cm. Two chages are needed. Edit routine nextdr.f and change line 1057 to become the following three statements:
if((colnut-chi-real(drad)*coleff).lt.coleff*drnxt ) then drnxt = (colnut-chi-real(drad)*coleff)*z/coleff endif
Next edit routine frstdr and change the block of code between lines 105 and 114 to read
if( colend.lt.5e29 ) then * this is useful for very thin columns, normally larger than 1st DR drcol = log10(colend ) - log10( hden * filfac * 20. ) else if( colpls.lt.5e29 ) then drcol = log10(colpls ) - log10( hii * filfac * 20. ) else if( colnut.lt.5e29 ) then drcol = log10(colnut ) - log10( hi * filfac * 50. ) else drcol = 30. endif drcol = 10.**min( 35., drcol )
Thanks to Ravi Sankrit for uncovering this problem. 1997 Oct 30.
Energy not conserved when large turbulences were assumed. When the helium line optical depth scale was overrun by a bad initial guess for the total optical depths, the continuum pumping rates for Lya were not reevaluated. This resulted in the pumping rates remaining at the values they had when the optical depth scale was first overrun. The solution is to continue to reevaluate the continuum pumping rates. Edit routine hetrans.f and add the following after line 238:
he1dst(2,1) = (1.-DstWght)*he1dst(2,1) + DstWght*he1as(2,1) * 1 eovrlp( heii*he1n(1) ,he1tau(2,1),he1lim(2,1) , damp , 2 he1opc(2,1), iphe1l(2,1) ,hevel,he1esc(2,1),'INCO' ) else damp = eshe2l( esin ) he1con(2,1) = he1jbr(2,1) * esin *he1as(2,1) endif
(add the else, followed by the evaluation of damp and he1con). After doing that go to the new line number 254 and add a similar reevaluation:
he2dst(2,1) = (1.-DstWght)*he2dst(2,1) + DstWght*he2as(2,1)* 1 eovrlp( heiii*he2n(1) ,he2tau(2,1),he2lim(2,1) , damp , 2 he2opc(2,1), iphe2l(2,1) , hevel , he2esc(2,1) ,'INCO' ) else damp = eshe2l( esin ) he2con(2,1) = he2jbr(2,1) * esin *he2as(2,1) endif
Edit routine hydropesc.f and move line 183 so that the evaluation of the continuum pumping rate follows the endif:
hfrcin(2,1) = min(1., max(0.,fracin ) ) endif * >>chng 97 sep 13, move out of it hcont(2,1) = hjbar(2,1) * esa0k2(htnext(2,1)) * hyas(2,1) *
Thanks to Fred Hamann for finding this problem, and to Rafael Reyes for finding the problem with the correction to the problem. 1997 Sept 13.
more problems on alphas:
Do not compile routines linefit.f, hcolst.f, he1col.f, and he2col.f, using the fast option. Use f77 -fpe0 -c name.f instead. There seem to be problems with the fast version of the Dec math floating point library. This only affects Alphas and these routines can safely be compiled with the fast option on other platforms.
Thanks to Roderick Johnstone for finding this problem. 1997 July 30.
The following list of changes define C90.03b. After making them in your copy, change line 7 of version.f to 90.03b.
CO matrix inversion routine could fail in fully molecular limit, if the OH abundance underflows first. Edit routine codriv.f and change line 41 from
if( oh.eq.0. ) return to if( oh.eq.0. ) then if( co/acarb.gt.colimt ) lotsco = .true. return endif
This had no effect on results but could cause the code to terminate and declare a failed calculation. Thanks to Bhaswati Mookerjea for finding this problem. 1997 July 29.
Trouble converging a 2000 K cloud exposed to hard x-rays. Edit tfirst.f and change line 372 from
looplimit = 8 to looplimit = 12
This had no effects on results but could prevent the code from finding an initial solution. Thanks to Jane Charlton for finding this problem. 1997 July 29
Incorrect type for argument in call to tabden. Edit routine frstdr.f and change line 159 to read
factor = hden/tabden( radius+drTabDen, dble(drTabDen) )
(pass the second argument to tabden as double precision.) Thanks to Karen Vanlandingham for finding this problem. 1997 July 24.
Incorrect variable for three-body recombination of hydrogen. This prevented the atom from going to the exact LTE limit at very high densities. Edit routine hlevel and change line 404 from
1 dble(rinduc(i)))/hlte(i) + dble(hcolnc(i)*eden) to 1 dble(rinduc(i)))/hlte(i) + dble(hcolnc(i)*edhi)
This had no effects on results except at densities greater than 1e14 cm^-3, where the slight departure from LTE it introduced caused three-body heating to be significant. 1997 July 25
Cooling due to 1s-2s collisions in HeI and HeII was not included. This had no effect on clouds with nebular temperatures, but clouds that equilibrated between 1e5K and 1e6K are about 10% cooler when this is included. The correction is to download a revised copy of coolr.f and replace the existing file. 1997 July 24.
Punch lines intensity command produced output on the wrong file when more than one punch command was entered. Edit file dopunch.f and changes lines 864 and 869 so that they read (change npunch to i): 1997 July 24
* punch line intensities - but do not do last zone twice call plinin( ipPnunit(i) ) else if( chTime.ne.'LAST'.and. 1 (nzone.eq.1.or.(nzone/levery)*levery.eq.nzone ) ) then * this is middle of calculation * punch line intensities call plinin( ipPnunit(i) )
Adaptive logic for zoning failed in PDRs ionized by cool continua. This inhibited convergence of certain PDR models. Change line 381 of nextdr.f from
else if(hii.gt.0. .and. hi.gt.0. .and. OHIIoHI.gt.0. ) then to else if(hii.gt.0. .and. hi.gt.0. .and. OHIIoHI.gt.1e-15 ) then
1997 July 9.
Incorrect logic for capping changes in electron density. This inhibited convergence of models that went to the fully molecular limit. Make two changes to ionize.f. On line 458 change
if( enew.gt.eden*fac ) then to if( enew.gt.eold*fac ) then
and on line 467 change
else if( enew.lt.eden/fac ) then to else if( enew.lt.eold/fac ) then
The problem did not affect results but inhibited convergence in some extreme situations. the code's internal self checks detected these problems, resulting in an error condition. Thanks to Bhaswati Mookerjea for finding this problem. 1997 July 9.
Convergence of hydrogen Lya destruction -- HeI 2^3S photoionization problem oscillated. Change line 772 of routine addopac.f from
if( opac(i) .gt.1e-30 ) then to if( opac(i) .gt.1e-28 ) then
Thanks to Michele Thornley for finding this problem. 1997 July 7.
Data statement out of order in hypho.f. Data statements must be the last item among the declarations in FORTRAN 77, although most compilers do not insist on this. If yours does, move the data statement setting the variable "first" on line 43 to the line after the real declaration. The region in question might look like the following after this correction:
real zero,one,two,four,ten,con1 data first /.true./ data zero/0./ , one/1./ , two/2./,four/4./,ten/10./, 1 con1/0.8559/
Thanks to George Jacoby for finding this problem. 1997 June 30.
Peter Martin's extra lines, turned on with print all command, were not rescaled by the normalize command when the scale factor was not unity. Edit addln2, adding the common block norm after line 26:
integer norm real sclnrm character chNormLab*4 * norm is integer wavelength of emission line on normalize command * sclnrm is the scale factor for its appearance * chNormLab is optional label1 common/norm/ norm , sclnrm , chNormLab
next the (new) change line 40 from
scaldu = sumnew / snorm to *change 97 june 30, did not have sclnrm so only worked when lines not rescaled scaldu = sumnew / snorm * sclnrm
Thanks to Brian Espey for finding this problem. 1997 June 30.
C90.03 to C90.03a
Typo in setting nWerner. Change line 392 of scalar.f from 1958 to 1957. It only affected the table stars werner command on some platforms. This had no effect here in Lexington, but did blow up one of Brian Espey's models. Change the character variable at line 7 of version.f from ' 90.03 ' to '90.03a' (simple replace the space after the with with an 'a', being careful to keep the character string the same length. Versions of the code marked C90.03a and later have this correct already made.
Thanks to Brian Espey for finding this problem. 1997 May 13.