examples/smpi/NAS/LU/jacld.f

   1
   2 c---------------------------------------------------------------------
   3 c---------------------------------------------------------------------
   4
   5       subroutine jacld(k)
   6
   7 c---------------------------------------------------------------------
   8 c---------------------------------------------------------------------
   9
  10
  11 c---------------------------------------------------------------------
  12 c   compute the lower triangular part of the jacobian matrix
  13 c---------------------------------------------------------------------
  14
  15       implicit none
  16
  17       include 'applu.incl'
  18
  19 c---------------------------------------------------------------------
  20 c  input parameters
  21 c---------------------------------------------------------------------
  22       integer k
  23
  24 c---------------------------------------------------------------------
  25 c  local variables
  26 c---------------------------------------------------------------------
  27       integer i, j
  28       double precision  r43
  29       double precision  c1345
  30       double precision  c34
  31       double precision  tmp1, tmp2, tmp3
  32
  33
  34
  35       r43 = ( 4.0d+00 / 3.0d+00 )
  36       c1345 = c1 * c3 * c4 * c5
  37       c34 = c3 * c4
  38
  39          do j = jst, jend
  40             do i = ist, iend
  41
  42 c---------------------------------------------------------------------
  43 c   form the block daigonal
  44 c---------------------------------------------------------------------
  45                tmp1 = 1.0d+00 / u(1,i,j,k)
  46                tmp2 = tmp1 * tmp1
  47                tmp3 = tmp1 * tmp2
  48
  49                d(1,1,i,j) =  1.0d+00
  50      >                       + dt * 2.0d+00 * (   tx1 * dx1
  51      >                                          + ty1 * dy1
  52      >                                          + tz1 * dz1 )
  53                d(1,2,i,j) =  0.0d+00
  54                d(1,3,i,j) =  0.0d+00
  55                d(1,4,i,j) =  0.0d+00
  56                d(1,5,i,j) =  0.0d+00
  57
  58                d(2,1,i,j) =  dt * 2.0d+00
  59      >          * (  tx1 * ( - r43 * c34 * tmp2 * u(2,i,j,k) )
  60      >             + ty1 * ( -       c34 * tmp2 * u(2,i,j,k) )
  61      >             + tz1 * ( -       c34 * tmp2 * u(2,i,j,k) ) )
  62                d(2,2,i,j) =  1.0d+00
  63      >          + dt * 2.0d+00
  64      >          * (  tx1 * r43 * c34 * tmp1
  65      >             + ty1 *       c34 * tmp1
  66      >             + tz1 *       c34 * tmp1 )
  67      >          + dt * 2.0d+00 * (   tx1 * dx2
  68      >                             + ty1 * dy2
  69      >                             + tz1 * dz2  )
  70                d(2,3,i,j) = 0.0d+00
  71                d(2,4,i,j) = 0.0d+00
  72                d(2,5,i,j) = 0.0d+00
  73
  74                d(3,1,i,j) = dt * 2.0d+00
  75      >      * (  tx1 * ( -       c34 * tmp2 * u(3,i,j,k) )
  76      >         + ty1 * ( - r43 * c34 * tmp2 * u(3,i,j,k) )
  77      >         + tz1 * ( -       c34 * tmp2 * u(3,i,j,k) ) )
  78                d(3,2,i,j) = 0.0d+00
  79                d(3,3,i,j) = 1.0d+00
  80      >         + dt * 2.0d+00
  81      >              * (  tx1 *       c34 * tmp1
  82      >                 + ty1 * r43 * c34 * tmp1
  83      >                 + tz1 *       c34 * tmp1 )
  84      >         + dt * 2.0d+00 * (  tx1 * dx3
  85      >                           + ty1 * dy3
  86      >                           + tz1 * dz3 )
  87                d(3,4,i,j) = 0.0d+00
  88                d(3,5,i,j) = 0.0d+00
  89
  90                d(4,1,i,j) = dt * 2.0d+00
  91      >      * (  tx1 * ( -       c34 * tmp2 * u(4,i,j,k) )
  92      >         + ty1 * ( -       c34 * tmp2 * u(4,i,j,k) )
  93      >         + tz1 * ( - r43 * c34 * tmp2 * u(4,i,j,k) ) )
  94                d(4,2,i,j) = 0.0d+00
  95                d(4,3,i,j) = 0.0d+00
  96                d(4,4,i,j) = 1.0d+00
  97      >         + dt * 2.0d+00
  98      >              * (  tx1 *       c34 * tmp1
  99      >                 + ty1 *       c34 * tmp1
 100      >                 + tz1 * r43 * c34 * tmp1 )
 101      >         + dt * 2.0d+00 * (  tx1 * dx4
 102      >                           + ty1 * dy4
 103      >                           + tz1 * dz4 )
 104                d(4,5,i,j) = 0.0d+00
 105
 106                d(5,1,i,j) = dt * 2.0d+00
 107      > * ( tx1 * ( - ( r43*c34 - c1345 ) * tmp3 * ( u(2,i,j,k) ** 2 )
 108      >             - ( c34 - c1345 ) * tmp3 * ( u(3,i,j,k) ** 2 )
 109      >             - ( c34 - c1345 ) * tmp3 * ( u(4,i,j,k) ** 2 )
 110      >             - ( c1345 ) * tmp2 * u(5,i,j,k) )
 111      >   + ty1 * ( - ( c34 - c1345 ) * tmp3 * ( u(2,i,j,k) ** 2 )
 112      >             - ( r43*c34 - c1345 ) * tmp3 * ( u(3,i,j,k) ** 2 )
 113      >             - ( c34 - c1345 ) * tmp3 * ( u(4,i,j,k) ** 2 )
 114      >             - ( c1345 ) * tmp2 * u(5,i,j,k) )
 115      >   + tz1 * ( - ( c34 - c1345 ) * tmp3 * ( u(2,i,j,k) ** 2 )
 116      >             - ( c34 - c1345 ) * tmp3 * ( u(3,i,j,k) ** 2 )
 117      >             - ( r43*c34 - c1345 ) * tmp3 * ( u(4,i,j,k) ** 2 )
 118      >             - ( c1345 ) * tmp2 * u(5,i,j,k) ) )
 119                d(5,2,i,j) = dt * 2.0d+00
 120      > * ( tx1 * ( r43*c34 - c1345 ) * tmp2 * u(2,i,j,k)
 121      >   + ty1 * (     c34 - c1345 ) * tmp2 * u(2,i,j,k)
 122      >   + tz1 * (     c34 - c1345 ) * tmp2 * u(2,i,j,k) )
 123                d(5,3,i,j) = dt * 2.0d+00
 124      > * ( tx1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k)
 125      >   + ty1 * ( r43*c34 -c1345 ) * tmp2 * u(3,i,j,k)
 126      >   + tz1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k) )
 127                d(5,4,i,j) = dt * 2.0d+00
 128      > * ( tx1 * ( c34 - c1345 ) * tmp2 * u(4,i,j,k)
 129      >   + ty1 * ( c34 - c1345 ) * tmp2 * u(4,i,j,k)
 130      >   + tz1 * ( r43*c34 - c1345 ) * tmp2 * u(4,i,j,k) )
 131                d(5,5,i,j) = 1.0d+00
 132      >   + dt * 2.0d+00 * ( tx1 * c1345 * tmp1
 133      >                    + ty1 * c1345 * tmp1
 134      >                    + tz1 * c1345 * tmp1 )
 135      >   + dt * 2.0d+00 * (  tx1 * dx5
 136      >                    +  ty1 * dy5
 137      >                    +  tz1 * dz5 )
 138
 139 c---------------------------------------------------------------------
 140 c   form the first block sub-diagonal
 141 c---------------------------------------------------------------------
 142                tmp1 = 1.0d+00 / u(1,i,j,k-1)
 143                tmp2 = tmp1 * tmp1
 144                tmp3 = tmp1 * tmp2
 145
 146                a(1,1,i,j) = - dt * tz1 * dz1
 147                a(1,2,i,j) =   0.0d+00
 148                a(1,3,i,j) =   0.0d+00
 149                a(1,4,i,j) = - dt * tz2
 150                a(1,5,i,j) =   0.0d+00
 151
 152                a(2,1,i,j) = - dt * tz2
 153      >           * ( - ( u(2,i,j,k-1)*u(4,i,j,k-1) ) * tmp2 )
 154      >           - dt * tz1 * ( - c34 * tmp2 * u(2,i,j,k-1) )
 155                a(2,2,i,j) = - dt * tz2 * ( u(4,i,j,k-1) * tmp1 )
 156      >           - dt * tz1 * c34 * tmp1
 157      >           - dt * tz1 * dz2
 158                a(2,3,i,j) = 0.0d+00
 159                a(2,4,i,j) = - dt * tz2 * ( u(2,i,j,k-1) * tmp1 )
 160                a(2,5,i,j) = 0.0d+00
 161
 162                a(3,1,i,j) = - dt * tz2
 163      >           * ( - ( u(3,i,j,k-1)*u(4,i,j,k-1) ) * tmp2 )
 164      >           - dt * tz1 * ( - c34 * tmp2 * u(3,i,j,k-1) )
 165                a(3,2,i,j) = 0.0d+00
 166                a(3,3,i,j) = - dt * tz2 * ( u(4,i,j,k-1) * tmp1 )
 167      >           - dt * tz1 * ( c34 * tmp1 )
 168      >           - dt * tz1 * dz3
 169                a(3,4,i,j) = - dt * tz2 * ( u(3,i,j,k-1) * tmp1 )
 170                a(3,5,i,j) = 0.0d+00
 171
 172                a(4,1,i,j) = - dt * tz2
 173      >        * ( - ( u(4,i,j,k-1) * tmp1 ) ** 2
 174      >            + 0.50d+00 * c2
 175      >            * ( ( u(2,i,j,k-1) * u(2,i,j,k-1)
 176      >                + u(3,i,j,k-1) * u(3,i,j,k-1)
 177      >                + u(4,i,j,k-1) * u(4,i,j,k-1) ) * tmp2 ) )
 178      >        - dt * tz1 * ( - r43 * c34 * tmp2 * u(4,i,j,k-1) )
 179                a(4,2,i,j) = - dt * tz2
 180      >             * ( - c2 * ( u(2,i,j,k-1) * tmp1 ) )
 181                a(4,3,i,j) = - dt * tz2
 182      >             * ( - c2 * ( u(3,i,j,k-1) * tmp1 ) )
 183                a(4,4,i,j) = - dt * tz2 * ( 2.0d+00 - c2 )
 184      >             * ( u(4,i,j,k-1) * tmp1 )
 185      >             - dt * tz1 * ( r43 * c34 * tmp1 )
 186      >             - dt * tz1 * dz4
 187                a(4,5,i,j) = - dt * tz2 * c2
 188
 189                a(5,1,i,j) = - dt * tz2
 190      >     * ( ( c2 * (  u(2,i,j,k-1) * u(2,i,j,k-1)
 191      >                 + u(3,i,j,k-1) * u(3,i,j,k-1)
 192      >                 + u(4,i,j,k-1) * u(4,i,j,k-1) ) * tmp2
 193      >       - c1 * ( u(5,i,j,k-1) * tmp1 ) )
 194      >            * ( u(4,i,j,k-1) * tmp1 ) )
 195      >       - dt * tz1
 196      >       * ( - ( c34 - c1345 ) * tmp3 * (u(2,i,j,k-1)**2)
 197      >           - ( c34 - c1345 ) * tmp3 * (u(3,i,j,k-1)**2)
 198      >           - ( r43*c34 - c1345 )* tmp3 * (u(4,i,j,k-1)**2)
 199      >          - c1345 * tmp2 * u(5,i,j,k-1) )
 200                a(5,2,i,j) = - dt * tz2
 201      >       * ( - c2 * ( u(2,i,j,k-1)*u(4,i,j,k-1) ) * tmp2 )
 202      >       - dt * tz1 * ( c34 - c1345 ) * tmp2 * u(2,i,j,k-1)
 203                a(5,3,i,j) = - dt * tz2
 204      >       * ( - c2 * ( u(3,i,j,k-1)*u(4,i,j,k-1) ) * tmp2 )
 205      >       - dt * tz1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k-1)
 206                a(5,4,i,j) = - dt * tz2
 207      >       * ( c1 * ( u(5,i,j,k-1) * tmp1 )
 208      >       - 0.50d+00 * c2
 209      >       * ( (  u(2,i,j,k-1)*u(2,i,j,k-1)
 210      >            + u(3,i,j,k-1)*u(3,i,j,k-1)
 211      >            + 3.0d+00*u(4,i,j,k-1)*u(4,i,j,k-1) ) * tmp2 ) )
 212      >       - dt * tz1 * ( r43*c34 - c1345 ) * tmp2 * u(4,i,j,k-1)
 213                a(5,5,i,j) = - dt * tz2
 214      >       * ( c1 * ( u(4,i,j,k-1) * tmp1 ) )
 215      >       - dt * tz1 * c1345 * tmp1
 216      >       - dt * tz1 * dz5
 217
 218 c---------------------------------------------------------------------
 219 c   form the second block sub-diagonal
 220 c---------------------------------------------------------------------
 221                tmp1 = 1.0d+00 / u(1,i,j-1,k)
 222                tmp2 = tmp1 * tmp1
 223                tmp3 = tmp1 * tmp2
 224
 225                b(1,1,i,j) = - dt * ty1 * dy1
 226                b(1,2,i,j) =   0.0d+00
 227                b(1,3,i,j) = - dt * ty2
 228                b(1,4,i,j) =   0.0d+00
 229                b(1,5,i,j) =   0.0d+00
 230
 231                b(2,1,i,j) = - dt * ty2
 232      >           * ( - ( u(2,i,j-1,k)*u(3,i,j-1,k) ) * tmp2 )
 233      >           - dt * ty1 * ( - c34 * tmp2 * u(2,i,j-1,k) )
 234                b(2,2,i,j) = - dt * ty2 * ( u(3,i,j-1,k) * tmp1 )
 235      >          - dt * ty1 * ( c34 * tmp1 )
 236      >          - dt * ty1 * dy2
 237                b(2,3,i,j) = - dt * ty2 * ( u(2,i,j-1,k) * tmp1 )
 238                b(2,4,i,j) = 0.0d+00
 239                b(2,5,i,j) = 0.0d+00
 240
 241                b(3,1,i,j) = - dt * ty2
 242      >           * ( - ( u(3,i,j-1,k) * tmp1 ) ** 2
 243      >      + 0.50d+00 * c2 * ( (  u(2,i,j-1,k) * u(2,i,j-1,k)
 244      >                           + u(3,i,j-1,k) * u(3,i,j-1,k)
 245      >                           + u(4,i,j-1,k) * u(4,i,j-1,k) )
 246      >                          * tmp2 ) )
 247      >       - dt * ty1 * ( - r43 * c34 * tmp2 * u(3,i,j-1,k) )
 248                b(3,2,i,j) = - dt * ty2
 249      >                   * ( - c2 * ( u(2,i,j-1,k) * tmp1 ) )
 250                b(3,3,i,j) = - dt * ty2 * ( ( 2.0d+00 - c2 )
 251      >                   * ( u(3,i,j-1,k) * tmp1 ) )
 252      >       - dt * ty1 * ( r43 * c34 * tmp1 )
 253      >       - dt * ty1 * dy3
 254                b(3,4,i,j) = - dt * ty2
 255      >                   * ( - c2 * ( u(4,i,j-1,k) * tmp1 ) )
 256                b(3,5,i,j) = - dt * ty2 * c2
 257
 258                b(4,1,i,j) = - dt * ty2
 259      >              * ( - ( u(3,i,j-1,k)*u(4,i,j-1,k) ) * tmp2 )
 260      >       - dt * ty1 * ( - c34 * tmp2 * u(4,i,j-1,k) )
 261                b(4,2,i,j) = 0.0d+00
 262                b(4,3,i,j) = - dt * ty2 * ( u(4,i,j-1,k) * tmp1 )
 263                b(4,4,i,j) = - dt * ty2 * ( u(3,i,j-1,k) * tmp1 )
 264      >                        - dt * ty1 * ( c34 * tmp1 )
 265      >                        - dt * ty1 * dy4
 266                b(4,5,i,j) = 0.0d+00
 267
 268                b(5,1,i,j) = - dt * ty2
 269      >          * ( ( c2 * (  u(2,i,j-1,k) * u(2,i,j-1,k)
 270      >                      + u(3,i,j-1,k) * u(3,i,j-1,k)
 271      >                      + u(4,i,j-1,k) * u(4,i,j-1,k) ) * tmp2
 272      >               - c1 * ( u(5,i,j-1,k) * tmp1 ) )
 273      >          * ( u(3,i,j-1,k) * tmp1 ) )
 274      >          - dt * ty1
 275      >          * ( - (     c34 - c1345 )*tmp3*(u(2,i,j-1,k)**2)
 276      >              - ( r43*c34 - c1345 )*tmp3*(u(3,i,j-1,k)**2)
 277      >              - (     c34 - c1345 )*tmp3*(u(4,i,j-1,k)**2)
 278      >              - c1345*tmp2*u(5,i,j-1,k) )
 279                b(5,2,i,j) = - dt * ty2
 280      >          * ( - c2 * ( u(2,i,j-1,k)*u(3,i,j-1,k) ) * tmp2 )
 281      >          - dt * ty1
 282      >          * ( c34 - c1345 ) * tmp2 * u(2,i,j-1,k)
 283                b(5,3,i,j) = - dt * ty2
 284      >          * ( c1 * ( u(5,i,j-1,k) * tmp1 )
 285      >          - 0.50d+00 * c2
 286      >          * ( (  u(2,i,j-1,k)*u(2,i,j-1,k)
 287      >               + 3.0d+00 * u(3,i,j-1,k)*u(3,i,j-1,k)
 288      >               + u(4,i,j-1,k)*u(4,i,j-1,k) ) * tmp2 ) )
 289      >          - dt * ty1
 290      >          * ( r43*c34 - c1345 ) * tmp2 * u(3,i,j-1,k)
 291                b(5,4,i,j) = - dt * ty2
 292      >          * ( - c2 * ( u(3,i,j-1,k)*u(4,i,j-1,k) ) * tmp2 )
 293      >          - dt * ty1 * ( c34 - c1345 ) * tmp2 * u(4,i,j-1,k)
 294                b(5,5,i,j) = - dt * ty2
 295      >          * ( c1 * ( u(3,i,j-1,k) * tmp1 ) )
 296      >          - dt * ty1 * c1345 * tmp1
 297      >          - dt * ty1 * dy5
 298
 299 c---------------------------------------------------------------------
 300 c   form the third block sub-diagonal
 301 c---------------------------------------------------------------------
 302                tmp1 = 1.0d+00 / u(1,i-1,j,k)
 303                tmp2 = tmp1 * tmp1
 304                tmp3 = tmp1 * tmp2
 305
 306                c(1,1,i,j) = - dt * tx1 * dx1
 307                c(1,2,i,j) = - dt * tx2
 308                c(1,3,i,j) =   0.0d+00
 309                c(1,4,i,j) =   0.0d+00
 310                c(1,5,i,j) =   0.0d+00
 311
 312                c(2,1,i,j) = - dt * tx2
 313      >          * ( - ( u(2,i-1,j,k) * tmp1 ) ** 2
 314      >     + c2 * 0.50d+00 * (  u(2,i-1,j,k) * u(2,i-1,j,k)
 315      >                        + u(3,i-1,j,k) * u(3,i-1,j,k)
 316      >                        + u(4,i-1,j,k) * u(4,i-1,j,k) ) * tmp2 )
 317      >          - dt * tx1 * ( - r43 * c34 * tmp2 * u(2,i-1,j,k) )
 318                c(2,2,i,j) = - dt * tx2
 319      >          * ( ( 2.0d+00 - c2 ) * ( u(2,i-1,j,k) * tmp1 ) )
 320      >          - dt * tx1 * ( r43 * c34 * tmp1 )
 321      >          - dt * tx1 * dx2
 322                c(2,3,i,j) = - dt * tx2
 323      >              * ( - c2 * ( u(3,i-1,j,k) * tmp1 ) )
 324                c(2,4,i,j) = - dt * tx2
 325      >              * ( - c2 * ( u(4,i-1,j,k) * tmp1 ) )
 326                c(2,5,i,j) = - dt * tx2 * c2
 327
 328                c(3,1,i,j) = - dt * tx2
 329      >              * ( - ( u(2,i-1,j,k) * u(3,i-1,j,k) ) * tmp2 )
 330      >         - dt * tx1 * ( - c34 * tmp2 * u(3,i-1,j,k) )
 331                c(3,2,i,j) = - dt * tx2 * ( u(3,i-1,j,k) * tmp1 )
 332                c(3,3,i,j) = - dt * tx2 * ( u(2,i-1,j,k) * tmp1 )
 333      >          - dt * tx1 * ( c34 * tmp1 )
 334      >          - dt * tx1 * dx3
 335                c(3,4,i,j) = 0.0d+00
 336                c(3,5,i,j) = 0.0d+00
 337
 338                c(4,1,i,j) = - dt * tx2
 339      >          * ( - ( u(2,i-1,j,k)*u(4,i-1,j,k) ) * tmp2 )
 340      >          - dt * tx1 * ( - c34 * tmp2 * u(4,i-1,j,k) )
 341                c(4,2,i,j) = - dt * tx2 * ( u(4,i-1,j,k) * tmp1 )
 342                c(4,3,i,j) = 0.0d+00
 343                c(4,4,i,j) = - dt * tx2 * ( u(2,i-1,j,k) * tmp1 )
 344      >          - dt * tx1 * ( c34 * tmp1 )
 345      >          - dt * tx1 * dx4
 346                c(4,5,i,j) = 0.0d+00
 347
 348                c(5,1,i,j) = - dt * tx2
 349      >          * ( ( c2 * (  u(2,i-1,j,k) * u(2,i-1,j,k)
 350      >                      + u(3,i-1,j,k) * u(3,i-1,j,k)
 351      >                      + u(4,i-1,j,k) * u(4,i-1,j,k) ) * tmp2
 352      >              - c1 * ( u(5,i-1,j,k) * tmp1 ) )
 353      >          * ( u(2,i-1,j,k) * tmp1 ) )
 354      >          - dt * tx1
 355      >          * ( - ( r43*c34 - c1345 ) * tmp3 * ( u(2,i-1,j,k)**2 )
 356      >              - (     c34 - c1345 ) * tmp3 * ( u(3,i-1,j,k)**2 )
 357      >              - (     c34 - c1345 ) * tmp3 * ( u(4,i-1,j,k)**2 )
 358      >              - c1345 * tmp2 * u(5,i-1,j,k) )
 359                c(5,2,i,j) = - dt * tx2
 360      >          * ( c1 * ( u(5,i-1,j,k) * tmp1 )
 361      >             - 0.50d+00 * c2
 362      >             * ( (  3.0d+00*u(2,i-1,j,k)*u(2,i-1,j,k)
 363      >                  + u(3,i-1,j,k)*u(3,i-1,j,k)
 364      >                  + u(4,i-1,j,k)*u(4,i-1,j,k) ) * tmp2 ) )
 365      >           - dt * tx1
 366      >           * ( r43*c34 - c1345 ) * tmp2 * u(2,i-1,j,k)
 367                c(5,3,i,j) = - dt * tx2
 368      >           * ( - c2 * ( u(3,i-1,j,k)*u(2,i-1,j,k) ) * tmp2 )
 369      >           - dt * tx1
 370      >           * (  c34 - c1345 ) * tmp2 * u(3,i-1,j,k)
 371                c(5,4,i,j) = - dt * tx2
 372      >           * ( - c2 * ( u(4,i-1,j,k)*u(2,i-1,j,k) ) * tmp2 )
 373      >           - dt * tx1
 374      >           * (  c34 - c1345 ) * tmp2 * u(4,i-1,j,k)
 375                c(5,5,i,j) = - dt * tx2
 376      >           * ( c1 * ( u(2,i-1,j,k) * tmp1 ) )
 377      >           - dt * tx1 * c1345 * tmp1
 378      >           - dt * tx1 * dx5
 379
 380             end do
 381          end do
 382
 383       return
 384       end