Changeset 2696 for palm/trunk/SOURCE/plant_canopy_model_mod.f90
- Timestamp:
- Dec 14, 2017 5:12:51 PM (6 years ago)
- Location:
- palm/trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk
-
palm/trunk/SOURCE
-
palm/trunk/SOURCE/plant_canopy_model_mod.f90
r2669 r2696 1 1 !> @file plant_canopy_model_mod.f90 2 2 !------------------------------------------------------------------------------! 3 ! This file is part of PALM.3 ! This file is part of the PALM model system. 4 4 ! 5 5 ! PALM is free software: you can redistribute it and/or modify it under the … … 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Bugfix for vertical loop index pch_index in case of Netcdf input 28 ! Introduce 2D index array incorporate canopy top index 29 ! Check if canopy on top of topography do not exceed vertical dimension 30 ! Add check for canopy_mode in case of Netcdf input. 31 ! Enable _FillValue output for 3d quantities 32 ! Bugfix in reading of PIDS leaf area density (MS) 33 ! 34 ! 2669 2017-12-06 16:03:27Z raasch 27 35 ! coupling_char removed 28 36 ! … … 151 159 CHARACTER (LEN=20) :: canopy_mode = 'block' !< canopy coverage 152 160 153 INTEGER(iwp) :: pch_index = 0 !< plant canopy height/top index 154 INTEGER(iwp) :: & 155 lad_vertical_gradient_level_ind(10) = -9999 !< lad-profile levels (index) 161 INTEGER(iwp) :: pch_index = 0 !< plant canopy height/top index 162 INTEGER(iwp) :: lad_vertical_gradient_level_ind(10) = -9999 !< lad-profile levels (index) 163 164 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pch_index_ji !< local plant canopy top 156 165 157 166 LOGICAL :: calc_beta_lad_profile = .FALSE. !< switch for calc. of lad from beta func. 158 167 LOGICAL :: plant_canopy = .FALSE. !< switch for use of canopy model 159 160 REAL(wp) :: alpha_lad = 9999999.9_wp !< coefficient for lad calculation 161 REAL(wp) :: beta_lad = 9999999.9_wp !< coefficient for lad calculation 162 REAL(wp) :: canopy_drag_coeff = 0.0_wp !< canopy drag coefficient (parameter) 163 REAL(wp) :: cdc = 0.0_wp !< canopy drag coeff. (abbreviation used in equations) 164 REAL(wp) :: cthf = 0.0_wp !< canopy top heat flux 165 REAL(wp) :: dt_plant_canopy = 0.0_wp !< timestep account. for canopy drag 166 REAL(wp) :: ext_coef = 0.6_wp !< extinction coefficient 167 REAL(wp) :: lad_surface = 0.0_wp !< lad surface value 168 REAL(wp) :: lai_beta = 0.0_wp !< leaf area index (lai) for lad calc. 169 REAL(wp) :: & 170 leaf_scalar_exch_coeff = 0.0_wp !< canopy scalar exchange coeff. 171 REAL(wp) :: & 172 leaf_surface_conc = 0.0_wp !< leaf surface concentration 173 REAL(wp) :: lsec = 0.0_wp !< leaf scalar exchange coeff. 174 REAL(wp) :: lsc = 0.0_wp !< leaf surface concentration 175 176 REAL(wp) :: & 177 lad_vertical_gradient(10) = 0.0_wp !< lad gradient 178 REAL(wp) :: & 179 lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m) 168 LOGICAL :: usm_lad_rma = .TRUE. !< use MPI RMA to access LAD for raytracing (instead of global array) 169 170 REAL(wp) :: alpha_lad = 9999999.9_wp !< coefficient for lad calculation 171 REAL(wp) :: beta_lad = 9999999.9_wp !< coefficient for lad calculation 172 REAL(wp) :: canopy_drag_coeff = 0.0_wp !< canopy drag coefficient (parameter) 173 REAL(wp) :: cdc = 0.0_wp !< canopy drag coeff. (abbreviation used in equations) 174 REAL(wp) :: cthf = 0.0_wp !< canopy top heat flux 175 REAL(wp) :: dt_plant_canopy = 0.0_wp !< timestep account. for canopy drag 176 REAL(wp) :: ext_coef = 0.6_wp !< extinction coefficient 177 REAL(wp) :: lad_surface = 0.0_wp !< lad surface value 178 REAL(wp) :: lai_beta = 0.0_wp !< leaf area index (lai) for lad calc. 179 REAL(wp) :: leaf_scalar_exch_coeff = 0.0_wp !< canopy scalar exchange coeff. 180 REAL(wp) :: leaf_surface_conc = 0.0_wp !< leaf surface concentration 181 REAL(wp) :: lsec = 0.0_wp !< leaf scalar exchange coeff. 182 REAL(wp) :: lsc = 0.0_wp !< leaf surface concentration 183 REAL(wp) :: prototype_lad !< prototype leaf area density for computing effective optical depth 184 185 REAL(wp) :: lad_vertical_gradient(10) = 0.0_wp !< lad gradient 186 REAL(wp) :: lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m) 180 187 181 188 REAL(wp), DIMENSION(:), ALLOCATABLE :: lad !< leaf area density … … 201 208 !-- Public variables and constants 202 209 PUBLIC pc_heating_rate, canopy_mode, cthf, dt_plant_canopy, lad, lad_s, & 203 pch_index, plant_canopy 210 pch_index, plant_canopy, prototype_lad, usm_lad_rma 204 211 205 212 … … 288 295 289 296 USE control_parameters, & 290 ONLY: cloud_physics, message_string, microphysics_seifert 297 ONLY: cloud_physics, coupling_char, message_string, & 298 microphysics_seifert 299 300 USE netcdf_data_input_mod, & 301 ONLY: input_file_static, input_pids_static 291 302 292 303 … … 331 342 CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 ) 332 343 ENDIF 344 ! 345 !-- If dynamic input file is used, canopy need to be read from file 346 IF ( input_pids_static .AND. & 347 TRIM( canopy_mode ) /= 'read_from_file_3d' ) THEN 348 message_string = 'Usage of dynamic input file ' // & 349 TRIM( input_file_static ) // & 350 TRIM( coupling_char ) // ' requires ' // & 351 'canopy_mode = read_from_file_3d' 352 CALL message( 'check_parameters', 'PA0999', 1, 2, 0, 6, 0 ) 353 ENDIF 333 354 334 355 … … 342 363 !> Subroutine defining 3D output variables 343 364 !------------------------------------------------------------------------------! 344 SUBROUTINE pcm_data_output_3d( av, variable, found, local_pf )365 SUBROUTINE pcm_data_output_3d( av, variable, found, local_pf, fill_value ) 345 366 346 367 USE control_parameters, & … … 356 377 CHARACTER (LEN=*) :: variable !< 357 378 358 INTEGER(iwp) :: av !< 359 INTEGER(iwp) :: i !< 360 INTEGER(iwp) :: j !< 361 INTEGER(iwp) :: k !< 379 INTEGER(iwp) :: av !< 380 INTEGER(iwp) :: i !< 381 INTEGER(iwp) :: j !< 382 INTEGER(iwp) :: k !< 383 INTEGER(iwp) :: k_topo !< topography top index 362 384 363 385 LOGICAL :: found !< 364 386 387 REAL(wp) :: fill_value 365 388 REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb:nz_do3d) :: local_pf !< 366 389 … … 368 391 found = .TRUE. 369 392 393 local_pf = REAL( fill_value, KIND = 4 ) 370 394 371 395 SELECT CASE ( TRIM( variable ) ) … … 375 399 DO i = nxl, nxr 376 400 DO j = nys, nyn 377 DO k = nzb_s_inner(j,i), nz_do3d 378 local_pf(i,j,k) = pc_heating_rate(k-nzb_s_inner(j,i),j,i) 379 ENDDO 401 IF ( pch_index_ji(j,i) /= 0 ) THEN 402 k_topo = get_topography_top_index( j, i, 's' ) 403 DO k = k_topo, k_topo + pch_index_ji(j,i) 404 local_pf(i,j,k) = pc_heating_rate(k-k_topo,j,i) 405 ENDDO 406 ENDIF 380 407 ENDDO 381 408 ENDDO … … 387 414 DO i = nxl, nxr 388 415 DO j = nys, nyn 389 DO k = nzb_s_inner(j,i), nz_do3d 390 local_pf(i,j,k) = lad_s(k-nzb_s_inner(j,i),j,i) 391 ENDDO 416 IF ( pch_index_ji(j,i) /= 0 ) THEN 417 k_topo = get_topography_top_index( j, i, 's' ) 418 DO k = k_topo, k_topo + pch_index_ji(j,i) 419 local_pf(i,j,k) = lad_s(k-k_topo,j,i) 420 ENDDO 421 ENDIF 392 422 ENDDO 393 423 ENDDO … … 570 600 passive_scalar, urban_surface 571 601 602 USE netcdf_data_input_mod, & 603 ONLY: leaf_area_density_f 604 572 605 USE surface_mod, & 573 606 ONLY: surf_def_h, surf_lsm_h, surf_usm_h … … 716 749 CASE ( 'read_from_file_3d' ) 717 750 ! 751 !-- Initialize LAD with data from file. If LAD is given in NetCDF file, 752 !-- use these values, else take LAD profiles from ASCII file. 753 !-- Please note, in NetCDF file LAD is only given up to the maximum 754 !-- canopy top, indicated by leaf_area_density_f%nz. 755 lad_s = 0.0_wp 756 IF ( leaf_area_density_f%from_file ) THEN 757 ! 758 !-- Set also pch_index, used to be the upper bound of the vertical 759 !-- loops. Therefore, use the global top of the canopy layer. 760 pch_index = leaf_area_density_f%nz - 1 761 762 DO i = nxl, nxr 763 DO j = nys, nyn 764 DO k = 0, leaf_area_density_f%nz - 1 765 IF ( leaf_area_density_f%var(k,j,i) /= & 766 leaf_area_density_f%fill ) & 767 lad_s(k,j,i) = leaf_area_density_f%var(k,j,i) 768 ENDDO 769 ENDDO 770 ENDDO 771 CALL exchange_horiz( lad_s, nbgp ) 772 ! 773 ! ASCII file 718 774 !-- Initialize canopy parameters cdc (canopy drag coefficient), 719 775 !-- lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration) 720 776 !-- from file which contains complete 3D data (separate vertical profiles for 721 777 !-- each location). 722 CALL pcm_read_plant_canopy_3d 778 ELSE 779 CALL pcm_read_plant_canopy_3d 780 ENDIF 723 781 724 782 CASE DEFAULT … … 732 790 733 791 END SELECT 792 ! 793 !-- Initialize 2D index array indicating canopy top index. 794 ALLOCATE( pch_index_ji(nysg:nyng,nxlg:nxrg) ) 795 pch_index_ji = 0 796 797 DO i = nxl, nxr 798 DO j = nys, nyn 799 DO k = 0, pch_index 800 IF ( lad_s(k,j,i) /= 0 ) pch_index_ji(j,i) = k 801 ENDDO 802 ! 803 !-- Check whether topography and local vegetation on top exceed 804 !-- height of the model domain. 805 k = get_topography_top_index( j, i, 's' ) 806 IF ( k + pch_index_ji(j,i) >= nzt + 1 ) THEN 807 message_string = 'Local vegetation height on top of ' // & 808 'topography exceeds height of model domain.' 809 CALL message( 'pcm_init', 'PA0999', 2, 2, 0, 6, 0 ) 810 ENDIF 811 812 ENDDO 813 ENDDO 814 815 CALL exchange_horiz_2d_int( pch_index_ji, nys, nyn, nxl, nxr, nbgp ) 734 816 735 817 ! … … 757 839 DO i = nxlg, nxrg 758 840 DO j = nysg, nyng 759 DO k = pch_index -1, 0, -1760 IF ( k == pch_index -1 ) THEN841 DO k = pch_index_ji(j,i)-1, 0, -1 842 IF ( k == pch_index_ji(j,i)-1 ) THEN 761 843 cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) + & 762 844 ( 0.5_wp * lad_s(k+1,j,i) * & … … 815 897 DO i = nxlg, nxrg 816 898 DO j = nysg, nyng 817 DO k = 1, pch_index 899 DO k = 1, pch_index_ji(j,i) 818 900 IF ( cum_lai_hf(0,j,i) /= 0.0_wp ) THEN 819 901 pc_heating_rate(k,j,i) = cthf * & … … 945 1027 946 1028 CASE DEFAULT 947 write(message_string, '(a,i2,a)') &1029 WRITE(message_string, '(a,i2,a)') & 948 1030 'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"' 949 1031 CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 ) … … 1029 1111 !-- Determine topography-top index on u-grid 1030 1112 k_wall = get_topography_top_index( j, i, 'u' ) 1031 DO k = k_wall+1, k_wall +pch_index1113 DO k = k_wall+1, k_wall + pch_index_ji(j,i) 1032 1114 1033 1115 kk = k - k_wall !- lad arrays are defined flat … … 1094 1176 k_wall = get_topography_top_index( j, i, 'v' ) 1095 1177 1096 DO k = k_wall+1, k_wall +pch_index1178 DO k = k_wall+1, k_wall + pch_index_ji(j,i) 1097 1179 1098 1180 kk = k - k_wall !- lad arrays are defined flat … … 1159 1241 k_wall = get_topography_top_index( j, i, 'w' ) 1160 1242 1161 DO k = k_wall+1, k_wall +pch_index-11243 DO k = k_wall+1, k_wall + pch_index_ji(j,i) - 1 1162 1244 1163 1245 kk = k - k_wall !- lad arrays are defined flat … … 1211 1293 k_wall = get_topography_top_index( j, i, 's' ) 1212 1294 1213 DO k = k_wall+1, k_wall +pch_index1295 DO k = k_wall+1, k_wall + pch_index_ji(j,i) 1214 1296 1215 1297 kk = k - k_wall !- lad arrays are defined flat … … 1228 1310 k_wall = get_topography_top_index( j, i, 's' ) 1229 1311 1230 DO k = k_wall+1, k_wall +pch_index1312 DO k = k_wall+1, k_wall + pch_index_ji(j,i) 1231 1313 1232 1314 kk = k - k_wall !- lad arrays are defined flat … … 1258 1340 k_wall = get_topography_top_index( j, i, 's' ) 1259 1341 1260 DO k = k_wall+1, k_wall +pch_index1342 DO k = k_wall+1, k_wall + pch_index_ji(j,i) 1261 1343 1262 1344 kk = k - k_wall !- lad arrays are defined flat … … 1287 1369 k_wall = get_topography_top_index( j, i, 's' ) 1288 1370 1289 DO k = k_wall+1, k_wall +pch_index1371 DO k = k_wall+1, k_wall + pch_index_ji(j,i) 1290 1372 1291 1373 kk = k - k_wall !- lad arrays are defined flat … … 1370 1452 1371 1453 ddt_3d = 1.0_wp / dt_3d 1372 1373 1454 ! 1374 1455 !-- Compute drag for the three velocity components and the SGS-TKE … … 1381 1462 !-- Determine topography-top index on u-grid 1382 1463 k_wall = get_topography_top_index( j, i, 'u' ) 1383 1384 DO k = k_wall+1, k_wall+pch_index 1464 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) 1385 1465 1386 1466 kk = k - k_wall !- lad arrays are defined flat 1467 1387 1468 ! 1388 1469 !-- In order to create sharp boundaries of the plant canopy, … … 1442 1523 k_wall = get_topography_top_index( j, i, 'v' ) 1443 1524 1444 DO k = k_wall +1, k_wall+pch_index1525 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) 1445 1526 1446 1527 kk = k - k_wall !- lad arrays are defined flat … … 1502 1583 k_wall = get_topography_top_index( j, i, 'w' ) 1503 1584 1504 DO k = k_wall +1, k_wall+pch_index-11585 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) - 1 1505 1586 1506 1587 kk = k - k_wall !- lad arrays are defined flat … … 1549 1630 k_wall = get_topography_top_index( j, i, 's' ) 1550 1631 1551 DO k = k_wall +1, k_wall+pch_index1632 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) 1552 1633 kk = k - k_wall !- lad arrays are defined flat 1553 1634 tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) … … 1562 1643 k_wall = get_topography_top_index( j, i, 's' ) 1563 1644 1564 DO k = k_wall +1, k_wall+pch_index1645 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) 1565 1646 1566 1647 kk = k - k_wall … … 1588 1669 k_wall = get_topography_top_index( j, i, 's' ) 1589 1670 1590 DO k = k_wall +1, k_wall+pch_index1671 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) 1591 1672 1592 1673 kk = k - k_wall … … 1614 1695 k_wall = get_topography_top_index( j, i, 's' ) 1615 1696 1616 DO k = k_wall +1, k_wall+pch_index1697 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) 1617 1698 1618 1699 kk = k - k_wall
Note: See TracChangeset
for help on using the changeset viewer.
