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 ------------- 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) 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 IMPLICIT REAL*8 (a-h,o-z) CHARACTER*7 sspecfile(0:2) REAL(8), DIMENSION(specgrid) :: sum PLANCKEVS=4.1356e-15 !eV*sec PLANCKERGS=6.626075e-27 !ergs*sec sspecfile(0)='sLi.ems' sspecfile(1)='sAg.ems' sspecfile(2)='sPt.ems' ! OPEN(unit=87,file='/home/manolo/SRMAC/'//sspecfile(me)) Pi= ACOS(-1.0) sum=0.0d0 DO i=1,iradat j=AINT(radat(i,7)) aux=src_rNa*src_fp(j)*radat(i,3)*(radat(i,2)/(4.0*Pi)) DO isgrid=1,specgrid SUM(isgrid)=aux*(6.626075e-27/4.1356e-15)*src_lineshape(i,isgrid)+SUM(isgrid) ENDDO ENDDO DO isgrid=1,specgrid j_emissivity(isgrid)=SUM(isgrid) ! WRITE(87,*) specwavlen(isgrid),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 IMPLICIT REAL*8 (a-h,o-z) CHARACTER*7 sspecfile(0:2) REAL(8), DIMENSION(specgrid) :: sum PLANCKEVS=4.1356e-15 !eV*sec PLANCKERGS=6.626075e-27 !ergs*sec sspecfile(0)='sLi.opc' sspecfile(1)='sAg.opc' sspecfile(2)='sPt.opc' ! OPEN(unit=89,file='/home/manolo/SRMAC/'//sspecfile(me)) Pi= ACOS(-1.0) sum=0.0d0 DO j=1,iradat 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(89,'(a1,3(1x,i4),2(2x,f14.7),4(2x,e12.5))') '#',j,k,i,radat(j,2),12398.0/radat(j,2),src_rNa,src_fp(i),src_fp(k),aux DO isgrid=1,specgrid SUM(isgrid)=aux*src_lineshape(j,isgrid)+SUM(isgrid) ENDDO ENDDO DO isgrid=1,specgrid k_opacity(isgrid)=SUM(isgrid) ! WRITE(89,*) specwavlen(isgrid),12398.0/specwavlen(isgrid),k_opacity(isgrid) ENDDO ! CLOSE(89) END SUBROUTINE OPACITY ! ********************************************************************** SUBROUTINE spectra ! ********************************************************************** USE gms use communicate IMPLICIT REAL*8 (a-h,o-z) CHARACTER*7 sspecfile(0:2) CHARACTER*4 charfoo REAL(8), DIMENSION(nlines) :: spec_aux1 REAL(8), DIMENSION(nlines) :: wavelen REAL(8), DIMENSION(nlines) :: avoigt REAL(8), DIMENSION(nlines) :: tint ! Non-Broaden line intensity REAL(8), DIMENSION(4001) :: profile character(len=70) :: envstring DATA iprtflag/0/ SAVE iprtflag ! --------------------------------- ! -- Diagnostics INTEGER :: unitnum CHARACTER(len=30) :: ctemp CHARACTER(len=3) :: cagain DATA iagain /0/ SAVE iagain iagain=iagain+1 write(cagain,'(i3.3)') iagain ! --------------------------------- if(me.ne.0)then ! --------------------------------------------------------------- ! -- Diagnostics envstring = 'WORKDIR' unitnum=89 ctemp = 'diag'//trim(adjustl(get_cmytid()))//"."//cagain call openf(unitnum,ctemp,envstring) ! --------------------------------------------------------------- endif hbar=6.5822e-16 sqrpi=SQRT(ACOS(-1.0)) spec_aux0=4.617e-5*SQRT(src_Te_eV/src_amass) sspecfile(0)='sspecLi' sspecfile(1)='sspecAg' sspecfile(2)='sspecPt' IF(iprtflag.EQ.1)THEN OPEN(unit=88,file='/home/manolo/mac/'//sspecfile(me)) ENDIF open(90,file='/tmp/mac/norm') DO i=1,iradat j=INT(radat(i,7),8) spec_aux1(i)=spec_aux0*radat(i,2) + fwhm tint(i)=src_rNa*src_fp(j)*radat(i,3)*radat(i,2) avoigt(i)= hbar*srdsum(j)/(2.0*spec_aux1(i)) ENDDO stephnu=(hnu2-hnu1)/4000.0 ilineshape=1 if(ilineshape.eq.0)then DO k=1,4001 hnup=hnu1+stephnu*(k-1) sspec=0.0d0 DO i=1,iradat aa=avoigt(i) vv=ABS(hnup-radat(i,2))/spec_aux1(i) hav=voigt(vv,aa) rint=tint(i)*hav/(spec_aux1(i)*sqrpi) sspec=sspec+rint ENDDO ! WRITE(88,*) HNUP,SSPEC wl=12398.0/hnup ! WRITE(88,*) WL,SSPEC specwavlen(k)=wl spec(k)=sspec ENDDO ! WRITE(88,*) else DO i=1,iradat IF(me.eq.0)then ! Li only if((i.eq.1).or.(i.eq.4).or.(i.eq.7).or.(i.eq.10.).or.(i.eq.13))then do k=1,4001 hnup=hnu1+stephnu*(k-1) levnum=i ! if(i.eq.1) then ! write(charfoo,'(i4.4)') k ! open(93,file='/home/manolo/SRMAC/hereb'//charfoo) ! write(93,*) levnum,' ',src_rNe,' ',hnup,' ',val ! close(93) ! endif call get_starkprofiles(levnum,src_rNe,hnup,val,workdir) src_lineshape(i,k)=val ! if(i.eq.1) then ! write(charfoo,'(i4.4)') k ! open(93,file='/home/manolo/SRMAC/herea'//charfoo) ! write(93,*) levnum,' ',src_rNe,' ',hnup,' ',val ! close(93) ! endif enddo else ! Li but not levels 1,4,7,10,13 aa=avoigt(i) DO k=1,4001 hnup=hnu1+stephnu*(k-1) specwavlen(k)=hnup ! specwavlen(k)=12398.0/hnup vv=ABS(hnup-radat(i,2))/spec_aux1(i) hav=voigt(vv,aa) src_lineshape(i,k)=hav/(spec_aux1(i)*sqrpi) ENDDO endif ELSE ! .NE. Li aa=avoigt(i) write(89,*) 'spec_aux0=',spec_aux0 write(89,*) 'radat(',i,'2)=',radat(i,2) write(89,*) 'avoigt(',i,')=',avoigt(i) write(89,*) 'spec_aux1(',i,')=',spec_aux1(i) DO k=1,4001 hnup=hnu1+stephnu*(k-1) specwavlen(k)=hnup ! specwavlen(k)=12398.0/hnup vv=ABS(hnup-radat(i,2))/spec_aux1(i) hav=voigt(vv,aa) src_lineshape(i,k)=hav/(spec_aux1(i)*sqrpi) profile(k)=src_lineshape(i,k) ENDDO aux=trap(hnu1,hnu2,4001,profile) write(90,*) i,' ',aux ENDIF ENDDO IF(iprtflag.EQ.1)THEN DO i=1,iradat DO k=1,4001 profile(k)=src_lineshape(i,k) ENDDO aux=trap(hnu1,hnu2,4001,profile) write(90,*) i,aux ENDDO DO k=1,specgrid WRITE(88,900) specwavlen(k),(src_lineshape(i,k), i=1,iradat) ENDDO iprtflag=0 CLOSE(88) close(89) ENDIF endif 900 FORMAT((1x,e12.5),500(2x,e12.5)) END SUBROUTINE SPECTRA ! *************************************************************************** REAL(8) FUNCTION TRAP(a,b,n,y) ! *************************************************************************** ! Integrator: Trapizoital Rule IMPLICIT REAL*8 (a-h,o-z) REAL(8),DIMENSION(5000) :: y sum=0.0d0 DO i=2,n-1 sum=2.0d0*y(i)+sum ENDDO trap=(sum+y(1)+y(n))*((b-a)/(2.0*n)) END FUNCTION TRAP ! **************************************************************************** REAL(8) FUNCTION VOIGT(x,y) ! **************************************************************************** IMPLICIT REAL*8 (a-h,o-z) SAVE ! *** routine computes the voigt function y/pi*integral *** ! *** - to + infinity of exp(-t*t)/(y*y+(x-t)*(x-t)) dt *** ! *** s.r. drayson j.q.s.r.t. vol. 16 , 611 (1976) *** DIMENSION b(22),ri(15),xn(15),yn(15),d0(25),d1(25),d2(25),d3(25), & &d4(25),hn(25),xx(3),hh(3),nby2(19),c(21) DATA b(1)/0./,b(2)/.7093602e-7/,xn/10.,9.,2*8.,7.,6.,5.,4.,7 & &*3./,yn/3*.6,.5,2*.4,4*.3,1.,.9,.8,2*.7/, & &h/.201/,xx/.5246476,1.65068,.7071068/, & &hh/.2562121,.2588268e-1,.2820948/,nby2/9.5,9.,8.5,8.,7.5,7.,6.5 & &,6.,5.5,5.,4.5,4.,3.5,3.,2.5,2.,1.5,1.,.5/,c/.7093602e-7, & &-.2518434e-6,.8566874e-6,-.2787638e-5,.866074e-5,-.2565551e-4,.7228775e-4, & &-.1933631e-3,.4899520e-3,-.1173267e-2,.2648762e-2,-.5623190e-2, & &.1119601e-1,-.2084976e-1,.3621573e-1,-.5851412e-1,.8770816e-1,-.121664, & &.15584,-.184,.2/ LOGICAL tru/.FALSE./ IF (tru) go to 104 ! *** region i. compute dawson's function at mesh points *** tru=.TRUE. DO 101 i=1,15 101 ri(i)=-i/2. DO 103 i=1,25 hn(i)=h*(i-.5) c0=4.*hn(i)*hn(i)/25.-2. DO 102 j=2,21 102 b(j+1)=c0*b(j)-b(j-1)+c(j) d0(i)=hn(i)*(b(22)-b(21))/5. d1(i)=1.-2.*hn(i)*d0(i) d2(i)=(hn(i)*d1(i)+d0(i))/ri(2) d3(i)=(hn(i)*d2(i)+d1(i))/ri(3) 103 d4(i)=(hn(i)*d3(i)+d2(i))/ri(4) 104 IF (x-5.) 105,112,112 105 IF (y-1.) 110,110,106 106 IF (x.GT.1.85*(3.6-y)) go to 112 ! *** region ii: continued fraction. compute number of terms needed. *** IF (y.LT.1.45) go to 107 i=y+y go to 108 107 i=11.*y 108 j=x+x+1.85 max=xn(j)*yn(i)+.46 min=min0(16,21-2*max) ! *** evaluate continued fraction *** uu=y vv=x DO 109 j=min,19 u=nby2(j)/(uu*uu+vv*vv) uu=y+u*uu 109 vv=x-u*vv voigt=uu/(uu*uu+vv*vv)/1.772454 RETURN 110 y2=y*y IF (x+y.GE.5.) go to 113 ! *** region i. compute dawson's function at x from taylor series. *** n=x/h dx=x-hn(n+1) u=(((d4(n+1)*dx+d3(n+1))*dx+d2(n+1))*dx+d1(n+1))*dx+d0(n+1) v=1.-2.*x*u ! *** taylor series expansion about y=0.0 *** vv=EXP(y2-x*x)*COS(2.*x*y)/1.128379-y*v uu=-y max=5.+(12.5-x)*.8*y DO 111 i=2,max,2 u=(x*v+u)/ri(i) v=(x*u+v)/ri(i+1) uu=-uu*y2 111 vv=vv+v*uu voigt=1.128379*vv RETURN 112 y2=y*y IF (y.LT.11.-.6875*x) go to 113 ! *** region iiib: 2-point gauss-hermite quadrature. *** u=x-xx(3) v=x+xx(3) voigt=y*(hh(3)/(y2+u*u)+hh(3)/(y2+v*v)) RETURN ! *** region iiia: 4-point gauss-hermite quadrature. *** 113 u=x-xx(1) v=x+xx(1) uu=x-xx(2) vv=x+xx(2) voigt=y*(hh(1)/(y2+u*u)+hh(1)/(y2+v*v)+hh(2)/(y2+uu*uu)+hh(2)/(y2+vv*vv)) RETURN END FUNCTION VOIGT ! ******************************************************* subroutine get_starkprofiles(lev,dne,eng,val,workdir) ! ******************************************************* ! ! -- Li Stark Broadened Profiles ! ! lev: (input) Level number ! dne: (input) Electron density ! eng: (input) Energy value to evaluate spectrum ! val: (output) Value of intensity of spectrum ! ! ******************************************************* IMPLICIT REAL*8 (a-h,o-z) PARAMETER(MZ=11,ML=200, MT=200, MF=2001, MA=10) TYPE profs REAL(8), DIMENSION(MF) :: x REAL(8), DIMENSION(MF) :: prof1 REAL(8), DIMENSION(MF) :: prof2 REAL(8), DIMENSION(MF) :: prof3 REAL(8), DIMENSION(MF) :: prof8 REAL(8), DIMENSION(MF) :: prof9 END TYPE profs TYPE(profs), TARGET :: p DATA pdne /-1.0/ DATA xp /0/ SAVE EXTERNAL fbuff if(pdne.ne.dne)then call load(dne,p,workdir) pdne=dne endif ! --- Converting indexies if(lev.eq.1) then lev = 1 elseif(lev.eq.4) then lev = 2 elseif(lev.eq.7) then lev = 3 elseif(lev.eq.10) then lev = 8 elseif(lev.eq.13) then lev = 9 endif if(lev .eq. 1)then call linlocate(eng,p%x,MF,iloc) call lininterp(iloc,p%x,p%prof1,MF,eng,val) elseif(lev.eq.2)then call linlocate(eng,p%x,MF,iloc) call lininterp(iloc,p%x,p%prof2,MF,eng,val) elseif(lev.eq.3)then call linlocate(eng,p%x,MF,iloc) call lininterp(iloc,p%x,p%prof3,MF,eng,val) elseif(lev.eq.8)then call linlocate(eng,p%x,MF,iloc) call lininterp(iloc,p%x,p%prof8,MF,eng,val) elseif(lev.eq.9)then call linlocate(eng,p%x,MF,iloc) call lininterp(iloc,p%x,p%prof9,MF,eng,val) else write(6,*) 'ERROR TRAP: get_starkprofiles ' endif END ! ******************************************************* SUBROUTINE load(dne,p,workdir) ! ******************************************************* ! ! dne: (input) Electron Density (cm^-3) ! p: (output) Structure of type profs containing profiles ! ! ******************************************************* IMPLICIT REAL*8 (a-h,o-z) PARAMETER(MZ=11,ML=200, MT=200, MF=2001, MA=10) DIMENSION dne_array(0:7),dne_array2(0:12) ! CHARACTER*34 fLi1(6), fLi2(11) CHARACTER*100 fLi1(6), fLi2(11) CHARACTER*100 c_aux CHARACTER*70 workdir CHARACTER(len=100), EXTERNAL :: dir TYPE profs REAL(8), DIMENSION(MF) :: x REAL(8), DIMENSION(MF) :: prof1 REAL(8), DIMENSION(MF) :: prof2 REAL(8), DIMENSION(MF) :: prof3 REAL(8), DIMENSION(MF) :: prof8 REAL(8), DIMENSION(MF) :: prof9 END TYPE profs TYPE(profs), TARGET :: p SAVE NPDF = MF dne_array(0)=0.0d0 dne_array(1)=1.0d17 dne_array(2)=2.0d17 dne_array(3)=5.0d17 dne_array(4)=1.0d18 dne_array(5)=2.0d18 dne_array(6)=5.0d18 dne_array(7)=1.0e28 c_aux = 'trn/Li2s2p.6' fLi1(1) = dir(workdir,c_aux) c_aux = 'trn/Li2s2p.5' fLi1(2) = dir(workdir,c_aux) c_aux = 'trn/Li2s2p.4' fLi1(3) = dir(workdir,c_aux) c_aux = 'trn/Li2s2p.3' fLi1(4) = dir(workdir,c_aux) c_aux = 'trn/Li2s2p.2' fLi1(5) = dir(workdir,c_aux) c_aux = 'trn/Li2s2p.1' fLi1(6) = dir(workdir,c_aux) ! fLi1(1)='/home/manolo/SRMAC/trn/Li2s2p.6' ! fLi1(2)='/home/manolo/SRMAC/trn/Li2s2p.5' ! fLi1(3)='/home/manolo/SRMAC/trn/Li2s2p.4' ! fLi1(4)='/home/manolo/SRMAC/trn/Li2s2p.3' ! fLi1(5)='/home/manolo/SRMAC/trn/Li2s2p.2' ! fLi1(6)='/home/manolo/SRMAC/trn/Li2s2p.1' iless=0 do j=1,6 ! j will alway over shoot if( dne .lt. dne_array(j) ) iless = 1 if( dne .eq. dne_array(j) ) iless = 2 if( iless .gt.0 ) goto 100 enddo 100 continue iflag=0 if((iless.eq.1).and.(j.eq.1)) then iflag=1 ! Case dne < profile density elseif(iless.eq.2) then iflag=2 ! Case dne = profile density else iflag=3 ! Case prof_min < dne < prof_max jupper=j jlower=j-1 iflag=3 print*,dne_array(jupper) print*,dne_array(jlower) endif if((iflag.eq.1).or.(iflag.eq.2)) then rlognelow = LOG10(dne_array(j)) write(1,*) 'j=',j write(1,*) 'iflag=',iflag open(unit=10,file=fLi1(j)) do i=1,NPDF read(10,*) p%x(i),p%prof1(i) p%prof2(i) = p%prof1(i) enddo close(10) elseif(iflag.eq.3) then rlogdne = LOG10(dne) rlognelow = LOG10(dne_array(jlower)) rlogneup = LOG10(dne_array(jupper)) deltalogne = rlogneup-rlognelow open(unit=10,file=fLi1(jlower)) open(unit=24,file=fLi1(jupper)) do i=1,NPDF read(10,*) p%x(i),prof_low read(24,*) engaux,prof_high deltaprof = prof_high-prof_low slope = deltaprof/deltalogne p%prof1(i) = prof_high+slope*(rlogdne-rlogneup) p%prof2(i) = p%prof1(i) enddo close(10) close(24) else print*,'error: else trap taken' stop endif write(1,*) 'finished Li2s2p' ! ---------------------------------------------- ! open(unit=10,file=fLi1(j)) ! do i=1,NPDF ! read(10,*) aux,rLi2s2prof ! p%x(i) = aux ! p%prof1(i) = rLi2s2prof ! p%prof2(i) = rLi2s2prof ! enddo ! close(10) ! --------------------------------------------- dne_array2(0) = 0.0d0 dne_array2(1) = 2.0e16 ! .8 dne_array2(2) = 3.0e16 ! .13 dne_array2(3) = 5.0e16 ! .7 dne_array2(4) = 7.0d16 ! .12 dne_array2(5) = 1.0d17 ! .6 dne_array2(6) = 1.5d17 ! .9 dne_array2(7) = 2.0d17 ! .5 dne_array2(8) = 3.0d17 ! .10 dne_array2(9) = 4.0d17 ! .11 dne_array2(10) = 5.0d17 ! .4 dne_array2(11) = 1.0d18 ! .3 dne_array2(12) = 1.0e28 c_aux = 'trn/Li3d2p.8' fLi2(1) = dir(workdir,c_aux) c_aux = 'trn/Li3d2p.13' fLi2(2) = dir(workdir,c_aux) c_aux = 'trn/Li3d2p.7' fLi2(3) = dir(workdir,c_aux) c_aux = 'trn/Li3d2p.12' fLi2(4) = dir(workdir,c_aux) c_aux = 'trn/Li3d2p.6 ' fLi2(5) = dir(workdir,c_aux) c_aux = 'trn/Li3d2p.9' fLi2(6) = dir(workdir,c_aux) c_aux = 'trn/Li3d2p.5' fLi2(7) = dir(workdir,c_aux) c_aux = 'trn/Li3d2p.10' fLi2(8) = dir(workdir,c_aux) c_aux = 'trn/Li3d2p.11' fLi2(9) = dir(workdir,c_aux) c_aux = 'trn/Li3d2p.4' fLi2(10)= dir(workdir,c_aux) c_aux = 'trn/Li3d2p.3' fLi2(11)= dir(workdir,c_aux) ! fLi2(1) ='/home/manolo/SRMAC/trn/Li3d2p.8' ! fLi2(2) ='/home/manolo/SRMAC/trn/Li3d2p.13' ! fLi2(3) ='/home/manolo/SRMAC/trn/Li3d2p.7' ! fLi2(4) ='/home/manolo/SRMAC/trn/Li3d2p.12' ! fLi2(5) ='/home/manolo/SRMAC/trn/Li3d2p.6 ' ! fLi2(6) ='/home/manolo/SRMAC/trn/Li3d2p.9' ! fLi2(7) ='/home/manolo/SRMAC/trn/Li3d2p.5' ! fLi2(8) ='/home/manolo/SRMAC/trn/Li3d2p.10' ! fLi2(9) ='/home/manolo/SRMAC/trn/Li3d2p.11' ! fLi2(10)='/home/manolo/SRMAC/trn/Li3d2p.4' ! fLi2(11)='/home/manolo/SRMAC/trn/Li3d2p.3' !----------------------------------------------- iless=0 do j=1,11 ! j will alway over shoot if( dne .lt. dne_array2(j) ) iless = 1 if( dne .eq. dne_array2(j) ) iless = 2 if( iless .gt.0 ) goto 200 enddo 200 continue iflag=0 if((iless.eq.1).and.(j.eq.1)) then iflag=1 ! Case dne < profile density elseif(iless.eq.2) then iflag=2 ! Case dne = profile density else iflag=3 ! Case prof_min < dne < prof_max jupper=j jlower=j-1 iflag=3 print*,dne_array2(jupper) print*,dne_array2(jlower) endif if((iflag.eq.1).or.(iflag.eq.2)) then rlognelow = LOG10(dne_array2(j)) print*,'j=',j print*,'iflag=',iflag open(unit=10,file=fLi2(j)) do i=1,NPDF read(10,*) p%x(i),p%prof3(i) p%prof8(i) = p%prof3(i) p%prof9(i) = p%prof8(i) enddo close(10) elseif(iflag.eq.3) then rlogdne = LOG10(dne) rlognelow = LOG10(dne_array2(jlower)) rlogneup = LOG10(dne_array2(jupper)) deltalogne = rlogneup-rlognelow open(unit=10,file=fLi2(jlower)) open(unit=24,file=fLi2(jupper)) do i=1,NPDF read(10,*) p%x(i),prof_low read(24,*) engaux,prof_high deltaprof = prof_high-prof_low slope = deltaprof/deltalogne p%prof3(i) = prof_high+slope*(rlogdne-rlogneup) p%prof8(i) = p%prof3(i) p%prof9(i) = p%prof8(i) enddo close(10) close(24) else print*,'error: else trap taken' stop endif end ! ************************************************** SUBROUTINE linlocate(eng,x,MF,iloc) ! ************************************************** implicit real*8 (a-h,o-z) real(8), dimension(MF) :: x deltay = MF-1 deltax = x(MF)-x(1) y = (deltay/deltax)*(eng-x(1))+ 1 iloc = int(y,8) END ! ************************************************** SUBROUTINE lininterp(iloc,x,prof,MF,eng,val) ! ************************************************** implicit real*8 (a-h,o-z) real(8), dimension(MF) :: x,prof slope = (prof(iloc+1)-prof(iloc))/(x(iloc+1)-x(iloc)) deltae= eng-x(iloc) val = slope*deltae+prof(iloc) END