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 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. ! mm: 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 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) 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