examples/smpi/NAS/LU/pintgr.f

   1
   2 c---------------------------------------------------------------------
   3 c---------------------------------------------------------------------
   4
   5       subroutine pintgr
   6
   7 c---------------------------------------------------------------------
   8 c---------------------------------------------------------------------
   9
  10       implicit none
  11
  12       include 'mpinpb.h'
  13       include 'applu.incl'
  14
  15 c---------------------------------------------------------------------
  16 c  local variables
  17 c---------------------------------------------------------------------
  18       integer i, j, k
  19       integer ibeg, ifin, ifin1
  20       integer jbeg, jfin, jfin1
  21       integer iglob, iglob1, iglob2
  22       integer jglob, jglob1, jglob2
  23       integer ind1, ind2
  24       double precision  phi1(0:isiz2+1,0:isiz3+1),
  25      >                  phi2(0:isiz2+1,0:isiz3+1)
  26       double precision  frc1, frc2, frc3
  27       double precision  dummy
  28
  29       integer IERROR
  30
  31
  32 c---------------------------------------------------------------------
  33 c   set up the sub-domains for integeration in each processor
  34 c---------------------------------------------------------------------
  35       ibeg = nx + 1
  36       ifin = 0
  37       iglob1 = ipt + 1
  38       iglob2 = ipt + nx
  39       if (iglob1.ge.ii1.and.iglob2.lt.ii2+nx) ibeg = 1
  40       if (iglob1.gt.ii1-nx.and.iglob2.le.ii2) ifin = nx
  41       if (ii1.ge.iglob1.and.ii1.le.iglob2) ibeg = ii1 - ipt
  42       if (ii2.ge.iglob1.and.ii2.le.iglob2) ifin = ii2 - ipt
  43       jbeg = ny + 1
  44       jfin = 0
  45       jglob1 = jpt + 1
  46       jglob2 = jpt + ny
  47       if (jglob1.ge.ji1.and.jglob2.lt.ji2+ny) jbeg = 1
  48       if (jglob1.gt.ji1-ny.and.jglob2.le.ji2) jfin = ny
  49       if (ji1.ge.jglob1.and.ji1.le.jglob2) jbeg = ji1 - jpt
  50       if (ji2.ge.jglob1.and.ji2.le.jglob2) jfin = ji2 - jpt
  51       ifin1 = ifin
  52       jfin1 = jfin
  53       if (ipt + ifin1.eq.ii2) ifin1 = ifin -1
  54       if (jpt + jfin1.eq.ji2) jfin1 = jfin -1
  55
  56 c---------------------------------------------------------------------
  57 c   initialize
  58 c---------------------------------------------------------------------
  59       do i = 0,isiz2+1
  60         do k = 0,isiz3+1
  61           phi1(i,k) = 0.
  62           phi2(i,k) = 0.
  63         end do
  64       end do
  65
  66       do j = jbeg,jfin
  67          jglob = jpt + j
  68          do i = ibeg,ifin
  69             iglob = ipt + i
  70
  71             k = ki1
  72
  73             phi1(i,j) = c2*(  u(5,i,j,k)
  74      >           - 0.50d+00 * (  u(2,i,j,k) ** 2
  75      >                         + u(3,i,j,k) ** 2
  76      >                         + u(4,i,j,k) ** 2 )
  77      >                        / u(1,i,j,k) )
  78
  79             k = ki2
  80
  81             phi2(i,j) = c2*(  u(5,i,j,k)
  82      >           - 0.50d+00 * (  u(2,i,j,k) ** 2
  83      >                         + u(3,i,j,k) ** 2
  84      >                         + u(4,i,j,k) ** 2 )
  85      >                        / u(1,i,j,k) )
  86          end do
  87       end do
  88
  89 c---------------------------------------------------------------------
  90 c  communicate in i and j directions
  91 c---------------------------------------------------------------------
  92       call exchange_4(phi1,phi2,ibeg,ifin1,jbeg,jfin1)
  93
  94       frc1 = 0.0d+00
  95
  96       do j = jbeg,jfin1
  97          do i = ibeg, ifin1
  98             frc1 = frc1 + (  phi1(i,j)
  99      >                     + phi1(i+1,j)
 100      >                     + phi1(i,j+1)
 101      >                     + phi1(i+1,j+1)
 102      >                     + phi2(i,j)
 103      >                     + phi2(i+1,j)
 104      >                     + phi2(i,j+1)
 105      >                     + phi2(i+1,j+1) )
 106          end do
 107       end do
 108
 109 c---------------------------------------------------------------------
 110 c  compute the global sum of individual contributions to frc1
 111 c---------------------------------------------------------------------
 112       dummy = frc1
 113       call MPI_ALLREDUCE( dummy,
 114      >                    frc1,
 115      >                    1,
 116      >                    dp_type,
 117      >                    MPI_SUM,
 118      >                    MPI_COMM_WORLD,
 119      >                    IERROR )
 120
 121       frc1 = dxi * deta * frc1
 122
 123 c---------------------------------------------------------------------
 124 c   initialize
 125 c---------------------------------------------------------------------
 126       do i = 0,isiz2+1
 127         do k = 0,isiz3+1
 128           phi1(i,k) = 0.
 129           phi2(i,k) = 0.
 130         end do
 131       end do
 132       jglob = jpt + jbeg
 133       ind1 = 0
 134       if (jglob.eq.ji1) then
 135         ind1 = 1
 136         do k = ki1, ki2
 137            do i = ibeg, ifin
 138               iglob = ipt + i
 139               phi1(i,k) = c2*(  u(5,i,jbeg,k)
 140      >             - 0.50d+00 * (  u(2,i,jbeg,k) ** 2
 141      >                           + u(3,i,jbeg,k) ** 2
 142      >                           + u(4,i,jbeg,k) ** 2 )
 143      >                          / u(1,i,jbeg,k) )
 144            end do
 145         end do
 146       end if
 147
 148       jglob = jpt + jfin
 149       ind2 = 0
 150       if (jglob.eq.ji2) then
 151         ind2 = 1
 152         do k = ki1, ki2
 153            do i = ibeg, ifin
 154               iglob = ipt + i
 155               phi2(i,k) = c2*(  u(5,i,jfin,k)
 156      >             - 0.50d+00 * (  u(2,i,jfin,k) ** 2
 157      >                           + u(3,i,jfin,k) ** 2
 158      >                           + u(4,i,jfin,k) ** 2 )
 159      >                          / u(1,i,jfin,k) )
 160            end do
 161         end do
 162       end if
 163
 164 c---------------------------------------------------------------------
 165 c  communicate in i direction
 166 c---------------------------------------------------------------------
 167       if (ind1.eq.1) then
 168         call exchange_5(phi1,ibeg,ifin1)
 169       end if
 170       if (ind2.eq.1) then
 171         call exchange_5 (phi2,ibeg,ifin1)
 172       end if
 173
 174       frc2 = 0.0d+00
 175       do k = ki1, ki2-1
 176          do i = ibeg, ifin1
 177             frc2 = frc2 + (  phi1(i,k)
 178      >                     + phi1(i+1,k)
 179      >                     + phi1(i,k+1)
 180      >                     + phi1(i+1,k+1)
 181      >                     + phi2(i,k)
 182      >                     + phi2(i+1,k)
 183      >                     + phi2(i,k+1)
 184      >                     + phi2(i+1,k+1) )
 185          end do
 186       end do
 187
 188 c---------------------------------------------------------------------
 189 c  compute the global sum of individual contributions to frc2
 190 c---------------------------------------------------------------------
 191       dummy = frc2
 192       call MPI_ALLREDUCE( dummy,
 193      >                    frc2,
 194      >                    1,
 195      >                    dp_type,
 196      >                    MPI_SUM,
 197      >                    MPI_COMM_WORLD,
 198      >                    IERROR )
 199
 200       frc2 = dxi * dzeta * frc2
 201
 202 c---------------------------------------------------------------------
 203 c   initialize
 204 c---------------------------------------------------------------------
 205       do i = 0,isiz2+1
 206         do k = 0,isiz3+1
 207           phi1(i,k) = 0.
 208           phi2(i,k) = 0.
 209         end do
 210       end do
 211       iglob = ipt + ibeg
 212       ind1 = 0
 213       if (iglob.eq.ii1) then
 214         ind1 = 1
 215         do k = ki1, ki2
 216            do j = jbeg, jfin
 217               jglob = jpt + j
 218               phi1(j,k) = c2*(  u(5,ibeg,j,k)
 219      >             - 0.50d+00 * (  u(2,ibeg,j,k) ** 2
 220      >                           + u(3,ibeg,j,k) ** 2
 221      >                           + u(4,ibeg,j,k) ** 2 )
 222      >                          / u(1,ibeg,j,k) )
 223            end do
 224         end do
 225       end if
 226
 227       iglob = ipt + ifin
 228       ind2 = 0
 229       if (iglob.eq.ii2) then
 230         ind2 = 1
 231         do k = ki1, ki2
 232            do j = jbeg, jfin
 233               jglob = jpt + j
 234               phi2(j,k) = c2*(  u(5,ifin,j,k)
 235      >             - 0.50d+00 * (  u(2,ifin,j,k) ** 2
 236      >                           + u(3,ifin,j,k) ** 2
 237      >                           + u(4,ifin,j,k) ** 2 )
 238      >                          / u(1,ifin,j,k) )
 239            end do
 240         end do
 241       end if
 242
 243 c---------------------------------------------------------------------
 244 c  communicate in j direction
 245 c---------------------------------------------------------------------
 246       if (ind1.eq.1) then
 247         call exchange_6(phi1,jbeg,jfin1)
 248       end if
 249       if (ind2.eq.1) then
 250         call exchange_6(phi2,jbeg,jfin1)
 251       end if
 252
 253       frc3 = 0.0d+00
 254
 255       do k = ki1, ki2-1
 256          do j = jbeg, jfin1
 257             frc3 = frc3 + (  phi1(j,k)
 258      >                     + phi1(j+1,k)
 259      >                     + phi1(j,k+1)
 260      >                     + phi1(j+1,k+1)
 261      >                     + phi2(j,k)
 262      >                     + phi2(j+1,k)
 263      >                     + phi2(j,k+1)
 264      >                     + phi2(j+1,k+1) )
 265          end do
 266       end do
 267
 268 c---------------------------------------------------------------------
 269 c  compute the global sum of individual contributions to frc3
 270 c---------------------------------------------------------------------
 271       dummy = frc3
 272       call MPI_ALLREDUCE( dummy,
 273      >                    frc3,
 274      >                    1,
 275      >                    dp_type,
 276      >                    MPI_SUM,
 277      >                    MPI_COMM_WORLD,
 278      >                    IERROR )
 279
 280       frc3 = deta * dzeta * frc3
 281       frc = 0.25d+00 * ( frc1 + frc2 + frc3 )
 282 c      if (id.eq.0) write (*,1001) frc
 283
 284       return
 285
 286  1001 format (//5x,'surface integral = ',1pe12.5//)
 287
 288       end