examples/smpi/NAS/SP/rhs.f

   1 c---------------------------------------------------------------------
   2 c---------------------------------------------------------------------
   3
   4        subroutine compute_rhs
   5
   6 c---------------------------------------------------------------------
   7 c---------------------------------------------------------------------
   8
   9        include 'header.h'
  10
  11        integer c, i, j, k, m
  12        double precision aux, rho_inv, uijk, up1, um1, vijk, vp1, vm1,
  13      >                  wijk, wp1, wm1
  14
  15
  16 c---------------------------------------------------------------------
  17 c loop over all cells owned by this node
  18 c---------------------------------------------------------------------
  19        do    c = 1, ncells
  20
  21 c---------------------------------------------------------------------
  22 c         compute the reciprocal of density, and the kinetic energy,
  23 c         and the speed of sound.
  24 c---------------------------------------------------------------------
  25
  26           do    k = -1, cell_size(3,c)
  27              do    j = -1, cell_size(2,c)
  28                 do    i = -1, cell_size(1,c)
  29                    rho_inv = 1.0d0/u(i,j,k,1,c)
  30                    rho_i(i,j,k,c) = rho_inv
  31                    us(i,j,k,c) = u(i,j,k,2,c) * rho_inv
  32                    vs(i,j,k,c) = u(i,j,k,3,c) * rho_inv
  33                    ws(i,j,k,c) = u(i,j,k,4,c) * rho_inv
  34                    square(i,j,k,c)     = 0.5d0* (
  35      >                        u(i,j,k,2,c)*u(i,j,k,2,c) +
  36      >                        u(i,j,k,3,c)*u(i,j,k,3,c) +
  37      >                        u(i,j,k,4,c)*u(i,j,k,4,c) ) * rho_inv
  38                    qs(i,j,k,c) = square(i,j,k,c) * rho_inv
  39 c---------------------------------------------------------------------
  40 c                  (don't need speed and ainx until the lhs computation)
  41 c---------------------------------------------------------------------
  42                    aux = c1c2*rho_inv* (u(i,j,k,5,c) - square(i,j,k,c))
  43                    aux = dsqrt(aux)
  44                    speed(i,j,k,c) = aux
  45                    ainv(i,j,k,c)  = 1.0d0/aux
  46                 end do
  47              end do
  48           end do
  49
  50 c---------------------------------------------------------------------
  51 c copy the exact forcing term to the right hand side;  because
  52 c this forcing term is known, we can store it on the whole of every
  53 c cell,  including the boundary
  54 c---------------------------------------------------------------------
  55
  56           do   m = 1, 5
  57              do   k = 0, cell_size(3,c)-1
  58                 do   j = 0, cell_size(2,c)-1
  59                    do   i = 0, cell_size(1,c)-1
  60                       rhs(i,j,k,m,c) = forcing(i,j,k,m,c)
  61                    end do
  62                 end do
  63              end do
  64           end do
  65
  66
  67 c---------------------------------------------------------------------
  68 c         compute xi-direction fluxes
  69 c---------------------------------------------------------------------
  70           do    k = start(3,c), cell_size(3,c)-end(3,c)-1
  71              do    j = start(2,c), cell_size(2,c)-end(2,c)-1
  72                 do    i = start(1,c), cell_size(1,c)-end(1,c)-1
  73                    uijk = us(i,j,k,c)
  74                    up1  = us(i+1,j,k,c)
  75                    um1  = us(i-1,j,k,c)
  76
  77                    rhs(i,j,k,1,c) = rhs(i,j,k,1,c) + dx1tx1 *
  78      >                    (u(i+1,j,k,1,c) - 2.0d0*u(i,j,k,1,c) +
  79      >                     u(i-1,j,k,1,c)) -
  80      >                    tx2 * (u(i+1,j,k,2,c) - u(i-1,j,k,2,c))
  81
  82                    rhs(i,j,k,2,c) = rhs(i,j,k,2,c) + dx2tx1 *
  83      >                    (u(i+1,j,k,2,c) - 2.0d0*u(i,j,k,2,c) +
  84      >                     u(i-1,j,k,2,c)) +
  85      >                    xxcon2*con43 * (up1 - 2.0d0*uijk + um1) -
  86      >                    tx2 * (u(i+1,j,k,2,c)*up1 -
  87      >                           u(i-1,j,k,2,c)*um1 +
  88      >                           (u(i+1,j,k,5,c)- square(i+1,j,k,c)-
  89      >                            u(i-1,j,k,5,c)+ square(i-1,j,k,c))*
  90      >                            c2)
  91
  92                    rhs(i,j,k,3,c) = rhs(i,j,k,3,c) + dx3tx1 *
  93      >                    (u(i+1,j,k,3,c) - 2.0d0*u(i,j,k,3,c) +
  94      >                     u(i-1,j,k,3,c)) +
  95      >                    xxcon2 * (vs(i+1,j,k,c) - 2.0d0*vs(i,j,k,c) +
  96      >                              vs(i-1,j,k,c)) -
  97      >                    tx2 * (u(i+1,j,k,3,c)*up1 -
  98      >                           u(i-1,j,k,3,c)*um1)
  99
 100                    rhs(i,j,k,4,c) = rhs(i,j,k,4,c) + dx4tx1 *
 101      >                    (u(i+1,j,k,4,c) - 2.0d0*u(i,j,k,4,c) +
 102      >                     u(i-1,j,k,4,c)) +
 103      >                    xxcon2 * (ws(i+1,j,k,c) - 2.0d0*ws(i,j,k,c) +
 104      >                              ws(i-1,j,k,c)) -
 105      >                    tx2 * (u(i+1,j,k,4,c)*up1 -
 106      >                           u(i-1,j,k,4,c)*um1)
 107
 108                    rhs(i,j,k,5,c) = rhs(i,j,k,5,c) + dx5tx1 *
 109      >                    (u(i+1,j,k,5,c) - 2.0d0*u(i,j,k,5,c) +
 110      >                     u(i-1,j,k,5,c)) +
 111      >                    xxcon3 * (qs(i+1,j,k,c) - 2.0d0*qs(i,j,k,c) +
 112      >                              qs(i-1,j,k,c)) +
 113      >                    xxcon4 * (up1*up1 -       2.0d0*uijk*uijk +
 114      >                              um1*um1) +
 115      >                    xxcon5 * (u(i+1,j,k,5,c)*rho_i(i+1,j,k,c) -
 116      >                              2.0d0*u(i,j,k,5,c)*rho_i(i,j,k,c) +
 117      >                              u(i-1,j,k,5,c)*rho_i(i-1,j,k,c)) -
 118      >                    tx2 * ( (c1*u(i+1,j,k,5,c) -
 119      >                             c2*square(i+1,j,k,c))*up1 -
 120      >                            (c1*u(i-1,j,k,5,c) -
 121      >                             c2*square(i-1,j,k,c))*um1 )
 122                 end do
 123              end do
 124           end do
 125
 126 c---------------------------------------------------------------------
 127 c         add fourth order xi-direction dissipation
 128 c---------------------------------------------------------------------
 129           if (start(1,c) .gt. 0) then
 130              i = 1
 131              do    m = 1, 5
 132                 do    k = start(3,c), cell_size(3,c)-end(3,c)-1
 133                    do    j = start(2,c), cell_size(2,c)-end(2,c)-1
 134                       rhs(i,j,k,m,c) = rhs(i,j,k,m,c)- dssp *
 135      >                    ( 5.0d0*u(i,j,k,m,c) - 4.0d0*u(i+1,j,k,m,c) +
 136      >                            u(i+2,j,k,m,c))
 137                    end do
 138                 end do
 139              end do
 140
 141              i = 2
 142              do    m = 1, 5
 143                 do    k = start(3,c), cell_size(3,c)-end(3,c)-1
 144                    do    j = start(2,c), cell_size(2,c)-end(2,c)-1
 145                       rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
 146      >                    (-4.0d0*u(i-1,j,k,m,c) + 6.0d0*u(i,j,k,m,c) -
 147      >                      4.0d0*u(i+1,j,k,m,c) + u(i+2,j,k,m,c))
 148                    end do
 149                 end do
 150              end do
 151           endif
 152
 153           do     m = 1, 5
 154              do     k = start(3,c), cell_size(3,c)-end(3,c)-1
 155                 do     j = start(2,c), cell_size(2,c)-end(2,c)-1
 156                    do  i = 3*start(1,c),cell_size(1,c)-3*end(1,c)-1
 157                       rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
 158      >                    (  u(i-2,j,k,m,c) - 4.0d0*u(i-1,j,k,m,c) +
 159      >                     6.0*u(i,j,k,m,c) - 4.0d0*u(i+1,j,k,m,c) +
 160      >                         u(i+2,j,k,m,c) )
 161                    end do
 162                 end do
 163              end do
 164           end do
 165
 166
 167           if (end(1,c) .gt. 0) then
 168              i = cell_size(1,c)-3
 169              do     m = 1, 5
 170                 do     k = start(3,c), cell_size(3,c)-end(3,c)-1
 171                    do     j = start(2,c), cell_size(2,c)-end(2,c)-1
 172                       rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
 173      >                    ( u(i-2,j,k,m,c) - 4.0d0*u(i-1,j,k,m,c) +
 174      >                      6.0d0*u(i,j,k,m,c) - 4.0d0*u(i+1,j,k,m,c) )
 175                    end do
 176                 end do
 177              end do
 178
 179              i = cell_size(1,c)-2
 180              do     m = 1, 5
 181                 do     k = start(3,c), cell_size(3,c)-end(3,c)-1
 182                    do     j = start(2,c), cell_size(2,c)-end(2,c)-1
 183                       rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
 184      >                    ( u(i-2,j,k,m,c) - 4.d0*u(i-1,j,k,m,c) +
 185      >                      5.d0*u(i,j,k,m,c) )
 186                    end do
 187                 end do
 188              end do
 189           endif
 190
 191 c---------------------------------------------------------------------
 192 c         compute eta-direction fluxes
 193 c---------------------------------------------------------------------
 194           do     k = start(3,c), cell_size(3,c)-end(3,c)-1
 195              do     j = start(2,c), cell_size(2,c)-end(2,c)-1
 196                 do     i = start(1,c), cell_size(1,c)-end(1,c)-1
 197                    vijk = vs(i,j,k,c)
 198                    vp1  = vs(i,j+1,k,c)
 199                    vm1  = vs(i,j-1,k,c)
 200                    rhs(i,j,k,1,c) = rhs(i,j,k,1,c) + dy1ty1 *
 201      >                   (u(i,j+1,k,1,c) - 2.0d0*u(i,j,k,1,c) +
 202      >                    u(i,j-1,k,1,c)) -
 203      >                   ty2 * (u(i,j+1,k,3,c) - u(i,j-1,k,3,c))
 204                    rhs(i,j,k,2,c) = rhs(i,j,k,2,c) + dy2ty1 *
 205      >                   (u(i,j+1,k,2,c) - 2.0d0*u(i,j,k,2,c) +
 206      >                    u(i,j-1,k,2,c)) +
 207      >                   yycon2 * (us(i,j+1,k,c) - 2.0d0*us(i,j,k,c) +
 208      >                             us(i,j-1,k,c)) -
 209      >                   ty2 * (u(i,j+1,k,2,c)*vp1 -
 210      >                          u(i,j-1,k,2,c)*vm1)
 211                    rhs(i,j,k,3,c) = rhs(i,j,k,3,c) + dy3ty1 *
 212      >                   (u(i,j+1,k,3,c) - 2.0d0*u(i,j,k,3,c) +
 213      >                    u(i,j-1,k,3,c)) +
 214      >                   yycon2*con43 * (vp1 - 2.0d0*vijk + vm1) -
 215      >                   ty2 * (u(i,j+1,k,3,c)*vp1 -
 216      >                          u(i,j-1,k,3,c)*vm1 +
 217      >                          (u(i,j+1,k,5,c) - square(i,j+1,k,c) -
 218      >                           u(i,j-1,k,5,c) + square(i,j-1,k,c))
 219      >                          *c2)
 220                    rhs(i,j,k,4,c) = rhs(i,j,k,4,c) + dy4ty1 *
 221      >                   (u(i,j+1,k,4,c) - 2.0d0*u(i,j,k,4,c) +
 222      >                    u(i,j-1,k,4,c)) +
 223      >                   yycon2 * (ws(i,j+1,k,c) - 2.0d0*ws(i,j,k,c) +
 224      >                             ws(i,j-1,k,c)) -
 225      >                   ty2 * (u(i,j+1,k,4,c)*vp1 -
 226      >                          u(i,j-1,k,4,c)*vm1)
 227                    rhs(i,j,k,5,c) = rhs(i,j,k,5,c) + dy5ty1 *
 228      >                   (u(i,j+1,k,5,c) - 2.0d0*u(i,j,k,5,c) +
 229      >                    u(i,j-1,k,5,c)) +
 230      >                   yycon3 * (qs(i,j+1,k,c) - 2.0d0*qs(i,j,k,c) +
 231      >                             qs(i,j-1,k,c)) +
 232      >                   yycon4 * (vp1*vp1       - 2.0d0*vijk*vijk +
 233      >                             vm1*vm1) +
 234      >                   yycon5 * (u(i,j+1,k,5,c)*rho_i(i,j+1,k,c) -
 235      >                             2.0d0*u(i,j,k,5,c)*rho_i(i,j,k,c) +
 236      >                             u(i,j-1,k,5,c)*rho_i(i,j-1,k,c)) -
 237      >                   ty2 * ((c1*u(i,j+1,k,5,c) -
 238      >                           c2*square(i,j+1,k,c)) * vp1 -
 239      >                          (c1*u(i,j-1,k,5,c) -
 240      >                           c2*square(i,j-1,k,c)) * vm1)
 241                 end do
 242              end do
 243           end do
 244
 245 c---------------------------------------------------------------------
 246 c         add fourth order eta-direction dissipation
 247 c---------------------------------------------------------------------
 248           if (start(2,c) .gt. 0) then
 249              j = 1
 250              do     m = 1, 5
 251                 do     k = start(3,c), cell_size(3,c)-end(3,c)-1
 252                    do     i = start(1,c), cell_size(1,c)-end(1,c)-1
 253                       rhs(i,j,k,m,c) = rhs(i,j,k,m,c)- dssp *
 254      >                    ( 5.0d0*u(i,j,k,m,c) - 4.0d0*u(i,j+1,k,m,c) +
 255      >                            u(i,j+2,k,m,c))
 256                    end do
 257                 end do
 258              end do
 259
 260              j = 2
 261              do     m = 1, 5
 262                 do     k = start(3,c), cell_size(3,c)-end(3,c)-1
 263                    do     i = start(1,c), cell_size(1,c)-end(1,c)-1
 264                       rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
 265      >                    (-4.0d0*u(i,j-1,k,m,c) + 6.0d0*u(i,j,k,m,c) -
 266      >                      4.0d0*u(i,j+1,k,m,c) + u(i,j+2,k,m,c))
 267                    end do
 268                 end do
 269              end do
 270           endif
 271
 272           do     m = 1, 5
 273              do     k = start(3,c), cell_size(3,c)-end(3,c)-1
 274                 do    j = 3*start(2,c), cell_size(2,c)-3*end(2,c)-1
 275                    do  i = start(1,c),cell_size(1,c)-end(1,c)-1
 276                       rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
 277      >                    (  u(i,j-2,k,m,c) - 4.0d0*u(i,j-1,k,m,c) +
 278      >                     6.0*u(i,j,k,m,c) - 4.0d0*u(i,j+1,k,m,c) +
 279      >                         u(i,j+2,k,m,c) )
 280                    end do
 281                 end do
 282              end do
 283           end do
 284
 285           if (end(2,c) .gt. 0) then
 286              j = cell_size(2,c)-3
 287              do     m = 1, 5
 288                 do     k = start(3,c), cell_size(3,c)-end(3,c)-1
 289                    do     i = start(1,c), cell_size(1,c)-end(1,c)-1
 290                       rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
 291      >                    ( u(i,j-2,k,m,c) - 4.0d0*u(i,j-1,k,m,c) +
 292      >                      6.0d0*u(i,j,k,m,c) - 4.0d0*u(i,j+1,k,m,c) )
 293                    end do
 294                 end do
 295              end do
 296
 297              j = cell_size(2,c)-2
 298              do     m = 1, 5
 299                 do     k = start(3,c), cell_size(3,c)-end(3,c)-1
 300                    do     i = start(1,c), cell_size(1,c)-end(1,c)-1
 301                       rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
 302      >                    ( u(i,j-2,k,m,c) - 4.d0*u(i,j-1,k,m,c) +
 303      >                      5.d0*u(i,j,k,m,c) )
 304                    end do
 305                 end do
 306              end do
 307           endif
 308
 309
 310 c---------------------------------------------------------------------
 311 c         compute zeta-direction fluxes
 312 c---------------------------------------------------------------------
 313           do    k = start(3,c), cell_size(3,c)-end(3,c)-1
 314              do     j = start(2,c), cell_size(2,c)-end(2,c)-1
 315                 do     i = start(1,c), cell_size(1,c)-end(1,c)-1
 316                    wijk = ws(i,j,k,c)
 317                    wp1  = ws(i,j,k+1,c)
 318                    wm1  = ws(i,j,k-1,c)
 319
 320                    rhs(i,j,k,1,c) = rhs(i,j,k,1,c) + dz1tz1 *
 321      >                   (u(i,j,k+1,1,c) - 2.0d0*u(i,j,k,1,c) +
 322      >                    u(i,j,k-1,1,c)) -
 323      >                   tz2 * (u(i,j,k+1,4,c) - u(i,j,k-1,4,c))
 324                    rhs(i,j,k,2,c) = rhs(i,j,k,2,c) + dz2tz1 *
 325      >                   (u(i,j,k+1,2,c) - 2.0d0*u(i,j,k,2,c) +
 326      >                    u(i,j,k-1,2,c)) +
 327      >                   zzcon2 * (us(i,j,k+1,c) - 2.0d0*us(i,j,k,c) +
 328      >                             us(i,j,k-1,c)) -
 329      >                   tz2 * (u(i,j,k+1,2,c)*wp1 -
 330      >                          u(i,j,k-1,2,c)*wm1)
 331                    rhs(i,j,k,3,c) = rhs(i,j,k,3,c) + dz3tz1 *
 332      >                   (u(i,j,k+1,3,c) - 2.0d0*u(i,j,k,3,c) +
 333      >                    u(i,j,k-1,3,c)) +
 334      >                   zzcon2 * (vs(i,j,k+1,c) - 2.0d0*vs(i,j,k,c) +
 335      >                             vs(i,j,k-1,c)) -
 336      >                   tz2 * (u(i,j,k+1,3,c)*wp1 -
 337      >                          u(i,j,k-1,3,c)*wm1)
 338                    rhs(i,j,k,4,c) = rhs(i,j,k,4,c) + dz4tz1 *
 339      >                   (u(i,j,k+1,4,c) - 2.0d0*u(i,j,k,4,c) +
 340      >                    u(i,j,k-1,4,c)) +
 341      >                   zzcon2*con43 * (wp1 - 2.0d0*wijk + wm1) -
 342      >                   tz2 * (u(i,j,k+1,4,c)*wp1 -
 343      >                          u(i,j,k-1,4,c)*wm1 +
 344      >                          (u(i,j,k+1,5,c) - square(i,j,k+1,c) -
 345      >                           u(i,j,k-1,5,c) + square(i,j,k-1,c))
 346      >                          *c2)
 347                    rhs(i,j,k,5,c) = rhs(i,j,k,5,c) + dz5tz1 *
 348      >                   (u(i,j,k+1,5,c) - 2.0d0*u(i,j,k,5,c) +
 349      >                    u(i,j,k-1,5,c)) +
 350      >                   zzcon3 * (qs(i,j,k+1,c) - 2.0d0*qs(i,j,k,c) +
 351      >                             qs(i,j,k-1,c)) +
 352      >                   zzcon4 * (wp1*wp1 - 2.0d0*wijk*wijk +
 353      >                             wm1*wm1) +
 354      >                   zzcon5 * (u(i,j,k+1,5,c)*rho_i(i,j,k+1,c) -
 355      >                             2.0d0*u(i,j,k,5,c)*rho_i(i,j,k,c) +
 356      >                             u(i,j,k-1,5,c)*rho_i(i,j,k-1,c)) -
 357      >                   tz2 * ( (c1*u(i,j,k+1,5,c) -
 358      >                            c2*square(i,j,k+1,c))*wp1 -
 359      >                           (c1*u(i,j,k-1,5,c) -
 360      >                            c2*square(i,j,k-1,c))*wm1)
 361                 end do
 362              end do
 363           end do
 364
 365 c---------------------------------------------------------------------
 366 c         add fourth order zeta-direction dissipation
 367 c---------------------------------------------------------------------
 368           if (start(3,c) .gt. 0) then
 369              k = 1
 370              do     m = 1, 5
 371                 do     j = start(2,c), cell_size(2,c)-end(2,c)-1
 372                    do     i = start(1,c), cell_size(1,c)-end(1,c)-1
 373                       rhs(i,j,k,m,c) = rhs(i,j,k,m,c)- dssp *
 374      >                    ( 5.0d0*u(i,j,k,m,c) - 4.0d0*u(i,j,k+1,m,c) +
 375      >                            u(i,j,k+2,m,c))
 376                    end do
 377                 end do
 378              end do
 379
 380              k = 2
 381              do     m = 1, 5
 382                 do     j = start(2,c), cell_size(2,c)-end(2,c)-1
 383                    do     i = start(1,c), cell_size(1,c)-end(1,c)-1
 384                       rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
 385      >                    (-4.0d0*u(i,j,k-1,m,c) + 6.0d0*u(i,j,k,m,c) -
 386      >                      4.0d0*u(i,j,k+1,m,c) + u(i,j,k+2,m,c))
 387                    end do
 388                 end do
 389              end do
 390           endif
 391
 392           do     m = 1, 5
 393              do     k = 3*start(3,c), cell_size(3,c)-3*end(3,c)-1
 394                 do     j = start(2,c), cell_size(2,c)-end(2,c)-1
 395                    do     i = start(1,c),cell_size(1,c)-end(1,c)-1
 396                       rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
 397      >                    (  u(i,j,k-2,m,c) - 4.0d0*u(i,j,k-1,m,c) +
 398      >                     6.0*u(i,j,k,m,c) - 4.0d0*u(i,j,k+1,m,c) +
 399      >                         u(i,j,k+2,m,c) )
 400                    end do
 401                 end do
 402              end do
 403           end do
 404
 405           if (end(3,c) .gt. 0) then
 406              k = cell_size(3,c)-3
 407              do     m = 1, 5
 408                 do     j = start(2,c), cell_size(2,c)-end(2,c)-1
 409                    do     i = start(1,c), cell_size(1,c)-end(1,c)-1
 410                       rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
 411      >                    ( u(i,j,k-2,m,c) - 4.0d0*u(i,j,k-1,m,c) +
 412      >                      6.0d0*u(i,j,k,m,c) - 4.0d0*u(i,j,k+1,m,c) )
 413                    end do
 414                 end do
 415              end do
 416
 417              k = cell_size(3,c)-2
 418              do     m = 1, 5
 419                 do     j = start(2,c), cell_size(2,c)-end(2,c)-1
 420                    do     i = start(1,c), cell_size(1,c)-end(1,c)-1
 421                       rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - dssp *
 422      >                    ( u(i,j,k-2,m,c) - 4.d0*u(i,j,k-1,m,c) +
 423      >                      5.d0*u(i,j,k,m,c) )
 424                    end do
 425                 end do
 426              end do
 427           endif
 428
 429           do     m = 1, 5
 430              do     k = start(3,c), cell_size(3,c)-end(3,c)-1
 431                 do     j = start(2,c), cell_size(2,c)-end(2,c)-1
 432                    do    i = start(1,c), cell_size(1,c)-end(1,c)-1
 433                       rhs(i,j,k,m,c) = rhs(i,j,k,m,c) * dt
 434                    end do
 435                 end do
 436              end do
 437           end do
 438
 439        end do
 440
 441        return
 442        end
 443
 444
 445
 446