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