Changeset 1359 for palm/trunk/SOURCE/lpm_set_attributes.f90
- Timestamp:
- Apr 11, 2014 5:15:14 PM (10 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/lpm_set_attributes.f90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm_set_attributes.f90
r1321 r1359 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! New particle structure integrated. 23 ! Kind definition added to all floating point numbers. 23 24 ! 24 25 ! Former revisions: … … 77 78 78 79 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 80 82 81 83 USE pegrid … … 87 89 88 90 INTEGER(iwp) :: i !: 91 INTEGER(iwp) :: ip !: 89 92 INTEGER(iwp) :: j !: 93 INTEGER(iwp) :: jp !: 90 94 INTEGER(iwp) :: k !: 95 INTEGER(iwp) :: kp !: 91 96 INTEGER(iwp) :: n !: 97 INTEGER(iwp) :: nb !: 98 99 INTEGER(iwp), DIMENSION(0:7) :: start_index !: 100 INTEGER(iwp), DIMENSION(0:7) :: end_index !: 92 101 93 102 REAL(wp) :: aa !: … … 101 110 REAL(wp) :: pt_int_l !: 102 111 REAL(wp) :: pt_int_u !: 103 REAL(wp) :: u_int !:104 112 REAL(wp) :: u_int_l !: 105 113 REAL(wp) :: u_int_u !: 106 REAL(wp) :: v_int !:107 114 REAL(wp) :: v_int_l !: 108 115 REAL(wp) :: v_int_u !: … … 113 120 REAL(wp) :: y !: 114 121 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 115 128 CALL cpu_log( log_point_s(49), 'lpm_set_attributes', 'start' ) 116 129 … … 122 135 !-- Set particle color depending on the absolute value of the horizontal 123 136 !-- 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 197 251 ENDDO 198 252 … … 204 258 !-- (This is also done in flow_statistics, but flow_statistics is called 205 259 !-- after this routine.) 206 sums_l(:,4,0) = 0.0 260 sums_l(:,4,0) = 0.0_wp 207 261 DO i = nxl, nxr 208 262 DO j = nys, nyn … … 214 268 215 269 #if defined( __parallel ) 216 217 270 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 218 271 CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, & 219 272 MPI_REAL, MPI_SUM, comm2d, ierr ) 220 221 273 #else 222 223 274 sums(:,4) = sums_l(:,4,0) 224 225 275 #endif 226 227 276 sums(:,4) = sums(:,4) / ngp_2dh(0) 228 277 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 268 347 ENDDO 269 348 … … 272 351 !-- Set particle color depending on the height above the bottom 273 352 !-- 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 289 381 ENDDO 290 382 … … 295 387 IF ( particle_dvrpsize == 'absw' ) THEN 296 388 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 337 459 ENDDO 338 460
Note: See TracChangeset
for help on using the changeset viewer.
