Qdiv.for

      SUBROUTINE QDIV
      INCLUDE 'Soldiv.fi'
C     COMPUTES PARTICLE AND HEAT FLUXES FROM CORE TO DIVERTOR

c  VP = 19.72*RMAJOR*(AMINOR**2)*ELONG
c  AP = 39.478*RMAJOR*AMINOR*SQRT((1.+ELONG**2)/2.)
      AXPT = DELXPT*6.2832*RMAJOR
      FUSION = PFUSION/(17.6*1.6E-13)
      wtxpt = axpt/ap
      xfrac = wtxpt
C     PARTICLE FLUX CONSISTENT WITH SOURCE TO PLASMA CORE

      XIONIN = YIONSOL*FIN*2.*6.28*RMAJOR
      XIONINXPT = YIONSOLXPT*FIN*2.*6.28*RMAJOR
      XNEUTIN = CURSEP*(1.-ALPHASEP)*(AP-AXPT)
      XNEUTINXPT = CURSEPXPT*(1.-ALPHASEPXPT)*AXPT
      RECYCLE = XNEUTIN + XNEUTINXPT + XIONIN + XIONINXPT
      fluxpartold = fluxpart
      if(jjoptped.eq.9) dn_dt = dln_dt*xnav*vp
      FLUXPART = (RECYCLE + SPELLET + FUELMP-2.*FUSION-dn_dt)/AP
      if(fluxpart.le.0.0.and.itern.lt.10) fluxpart = fluxpartold
C     CALCULATE AVERAGED INWARD NEUTRAL & ION FLUXES
      FLUXNEUTIN = (XNEUTIN + XNEUTINXPT)/AP
      FLUXIONIN  = (XIONIN + XIONINXPT)/AP

C  IONIZATION COOLING DUE TO INCIDENT NEUTRALS
      PIONCORE = (XNEUTIN+XNEUTINXPT)*XK*EIONPLAS*(1.-REABSORBCOR)
C  HEATING OF CORE BY INCIDENT NEUTRALS & IONS
      TXPT = TSOL
      TNSOLM = (TNSOLA - ALPHASEP*TBAR)/(1. - ALPHASEP)
      IF(TNSOLM.LE.0.0) TNSOLM = TNSOLA
      PINSEP = 1.5*XK*(TNSOLM*XNEUTIN + 2.*XIONIN*TSOL +
     2 TNSOLM*XNEUTINXPT + 2.*XIONINXPT*TXPT)
C     CALCULATES CORE HEAT LOSS FRACTIONS
      IF(IOPTOHM.EQ.1) POHMH = POHMIN
      HEATIN =  ((PFUSION/5.+PBEAM)*1.E6+POHMH +
     2            3.*XK*SPELLETREAL*TPEL + PINSEP)
      IF(IPFLUX.EQ.2) PRAD = FRACRAD*HEATIN - PIONCORE
      FRACRAD = (PRAD+PIONCORE)/HEATIN
      IF(FRACRAD.LE.0.95) GOTO 25
      PRAD = (0.95/FRACRAD)*PRAD
      PBREM = (0.95/FRACRAD)*PBREM
      PSYN = (0.95/FRACRAD)*PSYN
      PIMP = (0.95/FRACRAD)*PIMP
      PIONCORE = (0.95/FRACRAD)*PIONCORE
      FRACRAD = (PRAD+PIONCORE)/HEATIN
25    FRACSOL = 1.- FRACRAD
      powsol = fracsol*heatin
      if(jjoptped.eq.9) then
        dw_dt = dlnw_dt*wmhd
        pioncore = 0.0
        prad = pradcorex
      endif
50    FLUXHEAT =  ((PFUSION/5.+PBEAM+POHMH)*1.E6-PRAD - PIONCORE +
     2            3.*XK*SPELLETREAL*TPEL + PINSEP - dw_dt)/AP

C  CONDUCTIVE HEAT FRACTION
      FCOND = 1. - 3.*XK*TSEP*FLUXPART/FLUXHEAT
C  IF(FCOND.LE.0.0) FCOND = 0.1
C     INWARD ION PARTICLE AND HEAT LOSSES FROM SOL INTO CORE
      DNIN = (XIONIN+XIONINXPT)/(2.*6.2832*RMAJOR*BETAG)
      DQNIN = (PINSEP/AP)*XLPERP
C     CALCULATE CONFINEMENT FROM ITER 89-P
75    PHEAT = (PFUSION/5. + PBEAM) + (POHMH-PBREM-PSYN)*1.E-6
C     ***FOR DIII-D****DO NOT SUBTRACT PRAD FROM PHEAT FOR TAUE****
      PHEATTAU = (PFUSION/5. + PBEAM) + (POHMH)*1.E-6
C  *****************POHMH-PBREM-PSYN*********************************
      IF(PHEAT.EQ.0.0) PHEAT = 5.E5

C  CONFINEMENT DEGRADATION
      HN = HRAT*HCONF

      TAU89 =  0.048*(PLASMACUR**0.85)*(RMAJOR**1.2)*(B**0.2)*
     2       ((XNAV/1.E20)**0.1)*SQRT(AMASS*ELONG/PHEATTAU)
      TAU98 = 0.144*(PLASMACUR**0.93)*(RMAJOR**1.39)*(AMINOR**0.58)*
     2((XNAV/1.E20)**0.41)*(B**0.15)*(AMASS**0.19)*
     3(ELONG**0.78)/(PHEATTAU**0.69)
      TAUN89 = 0.05*(PLASMACUR**2)
      TAUN98 = 0.05*(PLASMACUR**2)
      HN89 = HRAT*H89
      HN98 = HRAT*H98
c  h89 & hn89 are the initial input values
      IF(jjoptped.eq.10) THEN
        OMEGTI = OMREWEAKI(14)
        OMEGTE = OMREWEAKE(14)
        IF(OMEGTI.LT.0.0) OMEGTI = 0.0
        IF(OMEGTE.LT.0.0) OMEGTE = 0.0
C     SATURATION
        IF(OMEGTI.GT.3.0e4) OMEGTI = 3.0e4
        IF(OMEGTE.GT.3.0e4) OMEGTE = 3.0e4
C  TURN OFF DEGRADATION
        IF(iopttran.eq.0) THEN
          OMEGTI = 0.0
          OMEGTE = 0.0
          OMEGTIN = 0.0
        END IF

        IF(IOPTTAUE.EQ.98) THEN
          HCONF = H98/(1. + H98*TAU98*C89*1.E-4*(OMEGTI+OMEGTE))
          HN = HN98/(1. + HN98*TAUN98*C89N*1.E-4*OMEGTI)
        ENDIF
        IF(IOPTTAUE.EQ.89) THEN
          HCONF = H89/(1. + H89*TAU89*C89*1.E-4*(OMEGTI+OMEGTE))
          HN = HN89/(1. + HN89*TAUN89*C89N*1.E-4*OMEGTI)
        ENDIF

C  POST-MARFE CALCULATION IF IOPTCONF = 1
        IF (IOPTCONF.EQ.1) THEN
          HCONF = 1.0
          HN = HRAT
        ENDIF

      ENDIF
69    CONTINUE
      taueold = taue
      IF(IOPTTAUE.EQ.89) TAUE = HCONF*TAU89
      IF(IOPTTAUE.EQ.98) TAUE = HCONF*TAU98
      if(itern.gt.10) taue = 0.5*(taueold+taue)
70    TAUP = TAURATIO*TAUE
C  TAUN = TAUP
      TAUN = 0.05*HN*(PLASMACUR**2)
C  taun = taue
      IF(FLUXHEAT.GT.0.) GOTO 100
C     TEMPORARY FIX IF PRAD > PHEAT, TO GET ON WITH ITERATION
      IFIX = IFIX + 1
      IF(IFIX.LT.5) GOTO 80
      GOTO 150
80    PRAD = 0.5*PRAD
      GOTO 50
C       CALCULATES NAV
100   XNAVOLD = XNAV
      XNAV = 0.5*(FLUXPART*TAUN*AP/VP + XNAVOLD)

C     CALCULATES TAV
      TAVOLD = TAV
      TAV = 0.5*(TAVMULT*PHEAT*1.E6*TAUE/(3.*XK*XNAV*VP)
     2          + TAVOLD)
      GOTO 200
150   WRITE (6,'(22A)') 'NEGATIVE HEAT FLUX.'
200   RETURN
      END