examples/smpi/NAS/SP/lhsz.f

   1
   2 c---------------------------------------------------------------------
   3 c---------------------------------------------------------------------
   4
   5        subroutine lhsz(c)
   6
   7 c---------------------------------------------------------------------
   8 c---------------------------------------------------------------------
   9
  10 c---------------------------------------------------------------------
  11 c This function computes the left hand side for the three z-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   j = start(2,c), cell_size(2,c)-end(2,c)-1
  27           do   i = start(1,c), cell_size(1,c)-end(1,c)-1
  28
  29              do   k = start(3,c)-1, cell_size(3,c)-end(3,c)
  30                 ru1 = c3c4*rho_i(i,j,k,c)
  31                 cv(k) = ws(i,j,k,c)
  32                 rhos(k) = dmax1(dz4 + con43 * ru1,
  33      >                          dz5 + c1c5 * ru1,
  34      >                          dzmax + ru1,
  35      >                          dz1)
  36              end do
  37
  38              do   k =  start(3,c), cell_size(3,c)-end(3,c)-1
  39                 lhs(i,j,k,1,c) =  0.0d0
  40                 lhs(i,j,k,2,c) = -dttz2 * cv(k-1) - dttz1 * rhos(k-1)
  41                 lhs(i,j,k,3,c) =  1.0 + c2dttz1 * rhos(k)
  42                 lhs(i,j,k,4,c) =  dttz2 * cv(k+1) - dttz1 * rhos(k+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(3,c) .gt. 0) then
  52           k = 1
  53           do    j = start(2,c), cell_size(2,c)-end(2,c)-1
  54              do    i = start(1,c), cell_size(1,c)-end(1,c)-1
  55                 lhs(i,j,k,3,c) = lhs(i,j,k,3,c) + comz5
  56                 lhs(i,j,k,4,c) = lhs(i,j,k,4,c) - comz4
  57                 lhs(i,j,k,5,c) = lhs(i,j,k,5,c) + comz1
  58
  59                 lhs(i,j,k+1,2,c) = lhs(i,j,k+1,2,c) - comz4
  60                 lhs(i,j,k+1,3,c) = lhs(i,j,k+1,3,c) + comz6
  61                 lhs(i,j,k+1,4,c) = lhs(i,j,k+1,4,c) - comz4
  62                 lhs(i,j,k+1,5,c) = lhs(i,j,k+1,5,c) + comz1
  63              end do
  64           end do
  65        endif
  66
  67        do    k = 3*start(3,c), cell_size(3,c)-3*end(3,c)-1
  68           do    j = start(2,c), cell_size(2,c)-end(2,c)-1
  69              do    i = start(1,c), cell_size(1,c)-end(1,c)-1
  70                 lhs(i,j,k,1,c) = lhs(i,j,k,1,c) + comz1
  71                 lhs(i,j,k,2,c) = lhs(i,j,k,2,c) - comz4
  72                 lhs(i,j,k,3,c) = lhs(i,j,k,3,c) + comz6
  73                 lhs(i,j,k,4,c) = lhs(i,j,k,4,c) - comz4
  74                 lhs(i,j,k,5,c) = lhs(i,j,k,5,c) + comz1
  75              end do
  76           end do
  77        end do
  78
  79        if (end(3,c) .gt. 0) then
  80           k = cell_size(3,c)-3
  81           do    j = start(2,c), cell_size(2,c)-end(2,c)-1
  82              do    i = start(1,c), cell_size(1,c)-end(1,c)-1
  83                 lhs(i,j,k,1,c) = lhs(i,j,k,1,c) + comz1
  84                 lhs(i,j,k,2,c) = lhs(i,j,k,2,c) - comz4
  85                 lhs(i,j,k,3,c) = lhs(i,j,k,3,c) + comz6
  86                 lhs(i,j,k,4,c) = lhs(i,j,k,4,c) - comz4
  87
  88                 lhs(i,j,k+1,1,c) = lhs(i,j,k+1,1,c) + comz1
  89                 lhs(i,j,k+1,2,c) = lhs(i,j,k+1,2,c) - comz4
  90                 lhs(i,j,k+1,3,c) = lhs(i,j,k+1,3,c) + comz5
  91              end do
  92           end do
  93        endif
  94
  95
  96 c---------------------------------------------------------------------
  97 c      subsequently, fill the other factors (u+c), (u-c)
  98 c---------------------------------------------------------------------
  99        do    k = start(3,c), cell_size(3,c)-end(3,c)-1
 100           do    j = start(2,c), cell_size(2,c)-end(2,c)-1
 101              do    i = start(1,c), cell_size(1,c)-end(1,c)-1
 102                 lhs(i,j,k,1+5,c)  = lhs(i,j,k,1,c)
 103                 lhs(i,j,k,2+5,c)  = lhs(i,j,k,2,c) -
 104      >                            dttz2 * speed(i,j,k-1,c)
 105                 lhs(i,j,k,3+5,c)  = lhs(i,j,k,3,c)
 106                 lhs(i,j,k,4+5,c)  = lhs(i,j,k,4,c) +
 107      >                            dttz2 * speed(i,j,k+1,c)
 108                 lhs(i,j,k,5+5,c) = lhs(i,j,k,5,c)
 109                 lhs(i,j,k,1+10,c) = lhs(i,j,k,1,c)
 110                 lhs(i,j,k,2+10,c) = lhs(i,j,k,2,c) +
 111      >                            dttz2 * speed(i,j,k-1,c)
 112                 lhs(i,j,k,3+10,c) = lhs(i,j,k,3,c)
 113                 lhs(i,j,k,4+10,c) = lhs(i,j,k,4,c) -
 114      >                            dttz2 * speed(i,j,k+1,c)
 115                 lhs(i,j,k,5+10,c) = lhs(i,j,k,5,c)
 116              end do
 117           end do
 118        end do
 119
 120        return
 121        end
 122
 123