Changeset 1738 for palm/trunk/SOURCE/flow_statistics.f90
- Timestamp:
- Dec 18, 2015 1:56:05 PM (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/flow_statistics.f90
r1710 r1738 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! bugfixes for calculations in statistical regions which do not contain grid 22 ! points in the lowest vertical levels, mean surface level height considered 23 ! in the calculation of the characteristic vertical velocity, 24 ! old upstream parts removed 22 25 ! 23 26 ! Former revisions: … … 222 225 INTEGER(iwp) :: j !< 223 226 INTEGER(iwp) :: k !< 227 INTEGER(iwp) :: k_surface_level !< 224 228 INTEGER(iwp) :: nt !< 225 229 INTEGER(iwp) :: omp_get_thread_num !< … … 1247 1251 !-- Final values are obtained by division by the total number of grid points 1248 1252 !-- used for summation. After that store profiles. 1253 !-- Check, if statistical regions do contain at least one grid point at the 1254 !-- respective k-level, otherwise division by zero will lead to undefined 1255 !-- values, which may cause e.g. problems with NetCDF output 1249 1256 !-- Profiles: 1250 1257 DO k = nzb, nzt+1 1251 sums(k,3) = sums(k,3) / ngp_2dh(sr) 1252 sums(k,8:11) = sums(k,8:11) / ngp_2dh_s_inner(k,sr) 1253 sums(k,12:22) = sums(k,12:22) / ngp_2dh(sr) 1254 sums(k,23:29) = sums(k,23:29) / ngp_2dh_s_inner(k,sr) 1255 sums(k,30:32) = sums(k,30:32) / ngp_2dh(sr) 1256 sums(k,33:34) = sums(k,33:34) / ngp_2dh_s_inner(k,sr) 1257 sums(k,35:39) = sums(k,35:39) / ngp_2dh(sr) 1258 sums(k,40) = sums(k,40) / ngp_2dh_s_inner(k,sr) 1259 sums(k,45:53) = sums(k,45:53) / ngp_2dh(sr) 1260 sums(k,54) = sums(k,54) / ngp_2dh_s_inner(k,sr) 1261 sums(k,55:63) = sums(k,55:63) / ngp_2dh(sr) 1262 sums(k,64) = sums(k,64) / ngp_2dh_s_inner(k,sr) 1263 sums(k,65:69) = sums(k,65:69) / ngp_2dh(sr) 1264 sums(k,70:80) = sums(k,70:80) / ngp_2dh_s_inner(k,sr) 1265 sums(k,81:88) = sums(k,81:88) / ngp_2dh(sr) 1266 sums(k,89:105) = sums(k,89:105) / ngp_2dh(sr) 1267 sums(k,106:pr_palm-2) = sums(k,106:pr_palm-2)/ ngp_2dh_s_inner(k,sr) 1258 sums(k,3) = sums(k,3) / ngp_2dh(sr) 1259 sums(k,12:22) = sums(k,12:22) / ngp_2dh(sr) 1260 sums(k,30:32) = sums(k,30:32) / ngp_2dh(sr) 1261 sums(k,35:39) = sums(k,35:39) / ngp_2dh(sr) 1262 sums(k,45:53) = sums(k,45:53) / ngp_2dh(sr) 1263 sums(k,55:63) = sums(k,55:63) / ngp_2dh(sr) 1264 sums(k,81:88) = sums(k,81:88) / ngp_2dh(sr) 1265 sums(k,89:114) = sums(k,89:114) / ngp_2dh(sr) 1266 IF ( ngp_2dh_s_inner(k,sr) /= 0 ) THEN 1267 sums(k,8:11) = sums(k,8:11) / ngp_2dh_s_inner(k,sr) 1268 sums(k,23:29) = sums(k,23:29) / ngp_2dh_s_inner(k,sr) 1269 sums(k,33:34) = sums(k,33:34) / ngp_2dh_s_inner(k,sr) 1270 sums(k,40) = sums(k,40) / ngp_2dh_s_inner(k,sr) 1271 sums(k,54) = sums(k,54) / ngp_2dh_s_inner(k,sr) 1272 sums(k,64) = sums(k,64) / ngp_2dh_s_inner(k,sr) 1273 sums(k,70:80) = sums(k,70:80) / ngp_2dh_s_inner(k,sr) 1274 sums(k,115:pr_palm-2) = sums(k,115:pr_palm-2) / ngp_2dh_s_inner(k,sr) 1275 ENDIF 1268 1276 ENDDO 1269 1277 1270 !-- Upstream-parts1271 sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)1272 1278 !-- u* and so on 1273 1279 !-- As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose … … 1422 1428 hom(:,1,114,sr) = sums(:,114) !: L 1423 1429 1424 hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)1425 ! upstream-parts u_x, u_y, u_z, v_x,1426 ! v_y, usw. (in last but one profile)1427 1430 hom(:,1,pr_palm,sr) = sums(:,pr_palm) 1428 1431 ! u*, w'u', w'v', t* (in last profile) … … 1445 1448 IF ( ocean ) THEN 1446 1449 DO k = nzt, nzb+1, -1 1447 IF ( first .AND. hom(k,1,18,sr) < 0.0_wp & 1448 .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp) THEN 1450 IF ( first .AND. hom(k,1,18,sr) < -1.0E-8_wp ) THEN 1449 1451 first = .FALSE. 1450 1452 height = zw(k) 1451 1453 ENDIF 1452 IF ( hom(k,1,18,sr) < 0.0_wp .AND. & 1453 abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND. & 1454 IF ( hom(k,1,18,sr) < -1.0E-8_wp .AND. & 1454 1455 hom(k-1,1,18,sr) > hom(k,1,18,sr) ) THEN 1455 1456 IF ( zw(k) < 1.5_wp * height ) THEN … … 1463 1464 ELSE 1464 1465 DO k = nzb, nzt-1 1465 IF ( first .AND. hom(k,1,18,sr) < 0.0_wp & 1466 .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp ) THEN 1466 IF ( first .AND. hom(k,1,18,sr) < -1.0E-8_wp ) THEN 1467 1467 first = .FALSE. 1468 1468 height = zw(k) 1469 1469 ENDIF 1470 IF ( hom(k,1,18,sr) < 0.0_wp .AND. & 1471 abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND. & 1470 IF ( hom(k,1,18,sr) < -1.0E-8_wp .AND. & 1472 1471 hom(k+1,1,18,sr) > hom(k,1,18,sr) ) THEN 1473 1472 IF ( zw(k) < 1.5_wp * height ) THEN … … 1519 1518 1520 1519 ! 1520 !-- Determine vertical index which is nearest to the mean surface level 1521 !-- height of the respective statistic region 1522 DO k = nzb, nzt 1523 IF ( zw(k) >= mean_surface_level_height(sr) ) THEN 1524 k_surface_level = k 1525 EXIT 1526 ENDIF 1527 ENDDO 1528 ! 1521 1529 !-- Computation of both the characteristic vertical velocity and 1522 1530 !-- the characteristic convective boundary layer temperature. 1523 !-- The horizontal average at nzb+1 is input for the average temperature. 1524 IF ( hom(nzb,1,18,sr) > 0.0_wp .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8_wp & 1525 .AND. z_i(1) /= 0.0_wp ) THEN 1526 hom(nzb+8,1,pr_palm,sr) = ( g / hom(nzb+1,1,4,sr) * & 1527 hom(nzb,1,18,sr) * & 1528 ABS( z_i(1) ) )**0.333333333_wp 1529 !-- so far this only works if Prandtl layer is used 1530 hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr) 1531 !-- The inversion height entering into the equation is defined with respect 1532 !-- to the mean surface level height of the respective statistic region. 1533 !-- The horizontal average at surface level index + 1 is input for the 1534 !-- average temperature. 1535 IF ( hom(k_surface_level,1,18,sr) > 1.0E-8_wp .AND. z_i(1) /= 0.0_wp )& 1536 THEN 1537 hom(nzb+8,1,pr_palm,sr) = & 1538 ( g / hom(k_surface_level+1,1,4,sr) * hom(k_surface_level,1,18,sr)& 1539 * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp 1531 1540 ELSE 1532 1541 hom(nzb+8,1,pr_palm,sr) = 0.0_wp 1533 hom(nzb+11,1,pr_palm,sr) = 0.0_wp1534 1542 ENDIF 1535 1543 … … 1689 1697 INTEGER(iwp) :: j !< 1690 1698 INTEGER(iwp) :: k !< 1699 INTEGER(iwp) :: k_surface_level !< 1691 1700 INTEGER(iwp) :: nt !< 1692 1701 INTEGER(iwp) :: omp_get_thread_num !< … … 3230 3239 !-- Final values are obtained by division by the total number of grid points 3231 3240 !-- used for summation. After that store profiles. 3241 !-- Check, if statistical regions do contain at least one grid point at the 3242 !-- respective k-level, otherwise division by zero will lead to undefined 3243 !-- values, which may cause e.g. problems with NetCDF output 3232 3244 !-- Profiles: 3233 3245 DO k = nzb, nzt+1 3234 sums(k,3) = sums(k,3) / ngp_2dh(sr) 3235 sums(k,8:11) = sums(k,8:11) / ngp_2dh_s_inner(k,sr) 3236 sums(k,12:22) = sums(k,12:22) / ngp_2dh(sr) 3237 sums(k,23:29) = sums(k,23:29) / ngp_2dh_s_inner(k,sr) 3238 sums(k,30:32) = sums(k,30:32) / ngp_2dh(sr) 3239 sums(k,33:34) = sums(k,33:34) / ngp_2dh_s_inner(k,sr) 3240 sums(k,35:39) = sums(k,35:39) / ngp_2dh(sr) 3241 sums(k,40) = sums(k,40) / ngp_2dh_s_inner(k,sr) 3242 sums(k,45:53) = sums(k,45:53) / ngp_2dh(sr) 3243 sums(k,54) = sums(k,54) / ngp_2dh_s_inner(k,sr) 3244 sums(k,55:63) = sums(k,55:63) / ngp_2dh(sr) 3245 sums(k,64) = sums(k,64) / ngp_2dh_s_inner(k,sr) 3246 sums(k,70:80) = sums(k,70:80) / ngp_2dh_s_inner(k,sr) 3247 sums(k,81:88) = sums(k,81:88) / ngp_2dh(sr) 3248 sums(k,89:114) = sums(k,89:114) / ngp_2dh(sr) 3249 sums(k,115:pr_palm-2) = sums(k,115:pr_palm-2)/ ngp_2dh_s_inner(k,sr) 3250 3246 sums(k,3) = sums(k,3) / ngp_2dh(sr) 3247 sums(k,12:22) = sums(k,12:22) / ngp_2dh(sr) 3248 sums(k,30:32) = sums(k,30:32) / ngp_2dh(sr) 3249 sums(k,35:39) = sums(k,35:39) / ngp_2dh(sr) 3250 sums(k,45:53) = sums(k,45:53) / ngp_2dh(sr) 3251 sums(k,55:63) = sums(k,55:63) / ngp_2dh(sr) 3252 sums(k,81:88) = sums(k,81:88) / ngp_2dh(sr) 3253 sums(k,89:114) = sums(k,89:114) / ngp_2dh(sr) 3254 IF ( ngp_2dh_s_inner(k,sr) /= 0 ) THEN 3255 sums(k,8:11) = sums(k,8:11) / ngp_2dh_s_inner(k,sr) 3256 sums(k,23:29) = sums(k,23:29) / ngp_2dh_s_inner(k,sr) 3257 sums(k,33:34) = sums(k,33:34) / ngp_2dh_s_inner(k,sr) 3258 sums(k,40) = sums(k,40) / ngp_2dh_s_inner(k,sr) 3259 sums(k,54) = sums(k,54) / ngp_2dh_s_inner(k,sr) 3260 sums(k,64) = sums(k,64) / ngp_2dh_s_inner(k,sr) 3261 sums(k,70:80) = sums(k,70:80) / ngp_2dh_s_inner(k,sr) 3262 sums(k,115:pr_palm-2) = sums(k,115:pr_palm-2) / ngp_2dh_s_inner(k,sr) 3263 ENDIF 3251 3264 ENDDO 3252 3265 3253 !-- Upstream-parts3254 sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)3255 3266 !-- u* and so on 3256 3267 !-- As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose … … 3271 3282 IF ( max_pr_user > 0 ) THEN 3272 3283 DO k = nzb, nzt+1 3273 sums(k,pr_palm+1:pr_palm+max_pr_user) = & 3274 sums(k,pr_palm+1:pr_palm+max_pr_user) / & 3275 ngp_2dh_s_inner(k,sr) 3284 IF ( ngp_2dh_s_inner(k,sr) /= 0 ) THEN 3285 sums(k,pr_palm+1:pr_palm+max_pr_user) = & 3286 sums(k,pr_palm+1:pr_palm+max_pr_user) / & 3287 ngp_2dh_s_inner(k,sr) 3288 ENDIF 3276 3289 ENDDO 3277 3290 ENDIF … … 3363 3376 END IF 3364 3377 3365 hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)3366 ! upstream-parts u_x, u_y, u_z, v_x,3367 ! v_y, usw. (in last but one profile)3368 3378 hom(:,1,pr_palm,sr) = sums(:,pr_palm) 3369 3379 ! u*, w'u', w'v', t* (in last profile) … … 3386 3396 IF ( ocean ) THEN 3387 3397 DO k = nzt, nzb+1, -1 3388 IF ( first .AND. hom(k,1,18,sr) < 0.0_wp & 3389 .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp ) THEN 3398 IF ( first .AND. hom(k,1,18,sr) < -1.0E-8_wp ) THEN 3390 3399 first = .FALSE. 3391 3400 height = zw(k) 3392 3401 ENDIF 3393 IF ( hom(k,1,18,sr) < 0.0_wp .AND. & 3394 abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND. & 3402 IF ( hom(k,1,18,sr) < -1.0E-8_wp .AND. & 3395 3403 hom(k-1,1,18,sr) > hom(k,1,18,sr) ) THEN 3396 3404 IF ( zw(k) < 1.5_wp * height ) THEN … … 3404 3412 ELSE 3405 3413 DO k = nzb, nzt-1 3406 IF ( first .AND. hom(k,1,18,sr) < 0.0_wp & 3407 .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp ) THEN 3414 IF ( first .AND. hom(k,1,18,sr) < -1.0E-8_wp ) THEN 3408 3415 first = .FALSE. 3409 3416 height = zw(k) 3410 3417 ENDIF 3411 IF ( hom(k,1,18,sr) < 0.0 .AND. & 3412 abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND. & 3418 IF ( hom(k,1,18,sr) < -1.0E-8_wp .AND. & 3413 3419 hom(k+1,1,18,sr) > hom(k,1,18,sr) ) THEN 3414 3420 IF ( zw(k) < 1.5_wp * height ) THEN … … 3460 3466 3461 3467 ! 3468 !-- Determine vertical index which is nearest to the mean surface level 3469 !-- height of the respective statistic region 3470 DO k = nzb, nzt 3471 IF ( zw(k) >= mean_surface_level_height(sr) ) THEN 3472 k_surface_level = k 3473 EXIT 3474 ENDIF 3475 ENDDO 3476 3477 ! 3462 3478 !-- Computation of both the characteristic vertical velocity and 3463 3479 !-- the characteristic convective boundary layer temperature. 3464 !-- The horizontal average at nzb+1 is input for the average temperature.3465 IF ( hom(nzb,1,18,sr) > 0.0_wp .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8_wp & 3466 .AND. z_i(1) /= 0.0_wp ) THEN 3467 hom(nzb+8,1,pr_palm,sr) = ( g / hom(nzb+1,1,4,sr) * & 3468 hom(nzb,1,18,sr) * &3469 ABS( z_i(1) ) )**0.333333333_wp3470 !-- so far this only works if Prandtl layer is used 3471 hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)3480 !-- The inversion height entering into the equation is defined with respect 3481 !-- to the mean surface level height of the respective statistic region. 3482 !-- The horizontal average at surface level index + 1 is input for the 3483 !-- average temperature. 3484 IF ( hom(nzb,1,18,sr) > 1.0E-8_wp .AND. z_i(1) /= 0.0_wp ) THEN 3485 hom(nzb+8,1,pr_palm,sr) = & 3486 ( g / hom(k_surface_level+1,1,4,sr) * hom(k_surface_level,1,18,sr)& 3487 * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp 3472 3488 ELSE 3473 3489 hom(nzb+8,1,pr_palm,sr) = 0.0_wp 3474 hom(nzb+11,1,pr_palm,sr) = 0.0_wp3475 3490 ENDIF 3476 3491
Note: See TracChangeset
for help on using the changeset viewer.