! Version Steady State 
! ------------------------------------------------------------------ !
      PROGRAM mm 
! ------------------------------------------------------------------ !
      USE stopwatch
      INCLUDE 'fpvm3.h'
      INTEGER :: istat        ! Spawn's kill notice (last calc)
      INTEGER mme, mptid, mmytid, mmsgtype, mmtid, mnproc, mtids(0:32)  
      INTEGER   igo, oldval
      INTEGER :: nin, nout
      REAL(8), DIMENSION(100) :: plistin
      REAL(8), DIMENSION(100) :: plistout
      CHARACTER*30 mycomp
      character*70 workdir
      call starttime
      CALL pvmfmytid( mmytid )
      CALL pvmfparent( mptid )  

      mmsgtype  = 1
      CALL pvmfrecv( mptid, mmsgtype, info )
      CALL pvmfunpack( INTEGER4, mnproc, 1, 1, info )
      CALL pvmfunpack( INTEGER4, mtids, mnproc, 1, info )
      CALL pvmfunpack( INTEGER4, nin, 1, 1, info )
      CALL pvmfunpack( INTEGER4, nout, 1, 1, info )
      CALL pvmfunpack( STRING, mycomp, 30, 1, info)
      CALL pvmfunpack( STRING, workdir, 70, 1, info ) 

!     Determine which slave I am (0 -- nproc-1)
      DO  i=0, mnproc-1
         IF( mtids(i) .EQ. mmytid ) mme = i
      ENDDO

      mmsgtype = 2
      call pvmfrecv( mptid, mmsgtype, info )
      call pvmfunpack( INTEGER4, igo, 1, 1, info )

      CALL subq1(mme,mycomp,workdir)

      call pvmfinitsend( PVMDEFAULT, info )
      call pvmfpack( INTEGER4, mme, 1, 1, info )
      call pvmfsend( mptid, mmsgtype, info )

      mmsgtype = mmsgtype + 1
      do
        call pvmfrecv( mptid, mmsgtype, info)
        call pvmfunpack( REAL8, plistin, nin, 1, info)
        call pvmfunpack( INTEGER4, istat, 1, 1, info)
        if(istat.eq.1) exit

          CALL subq2(nin,plistin,nout,plistout) 

        call pvmfinitsend( PVMDEFAULT, info )
        call pvmfpack( INTEGER4, mme, 1, 1, info)
        call pvmfpack( REAL8, plistout, nout, 1, info) 
        call pvmfsend( mptid, mmsgtype, info) 
      enddo

      call stoptime

      END
! ----------------------------------------------------------------------- !
      SUBROUTINE subq2(nin,plistin,nout,plistout)
! ----------------------------------------------------------------------- ! 
      REAL(8), DIMENSION(100) :: plistin
      REAL(8), DIMENSION(100) :: plistout

      CALL subq2ss(nin,plistin,nout,plistout)

      END
! ----------------------------------------------------------------------- !
      SUBROUTINE subq2ss(nin,plistin,nout,plistout)
! ----------------------------------------------------------------------- !
! Self-consistent interation section:
! ---------------------------------------------
      USE subs_time
      USE radtrans 
      USE stopwatch
      USE com_srmacm
      INCLUDE 'fpvm3.h'
      REAL(8), DIMENSION(100) :: plistin
      REAL(8), DIMENSION(100) :: plistout
      EXTERNAL cflusher
      ikill = 0         ! flag for killing slaves at final iter.
      iter = 1
      ilast = 0
      t =  plistin(1)
      dna = plistin(2)

      dne = 0.0d0
      do i=0,nproc
         dne = dne+(dna*fracspec(i)) * zbarinit(i)
         fdna(i)=dna*fracspec(i)
      enddo
      DO WHILE(iter .LE. maxiter) 
         WRITE(istdout,*)'  '
         WRITE(istdout,*)' --------------------------------------------'
         WRITE(istdout,*)'  '
         iter=iter+1
         msgtype = msgtype + 1
         CALL pvmfinitsend( PVMDEFAULT, info )
         CALL pvmfpack( REAL8, t, 1, 1, info )
         CALL pvmfpack( REAL8, dne, 1, 1, info ) 
         CALL pvmfpack( REAL8, fdna, nproc, 1, info )
         CALL pvmfmcast( nproc, tids, msgtype, info )

         WRITE(istdout,'(a3,1pe14.7,a8,e14.7)') 'Te=',t,' Ne_int=',dne
         WRITE(istdout,'(a3,1pe14.7)') 'Na=',dna
 
! -- Wait for results from nodes 
         msgtype  = msgtype + 1 
         DO i=1,nproc 
            CALL pvmfrecv( -1, msgtype, info ) 
            CALL pvmfunpack( INTEGER4, who, 1, 1, info )
            CALL pvmfunpack( INTEGER4,iconvrg(who),1,1,info )
            CALL pvmfunpack( INTEGER4,iconsum(who),1,1,info )
            CALL pvmfunpack( INTEGER4,nho(who),1,1,info )
            CALL pvmfunpack( REAL8, zbars(who), 1, 1, info )

            WRITE(istdout,*) 'Proc.#=',who,' ',elmt((who*2+1):(who*2+2)),  &
     &              ' zbar=',zbars(who)
            WRITE(istdout,*) ' Number of levels converged=',iconsum(who),  &
     &              ' out of',nho(who)
         ENDDO

         isumconvrg=0
         DO i=0,nproc-1  ! iconvrg = 0-> not converg.; 1-> has converg.
            IF( iconvrg(i).EQ. 1 ) isumconvrg = isumconvrg + 1
         ENDDO
         WRITE(istdout,*)' Iter=',iter-1,' #Species converted:',isumconvrg
         IF(ilast .EQ. 0) THEN
           IF(isumconvrg .EQ. nproc) ilast = 1
         ENDIF
         WRITE(istdout,*)'ilast=',ilast

! ---- CHECKING FOR LAST STEP OF A SELF-CONSITANT ITERATION: ilast=2

         IF( ilast .EQ. 2 ) THEN 
           IF( itd .EQ. 0) THEN ! Time Indep. Case  (Steady State)
             WRITE(istdout,*)' Final Interation Reached Steady State'
! --------  Sending Slave Termination Notice Time Indep. (Steady State)
              idiag = 1
              CALL send_killprt(nproc,tids,msgtype,0,idiag)
              IF(idiag.EQ.1) THEN 
                call cflusher ()
                CALL get_specinfo
                CALL radtransport(workdir)
              ENDIF
              GOTO 200
           ENDIF !(itd.eq.1).or.(itd.eq.0)
         ELSE ! (ilast .ne. 2)
            WRITE(istdout,*)'Slave Self-Consistent Continue Notice'
! ---       Sending Slave Self-Consistent Continue Notice 
            CALL send_killprt(nproc,tids,msgtype,0,0)
            dne = 0.0
            DO i=0,nproc-1
               dne= dne+ dna*fracspec(i)*zbars(i)
            ENDDO
            WRITE(istdout,'(a4,1pe14.7)') ' Ne=',dne 
            call cflusher ()
         ENDIF 
         IF(ilast .EQ. 1) ilast = 2
      ENDDO !do while
200   CONTINUE
      if(itd.eq.1) then
        CLOSE(25)
        CLOSE(26)
        CLOSE(21)
      endif
      WRITE(istdout,'(a4,1pe14.7)') ' Ne=',dne
! --------- End user program -------- 

!     program finished leave PVM before exiting 
!     CALL pvmfexit(info) 
      WRITE(istdout,*) 'ending the master1'
      END

! --------------------------------------------------- ! 
      SUBROUTINE shutdown( nproc, tids )
! --------------------------------------------------- ! 
! Purpose: Kill all tasks I spawned and then myself
! --------------------------------------------------- ! 
      INTEGER nproc, tids(*)

      DO 10 i=0, nproc
         CALL pvmfkill( tids(i), info )
  10  CONTINUE
      CALL pvmfexit( info )
      RETURN
      END

! ----------------------------------------------------------------------- !
      SUBROUTINE subq1(me,mycomp,tranworkdir)
! ----------------------------------------------------------------------- !
! Purpose: Reading in of data for multi-component calculations and
!          spawning of atomic kinetic elements
!
      USE subs_time
      USE radtrans 
      USE stopwatch
      USE com_srmacm
      INCLUDE 'fpvm3.h'
      EXTERNAL cflusher

     integer me
     character(len=30) mycomp
     character(len=20) ctemp
     character(len=100) cffoo
     character(70) tranworkdir

     workdir=tranworkdir
     mme = me
     WRITE(cme,'(i2.2)') mme
     istdout = 1
      i = 0
      ilast = 0
      nproc = 0
      sumfrac=0.0   
      dne = 0.0
      itimesteps = 0
      ctemp = 'srmacm.in'
      cffoo = dir(workdir,ctemp)
      OPEN(unit=20,file=cffoo)
!     OPEN(unit=20,file='/home/manolo/SRMAC/srmacm.in')
!     READ(20,*)
!     READ(20,'(a70)') workdir
      READ(20,*)
      ctemp = 'stdmaster.out'//cme
      cffoo = dir(workdir,ctemp)
      OPEN(unit=istdout,file=cffoo)
      WRITE(istdout,'(a70)') workdir
      WRITE(istdout,'(i2.2)') mme
      if(itd.eq.1)then 
        ctemp = 'srmacm.out'//cme
        cffoo = dir(workdir,ctemp)
        OPEN(unit=26,file=cffoo) 
        ctemp = 'srmac.col.out'//cme
        cffoo = dir(workdir,ctemp)
        OPEN(unit=21,file=cffoo)
        WRITE(21,'(a45)')'# itimesteps  dna  dne  te  Time(nsec)  ztot '
      endif 
      WRITE(istdout,*)'workdir=',workdir
      READ(20,'(I5)') maxiter
      WRITE(istdout,*) 'maxiter=',maxiter
      READ(20,'(I5)') idiag
      WRITE(istdout,*) 'idiag=',idiag
      READ(20,'(E14.7)') eps
      WRITE(istdout,*) 'eps=',eps
      READ(20,'(E14.7)') t
      WRITE(istdout,*) 't=',t
      READ(20,'(E14.7)') dna
      WRITE(istdout,*) 'dna=',dna
      READ(20,'(E14.7)') hnu1
      WRITE(istdout,*) 'hnu1=',hnu1
      READ(20,'(E14.7)') hnu2
      WRITE(istdout,*) 'hnu2=',hnu2
      READ(20,'(E14.7)') fwhm
      WRITE(istdout,*) 'fwhm=',fwhm
      READ(20,'(1x,i1)') itd 
      WRITE(istdout,*) 'itd=',itd 
10    CONTINUE      
         READ(20,*,END=20)
         READ(20,'(1x,a2)',END=20) auxelmt
         WRITE(istdout,*)'auxelmt= ',auxelmt
         elmt((i*2+1):(i*2+2)) = auxelmt
         READ(20,'(E14.7)',END=20) zbarinit(i)
         WRITE(istdout,*) 'zbarinit(',i,')=',zbarinit(i)
         READ(20,'(E14.7)',END=20) fracspec(i)
         WRITE(istdout,*) 'fracspec(',i,')=',fracspec(i)
!        dne = dne+(dna*fracspec(i)) * zbarinit(i)
!        fdna(i)=dna*fracspec(i)
!        sumfrac = sumfrac+fracspec(i)
         nproc = nproc + 1
         i=i+1
      GOTO 10
20    CONTINUE
      CLOSE(20)

! -- Time Dependent Input

      IF(itd .EQ. 1)THEN

          iprint = 1
          nnn = 1
          cffoo = dir(workdir,'srmacm_td.in')
          OPEN(unit=20,file=cffoo)
          READ(20,*)
          READ(20,'(1x,1pe14.7)') tinitial
          WRITE(istdout,*) 'tinitial=',tinitial 
          READ(20,'(1x,1pe14.7)') tfinal
          WRITE(istdout,*) 'tfinal=',tfinal
          READ(20,'(1x,i5)') nprints
          WRITE(istdout,*) 'nprints=',nprints 
          READ(20,'(1x,i5)') ntsteps
          WRITE(istdout,*)'ntsteps=',ntsteps
          READ(20,'(1x,1pe15.8)') dne
          WRITE(istdout,*)
          WRITE(istdout,*) '!!!!  WARN: CHECK DNE VALUE  !!!!'
          WRITE(istdout,'(6x,a5,1pe14.7)') ' dne= ',dne
          WRITE(istdout,*) '!!!!  WARN: CHECK DNE VALUE  !!!!'
          WRITE(istdout,*)'   '
          WRITE(istdout,*)'Reading srmacm_td.in complete'
          WRITE(istdout,*)'   '
          CLOSE(20)

          steptp=(tfinal-tinitial)/(nprints-1)
          DO i=1,nprints
             timep(i)=tinitial+(i-1)*steptp
             WRITE(istdout,*)'timep(',i,')=',timep(i)
          ENDDO

          cffoo = dir(workdir,'srmacm_thist.in')
          OPEN(unit=22,file=cffoo) ! read in sub timehist
          CALL timehist(tinitial,tfinal)
          WRITE(istdout,*)'    '
          WRITE(istdout,*)'Reading Time History Complete'
          WRITE(istdout,*)'    '
          CLOSE(22)

          deltat=(timep(iprint+1)-timep(iprint))/ntsteps
          time=timep(iprint)+nnn*deltat
      ENDIF
      WRITE(istdout,*)'deltat =',deltat
      WRITE(istdout,*)'dne    =',dne
      WRITE(istdout,*)'time   =',time
 
      call pvmfmytid( mytid )

      nodename = 'ms'
!     arch = '*'
!      arch = mycomp
!     CALL pvmfspawn( nodename, PVMHOST, mycomp, nproc, tids, numt )

      arch = '*'
      call pvmfspawn( nodename, PVMDEFAULT, arch, nproc, tids, numt )

!     Print out task IDs of spawned tasks and check for problems
      do i=0, nproc-1
         write(istdout,*) 'tid',i,tids(i)
      enddo
      IF( numt .LT. nproc ) THEN
         WRITE(istdout,*) 'trouble spawning ',nodename
         WRITE(istdout,*) ' Check tids for error code'
         CLOSE(istdout)
         CALL shutdown( numt, tids )
      ENDIF

! --- Broadcast tid information to all node programs 
      CALL pvmfinitsend( PVMDEFAULT, info )
      CALL pvmfpack( INTEGER4, nproc, 1, 1, info )
      CALL pvmfpack( INTEGER4, tids, nproc, 1, info )

      msgtype  = 1 
      CALL pvmfmcast( nproc, tids, msgtype, info )

      CALL pvmfinitsend( PVMDEFAULT, info )
      CALL pvmfpack( INTEGER4, idiag, 1, 1, info )
      CALL pvmfpack( INTEGER4, maxiter, 1, 1, info )
!     CALL pvmfpack( REAL8, fdna, nproc, 1, info )
      CALL pvmfpack( REAL8, hnu1, 1, 1, info ) ! see p.78 PVM manual
      CALL pvmfpack( REAL8, hnu2, 1, 1, info )
      CALL pvmfpack( REAL8, fwhm, 1, 1, info )
      CALL pvmfpack( REAL8,  eps, 1, 1, info )
      CALL pvmfpack( STRING, elmt, 66, 1, info )
      CALL pvmfpack( STRING, workdir,70, 1, info)
      CALL pvmfpack( INTEGER4, itd, 1, 1, info ) 
      IF(itd.EQ.1)THEN
         CALL pvmfpack( REAL8, time, 1, 1, info )
         ctemp = 'td_data.out'//cme
         cffoo = dir(workdir,ctemp)
         OPEN(25,file=cffoo)
      ENDIF
      msgtype  = 2 
      CALL pvmfmcast( nproc, tids, msgtype, info )
      END
! -------------------------------------------------------------------- !
      SUBROUTINE get_specinfo
! -------------------------------------------------------------------- ! 
      USE com_srmacm
      USE com_radtrans
      USE stopwatch
      INCLUDE 'fpvm3.h'
      CHARACTER(len=100) cffoo
      CHARACTER(Len=3) cagain
      CHARACTER(LEN=30) ctemp
      INTEGER :: nsize, itest,iagain
      INTEGER :: nlev(0:32)
      REAL(8), DIMENSION(MZ) :: yparab ! used to induce artifical grad.
      REAL(8), DIMENSION(SG) :: array
      REAL(8), DIMENSION(SG) :: emis1d
      REAL(8), DIMENSION(SG) :: opac1d
      REAL(8) :: p,h,k,bdry
      REAL(8), DIMENSION(SG) :: xaux1,xaux2,opc1,opc2,ems1,ems2
      EXTERNAL cflusher
      SAVE iagain
      DATA iagain /0/
      iagain = iagain + 1
      h=thick/2.0  ! Parbola vertex (h,k)
      k=1.0
      p=-(h**2)/(4.0*k)
         msgtype = msgtype + 1
         emis1d=0.0d0
         opac1d=0.0d0
         array=0.0d0
         DO i=0,nproc-1
           CALL pvmfrecv( tids(i), msgtype,info)
           CALL pvmfunpack( INTEGER4, nsize, 1, 1, info)
           CALL pvmfunpack( REAL8, xspec, nsize, 1, info )
           CALL pvmfunpack( REAL8, array, nsize, 1, info )
           emis1d(:)=array(:)+emis1d(:)
           CALL pvmfunpack( REAL8, array, nsize, 1, info )
           opac1d(:)=array(:)+opac1d(:)
         ENDDO
 
        bdry=thick/MZ            ! Length between zones       
        DO i=0,(MZ-1)            ! Filling zone center data
           opacity%zc(i+1)=i*bdry+bdry/2.0
           emissivity%zc(i+1)=i*bdry+bdry/2.0
        ENDDO

        goto 1000
        cffoo = dir(workdir,'opc3.98e16_0.8')
        OPEN(unit=30,file=cffoo)
        DO i=1,SG
           READ(30,*) xaux1(i), opc1(i)
        ENDDO
        CLOSE(30)
        cffoo = dir(workdir,'ems3.98e16_0.8')
        OPEN(unit=30,file=cffoo)
        DO i=1,SG
           READ(30,*) xaux1(i), ems1(i)
        ENDDO
        CLOSE(30)
        cffoo = dir(workdir,'opc1.58e17_0.8')
        OPEN(unit=31,file=cffoo)
        DO i=1,SG
           READ(31,*) xaux2(i), opc2(i)
        ENDDO
        CLOSE(31)
        cffoo = dir(workdir,'ems1.58e17_0.8')
        OPEN(unit=31,file=cffoo)
        DO i=1,SG
           READ(31,*) xaux2(i), ems2(i)
        ENDDO
        CLOSE(31)
1000    continue
!HERE 
!       DO i=1,MZ
!          if((i.lt.0).or.(i.gt.10)) then
!            emissivity%a2d(:,i)=ems1(:)
!            opacity%a2d(:,i)=opc1(:)
!          else
!            emissivity%a2d(:,i)=ems2(:)
!            opacity%a2d(:,i)=opc2(:)
!          endif
!       ENDDO

        DO i=1,MZ   ! Filling emis. and opac. arrays
           emissivity%a2d(:,i)=emis1d(:)
           opacity%a2d(:,i)=opac1d(:)
        ENDDO

        DO i=0,MZ
           hbdry(i)=i*bdry
        ENDDO

        write(cagain,'(i3.3)') iagain
        ctemp = '/home/manolo/srmacrt/opac'//cagain
        OPEN(30,file = ctemp)
        ctemp = '/home/manolo/srmacrt/emis'//cagain
        OPEN(31,file = ctemp)
        DO i=1,SG
           WRITE(30,*) xspec(i), opac1d(i)
           WRITE(31,*) xspec(i), emis1d(i)
        ENDDO
        WRITE(30,*) 
        WRITE(31,*)

!        OPEN(30,file='/home/manolo/SRMAC/opac2d')
!        OPEN(31,file='/home/manolo/SRMAC/ems2d')
!        DO i=1,SG
!          WRITE(30,900) xspec(i),(opacity%a2d(i,j),j=1,MZ)
!          WRITE(31,900) xspec(i),(emissivity%a2d(i,j),j=1,MZ)
!        ENDDO
!        CLOSE(30)
!        CLOSE(31)
         write(istdout,*) 'leaving get_spec'
         call cflusher ()      
900   FORMAT(1x,e12.5,5000(2x,e12.5))
      END SUBROUTINE get_specinfo

! ************************************************************* 
      SUBROUTINE send_killprt(nproc,tids,msgtype,ikill,iprt)
! *************************************************************
      INCLUDE 'fpvm3.h'
      INTEGER info, ikill, iprt, nproc, msgtype
      INTEGER tids(0:32)
            msgtype = msgtype + 1
            print*,'MASTER msgtype=',msgtype
            CALL pvmfinitsend( PVMDEFAULT, info )
            CALL pvmfpack( INTEGER4,  ikill, 1, 1, info )
            CALL pvmfpack( INTEGER4, iprt, 1, 1, info )
            CALL pvmfmcast( nproc, tids, msgtype, info )
      END SUBROUTINE send_killprt

! ******************************************************************
      SUBROUTINE print26_21(itimesteps,iterprt,dna,dne,t,time,ztot)
! ******************************************************************
         IMPLICIT REAL*8 (a-h,o-z)
         itimesteps = itimesteps + 1
         Timensec = Time*1.0e9
         WRITE(26,*) '------------------------------------'
         WRITE(26,'(a7,i4)')      ' #=    ',itimesteps
         WRITE(26,'(a7,i4)')      ' iter= ',iterprt
         WRITE(26,'(a7,1pe14.7)') ' Na=   ',dna
         WRITE(26,'(a7,1pe14.7)') ' Ne=   ',dne
         WRITE(26,'(a7,f14.7)')   ' Te=   ',t
         WRITE(26,'(a7,f14.7)')   ' Time= ',Timensec
         WRITE(26,'(a7,f14.7)')   ' Ztot= ',ztot
         WRITE(21,'(i4,5(2x,e16.9))') itimesteps,dna,dne,t,Timensec,ztot
      END SUBROUTINE print26_21
!***********************************************************************
      function dir(string1,string2)
!***********************************************************************
      implicit none
      character(len=*) :: dir
      character(len=*),intent(in)  :: string1,string2
      dir=trim(adjustl(string1))//adjustl(string2)
      end function dir
