examples/smpi/NAS/LU/erhs.f

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