diff --git a/CMakeLists.txt b/CMakeLists.txt index c63c0ffa..d6367e52 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -118,10 +118,7 @@ LIST(APPEND PROCESS_SRCS pulse_variables.f90 pf_power_variables.f90 dcll.f90 - read_and_get_atomic_data.f90 - plot_radiation.f90 rebco_variables.f90 - read_radiation.f90 primary_pumping_variables.f90 reinke_variables.f90 neoclassics_module.f90 diff --git a/source/fortran/commons.f90 b/source/fortran/commons.f90 deleted file mode 100755 index 2e33899b..00000000 --- a/source/fortran/commons.f90 +++ /dev/null @@ -1,76 +0,0 @@ -module mod_f90_kind - - implicit none - - ! Single byte integer - integer, parameter :: int_def = selected_int_kind(9) - - ! Long byte integer - integer, parameter :: int_long = selected_int_kind(18) - - ! Single precision - integer, parameter :: single = selected_real_kind(6,37) - - ! Double precision - integer, parameter :: double = selected_real_kind(15,300) - - ! logical - integer, parameter :: logical_def = kind(.true.) - - !----------------------------------------------------------------------- - ! Working precision - !----------------------------------------------------------------------- - integer, parameter :: rkind = double ! double - integer, parameter :: r4kind = single ! single - integer, parameter :: ikind = int_def ! integer - integer, parameter :: i8kind = int_long ! integer double - integer, parameter :: lkind = logical_def ! integer - integer, parameter :: skind = 512 ! string - -end module mod_f90_kind - -module param - use mod_f90_kind, only: rkind !EXTERNAL - implicit none - real(rkind), parameter :: pi = 3.14159265358979e0_rkind, & - c = 2.99792458e10_rkind -end module param - -module freq - use mod_f90_kind, only: rkind !EXTERNAL - implicit none - real(rkind), save :: om, xk0 -end module freq - -module mode - use mod_f90_kind, only: ikind !EXTERNAL - implicit none - integer(ikind), save :: imod -end module mode - -module machin - use mod_f90_kind, only: rkind !EXTERNAL - implicit none - real(rkind), save :: rmaj, rmin, sgnm, zsearch - character(len=160), save :: machname -end module machin - -module bsquar - use mod_f90_kind, only: rkind !EXTERNAL - implicit none - real(rkind), save :: conf, conf2, qb2 -end module bsquar - -module linliu - use mod_f90_kind, only: rkind !EXTERNAL - implicit none - save - complex(rkind) :: PCEXNORMO, PCEYNORMO, PCEZNORMO - complex(rkind) :: PCEXNORMX, PCEYNORMX, PCEZNORMX -end module linliu - -module abs_cd - use mod_f90_kind, only: ikind !EXTERNAL - implicit none - integer(ikind), save :: cdroutine, absroutine, iastra, profcalc -end module abs_cd diff --git a/source/fortran/init_module.f90 b/source/fortran/init_module.f90 index 604661dc..e46977b6 100644 --- a/source/fortran/init_module.f90 +++ b/source/fortran/init_module.f90 @@ -30,8 +30,6 @@ subroutine init_all_module_vars use pfcoil_module, only: init_pfcoil_module use physics_module, only: init_physics_module use physics_variables, only: init_physics_variables - use read_and_get_atomic_data, only: init_read_and_get_atomic_data - use read_radiation, only: init_read_radiation use scan_module, only: init_scan_module use sctfcoil_module, only: init_sctfcoil_module use tfcoil_variables, only: init_tfcoil_variables @@ -73,8 +71,6 @@ subroutine init_all_module_vars call init_pfcoil_module call init_physics_module call init_physics_variables - call init_read_and_get_atomic_data - call init_read_radiation call init_scan_module call init_sctfcoil_module call init_tfcoil_variables diff --git a/source/fortran/ode.f90 b/source/fortran/ode.f90 deleted file mode 100644 index 3181232d..00000000 --- a/source/fortran/ode.f90 +++ /dev/null @@ -1,1399 +0,0 @@ -module ode_mod - -#ifndef dp - use, intrinsic :: iso_fortran_env, only: dp=>real64 -#endif - - implicit none - -contains - -subroutine ode ( f, neqn, y, t, tout, relerr, abserr, iflag, work, iwork ) - -!***************************************************************************** -! -!! ODE is the user interface to an ordinary differential equation solver. -! -! Discussion: -! -! ODE integrates a system of NEQN first order ordinary differential -! equations of the form: -! dY(i)/dT = F(T,Y(1),Y(2),...,Y(NEQN)) -! Y(i) given at T. -! The subroutine integrates from T to TOUT. On return, the -! parameters in the call list are set for continuing the integration. -! The user has only to define a new value TOUT and call ODE again. -! -! The differential equations are actually solved by a suite of codes -! DE, STEP, and INTRP. ODE allocates virtual storage in the -! arrays WORK and IWORK and calls DE. DE is a supervisor which -! directs the solution. It calls the routines STEP and INTRP -! to advance the integration and to interpolate at output points. -! -! STEP uses a modified divided difference form of the Adams PECE -! formulas and local extrapolation. It adjusts the order and step -! size to control the local error per unit step in a generalized -! sense. Normally each call to STEP advances the solution one step -! in the direction of TOUT. For reasons of efficiency, DE integrates -! beyond TOUT internally, though never beyond T+10*(TOUT-T), and -! calls INTRP to interpolate the solution at TOUT. An option is -! provided to stop the integration at TOUT but it should be used -! only if it is impossible to continue the integration beyond TOUT. -! -! On the first call to ODE, the user must provide storage in the calling -! program for the arrays in the call list, -! Y(NEQN), WORK(100+21*NEQN), IWORK(5), -! declare F in an external statement, supply the double precision -! SUBROUTINE F ( T, Y, YP ) -! to evaluate dy(i)/dt = yp(i) = f(t,y(1),y(2),...,y(neqn)) -! and initialize the parameters: -! * NEQN, the number of equations to be integrated; -! * Y(1:NEQN), the vector of initial conditions; -! * T, the starting point of integration; -! * TOUT, the point at which a solution is desired; -! * RELERR, ABSERR, the relative and absolute local error tolerances; -! * IFLAG, an indicator to initialize the code. Normal input -! is +1. The user should set IFLAG = -1 only if it is -! impossible to continue the integration beyond TOUT. -! All parameters except F, NEQN and TOUT may be altered by the -! code on output, and so must be variables in the calling program. -! -! On normal return from ODE, IFLAG is 2, indicating that T has been -! set to TOUT, and Y has been set to the approximate solution at TOUT. -! -! If IFLAG is 3, then the program noticed that RELERR or ABSERR was -! too small; the output values of these variables are more appropriate, -! and integration can be resumed by setting IFLAG to 1. -! -! IFLAG is -2 if the user input IFLAG = -1, and the code was able to -! reach TOUT exactly. In that case, the output value of T is TOUT, -! and the output value of Y is the solution at TOUT, which was computed -! directly, and not by interpolation. -! -! Other values of IFLAG generally indicate an error. -! -! Normally, it is desirable to compute many points along the solution -! curve. After the first successful step, more steps may be taken -! simply by updating the value of TOUT and calling ODE again. -! -! Licensing: -! -! This code is distributed under the GNU LGPL license. -! -! Modified: -! -! 02 August 2009 -! -! Author: -! -! Original FORTRAN77 version by Lawrence Shampine, Marilyn Gordon. -! FORTRAN90 version by John Burkardt. -! -! Reference: -! -! Lawrence Shampine, Marilyn Gordon, -! Computer Solution of Ordinary Differential Equations: -! The Initial Value Problem, -! Freeman, 1975, -! ISBN: 0716704617, -! LC: QA372.S416. -! -! Parameters: -! -! Input, external F, the name of a user-supplied routine of the form -! subroutine f ( t, y, yp ) -! real ( kind = 8 ) t -! real ( kind = 8 ) y(neqn) -! real ( kind = 8 ) yp(neqn) -! which accepts input values T and Y(1:NEQN), evaluates the right hand -! sides of the ODE, and stores the result in YP(1:NEQN). -! -! Input, integer ( kind = 4 ) NEQN, the number of equations. -! -! Input/output, real ( kind = 8 ) Y(NEQN), the current solution. -! -! Input/output, real ( kind = 8 ) T, the current value of the independent -! variable. -! -! Input, real ( kind = 8 ) TOUT, the desired value of T on output. -! -! Input, real ( kind = 8 ) RELERR, ABSERR, the relative and absolute error -! tolerances. At each step, the code requires -! abs ( local error ) <= abs ( y ) * relerr + abserr -! for each component of the local error and solution vectors. -! -! Input/output, integer ( kind = 4 ) IFLAG, indicates the status of -! integration. On input, IFLAG is normally 1 (or -1 in the special case -! where TOUT is not to be exceeded.) On normal output, IFLAG is 2. Other -! output values are: -! * 3, integration did not reach TOUT because the error tolerances -! were too small. But RELERR and ABSERR were increased appropriately -! for continuing; -! * 4, integration did not reach TOUT because more than 500 steps were taken; -! * 5, integration did not reach TOUT because the equations appear to -! be stiff; -! * 6, invalid input parameters (fatal error). -! The value of IFLAG is returned negative when the input value is -! negative and the integration does not reach TOUT. -! -! Input/output, real ( kind = 8 ) WORK(100+21*NEQN), workspace. -! -! Input/output, integer ( kind = 4 ) IWORK(5), workspace. -! - implicit none - - integer ( kind = 4 ) neqn - - real ( kind = 8 ) abserr - external f - integer ( kind = 4 ), parameter :: ialpha = 1 - integer ( kind = 4 ), parameter :: ibeta = 13 - integer ( kind = 4 ), parameter :: idelsn = 93 - integer ( kind = 4 ) iflag - integer ( kind = 4 ), parameter :: ig = 62 - integer ( kind = 4 ), parameter :: ih = 89 - integer ( kind = 4 ), parameter :: ihold = 90 - integer ( kind = 4 ) ip - integer ( kind = 4 ), parameter :: iphase = 75 - integer ( kind = 4 ) iphi - integer ( kind = 4 ), parameter :: ipsi = 76 - integer ( kind = 4 ), parameter :: isig = 25 - integer ( kind = 4 ), parameter :: istart = 91 - integer ( kind = 4 ), parameter :: itold = 92 - integer ( kind = 4 ), parameter :: iv = 38 - integer ( kind = 4 ), parameter :: iw = 50 - integer ( kind = 4 ) iwt - integer ( kind = 4 ), parameter :: ix = 88 - integer ( kind = 4 ) iyp - integer ( kind = 4 ) iypout - integer ( kind = 4 ), parameter :: iyy = 100 - integer ( kind = 4 ) iwork(5) - logical nornd - logical phase1 - real ( kind = 8 ) relerr - logical start - real ( kind = 8 ) t - real ( kind = 8 ) tout - real ( kind = 8 ) work(100+21*neqn) - real ( kind = 8 ) y(neqn) - - iwt = iyy + neqn - ip = iwt + neqn - iyp = ip + neqn - iypout = iyp + neqn - iphi = iypout + neqn - - if ( abs ( iflag ) /= 1 ) then - start = ( 0.0D+00 < work(istart) ) - phase1 = ( 0.0D+00 < work(iphase) ) - nornd = ( iwork(2) /= -1 ) - end if - - call de ( f, neqn, y, t, tout, relerr, abserr, iflag, work(iyy), & - work(iwt), work(ip), work(iyp), work(iypout), work(iphi), & - work(ialpha), work(ibeta), work(isig), work(iv), work(iw), work(ig), & - phase1, work(ipsi), work(ix), work(ih), work(ihold), start, & - work(itold), work(idelsn), iwork(1), nornd, iwork(3), iwork(4), & - iwork(5) ) - - if ( start ) then - work(istart) = 1.0D+00 - else - work(istart) = -1.0D+00 - end if - - if ( phase1 ) then - work(iphase) = 1.0D+00 - else - work(iphase) = -1.0D+00 - end if - - if ( nornd ) then - iwork(2) = 1 - else - iwork(2) = -1 - end if - - return -! SJP Issue #835 -! For Intel compliance add "end subroutine" -end subroutine -subroutine de ( f, neqn, y, t, tout, relerr, abserr, iflag, yy, wt, p, yp, & - ypout, phi, alpha, beta, sig, v, w, g, phase1, psi, x, h, hold, start, & - told, delsgn, ns, nornd, k, kold, isnold ) - -!***************************************************************************** -! -!! DE carries out the ODE solution algorithm. -! -! Discussion: -! -! ODE merely allocates storage for DE, to relieve the user of the -! inconvenience of a long call list. Consequently, DE is used as -! described in the comments for ODE. -! -! Licensing: -! -! This code is distributed under the GNU LGPL license. -! -! Modified: -! -! 02 August 2009 -! -! Author: -! -! Original FORTRAN77 version by Lawrence Shampine, Marilyn Gordon. -! FORTRAN90 version by John Burkardt. -! -! Reference: -! -! Lawrence Shampine, Marilyn Gordon, -! Computer Solution of Ordinary Differential Equations: -! The Initial Value Problem, -! Freeman, 1975, -! ISBN: 0716704617, -! LC: QA372.S416. -! -! Parameters: -! -! Input, external F, the name of a user-supplied routine of the form -! subroutine f ( t, y, yp ) -! real ( kind = 8 ) t -! real ( kind = 8 ) y(neqn) -! real ( kind = 8 ) yp(neqn) -! which accepts input values T and Y(1:NEQN), evaluates the right hand -! sides of the ODE, and stores the result in YP(1:NEQN). -! -! Input, integer ( kind = 4 ) NEQN, the number of equations. -! -! Input/output, real ( kind = 8 ) Y(NEQN), the current solution. -! -! Input/output, real ( kind = 8 ) T, the current value of the independent -! variable. -! -! Input, real ( kind = 8 ) TOUT, the desired value of T on output. -! -! Input, real ( kind = 8 ) RELERR, ABSERR, the relative and absolute error -! tolerances. At each step, the code requires -! abs ( local error ) <= abs ( Y ) * RELERR + ABSERR -! for each component of the local error and solution vectors. -! -! Input/output, integer ( kind = 4 ) IFLAG, indicates the status of -! integration. On input, IFLAG is normally 1 (or -1 in the special case -! where TOUT is not to be exceeded.) On normal output, IFLAG is 2. Other -! output values are: -! * 3, integration did not reach TOUT because the error tolerances were -! too small. -! But RELERR and ABSERR were increased appropriately for continuing; -! * 4, integration did not reach TOUT because more than 500 steps were taken; -! * 5, integration did not reach TOUT because the equations appear to be -! stiff; -! * 6, invalid input parameters (fatal error). -! The value of IFLAG is returned negative when the input value is negative -! and the integration does not reach TOUT. -! -! Workspace, real ( kind = 8 ) YY(NEQN), used to hold old solution data. -! -! Input, real ( kind = 8 ) WT(NEQN), the error weight vector. -! -! Workspace, real ( kind = 8 ) P(NEQN). -! -! Workspace, real ( kind = 8 ) YP(NEQN), used to hold values of the -! solution derivative. -! -! Workspace, real ( kind = 8 ) YPOUT(NEQN), used to hold values of the -! solution derivative. -! -! Workspace, real ( kind = 8 ) PHI(NEQN,16), contains divided difference -! information about the polynomial interpolant to the solution. -! -! Workspace, real ( kind = 8 ) ALPHA(12), BETA(12), SIG(13). -! -! Workspace, real ( kind = 8 ) V(12), W(12), G(13). -! -! Input/output, logical PHASE1, indicates whether the program is in the -! first phase, when it always wants to increase the ODE method order. -! -! Workspace, real ( kind = 8 ) PSI(12), contains information about -! the polynomial interpolant to the solution. -! -! Input/output, real ( kind = 8 ) X, a "working copy" of T, the current value -! of the independent variable, which is adjusted as the code attempts -! to take a step. -! -! Input/output, real ( kind = 8 ) H, the current stepsize. -! -! Input/output, real ( kind = 8 ) HOLD, the last successful stepsize. -! -! Input/output, logical START, is TRUE on input for the first step. -! The program initializes data, and sets START to FALSE. -! -! Input/output, real ( kind = 8 ) TOLD, the previous value of T. -! -! Input/output, real ( kind = 8 ) DELSGN, the sign (+1 or -1) of -! TOUT - T. -! -! Input/output, integer ( kind = 4 ) NS, the number of steps taken with -! stepsize H. -! -! Input/output, logical NORND, ? -! -! Input, integer ( kind = 4 ) K, the order of the current ODE method. -! -! Input, integer ( kind = 4 ) KOLD, the order of the ODE method on the -! previous step. -! -! Input/output, integer ( kind = 4 ) ISNOLD, the previous value of ISN, the -! sign of IFLAG. -! -! Local parameters: -! -! Local, integer MAXNUM, the maximum number of steps allowed in one -! call to DE. -! - implicit none - - integer ( kind = 4 ) neqn - - real ( kind = 8 ) absdel - real ( kind = 8 ) abseps - real ( kind = 8 ) abserr - real ( kind = 8 ) alpha(12) - real ( kind = 8 ) beta(12) - logical crash - real ( kind = 8 ) del - real ( kind = 8 ) delsgn - real ( kind = 8 ) eps - external f - real ( kind = 8 ) fouru - real ( kind = 8 ) g(13) - real ( kind = 8 ) h - real ( kind = 8 ) hold - integer ( kind = 4 ) iflag - integer ( kind = 4 ) isn - integer ( kind = 4 ) isnold - integer ( kind = 4 ) k - integer ( kind = 4 ) kle4 - integer ( kind = 4 ) kold - ! integer ( kind = 4 ) l - integer ( kind = 4 ), parameter :: maxnum = 500 - logical nornd - integer ( kind = 4 ) nostep - integer ( kind = 4 ) ns - real ( kind = 8 ) p(neqn) - real ( kind = 8 ) phi(neqn,16) - logical phase1 - real ( kind = 8 ) psi(12) - real ( kind = 8 ) releps - real ( kind = 8 ) relerr - real ( kind = 8 ) sig(13) - logical start - logical stiff - real ( kind = 8 ) t - real ( kind = 8 ) tend - real ( kind = 8 ) told - real ( kind = 8 ) tout - real ( kind = 8 ) v(12) - real ( kind = 8 ) w(12) - real ( kind = 8 ) wt(neqn) - real ( kind = 8 ) x - real ( kind = 8 ) y(neqn) - real ( kind = 8 ) yp(neqn) - real ( kind = 8 ) ypout(neqn) - real ( kind = 8 ) yy(neqn) -! -! Test for improper parameters. -! - fouru = 4.0D+00 * epsilon ( fouru ) - - if ( neqn < 1 ) then - iflag = 6 - return - end if - - if ( t == tout ) then - iflag = 6 - return - end if - - if ( relerr < 0.0D+00 .or. abserr < 0.0D+00 ) then - iflag = 6 - print*, "relerr or abserr negative" - return - end if - - eps = max ( relerr, abserr ) - - if ( eps <= 0.0D+00 ) then - iflag = 6 - print*, "eps negative" - return - end if - - if ( iflag == 0 ) then - iflag = 6 - return - end if - - isn = sign ( 1, iflag ) - iflag = abs ( iflag ) - - if ( iflag /= 1 ) then - - if ( t /= told ) then - iflag = 6 - return - end if - - if ( iflag < 2 .or. 5 < iflag ) then - iflag = 6 - return - end if - end if -! -! On each call set interval of integration and counter for number of -! steps. Adjust input error tolerances to define weight vector for -! subroutine STEP. -! - del = tout - t - absdel = abs ( del ) - - if ( isn < 0 ) then - tend = tout - else - tend = t + 10.0D+00 * del - end if - - nostep = 0 - kle4 = 0 - stiff = .false. - releps = relerr / eps - abseps = abserr / eps -! -! On start and restart, also set work variables X and YY(*), store the -! direction of integration, and initialize the step size. -! - if ( iflag == 1 .or. isnold < 0 .or. delsgn * del <= 0.0D+00 ) then - - start = .true. - x = t - yy(1:neqn) = y(1:neqn) - delsgn = sign ( 1.0D+00, del ) - h = sign ( max ( abs ( tout - x ), fouru * abs ( x ) ), tout - x ) - - end if -! -! If already past the output point, then interpolate and return. -! - do - - if ( absdel <= abs ( x - t ) ) then - call intrp ( x, yy, tout, y, ypout, neqn, kold, phi, psi ) - iflag = 2 - t = tout - told = t - isnold = isn - exit - end if -! -! If we cannot go past the output point, and we are sufficiently -! close to it, then extrapolate and return. -! - if ( isn <= 0 .and. abs ( tout - x ) < fouru * abs ( x ) ) then - h = tout - x - call f ( x, yy, yp ) - y(1:neqn) = yy(1:neqn) + h * yp(1:neqn) - iflag = 2 - t = tout - told = t - isnold = isn - exit - end if -! -! Test for too many steps. -! - if ( maxnum <= nostep ) then - iflag = isn * 4 - if ( stiff ) then - iflag = isn * 5 - end if - y(1:neqn) = yy(1:neqn) - t = x - told = t - isnold = 1 - exit - end if -! -! Limit the step size, set the weight vector and take a step. -! - h = sign ( min ( abs ( h ), abs ( tend - x ) ), h ) - wt(1:neqn) = releps * abs ( yy(1:neqn) ) + abseps - - call step ( x, yy, f, neqn, h, eps, wt, start, & - hold, k, kold, crash, phi, p, yp, psi, & - alpha, beta, sig, v, w, g, phase1, ns, nornd ) -! -! Test for tolerances too small. -! - if ( crash ) then - iflag = isn * 3 - relerr = eps * releps - abserr = eps * abseps - y(1:neqn) = yy(1:neqn) - t = x - told = t - isnold = 1 - exit - end if -! -! Augment the step counter and test for stiffness. -! - nostep = nostep + 1 - kle4 = kle4 + 1 - - if ( 4 < kold ) then - kle4 = 0 - end if - - if ( 50 <= kle4 ) then - stiff = .true. - end if - - end do - - return -! SJP Issue #835 -! For Intel compliance add "end subroutine" -end subroutine -subroutine step ( x, y, f, neqn, h, eps, wt, start, hold, k, kold, crash, & - phi, p, yp, psi, alpha, beta, sig, v, w, g, phase1, ns, nornd ) - -!***************************************************************************** -! -!! STEP integrates the system of ODE's one step, from X to X+H. -! -! Discussion: -! -! This routine integrates a system of first order ordinary differential -! equations one step, normally from x to x+h, using a modified divided -! difference form of the Adams PECE formulas. Local extrapolation is -! used to improve absolute stability and accuracy. The code adjusts its -! order and step size to control the local error per unit step in a -! generalized sense. Special devices are included to control roundoff -! error and to detect when the user is requesting too much accuracy. -! -! STEP is normally not called directly by the user. However, it is -! possible to do so. -! -! On the first call to STEP, the user must pass in values for: -! * X, the initial value of the independent variable; -! * Y, the vector of initial values of dependent variables; -! * NEQN, the number of equations to be integrated; -! * H, the nominal step size indicating direction of integration -! and maximum size of step. H must be a variable, not a constant; -! * EPS, the local error tolerance per step. EPS must be variable; -! * WT, the vector of non-zero weights for error criterion; -! * START, set to TRUE. -! -! STEP requires the L2 norm of the vector with components -! local error(1:NEQN) / WT(1:NEQN) -! to be less than EPS for a successful step. The array WT allows the user -! to specify an error test appropriate for the problem. For example, -! if WT(L): -! = 1.0, specifies absolute error, -! = abs(Y(L)), specifies error relative to the most recent value of -! the L-th component of the solution, -! = abs(YP(L)), specifies error relative to the most recent value of -! the L-th component of the derivative, -! = max (WT(L),abs(Y(L))), specifies error relative to the largest -! magnitude of L-th component obtained so far, -! = abs(Y(L))*RELERR/EPS + ABSERR/EPS, specifies a mixed -! relative-absolute test where EPS = max ( RELERR, ABSERR ). -! -! On subsequent calls to STEP, the routine is designed so that all -! information needed to continue the integration, including the next step -! size H and the next order K, is returned with each step. With the -! exception of the step size, the error tolerance, and the weights, none -! of the parameters should be altered. The array WT must be updated after -! each step to maintain relative error tests like those above. -! -! Normally the integration is continued just beyond the desired endpoint -! and the solution interpolated there with subroutine INTRP. If it is -! impossible to integrate beyond the endpoint, the step size may be -! reduced to hit the endpoint since the code will not take a step -! larger than the H input. -! -! Changing the direction of integration, that is, the sign of H, requires -! the user to set START = TRUE before calling STEP again. This is the -! only situation in which START should be altered. -! -! A successful step: the subroutine returns after each successful step with -! START and CRASH both set to FALSE. X represents the independent variable -! advanced by one step of length HOLD from its input value; Y has been -! updated to the solution vector at the new value of X. All other parameters -! represent information corresponding to the new X needed to continue the -! integration. -! -! Unsuccessful steps: when the error tolerance is too small, the subroutine -! returns without taking a step and sets CRASH to TRUE. An appropriate step -! size and error tolerance for continuing are estimated and all other -! information is restored as upon input before returning. To continue -! with the larger tolerance, the user just calls the code again. A -! restart is neither required nor desirable. -! -! Licensing: -! -! This code is distributed under the GNU LGPL license. -! -! Modified: -! -! 02 August 2009 -! -! Author: -! -! Original FORTRAN77 version by Lawrence Shampine, Marilyn Gordon. -! FORTRAN90 version by John Burkardt. -! -! Reference: -! -! Lawrence Shampine, Marilyn Gordon, -! Computer Solution of Ordinary Differential Equations: -! The Initial Value Problem, -! Freeman, 1975, -! ISBN: 0716704617, -! LC: QA372.S416. -! -! Parameters: -! -! Input/output, real ( kind = 8 ) X, the value of the independent variable. -! -! Input/output, real ( kind = 8 ) Y(NEQN), the approximate solution at the -! current value of the independent variable. -! -! Input, external F, the name of a user-supplied routine of the form -! subroutine f ( t, y, yp ) -! real ( kind = 8 ) t -! real ( kind = 8 ) y(neqn) -! real ( kind = 8 ) yp(neqn) -! which accepts input values T and Y(1:NEQN), evaluates the right hand -! sides of the ODE, and stores the result in YP(1:NEQN). -! -! Input, integer ( kind = 4 ) NEQN, the number of equations. -! -! Input/output, real ( kind = 8 ) H, the suggested stepsize. -! -! Input/output, real ( kind = 8 ) EPS, the local error tolerance. -! -! Input, real ( kind = 8 ) WT(NEQN), the vector of error weights. -! -! Input/output, logical START, is set to TRUE before the first step. -! The program initializes data, and resets START to FALSE. -! -! Input/output, real ( kind = 8 ) HOLD, the step size used on the last -! successful step. -! -! Input/output, integer ( kind = 4 ) K, the appropriate order for the -! next step. -! -! Input/output, integer ( kind = 4 ) KOLD, the order used on the last -! successful step. -! -! Output, logical CRASH, is set to TRUE if no step can be taken. -! -! Workspace, real ( kind = 8 ) PHI(NEQN,16), contains divided difference -! information about the polynomial interpolant to the solution. -! -! Workspace, real ( kind = 8 ) P(NEQN). -! -! Workspace, real ( kind = 8 ) YP(NEQN), used to hold values of the -! solution derivative. -! -! Workspace, real ( kind = 8 ) PSI(12), contains information about -! the polynomial interpolant to the solution. -! -! Workspace, real ( kind = 8 ) ALPHA(12), BETA(12), SIG(13). -! -! Workspace, real ( kind = 8 ) V(12), W(12), G(13). -! -! Input/output, logical PHASE1, indicates whether the program is in the -! first phase, when it always wants to increase the ODE method order. -! -! Input/output, integer ( kind = 4 ) NS, the number of steps taken with -! stepsize H. -! -! Input/output, logical NORND, ? -! - implicit none - - integer ( kind = 4 ) neqn - - real ( kind = 8 ) absh - real ( kind = 8 ) alpha(12) - real ( kind = 8 ) beta(12) - logical crash - real ( kind = 8 ) eps - real ( kind = 8 ) erk - real ( kind = 8 ) erkm1 - real ( kind = 8 ) erkm2 - real ( kind = 8 ) erkp1 - real ( kind = 8 ) err - external f - real ( kind = 8 ) fouru - real ( kind = 8 ) g(13) - real ( kind = 8 ), dimension ( 13 ), parameter :: gstr = (/ & - 0.50D+00, 0.0833D+00, 0.0417D+00, 0.0264D+00, 0.0188D+00, & - 0.0143D+00, 0.0114D+00, 0.00936D+00, 0.00789D+00, 0.00679D+00, & - 0.00592D+00, 0.00524D+00, 0.00468D+00 /) - real ( kind = 8 ) h - real ( kind = 8 ) hnew - real ( kind = 8 ) hold - integer ( kind = 4 ) i - integer ( kind = 4 ) ifail - integer ( kind = 4 ) iq - integer ( kind = 4 ) j - integer ( kind = 4 ) k - integer ( kind = 4 ) km1 - integer ( kind = 4 ) km2 - integer ( kind = 4 ) knew - integer ( kind = 4 ) kold - integer ( kind = 4 ) kp1 - integer ( kind = 4 ) kp2 - integer ( kind = 4 ) l - logical nornd - integer ( kind = 4 ) ns - integer ( kind = 4 ) nsp1 - real ( kind = 8 ) p(neqn) - real ( kind = 8 ) p5eps - logical phase1 - real ( kind = 8 ) phi(neqn,16) - real ( kind = 8 ) psi(12) - real ( kind = 8 ) r - real ( kind = 8 ) rho - real ( kind = 8 ) round - real ( kind = 8 ) sig(13) - logical start - real ( kind = 8 ) total - real ( kind = 8 ) tau - real ( kind = 8 ) temp1 - real ( kind = 8 ) temp2 - real ( kind = 8 ), dimension ( 13 ), parameter :: two = (/ & - 2.0D+00, 4.0D+00, 8.0D+00, 16.0D+00, 32.0D+00, & - 64.0D+00, 128.0D+00, 256.0D+00, 512.0D+00, 1024.0D+00, & - 2048.0D+00, 4096.0D+00, 8192.0D+00/) - real ( kind = 8 ) twou - real ( kind = 8 ) v(12) - real ( kind = 8 ) w(12) - real ( kind = 8 ) wt(neqn) - real ( kind = 8 ) x - real ( kind = 8 ) xold - real ( kind = 8 ) y(neqn) - real ( kind = 8 ) yp(neqn) - - twou = 2.0D+00 * epsilon ( twou ) - fouru = 2.0D+00 * twou -! -! Check if the step size or error tolerance is too small. If this is the -! first step, initialize the PHI array and estimate a starting step size. -! -! If the step size is too small, determine an acceptable one. -! - crash = .true. - - if ( abs ( h ) < fouru * abs ( x ) ) then - h = sign ( fouru * abs ( x ), h ) - return - end if - - p5eps = 0.5D+00 * eps -! -! If the error tolerance is too small, increase it to an acceptable value. -! - round = twou * sqrt ( sum ( ( y(1:neqn) / wt(1:neqn) )**2 ) ) - - if ( p5eps < round ) then - eps = 2.0D+00 * round * ( 1.0D+00 + fouru ) - return - end if - - crash = .false. - g(1) = 1.0D+00 - g(2) = 0.5D+00 - sig(1) = 1.0D+00 -! -! Initialize. Compute an appropriate step size for the first step. -! - if ( start ) then - - call f ( x, y, yp ) - phi(1:neqn,1) = yp(1:neqn) - phi(1:neqn,2) = 0.0D+00 - total = sqrt ( sum ( ( yp(1:neqn) / wt(1:neqn) )**2 ) ) - absh = abs ( h ) - if ( eps < 16.0D+00 * total * h * h ) then - absh = 0.25D+00 * sqrt ( eps / total ) - end if - h = sign ( max ( absh, fouru * abs ( x ) ), h ) - hold = 0.0D+00 - k = 1 - kold = 0 - start = .false. - phase1 = .true. - nornd = .true. - - if ( p5eps <= 100.0D+00 * round ) then - nornd = .false. - phi(1:neqn,15) = 0.0D+00 - end if - - end if - - ifail = 0 -! -! Compute coefficients of formulas for this step. Avoid computing -! those quantities not changed when step size is not changed. -! - do - - kp1 = k + 1 - kp2 = k + 2 - km1 = k - 1 - km2 = k - 2 -! -! NS is the number of steps taken with size H, including the current -! one. When K < NS, no coefficients change. -! - if ( h /= hold ) then - ns = 0 - end if - - if ( ns <= kold ) then - ns = ns + 1 - end if - - nsp1 = ns + 1 -! -! Compute those components of ALPHA, BETA, PSI and SIG which change. -! - if ( ns <= k ) then - - beta(ns) = 1.0D+00 - alpha(ns) = 1.0D+00 / real ( ns, kind = 8 ) - temp1 = h * real ( ns, kind = 8 ) - sig(nsp1) = 1.0D+00 - - do i = nsp1, k - temp2 = psi(i-1) - psi(i-1) = temp1 - beta(i) = beta(i-1) * psi(i-1) / temp2 - temp1 = temp2 + h - alpha(i) = h / temp1 - sig(i+1) = real ( i, kind = 8 ) * alpha(i) * sig(i) - end do - - psi(k) = temp1 -! -! Compute coefficients G. -! -! Initialize V and set W. -! - if ( ns <= 1 ) then - - do iq = 1, k - v(iq) = 1.0D+00 / real ( iq * ( iq + 1 ), kind = 8 ) - w(iq) = v(iq) - end do -! -! If order was raised, update the diagonal part of V. -! - else - - if ( kold < k ) then - - v(k) = 1.0D+00 / real ( k * kp1, kind = 8 ) - - do j = 1, ns - 2 - i = k - j - v(i) = v(i) - alpha(j+1) * v(i+1) - end do - - end if -! -! Update V and set W. -! - do iq = 1, kp1 - ns - v(iq) = v(iq) - alpha(ns) * v(iq+1) - w(iq) = v(iq) - end do - g(nsp1) = w(1) - - end if -! -! Compute the G in the work vector W. -! - do i = ns + 2, kp1 - do iq = 1, kp2 - i - w(iq) = w(iq) - alpha(i-1) * w(iq+1) - end do - g(i) = w(1) - end do - - end if -! -! Predict a solution P, evaluate derivatives using predicted -! solution, estimate local error at order K and errors at orders K, -! K-1, K-2 as if a constant step size were used. -! -! Change PHI to PHI star. -! - do i = nsp1, k - phi(1:neqn,i) = beta(i) * phi(1:neqn,i) - end do -! -! Predict solution and differences. -! - phi(1:neqn,kp2) = phi(1:neqn,kp1) - phi(1:neqn,kp1) = 0.0D+00 - p(1:neqn) = 0.0D+00 - - do j = 1, k - i = kp1 - j - do l = 1, neqn - p(l) = p(l) + phi(l,i) * g(i) - phi(l,i) = phi(l,i) + phi(l,i+1) - end do - end do - - if ( .not. nornd ) then - do l = 1, neqn - tau = h * p(l) - phi(l,15) - p(l) = y(l) + tau - phi(l,16) = ( p(l) - y(l) ) - tau - end do - else - p(1:neqn) = y(1:neqn) + h * p(1:neqn) - end if - - xold = x - x = x + h - absh = abs ( h ) - call f ( x, p, yp ) -! -! Estimate the errors at orders K, K-1 and K-2. -! - erkm2 = 0.0D+00 - erkm1 = 0.0D+00 - erk = 0.0D+00 - - do l = 1, neqn - - if ( 0 < km2 ) then - erkm2 = erkm2 + ( ( phi(l,km1) + yp(l) - phi(l,1) ) / wt(l) )**2 - end if - - if ( 0 <= km2 ) then - erkm1 = erkm1 + ( ( phi(l,k) + yp(l) - phi(l,1) ) / wt(l) )**2 - end if - - erk = erk + ( ( yp(l) - phi(l,1) ) / wt(l) )**2 - - end do - - if ( 0 < km2 ) then - erkm2 = absh * sig(km1) * gstr(km2) * sqrt ( erkm2 ) - end if - - if ( 0 <= km2 ) then - erkm1 = absh * sig(k) * gstr(km1) * sqrt ( erkm1 ) - end if - - err = absh * sqrt ( erk ) * ( g(k) - g(kp1) ) - erk = absh * sqrt ( erk ) * sig(kp1) * gstr(k) - knew = k -! -! Test if the order should be lowered. -! - if ( 0 < km2 ) then - - if ( max ( erkm1, erkm2 ) <= erk ) then - knew = km1 - end if - - else if ( 0 == km2 ) then - - if ( erkm1 <= 0.5D+00 * erk ) then - knew = km1 - end if - - end if -! -! Test if the step was successful. -! - if ( err <= eps ) then - exit - end if -! -! The step is unsuccessful. Restore X, PHI and PSI. -! If third consecutive failure, set order to one. If the step fails more -! than three times, consider an optimal step size. Double the error -! tolerance and return if the estimated step size is too small for machine -! precision. -! -! Restore X, PHI and PSI. -! - phase1 = .false. - x = xold - do i = 1, k - phi(1:neqn,i) = ( phi(1:neqn,i) - phi(1:neqn,i+1) ) / beta(i) - end do - - do i = 2, k - psi(i-1) = psi(i) - h - end do -! -! On third failure, set the order to one. Thereafter, use optimal step size. -! - ifail = ifail + 1 - temp2 = 0.5D+00 - - if ( 3 < ifail ) then - - if ( p5eps < 0.25D+00 * erk ) then - temp2 = sqrt ( p5eps / erk ) - end if - end if - - if ( 3 <= ifail ) then - knew = 1 - end if - - h = temp2 * h - k = knew - - if ( abs ( h ) < fouru * abs ( x ) ) then - crash = .true. - h = sign ( fouru * abs ( x ), h ) - eps = eps + eps - return - end if - - end do -! -! The step is successful. Correct the predicted solution, evaluate -! the derivatives using the corrected solution and update the -! differences. Determine best order and step size for next step. -! - kold = k - hold = h -! -! Correct and evaluate. -! - if ( .not. nornd ) then - do l = 1, neqn - rho = h * g(kp1) * ( yp(l) - phi(l,1) ) - phi(l,16) - y(l) = p(l) + rho - phi(l,15) = ( y(l) - p(l) ) - rho - end do - else - y(1:neqn) = p(1:neqn) + h * g(kp1) * ( yp(1:neqn) - phi(1:neqn,1) ) - end if - - call f ( x, y, yp ) -! -! Update differences for the next step. -! - phi(1:neqn,kp1) = yp(1:neqn) - phi(1:neqn,1) - phi(1:neqn,kp2) = phi(1:neqn,kp1) - phi(1:neqn,kp2) - - do i = 1, k - phi(1:neqn,i) = phi(1:neqn,i) + phi(1:neqn,kp1) - end do -! -! Estimate error at order K+1 unless: -! * in first phase when always raise order, or, -! * already decided to lower order, or, -! * step size not constant so estimate unreliable. -! - erkp1 = 0.0D+00 - - if ( knew == km1 .or. k == 12 ) then - phase1 = .false. - end if - - if ( phase1 ) then - - k = kp1 - erk = erkp1 - - else if ( knew == km1 ) then - - k = km1 - erk = erkm1 - - else if ( kp1 <= ns ) then - - do l = 1, neqn - erkp1 = erkp1 + ( phi(l,kp2) / wt(l) )**2 - end do - erkp1 = absh * gstr(kp1) * sqrt ( erkp1 ) -! -! Using estimated error at order K+1, determine appropriate order -! for next step. -! - if ( k == 1 ) then - - if ( erkp1 < 0.5D+00 * erk ) then - k = kp1 - erk = erkp1 - end if - - else if ( erkm1 <= min ( erk, erkp1 ) ) then - - k = km1 - erk = erkm1 - - else if ( erkp1 < erk .and. k < 12 ) then - - k = kp1 - erk = erkp1 - - end if - - end if -! -! With the new order, determine appropriate step size for next step. -! - hnew = h + h - - if ( .not. phase1 ) then - - if ( p5eps < erk * two(k+1) ) then - - hnew = h - - if ( p5eps < erk ) then - temp2 = real ( k + 1, kind = 8 ) - r = ( p5eps / erk )**( 1.0D+00 / temp2 ) - hnew = absh * max ( 0.5D+00, min ( 0.9D+00, r ) ) - hnew = sign ( max ( hnew, fouru * abs ( x ) ), h ) - end if - - end if - - end if - - h = hnew - - return -! SJP Issue #835 -! For Intel compliance add "end subroutine" -end subroutine -subroutine intrp ( x, y, xout, yout, ypout, neqn, kold, phi, psi ) - -!***************************************************************************** -! -!! INTRP approximates the solution at XOUT by polynomial interpolation. -! -! Discussion: -! -! The methods in STEP approximate the solution near X by a polynomial. -! This routine approximates the solution at XOUT by evaluating the -! polynomial there. Information defining this polynomial is passed -! from STEP, so INTRP cannot be used alone. -! -! Licensing: -! -! This code is distributed under the GNU LGPL license. -! -! Modified: -! -! 02 August 2009 -! -! Author: -! -! Original FORTRAN77 version by Lawrence Shampine, Marilyn Gordon. -! FORTRAN90 version by John Burkardt. -! -! Reference: -! -! Lawrence Shampine, Marilyn Gordon, -! Computer Solution of Ordinary Differential Equations: -! The Initial Value Problem, -! Freeman, 1975, -! ISBN: 0716704617, -! LC: QA372.S416. -! -! Parameters: -! -! Input, real ( kind = 8 ) X, the point where the solution has been computed. -! -! Input, real ( kind = 8 ) Y(NEQN), the computed solution at X. -! -! Input, real ( kind = 8 ) XOUT, the point at which the solution is desired. -! -! Output, real ( kind = 8 ) YOUT(NEQN), the solution at XOUT. -! -! Output, real ( kind = 8 ) YPOUT(NEQN), the derivative of the solution -! at XOUT. -! -! Input, integer ( kind = 4 ) NEQN, the number of equations. -! -! Input, integer ( kind = 4 ) KOLD, the order used for the last -! successful step. -! -! Input, real ( kind = 8 ) PHI(NEQN,16), contains information about the -! interpolating polynomial. -! -! Input, real ( kind = 8 ) PSI(12), contains information about the -! interpolating polynomial. -! - implicit none - - integer ( kind = 4 ) neqn - - real ( kind = 8 ) eta - real ( kind = 8 ) g(13) - real ( kind = 8 ) gamma - real ( kind = 8 ) hi - integer ( kind = 4 ) i - integer ( kind = 4 ) j - integer ( kind = 4 ) ki - integer ( kind = 4 ) kold - ! integer ( kind = 4 ) l - real ( kind = 8 ) phi(neqn,16) - real ( kind = 8 ) psi(12) - real ( kind = 8 ) psijm1 - real ( kind = 8 ) rho(13) - real ( kind = 8 ) term - real ( kind = 8 ) w(13) - real ( kind = 8 ) x - real ( kind = 8 ) xout - real ( kind = 8 ) y(neqn) - real ( kind = 8 ) yout(neqn) - real ( kind = 8 ) ypout(neqn) - - hi = xout - x - ki = kold + 1 -! -! Initialize W for computing G. -! - do i = 1, ki - w(i) = 1.0D+00 / real ( i, kind = 8 ) - end do -! -! Compute G. -! - g(1) = 1.0D+00 - rho(1) = 1.0D+00 - term = 0.0D+00 - - do j = 2, ki - psijm1 = psi(j-1) - gamma = ( hi + term ) / psijm1 - eta = hi / psijm1 - do i = 1, ki + 1 - j - w(i) = gamma * w(i) - eta * w(i+1) - end do - g(j) = w(1) - rho(j) = gamma * rho(j-1) - term = psijm1 - end do -! -! Interpolate. -! - ypout(1:neqn) = 0.0D+00 - yout(1:neqn) = 0.0D+00 - - do j = 1, ki - i = ki + 1 - j - yout(1:neqn) = yout(1:neqn) + g(i) * phi(1:neqn,i) - ypout(1:neqn) = ypout(1:neqn) + rho(i) * phi(1:neqn,i) - end do - - yout(1:neqn) = y(1:neqn) + hi * yout(1:neqn) - - return -! SJP Issue #835 -! For Intel compliance add "end subroutine" -end subroutine -subroutine timestamp ( ) - -!***************************************************************************** -! -!! TIMESTAMP prints the current YMDHMS date as a time stamp. -! -! Example: -! -! 31 May 2001 9:45:54.872 AM -! -! Licensing: -! -! This code is distributed under the GNU LGPL license. -! -! Modified: -! -! 18 May 2013 -! -! Author: -! -! John Burkardt -! -! Parameters: -! -! None -! - implicit none - - character ( len = 8 ) ampm - integer ( kind = 4 ) d - integer ( kind = 4 ) h - integer ( kind = 4 ) m - integer ( kind = 4 ) mm - character ( len = 9 ), parameter, dimension(12) :: month = (/ & - 'January ', 'February ', 'March ', 'April ', & - 'May ', 'June ', 'July ', 'August ', & - 'September', 'October ', 'November ', 'December ' /) - integer ( kind = 4 ) n - integer ( kind = 4 ) s - integer ( kind = 4 ) values(8) - integer ( kind = 4 ) y - - call date_and_time ( values = values ) - - y = values(1) - m = values(2) - d = values(3) - h = values(5) - n = values(6) - s = values(7) - mm = values(8) - - if ( h < 12 ) then - ampm = 'AM' - else if ( h == 12 ) then - if ( n == 0 .and. s == 0 ) then - ampm = 'Noon' - else - ampm = 'PM' - end if - else - h = h - 12 - if ( h < 12 ) then - ampm = 'PM' - else if ( h == 12 ) then - if ( n == 0 .and. s == 0 ) then - ampm = 'Midnight' - else - ampm = 'AM' - end if - end if - end if - - write ( *, '(i2.2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & - d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) - - return -! SJP Issue #835 -! For Intel compliance add "end subroutine" -end subroutine - -end module ode_mod diff --git a/source/fortran/plot_radiation.f90 b/source/fortran/plot_radiation.f90 deleted file mode 100644 index 5925fe92..00000000 --- a/source/fortran/plot_radiation.f90 +++ /dev/null @@ -1,131 +0,0 @@ -module plot_radiation -#ifndef dp - use, intrinsic :: iso_fortran_env, only: dp=>real64 -#endif -implicit none - -contains - subroutine plot_Lz() - !! Write loss data to file for plotting - !! author: M Kovari, CCFE, Culham Science Centre - !! Write loss data to file for plotting - !! Compare to Figure 3 in Kallenbach 2016. - !! - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use impurity_radiation_module, only: nimp, imp_label - use read_radiation, only: read_lz - implicit none - - ! Subroutine declarations ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!! - - character(len=2) :: element - - real(dp) :: dummy - - integer :: i, j - - integer, parameter :: points = 27 - - ! Temperature plot points - real(dp), parameter :: te(points) = (/1., 1.2, 1.5, 2., 2.5, 3., 4., 5., 6., 7., 8., 9., & - 10., 12., 14., 16., 20., 30., 40., 50., 60., 70., 80., 90., 100., 150., 200./) - - real(dp) :: Lz_plot(nimp) - - real(dp), parameter :: netau = 0.5 - - open(unit=12,file='radiative_loss_functions.txt',status='replace') - write(12,'(30a11)')'Te (eV)', (imp_label(i), i=2,nimp) - - ! Just read data. Exclude hydrogen by starting at 2 - do i=2,nimp - element=imp_label(i) - dummy=read_lz(element,30.0d0,netau, mean_z=.false., mean_qz=.false., verbose=.false.) - enddo - - do i=1,points - do j=2,nimp - Lz_plot(j)=read_lz(imp_label(j),te(i),netau, mean_z=.false., mean_qz=.false., verbose=.false.) - enddo - write(12,'(30es11.3)')te(i), (Lz_plot(j), j=2,nimp) - enddo - close(unit=12) - - end subroutine plot_Lz - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine plot_z() - !! Write z and z^2 data to file for plotting - !! author: M Kovari, CCFE, Culham Science Centre - !! Write z and z^2 data to file for plotting - !! - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use impurity_radiation_module, only: nimp, imp_label - use read_radiation, only: read_lz - implicit none - - ! Subroutine declarations ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!! - - character(len=2) :: element - - real(dp) :: dummy - - integer :: i, j, k - - integer, parameter :: points = 27 - - real(dp), parameter :: te(points) = (/1., 1.2, 1.5, 2., 2.5, 3.,4., 5., 6., 7., & - 8., 9.,10.,12., 14., 16., 20., 30., 40., 50., 60., 70., 80., 90., 100., 150., 200./) - - real(dp) :: Z_plot(3,nimp) - - real(dp), parameter :: netau(3) = (/0.1, 1.0, 10.0/) - - open(unit=12,file='mean_Z.txt',status='replace') - - write(12,*)'Mean Z' - - write(12,'(a11, 42(3(a4, es7.1)))')'Te (ev)',((imp_label(i),netau(j), j=1,3),i=2,nimp) - - ! Just read data. Exclude hydrogen by starting at 2 - do i = 2, nimp - element = imp_label(i) - dummy = read_lz(element,30.0d0,1.0d0, mean_z=.true., mean_qz=.false., verbose=.true.) - dummy = read_lz(element,30.0d0,1.0d0, mean_z=.false., mean_qz=.true., verbose=.true.) - enddo - - do i = 1, points - do j = 1, 3 - do k = 2, nimp - element = imp_label(k) - Z_plot(j,k) = read_lz(element,te(i),netau(j), mean_z=.true., mean_qz=.false., verbose=.false.) - enddo - enddo - write(12,'(42es11.3)')te(i), ((Z_plot(j,k), j=1,3), k=2,nimp) - enddo - - write(12,*) - write(12,*)'Mean Z^2' - write(12,'(a11, 42(3(a4, es7.1)))')'Te (ev)', ((imp_label(i),netau(j), j=1,3),i=2,nimp) - do i = 1, points - do j = 1, 3 - do k = 2, nimp - element = imp_label(k) - Z_plot(j,k) = read_lz(element,te(i),netau(j), mean_z=.false., mean_qz=.true., verbose=.false.) - enddo - enddo - write(12,'(42es11.3)')te(i), ((Z_plot(j,k), j=1,3), k=2,nimp) - enddo - - close(unit=12) - - end subroutine plot_z - -end module plot_radiation diff --git a/source/fortran/precision_mod.f90 b/source/fortran/precision_mod.f90 deleted file mode 100755 index 15752494..00000000 --- a/source/fortran/precision_mod.f90 +++ /dev/null @@ -1,10 +0,0 @@ -module precision_mod - - integer, parameter :: real_8_ = selected_real_kind(12,100) - integer, parameter :: real_4_ = selected_real_kind(6,36) - - ! Precision can be adjusted below: - - integer, parameter :: p_ = real_8_ ! real_8_ - -end module precision_mod diff --git a/source/fortran/quanc8.f90 b/source/fortran/quanc8.f90 deleted file mode 100755 index 94493364..00000000 --- a/source/fortran/quanc8.f90 +++ /dev/null @@ -1,189 +0,0 @@ - subroutine quanc8(fun,a,b,abserr,relerr,result,errest,nofun,flag) -! -! double precision fun, a, b, abserr, relerr, result, errest, flag - double precision a, b, abserr, relerr, result, errest, flag - integer nofun - external fun -! -! estimate the integral of fun(x) from a to b -! to a user provided tolerance. -! an automatic adaptive routine based on -! the 8-panel newton-cotes rule. -! -! input .. -! -! fun the name of the integrand function subprogram fun(x). -! a the lower limit of integration. -! b the upper limit of integration.(b may be less than a.) -! relerr a relative error tolerance. (should be non-negative) -! abserr an absolute error tolerance. (should be non-negative) -! -! output .. -! -! result an approximation to the integral hopefully satisfying the -! least stringent of the two error tolerances. -! errest an estimate of the magnitude of the actual error. -! nofun the number of function values used in calculation of result. -! flag a reliability indicator. if flag is zero, then result -! probably satisfies the error tolerance. if flag is -! xxx.yyy , then xxx = the number of intervals which have -! not converged and 0.yyy = the fraction of the interval -! left to do when the limit on nofun was approached. -! - double precision w0,w1,w2,w3,w4,area,x0,f0,stone,step,cor11,temp - double precision qprev,qnow,qdiff,qleft,esterr,tolerr - double precision qright(31),f(16),x(16),fsave(8,30),xsave(8,30) - double precision dabs,dmax1 - integer levmin,levmax,levout,nomax,nofin,lev,nim,i,j -! -! *** stage 1 *** general initialization -! set constants. -! - levmin = 1 - levmax = 30 - levout = 6 - nomax = 5000 - nofin = nomax - 8*(levmax-levout+2**(levout+1)) -! -! trouble when nofun reaches nofin -! - w0 = 3956.0d0 / 14175.0d0 - w1 = 23552.0d0 / 14175.0d0 - w2 = -3712.0d0 / 14175.0d0 - w3 = 41984.0d0 / 14175.0d0 - w4 = -18160.0d0 / 14175.0d0 -! -! initialize running sums to zero. -! - flag = 0.0d0 - result = 0.0d0 - cor11 = 0.0d0 - errest = 0.0d0 - area = 0.0d0 - nofun = 0 - if (a .eq. b) return -! -! *** stage 2 *** initialization for first interval -! - lev = 0 - nim = 1 - x0 = a - x(16) = b - qprev = 0.0d0 - f0 = fun(x0) - stone = (b - a) / 16.0d0 - x(8) = (x0 + x(16)) / 2.0d0 - x(4) = (x0 + x(8)) / 2.0d0 - x(12) = (x(8) + x(16)) / 2.0d0 - x(2) = (x0 + x(4)) / 2.0d0 - x(6) = (x(4) + x(8)) / 2.0d0 - x(10) = (x(8) + x(12)) / 2.0d0 - x(14) = (x(12) + x(16)) / 2.0d0 - do 25 j = 2, 16, 2 - f(j) = fun(x(j)) - 25 continue - nofun = 9 -! -! *** stage 3 *** central calculation -! requires qprev,x0,x2,x4,...,x16,f0,f2,f4,...,f16. -! calculates x1,x3,...x15, f1,f3,...f15,qleft,qright,qnow,qdiff,area. -! - 30 x(1) = (x0 + x(2)) / 2.0d0 - f(1) = fun(x(1)) - do 35 j = 3, 15, 2 - x(j) = (x(j-1) + x(j+1)) / 2.0d0 - f(j) = fun(x(j)) - 35 continue - nofun = nofun + 8 - step = (x(16) - x0) / 16.0d0 - qleft = (w0*(f0 + f(8)) + w1*(f(1)+f(7)) + w2*(f(2)+f(6)) & - + w3*(f(3)+f(5)) + w4*f(4)) * step - qright(lev+1)=(w0*(f(8)+f(16))+w1*(f(9)+f(15))+w2*(f(10)+f(14)) & - + w3*(f(11)+f(13)) + w4*f(12)) * step - qnow = qleft + qright(lev+1) - qdiff = qnow - qprev - area = area + qdiff -! -! *** stage 4 *** interval convergence test -! - esterr = dabs(qdiff) / 1023.0d0 - tolerr = dmax1(abserr,relerr*dabs(area)) * (step/stone) - if (lev .lt. levmin) go to 50 - if (lev .ge. levmax) go to 62 - if (nofun .gt. nofin) go to 60 - if (esterr .le. tolerr) go to 70 -! -! *** stage 5 *** no convergence -! locate next interval. -! - 50 nim = 2*nim - lev = lev+1 -! -! store right hand elements for future use. -! - do 52 i = 1, 8 - fsave(i,lev) = f(i+8) - xsave(i,lev) = x(i+8) - 52 continue -! -! assemble left hand elements for immediate use. -! - qprev = qleft - do 55 i = 1, 8 - j = -i - f(2*j+18) = f(j+9) - x(2*j+18) = x(j+9) - 55 continue - go to 30 -! -! *** stage 6 *** trouble section -! number of function values is about to exceed limit. -! - 60 nofin = 2*nofin - levmax = levout - flag = flag + (b - x0) / (b - a) - go to 70 -! -! current level is levmax. -! - 62 flag = flag + 1.0d0 -! -! *** stage 7 *** interval converged -! add contributions into running sums. -! - 70 result = result + qnow - errest = errest + esterr - cor11 = cor11 + qdiff / 1023.0d0 -! -! locate next interval. -! - 72 if (nim .eq. 2*(nim/2)) go to 75 - nim = nim/2 - lev = lev-1 - go to 72 - 75 nim = nim + 1 - if (lev .le. 0) go to 80 -! -! assemble elements required for the next interval. -! - qprev = qright(lev) - x0 = x(16) - f0 = f(16) - do 78 i = 1, 8 - f(2*i) = fsave(i,lev) - x(2*i) = xsave(i,lev) - 78 continue - go to 30 -! -! *** stage 8 *** finalize and return -! - 80 result = result + cor11 -! -! make sure errest not less than roundoff level. -! - if (errest .eq. 0.0d0) return - 82 temp = dabs(result) + errest - if (temp .ne. dabs(result)) return - errest = 2.0d0*errest - go to 82 - end diff --git a/source/fortran/read_and_get_atomic_data.f90 b/source/fortran/read_and_get_atomic_data.f90 deleted file mode 100755 index 4957c3b6..00000000 --- a/source/fortran/read_and_get_atomic_data.f90 +++ /dev/null @@ -1,319 +0,0 @@ - -module read_and_get_atomic_data - !! Module for reading atomic data - !! author: M Kovari, CCFE, Culham Science Centre - !! N/A - !! This module reads atomic data for the PROCESS Kallenbach divertor model - !! - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Modules to import ! - ! !!!!!!!!!!!!!!!!!!!! - -#ifndef dp - use, intrinsic :: iso_fortran_env, only: dp=>real64 -#endif - implicit none - - ! Var for subroutine get_h_rates requiring re-initialisation on each new run - logical, save :: FirstCall - -contains - character(len=300) function hdatadir() - implicit none - character(len=200) :: process_dir - CALL get_environment_variable("PYTHON_PROCESS_ROOT", process_dir) - if (process_dir == "") then - hdatadir = INSTALLDIR//'/process/data/h_data/' - else - hdatadir = trim(process_dir)//'/data/h_data/' - end if - - end function hdatadir - - subroutine init_read_and_get_atomic_data - !! Initialise module variables - implicit none - - FirstCall = .true. - end subroutine init_read_and_get_atomic_data - - subroutine get_h_rates(density, temperature, s, al, Rcx, plt, prb, mass, verbose) - !! - !! author: M Kovari, CCFE, Culham Science Centre - !! density : input real : electron density (m^-3) - !! temperature : input real : electron temperature (eV) - !! s : output real : ionisation rate coefficient (m^3/s) - !! al : output real : recombination rate coefficient (m^3/s) - !! Rcx : output real : charge exchange rate coefficient (m^3/s) - !! plt : output real : line radiation power rate coefficient (Wm^3) - !! prb : output real : continuum radiation power rate coefficient (Wm^3) - !! mass : input real : relative atomic mass for CX rate coefficient - !! verbose : input logical : verbose switch - !! - !! - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use maths_library, only: interpolate - implicit none - - ! Subroutine declarations ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!! - - real(dp), intent(in) :: density, temperature, mass - - real(dp), intent(out) :: s, al, Rcx, plt, prb - - logical, intent(in) :: verbose - - integer :: m - - ! These arrays are read only once and then saved - ! Density values: "_d" - real(dp), save, dimension(24) :: scd_d, acd_d, ccd_d, plt_d,prb_d - - ! Temperature values: "_t" - real(dp), save, dimension(29) :: scd_t,acd_t,ccd_t, plt_t, prb_t - - ! Rate coefficients: "_r" - real(dp), save, dimension(24,29) :: scd_r, acd_r, ccd_r, plt_r,prb_r - - character(len=200) :: acd_file, scd_file, plt_file, prb_file, ccd_file - - real(dp) :: logdens, logtemp - - logical :: iexist - - integer :: ine, ite - - ! Maxima for log density and log temperature in each data file - real(dp), save :: max_scd_d, max_scd_t - real(dp), save :: max_acd_d, max_acd_t - real(dp), save :: max_ccd_d, max_ccd_t - real(dp), save :: max_plt_d, max_plt_t - real(dp), save :: max_prb_d, max_prb_t - - ine = 24 - ite = 29 - - logdens = log10(density/1.0D6) - logtemp = log10(temperature) - - ! Read the data - if(FirstCall) then - - FirstCall=.false. - - ! Add trailing / to hdatadir() if necessary - ! if (index(hdatadir(),'/',.true.) .ne. len(trim(hdatadir()))) hdatadir() = hdatadir()//'/' - ! if (index(hdatadir(),'\',.true.) .ne. len(trim(hdatadir()))) hdatadir() = hdatadir()//'\' - - acd_file = trim(hdatadir())//'acd12_h.dat' - scd_file = trim(hdatadir())//'scd12_h.dat' - plt_file = trim(hdatadir())//'plt12_h.dat' - prb_file = trim(hdatadir())//'prb12_h.dat' - - m = floor(mass+0.5) ! round to an integer - ! Select the correct atomic species for the CX rates: H, D or T - if (m == 1) then - ccd_file = trim(hdatadir())//'ccd96_h.dat' - elseif (m == 2) then - ccd_file = trim(hdatadir())//'ccd96_d.dat' - elseif (m == 3) then - ccd_file = trim(hdatadir())//'ccd96_t.dat' - else - write(*,*) 'The atomic mass is ', m, '. It must be in the range 1-3.' - end if - - inquire(file=trim(acd_file), exist=iexist) - if (.not.iexist) write(*,*) "ERROR File "// acd_file // ' does not exist' - inquire(file=trim(scd_file), exist=iexist) - if (.not.iexist) write(*,*) "ERROR File "// scd_file // ' does not exist' - inquire(file=trim(plt_file), exist=iexist) - if (.not.iexist) write(*,*) "ERROR File "// plt_file // ' does not exist' - inquire(file=trim(prb_file), exist=iexist) - if (.not.iexist) write(*,*) "ERROR File "// prb_file // ' does not exist' - inquire(file=trim(ccd_file), exist=iexist) - if (.not.iexist) write(*,*) "ERROR File "// ccd_file // ' does not exist' - - ! ionisation data - call read_atomdat(scd_d,scd_t,scd_r, ine=24, ite=29, filename=scd_file, verbose=verbose) - ! recombination data - call read_atomdat(acd_d,acd_t,acd_r, ine=24, ite=29, filename=acd_file, verbose=verbose) - ! CX data - call read_atomdat(ccd_d,ccd_t,ccd_r, ine=24, ite=29, filename=ccd_file, verbose=verbose) - ! line radiation data - call read_atomdat(plt_d,plt_t,plt_r, ine=24, ite=29, filename=plt_file, verbose=verbose) - ! continuum radiation data - call read_atomdat(prb_d,prb_t,prb_r, ine=24, ite=29, filename=prb_file, verbose=verbose) - - ! Store the maxima for log density and log temperature in each data file - ! Subtract a smidgen to ensure that the values submitted for interpolation are - ! the range. - max_scd_d=maxval(scd_d) - 0.00001 - max_scd_t=maxval(scd_t) - 0.00001 - - max_acd_d=maxval(acd_d) - 0.00001 - max_acd_t=maxval(acd_t) - 0.00001 - - max_ccd_d=maxval(ccd_d) - 0.00001 - max_ccd_t=maxval(ccd_t) - 0.00001 - - max_plt_d=maxval(plt_d) - 0.00001 - max_plt_t=maxval(plt_t) - 0.00001 - - max_prb_d=maxval(prb_d) - 0.00001 - max_prb_t=maxval(prb_t) - 0.00001 - - end if - - ! Using function interpolate(x_len, x_array, y_len, y_array, f, x, y, delta) - ! ionisation - s = 1.d-6*10.0**(interpolate(ine, scd_d, ite, scd_t, scd_r, & - min(logdens,max_scd_d), & - min(logtemp,max_scd_t))) - ! recombination - al = 1.d-6*10.0**(interpolate(ine, acd_d, ite, acd_t, acd_r, & - min(logdens,max_acd_d), & - min(logtemp,max_acd_t))) - ! Rcx - Rcx = 1.d-6*10.0**(interpolate(ine, ccd_d, ite, ccd_t, ccd_r, & - min(logdens,max_ccd_d), & - min(logtemp,max_ccd_t))) - ! line radiation - plt = 1.d-6*10.0**(interpolate(ine, plt_d, ite, plt_t, plt_r, & - min(logdens,max_plt_d), & - min(logtemp,max_plt_t))) - ! continuum radiation - prb = 1.d-6*10.0**(interpolate(ine, prb_d, ite, prb_t, prb_r, & - min(logdens,max_prb_d), & - min(logtemp,max_prb_t))) - - end subroutine get_h_rates - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine read_atomdat(density,temperature,rates, ine, ite, filename, verbose) - !! - !! author: M Kovari, CCFE, Culham Science Centre - !! density : input real : electron density (m^-3) - !! temperature : input real : electron temperature (eV) - !! s : input : ionisation rate coefficient (m^3/s) - !! al : input : recombination rate coefficient (m^3/s) - !! Rcx : input : charge exchange rate coefficient (m^3/s) - !! plt : input : line radiation power rate coefficient (Wm^3) - !! prb : input : continuum radiation power rate coefficient (Wm^3) - !! mass : input : relative atomic mass for CX rate coefficient - !! verbose : output : verbose switch - !! - !! - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - implicit none - - ! Subroutine declarations ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!! - - real(dp), dimension(ine), intent(out) :: density - real(dp), dimension(ite), intent(out) :: temperature - real(dp), dimension(ine,ite), intent(out) :: rates - character(len=200), intent(in) :: filename - logical, intent(in), optional::verbose - integer, intent(in)::ine,ite - integer::i, l - - if (verbose) write(*,*) filename - - ! Open the input/output external files - open(unit=8,file=filename,status='old') - read(8,*) - read(8,*) - read(8,'(8(f10.5))')(density(i), i= 1, ine) - read(8,'(8(f10.5))')(temperature(i), i= 1, ite) - - read(8,*) - read(8,'(8(f10.5))')((rates(I,L),I=1,ine),L=1,ite) - close(unit=8) - - end subroutine read_atomdat - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine unit_test_read() - ! Radiative cooling function Lz - ! To test the interpolation, use a point at the geometrical mean of the two - ! first temperatures and the two first values of ne.tau - use read_radiation, only: read_lz - implicit none - - real(dp):: s, al, Rcx, plt, prb, density, temperature, mass - real(dp):: te,netau,test_lz,estimate_lz - - te=sqrt(1.000 * 1.047) - netau=sqrt(0.1*1.0) - - test_lz= read_lz('Ar', te, netau, mean_z=.false., mean_qz=.false., verbose=.true.) - write(*,*) 'Compare Lz values from code with values interpolated by hand:' - write(*,*)'test_lz = ',test_lz - estimate_lz = sqrt(sqrt(1.772e-37))*sqrt(sqrt(2.916e-37))*sqrt(sqrt(1.772e-37))*sqrt(sqrt(2.916e-37)) - write(*,*)'estimate_lz = ',estimate_lz - - density = 10.0**((7.69897 + 8.00000)/2.0 +6.0) - temperature = 10.0**((-0.69897 -0.52288)/2.0 ) - mass=2.0 - write(*,'(5(e12.5),2x)')density, temperature, mass - - call get_h_rates(density, temperature, & - s, al, Rcx, plt, prb, & - mass, verbose=.true.) - - write(*,*) 'Compare atomic rate values from code with values interpolated by hand:' - write(*,'(5(e12.5),2x)')s, al, Rcx, plt, prb - - s = 10.0**((-37.87042 -37.83901 -27.99641 -27.97238)/4.0) / 1.0e6 - al = 10.0**((-11.85004 -11.83456 -11.99890 -11.98826)/4.0) / 1.0e6 - Rcx = 10.0**((-8.59260 -8.59260 -8.50455 -8.50455)/4.0) / 1.0e6 - plt= 10.0**((-47.19666d0 -47.19666d0 -39.88720d0 -39.88720d0)/4.0d0) / 1.0d6 - prb= 10.0**((-29.50732 -29.49208 -29.65257 -29.64218)/4.0) / 1.0e6 - write(*,'(5(e12.5),2x)')s, al, Rcx, plt, prb - return - end subroutine unit_test_read - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine plot_rates() - ! Reads rate coefficients for deuterium. - ! Compare to Figure 2 in Kallenbach 2016. - real(dp):: s(3), al(3), Rcx - real(dp):: plt(3), prb(3), mass - real(dp):: lz_deuterium(3) - real(dp):: dummy1, dummy2, dummy3, dummy4, dummy5 - integer::i,j - real(dp), parameter ::te(15)=(/1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,12.,14.,16.,18.,20./) - real(dp), parameter ::density(3)=(/1.e19,1.e20,1.e21/) - - open(unit=12,file='rate_coefficients.txt',status='replace') - write(12,'(30a11)')'te [eV]','Rcx', 'ionis19', 'recomb19', 'line rad19', 'cont rad19', 'tot rad19',& - 'ionis20', 'recomb20', 'line rad20', 'cont rad20', 'tot rad20',& - 'ionis21', 'recomb21', 'line rad21', 'cont rad21', 'tot rad21' - mass=2.0D0 - ! Just read data - call get_h_rates(1.d20, 1.0d0, dummy1, dummy2, dummy3, dummy4, dummy5, & - mass, verbose=.true.) - do i=1,15 - do j=1,3 - call get_h_rates(density(j), te(i), s(j), al(j), Rcx, plt(j), prb(j), & - mass, verbose=.false.) - lz_deuterium(j)=plt(j)+prb(j) - enddo - write(12,'(30es11.3)')te(i), Rcx, (s(j), al(j), plt(j), prb(j), lz_deuterium(j), j=1,3) - enddo - close(unit=12) - - end subroutine plot_rates - - -end module read_and_get_atomic_data diff --git a/source/fortran/read_radiation.f90 b/source/fortran/read_radiation.f90 deleted file mode 100755 index b48c2065..00000000 --- a/source/fortran/read_radiation.f90 +++ /dev/null @@ -1,344 +0,0 @@ -#ifndef INSTALLDIR -#error INSTALLDIR not defined! -#endif - -module read_radiation - !! Module for reading radiation data - !! author: M Kovari, CCFE, Culham Science Centre - !! N/A - !! This module contains the PROCESS Kallenbach divertor model - !! - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Modules to import ! - ! !!!!!!!!!!!!!!!!!!!! - -#ifndef dp - use, intrinsic :: iso_fortran_env, only: dp=>real64 -#endif - use impurity_radiation_module, only: nimp - implicit none - - ! List of impurities in the SOL/divertor model IS now same as the main plasma impurities - - ! Vars in function read_lz requiring re-initialisation on each new run - integer :: location - - ! Length of temperature and netau data - integer :: nt, nnetau - - ! Impurity data array - real(dp), dimension(200,5) :: impurity_data - - ! The values of the Lz, mean Z, and mean Z^2 stored in the data files - real(dp), dimension(nimp,200,5) :: data_lz, data_z, data_qz - - ! The values of ln(n.tau) at which the data is stored. - real(dp), dimension(nimp,5) :: lnetau_lz, lnetau_z, lnetau_qz - - ! The values of ln(Te) at which the data is stored. - real(dp), dimension(nimp, 200) :: logT_lz, logT_z, logT_qz - - ! First call boolean switch array - logical :: FirstCall(nimp) - -contains - character(len=300) function lzdir() - implicit none - character(len=200) :: process_dir - CALL get_environment_variable("PYTHON_PROCESS_ROOT", process_dir) - if (process_dir == "") then - lzdir = INSTALLDIR//'/process/data/lz_non_corona_14_elements/' - else - lzdir = trim(process_dir)//'/data/lz_non_corona_14_elements/' - end if - end function lzdir - - subroutine init_read_radiation - !! Initialise module variables - implicit none - - location = 0 - nt = 0 - nnetau = 0 - impurity_data = 0.0D0 - data_lz = 0.0D0 - data_z = 0.0D0 - data_qz = 0.0D0 - lnetau_lz = 0.0D0 - lnetau_z = 0.0D0 - lnetau_qz = 0.0D0 - logT_lz = 0.0D0 - logT_z = 0.0D0 - logT_qz = 0.0D0 - FirstCall = .true. - end subroutine init_read_radiation - - function read_lz(element, te, netau, mean_z, mean_qz, verbose) - !! Read Lz, mean_z or mean_qz data from database - !! author: M Koviar, CCFE, Culham Science Centre - !! element : input char : element name - !! te : input real : electron temperature [eV] - !! netau : input real : "non-coronal parameter" for radiative loss function [ms.1e20/m3] - !! mean_z : input logical : get mean charge instead of Lz - !! mean_qz : input logical : get mean quadratic charge instead of Lz - !! verbose : input logical : verbose boolean - !! read radiative loss data Lz calculated by M. O'Mullane from ADAS in 2016, - !! OR mean charge - !! A. Kallenbach, 11.1.2013, R. Dux (add readout of Z, Z^2, add a few elements) - !! Fortran: Michael Kovari 2/6/16 - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Modules to import ! - ! !!!!!!!!!!!!!!!!!!!! - - use maths_library, only: interpolate - use impurity_radiation_module, only: nimp, imp_label - implicit none - - ! Variable declarations ! - ! !!!!!!!!!!!!!!!!!!!!!!!! - - integer :: i, j - - ! "non-coronal parameter" for radiative loss function [ms.1e20/m3] - real(dp), intent(in) :: netau - - ! Electron temperature [eV] - real(dp), intent(in) :: te - - ! Verbose switch - logical, intent(in) :: verbose - - ! switches for geting mean charge or mean quadratic charge instead of Lz - logical, intent(in) :: mean_z, mean_qz - - ! Element name - character(len=*), intent(in) :: element - - ! Function output - real(dp) :: read_lz - - ! Natural logs of netau and te - real(dp) :: lnetau, lte - - ! Lz data filename - ! Beware: changing length can affect Kallenbach solutions / regression test - character(len=300) :: filename - - ! Find the index of the element. Exclude hydrogen by starting at 2 (Helium) - do i = 2, nimp - if (imp_label(i) .eq. element) then - location = i - exit - endif - enddo - - ! Check element was found - if (location.eq.0) then - write(*,*)'element '//element//' not supported' - end if - - ! Natural log netau and electron temperature - lnetau = log(netau) - lte = log(te) - - ! Read the data for that element - - ! FirstCall - if(FirstCall(location)) then - - FirstCall(location)=.false. - - ! Each data file may have its own values of T and n, stored in arrays logT_lz, lnetau_lz etc - ! Store data in logarithm form for interpolation - - ! Assign loss data filename - filename = trim(lzdir())//trim(element)//'_lz_tau.dat' - - ! Read the impuriy data - call read_impurity_data(filename, nt, nnetau, impurity_data, logT_lz(location,:), lnetau_lz(location,:)) - - ! Store log impurity data - data_lz(location,:,:) = log(impurity_data) - - if (verbose) then - write(*,*)'Lz' - write(*,*)'log(T) values' - write(*,'(8(f10.3))')(logT_lz(location,i), i= 1, nt) - do j=1,nnetau - write(*,*)'log(ne.tau)= ',lnetau_lz(location,j) - write(*,'(8(e10.3))')(impurity_data(i,j), i= 1, nt) - enddo - endif - - ! Assign z data filename - filename = trim(lzdir())//trim(element)//'_z_tau.dat' - - ! Read z data - call read_impurity_data(filename, nt, nnetau, impurity_data, logT_z(location,:), lnetau_z(location,:)) - - if (verbose) then - write(*,*)'Element = '// element//'./data/LZ_NON_CORONA/'//element//'_z_tau.dat' - write(*,*)'Mean Z' - write(*,*)'log(T) values' - write(*,'(8(f10.3))')(logT_z(location,i), i= 1, nt) - do j=1,nnetau - write(*,*)'log(ne.tau)= ',lnetau_z(location,j) - write(*,'(8(e10.3))')(impurity_data(i,j), i= 1, nt) - enddo - endif - - ! Store log z data - data_z(location,:,:) = log(impurity_data) - - ! Assign z^2 data filename - filename = trim(lzdir())//trim(element)//'_z2_tau.dat' - - ! Read root-mean square z data - call read_impurity_data(filename, nt, nnetau, impurity_data, logT_qz(location,:), lnetau_qz(location,:)) - - ! Square RMS z data - impurity_data = impurity_data**2 - - ! Store z^2 data - data_qz(location,:,:) = log(impurity_data) - - if (verbose) then - write(*,*)'Mean Z^2' - write(*,*)'log(T) values' - write(*,'(8(f10.3))')(logT_qz(location,i), i= 1, nt) - do j=1,nnetau - write(*,*)'log(ne.tau)= ',lnetau_qz(location,j) - write(*,'(8(e10.3))')(impurity_data(i,j), i= 1, nt) - enddo - endif - - endif - - ! Interpolate for the specified values of log(temperature) and log(n.tau) - ! Using function interpolate(x_len, x_array, y_len, y_array, f, x, y, delta) - if(mean_z) then - - read_lz = interpolate(nt, logT_z(location,:), nnetau, lnetau_z(location,:), data_z(location,:,:), lte, lnetau) - - elseif(mean_qz) then - - read_lz = interpolate(nt, logT_qz(location,:), nnetau, lnetau_qz(location,:), data_qz(location,:,:), lte, lnetau) - - else - - if( (minval(lnetau_lz(location,:)).gt.lnetau) .or. (maxval(lnetau_lz(location,:)).lt.lnetau) ) then - write(*,*)'lnetau is out of range, =', lnetau - endif - - read_lz = interpolate(nt, logT_lz(location,:), nnetau, lnetau_lz(location,:), data_lz(location,:,:), lte, lnetau) - - endif - - read_lz = exp(read_lz) - - end function read_lz - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine read_impurity_data(filename,nt,nnetau,impurity_data,logT_array,lnetau_array) - !! Read impurity data from ADAS database - !! author: M Kovari, CCFE, Culham Science Centre - !! filename : input character : filename of data - !! nt : output integer : length of temperature data - !! nnetau : output integer : length of netau data - !! impurity_data : output real array : impurity data array - !! logT_array : output real array : log temperature array - !! lnetau_array : output real array : log netau array - !! Read the impurity data from database - !! - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - implicit none - - ! Subroutine declarations ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Filename for data - character(len=*), intent(in) :: filename - - ! Length of temperature and netau arrays - integer, intent(out) :: nt, nnetau - - ! Impurity data array - real(dp), intent(out), dimension(200,5) :: impurity_data - - ! Log netau and temperature arrays - real(dp), intent(out), dimension(5) :: lnetau_array - real(dp), intent(out), dimension(200) :: logT_array - - ! netau data - real(dp) :: data_netau - - ! - integer :: pos, i, j, iostatus - - ! Electron temperature array [eV] - real(dp), dimension(200) :: T_array - - character(len=100) :: string, substring - - ! Open data file - open(unit=8,file=filename,status='old') - read(8,*) - read(8,*) - read(8,*) - read(8,*)nt, nnetau - - if ((nt.ne.200).or.(nnetau.ne.5))then - write(*,*)'read_impurity_data: There must be exactly 200 temperature values and 5 values of n.tau' - write(*,*)'nt = ',nt,' nnetau = ',nnetau - endif - - do - read(8,*,IOSTAT=iostatus)string - - ! The string 'Te[eV]' is present - if(index(string,'Te[eV]').ne.0) exit - - if(iostatus.ne.0)then - write(*,*)'Problem in reading impurity data from file ',filename - stop 1 - endif - - enddo - - read(8,'(8(f10.3))')(T_array(i), i= 1, nt) - - ! Create log temperature array [eV] - logT_array = log(T_array) - - ! Loop over netau - do j = 1, nnetau - - ! Read the non-coronal parameter ne.tau - read(8,'(a100)')string - - pos = index(string,'=') - - substring = string(pos+1:) - - read(substring,*)data_netau - - lnetau_array(j) = log(data_netau) - - ! Read the data - ! A fixed format read doesn't work for some reason - read(8,*)(impurity_data(i,j), i= 1, nt) - - enddo - - close(unit=8) - - end subroutine read_impurity_data - -end module read_radiation diff --git a/source/fortran/real_mod.f90 b/source/fortran/real_mod.f90 deleted file mode 100755 index 7070c371..00000000 --- a/source/fortran/real_mod.f90 +++ /dev/null @@ -1,75 +0,0 @@ -module real_mod - - interface real8 - - function real8_real_4_(x) - use precision_mod, only: real_4_, p_ - implicit none - real(p_) :: real8_real_4_ - real(real_4_), intent(in) :: x - !return real(x, real_8_) - end function real8_real_4_ - - function real8_real_8_(x) - use precision_mod, only: real_8_, p_ - implicit none - real(p_) :: real8_real_8_ - real(real_8_), intent(in) :: x - !return real(x, real_8_) - end function real8_real_8_ - - function real8_i(x) - use precision_mod, only: p_ - implicit none - real(p_) :: real8_i - integer, intent(in) :: x - !return real(x, real_8_) - end function real8_i - - function real8_c(x) - use precision_mod, only: real_4_, p_ - implicit none - real(p_) :: real8_c - complex(real_4_), intent(in) :: x - !return real(x, real_8_) - end function real8_c - - function real8_c8(x) - use precision_mod, only: real_8_, p_ - implicit none - real(p_) :: real8_c8 - complex(real_8_), intent(in) :: x - !return real(x, real_8_) - end function real8_c8 - - end interface - - interface cmplx8 - - function cmplx8_real_4_(x, y) - use precision_mod, only: real_4_, p_ - implicit none - complex(p_) :: cmplx8_real_4_ - real(real_4_), intent(in) :: x, y - !return cmplx(x, y, real_8_) - end function cmplx8_real_4_ - - function cmplx8_real_8_(x, y) - use precision_mod, only: real_8_, p_ - implicit none - complex(p_) :: cmplx8_real_8_ - real(real_8_), intent(in) :: x, y - !return cmplx(x, y, real_8_) - end function cmplx8_real_8_ - - function cmplx8_i(x, y) - use precision_mod, only: p_ - implicit none - complex(p_) :: cmplx8_i - integer, intent(in) :: x, y - !return cmplx(x, y, real_8_) - end function cmplx8_i - - end interface - -end module real_mod diff --git a/tests/unit/test_read_radiation.py b/tests/unit/test_read_radiation.py deleted file mode 100644 index 80390994..00000000 --- a/tests/unit/test_read_radiation.py +++ /dev/null @@ -1,40 +0,0 @@ -from typing import NamedTuple -import pytest - -from process.fortran import read_radiation, impurity_radiation_module - - -@pytest.fixture(scope="function") -def init_read_radiation(): - read_radiation.init_read_radiation() - impurity_radiation_module.init_impurity_radiation_module() - - -class ReadLzParam(NamedTuple): - element: str - te: float - netau: float - expected: float - mean_z: bool = False - mean_qz: bool = False - - -@pytest.mark.parametrize( - "read_lz_param", - ( - ReadLzParam("He", 5, 0.5, 5.958312853655463e-35), - ReadLzParam("He", 124.45, 0.5, 1.9852736835073994, mean_z=True), - ReadLzParam("He", 124.45, 0.5, 3.95858840328217, mean_qz=True), - ), -) -def test_read_lz(read_lz_param: ReadLzParam): - ret = read_radiation.read_lz( - read_lz_param.element, - read_lz_param.te, - read_lz_param.netau, - read_lz_param.mean_z, - read_lz_param.mean_qz, - False, - ) - - assert ret == pytest.approx(read_lz_param.expected)