ifephy Subroutine

  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: powfmw, 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)
        powfmw = 1.0D-6 * pdrive * gain

    else
        !  Driver Power
        reprat = rrin
        pdrive = reprat * edrive
        !  Gain
        powfmw = pfusife
        gain = powfmw / (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 = powfmw * 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 = powfmw * 0.5D0*sang / fwarea

    else
       wallmw = powfmw / 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)','(powfmw)',powfmw)
    call ovarre(outfile,'Neutron wall load (MW/m2)','(wallmw)',wallmw)

  end subroutine ifephy