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