- Timestamp:
- Dec 18, 2015 1:56:05 PM (9 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 4 edited
-
calc_mean_profile.f90 (modified) (2 diffs)
-
flow_statistics.f90 (modified) (14 diffs)
-
init_3d_model.f90 (modified) (10 diffs)
-
modules.f90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/calc_mean_profile.f90
r1683 r1738 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! bugfix: if a layer is completely filled with topography, no mean is calculated 22 22 ! 23 23 ! Former revisions: … … 126 126 #endif 127 127 128 hom(:,1,pr,0) = sums(:,pr) / ngp_2dh_s_inner(:,0) 128 DO k = nzb, nzt+1 129 IF ( ngp_2dh_s_inner(k,0) /= 0 ) THEN 130 hom(k,1,pr,0) = sums(k,pr) / ngp_2dh_s_inner(k,0) 131 ENDIF 132 ENDDO 129 133 130 134 ENDIF -
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 -
palm/trunk/SOURCE/init_3d_model.f90
r1735 r1738 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! calculate mean surface level height for each statistic region 22 22 ! 23 23 ! Former revisions: … … 273 273 274 274 USE statistics, & 275 ONLY: hom, hom_sum, pr_palm, rmask, spectrum_x, spectrum_y,&276 s tatistic_regions, sums, sums_divnew_l, sums_divold_l, sums_l,&277 sums_ l_l, sums_up_fraction_l, sums_wsts_bc_l, ts_value,&278 var_d, weight_pres, weight_substep275 ONLY: hom, hom_sum, mean_surface_level_height, pr_palm, rmask, & 276 spectrum_x, spectrum_y, statistic_regions, sums, sums_divnew_l, & 277 sums_divold_l, sums_l, sums_l_l, sums_up_fraction_l, & 278 sums_wsts_bc_l, ts_value, var_d, weight_pres, weight_substep 279 279 280 280 USE surface_layer_fluxes_mod, & … … 299 299 REAL(wp), DIMENSION(1:2) :: volume_flow_initial_l !< 300 300 301 REAL(wp), DIMENSION(:), ALLOCATABLE :: mean_surface_level_height_l !< 301 302 REAL(wp), DIMENSION(:), ALLOCATABLE :: ngp_3d_inner_l !< 302 303 REAL(wp), DIMENSION(:), ALLOCATABLE :: ngp_3d_inner_tmp !< … … 306 307 ! 307 308 !-- Allocate arrays 308 ALLOCATE( ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions), & 309 ALLOCATE( mean_surface_level_height(0:statistic_regions), & 310 mean_surface_level_height_l(0:statistic_regions), & 311 ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions), & 309 312 ngp_3d(0:statistic_regions), & 310 313 ngp_3d_inner(0:statistic_regions), & … … 1820 1823 ! 1821 1824 !-- Compute total sum of active mask grid points 1825 !-- and the mean surface level height for each statistic region 1822 1826 !-- ngp_2dh: number of grid points of a horizontal cross section through the 1823 1827 !-- total domain … … 1834 1838 ngp_sums = ( nz + 2 ) * ( pr_palm + max_pr_user ) 1835 1839 1840 mean_surface_level_height = 0.0_wp 1841 mean_surface_level_height_l = 0.0_wp 1842 1836 1843 DO sr = 0, statistic_regions 1837 1844 DO i = nxl, nxr … … 1841 1848 !-- All xy-grid points 1842 1849 ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1 1850 mean_surface_level_height_l(sr) = mean_surface_level_height_l(sr) + & 1851 zw(nzb_s_inner(j,i)) 1843 1852 ! 1844 1853 !-- xy-grid points above topography … … 1873 1882 MPI_SUM, comm2d, ierr ) 1874 1883 ngp_3d_inner = INT( ngp_3d_inner_tmp, KIND = SELECTED_INT_KIND( 18 ) ) 1884 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1885 CALL MPI_ALLREDUCE( mean_surface_level_height_l(0), & 1886 mean_surface_level_height(0), sr, MPI_REAL, & 1887 MPI_SUM, comm2d, ierr ) 1888 mean_surface_level_height = mean_surface_level_height / REAL( ngp_2dh ) 1875 1889 #else 1876 1890 ngp_2dh = ngp_2dh_l … … 1878 1892 ngp_2dh_s_inner = ngp_2dh_s_inner_l 1879 1893 ngp_3d_inner = INT( ngp_3d_inner_l, KIND = SELECTED_INT_KIND( 18 ) ) 1894 mean_surface_level_height = mean_surface_level_height_l / REAL( ngp_2dh_l ) 1880 1895 #endif 1881 1896 … … 1892 1907 ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) ) 1893 1908 1894 DEALLOCATE( ngp_2dh_l, ngp_2dh_outer_l, ngp_3d_inner_l, ngp_3d_inner_tmp ) 1909 DEALLOCATE( mean_surface_level_height_l, ngp_2dh_l, ngp_2dh_outer_l, & 1910 ngp_3d_inner_l, ngp_3d_inner_tmp ) 1895 1911 1896 1912 CALL location_message( 'leaving init_3d_model', .TRUE. ) -
palm/trunk/SOURCE/modules.f90
r1696 r1738 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! +mean_surface_level_height 22 22 ! 23 23 ! Former revisions: … … 1467 1467 LOGICAL :: flow_statistics_called = .FALSE. 1468 1468 REAL(wp) :: u_max, v_max, w_max 1469 REAL(wp), DIMENSION(:), ALLOCATABLE :: sums_divnew_l, sums_divold_l, & 1469 REAL(wp), DIMENSION(:), ALLOCATABLE :: mean_surface_level_height, & 1470 sums_divnew_l, sums_divold_l, & 1470 1471 var_d, weight_substep, & 1471 1472 weight_pres
Note: See TracChangeset
for help on using the changeset viewer.
