! $Id$ module tx_coefficients use tx_commons use tx_core_module implicit none private real(8), dimension(:,:,:,:), allocatable :: ELM real(8), dimension(:), allocatable :: & & rNuIN0, rNuCXN0, rNubeBE, rNubiBI, rNuTeiEI,& & rNuCXN1, rMueNe, rMuiNi, & !dPNeV, dPNiV, & & RUbthV, UethVR, UithVR, EthVR, RUerV, RUirV, UerVR, UirVR, & & FWpheBB, FWpheBB2, FWthpheB, FWpheB, FWphe, dAphV, & !!ion & FWphiBB, FWphiBB2, FWthphiB, FWphiB, & & BphBNi, BthBNi, Dbrpft, & & rNuei1EI, rNuei2BthEI, rNuei3EI, & & rNube1BE, rNube2BthBE, rNube3BE, & & Vbparaph, RVbparath, & & Chie1, Chie2, Chii1, Chii2, & & FVpchph, rGASPFA, RNeUer, RNiUir, & & rNue2NCN, rNui2NCN, rat_ei, & & rNueHLphthVR, rNuiHLphthVR ! for UHphSwitch == 1 miki_m 2010/8/26 !!rp_conv &, rNubLL real(8), dimension(:), allocatable :: UNITY real(8) :: DTt, DTf(1:NQM), invDT, BeamSW, RpplSW, ThntSW, FSVAHLT, FSVAHLE, & & fact = 1.d0 ! <= SOL loss accelerator integer(4), save :: ICALA = 0 public :: TXCALA contains !*************************************************************** ! ! Calculate ALC, BLC, CLC, PLC ! !*************************************************************** SUBROUTINE TXCALA INTEGER(4) :: NE, NR, NC, NQ, N, iHvsLC !*** Nodal Equation *****************************************! ! ALC : Coefficient matrix on NR+1 ! ! BLC : Coefficient matrix on NR ! ! CLC : Coefficient matrix on NR-1 ! ! NLCR : Matrix of identifying variables with NR index ! ! NLC : Matrix of identifying variables ! ! PLC : Coefficient matrix for non-variable terms ! ! NLCMAX : Number of right-hand-side terms ! !************************************************************! !*** Elemental Equation *************************************! ! ELM : Elemental matrix of one term, equation, element ! !************************************************************! allocate(ELM(1:NEMAX,1:4,0:NCM,1:NQMAX)) ELM(1:NEMAX,1:4,0:NCM,1:NQMAX) = 0.D0 ALC(0:NCM,1:NQMAX,0:NRMAX) = 0.D0 BLC(0:NCM,1:NQMAX,0:NRMAX) = 0.D0 CLC(0:NCM,1:NQMAX,0:NRMAX) = 0.D0 NLCR(0:NCM,1:NQMAX,0:NRMAX) = 0 NLC(0:NCM,1:NQMAX) = 0 PLC(1:NCM,1:NQMAX,0:NRMAX) = 0.D0 NLCMAX(1:NQMAX) = 0 ! Preconditioning invDT = 1.d0 / DT !!$ IF(DT <= 2.D-5) THEN !!$ DTt = 1.5d1 !!$ DTf(1) = 1.d0 !!$ DTf(2:NQMAX) = 1.5d1 !!$ ELSE DTt = 1.d0 DTf(1:NQMAX) = 1.d0 !!$ END IF ! In case of ohmic heating (i.e. no NBI), largeness of the term related to beam ! components comparable to that of other terms related to finite variables ! sometimes induces spurious values in beam components even if no NBI is activated. ! To avoid this, setting these terms, especially related to electromagnetic potentials, ! to zero is desirable during no NBI. IF(PNBH == 0.D0 .and. PNbV(0) == 0.D0) THEN BeamSW = 0.D0 ELSE BeamSW = 1.D0 END IF IF(FSRP == 1.D0) THEN RpplSW = 1.D0 ELSE RpplSW = 0.D0 END IF ! If ThntSW = 0, there is no particle source from PN02V. ThntSW = 1.D0 ! Annihilator of inherent particle pinch induced by FWth terms IF(MDVAHL == 1) THEN FSVAHLT = FSVAHL FSVAHLE = 0.D0 ELSE IF(MDVAHL >= 2) THEN IF(MDVAHL == 2) THEN FSVAHLT = FSVAHL ELSE FSVAHLT = 0.D0 END IF FSVAHLE = 1.D0 ELSE ! MDVAHL == 0 FSVAHLT = 1.D0 FSVAHLE = 0.D0 END IF ! Coefficients N = NRMAX allocate(rNuIN0(0:N), rNuCXN0(0:N), rNubeBE(0:N), rNubiBI(0:N), rNuTeiEI(0:N), & & rNuCXN1(0:N), rMueNe(0:N), rMuiNi(0:N), &!dPNeV(0:N), dPNiV(0:N), & & RUbthV(0:N), UethVR(0:N), UithVR(0:N), & & EthVR(0:N), RUerV(0:N), RUirV(0:N), UerVR(0:N), UirVR(0:N), & & FWpheBB(0:N), FWpheBB2(0:N), FWthpheB(0:N), FWpheB(0:N), FWphe(0:N), & & dAphV(0:N), BphBNi(0:N), BthBNi(0:N), Dbrpft(0:N), & & rNuei1EI(0:N), rNuei2BthEI(0:N), rNuei3EI(0:N), & & rNube1BE(0:N), rNube2BthBE(0:N), rNube3BE(0:N), & & Vbparaph(0:N), RVbparath(0:N), & & Chie1(0:N), Chie2(0:N), Chii1(0:N), Chii2(0:N), & & FVpchph(0:N), rGASPFA(0:N), RNeUer(0:N), RNiUir(0:N), & & rNueHLphthVR(0:N), rNuiHLphthVR(0:N))!, & ! & rNue2NCN(0:N), rNui2NCN(0:N), rat_ei(0:N)) !!ion allocate(FWphiBB(0:N), FWphiBB2(0:N), FWthphiB(0:N), FWphiB(0:N)) allocate(UNITY(0:N)) UNITY(0:N) = 1.D0 CALL LQCOEF ! Maxwell CALL LQm1CC CALL LQm2CC CALL LQm3CC CALL LQm4CC CALL LQm5CC ! Electron CALL LQe1CC CALL LQe2CC CALL LQe3CC CALL LQe4CC CALL LQe5CC ! Ion CALL LQi1CC CALL LQi2CC CALL LQi3CC CALL LQi4CC CALL LQi5CC ! Beam CALL LQb1CC CALL LQb3CC CALL LQb4CC ! Neutral CALL LQn1CC CALL LQn2CC CALL LQn3CC ! Ripple trapped beam CALL LQr1CC ! Elemental equations -> Nodal equations DO NE = 1, NEMAX NR = NE - 1 DO NQ = 1, NQMAX NC = 0 BLC(NC,NQ,NR) = BLC(NC,NQ,NR) + ELM(NE,1,NC,NQ) * DTt ALC(NC,NQ,NR) = ALC(NC,NQ,NR) + ELM(NE,2,NC,NQ) * DTt CLC(NC,NQ,NR+1) = CLC(NC,NQ,NR+1) + ELM(NE,3,NC,NQ) * DTt BLC(NC,NQ,NR+1) = BLC(NC,NQ,NR+1) + ELM(NE,4,NC,NQ) * DTt DO NC = 1, NLCMAX(NQ) ! iHvsLC : Heaviside function ! iHvsLC = 1 when NLC = 0 ; otherwise iHvsLC = 0. iHvsLC = 1 - ceiling(real(NLC(NC,NQ))/(NQMAX+1)) BLC(NC,NQ,NR) = BLC(NC,NQ,NR) + ELM(NE,1,NC,NQ) * DTf(NQ) ALC(NC,NQ,NR) = ALC(NC,NQ,NR) + ELM(NE,2,NC,NQ) * DTf(NQ) PLC(NC,NQ,NR) = PLC(NC,NQ,NR) + sum(ELM(NE,1:2,NC,NQ)) * DTf(NQ) * iHvsLC CLC(NC,NQ,NR+1) = CLC(NC,NQ,NR+1) + ELM(NE,3,NC,NQ) * DTf(NQ) BLC(NC,NQ,NR+1) = BLC(NC,NQ,NR+1) + ELM(NE,4,NC,NQ) * DTf(NQ) PLC(NC,NQ,NR+1) = PLC(NC,NQ,NR+1) + sum(ELM(NE,3:4,NC,NQ)) * DTf(NQ) * iHvsLC END DO END DO END DO DO NR = 0, NRMAX NLCR(0:NCM,1:NQMAX,NR) = NLC(0:NCM,1:NQMAX) END DO ! Dirichlet condition CALL BOUNDARY(NRMAX,LQm1,0) CALL BOUNDARY(0 ,LQm2,0) CALL BOUNDARY(0 ,LQe2,0) CALL BOUNDARY(NRMAX,LQe2,0) CALL BOUNDARY(0 ,LQe3,0) CALL BOUNDARY(NRMAX,LQe3,0) CALL BOUNDARY(NRMAX,LQe4,0) CALL BOUNDARY(0 ,LQi2,0) CALL BOUNDARY(NRMAX,LQi2,0) CALL BOUNDARY(0 ,LQi3,0) CALL BOUNDARY(NRMAX,LQi3,0) CALL BOUNDARY(NRMAX,LQi4,0) CALL BOUNDARY(NRMAX,LQn2,0) CALL BOUNDARY(NRMAX,LQn3,0) ! When ripple effect is on (FSRP /= 0), ripple diffusion term will be activated. ! Then we must impose two boundary conditions at each equation. CALL BOUNDARY(0 ,LQb3,0) CALL BOUNDARY(NRMAX,LQb3,0) CALL BOUNDARY(NRMAX,LQb4,0) ! IF(FSRP /= 0.D0) CALL BOUNDARY(0,LQr1,0) ! Neumann condition of the continuity equation at the boundary is naturally ! imposed through the diffusion term in the pressure equation. However, this ! condition is very weak because it comes from nondiagonal term. Then by setting ! the density constant in the last element, we can force the density gradient ! to be nought at the boundary explicitly. ! CALL BOUNDARY(NRMAX,LQe1,0,PNeV(NRMAX-1)) ! CALL BOUNDARY(NRMAX,LQi1,0,PNiV(NRMAX-1)) ! Integral term stemming from integration by parts in the diffusion term CALL BOUNDARY(NRMAX,LQm2,1, 2.D0*PSI(NRMAX)*BB/rMU0) CALL BOUNDARY(NRMAX,LQm3,1,-2.D0*R(NRMAX)*Bthb/rMUb2) ! CALL BOUNDARY(NRMAX,LQn1,1, 2.D0*R(NRMAX)*rGASPF) deallocate(ELM) deallocate(rNuIN0, rNuCXN0, rNubeBE, rNubiBI, rNuTeiEI,& & rNuCXN1, rMueNe, rMuiNi, & !dPNeV, dPNiV, & & RUbthV, UethVR, UithVR, EthVR, RUerV, RUirV, UerVR, UirVR, & & FWpheBB, FWpheBB2, FWthpheB, FWpheB, FWphe, dAphV, & & BphBNi, BthBNi, Dbrpft, & & rNuei1EI, rNuei2BthEI, rNuei3EI, & & rNube1BE, rNube2BthBE, rNube3BE, & & Vbparaph, RVbparath, & & Chie1, Chie2, Chii1, Chii2, & & FVpchph, rGASPFA, RNeUer, RNiUir, & & rNueHLphthVR, rNuiHLphthVR)!, & ! & rNue2NCN, rNui2NCN, rat_ei) !!ion deallocate(FWphiBB, FWphiBB2, FWthphiB, FWphiB) deallocate(UNITY) IF(ICALA ==0) ICALA = 1 RETURN END SUBROUTINE TXCALA !*************************************************************** ! ! Coefficients for Equations ! !************************************************************** SUBROUTINE LQCOEF use tx_interface, only : dfdx!, derivs INTEGER(4) :: NR REAL(8) :: BBL !!$ rNuIN0(0:NRMAX) = rNuION(0:NRMAX) * PNeV(0:NRMAX) & !!$ & / (PN01V(0:NRMAX) + PN02V(0:NRMAX) + PN03V(0:NRMAX)) !!$ rNuCXN0(0:NRMAX) = rNuiCX(0:NRMAX) * PNiV(0:NRMAX) & !!$ & / (PN01V(0:NRMAX) + PN02V(0:NRMAX) + PN03V(0:NRMAX)) rNuIN0(0:NRMAX) = FSION * SiVizA(0:NRMAX) * PNeV(0:NRMAX) * 1.D20 rNuCXN0(0:NRMAX) = FSCX * SiVcxA(0:NRMAX) * PNiV(0:NRMAX) * 1.D20 rNuCXN1(0:NRMAX) = rNuiCXT(0:NRMAX) * PNiV(0:NRMAX) * PT01V(0:NRMAX) rNuTeiEI(0:NRMAX) = rNuTei(0:NRMAX) * PNeV(0:NRMAX) / PNiV(0:NRMAX) !!oldNBI rNubeBE(0:NRMAX) = rNube(0:NRMAX) * PNbV(0:NRMAX) / PNeV(0:NRMAX) rNubiBI(0:NRMAX) = rNubi(0:NRMAX) * PNbV(0:NRMAX) / PNiV(0:NRMAX) ! rMueNe(0:NRMAX) = rMue(0:NRMAX) / PNeV(0:NRMAX) ! rMuiNi(0:NRMAX) = rMui(0:NRMAX) / PNiV(0:NRMAX) RUerV(0:NRMAX) = R(0:NRMAX) * UerV(0:NRMAX) RUirV(0:NRMAX) = R(0:NRMAX) * UirV(0:NRMAX) UerVR(1:NRMAX) = UerV(1:NRMAX) / R(1:NRMAX) UerVR(0) = 0.D0 ! Any value is OK. (Never affect the result.) UirVR(1:NRMAX) = UirV(1:NRMAX) / R(1:NRMAX) UirVR(0) = 0.D0 ! Any value is OK. (Never affect the result.) UethVR(1:NRMAX) = UethV(1:NRMAX) / R(1:NRMAX) UethVR(0) = 0.D0 ! Any value is OK. (Never affect the result.) UithVR(1:NRMAX) = UithV(1:NRMAX) / R(1:NRMAX) UithVR(0) = 0.D0 ! Any value is OK. (Never affect the result.) EthVR(1:NRMAX) = EthV(1:NRMAX) / R(1:NRMAX) EthVR(0) = 0.D0 ! Any value is OK. (Never affect the result.) ! CALL DERIVS(PSI,X,LQe1,NQMAX,NRMAX,dPNeV) ! CALL DERIVS(PSI,X,LQi1,NQMAX,NRMAX,dPNiV) dAphV(0:NRMAX) = dfdx(PSI,AphV,NRMAX,0) FWphe(0:NRMAX) =- 2.D0 * dAphV(0:NRMAX) * FWe(0:NRMAX) FWpheBB(0:NRMAX) =- 2.D0 * dAphV(0:NRMAX) * FWthphe(0:NRMAX) FWpheBB2(0:NRMAX) = FWe(0:NRMAX) * BthV(0:NRMAX)**2 !!ion FWphi(0:NRMAX) =- 2.D0 * dAphV(0:NRMAX) * FWi(0:NRMAX) !!ion FWphiBB(0:NRMAX) =- 2.D0 * dAphV(0:NRMAX) * FWthphi(0:NRMAX) !!ion FWphiBB2(0:NRMAX) = FWthi(0:NRMAX) * BthV(0:NRMAX)**2 rNueHLphthVR(1:NRMAX) = rNueHLphth(1:NRMAX) / R(1:NRMAX) rNueHLphthVR(0) = 0.D0 ! Any value is OK. (Never affect the result.) rNuiHLphthVR(1:NRMAX) = rNuiHLphth(1:NRMAX) / R(1:NRMAX) rNuiHLphthVR(0) = 0.D0 ! Any value is OK. (Never affect the result.) DO NR = 0, NRMAX BBL = SQRT(BphV(NR)**2 + BthV(NR)**2) BphBNi(NR) = BphV(NR) / (BBL * PNiV(NR)) IF(NR /= 0) BthBNi(NR) = BthV(NR) / (BBL * R(NR) * PNiV(NR)) if(MDLMOM == 0) then Vbparaph(NR) = Vbpara(NR) RVbparath(NR) = 0.d0 else Vbparaph(NR) = Vbpara(NR) * (BphV(NR) / BBL) RVbparath(NR) = Vbpara(NR) * (BthV(NR) / BBL) * R(NR) end if FWthpheB(NR) = FWthphe(NR) * BBL FWpheB (NR) = FWphe (NR) * BBL !!ion FWthphiB(NR) = FWthphi(NR) * BBL !!ion FWphiB (NR) = FWphi (NR) * BBL END DO BthBNi(0) = 0.D0 ! Any value is OK. (Never affect the result.) !!rp_conv rNubLL(0:NRMAX) = rNubL(0:NRMAX) * rip_rat(0:NRMAX) RUbthV(0:NRMAX) = R(0:NRMAX) * UbthV(0:NRMAX) Dbrpft(0:NRMAX) = Dbrp(0:NRMAX) * ft(0:NRMAX) rNuei1EI (0:NRMAX) = rNuei1 (0:NRMAX) * PNeV(0:NRMAX) / PNiV(0:NRMAX) rNuei2BthEI(0:NRMAX) = rNuei2Bth(0:NRMAX) * PNeV(0:NRMAX) / PNiV(0:NRMAX) rNuei3EI (0:NRMAX) = rNuei3 (0:NRMAX) * PNeV(0:NRMAX) / PNiV(0:NRMAX) rNube1BE (0:NRMAX) = rNube1 (0:NRMAX) * PNbV(0:NRMAX) / PNeV(0:NRMAX) rNube2BthBE(0:NRMAX) = rNube2Bth(0:NRMAX) * PNbV(0:NRMAX) / PNeV(0:NRMAX) rNube3BE (0:NRMAX) = rNube3 (0:NRMAX) * PNbV(0:NRMAX) / PNeV(0:NRMAX) Chie1(0:NRMAX) = Chie(0:NRMAX) + ChiNCTe(0:NRMAX) + ChiNCpe(0:NRMAX) Chie2(0:NRMAX) = Chie(0:NRMAX) + ChiNCTe(0:NRMAX) Chii1(0:NRMAX) = Chii(0:NRMAX) + ChiNCTi(0:NRMAX) + ChiNCpi(0:NRMAX) Chii2(0:NRMAX) = Chii(0:NRMAX) + ChiNCTi(0:NRMAX) FVpchph (0:NRMAX) = FVpch(0:NRMAX) / BphV(0:NRMAX) rGASPFA(0:NRMAX-1) = 0.D0 rGASPFA(NRMAX) = rGASPF RNeUer(0:NRMAX) = R(0:NRMAX) * PNeV(0:NRMAX) * UerV(0:NRMAX) RNiUir(0:NRMAX) = R(0:NRMAX) * PNiV(0:NRMAX) * UirV(0:NRMAX) ! rNue2NCN(0:NRMAX) = rNue2NC(0:NRMAX) / PNeV(0:NRMAX) ! rNui2NCN(0:NRMAX) = rNui2NC(0:NRMAX) / PNiV(0:NRMAX) ! rat_ei(0:NRMAX) = PNeV(0:NRMAX) / PNiV(0:NRMAX) END SUBROUTINE LQCOEF !*************************************************************** ! ! Poisson Equation: phi ! !************************************************************** SUBROUTINE LQm1CC integer(4) :: NEQ = LQm1 ! phi'(0) : 0 ELM(1:NEMAX,1:4,1,NEQ) = 4.D0 * sqeps0 / (AEE * 1.D20) * fem_int(18,UNITY) NLC(1,NEQ) = LQm1 ELM(1:NEMAX,1:4,2,NEQ) = fem_int(1) / sqeps0 NLC(2,NEQ) = LQe1 ELM(1:NEMAX,1:4,3,NEQ) = - PZ * fem_int(1) / sqeps0 NLC(3,NEQ) = LQi1 ELM(1:NEMAX,1:4,4,NEQ) = - PZ * fem_int(1) / sqeps0 * BeamSW NLC(4,NEQ) = LQb1 ELM(1:NEMAX,1:4,5,NEQ) = - PZ * fem_int(2,rip_rat) / sqeps0 * BeamSW * RpplSW NLC(5,NEQ) = LQr1 ! phi(b) : 0 NLCMAX(NEQ) = 5 RETURN END SUBROUTINE LQm1CC !*************************************************************** ! ! Ampere's Law: Atheta' ! !*************************************************************** SUBROUTINE LQm2CC integer(4) :: NEQ = LQm2 ! (r*Atheta)'(0) : 0 ELM(1:NEMAX,1:4,0,NEQ) = fem_int(1) * EPS0 * invDT NLC(0,NEQ) = LQm2 ! rot Bphi ELM(1:NEMAX,1:4,1,NEQ) = - 4.D0 * fem_int(18,UNITY) - 4.D0 * fem_int(4) NLC(1,NEQ) = LQm5 ! Electron current ELM(1:NEMAX,1:4,2,NEQ) = - AEE * 1.D20 * fem_int(1) NLC(2,NEQ) = LQe3 ! Ion current ELM(1:NEMAX,1:4,3,NEQ) = PZ * AEE * 1.D20 * fem_int(1) NLC(3,NEQ) = LQi3 ! Beam ion current ELM(1:NEMAX,1:4,4,NEQ) = PZ * AEE * 1.D20 * fem_int(1) * BeamSW NLC(4,NEQ) = LQb3 ! (1/r)*(r*Atheta)'(NRMAX) : BB NLCMAX(NEQ) = 4 RETURN END SUBROUTINE LQm2CC !*************************************************************** ! ! Ampere's Law: Aphi' ! !*************************************************************** SUBROUTINE LQm3CC integer(4) :: NEQ = LQm3 ! Aphi'(0) : 0 ELM(1:NEMAX,1:4,0,NEQ) = fem_int(1) * EPS0 * rMUb1 * invDT NLC(0,NEQ) = LQm3 ! rot Btheta ELM(1:NEMAX,1:4,1,NEQ) = - 4.D0 * fem_int(18,UNITY) NLC(1,NEQ) = LQm4 ! Electron current ELM(1:NEMAX,1:4,2,NEQ) = - rMUb1 * AEE * 1.D20 * fem_int(1) NLC(2,NEQ) = LQe4 ! Ion current ELM(1:NEMAX,1:4,3,NEQ) = rMUb1 * PZ * AEE * 1.D20 * fem_int(1) NLC(3,NEQ) = LQi4 ! Beam ion current ELM(1:NEMAX,1:4,4,NEQ) = rMUb1 * PZ * AEE * 1.D20 * fem_int(1) * BeamSW NLC(4,NEQ) = LQb4 ! ! Virtual current for helical system ELM(1:NEMAX,1:4,5,NEQ) = rMUb1 * fem_int(-1,AJV) NLC(5,NEQ) = 0 ! Aphi'(NRMAX) : -Bthb NLCMAX(NEQ) = 5 RETURN END SUBROUTINE LQm3CC !************************************************************** ! ! Faraday's Law : Aphi ! !*************************************************************** SUBROUTINE LQm4CC integer(4) :: NEQ = LQm4 ELM(1:NEMAX,1:4,0,NEQ) = fem_int(1) * invDT NLC(0,NEQ) = LQm4 ! Aphi' ELM(1:NEMAX,1:4,1,NEQ) = fem_int(1) / rMUb2 NLC(1,NEQ) = LQm3 NLCMAX(NEQ) = 1 RETURN END SUBROUTINE LQm4CC !*************************************************************** ! ! Faraday's Law : Atheta ! !*************************************************************** SUBROUTINE LQm5CC integer(4) :: NEQ = LQm5 ELM(1:NEMAX,1:4,0,NEQ) = fem_int(1) * invDT NLC(0,NEQ) = LQm5 ! r * Atheta' ELM(1:NEMAX,1:4,1,NEQ) = fem_int(1) / rMU0 NLC(1,NEQ) = LQm2 NLCMAX(NEQ) = 1 RETURN END SUBROUTINE LQm5CC !*************************************************************** ! ! Electron Density Equation ! !*************************************************************** SUBROUTINE LQe1CC integer(4) :: NEQ = LQe1 ELM(1:NEMAX,1:4,0,NEQ) = fem_int(1) * invDT NLC(0,NEQ) = LQe1 ! Divergence ELM(1:NEMAX,1:4,1,NEQ) = - 2.D0 * fem_int(4) NLC(1,NEQ) = LQe2 ! Ionization of n01, n02 and n03 ELM(1:NEMAX,1:4,2,NEQ) = fem_int(2,rNuIN0) NLC(2,NEQ) = LQn1 ELM(1:NEMAX,1:4,3,NEQ) = fem_int(2,rNuIN0) * ThntSW NLC(3,NEQ) = LQn2 ELM(1:NEMAX,1:4,4,NEQ) = fem_int(2,rNuIN0) * BeamSW NLC(4,NEQ) = LQn3 ! Loss to divertor ELM(1:NEMAX,1:4,5,NEQ) = - fem_int(2,rNuL) NLC(5,NEQ) = LQe1 ELM(1:NEMAX,1:4,6,NEQ) = PNeDIV * fem_int(-1,rNuL) NLC(6,NEQ) = 0 ! Generated by NBI (Ionization) ELM(1:NEMAX,1:4,7,NEQ) = (1.D0 - RatCX) * fem_int(-1,SNBe) NLC(7,NEQ) = 0 ! Diffusion of electrons (***AF 2008-06-08) ELM(1:NEMAX,1:4,8,NEQ) = - 4.D0 * fem_int(18,DMAGe) NLC(8,NEQ) = LQe1 NLCMAX(NEQ) = 8 RETURN END SUBROUTINE LQe1CC !*************************************************************** ! ! Electron Radial Flow (SUPG) ! !*************************************************************** SUBROUTINE LQe2CC integer(4) :: NEQ = LQe2 ! Ns*Usr(0) : fixed ELM(1:NEMAX,1:4,0,NEQ) = fem_int(1) * invDT & & + fem_int(8) * fem_int(0) * invDT NLC(0,NEQ) = LQe2 ! Advection ELM(1:NEMAX,1:4,1,NEQ) = - 2.D0 * fem_int( 3,RUerV) + fem_int(2,UerVR) & & +(- 2.D0 * fem_int(10,RUerV) + fem_int(9,UerVR)) * fem_int(0) ELM(1 ,1:4,1,NEQ) = - 2.D0 * fem_int_point( 3,0,RUerV) + fem_int_point(2,0,UerVR) & & +(- 2.D0 * fem_int_point(10,0,RUerV) + fem_int_point(9,0,UerVR)) & & * fem_int_point(0,1) NLC(1,NEQ) = LQe2 ! Centrifugal force ELM(1:NEMAX,1:4,2,NEQ) = fem_int(2,UethVR) & & + fem_int(9,UethVR) * fem_int(0) ELM(1 ,1:4,2,NEQ) = fem_int_point(2,0,UethVR) & & + fem_int_point(9,0,UethVR) * fem_int_point(0,1) NLC(2,NEQ) = LQe3 ! Pressure gradient force ELM(1:NEMAX,1:4,3,NEQ) = - 2.D0 * rKeV / AME * fem_int(17,UNITY) & & - 2.D0 * rKeV / AME * fem_int(18,UNITY) * fem_int(0) NLC(3,NEQ) = LQe5 ! Radial E force ELM(1:NEMAX,1:4,4,NEQ) = 2.D0 * (AEE / AME) * fem_int(17,PNeV) & & + 2.D0 * (AEE / AME) * fem_int(18,PNeV) * fem_int(0) NLC(4,NEQ) = LQm1 ! v x B force ELM(1:NEMAX,1:4,5,NEQ) = - (AEE / AME) * fem_int( 2,BphV) & & - (AEE / AME) * fem_int( 9,BphV) * fem_int(0) NLC(5,NEQ) = LQe3 ELM(1:NEMAX,1:4,6,NEQ) = - 2.D0 * (AEE / AME) * fem_int(16,AphV) & & - 2.D0 * (AEE / AME) * fem_int(20,UNITY,AphV) * fem_int(0) NLC(6,NEQ) = LQe4 NLCMAX(NEQ) = 6 IF(MDFIXT /= 0) THEN ELM(1:NEMAX,1:4,3,NEQ) = - 2.D0 * rKeV / AME & & *( fem_int(17,PTeV) + fem_int(18,PTeV) * fem_int(0) & & + fem_int(16,PTeV) + fem_int(34,PSI,PTeV) * fem_int(0)) NLC(3,NEQ) = LQe1 END IF RETURN END SUBROUTINE LQe2CC !*************************************************************** ! ! Electron Poloidal Flow ! !*************************************************************** SUBROUTINE LQe3CC integer(4) :: NEQ = LQe3 ! Ns*UsTheta(0) : 0 ELM(1:NEMAX,1:4, 0,NEQ) = fem_int(1) * invDT NLC( 0,NEQ) = LQe3 ! Advection ELM(1:NEMAX,1:4, 1,NEQ) = - 2.D0 * fem_int(3,RUerV) NLC( 1,NEQ) = LQe3 ! Viscosity force ! ELM(1:NEMAX,1:4, 2,NEQ) = - 4.D0 * fem_int(18,rMue) & ! & - 4.D0 * fem_int(39,rMueNe,dPNeV) & ! & - 4.D0 * fem_int( 3,rMue) ELM(1:NEMAX,1:4, 2,NEQ) = - 4.D0 * fem_int(18,rMue) & & - 4.D0 * fem_int( 3,rMue) NLC( 2,NEQ) = LQe3 ELM(1:NEMAX,1:4, 3,NEQ) = 4.D0 * fem_int(41,rMue,RUethV) NLC( 3,NEQ) = LQe1 ! Poloidal E force ELM(1:NEMAX,1:4, 4,NEQ) = (AEE / AME) * fem_int(2,PNeV) NLC( 4,NEQ) = LQm2 ! v x B force ELM(1:NEMAX,1:4, 5,NEQ) = (AEE / AME) * fem_int(2,BphV) NLC( 5,NEQ) = LQe2 ! Neoclassical viscosity force ELM(1:NEMAX,1:4, 6,NEQ) = - fem_int(2,rNueNC) NLC( 6,NEQ) = LQe3 ! Collisional friction force with ions ELM(1:NEMAX,1:4, 7,NEQ) = - fem_int(2,rNuei1) NLC( 7,NEQ) = LQe3 ELM(1:NEMAX,1:4, 8,NEQ) = fem_int(2,rNuei1EI) NLC( 8,NEQ) = LQi3 ELM(1:NEMAX,1:4, 9,NEQ) = 2.D0 * fem_int(37,rNuei2Bth,AphV) NLC( 9,NEQ) = LQe4 ELM(1:NEMAX,1:4,10,NEQ) = - 2.D0 * fem_int(37,rNuei2BthEI,AphV) NLC(10,NEQ) = LQi4 ! Collisional friction with beam ions ELM(1:NEMAX,1:4,11,NEQ) = - (AMB / AME) * fem_int(2,rNube1BE) NLC(11,NEQ) = LQe3 ELM(1:NEMAX,1:4,12,NEQ) = (AMB / AME) * fem_int(2,rNube1) * BeamSW NLC(12,NEQ) = LQb3 ELM(1:NEMAX,1:4,13,NEQ) = 2.D0 * (AMB / AME) * fem_int(37,rNube2BthBE,AphV) NLC(13,NEQ) = LQe4 ELM(1:NEMAX,1:4,14,NEQ) = - 2.D0 * (AMB / AME) * fem_int(37,rNube2Bth,AphV) * BeamSW NLC(14,NEQ) = LQb4 ! Turbulent particle transport driver !!$ ELM(1:NEMAX,1:4,15,NEQ) = 2.D0 * AEE / AME * fem_int(38,De,BphV) !!$ NLC(15,NEQ) = LQe1 !!$ !!$ ELM(1:NEMAX,1:4,16,NEQ) = 2.D0 * AEE / AME * fem_int(38,De,BphV/PTeV) * FSVAHLT !!$ NLC(16,NEQ) = LQe5 !!$ !!$ ELM(1:NEMAX,1:4,17,NEQ) = - 2.D0 * AEE / AME * fem_int(38,De,BphV) * FSVAHLT !!$ NLC(17,NEQ) = LQe1 !!$ ELM(1:NEMAX,1:4,15,NEQ) = - 1.D0 / AME * fem_int(2,FWthe) NLC(15,NEQ) = LQe3 ! ELM(1:NEMAX,1:4,16,NEQ) = 1.D0 / AME * fem_int(28,FWthe,rat_ei) ! NLC(16,NEQ) = LQi3 ELM(1:NEMAX,1:4,17,NEQ) = 1.D0 / AME * fem_int(15,FWpheBB) NLC(17,NEQ) = LQe4 ! ELM(1:NEMAX,1:4,18,NEQ) = - 1.D0 / AME * fem_int(44,FWpheBB,rat_ei) ! NLC(18,NEQ) = LQi4 ! Wave interaction force (electron driven) ELM(1:NEMAX,1:4,19,NEQ) = 1.D0 / AME * fem_int(44,FWthpheB,WPM) NLC(19,NEQ) = LQe1 ! Contribution of off-diagonal term due to heat flux ELM(1:NEMAX,1:4,20,NEQ) = - 2.D0 * (rKeV / (AME * AEE)) * fem_int(17,FWthphe) & & * (1.D0 - FSVAHLT) NLC(20,NEQ) = LQe5 ELM(1:NEMAX,1:4,21,NEQ) = + 2.D0 * (rKeV / (AME * AEE)) * fem_int(38,FWthphe,PTeV) & & * (1.D0 - FSVAHLT) NLC(21,NEQ) = LQe1 ELM(1:NEMAX,1:4,22,NEQ) = + 2.D0 * (1.D0 / AME ) * fem_int(37,FWthphe,PhiV) * FSVAHLE NLC(22,NEQ) = LQe1 ! Ad hoc turbulent pinch term ELM(1:NEMAX,1:4,23,NEQ) = - 1.D0 / AME * fem_int(15,FVpch) NLC(23,NEQ) = LQe1 !!ion ! Turbulent particle transport driver (ion driven) !!ion !!ion ELM(1:NEMAX,1:4,15,NEQ) = 1.D0 / AME * fem_int(2,FWthi) !!ion NLC(15,NEQ) = LQi3 !!ion !!ion ELM(1:NEMAX,1:4,16,NEQ) = - 1.D0 / AME * fem_int(2,FWthi) !!ion NLC(16,NEQ) = LQe3 !!ion !!ion ELM(1:NEMAX,1:4,17,NEQ) = - 1.D0 / AME * fem_int(15,FWphiBB) !!ion NLC(17,NEQ) = LQi4 !!ion !!ion ELM(1:NEMAX,1:4,18,NEQ) = 1.D0 / AME * fem_int(15,FWphiBB) !!ion NLC(18,NEQ) = LQe4 !!ion !!ion ELM(1:NEMAX,1:4,19,NEQ) = 2.D0 * (rKeV / (AME * AEE)) * fem_int(37,FWthphi,PTiV) & !!ion * (1.D0 - FSVAHLT) & !!ion & - 2.D0 * (1.D0 / AME ) * fem_int(37,FWthphi,PhiV) * FSVAHLE !!ion NLC(19,NEQ) = LQi1 !!ion !!ion ELM(1:NEMAX,1:4,20,NEQ) = 1.D0 / AME * fem_int(44,FWthphiB,WPM) !!ion NLC(20,NEQ) = LQi1 ! Loss to divertor ELM(1:NEMAX,1:4,24,NEQ) = - 2.D0 * fem_int(2,rNuL) * fact NLC(24,NEQ) = LQe3 ! Collisional friction force with neutrals ELM(1:NEMAX,1:4,25,NEQ) = - fem_int(2,rNu0e) NLC(25,NEQ) = LQe3 ! Helical neoclassical viscosity force (***AF 2008-06-08) ELM(1:NEMAX,1:4,26,NEQ) = - fem_int(2,rNueHLthth) NLC(26,NEQ) = LQe3 !---- 09/11/26 AF ! To deal with the existence of (m=0, n>0) component miki_m 2010/8/26 ! ELM(1:NEMAX,1:4,27,NEQ) = fem_int(15,rNueHLthph) IF(UHphSwitch == 0) THEN ELM(1:NEMAX,1:4,27,NEQ) = fem_int(15,rNueHLthph) else ELM(1:NEMAX,1:4,27,NEQ) = fem_int(22,rNueHLthph) endif NLC(27,NEQ) = LQe4 ! Diffusion of electrons (***AF 2008-06-08) ELM(1:NEMAX,1:4,28,NEQ) = - 4.D0 * fem_int(18,DMAGe) NLC(28,NEQ) = LQe3 ! Poloidal torque due to neoclassical heat flux ELM(1:NEMAX,1:4,29,NEQ) = FSNCPL / 1.D20 * fem_int(-2,R,FQeth) NLC(29,NEQ) = 0 ! ELM(1:NEMAX,1:4,29,NEQ) = FSNCPL * 2.D0 * rKeV / AEE * fem_int(17,rNue2NC) ! NLC(29,NEQ) = LQe5 ! ! ELM(1:NEMAX,1:4,30,NEQ) = - FSNCPL * 2.D0 * rKeV / AEE * fem_int(37,rNue2NCN,PNeV) ! NLC(30,NEQ) = LQe5 ! Ns*UsTheta(NRMAX) : 0 NLCMAX(NEQ) = 29 IF(MDFIXT /= 0) THEN ELM(1:NEMAX,1:4,20,NEQ) = - 2.D0 * (rKeV / (AME * AEE)) * fem_int(37,FWthphe,PTeV) & & * (1.D0 - FSVAHLT) NLC(20,NEQ) = LQe1 ELM(1:NEMAX,1:4,21,NEQ) = 0.D0 NLC(21,NEQ) = 0 ! ELM(1:NEMAX,1:4,29,NEQ) = FSNCPL * 2.D0 * rKeV / AEE * fem_int(37,rNue2NC,PTeV) ! NLC(29,NEQ) = LQe1 ! ! ELM(1:NEMAX,1:4,30,NEQ) = 0.D0 ! NLC(30,NEQ) = 0 END IF RETURN END SUBROUTINE LQe3CC !*************************************************************** ! ! Electron Toroidal Flow ! !*************************************************************** SUBROUTINE LQe4CC integer(4) :: NEQ = LQe4 ! Uephi(0)' : 0 ELM(1:NEMAX,1:4, 0,NEQ) = fem_int(1) * invDT NLC( 0,NEQ) = LQe4 ! Advection ELM(1:NEMAX,1:4, 1,NEQ) = - 2.D0 * fem_int(3,RUerV) NLC( 1,NEQ) = LQe4 ! Viscosity force ! ELM(1:NEMAX,1:4, 2,NEQ) = - 4.D0 *(fem_int(18,rMue) + fem_int(39,rMueNe,dPNeV)) ELM(1:NEMAX,1:4, 2,NEQ) = - 4.D0 * fem_int(18,rMue) NLC( 2,NEQ) = LQe4 ELM(1:NEMAX,1:4, 3,NEQ) = 4.D0 * fem_int(41,rMue,UephV) NLC( 3,NEQ) = LQe1 ! Toroidal E force ELM(1:NEMAX,1:4, 4,NEQ) = (AEE / AME) * fem_int(2,PNeV) NLC( 4,NEQ) = LQm3 ! v x B force ELM(1:NEMAX,1:4, 5,NEQ) = 2.D0 * (AEE / AME) * fem_int(6,AphV) NLC( 5,NEQ) = LQe2 ! Collisional friction with bulk ions ELM(1:NEMAX,1:4, 6,NEQ) = - fem_int(2,rNuei3) NLC( 6,NEQ) = LQe4 ELM(1:NEMAX,1:4, 7,NEQ) = fem_int(2,rNuei3EI) NLC( 7,NEQ) = LQi4 ELM(1:NEMAX,1:4, 8,NEQ) = 2.D0 * fem_int(29,rNuei2Bth,AphV) NLC( 8,NEQ) = LQe3 ELM(1:NEMAX,1:4, 9,NEQ) = - 2.D0 * fem_int(29,rNuei2BthEI,AphV) NLC( 9,NEQ) = LQi3 ! Collisional friction with beam ions ELM(1:NEMAX,1:4,10,NEQ) = - (AMB / AME) * fem_int(2,rNube3BE) NLC(10,NEQ) = LQe4 ELM(1:NEMAX,1:4,11,NEQ) = (AMB / AME) * fem_int(2,rNube3) * BeamSW NLC(11,NEQ) = LQb4 ELM(1:NEMAX,1:4,12,NEQ) = 2.D0 * (AMB / AME) * fem_int(29,rNube2BthBE,AphV) NLC(12,NEQ) = LQe3 ELM(1:NEMAX,1:4,13,NEQ) = - 2.D0 * (AMB / AME) * fem_int(29,rNube2Bth,AphV) * BeamSW NLC(13,NEQ) = LQb3 ! Turbulent particle transport driver !!$ ELM(1:NEMAX,1:4,14,NEQ) = - 2.D0 * AEE / AME * fem_int(38,De,BthV) !!$ NLC(14,NEQ) = LQe1 !!$ !!$ ELM(1:NEMAX,1:4,15,NEQ) = - 2.D0 * AEE / AME * fem_int(38,De,BthV/PTeV) * FSVAHLT !!$ NLC(15,NEQ) = LQe5 !!$ !!$ ELM(1:NEMAX,1:4,16,NEQ) = 2.D0 * AEE / AME * fem_int(38,De,BthV) * FSVAHLT !!$ NLC(16,NEQ) = LQe1 !!$ ELM(1:NEMAX,1:4,14,NEQ) = 1.D0 / AME * fem_int(2,FWpheBB) NLC(14,NEQ) = LQe3 ! ELM(1:NEMAX,1:4,15,NEQ) = - 1.D0 / AME * fem_int(28,FWpheBB,rat_ei) ! NLC(15,NEQ) = LQi3 ELM(1:NEMAX,1:4,16,NEQ) = - 1.D0 / AME * fem_int(2,FWpheBB2) NLC(16,NEQ) = LQe4 ! ELM(1:NEMAX,1:4,17,NEQ) = 1.D0 / AME * fem_int(28,FWpheBB2,rat_ei) ! NLC(17,NEQ) = LQi4 ! Wave interaction force (electron driven) ELM(1:NEMAX,1:4,18,NEQ) = - 1.D0 / AME * fem_int(44,FWpheB,WPM) NLC(18,NEQ) = LQe1 ! Contribution of off-diagonal term due to heat flux ELM(1:NEMAX,1:4,19,NEQ) = 2.D0 * (rKeV / (AME * AEE)) * fem_int(17,FWphe) & & * (1.D0 - FSVAHLT) NLC(19,NEQ) = LQe5 ELM(1:NEMAX,1:4,20,NEQ) = - 2.D0 * (rKeV / (AME * AEE)) * fem_int(38,FWphe,PTeV) & & * (1.D0 - FSVAHLT) NLC(20,NEQ) = LQe1 ELM(1:NEMAX,1:4,21,NEQ) = - 2.D0 * (1.D0 / AME ) * fem_int(37,FWphe,PhiV) * FSVAHLE NLC(21,NEQ) = LQe1 ! Ad hoc turbulent pinch term ELM(1:NEMAX,1:4,22,NEQ) = - 2.D0 / AME * fem_int(36,AphV,FVpchph) NLC(22,NEQ) = LQe1 !!ion ! Turbulent particle transport driver (ion driven) !!ion !!ion ELM(1:NEMAX,1:4,14,NEQ) = - 1.D0 / AME * fem_int(2,FWphiBB) !!ion NLC(14,NEQ) = LQi3 !!ion !!ion ELM(1:NEMAX,1:4,15,NEQ) = 1.D0 / AME * fem_int(2,FWphiBB) !!ion NLC(15,NEQ) = LQe3 !!ion !!ion ELM(1:NEMAX,1:4,16,NEQ) = 1.D0 / AME * fem_int(2,FWphiBB2) !!ion NLC(16,NEQ) = LQi4 !!ion !!ion ELM(1:NEMAX,1:4,17,NEQ) = - 1.D0 / AME * fem_int(2,FWphiBB2) !!ion NLC(17,NEQ) = LQe4 !!ion !!ion ELM(1:NEMAX,1:4,18,NEQ) = - 2.D0 * (rKeV / (AME * AEE)) * fem_int(37,FWphiBB,PTiV) & !!ion & * (1.D0 - FSVAHLT) & !!ion & + 2.D0 * (1.D0 / AME ) * fem_int(37,FWphiBB,PhiV) * FSVAHLE !!ion NLC(18,NEQ) = LQi1 !!ion !!ion ELM(1:NEMAX,1:4,19,NEQ) = - 1.D0 / AME * fem_int(44,FWphiB,WPM) !!ion NLC(19,NEQ) = LQi1 ! Loss to divertor ELM(1:NEMAX,1:4,23,NEQ) = - 2.D0 * fem_int(2,rNuL) * fact NLC(23,NEQ) = LQe4 ! Collisional friction force with neutrals ELM(1:NEMAX,1:4,24,NEQ) = - fem_int(2,rNu0e) NLC(24,NEQ) = LQe4 ! Helical neoclassical viscosity force (***AF 2008-06-08) ! To deal with the existence of (m=0, n>0) component miki_m 2010/8/26 ! ELM(1:NEMAX,1:4,25,NEQ) = fem_int(2,rNueHLphth) IF(UHphSwitch == 0) THEN ELM(1:NEMAX,1:4,25,NEQ) = fem_int(2,rNueHLphth) ELSE ELM(1:NEMAX,1:4,25,NEQ) = fem_int(2,rNueHLphthVR) ELM(1 ,1:4,25,NEQ) = fem_int_point(2,0,rNueHLphthVR) ENDIF NLC(25,NEQ) = LQe3 ! ELM(1:NEMAX,1:4,26,NEQ) = - fem_int(15,rNueHLphph) IF(UHphSwitch == 0) THEN ELM(1:NEMAX,1:4,26,NEQ) = - fem_int(15,rNueHLphph) else ELM(1:NEMAX,1:4,26,NEQ) = - fem_int(2,rNueHLphph) endif NLC(26,NEQ) = LQe4 ! Diffusion of electrons (***AF 2008-06-08) ELM(1:NEMAX,1:4,27,NEQ) = - 4.D0 * fem_int(18,DMAGe) NLC(27,NEQ) = LQe4 NLCMAX(NEQ) = 27 IF(MDFIXT /= 0) THEN ELM(1:NEMAX,1:4,19,NEQ) = 2.D0 * (rKeV / (AME * AEE)) * fem_int(37,FWphe,PTeV) & & * (1.D0 - FSVAHLT) NLC(19,NEQ) = LQe1 ELM(1:NEMAX,1:4,20,NEQ) = 0.D0 NLC(20,NEQ) = 0 END IF RETURN END SUBROUTINE LQe4CC !*************************************************************** ! ! Electron Energy Transport: Te ! !*************************************************************** SUBROUTINE LQe5CC integer(4) :: NEQ = LQe5 ! Temperature evolution IF(MDFIXT == 0) THEN ELM(1:NEMAX,1:4, 0,NEQ) = 1.5D0 * fem_int(1) * invDT NLC( 0,NEQ) = LQe5 ! Advection ELM(1:NEMAX,1:4, 1,NEQ) = - 5.D0 * fem_int(3,RUerV) NLC( 1,NEQ) = LQe5 ! Conduction transport ELM(1:NEMAX,1:4, 2,NEQ) = - 4.D0 * fem_int(18,Chie1) NLC( 2,NEQ) = LQe5 ELM(1:NEMAX,1:4, 3,NEQ) = 4.D0 * fem_int(41,Chie2,PTeV) NLC( 3,NEQ) = LQe1 ! Redundant heat advection ELM(1:NEMAX,1:4, 4,NEQ) = 2.D0 * fem_int(5,RUerV) NLC( 4,NEQ) = LQe5 ! Joule heating ELM(1:NEMAX,1:4, 5,NEQ) = - AEE / rKeV * fem_int(2,EthVR) ELM(1 ,1:4, 5,NEQ) = - AEE / rKeV * fem_int_point(2,0,EthVR) NLC( 5,NEQ) = LQe3 ELM(1:NEMAX,1:4, 6,NEQ) = - AEE / rKeV * fem_int(2,EphV) NLC( 6,NEQ) = LQe4 ! Collisional transfer with ions (Energy equilibration) ELM(1:NEMAX,1:4, 7,NEQ) = - 1.5d0 * fem_int(2,rNuTei) NLC( 7,NEQ) = LQe5 ELM(1:NEMAX,1:4, 8,NEQ) = 1.5d0 * fem_int(2,rNuTeiEI) NLC( 8,NEQ) = LQi5 ! Loss to diverter ELM(1:NEMAX,1:4, 9,NEQ) = - fem_int( 2,rNuL) NLC( 9,NEQ) = LQe5 ELM(1:NEMAX,1:4,10,NEQ) = PNeDIV * fem_int(-2,rNuL,PTeV) NLC(10,NEQ) = 0 ELM(1:NEMAX,1:4,11,NEQ) = - 1.5D0 * fem_int( 2,rNuLTe) NLC(11,NEQ) = LQe5 !!unstable ELM(1:NEMAX,1:4,12,NEQ) = 1.5D0 * PTeDIV * fem_int(-2,rNuLTe,PNeV) !!unstable NLC(12,NEQ) = 0 ELM(1:NEMAX,1:4,12,NEQ) = 1.5D0 * PTeDIV * fem_int( 2,rNuLTe) NLC(12,NEQ) = LQe1 ! Ionization loss of n01, n02 and n03 ELM(1:NEMAX,1:4,13,NEQ) = - (EION * 1.D-3) * fem_int(2,rNuIN0) NLC(13,NEQ) = LQn1 ELM(1:NEMAX,1:4,14,NEQ) = - (EION * 1.D-3) * fem_int(2,rNuIN0) NLC(14,NEQ) = LQn2 ELM(1:NEMAX,1:4,15,NEQ) = - (EION * 1.D-3) * fem_int(2,rNuIN0) * BeamSW NLC(15,NEQ) = LQn3 ! Collisional NBI heating (Perp + Tan) ELM(1:NEMAX,1:4,16,NEQ) = Eb * fem_int(-2,SNB,PNBcol_e) NLC(16,NEQ) = 0 ! Collisional heating with beam !!oldNBI ELM(1:NEMAX,1:4,21,NEQ) = - 0.5D0 * AMB * (PNBCD * Vb) / rKeV & !!oldNBI & * fem_int(2,rNubeBE) !!oldNBI NLC(21,NEQ) = LQe4 !!oldNBI !!oldNBI ELM(1:NEMAX,1:4,22,NEQ) = 0.5D0 * AMB * (PNBCD * Vb) / rKeV & !!oldNBI & * fem_int(2,rNube) * BeamSW !!oldNBI NLC(22,NEQ) = LQb4 ! Simpified Alpha heating ELM(1:NEMAX,1:4,17,NEQ) = 1.D0 / (1.D20 * rKeV) * fem_int(-1,PALFe) NLC(17,NEQ) = 0 ! Direct heating (RF) ELM(1:NEMAX,1:4,18,NEQ) = 1.D0 / (1.D20 * rKeV) * fem_int(-1,PRFe) NLC(18,NEQ) = 0 ! Radiation loss (Bremsstrahlung) ELM(1:NEMAX,1:4,19,NEQ) = - 1.D0 / (1.D20 * rKeV) * fem_int(-1,PBr) NLC(19,NEQ) = 0 ! Diffusion of electrons (***AF 2008-06-08) ELM(1:NEMAX,1:4,20,NEQ) = - 4.D0 * fem_int(18,DMAGe) NLC(20,NEQ) = LQe5 NLCMAX(NEQ) = 20 ELSE ! Fixed temperature profile ELM(1:NEMAX,1:4,0,NEQ) = fem_int(1) * invDT NLC(0,NEQ) = LQe5 NLCMAX(NEQ) = 0 END IF RETURN END SUBROUTINE LQe5CC !*************************************************************** ! ! Ion Density Equation ! !*************************************************************** SUBROUTINE LQi1CC integer(4) :: NEQ = LQi1 ELM(1:NEMAX,1:4,0,NEQ) = fem_int(1) * invDT NLC(0,NEQ) = LQi1 ! Divergence ELM(1:NEMAX,1:4,1,NEQ) = - 2.D0 * fem_int(4) NLC(1,NEQ) = LQi2 ! Ionization of n01, n02 and n03 ELM(1:NEMAX,1:4,2,NEQ) = 1.D0 / PZ * fem_int(2,rNuIN0) NLC(2,NEQ) = LQn1 ELM(1:NEMAX,1:4,3,NEQ) = 1.D0 / PZ * fem_int(2,rNuIN0) * ThntSW NLC(3,NEQ) = LQn2 ELM(1:NEMAX,1:4,4,NEQ) = 1.D0 / PZ * fem_int(2,rNuIN0) * BeamSW NLC(4,NEQ) = LQn3 ! Loss to divertor ELM(1:NEMAX,1:4,5,NEQ) = - 1.D0 / PZ * fem_int(2,rNuL) NLC(5,NEQ) = LQe1 ELM(1:NEMAX,1:4,6,NEQ) = PNeDIV / PZ * fem_int(-1,rNuL) NLC(6,NEQ) = 0 ! Particle source from beam ion ELM(1:NEMAX,1:4,7,NEQ) = fem_int(2,rNuB) NLC(7,NEQ) = LQb1 ! Particle source from ripple trapped beam ions ELM(1:NEMAX,1:4,8,NEQ) = fem_int(28,rNuB,rip_rat) * RpplSW NLC(8,NEQ) = LQr1 ! NBI kick up ions (Charge exchange) ELM(1:NEMAX,1:4,9,NEQ) = - RatCX * fem_int(-1,SNBi) NLC(9,NEQ) = 0 ! Parallel loss reduction due to the potential ! induced by the parallel loss of the beam ions ELM(1:NEMAX,1:4,10,NEQ) = fem_int(2,rNuLB) * BeamSW NLC(10,NEQ) = LQb1 ! Loss cone loss ELM(1:NEMAX,1:4,11,NEQ) = fem_int(-1,SiLC) NLC(11,NEQ) = 0 ! Ion orbit loss ELM(1:NEMAX,1:4,12,NEQ) = - fem_int(2,rNuOL) NLC(12,NEQ) = LQi1 ! Diffusion of ions (***AF 2008-06-08) ELM(1:NEMAX,1:4,13,NEQ) = - 4.D0 * fem_int(18,DMAGi) NLC(13,NEQ) = LQi1 NLCMAX(NEQ) = 13 RETURN END SUBROUTINE LQi1CC !*************************************************************** ! ! Ion Radial Flow (SUPG) ! !*************************************************************** SUBROUTINE LQi2CC integer(4) :: NEQ = LQi2 ! Ns*Usr(0) : fixed ELM(1:NEMAX,1:4,0,NEQ) = fem_int(1) * invDT & & + fem_int(8) * fem_int(0) * invDT NLC(0,NEQ) = LQi2 ! Advection ELM(1:NEMAX,1:4,1,NEQ) = - 2.D0 * fem_int( 3,RUirV) + fem_int(2,UirVR) & & +(- 2.D0 * fem_int(10,RUirV) + fem_int(9,UirVR)) * fem_int(0) ELM(1 ,1:4,1,NEQ) = - 2.D0 * fem_int_point( 3,0,RUirV) + fem_int_point(2,0,UirVR) & & +(- 2.D0 * fem_int_point(10,0,RUirV) + fem_int_point(9,0,UirVR)) & & * fem_int_point(0,1) NLC(1,NEQ) = LQi2 ! Centrifugal force ELM(1:NEMAX,1:4,2,NEQ) = fem_int(2,UithVR) & & + fem_int(9,UithVR) * fem_int(0) ELM(1 ,1:4,2,NEQ) = fem_int_point(2,0,UithVR) & & + fem_int_point(9,0,UithVR) * fem_int_point(0,1) NLC(2,NEQ) = LQi3 ! Pressure gradient force ELM(1:NEMAX,1:4,3,NEQ) = - 2.D0 * rKeV / AMI * fem_int(17,UNITY) & & - 2.D0 * rKeV / AMI * fem_int(18,UNITY) * fem_int(0) NLC(3,NEQ) = LQi5 ! Radial E force ELM(1:NEMAX,1:4,4,NEQ) = - 2.D0 * (PZ * AEE / AMI) * fem_int(17,PNiV) & & - 2.D0 * (PZ * AEE / AMI) * fem_int(18,PNiV) * fem_int(0) NLC(4,NEQ) = LQm1 ! v x B force ELM(1:NEMAX,1:4,5,NEQ) = (PZ * AEE / AMI) * fem_int( 2,BphV) & & + (PZ * AEE / AMI) * fem_int( 9,BphV) * fem_int(0) NLC(5,NEQ) = LQi3 ELM(1:NEMAX,1:4,6,NEQ) = 2.D0 * (PZ * AEE / AMI) * fem_int(16,AphV) & & + 2.D0 * (PZ * AEE / AMI) * fem_int(20,UNITY,AphV) * fem_int(0) NLC(6,NEQ) = LQi4 NLCMAX(NEQ) = 6 IF(MDFIXT /= 0) THEN ELM(1:NEMAX,1:4,3,NEQ) = - 2.D0 * rKeV / AMI & & *( fem_int(17,PTiV) + fem_int(18,PTiV) * fem_int(0) & & + fem_int(16,PTiV) + fem_int(34,PSI,PTiV) * fem_int(0)) NLC(3,NEQ) = LQi1 END IF RETURN END SUBROUTINE LQi2CC !*************************************************************** ! ! Ion Poloidal Flow ! !*************************************************************** SUBROUTINE LQi3CC integer(4) :: NEQ = LQi3 ! Ni*UiTheta(0) : 0 ELM(1:NEMAX,1:4, 0,NEQ) = fem_int(1) * invDT NLC( 0,NEQ) = LQi3 ! Advection ELM(1:NEMAX,1:4, 1,NEQ) = - 2.D0 * fem_int(3,RUirV) NLC( 1,NEQ) = LQi3 ! Viscosity force ! ELM(1:NEMAX,1:4, 2,NEQ) = - 4.D0 * fem_int(18,rMui) & ! & - 4.D0 * fem_int(39,rMuiNi,dPNiV) & ! & - 4.D0 * fem_int( 3,rMui) ELM(1:NEMAX,1:4, 2,NEQ) = - 4.D0 * fem_int(18,rMui) & & - 4.D0 * fem_int( 3,rMui) NLC( 2,NEQ) = LQi3 ELM(1:NEMAX,1:4, 3,NEQ) = 4.D0 * fem_int(41,rMui,RUithV) NLC( 3,NEQ) = LQi1 ! Poloidal E force ELM(1:NEMAX,1:4, 4,NEQ) = - (PZ * AEE / AMI) * fem_int(2,PNiV) NLC( 4,NEQ) = LQm2 ! v x B force ELM(1:NEMAX,1:4, 5,NEQ) = - (PZ * AEE / AMI) * fem_int(2,BphV) NLC( 5,NEQ) = LQi2 ! Neoclassical viscosity force ELM(1:NEMAX,1:4, 6,NEQ) = - fem_int(2,rNuiNC) NLC( 6,NEQ) = LQi3 ! Collisional friction force ELM(1:NEMAX,1:4, 7,NEQ) = - (AME / AMI) * fem_int(2,rNuei1EI) NLC( 7,NEQ) = LQi3 ELM(1:NEMAX,1:4, 8,NEQ) = (AME / AMI) * fem_int(2,rNuei1) NLC( 8,NEQ) = LQe3 ELM(1:NEMAX,1:4, 9,NEQ) = 2.D0 * (AME / AMI) * fem_int(37,rNuei2BthEI,AphV) NLC( 9,NEQ) = LQi4 ELM(1:NEMAX,1:4,10,NEQ) = - 2.D0 * (AME / AMI) * fem_int(37,rNuei2Bth,AphV) NLC(10,NEQ) = LQe4 ! Collisional friction with beam ions ELM(1:NEMAX,1:4,11,NEQ) = - (AMB / AMI) * fem_int(2,rNubiBI) NLC(11,NEQ) = LQi3 ELM(1:NEMAX,1:4,12,NEQ) = (AMB / AMI) * fem_int(2,rNubi) * BeamSW NLC(12,NEQ) = LQb3 ! Turbulent particle transport driver !!$ ELM(1:NEMAX,1:4,13,NEQ) = - 2.D0 * AEE / AMI * fem_int(38,De,BphV) !!$ NLC(13,NEQ) = LQe1 !!$ !!$ ELM(1:NEMAX,1:4,14,NEQ) = - 2.D0 * AEE / AMI * fem_int(38,De,BphV/PTeV) * FSVAHLT !!$ NLC(14,NEQ) = LQe5 !!$ !!$ ELM(1:NEMAX,1:4,15,NEQ) = 2.D0 * AEE / AMI * fem_int(38,De,BphV) * FSVAHLT !!$ NLC(15,NEQ) = LQe1 !!$ ELM(1:NEMAX,1:4,13,NEQ) = 1.D0 / AMI * fem_int(2,FWthe) NLC(13,NEQ) = LQe3 ! ELM(1:NEMAX,1:4,14,NEQ) = - 1.D0 / AMI * fem_int(28,FWthe,rat_ei) ! NLC(14,NEQ) = LQi3 ELM(1:NEMAX,1:4,15,NEQ) = - 1.D0 / AMI * fem_int(15,FWpheBB) NLC(15,NEQ) = LQe4 ! ELM(1:NEMAX,1:4,16,NEQ) = 1.D0 / AMI * fem_int(44,FWpheBB,rat_ei) ! NLC(16,NEQ) = LQi4 ! Wave interaction force (electron driven) ELM(1:NEMAX,1:4,17,NEQ) = - 1.D0 / AMI * fem_int(44,FWthpheB,WPM) NLC(17,NEQ) = LQe1 ! Contribution of off-diagonal term due to heat flux ELM(1:NEMAX,1:4,18,NEQ) = 2.D0 * (rKeV / (AMI * AEE)) * fem_int(17,FWthphe) & & * (1.D0 - FSVAHLT) NLC(18,NEQ) = LQe5 ELM(1:NEMAX,1:4,19,NEQ) = - 2.D0 * (rKeV / (AMI * AEE)) * fem_int(38,FWthphe,PTeV) & & * (1.D0 - FSVAHLT) NLC(19,NEQ) = LQe1 ELM(1:NEMAX,1:4,20,NEQ) = - 2.D0 * (1.D0 / AMI ) * fem_int(37,FWthphe,PhiV) * FSVAHLE NLC(20,NEQ) = LQe1 ! Ad hoc turbulent pinch term ELM(1:NEMAX,1:4,21,NEQ) = 1.D0 / AMI * fem_int(15,FVpch) NLC(21,NEQ) = LQe1 !!ion ! Turbulent particle transport driver (ion driven) !!ion !!ion ELM(1:NEMAX,1:4,13,NEQ) = - 1.D0 / AMI * fem_int(2,FWthi) !!ion NLC(13,NEQ) = LQi3 !!ion !!ion ELM(1:NEMAX,1:4,14,NEQ) = 1.D0 / AMI * fem_int(2,FWthi) !!ion NLC(14,NEQ) = LQe3 !!ion !!ion ELM(1:NEMAX,1:4,15,NEQ) = 1.D0 / AMI * fem_int(15,FWphiBB) !!ion NLC(15,NEQ) = LQi4 !!ion !!ion ELM(1:NEMAX,1:4,16,NEQ) = - 1.D0 / AMI * fem_int(15,FWphiBB) !!ion NLC(16,NEQ) = LQe4 !!ion !!ion ELM(1:NEMAX,1:4,17,NEQ) = - 2.D0 * (rKeV / (AMI * AEE)) * fem_int(37,FWthphi,PTiV) & !!ion & * (1.D0 - FSVAHLT) & !!ion & + 2.D0 * (1.D0 / AMI ) * fem_int(37,FWthphi,PhiV) * FSVAHLE !!ion NLC(17,NEQ) = LQi1 !!ion !!ion ELM(1:NEMAX,1:4,18,NEQ) = - 1.D0 / AMI * fem_int(44,FWthphiB,WPM) !!ion NLC(18,NEQ) = LQi1 ! Loss to divertor ELM(1:NEMAX,1:4,22,NEQ) = - 2.D0 * fem_int(2,rNuL) * fact NLC(22,NEQ) = LQi3 ! Collisional friction force with neutrals ELM(1:NEMAX,1:4,23,NEQ) = - fem_int(2,rNu0i) NLC(23,NEQ) = LQi3 ! Charge exchange force ELM(1:NEMAX,1:4,24,NEQ) = - fem_int(2,rNuiCX) NLC(24,NEQ) = LQi3 ! Loss cone loss ELM(1:NEMAX,1:4,25,NEQ) = fem_int(-1,SiLCth) NLC(25,NEQ) = 0 ! Ion orbit loss ELM(1:NEMAX,1:4,26,NEQ) = - fem_int(2,rNuOL) NLC(26,NEQ) = LQi3 ! Helical Neoclassical viscosity force ELM(1:NEMAX,1:4,27,NEQ) = - fem_int(2,rNuiHLthth) NLC(27,NEQ) = LQi3 ! To deal with the existence of (m=0, n>0) component miki_m 2010/8/26 ! ELM(1:NEMAX,1:4,28,NEQ) = fem_int(15,rNuiHLthph) IF(UHphSwitch == 0) THEN ELM(1:NEMAX,1:4,28,NEQ) = fem_int(15,rNuiHLthph) else ELM(1:NEMAX,1:4,28,NEQ) = fem_int(22,rNuiHLthph) endif NLC(28,NEQ) = LQi4 ! Diffusion of ions (***AF 2008-06-08) ELM(1:NEMAX,1:4,29,NEQ) = - 4.D0 * fem_int(18,DMAGi) NLC(29,NEQ) = LQi3 ! Poloidal torque due to neoclassical heat flux ELM(1:NEMAX,1:4,30,NEQ) = FSNCPL / 1.D20 * fem_int(-2,R,FQith) NLC(30,NEQ) = 0 ! ELM(1:NEMAX,1:4,30,NEQ) = - FSNCPL * 2.D0 * rKeV / (PZ * AEE) * fem_int(17,rNui2NC) ! NLC(30,NEQ) = LQi5 ! ! ELM(1:NEMAX,1:4,31,NEQ) = FSNCPL * 2.D0 * rKeV / (PZ * AEE) * fem_int(37,rNui2NCN,PNiV) ! NLC(31,NEQ) = LQi5 ! Additional torque input ELM(1:NEMAX,1:4,32,NEQ) = 1.D0 / (AMI * 1.D20) * fem_int(-1,Tqp) NLC(32,NEQ) = 0 ! Ns*UsTheta(NRMAX) : 0 NLCMAX(NEQ) = 32 IF(MDFIXT /= 0) THEN ELM(1:NEMAX,1:4,18,NEQ) = 2.D0 * (rKeV / (AMI * AEE)) * fem_int(37,FWthphe,PTeV) & & * (1.D0 - FSVAHLT) NLC(18,NEQ) = LQe1 ELM(1:NEMAX,1:4,19,NEQ) = 0.D0 NLC(19,NEQ) = 0 ! ELM(1:NEMAX,1:4,30,NEQ) = - FSNCPL * 2.D0 * rKeV / (PZ * AEE) * fem_int(37,rNui2NC,PTiV) ! NLC(30,NEQ) = LQi1 ! ! ELM(1:NEMAX,1:4,31,NEQ) = 0.D0 ! NLC(31,NEQ) = 0 END IF RETURN END SUBROUTINE LQi3CC !*************************************************************** ! ! Ion Toroidal Flow ! !*************************************************************** SUBROUTINE LQi4CC integer(4) :: NEQ = LQi4 ! Uiphi'(0) : 0 ELM(1:NEMAX,1:4, 0,NEQ) = fem_int(1) * invDT NLC( 0,NEQ) = LQi4 ! Advection ELM(1:NEMAX,1:4, 1,NEQ) = - 2.D0 * fem_int(3,RUirV) NLC( 1,NEQ) = LQi4 ! Viscosity force ! ELM(1:NEMAX,1:4, 2,NEQ) = - 4.D0 *(fem_int(18,rMui) + fem_int(39,rMuiNi,dPNiV)) ELM(1:NEMAX,1:4, 2,NEQ) = - 4.D0 * fem_int(18,rMui) NLC( 2,NEQ) = LQi4 ELM(1:NEMAX,1:4, 3,NEQ) = 4.D0 * fem_int(41,rMui,UiphV) NLC( 3,NEQ) = LQi1 ! Toroidal E force ELM(1:NEMAX,1:4, 4,NEQ) = - (PZ * AEE / AMI) * fem_int(2,PNiV) NLC( 4,NEQ) = LQm3 ! v x B force ELM(1:NEMAX,1:4, 5,NEQ) = - 2.D0 * (PZ * AEE / AMI) * fem_int(6,AphV) NLC( 5,NEQ) = LQi2 ! Collisional friction with bulk ions ELM(1:NEMAX,1:4, 6,NEQ) = - (AME / AMI) * fem_int(2,rNuei3EI) NLC( 6,NEQ) = LQi4 ELM(1:NEMAX,1:4, 7,NEQ) = (AME / AMI) * fem_int(2,rNuei3) NLC( 7,NEQ) = LQe4 ELM(1:NEMAX,1:4, 8,NEQ) = 2.D0 * (AME / AMI) * fem_int(29,rNuei2BthEI,AphV) NLC( 8,NEQ) = LQi3 ELM(1:NEMAX,1:4, 9,NEQ) = - 2.D0 * (AME / AMI) * fem_int(29,rNuei2Bth,AphV) NLC( 9,NEQ) = LQe3 ! Collisional friction with beam ions ELM(1:NEMAX,1:4,10,NEQ) = - (AMB / AMI) * fem_int(2,rNubiBI) NLC(10,NEQ) = LQi4 ELM(1:NEMAX,1:4,11,NEQ) = (AMB / AMI) * fem_int(2,rNubi) * BeamSW NLC(11,NEQ) = LQb4 ! Turbulent particle transport driver (electron driven) !!$ ELM(1:NEMAX,1:4,12,NEQ) = 2.D0 * AEE / AMI * fem_int(38,De,BthV) !!$ NLC(12,NEQ) = LQe1 !!$ !!$ ELM(1:NEMAX,1:4,13,NEQ) = 2.D0 * AEE / AMI * fem_int(38,De,BthV/PTeV) * FSVAHLT !!$ NLC(13,NEQ) = LQe5 !!$ !!$ ELM(1:NEMAX,1:4,14,NEQ) = - 2.D0 * AEE / AMI * fem_int(38,De,BthV) * FSVAHLT !!$ NLC(14,NEQ) = LQe1 !!$ ELM(1:NEMAX,1:4,12,NEQ) = - 1.D0 / AMI * fem_int(2,FWpheBB) NLC(12,NEQ) = LQe3 ! ELM(1:NEMAX,1:4,13,NEQ) = 1.D0 / AMI * fem_int(28,FWpheBB,rat_ei) ! NLC(13,NEQ) = LQi3 ELM(1:NEMAX,1:4,14,NEQ) = 1.D0 / AMI * fem_int(2,FWpheBB2) NLC(14,NEQ) = LQe4 ! ELM(1:NEMAX,1:4,15,NEQ) = - 1.D0 / AMI * fem_int(28,FWpheBB2,rat_ei) ! NLC(15,NEQ) = LQi4 ! Wave interaction force (electron driven) ELM(1:NEMAX,1:4,16,NEQ) = 1.D0 / AMI * fem_int(44,FWpheB,WPM) NLC(16,NEQ) = LQe1 ! Contribution of off-diagonal term due to heat flux ELM(1:NEMAX,1:4,17,NEQ) = - 2.D0 * (rKeV / (AMI * AEE)) * fem_int(17,FWphe) & & * (1.D0 - FSVAHLT) NLC(17,NEQ) = LQe5 ELM(1:NEMAX,1:4,18,NEQ) = + 2.D0 * (rKeV / (AMI * AEE)) * fem_int(38,FWphe,PTeV) & & * (1.D0 - FSVAHLT) NLC(18,NEQ) = LQe1 ELM(1:NEMAX,1:4,19,NEQ) = + 2.D0 * (1.D0 / AMI ) * fem_int(37,FWphe,PhiV) * FSVAHLE NLC(19,NEQ) = LQe1 ! Ad hoc turbulent pinch term ELM(1:NEMAX,1:4,20,NEQ) = 2.D0 / AMI * fem_int(36,AphV,FVpchph) NLC(20,NEQ) = LQe1 !!ion ! Turbulent particle transport driver (ion driven) !!ion !!ion ELM(1:NEMAX,1:4,12,NEQ) = 1.D0 / AMI * fem_int(2,FWphiBB) !!ion NLC(12,NEQ) = LQi3 !!ion !!ion ELM(1:NEMAX,1:4,13,NEQ) = - 1.D0 / AMI * fem_int(2,FWphiBB) !!ion NLC(13,NEQ) = LQe3 !!ion !!ion ELM(1:NEMAX,1:4,14,NEQ) = - 1.D0 / AMI * fem_int(2,FWphiBB2) !!ion NLC(14,NEQ) = LQi4 !!ion !!ion ELM(1:NEMAX,1:4,15,NEQ) = 1.D0 / AMI * fem_int(2,FWphiBB2) !!ion NLC(15,NEQ) = LQe4 !!ion !!ion ELM(1:NEMAX,1:4,16,NEQ) = 2.D0 * (rKeV / (AMI * AEE)) * fem_int(37,FWphiBB,PTiV) & !!ion & * (1.D0 - FSVAHLT) & !!ion & - 2.D0 * (1.D0 / AMI ) * fem_int(37,FWphiBB,PhiV) * FSVAHLE !!ion NLC(16,NEQ) = LQi1 !!ion !!ion ELM(1:NEMAX,1:4,17,NEQ) = 1.D0 / AMI * fem_int(44,FWphiB,WPM) !!ion NLC(17,NEQ) = LQi1 ! Loss to divertor ELM(1:NEMAX,1:4,21,NEQ) = - 2.D0 * fem_int(2,rNuL) * fact NLC(21,NEQ) = LQi4 ! Collisional friction force with neutrals ELM(1:NEMAX,1:4,22,NEQ) = - fem_int(2,rNu0i) NLC(22,NEQ) = LQi4 ! Charge exchange force ELM(1:NEMAX,1:4,23,NEQ) = - fem_int(2,rNuiCX) NLC(23,NEQ) = LQi4 ! Loss cone loss ELM(1:NEMAX,1:4,24,NEQ) = fem_int(-1,SiLCph) NLC(24,NEQ) = 0 ! Ion orbit loss ELM(1:NEMAX,1:4,25,NEQ) = - fem_int(2,rNuOL) NLC(25,NEQ) = LQi4 ! Helical Neoclassical viscosity force ! To deal with the existence of (m=0, n>0) component miki_m 2010/8/26 ! ELM(1:NEMAX,1:4,26,NEQ) = fem_int(2,rNuiHLphth) IF(UHphSwitch == 0) THEN ELM(1:NEMAX,1:4,26,NEQ) = fem_int(2,rNuiHLphth) ELSE ELM(1:NEMAX,1:4,26,NEQ) = fem_int(2,rNuiHLphthVR) ELM(1 ,1:4,26,NEQ) = fem_int_point(2,0,rNuiHLphthVR) ENDIF NLC(26,NEQ) = LQi3 ! ELM(1:NEMAX,1:4,27,NEQ) = - fem_int(15,rNuiHLphph) IF(UHphSwitch == 0) THEN ELM(1:NEMAX,1:4,27,NEQ) = - fem_int(15,rNuiHLphph) else ELM(1:NEMAX,1:4,27,NEQ) = - fem_int(2,rNuiHLphph) endif NLC(27,NEQ) = LQi4 ! Diffusion of ions (***AF 2008-06-08) ELM(1:NEMAX,1:4,28,NEQ) = - 4.D0 * fem_int(18,DMAGi) NLC(28,NEQ) = LQi4 ! Additional torque input ELM(1:NEMAX,1:4,29,NEQ) = 1.D0 / (AMI * RR * 1.D20) * fem_int(-1,Tqt) NLC(29,NEQ) = 0 NLCMAX(NEQ) = 29 IF(MDFIXT /= 0) THEN ELM(1:NEMAX,1:4,17,NEQ) = - 2.D0 * (rKeV / (AMI * AEE)) * fem_int(37,FWphe,PTeV) & & * (1.D0 - FSVAHLT) NLC(17,NEQ) = LQe1 ELM(1:NEMAX,1:4,18,NEQ) = 0.D0 NLC(18,NEQ) = 0 END IF RETURN END SUBROUTINE LQi4CC !*************************************************************** ! ! Ion Energy Transport: Ti ! !*************************************************************** SUBROUTINE LQi5CC integer(4) :: NEQ = LQi5 ! Temperature evolution IF(MDFIXT == 0) THEN ELM(1:NEMAX,1:4, 0,NEQ) = 1.5D0 * fem_int(1) * invDT NLC( 0,NEQ) = LQi5 ! Advection ELM(1:NEMAX,1:4, 1,NEQ) = - 5.D0 * fem_int(3,RUirV) NLC( 1,NEQ) = LQi5 ! Conduction transport ELM(1:NEMAX,1:4, 2,NEQ) = - 4.D0 * fem_int(18,Chii1) NLC( 2,NEQ) = LQi5 ELM(1:NEMAX,1:4, 3,NEQ) = 4.D0 * fem_int(41,Chii2,PTiV) NLC( 3,NEQ) = LQi1 ! Redundant heat convection term ELM(1:NEMAX,1:4,4,NEQ) = 2.D0 * fem_int(5,RUirV) NLC(4,NEQ) = LQi5 ! Joule heating ELM(1:NEMAX,1:4, 5,NEQ) = PZ * AEE / rKeV * fem_int(2,EthVR) ELM(1 ,1:4, 5,NEQ) = PZ * AEE / rKeV * fem_int_point(2,0,EthVR) NLC( 5,NEQ) = LQi3 ELM(1:NEMAX,1:4, 6,NEQ) = PZ * AEE / rKeV * fem_int(2,EphV) NLC( 6,NEQ) = LQi4 ! Collisional transfer with electrons (Energy equilibration) ELM(1:NEMAX,1:4, 7,NEQ) = - 1.5d0 * fem_int(2,rNuTeiEI) NLC( 7,NEQ) = LQi5 ELM(1:NEMAX,1:4, 8,NEQ) = 1.5d0 * fem_int(2,rNuTei) NLC( 8,NEQ) = LQe5 ! Heating due to beam momentum deposition !!oldNBI ELM(1:NEMAX,1:4, 9,NEQ) = - 0.5D0 * AMB * (PNBCD * Vb) / rKeV & !!oldNBI & * fem_int(2,rNubiBI) !!oldNBI NLC( 9,NEQ) = LQi4 !!oldNBI !!oldNBI ELM(1:NEMAX,1:4,10,NEQ) = 0.5D0 * AMB * (PNBCD * Vb) / rKeV & !!oldNBI & * fem_int(2,rNubi) * BeamSW !!oldNBI NLC(10,NEQ) = LQb4 ELM(1:NEMAX,1:4, 9,NEQ) = AMB * Vb / rKeV * fem_int(28,BthBNi,MNB) ELM(1 ,1:4, 9,NEQ) = AMB * Vb / rKeV * fem_int_point(28,0,BthBNi,MNB) NLC( 9,NEQ) = LQi3 ELM(1:NEMAX,1:4,10,NEQ) = AMB * Vb / rKeV * fem_int(28,BphBNi,MNB) NLC(10,NEQ) = LQi4 ! Loss to diverter ELM(1:NEMAX,1:4,11,NEQ) = - 1.D0 / PZ * fem_int( 2,rNuL) NLC(11,NEQ) = LQi5 ELM(1:NEMAX,1:4,12,NEQ) = PNeDIV / PZ * fem_int(-2,rNuL,PTiV) NLC(12,NEQ) = 0 ELM(1:NEMAX,1:4,13,NEQ) = - 1.5D0 * fem_int( 2,rNuLTi) NLC(13,NEQ) = LQi5 ELM(1:NEMAX,1:4,14,NEQ) = 1.5D0 * PTiDIV * fem_int(-2,rNuLTi,PNiV) NLC(14,NEQ) = 0 ! Ionization heating of n01, n02 and n03 ELM(1:NEMAX,1:4,15,NEQ) = 1.5D0 / PZ * fem_int(28,rNuIN0,PT01V) NLC(15,NEQ) = LQn1 ELM(1:NEMAX,1:4,16,NEQ) = 1.5D0 / PZ * fem_int(28,rNuIN0,PT02V) NLC(16,NEQ) = LQn2 ELM(1:NEMAX,1:4,17,NEQ) = 1.5D0 / PZ * fem_int(28,rNuIN0,PT03V) * BeamSW NLC(17,NEQ) = LQn3 ! Charge exchange loss due to slow neutrals ! (Thermal neutrals are assumed to have same temperature with ions.) ELM(1:NEMAX,1:4,18,NEQ) = - 1.5D0 * fem_int( 2,rNuiCXT) NLC(18,NEQ) = LQi5 ELM(1:NEMAX,1:4,19,NEQ) = 1.5D0 * fem_int(-1,rNuCXN1) NLC(19,NEQ) = 0 ! Collisional NBI heating (Perp + Tan) ELM(1:NEMAX,1:4,20,NEQ) = Eb * fem_int(-2,SNB,PNBcol_i) NLC(20,NEQ) = 0 ! Simplified Alpha heating ELM(1:NEMAX,1:4,21,NEQ) = 1.D0 / (1.D20 * rKeV) * fem_int(-1,PALFi) NLC(21,NEQ) = 0 ! Direct heating (RF) ELM(1:NEMAX,1:4,22,NEQ) = 1.D0 / (1.D20 * rKeV) * fem_int(-1,PRFi) NLC(22,NEQ) = 0 ! Diffusion of ions (***AF 2008-06-08) ELM(1:NEMAX,1:4,23,NEQ) = - 4.D0 * fem_int(18,DMAGi) NLC(23,NEQ) = LQi5 NLCMAX(NEQ) = 23 ELSE ! Fixed temperature profile ELM(1:NEMAX,1:4,0,NEQ) = fem_int(1) * invDT NLC(0,NEQ) = LQi5 NLCMAX(NEQ) = 0 END IF RETURN END SUBROUTINE LQi5CC !*************************************************************** ! ! Beam Ion Density ! !*************************************************************** SUBROUTINE LQb1CC integer(4) :: NEQ = LQb1 ELM(1:NEMAX,1:4,0,NEQ) = fem_int(1) * invDT NLC(0,NEQ) = LQb1 ! NBI particle source (Both charge exchange and ionization) ELM(1:NEMAX,1:4,1,NEQ) = fem_int(-1,SNBb) NLC(1,NEQ) = 0 ! Extracted NBI perpendicular component fallen into the ripple well region ELM(1:NEMAX,1:4,2,NEQ) = - fem_int(-2,SNBPDi,rip_rat) NLC(2,NEQ) = 0 ! Relaxation to thermal ions ELM(1:NEMAX,1:4,3,NEQ) = - fem_int(2,rNuB) NLC(3,NEQ) = LQb1 ! Loss to divertor ELM(1:NEMAX,1:4,4,NEQ) = - fem_int(2,rNuLB) * BeamSW NLC(4,NEQ) = LQb1 ! Ripple trapped beam ions collision with otherwise beam ions ELM(1:NEMAX,1:4,5,NEQ) = - fem_int(28,rNubrp2,rip_rat) * RpplSW NLC(5,NEQ) = LQb1 ELM(1:NEMAX,1:4,6,NEQ) = fem_int(28,rNubrp1,rip_rat) * RpplSW NLC(6,NEQ) = LQr1 !!rp_conv ELM(1:NEMAX,1:4,6,NEQ) = fem_int(-2,PNbrpLV,rNubLL) !!rp_conv NLC(6,NEQ) = 0 ! Ripple diffusion ELM(1:NEMAX,1:4,7,NEQ) = - 4.d0 * fem_int(18,Dbrpft) NLC(7,NEQ) = LQb1 NLCMAX(NEQ) = 7 RETURN END SUBROUTINE LQb1CC !*************************************************************** ! ! Beam Ion Poloidal Flow ! !*************************************************************** SUBROUTINE LQb3CC integer(4) :: NEQ = LQb3 ! Ubth(0) : 0 ELM(1:NEMAX,1:4,0,NEQ) = fem_int(1) * invDT NLC(0,NEQ) = LQb3 ! Poloidal E force ELM(1:NEMAX,1:4,1,NEQ) = - (PZ * AEE / AMB) * fem_int(2,PNbV) NLC(1,NEQ) = LQm2 ! Collisional friction force with electrons ELM(1:NEMAX,1:4,2,NEQ) = - fem_int(2,rNube1) NLC(2,NEQ) = LQb3 ELM(1:NEMAX,1:4,3,NEQ) = fem_int(2,rNube1BE) NLC(3,NEQ) = LQe3 ELM(1:NEMAX,1:4,4,NEQ) = 2.D0 * fem_int(37,rNube2Bth,AphV) * BeamSW NLC(4,NEQ) = LQb4 ELM(1:NEMAX,1:4,5,NEQ) = - 2.D0 * fem_int(37,rNube2BthBE,AphV) NLC(5,NEQ) = LQe4 ! Collisional friction force with ions ELM(1:NEMAX,1:4,6,NEQ) = - fem_int(2,rNubi) NLC(6,NEQ) = LQb3 ELM(1:NEMAX,1:4,7,NEQ) = fem_int(2,rNubiBI) NLC(7,NEQ) = LQi3 ! Collisional friction force with neutrals ELM(1:NEMAX,1:4,8,NEQ) = - fem_int(2,rNu0b) * BeamSW NLC(8,NEQ) = LQb3 ! Charge exchange force ELM(1:NEMAX,1:4,9,NEQ) = - fem_int(2,rNubCX) * BeamSW NLC(9,NEQ) = LQb3 ! NBI momentum source ELM(1:NEMAX,1:4,10,NEQ) = fem_int(-2,RVbparath,MNB) NLC(10,NEQ) = 0 ! Loss to divertor ELM(1:NEMAX,1:4,11,NEQ) = - fem_int(2,rNuLB) * BeamSW * fact NLC(11,NEQ) = LQb3 ! Momentum loss due to collisional ripple trapping ELM(1:NEMAX,1:4,12,NEQ) = - fem_int(28,rNubrp2,rip_rat) * RpplSW NLC(12,NEQ) = LQb3 ! Momentum diffusion arising from beam ion convective flux due to ripple ELM(1:NEMAX,1:4,13,NEQ) = - 4.D0 * fem_int(18,Dbrpft) NLC(13,NEQ) = LQb3 ELM(1:NEMAX,1:4,14,NEQ) = - 4.D0 * fem_int(45,Dbrpft,RUbthV) NLC(14,NEQ) = LQb1 !!neo ! Neoclassical viscosity force !!neo !!neo ELM(1:NEMAX,1:4,15,NEQ) = - fem_int(2,rNuiNC) !!neo NLC(15,NEQ) = LQb3 ! Ubth(NRMAX) : 0 NLCMAX(NEQ) = 14 RETURN END SUBROUTINE LQb3CC !*************************************************************** ! ! Beam Ion Toroidal Flow ! !*************************************************************** SUBROUTINE LQb4CC integer(4) :: NEQ = LQb4 ! - UbPhi(0)' : 0 ELM(1:NEMAX,1:4,0,NEQ) = fem_int(1) * invDT NLC(0,NEQ) = LQb4 ! Toroidal E force ELM(1:NEMAX,1:4,1,NEQ) = - (PZ * AEE / AMB) * fem_int(2,PNbV) NLC(1,NEQ) = LQm3 ! Collisional friction with electrons ELM(1:NEMAX,1:4,2,NEQ) = - fem_int(2,rNube3) NLC(2,NEQ) = LQb4 ELM(1:NEMAX,1:4,3,NEQ) = fem_int(2,rNube3BE) NLC(3,NEQ) = LQe4 ELM(1:NEMAX,1:4,4,NEQ) = 2.D0 * fem_int(29,rNube2Bth,AphV) NLC(4,NEQ) = LQb3 ELM(1:NEMAX,1:4,5,NEQ) = - 2.D0 * fem_int(29,rNube2BthBE,AphV) NLC(5,NEQ) = LQe3 ! Collisional friction with ions ELM(1:NEMAX,1:4,6,NEQ) = - fem_int(2,rNubi) NLC(6,NEQ) = LQb4 ELM(1:NEMAX,1:4,7,NEQ) = fem_int(2,rNubiBI) NLC(7,NEQ) = LQi4 ! Collisional friction force with neutrals ELM(1:NEMAX,1:4,8,NEQ) = - fem_int(2,rNu0b) * BeamSW NLC(8,NEQ) = LQb4 ! Charge exchange force ELM(1:NEMAX,1:4,9,NEQ) = - fem_int(2,rNubCX) * BeamSW NLC(9,NEQ) = LQb4 ! NBI momentum source ELM(1:NEMAX,1:4,10,NEQ) = fem_int(-2,Vbparaph,MNB) NLC(10,NEQ) = 0 ! Loss to divertor ELM(1:NEMAX,1:4,11,NEQ) = - fem_int(2,rNuLB) * BeamSW * fact NLC(11,NEQ) = LQb4 ! Momentum loss due to collisional ripple trapping ELM(1:NEMAX,1:4,12,NEQ) = - fem_int(28,rNubrp2,rip_rat) * RpplSW NLC(12,NEQ) = LQb4 ! Momentum diffusion arising from beam ion convective flux due to ripple ELM(1:NEMAX,1:4,13,NEQ) = - 4.D0 * fem_int(18,Dbrpft) NLC(13,NEQ) = LQb4 ELM(1:NEMAX,1:4,14,NEQ) = - 4.D0 * fem_int(45,Dbrpft,UbphV) NLC(14,NEQ) = LQb1 ! Ubphi(NRMAX) : 0 NLCMAX(NEQ) = 14 RETURN END SUBROUTINE LQb4CC !*************************************************************** ! ! Slow Neutral Transport: n01 ! !*************************************************************** SUBROUTINE LQn1CC integer(4) :: NEQ = LQn1 ELM(1:NEMAX,1:4,0,NEQ) = fem_int(1) * invDT NLC(0,NEQ) = LQn1 ! Diffusion of neutrals ELM(1:NEMAX,1:4,1,NEQ) = - 4.D0 * fem_int(18,D01) NLC(1,NEQ) = LQn1 ! Ionization ELM(1:NEMAX,1:4,2,NEQ) = - 1.D0 / PZ * fem_int(2,rNuIN0) NLC(2,NEQ) = LQn1 ! Generation of thermal neutrals by charge exchange ELM(1:NEMAX,1:4,3,NEQ) = - fem_int(2,rNuCXN0) NLC(3,NEQ) = LQn1 ! Recycling from divertor ELM(1:NEMAX,1:4,4,NEQ) = rGamm0 / PZ * fem_int(2,rNuL) NLC(4,NEQ) = LQe1 ELM(1:NEMAX,1:4,5,NEQ) = - rGamm0 * PNeDIV / PZ * fem_int(-1,rNuL) NLC(5,NEQ) = 0 ! Gas puff ELM(1:NEMAX,1:4,6,NEQ) = 2.D0 * fem_int(-3,R,rGASPFA) NLC(6,NEQ) = 0 NLCMAX(NEQ) = 6 RETURN END SUBROUTINE LQn1CC !*************************************************************** ! ! Thermal Neutral Transport: n02 ! !*************************************************************** SUBROUTINE LQn2CC integer(4) :: NEQ = LQn2 ELM(1:NEMAX,1:4,0,NEQ) = fem_int(1) * invDT NLC(0,NEQ) = LQn2 ! Diffusion of neutrals ELM(1:NEMAX,1:4,1,NEQ) = - 4.D0 * fem_int(18,D02) NLC(1,NEQ) = LQn2 ! Ionization ELM(1:NEMAX,1:4,2,NEQ) = - 1.D0 / PZ * fem_int(2,rNuIN0) NLC(2,NEQ) = LQn2 ! Generation of thermal neutrals by charge exchange ELM(1:NEMAX,1:4,3,NEQ) = fem_int(2,rNuCXN0) NLC(3,NEQ) = LQn1 ! ! NBI particle source (Charge exchange) ! ! ELM(1:NEMAX,1:4,4,NEQ) = RatCX * fem_int(-1,SNBi) ! NLC(4,NEQ) = 0 NLCMAX(NEQ) = 3 RETURN END SUBROUTINE LQn2CC !*************************************************************** ! ! Halo Neutral Transport: n03 ! !*************************************************************** SUBROUTINE LQn3CC integer(4) :: NEQ = LQn3 ELM(1:NEMAX,1:4,0,NEQ) = fem_int(1) * invDT NLC(0,NEQ) = LQn3 ! Diffusion of neutrals ELM(1:NEMAX,1:4,1,NEQ) = - 4.D0 * fem_int(18,D03) NLC(1,NEQ) = LQn3 ! Ionization ELM(1:NEMAX,1:4,2,NEQ) = - 1.D0 / PZ * fem_int(2,rNuIN0) NLC(2,NEQ) = LQn3 ! NBI particle source (Charge exchange) ELM(1:NEMAX,1:4,3,NEQ) = RatCX * fem_int(-1,SNBi) NLC(3,NEQ) = 0 NLCMAX(NEQ) = 3 RETURN END SUBROUTINE LQn3CC !*************************************************************** ! ! Ripple Trapped Beam Ion Density (SUPG) ! !*************************************************************** SUBROUTINE LQr1CC integer(4) :: ne, NEQ = LQr1 real(8) :: RUbrpl, SDbrpl, peclet, coef do ne = 1, nemax RUbrpl = RUbrp(ne-1) + RUbrp(ne) SDbrpl = Dbrp(ne-1)*PSI(ne-1) + Dbrp(ne)*PSI(ne) if (RUbrpl == 0.d0) then coef = 0.d0 elseif (SDbrpl == 0.d0) then coef = 1.d0 / sqrt(15.d0) * hpsi(ne) else peclet = 0.5d0 * RUbrpl * hpsi(ne) / SDbrpl coef = 0.5d0 * langevin(peclet) * hpsi(ne) end if coef=0.d0 !!! no SUPG ELM(NE,1:4,0,NEQ) = fem_int_point(1,NE) * invDT & & + fem_int_point(8,NE) * coef * invDT NLC(0,NEQ) = LQr1 ! NBI perpendicular particle source (Both charge exchange and ionization) ! (Beam ions by perpendicular NBI have few parallel velocity, hence ! they can be easily trapped by ripple wells if they go into the ripple well region.) ! (M.H. Redi, et al., NF 35 (1995) 1191, p.1201 sixth line from the bottom ! at right-hand-side column) ELM(NE,1:4,1,NEQ) =( fem_int_point(-1,NE,SNBPDi) & & + fem_int_point(-8,NE,SNBPDi) * coef) * RpplSW NLC(1,NEQ) = 0 ! Ripple trapped beam ions collision with otherwise beam ions ELM(NE,1:4,2,NEQ) =( fem_int_point(2,NE,rNubrp2) & & + fem_int_point(9,NE,rNubrp2) * coef) * RpplSW NLC(2,NEQ) = LQb1 ELM(NE,1:4,3,NEQ) =(- fem_int_point(2,NE,rNubrp1) & & - fem_int_point(9,NE,rNubrp1) * coef) * RpplSW NLC(3,NEQ) = LQr1 ! Relaxation to thermal ions ELM(NE,1:4,4,NEQ) =(- fem_int_point(2,NE,rNuB) & & - fem_int_point(9,NE,rNuB) * coef) * RpplSW NLC(4,NEQ) = LQr1 ! Ripple loss transport (convective) ELM(NE,1:4,5,NEQ) =(- 2.d0 * fem_int_point( 3,NE,RUbrp) & & - 2.d0 * fem_int_point(10,NE,RUbrp) * coef) * RpplSW NLC(5,NEQ) = LQr1 !!rp_conv ELM(NE,1:4,5,NEQ) = - fem_int_point(2,NE,rNubL) !!rp_conv NLC(5,NEQ) = LQr1 !!rp_conv ELM(NE,1:4,5,NEQ) = - fem_int_point(-2,NE,rNubL,PNbrpV) !!rp_conv NLC(5,NEQ) = 0 ! Ripple loss transport (diffusive) ELM(NE,1:4,6,NEQ) = - 4.d0 * fem_int_point(18,NE,Dbrp) * RpplSW NLC(6,NEQ) = LQr1 end do NLCMAX(NEQ) = 6 RETURN END SUBROUTINE LQr1CC !*************************************************************** ! ! Dirichlet condition ! !*************************************************************** SUBROUTINE BOUNDARY(NR,LQ,ID,VAL) integer(4), intent(in) :: NR, LQ, ID real(8), intent(in), optional :: VAL integer(4) :: NQ, NC, I integer(4), save :: IMAX(1:NQM) type list integer(4) :: IDXNC integer(4) :: IDXNQ end type list type(list), save :: IDX(1:NQM,0:1,1:50) IF(ID == 0) THEN ! Initialize ALC, BLC and CLC at NR ALC(0:NCM,LQ,NR) = 0.D0 BLC(0:NCM,LQ,NR) = 0.D0 CLC(0:NCM,LQ,NR) = 0.D0 PLC(1:NCM,LQ,NR) = 0.D0 IF(ICALA == 0) THEN I = 0 DO NQ = 1, NQMAX DO NC = 1, NLCMAX(NQ) IF(NLCR(NC,NQ,NR) == LQ) THEN I = I + 1 IDX(LQ,MOD(NR,NRMAX),I)%IDXNC = NC IDX(LQ,MOD(NR,NRMAX),I)%IDXNQ = NQ END IF END DO END DO IMAX(LQ) = I ELSE IF(PRESENT(VAL)) THEN IF(NR == 0) THEN DO I = 1, IMAX(LQ) NC = IDX(LQ,MOD(NR,NRMAX),I)%IDXNC NQ = IDX(LQ,MOD(NR,NRMAX),I)%IDXNQ PLC(NC,NQ,NR) = PLC(NC,NQ,NR) + BLC(NC,NQ,NR) * VAL PLC(NC,NQ,NR+1) = PLC(NC,NQ,NR+1) + CLC(NC,NQ,NR+1) * VAL END DO ELSE DO I = 1, IMAX(LQ) NC = IDX(LQ,MOD(NR,NRMAX),I)%IDXNC NQ = IDX(LQ,MOD(NR,NRMAX),I)%IDXNQ PLC(NC,NQ,NR) = PLC(NC,NQ,NR) + BLC(NC,NQ,NR) * VAL PLC(NC,NQ,NR-1) = PLC(NC,NQ,NR-1) + ALC(NC,NQ,NR-1) * VAL END DO END IF END IF IF(NR == 0) THEN DO I = 1, IMAX(LQ) NC = IDX(LQ,MOD(NR,NRMAX),I)%IDXNC NQ = IDX(LQ,MOD(NR,NRMAX),I)%IDXNQ BLC(NC,NQ,NR) = 0.D0 CLC(NC,NQ,NR+1) = 0.D0 END DO ELSE DO I = 1, IMAX(LQ) NC = IDX(LQ,MOD(NR,NRMAX),I)%IDXNC NQ = IDX(LQ,MOD(NR,NRMAX),I)%IDXNQ BLC(NC,NQ,NR) = 0.D0 ALC(NC,NQ,NR-1) = 0.D0 END DO END IF END IF ! Diagonal term on handled variable BLC(1,LQ,NR) = 1.D0 NLCR(1,LQ,NR) = LQ IF(PRESENT(VAL)) PLC(1,LQ,NR) = PLC(1,LQ,NR) - VAL ELSE NLCMAX(LQ) = NLCMAX(LQ) + 1 PLC(NLCMAX(LQ),LQ,NR) = VAL * DTf(LQ) NLCR(NLCMAX(LQ),LQ,NR) = 0 END IF END SUBROUTINE BOUNDARY !*************************************************************** ! ! Approximate Langevin function ! !*************************************************************** real(8) function langevin(x) result(y) real(8), intent(in) :: x if (x < -3.d0) then y = - 1.d0 - 1.d0 / x elseif (x > 3.d0) then y = 1.d0 - 1.d0 / x else y = x / 3.d0 * (1.d0 - abs(x) / 9.d0) end if end function langevin end module tx_coefficients