amns_external_functions.F90

Go to the documentation of this file.

module amns_external_functions

  use itm_types
  implicit none

  type fun_err_t
     integer                      :: ierr
     character (len=128)          :: cerr
  end type fun_err_t

contains

  subroutine nuclear_data_1001(function_parameters, x, f, with_warnings, fun_err)
    implicit none
    real (R8),         intent(in)  :: function_parameters(:,:), x(:)
    real (R8),         intent(out) :: f(:)
    logical,           intent(in)  :: with_warnings
    type (fun_err_t),  intent(out) :: fun_err
    integer :: i, j
    real (R8) :: e, s

    fun_err%ierr = 0
    fun_err%cerr = ''

    if(size(function_parameters,1).ne.12) then
       write(fun_err%cerr,*) "nuclear_data_1001: Expected the length of 'function_parameters' to be 12 in 'nuclear_data_1001' and not ", &
            size(function_parameters)
       fun_err%ierr = -1
       return
    endif
    if(size(x).ne.size(f)) then
       write(fun_err%cerr,*) "nuclear_data_1001: Expected the input and output vectors to be the same size in 'nuclear_data_1001'"
       fun_err%ierr = -1
       return
    endif

    do i=1, size(x)
       e = x(i) * 1.0e-3_r8
       j=1
       if(e.lt.function_parameters(11,j) .and. with_warnings) then
          write(*,*) 'extrapolating  below the desired range'
       endif
       do while(e.gt.function_parameters(12,j).and. &
            j.lt.size(function_parameters,2))
          j=j+1
       enddo
       if(e.gt.function_parameters(12,j)) then
          if (with_warnings) then
             write(*,*) 'extrapolating  above the desired range'
             write(*,*) j, x(i), function_parameters(12,j)
             write(*,*) 'Taking the boundary value when calculating S!'
          endif
          e = function_parameters(12,j)
       endif
       s = ( function_parameters(2,j) &
            & + e*(function_parameters(3,j) &
            & + e*(function_parameters(4,j) &
            & + e*( function_parameters(5,j) &
            & +e*function_parameters(6,j) ) ) ) ) &
            & /( 1 + e*(function_parameters(7,j)&
            &  + e*(function_parameters(8,j) &
            & + e*( function_parameters(9,j) &
            & +e*function_parameters(10,j) ) ) ) )
       ! We clamp the energy for S (the nuclear fusion probability).
       ! However, there is no need to clamp the energy in denominator.
       ! The denominator roughly describes the probability of the reactants to
       ! overcome the Coulomb barrier.
       e = x(i) * 1.0e-3_r8
       f(i) = s / (e * exp(function_parameters(1,j)/sqrt(e))) * 1.0e-31_r8
    enddo
  end subroutine nuclear_data_1001

  subroutine nuclear_data_1002(function_parameters, Tin, f, with_warnings, fun_err)
    implicit none
    real (R8), intent(in)  :: function_parameters(:,:), tin(:)
    real (R8), intent(out) :: f(:)
    logical,           intent(in)  :: with_warnings
    type (fun_err_t),  intent(out) :: fun_err
    integer :: it ! index for the different energy vaulues
    integer :: ir ! index for possible different temperature ranges

    real (R8) :: t
    real (R8) :: theta,xi,mrc2,bg
    real (R8) :: c(1:7)

    fun_err%ierr = 0
    fun_err%cerr = ''

    if(size(function_parameters,1).ne.11) then
       write(fun_err%cerr,*) "nuclear_data_1002: Expected the length of 'function_parameters' to be 11 in 'nuclear_data_1002' and not ", &
            size(function_parameters)
       fun_err%ierr = -1
       return
    endif
    if(size(tin).ne.size(f)) then
       write(fun_err%cerr,*) "nuclear_data_1002: Expected the input and output vectors to be the same size in 'nuclear_data_1002'"
       fun_err%ierr = -1
       return
    endif




    do it=1, size(tin)
       t=tin(it)*1e-3 ! temperature in keV
       ir=1
       if(t.lt.function_parameters(10,ir) .and. with_warnings) then
          write(*,*) 'extrapolating  below the desired range'
       endif
       do while(t.gt.function_parameters(11,ir).and. &
            ir.lt.size(function_parameters,2))
          ir=ir+1
       enddo
       if(t.gt.function_parameters(11,ir) .and. with_warnings) then
          write(*,*) 'extrapolating  above the desired range'
! dpc fix: changed 12 to 11
          write(*,*) ir, tin(it), function_parameters(11,ir)
       endif

       bg  =function_parameters(1,   ir)
       mrc2=function_parameters(2,   ir)
       c   =function_parameters(3:9, ir)
       ! (/Tmin,Tmax/) = function_parameters(10:11, iR)

       theta=t/(1-t*(c(2)+t*(c(4)+t*c(6)))&
           /(1+t*(c(3)+t*(c(5)+t*c(7)))))
       xi=(bg**2/(4*theta))**(1.0_r8/3.0_r8)
       f(it) = c(1)*theta*sqrt(xi/(mrc2*t**3))&
           *exp(-3*xi)*1.e-6 ! from cm^3/s to m^3s/



    end do

  end subroutine nuclear_data_1002

  subroutine rct_data_1003(function_parameters, Tin, f, with_warnings, fun_err)
    implicit none
    real (R8), intent(in)  :: function_parameters(:), tin(:)
    real (R8), intent(out) :: f(:)
    logical,           intent(in)  :: with_warnings
    type (fun_err_t),  intent(out) :: fun_err
    integer                :: it

    fun_err%ierr = 0
    fun_err%cerr = ''

    if(size(function_parameters,1).ne.2) then
       write(fun_err%cerr,*) "rct_data_1003: Expected the length of 'function_parameters' to be 2 in 'rct_data_1003' and not ", &
            size(function_parameters)
       fun_err%ierr = -1
       return
    endif
    if(size(tin).ne.size(f)) then
       write(fun_err%cerr,*) "rct_data_1003: Expected the input and output vectors to be the same size in 'rct_data_1003'"
       fun_err%ierr = -1
       return
    endif

    do it=1, size(tin)
       f(it) = function_parameters(1)*(1.0_r8 + function_parameters(2)* log(1.0_r8/tin(it)))**2*1e-20_r8
    end do

  end subroutine rct_data_1003

  subroutine sputter_data_1004(function_parameters, energy_arr, angle_arr, yield_arr, with_warnings, fun_err)
    use eckstein_yields
    implicit none
    real (R8), intent(in)  :: function_parameters(:,:), energy_arr(:), angle_arr(:)
    real (R8), intent(out) :: yield_arr(:)
    logical,           intent(in)  :: with_warnings
    type (fun_err_t),  intent(out) :: fun_err

    real (R8)              :: energy, angle, yield
    real(R8)               :: matchanglee0, matchanglee1
    real(R8)               :: m1, m2, z1, z2, q, lambda, mu, eth, f,b,c,esp
    integer                :: numtab, ianglee0, numpars, boundianglee0, ipars
    integer                :: id, ia
    logical, parameter     :: debug=.true.

    integer                :: have_angle_data

    fun_err%ierr = 0
    fun_err%cerr = ''

    numtab = ubound(function_parameters, 2)
    if(debug .and. with_warnings) then
       write(*,*) numtab
    endif
    do id = lbound(energy_arr,1), ubound(energy_arr,1)

       energy = energy_arr(id)
       angle  = angle_arr(id)

       if(debug .and. with_warnings) then
          write(*,*) energy, angle
       endif

       matchanglee0 = 0.0
       boundianglee0 = -1
       have_angle_data = -1
       do ia = lbound(function_parameters, 2), ubound(function_parameters, 2)
          matchanglee1 = function_parameters(3, ia)
          if(debug .and. with_warnings) then
             write(*,*) matchanglee0, energy, matchanglee1
          endif
          if((matchanglee0.le.energy).and.(matchanglee1.ge.energy)) then
             boundianglee0 = ia
             have_angle_data = 1
             exit
          end if
          if (matchanglee1 < matchanglee0)  then
             if (matchanglee1 < matchanglee0)  then
                if(debug .and. with_warnings) then
                   write(*, *) 'WARNING::sputteryield-> no angular reflection yield data available, returning value for perpendicular impact'
                endif
                have_angle_data = -1
                exit
             else
                write(fun_err%cerr,*) 'sputter_data_1004: sputteryield-> angular dependence energies are not sorted, rebuild the database'
                fun_err%ierr = -1
                return
             endif
          end if
          matchanglee0 = matchanglee1
       end do
                
       if(boundianglee0.lt.0) then
          if (debug .and. with_warnings) then
             write(*,*) 'sputter_data_1004: No Bound found for: ', energy, ' angle: ', angle, 'using boundianglee0 = ', numtab
          end if
          boundianglee0 = ubound(function_parameters, 2)
       end if
                
       m1     = function_parameters( 4, boundianglee0)
       m2     = function_parameters( 5, boundianglee0)
       z1     = function_parameters( 1, boundianglee0)
       z2     = function_parameters( 2, boundianglee0)
       q      = function_parameters( 6, boundianglee0)
       lambda = function_parameters( 9, boundianglee0)
       eth    = function_parameters( 8, boundianglee0)
       mu     = function_parameters( 7, boundianglee0)
       f      = function_parameters(10, boundianglee0)
       b      = function_parameters(11, boundianglee0)
       c      = function_parameters(12, boundianglee0)
       esp    = function_parameters(13, boundianglee0)
                
       yield = 0.0
       if (energy .gt. eth) then
          if (have_angle_data .gt. 0) then
             yield = seyield(energy, m1, m2, z1, z2, q, lambda, mu, eth) * sayield(energy, angle, f, b, c, esp)
          else
             yield = seyield(energy, m1, m2, z1, z2, q, lambda, mu, eth)
          endif
       end if

       yield_arr(id) = yield

    enddo

  end subroutine sputter_data_1004

  subroutine reflect_data_1005(function_parameters, energy_arr, angle_arr, refl_arr, with_warnings, fun_err)
    use eckstein_yields
    implicit none
    real (R8), intent(in)  :: function_parameters(:,:), energy_arr(:), angle_arr(:)
    real (R8), intent(out) :: refl_arr(:)
    logical,           intent(in)  :: with_warnings
    type (fun_err_t),  intent(out) :: fun_err

    real (R8)              :: energy, angle, refl
    real(R8)               :: m1, m2, z1, z2, a1, a2, a3, a4, e1, esp, anglee0, anglee1, c1, c2, c3, c4
    real(R8)               :: c1from, c2from, c3from, c4from
    real(R8)               :: c1to, c2to, c3to, c4to
    real(R8)               :: matchanglee0, matchanglee1
    integer                :: numtab, ianglee0, numpars, boundianglee0, boundianglee1, ipars
    integer                :: id, ia
    logical, parameter     :: debug=.true.

    integer                :: have_angle_data

    fun_err%ierr = 0
    fun_err%cerr = ''

    numtab = ubound(function_parameters, 2)
    do id = lbound(energy_arr,1), ubound(energy_arr,1)

       energy = energy_arr(id)
       angle  = angle_arr(id)

       matchanglee0 = 0.0
       boundianglee0 = -1
       boundianglee1 = -1
       have_angle_data = -1
       do ia = lbound(function_parameters, 2), ubound(function_parameters, 2)
          matchanglee1 = function_parameters(11, ia)
          if(debug .and. with_warnings) then
             write(*,*) matchanglee0, energy, matchanglee1
          endif
          if((matchanglee0.le.energy).and.(matchanglee1.ge.energy)) then
             if (ia .gt. 1) then
                if(ia .lt. numtab) then
                   boundianglee0 = ia - 1
                   boundianglee1 = ia
                else
                   boundianglee0 = numtab
                   boundianglee1 = numtab
                end if
             else
                boundianglee0 = 1
                boundianglee1 = 1
             end if
             have_angle_data = 1
             exit
          end if
          if (matchanglee1 < matchanglee0)  then
             if (matchanglee1 .lt. 0) then
                if(debug .and. with_warnings) then
                   write(*,*) 'WARNING::reflyield-> no angular reflection yield data available, returning value for perpendicular impact'
                endif
                have_angle_data = -1
                exit
             else
                write(fun_err%cerr,*) 'reflect_data_1005: reflyield-> angular dependence energies are not sorted, rebuild the database'
                fun_err%ierr = -1
                return
             end if
          endif
          matchanglee0 = matchanglee1
       end do
                
       if((boundianglee0.lt.0).or.(boundianglee1.lt.0)) then
          if (matchanglee1.lt.energy) then
             boundianglee0 = numtab
             boundianglee1 = numtab
          else
             boundianglee0 = 1
             boundianglee1 = 1
          end if
          if (debug .and. with_warnings) then
             write(*,*) 'reflect_data_1005: No Bound found for: ', energy, ' angle: ', angle, 'using boundianglee0 = ', boundianglee0, 'and boundianglee1 = ',boundianglee1
          end if
       end if
                
       z1  = function_parameters( 1, boundianglee0)
       z2  = function_parameters( 2, boundianglee0)
       m1  = function_parameters( 3, boundianglee0)
       m2  = function_parameters( 4, boundianglee0)             
       a1  = function_parameters( 5, boundianglee0)
       a2  = function_parameters( 6, boundianglee0)
       a3  = function_parameters( 7, boundianglee0)
       a4  = function_parameters( 8, boundianglee0)
       e1  = function_parameters( 9, boundianglee0)
       esp = function_parameters(10, boundianglee0)
                
       if ((angle.lt.1.0).or.(have_angle_data.eq.-1)) then
          if (debug .and. with_warnings) then
             write(*,*) 'angle = ',angle, ' is small or no angular data available: using perpendicular impact model with energy dependence'
          end if
          refl = reyield(energy, m1, m2, z1, z2, a1, a2, a3, a4)
       else
          if (debug .and. with_warnings) then
             write(*,*) 'angle = ',angle, ' is large using angular dependent model and interpolation of energy dependence'
          end if
          ! Interplolate the bounding table entries
          if(boundianglee0.lt.boundianglee1) then
             anglee0 = function_parameters(11, boundianglee0)
             anglee1 = function_parameters(11, boundianglee1)
             if(debug .and. with_warnings) then
                write(*,*) 'Interpolating angular dependent reflection yield for energy' ,energy,' between values for energies ', anglee0, ' and ', anglee1
             end if
             c1from = function_parameters(12, boundianglee0)
             c2from = function_parameters(13, boundianglee0)
             c3from = function_parameters(14, boundianglee0)
             c4from = function_parameters(15, boundianglee0)

             c1to = function_parameters(12, boundianglee1)
             c2to = function_parameters(13, boundianglee1)
             c3to = function_parameters(14, boundianglee1)
             c4to = function_parameters(15, boundianglee1)

             c1 = c1from + ((energy - anglee0) * (c1to - c1from) / (anglee1 - anglee0))
             c2 = c2from + ((energy - anglee0) * (c2to - c2from) / (anglee1 - anglee0))
             c3 = c3from + ((energy - anglee0) * (c3to - c3from) / (anglee1 - anglee0))
             c4 = c4from + ((energy - anglee0) * (c4to - c4from) / (anglee1 - anglee0))
          else
             anglee0 = function_parameters(11, boundianglee0)
             if(debug .and. with_warnings) then
                write(*,*) 'Energy: ',energy,'  for requested angle: ',angle,' not within db bounds using anglular dependent value at energy = ', anglee0
             end if
             c1 = function_parameters(12, boundianglee0)
             c2 = function_parameters(13, boundianglee0)
             c3 = function_parameters(14, boundianglee0)
             c4 = function_parameters(15, boundianglee0)
          end if
          refl = rayield(angle, c1, c2, c3, c4)
       end if

       refl_arr(id) = min(1.0_r8, refl)

    enddo

  end subroutine reflect_data_1005

  subroutine nuclear_data_1006(function_parameters, x, f, with_warnings, fun_err)
    implicit none
    real (R8),         intent(in)  :: function_parameters(:,:), x(:)
    real (R8),         intent(out) :: f(:)
    logical,           intent(in)  :: with_warnings
    type (fun_err_t),  intent(out) :: fun_err
    integer :: i, j
    real (R8) :: e, s

    fun_err%ierr = 0
    fun_err%cerr = ''

    if(size(function_parameters,1).ne.14) then
       write(fun_err%cerr,*) "nuclear_data_1006: Expected the length of 'function_parameters' to be 14 in 'nuclear_data_1006' and not ", &
            size(function_parameters)
       fun_err%ierr = -1
       return
    endif
    if(size(x).ne.size(f)) then
       write(fun_err%cerr,*) "nuclear_data_1006: Expected the input and output vectors to be the same size in 'nuclear_data_1006'"
       fun_err%ierr = -1
       return
    endif

    do i=1, size(x)
       e = x(i) * 1.0e-3_r8
       j=1
       if(e.lt.function_parameters(13,j) .and. with_warnings) then
          write(*,*) 'extrapolating  below the desired range'
       endif
       do while(e.gt.function_parameters(14,j).and. &
            j.lt.size(function_parameters,2))
          j=j+1
       enddo
       if(e.gt.function_parameters(14,j)) then
          if (with_warnings) then
             write(*,*) 'extrapolating  above the desired range'
             write(*,*) j, x(i), function_parameters(12,j)
             write(*,*) 'Taking the boundary value when calculating S!'
          endif
          e = function_parameters(14,j)
       endif
       s = ( function_parameters(2,j) &
            & + e*(function_parameters(3,j) &
            & + e*(function_parameters(4,j) &
            & + e*(function_parameters(5,j) &
            & + e*(function_parameters(6,j) &
            & + e* function_parameters(7,j) ) ) ) ) ) &
            & /( 1.0_r8 &
            & + e*(function_parameters(8,j) &
            & + e*(function_parameters(9,j) &
            & + e*(function_parameters(10,j) &
            & + e*(function_parameters(11,j) &
            & + e* function_parameters(12,j) ) ) ) ) )
       ! We clamp the energy for S (the nuclear fusion probability).
       ! However, there is no need to clamp the energy in denominator.
       ! The denominator roughly describes the probability of the reactants to
       ! overcome the Coulomb barrier.
       e = x(i) * 1.0e-3_r8
       f(i) = s / (e * exp(function_parameters(1,j)/sqrt(e))) * 1.0e-31_r8
    enddo
  end subroutine nuclear_data_1006

end module amns_external_functions