V3FIT
r8akherm3.f
1  subroutine r8akherm3(x,nx,y,ny,z,nz,fherm,nf2,nf3,
2  > ilinx,iliny,ilinz,ier)
3 C
4 C create a data set for Hermite interpolation, based on Akima's method
5 C [Hiroshi Akima, Communications of the ACM, Jan 1974, Vol. 17 No. 1]
6 C
7 C 3d routine
8 C
9 C input:
10 C
11 C nf2.get.nx and nf3.ge.ny expected!
12 C
13 !============
14 ! idecl: explicitize implicit INTEGER declarations:
15  IMPLICIT NONE
16  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
17  INTEGER iliny,ilinz,ier,ilinx
18 !============
19  integer nx,ny,nz,nf2,nf3 ! array dimensions
20  real*8 x(nx) ! x coordinate array
21  real*8 y(ny) ! y coordinate array
22  real*8 z(nz) ! z coordinate array
23  real*8 fherm(0:7,nf2,nf3,nz) ! data/Hermite array
24 C
25 C fherm(0,i,j,k) = function value f at x(i),y(j),z(k) **on input**
26 C
27 C fherm(1,i,j,k) = derivative df/dx at x(i),y(j),z(k) **on output**
28 C fherm(2,i,j,k) = derivative df/dy at x(i),y(j),z(k) **on output**
29 C fherm(3,i,j,k) = derivative df/dz at x(i),y(j),z(k) **on output**
30 C fherm(4,i,j,k) = derivative d2f/dxdy at x(i),y(j),z(k) **on output**
31 C fherm(5,i,j,k) = derivative d2f/dxdz at x(i),y(j),z(k) **on output**
32 C fherm(6,i,j,k) = derivative d2f/dydz at x(i),y(j),z(k) **on output**
33 C fherm(7,i,j,k) = derivative d3f/dxdydz at x(i),y(j),z(k) **on output**
34 C
35 C additional outputs:
36 C
37 C ilinx = 1 if x is evenly spaced (approximately)
38 C iliny = 1 if y is evenly spaced (approximately)
39 C ilinz = 1 if z is evenly spaced (approximately)
40 C
41 C ier = 0 if there are no errors
42 C
43 C note possible errors: x y or z NOT strict ascending
44 C nf2.lt.nx .or. nf3.lt.ny
45 C
46 C a default boundary condition is used, based on divided differences
47 C in the edge region. For more control of BC, use akherm3p...
48 C
49  call r8akherm3p(x,nx,y,ny,z,nz,fherm,nf2,nf3,
50  > ilinx,iliny,ilinz,0,0,0,ier)
51 C
52  return
53  end
54 C----------------------------
55  subroutine r8akherm3p(x,nx,y,ny,z,nz,fherm,nf2,nf3,
56  > ilinx,iliny,ilinz,ipx,ipy,ipz,ier)
57 C
58 C create a data set for Hermite interpolation, based on Akima's method
59 C [Hiroshi Akima, Communications of the ACM, Jan 1974, Vol. 17 No. 1]
60 C
61 C 3d routine --
62 C with independently settable boundary condition options:
63 C ipx|ipy|ipz =0 -- default, boundary conditions from divided diffs
64 C ipx|ipy|ipz =1 -- periodic boundary condition
65 C ipx|ipy|ipz =2 -- user supplied df/dx|df/dy|df/dz
66 C
67 C input:
68 C
69 C nf2.get.nx and nf3.ge.ny expected!
70 C
71 !============
72 ! idecl: explicitize implicit INTEGER declarations:
73  IMPLICIT NONE
74  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
75  INTEGER iliny,ilinz,ier,ilinx,ierx,iery,ierz,ierbc,iz,iy,ix
76  INTEGER icx,incx,ix1,icy,incy,iy1,icz,incz,iz1,izm2,izm1
77  INTEGER izmm2,izmm1,izp2,izp1,izpp2,izpp1,iym2,iym1,iymm2
78  INTEGER iymm1,iyp2,iyp1,iypp2,iypp1,ixm2,ixm1,ixmm2,ixmm1
79  INTEGER iflagx,ixp2,ixp1,ixpp2,ixpp1,iflagy,iflagz,k
80 !============
81 ! idecl: explicitize implicit REAL declarations:
82  real*8 cxp,cxm,cxpp,cxmm,cxtrap0,cxtrapn,cyp,cym,cypp,cymm
83  real*8 cytrap0,cytrapn,czp,czm,czpp,czmm,cztrap0,cztrapn
84  real*8 zansr,cxm2,cxp2,cym2,cyp2,fzp,fzm
85 !============
86  integer nx,ny,nz,nf2,nf3 ! array dimensions
87  REAL*8 x(nx) ! x coordinate array
88  real*8 y(ny) ! y coordinate array
89  real*8 z(nz) ! z coordinate array
90  real*8 fherm(0:7,nf2,nf3,nz) ! data/Hermite array
91 C
92  integer ipx ! =0: dflt; =1: periodic df/dx
93  integer ipy ! =0: dflt; =1: periodic df/dy
94  integer ipz ! =0: dflt; =1: periodic df/dz
95 C
96 C if ipx=2: fherm(1,1,1:ny,1:nz) & fherm(1,nx,1:ny,1:nz) = INPUT df/dx BCs
97 C if ipy=2: fherm(2,1:nx,1,1:nz) & fherm(2,1:nx,ny,1:nz) = INPUT df/dy BCs
98 C if ipz=2: fherm(3,1:nx,1:ny,1) & fherm(3,1:nx,1:ny,nz) = INPUT df/dz BCs
99 C
100 C for the rest of the grid points...
101 C
102 C fherm(0,i,j,k) = function value f at x(i),y(j),z(k) **on input**
103 C
104 C fherm(1,i,j,k) = derivative df/dx at x(i),y(j),z(k) **on output**
105 C fherm(2,i,j,k) = derivative df/dy at x(i),y(j),z(k) **on output**
106 C fherm(3,i,j,k) = derivative df/dz at x(i),y(j),z(k) **on output**
107 C fherm(4,i,j,k) = derivative d2f/dxdy at x(i),y(j),z(k) **on output**
108 C fherm(5,i,j,k) = derivative d2f/dxdz at x(i),y(j),z(k) **on output**
109 C fherm(6,i,j,k) = derivative d2f/dydz at x(i),y(j),z(k) **on output**
110 C fherm(7,i,j,k) = derivative d3f/dxdydz at x(i),y(j),z(k) **on output**
111 C
112 C additional outputs:
113 C
114 C ilinx = 1 if x is evenly spaced (approximately)
115 C iliny = 1 if y is evenly spaced (approximately)
116 C ilinz = 1 if z is evenly spaced (approximately)
117 C
118 C---------------------
119 C local...
120 C
121  real*8 wx(2),wy(2),wz(2),exy(2,2),exz(2,2),eyz(2,2),fxyz(2,2,2)
122 C
123  real*8, dimension(:,:,:), allocatable :: ftmp
124 c
125  real*8 xx(0:nx+1)
126  real*8 yy(0:ny+1)
127  real*8 zz(0:nz+1)
128 C
129 C----------------------------
130 C
131 C error check
132 C
133  ier=0
134 c
135  call r8splinck(x,nx,ilinx,1.0e-3_r8,ierx)
136  if(ierx.ne.0) ier=ier+1
137 c
138  if(ierx.ne.0) then
139  write(6,'('' ?akherm3: x axis not strict ascending'')')
140  endif
141 c
142  call r8splinck(y,ny,iliny,1.0e-3_r8,iery)
143  if(iery.ne.0) ier=ier+1
144 c
145  if(iery.ne.0) then
146  write(6,'('' ?akherm3: y axis not strict ascending'')')
147  endif
148 c
149  call r8splinck(z,nz,ilinz,1.0e-3_r8,ierz)
150  if(ierz.ne.0) ier=ier+1
151 c
152  if(ierz.ne.0) then
153  write(6,'('' ?akherm3: z axis not strict ascending'')')
154  endif
155 c
156  if(nf2.lt.nx) then
157  ier=ier+1
158  write(6,*) '?akherm3: fherm (x) array dimension too small.'
159  endif
160 c
161  if(nf3.lt.ny) then
162  ier=ier+1
163  write(6,*) '?akherm3: fherm (y) array dimension too small.'
164  endif
165 C
166  ierbc=0
167  call ibc_ck(ipx,'akherm3','X Bdy Cond',0,2,ierbc)
168  ier=ier+ierbc
169 C
170  ierbc=0
171  call ibc_ck(ipy,'akherm3','Y Bdy Cond',0,2,ierbc)
172  ier=ier+ierbc
173 C
174  ierbc=0
175  call ibc_ck(ipz,'akherm3','Z Bdy Cond',0,2,ierbc)
176  ier=ier+ierbc
177 C
178  if(ier.ne.0) return
179 C
180 C error check OK
181 C
182 C---------------------------------------
183 C
184 C get a temporary array for f -- will extend out by 1 zone in
185 C each direction so that numerical derivative evaluation can be
186 C done without a lot of special case logic near the edges...
187 C
188  allocate(ftmp(0:nx+1,0:ny+1,0:nz+1))
189 C
190  do iz=1,nz
191  do iy=1,ny
192  do ix=1,nx
193  ftmp(ix,iy,iz)=fherm(0,ix,iy,iz)
194  enddo
195  enddo
196  enddo
197 C
198 C also create expanded axes grids...
199 C
200  xx(1:nx)=x
201  yy(1:ny)=y
202  zz(1:nz)=z
203 c
204  xx(0)=2*x(1)-x(2)
205  xx(nx+1)=2*x(nx)-x(nx-1)
206  yy(0)=2*y(1)-y(2)
207  yy(ny+1)=2*y(ny)-y(ny-1)
208  zz(0)=2*z(1)-z(2)
209  zz(nz+1)=2*z(nz)-z(nz-1)
210 C
211 C---------------------------------------
212 C
213 C handle boundary conditions and create rows of extrapolated points
214 C in ftmp. first do ftmp(0,1:ny), ftmp(nx+1,1:ny),
215 C then ftmp(1:nx,0), ftmp(1:nx,ny+1),
216 c then ... fill in the corners
217 c
218 c also, for ipx.le.1 fill in the bdy fherm(1,*) values;
219 c for ipy.le.1 fill in the bdy fherm(2,*) values
220 c
221 c x bc's
222 c
223  do iz=1,nz
224  do iy=1,ny
225 c
226  cxp=(ftmp(2,iy,iz)-ftmp(1,iy,iz))/(xx(2)-xx(1))
227  cxm=(ftmp(nx,iy,iz)-ftmp(nx-1,iy,iz))/(xx(nx)-xx(nx-1))
228 c
229  if(ipx.eq.1) then
230 c
231 c periodic BC
232 c
233  if(nx.gt.2) then
234  cxpp=(ftmp(3,iy,iz)-ftmp(2,iy,iz))/(xx(3)-xx(2))
235  cxmm=(ftmp(nx-1,iy,iz)-ftmp(nx-2,iy,iz))/
236  > (xx(nx-1)-xx(nx-2))
237 c
238  call r8akherm0(cxmm,cxm,cxp,cxpp,wx,fherm(1,1,iy,iz))
239  fherm(1,nx,iy,iz)=fherm(1,1,iy,iz)
240 c
241  else
242  fherm(1,1,iy,iz)=cxp ! =cxm, nx=2
243  fherm(1,nx,iy,iz)=fherm(1,1,iy,iz)
244  endif
245 c
246  cxtrap0=cxm
247  cxtrapn=cxp
248 c
249  else if(ipx.eq.0) then
250 C
251 C default BC -- standard numeric extrapolation
252 C
253  if(nx.gt.2) then
254  cxpp=(ftmp(3,iy,iz)-ftmp(2,iy,iz))/(xx(3)-xx(2))
255  fherm(1,1,iy,iz)=1.5_r8*cxp-0.5_r8*cxpp
256 C
257  cxmm=(ftmp(nx-1,iy,iz)-ftmp(nx-2,iy,iz))/
258  > (xx(nx-1)-xx(nx-2))
259  fherm(1,nx,iy,iz)=1.5_r8*cxm-0.5_r8*cxmm
260 c
261  else
262  fherm(1,1,iy,iz)=cxp ! =cxm, nx=2
263  fherm(1,nx,iy,iz)=fherm(1,1,iy,iz)
264  endif
265 C
266 C extrapolate to slope to ghost points just past bdy...
267 C
268  cxtrap0=2.0_r8*fherm(1,1,iy,iz)-cxp
269  cxtrapn=2.0_r8*fherm(1,nx,iy,iz)-cxm
270 C
271  else
272 C
273 C BC supplied by user. Also use this for extrapolation...
274 C
275  cxtrap0=2.0_r8*fherm(1,1,iy,iz)-cxp
276  cxtrapn=2.0_r8*fherm(1,nx,iy,iz)-cxm
277 C
278  endif
279 C
280  ftmp(0,iy,iz)=ftmp(1,iy,iz)-cxtrap0*(xx(1)-xx(0))
281  ftmp(nx+1,iy,iz)=ftmp(nx,iy,iz)+cxtrapn*(xx(nx+1)-xx(nx))
282 C
283  enddo
284  enddo
285 c
286 c y bc's
287 c
288  do iz=1,nz
289  do ix=1,nx
290 c
291  cyp=(ftmp(ix,2,iz)-ftmp(ix,1,iz))/(yy(2)-yy(1))
292  cym=(ftmp(ix,ny,iz)-ftmp(ix,ny-1,iz))/(yy(ny)-yy(ny-1))
293 c
294  if(ipy.eq.1) then
295 c
296 c periodic BC
297 c
298  if(ny.gt.2) then
299  cypp=(ftmp(ix,3,iz)-ftmp(ix,2,iz))/(yy(3)-yy(2))
300  cymm=(ftmp(ix,ny-1,iz)-ftmp(ix,ny-2,iz))/
301  > (yy(ny-1)-yy(ny-2))
302 c
303  call r8akherm0(cymm,cym,cyp,cypp,wy,fherm(2,ix,1,iz))
304  fherm(2,ix,ny,iz)=fherm(2,ix,1,iz)
305 c
306  else
307  fherm(2,ix,1,iz) = cyp ! =cym, ny=2
308  fherm(2,ix,ny,iz)=fherm(2,ix,1,iz)
309  endif
310  cytrap0=cym
311  cytrapn=cyp
312 c
313  else if(ipy.eq.0) then
314 C
315 C default BC -- standard numeric extrapolation
316 C
317  if(ny.gt.2) then
318  cypp=(ftmp(ix,3,iz)-ftmp(ix,2,iz))/(yy(3)-yy(2))
319  fherm(2,ix,1,iz)=1.5_r8*cyp-0.5_r8*cypp
320 C
321  cymm=(ftmp(ix,ny-1,iz)-ftmp(ix,ny-2,iz))/
322  > (yy(ny-1)-yy(ny-2))
323  fherm(2,ix,ny,iz)=1.5_r8*cym-0.5_r8*cymm
324 c
325  else
326  fherm(2,ix,1,iz) = cyp ! =cym, ny=2
327  fherm(2,ix,ny,iz)=fherm(2,ix,1,iz)
328  endif
329 C
330 C extrapolate to slope to ghost points just past bdy...
331 C
332  cytrap0=2.0_r8*fherm(2,ix,1,iz)-cyp
333  cytrapn=2.0_r8*fherm(2,ix,ny,iz)-cym
334 C
335  else
336 C
337 C BC supplied by user. Also use this for extrapolation...
338 C
339  cytrap0=2.0_r8*fherm(2,ix,1,iz)-cyp
340  cytrapn=2.0_r8*fherm(2,ix,ny,iz)-cym
341 C
342  endif
343 C
344  ftmp(ix,0,iz)=ftmp(ix,1,iz)-cytrap0*(yy(1)-yy(0))
345  ftmp(ix,ny+1,iz)=ftmp(ix,ny,iz)+cytrapn*(yy(ny+1)-yy(ny))
346 C
347  enddo
348  enddo
349 c
350 c z bc's
351 c
352  do iy=1,ny
353  do ix=1,nx
354 c
355  czp=(ftmp(ix,iy,2)-ftmp(ix,iy,1))/(zz(2)-zz(1))
356  czm=(ftmp(ix,iy,nz)-ftmp(ix,iy,nz-1))/(zz(nz)-zz(nz-1))
357 c
358  if(ipy.eq.1) then
359 c
360 c periodic BC
361 c
362  if(nz.gt.2) then
363  czpp=(ftmp(ix,iy,3)-ftmp(ix,iy,2))/(zz(3)-zz(2))
364  czmm=(ftmp(ix,iy,nz-1)-ftmp(ix,iy,nz-2))/
365  > (zz(nz-1)-zz(nz-2))
366 c
367  call r8akherm0(czmm,czm,czp,czpp,wz,fherm(3,ix,iy,1))
368  fherm(3,ix,iy,nz)=fherm(3,ix,iy,1)
369 c
370  else
371  fherm(3,ix,iy,1) = czp
372  fherm(3,ix,iy,nz)=fherm(3,ix,iy,1)
373  endif
374 c
375  cztrap0=czm
376  cztrapn=czp
377 c
378  else if(ipy.eq.0) then
379 C
380 C default BC -- standard numeric extrapolation
381 C
382  if(nz.gt.2) then
383  czpp=(ftmp(ix,iy,3)-ftmp(ix,iy,2))/(zz(3)-zz(2))
384  fherm(3,ix,iy,1)=1.5_r8*czp-0.5_r8*czpp
385 C
386  czmm=(ftmp(ix,iy,nz-1)-ftmp(ix,iy,nz-2))/
387  > (zz(nz-1)-zz(nz-2))
388  fherm(3,ix,iy,nz)=1.5_r8*czm-0.5_r8*czmm
389 c
390  else
391  fherm(3,ix,iy,1) = czp
392  fherm(3,ix,iy,nz)=fherm(3,ix,iy,1)
393  endif
394 C
395 C extrapolate to slope to ghost points just past bdy...
396 C
397  cztrap0=2.0_r8*fherm(3,ix,iy,1)-czp
398  cztrapn=2.0_r8*fherm(3,ix,iy,nz)-czm
399 C
400  else
401 C
402 C BC supplied by user. Also use this for extrapolation...
403 C
404  cztrap0=2.0_r8*fherm(3,ix,iy,1)-czp
405  cztrapn=2.0_r8*fherm(3,ix,iy,nz)-czm
406 C
407  endif
408 C
409  ftmp(ix,iy,0)=ftmp(ix,iy,1)-cztrap0*(zz(1)-zz(0))
410  ftmp(ix,iy,nz+1)=ftmp(ix,iy,nz)+cztrapn*(zz(nz+1)-zz(nz))
411 C
412  enddo
413  enddo
414 c
415 c and for the edges...
416 c
417  do iz=1,nz
418  do ix=0,1
419  icx=ix*(nx+1)
420  incx=1-2*ix
421  ix1=icx+incx
422  do iy=0,1
423  icy=iy*(ny+1)
424  incy=1-2*iy
425  iy1=icy+incy
426 C
427  ftmp(icx,icy,iz)=ftmp(icx,iy1,iz)+ftmp(ix1,icy,iz)-
428  > ftmp(ix1,iy1,iz)
429 C
430  enddo
431  enddo
432  enddo
433 c
434  do iy=1,ny
435  do iz=0,1
436  icz=iz*(nz+1)
437  incz=1-2*iz
438  iz1=icz+incz
439  do ix=0,1
440  icx=ix*(nx+1)
441  incx=1-2*ix
442  ix1=icx+incx
443 C
444  ftmp(icx,iy,icz)=ftmp(icx,iy,iz1)+ftmp(ix1,iy,icz)-
445  > ftmp(ix1,iy,iz1)
446 C
447  enddo
448  enddo
449  enddo
450 c
451  do ix=1,nx
452  do iy=0,1
453  icy=iy*(ny+1)
454  incy=1-2*iy
455  iy1=icy+incy
456  do iz=0,1
457  icz=iz*(nz+1)
458  incz=1-2*iz
459  iz1=icz+incz
460 C
461  ftmp(ix,icy,icz)=ftmp(ix,icy,iz1)+ftmp(ix,iy1,icz)-
462  > ftmp(ix,iy1,iz1)
463 C
464  enddo
465  enddo
466  enddo
467 c
468 c and the corners...
469 c
470  do iz=0,1
471  icz=iz*(nz+1)
472  incz=1-2*iz
473  iz1=icz+incz
474  do iy=0,1
475  icy=iy*(ny+1)
476  incy=1-2*iy
477  iy1=icy+incy
478  do ix=0,1
479  icx=ix*(nx+1)
480  incx=1-2*ix
481  ix1=icx+incx
482 C
483  ftmp(icx,icy,icz)=(
484  > 2*(ftmp(icx,icy,iz1)+ftmp(icx,iy1,icz)+
485  > ftmp(ix1,icy,icz)) -
486  > (ftmp(icx,iy1,iz1)+ftmp(ix1,icy,iz1)+
487  > ftmp(ix1,iy1,icz)))/3
488 C
489  enddo
490  enddo
491  enddo
492 c
493 c-----------------------------------------------------
494 c boundary done; now fill in the interior
495 c
496 c
497  do iz=1,nz
498  izm2=iz
499  izm1=izm2-1
500 c
501  izmm2=iz-1
502  izmm1=izmm2-1
503 c
504  izp2=iz+1
505  izp1=izp2-1
506 c
507  izpp2=iz+2
508  izpp1=izpp2-1
509 c
510  do iy=1,ny
511  iym2=iy
512  iym1=iym2-1
513 c
514  iymm2=iy-1
515  iymm1=iymm2-1
516 c
517  iyp2=iy+1
518  iyp1=iyp2-1
519 c
520  iypp2=iy+2
521  iypp1=iypp2-1
522 c
523  do ix=1,nx
524 c
525 c x div. diffs in vicinity
526 c
527  ixm2=ix
528  ixm1=ixm2-1
529 c
530  ixmm2=ix-1
531  ixmm1=ixmm2-1
532 c
533  iflagx=0
534  cxm=(ftmp(ixm2,iy,iz)-ftmp(ixm1,iy,iz))/
535  > (xx(ixm2)-xx(ixm1))
536  if(ix.gt.1) then
537  cxmm=(ftmp(ixmm2,iy,iz)-ftmp(ixmm1,iy,iz))/
538  > (xx(ixmm2)-xx(ixmm1))
539  else
540  if(ipx.eq.1) then
541  cxmm=(ftmp(nx-1,iy,iz)-ftmp(nx-2,iy,iz))/
542  > (xx(nx-1)-xx(nx-2))
543  else
544  iflagx=1
545  endif
546  endif
547 c
548  ixp2=ix+1
549  ixp1=ixp2-1
550 c
551  ixpp2=ix+2
552  ixpp1=ixpp2-1
553 c
554  cxp=(ftmp(ixp2,iy,iz)-ftmp(ixp1,iy,iz))/
555  > (xx(ixp2)-xx(ixp1))
556  if(ix.lt.nx) then
557  cxpp=(ftmp(ixpp2,iy,iz)-ftmp(ixpp1,iy,iz))/
558  > (xx(ixpp2)-xx(ixpp1))
559  else
560  if(ipx.eq.1) then
561  cxpp=(ftmp(3,iy,iz)-ftmp(2,iy,iz))/
562  > (xx(3)-xx(2))
563  else
564  cxpp=cxp+(cxm-cxmm)
565  endif
566  endif
567 c
568  if(iflagx.eq.1) then
569  cxmm=cxm+(cxp-cxpp)
570  endif
571 c
572 c Akima weightings + df/dx for interior pts
573 c
574  call r8akherm0(cxmm,cxm,cxp,cxpp,wx,zansr)
575  if((ix.gt.1).and.(ix.lt.nx)) fherm(1,ix,iy,iz)=zansr
576 c
577 c y div. diffs in vicinity
578 c
579  iflagy=0
580  cym=(ftmp(ix,iym2,iz)-ftmp(ix,iym1,iz))/
581  > (yy(iym2)-yy(iym1))
582  if(iy.gt.1) then
583  cymm=(ftmp(ix,iymm2,iz)-ftmp(ix,iymm1,iz))/
584  > (yy(iymm2)-yy(iymm1))
585  else
586  if(ipy.eq.1) then
587  cymm=(ftmp(ix,ny-1,iz)-ftmp(ix,ny-2,iz))/
588  > (yy(ny-1)-yy(ny-2))
589  else
590  iflagy=1
591  endif
592  endif
593 c
594  cyp=(ftmp(ix,iyp2,iz)-ftmp(ix,iyp1,iz))/
595  > (yy(iyp2)-yy(iyp1))
596  if(iy.lt.ny) then
597  cypp=(ftmp(ix,iypp2,iz)-ftmp(ix,iypp1,iz))/
598  > (yy(iypp2)-yy(iypp1))
599  else
600  if(ipy.eq.1) then
601  cypp=(ftmp(ix,3,iz)-ftmp(ix,2,iz))/(yy(3)-yy(2))
602  else
603  cypp=cyp+(cym-cymm)
604  endif
605  endif
606 c
607  if(iflagy.eq.1) then
608  cymm=cym+(cyp-cypp)
609  endif
610 c
611 c Akima weightings + df/dy for interior pts
612 c
613  call r8akherm0(cymm,cym,cyp,cypp,wy,zansr)
614  if((iy.gt.1).and.(iy.lt.ny)) fherm(2,ix,iy,iz)=zansr
615 c
616 c z div. diffs in vicinity
617 c
618  iflagz=0
619  czm=(ftmp(ix,iy,izm2)-ftmp(ix,iy,izm1))/
620  > (zz(izm2)-zz(izm1))
621  if(iz.gt.1) then
622  czmm=(ftmp(ix,iy,izmm2)-ftmp(ix,iy,izmm1))/
623  > (zz(izmm2)-zz(izmm1))
624  else
625  if(ipz.eq.1) then
626  czmm=(ftmp(ix,iy,nz-1)-ftmp(ix,iy,nz-2))/
627  > (zz(nz-1)-zz(nz-2))
628  else
629  iflagz=1
630  endif
631  endif
632 c
633  czp=(ftmp(ix,iy,izp2)-ftmp(ix,iy,izp1))/
634  > (zz(izp2)-zz(izp1))
635  if(iz.lt.nz) then
636  czpp=(ftmp(ix,iy,izpp2)-ftmp(ix,iy,izpp1))/
637  > (zz(izpp2)-zz(izpp1))
638  else
639  if(ipz.eq.1) then
640  czpp=(ftmp(ix,iy,3)-ftmp(ix,iy,2))/(zz(3)-zz(2))
641  else
642  czpp=czp+(czm-czmm)
643  endif
644  endif
645 c
646  if(iflagz.eq.1) then
647  czmm=czm+(czp-czpp)
648  endif
649 c
650 c Akima weightings + df/dy for interior pts
651 c
652  call r8akherm0(czmm,czm,czp,czpp,wz,zansr)
653  if((iz.gt.1).and.(iz.lt.nz)) fherm(3,ix,iy,iz)=zansr
654 c
655 c cross derivatives (2nd order divided differences)
656 c
657  cxm2=(ftmp(ixm2,iym1,iz)-ftmp(ixm1,iym1,iz))/
658  > (xx(ixm2)-xx(ixm1))
659  exy(1,1)=(cxm-cxm2)/(yy(iym2)-yy(iym1))
660 c
661  cxm2=(ftmp(ixm2,iyp2,iz)-ftmp(ixm1,iyp2,iz))/
662  > (xx(ixm2)-xx(ixm1))
663  exy(1,2)=(cxm2-cxm)/(yy(iyp2)-yy(iyp1))
664 c
665  cxp2=(ftmp(ixp2,iym1,iz)-ftmp(ixp1,iym1,iz))/
666  > (xx(ixp2)-xx(ixp1))
667  exy(2,1)=(cxp-cxp2)/(yy(iym2)-yy(iym1))
668 c
669  cxp2=(ftmp(ixp2,iyp2,iz)-ftmp(ixp1,iyp2,iz))/
670  > (xx(ixp2)-xx(ixp1))
671  exy(2,2)=(cxp2-cxp)/(yy(iyp2)-yy(iyp1))
672 c
673 c---------
674 c
675  cxm2=(ftmp(ixm2,iy,izm1)-ftmp(ixm1,iy,izm1))/
676  > (xx(ixm2)-xx(ixm1))
677  exz(1,1)=(cxm-cxm2)/(zz(izm2)-zz(izm1))
678 c
679  cxm2=(ftmp(ixm2,iy,izp2)-ftmp(ixm1,iy,izp2))/
680  > (xx(ixm2)-xx(ixm1))
681  exz(1,2)=(cxm2-cxm)/(zz(izp2)-zz(izp1))
682 c
683  cxp2=(ftmp(ixp2,iy,izm1)-ftmp(ixp1,iy,izm1))/
684  > (xx(ixp2)-xx(ixp1))
685  exz(2,1)=(cxp-cxp2)/(zz(izm2)-zz(izm1))
686 c
687  cxp2=(ftmp(ixp2,iy,izp2)-ftmp(ixp1,iy,izp2))/
688  > (xx(ixp2)-xx(ixp1))
689  exz(2,2)=(cxp2-cxp)/(zz(izp2)-zz(izp1))
690 c
691 c---------
692 c
693  cym2=(ftmp(ix,iym2,izm1)-ftmp(ix,iym1,izm1))/
694  > (yy(iym2)-yy(iym1))
695  eyz(1,1)=(cym-cym2)/(zz(izm2)-zz(izm1))
696 c
697  cym2=(ftmp(ix,iym2,izp2)-ftmp(ix,iym1,izp2))/
698  > (yy(iym2)-yy(iym1))
699  eyz(1,2)=(cym2-cym)/(zz(izp2)-zz(izp1))
700 c
701  cyp2=(ftmp(ix,iyp2,izm1)-ftmp(ix,iyp1,izm1))/
702  > (yy(iyp2)-yy(iyp1))
703  eyz(2,1)=(cyp-cyp2)/(zz(izm2)-zz(izm1))
704 c
705  cyp2=(ftmp(ix,iyp2,izp2)-ftmp(ix,iyp1,izp2))/
706  > (yy(iyp2)-yy(iyp1))
707  eyz(2,2)=(cyp2-cyp)/(zz(izp2)-zz(izp1))
708 c
709 c dxdydz 3rd order cross deriv.
710 c
711  k=1
712 c
713  fzp=(ftmp(ixm2,iym2,izm2)-ftmp(ixm2,iym1,izm2)-
714  > ftmp(ixm1,iym2,izm2)+ftmp(ixm1,iym1,izm2))/
715  > ((xx(ixm2)-xx(ixm1))*(yy(iym2)-yy(iym1)))
716  fzm=(ftmp(ixm2,iym2,izm1)-ftmp(ixm2,iym1,izm1)-
717  > ftmp(ixm1,iym2,izm1)+ftmp(ixm1,iym1,izm1))/
718  > ((xx(ixm2)-xx(ixm1))*(yy(iym2)-yy(iym1)))
719  fxyz(1,1,k)=(fzp-fzm)/(zz(izm2)-zz(izm1))
720 c
721  fzp=(ftmp(ixm2,iyp2,izm2)-ftmp(ixm2,iyp1,izm2)-
722  > ftmp(ixm1,iyp2,izm2)+ftmp(ixm1,iyp1,izm2))/
723  > ((xx(ixm2)-xx(ixm1))*(yy(iyp2)-yy(iyp1)))
724  fzm=(ftmp(ixm2,iyp2,izm1)-ftmp(ixm2,iyp1,izm1)-
725  > ftmp(ixm1,iyp2,izm1)+ftmp(ixm1,iyp1,izm1))/
726  > ((xx(ixm2)-xx(ixm1))*(yy(iyp2)-yy(iyp1)))
727  fxyz(1,2,k)=(fzp-fzm)/(zz(izm2)-zz(izm1))
728 c
729  fzp=(ftmp(ixp2,iym2,izm2)-ftmp(ixp2,iym1,izm2)-
730  > ftmp(ixp1,iym2,izm2)+ftmp(ixp1,iym1,izm2))/
731  > ((xx(ixp2)-xx(ixp1))*(yy(iym2)-yy(iym1)))
732  fzm=(ftmp(ixp2,iym2,izm1)-ftmp(ixp2,iym1,izm1)-
733  > ftmp(ixp1,iym2,izm1)+ftmp(ixp1,iym1,izm1))/
734  > ((xx(ixp2)-xx(ixp1))*(yy(iym2)-yy(iym1)))
735  fxyz(2,1,k)=(fzp-fzm)/(zz(izm2)-zz(izm1))
736 c
737  fzp=(ftmp(ixp2,iyp2,izm2)-ftmp(ixp2,iyp1,izm2)-
738  > ftmp(ixp1,iyp2,izm2)+ftmp(ixp1,iyp1,izm2))/
739  > ((xx(ixp2)-xx(ixp1))*(yy(iyp2)-yy(iyp1)))
740  fzm=(ftmp(ixp2,iyp2,izm1)-ftmp(ixp2,iyp1,izm1)-
741  > ftmp(ixp1,iyp2,izm1)+ftmp(ixp1,iyp1,izm1))/
742  > ((xx(ixp2)-xx(ixp1))*(yy(iyp2)-yy(iyp1)))
743  fxyz(2,2,k)=(fzp-fzm)/(zz(izm2)-zz(izm1))
744 c
745  k=2
746 c
747  fzp=(ftmp(ixm2,iym2,izp2)-ftmp(ixm2,iym1,izp2)-
748  > ftmp(ixm1,iym2,izp2)+ftmp(ixm1,iym1,izp2))/
749  > ((xx(ixm2)-xx(ixm1))*(yy(iym2)-yy(iym1)))
750  fzm=(ftmp(ixm2,iym2,izp1)-ftmp(ixm2,iym1,izp1)-
751  > ftmp(ixm1,iym2,izp1)+ftmp(ixm1,iym1,izp1))/
752  > ((xx(ixm2)-xx(ixm1))*(yy(iym2)-yy(iym1)))
753  fxyz(1,1,k)=(fzp-fzm)/(zz(izp2)-zz(izp1))
754 c
755  fzp=(ftmp(ixm2,iyp2,izp2)-ftmp(ixm2,iyp1,izp2)-
756  > ftmp(ixm1,iyp2,izp2)+ftmp(ixm1,iyp1,izp2))/
757  > ((xx(ixm2)-xx(ixm1))*(yy(iyp2)-yy(iyp1)))
758  fzm=(ftmp(ixm2,iyp2,izp1)-ftmp(ixm2,iyp1,izp1)-
759  > ftmp(ixm1,iyp2,izp1)+ftmp(ixm1,iyp1,izp1))/
760  > ((xx(ixm2)-xx(ixm1))*(yy(iyp2)-yy(iyp1)))
761  fxyz(1,2,k)=(fzp-fzm)/(zz(izp2)-zz(izp1))
762 c
763  fzp=(ftmp(ixp2,iym2,izp2)-ftmp(ixp2,iym1,izp2)-
764  > ftmp(ixp1,iym2,izp2)+ftmp(ixp1,iym1,izp2))/
765  > ((xx(ixp2)-xx(ixp1))*(yy(iym2)-yy(iym1)))
766  fzm=(ftmp(ixp2,iym2,izp1)-ftmp(ixp2,iym1,izp1)-
767  > ftmp(ixp1,iym2,izp1)+ftmp(ixp1,iym1,izp1))/
768  > ((xx(ixp2)-xx(ixp1))*(yy(iym2)-yy(iym1)))
769  fxyz(2,1,k)=(fzp-fzm)/(zz(izp2)-zz(izp1))
770 c
771  fzp=(ftmp(ixp2,iyp2,izp2)-ftmp(ixp2,iyp1,izp2)-
772  > ftmp(ixp1,iyp2,izp2)+ftmp(ixp1,iyp1,izp2))/
773  > ((xx(ixp2)-xx(ixp1))*(yy(iyp2)-yy(iyp1)))
774  fzm=(ftmp(ixp2,iyp2,izp1)-ftmp(ixp2,iyp1,izp1)-
775  > ftmp(ixp1,iyp2,izp1)+ftmp(ixp1,iyp1,izp1))/
776  > ((xx(ixp2)-xx(ixp1))*(yy(iyp2)-yy(iyp1)))
777  fxyz(2,2,k)=(fzp-fzm)/(zz(izp2)-zz(izp1))
778 c
779 c the values
780 c
781  fherm(4,ix,iy,iz)=
782  > (wx(1)*(wy(1)*exy(1,1)+wy(2)*exy(1,2))+
783  > wx(2)*(wy(1)*exy(2,1)+wy(2)*exy(2,2)))/
784  > ((wx(1)+wx(2))*(wy(1)+wy(2)))
785 c
786  fherm(5,ix,iy,iz)=
787  > (wx(1)*(wz(1)*exz(1,1)+wz(2)*exz(1,2))+
788  > wx(2)*(wz(1)*exz(2,1)+wz(2)*exz(2,2)))/
789  > ((wx(1)+wx(2))*(wz(1)+wz(2)))
790 c
791  fherm(6,ix,iy,iz)=
792  > (wy(1)*(wz(1)*eyz(1,1)+wz(2)*eyz(1,2))+
793  > wy(2)*(wz(1)*eyz(2,1)+wz(2)*eyz(2,2)))/
794  > ((wy(1)+wy(2))*(wz(1)+wz(2)))
795 c
796  fherm(7,ix,iy,iz)=(wz(1)*
797  > (wx(1)*(wy(1)*fxyz(1,1,1)+wy(2)*fxyz(1,2,1))+
798  > wx(2)*(wy(1)*fxyz(2,1,1)+wy(2)*fxyz(2,2,1))) + wz(2)*
799  > (wx(1)*(wy(1)*fxyz(1,1,2)+wy(2)*fxyz(1,2,2))+
800  > wx(2)*(wy(1)*fxyz(2,1,2)+wy(2)*fxyz(2,2,2))) )/
801  > ((wx(1)+wx(2))*(wy(1)+wy(2))*(wz(1)+wz(2)))
802 c
803  enddo
804  enddo
805  enddo
806 C
807  deallocate (ftmp)
808  return
809  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17