ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
solution1.f90
Go to the documentation of this file.
1 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
7 ! + + + + + + + + + NUMERICAL SOLUTION + + + + + + + + +
8 
9 SUBROUTINE solution1 (SOLVER, ifail)
10 
11 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
12 ! +
13 ! This subroutine is prepared to solve single +
14 ! transport equation in standardised form +
15 ! adopted by the ETS. +
16 ! +
17 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
18 ! Source: 1-D transport code RITM +
19 ! Developers: M.Tokar, D.Kalupin +
20 ! Kontacts: M.Tokar@fz-juelich.de +
21 ! D.Kalupin@fz-juelich.de +
22 ! +
23 ! Comments: equation is solved in +
24 ! DIFFERENTIAL form +
25 ! +
26 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
27 
28 
29  USE type_solver
30 
31  use itm_types
32 
33  IMPLICIT NONE
34 
35  INTEGER, INTENT (INOUT) :: ifail
36 
37 ! +++ Input/Output with ETS:
38  TYPE (numerics), INTENT (INOUT) :: solver !contains all I/O quantities to numerics part
39 
40 
41 ! +++ Internal input/output parameters:
42  INTEGER :: idim, ndim !equation index and total number of equations to be solved
43 
44  INTEGER :: irho, nrho !radius index, number of radial points
45 
46  INTEGER :: flag !flag for equation: 0 - interpretative (not solved), 1 - predictive (solved)
47 
48  REAL (R8) :: rho(solver%nrho) !radii
49 
50  REAL (R8) :: amix !fraction of new sollution mixed
51 
52  REAL (R8) :: y(solver%nrho), ym(solver%nrho) !function at the current amd previous time steps
53  REAL (R8) :: dy(solver%nrho) !derivative of function
54 
55  REAL (R8) :: a(solver%nrho), b(solver%nrho) !coefficients
56  REAL (R8) :: c(solver%nrho), d(solver%nrho) !coefficients
57  REAL (R8) :: e(solver%nrho), f(solver%nrho) !coefficients
58  REAL (R8) :: g(solver%nrho), h !coefficients
59  REAL (R8) :: dd(solver%nrho), de(solver%nrho) !derivatives of coefficients
60 
61  REAL (R8) :: v(2), u(2), w(2) !boundary conditions
62 
63 
64 ! +++ Coefficients used by internal numerical solver:
65  INTEGER :: i, n !radius index, number of radial points
66  INTEGER :: isp !spline flag
67 
68  REAL (R8) :: x(solver%nrho) !radii
69 
70  REAL (R8) :: sol(solver%nrho) !solution
71  REAL (R8) :: dsol(solver%nrho) !derivative of solution
72 
73  REAL (R8) :: as(solver%nrho), bs(solver%nrho) !coefficients
74  REAL (R8) :: cs(solver%nrho) !coefficients
75 
76  REAL (R8) :: vs(2), us(2), ws(2) !boundary conditions
77 
78 
79 
80 ! + + + + + + + + + INTERFACE PART + + + + + + + + + + +
81 
82 
83 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
84 ! +++ Set up local variables (input, obtained from ETS):
85 ! Control parameters:
86  ndim = solver%NDIM
87  nrho = solver%NRHO
88  amix = solver%AMIX
89 
90 
91 
92 ! +++ Solution of equations starting from 1 to NDIM:
93  equation_loop1: DO idim = 1, ndim
94 
95  flag = solver%EQ_FLAG(idim)
96 
97 
98  IF (flag.EQ.0) goto 20 !equation is not solved
99 
100 
101 
102 ! +++ Numerical coefficients obtained from physics part in form:
103 !
104 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y
105 
106  rho_loop1: DO irho=1,nrho
107  rho(irho) = solver%RHO(irho)
108 
109  y(irho) = solver%Y(idim,irho)
110  dy(irho) = solver%DY(idim,irho)
111  ym(irho) = solver%YM(idim,irho)
112 
113  a(irho) = solver%A(idim,irho)
114  b(irho) = solver%B(idim,irho)
115  c(irho) = solver%C(idim,irho)
116  d(irho) = solver%D(idim,irho)
117  e(irho) = solver%E(idim,irho)
118  f(irho) = solver%F(idim,irho)
119  g(irho) = solver%G(idim,irho)
120 
121  END DO rho_loop1
122 
123  h = solver%H
124 
125  CALL derivn1(nrho,rho,d,dd) !Derivation of coefficient D
126  CALL derivn1(nrho,rho,e,de) !Derivation of coefficient E
127 
128 ! +++ Boundary conditions for numerical solver in form:
129 !
130 ! V*Y' + U*Y = W
131 
132 ! +++ On axis:
133 
134  v(1) = solver%V(idim,1)
135  u(1) = solver%U(idim,1)
136  w(1) = solver%W(idim,1)
137 
138 ! +++ At the edge:
139 
140  v(2) = solver%V(idim,2)
141  u(2) = solver%U(idim,2)
142  w(2) = solver%W(idim,2)
143 
144 
145 
146 
147 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
148 ! +++ Set up numerical coefficients used in the solver
149 ! (translation is done for differential form
150 ! of transport equation):
151 
152 ! indexies:
153  IF(solver%type .EQ. 1) THEN !SOLVER1 (w/o spline)
154  isp = 0
155  ELSE IF(solver%type .EQ. 4) THEN !SOLVER4 (with spline)
156  isp = 1
157  ELSE
158  write(*,*) 'Invalid solver%type = ',solver%type
159  stop
160  END IF
161 
162  n = nrho
163 ! coefficients:
164  rho_loop2: DO i = 1,n
165  x(i) = rho(i)
166  as(i) = -e(i)/d(i) + dd(i)/d(i)
167  bs(i) = a(i)*c(i)/h/d(i) + c(i)*g(i)/d(i) + de(i)/d(i)
168  cs(i) = c(i)/d(i) * (f(i) + b(i)*ym(i)/h)
169  END DO rho_loop2
170 ! boundary conditions:
171  vs(1) = v(1)
172  us(1) = u(1)
173  ws(1) = w(1)
174  vs(2) = v(2)
175  us(2) = u(2)
176  ws(2) = w(2)
177 
178 
179 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
180 ! +++ Solution of transport equation:
181 
182  CALL solver_ritm(x,n, vs,us,ws, as,bs,cs, sol,dsol, isp)
183 ! radii b.c. coeff. solution spline
184 
185 
186  rho_loop3: DO irho = 1,nrho
187 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
188 ! +++ New function and derivative:
189  y(irho) = y(irho)*(1.e0_r8-amix) + sol(irho)*amix
190  dy(irho) = dy(irho)*(1.e0_r8-amix) + dsol(irho)*amix
191 
192 
193 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
194 ! +++ Return solution to ETS:
195  solver%Y(idim,irho) = y(irho)
196  solver%DY(idim,irho) = dy(irho)
197 
198  END DO rho_loop3
199 
200 
201 
202 20 CONTINUE
203 
204 
205  END DO equation_loop1
206 
207 
208 
209  RETURN
210 
211 
212 
213 END SUBROUTINE solution1
214 
215 
216 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
217 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
218 
219 
220 
221 
222 
223 
224 ! + + + + + + + + + NUMERICAL PART + + + + + + + + + + +
225 
226 
227 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
235 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
236 SUBROUTINE solver_ritm(X,N, V,U,W, A,B,C, SOL,DSOL, ISP)
237 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
238 ! This is interface for RITM solver, which allows to !
239 ! use the solver on grids starting from zero !
240 ! !
241 ! The subroutine interploates the function at X=0 !
242 ! by parabola: Y=E+F*X^2
243 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
244 
245  use itm_types
246 
247  IMPLICIT NONE
248 
249  INTEGER, INTENT (INOUT) :: isp
250 
251  INTEGER :: i, n, n1
252 
253  REAL (R8) :: x(n), x1(n-1) !grid starting at x=0 and starting at x=/=0
254 
255  REAL (R8) :: a(n), b(n), c(n) !coefficients on the grid starting at x=0
256  REAL (R8) :: a1(n), b1(n), c1(n) !coefficients on the grid starting at x=/=0
257 
258  REAL (R8) :: v(2), u(2), w(2) !boundary conditions on the grid starting at x=0
259  REAL (R8) :: v1(2), u1(2), w1(2) !boundary conditions on the grid starting at x=/=0
260 
261  REAL (R8) :: sol(n), dsol(n) !solution and its derivative on the grid starting at x=0
262  REAL (R8) :: sol1(n), dsol1(n) !solution and its derivative on the grid starting at x=/=0
263 
264  REAL (R8) :: e, f, k1, k2 !coefficients of parabola
265 
266 
267  IF (x(1).EQ.0.e0_r8.AND.u(1).EQ.0.e0_r8.AND.w(1).EQ.0.e0_r8) THEN
268 
269  n1 = n - 1
270 
271  DO i = 1,n1
272  x1(i) = x(i+1)
273  a1(i) = a(i+1)
274  b1(i) = b(i+1)
275  c1(i) = c(i+1)
276  END DO
277 
278  k1 = 2.e0_r8 + 2.0*a1(1)*x1(1) - b1(1)*x1(1)**2
279  k2 = 2.e0_r8 + 2.0*a1(2)*x1(2) - b1(2)*x1(2)**2
280 
281  f = - (c1(1)*b1(2)-c1(2)*b1(1)) / (k1*b1(2)-k2*b1(1))
282  e = (f*k1+ c1(1))/b1(1)
283 
284  IF (isp.EQ.0) THEN !SOLVER1
285  v1(1) = 1.e0_r8/x1(1) - a1(1)
286  u1(1) = b1(1)
287  w1(1) = c1(1)
288  ELSE IF (isp.EQ.1) THEN !SOLVER4
289  v1(1) = 1.e0_r8
290  u1(1) = 0.e0_r8
291  w1(1) = 1.e0_r8 * f * x1(1)
292  END IF
293 
294  v1(2) = v(2)
295  u1(2) = u(2)
296  w1(2) = w(2)
297 
298 
299  CALL solver_ritm1(x1,n1, v1,u1,w1, a1,b1,c1, sol1,dsol1, isp)
300 ! radii b.c. coeff. solution spline
301 
302  DO i = 2,n
303  sol(i) = sol1(i-1)
304  dsol(i) = dsol1(i-1)
305  END DO
306 
307 ! SOL(1) = E
308 
309 
310  IF (isp.EQ.0) THEN !SOLVER1
311  CALL axis1(n, x, sol)
312  ELSE IF (isp.EQ.1) THEN !SOLVER4
313  sol(1) = e
314  END IF
315 
316  dsol(1) = 0.e0_r8
317 
318  ELSE IF (x(1).NE.0.e0_r8) THEN
319 
320  CALL solver_ritm1(x,n, v,u,w, a,b,c, sol,dsol, isp)
321 ! radii b.c. coeff. solution spline
322 
323  END IF
324 
325 
326 
327  RETURN
328 
329 
330 
331 END SUBROUTINE solver_ritm
332 
333 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
334 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
335 
336 
337 
338 
339 
340 
341 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
350 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
351 SUBROUTINE solver_ritm1(x,n, v,u,w, a,b,f, y,d, isp)
352 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
353 ! This subroutine provides solution for the equation !
354 ! !
355 ! Y''(x)+a(x)*Y'(x)=b(x)*Y(x)-f(x) !
356 ! !
357 ! with boundary conditions: v(1)*Y'(1)+u(1)*Y(1)=w(1) !
358 ! v(2)*Y'(n)+u(2)*Y(n)=w(2) !
359 ! !
360 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
361 
362  use itm_types
363 
364  IMPLICIT NONE
365 
366 ! +++ Input:
367  INTEGER :: n !number of radial points
368  INTEGER :: isp !flag for use of spline, default: isp=0
369 
370  REAL (R8) :: x(n) !radial coordinate
371  REAL (R8) :: a(n),b(n),f(n) !coefficients
372  REAL (R8) :: u(2),v(2),w(2) !boundary conditions
373 
374 ! +++ Output:
375  REAL (R8) :: y(n),d(n) !function and its derivative
376 
377 ! +++ Internal variables:
378  INTEGER :: i,l,m,mn,k,j,jm,ii
379 
380  REAL (R8) :: h(0:n+1)
381  REAL (R8) :: lam
382  REAL (R8) :: y0(n,3),d0(n,3),y1(n,2,3),d1(n,2,3),s1(n,2,3)
383  REAL (R8) :: cs(n,2,2),zt(n,2),et(n,2,2,2),mu(n,2,2),c(n,2)
384  REAL (R8) :: asi,arg,ay,cy,by,dy,bf,df,s,si,de,xi,af,cf,det
385 
386 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
387  h(0)=0.e0_r8
388 ! h(1)=x(1)
389  DO i=1,n-1
390  h(i)=x(i+1)-x(i)
391  END DO
392  h(n)=0.e0_r8
393  h(n+1)=0.e0_r8
394 
395 
396 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
397 ! +++ partial solution
398  DO i=1,n
399  DO l=1,3
400  y0(i,l)=f(i)/b(i)
401  d0(i,l)=0.e0_r8
402  END DO
403  END DO
404  IF(isp.EQ.0) goto 10
405  DO i=2,n-1
406 
407 ! +++ hyperbolic-trigonometric spline for f(x)
408  CALL hyp_trig_spline1(i,n,x,f,af,bf,cf,df,si)
409  jm=1
410  IF(i.EQ.2.OR.i.EQ.n-1) jm=2
411  DO j=1,jm
412  ii=i
413  IF(i.EQ.2.AND.j.EQ.2) ii=1
414  IF(i.EQ.n-1.AND.j.EQ.2) ii=n
415  asi=a(ii)*si
416 
417 ! +++ weights of trigonometric functions
418  arg=b(ii)+si**2
419  ay=(af*arg+cf*asi)/(arg**2+asi**2)
420  cy=(cf*arg-af*asi)/(arg**2+asi**2)
421 
422 ! +++ weights of hyperbolic functions
423  arg=b(ii)-si**2
424  by=(bf*arg+df*asi)/(arg**2-asi**2)
425  dy=(df*arg+bf*asi)/(arg**2-asi**2)
426  DO l=1,2
427  s=(-1)**l*si*h(i+l-2)/2.e0_r8
428  IF(i.EQ.2.AND.j.EQ.2) THEN
429  IF(l.EQ.1) s=-si*h(1)
430  IF(l.EQ.2) s=-si*h(1)/2.e0_r8
431  END IF
432  IF(i.EQ.n-1.AND.j.EQ.2) THEN
433  IF(l.EQ.1) s=si*h(n-1)/2.e0_r8
434  IF(l.EQ.2) s=si*h(n-1)
435  END IF
436  y0(ii,l)=ay*cos(s)+cy*sin(s)+ &
437  by*cosh(s)+dy*sinh(s)
438  d0(ii,l)=si*(-ay*sin(s)+cy*cos(s)+ &
439  by*sinh(s)+dy*cosh(s))
440  END DO
441  IF(j.EQ.1) y0(i,3)=ay+by
442  END DO
443  IF(i.EQ.2) y0(1,3)=y0(1,1)
444  IF(i.EQ.n-1) y0(n,3)=y0(n,2)
445  END DO
446 
447 
448 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
449 ! +++ general solution
450 10 DO i=1,n
451  de=a(i)**2/4.e0_r8+b(i)
452  DO k=1,2
453  lam=(-1)**k*sqrt(max(0.e0_r8,de))-a(i)/2.e0_r8
454  IF(abs(b(i)).LT.1.e-10_r8*a(i)**2/4.e0_r8) &
455  lam=-a(i)*(2.e0_r8-k)+b(i)/a(i)*(k-1.e0_r8)
456  DO l=1,3
457  s=(-1)**l*h(i+l-2)/2.e0_r8
458  IF(i.EQ.1.AND.l.EQ.1.OR.i.EQ.n.AND.l.EQ.2.OR.l.EQ.3) &
459  s=0.e0_r8
460  y1(i,k,l)=exp(s*lam-(h(i-1)*(abs(lam)-lam)+ &
461  h(i)*(abs(lam)+lam))/4.e0_r8)/(1.e0_r8+(sqrt(max &
462  (0.e0_r8,de))+abs(a(i)/2.e0_r8))*(h(i-1)+h(i))/2.e0_r8)
463  s1(i,k,l)=lam
464  IF(de.EQ.0.e0_r8) THEN
465  xi=(k-1.e0_r8)*(x(i)+s)
466  y1(i,k,l)=y1(i,k,l)*(1.e0_r8+xi/(abs(x(1))+abs(x(n))))
467  s1(i,k,l)=lam+(k-1.e0_r8)/(abs(x(1))+abs(x(n))+xi)
468  END IF
469  IF(de.LT.0.e0_r8) THEN
470  xi=(-1)**k*tan(sqrt(-de)*s)
471  y1(i,k,l)=y1(i,k,l)*cos(sqrt(-de)*s)*(1.e0_r8+xi)
472  s1(i,k,l)=lam+sqrt(-de)*(-1)**k*(1.e0_r8-xi)/(1.e0_r8+xi)
473  END IF
474  d1(i,k,l)=y1(i,k,l)*s1(i,k,l)
475  END DO
476  END DO
477  END DO
478 
479 
480 
481 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
482 ! +++ total solution
483  DO m=1,2
484  mn=n*(m-1)+2-m
485  DO k=1,2
486  cs(mn,k,m)=u(m)+v(m)*s1(mn,k,m)
487  END DO
488  zt(mn,m)=(w(m)-u(m)*y0(mn,m)-v(m)*d0(mn,m))
489  det=max(1.e0_r8,abs(cs(mn,1,m)),abs(cs(mn,2,m)))
490  IF(det.GT.0.e0_r8) THEN
491  DO k=1,2
492  cs(mn,k,m)=cs(mn,k,m)/det
493  END DO
494  zt(mn,m)=zt(mn,m)/det
495  END IF
496  DO i=mn,(n-1)*(2-m)+2*(m-1),-(-1)**m
497  j=i-(-1)**m
498  DO l=1,2
499  DO k=1,2
500  et(i,k,l,m)=(-1)**l*(s1(j,k,m)- &
501  s1(i,3-l,3-m))/(s1(i,2,3-m)-s1(i,1,3-m))
502  END DO
503  mu(i,l,m)=(-1)**l*(d0(j,m)-d0(i,3-m)-(y0(j,m)-y0(i,3-m))* &
504  s1(i,3-l,3-m))/(s1(i,2,3-m)-s1(i,1,3-m))
505  END DO
506  DO k=1,2
507  cs(j,k,m)=et(i,k,1,m)*y1(i,1,m)*y1(i,2,3-m)* &
508  cs(i,1,m)+et(i,k,2,m)*y1(i,1,3-m)*y1(i,2,m)*cs(i,2,m)
509  END DO
510  zt(j,m)=(zt(i,m)*y1(i,1,3-m)*y1(i,2,3-m)- &
511  mu(i,1,m)*y1(i,1,m)*y1(i,2,3-m)*cs(i,1,m)- &
512  mu(i,2,m)*y1(i,1,3-m)*y1(i,2,m)*cs(i,2,m))
513  det=max(abs(cs(j,1,m)),abs(cs(j,2,m)))
514  IF(det.GT.0.e0_r8) THEN
515  DO k=1,2
516  cs(j,k,m)=cs(j,k,m)/det
517  END DO
518  zt(j,m)=zt(j,m)/det
519  END IF
520  END DO
521  END DO
522 
523  DO i=1,n
524  DO k=1,2
525  j=3-k
526  c(i,k)=(zt(i,k)*y1(i,j,j)*cs(i,j,j)- &
527  zt(i,j)*y1(i,j,k)*cs(i,j,k))/ &
528  (y1(i,1,1)*y1(i,2,2)*cs(i,1,1)*cs(i,2,2)- &
529  y1(i,1,2)*y1(i,2,1)*cs(i,1,2)*cs(i,2,1))
530  END DO
531  y(i)=c(i,1)*y1(i,1,3)+c(i,2)*y1(i,2,3)+y0(i,3)
532  d(i)=c(i,1)*d1(i,1,3)+c(i,2)*d1(i,2,3)+d0(i,3)
533  END DO
534 
535 
536 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
537 ! +++ splined derivative:
538  DO i=2,n-1
539  d(i)=((y(i)-y(i-1))*h(i)/h(i-1)+ &
540  (y(i+1)-y(i))*h(i-1)/h(i))/(h(i-1)+h(i))
541  END DO
542 ! d(1)=((y(2)-y(1))*(h(2)/h(1)+2.e0_R8)-(y(3)-y(2))*h(1)/h(2))/ &
543 ! (h(1)+h(2))
544  d(n)=((y(n)-y(n-1))*(h(n-2)/h(n-1)+2.e0_r8)- &
545  (y(n-1)-y(n-2))*h(n-1)/h(n-2))/(h(n-2)+h(n-1))
546 
547 
548  RETURN
549 END SUBROUTINE solver_ritm1
550 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
551 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
552 
553 
554 
555 
556 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
564 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
565 SUBROUTINE hyp_trig_spline1(i,n,x,f,a,b,c,d,l)
566 
567  use itm_types
568 
569  IMPLICIT NONE
570 
571  INTEGER :: i,n
572 
573  REAL (R8) :: x(n), f(n)
574  REAL (R8) :: a, b, c, d, l
575  REAL (R8) :: l1,s1,s2,cs,ga
576 
577  s1=2.e0_r8*(x(i)-x(i-1))/(x(i+1)-x(i-1))
578  s2=2.e0_r8-s1
579 
580  IF(f(i).EQ.0.e0_r8) goto 30
581  cs=(f(i+1)*s1+f(i-1)*s2)/f(i)/2.e0_r8
582  l=1.e-5_r8
583 
584  IF(1.e0_r8+1.e-7_r8.LT.cs.AND.cs.LE.2.e0_r8) THEN
585 10 l1=l
586  ga=cosh(l)*(f(i-1)/sinh(l*s1)+f(i+1)/sinh(l*s2))/ &
587  (f(i)/tanh(l*s1)+f(i)/tanh(l*s2))
588  l=log(ga+sqrt(ga*ga-1.e0_r8))
589  IF(abs(l).GE.1.e0_r8) goto 30
590  IF(s1.NE.s2.AND.abs(1.e0_r8-l1/l).GT.1.e-5_r8) goto 10
591  a=0.e0_r8
592  b=f(i)
593  c=0.e0_r8
594  d=(f(i+1)*cosh(l*s1)-f(i-1)*cosh(l*s2))/sinh(2.e0_r8*l)
595  goto 40
596  END IF
597 
598  IF(-.5e0_r8.LE.cs.AND.cs.LT.1.e0_r8-1.e-7_r8) THEN
599 20 l1=l
600  ga=(f(i-1)*sin(l*s2)+f(i+1)*sin(l*s1))/2.e0_r8/f(i)/sin(l)
601  l=acos(ga)
602  IF(abs(l).GE.1.e0_r8) goto 30
603  IF(s1.NE.s2.AND.abs(1.e0_r8-l1/l).GT.1.e-5_r8) goto 20
604  a=f(i)
605  b=0.e0_r8
606  c=(f(i+1)*cos(l*s1)-f(i-1)*cos(l*s2))/sin(2.e0_r8*l)
607  d=0.e0_r8
608  goto 40
609  END IF
610 
611 30 l=1.e0_r8
612  a=(f(i+1)*(sin(s1)+sinh(s1))+f(i-1)*(sin(s2)+sinh(s2))- &
613  f(i)*(sin(s1)*cosh(s2)+sin(s2)*cosh(s1)+sinh(2.e0_r8)))/ &
614  (sin(2.e0_r8)+cos(s2)*sinh(s1)+cos(s1)*sinh(s2)- &
615  sin(s1)*cosh(s2)-sin(s2)*cosh(s1)-sinh(2.e0_r8))
616  b=f(i)-a
617  c=(f(i+1)-f(i-1)+a*(cos(s1)-cos(s2))+b*(cosh(s1)- &
618  cosh(s2)))/(sin(s1)+sinh(s1)+sin(s2)+sinh(s2))
619  d=c
620 
621 40 l=2.e0_r8*l/(x(i+1)-x(i-1))
622 
623 
624  RETURN
625 END SUBROUTINE hyp_trig_spline1
626 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
627 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
628 
629 
630 
631 
632 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
639 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
640 SUBROUTINE derivn1(N,X,Y,DY1)
641 
642  use itm_types
643 
644  IMPLICIT NONE
645 
646  INTEGER :: n ! number of radial points (input)
647  INTEGER :: i
648 
649  REAL (R8) :: x(n), & ! argument array (input)
650  y(n), & ! function array (input)
651  dy1(n) ! function derivative array (output)
652  REAL (R8) :: h(n),dy2(n)
653 
654  DO i=1,n-1
655  h(i)=x(i+1)-x(i)
656  END DO
657 
658  DO i=2,n-1
659  dy1(i)=((y(i+1)-y(i))*h(i-1)/h(i)+(y(i)-y(i-1))*h(i)/h(i-1)) &
660  /(h(i)+h(i-1))
661  dy2(i)=2.e0_r8*((y(i-1)-y(i))/h(i-1)+(y(i+1)-y(i))/h(i)) &
662  /(h(i)+h(i-1))
663  END DO
664 
665  dy1(1)=dy1(2)-dy2(2)*h(1)
666  dy1(n)=dy1(n-1)+dy2(n-1)*h(n-1)
667 
668  RETURN
669 END SUBROUTINE derivn1
670 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
671 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
672 
673 
674 
675 
676 
677 
678 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
684 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
685  SUBROUTINE axis1(n, r, f)
686 !-------------------------------------------------------!
687 ! !
688 ! This subroutine finds !
689 ! f(r_1=0) from f(r_2), f(r_3) and f(r_4) !
690 ! !
691 !-------------------------------------------------------!
692 
693  USE itm_types
694 
695  IMPLICIT NONE
696 
697  INTEGER n, i
698  REAL *8 h(n), r(n), f(n)
699 
700 
701  h(1)=r(1)
702  DO i=2,n-1
703  h(i)=r(i+1)-r(i)
704  END DO
705  h(n)=0.e0_r8
706 
707 
708  f(1) = ((f(2)*r(4)/h(2)+f(4)*r(2)/h(3))*r(3) &
709  -f(3)*(r(2)/h(2)+r(4)/h(3))*r(2)*r(4)/r(3)) &
710  /(r(4)-r(2))
711 
712 
713  RETURN
714 
715 
716  END SUBROUTINE axis1
717 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
718 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
solver
Definition: solver.f90:1
type_solver
Definition: type_solver.f90:1
solution1
subroutine solution1(SOLVER, ifail)
This subroutine is prepared to solve single transport equation in standardised form adopted by the ET...
Definition: solution1.f90:9
derivn1
subroutine derivn1(N, X, Y, DY1)
These subroutines calculate first and second derivatives, DY1 and DY2, of function Y respect to argum...
Definition: solution1.f90:640
solver_ritm1
subroutine solver_ritm1(x, n, v, u, w, a, b, f, y, d, isp)
This subroutine provides solution for the equation Y''(x)+a(x)*Y'(x)=b(x)*Y(x)-f(x) with boundary con...
Definition: solution1.f90:351
hyp_trig_spline1
subroutine hyp_trig_spline1(i, n, x, f, a, b, c, d, l)
This subroutine interpolates the experimental data by using qubic spline so that k additional points ...
Definition: solution1.f90:565
type_solver::numerics
Definition: type_solver.f90:22
solver_ritm
subroutine solver_ritm(X, N, V, U, W, A, B, C, SOL, DSOL, ISP)
This is interface for RITM solver, which allows to use the solver on grids starting from zero...
Definition: solution1.f90:236
axis1
subroutine axis1(n, r, f)
This subroutine finds f(r_1=0) from f(r_2), f(r_3) and f(r_4)
Definition: solution1.f90:685