c--------------------------------------------------------------------- c--------------------------------------------------------------------- subroutine pintgr c--------------------------------------------------------------------- c--------------------------------------------------------------------- implicit none include 'mpinpb.h' include 'applu.incl' c--------------------------------------------------------------------- c local variables c--------------------------------------------------------------------- integer i, j, k integer ibeg, ifin, ifin1 integer jbeg, jfin, jfin1 integer iglob, iglob1, iglob2 integer jglob, jglob1, jglob2 integer ind1, ind2 double precision phi1(0:isiz2+1,0:isiz3+1), > phi2(0:isiz2+1,0:isiz3+1) double precision frc1, frc2, frc3 double precision dummy integer IERROR c--------------------------------------------------------------------- c set up the sub-domains for integeration in each processor c--------------------------------------------------------------------- ibeg = nx + 1 ifin = 0 iglob1 = ipt + 1 iglob2 = ipt + nx if (iglob1.ge.ii1.and.iglob2.lt.ii2+nx) ibeg = 1 if (iglob1.gt.ii1-nx.and.iglob2.le.ii2) ifin = nx if (ii1.ge.iglob1.and.ii1.le.iglob2) ibeg = ii1 - ipt if (ii2.ge.iglob1.and.ii2.le.iglob2) ifin = ii2 - ipt jbeg = ny + 1 jfin = 0 jglob1 = jpt + 1 jglob2 = jpt + ny if (jglob1.ge.ji1.and.jglob2.lt.ji2+ny) jbeg = 1 if (jglob1.gt.ji1-ny.and.jglob2.le.ji2) jfin = ny if (ji1.ge.jglob1.and.ji1.le.jglob2) jbeg = ji1 - jpt if (ji2.ge.jglob1.and.ji2.le.jglob2) jfin = ji2 - jpt ifin1 = ifin jfin1 = jfin if (ipt + ifin1.eq.ii2) ifin1 = ifin -1 if (jpt + jfin1.eq.ji2) jfin1 = jfin -1 c--------------------------------------------------------------------- c initialize c--------------------------------------------------------------------- do i = 0,isiz2+1 do k = 0,isiz3+1 phi1(i,k) = 0. phi2(i,k) = 0. end do end do do j = jbeg,jfin jglob = jpt + j do i = ibeg,ifin iglob = ipt + i k = ki1 phi1(i,j) = c2*( u(5,i,j,k) > - 0.50d+00 * ( u(2,i,j,k) ** 2 > + u(3,i,j,k) ** 2 > + u(4,i,j,k) ** 2 ) > / u(1,i,j,k) ) k = ki2 phi2(i,j) = c2*( u(5,i,j,k) > - 0.50d+00 * ( u(2,i,j,k) ** 2 > + u(3,i,j,k) ** 2 > + u(4,i,j,k) ** 2 ) > / u(1,i,j,k) ) end do end do c--------------------------------------------------------------------- c communicate in i and j directions c--------------------------------------------------------------------- call exchange_4(phi1,phi2,ibeg,ifin1,jbeg,jfin1) frc1 = 0.0d+00 do j = jbeg,jfin1 do i = ibeg, ifin1 frc1 = frc1 + ( phi1(i,j) > + phi1(i+1,j) > + phi1(i,j+1) > + phi1(i+1,j+1) > + phi2(i,j) > + phi2(i+1,j) > + phi2(i,j+1) > + phi2(i+1,j+1) ) end do end do c--------------------------------------------------------------------- c compute the global sum of individual contributions to frc1 c--------------------------------------------------------------------- dummy = frc1 call MPI_ALLREDUCE( dummy, > frc1, > 1, > dp_type, > MPI_SUM, > MPI_COMM_WORLD, > IERROR ) frc1 = dxi * deta * frc1 c--------------------------------------------------------------------- c initialize c--------------------------------------------------------------------- do i = 0,isiz2+1 do k = 0,isiz3+1 phi1(i,k) = 0. phi2(i,k) = 0. end do end do jglob = jpt + jbeg ind1 = 0 if (jglob.eq.ji1) then ind1 = 1 do k = ki1, ki2 do i = ibeg, ifin iglob = ipt + i phi1(i,k) = c2*( u(5,i,jbeg,k) > - 0.50d+00 * ( u(2,i,jbeg,k) ** 2 > + u(3,i,jbeg,k) ** 2 > + u(4,i,jbeg,k) ** 2 ) > / u(1,i,jbeg,k) ) end do end do end if jglob = jpt + jfin ind2 = 0 if (jglob.eq.ji2) then ind2 = 1 do k = ki1, ki2 do i = ibeg, ifin iglob = ipt + i phi2(i,k) = c2*( u(5,i,jfin,k) > - 0.50d+00 * ( u(2,i,jfin,k) ** 2 > + u(3,i,jfin,k) ** 2 > + u(4,i,jfin,k) ** 2 ) > / u(1,i,jfin,k) ) end do end do end if c--------------------------------------------------------------------- c communicate in i direction c--------------------------------------------------------------------- if (ind1.eq.1) then call exchange_5(phi1,ibeg,ifin1) end if if (ind2.eq.1) then call exchange_5 (phi2,ibeg,ifin1) end if frc2 = 0.0d+00 do k = ki1, ki2-1 do i = ibeg, ifin1 frc2 = frc2 + ( phi1(i,k) > + phi1(i+1,k) > + phi1(i,k+1) > + phi1(i+1,k+1) > + phi2(i,k) > + phi2(i+1,k) > + phi2(i,k+1) > + phi2(i+1,k+1) ) end do end do c--------------------------------------------------------------------- c compute the global sum of individual contributions to frc2 c--------------------------------------------------------------------- dummy = frc2 call MPI_ALLREDUCE( dummy, > frc2, > 1, > dp_type, > MPI_SUM, > MPI_COMM_WORLD, > IERROR ) frc2 = dxi * dzeta * frc2 c--------------------------------------------------------------------- c initialize c--------------------------------------------------------------------- do i = 0,isiz2+1 do k = 0,isiz3+1 phi1(i,k) = 0. phi2(i,k) = 0. end do end do iglob = ipt + ibeg ind1 = 0 if (iglob.eq.ii1) then ind1 = 1 do k = ki1, ki2 do j = jbeg, jfin jglob = jpt + j phi1(j,k) = c2*( u(5,ibeg,j,k) > - 0.50d+00 * ( u(2,ibeg,j,k) ** 2 > + u(3,ibeg,j,k) ** 2 > + u(4,ibeg,j,k) ** 2 ) > / u(1,ibeg,j,k) ) end do end do end if iglob = ipt + ifin ind2 = 0 if (iglob.eq.ii2) then ind2 = 1 do k = ki1, ki2 do j = jbeg, jfin jglob = jpt + j phi2(j,k) = c2*( u(5,ifin,j,k) > - 0.50d+00 * ( u(2,ifin,j,k) ** 2 > + u(3,ifin,j,k) ** 2 > + u(4,ifin,j,k) ** 2 ) > / u(1,ifin,j,k) ) end do end do end if c--------------------------------------------------------------------- c communicate in j direction c--------------------------------------------------------------------- if (ind1.eq.1) then call exchange_6(phi1,jbeg,jfin1) end if if (ind2.eq.1) then call exchange_6(phi2,jbeg,jfin1) end if frc3 = 0.0d+00 do k = ki1, ki2-1 do j = jbeg, jfin1 frc3 = frc3 + ( phi1(j,k) > + phi1(j+1,k) > + phi1(j,k+1) > + phi1(j+1,k+1) > + phi2(j,k) > + phi2(j+1,k) > + phi2(j,k+1) > + phi2(j+1,k+1) ) end do end do c--------------------------------------------------------------------- c compute the global sum of individual contributions to frc3 c--------------------------------------------------------------------- dummy = frc3 call MPI_ALLREDUCE( dummy, > frc3, > 1, > dp_type, > MPI_SUM, > MPI_COMM_WORLD, > IERROR ) frc3 = deta * dzeta * frc3 frc = 0.25d+00 * ( frc1 + frc2 + frc3 ) c if (id.eq.0) write (*,1001) frc return 1001 format (//5x,'surface integral = ',1pe12.5//) end