 
      MODULE  MAC_CODE 
      private
      public mac, initialize_mac
      CONTAINS
! ********************************************************************
      SUBROUTINE INITIALIZE_MAC
! ********************************************************************
      use gms
      use mac_commons
      character(len=100),external:: dir
      integer :: ios
      external fbuff

      istdout=49
      clocation=element//'/'
      cauxdr='DIAG'//clocation//'macstd'
      cauxdr=dir(workdir,cauxdr)
      open(unit=istdout,file=cauxdr,iostat=ios)
!     if (ios /= 0)then
!        write(1,*) "Error sub initialize_mac :opening ", cauxdr
!        call fbuff()
!        stop
!     endif
      call output_open         ! Activates open files for OUTPUT
      call fbuff()
      if(idiag(8).eq.0)goto 4
          write(8,*) 'set data style linespoints'
          write(8,*) 'set yrange[0:1]'
          write(8,*) 'set xrange[0:70]'
          write(8,*) 'set xlabel "(eV)"'
          write(8,*) 'set ylabel "(Frational Population)"'
    4 continue
      call mac_initial  ! Atomic Kin. Initial.
      call energy_diff
      call rad_transdat
      frac=1.0d0
      call srd
      call srd_sum

      END SUBROUTINE INITIALIZE_MAC

! ******************************************************************** 
      SUBROUTINE MAC
! ******************************************************************** 
       use gms
       use mac_commons 
       implicit real*8(a-h,o-z)
       character(len=100),external:: dir
       Te_eV = src_Te_eV
       rNa   = src_rNa
!      zbar  = src_zbar
       rNe   = src_rNe
       iter  = isrc_iter

!t       clocation=element//'/'
       Te    = Te_eV/8.617385d-5

! -- BEGIN READ ONCE ----------------------------------------------------

      if(iReadFlag.eq.1) then
         call initialize_mac
!t        istdout=49
!t        cauxdr='DIAG'//clocation//'stdout'
!t        cauxdr=dir(workdir,cauxdr)
!t        open(unit=istdout,file=cauxdr)

        write(istdout,*) 'Te_eV =', src_Te_eV
        write(istdout,*) 'rNa   =', src_rNa
        write(istdout,*) 'zbar  =', src_zbar
        write(istdout,*) 'rNe   =', src_rNe
        write(istdout,*) 'iter  =', isrc_iter
        WRITE(istdout,*) 'itd   =', itd
        do i=1,num
           f(i)= 0.0
           fp(i)= 0.0
        enddo
!t call output_open was here 
!t       call initialize_mac
! ================================
!      call mac_initial  ! Atomic Kin. Initial.
!      call energy_diff
!      call rad_transdat
!      frac=1.0d0
!      call srd
!      call srd_sum
! ================================
     
! --- Time Dep. Calc. 
          if(itd.eq. 1)then
!            write(istdout,*)' in tindepfp itd=',itd
            cauxdr=clocation//'tindepfp'
            cauxdr=dir(workdir,cauxdr)
            write(istdout,*) cauxdr
            open(28,file=cauxdr,status='old',iostat=ierror,err=1000)
            do i=1,4
               read(28,*)
            enddo
            i=1
10          READ(28,*,END=20) indx1, fpp(i)
               WRITE(istdout,*) 'fpp(',i,')=',fpp(i)
               i=i+1
               GOTO 10
20          CONTINUE
            close(28)
            nread=i-1
           IF(modeltot.ne.nread)then
            write(istdout,*) 'ERROR MAC:modeltot=',modeltot,'.ne.',         &
     &      'nread=',nread
            stop
           ENDIF
          endif
! --- 
      iReadFlag=0
      write(istdout,*) 'READ ONCE COMPLETE'
      
      endif
! -- END READ ONCE ------------------------------------------------
! -- Time Dependent Code
!         write(istdout,*) 'itime=',itime
         if((itd .eq. 1).and.(itime.eq.1)) then ! itime=1 -> new time
          ! This section of code is not executed during the first 
          ! time step since itime =0.
          do i=1,modeltot
             fpp(i)=fp(i)
           enddo
!           write(istdout,*) 'Copy fpp = fp'
         endif
! ------------------------------
      call copyary 
      IF(ilte.eq.0)THEN
         CALL Stewart_Pyatt                  !NLTE self-cons.calc.
         CALL mac_kernal
        
      ELSE
         CALL mac_lte                        ! LTE self-cons. calc.
      ENDIF
      call converge(iconvrg,iconsum)
      
! ------------------------------------------------------------------
      if(iconvrg.eq.1)then
      kill_sum = 0
      do i= 1,num
         kill_sum = kill_lev(i) + kill_sum
      enddo
      if(itd.eq.0)then 
      call mac_output(iTempstep,iNa_step)
      endif
 
!     if(idiag(89).eq.0) goto 89
!     write(89,'( )')
! 89  continue
      if(idiag(90).eq.0) goto 90
         write(90,'( )')
  90  continue
      endif
      src_zbar  = zbar
      src_ntotal = modeltot
      do i=1,num
         src_fp(i)=fp(i)
      enddo
      return
1000  print*, 'CAN NOT OPEN FILE:',cauxdr
      write(istdout,*) 'CAN NOT OPEN FILE:',cauxdr
      close(istdout)
      stop 'CAN NOT OPENFILE tindepfp'
      END SUBROUTINE MAC

! ********************************************************************* 
      SUBROUTINE MAC_KERNAL 
! This program consitutes the kernal of Manolo's Atomic kinetic Code    
! ********************************************************************* 
      use mac_commons
      use gms
      implicit real*8 (a-h,o-z) 
!     parameter(num=200,mum=99,nnum=15,lnum=99,kmum=99) 
      dimension a      (num,num),                 &
  &             a_befor(num,num),                 &
  &             b      (num)

      dimension  avector(num*num)          !MINV                         
! --------------------------------------------------------------------  
!***************** BEGINNING OF ATOMIC KIN. KERNAL ******************   
! --------------------------------------------------------------------  
      if(isrd.eq.1) then 
! ===========================                                           
         call srd 
! ===========================                                           
      else 
         do i=1,modeltot 
            do j=1,modeltot 
                 rsrd(i,j) = 0.0d0 
            enddo 
         enddo 
      endif 
! -------------------------------------------                           
                                                                        
! ===========================================                           
      call rcoef_eceecd 
! ===========================================                           
      iconv = 0 
!********************************************************************** 
!****************** BEGINNING LEVEL OCCUPATION CALCS.****************** 
!********************************************************************** 
!     do  190 iter = 1,niter 
!       zbarp = zbar 
                                                                        
!       write (istdout,*)'   iter=',iter,'  zbar=',zbar                          
!       rNe=zbar*rNa 
! ----------------------------------------------------------------      
! -------------------------------------------                           
      if(impi.eq.1) then 
! ================================================================      
         call mpi 
! ================================================================      
      else 
         do i=1,modeltot 
            do j=1,modeltot 
               rmpi(i,j) = 0.0d0 
            enddo 
         enddo 
      endif 
                                                                        
! -------------------------------------------                           
      if(idr.eq.1) then 
! ===============================                                       
          call dr 
! ===============================                                       
      else 
         do i=1,modeltot 
            do j=1,modeltot 
               rdr(i,j) = 0.0d0 
            enddo 
         enddo 
      endif 
! --------------------------------------------------------------------- 
      if(irr.eq.1) then 
!       write (istdout,*) 'mac_kernal.f befor rcoef_pirr'                          
! ====================================================                  
         call rcoef_pirr   
         call rate_rr 
! ====================================================                  
      else 
         do i=1,modeltot 
            do j=1,modeltot 
               rrr(i,j) = 0.0d0 
            enddo 
         enddo 
      endif 
! --------------------------------------------------------------------- 
                                        !  The photoionizaton rate is ca
      if(ipi.eq.0) then 
                                        !  in rcoef_pirr.               
         do i=1,modeltot 
            do j=1,modeltot 
               rpi(i,j) = 0.0d0 
            enddo 
         enddo 
      endif 
! --------------------------------------------------------------------- 
                                                                        
! --------------------------------------------------------------------- 
      if(ieceecd.eq.1) then 
! ===========================================                           
        call rate_eceecd 
! ===========================================                           
      else 
         do i=1,modeltot 
            do j=1,modeltot 
               receecd(i,j) = 0.0d0 
            enddo 
         enddo 
      endif 
! --------------------------------------------------------------------- 
      if(ieciecr.eq.1) then 
! ===========================================                           
        call rcoef_eciecr 
        call rate_eciecr 
! ===========================================                           
      else 
         do i=1,modeltot 
            do j=1,modeltot 
               reciecr(i,j) = 0.0d0 
            enddo 
         enddo 
      endif 
!-----------------------------------------------------------------------
!***********************************************************************
! This region of the code produces the matrix A of rate coeff.         !
! The convension for all matries containing atomic rates:              !
! r(i,j): i = inital level, j = finial level.                          !
!       If i<j -> excitation or ionization                             !
!       If i>j -> deexcition or recombination                          !
!       If n the level of interest: r(n,m)->depopulating level n       !
!                                   r(m,n)->populating level n         !
!***********************************************************************
        ntotal = modeltot 
! ======================================================================
        call pressure_ionization 
! ======================================================================
                                                                        
        do n=1,num 
           do m=1,num 
              a(n,m)=0.0d0 
           enddo 
        enddo 
                                      ! n - is the level of interest.   
        do n = 1,ntotal 
                                      ! Incrimenting n -> marching down 
                                      ! the rows of the matrix A.       
! --------------------------------------------------------------------- 
! Solving for a(n,m) for m < n                                          
                                      ! n - is the level of interest.   
                                      ! m<n: pop. n from below.         
                                      ! n=1: has no lower levels        
          if(n.eq.1) goto 20 
                                      ! m,n: rates=excit or ioniz.      
          do m = 1,n-1 
              aux = receecd(m,n) + reciecr(m,n) 
              aux = aux + rmpi(m,n) + rpi(m,n) 
              a(n,m) = aux 
          enddo 
! --------------------------------------------------------------------- 
! Solving for a(n,n)                                                    
                                      ! n - is the level of interest.   
   20     aux = 0.0d0 
                                      ! n>m: depop. n to lower lev.     
                                      ! n=1: has no lower levels.       
          if (n.eq.1) goto 40 
                                      ! n,m: rates=deexcit or recomb.   
          do m = 1,n-1 
              aux = aux + receecd(n,m) + reciecr(n,m) 
              aux = aux + rsrd(n,m) + rrr(n,m) + rdr(n,m) 
          enddo 
                                      ! n - is the level of interest.   
                                      ! n<m: depop. n to higher lev.    
                                      ! n=ntotal: has no higher lev.    
   40     if (n.eq.ntotal) goto 60 
                                      ! n,m: rates=excit or ioniz.      
          do m = n+1,ntotal 
              aux = aux + receecd(n,m) + reciecr(n,m) 
              aux = aux + rmpi(n,m) + rpi(n,m) 
          enddo 
                                                                        
   60     a(n,n) = -aux 
! --------------------------------------------------------------------- 
! Solving for a(n,m) m > n                                              
                                      ! n - is the level of interest.   
                                      ! m>n: pop. level n from above.   
                                      ! n=ntotal: has no higher lev.    
          if (n.eq.ntotal) goto 80 
                                      ! m,n: rates=deexcit or recomb.   
          do m = n+1, ntotal 
              aux = receecd(m,n) + reciecr(m,n) 
              aux = aux + rsrd(m,n) + rrr(m,n) + rdr(m,n) 
              a(n,m) = aux 
          enddo 
   80     continue 
! --------------------------------------------------------------------- 
          enddo 
                                                                        
! **********************************************************************
                                                                        
! **************************************************************        
! This tests the form of the rate equations in creating A      !        
! by summing the colums of the matrix A.  The sum of the col.  !        
! coef. of A should be 0.0d0.  This implies that for every sink!        
! there is a source.                                           !        
! **************************************************************        
          if(idiag(30).eq.0) goto 13 
          write(30,*) '*****************************','iter=',iter 
          do m = 1,ntotal 
            auxp = 0.0d0 
            auxn = 0.0d0 
            do n=1,ntotal 
               if(a(n,m).le.0.0d0)then 
                                        ! adding all the pos. num.      
                 auxn = auxn + a(n,m) 
               else 
                                        ! adding all the neg. num.      
                 auxp = auxp + a(n,m) 
               endif 
            enddo 
            b(m) = auxn + auxp 
            write(30,*) 'A_col_sum(',m,')=',b(m) 
          enddo 
   13     continue 
! **********************************************************************
!FOR PRINTING OUT OF A-MATRIX BEFORE NORMAL. COND.                                                                         
!       do n =1,ntotal
!           write(6,'14(1x,1pe10.3)') (a(n,m),m=1,ntotal)
!       enddo 
      if((idiag(37).eq.0).and.(idiag(38).eq.0))goto 14 
      write (istdout,*) 'test_kernal' 
! ==================================                                    
      call test_kernal(ntotal) 
! ==================================                                    
   14 continue 
   if(itd.eq.0)then     ! Steady state
     do m = 1,ntotal 
       if(kill_lev(m).eq.0) a(m,m)= 1.0d0 
         a(ntotal,m) = 1.0d0 
         b(m) = 0.0d0 
     enddo 
     b(ntotal) = 1.0d0 
   elseif(itd.eq.1)then ! Time Dep
      do m=1,ntotal
        a(ntotal,m)=1.0d0
        b(m)=fpp(m)
      enddo
      b(ntotal)=1.0d0

      do n=1,ntotal-1
        do m=1,ntotal
          a(n,m)=-deltat*a(n,m) 
          if(n.eq.m) a(n,m)=1.0+a(n,m)
        enddo
      enddo
   else
     write(istdout,*) 'itd does not equal appropreate val HALTING CODE'
     STOP 'itd does not equal appropreate val HALTING CODE'
   endif                                                                   
! ================================================================      
!         call matrix_compression(a,levels,n_kills,ntotal,num)          
! ================================================================      
          if(idiag(35) .eq. 0) goto 8 
           write(35,*) '                      ' 
            do ii=1,ntotal 
              do jj=1,ntotal 
                 write(35,*) 'a(',ii,',',jj,')=',a(ii,jj) 
              enddo 
            enddo 
    8     continue 
                                                                        
          if(idiag(31) .eq. 0) goto 6 
            do ii=1,ntotal 
              do jj=1,ntotal 
                 a_befor(ii,jj)=a(ii,jj) 
                 write(31,*) 'a(',ii,',',jj,')=',a(ii,jj) 
              enddo 
            enddo 
    6     continue 
! ==========================================================            
        call gaussj(a,ntotal,num,b,1,1) 
         do ii=1,ntotal 
            f(ii)=b(ii)
         enddo 
! ==========================================================            
                                                                        
! ===================================================================== 
!FOR USE WITH MINV                                                      
!       nn=0                                                            
!       do j=1,ntotal                                                   
!          do i=1,ntotal                                                
!             nn=nn+1                                                   
!             avector(nn)=a(i,j)                                        
!          enddo                                                        
!       enddo                                                           
! ==========================================================            
!       call m_minv(avector,ntotal,detlog)                              
! ==========================================================            
!      write (istdout,*) 'detlog=',detlog,' rNe=',rNe                              
!       nn=0                                                            
!       do j=1,ntotal                                                   
!         do i=1,ntotal                                                
!             nn=nn+1                                                   
!            write(43,*)'inv_avector(',nn,')=',avector(nn)              
!             a(i,j)=avector(nn)                                        
!            write(44,*)'inv_av_a(',i,',',j,')=',a(i,j)                 
!          enddo                                                        
!       enddo                                                           
!       do n=1,ntotal                                                   
!          aux = 0.0d0                                                  
!          do m = 1,ntotal                                              
!             aux = aux + a(n,m) * b(m)   !used with minv               
!          enddo                                                        
!          f(n) = aux                                                   
!       enddo                                                           
!FOR USE WITH MINV                                                      
! ===================================================================== 
       if(idiag(32).eq.0) goto 7 
         do ii=1,ntotal 
           do jj=1,ntotal 
             write(32,*) 'inv_a(',ii,',',jj,')=',a(ii,jj) 
           enddo 
         enddo 
    7  continue 
                                                                        
                                                        ! needs a_befor(
       if((idiag(31).eq.1).and.(idiag(36).eq. 1)) then 
           do jj=1,ntotal 
             do ll=1,ntotal 
               sum=0.0d0 
               do kk=1,ntotal 
                  xaux=a_befor(ll,kk)*a(kk,jj) 
                  sum=sum+xaux 
               enddo 
               write(36,*) 'a_befor*a(',ll,',',jj,')=',sum 
            enddo 
           enddo 
       endif 
! --------------------------------------------------------------------- 
          sumf = 0.0d0 
          zbar = 0.0d0 
          do n = 1,ntotal 
             sumf = sumf + f(n) 
             zbar = zbar + f(n) * zeff(n) 
          enddo 
! ----------------------------------------------------------------------
!                                     ! delta - represents a threshold v
!     if (iter.eq.1) goto 170         ! at which the comparision for the
!       icheck = 0                    ! not applied. Therefore f(i)<delt
!     do n =1,ntotal                  ! ignored. 
!       reld = dabs(f(n)-fp(n))/(fp(n)+delta) 
!       if(reld.lt.eps) goto 160 
!       icheck = icheck +1 
! 160 enddo 
!     if (icheck.eq.0) goto 200       ! goto 200 -> conv. criteria reach 
! 170 do n = 1,ntotal 
!        fp(n) = f(n) 
!     enddo 
!     zbar = (zbar + zbarp) * 0.5d0 
!     zbarp = zbar 
!                                     ! 190 concludes the self interatio
! 190 continue                        ! do loop 
!                                                                       
!*********************** END LEVEL OCCUPATION CALCS.********************
! 200     rNe = zbar*rNa 
                                                                        
!********************** END OF ATOMIC KIN. KERNAL ********************* 
!     if(iter.eq.101) iconv =1 

      END SUBROUTINE MAC_KERNAL                                          





!********************************************************************** 
      SUBROUTINE CHAR_NUM(char,int) 
!**********************************************************************
      IMPLICIT REAL*8 (a-h,o-z) 
        CHARACTER*1 ones(0:9),char 
        DATA ones /"0","1","2","3","4","5","6","7","8","9"/ 
        DO i=0,9 
          IF(char.eq.ones(i))THEN 
          int=i 
          GOTO 10 
           ENDIF 
        ENDDO 
   10 CONTINUE 
      END SUBROUTINE CHAR_NUM                                          
                                                                        
!*****************************************************
      SUBROUTINE gaussj(a,n,np,b,m,mp) 
!*****************************************************
      INTEGER m,mp,n,np,NMAX 
      DOUBLE PRECISION a(np,np),b(np,mp) 
      PARAMETER (NMAX=500) 
      INTEGER i,icol,irow,j,k,l,ll,indxc(NMAX),indxr(NMAX),ipiv(NMAX) 
      DOUBLE PRECISION big,dum,pivinv 
      do 11 j=1,n 
        ipiv(j)=0 
   11 continue 
      do 22 i=1,n 
        big=0.d0 
        do 13 j=1,n 
          if(ipiv(j).ne.1)then 
            do 12 k=1,n 
              if (ipiv(k).eq.0) then 
                if (abs(a(j,k)).ge.big)then 
                  big=abs(a(j,k)) 
                  irow=j 
                  icol=k 
                endif 
              else if (ipiv(k).gt.1) then 
                pause 'singular matrix in gaussj' 
              endif 
   12       continue 
          endif 
   13   continue 
        ipiv(icol)=ipiv(icol)+1 
        if (irow.ne.icol) then 
          do 14 l=1,n 
            dum=a(irow,l) 
            a(irow,l)=a(icol,l) 
            a(icol,l)=dum 
   14     continue 
          do 15 l=1,m 
            dum=b(irow,l) 
            b(irow,l)=b(icol,l) 
            b(icol,l)=dum 
   15     continue 
        endif 
        indxr(i)=irow 
        indxc(i)=icol 
        if (a(icol,icol).eq.0.d0) pause 'singular matrix in gaussj' 
        pivinv=1.d0/a(icol,icol) 
        a(icol,icol)=1.d0 
        do 16 l=1,n 
          a(icol,l)=a(icol,l)*pivinv 
   16   continue 
        do 17 l=1,m 
          b(icol,l)=b(icol,l)*pivinv 
   17   continue 
        do 21 ll=1,n 
          if(ll.ne.icol)then 
            dum=a(ll,icol) 
            a(ll,icol)=0.d0 
            do 18 l=1,n 
              a(ll,l)=a(ll,l)-a(icol,l)*dum 
   18       continue 
            do 19 l=1,m 
              b(ll,l)=b(ll,l)-b(icol,l)*dum 
   19       continue 
          endif 
   21   continue 
   22 continue 
      do 24 l=n,1,-1 
        if(indxr(l).ne.indxc(l))then 
          do 23 k=1,n 
            dum=a(k,indxr(l)) 
            a(k,indxr(l))=a(k,indxc(l)) 
            a(k,indxc(l))=dum 
   23     continue 
        endif 
   24 continue 
      return 
      END SUBROUTINE GAUSSJ

! ***************************************************************
      SUBROUTINE OPEN_DIAG 
! ***************************************************************
      use mac_commons
      use gms
      implicit real*8 (a-h,o-z) 
      character(len=100),external :: dir
! ---------------------------------------                               
      cauxdr=clocation//'diag.in'
      cauxdr=dir(workdir,cauxdr)
      open(unit=10,file=cauxdr) 
      read(10,*) 
! ----------------------------                                          
                       ! header                                         
      read(10,*) 
                       ! --- !                                          
      read(10,*) 
                       ! subheading                                     
      read(10,*) 
      do i = 30,39 
         read(10,900) idiag(i) 
      enddo 
! ----------------------------                                          
      read(10,*) 
      read(10,*) 
      read(10,*) 
      do i = 40,49 
         read(10,900) idiag(i) 
      enddo 
! ----------------------------                                          
      read(10,*) 
      read(10,*) 
      read(10,*) 
      do i = 50,59 
         read(10,900) idiag(i) 
      enddo 
! ----------------------------                                          
      read(10,*) 
      read(10,*) 
      read(10,*) 
      do i = 60,65 
         read(10,900) idiag(i) 
      enddo 
! ----------------------------                                          
      close(10) 
                                                                        
! **************************** DIRECTORY ************************* !    
! ---------------------------------------------------------------- !    
! ATOMIC-DATA:                                                          
!     open(unit=40,file='DIAG/read_warn')      ! all read cs routines   
!     open(unit=41,file='DIAG/read_ece')       ! read_ece.f             
!     open(unit=42,file='DIAG/read_eci')       ! read_eci.f             
!     open(unit=44,file='DIAG/read_rr')        ! read_rr.f              
!     open(unit=45,file='DIAG/atomicstru')     ! read_atom.f            
!     open(unit=46,file='DIAG/gf')             ! read_atom.f            
!     open(unit=47,file='DIAG/radat')          ! rad_trans
! ---------------------------------------------------------------- !    
! RATES:                                                                
!     open(unit=51,file='DIAG/receecd')        ! rate_ece.f             
!     open(unit=52,file='DIAG/reciecr')        ! rate_eci.f             
!     open(unit=53,file='DIAG/rpi')            ! rcoef_rr.f             
!     open(unit=54,file='DIAG/rrr')            ! rate_rr.f              
!     open(unit=55,file='DIAG/rdr')            ! dr.f                   
!     open(unit=56,file='DIAG/rsrd')           ! srd.f                  
!     open(unit=57,file='DIAG/rmpi')           ! mpi.f                  
! ---------------------------------------------------------------- !    
! COEF:                                                                 
!     open(unit=61,file='DIAG/coef_eced')      ! rcoef_ece.f            
!     open(unit=62,file='DIAG/coef_ecir')      ! rcoef_eci.f            
!     open(unit=64,file='DIAG/coef_rr')        ! rate_rr.f              
! ---------------------------------------------------------------- !    
! MATRIX A:                                                             
!     open(unit=30,file='DIAG/a_col_sum')      ! mac_kernal.f           
!     open(unit=31,file='DIAG/A')              ! mac_kernal.f           
!     open(unit=32,file='DIAG/INV_A')          ! mac_kernal.f           
!     open(unit=33,file='DIAG/INDT')           ! mac_kernal.f           
!     open(unit=34,file='DIAG/A_STATE')        ! mac_kernal.f           
!     open(unit=35,file='DIAG/A')              ! mac_kernal.f           
!     open(unit=36,file='DIAG/AxINV_A')        ! mac_kernal.f           
!     open(unit=37,file='DIAG/col_sum.final')  ! test_kernal.f          
!     open(unit=38,file='DIAG/col_sum.indiv')  ! test_kernal.f          
! **************************************************************** !    
                                                                        
! **************************************************************** !    
! ATOMIC-DATA:                                                          
      if(idiag(40).eq.0)goto 40 
                                                ! all cs read routines  
        cauxdr='DIAG'//clocation//'read_warn'
        cauxdr=dir(workdir,cauxdr)
       open(unit=40,file=cauxdr) 
   40 continue 
      if(idiag(41).eq.0)goto 41 
                                               ! read_ece.f             
        cauxdr='DIAG'//clocation//'read_ece'
        cauxdr=dir(workdir,cauxdr)
      open(unit=41,file=cauxdr) 
   41 continue 
      if(idiag(42).eq.0)goto 42 
                                               ! read_eci.f             
        cauxdr='DIAG'//clocation//'read_eci'
        cauxdr=dir(workdir,cauxdr)
      open(unit=42,file=cauxdr) 
   42 continue 
      if(idiag(44).eq.0)goto 44 
                                               ! read_rr.f              
        cauxdr='DIAG'//clocation//'read_rr'
        cauxdr=dir(workdir,cauxdr)
      open(unit=44,file=cauxdr) 
   44 continue 
      if(idiag(45).eq.0)goto 45 
                                               ! read_atom.f            
        cauxdr='DIAG'//clocation//'atomicstru'
        cauxdr=dir(workdir,cauxdr)
      open(unit=45,file=cauxdr) 
   45 continue 
      if(idiag(46).eq.0)goto 46 
                                               ! read_atom.f            
        cauxdr='DIAG'//clocation//'gf'
        cauxdr=dir(workdir,cauxdr)
      open(unit=46,file=cauxdr) 
   46 continue
      if(idiag(47).eq.0)goto 47
                                               ! rad_transdat 
        cauxdr='DIAG'//clocation//'radat'
        cauxdr=dir(workdir,cauxdr)
      open(unit=47,file=cauxdr)
   47 continue
 
! ----------------------------------------------------------------      
! RATES:                                                                
      if(idiag(51).eq.0)goto 51 
                                               ! rate_ece.f             
        cauxdr='DIAG'//clocation//'receecd'
        cauxdr=dir(workdir,cauxdr)
      open(unit=51,file=cauxdr) 
   51 continue 
      if(idiag(52).eq.0)goto 52 
                                               ! rate_eci.f             
        cauxdr='DIAG'//clocation//'reciecr'
        cauxdr=dir(workdir,cauxdr)
      open(unit=52,file=cauxdr) 
   52 continue 
      if(idiag(53).eq.0)goto 53 
                                               ! rcoef_rr.f             
        cauxdr='DIAG'//clocation//'rpi'
        cauxdr=dir(workdir,cauxdr)
      open(unit=53,file=cauxdr) 
   53 continue 
      if(idiag(54).eq.0)goto 54 
                                               ! rate_rr.f              
        cauxdr='DIAG'//clocation//'rrr'
        cauxdr=dir(workdir,cauxdr)
      open(unit=54,file=cauxdr) 
   54 continue 
      if(idiag(55).eq.0)goto 55 
                                               ! dr.f                   
        cauxdr='DIAG'//clocation//'rdr'
        cauxdr=dir(workdir,cauxdr)
      open(unit=55,file=cauxdr) 
   55 continue 
      if(idiag(56).eq.0)goto 56 
                                               ! srd.f                  
        cauxdr='DIAG'//clocation//'rsrd'
        cauxdr=dir(workdir,cauxdr)
      open(unit=56,file=cauxdr) 
   56 continue 
      if(idiag(57).eq.0)goto 57 
                                               ! mpi.f                  
        cauxdr='DIAG'//clocation//'rmpi'
        cauxdr=dir(workdir,cauxdr)
      open(unit=57,file=cauxdr) 
   57 continue 
! ----------------------------------------------------------------      
! COEF:                                                                 
      if(idiag(61).eq.0)goto 61 
                                               ! rcoef_ece.f            
        cauxdr='DIAG'//clocation//'coef_eced'
        cauxdr=dir(workdir,cauxdr)
      open(unit=61,file=cauxdr) 
   61 continue 
      if(idiag(62).eq.0)goto 62 
                                               ! rcoef_eci.f            
        cauxdr='DIAG'//clocation//'coef_ecir'
        cauxdr=dir(workdir,cauxdr)
      open(unit=62,file=cauxdr) 
   62 continue 
      if(idiag(64).eq.0)goto 64 
                                               ! rate_rr.f              
        cauxdr='DIAG'//clocation//'coef_rr'
        cauxdr=dir(workdir,cauxdr)
      open(unit=64,file=cauxdr) 
   64 continue 
! ----------------------------------------------------------------      
! MATRIX A:                                                             
      if(idiag(30).eq.0)goto 30 
                                               ! mac_kernal.f           
        cauxdr='DIAG'//clocation//'a_col_sum'
        cauxdr=dir(workdir,cauxdr)
      open(unit=30,file=cauxdr) 
   30 continue 
      if(idiag(31).eq.0)goto 31 
                                               ! mac_kernal.f           
        cauxdr='DIAG'//clocation//'a_before'
        cauxdr=dir(workdir,cauxdr)
      open(unit=31,file=cauxdr) 
   31 continue 
      if(idiag(32).eq.0)goto 32 
                                               ! mac_kernal.f           
        cauxdr='DIAG'//clocation//'a_inv'
        cauxdr=dir(workdir,cauxdr)
      open(unit=32,file=cauxdr) 
   32 continue 
      if(idiag(33).eq.0)goto 33 
                                               ! mac_kernal.f           
        cauxdr='DIAG'//clocation//'indt'
        cauxdr=dir(workdir,cauxdr)
      open(unit=33,file=cauxdr) 
   33 continue 
      if(idiag(34).eq.0)goto 34 
                                               ! mac_kernal.f           
        cauxdr='DIAG'//clocation//'istate'
        cauxdr=dir(workdir,cauxdr)
      open(unit=34,file=cauxdr) 
   34 continue 
      if(idiag(35).eq.0)goto 35 
                                               ! mac_kernal.f           
        cauxdr='DIAG'//clocation//'a_state'
        cauxdr=dir(workdir,cauxdr)
      open(unit=35,file=cauxdr) 
   35 continue 
      if(idiag(36).eq.0)goto 36 
                                               ! mac_kernal.f           
        cauxdr='DIAG'//clocation//'axinv_a'
        cauxdr=dir(workdir,cauxdr)
      open(unit=36,file=cauxdr) 
   36 continue 
      if(idiag(37).eq.0)goto 37 
                                               ! test_kernal.f          
        cauxdr='DIAG'//clocation//'col_sum.final'
        cauxdr=dir(workdir,cauxdr)
      open(unit=37,file=cauxdr) 
   37 continue 
      if(idiag(38).eq.0)goto 38 
                                               ! test_kernal.f          
        cauxdr='DIAG'//clocation//'col_sum.indiv'
        cauxdr=dir(workdir,cauxdr)
      open(unit=38,file=cauxdr) 
   38 continue 
                                                                        
  900 format(15x,i1)
                                                                        
      END SUBROUTINE OPEN_DIAG
!*****************************************************************      
      SUBROUTINE DR 
!*****************************************************************      
!* Te -> [K]                                                            
      use mac_commons
      IMPLICIT REAL*8 (a-h,o-z) 
!     PARAMETER(num=200,mum=99,lnum=99) 
                                 !Nvalshell is the value of the         
                                 !principle quantum number.             
      DIMENSION Zeffdr(mum),                                            &
     &          rdr1(num,num)                                           
                                                                        
      DO i=1,mum 
          zeffdr(i)=DFLOAT(i) 
       ENDDO 
      DO i=1,modeltot 
         DO j=1,modeltot 
           RDR(i,j) = 0.0d0 
         ENDDO 
      ENDDO 
                                                                        
      boltzmann=8.617385D-5 
      t=Te/boltzmann 
      t=t/1.0d6 
                                                                        
       DO II=1,nterms-2 
          k=igrounds(ii) 
         i=igrounds(ii+1) 
         rnui=DFLOAT(Nvalshell(i)) 
         z=zeffdr(ii) 
         Blocal=z**0.5d0*(z+1.0d0)**2.5d0*(z**2+13.4d0)**(-0.5d0) 
         alittle=1.0d0+0.015*(z**3.0d0)*(z+1.0d0)**(-2) 
         sum=0.0d0 
         DO j=i,(igrounds(ii+2)-1) 
            rnuj=DFLOAT(Nvalshell(j)) 
            epslocal=(rnui**(-2))-(rnuj**(-2)) 
            IF(epslocal.le.1.0d-20)GOTO 21 
            xlocal=(z+1.0d0)*epslocal 
            Alocal=xlocal**0.5d0/(1.0d0+0.105*xlocal+0.015*xlocal**2) 
            Arg=0.158d0*((z+1.0d0)**2)*epslocal*(t**(-1))*(alittle**(-1)) 
              Sum=(gf(i,j)/(frac(i)*ig(i)))*Alocal*DEXP(-Arg)+Sum 
   21       CONTINUE 
         ENDDO 
         drcoef=3.0d-12*(t**(-1.5))*Blocal*Sum 
         rdr(i,k)=drcoef*rNe 
       ENDDO 
                                                                        
                                                                        
                                                                        
! -----------------------------------------------------                 
! DIAGNOSTICS                                                           
                                                                        
      if(idiag(55).eq.0)goto 100 
      do i=1,modeltot 
        do j=1,modeltot 
           write(55,*) i,j,RDR(i,j) 
        enddo 
      enddo 
  100 continue 
      return 
      END SUBROUTINE DR
                                                                        
                                                                        
!*******************************************************************    
       SUBROUTINE ENERGY_DIFF 
!******************************************************************     
      use mac_commons
      implicit real*8 (a-h,o-z) 
      j=1 
      do i=1,modeltot 
        if(i.lt.igrounds(j+1)) then 
          eng_diff(i) = rioneng(j)-energy(i) 
        else 
          j=j+1 
          eng_diff(i) = rioneng(j)-energy(i) 
        endif 
      enddo 
      END SUBROUTINE ENERGY_DIFF                                          

!***********************************************************************
!  This subroutine produces the filenames of the output files          *
      SUBROUTINE FILENAMES(prefix,fnames,sufix) 
!***********************************************************************
!*  NOTE: THIS SUBROUTINE WILL ONLY PRODUCE 999 FILES!                  
       implicit real*8 (a-h,o-z) 
                                                                        
       character*1 ones(10) 
       character*1 tens(10) 
       character*1 hundereds(10) 
       character*8 prefix 
       character*3 sufix 
       character*14 fnames(0:999) 
                                                                        
       data ones     /"0","1","2","3","4","5","6","7","8","9"/ 
       data tens     /"0","1","2","3","4","5","6","7","8","9"/ 
       data hundereds/"0","1","2","3","4","5","6","7","8","9"/ 
                                                                        
       jj=-1 
                                                                        
       do 20 i=1,10 
       do 30 j=1,10 
       do 40 k=1,10 
       jj=jj+1 
                                                                        
       fnames(jj)=prefix//hundereds(i)//tens(j)//ones(k)//sufix 
                                                                        
   40  continue 
   30  continue 
   20  continue 
      END SUBROUTINE FILENAMES 
                                                                        
                                                                        
! ----------------------------------------------------------------------
         DOUBLE PRECISION FUNCTION GMZBAR(tev) 
! ----------------------------------------------------------------------
! UNIT NUMBERS: 30.                                                     
! gmzbar.input must contain den., temp, and zbar                        
! NOTE: THIS FUNCTION REQUIRES THE COMPILER TO RETAIN ITS VALUES        
! This function returns  the interpolated value of zbar when it is given
! the temperature (eV) and the total atom number density.               
! This fuction reads the MAC output file zbar.in                        
                                                                        
      use gms
      use mac_commons
      implicit real*8(a-h,o-z) 
      dimension zzbar(251,33) 
       CHARACTER*1 char1,char2 
      character(len=100), external::dir
      data iflag/0/ 
                                             ! Temperature epsilons     
      epstmp = 1.0d-4 
                                             ! density epsilons         
      deltag = 1.0d14 
                                                                        
      t_low = 0.3d0 
      t_high = 50.0d0 
      d_low = 1.0d16 
      d_high = 1.0d24 
                                                                        
! ----------------------------------------------------------------------
! Setting up data array (read once):                                    
                                                                        
                                             ! allows for inputfile to  
      if(iflag.eq.1) go to 30 
                                             ! be read once.            
        cauxdr='DIAG'//clocation//'gmzbar.lte.in'
        cauxdr=dir(workdir,cauxdr)
      open(unit=70,file=cauxdr) 
      iflag=1 
      READ(70,'(a1,a1)') char1,char2 
      DO WHILE ((CHAR1.EQ.'#').AND.(CHAR2.EQ.'#')) 
          READ(70,'(a1,a1)') char1,char2 
       ENDDO 
       backspace(70) 
                                                                        
                                            ! number of temperature poin
      read(70,900) n_t_pts 
                                            ! number of denisty points. 
      read(70,900) n_d_pts 
       write (istdout,*) 'n_t_pt=',n_t_pts 
       write (istdout,*) 'n_d_pt=',n_d_pts 
                                                                        
!       read(70,901) tmphigh                  ! highest temp need cont. 
!       read(70,901) rnahigh                  ! highest den. need cont. 
                                                                        
                                                  ! slope used in temp. 
       Slope_T = (DBLE(n_t_pts-1))/(t_high-t_low) 
       Slope_D = (DBLE(n_d_pts-1))/(DLOG10(d_high)-DLOG10(d_low)) 
                                        ! slope used in dens. int. map  
                                                                        
                                             ! j is the index assoc. den
       do j=1,n_d_pts 
                                             ! i is the index assoc. tem
          do i=1,n_t_pts 
             read(70,*,END=35) auxtev,auxden,zzbar(i,j) 
            enddo 
          read(70,*) 
       enddo 
                                                                        
   35  continue 
       close(70) 
   30  continue 
                                                                        
! ----------------------------------------------------------------------
! Solution for data out of range:                                       
                                                                        
!      if((rna.gt.rnahigh).and.(tev.gt.tmphigh)) then  ! used when mac  
!       WRITE(16,*)'warn: BEYOND ZBAR LIMIT',            ! stops prematu
!     &              rna,'>?',rnahigh,' and ',         ! due to continuu
!     &              tev,'>?',tmphigh                                   
!      endif                                                            
                                                                        
      tevp = tev 
      rnap = rna 
                                                                        
      if(tev.le.t_low) then 
!        WRITE(16,*) 'warn: tev below range',tev                        
        tev = t_low + epstmp  
      endif 
      if(tev.ge.t_high) then 
!        WRITE(16,*)'warn: tev above range',tev                         
        tev = t_high - epstmp 
      endif 
      if(rna.le.d_low) then 
!        WRITE(16,*)'warn: rNa below range',rna                         
        rNa = d_low + deltag 
      endif 
      if(rna.ge.d_high) then 
!        WRITE(16,*)'warn: rNa above range',rna                         
        rna = d_high - deltag
      endif 
                                                                        
! ----------------------------------------------------------------------
! Temperature mapping to interger space:                                
                                                                        
! Slope = 5.0505051                                                     
         Tev_map = Slope_T*(tev-t_low)+1.0d0 
                                               ! DDINT - trunctate      
         Tev_lt  = INT(Tev_map) 
         iTev_lt = INT(Tev_lt) 
         iTev_gt = iTev_lt + 1 
                                                                        
! ----------------------------------------------------------------------
! Density mapping to integer space:                                     
                                                                        
         rNa_log = DLOG10(rNa) 
          rNa_log_map = Slope_D*(rNa_log - DLOG10(d_low)) + 1.0d0 
         rNa_lt = INT(rNa_log_map) 
         iNa_lt = INT(rNa_lt) 
         iNa_gt = iNa_lt + 1 
                                                                        
! ----------------------------------------------------------------------
! Interpolating zbar between iTe_lt & iTe_gt on curve iNa_lt:           
                                                                        
         z_lt = zzbar(iTev_lt,iNa_lt) 
         z_gt = zzbar(iTev_gt,iNa_lt) 
         slope = z_gt - z_lt 
         zbar_lower = slope*(Tev_map - dfloat(iTev_lt)) + z_lt 
! ----------------------------------------------------------------------
! Interpolating zbar between iTe_lt & iTe_gt on curve iNa_gt:           
                                                                        
         if(iNa_gt.gt.n_d_pts)                                          &
     &     write (istdout,*) 'HAULT:gmzbarnew.f iNa_gt > n_d_pts'                  
         if(iTev_gt.gt.n_t_pts)                                         &
     &     write (istdout,*) 'HUALT:gmzbarnew.f iTev_gt > n_t_pts'                 
         z_lt = zzbar(iTev_lt,iNa_gt) 
         z_gt = zzbar(iTev_gt,iNa_gt) 
                                                                        
         slope = z_gt - z_lt 
         zbar_upper = slope*(Tev_map - dfloat(iTev_gt)) + z_gt 
! ----------------------------------------------------------------------
!                                                                       
                                                   ! where 1.00 is the d
         slope = (zbar_upper - zbar_lower) 
                                                   ! in zbar_log between
                                                   ! adjacent points.   
                                                                        
         gmzbar = slope*(rNa_log_map - dfloat(iNa_lt)) + zbar_lower 
                                                                        
        tev = tevp 
        rna = rnap 
                                                                        
         return 
  900    format(3x,I10) 
  901    format(1x,1pE11.4) 
      END FUNCTION GMZBAR                                          
! *********************************************                         
      SUBROUTINE KILLED_LEV 
! *********************************************                         
      use mac_commons
      implicit real*8 (a-h,o-z) 
!     parameter(num=200,mum=99,nnum=15,lnum=99,kmum=99) 
                                                                        
      do i=1,modeltot 
         kill_lev(i)=1 
         r_stat_wt(i)=1.0d0 
         frac(i)=1.0 
      enddo 
                                                                        
      do i=1,modeltot 
         if(Eng_lower.gt.eng_diff(i)) then 
           kill_lev(i)=0 
                                                                        
                                                                     !re
          endif 
                                                !remaining              
        if(kill_lev(modeltot-1).eq.0) then 
           write (istdout,*) 'warn: too many levels removed do to lowering' 
            stop 'killed_lev'
         endif 
      enddo 
! -----------------------------------------------                       
! Gradual removal of statistical weight of the nearest surviving level  
                                                                        
      do i = 1,modeltot-1 
         if((kill_lev(i).eq.1).and.(kill_lev(i+1).eq.0)) then 
         frac(i)=(eng_diff(i)-Eng_lower)/(energy(i+1)-energy(i)) 
         endif 
      enddo 
! -------------------------------------------------------               
      END SUBROUTINE KILLED_LEV                                          

!****************************************************************       
      FUNCTION  GAMMLN(xx) 
!****************************************************************       
      implicit real*8(a-h,o-z) 
      dimension cof(6) 
      SAVE cof,stp 
      DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,            &
     &24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,    &
     &-.5395239384953d-5,2.5066282746310005d0/                          
      x=xx 
      y=x 
      tmp=x+5.5d0 
      tmp=(x+0.5d0)*log(tmp)-tmp 
      ser=1.000000000190015d0 
      do 11 j=1,6 
        y=y+1.d0 
        ser=ser+cof(j)/y 
   11 continue 
      gammln=tmp+log(stp*ser/x) 
      return 
      END FUNCTION GAMMLN                                          
                                                                        
!****************************************************************       
      SUBROUTINE M_GAULAG(x,w,n,alf) 
!****************************************************************       
      IMPLICIT REAL*8(A-H,O-Z) 
      DIMENSION w(n),x(n) 
      PARAMETER (EPS=3.D-14,MAXIT=10) 
!U    USES m_gammln                                                     
      do 13 i=1,n 
        if(i.eq.1)then 
          z=(1.+alf)*(3.+.92*alf)/(1.+2.4*n+1.8*alf) 
        else if(i.eq.2)then 
          z=z+(15.+6.25*alf)/(1.+.9*alf+2.5*n) 
        else 
          ai=1.0*i-2.0 
          z=z+((1.+2.55*ai)/(1.9*ai)+1.26*ai*alf/(1.+3.5*ai))*          &
     &(z-x(int(1.0*i-2.0)))/(1.+.3*alf)                                 
        endif 
        do 12 its=1,MAXIT 
          p1=1.d0 
          p2=0.d0 
          do 11 j=1,n 
            p3=p2 
            p2=p1 
            p1=((2.0*j-1.0+alf-z)*p2-(1.0*j-1.0+alf)*p3)/(1.0*j) 
   11     continue 
          pp=(n*p1-(1.0*n+alf)*p2)/z 
          z1=z 
          z=z1-p1/pp 
          if(dabs(z-z1).le.EPS)goto 1 
   12   continue 
        pause 'too many iterations in m_gaulag' 
    1   x(i)=z 
        w(i)=-dexp(gammln(alf+n)-gammln(dfloat(n)))/(pp*n*p2) 
   13 continue 
      return 
      END SUBROUTINE M_GAULAG                                   
                                                                        
!********************************************************************** 
         SUBROUTINE M_MINV(a,n,d) 
!********************************************************************** 
         IMPLICIT REAL*8(A-H,O-Z) 
         parameter(num=200) 
!                                                                       
! ENTRADA:   A , MATRIZ CUADRADA DE ORDEN N                             
!            L,M , VECTORES AUXILIARES DE CALCULO DE ORDEN N            
!                                                                       
! SALIDA:    A , MATRIZ CUADRADA INVERSA DE A                           
!            D , LOG DEL MODULO DEL DETERMINANTE DE A                   
!                                                                       
!        DIMENSION A(1),L(400),M(400)                                   
         DIMENSION A(num*num),L(num*num),M(num*num) 
!                                                                       
! UBICACION EN A DEL ELEMENTO DE MODULO MAXIMO                          
!                                                                       
         D=0.0D0 
         NK=-N 
         DO 80 K=1,N 
         NK=NK+N 
         L(K)=K 
         M(K)=K 
         KK=NK+K 
         BIGA=A(KK) 
         DO 20 J=K,N 
         IZ=N*(J-1) 
         DO 20 I=K,N 
         IJ=IZ+I 
         IF(DABS(BIGA)-DABS(A(IJ))) 15,20,20 
   15    BIGA=A(IJ) 
         L(K)=I 
         M(K)=J 
   20    CONTINUE 
!                                                                       
! INTERCAMBIO DE FILAS                                                  
!                                                                       
         J=L(K) 
         IF(J-K) 35,35,25 
   25    KI=K-N 
         DO 30 I=1,N 
         KI=KI+N 
         HOLD=-A(KI) 
         JI=KI-K+J 
         A(KI)=A(JI) 
   30    A(JI)=HOLD 
!                                                                       
! INTERCAMBIO DE COLUMNAS                                               
!                                                                       
   35    I=M(K) 
         IF(I-K) 45,45,38 
   38    JP=N*(I-1) 
         DO 40 J=1,N 
         JK=NK+J 
         JI=JP+J 
         HOLD=-A(JK) 
         A(JK)=A(JI) 
   40    A(JI)=HOLD 
!                                                                       
! DIVISION DE LA COLUMNA POR EL PIVOTE = (-1)*BIGA                      
!                                                                       
   45    IF(BIGA) 48,46,48 
   46    D=0.0D0 
         RETURN 
   48    DO 55 I=1,N 
         IF(I-K) 50,55,50 
   50    IK=NK+I 
         A(IK)=A(IK)/(-BIGA) 
   55    CONTINUE 
!                                                                       
! REDUCCION DE LA MATRIZ                                                
!                                                                       
         DO 65 I=1,N 
         IK=NK+I 
         HOLD=A(IK) 
         IJ=I-N 
         DO 65 J=1,N 
         IJ=IJ+N 
         IF(I-K) 60,65,60 
   60    IF(J-K) 62,65,62 
   62    KJ=IJ-I+K 
         A(IJ)=HOLD*A(KJ)+A(IJ) 
   65    CONTINUE 
!                                                                       
! DIVISION DE LAS FILAS POR EL PIVOTE                                   
!                                                                       
         KJ=K-N 
         DO 75 J=1,N 
         KJ=KJ+N 
         IF(J-K) 70,75,70 
   70    A(KJ)=A(KJ)/BIGA 
   75    CONTINUE 
!                                                                       
! PRODUCTO DE LOS PIVOTES Y SUBSTITUCION DEL PIVOTE POR SU INVERSO      
! CALCULO DEL DETERMINANTE                                              
!                                                                       
         D=D+DLOG10(DABS(BIGA)) 
         A(KK)=1.0D0/BIGA 
   80    CONTINUE 
!                                                                       
! INTERCAMBIO FINAL DE FILAS Y COLUMNAS                                 
!                                                                       
         K=N 
  100    K=K-1 
         IF(K) 150,150,105 
  105    I=L(K) 
         IF(I-K) 120,120,108 
  108    JQ=N*(K-1) 
         JR=N*(I-1) 
         DO 110 J=1,N 
         JK=JQ+J 
         HOLD=A(JK) 
         JI=JR+J 
         A(JK)=-A(JI) 
  110    A(JI)=HOLD 
  120    J=M(K) 
         IF(J-K) 100,100,125 
  125    KI=K-N 
         DO 130 I=1,N 
         KI=KI+N 
         HOLD=A(KI) 
         JI=KI-K+J 
         A(KI)=-A(JI) 
  130    A(JI)=HOLD 
         GO TO 100 
  150    RETURN 
         END SUBROUTINE M_MINV                                          
                                                                        
!****************************************************************       
      SUBROUTINE M_SPLINE(x,y,n,yp1,ypn,y2) 
!****************************************************************       
      implicit real*8 (a-h,o-z) 
      PARAMETER (NMAX=500) 
      dimension x(n),y(n),y2(n),u(NMAX) 
      if (yp1.gt..99e30) then 
        y2(1)=0. 
        u(1)=0. 
      else 
        y2(1)=-0.5 
        u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) 
      endif 
      do 11 i=2,n-1 
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1)) 
        p=sig*y2(i-1)+2. 
        y2(i)=(sig-1.)/p 
        u(i)=(6.*((y(i+1)-y(i))/(x(i+                                   &
     &1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*        &
     &u(i-1))/p                                                         
   11 continue 
      if (ypn.gt..99e30) then 
        qn=0. 
        un=0. 
      else 
        qn=0.5 
        un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) 
      endif 
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.) 
      do 12 k=n-1,1,-1 
        y2(k)=y2(k)*y2(k+1)+u(k) 
   12 continue 
      return 
      END SUBROUTINE M_SPLINE                                      
                                                                        
!****************************************************************       
      SUBROUTINE M_SPLINT(xa,ya,y2a,n,x,y) 
!****************************************************************       
      implicit real*8 (a-h,o-z) 
      dimension xa(n),y2a(n),ya(n) 
      klo=1 
      khi=n 
    1 if (khi-klo.gt.1) then 
        k=(khi+klo)/2 
        if(xa(k).gt.x)then 
          khi=k 
        else 
          klo=k 
        endif 
      goto 1 
      endif 
      h=xa(khi)-xa(klo) 
      if (h.eq.0.) pause 'bad xa input in m_splint' 
      a=(xa(khi)-x)/h 
      b=(x-xa(klo))/h 
      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**  &
     &2)/6.                                                             
      return 
      END SUBROUTINE M_SPLINT                                          
! Last modified 6/21/97                                                 
! ********************************************************************  
      SUBROUTINE MAC_INITIAL 
! ********************************************************************
      use mac_commons
      use gms 
      implicit real*8 (a-h,o-z) 
!     parameter(num=200,mum=99,nnum=15,kmum=99,lnum=99) 
      character(len=100), external:: dir                                               
      character*2 tnames(0:99) 
!     character*14 file_in 
      write (istdout,*) 'sub mac_init' 
! =================================      ! This subroutine contains all 
                                         ! open statements of the DIAG f
      call open_diag 
! =================================                                     
                          
! =================================                                     
      call tailnames(tnames) 
! =================================                                     
                                                                        
      do i=1,mum+1 
         igrounds(i) = 1000 
      enddo 
                                                                        
       cauxdr=clocation//'mac_init8.in'
       cauxdr=dir(workdir,cauxdr)
       open(unit=10,file=cauxdr) 
       read(10,*) 
        read(10,913)  element 
        read(10,909)  amass
        src_amass = amass
        read(10,905)  ieceecd 
        read(10,905)  ieciecr 
        read(10,905)  ipi 
        read(10,905)  irr 
        read(10,905)  idr 
        read(10,905)  isrd 
        read(10,905)  impi 
        read(10,*) 
        read(10,905)  iPress_Ion 
         write (istdout,*) 'ipress_ion',ipress_ion 
        read(10,906)  niter 
        read(10,907)  eps 
        read(10,907)  delta 
        read(10,907)  rLASER_I 
        read(10,907)  Ephoton 
        read(10,*) 
                                                                        
        read(10,904) nterms 
        read(10,*) 
       do iread = 1,nterms   ! itype -> 1=levels,2=terms,3=configurations
          read(10,914) ntermindex(iread),itype(iread) 
           write (istdout,*) 'ntermindex(',iread,')=',ntermindex(iread), &
     &       ' itype(',iread,')=',itype(iread)
       enddo 
        read(10,*) 
       write (istdout,*) 'nterms=',nterms 
        do iread = 1,nterms 
           fname_cat(iread)=clocation//"terms"//element//               &
     &    "."//tnames(iread)                                            
           write (istdout,*) 'fname_cat(',iread,')=',fname_cat(iread) 
       enddo 
                                                                        
       read(10,904) ngf 
       do iread = 1,ngf 
           fname_gf(iread)=clocation//"gf"//element//                   &
     &        "."//tnames(iread)                                        
         write (istdout,*) 'fname_gf(',iread,')=',fname_gf(iread) 
       enddo 
                                                                        
       read(10,904) nece 
       do iread = 1,nece 
          fname_ece(iread)=clocation//"ece"//element//                  &
     &        "."//tnames(iread)                                        
         write (istdout,*) 'fname_ece(',iread,')=',fname_ece(iread) 
       enddo 
                                                                        
       read(10,904) neci 
       read(10,*)
       do iread = 1,neci 
           read(10,'(22x,I1)') iecitype(iread)
           fname_eci(iread)=clocation//"eci"//element//                 &
     &    "."//tnames(iread)                                            
         write (istdout,*) 'fname_eci(',iread,')=',fname_eci(iread),' ',&
     &       iecitype(iread)
       enddo 

       read(10,*)                                                                 
       read(10,904) npi 
       do iread = 1,npi 
           fname_pi(iread)=clocation//"pi"//element//                   &
     &    "."//tnames(iread)                                            
          write (istdout,*) 'fname_pi(',iread,')=',fname_pi(iread) 
       enddo 
                                                                        
       read(10,*) 
                                                                        
       read(10,912) nrion 
       do iread = 1,nrion 
          read(10,909) rioneng(iread) 
          write (istdout,*) rioneng(iread) 
       enddo 
                                                                        
      close(10) 

      nstages = nterms 
      rIntensity=rLASER_I/(Ephoton*1.602d-19) 
      igrounds(0)=0 
       rioneng(0)=0 
                                                                        
! =========================================================             
        call read_levels 
        write (istdout,*) 'read_levels' 
! ---------------------------------------------                         
      call read_eceecd(ierror_ece) 
      write (istdout,*) 'read_eceed' 
      call read_eciecr(ierror_eci) 
      write (istdout,*) 'read_eciecr' 
      call read_pirr(ierror_pi) 
      write (istdout,*) 'read_pirr' 
      call m_gaulag(x,w,nnum,alf) 
      write (istdout,*) 'm_gaulag' 
      if((ierror_pi.eq.1).or.(ierror_eci.eq.1).or.(ierror_ece.eq.1))then 
          write (istdout,*) 'warn: improper cs detected' 
!          pause 
                                                                        
       endif 
! =========================================================             
                                                                        
      return 
  904 format(18x,I3) 
  905 format(18x,I1) 
  906 format(18x,I5) 
  907 format(18x,1pE8.2) 
  908 format(18x,a17) 
  909 format(18x,E14.7) 
  910 format(a3) 
!910   format(18x,a3)                                                   
  911 format(18x,a8) 
  912 format(18x,I3) 
  913 format(18x,a2) 
  914 format(18x,I3,I3) 
      END SUBROUTINE MAC_INITIAL                                        
!***********************************************************************
      SUBROUTINE MAC_LTE 
!***********************************************************************
      use mac_commons
      implicit real*8 (a-h,o-z) 
!     parameter(num=200,mum=99,nnum=15,kmum=99,lnum=99) 
!     open(unit=45,file= 'OUTPUT/mac_lte.out')                          
      ntotal = modeltot 
      do i=1,ntotal 
         f(i)=0.0d0 
      enddo 
                                                                        
      rkeVK=8.617385d-5 
                                                                        
!****************************************                               
!***** BEGINNING LEVEL OCCUPATION CALCS.                                
!****************************************                               
      DO iter = 1,niter 
         f(1)=1.0d0 
         j=2 
         sum=1.0d0 
         zbarp = zbar 
         rNe=zbar*rNa 
!---------------------------------                                      
         DO i=2,ntotal 
            IF(i.eq.igrounds(j)) THEN 
               aux1= rNe*2.070647635d-16*(Te**(-1.5)) 
               aux2= ((1.0d0*ig(igrounds(j-1)))/(1.0d0*ig(i))) 
                aux4=       rioneng(j-1)/(rkeVK*Te) 
                IF(aux4.gt.60.0d0)THEN 
                   f(i)=0.0d0 
                ELSE 
                  aux3=  dexp(aux4) 
                  aux=aux1*aux2*aux3 
                  f(i)=(f(igrounds(j-1)))/aux 
                ENDIF 
               j=j+1 
            ELSE 
               auxx1=f(i-1) 
               auxx2=((1.0d0*ig(i))/(1.0d0*ig(i-1))) 
                auxx4=(energy(i)-energy(i-1))/(rkeVK*Te) 
                IF(auxx4.gt.60.0d0)THEN 
                   f(i)=0 
                ELSE 
                  auxx3=dexp(-auxx4) 
                  f(i)= auxx1*auxx2*auxx3 
                ENDIF 
            ENDIF 
            sum=f(i)+sum 
         ENDDO 
!------------------------------------                                   
         DO i=1,ntotal 
            f(i)=f(i)/sum 
         ENDDO 
         zbar = 0.0d0 
         DO n = 1,ntotal 
            zbar = zbar + f(n) * zeff(n) 
         ENDDO 
!        write(45,900) (f(n), n = 1,ntotal)                             
         rNe=zbar*rNa 
!        write(45,*) 'zbar=',zbar,'  rNe=',rNe,' Te=',Te,'  iter=',iter 
!          write(45,*)                                                  
                                                                        
         if (iter.eq.1) goto 170 
         icheck = 0 
         do n = 1,ntotal 
            reld = dabs(f(n)-fp(n))/(fp(n)+1.0d-5) 
            if(reld.lt.eps) goto 160 
            icheck = icheck +1 
  160    enddo 
                                   ! goto 200 -> conv. criteria reached 
         if (icheck.eq.0) goto 200 
  170    do n = 1,ntotal 
            fp(n) = f(n) 
         enddo 
         zbar = (zbar + zbarp) * 0.5d0 
                                   ! 190 concludes the self interation  
      ENDDO 
                                   ! do loop                            
!************************************************************           
!******* END LEVEL OCCUPATION CALCS.*************************           
!************************************************************           
                                                                        
  200     rNe = zbar*rNa 
                                                                        
!********************** END OF ATOMIC KIN. KERNAL ********************* 
      return 
  900 format(10E9.2) 
      END SUBROUTINE MAC_LTE                                          
                                                                        
! **********************************************************************
      SUBROUTINE MAC_OUTPUT(iTempstep,iNa_step) 
! **********************************************************************
! This subroutine is called by mac_driver.f                             
! NOTE: THIS SUBROUTINE MUST RETAIN ITS VALUES: iswitch89,iswitch90                 
      use mac_commons
      implicit real*8 (a-h,o-z) 
!     parameter (num=200) 
      dimension rltelevels(num),ist(num) 
       character*60 char1 
       save iswitch89,iswitch90                                                      
       data iswitch89/0/
       data iswitch90/0/
      if(idiag(89).eq.0)goto 89 
          if(iswitch89.ne.89)then 
            iswitch89=89 
!           open(unit=10,file='mac_driver.in')   !need to be generalized
!  15      continue 
!          read(10,'(a60)',end=16) char1 
!           write(89,'(a60)') '## '//char1 
!           goto 15 
!  16      continue 
!          close(10) 
                                                                        

!          open(unit=10,file='mac_initial')       !need to be generalized
!  25      continue 
!           read(10,'(a60)',end=26) char1 
!          write(89,'(a63)') '## '//char1 
!           goto 25 
!  26      continue 
!          close(10) 
            write(89,'(a2,i5)')'# ',iTempstep +1 
            write(89,'(a2,i5)')'# ',iNa_step +1 
          endif 
         write(89,950) Te_eV,rNa,zbar 
   89 continue 
! --------------------------------------                                
      if(idiag(80).eq.0)goto 80 
         ftotal1 = 0.0d0 
         do ii = 1,28 
            ftotal1 = f(ii) + ftotal1 
         enddo 
         write(80,*) Te_eV, ftotal1 
   80 continue 
! --------------------------------------                                
      if(idiag(81).eq.0)goto 81 
      ftotal2 = 0.0d0 
      do ii = 29,39 
         ftotal2 = f(ii) + ftotal2 
      enddo 
      write(81,*) Te_eV, ftotal2 
   81 continue 
! --------------------------------------                                
      if(idiag(82).eq.0)goto 82 
      ftotal3 = 0.0d0 
      do ii = 40,51 
         ftotal3 = f(ii) + ftotal3 
      enddo 
      write(82,*) Te_eV, ftotal3 
   82 continue 
! --------------------------------------                                
      if(idiag(83).eq.0)goto 83 
      ftotal4 = 0.0d0 
      ftotal4 = f(52) 
      write(83,*) Te_eV, ftotal4 
   83 continue 
! --------------------------------------                                
                                                                        
!     sumf=0                                                            
!     do ii=1,modeltot                                                  
!        sumf=sumf+f(ii)                                                
!     enddo                                                             
                                                                        
      if(idiag(88).eq.0)goto 88 
      rewind(88)
      write(88,*)  ' rNe is only the contribution of this element' 
      write(88,'(a6,E15.8,a4,E15.8,a5,f3.1)') 'Te_eV=',Te_eV,' rNa',rNa,' eps=',eps 
      write(88,'(a5,E15.8,a4,E15.8,a6,i3)') 'zbar=',zbar,' rNe=',rNe,' iter=',iter 
      write(88,*) 'Globalindex |      Frac. Pop.     | Zeff |' 
      do ii = 1,modeltot 
         write(88,'(2x,i4,3x,E24.17,2x,f4.1)') ii,f(ii),zeff(ii) 
      enddo 
   88 continue 
                                                                        
      if(idiag(84).eq.0)goto 84 
      ftotal1 = 0.0d0 
      do ii = 1,28 
         ftotal1 = rltelevels(ii) + ftotal1 
      enddo 
      write(84,*) Te_eV, ftotal1 
   84 continue 
                                                                        
      if(idiag(85).eq.0)goto 85 
      ftotal2 = 0.0d0 
      do ii = 29,39 
         ftotal2 = rltelevels(ii) + ftotal2 
      enddo 
      write(85,*) Te_eV, ftotal2 
   85 continue 
                                                                        
      if(idiag(86).eq.0)goto 86 
      ftotal3 = 0.0d0 
      do ii = 40,51 
         ftotal3 = rltelevels(ii) + ftotal3 
      enddo 
      write(86,*) Te_eV, ftotal3 
   86 continue 
                                                                        
      if(idiag(87).eq.0)goto 87 
      ftotal4 = 0.0d0 
      ftotal4 = rltelevels(52) 
      write(87,*) Te_eV, ftotal4 
   87 continue 
                                                                        
      if(idiag(90).eq.0)goto 90 
          if(iswitch90.ne.90)then 
            iswitch90=90 
            open(unit=10,file='mac_driver.in')   !needs to be generalized
    5      continue 
           read(10,'(a60)',end=6) char1 
            write(90,'(a63)') '## '//char1 
            goto 5 
    6      continue 
           close(10) 
            write(90,'(a2,i5)')'# ',iTempstep +1 
            write(90,'(a2,i5)')'# ',iNa_step +1 
          endif 
         write(90,950) Te_eV,rNa,zbar 
   90 continue 
                                                                        
      if(idiag(91).eq.0)goto 91 
      write(91,*)  '                       ' 
      write(91,*) 'Te_eV=',Te_eV,' rNa',rNa,' eps=',eps 
      write(91,*) 'zbar=',zbar,' rNe=',rNe,' iter=',iter 
      do ii = 1,modeltot 
         write(91,*) '(',ii,')= ', rltelevels(ii) 
      enddo 
   91 continue 
                                                                        
! -------------------------------------------------------               
      if(idiag(34) .eq. 0) goto 9 
       kill=0 
      do i=1,modeltot 
         ist(i)=kill_lev(i) 
          if(ist(i).eq.0) kill=kill+1 
      enddo 
      write(34,*) 
      write(34,*) 'Na= ',rNa 
       write(34,*) '...... Eng_lower zbar TF_zbar Te_eV iconv #killed' 
      write(34,900)   ist(1),ist(2),ist(3),ist(4),ist(5),               &
     &               ist(6),ist(7),ist(8),ist(9),ist(10),               &
     &           ist(11),ist(12),ist(13),ist(14),ist(15),               &
     &           ist(16),ist(17),ist(18),ist(19),ist(20),               &
     &           ist(21),ist(22),ist(23),ist(24),ist(25),               &
     &           ist(26),ist(27),ist(28),ist(29),ist(30),               &
     &           ist(31),ist(32),ist(33),ist(34),ist(35),               &
     &           ist(36),ist(37),ist(38),ist(39),ist(40),               &
     &           ist(41),ist(42),ist(43),ist(44),ist(45),               &
     &           ist(46),ist(47),ist(48),ist(49),ist(50),               &
     &           ist(51),ist(52),Eng_lower,zbar,TF_zbar,Te_eV,iconv,kill
    9     continue 
! -------------------------------------------------------               
  900 format(28i1,1x,11i1,1x,12i1,1x,i1,1x,f5.2,1x,f5.2,1x,f5.2,        &
     &       1x,f5.2,1x,i1,1x,i3)                                       
                                                                        
  950 format(3(3x,1pE11.4)) 
  951 format(1x,1pE11.4,3x,1pE11.4) 
      END SUBROUTINE MAC_OUTPUT                                          
                                                                        
!***************************************************************        
       SUBROUTINE MATRIX_COMPRESSION(a,levels,nlev_rm,n,num) 
!***************************************************************        
! NOTE: Array levels must have the levels listed in increasing          
!       order                                                           
!                                                                       
       implicit real*8 (a-h,o-z) 
       dimension a(num,num),levels(num) 
                                                                        
       j=1 
   10  continue 
                                                                        
       if(levels(j).eq.n)then 
         n=n-1 
         return 
       endif 
                                                                        
         kk=levels(j) 
         do k=kk,n 
            do i=1,kk-1 
              a(k,i)=a(k+1,i) 
              a(i,k)=a(i,k+1) 
             enddo 
          enddo 
         do i=(kk+1),n 
            do k=(kk+1),n 
              a(k-1,i-1)=a(k,i) 
             enddo 
          enddo 
         n=n-1 
         j=j+1 
         levels(j) = levels(j)-1 
       if(j.le.nlev_rm)goto 10 
                                                                        
       return 
      END SUBROUTINE MATRIX_COMPRESSION                                          
!******************************************************************     
      SUBROUTINE MPI                            
!******************************************************************     
!* Planck's constant eV*s                                               
!* AlphaMPI statcoul/(g*(cm/s))                                         
!* Note:1eV=1.6d-12ergs                                                 
      use mac_commons                                                          
      implicit real*8 (a-h,o-z) 
!     parameter(num=200,mum=99,lnum=99) 
      dimension sigma(num,num) 
      pi         = DASIN(1.0d0)*2.0d0 
      alphampi   = 0.02648d0 
      plancks    = 4.14d-15 
                                                                        
      do i=1,num 
        do j=1,num 
           RMPI(i,j) = 0.0d0 
        enddo 
      enddo 
                                                                        
      hbar=plancks/(2.0d0*pi) 
      romega=Ephoton/hbar 
        do j=2,nstages-1 
          Eng_ioniz = rioneng(j-1)-Eng_lower 
          photon_num = Eng_ioniz/Ephoton 
          if((photon_num .gt.1).and.(photon_num .lt.4))then 
            nphotons= INT(Eng_ioniz/Ephoton)+1 
                                                                        
            sigma(igrounds(j-1),igrounds(j))=(alphampi**nphotons)       &
     &          *romega*(Eng_ioniz/Ephoton)**1.5                        &
     &          *(1.0d0/(romega**2*Eng_ioniz/Ephoton))                  &
     &          **nphotons                                              
                                                                        
                                                                        
            rmpi(igrounds(j-1),igrounds(j))=                            &
     &         sigma(igrounds(j-1),igrounds(j))                         &
     &         *rIntensity**nphotons                                    
          endif 
      enddo 
! -------------------------------------------------------------------   
! DIAGNOSITCS                                                           
      if(idiag(57).eq.0)goto 85 
      write(57,*) '                       ' 
      do l=1,modeltot 
         do k=1,modeltot 
          write(57,*)'rmpi=(',l,',',k,')=', rmpi(l,k) 
        enddo 
      enddo 
   85 continue 
      return 
      END SUBROUTINE MPI                                           

! *******************************************************               
      SUBROUTINE OUTPUT_OPEN 
! *******************************************************

      use gms
      use mac_commons
      implicit real*8 (a-h,o-z) 
      character*10 c_dir 
      character*17 c_file 
      character*27 c_place 
      character(len=100),external :: dir
      do i=1,100 
        idiag(i) = 0 
      enddo 
! ---------------------------------------                               
      c_place =clocation//'output.in'
      cauxdr=dir(workdir,c_place)
      open(unit=10,file=cauxdr) 
                       ! header                                         
      read(10,*) 
                       ! --- !                                          
      read(10,*) 
                       ! subheading                                     
      read(10,*) 
      read(10,900) idiag(8) 
      do i = 80,91 
         read(10,900) idiag(i) 
      enddo 
! ---------------------------------------                               
      c_dir = "OUTPUT."//clocation 
                                                                        
      if(idiag(8).eq.0)goto 8 
      c_file ="gnuinput.out" 
      c_place= c_dir//c_file 
      cauxdr=dir(workdir,c_place)
      open(unit=8,file=cauxdr) 
    8 continue 
                                                                        
      if(idiag(80).eq.0)goto 80 
      c_file ="f1.nlte" 
      c_place= c_dir//c_file 
      cauxdr=dir(workdir,c_place)
      open(unit=80,file=cauxdr) 
   80 continue 
                                                                        
      if(idiag(81).eq.0)goto 81 
      c_file ="f2.nlte" 
      c_place= c_dir//c_file 
      cauxdr=dir(workdir,c_place)
      open(unit=81,file=cauxdr) 
   81 continue 
                                                                        
      if(idiag(82).eq.0)goto 82 
      c_file ="f3.nlte" 
      c_place= c_dir//c_file 
      cauxdr=dir(workdir,c_place)
      open(unit=82,file=cauxdr) 
   82 continue 
                                                                        
      if(idiag(83).eq.0)goto 83 
      c_file ="f4.nlte" 
      c_place= c_dir//c_file 
      cauxdr=dir(workdir,c_place)
      open(unit=83,file=cauxdr) 
   83 continue 
                                                                        
      if(idiag(84).eq.0)goto 84 
      c_file ="f1.lte" 
      c_place= c_dir//c_file 
      cauxdr=dir(workdir,c_place)
      open(unit=84,file=cauxdr) 
   84 continue 
                                                                        
      if(idiag(85).eq.0)goto 85 
      c_file ="f2.lte" 
      c_place= c_dir//c_file 
      cauxdr=dir(workdir,c_place)
      open(unit=85,file=cauxdr) 
   85 continue 
                                                                        
      if(idiag(86).eq.0)goto 86 
      c_file ="f3.lte" 
      c_place= c_dir//c_file 
      cauxdr=dir(workdir,c_place)
      open(unit=86,file=cauxdr) 
   86 continue 
                                                                        
      if(idiag(87).eq.0)goto 87 
      c_file ="f4.lte" 
      c_place= c_dir//c_file 
      cauxdr=dir(workdir,c_place)
      open(unit=87,file=cauxdr) 
   87 continue 
                                                                        
      if(idiag(88).eq.0)goto 88 
      c_file = "complete.nlte" 
      c_place= c_dir//c_file 
      cauxdr=dir(workdir,c_place)
      open(unit=88,file=cauxdr) 
   88 continue 
                                                                        
      if(idiag(89).eq.0)goto 89 
      c_file ="zbar.nlte" 
      c_place= c_dir//c_file 
      cauxdr=dir(workdir,c_place)
      open(unit=89,file=cauxdr) 
   89 continue 
                                                                        
      if(idiag(90).eq.0)goto 90 
      c_file ="zbar.lte" 
      c_place= c_dir//c_file 
      cauxdr=dir(workdir,c_place)
      open(unit=90,file=cauxdr) 
   90 continue 
                                                                        
      if(idiag(91).eq.0)goto 91 
      c_file ="complete.lte" 
      c_place= c_dir//c_file 
      cauxdr=dir(workdir,c_place)
      open(unit=91,file=cauxdr) 
   91 continue 
! ---------------------------------------                               
  900 format(15x,i1) 
                                                                        
      END SUBROUTINE OUTPUT_OPEN
! --------------------------------------------------------              
! ===================== PARSER ===========================              
! --------------------------------------------------------              
! Delarations to call SUBROUTINE PARSER                                 
!     parameter(isize=5)                                                
!     character*40 config                                               
!     dimension iPrincQNs(isize),iOrbitalOccupant(isize)                
!     character*5 cOrbitalQNs                                           
! --------------------------------------------------------              
      SUBROUTINE PARSER(iPrincQNs,cOrbitalQNs,iOrbitalOccupant,config,k)
       implicit real*8 (a-h,o-z) 
      parameter(isize=5) 
       character*5 cOrbitalQNs 
       dimension iPrincQNs(isize),iOrbitalOccupant(isize) 
      character*45 config 
                                                                        
       istdout = 49
       iString = len(config) 
       cOrbitalQNs='    ' 
       do ii=1,isize 
          iPrincQNs(ii)=0 
          iOrbitalOccupant(ii)=0 
      enddo 
! -- Survey config array for orbital labels                             
      k=1 
       do i=1,iString 
              if(iFuncType(config(i:i)).eq. 2) then 
                 cOrbitalQNs(k:k)= config(i:i) 
             k=k+1 
          endif 
       enddo 
              ! total number of configuration terms                     
       k=k-1 
                                                                        
       i=0 
       j=1 
       jj=1 
       ival=1 
   10 continue 
                                                                        
       i=i+1 
       if(i.gt.iString) return 
       itype=iFuncType(config(i:i)) 
                                                                        
                                                    ! found a blank     
       if((itype .eq. 0).or.(itype .eq. 2))then 
         goto 10 
                                                    ! found numeral     
       elseif(itype .eq. 1)then 
              call alphaDigitsToInt(i,intr,config) 
              if(MOD(ival,2).eq.1)then 
             iPrincQNs(j)=intr 
                 j=j+1 
             ival=ival+1 
          else 
                     iOrbitalOccupant(jj)=intr 
                     jj=jj+1 
             ival=ival+1 
              endif 
              goto 10 
       else 
          write (istdout,*) 'PARSER ERROR: TYPE NOT FOUND' 
       endif 
       goto 10 
      END SUBROUTINE PARSER                                          
! ----------------------------------------------------------------------
! ========================= alphaDigitsToInt ===========================
! ----------------------------------------------------------------------
! Purpose:  Converts sequential numerical characters into an integer.   
! Remarks:  Must be employed when a number is detected at               
!           config(istart:istart). This subroutine then updates istart  
!                     where is finished in the array.                   
!           This routine is limited to a 30 digit integer               
! Date:              7/27/99                                            
! Author:       Manolo Sherrill                                         
! ----------------------------------------------------------------------
       SUBROUTINE ALPHADIGITSTOINT(istart,intr,config) 
       implicit real*8 (a-h,o-z) 
       dimension intarray(30) 
       character*45 config 
       do i=1,30 
         intarray(i) = 0 
       enddo 
       iString=len(config) 
                                                        ! end of config 
      if(istart.eq.iString)then 
              intr=icharToInt(config(istart:istart)) 
          return 
       endif 
       i=0 
                                                              ! numeral 
       do while((iFuncType(config(istart:istart)).eq. 1).and.           &
     &         (istart .le.iString))                                    
          i=i+1 
              intarray(i)=icharToInt(config(istart:istart)) 
              istart=istart+1 
       enddo 
       intr=0 
       do j=i,1,-1 
          intr=intr+intarray(j)*(10**(i-j)) 
       enddo 
                                                                        
      END SUBROUTINE ALPHADIGITSTOINT                                          
! --------------------------------------------------------              
! ===================== icharToInt ========================             
! --------------------------------------------------------              
      INTEGER FUNCTION ICHARTOINT(cChar) 
      use mac_commons
      implicit real*8 (a-h,o-z) 
       character*1 cChar 
                                                                        
      m=ichar(cChar) 
       if((m .lt. 48).and.(m .gt. 57)) then 
         write (istdout,*) 'icharToInt Error: Character non-numeric' 
       endif 
      icharToInt = m-48 
      END FUNCTION ICHARTOINT                                          
! --------------------------------------------------------              
       INTEGER FUNCTION IFUNCTYPE(cChar) 
! --------------------------------------------------------
      implicit real*8 (a-h,o-z) 
       character*1 cChar 
       istdout = 49                                                         
       iFuncType= -1 
       m=ichar(cChar) 
       if(m.eq.32)then 
              iFuncType=0 
                                                    ! numeral           
       elseif((m .ge. 48).and.(m .le. 57))   then 
          iFuncType=1 
                                                    ! upper case letter 
       elseif(((m .ge. 97).and.(m .le.122)).OR.                         &
     &       ((m .ge. 65).and.(m .le. 90))) then                        
                                                   ! lower case letter  
              iFuncType=2 
       else 
         write (istdout,*) 'iFuncType Error: digit does not fit any of the 3 cases'
       endif 
      END FUNCTION IFUNCTYPE                                          

!***********************************************************************
      SUBROUTINE PRESSURE_IONIZATION 
!***********************************************************************
      use mac_commons
      implicit real*8 (a-h,o-z) 
!     parameter (num=200) 
      do n=1,modeltot 
         kn=kill_lev(n) 
         do m=1,modeltot 
            knm=kn * kill_lev(m) 
            rsrd(n,m)   = rsrd(n,m) * knm 
            rmpi(n,m)   = rmpi(n,m) * knm 
            rdr(n,m)    = rdr(n,m)  * knm 
            rpi(n,m)    = rpi(n,m)  * knm 
            rrr(n,m)    = rrr(n,m)  * knm 
            receecd(n,m) = receecd(n,m) * knm 
            reciecr(n,m) = reciecr(n,m) * knm 
         enddo 
      enddo 
      END SUBROUTINE PRESSURE_IONIZATION
                                           
!***********************************************************************
      SUBROUTINE PRINTFNAME(rlte_levels) 
!***********************************************************************
      use mac_commons
      implicit real*8 (a-h,o-z) 
!     parameter(num=200,mum=99,lnum=99) 
      dimension                                                         &
     &          rlte_levels(num),                                       &
     &          energ_temp(num)                                         
                                                                        
      character*8 prefix 
      character*3 sufix 
      SAVE iset                                                                        
      if(iset.eq.1) goto 5 
                                                                        
       sufix  = ".te" 
       prefix = "OUTPUT/f" 
       call filenames(prefix,cfname,sufix) 
                                                                        
       sufix  = ".ne" 
       prefix = "OUTPUT/F" 
       call filenames(prefix,dfname,sufix) 
                                                                        
       sufix  = ".ps" 
       prefix = "OUTPUT/P" 
       call filenames(prefix,pfname,sufix) 
                                                                        
       iset=1 
    5 continue 
                                                                        
!** HARDWIRED                                                           
       do 10 ikl=1,modeltot 
                                                                        
         if(ikl.lt.igrounds(2)) then 
              energ_temp(ikl)=energy(ikl) 
         else 
            if((ikl.ge.igrounds(2)).and.(ikl.lt.igrounds(3))) then 
                energ_temp(ikl)=energy(ikl) + rioneng(1) 
            else 
                energ_temp(ikl)=energy(ikl)+rioneng(1)+rioneng(2) 
            endif 
         endif 
   10    continue 
!* Test                                                                 
!         open(unit=19,file=cfname(iTe))                                
!         open(unit=18,file=dfname(iTe))                                
                                                                        
!         do 40 j=1,modeltot                                            
!           write(19,951) energ_temp(j),rlte_levels(j)                  
!           write(18,951) energ_temp(j),f(j)                            
! 40        continue                                                    
                                                                        
       close(19) 
       close(18) 
                                                                        
  951  format(1x,1pE11.4,3x,1pE11.4) 
       return 
      END SUBROUTINE PRINTFNAME                                          
                                                                        
!***********************************************************************
      SUBROUTINE RATE_ECEECD 
!***********************************************************************
      use mac_commons
      use gms
      implicit real*8 (a-h,o-z) 
!     parameter (num=200) 
      do i=1,modeltot 
         do j=1,modeltot 
              receecd(i,j)=0.0d0 
              receecd(i,j)= rNe*coef_eced(i,j) 
         enddo 
      enddo 
! -----------------------------------------------------------           
!  DIAGNOSTICS                                                          
      if(idiag(51).eq.0) goto 120 
      do k=1,modeltot 
         do l=1,modeltot 
            write(51,900) 'receecd(',k,',',l,')=',receecd(k,l) 
        enddo 
      enddo 
      close(51)
  120 continue 
! ------------------------------------------------------------          
      return 
  900 format(A8,I3,A1,I3,A2,1pE14.7) 
      END SUBROUTINE RATE_ECEECD
                                                                        
!***********************************************************************
      SUBROUTINE RATE_ECIECR 
!***********************************************************************
      use mac_commons
      implicit real*8 (a-h,o-z) 
!     parameter(num=200) 
      do 10 i=1,modeltot 
         do 20 j=i,1,-1 
            reciecr(i,j)=0.0d0 
            reciecr(j,i)=0.0d0 
            reciecr(i,j)=(rNe**2)*coef_ecir(i,j) 
            reciecr(j,i)=rNe*coef_ecir(j,i) 
   20    continue 
   10 continue 
                                                                        
! ------------------------------------------------------                
! DIAGNOSTICS                                                           
      if(idiag(52).eq.0) goto 120 
                                                                        
      write(52,*) '                     ' 
      do k=1,modeltot 
         do l=1,modeltot 
            write(52,900) 'reciecr(',k,',',l,')=',reciecr(k,l) 
         enddo 
      enddo 
  120 continue 
                                                                        
      return 
  900 format(A8,I3,A1,I3,A2,1pE14.7) 
      END SUBROUTINE RATE_ECIECR                                          
                                                                        
                                                                        
                                                                        
!****************************************************************       
      SUBROUTINE RATE_RR 
!****************************************************************       
      use mac_commons
      implicit real*8 (a-h,o-z) 
!     parameter(num=200) 
      do i=1,modeltot 
        do j=1,modeltot 
           rrr(i,j)=0.0 
           rrr(i,j)=rNe*coef_rr(i,j) 
        enddo 
      enddo 
! ---------------------------------------------                         
! DIAGNOSTICS                                                           
                                                                        
      if(idiag(54).eq.0) goto 88 
      write(54,*) '                 ' 
      do k=1,modeltot 
         do l=1,modeltot 
            write(54,900) 'rrr(',k,',',l,')=',rrr(k,l) 
        enddo 
      enddo 
   88 continue 
      return 
  900 format(A4,I3,A1,I3,A2,1pE14.7) 
      END SUBROUTINE RATE_RR                                          
                                                                        
!*********************************************************************  
      SUBROUTINE RCOEF_ECEECD 
! Called from mac_kernal                                                
!*********************************************************************  
!  This subroutine calculates the electron impact excitation rate coef. 
!  when given the election temp. and electron density.  The electron    
!  energy distribution is assumed to be maxwellian.                     
!  last revision 6/26/97 8:16 am                                        
!*********************************************************************  
!* NOTE:  Te [Kelvin]                                                   
!* NOTE:  iglob1 < iglob2                                               
                                                                        
!*      i       Counts the ionization stages                            
!*      nptsece    Number of c.s data points for each c.s.                 
                                                                        
      use mac_commons
      use gms

      implicit real*8 (a-h,o-z) 
!     parameter (num=200,mum=99,lnum=99) 
!     parameter (npts=15,numcs=5000,numstages=8) 
!     parameter (nnum=15) 
      dimension                                                         &
     &          eng(nptsece),                                           &
     &          cs(nptsece),                                            &
     &          y2(nptsece),                                            &
     &          csint(nnum)                                             
      character*14 ecefn(0:999)
      character*8 prefix
      character*3 sufix
      data iii/0/
      save iii

      rkeVK  = 8.617385d-5 
      rkJK   = 1.380658d-23 
      remass = 9.1093897d-28 
      CONST  = 5.465383916d-11 
                                                                        
      do i=1,modeltot 
         do j=1,modeltot 
              coef_eced(i,j)=0.0d0 
         enddo 
      enddo 
!      write (istdout,*) 'rcoef_eceecd'                                            
                                                                        
!******** BEGIN 1ST DO WHILE LOOP (COUNTING OFF ION-STAGES) *********** 
      i=0 
   25 continue 
      i=i+1 
      if(i.eq.(nece+1))goto 110 
                                                                        
!******** BEGIN 2ND DO WHILE LOOP (COUNTING OFF RATES) **************** 
                                                                        
      j=0 
   30 continue 
       j=j+1 
      if(j.gt.jmax(i)) goto 100 
      ratio = ece_p_ratio(j,i)/Te 
      do mm=1,nptsece 
          eng(mm) = ece_eng(mm,j,i)*ratio 
                                                !  ece_cs = DLOG10(cs)  
          cs(mm)  = ece_cs(mm,j,i) 
      enddo 
      yp1 = ece_p_yp1(j,i)/ratio 
      ypn = ece_p_ypn(j,i)/ratio 
      call m_spline(eng,cs,nptsece,yp1,ypn,y2) 
                                                                        
! ------------------- Gausslag integration do while ------------------  
                                                                        
      sum=0.0d0 
      ii = 0 
   40 ii = ii + 1 
      if (ii.gt.nnum) goto 50 
      if (x(ii).lt.eng(nptsece)) then 
          xint      = x(ii) 
          call m_splint(eng,cs,y2,nptsece,xint,csint(ii)) 
          csinterpo = 10.0**(csint(ii)) 
          sum=sum+((1.0d0+(1.0d0/ratio)*x(ii))*csinterpo)*w(ii) 
      else 
         csint(ii) = ypn*(x(ii) - eng(nptsece)) + cs(nptsece) 
         csinterpo = 10.0**(csint(ii)) 
         sum=sum+((1.0d0+(1.0d0/ratio)*x(ii))*csinterpo)*w(ii) 
      endif 
      goto 40 
   50 continue 
                                                                        
      iglob1 = i_ece_iglob1(j,i) 
      iglob2 = i_ece_iglob2(j,i) 
                                                                        
! -------------------- Excitation Coef. ------------------------------- 
                                                                        
      coef_eced(iglob1,iglob2) = CONST*(Te**0.5)*                       &
     &          ratio*dexp(-ratio)*sum                                  
                                                                        
! -------------------- Dexcitation Coef. ------------------------------ 
                                                                        
! This format is NOT good! dexp(ratio) may lead to NaN                  
!     coef_eced(iglob2,iglob1) =                                        
!    *       coef_eced(iglob1,iglob2)*                                  
!    *       (1.0d0*ig(iglob1)*frac(iglob1))/(1.0d0*ig(iglob2)*         
!    *       frac(iglob2))*DEXP(ratio)                                  
                                                                        
                                                                        
      coef_eced(iglob2,iglob1) = CONST*(Te**0.5)*ratio*sum*             &
     &       (1.0d0*ig(iglob1)*frac(iglob1))/(1.0d0*ig(iglob2)*         &
     &       frac(iglob2))                                              
                                                                        
      goto 30 
                                                                        
!********************  END OF 2ND DO WHILE LOOP (j) ******************* 
                                                                        
  100 continue 
      goto 25 
  110 continue 
                                                                        
!********************  END OF 1ST DO WHILE LOOP (i) ******************* 
                                                                        
! --------------------------- DIAGNOSTICS ----------------------------- 
!     if(iconvrg .eq. 1) then
!      sufix  = ".rt"
!      prefix = "recexxxx"
!      call filenames(prefix,ecefn,sufix)

!     iii=iii+1
!     open(unit=61,file=ecefn(iii))
!     write(61,*)'            '
!     write(61,*) 'Na=',rna
!     write(61,*) 'Te=',te_eV
                                                                        
      if(idiag(61).eq.0) goto 120 
      do k=1,modeltot 
         do l=1,modeltot 
            write(61,900) 'coef_eced(',k,',',l,')=',coef_eced(k,l) 
         enddo 
      enddo 
!     endif
  120 continue 
                                                                        
      return 
  900 format(A10,I3,A1,I3,A2,1pE14.7) 
  915 format(12X,A4,23X,A4,2x,i7,i7,1X,E14.7) 
      END SUBROUTINE RCOEF_ECEECD
                                                                        
!*********************************************************************  
      SUBROUTINE RCOEF_ECIECR 
! Called from mac_kernal                                                
!*********************************************************************  
!  This subroutine calculates the electron impact ionization rate coef. 
!  when given the election temp. and electron density.  The electron    
!  energy distribution is assumed to be maxwellian.                     
!  last revision 6/26/97 8:16 am                                        
!*********************************************************************  
!* NOTE:  Te [Kelvin]                                                   
!* NOTE:  iglob1 < iglob2                                               
                                                                        
!*      i       Counts the ionization stages                            
!*      nptseci    Number of c.s data points for each c.s.                 
                                                                        
      use mac_commons
      implicit real*8 (a-h,o-z) 
!     parameter (num=200,mum=99,lnum=99) 
!     parameter (npts=21,numcs=5000,numstages=8)
!     parameter (nnum=15) 
      dimension                                                         &
     &          eng(nptseci),                                           &
     &          cs(nptseci),                                            &
     &          y2(nptseci),                                            &
     &          csint(nnum)

! -- Constants ---------------------------------------------------------
      rkeVK   = 8.617385d-5 
      rkergsK = 1.380658d-16 
      CONST  = 5.465383916d-11 
                                                                        
      do i=1,num 
         do j=1,num 
            coef_ecir(i,j)=0.0d0 
         enddo 
      enddo 
!       write (istdout,*) 'rcoef_eci.f'                                            
!******** BEGIN 1ST DO WHILE LOOP (COUNTING OFF ION-STAGES) *********** 
                                                                        
      i=0 
   25 continue 
      i=i+1 
      if(i.eq.(neci+1))goto 110 
!******** BEGIN 2ND DO WHILE LOOP (COUNTING OFF RATES) **************** 
!         *          *          *          *          *          *12    
      j=0 
   30 continue 
      j=j+1 
      if(j.gt.jmax_eci(i)) goto 100 
      engdif = (glob_eng(i_eci_iglob2(j,i))                             &
     &         -glob_eng(i_eci_iglob1(j,i)))-Eng_lower                  
                                                                        
      if(engdif.le.0.0d0)then 
         coef_ecir(i_eci_iglob1(j,i),i_eci_iglob2(j,i))=0.0d0 
         coef_ecir(i_eci_iglob2(j,i),i_eci_iglob1(j,i))=0.0d0 
         goto 30 
      endif 
                                                                        
      ratio = engdif/(rkeVK*Te) 
      do mm=1,nptseci 
          eng(mm) = (eci_eng(mm,j,i)-1.0d0)*ratio 
          cs(mm)  = DLOG10(eci_cs(mm,j,i)) 
      enddo 
      yp1 = ((cs(2)-cs(1))/(eng(2)-eng(1))) 
      ypn = ((cs(nptseci)-cs(nptseci-1))/                               &
     &      (eng(nptseci)-eng(nptseci-1)))                                    
! =========================================                             
      call m_spline(eng,cs,nptseci,yp1,ypn,y2) 
! =========================================                             
                                                                        
! ------------------- Gausslag integration do while ------------------  
                                                                        
      sum=0.0d0 
      ii = 0 
   40 ii = ii + 1 
      if (ii.gt.nnum) goto 50 
      if (x(ii).lt.eng(nptseci)) then 
         xint = x(ii) 
! ====================================================                  
         call m_splint(eng,cs,y2,nptseci,xint,csint(ii)) 
! ====================================================                  
         csinterpo = 10.0**(csint(ii)) 
         sum = sum+((1.0d0+(1.0d0/ratio)*x(ii))*csinterpo)*w(ii) 
      else 
         csint(ii) = ypn*(x(ii) - eng(nptseci)) + cs(nptseci) 
         csinterpo = 10.0**(csint(ii)) 
         sum=sum+((1.0d0+(1.0d0/ratio)*x(ii))*csinterpo)*w(ii) 
      endif 
      goto 40 
   50 continue 
                                                                        
      iglob1 = i_eci_iglob1(j,i) 
      iglob2 = i_eci_iglob2(j,i) 
                                                                        
       if(ypn.gt.0.0d0)then 
        ypn = ((cs(nptseci)-cs(nptseci-2))/(eng(nptseci)-eng(nptseci-2))) 
       endif 
         if(ypn.gt.0.0d0)then 
           write (istdout,*) 'eci_ypn(',iglob1,',',iglob2,')=',ypn,'< 0.0' 
           stop 'RCOEF_ECIECR' 
         endif 
                                                                        
! -------------------- Ionization Coef. ------------------------------- 
                                                                        
        coef_ecir(iglob1,iglob2)=CONST*(Te**0.5)*ratio*                 &
     &                      dexp(-ratio)*sum                            
                                                                        
! -------------------- Recombination Coef. ---------------------------- 
                                                                        
      C1 = 2.0706325d-16 
                                                                        
! Not a good format dexp(ratio) may become uncontroleable               
                                                                        
!       coef_ecir(iglob2,iglob1) = coef_ecir(iglob1,iglob2)*            
!    *     C1*(Te**(-1.5))*(1.0d0*ig(iglob1)*frac(iglob1)               
!    *      /(1.0d0*ig(iglob2)*frac(iglob2)))*                          
!    *     DEXP(ratio)                                                  
                                                                        
        coef_ecir(iglob2,iglob1) =CONST*(Te**0.5)*ratio*sum*            &
     &       C1*(Te**(-1.5))*(1.0d0*ig(iglob1)*frac(iglob1)             &
     &      /(1.0d0*ig(iglob2)*frac(iglob2)))                           
      goto 30 
                                                                        
!********************  END OF 2ND DO WHILE LOOP (j) ******************* 
                                                                        
  100 continue 
      goto 25 
  110 continue 
                                                                        
!********************  END OF 1ST DO WHILE LOOP (i) ******************* 
                                                                        
! --------------------------- DIAGNOSTICS ----------------------------- 
                                                                        
      if(idiag(62).eq.0) goto 120 
      do k=1,modeltot 
         do l=1,modeltot 
            write(62,900) 'coef_ecir(',k,',',l,')=',coef_ecir(k,l) 
         enddo 
      enddo 
                                                                        
  120 continue 
! --------------------------------------------------------------------- 
      return 
  900 format(A10,I3,A1,I3,A2,1pE14.7) 
  915 format(12X,A4,23X,A4,2x,i7,i7,1X,E14.7) 
      END SUBROUTINE RCOEF_ECIECR
                                          
!*********************************************************************  
      SUBROUTINE RCOEF_PIRR 
! Called from mac_kernal                                                
!*********************************************************************  
!  This subroutine calculates the photoionization rate and radiative    
!  recombination rate coef. when given the election temp. and electron  
!  density.  The electron energy distribution is assumed to be          
!  maxwellian. Note: Since photoionization rate needs on the Te and not 
!  the electron density we can calculate it here to prevent redundent   
!  calculations in the self consistency loop of mac_kernal.             
!  last revision 7/9/97 12:36 pm                                        
!*********************************************************************  
!* NOTE:  Te [Kelvin]                                                   
!* NOTE:  iglob1 < iglob2                                               
                                                                        
!*      g_i       Counts the ionization stages                          
!*      g_npts    Number of c.s data points for each c.s.               
                                                                        
      use mac_commons

      implicit real*8 (a-h,o-z) 
!     parameter (num=200,mum=99,lnum=99) 
!     parameter (npts=21,numcs=5000,numstages=8) 
!     parameter (nnum=15) 
      dimension                                                         &
     &          eng(nptsrr),                                            &
     &          cs(nptsrr),                                             &
     &          y2(nptsrr),                                             &
     &          csint(nnum)                                             
! -- Constants ---------------------------------------------------------
      alf = 0.0d0 
      rkeVK=8.617385d-5 
      remass = 9.1093897d-31 
                                                                        
                                                                        
      iikk = 0 
      do i=1,modeltot 
        do j=1,modeltot 
           coef_rr(i,j)=0.0 
           rpi(i,j)=0.0 
        enddo 
      enddo 
!******** BEGIN 1ST DO WHILE LOOP (i) (COUNTING OFF ION-STAGES) ******* 
      i=0 
   25 continue 
      i=i+1 
      if(i.eq.(npi+1))goto 110 
                                                                        
!*********** BEGIN 2ND DO WHILE LOOP (j)(COUNTING OFF RATES) ********** 
                                                                        
      j=0 
   30 continue 
      j=j+1 
      if(j.gt.jmax_rr(i)) goto 100 
      engdif = glob_eng(i_rr_iglob2(j,i))-                              &
     &         glob_eng(i_rr_iglob1(j,i))                               &
     &         - Eng_lower                                              
                                                                        
      if(engdif .lt. 1.0d-5) then 
        coef_rr(i_rr_iglob2(j,i),i_rr_iglob1(j,i))=0.0d0 
        rpi(i_rr_iglob1(j,i),i_rr_iglob2(j,i)) = 0.0d0 
        goto 30 
      endif 
                                                                        
      Te_eV = rkeVK*Te 
                                                                        
      ratio =  engdif/Te_eV 
      XEphoton   = (Ephoton/engdif) 
      ergsengdif = engdif*1.600217733d-12 
      do mm=1,nptsrr 
        eng(mm) = (rr_eng(mm,j,i)-1.0d0)*ratio 
         cs(mm)  = DLOG10(rr_cs(mm,j,i)) 
      enddo 
                                                                        
      yp1 = ((cs(2)-cs(1))/(eng(2)-eng(1))) 
      ypn = ((cs(nptsrr)-cs(nptsrr-1))/(eng(nptsrr)-eng(nptsrr-1))) 
! =========================================                             
      call m_spline(eng,cs,nptsrr,yp1,ypn,y2) 
! =========================================                             
! ------------------- Gausslag integration do while ------------------  
                                                                        
      sum=0.0d0 
      ii = 0 
   40 ii = ii + 1 
      if (ii.gt.nnum) goto 50 
      if (x(ii).lt.eng(nptsrr)) then 
          xint = x(ii) 
! ====================================================                  
          call m_splint(eng,cs,y2,nptsrr,xint,csint(ii)) 
! ====================================================                  
          csinterpo = 10.0**(csint(ii)) 
          sum=sum+((1.0d0+(1.0d0/ratio)*x(ii))**2*csinterpo)*w(ii) 
          interpo = interpo + 1 
      else 
          csint(ii) = ypn*(x(ii) - eng(nptsrr)) + cs(nptsrr) 
          csinterpo = 10.0**(csint(ii)) 
          sum=sum+((1.0d0+(1.0d0/ratio)*x(ii))**2*csinterpo)*w(ii) 
      endif 
      goto 40 
   50 continue 
                                                                        
                                                                        
      iglob1 = i_rr_iglob1(j,i) 
      iglob2 = i_rr_iglob2(j,i) 
                                                                        
       if(ypn.gt.0.0d0)then 
        ypn = ((cs(nptsrr)-cs(nptsrr-2))/(eng(nptsrr)-eng(nptsrr-2))) 
       endif 
       if(ypn.gt.0.0d0)then 
         write (istdout,*) 'ion=',i,'#=',j,'rr_ypn(',iglob1,',',iglob2,')=',ypn,'< 0.0' 
         stop 'RCOEF_PIRR'
       endif 
       if(XEphoton .le. 1.0d0) goto 60 
          iikk = iikk + 1 
          zEphoton  = (XEphoton-1.0d0)*ratio 
!           write(99,*)'iglob1=',iglob1,'  ',iglob2                     
!           write(99,*)'y2=',y2                                         
!           write(99,*)                                                 
!           write(99,*)'cs=',cs                                         
!       write(99,*)'eng=',eng                                           
                                                                        
! =====================================================                 
                                                               !Eval. pi
          call m_splint(eng,cs,y2,nptsrr,zEphoton,picsint) 
! =====================================================                 
          pics = 10.0**(picsint) 
!---- Photoionization rate ----!                                        
       rpi(iglob1,iglob2) = rIntensity*8.797356697D-17*(pics) 
   60  continue 
! -------------------- Recombination Coef. ---------------------------- 
      aux = 5.537169448d25*(((ergsengdif)**3)/ratio) 
      coef_rr(iglob2,iglob1) = aux * (Te**(-1.5))*                      &
     &                         ((1.0d0*ig(iglob1)*frac(iglob1))/        &
     &                         (1.0d0*ig(iglob2)*frac(iglob2)))*        &
     &                         sum                                      
                                                                        
      goto 30 
                                                                        
!********************  END OF 2ND DO WHILE LOOP (j) ******************* 
                                                                        
  100 continue 
      goto 25 
  110 continue 
                                                                        
!********************  END OF 1ST DO WHILE LOOP (i) ******************* 
                                                                        
! -----------------------------------------------------------           
! DIAGNOSTICS                                                           
                                                                        
      if(idiag(53).eq.0) goto 130 
      do k=1,modeltot 
         do l=1,modeltot 
            write(53,900) 'rpi(',k,',',l,')=',rpi(k,l) 
         enddo 
      enddo 
  130 continue 
! ------------------------------------------------------------          
      if(idiag(64).eq.0) goto 120 
      do k=1,modeltot 
         do l=1,modeltot 
            write(64,900) 'coef_rr(',k,',',l,')=',coef_rr(k,l) 
         enddo 
      enddo 
      close(28) 
  120 continue 
! ------------------------------------------------------------          

      return 
  900 format(A8,I3,A1,I3,A2,1pE14.7) 
  915 format(12X,A4,23X,A4,2x,i7,i7,1X,E14.7) 
      END  SUBROUTINE RCOEF_PIRR
                                          
!***********************************************************************
      SUBROUTINE READ_ECEECD(ierror_cs) 
! called by mac_initial                                                 
!***********************************************************************
!  This subroutine reads ACE (electron impact excitation crossectional  
!  data) and loads arrays to be used during the self consistent protion 
!  of MAC.                                                              
!  last revision 6/23/97 3:56pm                                         
!***********************************************************************
!* NOTE:  iglob1 < iglob2                                               
                                                                        
!*      g_i       Counts the ionization stages.                         
!*      g_j       Counts the number of cs in each ionization stage.     
!*      g_npts    Number of c.s data points for each c.s.               
                                                                        
!*********************************************************************  
      use mac_commons                                            
      use gms
      implicit real*8 (a-h,o-z) 
!     parameter (num=200,mum=99,lnum=99,kmum=99) 
!     parameter (npts=15,numcs=5000,numstages=8) 
      dimension eng(nptsece),cs(nptsece) 
      dimension tmpcs(2) 
      character(len=100),external :: dir
! -- Constants ---------------------------------------------------------
                                                                        
      rkeVK  = 8.617385d-5 
      rkJK   = 1.380658d-23 
      remass = 9.1093897d-28 
      CONST  = 5.465383916d-11 
                                                                        
!********** BEGIN 1ST DO WHILE LOOP (COUNTING OFF ION-STAGES) I ******* 
       ierror_cs=0 
      i=0 
   25 continue 
                                 ! i counts off ionizations.            
      i=i+1 
      if(i.eq.(nece+1))goto 110 
       cauxdr=fname_ece(i)
       cauxdr=dir(workdir,cauxdr)

      open(unit=20,file=cauxdr) 
                                                                        
                                                                        
                                                                        
!*********** BEGIN 2ND DO WHILE LOOP (COUNTING OFF RATES) J *********** 
                                                                        
                                        !counts off number of cs per i. 
      j=0 
   30 continue 
      read(20,915,END=100) gridxx,engxx,icat1,icat2,deltaeng 
      read(20,*) 

      i1 = iquery(icat1,i)
      i2 = iquery(icat2,i)

!      if((i1.le.ntermindex(i)).and.(i2.le.ntermindex(i))) then 
       if((i1.gt.0).and.(i2.gt.0))then 
          j=j+1 
          jmax(i)= j 
          iglob1 = (igrounds(i)-1)+i1 
          iglob2 = (igrounds(i)-1)+i2 
                                                                        
          i_ece_iglob1(j,i) = iglob1 
          i_ece_iglob2(j,i) = iglob2 
                                                                        
          engdif = glob_eng(iglob2) - glob_eng(iglob1) 
         if(engdif.lt.0.0d0)then 
            write (istdout,*) 'warn: readece.f:energy diff.<0' 
         endif 
          ece_p_ratio(j,i) = (engdif/(rkeVK)) 
                                                                        
           iswitch=0 
           ier=0 
                                   ! mm counts off the number of cs pts.
          do mm=1,nptsece 
             read(20,*,END=130) eng(mm),cs(mm) 
              if (cs(mm).le.0.0d0) cs(mm)=dabs(cs(mm)) 
              ccs=cs(mm) 
              if(mm.eq.1)then 
                       tmpcs(2)=ccs 
              else 
                 tmpcs(1)=tmpcs(2) 
                 tmpcs(2)=ccs 
                 if(iswitch.eq.0)then 
                                                               !check: g
                    if(tmpcs(2).lt.tmpcs(1)) iswitch=1 
                 endif 
                    if((iswitch.eq.1).and.(tmpcs(2).gt.tmpcs(1)).and.   &
     &               (ier.eq.0))then                                    
                   ier=1 
                   ierror_cs=1 
                   diflog=DABS(DLOG10(tmpcs(2))-DLOG10(tmpcs(1))) 
                      write(40,902) 'warn: irreg. cs detected in ',     &
     &            fname_ece(i),i1,i2,diflog                             
                endif 
               endif 
              cs(mm)  = DLOG10(cs(mm)) 
              eng(mm) = (eng(mm)-1.0d0) 
              ece_cs(mm,j,i)  = cs(mm) 
              ece_eng(mm,j,i) = eng(mm) 
          enddo 
                                                  ! -- Slope @ pts. 1 & 
                                                                        
       ece_p_yp1(j,i) = ((cs(2)-cs(1))/(eng(2)-eng(1))) 
       ece_p_ypn(j,i) = ((cs(nptsece)-cs(nptsece-1))/(eng(nptsece)-eng(nptsece-1))) 
                                                                        
       if(ece_p_ypn(j,i).gt.0.0d0)then 
        ece_p_ypn(j,i) = ((cs(nptsece)-cs(nptsece-2))/(eng(nptsece)-eng(nptsece-2))) 
       endif 
       if(ece_p_ypn(j,i).gt.0.0d0)then 
         write (istdout,*) 'ion=',i,'icat1=',icat1,' icat2=',icat2, &
     &     ' ece_p_ypn(',j,',',i,')=',ece_p_ypn(j,i) 
         stop 'ECE DATA ERROR' 
       endif 
                                                                        
          goto 30 
      else 
          do mm=1,nptsece 
             read(20,*,END=140) eng(mm),cs(mm) 
          enddo 
          goto 30 
      endif 
                                                                        
!*********************  END OF 2ND DO WHILE LOOP ********************** 
       goto 100 
                                                                        
  130 write (istdout,*) 'read halted prematurely',mm,'130 indexies satisfy criter.' 
      stop 'READ_ECEECD:read halted prematurely mm 130 indexies'
  140 write (istdout,*) 'read halted prematurely',mm,'140 indexies not satisfy cr.' 
      stop 'READ_ECEECD:read halted prematurely mm 140 indexies' 
  100 continue 
      close(20) 
       write (istdout,*) '# CS-ECE(',i,')=',jmax(i) 
      goto 25 
  110 continue 
!*********************  END OF 1ST DO WHILE LOOP ********************** 
                                                                        
! --------------------------------------------------------------------- 
! DIAGNOSTICS                                                           
      if(idiag(41).eq.0) goto 120 
                                                                        
      write(41,*) '                        ' 
      write(41,*) 'jmax=',jmax 
      do i = 1,nece 
        do j=1,jmax(i) 
          write(41,901) 'ece_p_yp1(',j,',',i,')=',ece_p_yp1(j,i) 
          write(41,901) 'ece_p_ypn(',j,',',i,')=',ece_p_ypn(j,i) 
          write(41,903) 'ece_p_ratio(',j,',',i,')=',ece_p_ratio(j,i) 
          write(41,904) 'i_ece_iglob1(',j,',',i,')=',i_ece_iglob1(j,i) 
          write(41,904) 'i_ece_iglob2(',j,',',i,')=',i_ece_iglob2(j,i) 
          do mm=1,nptsece
            write(41,900) 'mm=',mm,ece_eng(mm,j,i),(10**ece_cs(mm,j,i)) 
          enddo 
        enddo 
      enddo 
                                                                        
  900 format(a3,i2,2x,1pE14.7,2x,1pE14.7) 
  901 format(a10,i3,a1,i1,a2,1pE14.7) 
  902     format(a29,a17,i3,3x,i3,4x,1pE11.4) 
  903 format(a12,i3,a1,i1,a2,1pE14.7) 
  904 format(a13,i3,a1,i1,a2,i3) 
  120 continue 
                                                                        
      return 
  915 format(12X,A4,23X,A4,2x,i7,i7,1X,E14.7) 
      END SUBROUTINE READ_ECEECD                                          
!***********************************************************************
      SUBROUTINE QUERYALL(ival,ionix,ityp,ivec,kk)
!*********************************************************************** 
! This subroutine determines the association between a less detailed
! representation and the (natural) rep. chosen for a particular ionization 
! stage. Local ordered index
! Remark: This subroutine assumes at least one index association exists.
!
      use mac_commons
      implicit real*8 (a-h,o-z)
      dimension ivec(900)
      do i=1,900
         ivec(i)=0
      enddo
      kk=0
      do ii=1,ntermindex(ionix)
         if(ival.eq.ioracle(ii,ionix,ityp)) then  
           kk=kk+1
           ivec(kk)=ii
         endif
      enddo 
      if(kk.eq.0) then 
         print*,'ERROR: QUERYALL',ival,' association not found'
         stop 'HALT: QUERYALL'
      endif  
        
      END SUBROUTINE QUERYALL

!***********************************************************************
      INTEGER FUNCTION ICHK(ival,ionix,ityp)
!***********************************************************************
! If association is not found for ival
! then ichk=-1
      use mac_commons
      implicit real*8 (a-h,o-z)
      ichk = -1
      do ii=1,ntermindex(ionix)
       io=ioracle(ii,ionix,ityp)
       if(io.eq.ival) then
         ichk=ii
         return
       endif
      enddo

      END FUNCTION ICHK

!***********************************************************************
      SUBROUTINE READ_ECIECR(ierror_cs) 
! called by mac_initial                                                 
!***********************************************************************
!  This subroutine reads ACE (electron impact ionization crossectional  
!  data) and loads arrays to be used during the self consistent protion 
!  of MAC.                                                              
!  last revision 6/30/97 9:51 pm                                        
!***********************************************************************
!* NOTE:  iglob1 < iglob2                                               
                                                                        
!*      g_i       Counts the ionization stages.                         
!*      g_j       Counts the number of cs in each ionization stage.     
!*      g_npts    Number of c.s data points for each c.s.               
                                                                        
!*********************************************************************  
      use mac_commons                                 
      use gms
      implicit real*8 (a-h,o-z) 
!     parameter (num=200,mum=99,lnum=99,kmum=99) 
!     parameter (npts=21,numcs=5000,numstages=8) 
      integer ivec1(900),ivec2(900),iv1(9000),iv2(9000)
       dimension tmpcs(2),v_cs(nptseci),v_eng(nptseci) 
      character(len=100),external :: dir
! -- Constants ---------------------------------------------------------
      rkeVK   = 8.617385d-5 
      rkergsK = 1.380658d-16 
      CONST  = 5.465383916d-11 
                                                                        
!*********** BEGIN 1ST DO WHILE LOOP (COUNTING OFF ION-STAGES)          
       ierror_cs=0 
      i=0 
   25 continue 
                         ! i counts off the lower ionization stage.     
      i=i+1 
      if(i.eq.(neci+1))goto 110 
      cauxdr=fname_eci(i)
      cauxdr=dir(workdir,cauxdr)
      open(unit=20,file=cauxdr) 
      read(20,*) 
      read(20,*) 
      read(20,*) 
                                                                        
!************ BEGIN 2ND DO WHILE LOOP (COUNTING OFF RATES)              
                                                                        
      j=0 
   30 continue 
      read(20,915,END=100) gridxx,engxx,icat1,icat2,deltaeng 
      i1aux=ichk(icat1,i,iecitype(i))
      i2aux=ichk(icat2,i+1,iecitype(i))                                                                    
      
      if((i1aux.gt.0).and.(i2aux.gt.0))then   ! Check indexies if the exist in ion. stage 
 
! -- Mapping ionization c.s. indexies

        if((iecitype(i).eq.itype(i)).and.(iecitype(i).eq.itype(i+1)))then
           iv1(1)=i1aux
           iv2(1)=i2aux
           ktot=1
        elseif((iecitype(i).ne.itype(i)).and.(iecitype(i).eq.itype(i+1)))then
           call queryall(icat1,i,iecitype(i),ivec1,k1)
           ktot=k1
           do ijj=1,k1
              iv1(ijj)=ivec1(ijj)
              iv2(ijj)=i2aux
           enddo 
        elseif((iecitype(i).eq.itype(i)).and.(iecitype(i).ne.itype(i+1)))then
           print*,'READ_ECI: upper index mappings req'
           call queryall(icat2,i+1,iecitype(i),ivec2,k2)
           ktot=k2
           do ijj=1,k2
              iv1(ijj)=i1aux
              iv2(ijj)=ivec2(ijj)
           enddo
        elseif((iecitype(i).ne.itype(i)).and.(iecitype(i).ne.itype(i+1)))then
           print*,'READ_ECI: 2 mappings req'
           call queryall(icat1,i,iecitype(i),ivec1,k1)
           call queryall(icat2,i+1,iecitype(i),ivec2,k2)
           ktot=k1*k2
           ict=0
           do ij1=1,k1
           do ij2=1,k2
              ict=ict+1
              iv1(ict)=ivec1(ij1)
              iv2(ict)=ivec2(ij2)
           enddo
           enddo
        else
           stop 'ERROR: READ_ECI INDEX CONDITIONAL TRAP'
        endif 

! -- Check c.s.

          ier=0 
          iswitch=0 
                                  ! mm counts off the number of cs pts. 
         do mm=1,nptseci 
            read(20,*) eng,cs 
             if (cs.le.0.0d0) cs=dabs(cs) 
             if(mm.eq.1)then 
                     tmpcs(2)=cs 
             else 
                   tmpcs(1)=tmpcs(2) 
               tmpcs(2)=cs 
                                                                        
               if(iswitch.eq.0)then 
                                                               !check: g
                 if(tmpcs(2).lt.tmpcs(1)) iswitch=1 
               endif 
                  if((iswitch.eq.1).and.(tmpcs(2).gt.tmpcs(1)).and.     &
     &             (ier.eq.0))then                                      
                 ier=1 
                 ierror_cs=1 
                 diflog=DABS(DLOG10(tmpcs(2))-DLOG10(tmpcs(1))) 
                  write(40,901) 'warn: irreg. cs detected in ',fname_eci(i),&
     &                   icat1,icat2,diflog                                   
                endif 
              endif 
!           eci_cs(mm,j,i)  = cs 
!           eci_eng(mm,j,i) = eng 
            v_cs(mm) = cs
            v_eng(mm)= eng
         enddo  !mm  loop 

        do ijj=1,ktot
          i1=iv1(ijj)
          i2=iv2(ijj)
          j=j+1
          jmax_eci(i)=j
          iglob1 = (igrounds(i)-1)+i1
          iglob2 = (igrounds(i+1)-1)+i2
          i_eci_iglob1(j,i)=iglob1
          i_eci_iglob2(j,i)=iglob2
          do mm=1,nptseci
             eci_cs(mm,j,i) = v_cs(mm)
             eci_eng(mm,j,i)= v_eng(mm)
          enddo
         enddo  !ijj loop

         goto 30 
      else 
         do mm=1,nptseci 
            read(20,*) eng,cs 
         enddo 
         goto 30 
      endif 
                                                                        
!*********************  END OF 2ND DO WHILE LOOP ********************** 
                                                                        
  100 continue 
      close(20) 
       write (istdout,*) '# CS-ECI(',i,')=',jmax_eci(i) 
      goto 25 
  110 continue 
                                                                        
!*********************  END OF 1ST DO WHILE LOOP ********************** 
! ----------------------------------------------------------------------
! DIAGNOSTICS                                                           
      if(idiag(42).eq.0) goto 120 
                                                                        
       write(42,*) '                  ' 
        write(42,*) 'jmax_eci=',jmax_eci 
        do i = 1,neci 
          do j=1,jmax_eci(i) 
           write(42,904) 'i_eci_iglob1(',j,',',i,')=',i_eci_iglob1(j,i) 
           write(42,904) 'i_eci_iglob2(',j,',',i,')=',i_eci_iglob2(j,i) 
           do mm=1,nptseci 
              write(42,900) 'mm=',mm,eci_cs(mm,j,i),eci_eng(mm,j,i) 
           enddo 
         enddo 
        enddo 
                                                                        
  900 format(a3,i2,2x,1pE14.7,2x,1pE14.7) 
  901     format(a29,a17,i3,3x,i3,4x,1pE11.4) 
  904 format(a13,i3,a1,i1,a2,i3) 
  120 continue 
      return 
  915 format(12X,A3,16X,A3,3x,i6,i6,1X,E11.4) 
      END SUBROUTINE READ_ECIECR                                           

!***********************************************************************
      SUBROUTINE READ_LEVELS
!***********************************************************************
! This program reads cats levels files and                                
! determines the statistical weight of each levels.                       
!*  iglobindx: is the global index beginning with the lowest energy leve
!*             in the lowest ionization stage.                          
                                                                        
      use mac_commons
      use gms
      implicit real*8 (a-h,o-z) 
!     parameter(num=200,mum=99,lnum=99) 
      character*1 orbital_aux,blk,char
      character*45 config_aux 
      character(len=100),external::dir
      character(len=9),dimension(4) ::gfile           !hardwired
      dimension totalj(num)
      integer Translate 
      external Translate 
      do i=1,num 
         totalj(i)=0.0
         do j=1,num 
            gf(i,j)=0.0d0 
         enddo 
      enddo 
      do i=1,num
         do j=1,mum
            do k=1,3
               ioracle(i,j,k)=0
            enddo
         enddo
      enddo

      gfile(1)=element//'/gAg.01'
      gfile(2)=element//'/gAg.02'
      gfile(3)=element//'/gAg.03'
      gfile(4)=element//'/gAg.04'

      igrounds(1)=1 
      iglobindx=1 
      sum_ip = 0.0d0 
      i=0 

!**** Begin 1st do while                                                
   25 continue 
      i=i+1              !i - counts the ionization
       p_eng=0.0d0 
         if(i.eq.(nterms+1)) goto 40 
         cauxdr=gfile(i)
         cauxdr=dir(workdir,gfile(i))
         open(unit=21,file=cauxdr)

         cauxdr=fname_cat(i)
         cauxdr=dir(workdir,cauxdr)
         open(unit=20,file=cauxdr) 
         read(20,*) 
         read(20,*)
         sum_ip = sum_ip + rioneng(i-1) 

! --------------------------------------------------------------------------
!  This section reads all atomic structure cats data files
!  for all ionization stages.
! --------------------------------------------------------------------------
! ntermsindex(i) carries the # of (levels/terms/configs) for each ionization stage.
! ioracle -> contains the mapping between the different representations:
!            levels/terms/configs  ioracle(i,j,k) i=energies, j=ions, k=rep.
 
      do ii=1,ntermindex(i) !ii - counts the # (levels/terms/configs) per ion  
       if(itype(i) .eq.1) then ! levels
         read(20,911,END=30) index_aux,                                 &
     &                       iconfig_aux,                               &
     &                       mult_aux,                                  &
     &                       orbital_aux,                               &
     &                       totalj_aux,                                &
     &                       energy_aux,                                &
     &                       config_aux                                 
         icptype_aux = 1       
         ioracle(ii,i,1) = index_aux
         ioracle(ii,i,3) = iconfig_aux
       elseif(itype(i) .eq. 2) then ! terms                                                       
         read(20,912,END=30) index_aux,                                 &
     &                       iconfig_aux,                               &
     &                       mult_aux,                                  &
     &                       orbital_aux,                               &
     &                       energy_aux,                                &
     &                       config_aux
         totalj_aux  = 0
         icptype_aux  = 2 
         ioracle(ii,i,2) = index_aux
         ioracle(ii,i,3) = iconfig_aux
       elseif(itype(i) .eq. 3) then ! configs
         read(20,913,END=30) index_aux,                                 &
     &                       energy_aux,                                &
     &                       config_aux
         iconfig_aux = index_aux
         orbital_aux = ' '
         mult_aux    = 0
         totalj_aux  = 0
         icptype_aux  = 3 
         ioracle(ii,i,3) = index_aux
       else
         print*,'ERROR READ_LEVELS: itype not an acceptable value'
         stop 'ERROR READ_LEVELS: itype not an acceptable value'
       endif
911 FORMAT(I7,I7,9X,I4,A1,F9.1,3X,1pE14.7,5X,A45)
912 FORMAT(4X,I3,5X,I2,7X,I1,A1,3X,1pE14.7,5X,A45)
913 FORMAT(I5,1x,1pe14.7,2x,A45)
! --------------------------------------------------------------------------
! ntermsindex(i) carries the # of (levels/terms/configs) for each ionization stage. 

! This section expects that the energy of the cats output be in ascending order

           index(iglobindx)    =  index_aux 
           iconfig(iglobindx)  =  iconfig_aux 
           mult(iglobindx)     =  mult_aux 
           orbital(iglobindx)  =  orbital_aux 
           totalj(iglobindx)   =  totalj_aux
           energy(iglobindx)   =  energy_aux 
           config(iglobindx)   =  config_aux 
           icptype(iglobindx)   =  icptype_aux
          if(energy_aux.le.p_eng) write(istdout,*) fname_cat(i),'(',index_aux,')',&
     &       'is less than previous value'                              
           p_eng=energy_aux 
           glob_eng(iglobindx) = energy_aux + sum_ip 
           zeff(iglobindx)=i-1 
                                                                        
            blk=' ' 
            DO icnt=45,1,-1 
              IF(config_aux(icnt:icnt).ne.blk)THEN 
                 char = config_aux(icnt-3:icnt-3) 
                 CALL char_num(char,intr) 
                 Nvalshell(iglobindx) = intr 
                 GOTO 26 
              ENDIF 
            ENDDO 
   26      CONTINUE 
!-----------------------------------------------------                  
         if(itype(i).eq.1) then    ! levels
            ig(iglobindx) = NINT(2*totalj(iglobindx) + 1)  
         elseif(itype(i).eq.2)then ! terms 
            L = Translate(orbital(iglobindx))
            ig(iglobindx) = (2*L+1)*mult(iglobindx)
         elseif(itype(i).eq.3)then ! configs
            read(21,*) ig(iglobindx) 
         else
            print*,'ERROR: READ_LEVELS itype(',i,')=',itype(i),   &
     &      'not appropreate val'
            stop 'ERROR: READ_LEVELS itype' 
         endif

         iglobindx=iglobindx+1 

      enddo                                                                        
         
      close(20) 
      close(21) 
      igrounds(i+1)=iglobindx 
       write (istdout,*) 'igrounds(',i+1,')=',igrounds(i+1) 
      goto 25 
                                                                        
   40 continue 
!**** End of do while                                               

         modeltot=iglobindx - 1 

      if(modeltot+1.ge.num) write (istdout,*) 'warn: (modeltot+1) .ge. num',       &
     &   'problems may arise in killed_lev.f'                           
      write (istdout,*) 'modeltot=',modeltot 
! **********************************************************************
! READING gf's indexies                                                 
      j=1 ! ionization stage 1=I,2=II,etc 
   42 continue 
      cauxdr=fname_gf(j)
      cauxdr=dir(workdir,cauxdr)
      open(unit=21,file=cauxdr) 
      read(21,*) 
      read(21,*) 
   45 continue 
!        read(21,902,end=50) tempgf,indx1,indx2 
         read(21,902,end=50) tempgf,icat1,icat2 
         indx1 = iquery(icat1,j) 
         indx2 = iquery(icat2,j)
         if((indx1.gt.0).and.(indx2.gt.0))then
             indx1glob = indx1+(igrounds(j)-1) 
             indx2glob = indx2+(igrounds(j)-1) 
             gf(indx1glob,indx2glob)=tempgf 
         else 
         goto 45 
         endif 
      goto 45 
   50 continue 
      close(21) 
      j=j+1 
      if(j.eq.(ngf+1)) goto 55 
      goto 42 
   55 continue 
                                                                        
! ----------------------------------------------------------------------
! DIAGNOSTICS                                                           
      if (idiag(46).eq.0) goto 67 
      do ii=1,modeltot 
        do jj=1,modeltot 
         write(46,904) 'gf(',ii,',',jj,')=',gf(ii,jj) 
        enddo 
      enddo 
   67 continue 
                                                                        
! -------------------------------------------                           
! DIAGNOSTICS                                                           
      if (idiag(45).eq.0) goto 89 
       do k=1,modeltot 
         auxeng=energy(k)*8066.188 
          write(45,901)    k,                                           &
     &                     index(k),                                    &
     &                     ig(k),                                       &
     &                     mult(k),                                     &
     &                     orbital(k),                                  &
     &                     totalj(k),                                   &
     &                     energy(k),                                   &
     &                     auxeng,                                      &
     &                     glob_eng(k),                                 &
     &                     config(k),                                   &
     &                     icptype(k) 
   
                                    !energy=[eV]                    
       enddo 
      tmax=45
      write(45,*) 'ORACLE 3-D ARRAY MAX # LEV ENTRIES',TMAX
      do ii=1,tmax
         write(45,'(i4,4x,5(2x,i3))') ii,(ioracle(ii,jj,1),jj=1,5)
      enddo
      write(45,*)
      do ii=1,tmax
         write(45,'(i4,4x,5(2x,i3))')  ii,(ioracle(ii,jj,2),jj=1,5)
      enddo
      write(45,*)
      do ii=1,tmax
         write(45,'(i4,4x,5(2x,i3))')  ii,(ioracle(ii,jj,3),jj=1,5)
      enddo

   89  continue 
! **********************************************************************
  901 FORMAT(4X,I3,2X,I3,2x,I4,7X,I1,A1,f9.1,3X,1pE14.7,3X,1pE14.7,3X,       &
     &       1pE14.7,5X,A45,3x,I1)                                            
  902 FORMAT(25x,1pE11.4,72x,I7,I7) 
  903 FORMAT(25x,1pE11.4) 
  904 FORMAT(a3,I4,a1,I4,a2,1pE11.4) 
  905 FORMAT(1pE11.4) 
      return 
30    continue
      print*,'ERROR: less entires in ',cauxdr,' then expected'
      stop 'ERROR: less entires in cauxdr then expected'
      END SUBROUTINE READ_LEVELS                                          
!***********************************************************************
      integer function iquery(ival,ionix)
!***********************************************************************
! If association is not found for ival 
! then iquery=-1 
      use mac_commons
      implicit real*8 (a-h,o-z)
      iquery = -1 
      do ii=1,ntermindex(ionix)
       io=ioracle(ii,ionix,itype(ionix)) 
       if(io.eq.ival) then
         iquery=ii
         return 
       endif
      enddo  

      end function iquery
                                                                        
!***********************************************************************
      SUBROUTINE READ_PIRR(ierror_cs) 
! Called by mac_initial                                                 
!***********************************************************************
!  This subroutine reads the photoionization atomic data               *
!***********************************************************************
!*    last revision 7/8/97 1:20pm                                       
!*  j > i                                                               
!*  NOTE: rIntensity -> [number of photons/(cm^2*s)]                    
!*        Ephoton    -> [eV]                                            
!*        Te         -> [K]                                             
! NOTE: PIRR ASSUMES ECIECR c.s REPRESENTATION (lev,terms,configs)
      use mac_commons
      use gms
      implicit real*8 (a-h,o-z) 
!     parameter(num=200,mum=99,lnum=99,kmum=99) 
!     parameter(npts=21,numcs=5000,numstages=8) 
      integer ivec1(900),ivec2(900),iv1(9000),iv2(9000)
      dimension tmpcs(2),v_cs(nptseci),v_eng(nptseci)
      character(len=100), external :: dir
! -- Constants -------------------------------------------------------- 
      alf = 0.0d0 
      rkeVK=8.617385d-5 
      remass = 9.1093897d-31 

!*********** BEGIN 1ST DO WHILE LOOP (COUNTING OFF ION-STAGES)          
      ierror_cs=0
      i=0 
   25 continue 
                              ! i counts off the lower ionization stage.
      i=i+1 
      if(i.eq.(npi+1))goto 110 
      cauxdr=fname_pi(i)
      cauxdr=dir(workdir,cauxdr)
      open(unit=21,file=cauxdr) 
      read(21,*) 
      read(21,*) 
      read(21,*) 
                                                                        
!************ BEGIN 2ND DO WHILE LOOP (COUNTING OFF RATES)              

      j=0 
   30 continue 
      read(21,915,end=100) gridxx,engxx,icat1,icat2,deltaeng 
      i1aux=ichk(icat1,i,iecitype(i))
      i2aux=ichk(icat2,i+1,iecitype(i))            

      if((i1aux.gt.0).and.(i2aux.gt.0))then   ! Check indexies if the exist in ion. stage

! -- Mapping ionization c.s. indexies

        if((iecitype(i).eq.itype(i)).and.(iecitype(i).eq.itype(i+1)))then
           iv1(1)=i1aux
           iv2(1)=i2aux
           ktot=1
        elseif((iecitype(i).ne.itype(i)).and.(iecitype(i).eq.itype(i+1)))then
!          print*,'READ_PIRR: lower index mappings req'
           call queryall(icat1,i,iecitype(i),ivec1,k1)
           ktot=k1
           do ijj=1,k1
              iv1(ijj)=ivec1(ijj)
              iv2(ijj)=i2aux
           enddo
        elseif((iecitype(i).eq.itype(i)).and.(iecitype(i).ne.itype(i+1)))then
!          print*,'READ_PIRR: upper index mappings req not TESTED'
           call queryall(icat2,i+1,iecitype(i),ivec2,k2)
           ktot=k2
           do ijj=1,k2
              iv1(ijj)=i1aux
              iv2(ijj)=ivec2(ijj)
           enddo
        elseif((iecitype(i).ne.itype(i)).and.(iecitype(i).ne.itype(i+1)))then
!          print*,'READ_PIRR: 2 mappings req not TESTED'
           call queryall(icat1,i,iecitype(i),ivec1,k1)
           call queryall(icat2,i+1,iecitype(i),ivec2,k2)
           ktot=k1*k2
           ict=0
           do ij1=1,k1
           do ij2=1,k2
              ict=ict+1
              iv1(ict)=ivec1(ij1)
              iv2(ict)=ivec2(ij2)
           enddo
           enddo
        else
           stop 'ERROR: READ_ECI INDEX CONDITIONAL TRAP'
        endif

! -- Check c.s.

         iswitch=0 
         ier=0 
         do mm=1,nptsrr 
           read(21,*) eng,cs 
           if (cs.le.0.0d0) cs=dabs(cs) 
             if(mm.eq.1)then 
               tmpcs(2)=cs 
             else 
               tmpcs(1)=tmpcs(2) 
               tmpcs(2)=cs 
               if(iswitch.eq.0)then 
                 if(tmpcs(2).lt.tmpcs(1)) iswitch=1 
               endif 
                 if((iswitch.eq.1).and.(tmpcs(2).gt.tmpcs(1)).and.     &
     &              (ier.eq.0))then                                     
                   ier=1 
                   ierror_cs=1 
                   diflog=DABS(DLOG10(tmpcs(2))-DLOG10(tmpcs(1))) 
                   write(40,901) 'warn: irreg. cs detected in ',fname_pi(i), &
     &                 icat1,icat2,diflog                                   
                                           ! End Checking the form of CS
               endif 
           endif 
!          rr_cs(mm,j,i)  = cs/8.797356697d-17 
!          rr_eng(mm,j,i) = eng 
           v_cs(mm) = cs/8.797356697d-17
           v_eng(mm)= eng
         enddo  !mm loop
      
        do ijj=1,ktot
          i1=iv1(ijj)
          i2=iv2(ijj)
          j=j+1
          jmax_rr(i)=j
          iglob1 = (igrounds(i)-1)+i1
          iglob2 = (igrounds(i+1)-1)+i2
          i_rr_iglob1(j,i) = iglob1
          i_rr_iglob2(j,i) = iglob2
          do mm=1,nptsrr
             rr_cs(mm,j,i) = v_cs(mm)
             rr_eng(mm,j,i)= v_eng(mm)
          enddo
        enddo  !ijj loop
                                                                  
        goto 30 
      else 
         do mm=1,nptsrr
            read(21,*) eng,cs 
         enddo 
         goto 30 
      endif 
!*********************  END OF 2ND DO WHILE LOOP ********************** 
  100 continue 
      close(21) 
       write (istdout,*) '# CS-PI(',i,')=',jmax_rr(i) 
      goto 25 
  110 continue 
                                                                        
!*********************  END OF 1ST DO WHILE LOOP ********************** 
                                                                        
! --------------------------------------------------------------------- 
! DIAGNOSTICS                                                           
      if(idiag(44).eq.0) goto 120 
      write(44,*) '                  '                                                                  
      write(44,*) 'jmax_rr=',jmax_rr 
      do i = 1,npi 
        do j=1,jmax_rr(i) 
          write(44,904) 'i_rr_iglob1(',j,',',i,')=',i_rr_iglob1(j,i) 
          write(44,904) 'i_rr_iglob2(',j,',',i,')=',i_rr_iglob2(j,i) 
          do mm=1,nptsrr 
           write(44,900) 'mm=',mm,rr_eng(mm,j,i),rr_cs(mm,j,i)*8.797356697d-17
          enddo 
         enddo 
      enddo 
                                                                        
  900 format(a3,i2,2x,1pE14.7,2x,1pE14.7) 
  901     format(a29,a17,i3,3x,i3,4x,1pE11.4) 
  904 format(a12,i3,a1,i1,a2,i3) 
  120 continue 
      return 
  915 format(12X,A3,16X,A3,3x,i7,i7,1X,1pE11.4) 
      END SUBROUTINE READ_PIRR                                          

!********************************************************************** 
      SUBROUTINE SRD 
!********************************************************************** 
!NOTE: This routine must called over and over due to change in
!      frac(j) values during continum lowering

      use mac_commons

      implicit real*8 (a-h,o-z) 
!     parameter(num=200,mum=99,lnum=99) 
      rsrd = 0.0d0 
      do j=1,modeltot 
         do i=1,modeltot 
            rsrd(j,i)=4.34d7*(energy(j)-energy(i))**2*                     &
     &             (gf(i,j)/(1.0d0*ig(j)*frac(j)))                      
         enddo 
      enddo 
! -------------------------------------------------------------         
! DIAGNOSTICS:                                                          
                                                                        
      if(idiag(56).eq.0)goto 90 
      do l=1,modeltot 
         do k=1,modeltot 
              if(rsrd(l,k).ne.0.0d0)then 
               write(56,901) l,k,rsrd(l,k) 
              endif 
        enddo 
      enddo 
   90 continue 
  901 format (i7,i7,1pE14.7) 
      return 
      END SUBROUTINE SRD
                                                                        
!***********************************************************************
      SUBROUTINE STEWART_PYATT 
!***********************************************************************
      use mac_commons
      implicit real*8 (a-h,o-z) 
       external gmzbar 
!     parameter(num=200) 
      data itag/ 1 / 
      save itag 
      pi = 3.1415927 
      e  = 4.8032d-10                      ! ergs/K 
      rk = 1.3807d-16 
                                           ! for T_F model              
      z_charg    = 6 
                                           ! for T_F model              
      atomass    = 12.011 
                                           ! 1.1331d23atoms in 2.26g/cc 
      rho = rNa*(atomass/6.022d23) 
                                                                        
      if(itag .eq. 1)then 
!     write(istdout,*) 'WARN:  STEWART_PYATT.F SET TO CARBON !!!!!!!!!!!!!' 
      itag=0 
      endif 
       if(iPress_Ion.eq.1)then 
! =======================================================               
              call tfermi(z_charg,atomass,rho,Te_eV,TF_zbar) 
! =======================================================               
              T = Te_eV/8.617383d-5 
              R = (3.0d0/(4.0d0*pi*rNa))**(1.0/3.0) 
              gamma = (TF_zbar**2 * e**2)/(R*rk*T) 
              auxx = ((3.0d0*gamma)**1.5 + 1.0d0)**0.66667 
              aux_lower = ((rk*T)/(2.0d0*TF_zbar))*(auxx-1.0d0) 
              Eng_lower = aux_lower/1.602177d-12 
                                                                        
! =======================================================               
!          zbarlte = gmzbar(Te_eV,rNa,istdout)                                  
! =======================================================               
!                                                                       
!          if(zbarlte .lt. 1.0d0)then                                   
!              zbarlte = 1.0d0                                          
!            endif                                                      
!           TF_zbar = zbarlte      ! remnant from T-F for printing      
!              rNe_lte = zbarlte*rNa                                    
!              T = Te_eV/8.617383d-5                                    
!              R = (3.0d0/(4.0d0*pi*rNa))**(1.0/3.0)                    
!           Debye = DSQRT(rk*T/(4.0*pi*(e**2)*(zbarlte+1)*rNe_lte))     
!         Eng_lower = Te_eV*(((R/Debye)**3.0 + 1.0d0)**0.6666-1.0d0)/   
!     &                (2.0d0*(zbarlte+1.0d0))                          
                                                                        
       else 
                              ! No contiume lowering                    
              Eng_lower = 0 
       endif 
!===============================================================        
      call killed_lev 
!===============================================================        
!     write(40,900)rNa,Te_eV,gamma,Eng_lower                            
  900 format(1pE12.5,4x,1pE10.3,4x,1pE10.3,4x,1pE10.3) 
      END SUBROUTINE STEWART_PYATT

!***********************************************************************
!  This subroutine produces the filenames of the output files          *
      SUBROUTINE TAILNAMES(tnames) 
!***********************************************************************
!*  NOTE: THIS SUBROUTINE WILL ONLY PRODUCE 999 FILES!                  
       implicit real*8 (a-h,o-z) 
                                                                        
       character*1 ones(10) 
       character*1 tens(10) 
!      character*1 hundereds(10)                                        
!      character*8 prefix                                               
!      character*3 sufix                                                
       character*2 tnames(0:99) 
                                                                        
       data ones     /"0","1","2","3","4","5","6","7","8","9"/ 
       data tens     /"0","1","2","3","4","5","6","7","8","9"/ 
!      data hundereds/"0","1","2","3","4","5","6","7","8","9"/          
                                                                        
       jj=-1 
                                                                        
!      do 20 i=1,10                                                     
       do 30 j=1,10 
       do 40 k=1,10 
       jj=jj+1 
                                                                        
       tnames(jj)=tens(j)//ones(k) 
                                                                        
   40  continue 
   30  continue 
!0     continue                                                         
      END SUBROUTINE TAILNAMES                                           
                                                                        
! **********************************************************************
      SUBROUTINE TEST_KERNAL(ntotal) 
! **********************************************************************
! This program is designed to sum the col. of the matrix A by summing th
! contributions of the individual atomic processes first and then perfor
! a final sum of thoses components.  This test was developed after non-z
! elements were detected after a conventional summing of the components 
! the matrix A.                                                         
! NOTE: At this time the different atomic processes must be hardwired.  
!       This routine is not automated.                                  
! **********************************************************************
      use mac_commons                                                                        
      implicit real*8 (a-h,o-z) 
!     parameter(num=200,mum=99,nnum=15,lnum=99,kmum=99) 
      dimension a_eced(num,num),                                        &
     &          a_ecir(num,num),                                        &
     &          a_rate(num,num),                                        &
     &          auxed(num),                                             &
     &          auxir(num),                                             &
     &          auxrate(num)                                            
! --------------------------------------------------------------------- 
        do n = 1,ntotal 
! --------------------------------------------------------------------- 
! Solving for a(n,m) for m < n                                          
                                                                        
          if(n.eq.1) goto 120 
                                                                        
          do m = 1,n-1 
              a_eced(n,m) = receecd(m,n) 
              a_ecir(n,m) = reciecr(m,n) 
              a_rate(n,m) = rpi(m,n) 
!             a_rate(n,m) = rmpi(m,n)                                   
          enddo 
! --------------------------------------------------------------------- 
! Solving for a(n,n)                                                    
                                                                        
  120     aux1 = 0.0d0 
          aux2 = 0.0d0 
          aux3 = 0.0d0 
                                                                        
          if (n.eq.1) goto 140 
                                                                        
          do m = 1,n-1 
              aux1 =  receecd(n,m) + aux1 
              aux2 =  reciecr(n,m) + aux2 
!             aux3 =  0.0d0 + aux3           ! pi,rmpi                  
!             aux3 =  rsrd(n,m)                                         
!             aux3 =  rrr(n,m)                                          
!             aux3 =  rdr(n,m)                                          
          enddo 
                                                                        
  140     if (n.eq.ntotal) goto 165 
          aux11 = 0.0d0 
          aux12 = 0.0d0 
          aux13 = 0.0d0 
                                                                        
          do m = n+1,ntotal 
              aux11 = receecd(n,m) + aux11 
              aux12 = reciecr(n,m) + aux12 
              aux13 = rpi(n,m)     + aux13 
!             aux13 = rmpi(n,m)    + aux13                              
!             aux13 = 0.0d0 + aux13           ! rsrd,rrr,rdr            
          enddo 
                                                                        
          auxed(n)   = aux1 + aux11 
          auxir(n)   = aux2 + aux12 
          auxrate(n) = aux13 
                                                                        
  165     auxir(ntotal)=aux2 
! --------------------------------------------------------------------- 
! Solving for a(n,m) m > n                                              
                                                                        
          if (n.eq.ntotal) goto 180 
                                                                        
          do m = n+1, ntotal 
              a_eced(n,m) = receecd(m,n) 
              a_ecir(n,m) = reciecr(m,n) 
                                               !pi,mpi                  
              a_rate(n,m) = 0.0d0 
!             a_rate(n,m) = rsrd(m,n)                                   
!             a_rate(n,m) = rrr(m,n)                                    
!             a_rate(n,m) = rdr(m,n)                                    
          enddo 
  180     continue 
! --------------------------------------------------------------------- 
          enddo 
                                                                        
          write(37,*)'                      ' 
          write(37,*)'col# ! sum ! s_ed ! s_ir ! s_rate !' 
                                                   ! col.               
          do m = 1,ntotal 
             sum_ed = 0.0d0 
             sum_ir = 0.0d0 
             sum_rate = 0.0d0 
                                                   ! row                
             do n = 1,ntotal 
                sum_ed = a_eced(n,m) + sum_ed 
                sum_ir = a_ecir(n,m) + sum_ir 
                sum_rate = a_rate(n,m) + sum_rate 
             enddo 
                                                                        
      write(38,*)'            ' 
      write(38,*)'m=',m 
      write(38,*)'sum_ed=',sum_ed,'  auxed=',auxed(m) 
      write(38,*)'sum_ir=',sum_ir,'  auxir=',auxir(m) 
      write(38,*)'sum_rate=',sum_rate,'   auxrate=',auxrate(m) 
                                                                        
             s_ed = sum_ed - auxed(m) 
             s_ir = sum_ir - auxir(m) 
             s_rate = sum_rate - auxrate(m) 
             sum = s_ed + s_ir + s_rate 
                                                                        
          if(idiag(37).eq.0)goto 206 
             write(37,*) 'col=',m,'  ',sum,'  ',s_ed,                   &
     &          '   ', s_ir,'   ',s_rate                                
  206     continue 
          enddo 
      END SUBROUTINE TEST_KERNAL                          

!***********************************************************************
         SUBROUTINE TFERMI(zn,am,rho,t,zbar) 
!***********************************************************************
! CALCULATES AVERAGE IONIZATION ZBAR USING THE THOMAS-FERMI ATOMIC MODEL
! SCALING AND FITTING PARAMETERS FROM R.M. MORE, LLNL UCRL-84991 (1981) 
!                                                                       
         IMPLICIT REAL*8(A-H,O-Z) 
!        DIMENSION TEMP(10),DEN(10),BIGZ(10,10)                         
!        ZN=92.0D0                                                      
!        AM=238.0D0                                                     
!        WRITE(6,*) 'ENTER RHO (G/CC) AND T (EV):'                      
!        READ(5,*) RHO,T                                                
!CCC         DO 10 J=1,10                                               
!CCC         T=100.0D0+(J-1)*100.0D0                                    
!CCC         TEMP(J)=T                                                  
!CCC         DO 10 I=1,5                                                
!CCC         RHO=0.000010D0*10**I                                       
!CCC         DEN(I)=RHO                                                 
         ZBAR=TFZBAR(RHO,T,ZN,AM) 
!        WRITE(26,*) t,ZBAR                                             
!CCC         BIGZ(I,J)=ZBAR                                             
!CCC 10      CONTINUE                                                   
!CCC         WRITE(8,800) (TEMP(J),J=1,10)                              
!CCC         DO 20 I=1,5                                                
!CCC         WRITE(8,810) (BIGZ(I,J),J=1,10)                            
!CCC 20      CONTINUE                                                   
!CCC 800     FORMAT(8X,10E8.1,/)                                        
!CCC 810     FORMAT(E8.1,10F8.2,/)                                      
         return 
      END SUBROUTINE TFERMI                                          

!*****************************************************                  
      FUNCTION TFZBAR(RHO,T,ZN,AM) 
!*****************************************************                  
         IMPLICIT REAL*8(A-H,O-Z) 
!                                                                       
         A1=0.0033230D0 
         A2=0.9718320D0 
         A3=9.261480D-5 
         A4=3.101650D0 
         B0=-1.7630D0 
         B1=1.431750D0 
         B2=0.3154630D0 
         C1=-0.3666670D0 
         C2=0.9833330D0 
         ALPHA=14.31390D0 
         BETA=0.66240D0 
!                                                                       
         R=RHO/(ZN*AM) 
         T0=T/ZN**(4.0D0/3.0D0) 
         TF=T0/(1.0D0+T0) 
         A=A1*T0**A2+A3*T0**A4 
         B=-DEXP(B0+B1*TF+B2*TF**7) 
         C=C1*TF+C2 
         Q1=A*R**B 
         Q=(R**C+Q1**C)**(1.0D0/C) 
         X=ALPHA*Q**BETA 
         FX=X/(1.0D0+X+DSQRT(1.0D0+2.0D0*X)) 
         TFZBAR=ZN*FX 
         RETURN 
      END FUNCTION TFZBAR                                          
! **********************************************************************
      SUBROUTINE rad_transdat
! **********************************************************************
! Purpose: Reread gfXX.?? files to obtain radiation transport data
! 
      use mac_commons
      use gms

      implicit real*8 (a-h,o-z)
      integer(4),dimension(10000) ::istatwtup,istatwtlow
      character(len=100),external::dir
      icnt=0
      j=1 ! ionization stage 1=I,2=II,etc
   42 continue
      cauxdr=fname_gf(j)
      cauxdr=dir(workdir,cauxdr)
      open(unit=21,file=cauxdr)
      read(21,*)
      read(21,*)
   45 continue
!  icat2 > icat1
         read(21,902,end=50) rl,e,ggf,rloggf,fa,fe,aa,ae,rlifet,icat1,icat2
         indx1 = iquery(icat1,j)
         indx2 = iquery(icat2,j)
         if((indx1.gt.0).and.(indx2.gt.0))then
             if(me.eq.0)then
               e=(energy(indx2+(igrounds(j)-1))) - (energy(indx1+(igrounds(j)-1)))
             endif 
             if((e.lt.hnu1).or.(e.gt.hnu2)) goto 45
             indx1glob = indx1+(igrounds(j)-1)
             indx2glob = indx2+(igrounds(j)-1)
             icnt=icnt+1
             call einsteinb(aa,e,indx1glob,indx2glob,b_stmemm,b_phtabs)
             radat(icnt,1)=rl
             radat(icnt,2)=e
             radat(icnt,3)=aa
             radat(icnt,4)=b_stmemm
             radat(icnt,5)=b_phtabs
             radat(icnt,6)=indx1glob
             radat(icnt,7)=indx2glob           
             radat(icnt,8)=ig(indx1glob) !lower
             radat(icnt,9)=ig(indx2glob)  !upper
         else
         goto 45
         endif
      goto 45
   50 continue
      close(21)
      j=j+1
      if(j.eq.(ngf+1)) goto 55
      goto 42
   55 continue
      iradat=icnt
! -------------------------------------------------------------
! DIAGNOSTICS:

      if(idiag(47).eq.0)goto 90
      do ii=1,icnt
       write(47,901)(radat(ii,j), j=1,9)
      enddo  
      close(47)      
90    continue 
902   FORMAT( 9(1x,E11.5), 2(i7))
901   FORMAT( 2x,1e11.5,2x,F11.5,2x,3(2x,e11.5),2(2x,f7.1),2x,3(1e11.5))
      END SUBROUTINE RAD_TRANSDAT

! **********************************************************************
      SUBROUTINE EINSTEINB(aa,e,ilow,iup,b_stmemm,b_phtabs)
! **********************************************************************
! A(i,j)=rsrd(i,j)
! A(initial,final) = A(upper,lower)
! B(upper,lower)
! pi = 2.0*ASIN(1.0)
! B(j,i)=A(j,i)*c^2/(2h*v^3) v=[sec^-1]   NOTE: DIFFERENT FROM COWAN BY 4PI
      use mac_commons
      use gms

      implicit real*8 (a-h,o-z)

      h=6.6260755e-27 ! ergs*s
      c=2.99792458d10  ! cm/s

      const = c**2/(2*h) != 6.7818560099e+46 cm^2/(ergs*s^3) = 6.7819569724d46

      if(srcnum.ne.num) then
        write(istdout,*) 'ERROR EINSTEINB : SRCNUM.NE.NUM'
        STOP
      endif

      ithin=0  ! ithin = 1 -> opt. thin; ithin = 0 -> opt. thick
!     write(istdout,*) 'WARN: ithin=',ithin,' 1 -> opt. thin;',         &
!     &       ' ithin = 0 -> opt. thick'
      if(ithin.eq.1) goto 80
          v=e*2.41799e14 ! eV -> Hz
          b_stmemm = aa*const/(v**3)
          b_phtabs = (dble(ig(iup))/dble(ig(ilow)))*b_stmemm
80    continue

      END SUBROUTINE EINSTEINB 

! **********************************************************************
      SUBROUTINE SRD_SUM
! **********************************************************************
! NOTE: j>i

      use mac_commons
      use gms
      implicit real*8 (a-h,o-z)

      do j=2,modeltot
        sum=0
        do i=1,(j-1)
          sum=sum+rsrd(j,i)
        enddo
        srdsum(j)=sum
        write(istdout,*) 'srdsum(',j,')=',srdsum(j)
      enddo

      END SUBROUTINE SRD_SUM
      END MODULE MAC_CODE

!***********************************************************************
!  This subroutine translates s->0,p->1,d->2,...etc                     
      INTEGER FUNCTION TRANSLATE(L) 
!***********************************************************************
      use mac_commons
      implicit real*8(a-h,o-z) 
      character*1 L 
      if (L.eq.'s') then 
         Translate=0 
      elseif (L.eq.'p') then 
         Translate=1 
      elseif (L.eq.'d') then 
         Translate=2 
      elseif (L.eq.'f') then 
         Translate=3 
      elseif (L.eq.'g') then 
         Translate=4 
       elseif (L.eq.'h') then 
          Translate=5 
      else 
         write (istdout,*) 'Problem in Function Translate' 
      endif 
      return 
      END FUNCTION TRANSLATE    
!*********************************************************************
      subroutine copyary
!*********************************************************************
      use mac_commons
      do i=1,num
         fp(i) = f(i)
      enddo
      end
!***********************************************************************
      subroutine converge(iconvrg,iconsum)
!***********************************************************************
      use mac_commons
      double precision reldiff(num)
      integer iconvrg,iconsum
      iconsum = modeltot 
      iconvrg = 1  ! 0-> not converged; 1-> has converged
      do i=1,modeltot
!        reldiff(i)= ( abs( f(i) - fp(i) ) )/((f(i)+fp(i))/2.0)
         reldiff(i)  = dabs(f(i)-fp(i))/(fp(i)+delta)
         if( reldiff(i) .gt. eps ) then
            iconsum = iconsum - 1  !# lev not converged
            iconvrg = 0
         endif
      enddo
      end
!***********************************************************************
  character(len=*) function dir(string1,string2)
!***********************************************************************
  implicit none
  character(len=*),intent(in)  :: string1,string2
  dir=trim(adjustl(string1))//adjustl(string2)
  end function dir

