! $Id$
! ***********************************************************
! CALCULATE TRANSPORT COEFFICIENTS AND SOURCE
! ***********************************************************
SUBROUTINE TRCALC(IERR)
USE TRCOMM, ONLY : AJBS, AJRF, AR1RHOG, ARRHOG, BP, DT, DVRHOG, MDLEQ0, &
MDLEQB, MDLJBS, MDLNF, MDLUF, MDLPR, MDNCLS, NRAMAX, NRMAX, &
NROMAX, NSM, NSMAX, PBCL, PBIN, PCX, PELTIM, PEX, PFCL, PFIN, &
PI, PIE, PIN, PN, PNB, PNF, POH, PRB, PRC, PRF, PRL, PRSUM, &
Q0, QP, RDP, RG, RHOA, RR, SCX, SEX, SIE, SNB, SNF, SPE, SSIN, &
T, TAUF, TTRHOG, RDPVRHOG
USE tr_cytran_mod
IMPLICIT NONE
INTEGER(4),INTENT(OUT) :: IERR
INTEGER(4) :: NR
REAL(8) :: FCTR
IF(RHOA.NE.1.D0) NRMAX=NROMAX
IERR=0
SIE(1:NRMAX)=0.D0
SNF(1:NRMAX)=0.D0
SNB(1:NRMAX)=0.D0
POH(1:NRMAX)=0.D0
PIE(1:NRMAX)=0.D0
PCX(1:NRMAX)=0.D0
PRB(1:NRMAX)=0.D0
PRC(1:NRMAX)=0.D0
PRL(1:NRMAX)=0.D0
PRSUM(1:NRMAX)=0.D0
PNB(1:NRMAX)=0.D0
PNF(1:NRMAX)=0.D0
PBIN(1:NRMAX)=0.D0
PFIN(1:NRMAX)=0.D0
! AJNB(1:NRMAX)=0.D0
AJRF(1:NRMAX)=0.D0
AJBS(1:NRMAX)=0.D0
SPE(1:NRMAX,1:NSMAX)=0.D0
PBCL(1:NRMAX,1:NSMAX)=0.D0
PFCL(1:NRMAX,1:NSMAX)=0.D0
IF(MDLUF.NE.0) THEN
PRF(1:NRMAX,3:NSM)=0.D0
ELSE
PRF(1:NRMAX,1:NSM)=0.D0
ENDIF
BP(1:NRMAX)=AR1RHOG(1:NRMAX)*RDP(1:NRMAX)/RR
QP(1:NRMAX)=TTRHOG(1:NRMAX)*ARRHOG(1:NRMAX)/(4.D0*PI**2*RDPVRHOG(1:NRMAX))
Q0=FCTR(RG(1),RG(2),QP(1),QP(2))
! *** RADIAL ELECTRIC FIELD ***
CALL TRERAD
IF(T.LT.PELTIM+0.5D0*DT.AND. T.GE.PELTIM-0.5D0*DT) CALL TRPELT
CALL TRZEFF
IF(MDLPR.GT.0) CALL TR_CYTRAN
IF(MDNCLS.NE.0) THEN
CALL TR_NCLASS(IERR)
IF(IERR.NE.0) RETURN
ENDIF
CALL TRCOEF
CALL TRLOSS
IF(MDLUF.NE.1.AND.MDLUF.NE.3) CALL TRPWRF
CALL TRPWNB
IF(MDNCLS.NE.0) THEN
CALL TRAJBS_NCLASS
ELSE
select case(MDLJBS)
case(1)
CALL TRAJBS
case(2)
CALL TRAJBS
case(3)
CALL TRAJBS
case(4)
CALL TRAJBSNEW
case(5)
CALL TRAJBSSAUTER
case default
CALL TRAJBS
end select
ENDIF
SELECT CASE(MDLNF)
CASE(0)
TAUF(1:NRMAX)=1.D0
CASE(1:4)
CALL TRNFDT
CASE(5:6)
CALL TRNFDHE3
END SELECT
CALL TRAJOH
DO NR=1,NRMAX
IF(MDLEQ0.EQ.0) THEN
SSIN(NR,1)=SIE(NR) +SNB(NR)+SEX(NR,1)
SSIN(NR,2)=PN(2)*SIE(NR)/(PN(2)+PN(3))-SNF(NR)+SNB(NR)+SEX(NR,2)
SSIN(NR,3)=PN(3)*SIE(NR)/(PN(2)+PN(3))-SNF(NR) +SEX(NR,3)
SSIN(NR,4)= SNF(NR) +SEX(NR,4)
SSIN(NR,7)=-SIE(NR) -SCX(NR)
SSIN(NR,8)= SNB(NR)+SCX(NR)
ELSEIF(MDLEQ0.EQ.1) THEN
SSIN(NR,1)= SNB(NR)+SEX(NR,1)
SSIN(NR,2)= -SNF(NR)+SNB(NR)+SEX(NR,2)
SSIN(NR,3)= -SNF(NR) +SEX(NR,3)
SSIN(NR,4)= SNF(NR) +SEX(NR,4)
SSIN(NR,7)=0.D0
SSIN(NR,8)= SNB(NR)
ENDIF
PIN(NR,1)=PBCL(NR,1)+PFCL(NR,1)+PRF(NR,1) &
& +POH(NR)-PRSUM(NR)-PIE(NR)+PEX(NR,1)
PIN(NR,2)=PBCL(NR,2)+PFCL(NR,2)+PRF(NR,2) &
& -PN(2)*PCX(NR)/(PN(2)+PN(3))+PEX(NR,2)
PIN(NR,3)=PBCL(NR,3)+PFCL(NR,3)+PRF(NR,3) &
& -PN(3)*PCX(NR)/(PN(2)+PN(3))+PEX(NR,3)
PIN(NR,4)=PBCL(NR,4)+PFCL(NR,4)+PRF(NR,4)+PEX(NR,4)
ENDDO
IF(RHOA.NE.1.D0) NRMAX=NRAMAX
RETURN
END SUBROUTINE TRCALC
! ***********************************************************
! RADIAL ELECTRIC FIELD
! **********************************************************
SUBROUTINE TRERAD
USE TRCOMM, ONLY : AEE, AMM, BB, BP, DR, EPSRHO, ER, MDLER, NRMAX, &
& PA, PADD, PBM, PNSS, PTS, PZ, QP, RHOG, RHOM, &
& RJCB, RKEV, RN, RNF, RT, SUMPBM, VPOL, VTOR
IMPLICIT NONE
INTEGER(4):: NR
REAL(8) :: ALPHA_NEO, CS, DPD, DRL, EPS, F_UNTRAP, RHO_S, RLNI, &
& RLTI, TEL, TERM_DP, TIL
REAL(8) :: DERIV3P
IF(SUMPBM.EQ.0.D0) THEN
PADD(1:NRMAX)=0.D0
ELSE
PADD(1:NRMAX)=PBM(1:NRMAX)*1.D-20/RKEV-RNF(1:NRMAX,1)*RT(1:NRMAX,2)
ENDIF
DO NR=1,NRMAX
DRL=RJCB(NR)/DR
IF(NR.EQ.NRMAX) THEN
DPD = DERIV3P(PNSS(2)*PTS(2),RN(NR ,2)*RT(NR ,2)-PADD(NR ), &
& RN(NR-1,2)*RT(NR-1,2)-PADD(NR-1),RHOG(NR),RHOM(NR),RHOM(NR-1))
TERM_DP = DPD*RKEV/(PZ(2)*AEE*PNSS(2))
ELSE
DPD =( RN(NR+1,2)*RT(NR+1,2)-PADD(NR+1)-(RN(NR ,2)*RT(NR ,2)-PADD(NR )))*DRL
TERM_DP = DPD*RKEV/(PZ(2)*AEE*0.5D0*(RN(NR+1,2)+RN(NR,2)))
ENDIF
IF(MDLER.EQ.0) THEN
! pressure gradient only
ER(NR) = TERM_DP
ELSEIF(MDLER.EQ.1) THEN
! nabla p + toroidal rotation
ER(NR) = TERM_DP+VTOR(NR)*BP(NR)
ELSEIF(MDLER.EQ.2) THEN
! nabla p + V_tor + poloidal rotation
ER(NR) = TERM_DP+VTOR(NR)*BP(NR)-VPOL(NR)*BB
ELSEIF(MDLER.EQ.3) THEN
! Waltz definition
EPS = EPSRHO(NR)
F_UNTRAP = 1.D0-1.46D0*SQRT(EPS)+0.46D0*EPS**1.5D0
ALPHA_NEO = 1.D0-0.8839D0*F_UNTRAP/(0.3477D0+0.4058D0*F_UNTRAP)
IF(NR.EQ.NRMAX) THEN
TEL = PTS(1)
TIL = PTS(2)
RLNI = -DERIV3P(PNSS(2),RN(NR,2),RN(NR-1,2),RHOG(NR),RHOM(NR),RHOM(NR-1))/PNSS(2)
RLTI = -DERIV3P(PTS(2),RT(NR,2),RT(NR-1,2),RHOG(NR),RHOM(NR),RHOM(NR-1))/PTS(2)
ELSE
TEL = 0.5D0*(RT(NR,1)+RT(NR+1,1))
TIL = 0.5D0*(RT(NR,2)+RT(NR+1,2))
RLNI = -(LOG(RN(NR+1,2))-LOG(RN(NR,2)))*DRL
RLTI = -(LOG(RT(NR+1,2))-LOG(RT(NR,2)))*DRL
ENDIF
CS = SQRT(TEL*RKEV/(PA(2)*AMM))
RHO_S = CS*PA(2)*AMM/(PZ(2)*AEE*BB)
ER(NR) =-BB*( (TIL/TEL)*RHO_S*CS*(RLNI+ALPHA_NEO*RLTI)-EPS/QP(NR)*VTOR(NR))
ENDIF
ENDDO
RETURN
END SUBROUTINE TRERAD
! ***********************************************************
! CALCULATE Z-C
! **********************************************************
REAL(8) FUNCTION TRZEC(TE)
IMPLICIT NONE
REAL(8):: TE,TEL
REAL(8):: BC00=2.093736D+03, BC01=5.153766D+03,BC02=5.042105D+03, BC03=2.445345D+03, BC04=5.879207D+02, BC05=5.609128D+01
REAL(8):: BC10=7.286051D+00, BC11=4.506650D+01,BC12=1.595483D+02, BC13=2.112702D+02, BC14=1.186473D+02, BC15=2.400040D+01
REAL(8):: BC20=5.998782D+00, BC21=6.808206D-03,BC22=-3.096242D-02,BC23=-4.794194D-02,BC24=1.374732D-01, BC25=5.157571D-01
REAL(8):: BC30=5.988153D+00, BC31=7.335329D-02,BC32=-1.754858D-01,BC33=2.034126D-01, BC34=-1.145930D-01,BC35=2.517150D-02
REAL(8):: BC40=-3.412166D+01,BC41=1.250454D+02,BC42=-1.550822D+02,BC43=9.568297D+01, BC44=-2.937297D+01,BC45=3.589667D+00
TEL = LOG10(TE)
IF(TE.LE.3.D-3) THEN
TRZEC=0.D0
ELSEIF(TE.LE.2.D-2) THEN
TRZEC= BC00+(BC01*TEL)+(BC02*TEL**2)+(BC03*TEL**3)+(BC04*TEL**4)+(BC05*TEL**5)
ELSEIF(TE.LE.0.2D0) THEN
TRZEC= BC10+(BC11*TEL)+(BC12*TEL**2)+(BC13*TEL**3)+(BC14*TEL**4)+(BC15*TEL**5)
ELSEIF(TE.LE.2.D0) THEN
TRZEC= BC20+(BC21*TEL)+(BC22*TEL**2)+(BC23*TEL**3)+(BC24*TEL**4)+(BC25*TEL**5)
ELSEIF(TE.LE.20.D0) THEN
TRZEC= BC30+(BC31*TEL)+(BC32*TEL**2)+(BC33*TEL**3)+(BC34*TEL**4)+(BC35*TEL**5)
ELSEIF(TE.LE.100.D0) THEN
TRZEC= BC40+(BC41*TEL)+(BC42*TEL**2)+(BC43*TEL**3)+(BC44*TEL**4)+(BC45*TEL**5)
ELSE
TRZEC= 6.D0
ENDIF
RETURN
END FUNCTION TRZEC
! ***********************************************************
! CALCULATE Z-FE
! ***********************************************************
REAL(8) FUNCTION TRZEFE(TE)
IMPLICIT NONE
REAL(8):: TE, TEL
REAL(8):: BF10=8.778318D+00, BF11=-5.581412D+01,BF12=-1.225124D+02,BF13=-1.013985D+02,BF14=-3.914244D+01, BF15=-5.851445D+00
REAL(8):: BF20=1.959496D+01, BF21=1.997852D+01, BF22=9.593373D+00, BF23=-7.609962D+01,BF24=-1.007190D+02, BF25=-1.144046D+01
REAL(8):: BF30=1.478150D+01, BF31=6.945311D+01, BF32=-1.991202D+02,BF33=2.653804D+02, BF34=-1.631423D+02, BF35=3.770026D+01
REAL(8):: BF40=2.122081D+01, BF41=6.891607D+00, BF42=-4.076853D+00,BF43=1.577171D+00, BF44=-5.139468D-01, BF45=8.934176D-02
TEL = LOG10(TE)
IF(TE.LE.3.D-3) THEN
TRZEFE=0.D0
ELSEIF(TE.LE.0.2D0) THEN
TRZEFE = BF10+(BF11*TEL)+(BF12*TEL**2)+(BF13*TEL**3)+(BF14*TEL**4)+(BF15*TEL**5)
ELSEIF(TE.LE.2.D0) THEN
TRZEFE = BF20+(BF21*TEL)+(BF22*TEL**2)+(BF23*TEL**3)+(BF24*TEL**4)+(BF25*TEL**5)
ELSEIF(TE.LE.20.D0) THEN
TRZEFE = BF30+(BF31*TEL)+(BF32*TEL**2)+(BF33*TEL**3)+(BF34*TEL**4)+(BF35*TEL**5)
ELSEIF(TE.LE.100.D0) THEN
TRZEFE = BF40+(BF41*TEL)+(BF42*TEL**2)+(BF43*TEL**3)+(BF44*TEL**4)+(BF45*TEL**5)
ELSE
TRZEFE = 26.D0
ENDIF
RETURN
END FUNCTION TRZEFE
! ***********************************************************
! CALCULATE Z-EFFF
! ***********************************************************
SUBROUTINE TRZEFF
USE TRCOMM, ONLY : ANC, ANFE, MDLEQN, MDLUF, NRMAX, PZ, PZC, PZFE, RN, RNF, RT, ZEFF
IMPLICIT NONE
INTEGER(4):: NR
REAL(8) :: TE, TRZEC, TRZEFE
IF(MDLUF.EQ.0) THEN
DO NR=1,NRMAX
TE =RT(NR,1)
ZEFF(NR) =(PZ(2) *PZ(2) *RN(NR,2) +PZ(3) *PZ(3) *RN(NR,3) +PZ(4) *PZ(4) *RN(NR,4) &
& +TRZEC(TE)**2 *ANC (NR) +TRZEFE(TE)**2 *ANFE(NR))/RN(NR,1)
ENDDO
ELSE
IF(MDLEQN.EQ.0) THEN ! fixed density
DO NR=1,NRMAX
ZEFF(NR) =(PZ(2) *PZ(2) *RN(NR,2) +PZ(3) *PZ(3) *RN(NR,3) +PZ(2) *PZ(2) *RNF(NR,1))/RN(NR,1)
ENDDO
ENDIF
ENDIF
DO NR=1,NRMAX
TE=RT(NR,1)
PZC(NR)=TRZEC(TE)
PZFE(NR)=TRZEFE(TE)
ENDDO
!
RETURN
END SUBROUTINE TRZEFF
! ***********************************************************
! RADIATION POWER - C
! ***********************************************************
REAL(8) FUNCTION TRRPC(TE)
IMPLICIT NONE
REAL(8):: TE,TEL,ARG
REAL(8):: AC00=1.965300D+03, AC01=4.572039D+03, AC02=4.159590D+03, AC03=1.871560D+03, AC04=4.173889D+02, AC05=3.699382D+01
REAL(8):: AC10=7.467599D+01, AC11=4.549038D+02, AC12=8.372937D+02, AC13=7.402515D+02, AC14=3.147607D+02, AC15=5.164578D+01
REAL(8):: AC20=-2.120151D+01,AC21=-3.668933D-01,AC22=7.295099D-01, AC23=-1.944827D-01,AC24=-1.263576D-01,AC25=-1.491027D-01
REAL(8):: AC30=-2.121979D+01,AC31=-2.346986D-01,AC32=4.093794D-01, AC33=7.874548D-02, AC34=-1.841379D-01,AC35=5.590744D-02
REAL(8):: AC40=-2.476796D+01,AC41=9.408181D+00, AC42=-9.657446D+00,AC43=4.999161D+00, AC44=-1.237382D+00,AC45=1.160610D-01
TEL = LOG10(TE)
IF(TE.LE.3.D-3) THEN
TEL=LOG10(3.D-3)
ARG=AC00+(AC01*TEL)+(AC02*TEL**2)+(AC03*TEL**3)+(AC04*TEL**4)+(AC05*TEL**5)
ELSEIF(TE.LE.2.D-2) THEN
ARG=AC00+(AC01*TEL)+(AC02*TEL**2)+(AC03*TEL**3)+(AC04*TEL**4)+(AC05*TEL**5)
ELSEIF(TE.LE.0.2D0) THEN
ARG=AC10+(AC11*TEL)+(AC12*TEL**2)+(AC13*TEL**3)+(AC14*TEL**4)+(AC15*TEL**5)
ELSEIF(TE.LE.2.D0) THEN
ARG=AC20+(AC21*TEL)+(AC22*TEL**2)+(AC23*TEL**3)+(AC24*TEL**4)+(AC25*TEL**5)
ELSEIF(TE.LE.20.D0) THEN
ARG=AC30+(AC31*TEL)+(AC32*TEL**2)+(AC33*TEL**3)+(AC34*TEL**4)+(AC35*TEL**5)
ELSEIF(TE.LE.1.D2) THEN
ARG=AC40+(AC41*TEL)+(AC42*TEL**2)+(AC43*TEL**3)+(AC44*TEL**4)+(AC45*TEL**5)
ENDIF
TRRPC = 10.D0**(ARG-13.D0)
RETURN
END FUNCTION TRRPC
! ***********************************************************
! RADIATION POWER - FE
! ***********************************************************
REAL(8) FUNCTION TRRPFE(TE)
IMPLICIT NONE
REAL(8)::TE,TEL,ARG
REAL(8)::AF10=-2.752599D+01,AF11=-3.908228D+01,AF12=-6.469423D+01,AF13=-5.555048D+01,AF14=-2.405568D+01,AF15=-4.093160D+00
REAL(8)::AF20=-1.834973D+01,AF21=-1.252028D+00,AF22=-7.533115D+00,AF23=-3.289693D+00,AF24=2.866739D+01, AF25=2.830249D+01
REAL(8)::AF30=-1.671042D+01,AF31=-1.646143D+01,AF32=3.766238D+01, AF33=-3.944080D+01,AF34=1.918529D+01, AF35=-3.509238D+00
REAL(8)::AF40=-2.453957D+01,AF41=1.795222D+01, AF42=-2.356360D+01,AF43=1.484503D+01, AF44=-4.542323D+00,AF45=5.477462D-01
TEL = LOG10(TE)
IF(TE.LE.3.D-3) THEN
TEL=LOG10(3.D-3)
ARG=AF10+(AF11*TEL)+(AF12*TEL**2)+(AF13*TEL**3)+(AF14*TEL**4)+(AF15*TEL**5)
ELSEIF(TE.LE.0.2D0) THEN
ARG=AF10+(AF11*TEL)+(AF12*TEL**2)+(AF13*TEL**3)+(AF14*TEL**4)+(AF15*TEL**5)
ELSEIF(TE.LE.2.D0) THEN
ARG=AF20+(AF21*TEL)+(AF22*TEL**2)+(AF23*TEL**3)+(AF24*TEL**4)+(AF25*TEL**5)
ELSEIF(TE.LE.20.D0) THEN
ARG=AF30+(AF31*TEL)+(AF32*TEL**2)+(AF33*TEL**3)+(AF34*TEL**4)+(AF35*TEL**5)
ELSEIF(TE.LE.1.D2) THEN
ARG=AF40+(AF41*TEL)+(AF42*TEL**2)+(AF43*TEL**3)+(AF44*TEL**4)+(AF45*TEL**5)
ENDIF
TRRPFE = 10.D0**(ARG-13.D0)
RETURN
END FUNCTION TRRPFE
! ***********************************************************
! POWER LOSS
! ***********************************************************
SUBROUTINE TRLOSS
USE TRCOMM, ONLY : AEE, ANC, ANFE, ANNU, DT, KUFDEV, MDLUF, MDLPR, &
NRMAX, NT, NTUM, PCX, PIE, PRB, PRC, PRSUM, PRL, PRLU, RKEV, RN, &
RT, SCX, SIE, TSCX, TSIE
USE TRCOM1, ONLY : NTAMAX, NTXMAX, PNBI, TMU
USE tr_cytran_mod, ONLY: tr_cytran
IMPLICIT NONE
INTEGER(4):: IERR, NR
REAL(8):: ANDX, ANE, ANHE, ANT, EION, PLC, PLD, PLFE, PLHE, PLTT, PRLL, SCH, SION, TD, TE, TN, TNU, TRRPC, TRRPFE, TSL
IF(MDLUF.EQ.1.OR.MDLUF.EQ.3) THEN
IF(NT.EQ.0) THEN
TSL=DT*DBLE(1)
DO NR=1,NRMAX
PRSUM(NR)=PRLU(1,NR)
ENDDO
ELSE
TSL=DT*DBLE(NT)
IF(KUFDEV.EQ.'X') THEN
DO NR=1,NRMAX
CALL TIMESPL(TSL,PRLL,TMU,PRLU(1,NR),NTXMAX,NTUM,IERR)
PRSUM(NR)=PRLL
ENDDO
ELSE
DO NR=1,NRMAX
IF(PNBI.LT.12.D6) THEN
PRSUM(NR)=PRLU(1,NR)
ELSE
PRSUM(NR)=PRLU(2,NR)
IF(NT.EQ.NTAMAX) PRSUM(NR)=PRLU(3,NR)
ENDIF
ENDDO
ENDIF
ENDIF
ELSEIF(MDLUF.EQ.2) THEN
IF(NT.EQ.0) THEN
TSL=DT*DBLE(1)
DO NR=1,NRMAX
PRSUM(NR)=PRLU(1,NR)
ENDDO
ELSE
TSL=DT*DBLE(NT)
DO NR=1,NRMAX
PRSUM(NR)=PRLU(1,NR)
ENDDO
ENDIF
ELSE
DO NR=1,NRMAX
ANE =RN(NR,1)
ANDX=RN(NR,2)
ANT =RN(NR,3)
ANHE=RN(NR,4)
TE =RT(NR,1)
! Radiation loss caused by impurities
PLFE = ANE*ANFE(NR)*TRRPFE(TE)*1.D40
PLC = ANE*ANC (NR)*TRRPC (TE)*1.D40
PRL(NR)=PLFE+PLC
! Bremsstrahlung
PLD = ANE*ANDX*5.35D-37*1.D0**2*SQRT(ABS(TE))*1.D40
PLTT = ANE*ANT *5.35D-37*1.D0**2*SQRT(ABS(TE))*1.D40
PLHE = ANE*ANHE*5.35D-37*2.D0**2*SQRT(ABS(TE))*1.D40
PRB(NR)= PLD+PLTT+PLHE
! Sumup
SELECT CASE(MDLPR)
CASE(0)
PRSUM(NR)=PRL(NR)+PRB(NR)
CASE(1,2)
PRSUM(NR)=PRL(NR)+PRB(NR)+PRC(NR)
END SELECT
ENDDO
ENDIF
! ****** IONIZATION LOSS ******
DO NR=1,NRMAX
ANE=RN(NR,1)
TE =RT(NR,1)
EION = 13.64D0
! SIONS = 1.76D-13/SQRT(TE*1.D3)*FKN(EION/(TE*1.D3))
TN = MAX(TE*1.D3/EION,1.D-2)
SION = 1.D-11*SQRT(TN)*EXP(-1.D0/TN)/(EION**1.5D0*(6.D0+TN))
PIE(NR) = ANE*ANNU(NR)*SION*1.D40*EION*AEE
SIE(NR) = ANE*ANNU(NR)*SION*1.D20
TSIE(NR)= ANE *SION*1.D20
ENDDO
! ****** CHARGE EXCHANGE LOSS ******
DO NR=1,NRMAX
ANE =RN(NR,1)
ANDX=RN(NR,2)
TE =RT(NR,1)
TD =RT(NR,2)
TNU = 0.D0
TN = LOG10(MAX(TD*1.D3,50.D0))
SCH= 1.57D-16*SQRT(ABS(TD)*1.D3)*(TN*TN-14.63D0*TN+53.65D0)
SCX(NR) =ANDX*ANNU(NR)*SCH*1.D20
TSCX(NR)=ANDX *SCH*1.D20
PCX(NR)=1.5D0*ANDX*ANNU(NR)*SCH*(TD-TNU)*RKEV*1.D40
ENDDO
RETURN
END SUBROUTINE TRLOSS
REAL(8) FUNCTION FKN(X)
IMPLICIT NONE
REAL(8) X,GAMMA
GAMMA=0.577215664901532D0
FKN=-(LOG10(X)+GAMMA-X+X**2/4.D0-X**3/1.8D1)
RETURN
END FUNCTION FKN
! ***********************************************
! BOOTSTRAP CURRENT (NCLASS)
! ***********************************************
SUBROUTINE TRAJBS_NCLASS
USE TRCOMM, ONLY : AJBS, AJBSNC, BB, CJBSP, CJBST, DR, NRMAX, NSMAX, PBSCD, PNSS, PTS, RG, RM, RN, RT
IMPLICIT NONE
INTEGER(4):: NR, NS, NSW
REAL(8) :: DRPNW, DRTNW, RPNW, RTNW, SUML
REAL(8),DIMENSION(NRMAX):: AJBSL
REAL(8) :: DERIV3P
IF(PBSCD.LE.0.D0) RETURN
NSW=1
IF(NSW.EQ.0) THEN
AJBSL(1:NRMAX)=PBSCD*AJBSNC(1:NRMAX)
ELSE
DO NR=1,NRMAX-1
SUML=0.D0
DO NS=1,NSMAX
RTNW=0.5D0*(RT(NR+1,NS)+RT(NR,NS))
RPNW=0.5D0*(RN(NR+1,NS)*RT(NR+1,NS)+RN(NR ,NS)*RT(NR ,NS))
DRTNW=(RT(NR+1,NS)-RT(NR,NS))/DR
DRPNW=(RN(NR+1,NS)*RT(NR+1,NS)-RN(NR,NS)*RT(NR,NS))/DR
SUML=SUML+CJBST(NR,NS)*DRTNW/RTNW +CJBSP(NR,NS)*DRPNW/RPNW
ENDDO
AJBSL(NR)=-PBSCD*SUML/BB
ENDDO
NR=NRMAX
SUML=0.D0
DO NS=1,NSMAX
RTNW=PTS(NS)
RPNW=PNSS(NS)*PTS(NS)
DRTNW=DERIV3P(PTS(NS),RT(NR,NS),RT(NR-1,NS),RG(NR),RM(NR),RM(NR-1))
DRPNW=DERIV3P(PNSS(NS)*PTS(NS),RN(NR ,NS)*RT(NR ,NS),RN(NR-1,NS)*RT(NR-1,NS),RG(NR),RM(NR),RM(NR-1))
SUML=SUML+CJBST(NR,NS)*DRTNW/RTNW +CJBSP(NR,NS)*DRPNW/RPNW
!CC if(ns.eq.3) write(6,*) NR,PNSS(NS),PTS(NS)
ENDDO
AJBSL(NR)=-PBSCD*SUML/BB
ENDIF
!
AJBS(1)=0.5D0*AJBSL(1)
DO NR=2,NRMAX
AJBS(NR)=0.5D0*(AJBSL(NR)+AJBSL(NR-1))
ENDDO
RETURN
END SUBROUTINE TRAJBS_NCLASS
! ***********************************************
! BOOTSTRAP CURRENT (O. SAUTER)
! ***********************************************
SUBROUTINE TRAJBSSAUTER
USE TRCOMM, ONLY : AJBS, BB, DR, EPSRHO, MDLTPF, NRMAX, NSM, PADD, PBSCD, PNSS, PTS, PZ, QP, RDP, RHOG, RHOM, RKEV, &
& RN, RPE, RR, RT, RW, TTRHOG, ZEFF
IMPLICIT NONE
INTEGER(4):: NR, NS
REAL(8) :: ANE, DPE, DPI, DRL, DTE, DTI, EPS, EPSS, F31TEFF, F32EETEFF, F32EITEFF, F34TEFF, FT, FTPF, PE, PPI, QL, RL31,&
& RL32, RL34, RLNLAME, RLNLAMII, RNM, RNP, RNTM, RNTP, RNUE, RNUI, RPIM, RPIP, SALFA, SALFA0, TE, TI, ZEFFL
REAL(8),DIMENSION(NRMAX):: AJBSL, ANI
REAL(8) :: DERIV3P, F31, F32EE, F32EI
IF(PBSCD.LE.0.D0) RETURN
DO NR=1,NRMAX-1
EPS=EPSRHO(NR)
EPSS=SQRT(EPS)**3
QL=ABS(QP(NR))
ZEFFL=0.5D0*(ZEFF(NR)+ZEFF(NR+1))
DRL=1.D0/DR
RNTP=0.D0
RNP =0.D0
RNTM=0.D0
RNM =0.D0
DO NS=2,NSM
RNTP=RNTP+RN(NR+1,NS)*RT(NR+1,NS)
RNP =RNP +RN(NR+1,NS)
RNTM=RNTM+RN(NR ,NS)*RT(NR ,NS)
RNM =RNM +RN(NR ,NS)
ENDDO
RNTP=RNTP+RW(NR+1,1)+RW(NR+1,2)
RNTM=RNTM+RW(NR ,1)+RW(NR ,2)
RPIP=RNTP+PADD(NR+1)
RPIM=RNTM+PADD(NR )
! ****** ION PARAMETER ******
! *** ANI is the the ion density (ni) ***
! *** TI is the ion temperature (Ti) ***
! *** DTI is the derivative of ion temperature (dTi/dr) ***
! *** PPI is the ion pressure (Pi) ***
! *** DPI is the derivative of ion pressure (dPi/dr) ***
ANI(NR)=0.D0
DO NS=2,NSM
ANI(NR)=ANI(NR)+0.5D0*PZ(NS)*(RN(NR+1,NS)+RN(NR,NS))
ENDDO
PPI=0.5D0*(RPIP+RPIM)
DPI=(RPIP-RPIM)*DRL
TI =0.5D0*(RNTP/RNP+RNTM/RNM)
DTI=(RNTP/RNP-RNTM/RNM)*DRL
rLnLamii=30.D0-LOG(PZ(2)**3*SQRT(ANI(NR)*1.D20)/(ABS(TI*1.D3)**1.5D0))
RNUI=4.90D-18*QL*RR*ANI(NR)*1.D20*PZ(2)**4*rLnLamii /(ABS(TI*1.D3)**2*EPSS)
!
! ****** ELECTORON PARAMETER ******
! *** ANE is the the electron density (ne) ***
! *** TE is the electron temperature (Te) ***
! *** PE is the electron pressure (Pe) ***
! *** DTE is the derivative of electron temperature (dTe/dr) ***
! *** DPE is the derivative of electron pressure (dPe/dr) ***
ANE=0.5D0*(RN(NR+1,1)+RN(NR,1))
TE= 0.5D0*(RT(NR+1,1)+RT(NR,1))
PE= 0.5D0*(RN(NR+1,1)*RT(NR+1,1)+RN(NR,1)*RT(NR,1))
DTE=(RT(NR+1,1)-RT(NR,1))*DRL
DPE=(RN(NR+1,1)*RT(NR+1,1)-RN(NR,1)*RT(NR,1))*DRL
rLnLame=31.3D0-LOG(SQRT(ANE*1.D20)/ABS(TE*1.D3))
RNUE=6.921D-18*QL*RR*ANE*1.D20*ZEFFL*rLnLame /(ABS(TE*1.D3)**2*EPSS)
RPE=PE/(PE+PPI)
FT=FTPF(MDLTPF,EPS)
! F33TEFF=FT/(1.D0+(0.55D0-0.1D0*FT)*SQRT(RNUE)+0.45D0*(1.D0-FT)*RNUE/ZEFFL**1.5)
F31TEFF=FT/(1.D0+(1.D0-0.1D0*FT)*SQRT(RNUE) +0.5D0*(1.D0-FT)*RNUE/ZEFFL)
F32EETEFF=FT/(1.D0+0.26D0*(1.D0-FT)*SQRT(RNUE)+0.18D0*(1.D0-0.37D0*FT)*RNUE/SQRT(ZEFFL))
F32EITEFF=FT/(1.D0+(1.D0+0.6D0*FT)*SQRT(RNUE) +0.85D0*(1.D0-0.37D0*FT)*RNUE*(1.D0+ZEFFL))
F34TEFF=FT/(1.D0+(1.D0-0.1D0*FT)*SQRT(RNUE) +0.5D0*(1.D0-0.5D0*FT)*RNUE/ZEFFL)
SALFA0=-1.17D0*(1.D0-FT)/(1.D0-0.22D0*FT-0.19D0*FT**2)
SALFA=((SALFA0+0.25D0*(1.D0-FT**2)*SQRT(RNUI))/(1.D0+0.5D0*SQRT(RNUI))+0.315D0*RNUI**2*FT**6) &
& /(1.D0+0.15D0*RNUI**2*FT**6)
! RNZ=0.58D0+0.74D0/(0.76D0+ZEFFL)
! SGMSPTZ=1.9012D4*(TE*1.D3)**1.5/(ZEFFL*RNZ*rLnLame)
! SGMNEO=SGMSPTZ*F33(F33TEFF,ZEFFL)
RL31=F31(F31TEFF,ZEFFL)
RL32=F32EE(F32EETEFF,ZEFFL)+F32EI(F32EITEFF,ZEFFL)
RL34=F31(F34TEFF,ZEFFL)
AJBSL(NR)=-PBSCD*TTRHOG(NR)*PE*1.D20*RKEV *( RL31*(DPE/PE+DPI/PE)+RL32*DTE/TE &
& +RL34*SALFA*(1.D0-RPE)/RPE*DTI/TI)/RDP(NR)/BB
ENDDO
NR=NRMAX
EPS=EPSRHO(NR)
EPSS=SQRT(EPS)**3
QL=ABS(QP(NR))
ZEFFL=2.D0*ZEFF(NR-1)-ZEFF(NR-2)
DRL=1.D0/DR
! In the following, we assume that
! 1. pressures of beam and fusion at rho=1 are negligible,
! 2. calibration of Pbeam (from UFILE) is not necessary at rho=1.
RNTP=0.D0
RNP =0.D0
RNTM=0.D0
RNM =0.D0
DO NS=2,NSM
RNTP=RNTP+PNSS(NS)*PTS(NS)
RNP =RNP +PNSS(NS)
RNTM=RNTM+RN(NR-1,NS)*RT(NR-1,NS)+RN(NR ,NS)*RT(NR ,NS)
RNM =RNM +RN(NR-1,NS)+RN(NR ,NS)
ENDDO
RNTM=RNTM+RW(NR-1,1)+RW(NR-1,2)+RW(NR ,1)+RW(NR ,2)
RPIP=RNTP
RPIM=RNTM+PADD(NR-1)+PADD(NR )
RNTM=0.5D0*RNTM
RNM =0.5D0*RNM
RPIM=0.5D0*RPIM
! ****** ION PARAMETER ******
! *** ANI is the the ion density (ni) ***
! *** TI is the ion temperature (Ti) ***
! *** DTI is the derivative of ion temperature (dTi/dr) ***
! *** PPI is the ion pressure (Pi) ***
! *** DPI is the derivative of ion pressure (dPi/dr) ***
ANI(NR)=0.D0
DO NS=2,NSM
ANI(NR)=ANI(NR)+PZ(NS)*PNSS(NS)
ENDDO
PPI=RPIP
DPI=(RPIP-RPIM)*DRL
TI =RNTP/RNP
DTI=(RNTP/RNP-RNTM/RNM)*DRL
rLnLamii=30.D0-LOG(PZ(2)**3*SQRT(ANI(NR)*1.D20)/(ABS(TI*1.D3)**1.5D0))
RNUI=4.90D-18*QL*RR*ANI(NR)*1.D20*PZ(2)**4*rLnLamii /(ABS(TI*1.D3)**2*EPSS)
!
! ****** ELECTORON PARAMETER ******
! *** ANE is the the electron density (ne) ***
! *** TE is the electron temperature (Te) ***
! *** PE is the electron pressure (Pe) ***
! *** DTE is the derivative of electron temperature (dTe/dr) ***
! *** DPE is the derivative of electron pressure (dPe/dr) ***
ANE=PNSS(1)
TE =PTS(1)
PE =PNSS(1)*PTS(1)
DTE=DERIV3P(PTS(1),RT(NR,1),RT(NR-1,1),RHOG(NR),RHOM(NR),RHOM(NR-1))
DPE=DERIV3P(PNSS(1)*PTS(1),RN(NR,1)*RT(NR,1),RN(NR-1,1)*RT(NR-1,1),RHOG(NR),RHOM(NR),RHOM(NR-1))
!
rLnLame=31.3D0-LOG(SQRT(ANE*1.D20)/ABS(TE*1.D3))
RNUE=6.921D-18*QL*RR*ANE*1.D20*ZEFFL*rLnLame /(ABS(TE*1.D3)**2*EPSS)
!
RPE=PE/(PE+PPI)
FT=FTPF(MDLTPF,EPS)
! F33TEFF=FT/(1.D0+(0.55D0-0.1D0*FT)*SQRT(RNUE)+0.45D0*(1.D0-FT)*RNUE/ZEFFL**1.5)
F31TEFF=FT/(1.D0+(1.D0-0.1D0*FT)*SQRT(RNUE) +0.5D0*(1.D0-FT)*RNUE/ZEFFL)
F32EETEFF=FT/(1.D0+0.26D0*(1.D0-FT)*SQRT(RNUE)+0.18D0*(1.D0-0.37D0*FT)*RNUE/SQRT(ZEFFL))
F32EITEFF=FT/(1.D0+(1.D0+0.6D0*FT)*SQRT(RNUE) +0.85D0*(1.D0-0.37D0*FT)*RNUE*(1.D0+ZEFFL))
F34TEFF=FT/(1.D0+(1.D0-0.1D0*FT)*SQRT(RNUE) +0.5D0*(1.D0-0.5D0*FT)*RNUE/ZEFFL)
SALFA0=-1.17D0*(1.D0-FT)/(1.D0-0.22D0*FT-0.19D0*FT**2)
SALFA=((SALFA0+0.25D0*(1.D0-FT**2)*SQRT(RNUI)) /(1.D0+0.5D0*SQRT(RNUI))+0.315D0*RNUI**2*FT**6) &
& /(1.D0+0.15D0*RNUI**2*FT**6)
! RNZ=0.58D0+0.74D0/(0.76D0+ZEFFL)
! SGMSPTZ=1.9012D4*(TE*1.D3)**1.5/(ZEFFL*RNZ*rLnLame)
! SGMNEO=SGMSPTZ*F33(F33TEFF,ZEFFL)
!
RL31=F31(F31TEFF,ZEFFL)
RL32=F32EE(F32EETEFF,ZEFFL)+F32EI(F32EITEFF,ZEFFL)
RL34=F31(F34TEFF,ZEFFL)
AJBSL(NR)=-PBSCD*TTRHOG(NR)*PE*1.D20*RKEV*( RL31*(DPE/PE+DPI/PE)+RL32*DTE/TE &
& +RL34*SALFA*(1.D0-RPE)/RPE*DTI/TI)/RDP(NR)/BB
AJBS(1)=0.5D0*AJBSL(1)
DO NR=2,NRMAX
AJBS(NR)=0.5D0*(AJBSL(NR)+AJBSL(NR-1))
ENDDO
RETURN
END SUBROUTINE TRAJBSSAUTER
! *********************
! * Fitting Function *
! *********************
REAL(8) FUNCTION F33(X,Z)
IMPLICIT NONE
REAL(8) X,Z
F33 = 1.D0-(1.D0+0.36D0/Z)*X+0.59D0/Z*X**2-0.23D0/Z*X**3
RETURN
END FUNCTION F33
REAL(8) FUNCTION F31(X,Z)
IMPLICIT NONE
REAL(8) X,Z
F31 = (1.D0+1.4D0/(Z+1.D0))*X-1.9D0/(Z+1.D0)*X**2 +0.3D0/(Z+1.D0)*X**3+0.2D0/(Z+1.D0)*X**4
RETURN
END FUNCTION F31
REAL(8) FUNCTION F32EE(X,Z)
IMPLICIT NONE
REAL(8) X,Z
F32EE = (0.05D0+0.62D0*Z)/(Z*(1.D0+0.44D0*Z))*(X-X**4)&
& +1.D0/(1.D0+0.22D0*Z)*(X**2-X**4-1.2D0*(X**3-X**4))+1.2D0/(1.D0+0.5D0*Z)*X**4
RETURN
END FUNCTION F32EE
REAL(8) FUNCTION F32EI(X,Z)
IMPLICIT NONE
REAL(8) X,Z
F32EI =-(0.56D0+1.93D0*Z)/(Z*(1.D0+0.44D0*Z))*(X-X**4) &
& +4.95D0/(1.D0+2.48D0*Z)*(X**2-X**4-0.55D0*(X**3-X**4))-1.2D0/(1.D0+0.5D0*Z)*X**4
RETURN
END FUNCTION F32EI
! ************************************************
! BOOTSTRAP CURRENT (Hirshman (Wilson))
! ************************************************
SUBROUTINE TRAJBSNEW
USE TRCOMM, ONLY : AJBS, BB, DR, EPSRHO, NRMAX, NSM, PADD, PBSCD, PNSS, PTS, PZ, RDP, RHOG, RHOM, RKEV, RN, RT, &
& RW, TTRHOG
IMPLICIT NONE
INTEGER(4):: NR, NS
REAL(8) :: DDD, DDX, DPE, DPI, DRL, DTE, DTI, EPS, FT, PE, PPI, RL31, RL32, RNM, RNP, RNTM, RNTP, RPIM, RPIP, TE, TI
REAL(8),DIMENSION(NRMAX):: AJBSL, ANI
REAL(8) :: DERIV3P
IF(PBSCD.LE.0.D0) RETURN
DO NR=1,NRMAX-1
EPS=EPSRHO(NR)
! EPSS=SQRT(EPS)**3
! QL=ABS(QP(NR))
! ZEFFL=0.5D0*(ZEFF(NR)+ZEFF(NR+1))
DRL=1.D0/DR
RNTP=0.D0
RNP= 0.D0
RNTM=0.D0
RNM =0.D0
DO NS=2,NSM
RNTP=RNTP+RN(NR+1,NS)*RT(NR+1,NS)
RNP =RNP +RN(NR+1,NS)
RNTM=RNTM+RN(NR ,NS)*RT(NR ,NS)
RNM =RNM +RN(NR ,NS)
ENDDO
RNTP=RNTP+RW(NR+1,1)+RW(NR+1,2)
RNTM=RNTM+RW(NR ,1)+RW(NR ,2)
RPIP=RNTP+PADD(NR+1)
RPIM=RNTM+PADD(NR )
! ****** ION PARAMETER ******
! ***** ANI is the the ion density (ni) *****
! ***** TI is the ion temperature (Ti) *****
! ***** DTI is the derivative of ion temperature (dTi/dr) *****
! ***** PPI is the ion pressure (Pi) *****
! ***** DPI is the derivative of ion pressure (dPi/dr) *****
! ***** VTI is the ion velocity (VTi) *****
ANI(NR)=0.D0
DO NS=2,NSM
ANI(NR)=ANI(NR)+0.5D0*PZ(NS)*(RN(NR+1,NS)+RN(NR,NS))
ENDDO
PPI=0.5D0*(RPIP+RPIM)
DPI=(RPIP-RPIM)*DRL
TI =0.5D0*(RNTP/RNP+RNTM/RNM)
DTI=(RNTP/RNP-RNTM/RNM)*DRL
! VTI=SQRT(ABS(TI)*RKEV/AMM)
! ANE=0.5D0*(RN(NR+1,1)+RN(NR,1))
! rLnLam=17.3D0-DLOG(ANE)*0.5D0+DLOG(ABS(TI))*1.5D0
! TAUI=12.D0*PI*SQRT(PI)*EPS0**2*SQRT(AMM)
! & *(ABS(TI)*RKEV)**1.5D0/(ANI(NR)*1.D20
! & *ZEFFL**4*AEE**4*rLnLam)
! RNUI=QL*RR/(TAUI*VTI*EPSS)
! ****** ELECTORON PARAMETER ******
! ***** ANE is the the electron density (ne) *****
! ***** TE is the electron temperature (Te) *****
! ***** DTE is the derivative of electron temperature (dTe/dr) ****
! ***** PE is the electron pressure (Pe) *****
! ***** DPE is the derivative of electron pressure (dPe/dr) *****
! ***** VTE is the electron velocity (VTe) *****
TE= 0.5D0*(RT(NR+1,1)+RT(NR,1))
PE= 0.5D0*(RN(NR+1,1)*RT(NR+1,1)+RN(NR,1)*RT(NR,1))
DTE=(RT(NR+1,1)-RT(NR,1))*DRL
DPE=(RN(NR+1,1)*RT(NR+1,1)-RN(NR,1)*RT(NR,1))*DRL
! VTE=SQRT(ABS(TE)*RKEV/AME)
! rLnLam=15.2D0-DLOG(ANE)*0.5D0+DLOG(ABS(TE))
! TAUE=6.D0*PI*SQRT(2.D0*PI)*EPS0**2*SQRT(AME)*(ABS(TE)*RKEV)**1.5D0/(ANI*1.D20*ZEFFL**2*AEE**4*rLnLam)
! RNUE=QL*RR/(TAUE*VTE*EPSS)
! FT=FTPF(MDLTPF,EPS)
FT=(1.46D0*SQRT(EPS)+2.4D0*EPS)/(1.D0-EPS)**1.5D0
! DDX=2.4D0+5.4D0*FT+2.6D0*FT**2
! DDD=-1.17D0/(1.D0+0.46D0*FT)
! C1=(4.D0+2.6D0*FT)/((1.D0+1.02D0*DSQRT(RNUE)+1.07D0*RNUE)*(1.D0+1.07D0*EPSS*RNUE))
! C2=C1*TI/TE
! C3=(7.D0+6.5D0*FT)/((1.D0+0.57D0*DSQRT(RNUE)+0.61D0*RNUE)*(1.D0+0.61D0*EPSS*RNUE))-2.5D0*C1
! C4=((DDD+0.35D0*DSQRT(RNUI))/(1.D0+0.7D0*DSQRT(RNUI))
! & +2.1D0*EPS**3*RNUI**2)*C2/((1.D0-EPS**3*RNUI**2)*(1.D0+EPS**3*RNUE**2))
! AJBSL(NR)=-PBSCD*FT*PE*1.D20*RKEV/(DDX*BP(NR))*(C1*(DPE/PE)
! & +C2*(DPI/PPI)+C3*(DTE/TE)+C4*(DTI/TI))
! *** S. P. Hirshman, Phys Fluids 31, 1988 3150 ***
! *** cited (H.R. WILSON, Nucl.Fusion 32,no.2,1992 259-263) ***
DDX=1.414D0*PZ(2)+PZ(2)**2+FT*(0.754D0+2.657D0*PZ(2) &
& +2.D0*PZ(2)**2)+FT**2*(0.348D0+1.243D0*PZ(2)+PZ(2)**2)
RL31= FT*(0.754D0+2.210D0*PZ(2)+PZ(2)**2 +FT*(0.348D0+1.243D0*PZ(2)+PZ(2)**2))/DDX
RL32=-FT*(0.884D0+2.074D0*PZ(2))/DDX
DDD=-1.172D0/(1.D0+0.462D0*FT)
AJBSL(NR)=-PBSCD*TTRHOG(NR)*PE*1.D20*RKEV*(RL31*((DPE/PE)+(TI/(PZ(2)*TE)) &
& *((DPI/PPI)+DDD*(DTI/TI)))+RL32*(DTE/TE))/RDP(NR)/BB
ENDDO
NR=NRMAX
EPS=EPSRHO(NR)
! EPSS=SQRT(EPS)**3
! QL=ABS(QP(NR))
! ZEFFL=2.D0*ZEFF(NR-1)-ZEFF(NR-2)
DRL=1.D0/DR
! In the following, we assume that
! 1. pressures of beam and fusion at rho=1 are negligible,
! 2. calibration of Pbeam (from UFILE) is not necessary at rho=1.
RNTP=0.D0
RNP= 0.D0
RNTM=0.D0
RNM =0.D0
DO NS=2,NSM
RNTP=RNTP+PNSS(NS)*PTS(NS)
RNP =RNP +PNSS(NS)
RNTM=RNTM+RN(NR-1,NS)*RT(NR-1,NS)+RN(NR ,NS)*RT(NR ,NS)
RNM =RNM +RN(NR-1,NS)+RN(NR ,NS)
ENDDO
RNTM=RNTP+RW(NR-1,1)+RW(NR-1,2)+RW(NR ,1)+RW(NR ,2)
RPIP=RNTP
RPIM=RNTM+PADD(NR-1)+PADD(NR )
RNTM=0.5D0*RNTM
RNM =0.5D0*RNM
RPIM=0.5D0*RPIM
! ****** ION PARAMETER ******
! ***** ANI is the the ion density (ni) *****
! ***** TI is the ion temperature (Ti) *****
! ***** DTI is the derivative of ion temperature (dTi/dr) *****
! ***** PPI is the ion pressure (Pi) *****
! ***** DPI is the derivative of ion pressure (dPi/dr) *****
! ***** VTI is the ion velocity (VTi) *****
ANI(NR)=0.D0
DO NS=2,NSM
ANI(NR)=ANI(NR)+PZ(NS)*PNSS(NS)
ENDDO
PPI=0.5D0*(RPIP+RPIM)
DPI=(RPIP-RPIM)*DRL
TI =0.5D0*(RNTP/RNP+RNTM/RNM)
DTI=(RNTP/RNP-RNTM/RNM)*DRL
! VTI=SQRT(ABS(TI)*RKEV/AMM)
! ANE=PNSS(1)
! rLnLam=17.3D0-DLOG(ANE)*0.5D0+DLOG(ABS(TI))*1.5D0
! TAUI=12.D0*PI*SQRT(PI)*EPS0**2*SQRT(AMM)
! & *(ABS(TI)*RKEV)**1.5D0/(ANI(NR)*1.D20
! & *ZEFFL**4*AEE**4*rLnLam)
! RNUI=QL*RR/(TAUI*VTI*EPSS)
! ****** ELECTORON PARAMETER ******
! ***** ANE is the the electron density (ne) *****
! ***** TE is the electron temperature (Te) *****
! ***** DTE is the derivative of electron temperature (dTe/dr) ****
! ***** PE is the electron pressure (Pe) *****
! ***** DPE is the derivative of electron pressure (dPe/dr) *****
! ***** VTE is the electron velocity (VTe) *****
TE =PTS(1)
PE =PNSS(1)*PTS(1)
DTE=DERIV3P(PTS(1),RT(NR,1),RT(NR-1,1),RHOG(NR),RHOM(NR),RHOM(NR-1))
DPE=DERIV3P(PNSS(1)*PTS(1),RN(NR,1)*RT(NR,1),RN(NR-1,1)*RT(NR-1,1),RHOG(NR),RHOM(NR),RHOM(NR-1))
! VTE=SQRT(ABS(TE)*RKEV/AME)
! rLnLam=15.2D0-DLOG(ANE)*0.5D0+DLOG(ABS(TE))
! TAUE=6.D0*PI*SQRT(2.D0*PI)*EPS0**2*SQRT(AME)
! & *(ABS(TE)*RKEV)**1.5D0/(ANI*1.D20*ZEFFL**2*AEE**4*rLnLam)
! RNUE=QL*RR/(TAUE*VTE*EPSS)
! FT=FTPF(MDLTPF,EPS)
FT=(1.46D0*SQRT(EPS)+2.4D0*EPS)/(1.D0-EPS)**1.5D0
! DDX=2.4D0+5.4D0*FT+2.6D0*FT**2
! DDD=-1.17D0/(1.D0+0.46D0*FT)
! C1=(4.D0+2.6D0*FT)/((1.D0+1.02D0*DSQRT(RNUE)+1.07D0*RNUE)*(1.D0+1.07D0*EPSS*RNUE))
! C2=C1*TI/TE
! C3=(7.D0+6.5D0*FT)/((1.D0+0.57D0*DSQRT(RNUE)+0.61D0*RNUE)*(1.D0+0.61D0*EPSS*RNUE))
! & -2.5D0*C1
! C4=((DDD+0.35D0*DSQRT(RNUI))/(1.D0+0.7D0*DSQRT(RNUI))+2.1D0*EPS**3*RNUI**2)*C2
! & /((1.D0-EPS**3*RNUI**2)*(1.D0+EPS**3*RNUE**2))
! AJBSL(NR)=-PBSCD*FT*PE*1.D20*RKEV/(DDX*BP(NR))*(C1*(DPE/PE)
! & +C2*(DPI/PPI)+C3*(DTE/TE)+C4*(DTI/TI))
! *** S. P. Hirshman, Phys Fluids 31, 1988 3150 ***
! *** cited (H.R. WILSON, Nucl.Fusion 32,no.2,1992 259-263) ***
DDX=1.414D0*PZ(2)+PZ(2)**2+FT*(0.754D0+2.657D0*PZ(2) &
& +2.D0*PZ(2)**2)+FT**2*(0.348D0+1.243D0*PZ(2)+PZ(2)**2)
RL31=FT*( 0.754D0+2.21D0*PZ(2)+PZ(2)**2+FT*(0.348D0+1.243D0 &
& *PZ(2)+PZ(2)**2))/DDX
RL32=-FT*(0.884D0+2.074D0*PZ(2))/DDX
DDD=-1.172D0/(1.D0+0.462D0*FT)
AJBSL(NR)=-PBSCD*TTRHOG(NR)*PE*1.D20*RKEV &
& *(RL31*((DPE/PE)+(TI/(PZ(2)*TE))*((DPI/PPI)+DDD*(DTI/TI)))+RL32*(DTE/TE))/RDP(NR)/BB
AJBS(1)=0.5D0*AJBSL(1)
DO NR=2,NRMAX
AJBS(NR)=0.5D0*(AJBSL(NR)+AJBSL(NR-1))
ENDDO
RETURN
END SUBROUTINE TRAJBSNEW
! ***********************************************************
! BOOTSTRAP CURRENT (Hinton & Hazeltine)
! ***********************************************************
SUBROUTINE TRAJBS
USE TRCOMM, ONLY : AJBS, AME, AMM, BB, BP, DR, EPSRHO, NRMAX, NSMAX, PA, PBSCD, PNSS, PTS, PZ, QP, RHOG, RHOM, &
& RJCB, RKEV, RN, RR, RT, ZEFF
IMPLICIT NONE
INTEGER(4):: NR
REAL(8) :: A, AMA, AMD, AMT, ANA, ANDX, ANE, ANT, BPL, DPA, DPD, DPE, DPT, DRL, DTA, DTD, DTE, DTT, EPS, EPSS, FACT, &
& FTAUE, FTAUI, H, PAL, PDL, PEL, PTL, RK13E, RK23E, RK3A, RK3D, RK3T, RNUA, RNUD, RNUE, RNUT, TAL, TAUA, &
& TAUD, TAUE, TAUT, TDL, TEL, TTL, VTA, VTD, VTE, VTT, ZEFFL
REAL(8) :: DERIV3P
REAL(8),DIMENSION(NRMAX):: AJBSL
REAL(8):: RK13=2.30D0, RA13=1.02D0, RB13=1.07D0, RC13=1.07D0, RK23=4.19D0, RA23=0.57D0, RB23=0.61D0, RC23=0.61D0
! ZEFF=1
! DATA RK11,RA11,RB11,RC11/1.04D0,2.01D0,1.53D0,0.89D0/
! DATA RK12,RA12,RB12,RC12/1.20D0,0.76D0,0.67D0,0.56D0/
! DATA RK22,RA22,RB22,RC22/2.55D0,0.45D0,0.43D0,0.43D0/
!! DATA RK13,RA13,RB13,RC13/2.30D0,1.02D0,1.07D0,1.07D0/
!! DATA RK23,RA23,RB23,RC23/4.19D0,0.57D0,0.61D0,0.61D0/
! DATA RK33,RA33,RB33,RC33/1.83D0,0.68D0,0.32D0,0.66D0/
! DATA RK2 ,RA2 ,RB2 ,RC2 /0.66D0,1.03D0,0.31D0,0.74D0/
IF(PBSCD.LE.0.D0) RETURN
AMD=PA(2)*AMM
AMT=PA(3)*AMM
AMA=PA(4)*AMM
DO NR=1,NRMAX-1
EPS=EPSRHO(NR)
EPSS=SQRT(EPS)**3
ANE=0.5D0*(RN(NR+1,1)+RN(NR,1))
ANDX=0.5D0*(RN(NR+1,2)+RN(NR,2))
ANT =0.5D0*(RN(NR+1,3)+RN(NR,3))
ANA =0.5D0*(RN(NR+1,4)+RN(NR,4))
TEL=ABS(0.5D0*(RT(NR+1,1)+RT(NR,1)))
TDL=ABS(0.5D0*(RT(NR+1,2)+RT(NR,2)))
TTL=ABS(0.5D0*(RT(NR+1,3)+RT(NR,3)))
TAL=ABS(0.5D0*(RT(NR+1,4)+RT(NR,4)))
PEL=0.5D0*(RN(NR+1,1)*RT(NR+1,1)+RN(NR,1)*RT(NR,1))
PDL=0.5D0*(RN(NR+1,2)*RT(NR+1,2)+RN(NR,2)*RT(NR,2))
PTL=0.5D0*(RN(NR+1,3)*RT(NR+1,3)+RN(NR,3)*RT(NR,3))
PAL=0.5D0*(RN(NR+1,4)*RT(NR+1,4)+RN(NR,4)*RT(NR,4))
ZEFFL=0.5D0*(ZEFF(NR+1)+ZEFF(NR))
TAUE = FTAUE(ANE,ANDX,TEL,ZEFFL)
TAUD = FTAUI(ANE,ANDX,TDL,PZ(2),PA(2))
TAUT = FTAUI(ANE,ANT ,TTL,PZ(3),PA(3))
TAUA = FTAUI(ANE,ANA ,TAL,PZ(4),PA(4))
VTE=SQRT(TEL*RKEV/AME)
VTD=SQRT(TDL*RKEV/AMD)
VTT=SQRT(TTL*RKEV/AMT)
VTA=SQRT(TAL*RKEV/AMA)
RNUE=ABS(QP(NR))*RR/(TAUE*VTE*EPSS)
RNUD=ABS(QP(NR))*RR/(TAUD*VTD*EPSS)
RNUT=ABS(QP(NR))*RR/(TAUT*VTT*EPSS)
RNUA=ABS(QP(NR))*RR/(TAUA*VTA*EPSS)
! RK11E=RK11*(1.D0/(1.D0+RA11*SQRT(RNUE)+RB11*RNUE)+(EPSS*RC11)**2/RB11*RNUE/(1.D0+RC11*RNUE*EPSS))
! RK12E=RK12*(1.D0/(1.D0+RA12*SQRT(RNUE)+RB12*RNUE)+(EPSS*RC12)**2/RB12*RNUE/(1.D0+RC12*RNUE*EPSS))
! RK22E=RK22*(1.D0/(1.D0+RA22*SQRT(RNUE)+RB22*RNUE)+(EPSS*RC22)**2/RB22*RNUE/(1.D0+RC22*RNUE*EPSS))
RK13E=RK13/(1.D0+RA13*SQRT(RNUE)+RB13*RNUE)/(1.D0+RC13*RNUE*EPSS)
RK23E=RK23/(1.D0+RA23*SQRT(RNUE)+RB23*RNUE)/(1.D0+RC23*RNUE*EPSS)
! RK33E=RK33/(1.D0+RA33*SQRT(RNUE)+RB33*RNUE)/(1.D0+RC33*RNUE*EPSS)
! RK2D =RK2 *(1.D0/(1.D0+RA2 *SQRT(RNUD)+RB2 *RNUD)+(EPSS*RC2 )**2/RB2 *RNUD/(1.D0+RC2 *RNUD*EPSS))
! RK2T =RK2 *(1.D0/(1.D0+RA2 *SQRT(RNUT)+RB2 *RNUT)+(EPSS*RC2 )**2/RB2 *RNUT/(1.D0+RC2 *RNUT*EPSS))
! RK2A =RK2 *(1.D0/(1.D0+RA2 *SQRT(RNUA)+RB2 *RNUA)+(EPSS*RC2 )**2/RB2 *RNUA/(1.D0+RC2 *RNUA*EPSS))
RK3D=((1.17D0-0.35D0*SQRT(RNUD))/(1.D0+0.7D0*SQRT(RNUD))-2.1D0*(RNUD*EPSS)**2)/(1.D0+(RNUD*EPSS)**2)
RK3T=((1.17D0-0.35D0*SQRT(RNUT))/(1.D0+0.7D0*SQRT(RNUT))-2.1D0*(RNUT*EPSS)**2)/(1.D0+(RNUT*EPSS)**2)
RK3A=((1.17D0-0.35D0*SQRT(RNUA))/(1.D0+0.7D0*SQRT(RNUA))-2.1D0*(RNUA*EPSS)**2)/(1.D0+(RNUA*EPSS)**2)
DRL=RJCB(NR)/DR
DTE=(RT(NR+1,1)-RT(NR,1))*DRL
DTD=(RT(NR+1,2)-RT(NR,2))*DRL
DTT=(RT(NR+1,3)-RT(NR,3))*DRL
DTA=(RT(NR+1,4)-RT(NR,4))*DRL
DPE=(RN(NR+1,1)*RT(NR+1,1)-RN(NR,1)*RT(NR,1))*DRL
DPD=(RN(NR+1,2)*RT(NR+1,2)-RN(NR,2)*RT(NR,2))*DRL
DPT=(RN(NR+1,3)*RT(NR+1,3)-RN(NR,3)*RT(NR,3))*DRL
DPA=(RN(NR+1,4)*RT(NR+1,4)-RN(NR,4)*RT(NR,4))*DRL
BPL=BP(NR)
FACT=1.D0/(1.D0+(RNUE*EPS)**2)
IF(NSMAX.EQ.2) THEN
A= DPE/PEL- 2.5D0*DTE/TEL +(PDL/PEL)*(DPD/PDL-RK3D*FACT*DTD/TDL)
ELSEIF(NSMAX.EQ.3) THEN
A= DPE/PEL- 2.5D0*DTE/TEL &
& +(PDL/PEL)*(DPD/PDL-RK3D*FACT*DTD/TDL)+(PTL/PEL)*(DPT/PTL-RK3T*FACT*DTT/TTL)
ELSEIF(NSMAX.EQ.4) THEN
A= DPE/PEL- 2.5D0*DTE/TEL +(PDL/PEL)*(DPD/PDL-RK3D*FACT*DTD/TDL) &
& +(PTL/PEL)*(DPT/PTL-RK3T*FACT*DTT/TTL)+(PAL/PEL)*(DPA/PAL-RK3A*FACT*DTA/TAL)
ENDIF
H=BB/(BB+BP(NR))
AJBSL(NR)=-PBSCD*H*SQRT(EPS)*ANE*1.D20*TEL*RKEV/BPL *(RK13E*A+RK23E*DTE/TEL)
!
! WRITE(6,'(I3,1P6E12.4)') NR,RNUE,RK13E,RK23E,A,BPL,AJBSL(NR)
! IF(NR.EQ.NRMAX) THEN
! WRITE(6,'(1P6E12.4)') DN,DTE,ANE,PNSS(1),TE,PTS(1)
! WRITE(6,'(1P6E12.4)') RK13E*(RN(NR,1)*RT(NR,1)
! & +RN(NR,2)*RT(NR,2)
! & +RN(NR,3)*RT(NR,3)
! & +RN(NR,4)*RT(NR,4)
! & +RW(NR,1)
! & +RW(NR,2) )*DN/ANE,
! & +(RK23E-1.5D0*RK13E)*RN(NR,1)*DTE,
! & +RK13D*(RK23D-1.5D0)*RN(NR,2)*DTD,
! & BPL,AJBS(NR-1),AJBS(NR)
! ENDIF
ENDDO
NR=NRMAX
EPS=EPSRHO(NR)
EPSS=SQRT(EPS)**3
ANE=PNSS(1)
ANDX=PNSS(2)
! ANT =PNSS(3)
! ANA =PNSS(4)
TEL=ABS(PTS(1))
TDL=ABS(PTS(2))
TTL=ABS(PTS(3))
TAL=ABS(PTS(4))
PEL=PNSS(1)*PTS(1)
PDL=PNSS(2)*PTS(2)
PTL=PNSS(3)*PTS(3)
PAL=PNSS(4)*PTS(4)
ZEFFL=2.D0*ZEFF(NR-1)-ZEFF(NR-2)
TAUE = FTAUE(ANE,ANDX,TEL,ZEFFL)
TAUD = FTAUI(ANE,ANDX,TDL,PZ(2),PA(2))
TAUT = FTAUI(ANE,ANT ,TTL,PZ(3),PA(3))
TAUA = FTAUI(ANE,ANA ,TAL,PZ(4),PA(4))
VTE=SQRT(TEL*RKEV/AME)
VTD=SQRT(TDL*RKEV/AMD)
VTT=SQRT(TTL*RKEV/AMT)
VTA=SQRT(TAL*RKEV/AMA)
RNUE=ABS(QP(NR))*RR/(TAUE*VTE*EPSS)
RNUD=ABS(QP(NR))*RR/(TAUD*VTD*EPSS)
RNUT=ABS(QP(NR))*RR/(TAUT*VTT*EPSS)
RNUA=ABS(QP(NR))*RR/(TAUA*VTA*EPSS)
! RK11E=RK11*(1.D0/(1.D0+RA11*SQRT(RNUE)+RB11*RNUE)+(EPSS*RC11)**2/RB11*RNUE/(1.D0+RC11*RNUE*EPSS))
! RK12E=RK12*(1.D0/(1.D0+RA12*SQRT(RNUE)+RB12*RNUE)+(EPSS*RC12)**2/RB12*RNUE/(1.D0+RC12*RNUE*EPSS))
! RK22E=RK22*(1.D0/(1.D0+RA22*SQRT(RNUE)+RB22*RNUE)+(EPSS*RC22)**2/RB22*RNUE/(1.D0+RC22*RNUE*EPSS))
RK13E=RK13/(1.D0+RA13*SQRT(RNUE)+RB13*RNUE) /(1.D0+RC13*RNUE*EPSS)
RK23E=RK23/(1.D0+RA23*SQRT(RNUE)+RB23*RNUE) /(1.D0+RC23*RNUE*EPSS)
! RK33E=RK33/(1.D0+RA33*SQRT(RNUE)+RB33*RNUE)/(1.D0+RC33*RNUE*EPSS)
! RK2D =RK2 *(1.D0/(1.D0+RA2 *SQRT(RNUD)+RB2 *RNUD)+(EPSS*RC2 )**2/RB2 *RNUD/(1.D0+RC2 *RNUD*EPSS))
! RK2T =RK2 *(1.D0/(1.D0+RA2 *SQRT(RNUT)+RB2 *RNUT)+(EPSS*RC2 )**2/RB2 *RNUT/(1.D0+RC2 *RNUT*EPSS))
! RK2A =RK2 *(1.D0/(1.D0+RA2 *SQRT(RNUA)+RB2 *RNUA)+(EPSS*RC2 )**2/RB2 *RNUA/(1.D0+RC2 *RNUA*EPSS))
RK3D=((1.17D0-0.35D0*SQRT(RNUD))/(1.D0+0.7D0*SQRT(RNUD))-2.1D0*(RNUD*EPSS)**2)/(1.D0+(RNUD*EPSS)**2)
RK3T=((1.17D0-0.35D0*SQRT(RNUT))/(1.D0+0.7D0*SQRT(RNUT))-2.1D0*(RNUT*EPSS)**2)/(1.D0+(RNUT*EPSS)**2)
RK3A=((1.17D0-0.35D0*SQRT(RNUA))/(1.D0+0.7D0*SQRT(RNUA))-2.1D0*(RNUA*EPSS)**2)/(1.D0+(RNUA*EPSS)**2)
DRL=RJCB(NR)/DR
DTE=DERIV3P(PTS(1),RT(NR,1),RT(NR-1,1),RHOG(NR),RHOM(NR),RHOM(NR-1))
DTD=DERIV3P(PTS(2),RT(NR,2),RT(NR-1,2),RHOG(NR),RHOM(NR),RHOM(NR-1))
DTT=DERIV3P(PTS(3),RT(NR,3),RT(NR-1,3),RHOG(NR),RHOM(NR),RHOM(NR-1))
DTA=DERIV3P(PTS(4),RT(NR,4),RT(NR-1,4),RHOG(NR),RHOM(NR),RHOM(NR-1))
DPE=DERIV3P(PNSS(1)*PTS(1),RN(NR,1)*RT(NR,1),RN(NR-1,1)*RT(NR-1,1),RHOG(NR),RHOM(NR),RHOM(NR-1))
DPD=DERIV3P(PNSS(2)*PTS(2),RN(NR,2)*RT(NR,2),RN(NR-1,2)*RT(NR-1,2),RHOG(NR),RHOM(NR),RHOM(NR-1))
DPT=DERIV3P(PNSS(3)*PTS(3),RN(NR,3)*RT(NR,3),RN(NR-1,3)*RT(NR-1,3),RHOG(NR),RHOM(NR),RHOM(NR-1))
DPA=DERIV3P(PNSS(4)*PTS(4),RN(NR,4)*RT(NR,4),RN(NR-1,4)*RT(NR-1,4),RHOG(NR),RHOM(NR),RHOM(NR-1))
BPL=BP(NR)
FACT=1.D0/(1.D0+(RNUE*EPS)**2)
IF(NSMAX.EQ.2) THEN
A= DPE/PEL- 2.5D0*DTE/TEL +(PDL/PEL)*(DPD/PDL-RK3D*FACT*DTD/TDL)
ELSEIF(NSMAX.EQ.3) THEN
A= DPE/PEL- 2.5D0*DTE/TEL &
& +(PDL/PEL)*(DPD/PDL-RK3D*FACT*DTD/TDL)+(PTL/PEL)*(DPT/PTL-RK3T*FACT*DTT/TTL)
ELSEIF(NSMAX.EQ.4) THEN
A= DPE/PEL- 2.5D0*DTE/TEL +(PDL/PEL)*(DPD/PDL-RK3D*FACT*DTD/TDL) &
& +(PTL/PEL)*(DPT/PTL-RK3T*FACT*DTT/TTL)+(PAL/PEL)*(DPA/PAL-RK3A*FACT*DTA/TAL)
ENDIF
H=BB/(BB+BP(NR))
AJBSL(NR)=-PBSCD*H*SQRT(EPS)*ANE*1.D20*TEL*RKEV/BPL*(RK13E*A+RK23E*DTE/TEL)
!
! WRITE(6,'(I3,1P6E12.4)') NR,RNUE,RK13E,RK23E,A,BPL,AJBSL(NR)
! IF(NR.EQ.NRMAX) THEN
! WRITE(6,'(1P6E12.4)') DN,DTE,ANE,PNSS(1),TE,PTS(1)
! WRITE(6,'(1P6E12.4)') RK13E*(RN(NR,1)*RT(NR,1)
! & +RN(NR,2)*RT(NR,2)
! & +RN(NR,3)*RT(NR,3)
! & +RN(NR,4)*RT(NR,4)
! & +RW(NR,1)
! & +RW(NR,2) )*DN/ANE,
! & +(RK23E-1.5D0*RK13E)*RN(NR,1)*DTE,
! & +RK13D*(RK23D-1.5D0)*RN(NR,2)*DTD,
! & BPL,AJBS(NR-1),AJBS(NR)
! ENDIF
AJBS(1)=0.5D0*AJBSL(1)
DO NR=2,NRMAX
AJBS(NR)=0.5D0*(AJBSL(NR)+AJBSL(NR-1))
ENDDO
RETURN
END SUBROUTINE TRAJBS
! ***********************************************************
! CALCULATE AJ, AJOH, POH, EZOH
! ***********************************************************
SUBROUTINE TRAJOH
USE TRCOMM, ONLY : ABRHOG, AJ, AJBS, AJNB, AJOH, AJRF, AJTOR, BB, DR, DVRHO, DVRHOG, ETA, EZOH, KUFDEV, MDLEQB, &
& MDLJQ, MDLUF, NRMAX, POH, RDP, RMU0, RR, TTRHO, TTRHOG, rm, ABVRHOG, RDPVRHOG, PI, rm, rg
IMPLICIT NONE
INTEGER(4):: NR
REAL(8) :: FACTOR0, FACTORM, FACTORP
IF(MDLEQB.EQ.1.OR.MDLJQ.EQ.1.OR.(MDLUF.EQ.0.OR.MDLUF.EQ.3)) THEN
NR=1
FACTOR0=TTRHO(NR)**2/(RMU0*BB*DVRHO(NR))
FACTORP=ABVRHOG(NR )/TTRHOG(NR )
AJ(NR) =FACTOR0*FACTORP*RDPVRHOG(NR)/DR
DO NR=2,NRMAX
FACTOR0=TTRHO(NR)**2/(RMU0*BB*DVRHO(NR))
FACTORM=ABVRHOG(NR-1)/TTRHOG(NR-1)
FACTORP=ABVRHOG(NR )/TTRHOG(NR )
AJ(NR) =FACTOR0*(FACTORP*RDPVRHOG(NR)-FACTORM*RDPVRHOG(NR-1))/DR
ENDDO
NR=1
FACTOR0=RR/(RMU0*DVRHO(NR))
FACTORP=ABVRHOG(NR )
AJTOR(NR) =FACTOR0*FACTORP*RDPVRHOG(NR)/DR
DO NR=2,NRMAX
FACTOR0=RR/(RMU0*DVRHO(NR))
FACTORM=ABVRHOG(NR-1)
FACTORP=ABVRHOG(NR )
AJTOR(NR) =FACTOR0*(FACTORP*RDPVRHOG(NR)-FACTORM*RDPVRHOG(NR-1))/DR
ENDDO
ENDIF
AJOH(1:NRMAX) = AJ(1:NRMAX)-(AJNB(1:NRMAX)+AJRF(1:NRMAX)+AJBS(1:NRMAX))
EZOH(1:NRMAX) = ETA(1:NRMAX)*AJOH(1:NRMAX)
!!! IF(KUFDEV.EQ.'lhd') THEN
!!! POH(1:NRMAX) = 0.D0
!!! ELSE
POH(1:NRMAX) = EZOH(1:NRMAX)*AJOH(1:NRMAX)
!!! ENDIF
RETURN
END SUBROUTINE TRAJOH
! ***********************************************************
! SAWTOOTH OSCILLATION
! ***********************************************************
SUBROUTINE TRSAWT
USE TRCOMM, ONLY : AR1RHOG, ARRHOG, BP, DR, DVRHO, DVRHOG, MDLST, NRMAX, NSM, PI, QP, RDP, RG, RM, RN, RR, RT, &
& T, TTRHOG, RDPVRHOG
IMPLICIT NONE
INTEGER(4):: IONE, IZEROX, LN, LQ, LT, NR, NS
REAL(8) :: RNN, RTN, SUML, SUML1, SUML2
REAL(8),DIMENSION(NRMAX):: QONE
IF(MDLST.EQ.0) RETURN
IF(MOD(MDLST,2).EQ.1) THEN
LT = 1
ELSE
LT = 0
ENDIF
IF(MOD(MDLST/2,2).EQ.1) THEN
LN = 1
ELSE
LN = 0
ENDIF
IF(MOD(MDLST/4,2).EQ.1) THEN
LQ = 1
ELSE
LQ = 0
ENDIF
WRITE(6,601) MDLST,T
601 FORMAT(' ','# SAWTOOTH OSCILLATION -TYPE ',I1, ' AT ',F7.3,' SEC')
SUML=0.D0
DO NR=1,NRMAX
SUML = SUML+(1.D0/ABS(QP(NR))-1.D0)*(DVRHOG(NR)/2.D0*PI*RR)*DR
IF(SUML.LT.0.D0) GOTO 1000
ENDDO
NR=NRMAX
1000 IZEROX=NR
DO NR=1,NRMAX
IF(QP(NR).GE.1.D0) GOTO 2000
ENDDO
NR=NRMAX
2000 IONE=NR
SUML = SUM((1.D0/ABS(QP(IONE:IZEROX))-1.D0)*DVRHOG(IONE:IZEROX))/(2.D0*PI*RR)*DR
QONE(1:IZEROX) = 1.D0+SUML*4.D0*RG(1:IZEROX)**2/RG(IZEROX)**4
IF(LT.EQ.1) THEN
DO NS=1,NSM
SUML1 = SUM(RN(1:IZEROX,NS) *DVRHO(1:IZEROX))
SUML2 = SUM(RN(1:IZEROX,NS)*RT(1:IZEROX,NS)*DVRHO(1:IZEROX))
RTN = SUML2/SUML1
RT(1:IZEROX,NS) = RTN
ENDDO
ENDIF
IF(LN.EQ.1) THEN
DO NS=1,NSM
SUML1 = SUM( DVRHO(1:IZEROX))
SUML2 = SUM(RN(1:IZEROX,NS)*DVRHO(1:IZEROX))
RNN = SUML2/SUML1
RN(1:IZEROX,NS) = RNN
ENDDO
ENDIF
IF(LQ.EQ.1) THEN
QP(1:IZEROX) = 1.D0/QONE(1:IZEROX)
RDPVRHOG(1:IZEROX)=TTRHOG(1:IZEROX)*ARRHOG(1:IZEROX) &
& /(4.D0*PI**2*QP(1:IZEROX))
RDP(1:IZEROX)=RDPVRHOG(1:IZEROX)*DVRHOG(1:IZEROX)
BP(1:IZEROX) =AR1RHOG(1:IZEROX)*RDP(1:IZEROX)/RR
ENDIF
WRITE(6,602) RM(IONE),RM(IZEROX),RTN,RNN
602 FORMAT(' ',' R-ONE,R-ZERO,RTN,RNN = ',4F8.3)
DO NR=1,IZEROX+2
WRITE(6,'(A,I5,1P5E12.4)') 'NR,R,Q,B,Q,S=', &
NR,RM(NR),QP(NR),BP(NR),QONE(NR),SUML
ENDDO
RETURN
END SUBROUTINE TRSAWT
! ***********************************************************
! TRAPPED PARTICLE FRACTION
! ***********************************************************
REAL(8) FUNCTION FTPF(ID,EPS)
IMPLICIT NONE
INTEGER(4) :: ID, IERR, N, I
INTEGER(4),PARAMETER :: IMAX=20
REAL(8) :: EPS, EPSC, FTLL, FTUL, OMEGA, PI, S
REAL(8),DIMENSION(IMAX,IMAX):: TABLE
REAL(8) :: FTL, FTU
EXTERNAL FTL, FTU
IF(ID.EQ.1) THEN
! Y. R. Lin-Liu and R. L. Miller, PoP 2 1666 (1995), eqs(7)(13)(18)(19)
PI=3.14159265358979323846D0
EPSC=1.D-9
I=IMAX
FTUL=1.D0-(1.D0-1.5D0*SQRT(EPS)+0.5D0*EPS**1.5D0)/SQRT(1-EPS**2)
CALL RMBRG(0.D0,2.D0*PI,EPSC,S,I,N,IERR,TABLE,EPS,FTL)
FTLL=1.D0-(1.D0-EPS)**1.5D0/SQRT(1.D0+EPS)*(S/(2.D0*PI))
OMEGA=(3.D0*SQRT(2.D0)/2.D0*0.69D0-3.D0*SQRT(2.D0)/PI)/(1.5D0-3.D0*SQRT(2.D0)/PI)
FTPF=OMEGA*FTUL+(1.D0-OMEGA)*FTLL
ELSEIF(ID.EQ.2) THEN
! S. P. Hirshman et al., NF 17 611 (1977)
FTPF=1.D0-(1.D0-EPS)**2.D0/(DSQRT(1.D0-EPS**2)*(1.D0+1.46D0*DSQRT(EPS)))
ELSEIF(ID.EQ.3) THEN
! Y. R. Lin-Liu and R. L. Miller, PoP 2 1666 (1995), eqs(16)(17)(18)
PI=3.14159265358979323846D0
FTUL=1.5D0*SQRT(EPS)
FTLL=3.D0*SQRT(2.D0)/PI*SQRT(EPS)
FTPF=0.75D0*FTUL+0.25D0*FTLL
ELSEIF(ID.EQ.4) THEN
! M. N. Rosenbluth et al., PoF 15 116 (1972)
FTPF=1.46D0*SQRT(EPS)
ELSE
! Y. B. Kim et al., PoF B 3 2050 (1991) eq(C18), default
FTPF=1.46D0*SQRT(EPS)-0.46D0*(EPS)**1.5D0
ENDIF
RETURN
END FUNCTION FTPF
! *********************************************
! FUNCTION FOR ROMBERG INTEGRATION
! *********************************************
REAL(8) FUNCTION FTU(X,EPS)
IMPLICIT NONE
REAL(8):: EPS, X
FTU = X/SQRT(1.D0-X*(1.D0-EPS))
RETURN
END FUNCTION FTU
REAL(8) FUNCTION FTL(X,EPS)
!
IMPLICIT NONE
REAL(8):: X, EPS
REAL(8):: H
H = (1.D0 - EPS) / (1.D0 + EPS * COS(X))
FTL = (1.D0 - SQRT(1.D0 - H) * (1.D0 + 0.5D0 * H)) / H**2
RETURN
END FUNCTION FTL
!
! *********************************************
! ROMBERG INTEGRATION METHOD
! *********************************************
SUBROUTINE RMBRG(A,B,EPS,S,IMAX,N,IERR,T,ARG,F)
!
! A : lower bound
! B : upper bound
! EPS : stopping criterion
! IMAX : maximum division number
! F : formula of integrand
!