1 subroutine akherm2(x,nx,y,ny,fherm,nf2,ilinx,iliny,ier)
11 real fherm(0:3,nf2,ny)
29 call akherm2p(x,nx,y,ny,fherm,nf2,ilinx,iliny,0,0,ier)
34 subroutine akherm2p(x,nx,y,ny,fherm,nf2,ilinx,iliny,ipx,ipy,ier)
48 real fherm(0:3,nf2,ny)
67 real wx(2),wy(2),e(2,2)
69 real,
dimension(:,:),
allocatable :: ftmp
80 call splinck(x,nx,ilinx,1.0e-3,ierx)
81 if(ierx.ne.0) ier=ier+1
84 write(6,
'('' ?akherm2: x axis not strict ascending'')')
87 call splinck(y,ny,iliny,1.0e-3,iery)
88 if(iery.ne.0) ier=ier+1
91 write(6,
'('' ?akherm2: y axis not strict ascending'')')
96 write(6,*)
'?akherm2: fherm array dimension too small.'
100 call ibc_ck(ipx,
'akherm2',
'X Bdy Cond',0,2,ierbc)
104 call ibc_ck(ipy,
'akherm2',
'Y Bdy Cond',0,2,ierbc)
115 allocate(ftmp(0:nx+1,0:ny+1))
119 ftmp(ix,iy)=fherm(0,ix,iy)
128 xx(nx+1)=2*x(nx)-x(nx-1)
130 yy(ny+1)=2*y(ny)-y(ny-1)
146 cxp=(ftmp(2,iy)-ftmp(1,iy))/(xx(2)-xx(1))
147 cxm=(ftmp(nx,iy)-ftmp(nx-1,iy))/(xx(nx)-xx(nx-1))
154 cxpp=(ftmp(3,iy)-ftmp(2,iy))/(xx(3)-xx(2))
155 cxmm=(ftmp(nx-1,iy)-ftmp(nx-2,iy))/(xx(nx-1)-xx(nx-2))
157 call akherm0(cxmm,cxm,cxp,cxpp,wx,fherm(1,1,iy))
158 fherm(1,nx,iy)=fherm(1,1,iy)
161 fherm(1,nx,iy) = fherm(1,1,iy)
167 else if(ipx.eq.0)
then
172 cxpp=(ftmp(3,iy)-ftmp(2,iy))/(xx(3)-xx(2))
173 fherm(1,1,iy)=1.5*cxp-0.5*cxpp
175 cxmm=(ftmp(nx-1,iy)-ftmp(nx-2,iy))/(xx(nx-1)-xx(nx-2))
176 fherm(1,nx,iy)=1.5*cxm-0.5*cxmm
180 fherm(1,nx,iy) = fherm(1,1,iy)
185 cxtrap0=2.0*fherm(1,1,iy)-cxp
186 cxtrapn=2.0*fherm(1,nx,iy)-cxm
192 cxtrap0=2.0*fherm(1,1,iy)-cxp
193 cxtrapn=2.0*fherm(1,nx,iy)-cxm
197 ftmp(0,iy)=ftmp(1,iy)-cxtrap0*(xx(1)-xx(0))
198 ftmp(nx+1,iy)=ftmp(nx,iy)+cxtrapn*(xx(nx+1)-xx(nx))
206 cyp=(ftmp(ix,2)-ftmp(ix,1))/(yy(2)-yy(1))
207 cym=(ftmp(ix,ny)-ftmp(ix,ny-1))/(yy(ny)-yy(ny-1))
214 cypp=(ftmp(ix,3)-ftmp(ix,2))/(yy(3)-yy(2))
215 cymm=(ftmp(ix,ny-1)-ftmp(ix,ny-2))/(yy(ny-1)-yy(ny-2))
217 call akherm0(cymm,cym,cyp,cypp,wy,fherm(2,ix,1))
218 fherm(2,ix,ny)=fherm(2,ix,1)
222 fherm(2,ix,ny)=fherm(2,ix,1)
228 else if(ipy.eq.0)
then
233 cypp=(ftmp(ix,3)-ftmp(ix,2))/(yy(3)-yy(2))
234 fherm(2,ix,1)=1.5*cyp-0.5*cypp
236 cymm=(ftmp(ix,ny-1)-ftmp(ix,ny-2))/(yy(ny-1)-yy(ny-2))
237 fherm(2,ix,ny)=1.5*cym-0.5*cymm
241 fherm(2,ix,ny)=fherm(2,ix,1)
246 cytrap0=2.0*fherm(2,ix,1)-cyp
247 cytrapn=2.0*fherm(2,ix,ny)-cym
253 cytrap0=2.0*fherm(2,ix,1)-cyp
254 cytrapn=2.0*fherm(2,ix,ny)-cym
258 ftmp(ix,0)=ftmp(ix,1)-cytrap0*(yy(1)-yy(0))
259 ftmp(ix,ny+1)=ftmp(ix,ny)+cytrapn*(yy(ny+1)-yy(ny))
275 ftmp(icx,icy)=ftmp(icx,iy1)+ftmp(ix1,icy)-ftmp(ix1,iy1)
308 cxm=(ftmp(ixm2,iy)-ftmp(ixm1,iy))/(xx(ixm2)-xx(ixm1))
310 cxmm=(ftmp(ixmm2,iy)-ftmp(ixmm1,iy))/
311 > (xx(ixmm2)-xx(ixmm1))
314 cxmm=(ftmp(nx-1,iy)-ftmp(nx-2,iy))/
315 > (xx(nx-1)-xx(nx-2))
327 cxp=(ftmp(ixp2,iy)-ftmp(ixp1,iy))/(xx(ixp2)-xx(ixp1))
329 cxpp=(ftmp(ixpp2,iy)-ftmp(ixpp1,iy))/
330 > (xx(ixpp2)-xx(ixpp1))
333 cxpp=(ftmp(3,iy)-ftmp(2,iy))/(xx(3)-xx(2))
345 call akherm0(cxmm,cxm,cxp,cxpp,wx,zansr)
346 if((ix.gt.1).and.(ix.lt.nx)) fherm(1,ix,iy)=zansr
351 cym=(ftmp(ix,iym2)-ftmp(ix,iym1))/(yy(iym2)-yy(iym1))
353 cymm=(ftmp(ix,iymm2)-ftmp(ix,iymm1))/
354 > (yy(iymm2)-yy(iymm1))
357 cymm=(ftmp(ix,ny-1)-ftmp(ix,ny-2))/
358 > (yy(ny-1)-yy(ny-2))
364 cyp=(ftmp(ix,iyp2)-ftmp(ix,iyp1))/(yy(iyp2)-yy(iyp1))
366 cypp=(ftmp(ix,iypp2)-ftmp(ix,iypp1))/
367 > (yy(iypp2)-yy(iypp1))
370 cypp=(ftmp(ix,3)-ftmp(ix,2))/(yy(3)-yy(2))
382 call akherm0(cymm,cym,cyp,cypp,wy,zansr)
383 if((iy.gt.1).and.(iy.lt.ny)) fherm(2,ix,iy)=zansr
387 cxm2=(ftmp(ixm2,iym1)-ftmp(ixm1,iym1))/
388 > (xx(ixm2)-xx(ixm1))
389 e(1,1)=(cxm-cxm2)/(yy(iym2)-yy(iym1))
391 cxm2=(ftmp(ixm2,iyp2)-ftmp(ixm1,iyp2))/
392 > (xx(ixm2)-xx(ixm1))
393 e(1,2)=(cxm2-cxm)/(yy(iyp2)-yy(iyp1))
395 cxp2=(ftmp(ixp2,iym1)-ftmp(ixp1,iym1))/
396 > (xx(ixp2)-xx(ixp1))
397 e(2,1)=(cxp-cxp2)/(yy(iym2)-yy(iym1))
399 cxp2=(ftmp(ixp2,iyp2)-ftmp(ixp1,iyp2))/
400 > (xx(ixp2)-xx(ixp1))
401 e(2,2)=(cxp2-cxp)/(yy(iyp2)-yy(iyp1))
405 fherm(3,ix,iy)=(wx(1)*(wy(1)*e(1,1)+wy(2)*e(1,2))+
406 > wx(2)*(wy(1)*e(2,1)+wy(2)*e(2,2)))/
407 > ((wx(1)+wx(2))*(wy(1)+wy(2)))