subroutine interpolate(gn, q, p, nx, nd, dd, i) integer nd, nx, i,nd1 parameter (nd1=8) double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1) double precision q(-(nd+1):(nd+1)) double precision p(-(nd+1):(nd+1)) double precision dd double precision z integer b1,b2,b3,b4 parameter (z=1.) parameter (b1=-(nd1+1),b2=(nd1+1),b3=-(nd1+1)/2, & b4=(nd1+1)/2) integer D, U, k, l integer jg(-(nd+1):(nd+1),-(nd+1):(nd+1)) integer ig(-(nd+1):(nd+1),-(nd+1):(nd+1)) integer ig2(-(nd+1):(nd+1),-(nd+1):(nd+1)) integer jg2(-(nd+1):(nd+1),-(nd+1):(nd+1)) integer k1(-(nd+1):(nd+1),-(nd+1):(nd+1)) integer k2(-(nd+1):(nd+1),-(nd+1):(nd+1)) integer l1(-(nd+1):(nd+1),-(nd+1):(nd+1)) integer l2(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision qf(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision pf(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision cp(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision cq(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision cp1(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision cp2(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision cq1(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision cq2(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision q1(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision q2(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision p1(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision p2(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision v1(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision v2(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision q1p(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision q2p(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision p1p(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision p2p(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision r1(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision r2(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision r3(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision r4(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision rx(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision qd1(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision qd2(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision pd1(-(nd+1):(nd+1),-(nd+1):(nd+1)) double precision pd2(-(nd+1):(nd+1),-(nd+1):(nd+1)) c double precision :: q_s2n, p_s2n, cubic c double precision :: q_s2n, p_s2n, cubic INTERFACE PURE DOUBLE PRECISION FUNCTION Q_S2N(QS,PS) DOUBLE PRECISION, INTENT(IN) :: QS,PS END FUNCTION Q_S2N PURE DOUBLE PRECISION FUNCTION P_S2N(QS,PS) DOUBLE PRECISION, INTENT(IN) :: QS,PS END FUNCTION P_S2N PURE DOUBLE PRECISION FUNCTION CUBIC(yp,yp1,yp2,yp3,yp4,f1 &,f2,f3,f4) DOUBLE PRECISION, INTENT(IN) ::yp,yp1,yp2,yp3,yp4,f1,f2 & ,f3,f4 END FUNCTION CUBIC END INTERFACE c Fourth order interpolation at the patch boundary c Upper is the grid we are looking at (rectangular) c Lower is the other patch, with circular arc coordinate lines c do both grids do D = 0,1 if (D.eq.0) U = 1 if (D.eq.1) U = 0 c Find the coordinate of the fictitious point in the lower grid forall (l=b1:b2) qf(b1,l)=q_s2n(q(b1),p(l)) forall (l=b1:b2) pf(b1,l)=p_s2n(q(b1),p(l)) forall (l=b1:b2) qf(b2,l)=q_s2n(q(b2),p(l)) forall (l=b1:b2) pf(b2,l)=p_s2n(q(b2),p(l)) forall (k=-nd:nd)qf(k,b1)=q_s2n(q(k),p(b1)) forall (k=-nd:nd)pf(k,b1)=p_s2n(q(k),p(b1)) forall (k=-nd:nd)qf(k,b2)=q_s2n(q(k),p(b2)) forall (k=-nd:nd)pf(k,b2)=p_s2n(q(k),p(b2)) c All the operations for mask(k,l)=11 c Guard point index and outer neighbor index forall(l=b3:b4) ig(b1,l)=int(qf(b1,l)/dd) forall(l=b3:b4) jg(b1,l)=int(pf(b1,l)/dd) & + int((sign(z,pf(b1,l))-1)/2) forall(l=b3:b4) ig(b2,l)=int(qf(b2,l)/dd) forall(l=b3:b4) jg(b2,l)=int(pf(b2,l)/dd) & + int((sign(z,pf(b2,l))-1)/2) forall(l=b3:b4) ig2(b1,l)=ig(b1,l) & - int(sign(z,q(ig(b1,l)))) forall(l=b3:b4) jg2(b1,l)=jg(b1,l) & - int(sign(z,p(jg(b1,l)))) forall(l=b3:b4) ig2(b2,l)=ig(b2,l) & - int(sign(z,q(ig(b2,l)))) forall(l=b3:b4) jg2(b2,l)=jg(b2,l) & - int(sign(z,p(jg(b2,l)))) c find the upper grid points used in the cubic interpolation forall(l=b3:b4) k1(b1,l) = b1 - 1 * int(sign(z,q(b1))) forall(l=b3:b4) k2(b1,l) = b1 - 2 * int(sign(z,q(b1))) forall(l=b3:b4) l1(b1,l) = l - 1 * int(sign(z,p(l))) forall(l=b3:b4) l2(b1,l) = l - 2 * int(sign(z,p(l))) forall(l=b3:b4) k1(b2,l) = b2 - 1 * int(sign(z,q(b2))) forall(l=b3:b4) k2(b2,l) = b2 - 2 * int(sign(z,q(b2))) forall(l=b3:b4) l1(b2,l) = l - 1 * int(sign(z,p(l))) forall(l=b3:b4) l2(b2,l) = l - 2 * int(sign(z,p(l))) c calculate the intersection points depending on which combination of c crossing lines we use forall(l=b3:b4) q1(b1,l) = (1. + sqrt(1. - 4. * p(l)**2 & * q(ig(b1,l))**2))/ & (2.* q(ig(b1,l))) forall(l=b3:b4) p1p(b1,l) = -2.* p(l) * q(ig(b1,l))**2 & /(1.+ sqrt(1.- 4. * p(l)**2 & * q(ig(b1,l))**2)) forall(l=b3:b4) q2(b1,l) = (1. + sqrt(1. - 4. * p(l)**2 & * q(ig2(b1,l))**2))/ & (2.* q(ig2(b1,l))) forall(l=b3:b4) p2p(b1,l) = -2.* p(l) * q(ig2(b1,l))**2 & /(1. + sqrt(1. - 4.* p(l)**2 & * q(ig2(b1,l))**2)) forall(l=b3:b4) q1(b2,l) = (1. + sqrt(1. - 4. * p(l)**2 & * q(ig(b2,l))**2))/ & (2.* q(ig(b2,l))) forall(l=b3:b4) p1p(b2,l) = -2.* p(l) * q(ig(b2,l))**2 & /(1.+ sqrt(1.- 4. * p(l)**2 & * q(ig(b2,l))**2)) forall(l=b3:b4) q2(b2,l) = (1. + sqrt(1. - 4. * p(l)**2 & * q(ig2(b2,l))**2))/ & (2.* q(ig2(b2,l))) forall(l=b3:b4) p2p(b2,l) = -2.* p(l) * q(ig2(b2,l))**2 & /(1. + sqrt(1. - 4.* p(l)**2 & * q(ig2(b2,l))**2)) c Do the lower interpolation depending on which direction is used forall(l=b3:b4) v1(b1,l) = cubic(p1p(b1,l),p(jg(b1,l)-1), & p(jg(b1,l)),p(jg(b1,l)+1), & p(jg(b1,l)+2),gn(ig(b1,l),jg(b1,l)-1,i,U), & gn(ig(b1,l),jg(b1,l),i,U), & gn(ig(b1,l),jg(b1,l)+1,i,U) & ,gn(ig(b1,l),jg(b1,l)+2,i,U)) forall(l=b3:b4) v2(b1,l) = cubic(p2p(b1,l),p(jg(b1,l)-1), & p(jg(b1,l)),p(jg(b1,l)+1) & ,p(jg(b1,l)+2),gn(ig2(b1,l), & jg(b1,l)-1,i,U), & gn(ig2(b1,l),jg(b1,l),i,U), & gn(ig2(b1,l),jg(b1,l)+1,i,U) & ,gn(ig2(b1,l),jg(b1,l)+2,i,U)) forall(l=b3:b4) v1(b2,l) = cubic(p1p(b2,l),p(jg(b2,l)-1) & ,p(jg(b2,l)),p(jg(b2,l)+1), & p(jg(b2,l)+2),gn(ig(b2,l),jg(b2,l)-1,i,U), & gn(ig(b2,l),jg(b2,l),i,U), & gn(ig(b2,l),jg(b2,l)+1,i,U) & ,gn(ig(b2,l),jg(b2,l)+2,i,U)) forall(l=b3:b4) v2(b2,l) = cubic(p2p(b2,l),p(jg(b2,l)-1) & ,p(jg(b2,l)),p(jg(b2,l)+1) & ,p(jg(b2,l)+2),gn(ig2(b2,l),jg(b2,l)-1,i,U), & gn(ig2(b2,l),jg(b2,l),i,U), & gn(ig2(b2,l),jg(b2,l)+1,i,U) & ,gn(ig2(b2,l),jg(b2,l)+2,i,U)) c Do the upper interpolation depending on which direction is used forall(l=b3:b4) gn(b1,l,i,D) = cubic(q(b1),q(k2(b1,l)), & q(k1(b1,l)),q1(b1,l),q2(b1,l) & ,gn(k2(b1,l),l,i,D), & gn(k1(b1,l),l,i,D),v1(b1,l),v2(b1,l)) forall(l=b3:b4) gn(b2,l,i,D) = cubic(q(b2),q(k2(b2,l)), & q(k1(b2,l)), & q1(b2,l),q2(b2,l) & ,gn(k2(b2,l),l,i,D),gn(k1(b2,l),l,i,D), & v1(b2,l),v2(b2,l)) c All the operations for mask(k,l)=12 c The guard point index and the outer neighbor index forall(l=(b4+1):nd) ig(b1,l) = int(qf(b1,l)/dd) & + int((sign(z,qf(b1,l))-1)/2) forall(l=(b4+1):nd) jg(b1,l) = int(pf(b1,l)/dd) forall(l=(b4+1):nd) ig(b2,l) = int(qf(b2,l)/dd) & + int((sign(z,qf(b2,l))-1)/2) forall(l=(b4+1):nd) jg(b2,l) = int(pf(b2,l)/dd) forall(l=-nd:(b3-1)) ig(b1,l) = int(qf(b1,l)/dd) & + int((sign(z,qf(b1,l))-1)/2) forall(l=-nd:(b3-1)) jg(b1,l) = int(pf(b1,l)/dd) forall(l=-nd:(b3-1)) ig(b2,l) = int(qf(b2,l)/dd) & + int((sign(z,qf(b2,l))-1)/2) forall(l=-nd:(b3-1)) jg(b2,l) = int(pf(b2,l)/dd) forall(l=(b4+1):nd) ig2(b1,l) = ig(b1,l) & - int(sign(z,q(ig(b1,l)))) forall(l=(b4+1):nd) jg2(b1,l) = jg(b1,l) & - int(sign(z,p(jg(b1,l)))) forall(l=(b4+1):nd) ig2(b2,l) = ig(b2,l) & - int(sign(z,q(ig(b2,l)))) forall(l=(b4+1):nd) jg2(b2,l) = jg(b2,l) & - int(sign(z,p(jg(b2,l)))) forall(l=-nd:(b3-1)) ig2(b1,l) = ig(b1,l) & - int(sign(z,q(ig(b1,l)))) forall(l=-nd:(b3-1)) jg2(b1,l) = jg(b1,l) & - int(sign(z,p(jg(b1,l)))) forall(l=-nd:(b3-1)) ig2(b2,l) = ig(b2,l) & - int(sign(z,q(ig(b2,l)))) forall(l=-nd:(b3-1)) jg2(b2,l) = jg(b2,l) & - int(sign(z,p(jg(b2,l)))) c Find the upper grid grid points used in the cubic interpolation forall(l=(b4+1):nd) k1(b1,l) = b1 - 1 * int(sign(z,q(b1))) forall(l=(b4+1):nd) k2(b1,l) = b1 - 2 * int(sign(z,q(b1))) forall(l=(b4+1):nd) k1(b2,l) = b2 - 1 * int(sign(z,q(b2))) forall(l=(b4+1):nd) k2(b2,l) = b2 - 2 * int(sign(z,q(b2))) forall(l=-nd:(b3-1)) k1(b1,l) = b1 - 1 * int(sign(z,q(b1))) forall(l=-nd:(b3-1)) k2(b1,l) = b1 - 2 * int(sign(z,q(b1))) forall(l=-nd:(b3-1)) k1(b2,l) = b2 - 1 * int(sign(z,q(b2))) forall(l=-nd:(b3-1)) k2(b2,l) = b2 - 2 * int(sign(z,q(b2))) c Calculate the intersection points depending on which combination of c crossing lines we use forall(l=(b4+1):nd) q1(b1,l) = sign(z,q(b1)) & * sqrt(-p(l)/p(jg(b1,l)) - p(l)**2) forall(l=(b4+1):nd) q1p(b1,l) = - (p(jg(b1,l)) & /p(l)) * q1(b1,l) forall(l=(b4+1):nd) q2(b1,l) = sign(z,q(b1)) & * sqrt(-p(l)/p(jg2(b1,l)) - p(l)**2) forall(l=(b4+1):nd) q2p(b1,l) = - (p(jg2(b1,l)) & /p(l)) * q2(b1,l) forall(l=(b4+1):nd) q1(b2,l) = sign(z,q(b2)) & * sqrt(-p(l)/p(jg(b2,l)) - p(l)**2) forall(l=(b4+1):nd) q1p(b2,l) = - (p(jg(b2,l)) & /p(l)) * q1(b2,l) forall(l=(b4+1):nd) q2(b2,l) = sign(z,q(b2)) & * sqrt(-p(l)/p(jg2(b2,l)) - p(l)**2) forall(l=(b4+1):nd) q2p(b2,l) = - (p(jg2(b2,l)) & /p(l)) * q2(b2,l) forall(l=-nd:(b3-1)) q1(b1,l) = sign(z,q(b1)) & * sqrt(-p(l)/p(jg(b1,l)) - p(l)**2) forall(l=-nd:(b3-1)) q1p(b1,l) = - (p(jg(b1,l)) & /p(l)) * q1(b1,l) forall(l=-nd:(b3-1)) q2(b1,l) = sign(z,q(b1)) & * sqrt(-p(l)/p(jg2(b1,l)) - p(l)**2) forall(l=-nd:(b3-1)) q2p(b1,l) = - (p(jg2(b1,l)) & /p(l)) * q2(b1,l) forall(l=-nd:(b3-1)) q1(b2,l) = sign(z,q(b2)) & * sqrt(-p(l)/p(jg(b2,l)) - p(l)**2) forall(l=-nd:(b3-1)) q1p(b2,l) = - (p(jg(b2,l)) & /p(l)) * q1(b2,l) forall(l=-nd:(b3-1)) q2(b2,l) = sign(z,q(b2)) * sqrt(-p(l) & /p(jg2(b2,l)) - p(l)**2) forall(l=-nd:(b3-1)) q2p(b2,l) = - (p(jg2(b2,l)) & /p(l)) * q2(b2,l) c Do the lower interpolation depending on which direction is used forall(l=(b4+1):nd) v1(b1,l) = cubic (q1p(b1,l),q(ig(b1,l)-1), & q(ig(b1,l)),q(ig(b1,l)+1) & ,q(ig(b1,l)+2), & gn(ig(b1,l)-1,jg(b1,l),i,U) & ,gn(ig(b1,l),jg(b1,l),i,U), & gn(ig(b1,l)+1,jg(b1,l),i,U) & ,gn(ig(b1,l)+2,jg(b1,l),i,U)) forall(l=(b4+1):nd) v2(b1,l) = cubic (q2p(b1,l),q(ig(b1,l)-1), & q(ig(b1,l)),q(ig(b1,l)+1) & ,q(ig(b1,l)+2),gn(ig(b1,l)-1,jg2(b1,l),i,U) & ,gn(ig(b1,l),jg2(b1,l),i,U), & gn(ig(b1,l)+1,jg2(b1,l),i,U), & gn(ig(b1,l)+2,jg2(b1,l),i,U)) forall(l=-nd:(b3-1)) v1(b1,l) = cubic (q1p(b1,l),q(ig(b1,l)-1), & q(ig(b1,l)),q(ig(b1,l)+1) & ,q(ig(b1,l)+2),gn(ig(b1,l)-1,jg(b1,l),i,U) & ,gn(ig(b1,l),jg(b1,l),i,U), & gn(ig(b1,l)+1,jg(b1,l),i,U) & ,gn(ig(b1,l)+2,jg(b1,l),i,U)) forall(l=-nd:(b3-1)) v2(b1,l) = cubic (q2p(b1,l),q(ig(b1,l)-1), & q(ig(b1,l)),q(ig(b1,l)+1) & ,q(ig(b1,l)+2),gn(ig(b1,l)-1,jg2(b1,l),i,U) & ,gn(ig(b1,l),jg2(b1,l),i,U), & gn(ig(b1,l)+1,jg2(b1,l),i,U), & gn(ig(b1,l)+2,jg2(b1,l),i,U)) forall(l=(b4+1):nd) v1(b2,l) = cubic (q1p(b2,l),q(ig(b2,l)-1), & q(ig(b2,l)),q(ig(b2,l)+1) & ,q(ig(b2,l)+2), & gn(ig(b2,l)-1,jg(b2,l),i,U) & ,gn(ig(b2,l),jg(b2,l),i,U), & gn(ig(b2,l)+1,jg(b2,l),i,U) & ,gn(ig(b2,l)+2,jg(b2,l),i,U)) forall(l=(b4+1):nd) v2(b2,l) = cubic (q2p(b2,l),q(ig(b2,l)-1), & q(ig(b2,l)),q(ig(b2,l)+1) & ,q(ig(b2,l)+2), & gn(ig(b2,l)-1,jg2(b2,l),i,U) & ,gn(ig(b2,l),jg2(b2,l),i,U), & gn(ig(b2,l)+1,jg2(b2,l),i,U), & gn(ig(b2,l)+2,jg2(b2,l),i,U)) forall(l=-nd:(b3-1)) v1(b2,l) = cubic (q1p(b2,l),q(ig(b2,l)-1), & q(ig(b2,l)),q(ig(b2,l)+1) & ,q(ig(b2,l)+2), & gn(ig(b2,l)-1,jg(b2,l),i,U) & ,gn(ig(b2,l),jg(b2,l),i,U), & gn(ig(b2,l)+1,jg(b2,l),i,U) & ,gn(ig(b2,l)+2,jg(b2,l),i,U)) forall(l=-nd:(b3-1)) v2(b2,l) = cubic (q2p(b2,l),q(ig(b2,l)-1), & q(ig(b2,l)),q(ig(b2,l)+1) & ,q(ig(b2,l)+2),gn(ig(b2,l)-1,jg2(b2,l),i,U) & ,gn(ig(b2,l),jg2(b2,l),i,U), & gn(ig(b2,l)+1,jg2(b2,l),i,U), & gn(ig(b2,l)+2,jg2(b2,l),i,U)) c Do the upper interpolation depending on which direction is used forall(l=(b4+1):nd) gn(b1,l,i,D) = cubic(q(b1),q(k2(b1,l)), & q(k1(b1,l)),q1(b1,l),q2(b1,l) & ,gn(k2(b1,l),l,i,D), & gn(k1(b1,l),l,i,D),v1(b1,l),v2(b1,l)) forall(l=(b4+1):nd) gn(b2,l,i,D) = cubic(q(b2),q(k2(b2,l)), & q(k1(b2,l)),q1(b2,l),q2(b2,l) & ,gn(k2(b2,l),l,i,D), & gn(k1(b2,l),l,i,D),v1(b2,l),v2(b2,l)) forall(l=-nd:(b3-1)) gn(b1,l,i,D) = cubic(q(b1),q(k2(b1,l)), & q(k1(b1,l)),q1(b1,l),q2(b1,l) & ,gn(k2(b1,l),l,i,D), & gn(k1(b1,l),l,i,D),v1(b1,l),v2(b1,l)) forall(l=-nd:(b3-1)) gn(b2,l,i,D) = cubic(q(b2),q(k2(b2,l)), & q(k1(b2,l)),q1(b2,l),q2(b2,l) & ,gn(k2(b2,l),l,i,D), & gn(k1(b2,l),l,i,D),v1(b2,l),v2(b2,l)) c All the operations for mask(k,l)=21 c The guard point index and the outer neighbor index forall(k=(b4+1):nd) ig(k,b1) = int(qf(k,b1)/dd) forall(k=(b4+1):nd) jg(k,b1) = int(pf(k,b1)/dd) & + int((sign(z,pf(k,b1))-1)/2) forall(k=(b4+1):nd) ig(k,b2) = int(qf(k,b2)/dd) forall(k=(b4+1):nd) jg(k,b2) = int(pf(k,b2)/dd) & + int((sign(z,pf(k,b2))-1)/2) forall(k=-nd:(b3-1)) ig(k,b1) = int(qf(k,b1)/dd) forall(k=-nd:(b3-1)) jg(k,b1) = int(pf(k,b1)/dd) & + int((sign(z,pf(k,b1))-1)/2) forall(k=-nd:(b3-1)) ig(k,b2) = int(qf(k,b2)/dd) forall(k=-nd:(b3-1)) jg(k,b2) = int(pf(k,b2)/dd) & + int((sign(z,pf(k,b2))-1)/2) forall(k=(b4+1):nd) ig2(k,b1) = ig(k,b1) & - int(sign(z,q(ig(k,b1)))) forall(k=(b4+1):nd) jg2(k,b1) = jg(k,b1) & - int(sign(z,p(jg(k,b1)))) forall(k=(b4+1):nd) ig2(k,b2) = ig(k,b2) & - int(sign(z,q(ig(k,b2)))) forall(k=(b4+1):nd) jg2(k,b2) = jg(k,b2) & - int(sign(z,p(jg(k,b2)))) forall(k=-nd:(b3-1)) ig2(k,b1) = ig(k,b1) & - int(sign(z,q(ig(k,b1)))) forall(k=-nd:(b3-1)) jg2(k,b1) = jg(k,b1) & - int(sign(z,p(jg(k,b1)))) forall(k=-nd:(b3-1)) ig2(k,b2) = ig(k,b2) & - int(sign(z,q(ig(k,b2)))) forall(k=-nd:(b3-1)) jg2(k,b2) = jg(k,b2) & - int(sign(z,p(jg(k,b2)))) c Find the upper grid points used in the cubic interpolation forall(k=(b4+1):nd) l1(k,b1) = b1 - 1 * int(sign(z,p(b1))) forall(k=(b4+1):nd) l2(k,b1) = b1 - 2 * int(sign(z,p(b1))) forall(k=(b4+1):nd) l1(k,b2) = b2 - 1 * int(sign(z,p(b2))) forall(k=(b4+1):nd) l2(k,b2) = b2 - 2 * int(sign(z,p(b2))) forall(k=-nd:(b3-1)) l1(k,b1) = b1 - 1 * int(sign(z,p(b1))) forall(k=-nd:(b3-1)) l2(k,b1) = b1 - 2 * int(sign(z,p(b1))) forall(k=-nd:(b3-1)) l1(k,b2) = b2 - 1 * int(sign(z,p(b2))) forall(k=-nd:(b3-1)) l2(k,b2) = b2 - 2 * int(sign(z,p(b2))) c Calculate the intersection points depending on which combination of crossing c lines we use forall(k=(b4+1):nd) p1(k,b1) = sign(z,p(b1)) * sqrt(q(k) & /q(ig(k,b1)) - q(k)**2) forall(k=(b4+1):nd) p1p(k,b1) = -(q(ig(k,b1)) & /q(k)) * p1(k,b1) forall(k=(b4+1):nd) p2(k,b1) = sign(z,p(b1)) * sqrt(q(k) & /q(ig2(k,b1)) - q(k)**2) forall(k=(b4+1):nd) p2p(k,b1) = -(q(ig2(k,b1))/q(k)) * p2(k,b1) forall(k=(b4+1):nd) p1(k,b2) = sign(z,p(b2)) * sqrt(q(k) & /q(ig(k,b2)) - q(k)**2) forall(k=(b4+1):nd) p1p(k,b2) = -(q(ig(k,b2))/q(k)) & * p1(k,b2) forall(k=(b4+1):nd) p2(k,b2) = sign(z,p(b2)) & * sqrt(q(k)/q(ig2(k,b2)) - q(k)**2) forall(k=(b4+1):nd) p2p(k,b2) = -(q(ig2(k,b2))/q(k)) * p2(k,b2) forall(k=-nd:(b3-1)) p1(k,b1) = sign(z,p(b1)) & * sqrt(q(k)/q(ig(k,b1)) - q(k)**2) forall(k=-nd:(b3-1)) p1p(k,b1) = -(q(ig(k,b1))/q(k)) * p1(k,b1) forall(k=-nd:(b3-1)) p2(k,b1) = sign(z,p(b1)) & * sqrt(q(k)/q(ig2(k,b1)) - q(k)**2) forall(k=-nd:(b3-1)) p2p(k,b1) = -(q(ig2(k,b1))/q(k)) * p2(k,b1) forall(k=-nd:(b3-1)) p1(k,b2) = sign(z,p(b2)) & * sqrt(q(k)/q(ig(k,b2)) - q(k)**2) forall(k=-nd:(b3-1)) p1p(k,b2) = -(q(ig(k,b2))/q(k)) * p1(k,b2) forall(k=-nd:(b3-1)) p2(k,b2) = sign(z,p(b2)) & * sqrt(q(k)/q(ig2(k,b2)) - q(k)**2) forall(k=-nd:(b3-1)) p2p(k,b2) = -(q(ig2(k,b2))/q(k)) * p2(k,b2) c Do the lower interpolation depending on which direction is used forall(k=(b4+1):nd) v1(k,b1) = cubic(p1p(k,b1),p(jg(k,b1)-1), & p(jg(k,b1)),p(jg(k,b1)+1) & ,p(jg(k,b1)+2), & gn(ig(k,b1),jg(k,b1)-1,i,U), & gn(ig(k,b1),jg(k,b1),i,U), & gn(ig(k,b1),jg(k,b1)+1,i,U), & gn(ig(k,b1),jg(k,b1)+2,i,U)) forall(k=(b4+1):nd) v2(k,b1) = cubic(p2p(k,b1),p(jg(k,b1)-1), & p(jg(k,b1)),p(jg(k,b1)+1) & ,p(jg(k,b1)+2), & gn(ig2(k,b1),jg(k,b1)-1,i,U), & gn(ig2(k,b1),jg(k,b1),i,U), & gn(ig2(k,b1),jg(k,b1)+1,i,U), & gn(ig2(k,b1),jg(k,b1)+2,i,U)) forall(k=(b4+1):nd) v1(k,b2) = cubic(p1p(k,b2),p(jg(k,b2)-1), & p(jg(k,b2)),p(jg(k,b2)+1) & ,p(jg(k,b2)+2), & gn(ig(k,b2),jg(k,b2)-1,i,U), & gn(ig(k,b2),jg(k,b2),i,U), & gn(ig(k,b2),jg(k,b2)+1,i,U), & gn(ig(k,b2),jg(k,b2)+2,i,U)) forall(k=(b4+1):nd) v2(k,b2) = cubic(p2p(k,b2),p(jg(k,b2)-1), & p(jg(k,b2)),p(jg(k,b2)+1) & ,p(jg(k,b2)+2), & gn(ig2(k,b2),jg(k,b2)-1,i,U), & gn(ig2(k,b2),jg(k,b2),i,U), & gn(ig2(k,b2),jg(k,b2)+1,i,U), & gn(ig2(k,b2),jg(k,b2)+2,i,U)) forall(k=-nd:(b3-1)) v1(k,b1) = cubic(p1p(k,b1),p(jg(k,b1)-1), & p(jg(k,b1)),p(jg(k,b1)+1) & ,p(jg(k,b1)+2), & gn(ig(k,b1),jg(k,b1)-1,i,U), & gn(ig(k,b1),jg(k,b1),i,U), & gn(ig(k,b1),jg(k,b1)+1,i,U), & gn(ig(k,b1),jg(k,b1)+2,i,U)) forall(k=-nd:(b3-1)) v2(k,b1) = cubic(p2p(k,b1),p(jg(k,b1)-1), & p(jg(k,b1)),p(jg(k,b1)+1) & ,p(jg(k,b1)+2), & gn(ig2(k,b1),jg(k,b1)-1,i,U), & gn(ig2(k,b1),jg(k,b1),i,U), & gn(ig2(k,b1),jg(k,b1)+1,i,U), & gn(ig2(k,b1),jg(k,b1)+2,i,U)) forall(k=-nd:(b3-1)) v1(k,b2) = cubic(p1p(k,b2),p(jg(k,b2)-1), & p(jg(k,b2)),p(jg(k,b2)+1) & ,p(jg(k,b2)+2), & gn(ig(k,b2),jg(k,b2)-1,i,U), & gn(ig(k,b2),jg(k,b2),i,U), & gn(ig(k,b2),jg(k,b2)+1,i,U), & gn(ig(k,b2),jg(k,b2)+2,i,U)) forall(k=-nd:(b3-1)) v2(k,b2) = cubic(p2p(k,b2),p(jg(k,b2)-1), & p(jg(k,b2)),p(jg(k,b2)+1) & ,p(jg(k,b2)+2), & gn(ig2(k,b2),jg(k,b2)-1,i,U), & gn(ig2(k,b2),jg(k,b2),i,U), & gn(ig2(k,b2),jg(k,b2)+1,i,U), & gn(ig2(k,b2),jg(k,b2)+2,i,U)) c Do the upper interpolation depending on which direction is used forall(k=(b4+1):nd) gn(k,b1,i,D) = cubic(p(b1),p(l2(k,b1)), & p(l1(k,b1)),p1(k,b1) & ,p2(k,b1),gn(k,l2(k,b1),i,D) & ,gn(k,l1(k,b1),i,D) & ,v1(k,b1),v2(k,b1)) forall(k=(b4+1):nd) gn(k,b2,i,D) = cubic(p(b2),p(l2(k,b2)) & ,p(l1(k,b2)),p1(k,b2) & ,p2(k,b2),gn(k,l2(k,b2),i,D) & ,gn(k,l1(k,b2),i,D) & ,v1(k,b2),v2(k,b2)) forall(k=-nd:(b3-1)) gn(k,b1,i,D) = cubic(p(b1),p(l2(k,b1)) & ,p(l1(k,b1)),p1(k,b1) & ,p2(k,b1),gn(k,l2(k,b1),i,D) & ,gn(k,l1(k,b1),i,D) & ,v1(k,b1),v2(k,b1)) forall(k=-nd:(b3-1)) gn(k,b2,i,D) = cubic(p(b2),p(l2(k,b2)) & ,p(l1(k,b2)),p1(k,b2) & ,p2(k,b2),gn(k,l2(k,b2),i,D) & ,gn(k,l1(k,b2),i,D) & ,v1(k,b2),v2(k,b2)) c All the operations for mask(k,l)=22 c The guard point index and the outer neighbor index forall(k=b3:b4) ig(k,b1) = int(qf(k,b1)/dd) & + int((sign(z,qf(k,b1)) - 1)/2) forall(k=b3:b4) jg(k,b1) = int(pf(k,b1)/dd) forall(k=b3:b4) ig(k,b2) = int(qf(k,b2)/dd) & + int((sign(z,qf(k,b2)) - 1)/2) forall(k=b3:b4) jg(k,b2) = int(pf(k,b2)/dd) forall(k=b3:b4) ig2(k,b1) = ig(k,b1) - int(sign(z,q(ig(k,b1)))) forall(k=b3:b4) jg2(k,b1) = jg(k,b1) - int(sign(z,p(jg(k,b1)))) forall(k=b3:b4) ig2(k,b2) = ig(k,b2) - int(sign(z,q(ig(k,b2)))) forall(k=b3:b4) jg2(k,b2) = jg(k,b2) - int(sign(z,p(jg(k,b2)))) c Find the upper grid points used in the cubic interpolation forall(k=b3:b4) l1(k,b1) = b1 - 1 * int(sign(z,p(b1))) forall(k=b3:b4) l2(k,b1) = b1 - 2 * int(sign(z,p(b1))) forall(k=b3:b4) l1(k,b2) = b2 - 1 * int(sign(z,p(b2))) forall(k=b3:b4) l2(k,b2) = b2 - 2 * int(sign(z,p(b2))) c Calculate the intersection points depending on which combination of c crossing lines we use forall(k=b3:b4) p1(k,b1) = (-1. - sqrt(1. - 4. * q(k)**2 & * (p(jg(k,b1)))**2))/ & (2.* p(jg(k,b1))) forall(k=b3:b4) q1p(k,b1) = 2.*q(k)*(p(jg(k,b1)))**2 & /(1. + sqrt(1. - 4. * q(k) & **2 * (p(jg(k,b1)))**2)) forall(k=b3:b4) p2(k,b1) = (-1. - sqrt(1. - 4. * q(k)**2 & * p(jg2(k,b1))**2))/ & (2.*p(jg2(k,b1))) forall(k=b3:b4) q2p(k,b1) = 2.*q(k)*(p(jg2(k,b1)))**2 & /(1. +sqrt(1. - 4. * & q(k)**2 * (p(jg2(k,b1)))**2)) forall(k=b3:b4) p1(k,b2) = (-1. - sqrt(1. - 4. * q(k)**2 & * p(jg(k,b2))**2))/ & (2.* p(jg(k,b2))) forall(k=b3:b4) q1p(k,b2) = 2.*q(k)*(p(jg(k,b2)))**2 & /(1. + sqrt(1. - 4. * q(k) & **2 * p(jg(k,b2))**2)) forall(k=b3:b4) p2(k,b2) = (-1. - sqrt(1. - 4. * q(k)**2 & * p(jg2(k,b2))**2))/ & (2.*p(jg2(k,b2))) forall(k=b3:b4) q2p(k,b2) = 2.*q(k)*p(jg2(k,b2))**2 & /(1. +sqrt(1. - 4. * & q(k)**2 * p(jg2(k,b2))**2)) c Do the lower interpolation depending on which direction we use forall(k=b3:b4) v1(k,b1) = cubic(q1p(k,b1), q(ig(k,b1)-1), & q(ig(k,b1)) & , q(ig(k,b1)+1),q(ig(k,b1)+2) & ,gn(ig(k,b1)-1,jg(k,b1),i,U), & gn(ig(k,b1),jg(k,b1),i,U) & ,gn(ig(k,b1)+1,jg(k,b1),i,U), & gn(ig(k,b1)+2,jg(k,b1),i,U)) forall(k=b3:b4) v1(k,b2) = cubic(q1p(k,b2), q(ig(k,b2)-1), & q(ig(k,b2)) & , q(ig(k,b2)+1),q(ig(k,b2)+2) & ,gn(ig(k,b2)-1,jg(k,b2),i,U), & gn(ig(k,b2),jg(k,b2),i,U) & ,gn(ig(k,b2)+1,jg(k,b2),i,U) & ,gn(ig(k,b2)+2,jg(k,b2),i,U)) forall(k=b3:b4) v2(k,b1) = cubic(q2p(k,b1), q(ig(k,b1)-1), & q(ig(k,b1)) & , q(ig(k,b1)+1),q(ig(k,b1)+2) & ,gn(ig(k,b1)-1,jg2(k,b1),i,U), & gn(ig(k,b1),jg2(k,b1),i,U) & ,gn(ig(k,b1)+1,jg2(k,b1),i,U), & gn(ig(k,b1)+2,jg2(k,b1),i,U)) forall(k=b3:b4) v2(k,b2) = cubic(q2p(k,b2), q(ig(k,b2)-1), & q(ig(k,b2)) & , q(ig(k,b2)+1),q(ig(k,b2)+2) & ,gn(ig(k,b2)-1,jg2(k,b2),i,U), & gn(ig(k,b2),jg2(k,b2),i,U) & ,gn(ig(k,b2)+1,jg2(k,b2),i,U), & gn(ig(k,b2)+2,jg2(k,b2),i,U)) c Do the upper interpolation depending on which direction is used forall(k=b3:b4) gn(k,b1,i,D) = cubic(p(b1), p(l2(k,b1)) & , p(l1(k,b1)),p1(k,b1) & ,p2(k,b1),gn(k,l2(k,b1),i,D), & gn(k,l1(k,b1),i,D) & ,v1(k,b1),v2(k,b1)) forall(k=b3:b4) gn(k,b2,i,D) = cubic(p(b2), p(l2(k,b2)), & p(l1(k,b2)),p1(k,b2) & ,p2(k,b2),gn(k,l2(k,b2),i,D), & gn(k,l1(k,b2),i,D) & ,v1(k,b2),v2(k,b2)) c All the operations for mask(k,l)=3 c guard point index and the outer neighbor index ig(b1,b1)=int(qf(b1,b1)/dd) ig(b1,b2)=int(qf(b1,b2)/dd) ig(b2,b1)=int(qf(b2,b1)/dd) ig(b2,b2)=int(qf(b2,b2)/dd) jg(b1,b1)=int(pf(b1,b1)/dd) jg(b1,b2)=int(pf(b1,b2)/dd) jg(b2,b1)=int(pf(b2,b1)/dd) jg(b2,b2)=int(pf(b2,b2)/dd) ig2(b1,b1)= ig(b1,b1) - int(sign (z,q(ig(b1,b1)))) ig2(b1,b2)= ig(b1,b2) - int(sign (z,q(ig(b1,b2)))) ig2(b2,b1)= ig(b2,b1) - int(sign (z,q(ig(b2,b1)))) ig2(b2,b2)= ig(b2,b2) - int(sign (z,q(ig(b2,b2)))) jg2(b1,b1)= jg(b1,b1) - int(sign (z,p(jg(b1,b1)))) jg2(b1,b2)= jg(b1,b2) - int(sign (z,p(jg(b1,b2)))) jg2(b2,b1)= jg(b2,b1) - int(sign (z,p(jg(b2,b1)))) jg2(b2,b2)= jg(b2,b2) - int(sign (z,p(jg(b2,b2)))) c Find the upper grid points used in the cubic interpolation k1(b1,b1) = b1 - 1 * int(sign(z,q(b1))) k1(b1,b2) = b1 - 1 * int(sign(z,q(b1))) k1(b2,b1) = b2 - 1 * int(sign(z,q(b2))) k1(b2,b2) = b2 - 1 * int(sign(z,q(b2))) k2(b1,b1) = b1 - 2 * int(sign(z,q(b1))) k2(b1,b2) = b1 - 2 * int(sign(z,q(b1))) k2(b2,b1) = b2 - 2 * int(sign(z,q(b2))) k2(b2,b2) = b2 - 2 * int(sign(z,q(b2))) l1(b1,b1) = b1 - 1 * int(sign(z,p(b1))) l1(b1,b2) = b2 - 1 * int(sign(z,p(b2))) l1(b2,b1) = b1 - 1 * int(sign(z,p(b1))) l1(b2,b2) = b2 - 1 * int(sign(z,p(b2))) l2(b1,b1) = b1 - 2 * int(sign(z,p(b1))) l2(b1,b2) = b2 - 2 * int(sign(z,p(b2))) l2(b2,b1) = b1 - 2 * int(sign(z,p(b1))) l2(b2,b2) = b2 - 2 * int(sign(z,p(b2))) c Do the lower interpolation depending on which direction is used v1(b1,b1) = gn(ig(b1,b1),jg(b1,b1),i,U) v1(b1,b2) = gn(ig(b1,b2),jg(b1,b2),i,U) v1(b2,b1) = gn(ig(b2,b1),jg(b2,b1),i,U) v1(b2,b2) = gn(ig(b2,b2),jg(b2,b2),i,U) v2(b1,b1) = gn(ig2(b1,b1),jg2(b1,b1),i,U) v2(b1,b2) = gn(ig2(b1,b2),jg2(b1,b2),i,U) v2(b2,b1) = gn(ig2(b2,b1),jg2(b2,b1),i,U) v2(b2,b2) = gn(ig2(b2,b2),jg2(b2,b2),i,U) c Do the upper interpolation depending on which direction is used qd1(b1,b1) = q_s2n(q(ig(b1,b1)),p(jg(b1,b1))) qd1(b1,b2) = q_s2n(q(ig(b1,b2)),p(jg(b1,b2))) qd1(b2,b1) = q_s2n(q(ig(b2,b1)),p(jg(b2,b1))) qd1(b2,b2) = q_s2n(q(ig(b2,b2)),p(jg(b2,b2))) pd1(b1,b1) = p_s2n(q(ig(b1,b1)),p(jg(b1,b1))) pd1(b1,b2) = p_s2n(q(ig(b1,b2)),p(jg(b1,b2))) pd1(b2,b1) = p_s2n(q(ig(b2,b1)),p(jg(b2,b1))) pd1(b2,b2) = p_s2n(q(ig(b2,b2)),p(jg(b2,b2))) qd2(b1,b1) = q_s2n(q(ig2(b1,b1)),p(jg2(b1,b1))) qd2(b1,b2) = q_s2n(q(ig2(b1,b2)),p(jg2(b1,b2))) qd2(b2,b1) = q_s2n(q(ig2(b2,b1)),p(jg2(b2,b1))) qd2(b2,b2) = q_s2n(q(ig2(b2,b2)),p(jg2(b2,b2))) pd2(b1,b1) = p_s2n(q(ig2(b1,b1)),p(jg2(b1,b1))) pd2(b1,b2) = p_s2n(q(ig2(b1,b2)),p(jg2(b1,b2))) pd2(b2,b1) = p_s2n(q(ig2(b2,b1)),p(jg2(b2,b1))) pd2(b2,b2) = p_s2n(q(ig2(b2,b2)),p(jg2(b2,b2))) r1(b1,b1) = sqrt(q(k1(b1,b1))**2+p(l1(b1,b1))**2) r1(b1,b2) = sqrt(q(k1(b1,b2))**2+p(l1(b1,b2))**2) r1(b2,b1) = sqrt(q(k1(b2,b1))**2+p(l1(b2,b1))**2) r1(b2,b2) = sqrt(q(k1(b2,b2))**2+p(l1(b2,b2))**2) r2(b1,b1) = sqrt(q(k2(b1,b1))**2+p(l2(b1,b1))**2) r2(b1,b2) = sqrt(q(k2(b1,b2))**2+p(l2(b1,b2))**2) r2(b2,b1) = sqrt(q(k2(b2,b1))**2+p(l2(b2,b1))**2) r2(b2,b2) = sqrt(q(k2(b2,b2))**2+p(l2(b2,b2))**2) rx(b1,b1) = sqrt(q(b1)**2+p(b1)**2) rx(b1,b2) = sqrt(q(b1)**2+p(b2)**2) rx(b2,b1) = sqrt(q(b2)**2+p(b1)**2) rx(b2,b2) = sqrt(q(b2)**2+p(b2)**2) r3(b1,b1) = sqrt(qd1(b1,b1)**2 + pd1(b1,b1)**2) r3(b1,b2) = sqrt(qd1(b1,b2)**2 + pd1(b1,b2)**2) r3(b2,b1) = sqrt(qd1(b2,b1)**2 + pd1(b2,b1)**2) r3(b2,b2) = sqrt(qd1(b2,b2)**2 + pd1(b2,b2)**2) r4(b1,b1) = sqrt(qd2(b1,b1)**2 + pd2(b1,b1)**2) r4(b1,b2) = sqrt(qd2(b1,b2)**2 + pd2(b1,b2)**2) r4(b2,b1) = sqrt(qd2(b2,b1)**2 + pd2(b2,b1)**2) r4(b2,b2) = sqrt(qd2(b2,b2)**2 + pd2(b2,b2)**2) gn(b1,b1,i,D) = cubic(rx(b1,b1),r1(b1,b1), & r2(b1,b1),r3(b1,b1),r4(b1,b1) & ,gn(k1(b1,b1),l1(b1,b1),i,D), & gn(k2(b1,b1),l2(b1,b1),i,D) & ,v1(b1,b1),v2(b1,b1)) gn(b1,b2,i,D) = cubic(rx(b1,b2),r1(b1,b2), & r2(b1,b2),r3(b1,b2),r4(b1,b2) & ,gn(k1(b1,b2),l1(b1,b2),i,D), & gn(k2(b1,b2),l2(b1,b2),i,D) & ,v1(b1,b2),v2(b1,b2)) gn(b2,b1,i,D) = cubic(rx(b2,b1),r1(b2,b1), & r2(b2,b1),r3(b2,b1),r4(b2,b1) & ,gn(k1(b2,b1),l1(b2,b1),i,D), & gn(k2(b2,b1),l2(b2,b1),i,D) & ,v1(b2,b1),v2(b2,b1)) gn(b2,b2,i,D) = cubic(rx(b2,b2),r1(b2,b2), & r2(b2,b2),r3(b2,b2),r4(b2,b2) & ,gn(k1(b2,b2),l1(b2,b2),i,D), & gn(k2(b2,b2),l2(b2,b2),i,D) & ,v1(b2,b2),v2(b2,b2)) if(i.eq.1) then l=b2 do k=-nd,b3-1 write(24,*)k,l,i,D,v1(k,l),v2(k,l) enddo endif end do return end