Changeset 4770
- Timestamp:
- Nov 3, 2020 2:04:53 PM (3 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/plant_canopy_model_mod.f90 (modified) (29 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/plant_canopy_model_mod.f90
r4768 r4770 27 27 ! ----------------- 28 28 ! $Id$ 29 ! Consider basal-area density as an additional sink for momentum in the 30 ! prognostic equations for momentum and SGS TKE 31 ! 32 ! 4768 2020-11-02 19:11:23Z suehring 29 33 ! Enable 3D data output also with 64-bit precision 30 34 ! … … 293 297 REAL(wp), DIMENSION(:), ALLOCATABLE :: lad !< leaf area density 294 298 REAL(wp), DIMENSION(:), ALLOCATABLE :: pre_lad !< preliminary lad 295 299 300 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: bad_s !< basal-area density on scalar-grid 296 301 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: cum_lai_hf !< cumulative lai for heatflux calc. 297 302 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: lad_s !< lad on scalar-grid … … 546 551 unit = 'K s-1' 547 552 548 CASE ( 'pcm_ lad' )553 CASE ( 'pcm_bad', 'pcm_lad' ) 549 554 unit = 'm2 m-3' 550 555 … … 870 875 ENDIF 871 876 877 CASE ( 'pcm_bad' ) 878 IF ( av == 0 ) THEN 879 DO i = nxl, nxr 880 DO j = nys, nyn 881 DO k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do ) 882 local_pf(i,j,k) = bad_s(k,j,i) 883 ENDDO 884 ENDDO 885 ENDDO 886 ENDIF 887 872 888 CASE DEFAULT 873 889 found = .FALSE. … … 898 914 SELECT CASE ( TRIM( var ) ) 899 915 900 CASE ( 'pcm_heatrate', 'pcm_ lad', 'pcm_transpirationrate', 'pcm_latentrate')916 CASE ( 'pcm_heatrate', 'pcm_bad', 'pcm_lad', 'pcm_transpirationrate', 'pcm_latentrate') 901 917 grid_x = 'x' 902 918 grid_y = 'y' … … 1039 1055 1040 1056 LOGICAL :: lad_on_top = .FALSE. !< dummy flag to indicate that LAD is defined on a building roof 1057 LOGICAL :: bad_on_top = .FALSE. !< dummy flag to indicate that BAD is defined on a building roof 1041 1058 1042 1059 REAL(wp) :: canopy_height !< canopy height for lad-profile construction … … 1142 1159 1143 1160 ! 1144 !-- Allocate 3D-array for the leaf area density (lad_s). 1161 !-- Allocate 3D-array for the leaf-area density (lad_s) as well as 1162 !-- for basal-area densitiy (bad_s). Note, by default bad_s is zero. 1145 1163 ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1146 1164 ALLOCATE( bad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1165 bad_s = 0.0_wp 1147 1166 ! 1148 1167 !-- Initialization of the canopy coverage in the model domain: … … 1303 1322 #endif 1304 1323 IF ( lad_on_top ) THEN 1305 WRITE( message_string, * ) & 1306 'Resolved plant-canopy is ' // & 1307 'defined on top of an artificially '// & 1308 'created building grid point(s) ' // & 1309 '(emerged from the filtering) - ' // & 1310 'LAD profile is omitted at this / ' // & 1324 WRITE( message_string, * ) & 1325 'Resolved plant-canopy is defined on top of an ' // & 1326 'artificially created building grid point(s) '// & 1327 'the filtering) - LAD/BAD profile is omitted at this / ' //& 1311 1328 'these grid point(s).' 1312 1329 CALL message( 'pcm_init', 'PA0313', 0, 0, 0, 6, 0 ) … … 1321 1338 ELSE 1322 1339 CALL pcm_read_plant_canopy_3d 1340 ENDIF 1341 ! 1342 !-- Initialize LAD with data from file. If LAD is given in NetCDF file, 1343 !-- use these values, else take LAD profiles from ASCII file. 1344 !-- Please note, in NetCDF file LAD is only given up to the maximum 1345 !-- canopy top, indicated by basal_area_density_f%nz. 1346 bad_s = 0.0_wp 1347 IF ( basal_area_density_f%from_file ) THEN 1348 ! 1349 !-- Set also pch_index, used to be the upper bound of the vertical 1350 !-- loops. Therefore, use the global top of the canopy layer. 1351 pch_index = basal_area_density_f%nz - 1 1352 1353 DO i = nxl, nxr 1354 DO j = nys, nyn 1355 DO k = 0, basal_area_density_f%nz - 1 1356 IF ( basal_area_density_f%var(k,j,i) /= basal_area_density_f%fill ) & 1357 bad_s(k,j,i) = basal_area_density_f%var(k,j,i) 1358 ENDDO 1359 ! 1360 !-- Check if resolved vegetation is mapped onto buildings. 1361 !-- Please see comment for leaf_area density 1362 IF ( ANY( bad_s(:,j,i) /= 0.0_wp ) .AND. & 1363 ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) ) .AND. & 1364 ANY( BTEST( wall_flags_total_0(:,j,i), 4 ) ) ) THEN 1365 bad_s(:,j,i) = 0.0_wp 1366 bad_on_top = .TRUE. 1367 ENDIF 1368 ENDDO 1369 ENDDO 1370 #if defined( __parallel ) 1371 CALL MPI_ALLREDUCE( MPI_IN_PLACE, bad_on_top, 1, MPI_LOGICAL, & 1372 MPI_LOR, comm2d, ierr) 1373 #endif 1374 IF ( bad_on_top ) THEN 1375 WRITE( message_string, * ) & 1376 'Resolved plant-canopy is defined on top of an ' // & 1377 'artificially created building grid point(s) '// & 1378 'the filtering) - LAD/BAD profile is omitted at this / ' //& 1379 'these grid point(s).' 1380 CALL message( 'pcm_init', 'PA0313', 0, 0, 0, 6, 0 ) 1381 ENDIF 1382 CALL exchange_horiz( bad_s, nbgp ) 1323 1383 ENDIF 1324 1384 … … 1355 1415 DO j = nysg, nyng 1356 1416 DO k = 0, pch_index 1357 IF ( lad_s(k,j,i) /= 0 ) pch_index_ji(j,i) = k1417 IF ( lad_s(k,j,i) /= 0.0_wp .OR. bad_s(k,j,i) /= 0.0_wp ) pch_index_ji(j,i) = k 1358 1418 ENDDO 1359 1419 ! … … 1862 1922 LOGICAL :: building_edge_w !< control flag indicating a westward-facing building edge 1863 1923 1924 REAL(wp) :: bad_local !< local bad value 1864 1925 REAL(wp) :: ddt_3d !< inverse of the LES timestep (dt_3d) 1865 1926 REAL(wp) :: lad_local !< local lad value … … 1906 1967 !-- this is not applied, here the LAD is equals the LAD at grid 1907 1968 !-- point (k,j,i), in order to avoid that LAD is mistakenly mapped 1908 !-- on top of a roof where (usually) is no LAD is defined. 1969 !-- on top of a roof where (usually) is no LAD is defined. The same is 1970 !-- also valid for bad_s. 1909 1971 lad_local = lad_s(kk,j,i) 1910 IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp&1972 IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp & 1911 1973 .AND. .NOT. building_edge_w ) lad_local = lad_s(kk,j,i-1) 1974 1975 bad_local = bad_s(kk,j,i) 1976 IF ( bad_local == 0.0_wp .AND. bad_s(kk,j,i-1) > 0.0_wp & 1977 .AND. .NOT. building_edge_w ) bad_local = bad_s(kk,j,i-1) 1912 1978 ! 1913 1979 !-- In order to avoid that LAD is mistakenly considered at right- … … 1915 1981 !-- u-component at index j,i is still on the building while the 1916 1982 !-- topography top for the scalar isn't), LAD is taken from grid 1917 !-- point (j,i-1). 1918 IF ( lad_local > 0.0_wp .AND. lad_s(kk,j,i-1) == 0.0_wp&1983 !-- point (j,i-1). The same is also valid for bad_s. 1984 IF ( lad_local > 0.0_wp .AND. lad_s(kk,j,i-1) == 0.0_wp & 1919 1985 .AND. building_edge_e ) lad_local = lad_s(kk,j,i-1) 1986 IF ( bad_local > 0.0_wp .AND. bad_s(kk,j,i-1) == 0.0_wp & 1987 .AND. building_edge_e ) bad_local = bad_s(kk,j,i-1) 1920 1988 1921 1989 pre_tend = 0.0_wp … … 1924 1992 !-- Calculate preliminary value (pre_tend) of the tendency 1925 1993 pre_tend = - canopy_drag_coeff * & 1926 lad_local *&1994 ( lad_local + bad_local ) * & 1927 1995 SQRT( u(k,j,i)**2 + & 1928 1996 ( 0.25_wp * ( v(k,j,i-1) + & … … 1987 2055 !-- point (k,j,i), in order to avoid that LAD is mistakenly mapped 1988 2056 !-- on top of a roof where (usually) is no LAD is defined. 2057 !-- The same is also valid for bad_s. 1989 2058 lad_local = lad_s(kk,j,i) 1990 IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp&2059 IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp & 1991 2060 .AND. .NOT. building_edge_s ) lad_local = lad_s(kk,j-1,i) 2061 2062 bad_local = bad_s(kk,j,i) 2063 IF ( bad_local == 0.0_wp .AND. bad_s(kk,j-1,i) > 0.0_wp & 2064 .AND. .NOT. building_edge_s ) bad_local = bad_s(kk,j-1,i) 1992 2065 ! 1993 2066 !-- In order to avoid that LAD is mistakenly considered at right- … … 1995 2068 !-- u-component at index j,i is still on the building while the 1996 2069 !-- topography top for the scalar isn't), LAD is taken from grid 1997 !-- point (j,i-1). 1998 IF ( lad_local > 0.0_wp .AND. lad_s(kk,j-1,i) == 0.0_wp&2070 !-- point (j,i-1). The same is also valid for bad_s. 2071 IF ( lad_local > 0.0_wp .AND. lad_s(kk,j-1,i) == 0.0_wp & 1999 2072 .AND. building_edge_n ) lad_local = lad_s(kk,j-1,i) 2073 2074 IF ( bad_local > 0.0_wp .AND. bad_s(kk,j-1,i) == 0.0_wp & 2075 .AND. building_edge_n ) bad_local = bad_s(kk,j-1,i) 2000 2076 2001 2077 pre_tend = 0.0_wp … … 2004 2080 !-- Calculate preliminary value (pre_tend) of the tendency 2005 2081 pre_tend = - canopy_drag_coeff * & 2006 lad_local *&2082 ( lad_local + bad_local ) * & 2007 2083 SQRT( ( 0.25_wp * ( u(k,j-1,i) + & 2008 2084 u(k,j-1,i+1) + & … … 2054 2130 !-- Calculate preliminary value (pre_tend) of the tendency 2055 2131 pre_tend = - canopy_drag_coeff * & 2056 (0.5_wp * & 2057 ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) * & 2132 ( 0.5_wp * & 2133 ( lad_s(kk+1,j,i) + lad_s(kk,j,i) ) + & 2134 0.5_wp * & 2135 ( bad_s(kk+1,j,i) + bad_s(kk,j,i) ) ) * & 2058 2136 SQRT( ( 0.25_wp * ( u(k,j,i) + & 2059 2137 u(k,j,i+1) + & … … 2161 2239 tend(k,j,i) = tend(k,j,i) - & 2162 2240 2.0_wp * canopy_drag_coeff * & 2163 lad_s(kk,j,i) *&2241 ( lad_s(kk,j,i) + bad_s(kk,j,i) ) * & 2164 2242 SQRT( ( 0.5_wp * ( u(k,j,i) + & 2165 2243 u(k,j,i+1) ) & … … 2253 2331 LOGICAL :: building_edge_w !< control flag indicating a westward-facing building edge 2254 2332 2333 REAL(wp) :: bad_local !< local lad value 2255 2334 REAL(wp) :: ddt_3d !< inverse of the LES timestep (dt_3d) 2256 2335 REAL(wp) :: lad_local !< local lad value … … 2296 2375 !-- point (k,j,i), in order to avoid that LAD is mistakenly mapped 2297 2376 !-- on top of a roof where (usually) is no LAD is defined. 2377 !-- The same is also valid for bad_s. 2298 2378 lad_local = lad_s(kk,j,i) 2299 IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp .AND.&2379 IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp .AND. & 2300 2380 .NOT. building_edge_w ) lad_local = lad_s(kk,j,i-1) 2381 2382 bad_local = bad_s(kk,j,i) 2383 IF ( bad_local == 0.0_wp .AND. bad_s(kk,j,i-1) > 0.0_wp .AND. & 2384 .NOT. building_edge_w ) bad_local = bad_s(kk,j,i-1) 2301 2385 ! 2302 2386 !-- In order to avoid that LAD is mistakenly considered at right- … … 2304 2388 !-- u-component at index j,i is still on the building while the 2305 2389 !-- topography top for the scalar isn't), LAD is taken from grid 2306 !-- point (j,i-1). 2307 IF ( lad_local > 0.0_wp .AND. lad_s(kk,j,i-1) == 0.0_wp .AND.&2390 !-- point (j,i-1). The same is also valid for bad_s. 2391 IF ( lad_local > 0.0_wp .AND. lad_s(kk,j,i-1) == 0.0_wp .AND. & 2308 2392 building_edge_e ) lad_local = lad_s(kk,j,i-1) 2393 IF ( bad_local > 0.0_wp .AND. bad_s(kk,j,i-1) == 0.0_wp .AND. & 2394 building_edge_e ) bad_local = bad_s(kk,j,i-1) 2309 2395 2310 2396 pre_tend = 0.0_wp … … 2313 2399 !-- Calculate preliminary value (pre_tend) of the tendency 2314 2400 pre_tend = - canopy_drag_coeff * & 2315 lad_local * &2401 ( lad_local + bad_local ) * & 2316 2402 SQRT( u(k,j,i)**2 + & 2317 2403 ( 0.25_wp * ( v(k,j,i-1) + & … … 2374 2460 !-- point (k,j,i), in order to avoid that LAD is mistakenly mapped 2375 2461 !-- on top of a roof where (usually) is no LAD is defined. 2462 !-- The same is also valid for bad_s. 2376 2463 lad_local = lad_s(kk,j,i) 2377 IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp .AND.&2464 IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp .AND. & 2378 2465 .NOT. building_edge_s ) lad_local = lad_s(kk,j-1,i) 2466 2467 bad_local = bad_s(kk,j,i) 2468 IF ( bad_local == 0.0_wp .AND. bad_s(kk,j-1,i) > 0.0_wp .AND. & 2469 .NOT. building_edge_s ) bad_local = bad_s(kk,j-1,i) 2379 2470 ! 2380 2471 !-- In order to avoid that LAD is mistakenly considered at right- … … 2382 2473 !-- u-component at index j,i is still on the building while the 2383 2474 !-- topography top for the scalar isn't), LAD is taken from grid 2384 !-- point (j,i-1). 2385 IF ( lad_local > 0.0_wp .AND. lad_s(kk,j-1,i) == 0.0_wp .AND.&2475 !-- point (j,i-1). The same is also valid for bad_s. 2476 IF ( lad_local > 0.0_wp .AND. lad_s(kk,j-1,i) == 0.0_wp .AND. & 2386 2477 building_edge_n ) lad_local = lad_s(kk,j-1,i) 2478 IF ( bad_local > 0.0_wp .AND. bad_s(kk,j-1,i) == 0.0_wp .AND. & 2479 building_edge_n ) bad_local = bad_s(kk,j-1,i) 2387 2480 2388 2481 pre_tend = 0.0_wp … … 2391 2484 !-- Calculate preliminary value (pre_tend) of the tendency 2392 2485 pre_tend = - canopy_drag_coeff * & 2393 lad_local *&2486 ( lad_local + bad_local ) * & 2394 2487 SQRT( ( 0.25_wp * ( u(k,j-1,i) + & 2395 2488 u(k,j-1,i+1) + & … … 2437 2530 !-- Calculate preliminary value (pre_tend) of the tendency 2438 2531 pre_tend = - canopy_drag_coeff * & 2439 (0.5_wp * & 2440 ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) * & 2441 SQRT( ( 0.25_wp * ( u(k,j,i) + & 2532 ( 0.5_wp * & 2533 ( lad_s(kk+1,j,i) + lad_s(kk,j,i) ) + & 2534 0.5_wp * & 2535 ( bad_s(kk+1,j,i) + bad_s(kk,j,i) ) ) * & 2536 SQRT( ( 0.25_wp * ( u(k,j,i) + & 2442 2537 u(k,j,i+1) + & 2443 2538 u(k+1,j,i) + & … … 2527 2622 tend(k,j,i) = tend(k,j,i) - & 2528 2623 2.0_wp * canopy_drag_coeff * & 2529 lad_s(kk,j,i) *&2624 ( lad_s(kk,j,i) + bad_s(kk,j,i) ) * & 2530 2625 SQRT( ( 0.5_wp * ( u(k,j,i) + & 2531 2626 u(k,j,i+1) ) & … … 2540 2635 e(k,j,i) 2541 2636 ENDDO 2542 2543 2637 ! 2544 2638 !-- scalar concentration
Note: See TracChangeset
for help on using the changeset viewer.
