examples/smpi/NAS/LU/rhs.f

   1
   2 c---------------------------------------------------------------------
   3 c---------------------------------------------------------------------
   4
   5       subroutine rhs
   6
   7 c---------------------------------------------------------------------
   8 c---------------------------------------------------------------------
   9
  10 c---------------------------------------------------------------------
  11 c   compute the right hand sides
  12 c---------------------------------------------------------------------
  13
  14       implicit none
  15
  16       include 'applu.incl'
  17
  18 c---------------------------------------------------------------------
  19 c  local variables
  20 c---------------------------------------------------------------------
  21       integer i, j, k, m
  22       integer iex
  23       integer L1, L2
  24       integer ist1, iend1
  25       integer jst1, jend1
  26       double precision  q
  27       double precision  u21, u31, u41
  28       double precision  tmp
  29       double precision  u21i, u31i, u41i, u51i
  30       double precision  u21j, u31j, u41j, u51j
  31       double precision  u21k, u31k, u41k, u51k
  32       double precision  u21im1, u31im1, u41im1, u51im1
  33       double precision  u21jm1, u31jm1, u41jm1, u51jm1
  34       double precision  u21km1, u31km1, u41km1, u51km1
  35
  36
  37
  38       do k = 1, nz
  39          do j = 1, ny
  40             do i = 1, nx
  41                do m = 1, 5
  42                   rsd(m,i,j,k) = - frct(m,i,j,k)
  43                end do
  44             end do
  45          end do
  46       end do
  47
  48 c---------------------------------------------------------------------
  49 c   xi-direction flux differences
  50 c---------------------------------------------------------------------
  51
  52 c---------------------------------------------------------------------
  53 c   iex = flag : iex = 0  north/south communication
  54 c              : iex = 1  east/west communication
  55 c---------------------------------------------------------------------
  56       iex   = 0
  57
  58 c---------------------------------------------------------------------
  59 c   communicate and receive/send two rows of data
  60 c---------------------------------------------------------------------
  61       call exchange_3(u,iex)
  62
  63       L1 = 0
  64       if (north.eq.-1) L1 = 1
  65       L2 = nx + 1
  66       if (south.eq.-1) L2 = nx
  67
  68       ist1 = 1
  69       iend1 = nx
  70       if (north.eq.-1) ist1 = 4
  71       if (south.eq.-1) iend1 = nx - 3
  72
  73       do k = 2, nz - 1
  74          do j = jst, jend
  75             do i = L1, L2
  76                flux(1,i,j,k) = u(2,i,j,k)
  77                u21 = u(2,i,j,k) / u(1,i,j,k)
  78
  79                q = 0.50d+00 * (  u(2,i,j,k) * u(2,i,j,k)
  80      >                         + u(3,i,j,k) * u(3,i,j,k)
  81      >                         + u(4,i,j,k) * u(4,i,j,k) )
  82      >                      / u(1,i,j,k)
  83
  84                flux(2,i,j,k) = u(2,i,j,k) * u21 + c2 *
  85      >                        ( u(5,i,j,k) - q )
  86                flux(3,i,j,k) = u(3,i,j,k) * u21
  87                flux(4,i,j,k) = u(4,i,j,k) * u21
  88                flux(5,i,j,k) = ( c1 * u(5,i,j,k) - c2 * q ) * u21
  89             end do
  90
  91             do i = ist, iend
  92                do m = 1, 5
  93                   rsd(m,i,j,k) =  rsd(m,i,j,k)
  94      >                 - tx2 * ( flux(m,i+1,j,k) - flux(m,i-1,j,k) )
  95                end do
  96             end do
  97
  98             do i = ist, L2
  99                tmp = 1.0d+00 / u(1,i,j,k)
 100
 101                u21i = tmp * u(2,i,j,k)
 102                u31i = tmp * u(3,i,j,k)
 103                u41i = tmp * u(4,i,j,k)
 104                u51i = tmp * u(5,i,j,k)
 105
 106                tmp = 1.0d+00 / u(1,i-1,j,k)
 107
 108                u21im1 = tmp * u(2,i-1,j,k)
 109                u31im1 = tmp * u(3,i-1,j,k)
 110                u41im1 = tmp * u(4,i-1,j,k)
 111                u51im1 = tmp * u(5,i-1,j,k)
 112
 113                flux(2,i,j,k) = (4.0d+00/3.0d+00) * tx3 * (u21i-u21im1)
 114                flux(3,i,j,k) = tx3 * ( u31i - u31im1 )
 115                flux(4,i,j,k) = tx3 * ( u41i - u41im1 )
 116                flux(5,i,j,k) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
 117      >              * tx3 * ( ( u21i  **2 + u31i  **2 + u41i  **2 )
 118      >                      - ( u21im1**2 + u31im1**2 + u41im1**2 ) )
 119      >              + (1.0d+00/6.0d+00)
 120      >              * tx3 * ( u21i**2 - u21im1**2 )
 121      >              + c1 * c5 * tx3 * ( u51i - u51im1 )
 122             end do
 123
 124             do i = ist, iend
 125                rsd(1,i,j,k) = rsd(1,i,j,k)
 126      >              + dx1 * tx1 * (            u(1,i-1,j,k)
 127      >                             - 2.0d+00 * u(1,i,j,k)
 128      >                             +           u(1,i+1,j,k) )
 129                rsd(2,i,j,k) = rsd(2,i,j,k)
 130      >          + tx3 * c3 * c4 * ( flux(2,i+1,j,k) - flux(2,i,j,k) )
 131      >              + dx2 * tx1 * (            u(2,i-1,j,k)
 132      >                             - 2.0d+00 * u(2,i,j,k)
 133      >                             +           u(2,i+1,j,k) )
 134                rsd(3,i,j,k) = rsd(3,i,j,k)
 135      >          + tx3 * c3 * c4 * ( flux(3,i+1,j,k) - flux(3,i,j,k) )
 136      >              + dx3 * tx1 * (            u(3,i-1,j,k)
 137      >                             - 2.0d+00 * u(3,i,j,k)
 138      >                             +           u(3,i+1,j,k) )
 139                rsd(4,i,j,k) = rsd(4,i,j,k)
 140      >          + tx3 * c3 * c4 * ( flux(4,i+1,j,k) - flux(4,i,j,k) )
 141      >              + dx4 * tx1 * (            u(4,i-1,j,k)
 142      >                             - 2.0d+00 * u(4,i,j,k)
 143      >                             +           u(4,i+1,j,k) )
 144                rsd(5,i,j,k) = rsd(5,i,j,k)
 145      >          + tx3 * c3 * c4 * ( flux(5,i+1,j,k) - flux(5,i,j,k) )
 146      >              + dx5 * tx1 * (            u(5,i-1,j,k)
 147      >                             - 2.0d+00 * u(5,i,j,k)
 148      >                             +           u(5,i+1,j,k) )
 149             end do
 150
 151 c---------------------------------------------------------------------
 152 c   Fourth-order dissipation
 153 c---------------------------------------------------------------------
 154             IF (north.eq.-1) then
 155              do m = 1, 5
 156                rsd(m,2,j,k) = rsd(m,2,j,k)
 157      >           - dssp * ( + 5.0d+00 * u(m,2,j,k)
 158      >                      - 4.0d+00 * u(m,3,j,k)
 159      >                      +           u(m,4,j,k) )
 160                rsd(m,3,j,k) = rsd(m,3,j,k)
 161      >           - dssp * ( - 4.0d+00 * u(m,2,j,k)
 162      >                      + 6.0d+00 * u(m,3,j,k)
 163      >                      - 4.0d+00 * u(m,4,j,k)
 164      >                      +           u(m,5,j,k) )
 165              end do
 166             END IF
 167
 168             do i = ist1,iend1
 169                do m = 1, 5
 170                   rsd(m,i,j,k) = rsd(m,i,j,k)
 171      >              - dssp * (            u(m,i-2,j,k)
 172      >                        - 4.0d+00 * u(m,i-1,j,k)
 173      >                        + 6.0d+00 * u(m,i,j,k)
 174      >                        - 4.0d+00 * u(m,i+1,j,k)
 175      >                        +           u(m,i+2,j,k) )
 176                end do
 177             end do
 178
 179             IF (south.eq.-1) then
 180              do m = 1, 5
 181                rsd(m,nx-2,j,k) = rsd(m,nx-2,j,k)
 182      >           - dssp * (             u(m,nx-4,j,k)
 183      >                      - 4.0d+00 * u(m,nx-3,j,k)
 184      >                      + 6.0d+00 * u(m,nx-2,j,k)
 185      >                      - 4.0d+00 * u(m,nx-1,j,k)  )
 186                rsd(m,nx-1,j,k) = rsd(m,nx-1,j,k)
 187      >           - dssp * (             u(m,nx-3,j,k)
 188      >                      - 4.0d+00 * u(m,nx-2,j,k)
 189      >                      + 5.0d+00 * u(m,nx-1,j,k) )
 190              end do
 191             END IF
 192
 193          end do
 194       end do
 195
 196 c---------------------------------------------------------------------
 197 c   eta-direction flux differences
 198 c---------------------------------------------------------------------
 199
 200 c---------------------------------------------------------------------
 201 c   iex = flag : iex = 0  north/south communication
 202 c---------------------------------------------------------------------
 203       iex   = 1
 204
 205 c---------------------------------------------------------------------
 206 c   communicate and receive/send two rows of data
 207 c---------------------------------------------------------------------
 208       call exchange_3(u,iex)
 209
 210       L1 = 0
 211       if (west.eq.-1) L1 = 1
 212       L2 = ny + 1
 213       if (east.eq.-1) L2 = ny
 214
 215       jst1 = 1
 216       jend1 = ny
 217       if (west.eq.-1) jst1 = 4
 218       if (east.eq.-1) jend1 = ny - 3
 219
 220       do k = 2, nz - 1
 221          do j = L1, L2
 222             do i = ist, iend
 223                flux(1,i,j,k) = u(3,i,j,k)
 224                u31 = u(3,i,j,k) / u(1,i,j,k)
 225
 226                q = 0.50d+00 * (  u(2,i,j,k) * u(2,i,j,k)
 227      >                         + u(3,i,j,k) * u(3,i,j,k)
 228      >                         + u(4,i,j,k) * u(4,i,j,k) )
 229      >                      / u(1,i,j,k)
 230
 231                flux(2,i,j,k) = u(2,i,j,k) * u31
 232                flux(3,i,j,k) = u(3,i,j,k) * u31 + c2 * (u(5,i,j,k)-q)
 233                flux(4,i,j,k) = u(4,i,j,k) * u31
 234                flux(5,i,j,k) = ( c1 * u(5,i,j,k) - c2 * q ) * u31
 235             end do
 236          end do
 237
 238          do j = jst, jend
 239             do i = ist, iend
 240                do m = 1, 5
 241                   rsd(m,i,j,k) =  rsd(m,i,j,k)
 242      >                   - ty2 * ( flux(m,i,j+1,k) - flux(m,i,j-1,k) )
 243                end do
 244             end do
 245          end do
 246
 247          do j = jst, L2
 248             do i = ist, iend
 249                tmp = 1.0d+00 / u(1,i,j,k)
 250
 251                u21j = tmp * u(2,i,j,k)
 252                u31j = tmp * u(3,i,j,k)
 253                u41j = tmp * u(4,i,j,k)
 254                u51j = tmp * u(5,i,j,k)
 255
 256                tmp = 1.0d+00 / u(1,i,j-1,k)
 257                u21jm1 = tmp * u(2,i,j-1,k)
 258                u31jm1 = tmp * u(3,i,j-1,k)
 259                u41jm1 = tmp * u(4,i,j-1,k)
 260                u51jm1 = tmp * u(5,i,j-1,k)
 261
 262                flux(2,i,j,k) = ty3 * ( u21j - u21jm1 )
 263                flux(3,i,j,k) = (4.0d+00/3.0d+00) * ty3 * (u31j-u31jm1)
 264                flux(4,i,j,k) = ty3 * ( u41j - u41jm1 )
 265                flux(5,i,j,k) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
 266      >              * ty3 * ( ( u21j  **2 + u31j  **2 + u41j  **2 )
 267      >                      - ( u21jm1**2 + u31jm1**2 + u41jm1**2 ) )
 268      >              + (1.0d+00/6.0d+00)
 269      >              * ty3 * ( u31j**2 - u31jm1**2 )
 270      >              + c1 * c5 * ty3 * ( u51j - u51jm1 )
 271             end do
 272          end do
 273
 274          do j = jst, jend
 275             do i = ist, iend
 276
 277                rsd(1,i,j,k) = rsd(1,i,j,k)
 278      >              + dy1 * ty1 * (            u(1,i,j-1,k)
 279      >                             - 2.0d+00 * u(1,i,j,k)
 280      >                             +           u(1,i,j+1,k) )
 281
 282                rsd(2,i,j,k) = rsd(2,i,j,k)
 283      >          + ty3 * c3 * c4 * ( flux(2,i,j+1,k) - flux(2,i,j,k) )
 284      >              + dy2 * ty1 * (            u(2,i,j-1,k)
 285      >                             - 2.0d+00 * u(2,i,j,k)
 286      >                             +           u(2,i,j+1,k) )
 287
 288                rsd(3,i,j,k) = rsd(3,i,j,k)
 289      >          + ty3 * c3 * c4 * ( flux(3,i,j+1,k) - flux(3,i,j,k) )
 290      >              + dy3 * ty1 * (            u(3,i,j-1,k)
 291      >                             - 2.0d+00 * u(3,i,j,k)
 292      >                             +           u(3,i,j+1,k) )
 293
 294                rsd(4,i,j,k) = rsd(4,i,j,k)
 295      >          + ty3 * c3 * c4 * ( flux(4,i,j+1,k) - flux(4,i,j,k) )
 296      >              + dy4 * ty1 * (            u(4,i,j-1,k)
 297      >                             - 2.0d+00 * u(4,i,j,k)
 298      >                             +           u(4,i,j+1,k) )
 299
 300                rsd(5,i,j,k) = rsd(5,i,j,k)
 301      >          + ty3 * c3 * c4 * ( flux(5,i,j+1,k) - flux(5,i,j,k) )
 302      >              + dy5 * ty1 * (            u(5,i,j-1,k)
 303      >                             - 2.0d+00 * u(5,i,j,k)
 304      >                             +           u(5,i,j+1,k) )
 305
 306             end do
 307          end do
 308
 309 c---------------------------------------------------------------------
 310 c   fourth-order dissipation
 311 c---------------------------------------------------------------------
 312          IF (west.eq.-1) then
 313             do i = ist, iend
 314              do m = 1, 5
 315                rsd(m,i,2,k) = rsd(m,i,2,k)
 316      >           - dssp * ( + 5.0d+00 * u(m,i,2,k)
 317      >                      - 4.0d+00 * u(m,i,3,k)
 318      >                      +           u(m,i,4,k) )
 319                rsd(m,i,3,k) = rsd(m,i,3,k)
 320      >           - dssp * ( - 4.0d+00 * u(m,i,2,k)
 321      >                      + 6.0d+00 * u(m,i,3,k)
 322      >                      - 4.0d+00 * u(m,i,4,k)
 323      >                      +           u(m,i,5,k) )
 324              end do
 325             end do
 326          END IF
 327
 328          do j = jst1, jend1
 329             do i = ist, iend
 330                do m = 1, 5
 331                   rsd(m,i,j,k) = rsd(m,i,j,k)
 332      >              - dssp * (            u(m,i,j-2,k)
 333      >                        - 4.0d+00 * u(m,i,j-1,k)
 334      >                        + 6.0d+00 * u(m,i,j,k)
 335      >                        - 4.0d+00 * u(m,i,j+1,k)
 336      >                        +           u(m,i,j+2,k) )
 337                end do
 338             end do
 339          end do
 340
 341          IF (east.eq.-1) then
 342             do i = ist, iend
 343              do m = 1, 5
 344                rsd(m,i,ny-2,k) = rsd(m,i,ny-2,k)
 345      >           - dssp * (             u(m,i,ny-4,k)
 346      >                      - 4.0d+00 * u(m,i,ny-3,k)
 347      >                      + 6.0d+00 * u(m,i,ny-2,k)
 348      >                      - 4.0d+00 * u(m,i,ny-1,k)  )
 349                rsd(m,i,ny-1,k) = rsd(m,i,ny-1,k)
 350      >           - dssp * (             u(m,i,ny-3,k)
 351      >                      - 4.0d+00 * u(m,i,ny-2,k)
 352      >                      + 5.0d+00 * u(m,i,ny-1,k) )
 353              end do
 354             end do
 355          END IF
 356
 357       end do
 358
 359 c---------------------------------------------------------------------
 360 c   zeta-direction flux differences
 361 c---------------------------------------------------------------------
 362       do k = 1, nz
 363          do j = jst, jend
 364             do i = ist, iend
 365                flux(1,i,j,k) = u(4,i,j,k)
 366                u41 = u(4,i,j,k) / u(1,i,j,k)
 367
 368                q = 0.50d+00 * (  u(2,i,j,k) * u(2,i,j,k)
 369      >                         + u(3,i,j,k) * u(3,i,j,k)
 370      >                         + u(4,i,j,k) * u(4,i,j,k) )
 371      >                      / u(1,i,j,k)
 372
 373                flux(2,i,j,k) = u(2,i,j,k) * u41
 374                flux(3,i,j,k) = u(3,i,j,k) * u41
 375                flux(4,i,j,k) = u(4,i,j,k) * u41 + c2 * (u(5,i,j,k)-q)
 376                flux(5,i,j,k) = ( c1 * u(5,i,j,k) - c2 * q ) * u41
 377             end do
 378          end do
 379       end do
 380
 381       do k = 2, nz - 1
 382          do j = jst, jend
 383             do i = ist, iend
 384                do m = 1, 5
 385                   rsd(m,i,j,k) =  rsd(m,i,j,k)
 386      >                - tz2 * ( flux(m,i,j,k+1) - flux(m,i,j,k-1) )
 387                end do
 388             end do
 389          end do
 390       end do
 391
 392       do k = 2, nz
 393          do j = jst, jend
 394             do i = ist, iend
 395                tmp = 1.0d+00 / u(1,i,j,k)
 396
 397                u21k = tmp * u(2,i,j,k)
 398                u31k = tmp * u(3,i,j,k)
 399                u41k = tmp * u(4,i,j,k)
 400                u51k = tmp * u(5,i,j,k)
 401
 402                tmp = 1.0d+00 / u(1,i,j,k-1)
 403
 404                u21km1 = tmp * u(2,i,j,k-1)
 405                u31km1 = tmp * u(3,i,j,k-1)
 406                u41km1 = tmp * u(4,i,j,k-1)
 407                u51km1 = tmp * u(5,i,j,k-1)
 408
 409                flux(2,i,j,k) = tz3 * ( u21k - u21km1 )
 410                flux(3,i,j,k) = tz3 * ( u31k - u31km1 )
 411                flux(4,i,j,k) = (4.0d+00/3.0d+00) * tz3 * (u41k-u41km1)
 412                flux(5,i,j,k) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
 413      >              * tz3 * ( ( u21k  **2 + u31k  **2 + u41k  **2 )
 414      >                      - ( u21km1**2 + u31km1**2 + u41km1**2 ) )
 415      >              + (1.0d+00/6.0d+00)
 416      >              * tz3 * ( u41k**2 - u41km1**2 )
 417      >              + c1 * c5 * tz3 * ( u51k - u51km1 )
 418             end do
 419          end do
 420       end do
 421
 422       do k = 2, nz - 1
 423          do j = jst, jend
 424             do i = ist, iend
 425                rsd(1,i,j,k) = rsd(1,i,j,k)
 426      >              + dz1 * tz1 * (            u(1,i,j,k-1)
 427      >                             - 2.0d+00 * u(1,i,j,k)
 428      >                             +           u(1,i,j,k+1) )
 429                rsd(2,i,j,k) = rsd(2,i,j,k)
 430      >          + tz3 * c3 * c4 * ( flux(2,i,j,k+1) - flux(2,i,j,k) )
 431      >              + dz2 * tz1 * (            u(2,i,j,k-1)
 432      >                             - 2.0d+00 * u(2,i,j,k)
 433      >                             +           u(2,i,j,k+1) )
 434                rsd(3,i,j,k) = rsd(3,i,j,k)
 435      >          + tz3 * c3 * c4 * ( flux(3,i,j,k+1) - flux(3,i,j,k) )
 436      >              + dz3 * tz1 * (            u(3,i,j,k-1)
 437      >                             - 2.0d+00 * u(3,i,j,k)
 438      >                             +           u(3,i,j,k+1) )
 439                rsd(4,i,j,k) = rsd(4,i,j,k)
 440      >          + tz3 * c3 * c4 * ( flux(4,i,j,k+1) - flux(4,i,j,k) )
 441      >              + dz4 * tz1 * (            u(4,i,j,k-1)
 442      >                             - 2.0d+00 * u(4,i,j,k)
 443      >                             +           u(4,i,j,k+1) )
 444                rsd(5,i,j,k) = rsd(5,i,j,k)
 445      >          + tz3 * c3 * c4 * ( flux(5,i,j,k+1) - flux(5,i,j,k) )
 446      >              + dz5 * tz1 * (            u(5,i,j,k-1)
 447      >                             - 2.0d+00 * u(5,i,j,k)
 448      >                             +           u(5,i,j,k+1) )
 449             end do
 450          end do
 451       end do
 452
 453 c---------------------------------------------------------------------
 454 c   fourth-order dissipation
 455 c---------------------------------------------------------------------
 456       do j = jst, jend
 457          do i = ist, iend
 458             do m = 1, 5
 459                rsd(m,i,j,2) = rsd(m,i,j,2)
 460      >           - dssp * ( + 5.0d+00 * u(m,i,j,2)
 461      >                      - 4.0d+00 * u(m,i,j,3)
 462      >                      +           u(m,i,j,4) )
 463                rsd(m,i,j,3) = rsd(m,i,j,3)
 464      >           - dssp * ( - 4.0d+00 * u(m,i,j,2)
 465      >                      + 6.0d+00 * u(m,i,j,3)
 466      >                      - 4.0d+00 * u(m,i,j,4)
 467      >                      +           u(m,i,j,5) )
 468             end do
 469          end do
 470       end do
 471
 472       do k = 4, nz - 3
 473          do j = jst, jend
 474             do i = ist, iend
 475                do m = 1, 5
 476                   rsd(m,i,j,k) = rsd(m,i,j,k)
 477      >              - dssp * (            u(m,i,j,k-2)
 478      >                        - 4.0d+00 * u(m,i,j,k-1)
 479      >                        + 6.0d+00 * u(m,i,j,k)
 480      >                        - 4.0d+00 * u(m,i,j,k+1)
 481      >                        +           u(m,i,j,k+2) )
 482                end do
 483             end do
 484          end do
 485       end do
 486
 487       do j = jst, jend
 488          do i = ist, iend
 489             do m = 1, 5
 490                rsd(m,i,j,nz-2) = rsd(m,i,j,nz-2)
 491      >           - dssp * (             u(m,i,j,nz-4)
 492      >                      - 4.0d+00 * u(m,i,j,nz-3)
 493      >                      + 6.0d+00 * u(m,i,j,nz-2)
 494      >                      - 4.0d+00 * u(m,i,j,nz-1)  )
 495                rsd(m,i,j,nz-1) = rsd(m,i,j,nz-1)
 496      >           - dssp * (             u(m,i,j,nz-3)
 497      >                      - 4.0d+00 * u(m,i,j,nz-2)
 498      >                      + 5.0d+00 * u(m,i,j,nz-1) )
 499             end do
 500          end do
 501       end do
 502
 503       return
 504       end