c ***********************786*********************************** c
c
c     program lifetmss.f
c
c * A revised version. - Sultana N. Nahar, OSU, June 2020
c
c ************************************************************* c
c
c *  Program to read file "ion.f.bp.txt" for energy and oscillator 
c strengths/A-values etc. It then ccalculates lifetimes from A-values
c
c  It reads the A-values from file "efbpss.ion" computed by SUPERSTRUCTURE
c  It reads two input files:
c * 1. lftmin = input parameter file (unit 10). It contains:
c nlvi=number of levels for which lifetimes are to be calculated. But
c current option is free - any negative value for all levels, any
c positive value for specific levels
c isi,ili,ipi,aj,ici=spin multiplicity, angular momentum, parity, j-value
c  and configuration number of the interested level
c * 2. "efbpss.ion"/"ion.f.bp.txt" file (unit 20) for energy and 
c transition parameters
c
      parameter(nl=50,ns=1000)
      character*1 prm*3,l0(ns),p0(ns),sr0(ns),ili(ns),pri(ns),sri(ns),
     1trsn(6)*3,trs*3,p1*1,p2*1,l1*1,l2*1,pr(2)*1,cfg(ns)*14,
     1etyp*3,sr(nl)*1,cnfg*70,lftmin*7,efion*20
      dimension is0(ns),j20(ns),ic0(ns),en0(ns),lv0(ns),isi(ns),ici(ns),
     1j2i(ns),engi(ns),lvi(ns),aa(2),ss(2),ne(nl),eni(nl)
      data pr/'e','o'/,trsn/'E1d','E1i','E3 ','E2 ','M2 ','M1 '/
     1cm/109737.31/
  101 format(i5,2x,i1,2a1,1x,i2,1x,a1,1x,i2,24x,a14,7x,1pe11.4)
  103 format(/,5x,'LS cf_i',3x,'gi  lvi',2x,' <-',3x,'LS cf_j',2x,
     1'gj  lvj',1x,'f(E1)/S(E2',3x,'Aji',6x,'Eij',/,47x,'E3,M1,M2)',
     13x,'(sec-1)',3x,'(A)',/)
  105 format(i5,2x,i1,2a1,1x,i2,1x,a1,1x,i2,13x,1pe12.5)
  111 format(' Radiative decay rates of level j to various levels',
     1'(j -> i):')
c
c * print names of input/output files
c
      print 100
  100 format(/,' print names of input files: lftmin and',
     1'ion.f in a7,a20')
      read(5,'(a7,a20)')lftmin,efion
c
      open(unit=10,file=lftmin,status='old')
      open(unit=20,file=efion,status='old')
      iu=20
c
      print 116
  116 format(/,' output files: fort.7, fort.9',/)
c
c * read and write ion and configurations
c
      call efbpsscf(iu,etyp)
c
      write(7,'(a1)')
      write(9,'(a1)')
c
c * get fine structure energy levels of the ion
c
      call efbpsslv(iu,etyp,lv0,is0,l0,p0,ic0,sr0,j20,en0,ien,cfg)
c
c * read the levels for which lifetimes are to be calculated: nlvi is
c negative for all levels, positive for specific levels 
c
      read(10,*)nlvi
c
      call selectlv(isi,ili,pri,j2i,sri,ici,lvi,engi,nlvi,is0,l0,
     1p0,j20,sr0,ic0,lv0,en0,ien)
c
      print 150,nlvi
  150 format(' total levels=',i5)
c
c * start on calculations of lifetimes
c
      write(7,111)
      write(7,103)
c
c * find the level for which lifetime is to be calculated
c
      do 15 ie=1,nlvi 
c
c * print level on screen and in files fort.7 and fort.9
c
      print 151,lvi(ie)
  151 format(' ss-level=',i5)
c
      ji=j2i(ie)-1
      ecm=engi(ie)*cm
      write(7,127)lvi(ie),ici(ie),sri(ie),isi(ie),ili(ie),pri(ie),
     1ji,engi(ie),ecm
      write(9,127)lvi(ie),ici(ie),sri(ie),isi(ie),ili(ie),pri(ie),
     1ji,engi(ie),ecm
c
  127 format(' lifetime: sslevel j=',i4,', cf=',i2,',',1x,a1,i1,2a1,
     1i2,'/2[E=',1pe10.3,' Ry=',1pe14.7,' /cm]')
c
c * initialize to start lifetime computation
c
c  11 igi=j2i(ie)+1
      igi=j2i(ie)
c
c * read transition tables of type E1d, E1i, and get sum of a (asm)
c
      asm=0.
      ism=0
c
      do 1 ii=1,2
      call trnse1(iu,ii,igi,pri(ie),lvi(ie),asm,ism)
    1 continue
c
c * read forbidden transitions and get sum of A-values
c
      smfa=0.
c
      do 2 ii=3,4
      call trnsfrbd(iu,ii,igi,pri(ie),lvi(ie),smfa,ism) 
    2 continue
c
c * get the total sum of A-values
c
      sma=asm+smfa
c
c * write the sum of A-values and lifetimes
c
      write(7,115)smfa,asm
      write(9,115)smfa,asm
      write(7,104)ism,sma
      write(9,104)ism,sma
c
      if (sma.gt.0.) then
      alftm=1./sma
      write(7,119)alftm
c     write(9,107)sri,isi,ili,ipi,ji2,engi,lvi,ici,alftm
      write(9,119)alftm
      else
      write(7,102)sri(ie),isi(ie),ili(ie),pri(ie),ji,engi(ie),
     1ici(ie),lvi(ie)
      write(9,102)sri(ie),isi(ie),ili(ie),pri(ie),ji,engi(ie),
     1ici(ie),lvi(ie)
      endif
c
  104 format(1x,'TotalSum(Af+Aa)= Aji(',i4,' transitions) to the ',
     1'level=',5x,1pe10.3,' s-1')
c 107 format(/,' lifetime:level ',a1,i1,2a1,i2,'/2 [E=',1pe10.3,
cc   1' Ry,cf=',i2,',sslvl=',i4,']=',1pe10.3,' s',/)
c    1' Ry,sslvl=',i4,',cf=',i2,'] =',1pe10.3,' s',/)
  102 format(' lifetime:level ',a1,i1,2a1,i2,'/2 [E=',1pe10.3,
     1' Ry,cf=',i2,',sslvl=',i4,'] =','Infy:No Transition',/)
  115 format(' Summed A-values: Af (forbidden)=',1pe10.3,', Aa ',
     1'(allowed)=',1pe10.3,' s-1')
c 117 format(/,' lifetime:level ',a1,i1,2a1,i2,'/2 [E=',1pe10.3,
c    1' Ry=',1pe14.7,' cm-1, sslvl=',i4,',cf=',i2,']',/)
  119 format(' Lifetime (=1/Aji)=',39x,1pe10.3,' s',/)
c
c * rewind file to read and compute lifetime of the next level
c
      rewind iu
c
   15 continue
c
      stop
      end
c
c ************************************************************* c
c
      subroutine efbpsscf(iu,etyp)
c
c ************************************************************* c
c
c  * routine to read ion information and configurations from file 
c "efbpss.ion". 
c
c * etyp="cal" -> calculated energies, ="exp" -> observed energies
c
      parameter(nl=50,ns=1000)
      character prm*3,cnfg*70,etyp*3
c
  100 format(/,' nz =',i4,2x,', No. of core electrons=',i3)
c
c * from "efbpss.ion" read and write ion information
c
    9 read(iu,'(5x,a3)')prm
      if (prm.ne.' nz') goto 9
c
      backspace iu
      read(iu,'(10x,i3,9x,i2)')nz,nelc
c
      write(7,100)nz,nelc
      write(9,100)nz,nelc
c
c * read and write number of configurations
c
   15 read(iu,'(5x,a3)')prm
      if (prm.ne.'con') goto 15
c
      backspace iu
      read(iu,'(i5)')ncf
      read(iu,'(a1)')
c
      write(7,120)ncf
      write(9,120)ncf
  120 format(/,' Number of configurations:',i5,/)
c
c * read and write out the configurations
c
      do 7 i=1,ncf
      read(iu,'(i3,1x,a70)')ic,cnfg
      write(7,'(i3,1x,a70)')ic,cnfg
      write(9,'(i3,1x,a70)')ic,cnfg
    7 continue
c
      read(iu,'(a1)')
c
c * find energies are from observed/calculated set from NOTE statement
c
      read(iu,'(1x,a3)')prm
c
      if (prm.eq.'NOT') etyp='exp'
      if (prm.ne.'NOT') etyp='cal'
c
      return
      end
c
c ************************************************************** c
c                786
      subroutine efbpsslv(iu,etyp,lv0,is0,l0,p0,ic0,sr0,j20,en0,ien,
     1cfg)
c
c ************************************************************** c
c    
c * routine to read the energy levels of the ion from efbpss.ion file 
c * isi,ili,ipi,aj,ici=spin multiplicity, angular momentum, parity, j-value
c  and configuration number of the interested level
c
      parameter(ns=1000)
      character*1 prm*3,etyp*3,l0(ns),sr0(ns),p0(ns),cfg(ns)*14
      dimension lv0(ns),is0(ns),ic0(ns),j20(ns),en0(ns)
c
c * pass lines to the fine structure energy table
c
    1 read(iu,'(11x,a3)')prm
      if (prm.ne.'fin') goto 1
c
c * ien=total number of calculated energy levels
c
      backspace iu
      read(iu,'(34x,i5)')ien
c
      write(7,106)ien
      write(9,106)ien
  106 format(/,' Total number of calculated levels=',i5,/)
c
    2 read(iu,'(2x,a3)')prm
      if (prm.ne.' ie') goto 2
      read(iu,'(a1)')
c
c * start reading levels and energies from efbpss.ion, and save them 
c
      do 3 i1=1,ien
c
      if (etyp.eq.'cal')
     1read(iu,105)lv0(i1),is0(i1),l0(i1),p0(i1),ic0(i1),sr0(i1),
     1j20(i1),en0(i1)
c
      if (etyp.eq.'exp')
     1read(iu,101)lv0(i1),is0(i1),l0(i1),p0(i1),ic0(i1),sr0(i1),
     1j20(i1),cfg(i1),en0(i1)
c
c     write(7,105)lv0(i1),is0(i1),l0(i1),p0(i1),ic0(i1),sr0(i1),
c    1j20(i1),en0(i1)
c
    3 continue
c
  101 format(i5,2x,i1,2a1,1x,i2,1x,a1,1x,i2,24x,a14,7x,1pe11.4)
  105 format(i5,2x,i1,2a1,1x,i2,1x,a1,1x,i2,14x,1pe12.5)
c 105 format(i5,2x,i1,2a1,1x,i3,1x,a1,1x,i2,13x,1x,1pe12.5)
c
      return
      end
c
c ************************************************************** c
c                786
      subroutine selectlv(isi,ili,pri,j2i,sri,ici,lvi,engi,nlvi,is0,l0,
     1p0,j20,sr0,ic0,lv0,en0,ien)
c 
c ************************************************************** c
c
c * routine to define the levels, all or selected, for which lifetimes 
c are to be calculated, It also identifies the location of the selected
c levels in the complete set of levels of the ion 
c
      parameter(ns=1000)
c     real*8 a-h,j2i,j20,o-z
      character*1 ili(ns),pri(ns),l0(ns),p0(ns),sri(ns),sr0(ns)
      dimension isi(ns),j2i(ns),ici(ns),lvi(ns),is0(ns),ic0(ns),
     1j20(ns),lv0(ns),en0(ns),engi(ns)
c
c * use initial sign of nlvi for selected/all levels for lifetime calculations
c
      if (nlvi.lt.0) then
c
      do 1 i=1,ien
      isi(i)=is0(i)
      ili(i)=l0(i)
      pri(i)= p0(i)
      j2i(i)=j20(i)
      sri(i)=sr0(i)
      ici(i)=ic0(i)
      lvi(i)=lv0(i)
      engi(i)=en0(i)
    1 continue
c
      nlvi=ien
      return
      endif
c
c * for calculation of lifetimes for selected levels
c
      if (nlvi.ge.0) then
c
c * ie=counter for the selected levels
c
      i=0
    3 i=i+1
c
c * read the symmetry of the interested level for lifetime
c
      read(10,100)isi(i),ili(i),pri(i),aj,ici(i)
c
      if (isi(i).le.0) goto 2
c
      j2i(i)=2*aj+1
c
c     print 100,isi(i),ili(i),pri(i),aj,j2i(i),ici(i)
c 100 format(i2,2a1,f4.1,i3,i5)
      goto 3
c
    2 nlvi=i-1
c
      print 150,nlvi
  150 format(' total levels=',i5)
c
c * identify the selected lelves in the set of SS energy levles
c
      do 5 ie=1,nlvi
      do 10 i1=1,ien
c
      if (is0(i1).eq.isi(ie).and.l0(i1).eq.ili(ie).and.p0(i1).eq.
     1pri(ie).and.ic0(i1).eq.ici(ie).and.j20(i1).eq.j2i(ie)) then
c
c * skip for duplicate level with same symmetry and configuration
c
      if (nlvi.gt.1) then
      do 4 i=1,ie
      if (engi(i).eq.en0(i1)) go to 10
    4 continue
      endif
c
      sri(ie)=sr0(i1)
      engi(ie)=en0(i1)
      lvi(ie)=lv0(i1)
      goto 5
      endif
c
   10 continue
    5 continue
c
      endif
c
      return
      end
c
c ************************************************************** c
c
      subroutine trnse1(iu,ii,igi,pri,lvi,asm,ism) 
c
c ************************************************************** c
c
c * routine to read the E1d and E1i transition probabilities, add the
c A-values decaying from level igi,pri,lvi to "asm" to calculate the
c lifetimes.
c
      parameter(ns=1000)
      character*1 prm*3,trsn(2)*3,trs*3,pri,l1,l2,p1,p2 
      data trsn/'E1d','E1i'/
c
c 112 format(1x,a3,1x,i1,2a1,1x,i2,4x,i2,1x,a1,1x,i4,7x,i1,2a1,1x,i2,3x,
c    1i2,1x,a1,1x,i4,1x,1pe9.2,1pe10.3,1pe11.4)
  112 format(1x,a3,1x,i1,2a1,1x,i2,4x,i2,1x,i4,7x,i1,2a1,1x,i2,3x,
     1i2,1x,i4,1x,1pe9.2,1pe11.3,1pe12.4)
c
c * check the transition tables of type E1d, E1i
c
      trs=trsn(ii)
    8 read(iu,'(16x,a3)')prm
      if (prm.ne.trs) goto 8
c
      do 12 i=1,3
   12 read(iu,'(a1)')
c
c * read transitions of E1d or E1i types
c
   14 read(iu,'(1x,a3)')prm
c
c * if not end of set, read transitions
c
      if (prm.ne.'   ') then
      backspace iu
c
c * E1d or E1i transitions
c
      read(iu,118)
     1i1,i2,is1,l1,p1,ic1,is2,l2,p2,ic2,ig1,ig2,wl,e1,e2,fl,s,a
c
  118  format(i3,1x,i3,2(1x,i1,2a1,i2),2i3,f10.2,2f8.2,1pe9.2,
     11pe11.3,1pe9.2)
c 118  format(2i4,2(1x,i1,2a1,i2),2i3,f10.2,2f8.2,1pe9.2,1pe10.3,
c    11pe9.2)
c
c * add A-values if the upper level is matched to the interested one
c
      if (igi.eq.ig2.and.pri.eq.p2.and.lvi.eq.i2) then
c
      asm=asm+a
      ism=ism+1
c
      write(7,112)trs,is1,l1,p1,ic1,ig1,i1,is2,l2,p2,ic2,
     1ig2,i2,fl,a,wl
c
      endif
c
      goto 14
      endif
c
      return
      end
c
c ************************************************************** c
c
      subroutine trnsfrbd(iu,ii,igi,pri,lvi,asm,ism) 
c
c ************************************************************** c
c
c * routine to read table of E3,M2 or E1,M1 transitions and sum the
c A-values for contributions to a-sum for forbidden transitions
c
      parameter(ns=1000)
      character*1 prm*3,trsn(6)*3,trs*3,pri,l1,l2,p1,p2 
      dimension ss(2),aa(2)
      data trsn/'E1d','E1i','E3 ','E2 ','M2 ','M1 '/
c
c * check the transition tables for type E2, E3, M1, M2
c
      trs=trsn(ii)
    8 read(iu,'(16x,a3)')prm
      if (prm.ne.trs) goto 8
c
      do 12 i=1,3
   12 read(iu,'(a1)')
c
c * read transitions of various types
c
   14 read(iu,'(1x,a3)')prm
c
c * if not end of the set of the set, read transitions
c
      if (prm.ne.'   ') then
      backspace iu
c
c * for forbidden transitions
c
      read(iu,114)
     1i1,i2,is1,l1,p1,ic1,is2,l2,p2,ic2,ig1,ig2,wl,ei,ef,ss(1),aa(1),
     1ss(2),aa(2)
c
  114 format(i3,i4,2(1x,i1,2a1,i2),2i3,f10.2,1p2e10.2,1p4e9.2)
c 114 format(2i4,2(1x,i1,2a1,i2),2i3,f10.2,1p2e10.3,1p4e9.2)
c
c * process to add A-values if the level is matched to the interested one
c
      if (igi.eq.ig2.and.pri.eq.p2.and.lvi.eq.i2) then
c
c * add contributions of two types of transitions, e.g. E2,M1 and E3,M2
c
      do 18 k=1,2
c
c * add contributions of two types of transitions, e.g. E2,M1 and E3,M2
c
      if (aa(k).gt.0.) then
      if (k.eq.1) trs=trsn(ii)
      if (k.eq.2) trs=trsn(ii+2)
      a=aa(k)
      fl=ss(k)
      endif
c
      asm=asm+a
      ism=ism+1
c
      write(7,112)trs,is1,l1,p1,ic1,ig1,i1,is2,l2,p2,ic2,
     1ig2,i2,fl,a,wl
c
   18 continue 
c
      endif
c
c 112 format(1x,a3,1x,i1,2a1,1x,i2,4x,i2,1x,a1,1x,i4,7x,i1,2a1,1x,i2,3x,
c    1i2,1x,a1,1x,i4,1x,1pe9.2,1pe10.3,1pe11.4)
  112 format(1x,a3,1x,i1,2a1,1x,i2,4x,i2,1x,i4,7x,i1,2a1,1x,i2,3x,
     1i2,1x,i4,1x,1pe9.2,1pe11.3,1pe12.4)
c
      goto 14
      endif
c
    3 return
      end