examples/smpi/NAS/SP/lhsy.f

   1
   2 c---------------------------------------------------------------------
   3 c---------------------------------------------------------------------
   4
   5        subroutine lhsy(c)
   6
   7 c---------------------------------------------------------------------
   8 c---------------------------------------------------------------------
   9
  10 c---------------------------------------------------------------------
  11 c This function computes the left hand side for the three y-factors
  12 c---------------------------------------------------------------------
  13
  14        include 'header.h'
  15
  16        double precision ru1
  17        integer          i, j, k, c
  18
  19 c---------------------------------------------------------------------
  20 c      treat only cell c
  21 c---------------------------------------------------------------------
  22
  23 c---------------------------------------------------------------------
  24 c      first fill the lhs for the u-eigenvalue
  25 c---------------------------------------------------------------------
  26        do  k = start(3,c), cell_size(3,c)-end(3,c)-1
  27           do  i = start(1,c), cell_size(1,c)-end(1,c)-1
  28
  29              do  j = start(2,c)-1, cell_size(2,c)-end(2,c)
  30                 ru1 = c3c4*rho_i(i,j,k,c)
  31                 cv(j) = vs(i,j,k,c)
  32                 rhoq(j) = dmax1( dy3 + con43 * ru1,
  33      >                           dy5 + c1c5*ru1,
  34      >                           dymax + ru1,
  35      >                           dy1)
  36              end do
  37
  38              do  j = start(2,c), cell_size(2,c)-end(2,c)-1
  39                 lhs(i,j,k,1,c) =  0.0d0
  40                 lhs(i,j,k,2,c) = -dtty2 * cv(j-1) - dtty1 * rhoq(j-1)
  41                 lhs(i,j,k,3,c) =  1.0 + c2dtty1 * rhoq(j)
  42                 lhs(i,j,k,4,c) =  dtty2 * cv(j+1) - dtty1 * rhoq(j+1)
  43                 lhs(i,j,k,5,c) =  0.0d0
  44              end do
  45           end do
  46        end do
  47
  48 c---------------------------------------------------------------------
  49 c      add fourth order dissipation
  50 c---------------------------------------------------------------------
  51        if (start(2,c) .gt. 0) then
  52           j = 1
  53           do   k = start(3,c), cell_size(3,c)-end(3,c)-1
  54              do   i = start(1,c), cell_size(1,c)-end(1,c)-1
  55
  56                 lhs(i,j,k,3,c) = lhs(i,j,k,3,c) + comz5
  57                 lhs(i,j,k,4,c) = lhs(i,j,k,4,c) - comz4
  58                 lhs(i,j,k,5,c) = lhs(i,j,k,5,c) + comz1
  59
  60                 lhs(i,j+1,k,2,c) = lhs(i,j+1,k,2,c) - comz4
  61                 lhs(i,j+1,k,3,c) = lhs(i,j+1,k,3,c) + comz6
  62                 lhs(i,j+1,k,4,c) = lhs(i,j+1,k,4,c) - comz4
  63                 lhs(i,j+1,k,5,c) = lhs(i,j+1,k,5,c) + comz1
  64              end do
  65           end do
  66        endif
  67
  68        do   k = start(3,c), cell_size(3,c)-end(3,c)-1
  69           do   j=3*start(2,c), cell_size(2,c)-3*end(2,c)-1
  70              do   i = start(1,c), cell_size(1,c)-end(1,c)-1
  71
  72                 lhs(i,j,k,1,c) = lhs(i,j,k,1,c) + comz1
  73                 lhs(i,j,k,2,c) = lhs(i,j,k,2,c) - comz4
  74                 lhs(i,j,k,3,c) = lhs(i,j,k,3,c) + comz6
  75                 lhs(i,j,k,4,c) = lhs(i,j,k,4,c) - comz4
  76                 lhs(i,j,k,5,c) = lhs(i,j,k,5,c) + comz1
  77              end do
  78           end do
  79        end do
  80
  81        if (end(2,c) .gt. 0) then
  82           j = cell_size(2,c)-3
  83           do   k = start(3,c), cell_size(3,c)-end(3,c)-1
  84              do   i = start(1,c), cell_size(1,c)-end(1,c)-1
  85                 lhs(i,j,k,1,c) = lhs(i,j,k,1,c) + comz1
  86                 lhs(i,j,k,2,c) = lhs(i,j,k,2,c) - comz4
  87                 lhs(i,j,k,3,c) = lhs(i,j,k,3,c) + comz6
  88                 lhs(i,j,k,4,c) = lhs(i,j,k,4,c) - comz4
  89
  90                 lhs(i,j+1,k,1,c) = lhs(i,j+1,k,1,c) + comz1
  91                 lhs(i,j+1,k,2,c) = lhs(i,j+1,k,2,c) - comz4
  92                 lhs(i,j+1,k,3,c) = lhs(i,j+1,k,3,c) + comz5
  93              end do
  94           end do
  95        endif
  96
  97 c---------------------------------------------------------------------
  98 c      subsequently, do the other two factors
  99 c---------------------------------------------------------------------
 100        do    k = start(3,c), cell_size(3,c)-end(3,c)-1
 101           do    j = start(2,c), cell_size(2,c)-end(2,c)-1
 102              do    i = start(1,c), cell_size(1,c)-end(1,c)-1
 103                 lhs(i,j,k,1+5,c)  = lhs(i,j,k,1,c)
 104                 lhs(i,j,k,2+5,c)  = lhs(i,j,k,2,c) -
 105      >                            dtty2 * speed(i,j-1,k,c)
 106                 lhs(i,j,k,3+5,c)  = lhs(i,j,k,3,c)
 107                 lhs(i,j,k,4+5,c)  = lhs(i,j,k,4,c) +
 108      >                            dtty2 * speed(i,j+1,k,c)
 109                 lhs(i,j,k,5+5,c) = lhs(i,j,k,5,c)
 110                 lhs(i,j,k,1+10,c) = lhs(i,j,k,1,c)
 111                 lhs(i,j,k,2+10,c) = lhs(i,j,k,2,c) +
 112      >                            dtty2 * speed(i,j-1,k,c)
 113                 lhs(i,j,k,3+10,c) = lhs(i,j,k,3,c)
 114                 lhs(i,j,k,4+10,c) = lhs(i,j,k,4,c) -
 115      >                            dtty2 * speed(i,j+1,k,c)
 116                 lhs(i,j,k,5+10,c) = lhs(i,j,k,5,c)
 117              end do
 118           end do
 119        end do
 120
 121        return
 122        end
 123
 124
 125