Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2011-2019 Bill Paxton & The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module ion_tables_eval
21 :
22 : use const_def, only: dp, one_sixth, arg_not_provided
23 : use ionization_def
24 : use math_lib
25 : use utils_lib, only: mesa_error
26 :
27 : implicit none
28 :
29 : logical, parameter :: dbg = .false.
30 :
31 : contains
32 :
33 4 : subroutine Get_ion_Results(&
34 : Z_in, X_in, arho, alogrho, atemp, alogtemp, &
35 : res, ierr)
36 : use const_def, only: dp
37 : use utils_lib, only: is_bad
38 : use ion_tables_load, only: Load_ion_Table
39 :
40 : ! INPUT
41 :
42 : real(dp), intent(in) :: Z_in ! the desired Z
43 : real(dp), intent(in) :: X_in ! the desired X
44 :
45 : real(dp), intent(in) :: arho, alogrho ! the density
46 : ! provide both if you have them. else pass one and set the other to = arg_not_provided
47 :
48 : real(dp), intent(in) :: atemp, alogtemp ! the temperature
49 : ! provide both if you have them. else pass one and set the other to = arg_not_provided
50 :
51 :
52 : ! OUTPUT
53 :
54 : real(dp), intent(inout) :: res(num_ion_vals)
55 : integer, intent(out) :: ierr
56 :
57 : real(dp), dimension(num_ion_vals) :: res1, d_dlnRho_c_T1, d_dlnT_c_Rho1
58 : real(dp), dimension(num_ion_vals) :: res2, d_dlnRho_c_T2, d_dlnT_c_Rho2
59 2 : real(dp) :: X, Z, Rho, logRho, T, logT
60 : integer :: iregion
61 : character (len=256) :: message
62 :
63 : real(dp) :: alfa, beta, c_dx, c_dy, logRho_lo, logRho_hi, &
64 : logT1, logT2, logT7, logT8, logRho3, logRho4
65 : real(dp) :: logQ, A, B, sin_pi_alfa, dA_dlnT, dA_dlnRho, dB_dlnT, dB_dlnRho
66 : real(dp) :: d_dx_dlogT, d_dx_dlogRho, d_dy_dlogT, d_dy_dlogRho
67 : real(dp) ::&
68 : lnPgas, Pgas, dlnPgas_dlnRho, dlnPgas_dlnT, dPgas_dlnRho, dPgas_dlnT, &
69 : Prad, dPrad_dlnT, P, dP_dlnRho, dP_dlnT, dlnP_dlnRho, dlnP_dlnT, &
70 : lnE, energy, dlnE_dlnRho, dlnE_dlnT, &
71 : lnS, entropy, dlnS_dlnRho, dlnS_dlnT, &
72 : d_alfa_dlogT, d_alfa_dlogRho, d_beta_dlogT, d_beta_dlogRho
73 : real(dp), parameter :: dZ_transition = 0.01d0 ! must be > 0
74 : real(dp), parameter :: logT_margin = 0.1d0
75 : real(dp), parameter :: tiny = 1d-20
76 :
77 : logical :: debug
78 :
79 : include 'formats'
80 :
81 2 : ierr = 0
82 2 : debug = dbg
83 :
84 2 : if (.not. ion_is_initialized) then
85 2 : call Load_ion_Table(ierr)
86 2 : if (ierr /= 0) then
87 : if (dbg) write(*,*) 'Load_ion_Table ierr', ierr
88 0 : return
89 : end if
90 : end if
91 :
92 2 : if (is_bad(X_in) .or. is_bad(Z_in)) then
93 0 : ierr = -1
94 0 : return
95 : end if
96 :
97 2 : X = X_in; Z = Z_in
98 2 : if (X < tiny) X = 0
99 2 : if (Z < tiny) Z = 0
100 :
101 : !..get temp and rho args
102 2 : T = atemp; logT = alogtemp
103 2 : if (atemp == arg_not_provided .and. alogtemp == arg_not_provided) then
104 0 : ierr = -2; return
105 : end if
106 2 : if (alogtemp == arg_not_provided) logT = log10(T)
107 2 : if (atemp == arg_not_provided) T = exp10(logT)
108 :
109 2 : if (T <= 0) then
110 0 : ierr = -1
111 0 : return
112 : end if
113 :
114 2 : Rho = arho; logrho = alogrho
115 2 : if (arho == arg_not_provided .and. alogrho == arg_not_provided) then
116 0 : ierr = -3; return
117 : end if
118 2 : if (alogrho == arg_not_provided) logRho = log10(Rho)
119 2 : if (arho == arg_not_provided) Rho = exp10(logRho)
120 :
121 2 : if (Rho <= 0) then
122 0 : ierr = -1
123 0 : return
124 : end if
125 :
126 2 : if (is_bad(Rho) .or. is_bad(T)) then
127 0 : ierr = -1
128 0 : return
129 : end if
130 2 : call Get_ion_ZResults(Z, X, Rho, logRho, T, logT, res, ierr)
131 :
132 : end subroutine Get_ion_Results
133 :
134 :
135 2 : subroutine Get_ion_ZResults(&
136 : Z, X, Rho, logRho, T, logT, &
137 : res, ierr)
138 : use chem_def
139 : real(dp), intent(in) :: Z, X, Rho, logRho, T, logT
140 : real(dp), intent(inout) :: res(num_ion_vals)
141 : integer, intent(out) :: ierr
142 :
143 : real(dp), dimension(num_ion_vals, num_ion_Zs) :: &
144 292 : res_zx, d_dlnRho_c_T_zx, d_dlnT_c_Rho_zx, d_res_dX
145 : real(dp), dimension(num_ion_vals) :: d_dX, d_dZ
146 8 : real(dp) :: d_res_dZ(num_ion_vals), dZ, denom, c(3),&
147 : dP_dZ, dlnP_dZ
148 : real(dp), parameter :: tiny = 1d-6
149 :
150 : character (len=256) :: message
151 : integer :: iz, j, ci
152 : logical, parameter :: ion_dbg = dbg
153 :
154 2 : ierr = 0
155 :
156 : if (num_ion_Zs < 3) then
157 : write(*, *) 'error: Get_ion_ZResults assumes num_ion_Zs >= 3'
158 : call mesa_error(__FILE__,__LINE__)
159 : end if
160 :
161 : if (ion_Zs(1) /= 0) then
162 : write(*, *) 'error: Get_ion_ZResults assumes ion_Zs(1) == 0'
163 : call mesa_error(__FILE__,__LINE__)
164 : end if
165 :
166 : if (abs(ion_Zs(1) - 2*ion_Zs(2) + ion_Zs(3)) > tiny) then
167 : write(*, *) 'error: Get_ion_ZResults assumes equal spaced Zs(1:3)'
168 : call mesa_error(__FILE__,__LINE__)
169 : end if
170 :
171 2 : if (Z < tiny) then
172 0 : call Get_ion_for_X(1, X, Rho, logRho, T, logT, res, ierr)
173 0 : return
174 : end if
175 :
176 2 : if (Z > ion_Zs(3)) then
177 :
178 0 : if (Z <= ion_Zs(4)) then
179 0 : call do_interp2(3,4,ierr)
180 0 : if (ierr /= 0) return
181 0 : else if (Z <= ion_Zs(5) - tiny) then
182 0 : call do_interp2(4,5,ierr)
183 0 : if (ierr /= 0) return
184 : else
185 0 : call Get_ion_for_X(5, X, Rho, logRho, T, logT, res, ierr)
186 0 : return
187 : end if
188 :
189 : else
190 :
191 8 : do iz = 1, 3
192 6 : call Get_ion_for_X(iz, X, Rho, logRho, T, logT, res_zx(:, iz), ierr)
193 8 : if (ierr /= 0) return
194 : end do
195 :
196 2 : dZ = ion_Zs(2) - ion_Zs(1)
197 2 : denom = 2*dZ*dZ
198 2 : c(1) = (2*dZ*dZ - 3*dZ*Z + Z*Z)/denom
199 2 : c(2) = 2*(2*dZ-Z)*Z/denom
200 2 : c(3) = Z*(Z-dZ)/denom
201 :
202 : end if
203 :
204 58 : res(:) = c(1)*res_zx(:, 1) + c(2)*res_zx(:, 2) + c(3)*res_zx(:, 3)
205 :
206 :
207 : contains
208 :
209 0 : subroutine do_interp2(iz1, iz2, ierr)
210 2 : use const_def, only: pi
211 : integer, intent(in) :: iz1, iz2
212 : integer, intent(out) :: ierr
213 0 : real(dp) :: Z1, Z2, alfa, beta
214 0 : ierr = 0
215 0 : Z1 = ion_Zs(iz1)
216 0 : Z2 = ion_Zs(iz2)
217 0 : alfa = (Z - Z1) / (Z2 - Z1)
218 0 : alfa = (1 - cospi(alfa))/2 ! smooth the transitions
219 0 : beta = 1 - alfa
220 0 : call Get_ion_for_X(iz1, X, Rho, logRho, T, logT, res_zx(:, 1), ierr)
221 0 : if (ierr /= 0) return
222 0 : call Get_ion_for_X(iz2, X, Rho, logRho, T, logT, res_zx(:, 2), ierr)
223 0 : if (ierr /= 0) return
224 0 : c(1) = beta
225 0 : c(2) = alfa
226 0 : c(3) = 0
227 0 : res_zx(:,3) = 0
228 : end subroutine do_interp2
229 :
230 : end subroutine Get_ion_ZResults
231 :
232 :
233 6 : subroutine Get_ion_for_X(iz, X, Rho, logRho, T, logT,res, ierr)
234 : integer, intent(in) :: iz ! the index in ion_Zs
235 : real(dp), intent(in) :: X, Rho, logRho, T, logT
236 : real(dp), intent(inout) :: res(num_ion_vals)
237 : integer, intent(out) :: ierr
238 :
239 702 : real(dp), dimension(num_ion_vals, 4) :: res_zx
240 6 : real(dp) :: dX, c(4), denom, delX, coef
241 : character (len=256) :: message
242 : integer :: ix, ix_lo, ix_hi, j, num_Xs
243 : real(dp), parameter :: tiny = 1d-6
244 : logical, parameter :: dbg_for_X = dbg ! .or. .true.
245 :
246 6 : ierr = 0
247 :
248 : if (num_ion_Xs /= 6) then
249 : write(*, *) 'error: Get_ion_for_X assumes num_ion_Xs == 6'
250 : call mesa_error(__FILE__,__LINE__)
251 : end if
252 :
253 : if (ion_Xs(1) /= 0) then
254 : write(*, *) 'error: Get_ion_for_X assumes ion_Xs(1) == 0'
255 : call mesa_error(__FILE__,__LINE__)
256 : end if
257 :
258 6 : num_Xs = num_ion_Xs_for_Z(iz)
259 :
260 6 : if (X < tiny .or. num_Xs == 1) then
261 0 : call Get_ion_XTable_Results(1, iz, Rho, logRho, T, logT, res, ierr)
262 0 : return
263 : end if
264 :
265 6 : dX = ion_Xs(2)-ion_Xs(1)
266 :
267 26 : do ix = 3, num_Xs
268 26 : if (abs(dX - (ion_Xs(ix) - ion_Xs(ix-1))) > tiny) then
269 0 : write(*, *) 'error: Get_ion_for_X assumes equal spaced Xs'
270 0 : call mesa_error(__FILE__,__LINE__)
271 : end if
272 : end do
273 :
274 6 : ix_hi = -1
275 6 : if (X <= ion_Xs(2)) then
276 0 : ix_lo = 1; ix_hi = 3
277 6 : else if (X >= ion_Xs(num_Xs-1)) then
278 4 : ix_lo = num_Xs-2; ix_hi = num_Xs
279 : else
280 6 : do ix = 3, num_Xs-1
281 6 : if (X <= ion_Xs(ix)) then
282 2 : ix_lo = ix-2; ix_hi = ix+1; exit
283 : end if
284 : end do
285 : end if
286 :
287 6 : if (ix_hi < 0) then
288 0 : write(*, *) 'X', X
289 0 : write(*, *) 'ix_lo', ix_lo
290 0 : write(*, *) 'ix_hi', ix_hi
291 0 : write(*, *) 'error: Get_ion_for_X logic bug'
292 0 : call mesa_error(__FILE__,__LINE__)
293 : end if
294 :
295 : if (dbg_for_X) then
296 : write(*, *) 'X', X
297 : write(*, *) 'ix_lo', ix_lo
298 : write(*, *) 'ix_hi', ix_hi
299 : end if
300 :
301 26 : do ix=ix_lo, ix_hi
302 20 : j = ix-ix_lo+1
303 : call Get_ion_XTable_Results(ix, iz, Rho, logRho, T, logT, &
304 20 : res_zx(:, j), ierr)
305 26 : if (ierr /= 0) return
306 : end do
307 :
308 6 : delX = X - ion_Xs(ix_lo)
309 6 : if (ix_hi-ix_lo==2) then
310 :
311 4 : denom = 2*dX*dX
312 4 : c(1) = (2*dX*dX - 3*dX*delX + delX*delX)/denom
313 4 : c(2) = 2*(2*dX-delX)*delX/denom
314 4 : c(3) = delX*(delX-dX)/denom
315 116 : res(:) = c(1)*res_zx(:, 1) + c(2)*res_zx(:, 2) + c(3)*res_zx(:, 3)
316 :
317 : 1 format(a30, e25.15)
318 : if (dbg_for_X) then
319 : end if
320 :
321 : else
322 :
323 2 : coef = (X-ion_Xs(ix_lo+1))/dX
324 : ! coef = fractional location of X between 2nd and 3rd X's for fit.
325 : ! coef is weight for the quadratic based on points 2, 3, 4 of fit.
326 : ! (1-coef) is weight for quadratic based on points 1, 2, 3 of fit.
327 2 : if (coef < 0 .or. coef > 1) then
328 0 : write(*,*) 'logic bug in Get_ion_for_X'
329 0 : call mesa_error(__FILE__,__LINE__)
330 : end if
331 2 : c(1) = -coef*(coef-1)*(coef-1)/2
332 2 : c(2) = (2 + coef*coef*(-5 + 3*coef))/2
333 2 : c(3) = (coef + coef*coef*(4 - 3*coef))/2
334 2 : c(4) = coef*coef*(coef-1)/2
335 : res(:) = c(1)*res_zx(:, 1) + c(2)*res_zx(:, 2) &
336 58 : + c(3)*res_zx(:, 3) + c(4)*res_zx(:, 4)
337 :
338 : end if
339 :
340 :
341 : end subroutine Get_ion_for_X
342 :
343 :
344 40 : subroutine Get_ion_XTable_Results(ix, iz, Rho, logRho_in, T, logT_in, &
345 : res, ierr)
346 : integer, intent(in) :: ix, iz
347 : real(dp), intent(in) :: Rho, logRho_in, T, logT_in
348 : real(dp), intent(inout) :: res(num_ion_vals)
349 : integer, intent(out) :: ierr
350 :
351 20 : real(dp) :: logQ0, logQ1, logT0, logT1
352 : integer :: iQ, jtemp
353 :
354 20 : real(dp) :: logRho, logT, logQ
355 :
356 : include 'formats'
357 :
358 20 : logRho = logRho_in
359 20 : logT = logT_in
360 20 : logQ = logRho - 2*logT + 12
361 :
362 : ierr = 0
363 :
364 20 : call Locate_logQ(logQ, iQ, logQ0, logQ1, ierr)
365 20 : if (ierr /= 0) return
366 :
367 20 : call Locate_logT(logT, jtemp, logT0, logT1, ierr)
368 20 : if (ierr /= 0) return
369 :
370 : call Do_ion_Interpolations(&
371 : ion_num_logQs, ion_logQs, ion_num_logTs, ion_logTs, &
372 : ion_tbl(:, :, :, :, ix, iz), &
373 : iQ, jtemp, logQ0, logQ, logQ1, logT0, logT, logT1, &
374 20 : res, ierr)
375 :
376 :
377 : end subroutine Get_ion_XTable_Results
378 :
379 :
380 20 : subroutine Locate_logQ(logQ, iQ, logQ0, logQ1, ierr)
381 : real(dp), intent(inout) :: logQ
382 : integer, intent(out) :: iQ
383 : real(dp), intent(out) :: logQ0, logQ1
384 : integer, intent(out) :: ierr
385 :
386 20 : ierr = 0
387 20 : iQ = int((logQ - ion_logQ_min) / ion_del_logQ + 1d-4) + 1
388 :
389 20 : if (iQ < 1 .or. iQ >= ion_num_logQs) then
390 :
391 0 : if (iQ < 1) then
392 0 : iQ = 1
393 0 : logQ0 = ion_logQ_min
394 0 : logQ1 = logQ0 + ion_del_logQ
395 0 : logQ = logQ0
396 : else
397 0 : iQ = ion_num_logQs-1
398 0 : logQ0 = ion_logQ_min + (iQ-1) * ion_del_logQ
399 0 : logQ1 = logQ0 + ion_del_logQ
400 0 : logQ = logQ1
401 : end if
402 :
403 : else
404 :
405 20 : logQ0 = ion_logQ_min + (iQ-1) * ion_del_logQ
406 20 : logQ1 = logQ0 + ion_del_logQ
407 :
408 : end if
409 :
410 20 : end subroutine Locate_logQ
411 :
412 :
413 20 : subroutine Locate_logT(logT, iT, logT0, logT1, ierr)
414 : real(dp), intent(inout) :: logT
415 : integer, intent(out) :: iT
416 : real(dp), intent(out) :: logT0, logT1
417 : integer, intent(out) :: ierr
418 :
419 20 : ierr = 0
420 20 : iT = int((logT - ion_logT_min) / ion_del_logT + 1d-4) + 1
421 :
422 20 : if (iT < 1 .or. iT >= ion_num_logTs) then
423 :
424 0 : if (iT < 1) then
425 0 : iT = 1
426 0 : logT0 = ion_logT_min
427 0 : logT1 = logT0 + ion_del_logT
428 0 : logT = logT0
429 : else
430 0 : iT = ion_num_logTs-1
431 0 : logT0 = ion_logT_min + (iT-1) * ion_del_logT
432 0 : logT1 = logT0 + ion_del_logT
433 0 : logT = logT1
434 : end if
435 :
436 : else
437 :
438 20 : logT0 = ion_logT_min + (iT-1) * ion_del_logT
439 20 : logT1 = logT0 + ion_del_logT
440 :
441 : end if
442 :
443 20 : end subroutine Locate_logT
444 :
445 :
446 20 : subroutine Do_ion_Interpolations(&
447 20 : nx, x, ny, y, fin, i, j, &
448 : x0, xget, x1, y0, yget, y1, &
449 : res, ierr)
450 : ! > fval, df_dx, df_dy, ierr)
451 : use interp_2d_lib_sg, only: interp_evbipm_sg
452 : integer, intent(in) :: nx, ny
453 : real(dp), intent(in) :: x(nx), y(ny), fin(4,num_ion_vals,nx,ny)
454 : integer, intent(in) :: i, j ! target cell in f
455 : real(dp), intent(in) :: x0, xget, x1 ! x0 <= xget <= x1; x0 = xs(i), x1 = xs(i+1)
456 : real(dp), intent(in) :: y0, yget, y1 ! y0 <= yget <= y1; y0 = ys(j), y1 = ys(j+1)
457 : real(dp), intent(inout) :: res(num_ion_vals)
458 : integer, intent(out) :: ierr
459 :
460 : real(dp), parameter :: z36th = 1d0/36d0
461 20 : real(dp) :: xp, xpi, xp2, xpi2, ax, axbar, bx, bxbar, cx, cxi, hx2, cxd, cxdi, hx, hxi
462 20 : real(dp) :: yp, ypi, yp2, ypi2, ay, aybar, by, bybar, cy, cyi, hy2, cyd, cydi, hy, hyi
463 20 : real(dp) :: sixth_hx2, sixth_hy2, z36th_hx2_hy2
464 20 : real(dp) :: sixth_hx, sixth_hxi_hy2, z36th_hx_hy2
465 20 : real(dp) :: sixth_hx2_hyi, sixth_hy, z36th_hx2_hy
466 : integer :: k, ip1, jp1
467 : real(dp) :: f(4,nx,ny)
468 :
469 : include 'formats'
470 :
471 20 : ierr = 0
472 20 : hx=x1-x0
473 20 : hxi=1.0d0/hx
474 20 : hx2=hx*hx
475 :
476 20 : xp=(xget-x0)*hxi
477 20 : xpi=1.0d0-xp
478 20 : xp2=xp*xp
479 20 : xpi2=xpi*xpi
480 :
481 20 : ax=xp2*(3.0d0-2.0d0*xp)
482 20 : axbar=1.0d0-ax
483 :
484 20 : bx=-xp2*xpi
485 20 : bxbar=xpi2*xp
486 :
487 20 : cx=xp*(xp2-1.0d0)
488 20 : cxi=xpi*(xpi2-1.0d0)
489 20 : cxd=3.0d0*xp2-1.0d0
490 20 : cxdi=-3.0d0*xpi2+1.0d0
491 :
492 20 : hy=y1-y0
493 20 : hyi=1.0d0/hy
494 20 : hy2=hy*hy
495 :
496 20 : yp=(yget-y0)*hyi
497 20 : ypi=1.0d0-yp
498 20 : yp2=yp*yp
499 20 : ypi2=ypi*ypi
500 :
501 20 : ay=yp2*(3.0d0-2.0d0*yp)
502 20 : aybar=1.0d0-ay
503 :
504 20 : by=-yp2*ypi
505 20 : bybar=ypi2*yp
506 :
507 20 : cy=yp*(yp2-1.0d0)
508 20 : cyi=ypi*(ypi2-1.0d0)
509 20 : cyd=3.0d0*yp2-1.0d0
510 20 : cydi=-3.0d0*ypi2+1.0d0
511 :
512 20 : sixth_hx2 = one_sixth*hx2
513 20 : sixth_hy2 = one_sixth*hy2
514 20 : z36th_hx2_hy2 = z36th*hx2*hy2
515 :
516 20 : sixth_hx = one_sixth*hx
517 20 : sixth_hxi_hy2 = one_sixth*hxi*hy2
518 20 : z36th_hx_hy2 = z36th*hx*hy2
519 :
520 20 : sixth_hx2_hyi = one_sixth*hx2*hyi
521 20 : sixth_hy = one_sixth*hy
522 20 : z36th_hx2_hy = z36th*hx2*hy
523 :
524 20 : ip1 = i+1
525 20 : jp1 = j+1
526 :
527 580 : do k=1,num_ion_vals
528 : ! bicubic spline interpolation
529 : res(k) = dble(&
530 : xpi*(&
531 : ypi*fin(1,k,i,j) +yp*fin(1,k,i,jp1))&
532 : +xp*(ypi*fin(1,k,ip1,j)+yp*fin(1,k,ip1,jp1))&
533 : +sixth_hx2*(&
534 : cxi*(ypi*fin(2,k,i,j) +yp*fin(2,k,i,jp1))+&
535 : cx*(ypi*fin(2,k,ip1,j)+yp*fin(2,k,ip1,jp1)))&
536 : +sixth_hy2*(&
537 : xpi*(cyi*fin(3,k,i,j) +cy*fin(3,k,i,jp1))+&
538 : xp*(cyi*fin(3,k,ip1,j)+cy*fin(3,k,ip1,jp1)))&
539 : +z36th_hx2_hy2*(&
540 : cxi*(cyi*fin(4,k,i,j) +cy*fin(4,k,i,jp1))+&
541 580 : cx*(cyi*fin(4,k,ip1,j)+cy*fin(4,k,ip1,jp1))))
542 :
543 : ! derivatives of bicubic splines
544 : ! df_dx(k) =
545 : ! > hxi*(
546 : ! > -(ypi*fin(1,k,i,j) +yp*fin(1,k,i,jp1))
547 : ! > +(ypi*fin(1,k,ip1,j)+yp*fin(1,k,ip1,jp1)))
548 : ! > +sixth_hx*(
549 : ! > cxdi*(ypi*fin(2,k,i,j) +yp*fin(2,k,i,jp1))+
550 : ! > cxd*(ypi*fin(2,k,ip1,j)+yp*fin(2,k,ip1,jp1)))
551 : ! > +sixth_hxi_hy2*(
552 : ! > -(cyi*fin(3,k,i,j) +cy*fin(3,k,i,jp1))
553 : ! > +(cyi*fin(3,k,ip1,j)+cy*fin(3,k,ip1,jp1)))
554 : ! > +z36th_hx_hy2*(
555 : ! > cxdi*(cyi*fin(4,k,i,j) +cy*fin(4,k,i,jp1))+
556 : ! > cxd*(cyi*fin(4,k,ip1,j)+cy*fin(4,k,ip1,jp1)))
557 : !
558 : ! df_dy(k) =
559 : ! > hyi*(
560 : ! > xpi*(-fin(1,k,i,j) +fin(1,k,i,jp1))+
561 : ! > xp*(-fin(1,k,ip1,j)+fin(1,k,ip1,jp1)))
562 : ! > +sixth_hx2_hyi*(
563 : ! > cxi*(-fin(2,k,i,j) +fin(2,k,i,jp1))+
564 : ! > cx*(-fin(2,k,ip1,j)+fin(2,k,ip1,jp1)))
565 : ! > +sixth_hy*(
566 : ! > xpi*(cydi*fin(3,k,i,j) +cyd*fin(3,k,i,jp1))+
567 : ! > xp*(cydi*fin(3,k,ip1,j)+cyd*fin(3,k,ip1,jp1)))
568 : ! > +z36th_hx2_hy*(
569 : ! > cxi*(cydi*fin(4,k,i,j) +cyd*fin(4,k,i,jp1))+
570 : ! > cx*(cydi*fin(4,k,ip1,j)+cyd*fin(4,k,ip1,jp1)))
571 :
572 : end do
573 :
574 20 : end subroutine Do_ion_Interpolations
575 :
576 : end module ion_tables_eval
|
