ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
itm_toolbox.F90
Go to the documentation of this file.
1 
10 
11 module itm_toolbox
12 
13  use itm_constants
14 
15  implicit none
16 
17  integer, parameter :: BUFLEN = 256
18 
19 contains
20 
28 
29  function interpolate(x_in, y_in, x_out)
30 
31  implicit none
32 
33  REAL(R8), INTENT(IN) :: x_in(:), y_in(:), x_out(:)
34  REAL(R8) :: interpolate(size(x_out)), y_out(size(x_out))
35 
36  call l3interp(y_in, x_in, size(x_in), y_out, x_out, size(x_out))
37  interpolate = y_out
38 
39  return
40  end function interpolate
41 
42 
43 
54 
55  function rho_with_accumulation(function_string, x)
56 
57  implicit none
58 
59  real(R8) :: x(:), rho_with_accumulation(1:size(x))
60  character (len=BUFLEN) :: function_string
61  real(R8), ALLOCATABLE :: rho_1(:), rho_2(:), rho_3(:), rho_4(:), rho_5(:)
62  real(R8) :: dummy1
63  integer :: nrho, irho
64 
65  nrho=size(x)
66  allocate(rho_1(0:10*nrho), rho_2(0:10*nrho), rho_3(0:10*nrho), rho_4(nrho), rho_5(nrho))
67  rho_1 = (/ (1.0_r8/(10*nrho) * (irho), irho=0,nrho*10) /)
68  rho_2 = profile(function_string, rho_1)
69  rho_3(0)=0
70  do irho=1,10*nrho
71  call cubint(10*nrho+1, rho_1, rho_2, 1, irho+1, rho_3(irho), dummy1)
72  enddo
73  rho_3=rho_3/rho_3(10*nrho)
74  rho_4 = (/ (1.0_r8/(nrho-1) * (irho-1), irho=1,nrho) /)
75  call l3interp(rho_1, rho_3, 10*nrho+1, rho_5, rho_4, nrho)
76  rho_with_accumulation = rho_5
77  deallocate(rho_1, rho_2, rho_3, rho_4, rho_5)
78 
79  end function rho_with_accumulation
80 
81 
82 
91 
92  function profile(function_string, x)
93 
94  use fortranparser, only : equationparser
95 
96  implicit none
97 
98  real(R8) :: x(:), profile(1:size(x))
99  character (len=BUFLEN) :: function_string
100 
101  type(equationparser) :: function_descriptor
102  character(len=10) :: variables(1) = ['x']
103 
104  integer :: i
105 
106  function_descriptor = equationparser(trim(function_string), variables)
107 ! if(.not.c_associated(function_descriptor)) then
108 ! write(*,*) 'Invalid function ', trim(function_string)
109 ! stop
110 ! endif
111 
112  do i = 1, size(x)
113  profile(i) = function_descriptor%evaluate([x(i)])
114  enddo
115 
116 ! call evaluator_destroy(function_descriptor)
117 
118  end function profile
119 
120  subroutine cubint ( ntab, xtab, ftab, ia_in, ib_in, result, error )
121  USE itm_types, only : r8
122 
123  ! obtained from
124  ! http://people.sc.fsu.edu/~jburkardt/f_src/intlib/intlib.html
125 
126  !*****************************************************************************80
127  !
128  !! CUBINT approximates an integral using cubic interpolation of data.
129  !
130  ! Discussion:
131  !
132  ! The integral to be approximated is
133  !
134  ! Integral ( XTAB(IB) <= X <= XTAB(IA) ) F(X) DX
135  !
136  ! The routine estimates the error in integration.
137  !
138  ! Modified:
139  !
140  ! 10 February 2006
141  !
142  ! Reference:
143  !
144  ! Philip Davis, Philip Rabinowitz,
145  ! Methods of Numerical Integration,
146  ! Second Edition,
147  ! Dover, 2007,
148  ! ISBN: 0486453391,
149  ! LC: QA299.3.D28.
150  !
151  ! Philip Gill, GF Miller,
152  ! An algorithm for the integration of unequally spaced data,
153  ! The Computer Journal,
154  ! Number 15, Number 1, 1972, pages 80-83.
155  !
156  ! Parameters:
157  !
158  ! Input, integer NTAB, the number of tabulated points.
159  ! NTAB must be at least 4.
160  !
161  ! Input, real (R8) XTAB(NTAB), contains the points at which the
162  ! function was tabulated. XTAB should contain distinct
163  ! values, given in ascending order.
164  !
165  ! Input, real (R8) FTAB(NTAB), contains the tabulated function
166  ! values, FTAB(I) = F(XTAB(I)).
167  !
168  ! Input, integer IA, the entry of XTAB at which integration
169  ! is to begin. IA must be no less than 1 and no greater
170  ! than NTAB.
171  !
172  ! Input, integer IB, the entry of XTAB at which integration
173  ! is to end. IB must be no less than 1 and no greater than
174  ! NTAB.
175  !
176  ! Output, real (R8) RESULT, the approximate value of the
177  ! integral from XTAB(IA) to XTAB(IB) of the function.
178  !
179  ! Output, real (R8) ERROR, an estimate of the error in
180  ! integration.
181  !
182  implicit none
183 
184  integer ntab
185 
186  real (R8) c
187  real (R8) d1
188  real (R8) d2
189  real (R8) d3
190  real (R8) error
191  real (R8) ftab(ntab)
192  real (R8) h1
193  real (R8) h2
194  real (R8) h3
195  real (R8) h4
196  integer i
197  integer ia, ia_in
198  integer ib, ib_in
199  integer ind
200  integer it
201  integer j
202  integer k
203  real (R8) r1
204  real (R8) r2
205  real (R8) r3
206  real (R8) r4
207  real (R8) result
208  real (R8) s
209  real (R8) term
210  real (R8) xtab(ntab)
211 
212  ia = ia_in
213  ib = ib_in
214  result = 0.0_r8
215  error = 0.0_r8
216 
217  if ( ia == ib ) then
218  return
219  end if
220 
221  if ( ntab < 4 ) then
222  write ( *, '(a)' ) ' '
223  write ( *, '(a)' ) 'CUBINT - Fatal error!'
224  write ( *, '(a,i8)' ) ' NTAB must be at least 4, but input NTAB = ', ntab
225  stop
226  end if
227 
228  if ( ia < 1 ) then
229  write ( *, '(a)' ) ' '
230  write ( *, '(a)' ) 'CUBINT - Fatal error!'
231  write ( *, '(a,i8)' ) ' IA must be at least 1, but input IA = ', ia
232  stop
233  end if
234 
235  if ( ntab < ia ) then
236  write ( *, '(a)' ) ' '
237  write ( *, '(a)' ) 'CUBINT - Fatal error!'
238  write ( *, '(a,i8)' ) ' IA must be <= NTAB, but input IA = ', ia
239  stop
240  end if
241 
242  if ( ib < 1 ) then
243  write ( *, '(a)' ) ' '
244  write ( *, '(a)' ) 'CUBINT - Fatal error!'
245  write ( *, '(a,i8)' ) ' IB must be at least 1, but input IB = ', ib
246  stop
247  end if
248 
249  if ( ntab < ib ) then
250  write ( *, '(a)' ) ' '
251  write ( *, '(a)' ) 'CUBINT - Fatal error!'
252  write ( *, '(a,i8)' ) ' IB must be <= NTAB, but input IB = ', ib
253  stop
254  end if
255  !
256  ! Temporarily switch IA and IB, and store minus sign in IND
257  ! so that, while integration is carried out from low X's
258  ! to high ones, the sense of the integral is preserved.
259  !
260  if ( ib < ia ) then
261  ind = -1
262  it = ib
263  ib = ia
264  ia = it
265  else
266  ind = 1
267  end if
268 
269  s = 0.0_r8
270  c = 0.0_r8
271  r4 = 0.0_r8
272  j = ntab-2
273  if ( ia < ntab-1 .or. ntab == 4 ) then
274  j = max( 3, ia )
275  end if
276 
277  k = 4
278  if ( 2 < ib .or. ntab == 4 ) then
279  k = min( ntab, ib + 2 ) - 1
280  end if
281 
282  do i = j, k
283 
284  if ( i <= j ) then
285 
286  h2 = xtab(j-1) - xtab(j-2)
287  d3 = ( ftab(j-1) - ftab(j-2) ) / h2
288  h3 = xtab(j) - xtab(j-1)
289  d1 = ( ftab(j) - ftab(j-1) ) / h3
290  h1 = h2 + h3
291  d2 = ( d1 - d3 ) / h1
292  h4 = xtab(j+1) - xtab(j)
293  r1 = ( ftab(j+1) - ftab(j) ) / h4
294  r2 = ( r1 - d1 ) / ( h4 + h3 )
295  h1 = h1 + h4
296  r3 = (r2-d2) / h1
297 
298  if ( ia <= 1 ) then
299  result = h2 * ( ftab(1) + h2 * ( 0.5_r8 * d3 - h2 &
300  * ( d2 / 6.0_r8 -(h2+h3+h3)*r3/12.0_r8)))
301  s = -h2**3 * (h2*(3.0_r8*h2+5.0_r8*h4)+10.0_r8*h3*h1) / 60.0_r8
302  end if
303 
304  else
305 
306  h4 = xtab(i+1) - xtab(i)
307  r1 = ( ftab(i+1) - ftab(i) ) / h4
308  r4 = h4 + h3
309  r2 = ( r1 - d1 ) / r4
310  r4 = r4 + h2
311  r3 = ( r2 - d2 ) / r4
312  r4 = ( r3 - d3 ) / ( r4 + h1 )
313 
314  end if
315 
316  if ( ia < i .and. i <= ib ) then
317 
318  term = h3 * ( ( ftab(i) + ftab(i-1) ) * 0.5_r8 &
319  -h3 * h3 * ( d2 + r2 + ( h2 - h4 ) * r3 ) / 12.0_r8 )
320  result = result + term
321  c = h3**3 * ( 2.0_r8 * h3 * h3 &
322  + 5.0_r8 * ( h3 * ( h4 + h2 ) + 2.0_r8 * h2 * h4 ) ) / 120.0_r8
323  error = error + (c+s)*r4
324 
325  if ( i /= j ) then
326  s = c
327  else
328  s = s + c + c
329  end if
330 
331  else
332 
333  error = error + r4 * s
334 
335  end if
336 
337  if ( k <= i ) then
338 
339  if ( ntab <= ib ) then
340  term = h4 * ( ftab(ntab) - h4 * ( 0.5 * r1 &
341  + h4 * ( r2 / 6.0_r8 + ( h3 + h3 + h4 ) * r3 / 12.0_r8 )))
342  result = result + term
343  error = error - h4**3 * r4 * &
344  ( h4 * ( 3.0_r8 * h4 + 5.0_r8 * h2 ) &
345  + 10.0_r8 * h3 * ( h2 + h3 + h4 ) ) / 60.0_r8
346  end if
347 
348  if ( ntab-1 <= ib ) then
349  error = error + s * r4
350  end if
351 
352  else
353 
354  h1 = h2
355  h2 = h3
356  h3 = h4
357  d1 = r1
358  d2 = r2
359  d3 = r3
360  end if
361 
362  end do
363  !
364  ! Restore original values of IA and IB, reverse signs
365  ! of RESULT and ERROR, to account for integration
366  ! that proceeded from high X to low X.
367  !
368  if ( ind /= 1 ) then
369  it = ib
370  ib = ia
371  ia = it
372  result = -result
373  error = -error
374  end if
375 
376  return
377  end subroutine cubint
378 
379  subroutine avint ( ntab, xtab, ytab, a, b, result )
380  USE itm_types, only : r8
381 
382  ! obtained from
383  ! http://people.sc.fsu.edu/~jburkardt/f_src/intlib/intlib.html
384 
385  !*****************************************************************************80
386  !
387  !! AVINT estimates the integral of unevenly spaced data.
388  !
389  ! Discussion:
390  !
391  ! The data is given as NTAB pairs of values
392  ! ( XTAB(1:NTAB), YTAB(1:NTAB) ).
393  !
394  ! The quadrature method uses overlapping parabolas and smoothing.
395  !
396  ! Modified:
397  !
398  ! 10 February 2006
399  !
400  ! Reference:
401  !
402  ! Philip Davis, Philip Rabinowitz,
403  ! Methods of Numerical Integration,
404  ! Second Edition,
405  ! Dover, 2007,
406  ! ISBN: 0486453391,
407  ! LC: QA299.3.D28.
408  !
409  ! Paul Hennion,
410  ! Algorithm 77:
411  ! Interpolation, Differentiation and Integration,
412  ! Communications of the ACM,
413  ! Volume 5, page 96, 1962.
414  !
415  ! Parameters:
416  !
417  ! Input, integer NTAB, the number of entries in XTAB and
418  ! YTAB. NTAB must be at least 2.
419  !
420  ! Input, real (R8) XTAB(NTAB), the abscissas at which the
421  ! function values are given. The XTAB's must be distinct
422  ! and in ascending order.
423  !
424  ! Input, real (R8) YTAB(NTAB), the function values,
425  ! YTAB(I) = F(XTAB(I)).
426  !
427  ! Input, real (R8) A, the lower limit of integration. A should
428  ! be, but need not be, near one endpoint of the interval
429  ! (X(1), X(NTAB)).
430  !
431  ! Input, real (R8) B, the upper limit of integration. B should
432  ! be, but need not be, near one endpoint of the interval
433  ! (X(1), X(NTAB)).
434  !
435  ! Output, real (R8) RESULT, the approximate value of the integral.
436  !
437  implicit none
438 
439  integer ntab
440 
441  real (R8) a
442  real (R8) b
443  real (R8) ba
444  real (R8) bb
445  real (R8) bc
446  real (R8) ca
447  real (R8) cb
448  real (R8) cc
449  real (R8) fa
450  real (R8) fb
451  integer i
452  integer inlft
453  integer inrt
454  integer istart
455  integer istop
456  real (R8) result
457  real (R8) slope
458  real (R8) syl
459  real (R8) syl2
460  real (R8) syl3
461  real (R8) syu
462  real (R8) syu2
463  real (R8) syu3
464  real (R8) term1
465  real (R8) term2
466  real (R8) term3
467  real (R8) total
468  real (R8) x1
469  real (R8) x12
470  real (R8) x13
471  real (R8) x2
472  real (R8) x23
473  real (R8) x3
474  real (R8) xtab(ntab)
475  real (R8) ytab(ntab)
476 
477  result = 0.0_r8
478 
479  if ( a == b ) then
480  return
481  end if
482 
483  if ( b < a ) then
484  end if
485 
486  if ( ntab < 2 ) then
487  write ( *, '(a)' ) ' '
488  write ( *, '(a)' ) 'AVINT - Fatal error!'
489  write ( *, '(a,i8)' ) ' NTAB is less than 3. NTAB = ', ntab
490  stop
491  end if
492 
493  do i = 2, ntab
494 
495  if ( xtab(i) <= xtab(i-1) ) then
496  write ( *, '(a)' ) ' '
497  write ( *, '(a)' ) 'AVINT - Fatal error!'
498  write ( *, '(a)' ) ' XTAB(I) is not greater than XTAB(I-1).'
499  write ( *, '(a,i8)' ) ' Here, I = ', i
500  write ( *, '(a,g14.6)' ) ' XTAB(I-1) = ', xtab(i-1)
501  write ( *, '(a,g14.6)' ) ' XTAB(I) = ', xtab(i)
502  stop
503  end if
504 
505  end do
506  !
507  ! Special case for NTAB = 2.
508  !
509  if ( ntab == 2 ) then
510  slope = ( ytab(2) - ytab(1) ) / ( xtab(2) - xtab(1) )
511  fa = ytab(1) + slope * ( a - xtab(1) )
512  fb = ytab(2) + slope * ( b - xtab(2) )
513  result = 0.5_r8 * ( fa + fb ) * ( b - a )
514  return
515  end if
516 
517  if ( xtab(ntab-2) < a .or. b < xtab(3) ) then
518  write ( *, '(a)' ) ' '
519  write ( *, '(a)' ) 'AVINT - Fatal error!'
520  write ( *, '(a)' ) ' There were less than 3 function values'
521  write ( *, '(a)' ) ' between the limits of integration.'
522  stop
523  end if
524 
525  i = 1
526  do
527 
528  if ( a <= xtab(i) ) then
529  exit
530  end if
531 
532  i = i + 1
533 
534  end do
535 
536  inlft = i
537 
538  i = ntab
539 
540  do
541 
542  if ( xtab(i) <= b ) then
543  exit
544  end if
545 
546  i = i - 1
547 
548  end do
549 
550  inrt = i
551 
552  if ( inrt - inlft < 2 ) then
553  write ( *, '(a)' ) ' '
554  write ( *, '(a)' ) 'AVINT - Fatal error!'
555  write ( *, '(a)' ) ' There were less than 3 function values'
556  write ( *, '(a)' ) ' between the limits of integration.'
557  stop
558  end if
559 
560  if ( inlft == 1 ) then
561  istart = 2
562  else
563  istart = inlft
564  end if
565 
566  if ( inrt == ntab ) then
567  istop = ntab - 1
568  else
569  istop = inrt
570  end if
571 
572  total = 0.0_r8
573 
574  syl = a
575  syl2 = syl * syl
576  syl3 = syl2 * syl
577 
578  do i = istart, istop
579 
580  x1 = xtab(i-1)
581  x2 = xtab(i)
582  x3 = xtab(i+1)
583 
584  x12 = x1 - x2
585  x13 = x1 - x3
586  x23 = x2 - x3
587 
588  term1 = ( ytab(i-1) ) / ( x12 * x13 )
589  term2 = - ( ytab(i) ) / ( x12 * x23 )
590  term3 = ( ytab(i+1) ) / ( x13 * x23 )
591 
592  ba = term1 + term2 + term3
593  bb = - ( x2 + x3 ) * term1 - ( x1 + x3 ) * term2 - ( x1 + x2 ) * term3
594  bc = x2 * x3 * term1 + x1 * x3 * term2 + x1 * x2 * term3
595 
596  if ( i == istart ) then
597  ca = ba
598  cb = bb
599  cc = bc
600  else
601  ca = 0.5_r8 * ( ba + ca )
602  cb = 0.5_r8 * ( bb + cb )
603  cc = 0.5_r8 * ( bc + cc )
604  end if
605 
606  syu = x2
607  syu2 = syu * syu
608  syu3 = syu2 * syu
609 
610  total = total + ca * ( syu3 - syl3 ) / 3.0_r8 &
611  + cb * ( syu2 - syl2 ) / 2.0_r8 &
612  + cc * ( syu - syl )
613  ca = ba
614  cb = bb
615  cc = bc
616 
617  syl = syu
618  syl2 = syu2
619  syl3 = syu3
620 
621  end do
622 
623  syu = b
624  syu2 = syu * syu
625  syu3 = syu2 * syu
626 
627  result = total + ca * ( syu3 - syl3 ) / 3.0_r8 &
628  + cb * ( syu2 - syl2 ) / 2.0_r8 &
629  + cc * ( syu - syl )
630 
631  return
632  end subroutine avint
633 
639 
640  SUBROUTINE l3interp(y_in,x_in,nr_in,y_out,x_out,nr_out)
641 
642  USE itm_types
643 
644  IMPLICIT NONE
645 
646  INTEGER(ITM_I4) :: nr_in,nr_out
647  REAL(R8), DIMENSION(nr_in), INTENT(IN) :: y_in,x_in
648  REAL(R8), DIMENSION(nr_out), INTENT(IN) :: x_out
649  REAL(R8), DIMENSION(nr_out), INTENT(OUT) :: y_out
650 
651  REAL(R8) :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
652  INTEGER(ITM_I4) :: j,jm,j0,j1,j2
653  INTEGER :: jstart,jfirst,jlast,jstep
654 
655  IF (x_in(nr_in) > x_in(1)) THEN
656  jstart=3
657  jfirst=1
658  jlast=nr_out
659  jstep=1
660  ELSE
661  jstart=nr_out-2
662  jfirst=nr_out
663  jlast=1
664  jstep=-1
665  END IF
666 
667  j1=jstart
668  DO j=jfirst,jlast,jstep
669  x=x_out(j)
670  DO WHILE (x >= x_in(j1) .AND. j1 < nr_in-1 .AND. j1 > 2)
671  j1=j1+jstep
672  END DO
673 
674  j2=j1+jstep
675  j0=j1-jstep
676  jm=j1-2*jstep
677 
678  !... extrapolate inside or outside
679 
680  x2=x_in(j2)
681  x1=x_in(j1)
682  x0=x_in(j0)
683  xm=x_in(jm)
684 
685  aintm=(x-x0)*(x-x1)*(x-x2)/((xm-x0)*(xm-x1)*(xm-x2))
686  aint0=(x-xm)*(x-x1)*(x-x2)/((x0-xm)*(x0-x1)*(x0-x2))
687  aint1=(x-xm)*(x-x0)*(x-x2)/((x1-xm)*(x1-x0)*(x1-x2))
688  aint2=(x-xm)*(x-x0)*(x-x1)/((x2-xm)*(x2-x0)*(x2-x1))
689 
690  y_out(j)=aintm*y_in(jm)+aint0*y_in(j0) &
691  +aint1*y_in(j1)+aint2*y_in(j2)
692 
693  END DO
694 
695  END SUBROUTINE l3interp
696 
697 
703 
704  SUBROUTINE l3deriv(y_in,x_in,nr_in,dydx_out,x_out,nr_out)
705 
706  USE itm_types
707 
708  IMPLICIT NONE
709 
710  INTEGER(ITM_I4) :: nr_in,nr_out
711  REAL(R8), DIMENSION(nr_in), INTENT(IN) :: y_in,x_in
712  REAL(R8), DIMENSION(nr_out), INTENT(IN) :: x_out
713  REAL(R8), DIMENSION(nr_out), INTENT(OUT) :: dydx_out
714 
715  REAL(R8) :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
716  INTEGER(ITM_I4) :: j,jm,j0,j1,j2
717  INTEGER(ITM_I4) :: jstart,jfirst,jlast,jstep
718 
719  IF (x_in(nr_in) > x_in(1)) THEN
720  jstart=3
721  jfirst=1
722  jlast=nr_out
723  jstep=1
724  ELSE
725  jstart=nr_out-2
726  jfirst=nr_out
727  jlast=1
728  jstep=-1
729  END IF
730 
731  j1=jstart
732  DO j=jfirst,jlast,jstep
733  x=x_out(j)
734  DO WHILE (x >= x_in(j1) .AND. j1 < nr_in-1 .AND. j1 > 2)
735  j1=j1+jstep
736  END DO
737 
738  j2=j1+jstep
739  j0=j1-jstep
740  jm=j1-2*jstep
741 
742  !... extrapolate inside or outside
743 
744  x2=x_in(j2)
745  x1=x_in(j1)
746  x0=x_in(j0)
747  xm=x_in(jm)
748 
749  aintm=((x-x1)*(x-x2)+(x-x0)*(x-x2)+(x-x0)*(x-x1)) &
750  /((xm-x0)*(xm-x1)*(xm-x2))
751  aint0=((x-x1)*(x-x2)+(x-xm)*(x-x2)+(x-xm)*(x-x1)) &
752  /((x0-xm)*(x0-x1)*(x0-x2))
753  aint1=((x-x0)*(x-x2)+(x-xm)*(x-x2)+(x-xm)*(x-x0)) &
754  /((x1-xm)*(x1-x0)*(x1-x2))
755  aint2=((x-x0)*(x-x1)+(x-xm)*(x-x1)+(x-xm)*(x-x0)) &
756  /((x2-xm)*(x2-x0)*(x2-x1))
757 
758  dydx_out(j)=aintm*y_in(jm)+aint0*y_in(j0) &
759  +aint1*y_in(j1)+aint2*y_in(j2)
760 
761  END DO
762 
763  END SUBROUTINE l3deriv
764 
765 
771 
772  SUBROUTINE lderiv(y,x,dydx,nx)
773 
774  USE itm_types
775 
776  IMPLICIT NONE
777 
778  INTEGER(ITM_I4) :: nx
779  REAL(R8) :: dxp,dxm,dx0
780  REAL(R8), DIMENSION(nx), INTENT(IN) :: x,y
781  REAL(R8), DIMENSION(nx), INTENT(OUT) :: dydx
782 
783  INTEGER(ITM_I4) :: i,ip,im
784 
785  DO i=2,nx-1
786  ip=min(i+1,nx)
787  im=max(i-1,1)
788  dxp=x(ip)-x(i)
789  dxm=x(i)-x(im)
790  dx0=x(ip)-x(im)
791  dydx(i)=(dxm*dxm*(y(ip)-y(i))+dxp*dxp*(y(i)-y(im)))/(dxm*dxp*dx0)
792  END DO
793  dydx(1)=2.*dydx(2) - dydx(3)
794  dydx(nx)=2.*dydx(nx-1) - dydx(nx-2)
795 
796  END SUBROUTINE lderiv
797 
798 end module itm_toolbox
l3deriv
subroutine l3deriv(y_in, x_in, nr_in, dydx_out, x_out, nr_out)
Definition: l3interp.f90:59
avint
subroutine avint(ntab, xtab, ytab, a, b, result)
Definition: avint.f90:1
itm_toolbox::rho_with_accumulation
real(r8) function, dimension(1:size(x)) rho_with_accumulation(function_string, x)
itm_toolbox:rho_with_accumulation
Definition: itm_toolbox.F90:55
l3interp
subroutine l3interp(y_in, x_in, nr_in, y_out, x_out, nr_out)
Definition: l3interp.f90:1
itm_toolbox
prototype ITM_TOOLBOX
Definition: itm_toolbox.F90:11
lderiv
subroutine lderiv(y, x, dydx, nx)
Definition: l3interp.f90:121
profile
real(r8) function, dimension(1:size(x)) profile(function_string, x)
Definition: ets_workflow.F90:1934
phys_functions::r2
real(r8) function r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:121
cubint
subroutine cubint(ntab, xtab, ftab, ia_in, ib_in, result, error)
Definition: cubint.f90:1
itm_toolbox::interpolate
real(r8) function, dimension(size(x_out)) interpolate(x_in, y_in, x_out)
&quot;generic&quot; interpolation routine (only calls l3interp at the moment)
Definition: itm_toolbox.F90:29