Ignore:
Timestamp:
Apr 11, 2014 5:15:14 PM (7 years ago)
Author:
hoffmann
Message:

new Lagrangian particle structure integrated

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/lpm_set_attributes.f90

    r1321 r1359  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! New particle structure integrated.
     23! Kind definition added to all floating point numbers.
    2324!
    2425! Former revisions:
     
    7778
    7879    USE particle_attributes,                                                   &
    79         ONLY:  number_of_particles, offset_ocean_nzt, particles
     80        ONLY:  block_offset, grid_particles, number_of_particles,              &
     81               offset_ocean_nzt, particles, prt_count
    8082
    8183    USE pegrid
     
    8789
    8890    INTEGER(iwp) ::  i        !:
     91    INTEGER(iwp) ::  ip       !:
    8992    INTEGER(iwp) ::  j        !:
     93    INTEGER(iwp) ::  jp       !:
    9094    INTEGER(iwp) ::  k        !:
     95    INTEGER(iwp) ::  kp       !:
    9196    INTEGER(iwp) ::  n        !:
     97    INTEGER(iwp) ::  nb       !:
     98
     99    INTEGER(iwp), DIMENSION(0:7) ::  start_index !:
     100    INTEGER(iwp), DIMENSION(0:7) ::  end_index   !:
    92101
    93102    REAL(wp)    ::  aa        !:
     
    101110    REAL(wp)    ::  pt_int_l  !:
    102111    REAL(wp)    ::  pt_int_u  !:
    103     REAL(wp)    ::  u_int     !:
    104112    REAL(wp)    ::  u_int_l   !:
    105113    REAL(wp)    ::  u_int_u   !:
    106     REAL(wp)    ::  v_int     !:
    107114    REAL(wp)    ::  v_int_l   !:
    108115    REAL(wp)    ::  v_int_u   !:
     
    113120    REAL(wp)    ::  y         !:
    114121
     122    REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_int                  !:
     123    REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_int                  !:
     124    REAL(wp), DIMENSION(:), ALLOCATABLE ::  xv                     !:
     125    REAL(wp), DIMENSION(:), ALLOCATABLE ::  yv                     !:
     126    REAL(wp), DIMENSION(:), ALLOCATABLE ::  zv                     !:
     127
    115128    CALL cpu_log( log_point_s(49), 'lpm_set_attributes', 'start' )
    116129
     
    122135!--    Set particle color depending on the absolute value of the horizontal
    123136!--    velocity
    124        DO  n = 1, number_of_particles
    125 !
    126 !--       Interpolate u velocity-component, determine left, front, bottom
    127 !--       index of u-array
    128           i = ( particles(n)%x + 0.5 * dx ) * ddx
    129           j =   particles(n)%y * ddy
    130           k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
    131               + offset_ocean_nzt                     ! only exact if equidistant
    132 
    133 !
    134 !--       Interpolation of the velocity components in the xy-plane
    135           x  = particles(n)%x + ( 0.5 - i ) * dx
    136           y  = particles(n)%y - j * dy
    137           aa = x**2          + y**2
    138           bb = ( dx - x )**2 + y**2
    139           cc = x**2          + ( dy - y )**2
    140           dd = ( dx - x )**2 + ( dy - y )**2
    141           gg = aa + bb + cc + dd
    142 
    143           u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
    144                     + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1) &
    145                     ) / ( 3.0 * gg ) - u_gtrans
    146           IF ( k+1 == nzt+1 )  THEN
    147              u_int = u_int_l
    148           ELSE
    149              u_int_u = ( ( gg-aa ) * u(k+1,j,i)   + ( gg-bb ) * u(k+1,j,i+1)   &
    150                        + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1) &
    151                        ) / ( 3.0 * gg ) - u_gtrans
    152              u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * &
    153                                ( u_int_u - u_int_l )
    154           ENDIF
    155 
    156 !
    157 !--       Same procedure for interpolation of the v velocity-component (adopt
    158 !--       index k from u velocity-component)
    159           i =   particles(n)%x * ddx
    160           j = ( particles(n)%y + 0.5 * dy ) * ddy
    161 
    162           x  = particles(n)%x - i * dx
    163           y  = particles(n)%y + ( 0.5 - j ) * dy
    164           aa = x**2          + y**2
    165           bb = ( dx - x )**2 + y**2
    166           cc = x**2          + ( dy - y )**2
    167           dd = ( dx - x )**2 + ( dy - y )**2
    168           gg = aa + bb + cc + dd
    169 
    170           v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
    171                  + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
    172                  ) / ( 3.0 * gg ) - v_gtrans
    173           IF ( k+1 == nzt+1 )  THEN
    174              v_int = v_int_l
    175           ELSE
    176              v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
    177                        + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
    178                        ) / ( 3.0 * gg ) - v_gtrans
    179              v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * &
    180                                ( v_int_u - v_int_l )
    181           ENDIF
    182 
    183           absuv = SQRT( u_int**2 + v_int**2 )
    184 
    185 !
    186 !--       Limit values by the given interval and normalize to interval [0,1]
    187           absuv = MIN( absuv, color_interval(2) )
    188           absuv = MAX( absuv, color_interval(1) )
    189 
    190           absuv = ( absuv - color_interval(1) ) / &
    191                   ( color_interval(2) - color_interval(1) )
    192 
    193 !
    194 !--       Number of available colors is defined in init_dvrp
    195           particles(n)%class = 1 + absuv * ( dvrp_colortable_entries_prt - 1 )
    196 
     137       DO  ip = nxl, nxr
     138          DO  jp = nys, nyn
     139             DO  kp = nzb+1, nzt
     140
     141                number_of_particles = prt_count(kp,jp,ip)
     142                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
     143                IF ( number_of_particles <= 0 )  CYCLE
     144                start_index = grid_particles(kp,jp,ip)%start_index
     145                end_index   = grid_particles(kp,jp,ip)%end_index
     146
     147                ALLOCATE( u_int(1:number_of_particles), &
     148                          v_int(1:number_of_particles), &
     149                          xv(1:number_of_particles),    &
     150                          yv(1:number_of_particles),    &
     151                          zv(1:number_of_particles) )
     152
     153                xv = particles(1:number_of_particles)%x
     154                yv = particles(1:number_of_particles)%y
     155                zv = particles(1:number_of_particles)%z
     156
     157                DO  nb = 0,7
     158
     159                   i = ip
     160                   j = jp + block_offset(nb)%j_off
     161                   k = kp + block_offset(nb)%k_off
     162
     163                   DO  n = start_index(nb), end_index(nb)
     164!
     165!--                   Interpolation of the velocity components in the xy-plane
     166                      x  = xv(n) + ( 0.5_wp - i ) * dx
     167                      y  = yv(n) - j * dy
     168                      aa = x**2          + y**2
     169                      bb = ( dx - x )**2 + y**2
     170                      cc = x**2          + ( dy - y )**2
     171                      dd = ( dx - x )**2 + ( dy - y )**2
     172                      gg = aa + bb + cc + dd
     173
     174                      u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) *     &
     175                                  u(k,j,i+1) + ( gg - cc ) * u(k,j+1,i) +      &
     176                                  ( gg - dd ) * u(k,j+1,i+1)                   &
     177                                ) / ( 3.0_wp * gg ) - u_gtrans
     178
     179                      IF ( k+1 == nzt+1 )  THEN
     180                         u_int(n) = u_int_l
     181                      ELSE
     182                         u_int_u = ( ( gg - aa ) * u(k+1,j,i)   + ( gg - bb ) *  &
     183                                     u(k+1,j,i+1) + ( gg - cc ) * u(k+1,j+1,i) + &
     184                                     ( gg - dd ) * u(k+1,j+1,i+1)                &
     185                                   ) / ( 3.0_wp * gg ) - u_gtrans
     186                         u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dz *         &
     187                                           ( u_int_u - u_int_l )
     188                      ENDIF
     189
     190                   ENDDO
     191
     192                   i = ip + block_offset(nb)%i_off
     193                   j = jp
     194                   k = kp + block_offset(nb)%k_off
     195
     196                   DO  n = start_index(nb), end_index(nb)
     197!
     198!--                   Same procedure for interpolation of the v velocity-component
     199                      x  = xv(n) - i * dx
     200                      y  = yv(n) + ( 0.5_wp - j ) * dy
     201                      aa = x**2          + y**2
     202                      bb = ( dx - x )**2 + y**2
     203                      cc = x**2          + ( dy - y )**2
     204                      dd = ( dx - x )**2 + ( dy - y )**2
     205                      gg = aa + bb + cc + dd
     206
     207                      v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) *     &
     208                                  v(k,j,i+1) + ( gg - cc ) * v(k,j+1,i) +      &
     209                                  ( gg - dd ) * v(k,j+1,i+1)                   &
     210                                ) / ( 3.0_wp * gg ) - v_gtrans
     211
     212                      IF ( k+1 == nzt+1 )  THEN
     213                         v_int(n) = v_int_l
     214                      ELSE
     215                         v_int_u  = ( ( gg - aa ) * v(k+1,j,i) + ( gg - bb ) *    &
     216                                      v(k+1,j,i+1) + ( gg - cc ) * v(k+1,j+1,i) + &
     217                                      ( gg - dd ) * v(k+1,j+1,i+1)                &
     218                                    ) / ( 3.0_wp * gg ) - v_gtrans
     219                         v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dz *         &
     220                                           ( v_int_u - v_int_l )
     221                      ENDIF
     222
     223                   ENDDO
     224
     225                ENDDO
     226
     227                DO  n = 1, number_of_particles
     228
     229                   absuv = SQRT( u_int(n)**2 + v_int(n)**2 )
     230
     231!
     232!--                Limit values by the given interval and normalize to
     233!--                interval [0,1]
     234                   absuv = MIN( absuv, color_interval(2) )
     235                   absuv = MAX( absuv, color_interval(1) )
     236
     237                   absuv = ( absuv - color_interval(1) ) / &
     238                           ( color_interval(2) - color_interval(1) )
     239
     240!
     241!--                Number of available colors is defined in init_dvrp
     242                   particles(n)%class = 1 + absuv *                            &
     243                                            ( dvrp_colortable_entries_prt - 1 )
     244
     245                ENDDO
     246
     247                DEALLOCATE( u_int, v_int, xv, yv, zv )
     248
     249             ENDDO
     250          ENDDO
    197251       ENDDO
    198252
     
    204258!--    (This is also done in flow_statistics, but flow_statistics is called
    205259!--    after this routine.)
    206        sums_l(:,4,0) = 0.0
     260       sums_l(:,4,0) = 0.0_wp
    207261       DO  i = nxl, nxr
    208262          DO  j =  nys, nyn
     
    214268
    215269#if defined( __parallel )
    216 
    217270       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    218271       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, &
    219272                           MPI_REAL, MPI_SUM, comm2d, ierr )
    220 
    221273#else
    222 
    223274       sums(:,4) = sums_l(:,4,0)
    224 
    225275#endif
    226 
    227276       sums(:,4) = sums(:,4) / ngp_2dh(0)
    228277
    229        DO  n = 1, number_of_particles
    230 !
    231 !--       Interpolate temperature to the current particle position
    232           i = particles(n)%x * ddx
    233           j = particles(n)%y * ddy
    234           k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
    235               + offset_ocean_nzt                     ! only exact if equidistant
    236 
    237           x  = particles(n)%x - i * dx
    238           y  = particles(n)%y - j * dy
    239           aa = x**2          + y**2
    240           bb = ( dx - x )**2 + y**2
    241           cc = x**2          + ( dy - y )**2
    242           dd = ( dx - x )**2 + ( dy - y )**2
    243           gg = aa + bb + cc + dd
    244 
    245           pt_int_l = ( ( gg - aa ) * pt(k,j,i)   + ( gg - bb ) * pt(k,j,i+1)   &
    246                      + ( gg - cc ) * pt(k,j+1,i) + ( gg - dd ) * pt(k,j+1,i+1) &
    247                      ) / ( 3.0 * gg ) - sums(k,4)
    248 
    249           pt_int_u = ( ( gg-aa ) * pt(k+1,j,i)   + ( gg-bb ) * pt(k+1,j,i+1)   &
    250                      + ( gg-cc ) * pt(k+1,j+1,i) + ( gg-dd ) * pt(k+1,j+1,i+1) &
    251                      ) / ( 3.0 * gg ) - sums(k,4)
    252 
    253           pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz * &
    254                               ( pt_int_u - pt_int_l )
    255 
    256 !
    257 !--       Limit values by the given interval and normalize to interval [0,1]
    258           pt_int = MIN( pt_int, color_interval(2) )
    259           pt_int = MAX( pt_int, color_interval(1) )
    260 
    261           pt_int = ( pt_int - color_interval(1) ) / &
    262                    ( color_interval(2) - color_interval(1) )
    263 
    264 !
    265 !--       Number of available colors is defined in init_dvrp
    266           particles(n)%class = 1 + pt_int * ( dvrp_colortable_entries_prt - 1 )
    267 
     278       DO  ip = nxl, nxr
     279          DO  jp = nys, nyn
     280             DO  kp = nzb+1, nzt
     281
     282                number_of_particles = prt_count(kp,jp,ip)
     283                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
     284                IF ( number_of_particles <= 0 )  CYCLE
     285                start_index = grid_particles(kp,jp,ip)%start_index
     286                end_index   = grid_particles(kp,jp,ip)%end_index
     287
     288                ALLOCATE( xv(1:number_of_particles), &
     289                          yv(1:number_of_particles), &
     290                          zv(1:number_of_particles) )
     291
     292                xv = particles(1:number_of_particles)%x
     293                yv = particles(1:number_of_particles)%y
     294                zv = particles(1:number_of_particles)%z
     295
     296                DO  nb = 0,7
     297
     298                   i = ip + block_offset(nb)%i_off
     299                   j = jp + block_offset(nb)%j_off
     300                   k = kp + block_offset(nb)%k_off
     301
     302                   DO  n = start_index(nb), end_index(nb)
     303!
     304!--                   Interpolate temperature to the current particle position
     305                      x  = xv(n) - i * dx
     306                      y  = yv(n) - j * dy
     307                      aa = x**2          + y**2
     308                      bb = ( dx - x )**2 + y**2
     309                      cc = x**2          + ( dy - y )**2
     310                      dd = ( dx - x )**2 + ( dy - y )**2
     311                      gg = aa + bb + cc + dd
     312
     313                      pt_int_l = ( ( gg - aa ) * pt(k,j,i)   + ( gg - bb ) *   &
     314                                   pt(k,j,i+1) + ( gg - cc ) * pt(k,j+1,i) +   &
     315                                   ( gg - dd ) * pt(k,j+1,i+1)                 &
     316                                 ) / ( 3.0_wp * gg ) - sums(k,4)
     317
     318                      pt_int_u = ( ( gg - aa ) * pt(k+1,j,i) + ( gg - bb ) *     &
     319                                   pt(k+1,j,i+1) + ( gg - cc ) * pt(k+1,j+1,i) + &
     320                                   ( gg - dd ) * pt(k+1,j+1,i+1)                 &
     321                                 ) / ( 3.0_wp * gg ) - sums(k,4)
     322
     323                      pt_int = pt_int_l + ( zv(n) - zu(k) ) / dz *    &
     324                                          ( pt_int_u - pt_int_l )
     325
     326!
     327!--                   Limit values by the given interval and normalize to
     328!--                   interval [0,1]
     329                      pt_int = MIN( pt_int, color_interval(2) )
     330                      pt_int = MAX( pt_int, color_interval(1) )
     331
     332                      pt_int = ( pt_int - color_interval(1) ) /                &
     333                               ( color_interval(2) - color_interval(1) )
     334
     335!
     336!--                   Number of available colors is defined in init_dvrp
     337                      particles(n)%class = 1 + pt_int *                        &
     338                                           ( dvrp_colortable_entries_prt - 1 )
     339
     340                   ENDDO
     341                ENDDO
     342
     343                DEALLOCATE( xv, yv, zv )
     344
     345             ENDDO
     346          ENDDO
    268347       ENDDO
    269348
     
    272351!--    Set particle color depending on the height above the bottom
    273352!--    boundary (z=0)
    274        DO  n = 1, number_of_particles
    275 
    276           height = particles(n)%z
    277 !
    278 !--       Limit values by the given interval and normalize to interval [0,1]
    279           height = MIN( height, color_interval(2) )
    280           height = MAX( height, color_interval(1) )
    281 
    282           height = ( height - color_interval(1) ) / &
    283                    ( color_interval(2) - color_interval(1) )
    284 
    285 !
    286 !--       Number of available colors is defined in init_dvrp
    287           particles(n)%class = 1 + height * ( dvrp_colortable_entries_prt - 1 )
    288 
     353       DO  ip = nxl, nxr
     354          DO  jp = nys, nyn
     355             DO  kp = nzb+1, nzt
     356
     357                number_of_particles = prt_count(kp,jp,ip)
     358                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
     359                IF ( number_of_particles <= 0 )  CYCLE
     360                DO  n = 1, number_of_particles
     361
     362                   height = particles(n)%z
     363!
     364!--                Limit values by the given interval and normalize to
     365!--                interval [0,1]
     366                   height = MIN( height, color_interval(2) )
     367                   height = MAX( height, color_interval(1) )
     368
     369                   height = ( height - color_interval(1) ) / &
     370                            ( color_interval(2) - color_interval(1) )
     371
     372!
     373!--                Number of available colors is defined in init_dvrp
     374                   particles(n)%class = 1 + height *                           &
     375                                            ( dvrp_colortable_entries_prt - 1 )
     376
     377                ENDDO
     378
     379             ENDDO
     380          ENDDO
    289381       ENDDO
    290382
     
    295387    IF ( particle_dvrpsize == 'absw' )  THEN
    296388
    297        DO  n = 1, number_of_particles
    298 !
    299 !--       Interpolate w-component to the current particle position
    300           i = particles(n)%x * ddx
    301           j = particles(n)%y * ddy
    302           k = particles(n)%z / dz
    303 
    304           x  = particles(n)%x - i * dx
    305           y  = particles(n)%y - j * dy
    306           aa = x**2          + y**2
    307           bb = ( dx - x )**2 + y**2
    308           cc = x**2          + ( dy - y )**2
    309           dd = ( dx - x )**2 + ( dy - y )**2
    310           gg = aa + bb + cc + dd
    311 
    312           w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)   &
    313                     + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
    314                     ) / ( 3.0 * gg )
    315 
    316           IF ( k+1 == nzt+1 )  THEN
    317              w_int = w_int_l
    318           ELSE
    319              w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + ( gg-bb ) * w(k+1,j,i+1)   &
    320                        + ( gg-cc ) * w(k+1,j+1,i) + ( gg-dd ) * w(k+1,j+1,i+1) &
    321                        ) / ( 3.0 * gg )
    322              w_int = w_int_l + ( particles(n)%z - zw(k) ) / dz * &
    323                                ( w_int_u - w_int_l )
    324           ENDIF
    325 
    326 !
    327 !--       Limit values by the given interval and normalize to interval [0,1]
    328           w_int = ABS( w_int )
    329           w_int = MIN( w_int, dvrpsize_interval(2) )
    330           w_int = MAX( w_int, dvrpsize_interval(1) )
    331 
    332           w_int = ( w_int - dvrpsize_interval(1) ) / &
    333                   ( dvrpsize_interval(2) - dvrpsize_interval(1) )
    334 
    335           particles(n)%dvrp_psize = ( 0.25 + w_int * 0.6 ) * dx
    336 
     389       DO  ip = nxl, nxr
     390          DO  jp = nys, nyn
     391             DO  kp = nzb+1, nzt
     392
     393                number_of_particles = prt_count(kp,jp,ip)
     394                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
     395                IF ( number_of_particles <= 0 )  CYCLE
     396                start_index = grid_particles(kp,jp,ip)%start_index
     397                end_index   = grid_particles(kp,jp,ip)%end_index
     398
     399                ALLOCATE( xv(1:number_of_particles), &
     400                          yv(1:number_of_particles) )
     401
     402                xv = particles(1:number_of_particles)%x
     403                yv = particles(1:number_of_particles)%y
     404                zv = particles(1:number_of_particles)%z
     405
     406                DO  nb = 0,7
     407
     408                   i = ip + block_offset(nb)%i_off
     409                   j = jp + block_offset(nb)%j_off
     410                   k = kp-1
     411
     412                   DO  n = start_index(nb), end_index(nb)
     413!
     414!--                   Interpolate w-component to the current particle position
     415                      x  = xv(n) - i * dx
     416                      y  = yv(n) - j * dy
     417                      aa = x**2          + y**2
     418                      bb = ( dx - x )**2 + y**2
     419                      cc = x**2          + ( dy - y )**2
     420                      dd = ( dx - x )**2 + ( dy - y )**2
     421                      gg = aa + bb + cc + dd
     422
     423                      w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) *     &
     424                                  w(k,j,i+1) + ( gg - cc ) * w(k,j+1,i) +      &
     425                                  ( gg - dd ) * w(k,j+1,i+1)                   &
     426                                ) / ( 3.0_wp * gg )
     427
     428                      IF ( k+1 == nzt+1 )  THEN
     429                         w_int = w_int_l
     430                      ELSE
     431                         w_int_u = ( ( gg - aa ) * w(k+1,j,i)   + ( gg - bb ) *  &
     432                                     w(k+1,j,i+1) + ( gg - cc ) * w(k+1,j+1,i) + &
     433                                     ( gg - dd ) * w(k+1,j+1,i+1)                &
     434                                   ) / ( 3.0_wp * gg )
     435                         w_int = w_int_l + ( zv(n) - zw(k) ) / dz *     &
     436                                           ( w_int_u - w_int_l )
     437                      ENDIF
     438
     439!
     440!--                   Limit values by the given interval and normalize to
     441!--                   interval [0,1]
     442                      w_int = ABS( w_int )
     443                      w_int = MIN( w_int, dvrpsize_interval(2) )
     444                      w_int = MAX( w_int, dvrpsize_interval(1) )
     445
     446                      w_int = ( w_int - dvrpsize_interval(1) ) / &
     447                              ( dvrpsize_interval(2) - dvrpsize_interval(1) )
     448
     449                      particles(n)%dvrp_psize = ( 0.25_wp + w_int * 0.6_wp ) * &
     450                                                dx
     451
     452                   ENDDO
     453                ENDDO
     454
     455                DEALLOCATE( xv, yv, zv )
     456
     457             ENDDO
     458          ENDDO
    337459       ENDDO
    338460
Note: See TracChangeset for help on using the changeset viewer.