      PROGRAM ms 
      USE gms
      use communicate
      INCLUDE 'fpvm3.h'
      character(len=30)  :: myname = 'ms'
      character(len=30)  :: myparnt= 'mm'
      character(len=100) :: c_aux
      character(len=70)  :: envstring
      character(len=100), EXTERNAL :: dir 
      integer :: ikill, iter
      integer :: info
      integer :: msg
      integer :: nints
      integer, dimension(0:32) :: ints
      real(8) :: dne
      real(8) :: work
      real(8) :: t, radat1d(nlines*7)
      external fbuff

      iunit = 1
      iReadFlag = 1

      envstring="STDOUT"
      call myinit( myname, myparnt, envstring)
!     write(1,*) 'aft myinit'
      call fbuff
      call x_rtids_p( 80, myname )
!     write(1,*) 'aft x_rtids_p'
      call fbuff

      msg = 81
      call ms_ridat_mm(msg)
      envstring="WORKDIR"
      call envarf(envstring,workdir)
      call AtomicKinetics(1)

! --------- Begin Iterative Secton -------------
      write(1,*) 'element=',element
      do
         msg = 2
         call x_rflg_p(msg,istat)
         if(istat.eq.0) exit 
         iter = 1
         ieinstein = 1
         msg = 82
         do 
            if(iter.eq.maxiter)then
               write(iunit,*) 'Error ms.f90:iter = maxiter'; call fbuff()
               stop
            endif
            iter=iter+1
            msg=msg+1
            call ms_rdat_mm(msg,iter)
            call atomickinetics(0)
            msg=msg+1
            call ms_sdat_mm(msg) 
            msg=msg+1
            call x_rflg_p(msg,istat)
            if(istat.eq.1) exit 
         enddo
     
         call spectra 
         call emissivity
         call opacity
         call ms_srt_mm(1)
! Printing 
         write(1,*) '#D Te_eV=',src_Te_eV,' Na=',src_rNa
         write(1,*) '#D Ne=',src_rNe
         write(1,*) '#D Zbar=',src_zbar
!        do i=1,src_ntotal
!           write(1,*) '#fp ',i,' ',src_fp(i)
!        enddo
    
      enddo
! ------ Printing Time Independent Pop.
!
!     IF(itd.EQ.0)THEN  !Steady state
!       c_aux = 'DIAG'//element//'/tindepfp'
!       c_aux = dir(workdir,c_aux)
!       OPEN(28,file=c_aux)
!       WRITE(28,*) 'dne reported below is only this elements ', &
!    &                  'contribution to the total-dne'
!           WRITE(28,*)  'Te_eV=',src_Te_eV, ' Ne=',src_rNe
!           WRITE(28,*)  'Zbar=',src_zbar
!           WRITE(28,*)
!           DO i=1,src_ntotal
!              WRITE(28,*) i,src_fp(i)
!           ENDDO
!           CLOSE(28)
!        ENDIF
!       
!     close(iunit)

! --------- End user program -------- 

      CALL pvmfexit(info) 
      END

!**********************************************************************
      SUBROUTINE AtomicKinetics(iset)
!**********************************************************************
      USE mac_code
      USE gms
!     OPEN(2,file='/home/manolo/SRMAC/foo')
!       WRITE(1,*) 'me=',me,' nproc=',nproc
!       WRITE(1,*) 'workdir=',workdir
!       WRITE(1,*) 'Te_eV =', src_Te_eV
!       WRITE(1,*) 'rNa   =', src_rNa
!       WRITE(1,*) 'zbar  =', src_zbar
!       WRITE(1,*) 'rNe   =', src_rNe
!       WRITE(1,*) 'iter  =', isrc_iter
!       WRITE(1,*) 'iset  =', iset
!     CLOSE(2)

      IF(iset.EQ. 1) THEN
        WRITE(iunit,*) 'initializing mac'
        call initialize_mac
        write(iunit,*) 'after init mac'
      ELSE
        IF(itd .EQ. 0) THEN  ! Time Indep.
!         WRITE(iunit,*) 'Time Independent Calc'
!         WRITE(iunit,*) 'iReadFlag=',iReadFlag
          CALL mac
        ELSE  ! Time Dep.
          WRITE(iunit,*) 'Time Dependent Calc'
          WRITE(iunit,*) 'iReadFlag=',iReadFlag
          CALL mac
        ENDIF
      ENDIF
      END
! **********************************************************************
      subroutine emissivity
! **********************************************************************
! radat(i,3)=A(j,k) 
! radat(i,2)=DeltaE(eV)
! radat(i,7)=Upper global index 
      use gms
      use communicate
      implicit real*8 (a-h,o-z)
      integer :: ii, jj
      real(8), dimension(specgrid) :: sum  
! ---------------------------------
! -- Diagnostics
      integer           :: iagain
      integer           :: unitnum
      character(len=30) :: ctemp
      character(len=3)  :: cagain
      character(len=70) :: envstring
      data iagain /0/
      save iagain
      external fbuff
      iagain=iagain+1
      write(cagain,'(i3.3)') iagain
! ---------------------------------
! ---------------------------------------------------------------
! -- Diagnostics
!        envstring = 'STDOUT'
!        unitnum=87
!        ctemp = 'ems'//trim(adjustl(get_cmytid()))//"."//cagain
!        call openf(unitnum,ctemp,envstring)
! ---------------------------------------------------------------

      PLANCKEVS=4.1356e-15    !eV*sec
      PLANCKERGS=6.626075e-27 !ergs*sec

      Pi= ACOS(-1.0)
      sum=0.0d0
      do kk=1,lineshape%ntrans
        i=lineshape%itrans(kk)
        j=AINT(radat(i,7)) 
        aux=src_rNa*src_fp(j)*radat(i,3)*(radat(i,2)/(4.0*Pi))
!       write(1,*) 'emissivity: radat trans #',i,' up=',j,' low=',AINT(radat(i,6))
        DO isgrid=1,gsize
          SUM(isgrid)=aux*(6.626075e-27/4.1356e-15)*&
     &     lineshape%prfdat(kk)%prf(isgrid)%y+SUM(isgrid)
        ENDDO
      ENDDO 
      DO isgrid=1,gsize
        j_emissivity(isgrid)=SUM(isgrid)
!       WRITE(87,*) j_emissivity(isgrid)
      ENDDO
!     CLOSE(87)
      end subroutine emissivity
! **********************************************************************
      SUBROUTINE OPACITY
! **********************************************************************
! radat(i,2)=DeltaE(eV) 
! radat(icnt,4)=B(k,i) stimulated emission
! radat(icnt,5)=B(i,k) photo absorption
      use gms
      use communicate
      implicit real*8 (a-h,o-z)
      integer :: ii,jj
      real(8), dimension(specgrid) :: sum
! ---------------------------------
! -- Diagnostics
      integer           :: iagain
      integer           :: unitnum
      real(8)           :: prtaux   ! LTE nj/ni
      character(len=30) :: ctemp
      character(len=3)  :: cagain
      character(len=70) :: envstring
      data iagain /0/
      save iagain
      external fbuff
      iagain=iagain+1
      write(cagain,'(i3.3)') iagain
! ---------------------------------
! ---------------------------------------------------------------
! -- Diagnostics
!        envstring = 'STDOUT'
!        unitnum=89
!        ctemp = 'opc'//trim(adjustl(get_cmytid()))//"."//cagain
!        call openf(unitnum,ctemp,envstring)
! ---------------------------------------------------------------
      PLANCKEVS=4.1356e-15    !eV*sec
      PLANCKERGS=6.626075e-27 !ergs*sec
      Pi= ACOS(-1.0)
      sum=0.0d0

      if(element.eq.'LI')then
        write(1,'(a3,a2,2(1x,1pe14.7))') '#ML',element,src_rNa/0.43d0, src_Te_eV
      else
        write(1,'(a3,a2,2(1x,1pe14.7))') '#ML',element,src_rNa/0.57d0, src_Te_eV
      endif

      write(1,901) element,src_rNa, src_Te_eV, src_zbar, src_rNe
      do kk=1,lineshape%ntrans
         j=lineshape%itrans(kk) 
         k=radat(j,7) ! upper global index     
         i=radat(j,6) ! lower global index 
         aux=src_rNa*(src_fp(i)*radat(j,5)-src_fp(k)*radat(j,4))*&
     &       (radat(j,2)*PLANCKERGS/4.0*Pi)          
!        write(1,*) 'opacity: radat trans #',j,' up=',k,' low=',i 

         prtaux=(radat(j,9)/radat(j,8))*exp(-radat(j,2)/src_Te_eV)  !LTE nj/ni
!        write(1,900) element,((src_fp(k)/src_fp(i))/prtaux),src_fp(i),src_fp(k),(radat(j,iprt),iprt=1,9)
         
         DO isgrid=1,gsize
           SUM(isgrid)=aux*lineshape%prfdat(kk)%prf(isgrid)%y+SUM(isgrid)
         ENDDO 
      enddo
      DO isgrid=1,gsize
        k_opacity(isgrid)=SUM(isgrid)
!       WRITE(89,*) k_opacity(isgrid)
!       call fbuff()
      ENDDO
!     CLOSE(89) 
900   format('#',a2,1x,8(1x,1e14.7),4(1x,f4.0))
901   format('#',a2,'NaTeZNe',1x,10(1x,1pe14.7))
      END SUBROUTINE OPACITY
! **********************************************************************
      subroutine spectra
! **********************************************************************
      use gms
      use communicate
      implicit none
      character(len=7)  :: sspecfile(0:2) 
!     character(len=70) :: envstring
      integer :: kk
      integer :: nn
      integer :: i
      integer :: ict
      integer :: ionce
! ---------------------------------
! -- Diagnostics
      integer           :: iagain
      integer           :: unitnum
      integer           :: idg
      integer           :: iflg,dinx !debugging
      character(len=30) :: ctemp
      character(len=3)  :: cagain
      character(len=70) :: envstring
      real(8), dimension(0:10) :: diag_rNe,idiag_rNe
      data iagain /0/
      save iagain
      data ionce /0/
      save ionce
      external fbuff     
      
      iagain=iagain+1
      write(cagain,'(i3.3)') iagain
! ---------------------------------
! ---------------------------------------------------------------
! -- Diagnostics
!        envstring = 'STDOUT'
!        unitnum=87
!        ctemp = 'spd'//trim(adjustl(get_cmytid()))//"."//cagain
!        call openf(unitnum,ctemp,envstring)
! ---------------------------------------------------------------
      if(ionce.eq.0)then
         call load
         ionce=1
      endif
!     write(1,*) 'after load'
      call fbuff()
      ict=0
      do kk=1, prfnames%nfnames   ! indx prf dat files
         do nn=1,trns(kk)%ntrans  ! indx trans in each dat file 
            ict=ict+1
            lineshape%ntrans=ict
            lineshape%itrans(ict)=trns(kk)%itrans(nn)
            call get_profile(kk,src_rNe,iflg)
            lineshape%prfdat(ict)=profile
!           write(87,*) '#',' rNe=',src_rNe
!           write(87,*) '#','kk=',kk,' iflg=',iflg
!           write(87,*) '#  dne=',lineshape%prfdat(ict)%dne
!           do idg=1,lineshape%prfdat(ict)%npts
!              write(87,*)lineshape%prfdat(ict)%prf(idg)%x, &
!    &                    lineshape%prfdat(ict)%prf(idg)%y  
!           enddo
!           write(87,*)
!           write(1,*)'transition=',lineshape%itrans(ict),&
!    &     ' dat file=',kk
        enddo                  
      enddo
!     close(87)
      end subroutine spectra
! *******************************************************
      subroutine get_profile(kk,dne,iflg)
! *******************************************************
! -- Li Stark Broadened Profiles  
!
!     kk : (input)    Index for profile data file
!     dne: (input)    Electron density
!     profile: (output-global gms) Then interpolated data file 
!
! *******************************************************
      use gms
      integer :: ix        ! aux index
      integer :: kk        ! profl. dat-file index 
      integer :: pkk       ! Prev. profl. dat-file index
      integer :: j         ! Density index for prof. dat-file
      integer :: iflg      ! Specifies class of location of dne in p-dat-file
      real(8) :: dne       ! Electron denisty
      real(8) :: pdne      ! Prev. electron density
      real(8), dimension(specgrid) :: slope
      real(8), dimension(specgrid) :: deltaprof
      data pkk /-1/
      data pdne /-1.0/
!      if((pkk.eq.kk).and.(dne.eq.pdne)) return
      iflg=0           ! Case prof_min < dne < prof_max
! -- Checking if dne profile already exits in prof-dat file 
      do j=1,trns(kk)%ntrans
         if(dne.eq.(trns(kk)%prfdat(j)%dne))then
            iflg = 1   ! Case dne = profile density
            goto 200
         endif
      enddo
! -- Checking if dne is less than lowest stored profile
      if(dne .lt. (trns(kk)%prfdat(1)%dne)) then
         iflg = 2      ! Case dne < profile density
         goto 200
      endif
! -- Checking if dne is greater than the largest stored profile
      ix = trns(kk)%nprof
!     write(1,*) 'ix=',ix
      if(dne.gt.(trns(kk)%prfdat(ix)%dne)) then
         iflg = 3      
         goto 200
      endif  
! -- 
200   continue
      if(iflg.eq.0) then 
        do j=1,trns(kk)%nprof-1
           if(dne.gt.(trns(kk)%prfdat(j)%dne)) then
              jupr=j+1
              jlwr=j
           endif
        enddo
      endif

      if(iflg.eq.1) then
        profile=trns(kk)%prfdat(j)
      elseif(iflg.eq.2) then
        profile=trns(kk)%prfdat(1)
      elseif(iflg.eq.3) then
        profile=trns(kk)%prfdat(ix)
      elseif(iflg.eq.0) then
        rlogdne = LOG10(dne)
        rlognelow = LOG10(trns(kk)%prfdat(jlwr)%dne)
        rlogneup  = LOG10(trns(kk)%prfdat(jupr)%dne)
        deltalogne = rlogneup-rlognelow
        deltaprof(:) = LOG10(trns(kk)%prfdat(jupr)%prf(:)%y)-&
     &                 LOG10(trns(kk)%prfdat(jlwr)%prf(:)%y)
        slope = deltaprof/deltalogne
        profile%prf(:)%y = LOG10(trns(kk)%prfdat(jupr)%prf(:)%y)+&
     &                     slope(:)*(rlogdne-rlogneup)
        profile%prf(:)%y= 10**profile%prf(:)%y
        profile%dne = trns(kk)%prfdat(jupr)%dne
      else
        print*,'error: else trap taken'
        stop
      endif
      pkk=kk
      pdne=dne
      end subroutine get_profile
    
! *******************************************************
      subroutine load
! ******************************************************* 
! Purpose: Read and loads trns data structure containing 
!     p: (output)   Structure of type profs containing profiles
! *******************************************************
      use gms
      use communicate
      implicit none
      integer :: i, ii ! profile density index
      integer :: j     ! profile mesh index
      integer :: k     ! dat file index
      integer :: nn
      integer :: ipk
      character(len=1)   :: chr
      character(len=30)  :: proflist
      character(len=70)  :: envstring
!     type(profdat) :: profile
! ---------------------------------
! -- Diagnostics
      integer           :: iagain
      integer           :: unitnum
      character(len=30) :: ctemp
      character(len=3)  :: cagain
      data iagain /0/
      save iagain
      external fbuff
      iagain=iagain+1
      write(cagain,'(i3.3)') iagain
! ---------------------------------
! ---------------------------------------------------------------
! -- Diagnostics
!        envstring = 'STDOUT'
!        unitnum=89
!        ctemp = 'load'//trim(adjustl(get_cmytid()))//"."//cagain
!        call openf(unitnum,ctemp,envstring)
! ---------------------------------------------------------------
      envstring = 'WORKDIR'
      unitnum=10
      proflist='newprof/proflist'//element//'.in'
      call openf(unitnum,proflist,envstring)
      ipk=0
      do  
         ipk=ipk+1
         read(10,*,end=80) prfnames%fname(ipk)
      enddo
80    continue
      prfnames%nfnames=ipk-1
      close(10)
      envstring = 'WORKDIR'
      unitnum=10
      do k=1,prfnames%nfnames
        call openf(unitnum,prfnames%fname(k),envstring) 
        read(10,*) chr, trns(k)%ntrans
        read(10,*) chr, (trns(k)%itrans(nn),nn=1,trns(k)%ntrans)
        print*,'ntrans=',trns(k)%ntrans,' itrans=',trns(k)%itrans
        i=0
        do
          i=i+1
          read(10,*,end=180) chr, trns(k)%prfdat(i)%dne
!         write(89,*) trns(k)%prfdat(i)%dne,' ',i,k
          call fbuff()
          read(10,*)
          read(10,*)
          read(10,*) chr, trns(k)%prfdat(i)%npts
!         write(1,*) chr, trns(k)%prfdat(i)%npts
          read(10,*)
          do j=1,trns(k)%prfdat(i)%npts
             read(10,*) trns(k)%prfdat(i)%prf(j)
!            write(1,*) trns(k)%prfdat(i)%prf(j)
          enddo
          gsize = trns(k)%prfdat(i)%npts  ! all profiles for a particular
        enddo                             ! element must have the same grid
180     continue 
        close(10)
        trns(k)%nprof=i-1
      enddo
      do k=1,prfnames%nfnames
        i=0
        do
          i=i+1
          if(i.gt.trns(k)%nprof)goto 190
          profile = trns(k)%prfdat(i)
!         write(89,*) '#',profile%dne
!         write(89,*) '#',profile%npts
          specwavlen=profile%prf%x
!         do j=1,profile%npts
!            write(89,*) profile%prf(j)
!            call fbuff()
!         enddo
        enddo
190   continue
      enddo
!     close(89)
      end subroutine load
