1 subroutine r8akherm2(x,nx,y,ny,fherm,nf2,ilinx,iliny,ier)
11 INTEGER,
PARAMETER :: R8=selected_real_kind(12,100)
12 INTEGER ilinx,iliny,ier,nf2
17 real*8 fherm(0:3,nf2,ny)
35 call r8akherm2p(x,nx,y,ny,fherm,nf2,ilinx,iliny,0,0,ier)
40 subroutine r8akherm2p(x,nx,y,ny,fherm,nf2,ilinx,iliny,ipx,ipy,ier)
54 INTEGER,
PARAMETER :: R8=selected_real_kind(12,100)
55 INTEGER ilinx,iliny,ier,nf2,ierx,iery,ierbc,iy,ix,icx,icy
56 INTEGER incx,incy,ix1,iy1,iym2,iym1,iymm2,iymm1,iyp2,iyp1
57 INTEGER iypp2,iypp1,ixm2,ixm1,ixmm2,ixmm1,iflagx,ixp2,ixp1
58 INTEGER ixpp2,ixpp1,iflagy
61 real*8 cxp,cxm,cxpp,cxmm,cxtrap0,cxtrapn,cyp,cym,cypp,cymm
62 real*8 cytrap0,cytrapn,zansr,cxm2,cxp2
67 real*8 fherm(0:3,nf2,ny)
86 real*8 wx(2),wy(2),e(2,2)
88 real*8,
dimension(:,:),
allocatable :: ftmp
99 call r8splinck(x,nx,ilinx,1.0e-3_r8,ierx)
100 if(ierx.ne.0) ier=ier+1
103 write(6,
'('' ?akherm2: x axis not strict ascending'')')
106 call r8splinck(y,ny,iliny,1.0e-3_r8,iery)
107 if(iery.ne.0) ier=ier+1
110 write(6,
'('' ?akherm2: y axis not strict ascending'')')
115 write(6,*)
'?akherm2: fherm array dimension too small.'
119 call ibc_ck(ipx,
'akherm2',
'X Bdy Cond',0,2,ierbc)
123 call ibc_ck(ipy,
'akherm2',
'Y Bdy Cond',0,2,ierbc)
134 allocate(ftmp(0:nx+1,0:ny+1))
138 ftmp(ix,iy)=fherm(0,ix,iy)
147 xx(nx+1)=2*x(nx)-x(nx-1)
149 yy(ny+1)=2*y(ny)-y(ny-1)
165 cxp=(ftmp(2,iy)-ftmp(1,iy))/(xx(2)-xx(1))
166 cxm=(ftmp(nx,iy)-ftmp(nx-1,iy))/(xx(nx)-xx(nx-1))
173 cxpp=(ftmp(3,iy)-ftmp(2,iy))/(xx(3)-xx(2))
174 cxmm=(ftmp(nx-1,iy)-ftmp(nx-2,iy))/(xx(nx-1)-xx(nx-2))
176 call r8akherm0(cxmm,cxm,cxp,cxpp,wx,fherm(1,1,iy))
177 fherm(1,nx,iy)=fherm(1,1,iy)
180 fherm(1,nx,iy) = fherm(1,1,iy)
186 else if(ipx.eq.0)
then
191 cxpp=(ftmp(3,iy)-ftmp(2,iy))/(xx(3)-xx(2))
192 fherm(1,1,iy)=1.5_r8*cxp-0.5_r8*cxpp
194 cxmm=(ftmp(nx-1,iy)-ftmp(nx-2,iy))/(xx(nx-1)-xx(nx-2))
195 fherm(1,nx,iy)=1.5_r8*cxm-0.5_r8*cxmm
199 fherm(1,nx,iy) = fherm(1,1,iy)
204 cxtrap0=2.0_r8*fherm(1,1,iy)-cxp
205 cxtrapn=2.0_r8*fherm(1,nx,iy)-cxm
211 cxtrap0=2.0_r8*fherm(1,1,iy)-cxp
212 cxtrapn=2.0_r8*fherm(1,nx,iy)-cxm
216 ftmp(0,iy)=ftmp(1,iy)-cxtrap0*(xx(1)-xx(0))
217 ftmp(nx+1,iy)=ftmp(nx,iy)+cxtrapn*(xx(nx+1)-xx(nx))
225 cyp=(ftmp(ix,2)-ftmp(ix,1))/(yy(2)-yy(1))
226 cym=(ftmp(ix,ny)-ftmp(ix,ny-1))/(yy(ny)-yy(ny-1))
233 cypp=(ftmp(ix,3)-ftmp(ix,2))/(yy(3)-yy(2))
234 cymm=(ftmp(ix,ny-1)-ftmp(ix,ny-2))/(yy(ny-1)-yy(ny-2))
236 call r8akherm0(cymm,cym,cyp,cypp,wy,fherm(2,ix,1))
237 fherm(2,ix,ny)=fherm(2,ix,1)
241 fherm(2,ix,ny)=fherm(2,ix,1)
247 else if(ipy.eq.0)
then
252 cypp=(ftmp(ix,3)-ftmp(ix,2))/(yy(3)-yy(2))
253 fherm(2,ix,1)=1.5_r8*cyp-0.5_r8*cypp
255 cymm=(ftmp(ix,ny-1)-ftmp(ix,ny-2))/(yy(ny-1)-yy(ny-2))
256 fherm(2,ix,ny)=1.5_r8*cym-0.5_r8*cymm
260 fherm(2,ix,ny)=fherm(2,ix,1)
265 cytrap0=2.0_r8*fherm(2,ix,1)-cyp
266 cytrapn=2.0_r8*fherm(2,ix,ny)-cym
272 cytrap0=2.0_r8*fherm(2,ix,1)-cyp
273 cytrapn=2.0_r8*fherm(2,ix,ny)-cym
277 ftmp(ix,0)=ftmp(ix,1)-cytrap0*(yy(1)-yy(0))
278 ftmp(ix,ny+1)=ftmp(ix,ny)+cytrapn*(yy(ny+1)-yy(ny))
294 ftmp(icx,icy)=ftmp(icx,iy1)+ftmp(ix1,icy)-ftmp(ix1,iy1)
327 cxm=(ftmp(ixm2,iy)-ftmp(ixm1,iy))/(xx(ixm2)-xx(ixm1))
329 cxmm=(ftmp(ixmm2,iy)-ftmp(ixmm1,iy))/
330 > (xx(ixmm2)-xx(ixmm1))
333 cxmm=(ftmp(nx-1,iy)-ftmp(nx-2,iy))/
334 > (xx(nx-1)-xx(nx-2))
346 cxp=(ftmp(ixp2,iy)-ftmp(ixp1,iy))/(xx(ixp2)-xx(ixp1))
348 cxpp=(ftmp(ixpp2,iy)-ftmp(ixpp1,iy))/
349 > (xx(ixpp2)-xx(ixpp1))
352 cxpp=(ftmp(3,iy)-ftmp(2,iy))/(xx(3)-xx(2))
364 call r8akherm0(cxmm,cxm,cxp,cxpp,wx,zansr)
365 if((ix.gt.1).and.(ix.lt.nx)) fherm(1,ix,iy)=zansr
370 cym=(ftmp(ix,iym2)-ftmp(ix,iym1))/(yy(iym2)-yy(iym1))
372 cymm=(ftmp(ix,iymm2)-ftmp(ix,iymm1))/
373 > (yy(iymm2)-yy(iymm1))
376 cymm=(ftmp(ix,ny-1)-ftmp(ix,ny-2))/
377 > (yy(ny-1)-yy(ny-2))
383 cyp=(ftmp(ix,iyp2)-ftmp(ix,iyp1))/(yy(iyp2)-yy(iyp1))
385 cypp=(ftmp(ix,iypp2)-ftmp(ix,iypp1))/
386 > (yy(iypp2)-yy(iypp1))
389 cypp=(ftmp(ix,3)-ftmp(ix,2))/(yy(3)-yy(2))
401 call r8akherm0(cymm,cym,cyp,cypp,wy,zansr)
402 if((iy.gt.1).and.(iy.lt.ny)) fherm(2,ix,iy)=zansr
406 cxm2=(ftmp(ixm2,iym1)-ftmp(ixm1,iym1))/
407 > (xx(ixm2)-xx(ixm1))
408 e(1,1)=(cxm-cxm2)/(yy(iym2)-yy(iym1))
410 cxm2=(ftmp(ixm2,iyp2)-ftmp(ixm1,iyp2))/
411 > (xx(ixm2)-xx(ixm1))
412 e(1,2)=(cxm2-cxm)/(yy(iyp2)-yy(iyp1))
414 cxp2=(ftmp(ixp2,iym1)-ftmp(ixp1,iym1))/
415 > (xx(ixp2)-xx(ixp1))
416 e(2,1)=(cxp-cxp2)/(yy(iym2)-yy(iym1))
418 cxp2=(ftmp(ixp2,iyp2)-ftmp(ixp1,iyp2))/
419 > (xx(ixp2)-xx(ixp1))
420 e(2,2)=(cxp2-cxp)/(yy(iyp2)-yy(iyp1))
424 fherm(3,ix,iy)=(wx(1)*(wy(1)*e(1,1)+wy(2)*e(1,2))+
425 > wx(2)*(wy(1)*e(2,1)+wy(2)*e(2,2)))/
426 > ((wx(1)+wx(2))*(wy(1)+wy(2)))