ife.f90 Source File


Contents

Source Code


Source Code

! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

module ife_module
  !! Module containing Inertial Fusion Energy device routines
  !! author: P J Knight, CCFE, Culham Science Centre
  !! N/A
  !! This module contains routines for calculating the
  !! parameters of an Inertial Fusion Energy power plant.
  !!
  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#ifndef dp
  use, intrinsic :: iso_fortran_env, only: dp=>real64
#endif
  implicit none

contains

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine ifephy(outfile,iprint)

    !! Routine to calculate the physics parameters of an Inertial Fusion
    !! Energy power plant
    !! author: P J Knight, CCFE, Culham Science Centre
    !! outfile : input integer : output file unit
    !! iprint : input integer : switch for writing to output file (1=yes)
    !! This routine calculates the physics parameters of an Inertial Fusion
    !! Energy power plant.
    !! F/MI/PJK/LOGBOOK12, pp.68,85
    !!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    use constants, only: pi
    use build_variables, only: fwarea
    use error_handling, only: idiags, report_error
    use process_output, only: oheadr, oblnkl, ocmmnt, ovarre
    use ife_variables, only: ifedrv, edrive, gainve, etave, gain, etadrv, tgain, &
      drveff, reprat, pdrive, rrin, pfusife, ifetyp, zl1, r1, zu1, flirad
    use physics_variables, only: fusion_power, wallmw

		use constants, only: pi
    implicit none

    !  Arguments
    integer, intent(in) :: outfile,iprint

    !  Local variables

    real(dp) :: aaion,bmax,dpp,dtheta,emitt,etai,lf,phi,qion, &
         sang,sigma,sigma0,tauf,theta,vi
    integer :: nbeams

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !  Driver calculations

    select case (ifedrv)

    case (-1)
       !  Target gain and driver efficiency dependencies on
       !  driver energy are input

       call driver(edrive,gainve,etave,gain,etadrv)

    case (0)  !  Target gain and driver efficiency are input

       gain = tgain
       etadrv = drveff

    case (1)  !  Laser driver based on SOMBRERO design

       call lasdrv(edrive,gain,etadrv)

    case (2)  !  Heavy-ion beam driver based on OSIRIS design

       aaion = 131.0D0
       bmax = 10.0D0
       dpp = 1.25D-3
       dtheta = 1.3D-4
       emitt = 1.0D-5
       etai = 0.8D0
       lf = 5.0D0
       nbeams = 12
       qion = 1.0D0
       sigma = 8.0D0
       sigma0 = 80.0D0
       tauf = 1.0D-7
       theta = 30.0D-3
       vi = 3.0D6

       call iondrv(aaion,bmax,dpp,dtheta,edrive,emitt,etai,lf, &
            nbeams,qion,sigma,sigma0,tauf,theta,vi,gain,etadrv)

    case(3)
       etadrv = drveff

    case default
       idiags(1) = ifedrv
       call report_error(127)

    end select

    if (ifedrv /= 3) then
        !  Repetition rate (Hz)
        reprat = pdrive / edrive
        !  Fusion power (MW)
        fusion_power = 1.0D-6 * pdrive * gain

    else
        !  Driver Power
        reprat = rrin
        pdrive = reprat * edrive
        !  Gain
        fusion_power = pfusife
        gain = fusion_power / (1.0D-6 * pdrive)
    end if

    !  Wall load (assume total fusion power applies)

    if (ifetyp == 1) then

       !  OSIRIS-type build: First wall subtends a solid angle of 2 pi * SANG

       phi = 0.5D0*pi + atan(zl1/r1)
       sang = 1.0D0 - cos(phi)
       wallmw = fusion_power * 0.5D0*sang / fwarea

    else if (ifetyp == 4) then

       ! 2019 build only has first wall at the top which has a tube at
       ! its centre.  This calculates solid angle and removes tube.

       phi = atan(r1/zu1)
       sang = 1.0D0 - cos(phi)
       phi = atan(flirad/zu1)
       sang = sang - (1.0D0 - cos(phi))
       wallmw = fusion_power * 0.5D0*sang / fwarea

    else
       wallmw = fusion_power / fwarea
    end if

    if (iprint == 0) return

    !  Output section

    call oheadr(outfile,'Physics / Driver Issues')

    select case (ifedrv)

    case (-1, 0)
       call ocmmnt(outfile,'Driver type : generic')
    case (1)
       call ocmmnt(outfile,'Driver type : laser')
    case (2)
       call ocmmnt(outfile,'Driver type : heavy ion beam')
    end select
    call oblnkl(outfile)

    call ovarre(outfile,'Driver energy (J)','(edrive)',edrive)
    call ovarre(outfile,'Driver efficiency','(etadrv)',etadrv)
    call ovarre(outfile,'Driver power reaching target (W)','(pdrive)', &
         pdrive)
    call ovarre(outfile,'Driver repetition rate (Hz)','(reprat)',reprat)
    call ovarre(outfile,'Target gain','(gain)',gain)
    call ovarre(outfile,'Fusion power (MW)','(fusion_power)',fusion_power)
    call ovarre(outfile,'Neutron wall load (MW/m2)','(wallmw)',wallmw)

  end subroutine ifephy

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine iondrv(aaion,bmax,dpp,dtheta,edrive,emitt,etai,lf, &
       nbeams,qion,sigma,sigma0,tauf,theta,vi,gain,etadrv)

    !! Routine to calculate parameters of a heavy ion driver
    !! suitable for inertial fusion energy
    !! author: P J Knight, CCFE, Culham Science Centre
    !! aaion  : input real : Ion mass (amu)
    !! bmax   : input real : Maximum field at the superconductor (T)
    !! dpp    : input real : Beam momentum spread
    !! dtheta : input real : Pointing error (rad)
    !! edrive : input real : Driver energy (J)
    !! emitt  : input real : Normalised beam emittance (metre-rad)
    !! etai   : input real : Axial quadrupole packing fraction at injection
    !! lf     : input real : Distance from final focussing quad to target (m)
    !! nbeams : input integer : Number of beams
    !! qion   : input real : Ion charge state
    !! sigma  : input real : Depressed tune (incl. space charge effects) (deg)
    !! sigma0 : input real : Phase advance per lattice period (tune) (deg)
    !! tauf   : input real : Post-acceleration pulse duration (s)
    !! theta  : input real : Final focussing half angle (rad)
    !! vi     : input real : Injection voltage (V)
    !! gain   : output real : Target gain
    !! etadrv : output real : Driver efficiency
    !! This routine calculates the parameters of a heavy ion driver
    !! suitable for inertial fusion energy.
    !! <P>Currently, the complicated model taken from the reference is not
    !! complete, so it is recommended that the simple model is used
    !! (set <CODE>ISIMP=1</CODE>)
    !! Heavy-ion Driver Design and Scaling, R. Bieri et al.,
    !! Fusion Technology, vol.21 (1992) 1583
    !! Meier and Bieri, Fusion Technology, vol.21 (1992) 1547
    !!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    use constants, only: electron_charge, proton_mass, degrad, pi
    implicit none

    !  Arguments

    real(dp), intent(in) :: aaion,bmax,dpp,dtheta,edrive,emitt,etai,lf,qion, &
         sigma,sigma0,tauf,theta,vi
    integer, intent(in) :: nbeams
    real(dp), intent(out) :: gain,etadrv

    !  Local variables

    real(dp), parameter :: c2 = 8.98755178737D16
    integer, parameter :: isimp = 1  !  Switch for simple model (1=yes)

    real(dp) :: ci,de,dgap,dlcore,drcore,e,eomc2,fins,floss, &
         ibf,ibfo,ibi,ibpc,lfocus,lpf,lpfo,lpi,lppc,lq,phif,phifo, &
         phii,phipc,rion,rs,rs1,rs2,rs3,rs4,sig,sig0,taufo,taui, &
         taupc,tbrad,vf,vfo,vpc,vs,vscore,vshe,vsle,vspc,xhe,xle,xpc
    real(dp), dimension(10) :: eve,gve
    integer :: i,ie,j,nche,ncle,ncpc,nqarrs,nqle,nqpche,nquads

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    if (isimp /= 1) then

       !  Complex IONDRV model (incomplete)...
       !  ====================

       !  Electron charge / (proton mass * c**2)

       eomc2 = electron_charge / (proton_mass*c2)

       !  Degrees to radians

       sig = sigma * degrad
       sig0 = sigma0 * degrad

       !  Beam current coefficient

       ci = cbeam(etai,bmax,aaion,qion,sig,sig0,emitt,vi,eomc2)

       !  Final voltage (V)

       vf = (edrive/(dble(nbeams)*tauf*ci))**(0.666666D0)

       !  Final beam current (A)

       ibf = ci*sqrt(vf)

       !  Ion range (g/cm^2)

       rion = (3.04D-5 * aaion + 349.0D0*aaion**(-2.21D0)) * &
            (qion*vf/1.0D9)**(1.72D0-0.00275D0*aaion)

       !  Convert to kg/m^2

       rion = rion*10.0D0

       !  Total spot size (rms of four terms):
       !  (1) Emittance term

       rs1 = emitt / (betgam(aaion,qion,vf)*theta)

       !  (2) Dispersion (momentum spread) term

       rs2 = 8.0D0*dpp*lf*theta

       !  (3) Alignment and jitter term

       rs3 = dtheta*lf

       !  (4) Space charge term (inaccurate at present)

       rs4 = 1.0D-3

       rs = sqrt(rs1**2 + rs2**2 + rs3**2 + rs4**2)

       !  Target gain (inaccurate at present)

       gain = 62.6D0 - 142.3D0*rs**1.25D0 * rion**0.633D0 + &
            (32.8D0+82.1D0*rs)**0.783D0 * rion**0.647D0 * &
            log(1.0D-6*edrive)

       !  Source/Injector characteristics:
       !  Transport beam radius (m)

       tbrad = ( (emitt**2 * sig0 * c2)/(sig**2 * bmax/1.5D0 * etai * &
            sqrt(vi)) )**(0.333333D0) * sqrt(2.0D0 * aaion * proton_mass / &
            (qion * electron_charge))

       !  Extractor voltage (V)

       vs = 4.0D0/3.0D0 * vi

       !  Source extraction gap width (m)

       dgap = 0.8D0

       !  Injected current per beam (A)

       ibi = 4.0D0 * pi * tbrad**2 * 5.46D-8 * &
            sqrt( qion*vs**3/(aaion*dgap**4) )

       !  Beam loss fraction between injector and target

       floss = 0.08D0

       !  Pulse duration at injection (s)

       taui = 1.0D0/(1.0D0-floss) * tauf * ibf / ibi

       !  Pulse length at injection (m)

       lpi = taui * sqrt( 2.0D0 * qion * electron_charge * vi / (aaion * proton_mass) )

       !  Initial voltage gradient (V/m)

       phii = 0.3D0 * vi / lpi

       !  End of Low Energy Transport Stage:
       !  Insulated fraction of driver length

       fins = 0.85D0

       !  Maximum voltage gradient, limited by insulator flash-over (V/m)

       phifo = 1.0D6 * fins

       !  Voltage at which this limit is reached (V)
       !  The end of the low energy stage is set by the point at which
       !  the limit is reached.

       vfo = (phifo/phii)**(2.0D0/3.0D0) * vi

       !  Beam current (A)

       ibfo = ibi * (vfo/vi)

       !  Pulse duration (s)

       taufo = taui * (vi/vfo)

       !  Pulse length (m)

       lpfo = taui * vi / vfo * sqrt( 2.0D0 * qion * electron_charge * vfo / &
            (aaion * proton_mass) )

       !  Length of the low energy transport stage (m)
       !  (rearrangement of integral of PHI(V).dl)

       xle = 2.0D0*vi**1.5D0 / phii * (1.0D0/sqrt(vi) - 1.0D0/sqrt(vfo))

       !  End of Pulse Compression Stage:
       !  Pulse length = final pulse length (m)

       lpf = tauf * sqrt( 2.0D0 * qion * electron_charge * vf / (aaion * proton_mass) )
       lppc = lpf

       !  Length of the pulse compression region (m)

       xpc = (lpfo - lpf)/0.3D0

       !  Voltage gradient (V/m) - limit already reached

       phipc = phifo

       !  Voltage (V)

       vpc = vfo + (xpc*phipc)

       !  Beam current (A)

       ibpc = ibfo * (lpfo/lppc) * sqrt(vpc/vfo)

       !  Pulse duration (s)

       taupc = taufo * (lppc/lpfo) * sqrt(vfo/vpc)

       !  End of High Energy Transport Stage:
       !  Voltage gradient (V/m) - limit already reached

       phif = phifo

       !  Length of the high energy transport stage (m)

       xhe = (vf - vfo)/phifo

       write(*,*) '                 Injector:'
       write(*,*) '   vi = ',vi
       write(*,*) '  ibi = ',ibi
       write(*,*) ' taui = ',taui
       write(*,*) ' phii = ',phii
       write(*,*) '  lpi = ',lpi
       write(*,*) '                 End of LET stage:'
       write(*,*) '  vfo = ',vfo
       write(*,*) ' ibfo = ',ibfo
       write(*,*) 'taufo = ',taufo
       write(*,*) 'phifo = ',phifo
       write(*,*) ' lpfo = ',lpfo
       write(*,*) '  xle = ',xle
       write(*,*) '                 End of PC stage:'
       write(*,*) '  vpc = ',vpc
       write(*,*) ' ibpc = ',ibpc
       write(*,*) 'taupc = ',taupc
       write(*,*) 'phipc = ',phipc
       write(*,*) ' lppc = ',lppc
       write(*,*) '  xpc = ',xpc
       write(*,*) '                 End of HET stage:'
       write(*,*) '   vf = ',vf
       write(*,*) '  ibf = ',ibf
       write(*,*) ' tauf = ',tauf
       write(*,*) ' phif = ',phif
       write(*,*) '  lpf = ',lpf
       write(*,*) '  xhe = ',xhe

       !  Volt-second requirements and number of cores

       !  Volt-seconds per core (2.4T is field swing of the core
       !  material metglas)

       drcore = 0.8D0
       dlcore = 0.2D0
       vscore = 2.4D0 * drcore * dlcore

       !  Low energy stage

       vsle = taui * vi * log(vfo/vi)
       ncle = 1 + int(vsle/vscore)

       !  Pulse compression stage

       vspc = 2.0D0 * tauf * sqrt(vf) * (sqrt(vpc)-sqrt(vfo))
       ncpc = 1 + int(vspc/vscore)

       !  High energy stage

       drcore = 0.4D0
       dlcore = 0.1D0
       vscore = 2.4D0 * DRCORE * DLCORE

       vshe = 2.0D0 * tauf * vf * (1.0D0 - sqrt(vpc/vf))
       nche = 1 + int(vshe/vscore)

       write(*,*) '                 Number of cores: '
       write(*,*) ' ncle = ',ncle
       write(*,*) ' ncpc = ',ncpc
       write(*,*) ' nche = ',nche

       !  Quadrupole requirements

       !  Effective quadrupole length (m)

       lq = ( emitt * etai * sig0**2 * sqrt(c2*vi) / &
            (sig*(bmax/1.5D0)**2) )**(0.333333D0) * &
            sqrt(2.0D0*aaion*proton_mass/(qion*electron_charge))

       write(*,*) '   lq = ',lq

       !  Actual quadrupole length (m)

       write(*,*) '  alq = ',lq/etai

       !  Length of quadrupole focussing fields: low energy stage

       lfocus = etai * vi/phii * (1.0D0 - vi/vfo)
       nqle = 1 + int(lfocus/lq)

       !  Pulse compression + high energy stages

       lfocus = 2.0D0 * etai * sqrt(vi)/phifo * (sqrt(vf) - sqrt(vfo))

       nqpche = 1 + int(lfocus/lq)

       !  Number of quadrupole arrays

       nqarrs = nqle + nqpche

       !  Number of quadrupoles

       nquads = nqarrs * nbeams

       write(*,*) 'nqarrs = ',nqarrs
       write(*,*) 'nquads = ',nquads

       etadrv = 0.28D0

    else

       !  Simple model
       !  ============

       !  gve(k): target gain at edrive = k MegaJoules

       gve(1)  = 25.0D0
       gve(2)  = 44.0D0
       gve(3)  = 62.0D0
       gve(4)  = 76.0D0
       gve(5)  = 87.0D0
       gve(6)  = 97.0D0
       gve(7)  = 107.0D0
       gve(8)  = 115.0D0
       gve(9)  = 125.0D0
       gve(10) = 132.0D0

       !  eve(k): driver efficiency at edrive = k MegaJoules

       eve(1)  = 0.232D0
       eve(2)  = 0.256D0
       eve(3)  = 0.269D0
       eve(4)  = 0.276D0
       eve(5)  = 0.282D0
       eve(6)  = 0.286D0
       eve(7)  = 0.290D0
       eve(8)  = 0.292D0
       eve(9)  = 0.294D0
       eve(10) = 0.296D0

       e = 1.0D-6 * edrive
       ie = int(e)
       de = e - dble(IE)

       !  Assume linear interpolations and extrapolations
       !  Would be better to prevent extrapolation

       if (ie <= 1) then

          gain   = gve(2) - 1.0D-6*(edrive-2.0D6)*(gve(1)-gve(2))
          etadrv = eve(2) - 1.0D-6*(edrive-2.0d6)*(eve(1)-eve(2))

       else if (ie >= 9) then

          gain   = gve(9) + 1.0D-6*(edrive-9.0D6)*(gve(10)-gve(9))
          etadrv = eve(9) + 1.0D-6*(edrive-9.0D6)*(eve(10)-eve(9))

       else

          gain   = gve(ie) + de*(gve(ie+1)-gve(ie))
          etadrv = eve(ie) + de*(eve(ie+1)-eve(ie))

       end if

       !  Ensure sensible values

       gain = max(0.01D0,gain)
       etadrv = max(0.01D0,etadrv)

    end if

  contains

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    function cbeam(etai,bmax,aaion,qion,sigma,sigma0,emitt,vi,eomc2)

      !! Routine to evaluate the beam current coefficient, CI
      !! author: P J Knight, CCFE, Culham Science Centre
      !! etai   : input real : Axial quadrupole packing fraction at injection
      !! bmax   : input real : Maximum field at the superconductor (T)
      !! aaion  : input real : Ion mass (amu)
      !! qion   : input real : Ion charge state
      !! sigma  : input real : Depressed tune (incl. space charge effects) (rad)
      !! sigma0 : input real : Phase advance per lattice period (tune) (rad)
      !! emitt  : input real : Normalised beam emittance (metre-rad)
      !! vi     : input real : Injection voltage (V)
      !! eomc2  : input real : Electron charge / (proton mass * light speed**2)
      !! This routine calculates the beam current coefficient
      !! of a heavy ion driver suitable for inertial fusion energy.
      !! Heavy-ion Driver Design and Scaling, R. Bieri et al.,
      !! Fusion Technology, vol.21 (1992) 1583
      !!       !
      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      implicit none

      real(dp) :: cbeam

      !  Arguments

      real(dp), intent(in) :: etai,bmax,aaion,qion,sigma,sigma0, &
           emitt,vi,eomc2

      !  Local variables

      real(dp) :: bbe

      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      !  Maximum field at the beam edge

      bbe = bmax/1.5D0

      !  Perform calculation

      cbeam = 2.89D6*(1.0D0-(sigma/sigma0)**2) * &
           ( sigma0**4 * etai**2 * vi * aaion/qion * (emitt/sigma)**2 &
           * bbe**2 )**0.333333D0 &
           * (2.0D0 * qion * eomc2 / aaion)**(5.0D0/6.0D0)

    end function cbeam

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    function betgam(aaion,qion,v)

      !! Routine to calculate the relativistic factor (beta*gamma) for the
      !! heavy ions in the beam driver
      !! author: P J Knight, CCFE, Culham Science Centre
      !! aaion  : input real : Ion mass (amu)
      !! qion   : input real : Ion charge state
      !! v      : input real : Acceleration voltage (V)
      !! This routine calculates the relativistic factor (beta*gamma) for the
      !! heavy ions in a beam driver suitable for inertial fusion energy.
      !! Heavy-ion Driver Design and Scaling, R. Bieri et al.,
      !! Fusion Technology, vol.21 (1992) 1583
      !!       !
      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      implicit none

      real(dp) :: betgam

      !  Arguments

      real(dp), intent(in) :: aaion,qion,v

      !  Local variables

      real(dp), parameter :: c2 = 8.98755178737D16

      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      betgam = sqrt( 2.0D0*qion*electron_charge*v / (aaion*proton_mass*c2) )

    end function betgam

  end subroutine iondrv

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine lasdrv(edrive,gain,etadrv)

    !! Routine to calculate parameters of a laser driver
    !! suitable for inertial fusion energy
    !! author: P J Knight, CCFE, Culham Science Centre
    !! edrive : input real : Driver energy (J)
    !! gain   : output real : Target gain
    !! etadrv : output real : Driver efficiency
    !! This routine calculates the parameters of a laser driver
    !! suitable for inertial fusion energy.
    !! Gain and driver efficiency data are taken from Figures 1 and 2 of
    !! Meier and Rosenberg.
    !! Meier and Rosenberg, Fusion Technology vol.21 (1992) p.1552
    !! F/MI/PJK/LOGBOOK12, p.86
    !!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    implicit none

    !  Arguments

    real(dp), intent(in) :: edrive
    real(dp), intent(out) :: etadrv,gain

    !  Local variables

    real(dp) :: e,de
    real(dp), dimension(10) :: gve,eve
    integer :: ie

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !  GVE(K): target gain at EDRIVE = K MegaJoules

    gve(1)  = 63.0D0
    gve(2)  = 95.0D0
    gve(3)  = 112.0D0
    gve(4)  = 125.0D0
    gve(5)  = 136.0D0
    gve(6)  = 144.0D0
    gve(7)  = 151.0D0
    gve(8)  = 157.0D0
    gve(9)  = 162.0D0
    gve(10) = 166.0D0

    !  EVE(K): driver efficiency at EDRIVE = K MegaJoules

    eve(1)  = 0.082D0
    eve(2)  = 0.079D0
    eve(3)  = 0.076D0
    eve(4)  = 0.072D0
    eve(5)  = 0.069D0
    eve(6)  = 0.064D0
    eve(7)  = 0.059D0
    eve(8)  = 0.054D0
    eve(9)  = 0.048D0
    eve(10) = 0.042D0

    e = 1.0D-6 * edrive
    ie = int(e)
    de = e - dble(ie)

    !  Assume linear interpolations and extrapolations
    !  Would be better to prevent extrapolation

    if (ie <= 1) then

       gain   = gve(2) - 1.0D-6*(edrive-2.0D6)*(gve(1)-gve(2))
       etadrv = eve(2) - 1.0D-6*(edrive-2.0D6)*(eve(1)-eve(2))

    else if (ie >= 9) then

       gain   = gve(9) + 1.0D-6*(edrive-9.0D6)*(gve(10)-gve(9))
       etadrv = eve(9) + 1.0D-6*(edrive-9.0D6)*(eve(10)-eve(9))

    else

       gain = gve(ie)   + de*(gve(ie+1)-gve(ie))
       etadrv = eve(ie) + de*(eve(ie+1)-eve(ie))

    end if

    !  Ensure sensible values

    gain = max(0.01D0,gain)
    etadrv = max(0.01D0,etadrv)

  end subroutine lasdrv

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine driver(edrive,gainve,etave,gain,etadrv)

    !! Routine to calculate parameters of a generic driver
    !! suitable for inertial fusion energy
    !! author: P J Knight, CCFE, Culham Science Centre
    !! edrive : input real : Driver energy (J)
    !! gainve(10) : input real array : Gain vs energy data
    !! etave(10) : input real array : Driver efficiency vs energy data
    !! gain   : output real : Target gain
    !! etadrv : output real : Driver efficiency
    !! This routine calculates the parameters of a generic driver
    !! suitable for inertial fusion energy.
    !! Gain and driver efficiency data are interpolated from input data.
    !! F/MI/PJK/LOGBOOK12, p.85
    !!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    implicit none

    !  Arguments

    real(dp), intent(in) :: edrive
    real(dp), dimension(10), intent(in) :: etave,gainve
    real(dp), intent(out) :: etadrv,gain

    !  Local variables

    real(dp) :: de,e
    integer :: ie

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !  The arrays contain data points for EDRIVE = 1MJ, 2MJ, ... , 10MJ

    e = 1.0D-6 * edrive
    ie = int(e)
    de = e - dble(ie)

    !  Assume linear interpolations and extrapolations

    if (ie <= 1) then

       gain =  gainve(2) - 1.0D-6*(edrive-2.0D6)*(gainve(1)-gainve(2))
       etadrv = etave(2) - 1.0D-6*(edrive-2.0D6)*(etave(1)-etave(2))

    else if (ie >= 9) then

       gain =  gainve(9) + 1.0D-6*(edrive-9.0D6)*(gainve(10)-gainve(9))
       etadrv = etave(9) + 1.0D-6*(edrive-9.0D6)*(etave(10)-etave(9))

    else

       gain =  gainve(ie) + de*(gainve(ie+1)-gainve(ie))
       etadrv = etave(ie) + de*(etave(ie+1)-etave(ie))

    end if

    !  Ensure sensible values

    gain = max(0.01D0,gain)
    etadrv = max(0.01D0,etadrv)

  end subroutine driver

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine ifebld(outfile,iprint)
    !! Routine to create the build of an inertial fusion energy device
    !! and to calculate the material volumes for the device core
    !! author: P J Knight, CCFE, Culham Science Centre
    !! outfile : input integer : output file unit
    !! iprint : input integer : switch for writing to output file (1=yes)
    !! This routine constructs the build of an inertial fusion energy device
    !! and calculates the material volumes for the device core.
    !! F/MI/PJK/LOGBOOK12, p.52
    !!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    use process_output, only: oheadr, obuild, ovarre
    use ife_variables, only: ifetyp, chrad, r1, fwdr, r2, v1dr, r3, bldr, r4, &
      v2dr, r5, shdr, r6, v3dr, r7, zl7, v3dzl, zl6, shdzl, zl5, v2dzl, zl4, &
      bldzl, zl3, v1dzl, zl2, fwdzl, zl1, chdzl, chdzu, zu1, fwdzu, zu2, &
      v1dzu, zu3, bldzu, zu4, v2dzu, zu5, shdzu, zu6, v3dzu, zu7, fwmatv, &
      v1matv, blmatv, v2matv, shmatv, v3matv, chmatv
		use constants, only: mfile
    implicit none

    !  Arguments
    integer, intent(in) :: outfile,iprint

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    select case (ifetyp)
    case (1)  !  OSIRIS-type device
       call osibld
    case (2)  !  SOMBRERO-type device
       call sombld
    case (3)  !  HYLIFE-II-type device
       call hylbld
    case (4)  !  2019 device
       call bld2019
    case default  !  Generic device
       call genbld
    end select

    if (iprint == 0) return

    !  Output section

    call oheadr(outfile,'Radial Build')
    write(outfile,20)
    20  format(T43,'Thickness (m)',T60,'Radius (m)')

    call obuild(outfile,'Device centreline',0.0D0,0.0D0)
    call obuild(outfile,'Chamber',chrad,r1)
    call ovarre(mfile,'Chamber (m)','(chrad)',chrad)
    call obuild(outfile,'First Wall',fwdr,r2)
    call ovarre(mfile,'First Wall (m)','(fwdr)',fwdr)
    call obuild(outfile,'Void 1',v1dr,r3)
    call ovarre(mfile,'Void 1 (m)','(v1dr)',v1dr)
    call obuild(outfile,'Blanket',bldr,r4)
    call ovarre(mfile,'Blanket (m)','(bldr)',bldr)
    call obuild(outfile,'Void 2',v2dr,r5)
    call ovarre(mfile,'Void 2 (m)','(v2dr)',v2dr)
    call obuild(outfile,'Shield',shdr,r6)
    call ovarre(mfile,'Shield (m)','(shdr)',shdr)
    call obuild(outfile,'Void 3',v3dr,r7)
    call ovarre(mfile,'Void 3 (m)','(v3dr)',v3dr)

    30  format(T43,'Thickness (m)',T60,'Height (m)')

    if (ifetyp /= 4) then

        call oheadr(outfile,'Vertical Build')
        write(outfile,30)

        call obuild(outfile,'Base of device',0.0D0,-zl7)
        call obuild(outfile,'Void 3',v3dzl,-zl6)
        call obuild(outfile,'Shield',shdzl,-zl5)
        call obuild(outfile,'Void 2',v2dzl,-zl4)
        call obuild(outfile,'Blanket',bldzl,-zl3)
        call obuild(outfile,'Void 1',v1dzl,-zl2)
        call obuild(outfile,'First Wall',fwdzl,-zl1)
        call obuild(outfile,'Chamber',chdzl,0.0D0)
        call obuild(outfile,'Chamber',chdzu,zu1)
        call obuild(outfile,'First Wall',fwdzu,zu2)
        call obuild(outfile,'Void 1',v1dzu,zu3)
        call obuild(outfile,'Blanket',bldzu,zu4)
        call obuild(outfile,'Void 2',v2dzu,zu5)
        call obuild(outfile,'Shield',shdzu,zu6)
        call obuild(outfile,'Void 3',v3dzu,zu7)

    else

        call oheadr(outfile,'Vertical Build - Midplane')
        write(outfile,30)

        call obuild(outfile,'Base of device',0.0D0,-zl7)
        call obuild(outfile,'Void 3',v3dzl,-zl6)
        call obuild(outfile,'Shield',shdzl,-zl5)
        call obuild(outfile,'Void 2',v2dzl,-zl4)
        call obuild(outfile,'Blanket',bldzl,-zl3)
        call obuild(outfile,'Void 1',v1dzl,-zl2)
        call obuild(outfile,'First Wall',fwdzl,-zl1)
        call obuild(outfile,'Chamber',chdzl,0.0D0)
        call obuild(outfile,'Chamber',chdzu,zu1)
        call obuild(outfile,'First Wall',fwdzu,zu2)
        call obuild(outfile,'Void 1',v1dzu,zu3)
        call obuild(outfile,'Blanket',bldzu-bldzu,zu4-bldzu)
        call obuild(outfile,'Void 2',v2dzu,zu5-bldzu)
        call obuild(outfile,'Shield',shdzu,zu6-bldzu)
        call obuild(outfile,'Void 3',v3dzu+bldzu,zu7)

        call oheadr(outfile,'Vertical Build - Edge')
        write(outfile,30)

        call obuild(outfile,'Base of device',0.0D0,-zl7)
        call obuild(outfile,'Void 3',v3dzl,-zl6)
        call ovarre(mfile,'Void 3 lower (m)','(v3dzl)',v3dzl)
        call obuild(outfile,'Shield',shdzl,-zl5)
        call ovarre(mfile,'Shield lower (m)','(shdzl)',shdzl)
        call obuild(outfile,'Void 2',v2dzl,-zl4)
        call ovarre(mfile,'Void 2 lower (m)','(v2dzl)',v2dzl)
        call obuild(outfile,'Blanket',bldzl,-zl3)
        call ovarre(mfile,'Blanket lower (m)','(bldzl)',bldzl)
        call obuild(outfile,'Void 1',v1dzl,-zl2)
        call ovarre(mfile,'Void 1 lower (m)','(v1dzl)',v1dzl)
        call obuild(outfile,'First Wall',fwdzl,-zl1)
        call ovarre(mfile,'First Wall lower (m)','(fwdzl)',fwdzl)
        call obuild(outfile,'Chamber',chdzl,0.0D0)
        call ovarre(mfile,'Chamber lower (m)','(chdzl)',chdzl)
        call obuild(outfile,'Chamber',chdzu,zu1)
        call ovarre(mfile,'Chamber upper (m)','(chdzu)',chdzu)
        call obuild(outfile,'First Wall',fwdzu,zu2)
        call ovarre(mfile,'First Wall upper (m)','(fwdzu)',fwdzu)
        call obuild(outfile,'Void 1',v1dzu,zu3)
        call ovarre(mfile,'Void 1 upper (m)','(v1dzu)',v1dzu)
        call obuild(outfile,'Blanket',bldzu,zu4)
        call ovarre(mfile,'Blanket upper (m)','(bldzu)',bldzu)
        call obuild(outfile,'Void 2',v2dzu,zu5)
        call ovarre(mfile,'Void 2 upper (m)','(v2dzu)',v2dzu)
        call obuild(outfile,'Shield',shdzu,zu6)
        call ovarre(mfile,'Shield upper (m)','(shdzu)',shdzu)
        call obuild(outfile,'Void 3',v3dzu,zu7)
        call ovarre(mfile,'Void 3 upper (m)','(v3dzu)',v3dzu)

    end if

    !  Print matrix of material volumes

    call oheadr(outfile,'Material volumes')

    write(outfile,*) '         Chamber  1st wall  Void 1  Blanket  ' &
         //' Void 2   Shield   Void 3'
    write(outfile,'(A9,7(1pe9.2))') 'void     ', &
         chmatv(0), &
         (fwmatv(1,0)+fwmatv(2,0)+fwmatv(3,0)), &
         (v1matv(1,0)+v1matv(2,0)+v1matv(3,0)), &
         (blmatv(1,0)+blmatv(2,0)+blmatv(3,0)), &
         (v2matv(1,0)+v2matv(2,0)+v2matv(3,0)), &
         (shmatv(1,0)+shmatv(2,0)+shmatv(3,0)), &
         (v3matv(1,0)+v3matv(2,0)+v3matv(3,0))
    write(outfile,'(A9,7(1pe9.2))') 'steel    ', &
         chmatv(1), &
         (fwmatv(1,1)+fwmatv(2,1)+fwmatv(3,1)), &
         (v1matv(1,1)+v1matv(2,1)+v1matv(3,1)), &
         (blmatv(1,1)+blmatv(2,1)+blmatv(3,1)), &
         (v2matv(1,1)+v2matv(2,1)+v2matv(3,1)), &
         (shmatv(1,1)+shmatv(2,1)+shmatv(3,1)), &
         (v3matv(1,1)+v3matv(2,1)+v3matv(3,1))

    write(outfile,'(A9,7(1pe9.2))') 'carbon   ', &
         chmatv(2), &
         (fwmatv(1,2)+fwmatv(2,2)+fwmatv(3,2)), &
         (v1matv(1,2)+v1matv(2,2)+v1matv(3,2)), &
         (blmatv(1,2)+blmatv(2,2)+blmatv(3,2)), &
         (v2matv(1,2)+v2matv(2,2)+v2matv(3,2)), &
         (shmatv(1,2)+shmatv(2,2)+shmatv(3,2)), &
         (v3matv(1,2)+v3matv(2,2)+v3matv(3,2))

    write(outfile,'(A9,7(1pe9.2))') 'FLiBe    ', &
         chmatv(3), &
         (fwmatv(1,3)+fwmatv(2,3)+fwmatv(3,3)), &
         (v1matv(1,3)+v1matv(2,3)+v1matv(3,3)), &
         (blmatv(1,3)+blmatv(2,3)+blmatv(3,3)), &
         (v2matv(1,3)+v2matv(2,3)+v2matv(3,3)), &
         (shmatv(1,3)+shmatv(2,3)+shmatv(3,3)), &
         (v3matv(1,3)+v3matv(2,3)+v3matv(3,3))

    write(outfile,'(A9,7(1pe9.2))') 'Li2O     ', &
         chmatv(4), &
         (fwmatv(1,4)+fwmatv(2,4)+fwmatv(3,4)), &
         (v1matv(1,4)+v1matv(2,4)+v1matv(3,4)), &
         (blmatv(1,4)+blmatv(2,4)+blmatv(3,4)), &
         (v2matv(1,4)+v2matv(2,4)+v2matv(3,4)), &
         (shmatv(1,4)+shmatv(2,4)+shmatv(3,4)), &
         (v3matv(1,4)+v3matv(2,4)+v3matv(3,4))

    write(outfile,'(A9,7(1pe9.2))') 'concrete ', &
         chmatv(5), &
         (fwmatv(1,5)+fwmatv(2,5)+fwmatv(3,5)), &
         (v1matv(1,5)+v1matv(2,5)+v1matv(3,5)), &
         (blmatv(1,5)+blmatv(2,5)+blmatv(3,5)), &
         (v2matv(1,5)+v2matv(2,5)+v2matv(3,5)), &
         (shmatv(1,5)+shmatv(2,5)+shmatv(3,5)), &
         (v3matv(1,5)+v3matv(2,5)+v3matv(3,5))

    write(outfile,'(A9,7(1pe9.2))') 'helium   ', &
         chmatv(6), &
         (fwmatv(1,6)+fwmatv(2,6)+fwmatv(3,6)), &
         (v1matv(1,6)+v1matv(2,6)+v1matv(3,6)), &
         (blmatv(1,6)+blmatv(2,6)+blmatv(3,6)), &
         (v2matv(1,6)+v2matv(2,6)+v2matv(3,6)), &
         (shmatv(1,6)+shmatv(2,6)+shmatv(3,6)), &
         (v3matv(1,6)+v3matv(2,6)+v3matv(3,6))

    write(outfile,'(A9,7(1pe9.2))') 'xenon    ', &
         chmatv(7), &
         (fwmatv(1,7)+fwmatv(2,7)+fwmatv(3,7)), &
         (v1matv(1,7)+v1matv(2,7)+v1matv(3,7)), &
         (blmatv(1,7)+blmatv(2,7)+blmatv(3,7)), &
         (v2matv(1,7)+v2matv(2,7)+v2matv(3,7)), &
         (shmatv(1,7)+shmatv(2,7)+shmatv(3,7)), &
         (v3matv(1,7)+v3matv(2,7)+v3matv(3,7))

    write(outfile,'(A9,7(1pe9.2))') 'lithium  ', &
         chmatv(8), &
         (fwmatv(1,8)+fwmatv(2,8)+fwmatv(3,8)), &
         (v1matv(1,8)+v1matv(2,8)+v1matv(3,8)), &
         (blmatv(1,8)+blmatv(2,8)+blmatv(3,8)), &
         (v2matv(1,8)+v2matv(2,8)+v2matv(3,8)), &
         (shmatv(1,8)+shmatv(2,8)+shmatv(3,8)), &
         (v3matv(1,8)+v3matv(2,8)+v3matv(3,8))

  contains

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine bld2019
     !! Routine to create the build of a 2019 inertial fusion energy
     !! device, and to calculate the material volumes for the device core
     !! author: S I Muldrew, CCFE, Culham Science Centre
     !! None
     !! This routine constructs the build of a modern inertial fusion energy
     !! device, assumed to be cylindrically-symmetric, with a pool at bottom
     !! and top corners and with a lower shield at the centre.  See diagram
     !! attached to Issue #907.
     !! Issue #907
     !!      !
     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    use constants, only: pi
    use build_variables, only: fwarea
    use error_handling, only: report_error
    use buildings_variables, only: trcl, stcl
    use fwbs_variables, only: tbr, emult
    use ife_variables, only: blmatf, bldrc, lipmw, etali, taufall, rrmax, chvol, &
      fwvol, v1vol, blvol, v2vol, shvol, flirad, v3vol, maxmat, chmatf, fwmatf, &
      v1matf, v2matf, shmatf, v3matf
    implicit none

     !  Arguments

     !  Local variables

     integer :: i,j
     real(dp), save :: g, vel, acurt, mdot, phi, sang, li_frac

     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     ! Check input
     if ((fwdr.gt.1.0D-9).or.(v1dr.gt.1.0d-9)) then
          call report_error(230)
     else if ((fwdzu.gt.1.0D-9).or.(v1dzu.gt.1.0D-9).or.(v2dzu.gt.1.0D-9)) then
          call report_error(231)
     else if ((fwdzl.gt.1.0D-9).or.(v1dzl.gt.1.0d-9).or.(v2dzu.gt.1.0D-9)) then
          call report_error(232)
     end if

     ! Lithium Pump

     ! Acceleration due to gravity (m/s^2)
     g = 9.81D0

     ! Velocity
     vel = sqrt(2.0D0*g*(chdzu+bldzu))

     ! Lithium Fraction
     blmatf(1,8) = 0.91*sqrt(bldzu/(chdzu+bldzu))
     blmatf(1,0) = 1.0D0 - blmatf(1,8)

     ! Spatial Thickness
     bldr = bldrc / blmatf(1,8)

     ! Area
     acurt = pi * ((chrad+bldr)**2.0D0-chrad**2.0D0)

     ! Mass Flow
     mdot = 512.0D0 * vel * blmatf(1,8) * acurt

     ! Pump Power (MW)
     lipmw = 1.0D-6 * mdot * g * (chdzl+chdzu+bldzu+bldzl) / etali

     ! Fall Time
     taufall = 2.0D0 * (chdzl+chdzu+bldzu) / vel

     rrmax = 1.0D0 / taufall

     ! TBR and Emult model was for spherical lithium
     ! Remove reactor head
     phi = atan(r1/zu1)
     sang = 1.0D0 - cos(phi)
     li_frac = 1.0D0 - 0.5D0 * sang

     ! TBR
     tbr = 3.7418D0 * (1.0D0 / (1.0D0 + exp(-2.6366D0 * bldrc)) - 0.5D0)
     tbr = tbr * li_frac

     ! Energy Multiplication
     emult = 2.2414D0 * (1.0D0 / (1.0D0 + exp(-3.0038D0 * bldrc)) - 0.5D0)
     emult = emult * li_frac

     !  Radial build

     r1 = chrad
     r2 = r1 + fwdr
     r3 = r2 + v1dr
     r4 = r3 + bldr
     r5 = r4 + v2dr
     r6 = r5 + shdr
     r7 = r6 + v3dr

     !  Vertical build (below midplane)

     zl1 = chdzl
     zl2 = zl1 + fwdzl
     zl3 = zl2 + v1dzl
     zl4 = zl3 + bldzl
     zl5 = zl4 + v2dzl
     zl6 = zl5 + shdzl
     zl7 = zl6 + v3dzl

     !  Vertical build (above midplane)

     zu1 = chdzu
     zu2 = zu1 + fwdzu
     zu3 = zu2 + v1dzu
     zu4 = zu3 + bldzu
     zu5 = zu4 + v2dzu
     zu6 = zu5 + shdzu

     v3dzu = (zu6+zl6) + trcl + stcl + 5.1D0 + 9.41D-6 * 1.0D5

     zu7 = zu6 + v3dzu

     !  Component volumes
     !  The following notation applies below:
     !  J=1 : side part
     !  J=2 : top part
     !  J=3 : bottom part

     !  Chamber

     chvol = pi * r1*r1 * (zu1 + zl1)

     !  First wall

     fwvol(1) = pi * (r2*r2 - r1*r1) * (zu1 + zl1)
     fwvol(2) = pi * r2*r2 * (zu2 - zu1)
     fwvol(3) = pi * r2*r2 * (zl2 - zl1)

     !  First void

     v1vol(1) = pi * (r3*r3 - r2*r2) * (zu2 + zl2)
     v1vol(2) = pi * r3*r3 * (zu3 - zu2)
     v1vol(3) = pi * r3*r3 * (zl3 - zl2)

     !  Blanket
     !  Radial Blanket - between void 2 and chamber
     blvol(1) = pi * (r4*r4 - r3*r3) * (zu3 + zl3)
     !  Upper Blanket - Pool radially between shield and
     !  chamber of input height.
     blvol(2) = pi * (r5*r5 - r3*r3) * bldzu
     !  Lower Blanket - Pool filling base of device
     blvol(3) = pi * r5*r5 * (zl4 - zl3)

     !  Second void

    v2vol(1) = pi * (r5*r5 - r4*r4) * (chdzl+chdzu)
    v2vol(2) = 0.0D0
    v2vol(3) = 0.0D0

    !  Shield
    shvol(1) = pi * (r6*r6 - r5*r5) * (zu5 + zl5)
    ! Top Section is in three parts to account for the dip at
    ! the centre.  The first is the horizontal top, the second is the
    ! horizontal
    shvol(2) = pi * (((r6*r6 - (chrad-shdr)*(chrad-shdr)) * shdzu) &
      + ((r1*r1 - flirad*flirad) * shdzu) &
      + ( (r1*r1 - (r1-shdzu)*(r1-shdzu)) *(bldzu-shdzu) ))
    shvol(3) = pi * r6*r6 * (zl6 - zl5)

    !  Third void

    v3vol(1) = pi * (r7*r7 - r6*r6) * (zu6 + zl6)
    v3vol(2) = pi * r7*r7 * (zu7 - zu6) + pi * ((r1-shdzu)*(r1-shdzu) - flirad*flirad) * bldzu
    v3vol(3) = pi * r7*r7 * (zl7 - zl6)

     !  Material volumes

     do i = 0,maxmat
        chmatv(i) = max(0.0D0, chvol * chmatf(i))
        do j = 1,3
           fwmatv(j,i) = max(0.0D0, fwvol(j) * fwmatf(j,i))
           v1matv(j,i) = max(0.0D0, v1vol(j) * v1matf(j,i))
           blmatv(j,i) = max(0.0D0, blvol(j) * blmatf(j,i))
           v2matv(j,i) = max(0.0D0, v2vol(j) * v2matf(j,i))
           shmatv(j,i) = max(0.0D0, shvol(j) * shmatf(j,i))
           v3matv(j,i) = max(0.0D0, v3vol(j) * v3matf(j,i))
        end do
     end do

     !  First wall area
     !  The chamber is surrounded by liquid on three sides
     !  with only the top being solid.  This is considered part
     !  of the shield. There is a target injector tube at the
     !  centre of this area.
     fwarea = pi * (r1*r1 - flirad*flirad)

   end subroutine bld2019

   ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    subroutine genbld

      !! Routine to create the build of a generic inertial fusion energy
      !! device, and to calculate the material volumes for the device core
      !! author: P J Knight, CCFE, Culham Science Centre
      !! None
      !! This routine constructs the build of a generic inertial fusion energy
      !! device, assumed to be cylindrically-symmetric, and to calculate
      !! the material volumes for the device core.
      !! F/MI/PJK/LOGBOOK12, p.52
      !!       !
      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      use constants, only: pi
      use build_variables, only: fwarea
      use ife_variables, only: chvol, fwvol, v1vol, blvol, v2vol, shvol, v3vol, &
        maxmat, chmatf, fwmatf, v1matf, blmatf, v2matf, shmatf, v3matf

      implicit none

      !  Local variables
      integer :: i,j

      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      !  Radial build

      r1 = chrad
      r2 = r1 + fwdr
      r3 = r2 + v1dr
      r4 = r3 + bldr
      r5 = r4 + v2dr
      r6 = r5 + shdr
      r7 = r6 + v3dr

      !  Vertical build (below midplane)

      zl1 = chdzl
      zl2 = zl1 + fwdzl
      zl3 = zl2 + v1dzl
      zl4 = zl3 + bldzl
      zl5 = zl4 + v2dzl
      zl6 = zl5 + shdzl
      zl7 = zl6 + v3dzl

      !  Vertical build (above midplane)

      zu1 = chdzu
      zu2 = zu1 + fwdzu
      zu3 = zu2 + v1dzu
      zu4 = zu3 + bldzu
      zu5 = zu4 + v2dzu
      zu6 = zu5 + shdzu
      zu7 = zu6 + v3dzu

      !  Component volumes
      !  The following notation applies below:
      !  J=1 : side part
      !  J=2 : top part
      !  J=3 : bottom part

      !  Chamber

      chvol = pi * r1*r1 * (zu1 + zl1)

      !  First wall

      fwvol(1) = pi * (r2*r2 - r1*r1) * (zu1 + zl1)
      fwvol(2) = pi * r2*r2 * (zu2 - zu1)
      fwvol(3) = pi * r2*r2 * (zl2 - zl1)

      !  First void

      v1vol(1) = pi * (r3*r3 - r2*r2) * (zu2 + zl2)
      v1vol(2) = pi * r3*r3 * (zu3 - zu2)
      v1vol(3) = pi * r3*r3 * (zl3 - zl2)

      !  Blanket

      blvol(1) = pi * (r4*r4 - r3*r3) * (zu3 + zl3)
      blvol(2) = pi * r4*r4 * (zu4 - zu3)
      blvol(3) = pi * r4*r4 * (zl4 - zl3)

      !  Second void

      v2vol(1) = pi * (r5*r5 - r4*r4) * (zu4 + zl4)
      v2vol(2) = pi * r5*r5 * (zu5 - zu4)
      v2vol(3) = pi * r5*r5 * (zl5 - zl4)

      !  Shield

      shvol(1) = pi * (r6*r6 - r5*r5) * (zu5 + zl5)
      shvol(2) = pi * r6*r6 * (zu6 - zu5)
      shvol(3) = pi * r6*r6 * (zl6 - zl5)

      !  Third void

      v3vol(1) = pi * (r7*r7 - r6*r6) * (zu6 + zl6)
      v3vol(2) = pi * r7*r7 * (zu7 - zu6)
      v3vol(3) = pi * r7*r7 * (zl7 - zl6)

      !  Material volumes

      do i = 0,maxmat
         chmatv(i) = max(0.0D0, chvol * chmatf(i))
         do j = 1,3
            fwmatv(j,i) = max(0.0D0, fwvol(j) * fwmatf(j,i))
            v1matv(j,i) = max(0.0D0, v1vol(j) * v1matf(j,i))
            blmatv(j,i) = max(0.0D0, blvol(j) * blmatf(j,i))
            v2matv(j,i) = max(0.0D0, v2vol(j) * v2matf(j,i))
            shmatv(j,i) = max(0.0D0, shvol(j) * shmatf(j,i))
            v3matv(j,i) = max(0.0D0, v3vol(j) * v3matf(j,i))
         end do
      end do

      !  First wall area

      fwarea = 2.0D0*pi*r1 * ((zu1 + zl1) + r1)

    end subroutine genbld

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    subroutine osibld
      !! Routine to create the build of an inertial fusion energy
      !! device, based on the design of the OSIRIS study,
      !! and to calculate the material volumes for the device core
      !! author: P J Knight, CCFE, Culham Science Centre
      !! None
      !! This routine constructs the build of an inertial fusion energy
      !! device, based on the design of the OSIRIS study, and to calculate
      !! the material volumes for the device core.
      !! F/MI/PJK/LOGBOOK12, p.56
      !!       !
      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      use constants, only: pi
      use build_variables, only: fwarea

      implicit none

      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      !  Careful choice of thicknesses, and assuming that the FLiBe
      !  inlet radius is small, allows the generic build calculation
      !  to be roughly applicable.
      call genbld

      !  First wall area: no true first wall at bottom of chamber
      fwarea = 2.0D0*pi*r1*(zu1 + zl1) + pi*r1*r1

    end subroutine osibld

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    subroutine sombld
      !! Routine to create the build of an inertial fusion energy
      !! device, based on the design of the SOMBRERO study,
      !! and to calculate the material volumes for the device core
      !! author: P J Knight, CCFE, Culham Science Centre
      !! None
      !! This routine constructs the build of an inertial fusion energy
      !! device, based on the design of the SOMBRERO study, and to calculate
      !! the material volumes for the device core.
      !! Sviatoslavsky et al, Fusion Technology vol.21 (1992) 1470
      !! F/MI/PJK/LOGBOOK12, p.53
      !!       !
      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      use constants, only: pi
      use build_variables, only: fwarea
      use ife_variables, only: chvol, fwvol, v1vol, blvol, v2vol, somtdr, sombdr, &
        shvol, v3vol, maxmat, chmatf, fwmatf, blmatf, v1matf, v2matf, shmatf, v3matf

      implicit none

      !  Local variables

      real(dp), parameter :: third  = 1.0D0/3.0D0
      real(dp) :: chcylh,ddz,dvol
      integer :: i,j

      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      !  Radial build
      r1 = chrad
      r2 = r1 + fwdr
      r3 = r2 + v1dr
      r4 = r3 + bldr
      r5 = r4 + v2dr
      r6 = r5 + shdr
      r7 = r6 + v3dr

      !  Vertical build (below midplane)
      zl1 = chdzl
      zl2 = zl1 + fwdzl
      zl3 = zl2 + v1dzl
      zl4 = zl3 + bldzl
      zl5 = zl4 + v2dzl
      zl6 = zl5 + shdzl
      zl7 = zl6 + v3dzl

      !  Vertical build (above midplane)
      zu1 = chdzu
      zu2 = zu1 + fwdzu
      zu3 = zu2 + v1dzu
      zu4 = zu3 + bldzu
      zu5 = zu4 + v2dzu
      zu6 = zu5 + shdzu
      zu7 = zu6 + v3dzu

      !  The SOMBRERO chamber is made up of a cylindrical first wall/
      !  blanket, with conical regions above and below. Outside this is
      !  a cylindrical shield.

      !  Component volumes
      !  The following notation applies below:
      !  J=1 : side part
      !  J=2 : top part
      !  J=3 : bottom part

      !  Chamber : CHCYLH is the height of the cylindrical part
      chcylh = chdzu + chdzl - 2.0D0*chrad

      chvol = pi * r1*r1 * (chcylh + 2.0D0*third * chrad)

      !  First wall
      fwvol(1) = pi * (r2*r2 - r1*r1) * chcylh
      fwvol(2) = third * pi * ( r2*r2*(chrad+fwdzu) - r1*r1*chrad )
      fwvol(3) = third * pi * ( r2*r2*(chrad+fwdzl) - r1*r1*chrad )

      !  First void
      v1vol(1) = pi * (r3*r3 - r2*r2) * chcylh
      v1vol(2) = third * pi * ( r3*r3*(chrad+fwdzu+v1dzu) -  &
           r2*r2*(chrad+fwdzu) )
      v1vol(3) = third * pi * ( r3*r3*(chrad+fwdzl+v1dzl) -  &
           r2*r2*(chrad+fwdzl) )

      !  Blanket:  SOMTDR and SOMBDR are the radii of the cylindrical
      !  sections at the top/bottom of the blanket
      !  DDZ = Height of top cylindrical section (by similar triangles)
      !  DVOL = Volume of top cylindrical section, less the internal cone

      blvol(1) = pi * (r4*r4 - r3*r3) * chcylh

      blvol(2) = third * pi * ( r4*r4*(chrad+fwdzu+v1dzu+bldzu) -  &
           r3*r3*(chrad+fwdzu+v1dzu) )
      ddz = (chrad+fwdzu+v1dzu+bldzu)/(chrad+fwdr+v1dr+bldr)*somtdr
      dvol = 2.0D0*third * pi * somtdr*somtdr * ddz

      blvol(2) = blvol(2) + dvol

      !  Ditto for bottom region...

      blvol(3) = third * pi * ( r4*r4*(chrad+fwdzl+v1dzl+bldzl) -  &
           r3*r3*(chrad+fwdzl+v1dzl) )
      ddz = (chrad+fwdzl+v1dzl+bldzl)/(chrad+fwdr+v1dr+bldr)*sombdr
      dvol = 2.0D0*third * pi * sombdr*sombdr * ddz

      blvol(3) = blvol(3) + dvol

      !  Second void
      v2vol(1) = pi * (r5*r5 - r4*r4) * chcylh
      v2vol(2) = pi * r5*r5 * (zu5 - chdzu + chrad) - ( &
           fwvol(2) + v1vol(2) + blvol(2) + (third*pi*r1*r1*chrad) )
      v2vol(3) = pi * r5*r5 * (zl5 - chdzl + chrad) - ( &
           fwvol(3) + v1vol(3) + blvol(3) + (third*pi*r1*r1*chrad) )

      !  Shield
      shvol(1) = pi * (r6*r6 - r5*r5) * (zu6 + zl6)
      shvol(2) = pi * r5*r5 * (zu6 - zu5)
      shvol(3) = pi * r5*r5 * (zl6 - zl5)

      !  Third void
      v3vol(1) = pi * (r7*r7 - r6*r6) * (zu7 + zl7)
      v3vol(2) = pi * r6*r6 * (zu7 - zu6)
      v3vol(3) = pi * r6*r6 * (zl7 - zl6)

      !  Material volumes
      do i = 0,maxmat
        chmatv(i) = max(0.0D0, chvol * chmatf(i))
        do j = 1,3
          fwmatv(j,i) = max(0.0D0, fwvol(j) * fwmatf(j,i))
          v1matv(j,i) = max(0.0D0, v1vol(j) * v1matf(j,i))
          blmatv(j,i) = max(0.0D0, blvol(j) * blmatf(j,i))
          v2matv(j,i) = max(0.0D0, v2vol(j) * v2matf(j,i))
          shmatv(j,i) = max(0.0D0, shvol(j) * shmatf(j,i))
          v3matv(j,i) = max(0.0D0, v3vol(j) * v3matf(j,i))
        end do
      end do

      !  First wall area
      fwarea = 2.0D0*pi*r1*( (zu1 + zl1) + r1*sqrt(2.0D0) )

    end subroutine sombld

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    subroutine hylbld
      !! Routine to create the build of an inertial fusion energy
      !! device, based on the design of the HYLIFE-II study,
      !! and to calculate the material volumes for the device core
      !! author: P J Knight, CCFE, Culham Science Centre
      !! None
      !! This routine constructs the build of an inertial fusion energy
      !! device, based on the design of the HYLIFE-II study, and to calculate
      !! the material volumes for the device core.
      !! Moir, Fusion Technology vol.21 (1992) 1475
      !! F/MI/PJK/LOGBOOK12, p.57
      !!       !
      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      use constants, only: pi
      use build_variables, only: fwarea
      use ife_variables, only: chvol, fwvol, v1vol, blvol, v2vol, shvol, v3vol, &
        maxmat, chmatf, fwmatf, blmatf, v1matf, v2matf, shmatf, v3matf, flirad

      implicit none

      !  Local variables

      real(dp), parameter :: third = 1.0D0/3.0D0
      real(dp) :: chcylh,ddz,dvol
      integer :: i,j

      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      !  Radial build
      r1 = chrad
      r2 = r1 + fwdr
      r3 = r2 + v1dr
      r4 = r3 + bldr
      r5 = r4 + v2dr
      r6 = r5 + shdr
      r7 = r6 + v3dr

      !  Vertical build (below midplane)
      zl1 = chdzl
      zl2 = zl1 + fwdzl
      zl3 = zl2 + v1dzl
      zl4 = zl3 + bldzl
      zl5 = zl4 + v2dzl
      zl6 = zl5 + shdzl
      zl7 = zl6 + v3dzl

      !  Vertical build (above midplane)
      zu1 = chdzu
      zu2 = zu1 + fwdzu
      zu3 = zu2 + v1dzu
      zu4 = zu3 + bldzu
      zu5 = zu4 + v2dzu
      zu6 = zu5 + shdzu
      zu7 = zu6 + v3dzu

      !  The HYLIFE-II chamber is assumed to be mostly cylindrical, but
      !  with a conical region below the midplane that causes the Flibe
      !  to flow downwards and outwards towards the outlet.

      !  Component volumes
      !  The following notation applies below:
      !  J=1 : side part
      !  J=2 : top part
      !  J=3 : bottom part
      chvol = pi * r1*r1 * ( (zu1 + zl5) - third*(zl5 - zl1) )

      !  First wall
      !  FLIRAD is the radius of the Flibe inlet
      fwvol(1) = pi * (r2*r2 - r1*r1) * (zu2 + zl5)
      fwvol(2) = pi * (r1*r1 - flirad*flirad) * (zu2 - zu1)
      fwvol(3) = third * pi * ( r2*r2*(zl5-zl1) - r1*r1*(zl5-zl2) )

      !  First void
      v1vol(1) = pi * (r3*r3 - r2*r2) * (zu2 + zl3)
      v1vol(2) = pi * (r4*r4 - flirad*flirad) * (zu3 - zu2)
      v1vol(3) = third * pi * r1*r1 * (zl3 - zl2)

      !  Blanket
      blvol(1) = pi * (r4*r4 - r3*r3) * (zu2 + zl3)
      blvol(2) = pi * (r4*r4 - flirad*flirad) * (zu4 - zu3)
      blvol(3) = pi * r4*r4 * (zl4 - zl3)

      !  Second void
      v2vol(1) = pi * (r5*r5 - r4*r4) * (zu4 + zl4)
      v2vol(2) = pi * (r5*r5 - flirad*flirad) * (zu5 - zu4)
      v2vol(3) = pi * r5*r5 * (zl5 - zl4)

      !  Shield
      shvol(1) = pi * (r6*r6 - r5*r5) * (zu5 + zl5)
      shvol(2) = pi * r6*r6 * (zu6 - zu5)
      shvol(3) = pi * r6*r6 * (zl6 - zl5)

      !  Third void
      v3vol(1) = pi * (r7*r7 - r6*r6) * (zu6 + zl6)
      v3vol(2) = pi * r7*r7 * (zu7 - zu6)
      v3vol(3) = pi * r7*r7 * (zl7 - zl6)

      !  Material volumes
      do i = 0,maxmat
        chmatv(i) = max(0.0D0, chvol * chmatf(i))
        do j = 1,3
          fwmatv(j,i) = max(0.0D0, fwvol(j) * fwmatf(j,i))
          v1matv(j,i) = max(0.0D0, v1vol(j) * v1matf(j,i))
          blmatv(j,i) = max(0.0D0, blvol(j) * blmatf(j,i))
          v2matv(j,i) = max(0.0D0, v2vol(j) * v2matf(j,i))
          shmatv(j,i) = max(0.0D0, shvol(j) * shmatf(j,i))
          v3matv(j,i) = max(0.0D0, v3vol(j) * v3matf(j,i))
        end do
      end do

      !  First wall area
      fwarea = 2.0D0 * pi * r1 * (zu1 + zl5)
      fwarea = fwarea + pi * (r1*r1 - flirad*flirad)
      fwarea = fwarea + pi * r1 * sqrt(r1*r1 + (zl3-zl1)**2)

    end subroutine hylbld

  end subroutine ifebld

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine ifestr
    !! Routine to calculate the support structural masses for the core of
    !! an Inertial Fusion Energy power plant
    !! author: P J Knight, CCFE, Culham Science Centre
    !! This routine calculates the support structural masses for the core of
    !! an Inertial Fusion Energy power plant.
    !! <P>In fact, the output masses are all trivially zero, as they are
    !! magnetic fusion specific. F/MI/PJK/LOGBOOK12, p.87
    !!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    use structure_variables, only: aintmass, clgsmass, coldmass, fncmass, gsmass

    implicit none

    !  Set all outputs to zero, as they are magnetic fusion specific
    aintmass = 0.0D0
    clgsmass = 0.0D0
    coldmass = 0.0D0
    fncmass  = 0.0D0
    gsmass   = 0.0D0

  end subroutine ifestr

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine ifetgt()
    !! Routine to calculate the power requirements of the target
    !! delivery system and the target factory
    !! author: P J Knight, CCFE, Culham Science Centre
    !! This routine calculates the power requirements of the target
    !! delivery system and the target factory, for an Inertial
    !! Fusion Energy power plant.
    !! F/MI/PJK/LOGBOOK12, pp.87-88
    !!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    use ife_variables, only: tfacmw, ptargf, reprat

    implicit none

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !  Target delivery system power (MWe) - effectively negligible

    !  tdspmw = 1.0D-2

    !  Target factory power (MWe)
    !  Assumed to scale with repetition rate (not quite linearly)
    tfacmw = ptargf * (reprat/6.0D0)**0.7D0

  end subroutine ifetgt

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine ifefbs(outfile,iprint)
    !! Routine to calculate the first wall, blanket and shield volumes,
    !! masses and other parameters, for an Inertial Fusion Energy device
    !! author: P J Knight, CCFE, Culham Science Centre
    !! outfile : input integer : output file unit
    !! iprint : input integer : switch for writing to output file (1=yes)
    !! This routine calculates the first wall, blanket and shield volumes,
    !! masses and other parameters, for an Inertial Fusion Energy device.
    !! F/MI/PJK/LOGBOOK12, p.86
    !! Moir et al., Fusion Technology, vol.25 (1994) p.5
    !!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    use build_variables, only: fwarea
    use error_handling, only: report_error, fdiags
    use process_output, only: ovarre, oheadr
    use cost_variables, only: tlife, abktflnc, cfactr
    use fwbs_variables, only: denstl, fwmass, whtblkt, whtshld, whtblbe, whtblvd, &
      whtblss, wtblli2o, whtblli, bktlife, fwlife
    use ife_variables, only: maxmat, chmatm, chmatv, fwmatm, fwmatv, v1matm, &
      v1matv, blmatm, blmatv, v2matm, v2matv, shmatm, shmatv, v3matm, v3matv, &
      mflibe, fbreed, ifetyp
    use physics_variables, only: wallmw

    implicit none

    !  Arguments
    integer, intent(in) :: outfile,iprint

    !  Local variables

    real(dp) :: den,life
    real(dp), dimension(0:maxmat) :: matden
    integer :: i,j

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !  Material densities
    !  0 = void
    !  1 = steel
    !  2 = carbon
    !  3 = FLiBe (inferred from Moir et al)
    !  4 = Li2O
    !  5 = concrete
    !  6 = helium (at typical coolant temperatures)
    !  7 = xenon (taken as ten times the normal tabulated value)
    !  8 = lithium (liquid, source Wikipedia)

    matden(0) = 0.0D0
    matden(1) = denstl
    matden(2) = 2300.0D0
    matden(3) = 2020.0D0
    matden(4) = 2010.0D0
    matden(5) = 2400.0D0
    matden(6) = 1.517D0
    matden(7) = 55.0D0
    matden(8) = 512.0D0

    !  Material masses
    do i = 0,maxmat
      den = matden(i)
      chmatm(i) = chmatv(i) * den
      do j = 1,3
        fwmatm(j,i) = fwmatv(j,i) * den
        v1matm(j,i) = v1matv(j,i) * den
        blmatm(j,i) = blmatv(j,i) * den
        v2matm(j,i) = v2matv(j,i) * den
        shmatm(j,i) = shmatv(j,i) * den
        v3matm(j,i) = v3matv(j,i) * den
      end do
    end do

    !  Total masses of components (excluding coolant)
    fwmass = 0.0D0
    whtblkt = 0.0D0
    whtshld = 0.0D0
    do i = 1,5
      do j = 1,3
        fwmass = fwmass + fwmatm(j,i)
        whtblkt = whtblkt + blmatm(j,i)
        whtshld = whtshld + shmatm(j,i)
      end do
    end do

    !  Other masses
    whtblbe = 0.0D0
    whtblvd = 0.0D0
    whtblss = 0.0D0
    wtblli2o = 0.0D0
    whtblli = 0.0D0
    do j = 1,3
      whtblss = whtblss + blmatm(j,1)
      wtblli2o = wtblli2o + blmatm(j,4)
      whtblli = whtblli + blmatm(j,8)
    end do

    !  Total mass of FLiBe
    mflibe = chmatm(3)
    do j = 1,3
      mflibe = mflibe + fwmatm(j,3) + v1matm(j,3) + blmatm(j,3) + &
        v2matm(j,3) + shmatm(j,3) + v3matm(j,3)
    end do

    !  A fraction FBREED of the total breeder inventory is outside the
    !  core region, i.e. is in the rest of the heat transport system
    if ((fbreed < 0.0D0).or.(fbreed > 0.999D0)) then
      fdiags(1) = fbreed ; call report_error(26)
    end if

    !  Following assumes that use of FLiBe and Li2O are
    !  mutually exclusive
    mflibe = mflibe / (1.0D0 - fbreed)
    wtblli2o = wtblli2o / (1.0D0 - fbreed)
    whtblli = whtblli / (1.0D0 - fbreed)

    !  Blanket and first wall lifetimes (HYLIFE-II: = plant life)
    if ((ifetyp == 3).or.(ifetyp == 4)) THEN
      life = tlife
    else
      life = min( tlife, abktflnc/(wallmw*cfactr) )
    end if

    bktlife = life
    fwlife = life

    !  Cryostat mass (=zero)
    !cryomass = 0.0D0

    if (iprint == 0) return

    !  Output section
    call oheadr(outfile,'First Wall, Blanket, Shield')
    call ovarre(outfile,'First wall area (m2)','(fwarea)',fwarea)
    call ovarre(outfile,'First wall mass (kg)','(fwmass)',fwmass)
    call ovarre(outfile,'Blanket mass (kg)','(whtblkt)',whtblkt)
    call ovarre(outfile,'Blanket lithium mass (kg)','(whtblli)',whtblli)
    call ovarre(outfile,'Total mass of FLiBe (kg)','(mflibe)',mflibe)
    call ovarre(outfile,'Shield mass (kg)','(whtshld)',whtshld)

  end subroutine ifefbs

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine ifepw1

    !! Routine to calculate the first part of the heat transport
    !! and plant power balance constituents, for an IFE power plant
    !! author: P J Knight, CCFE, Culham Science Centre
    !! None
    !! This routine calculates the first part of the heat transport
    !! and plant power balance constituents, for an IFE power plant.
    !! F/MI/PJK/LOGBOOK12, pp.67,89
    !! Bourque et al., Fusion Technology vol.21 (1992) 1465
    !!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    use fwbs_variables, only: emult, fhole, pnucblkt, pnucshld, pnucloss
    use heat_transport_variables, only: priheat, pthermmw, pfwdiv, nphx, pinjwp, &
      pinjht, crypmw, helpow
    use ife_variables, only: pdrive, ifetyp, etadrv, pifecr
    use physics_variables, only: fusion_power

    implicit none

    !  Local variables

    real(dp) :: pdrvmw

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !  Driver power reaching target (MW)

    pdrvmw = 1.0D-6 * pdrive

    !  Primary nuclear heating (MW)
    !  Total thermal power removed from fusion core

    priheat = emult * fusion_power

    !  Useful (high-grade) thermal power (MW)

    pthermmw = priheat * (1.0D0-fhole)

    !  Assume 0.24 of thermal power is intercepted by the first wall
    !  (Bourque et al)
    !  HYLIFE-II case: Assume FLiBe flows intercept all fusion power
    !  and provide the energy multiplication as though it were a
    !  conventional blanket

    if ((ifetyp /= 3).and.(ifetyp /= 4)) then
       pfwdiv = 0.24D0 * pthermmw
       pnucblkt = pthermmw - pfwdiv
    else
       pfwdiv = 0.0D0
       pnucblkt = pthermmw
    end if
    pnucshld = 0.0D0

    !  Lost fusion power (MW)

    pnucloss = priheat - pthermmw  !  = priheat*fhole

    !  Number of primary heat exchangers

    !rnphx = max(2.0D0, (pthermmw/400.0D0 + 0.8D0) )
    nphx = ceiling(pthermmw/1000.0D0)

    !  Secondary heat (some of it... rest calculated in IFEPW2)

    !  Wall plug driver power (MW)

    pinjwp = pdrvmw/etadrv

    !  Waste driver power (MW)

    pinjht = pinjwp - pdrvmw

    !  Cryogenic power (MW)
    !  Cryogenic temperature is assumed to be 4.5K

    crypmw = pifecr
    helpow = 1.0D6 * crypmw * (0.13D0 * 4.5D0)/(293.0D0 - 4.5D0)

  end subroutine ifepw1

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine ifeacp(outfile,iprint)

    !! Routine to calculate AC power requirements for an IFE power plant
    !! author: P J Knight, CCFE, Culham Science Centre
    !! outfile : input integer : output file unit
    !! iprint : input integer : switch for writing to output file (1=yes)
    !! This routine calculates the AC power requirements for an IFE power plant.
    !! F/MI/PJK/LOGBOOK12, p.68
    !!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    use buildings_variables, only: efloor
    use process_output, only: oheadr, ovarre, oblnkl
    use heat_transport_variables, only: baseel, pwpm2, pacpmw, crypmw, vachtmw, &
      trithtmw, pinjwp, fcsht, tlvpmw
    use ife_variables, only: tdspmw, tfacmw, htpmw_ife, reprat, lipmw, ifetyp

    implicit none

    ! Arguments
    integer, intent(in) :: iprint,outfile

    !  Local variables

    real(dp), save :: basemw,pmwpm2

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !  Facility base load, MW (loads not dependent on floor area)

    basemw = baseel * 1.0D-6

    !  Power needed per floor area, MW/m2

    pmwpm2 = pwpm2 * 1.0D-6

    !  Total pulsed power system load, MW

    pacpmw = crypmw + vachtmw + tdspmw + tfacmw + &
         (htpmw_ife*reprat/6.0D0) + trithtmw + pinjwp + basemw + &
         (efloor*pmwpm2) + lipmw

    !  Total baseline power to facility loads, MW

    fcsht  = basemw + (efloor*pmwpm2)

    !  Estimate of the total low voltage power, MW

    tlvpmw = fcsht + trithtmw + (htpmw_ife*reprat/6.0D0) + vachtmw + &
         0.5D0*crypmw + tfacmw

    if (iprint == 0) return

    !  Output section
    call oheadr(outfile,'AC Power')

    call ovarre(outfile,'Facility base load (MW)','(basemw)',basemw)
    call ovarre(outfile,'Total floor space (m2)','(efloor)',efloor)
    call ovarre(outfile,'Power/floor area (MW/m2)','(pmwpm2)',pmwpm2)
    call ovarre(outfile,'Driver power supplies (MW)','(pinjwp)', &
         pinjwp)
    call ovarre(outfile,'Target delivery system (MW)','(tdspmw)', &
         tdspmw)
    call ovarre(outfile,'Target factory (MW)','(tfacmw)', &
         tfacmw)
    call ovarre(outfile,'Tritium processing plant (MW)','(trithtmw)', &
         trithtmw)
    call ovarre(outfile,'Vacuum pump motors (MW)','(vachtmw)',vachtmw)
    call ovarre(outfile,'Cryogenic comp motors (MW)','(crypmw)',crypmw)
    call ovarre(outfile,'Heat transport system pump motors (MW)', &
         '(htpmw_ife*reprat/6)',htpmw_ife*reprat/6.0D0)
    if (ifetyp.eq.4) then
        call ovarre(outfile,'Lithium Pump Power (MW)','(lipmw)',lipmw)
    end if
    call oblnkl(outfile)
    call ovarre(outfile,'Total pulsed power (MW)','(pacpmw)',pacpmw)
    call ovarre(outfile,'Total base power reqd at all times (MW)', &
         '(fcsht)',fcsht)
    call ovarre(outfile,'Total low voltage power (MW)','(tlvpmw)',tlvpmw)

  end subroutine ifeacp

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine ifepw2(outfile,iprint)

    !! Routine to calculate the rest of the IFE heat transport
    !! and plant power balance constituents, not already calculated in
    !! IFEPW1 or IFEACP
    !! author: P J Knight, CCFE, Culham Science Centre
    !! outfile : input integer : output file unit
    !! iprint : input integer : switch for writing to output file (1=yes)
    !! This routine calculates the rest of the IFE heat transport
    !! and plant power balance constituents, not already calculated in
    !! routines <A HREF="ifepw1.html">IFEPW1</A> or
    !! <A HREF="ifeacp.html">IFEACP</A>.
    !! F/MI/PJK/LOGBOOK12, p.67
    !!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    use cost_variables, only: ireactor, ipnet
    use process_output, only: oheadr, osubhd, ovarre, oblnkl, ovarin
    use fwbs_variables, only: pnucloss, emult, tbr, pnucblkt
    use heat_transport_variables, only: fachtmw, fcsht, psechtmw, pinjht, &
      vachtmw, trithtmw, crypmw, pgrossmw, pthermmw, etath, fgrosbop, precircmw, &
      pacpmw, pnetelmw, pinjwp, pfwdiv, nphx
    use ife_variables, only: tdspmw, tfacmw, htpmw_ife, fauxbop, ifetyp, taufall

    implicit none

    !  Arguments
    integer, intent(in) :: outfile,iprint

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !  Facility heat removal (fcsht calculated in IFEACP)
    fachtmw = fcsht

    !  Total secondary heat
    psechtmw = pinjht + pnucloss + fachtmw + vachtmw + trithtmw + &
      tdspmw + tfacmw + crypmw + htpmw_ife

    !  Calculate powers relevant to a power-producing plant
    if (ireactor == 1) then

      !  Gross electric power
      pgrossmw = pthermmw * etath

      !  Balance of plant recirculating power fraction
      fgrosbop = min( 0.5D0, ( fauxbop/(pgrossmw/1000.0D0)**0.6D0) )

      !  Total recirculating power
      precircmw = (fgrosbop*pgrossmw) + pacpmw

      !  Net electric power
      pnetelmw = pgrossmw - precircmw

    end if

    if (iprint == 0) return

    !  Output section
    call oheadr(outfile,'Power / Heat Transport')
    call ovarre(outfile,'Fusion power escaping via holes (MW)', '(pnucloss)', pnucloss)
    call ovarre(outfile,'Power multiplication factor','(emult)',emult)
    if (ifetyp.eq.4) then
        call ovarre(outfile,'Tritium Breeding Ratio','(tbr)',tbr)
        call ovarre(outfile,'Lithium Fall Time (s)','(taufall)',taufall)
    end if
    call ovarre(outfile,'Driver wall plug power (MW)','(pinjwp)', pinjwp)
    call ovarre(outfile,'First wall nuclear heating (MW)','(pfwdiv)', pfwdiv)
    call ovarre(outfile,'Blanket nuclear heating (MW)','(pnucblkt)', pnucblkt)
    call ovarre(outfile,'Primary heat (MW)','(pthermmw)', pthermmw)
    call ovarre(outfile,'Secondary heat (MW)','(psechtmw)', psechtmw)
    call oblnkl(outfile)
    call ovarre(outfile,'Heat removal from driver power (MW)', '(pinjht)', pinjht)
    call ovarre(outfile,'Heat removal from cryogenic plant (MW)', '(crypmw)', crypmw)
    call ovarre(outfile,'Heat removal from vacuum pumps (MW)', '(vachtmw)', vachtmw)
    call ovarre(outfile,'Heat removal from target factory (MW)', '(tfacmw)', tfacmw)
    call ovarre(outfile,'Heat removal from delivery system (MW)', '(tdspmw)', tdspmw)
    call ovarre(outfile,'Heat removal from tritium plant (MW)', '(trithtmw)', trithtmw)
    call ovarre(outfile,'Heat removal from facilities (MW)', '(fachtmw)', fachtmw)
    call ovarin(outfile,'Number of primary heat exchangers', '(nphx)', nphx)

    if (ireactor /= 1) return
    call osubhd(outfile,'Reactor powers :')
    call ovarre(outfile,'Gross electric power (MW)', '(pgrossmw)', pgrossmw)
    call ovarre(outfile,'Net electric power (MW)', '(pnetelmw)', pnetelmw)
    call ovarre(outfile,'Balance of plant aux. power fraction', '(fgrosbop)', fgrosbop)

  end subroutine ifepw2

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine ifevac
    !! Routine to calculate parameters of the vacuum system for an
    !! Inertial Fusion Energy power plant
    !! author: P J Knight, CCFE, Culham Science Centre
    !! This routine calculates the parameters of the vacuum system for an
    !! Inertial Fusion Energy power plant.
    !! <P>The calculated values are hard-wired; they are based loosely
    !! on those for a tokamak of 6m major radius. F/MI/PJK/LOGBOOK12, p.87
    !!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    use vacuum_variables, only: dlscal, nvduct, vacdshm, vcdimax, vpumpn

    implicit none

    dlscal = 2.0D0
    nvduct = 16
    vacdshm = 0.0D0
    vcdimax = 0.3D0
    vpumpn = 32

  end subroutine ifevac

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine ifebdg(outfile,iprint)
    !! Routine to calculate the volumes of the buildings required for
    !! an Inertial Fusion Energy power plant
    !! author: P J Knight, CCFE, Culham Science Centre
    !! outfile : input integer : output file unit
    !! iprint : input integer : switch for writing to output file (1=yes)
    !! This routine calculates the volumes of the buildings required for
    !! an Inertial Fusion Energy power plant. The method is based
    !! closely on that for tokamaks etc. in routine
    !! <A HREF="bldgs.html">BLDGS</A>.
    !! F/MI/PJK/LOGBOOK12, p.87
    !!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    use buildings_variables, only: wrbi, rbwt, rbrt, fndt, trcl, hcwt, hccl, &
      wgt2, stcl, pibv, efloor, triv, conv, admv, shov, admvol, convol, elevol, &
      rbvol, rmbvol, shovol, volrci, wsvol, volnucb
    use process_output, only: oheadr, ovarre
    use fwbs_variables, only: whtshld
    use heat_transport_variables, only: helpow
    use ife_variables, only: zl7, zu7, r7, zl6, zu6, r6

    implicit none

    !  Arguments
    integer, intent(in) :: outfile,iprint

    !  Local variables

    real(dp) :: cran,dcl,dcw,fac2,fac3,hcl,hcw,hrbi,rbh,rbl,rbw, &
         rmbh,rmbl,rmbw,rwl,rww,shh,tch,tcl,tcw,wgts,wsa
    real(dp), save :: cryv,elev,rbv,rmbv,vrci,wsv

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !  Reactor building
    !  ================

    !  Total internal height
    hrbi = zl7 + zu7

    !  Distance from centre of device to wall
    wrbi = r7

    !  Internal volume (square floor)
    vrci = (2.0D0 * wrbi)**2 * hrbi

    !  External dimensions
    !  RBWT = wall thickness
    !  RBRT = roof thickness
    !  FNDT = foundation thickness
    rbw = 2.0D0 * (r7 + rbwt)
    rbl = rbw
    rbh = hrbi + rbrt + fndt

    !  External volume
    rbv = rbw * rbl * rbh

    !  Maintenance building
    !  ====================

    !  The reactor maintenance building includes the hot cells, the
    !  decontamination chamber, the transfer corridors, and the waste
    !  treatment building.  The dimensions of these areas are scaled
    !  from a reference (tokamak) design based on the shield sector size.

    !  Shield height

    shh = zl6 + zu6

    !  Transport corridor size

    tcw = r6 + 4.0D0*trcl
    tcl = 5.0D0*tcw + 2.0D0*hcwt

    !  Decontamination cell size

    dcw = 2.0D0*tcw + 1.0D0
    dcl = 2.0D0*tcw + 1.0D0

    !  Hot cell size

    hcw = r6 + 3.0D0*hccl + 2.0D0
    hcl = 3.0D0*r6 + 4.0D0*hccl + tcw

    !  Radioactive waste treatment

    rww = dcw
    rwl = hcl-dcl-hcwt

    !  Maintenance building dimensions

    rmbw = hcw + dcw + 3.0D0*hcwt
    rmbl = hcl + 2.0D0*hcwt

    !  Height

    if (wgt2 > 1.0D0) then
       wgts = wgt2
    else
       wgts = whtshld
    end if
    cran = 9.41D-6*wgts + 5.1D0
    rmbh = 10.0D0 + (zl6+zu6) + trcl + cran + 5.1D0 + stcl + fndt
    tch = shh + stcl + fndt

    !  Volume

    fac2 = 2.8D0
    rmbv = fac2*rmbw*rmbl*rmbh + tcw*tcl*tch

    !  Warm shop and hot cell gallery

    wsa = (rmbw+7.0D0)*20.0D0 + rmbl*7.0D0
    fac3 = 1.9D0
    wsv = fac3 * wsa*rmbh

    !  Cryogenic building volume

    cryv = 55.0D0 * sqrt(helpow)

    !  Electrical building volume
    !  (set equal to power injection (i.e. driver) building volume)

    elev = pibv

    !  Calculate effective floor area for ac power module

    efloor = (rbv+rmbv+wsv+triv+elev+conv+cryv+admv+shov)/6.0D0

    !  Convert local into global variables

    admvol = admv
    convol = conv
    elevol = elev
    rbvol  = rbv
    rmbvol = rmbv
    shovol = shov
    volrci = vrci
    wsvol  = wsv

    !  Total volume of nuclear buildings

    volnucb = ( vrci + rmbv + wsv + triv + cryv )

    if (iprint == 0) return

    !  Output section

    call oheadr(outfile,'Plant Buildings System')
    call ovarre(outfile,'Internal volume of reactor building (m3)', &
         '(vrci)',vrci)
    call ovarre(outfile,'Dist from device centre to bldg wall (m)', &
         '(wrbi)',wrbi)
    call ovarre(outfile,'Effective floor area (m2)','(efloor)',efloor)
    call ovarre(outfile,'Reactor building volume (m3)','(rbv)',rbv)
    call ovarre(outfile,'Reactor maintenance building volume (m3)', &
         '(rmbv)',rmbv)
    call ovarre(outfile,'Warmshop volume (m3)','(wsv)',wsv)
    call ovarre(outfile,'Tritium building volume (m3)','(triv)',triv)
    call ovarre(outfile,'Electrical building volume (m3)','(elev)',elev)
    call ovarre(outfile,'Control building volume (m3)','(conv)',conv)
    call ovarre(outfile,'Cryogenics building volume (m3)','(cryv)',cryv)
    call ovarre(outfile,'Administration building volume (m3)','(admv)', &
         admv)
    call ovarre(outfile,'Shops volume (m3)','(shov)',shov)
    call ovarre(outfile,'Total volume of nuclear buildings (m3)', &
         '(volnucb)',volnucb)

  end subroutine ifebdg

end module ife_module