older/c9003_hotfix.htm
Home Up C94.00 Other versions Revision history The future etcetera Acknowledgements Contacts, Mailing list Links Site map, search

 

Home
Up

 

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.

C90.03b

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.

Hit Counter
Last changed 01/27/04.
Return to the Cloudy Home Page.
Copyright 1978-2003 Gary J. Ferland