Changeset 1484
- Timestamp:
- Oct 21, 2014 10:53:05 AM (10 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 13 edited
-
Makefile (modified) (7 diffs)
-
check_parameters.f90 (modified) (7 diffs)
-
header.f90 (modified) (8 diffs)
-
init_3d_model.f90 (modified) (4 diffs)
-
modules.f90 (modified) (12 diffs)
-
package_parin.f90 (modified) (9 diffs)
-
parin.f90 (modified) (11 diffs)
-
plant_canopy_model.f90 (modified) (14 diffs)
-
prognostic_equations.f90 (modified) (3 diffs)
-
read_var_list.f90 (modified) (10 diffs)
-
timestep.f90 (modified) (7 diffs)
-
user_init_plant_canopy.f90 (modified) (5 diffs)
-
write_var_list.f90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/Makefile
r1445 r1484 20 20 # Current revisions: 21 21 # ------------------ 22 # 22 # plant_canopy_model-dependency added for check_parameters, header, init_3d_model, 23 # package_parin, read_var_list, user_init_plant_canopy, write_var_list 23 24 # 24 25 # Former revisions: … … 288 289 check_for_restart.o: modules.o mod_kinds.o 289 290 check_open.o: modules.o mod_kinds.o mod_particle_attributes.o 290 check_parameters.o: modules.o mod_kinds.o subsidence.o 291 check_parameters.o: modules.o mod_kinds.o subsidence.o plant_canopy_model.o 291 292 close_file.o: modules.o mod_kinds.o 292 293 compute_vpt.o: modules.o mod_kinds.o … … 318 319 flow_statistics.o: modules.o cpulog.o mod_kinds.o 319 320 global_min_max.o: modules.o mod_kinds.o 320 header.o: modules.o cpulog.o mod_kinds.o subsidence.o321 header.o: modules.o cpulog.o mod_kinds.o plant_canopy_model.o subsidence.o 321 322 impact_of_latent_heat.o: modules.o mod_kinds.o 322 323 inflow_turbulence.o: modules.o cpulog.o mod_kinds.o 323 324 init_1d_model.o: modules.o mod_kinds.o 324 325 init_3d_model.o: modules.o cpulog.o mod_kinds.o random_function.o advec_ws.o \ 325 ls_forcing.o lpm_init.o random_generator_parallel.o326 ls_forcing.o lpm_init.o plant_canopy_model.o random_generator_parallel.o 326 327 init_advec.o: modules.o mod_kinds.o 327 328 init_cloud_physics.o: modules.o mod_kinds.o … … 379 380 netcdf.o: modules.o mod_kinds.o 380 381 nudging.o: modules.o cpulog.o mod_kinds.o 381 package_parin.o: modules.o mod_kinds.o 382 package_parin.o: modules.o mod_kinds.o plant_canopy_model.o 382 383 palm.o: modules.o cpulog.o ls_forcing.o mod_kinds.o nudging.o 383 384 parin.o: modules.o cpulog.o mod_kinds.o progress_bar.o … … 400 401 random_generator_parallel.o: mod_kinds.o 401 402 read_3d_binary.o: modules.o cpulog.o mod_kinds.o random_function.o random_generator_parallel.o 402 read_var_list.o: modules.o mod_kinds.o 403 read_var_list.o: modules.o mod_kinds.o plant_canopy_model.o 403 404 run_control.o: modules.o cpulog.o mod_kinds.o 404 405 set_slicer_attributes_dvrp.o: modules.o mod_kinds.o … … 434 435 user_init_3d_model.o: modules.o mod_kinds.o user_module.o 435 436 user_init_grid.o: modules.o mod_kinds.o user_module.o 436 user_init_plant_canopy.o: modules.o mod_kinds.o user_module.o 437 user_init_plant_canopy.o: modules.o mod_kinds.o user_module.o plant_canopy_model.o 437 438 user_last_actions.o: modules.o mod_kinds.o user_module.o 438 439 user_lpm_advec.o: modules.o mod_kinds.o user_module.o … … 446 447 wall_fluxes.o: modules.o mod_kinds.o 447 448 write_3d_binary.o: modules.o cpulog.o mod_kinds.o random_function.o random_generator_parallel.o 448 write_var_list.o: modules.o mod_kinds.o 449 write_var_list.o: modules.o mod_kinds.o plant_canopy_model.o -
palm/trunk/SOURCE/check_parameters.f90
r1456 r1484 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Changes due to new module structure of the plant canopy model: 23 ! module plant_canopy_model_mod added, 24 ! checks regarding usage of new method for leaf area density profile 25 ! construction added, 26 ! lad-profile construction moved to new subroutine init_plant_canopy within 27 ! the module plant_canopy_model_mod, 28 ! drag_coefficient renamed to canopy_drag_coeff. 29 ! Missing KIND-attribute for REAL constant added 23 30 ! 24 31 ! Former revisions: … … 246 253 USE particle_attributes 247 254 USE pegrid 255 USE plant_canopy_model_mod 248 256 USE profil_parameter 249 257 USE spectrum … … 891 899 ENDIF 892 900 893 IF ( plant_canopy .AND. ( drag_coefficient == 0.0_wp ) ) THEN 894 message_string = 'plant_canopy = .TRUE. requires a non-zero drag ' // & 895 'coefficient & given value is drag_coefficient = 0.0' 896 CALL message( 'check_parameters', 'PA0041', 1, 2, 0, 6, 0 ) 897 ENDIF 898 899 IF ( plant_canopy .AND. cloud_physics .AND. icloud_scheme == 0 ) THEN 900 message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' // & 901 ' seifert_beheng' 902 CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 ) 903 ENDIF 901 IF ( plant_canopy ) THEN 902 903 IF ( canopy_drag_coeff == 0.0_wp ) THEN 904 message_string = 'plant_canopy = .TRUE. requires a non-zero drag '// & 905 'coefficient & given value is canopy_drag_coeff = 0.0' 906 CALL message( 'check_parameters', 'PA0041', 1, 2, 0, 6, 0 ) 907 ENDIF 908 909 IF ( ( alpha_lad /= 9999999.9_wp .AND. beta_lad == 9999999.9_wp ) .OR.& 910 beta_lad /= 9999999.9_wp .AND. alpha_lad == 9999999.9_wp ) THEN 911 message_string = 'using the beta function for the construction ' // & 912 'of the leaf area density profile requires ' // & 913 'both alpha_lad and beta_lad to be /= 9999999.9' 914 CALL message( 'check_parameters', 'PA0118', 1, 2, 0, 6, 0 ) 915 ENDIF 916 917 IF ( calc_beta_lad_profile .AND. lai_beta == 0.0_wp ) THEN 918 message_string = 'using the beta function for the construction ' // & 919 'of the leaf area density profile requires ' // & 920 'a non-zero lai_beta, but given value is ' // & 921 'lai_beta = 0.0' 922 CALL message( 'check_parameters', 'PA0119', 1, 2, 0, 6, 0 ) 923 ENDIF 924 925 IF ( calc_beta_lad_profile .AND. lad_surface /= 0.0_wp ) THEN 926 message_string = 'simultaneous setting of alpha_lad /= 9999999.9' // & 927 'and lad_surface /= 0.0 is not possible, ' // & 928 'use either vertical gradients or the beta ' // & 929 'function for the construction of the leaf area '// & 930 'density profile' 931 CALL message( 'check_parameters', 'PA0120', 1, 2, 0, 6, 0 ) 932 ENDIF 933 934 IF ( cloud_physics .AND. icloud_scheme == 0 ) THEN 935 message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' // & 936 ' seifert_beheng' 937 CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 ) 938 ENDIF 939 940 ENDIF 904 941 905 942 IF ( .NOT. ( loop_optimization == 'cache' .OR. & … … 924 961 IF ( ocean ) sa_init = sa_surface 925 962 IF ( passive_scalar ) q_init = s_surface 926 IF ( plant_canopy ) lad = 0.0_wp927 963 928 964 ! … … 1264 1300 ENDIF 1265 1301 1266 !1267 !-- If required compute the profile of leaf area density used in the plant1268 !-- canopy model1269 IF ( plant_canopy ) THEN1270 1271 i = 11272 gradient = 0.0_wp1273 1274 IF ( .NOT. ocean ) THEN1275 1276 lad(0) = lad_surface1277 1278 lad_vertical_gradient_level_ind(1) = 01279 DO k = 1, pch_index1280 IF ( i < 11 ) THEN1281 IF ( lad_vertical_gradient_level(i) < zu(k) .AND. &1282 lad_vertical_gradient_level(i) >= 0.0_wp ) THEN1283 gradient = lad_vertical_gradient(i)1284 lad_vertical_gradient_level_ind(i) = k - 11285 i = i + 11286 ENDIF1287 ENDIF1288 IF ( gradient /= 0.0_wp ) THEN1289 IF ( k /= 1 ) THEN1290 lad(k) = lad(k-1) + dzu(k) * gradient1291 ELSE1292 lad(k) = lad_surface + dzu(k) *gradient1293 ENDIF1294 ELSE1295 lad(k) = lad(k-1)1296 ENDIF1297 ENDDO1298 1299 ENDIF1300 1301 !1302 !-- In case of no given leaf area density gradients, choose a vanishing1303 !-- gradient1304 IF ( lad_vertical_gradient_level(1) == -9999999.9_wp ) THEN1305 lad_vertical_gradient_level(1) = 0.0_wp1306 ENDIF1307 1308 ENDIF1309 1302 1310 1303 ENDIF … … 1631 1624 ENDIF 1632 1625 1633 IF ( top_heatflux /= 0.0 .AND. top_heatflux /= 9999999.9_wp ) &1626 IF ( top_heatflux /= 0.0_wp .AND. top_heatflux /= 9999999.9_wp ) & 1634 1627 THEN 1635 1628 message_string = 'heatflux must not be set for pure neutral flow' … … 1828 1821 rayleigh_damping_factor = 0.0_wp 1829 1822 ELSE 1830 IF ( rayleigh_damping_factor < 0.0 .OR. rayleigh_damping_factor > 1.0_wp ) &1823 IF ( rayleigh_damping_factor < 0.0_wp .OR. rayleigh_damping_factor > 1.0_wp ) & 1831 1824 THEN 1832 1825 WRITE( message_string, * ) 'rayleigh_damping_factor = ', & -
palm/trunk/SOURCE/header.f90
r1483 r1484 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Changes due to new module structure of the plant canopy model: 23 ! module plant_canopy_model_mod and output for new canopy model parameters 24 ! (alpha_lad, beta_lad, lai_beta,...) added, 25 ! drag_coefficient, leaf_surface_concentration and scalar_exchange_coefficient 26 ! renamed to canopy_drag_coeff, leaf_surface_conc and leaf_scalar_exch_coeff, 27 ! learde renamed leaf_area_density. 28 ! Bugfix: DO-WHILE-loop for lad header information additionally restricted 29 ! by maximum number of gradient levels (currently 10) 23 30 ! 24 31 ! Former revisions: … … 165 172 166 173 USE arrays_3d, & 167 ONLY: lad,pt_init, qsws, q_init, sa_init, shf, ug, vg, w_subs, zu174 ONLY: pt_init, qsws, q_init, sa_init, shf, ug, vg, w_subs, zu 168 175 169 176 USE control_parameters … … 205 212 206 213 USE pegrid 214 215 USE plant_canopy_model_mod, & 216 ONLY: alpha_lad, beta_lad, calc_beta_lad_profile, canopy_drag_coeff, & 217 canopy_mode, cthf, lad, lad_surface, lad_vertical_gradient, & 218 lad_vertical_gradient_level, lad_vertical_gradient_level_ind, & 219 lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, pch_index, & 220 plant_canopy 207 221 208 222 USE spectrum, & … … 242 256 CHARACTER (LEN=86) :: coordinates !: 243 257 CHARACTER (LEN=86) :: gradients !: 244 CHARACTER (LEN=86) :: lea rde!:258 CHARACTER (LEN=86) :: leaf_area_density !: 245 259 CHARACTER (LEN=86) :: slices !: 246 260 CHARACTER (LEN=86) :: temperatures !: … … 270 284 INTEGER(iwp) :: io !: 271 285 INTEGER(iwp) :: j !: 286 INTEGER(iwp) :: k !: 272 287 INTEGER(iwp) :: l !: 273 288 INTEGER(iwp) :: ll !: 274 289 INTEGER(iwp) :: mpi_type !: 275 290 291 REAL(wp) :: canopy_height !: canopy height (in m) 276 292 REAL(wp) :: cpuseconds_per_simulated_second !: 277 293 … … 760 776 761 777 IF ( plant_canopy ) THEN 762 763 WRITE ( io, 280 ) canopy_mode, pch_index, drag_coefficient 778 779 canopy_height = pch_index * dz 780 781 WRITE ( io, 280 ) canopy_mode, canopy_height, pch_index, & 782 canopy_drag_coeff 764 783 IF ( passive_scalar ) THEN 765 WRITE ( io, 281 ) scalar_exchange_coefficient,&766 leaf_surface_concentration784 WRITE ( io, 281 ) leaf_scalar_exch_coeff, & 785 leaf_surface_conc 767 786 ENDIF 768 787 769 788 ! 770 789 !-- Heat flux at the top of vegetation 771 WRITE ( io, 282 ) cthf 772 773 ! 774 !-- Leaf area density profile 775 !-- Building output strings, starting with surface value 776 WRITE ( learde, '(F6.4)' ) lad_surface 777 gradients = '------' 778 slices = ' 0' 779 coordinates = ' 0.0' 780 i = 1 781 DO WHILE ( lad_vertical_gradient_level_ind(i) /= -9999 ) 782 783 WRITE (coor_chr,'(F7.2)') lad(lad_vertical_gradient_level_ind(i)) 784 learde = TRIM( learde ) // ' ' // TRIM( coor_chr ) 785 786 WRITE (coor_chr,'(F7.2)') lad_vertical_gradient(i) 787 gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr ) 788 789 WRITE (coor_chr,'(I7)') lad_vertical_gradient_level_ind(i) 790 slices = TRIM( slices ) // ' ' // TRIM( coor_chr ) 791 792 WRITE (coor_chr,'(F7.1)') lad_vertical_gradient_level(i) 793 coordinates = TRIM( coordinates ) // ' ' // TRIM( coor_chr ) 794 795 i = i + 1 796 ENDDO 797 798 WRITE ( io, 283 ) TRIM( coordinates ), TRIM( learde ), & 799 TRIM( gradients ), TRIM( slices ) 800 801 ENDIF 790 WRITE ( io, 282 ) cthf 791 792 ! 793 !-- Leaf area density profile, calculated either from given vertical 794 !-- gradients or from beta probability density function. 795 IF ( .NOT. calc_beta_lad_profile ) THEN 796 797 !-- Building output strings, starting with surface value 798 WRITE ( leaf_area_density, '(F6.4)' ) lad_surface 799 gradients = '------' 800 slices = ' 0' 801 coordinates = ' 0.0' 802 i = 1 803 DO WHILE ( i < 11 .AND. lad_vertical_gradient_level_ind(i) /= -9999 ) 804 805 WRITE (coor_chr,'(F7.2)') lad(lad_vertical_gradient_level_ind(i)) 806 leaf_area_density = TRIM( leaf_area_density ) // ' ' // TRIM( coor_chr ) 807 808 WRITE (coor_chr,'(F7.2)') lad_vertical_gradient(i) 809 gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr ) 810 811 WRITE (coor_chr,'(I7)') lad_vertical_gradient_level_ind(i) 812 slices = TRIM( slices ) // ' ' // TRIM( coor_chr ) 813 814 WRITE (coor_chr,'(F7.1)') lad_vertical_gradient_level(i) 815 coordinates = TRIM( coordinates ) // ' ' // TRIM( coor_chr ) 816 817 i = i + 1 818 ENDDO 819 820 WRITE ( io, 283 ) TRIM( coordinates ), TRIM( leaf_area_density ), & 821 TRIM( gradients ), TRIM( slices ) 822 823 ELSE 824 825 WRITE ( leaf_area_density, '(F6.4)' ) lad_surface 826 coordinates = ' 0.0' 827 828 DO k = 1, pch_index 829 830 WRITE (coor_chr,'(F7.2)') lad(k) 831 leaf_area_density = TRIM( leaf_area_density ) // ' ' // TRIM( coor_chr ) 832 833 WRITE (coor_chr,'(F7.1)') zu(k) 834 coordinates = TRIM( coordinates ) // ' ' // TRIM( coor_chr ) 835 836 ENDDO 837 838 WRITE ( io, 284 ) TRIM( coordinates ), TRIM( leaf_area_density ), alpha_lad, & 839 beta_lad, lai_beta 840 841 ENDIF 842 843 ENDIF 844 802 845 803 846 ! … … 1836 1879 ' ------------------------------'// & 1837 1880 ' Canopy mode: ', A / & 1838 ' Canopy top: ',I4/ &1881 ' Canopy height: ',F6.2,'m (',I4,' grid points)' / & 1839 1882 ' Leaf drag coefficient: ',F6.2 /) 1840 281 FORMAT (/ ' Scalar _exchange_coefficient: ',F6.2 / &1883 281 FORMAT (/ ' Scalar exchange coefficient: ',F6.2 / & 1841 1884 ' Scalar concentration at leaf surfaces in kg/m**3: ',F6.2 /) 1842 1885 282 FORMAT (' Predefined constant heatflux at the top of the vegetation: ',F6.2,' K m/s') … … 1846 1889 ' Gradient: ',A,' m**2/m**4'/ & 1847 1890 ' Gridpoint: ',A) 1891 284 FORMAT (//' Characteristic levels of the leaf area density and coefficients:'// & 1892 ' Height: ',A,' m'/ & 1893 ' Leaf area density: ',A,' m**2/m**3'/ & 1894 ' Coefficient alpha: ',F6.2 / & 1895 ' Coefficient beta: ',F6.2 / & 1896 ' Leaf area index: ',F6.2,' m**2/m**2' /) 1848 1897 1849 1898 300 FORMAT (//' Boundary conditions:'/ & -
palm/trunk/SOURCE/init_3d_model.f90
r1432 r1484 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Changes due to new module structure of the plant canopy model: 23 ! canopy-related initialization (e.g. lad and canopy_heat_flux) moved to new 24 ! subroutine init_plant_canopy within the module plant_canopy_model_mod, 25 ! call of subroutine init_plant_canopy added. 23 26 ! 24 27 ! Former revisions: … … 217 220 USE pegrid 218 221 222 USE plant_canopy_model_mod, & 223 ONLY: init_plant_canopy, plant_canopy 224 219 225 USE random_function_mod 220 226 … … 517 523 518 524 ! 519 !-- 3D-arrays for the leaf area density and the canopy drag coefficient520 IF ( plant_canopy ) THEN521 ALLOCATE ( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &522 lad_u(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &523 lad_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &524 lad_w(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &525 cdc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )526 527 IF ( passive_scalar ) THEN528 ALLOCATE ( sls(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &529 sec(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )530 ENDIF531 532 IF ( cthf /= 0.0_wp ) THEN533 ALLOCATE ( lai(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &534 canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )535 ENDIF536 537 ENDIF538 539 !540 525 !-- 4D-array for storing the Rif-values at vertical walls 541 526 IF ( topography /= 'flat' ) THEN … … 1571 1556 1572 1557 ! 1573 !-- Initialization of the leaf area density 1574 IF ( plant_canopy ) THEN 1575 1576 SELECT CASE ( TRIM( canopy_mode ) ) 1577 1578 CASE( 'block' ) 1579 1580 DO i = nxlg, nxrg 1581 DO j = nysg, nyng 1582 lad_s(:,j,i) = lad(:) 1583 cdc(:,j,i) = drag_coefficient 1584 IF ( passive_scalar ) THEN 1585 sls(:,j,i) = leaf_surface_concentration 1586 sec(:,j,i) = scalar_exchange_coefficient 1587 ENDIF 1588 ENDDO 1589 ENDDO 1590 1591 CASE DEFAULT 1592 1593 ! 1594 !-- The DEFAULT case is reached either if the parameter 1595 !-- canopy mode contains a wrong character string or if the 1596 !-- user has coded a special case in the user interface. 1597 !-- There, the subroutine user_init_plant_canopy checks 1598 !-- which of these two conditions applies. 1599 CALL user_init_plant_canopy 1600 1601 END SELECT 1602 1603 CALL exchange_horiz( lad_s, nbgp ) 1604 CALL exchange_horiz( cdc, nbgp ) 1605 1606 IF ( passive_scalar ) THEN 1607 CALL exchange_horiz( sls, nbgp ) 1608 CALL exchange_horiz( sec, nbgp ) 1609 ENDIF 1610 1611 ! 1612 !-- Sharp boundaries of the plant canopy in horizontal directions 1613 !-- In vertical direction the interpolation is retained, as the leaf 1614 !-- area density is initialised by prescribing a vertical profile 1615 !-- consisting of piecewise linear segments. The upper boundary 1616 !-- of the plant canopy is now defined by lad_w(pch_index,:,:) = 0.0. 1617 1618 DO i = nxl, nxr 1619 DO j = nys, nyn 1620 DO k = nzb, nzt+1 1621 IF ( lad_s(k,j,i) > 0.0_wp ) THEN 1622 lad_u(k,j,i) = lad_s(k,j,i) 1623 lad_u(k,j,i+1) = lad_s(k,j,i) 1624 lad_v(k,j,i) = lad_s(k,j,i) 1625 lad_v(k,j+1,i) = lad_s(k,j,i) 1626 ENDIF 1627 ENDDO 1628 DO k = nzb, nzt 1629 lad_w(k,j,i) = 0.5_wp * ( lad_s(k+1,j,i) + lad_s(k,j,i) ) 1630 ENDDO 1631 ENDDO 1632 ENDDO 1633 1634 lad_w(pch_index,:,:) = 0.0_wp 1635 lad_w(nzt+1,:,:) = lad_w(nzt,:,:) 1636 1637 CALL exchange_horiz( lad_u, nbgp ) 1638 CALL exchange_horiz( lad_v, nbgp ) 1639 CALL exchange_horiz( lad_w, nbgp ) 1640 1641 ! 1642 !-- Initialisation of the canopy heat source distribution 1643 IF ( cthf /= 0.0_wp ) THEN 1644 ! 1645 !-- Piecewise evaluation of the leaf area index by 1646 !-- integration of the leaf area density 1647 lai(:,:,:) = 0.0_wp 1648 DO i = nxlg, nxrg 1649 DO j = nysg, nyng 1650 DO k = pch_index-1, 0, -1 1651 lai(k,j,i) = lai(k+1,j,i) + & 1652 ( 0.5_wp * ( lad_w(k+1,j,i) + & 1653 lad_s(k+1,j,i) ) * & 1654 ( zw(k+1) - zu(k+1) ) ) + & 1655 ( 0.5_wp * ( lad_w(k,j,i) + & 1656 lad_s(k+1,j,i) ) * & 1657 ( zu(k+1) - zw(k) ) ) 1658 ENDDO 1659 ENDDO 1660 ENDDO 1661 1662 ! 1663 !-- Evaluation of the upward kinematic vertical heat flux within the 1664 !-- canopy 1665 DO i = nxlg, nxrg 1666 DO j = nysg, nyng 1667 DO k = 0, pch_index 1668 canopy_heat_flux(k,j,i) = cthf * & 1669 exp( -0.6_wp * lai(k,j,i) ) 1670 ENDDO 1671 ENDDO 1672 ENDDO 1673 1674 ! 1675 !-- The near surface heat flux is derived from the heat flux 1676 !-- distribution within the canopy 1677 shf(:,:) = canopy_heat_flux(0,:,:) 1678 1679 ENDIF 1680 1681 ENDIF 1558 !-- If required, initialize quantities needed for the plant canopy model 1559 IF ( plant_canopy ) CALL init_plant_canopy 1682 1560 1683 1561 ! -
palm/trunk/SOURCE/modules.f90
r1469 r1484 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Changes due to new module structure of the plant canopy model: 23 ! canopy-model related parameters/variables moved to module 24 ! plant_canopy_model_mod 23 25 ! 24 26 ! Former revisions: … … 38 40 ! 1429 2014-07-15 12:53:45Z knoop 39 41 ! +ensemble_member_nr to prepare the random_generator for ensemble runs 40 ! 42 ! 41 43 ! 1382 2014-04-30 12:15:41Z boeske 42 44 ! Renamed variables which store large scale forcing tendencies … … 294 296 REAL(wp), DIMENSION(:), ALLOCATABLE :: & 295 297 c_u_m, c_u_m_l, c_v_m, c_v_m_l, c_w_m, c_w_m_l, ddzu, ddzu_pres, & 296 dd2zu, dzu, ddzw, dzw, hyp, inflow_damping_factor, l ad, l_grid,&298 dd2zu, dzu, ddzw, dzw, hyp, inflow_damping_factor, l_grid, & 297 299 nc_1d, nr_1d, ptdf_x, ptdf_y, p_surf, pt_surf, pt_1d, pt_init, & 298 300 qsws_surf, q_1d, q_init, q_surf, qc_1d, qr_1d, rdf, rdf_sc, & … … 313 315 314 316 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 315 canopy_heat_flux, cdc, d, de_dx, de_dy, de_dz, diss, diss_l_e,&317 d, de_dx, de_dy, de_dz, diss, diss_l_e, & 316 318 diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr, diss_l_sa, diss_l_u, & 317 319 diss_l_v, diss_l_w, flux_l_e, flux_l_nr, flux_l_pt, flux_l_q, & 318 flux_l_qr, flux_l_sa, flux_l_u, flux_l_v, flux_l_w, kh, km, lad_s,&319 l ad_u, lad_v, lad_w, lai, l_wall, p_loc, sec, sls, tend, tric,&320 flux_l_qr, flux_l_sa, flux_l_u, flux_l_v, flux_l_w, kh, km, & 321 l_wall, p_loc, tend, tric, & 320 322 u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n, v_m_r, v_m_s, w_m_l, & 321 323 w_m_n, w_m_r, w_m_s … … 533 535 bc_sa_t = 'neumann', & 534 536 bc_uv_b = 'dirichlet', bc_uv_t = 'dirichlet', & 535 canopy_mode = 'block', &536 537 cloud_scheme = 'seifert_beheng', & 537 538 coupling_mode = 'uncoupled', & … … 586 587 nr_timesteps_this_run = 0, & 587 588 nsor = 20, nsor_ini = 100, n_sor, normalizing_region = 0, & 588 nz_do3d = -9999, p ch_index = 0, prt_time_count = 0, &589 nz_do3d = -9999, prt_time_count = 0, & 589 590 recycling_plane, runnr = 0, & 590 591 skip_do_avs = 0, subdomain_size, terminate_coupled = 0, & … … 596 597 do3d_no(0:1) = 0, do3d_time_count(0:1), & 597 598 domask_no(max_masks,0:1) = 0, domask_time_count(max_masks,0:1),& 598 lad_vertical_gradient_level_ind(10) = -9999, &599 599 mask_size(max_masks,3) = -1, mask_size_l(max_masks,3) = -1, & 600 600 mask_start_l(max_masks,3) = -1, & … … 652 652 outflow_l = .FALSE., outflow_n = .FALSE., outflow_r = .FALSE., & 653 653 outflow_s = .FALSE., passive_scalar = .FALSE., & 654 p lant_canopy = .FALSE., prandtl_layer = .TRUE., &654 prandtl_layer = .TRUE., & 655 655 precipitation = .FALSE., & 656 656 radiation = .FALSE., random_heatflux = .FALSE., & … … 682 682 canyon_width_x = 9999999.9_wp, canyon_width_y = 9999999.9_wp, & 683 683 canyon_wall_left = 9999999.9_wp, canyon_wall_south = 9999999.9_wp, & 684 c thf = 0.0_wp, cfl_factor = -1.0_wp, cos_alpha_surface, &684 cfl_factor = -1.0_wp, cos_alpha_surface, & 685 685 coupling_start_time = 0.0_wp, disturbance_amplitude = 0.25_wp, & 686 686 disturbance_energy_limit = 0.01_wp, & 687 687 disturbance_level_b = -9999999.9_wp, & 688 688 disturbance_level_t = -9999999.9_wp, & 689 dp_level_b = 0.0_wp, drag_coefficient = 0.0_wp,&689 dp_level_b = 0.0_wp, & 690 690 dt = -1.0_wp, dt_averaging_input = 0.0_wp, & 691 691 dt_averaging_input_pr = 9999999.9_wp, dt_coupling = 9999999.9_wp, & … … 703 703 f = 0.0_wp, fs = 0.0_wp, g = 9.81_wp, inflow_damping_height = 9999999.9_wp, & 704 704 inflow_damping_width = 9999999.9_wp, kappa = 0.4_wp, km_constant = -1.0_wp,& 705 lad_surface = 0.0_wp, leaf_surface_concentration = 0.0_wp, &706 705 mask_scale_x = 1.0_wp, mask_scale_y = 1.0_wp, mask_scale_z = 1.0_wp, & 707 706 maximum_cpu_time_allowed = 0.0_wp, & … … 719 718 restart_time = 9999999.9_wp, rho_reference, rho_surface, & 720 719 rif_max = 1.0_wp, rif_min = -5.0_wp, roughness_length = 0.1_wp, & 721 sa_surface = 35.0_wp, scalar_exchange_coefficient = 0.0_wp,&720 sa_surface = 35.0_wp, & 722 721 simulated_time = 0.0_wp, simulated_time_at_begin, sin_alpha_surface, & 723 722 skip_time_data_output = 0.0_wp, skip_time_data_output_av = 9999999.9_wp,& … … 745 744 REAL(wp) :: do2d_xy_last_time(0:1) = -1.0_wp, do2d_xz_last_time(0:1) = -1.0_wp, & 746 745 do2d_yz_last_time(0:1) = -1.0_wp, dpdxy(1:2) = 0.0_wp, & 747 dt_domask(max_masks) = 9999999.9_wp, lad_vertical_gradient(10) = 0.0_wp,& 748 lad_vertical_gradient_level(10) = -9999999.9_wp, & 746 dt_domask(max_masks) = 9999999.9_wp, & 749 747 mask_scale(3), & 750 748 pt_vertical_gradient(10) = 0.0_wp, & -
palm/trunk/SOURCE/package_parin.f90
r1368 r1484 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Changes due to new module structure of the plant canopy model: 23 ! module plant_canopy_model_mod added, 24 ! new package/namelist canopy_par added, i.e. the canopy model is no longer 25 ! steered over the inipar-namelist, 26 ! drag_coefficient, leaf_surface_concentration and scalar_exchange_coefficient 27 ! renamed to canopy_drag_coeff, leaf_surface_conc and leaf_scalar_exch_coeff. 28 ! Changed statement tags in CONTINUE-statement 23 29 ! 24 30 ! Former revisions: … … 103 109 write_particle_statistics 104 110 111 USE plant_canopy_model_mod, & 112 ONLY: alpha_lad, beta_lad, calc_beta_lad_profile, canopy_drag_coeff, & 113 canopy_mode, cthf, lad_surface, & 114 lad_vertical_gradient, lad_vertical_gradient_level, lai_beta, & 115 leaf_scalar_exch_coeff, leaf_surface_conc, pch_index, & 116 plant_canopy 117 105 118 USE spectrum, & 106 119 ONLY: comp_spectra_level, data_output_sp, plot_spectra_level, & … … 110 123 111 124 CHARACTER (LEN=80) :: line !: 125 126 NAMELIST /canopy_par/ alpha_lad, beta_lad, canopy_drag_coeff, & 127 canopy_mode, cthf, & 128 lad_surface, & 129 lad_vertical_gradient, & 130 lad_vertical_gradient_level, & 131 lai_beta, & 132 leaf_scalar_exch_coeff, & 133 leaf_surface_conc, pch_index 112 134 113 135 NAMELIST /dvrp_graphics_par/ clip_dvrp_l, clip_dvrp_n, clip_dvrp_r, & … … 159 181 line = ' ' 160 182 183 ! 184 !-- Try to find canopy package 185 REWIND ( 11 ) 186 line = ' ' 187 DO WHILE ( INDEX( line, '&canopy_par' ) == 0 ) 188 READ ( 11, '(A)', END=10 ) line 189 ENDDO 190 BACKSPACE ( 11 ) 191 192 ! 193 !-- Read user-defined namelist 194 READ ( 11, canopy_par ) 195 196 ! 197 !-- Set flag that indicates that canopy model is switched on 198 plant_canopy = .TRUE. 199 200 ! 201 !-- Set flag that indicates that the lad-profile shall be calculated by using 202 !-- a beta probability density function 203 IF ( alpha_lad /= 9999999.9_wp .AND. beta_lad /= 9999999.9_wp ) & 204 calc_beta_lad_profile = .TRUE. 205 206 10 CONTINUE 207 208 161 209 #if defined( __dvrp_graphics ) 162 210 REWIND ( 11 ) 163 211 line = ' ' 164 212 DO WHILE ( INDEX( line, '&dvrp_graphics_par' ) == 0 ) 165 READ ( 11, '(A)', END= 10 ) line213 READ ( 11, '(A)', END=20 ) line 166 214 ENDDO 167 215 BACKSPACE ( 11 ) … … 171 219 READ ( 11, dvrp_graphics_par ) 172 220 173 10 CONTINUE221 20 CONTINUE 174 222 #endif 175 223 … … 179 227 line = ' ' 180 228 DO WHILE ( INDEX( line, '&particles_par' ) == 0 ) 181 READ ( 11, '(A)', END= 20 ) line229 READ ( 11, '(A)', END=30 ) line 182 230 ENDDO 183 231 BACKSPACE ( 11 ) … … 191 239 particle_advection = .TRUE. 192 240 193 20 CONTINUE241 30 CONTINUE 194 242 195 243 … … 198 246 line = ' ' 199 247 DO WHILE ( INDEX( line, '&spectra_par' ) == 0 ) 200 READ ( 11, '(A)', END= 30 ) line248 READ ( 11, '(A)', END=40 ) line 201 249 ENDDO 202 250 BACKSPACE ( 11 ) … … 211 259 IF ( dt_dosp == 9999999.9_wp ) dt_dosp = dt_data_output 212 260 213 30 CONTINUE261 40 CONTINUE 214 262 #endif 215 263 -
palm/trunk/SOURCE/parin.f90
r1430 r1484 19 19 ! 20 20 ! Current revisions: 21 ! ------------------ 22 ! 21 ! ----------------- 22 ! Changes due to new module structure of the plant canopy model: 23 ! canopy-model related parameters moved to new package canopy_par in 24 ! subroutine package_parin 23 25 ! 24 26 ! Former revisions: … … 151 153 152 154 USE arrays_3d, & 153 ONLY: lad, pt_init, q_init, ref_state, sa_init, ug, u_init, v_init,&155 ONLY: pt_init, q_init, ref_state, sa_init, ug, u_init, v_init, & 154 156 vg 155 157 … … 165 167 building_length_y, building_wall_left, building_wall_south, & 166 168 call_microphysics_at_all_substeps, call_psolver_at_all_substeps,& 167 can opy_mode, canyon_height,&169 canyon_height, & 168 170 canyon_width_x, canyon_width_y, canyon_wall_left, & 169 171 canyon_wall_south, cfl_factor, & 170 172 cloud_droplets, cloud_physics, cloud_scheme, & 171 173 conserve_volume_flow, conserve_volume_flow_mode, & 172 coupling_start_time, create_disturbances, c thf, cycle_mg,&174 coupling_start_time, create_disturbances, cycle_mg, & 173 175 data_output, data_output_masks, & 174 176 data_output_pr, data_output_2d_on_each_pe, & … … 176 178 disturbance_level_b, disturbance_level_t, dissipation_1d, & 177 179 do2d_at_begin, do3d_at_begin, & 178 dp_external, dp_level_b, dp_smooth, dpdxy, drag_coefficient,&180 dp_external, dp_level_b, dp_smooth, dpdxy, & 179 181 drizzle, dt, dt_averaging_input, dt_averaging_input_pr, & 180 182 dt_coupling, dt_data_output, dt_data_output_av, dt_disturb, & … … 187 189 inflow_damping_width, inflow_disturbance_begin, & 188 190 inflow_disturbance_end, initializing_actions, io_blocks, & 189 io_group, km_constant, lad_surface, lad_vertical_gradient,&190 la d_vertical_gradient_level, large_scale_forcing,&191 large_scale_subsidence, leaf_surface_concentration,&191 io_group, km_constant, & 192 large_scale_forcing, & 193 large_scale_subsidence, & 192 194 loop_optimization, masking_method, mask_scale_x, mask_scale_y, & 193 195 mask_scale_z, mask_x, mask_y, mask_z, mask_x_loop, & … … 197 199 momentum_advec, netcdf_data_format, netcdf_precision, neutral, & 198 200 ngsrb, normalizing_region, nsor, nsor_ini, nudging, ocean, & 199 omega, omega_sor, passive_scalar, p ch_index, phi, nz_do3d,&200 p lant_canopy, prandtl_layer, prandtl_number, precipitation,&201 omega, omega_sor, passive_scalar, phi, nz_do3d, & 202 prandtl_layer, prandtl_number, precipitation, & 201 203 precipitation_amount_interval, psolver, pt_damping_factor, & 202 204 pt_damping_width, pt_reference, pt_surface, & … … 211 213 run_identifier, sa_surface, sa_vertical_gradient, & 212 214 sa_vertical_gradient_level, scalar_advec, & 213 scalar_ exchange_coefficient, scalar_rayleigh_damping,&215 scalar_rayleigh_damping, & 214 216 section_xy, section_xz, section_yz, skip_time_data_output, & 215 217 skip_time_data_output_av, skip_time_dopr, skip_time_do2d_xy, & … … 270 272 building_length_y, building_wall_left, building_wall_south, & 271 273 call_psolver_at_all_substeps, call_microphysics_at_all_substeps, & 272 can opy_mode, canyon_height,&274 canyon_height, & 273 275 canyon_width_x, canyon_width_y, canyon_wall_left, & 274 276 canyon_wall_south, c_sedimentation, cfl_factor, cloud_droplets, & 275 277 cloud_physics, cloud_scheme, collective_wait, & 276 278 conserve_volume_flow, & 277 conserve_volume_flow_mode, coupling_start_time, cthf,&279 conserve_volume_flow_mode, coupling_start_time, & 278 280 curvature_solution_effects, cycle_mg, damp_level_1d, & 279 281 dissipation_1d, & 280 dp_external, dp_level_b, dp_smooth, dpdxy, drag_coefficient,&282 dp_external, dp_level_b, dp_smooth, dpdxy, & 281 283 drizzle, dt, dt_pr_1d, dt_run_control_1d, dx, dy, dz, dz_max, & 282 284 dz_stretch_factor, dz_stretch_level, end_time_1d, & … … 285 287 inflow_damping_height, inflow_damping_width, & 286 288 inflow_disturbance_begin, inflow_disturbance_end, & 287 initializing_actions, km_constant, lad_surface, & 288 lad_vertical_gradient, lad_vertical_gradient_level, & 289 initializing_actions, km_constant, & 289 290 large_scale_forcing, large_scale_subsidence, & 290 l eaf_surface_concentration, limiter_sedimentation,&291 limiter_sedimentation, & 291 292 loop_optimization, masking_method, mg_cycles, & 292 293 mg_switch_to_pe0_level, mixing_length_1d, momentum_advec, & 293 294 nc_const, netcdf_precision, neutral, ngsrb, & 294 295 nsor, nsor_ini, nudging, nx, ny, nz, ocean, omega, omega_sor, & 295 passive_scalar, p ch_index, phi, plant_canopy, prandtl_layer,&296 passive_scalar, phi, prandtl_layer, & 296 297 prandtl_number, precipitation, psolver, pt_damping_factor, & 297 298 pt_damping_width, pt_reference, pt_surface, & … … 304 305 rif_max, rif_min, roughness_length, sa_surface, & 305 306 sa_vertical_gradient, sa_vertical_gradient_level, scalar_advec, & 306 scalar_ exchange_coefficient, scalar_rayleigh_damping,&307 scalar_rayleigh_damping, & 307 308 statistic_regions, subs_vertical_gradient, & 308 309 subs_vertical_gradient_level, surface_heatflux, surface_pressure, & … … 478 479 !-- ATTENTION: in case of changes to the following statement please 479 480 !-- also check the allocate statement in routine read_var_list 480 ALLOCATE( lad(0:nz+1), pt_init(0:nz+1), q_init(0:nz+1),&481 ALLOCATE( pt_init(0:nz+1), q_init(0:nz+1), & 481 482 ref_state(0:nz+1), sa_init(0:nz+1), ug(0:nz+1), & 482 483 u_init(0:nz+1), v_init(0:nz+1), vg(0:nz+1), & -
palm/trunk/SOURCE/plant_canopy_model.f90
r1341 r1484 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Changes due to new module structure of the plant canopy model: 23 ! module plant_canopy_model_mod now contains a subroutine for the 24 ! initialization of the canopy model (init_plant_canopy), 25 ! limitation of the canopy drag (previously accounted for by calculation of 26 ! a limiting canopy timestep for the determination of the maximum LES timestep 27 ! in subroutine timestep) is now realized by the calculation of pre-tendencies 28 ! and preliminary velocities in subroutine plant_canopy_model, 29 ! some redundant MPI communication removed in subroutine init_plant_canopy 30 ! (was previously in init_3d_model), 31 ! unnecessary 3d-arrays lad_u, lad_v, lad_w removed - lad information on the 32 ! respective grid is now provided only by lad_s (e.g. in the calculation of 33 ! the tendency terms or of cum_lai_hf), 34 ! drag_coefficient, lai, leaf_surface_concentration, 35 ! scalar_exchange_coefficient, sec and sls renamed to canopy_drag_coeff, 36 ! cum_lai_hf, leaf_surface_conc, leaf_scalar_exch_coeff, lsec and lsc, 37 ! respectively, 38 ! unnecessary 3d-arrays cdc, lsc and lsec now defined as single-value constants, 39 ! USE-statements and ONLY-lists modified accordingly 23 40 ! 24 41 ! Former revisions: … … 46 63 ! Description: 47 64 ! ------------ 48 ! Evaluation of sinks and sources of momentum, heat and scalar concentration 49 ! due to canopy elements 65 ! 1) Initialization of the canopy model, e.g. construction of leaf area density 66 ! profile (subroutine init_plant_canopy). 67 ! 2) Calculation of sinks and sources of momentum, heat and scalar concentration 68 ! due to canopy elements (subroutine plant_canopy_model). 50 69 !------------------------------------------------------------------------------! 70 USE arrays_3d, & 71 ONLY: dzu, dzw, e, q, shf, tend, u, v, w, zu, zw 72 73 USE indices, & 74 ONLY: nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv, & 75 nz, nzb, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt 76 77 USE kinds 78 79 80 IMPLICIT NONE 81 82 83 CHARACTER (LEN=20) :: canopy_mode = 'block' !: canopy coverage 84 85 INTEGER(iwp) :: pch_index = 0 !: plant canopy height/top index 86 INTEGER(iwp) :: & 87 lad_vertical_gradient_level_ind(10) = -9999 !: lad-profile levels (index) 88 89 LOGICAL :: calc_beta_lad_profile = .FALSE. !: switch for calc. of lad from beta func. 90 LOGICAL :: plant_canopy = .FALSE. !: switch for use of canopy model 91 92 REAL(wp) :: alpha_lad = 9999999.9_wp !: coefficient for lad calculation 93 REAL(wp) :: beta_lad = 9999999.9_wp !: coefficient for lad calculation 94 REAL(wp) :: canopy_drag_coeff = 0.0_wp !: canopy drag coefficient (parameter) 95 REAL(wp) :: cdc = 0.0_wp !: canopy drag coeff. (abbreviation used in equations) 96 REAL(wp) :: cthf = 0.0_wp !: canopy top heat flux 97 REAL(wp) :: dt_plant_canopy = 0.0_wp !: timestep account. for canopy drag 98 REAL(wp) :: lad_surface = 0.0_wp !: lad surface value 99 REAL(wp) :: lai_beta = 0.0_wp !: leaf area index (lai) for lad calc. 100 REAL(wp) :: & 101 leaf_scalar_exch_coeff = 0.0_wp !: canopy scalar exchange coeff. 102 REAL(wp) :: & 103 leaf_surface_conc = 0.0_wp !: leaf surface concentration 104 REAL(wp) :: lsec = 0.0_wp !: leaf scalar exchange coeff. 105 REAL(wp) :: lsc = 0.0_wp !: leaf surface concentration 106 107 REAL(wp) :: & 108 lad_vertical_gradient(10) = 0.0_wp !: lad gradient 109 REAL(wp) :: & 110 lad_vertical_gradient_level(10) = -9999999.9_wp !: lad-prof. levels (in m) 111 112 REAL(wp), DIMENSION(:), ALLOCATABLE :: lad !: leaf area density 113 REAL(wp), DIMENSION(:), ALLOCATABLE :: pre_lad !: preliminary lad 114 115 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 116 canopy_heat_flux !: canopy heat flux 117 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: cum_lai_hf !: cumulative lai for heatflux calc. 118 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: lad_s !: lad on scalar-grid 119 120 121 SAVE 122 51 123 52 124 PRIVATE 53 PUBLIC plant_canopy_model 125 PUBLIC alpha_lad, beta_lad, calc_beta_lad_profile, canopy_drag_coeff, & 126 canopy_mode, cdc, cthf, dt_plant_canopy, init_plant_canopy, lad, & 127 lad_s, lad_surface, lad_vertical_gradient, & 128 lad_vertical_gradient_level, lad_vertical_gradient_level_ind, & 129 lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, lsc, lsec, & 130 pch_index, plant_canopy, plant_canopy_model 131 132 133 INTERFACE init_plant_canopy 134 MODULE PROCEDURE init_plant_canopy 135 END INTERFACE init_plant_canopy 54 136 55 137 INTERFACE plant_canopy_model … … 58 140 END INTERFACE plant_canopy_model 59 141 142 143 144 60 145 CONTAINS 61 146 62 147 63 148 !------------------------------------------------------------------------------! 64 ! Call for all grid points 149 ! Description: 150 ! ------------ 151 !-- Initialization of the plant canopy model 152 !------------------------------------------------------------------------------! 153 SUBROUTINE init_plant_canopy 154 155 156 USE control_parameters, & 157 ONLY: dz, ocean, passive_scalar 158 159 160 IMPLICIT NONE 161 162 INTEGER(iwp) :: i !: running index 163 INTEGER(iwp) :: j !: running index 164 INTEGER(iwp) :: k !: running index 165 166 REAL(wp) :: int_bpdf !: vertical integral for lad-profile construction 167 REAL(wp) :: dzh !: vertical grid spacing in units of canopy height 168 REAL(wp) :: gradient !: gradient for lad-profile construction 169 REAL(wp) :: canopy_height !: canopy height for lad-profile construction 170 171 ! 172 !-- Allocate one-dimensional arrays for the computation of the 173 !-- leaf area density (lad) profile 174 ALLOCATE( lad(0:nz+1), pre_lad(0:nz+1) ) 175 lad = 0.0_wp 176 pre_lad = 0.0_wp 177 178 ! 179 !-- Compute the profile of leaf area density used in the plant 180 !-- canopy model. The profile can either be constructed from 181 !-- prescribed vertical gradients of the leaf area density or by 182 !-- using a beta probability density function (see e.g. Markkanen et al., 183 !-- 2003: Boundary-Layer Meteorology, 106, 437-459) 184 IF ( .NOT. calc_beta_lad_profile ) THEN 185 186 ! 187 !-- Use vertical gradients for lad-profile construction 188 i = 1 189 gradient = 0.0_wp 190 191 IF ( .NOT. ocean ) THEN 192 193 lad(0) = lad_surface 194 lad_vertical_gradient_level_ind(1) = 0 195 196 DO k = 1, pch_index 197 IF ( i < 11 ) THEN 198 IF ( lad_vertical_gradient_level(i) < zu(k) .AND. & 199 lad_vertical_gradient_level(i) >= 0.0_wp ) THEN 200 gradient = lad_vertical_gradient(i) 201 lad_vertical_gradient_level_ind(i) = k - 1 202 i = i + 1 203 ENDIF 204 ENDIF 205 IF ( gradient /= 0.0_wp ) THEN 206 IF ( k /= 1 ) THEN 207 lad(k) = lad(k-1) + dzu(k) * gradient 208 ELSE 209 lad(k) = lad_surface + dzu(k) * gradient 210 ENDIF 211 ELSE 212 lad(k) = lad(k-1) 213 ENDIF 214 ENDDO 215 216 ENDIF 217 218 ! 219 !-- In case of no given leaf area density gradients, choose a vanishing 220 !-- gradient. This information is used for the HEADER and the RUN_CONTROL 221 !-- file. 222 IF ( lad_vertical_gradient_level(1) == -9999999.9_wp ) THEN 223 lad_vertical_gradient_level(1) = 0.0_wp 224 ENDIF 225 226 ELSE 227 228 ! 229 !-- Use beta function for lad-profile construction 230 int_bpdf = 0.0_wp 231 canopy_height = pch_index * dz 232 233 DO k = nzb, pch_index 234 int_bpdf = int_bpdf + & 235 ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) * & 236 ( ( 1.0_wp - ( zw(k) / canopy_height ) )**( beta_lad-1.0_wp ) ) * & 237 ( ( zw(k+1)-zw(k) ) / canopy_height ) ) 238 ENDDO 239 240 ! 241 !-- Preliminary lad profile (defined on w-grid) 242 DO k = nzb, pch_index 243 pre_lad(k) = lai_beta * & 244 ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) * & 245 ( ( 1.0_wp - ( zw(k) / canopy_height ) )**( beta_lad-1.0_wp ) ) / & 246 int_bpdf & 247 ) / canopy_height 248 ENDDO 249 250 ! 251 !-- Final lad profile (defined on scalar-grid level, since most prognostic 252 !-- quantities are defined there, hence, less interpolation is required 253 !-- when calculating the canopy tendencies) 254 lad(0) = pre_lad(0) 255 DO k = nzb+1, pch_index 256 lad(k) = 0.5 * ( pre_lad(k-1) + pre_lad(k) ) 257 ENDDO 258 259 ENDIF 260 261 ! 262 !-- Allocate 3D-array for the leaf area density (lad_s). In case of a 263 !-- prescribed canopy-top heat flux (cthf), allocate 3D-arrays for 264 !-- the cumulative leaf area index (cum_lai_hf) and the canopy heat flux. 265 ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 266 267 IF ( cthf /= 0.0_wp ) THEN 268 ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 269 canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 270 ENDIF 271 272 ! 273 !-- Initialize canopy parameters cdc (canopy drag coefficient), 274 !-- lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration) 275 !-- with the prescribed values 276 cdc = canopy_drag_coeff 277 lsec = leaf_scalar_exch_coeff 278 lsc = leaf_surface_conc 279 280 ! 281 !-- Initialization of the canopy coverage in the model domain: 282 !-- Setting the parameter canopy_mode = 'block' initializes a canopy, which 283 !-- fully covers the domain surface 284 SELECT CASE ( TRIM( canopy_mode ) ) 285 286 CASE( 'block' ) 287 288 DO i = nxlg, nxrg 289 DO j = nysg, nyng 290 lad_s(:,j,i) = lad(:) 291 ENDDO 292 ENDDO 293 294 CASE DEFAULT 295 296 ! 297 !-- The DEFAULT case is reached either if the parameter 298 !-- canopy mode contains a wrong character string or if the 299 !-- user has coded a special case in the user interface. 300 !-- There, the subroutine user_init_plant_canopy checks 301 !-- which of these two conditions applies. 302 CALL user_init_plant_canopy 303 304 END SELECT 305 306 ! 307 !-- Initialization of the canopy heat source distribution 308 IF ( cthf /= 0.0_wp ) THEN 309 ! 310 !-- Piecewise calculation of the leaf area index by vertical 311 !-- integration of the leaf area density 312 cum_lai_hf(:,:,:) = 0.0_wp 313 DO i = nxlg, nxrg 314 DO j = nysg, nyng 315 DO k = pch_index-1, 0, -1 316 IF ( k == pch_index-1 ) THEN 317 cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) + & 318 ( 0.5_wp * lad_s(k+1,j,i) * & 319 ( zw(k+1) - zu(k+1) ) ) + & 320 ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) + & 321 lad_s(k,j,i) ) + & 322 lad_s(k+1,j,i) ) * & 323 ( zu(k+1) - zw(k) ) ) 324 ELSE 325 cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) + & 326 ( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) + & 327 lad_s(k+1,j,i) ) + & 328 lad_s(k+1,j,i) ) * & 329 ( zw(k+1) - zu(k+1) ) ) + & 330 ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) + & 331 lad_s(k,j,i) ) + & 332 lad_s(k+1,j,i) ) * & 333 ( zu(k+1) - zw(k) ) ) 334 ENDIF 335 ENDDO 336 ENDDO 337 ENDDO 338 339 ! 340 !-- Calculation of the upward kinematic vertical heat flux within the 341 !-- canopy 342 DO i = nxlg, nxrg 343 DO j = nysg, nyng 344 DO k = 0, pch_index 345 canopy_heat_flux(k,j,i) = cthf * & 346 exp( -0.6_wp * cum_lai_hf(k,j,i) ) 347 ENDDO 348 ENDDO 349 ENDDO 350 351 ! 352 !-- The surface heat flux is set to the surface value of the calculated 353 !-- in-canopy heat flux distribution 354 shf(:,:) = canopy_heat_flux(0,:,:) 355 356 ENDIF 357 358 359 360 END SUBROUTINE init_plant_canopy 361 362 363 364 !------------------------------------------------------------------------------! 365 ! Description: 366 ! ------------ 367 !-- Calculation of the tendency terms, accounting for the effect of the plant 368 !-- canopy on momentum and scalar quantities. 369 !-- 370 !-- The canopy is located where the leaf area density lad_s(k,j,i) > 0.0 371 !-- (defined on scalar grid), as initialized in subroutine init_plant_canopy. 372 !-- The lad on the w-grid is vertically interpolated from the surrounding 373 !-- lad_s. The upper boundary of the canopy is defined on the w-grid at 374 !-- k = pch_index. Here, the lad is zero. 375 !-- 376 !-- The canopy drag must be limited (previously accounted for by calculation of 377 !-- a limiting canopy timestep for the determination of the maximum LES timestep 378 !-- in subroutine timestep), since it is physically impossible that the canopy 379 !-- drag alone can locally change the sign of a velocity component. This 380 !-- limitation is realized by calculating preliminary tendencies and velocities. 381 !-- It is subsequently checked if the preliminary new velocity has a different 382 !-- sign than the current velocity. If so, the tendency is limited in a way that 383 !-- the velocity can at maximum be reduced to zero by the canopy drag. 384 !-- 385 !-- 386 !-- Call for all grid points 65 387 !------------------------------------------------------------------------------! 66 388 SUBROUTINE plant_canopy_model( component ) 67 389 68 USE arrays_3d, &69 ONLY: canopy_heat_flux, cdc, dzw, e, lad_s, lad_u, lad_v, lad_w, &70 q, sec, sls, tend, u, v, w71 390 72 391 USE control_parameters, & 73 ONLY: pch_index, message_string 74 75 USE indices, & 76 ONLY: nxl, nxlu, nxr, nys, nysv, nyn, nzb_s_inner, nzb_u_inner, & 77 nzb_v_inner, nzb_w_inner 392 ONLY: dt_3d, message_string 78 393 79 394 USE kinds … … 81 396 IMPLICIT NONE 82 397 83 INTEGER(iwp) :: component !: 84 INTEGER(iwp) :: i !: 85 INTEGER(iwp) :: j !: 86 INTEGER(iwp) :: k !: 398 INTEGER(iwp) :: component !: prognostic variable (u,v,w,pt,q,e) 399 INTEGER(iwp) :: i !: running index 400 INTEGER(iwp) :: j !: running index 401 INTEGER(iwp) :: k !: running index 402 403 REAL(wp) :: ddt_3d !: inverse of the LES timestep (dt_3d) 404 REAL(wp) :: lad_local !: local lad value 405 REAL(wp) :: pre_tend !: preliminary tendency 406 REAL(wp) :: pre_u !: preliminary u-value 407 REAL(wp) :: pre_v !: preliminary v-value 408 REAL(wp) :: pre_w !: preliminary w-value 409 410 411 ddt_3d = 1.0_wp / dt_3d 87 412 88 413 ! 89 !-- Compute drag for the three velocity components and the SGS-TKE 414 !-- Compute drag for the three velocity components and the SGS-TKE: 90 415 SELECT CASE ( component ) 91 416 … … 96 421 DO j = nys, nyn 97 422 DO k = nzb_u_inner(j,i)+1, pch_index 98 tend(k,j,i) = tend(k,j,i) - & 99 cdc(k,j,i) * lad_u(k,j,i) * & 100 SQRT( u(k,j,i)**2 + & 101 ( ( v(k,j,i-1) + & 102 v(k,j,i) + & 103 v(k,j+1,i) + & 104 v(k,j+1,i-1) ) & 105 / 4.0_wp )**2 + & 106 ( ( w(k-1,j,i-1) + & 107 w(k-1,j,i) + & 108 w(k,j,i-1) + & 109 w(k,j,i) ) & 110 / 4.0_wp )**2 ) * & 111 u(k,j,i) 423 424 ! 425 !-- In order to create sharp boundaries of the plant canopy, 426 !-- the lad on the u-grid at index (k,j,i) is equal to 427 !-- lad_s(k,j,i), rather than being interpolated from the 428 !-- surrounding lad_s, because this would yield smaller lad 429 !-- at the canopy boundaries than inside of the canopy. 430 !-- For the same reason, the lad at the rightmost(i+1)canopy 431 !-- boundary on the u-grid equals lad_s(k,j,i). 432 lad_local = lad_s(k,j,i) 433 IF ( lad_local == 0.0_wp .AND. & 434 lad_s(k,j,i-1) > 0.0_wp ) THEN 435 lad_local = lad_s(k,j,i-1) 436 ENDIF 437 438 pre_tend = 0.0_wp 439 pre_u = 0.0_wp 440 ! 441 !-- Calculate preliminary value (pre_tend) of the tendency 442 pre_tend = - cdc * & 443 lad_local * & 444 SQRT( u(k,j,i)**2 + & 445 ( 0.25_wp * ( v(k,j,i-1) + & 446 v(k,j,i) + & 447 v(k,j+1,i) + & 448 v(k,j+1,i-1) ) & 449 )**2 + & 450 ( 0.25_wp * ( w(k-1,j,i-1) + & 451 w(k-1,j,i) + & 452 w(k,j,i-1) + & 453 w(k,j,i) ) & 454 )**2 & 455 ) * & 456 u(k,j,i) 457 458 ! 459 !-- Calculate preliminary new velocity, based on pre_tend 460 pre_u = u(k,j,i) + dt_3d * pre_tend 461 ! 462 !-- Compare sign of old velocity and new preliminary velocity, 463 !-- and in case the signs are different, limit the tendency 464 IF ( SIGN(pre_u,u(k,j,i)) /= pre_u ) THEN 465 pre_tend = - u(k,j,i) * ddt_3d 466 ELSE 467 pre_tend = pre_tend 468 ENDIF 469 ! 470 !-- Calculate final tendency 471 tend(k,j,i) = tend(k,j,i) + pre_tend 472 112 473 ENDDO 113 474 ENDDO … … 120 481 DO j = nysv, nyn 121 482 DO k = nzb_v_inner(j,i)+1, pch_index 122 tend(k,j,i) = tend(k,j,i) - & 123 cdc(k,j,i) * lad_v(k,j,i) * & 124 SQRT( ( ( u(k,j-1,i) + & 125 u(k,j-1,i+1) + & 126 u(k,j,i) + & 127 u(k,j,i+1) ) & 128 / 4.0_wp )**2 + & 129 v(k,j,i)**2 + & 130 ( ( w(k-1,j-1,i) + & 131 w(k-1,j,i) + & 132 w(k,j-1,i) + & 133 w(k,j,i) ) & 134 / 4.0_wp )**2 ) * & 135 v(k,j,i) 483 484 ! 485 !-- In order to create sharp boundaries of the plant canopy, 486 !-- the lad on the v-grid at index (k,j,i) is equal to 487 !-- lad_s(k,j,i), rather than being interpolated from the 488 !-- surrounding lad_s, because this would yield smaller lad 489 !-- at the canopy boundaries than inside of the canopy. 490 !-- For the same reason, the lad at the northmost(j+1) canopy 491 !-- boundary on the v-grid equals lad_s(k,j,i). 492 lad_local = lad_s(k,j,i) 493 IF ( lad_local == 0.0_wp .AND. & 494 lad_s(k,j-1,i) > 0.0_wp ) THEN 495 lad_local = lad_s(k,j-1,i) 496 ENDIF 497 498 pre_tend = 0.0_wp 499 pre_v = 0.0_wp 500 ! 501 !-- Calculate preliminary value (pre_tend) of the tendency 502 pre_tend = - cdc * & 503 lad_local * & 504 SQRT( ( 0.25_wp * ( u(k,j-1,i) + & 505 u(k,j-1,i+1) + & 506 u(k,j,i) + & 507 u(k,j,i+1) ) & 508 )**2 + & 509 v(k,j,i)**2 + & 510 ( 0.25_wp * ( w(k-1,j-1,i) + & 511 w(k-1,j,i) + & 512 w(k,j-1,i) + & 513 w(k,j,i) ) & 514 )**2 & 515 ) * & 516 v(k,j,i) 517 518 ! 519 !-- Calculate preliminary new velocity, based on pre_tend 520 pre_v = v(k,j,i) + dt_3d * pre_tend 521 ! 522 !-- Compare sign of old velocity and new preliminary velocity, 523 !-- and in case the signs are different, limit the tendency 524 IF ( SIGN(pre_v,v(k,j,i)) /= pre_v ) THEN 525 pre_tend = - v(k,j,i) * ddt_3d 526 ELSE 527 pre_tend = pre_tend 528 ENDIF 529 ! 530 !-- Calculate final tendency 531 tend(k,j,i) = tend(k,j,i) + pre_tend 532 136 533 ENDDO 137 534 ENDDO … … 143 540 DO i = nxl, nxr 144 541 DO j = nys, nyn 145 DO k = nzb_w_inner(j,i)+1, pch_index 146 tend(k,j,i) = tend(k,j,i) - & 147 cdc(k,j,i) * lad_w(k,j,i) * & 148 SQRT( ( ( u(k,j,i) + & 149 u(k,j,i+1) + & 150 u(k+1,j,i) + & 151 u(k+1,j,i+1) ) & 152 / 4.0_wp )**2 + & 153 ( ( v(k,j,i) + & 154 v(k,j+1,i) + & 155 v(k+1,j,i) + & 156 v(k+1,j+1,i) ) & 157 / 4.0_wp )**2 + & 158 w(k,j,i)**2 ) * & 159 w(k,j,i) 542 DO k = nzb_w_inner(j,i)+1, pch_index-1 543 544 pre_tend = 0.0_wp 545 pre_w = 0.0_wp 546 ! 547 !-- Calculate preliminary value (pre_tend) of the tendency 548 pre_tend = - cdc * & 549 (0.5_wp * & 550 ( lad_s(k+1,j,i) + lad_s(k,j,i) )) * & 551 SQRT( ( 0.25_wp * ( u(k,j,i) + & 552 u(k,j,i+1) + & 553 u(k+1,j,i) + & 554 u(k+1,j,i+1) ) & 555 )**2 + & 556 ( 0.25_wp * ( v(k,j,i) + & 557 v(k,j+1,i) + & 558 v(k+1,j,i) + & 559 v(k+1,j+1,i) ) & 560 )**2 + & 561 w(k,j,i)**2 & 562 ) * & 563 w(k,j,i) 564 ! 565 !-- Calculate preliminary new velocity, based on pre_tend 566 pre_w = w(k,j,i) + dt_3d * pre_tend 567 ! 568 !-- Compare sign of old velocity and new preliminary velocity, 569 !-- and in case the signs are different, limit the tendency 570 IF ( SIGN(pre_w,w(k,j,i)) /= pre_w ) THEN 571 pre_tend = - w(k,j,i) * ddt_3d 572 ELSE 573 pre_tend = pre_tend 574 ENDIF 575 ! 576 !-- Calculate final tendency 577 tend(k,j,i) = tend(k,j,i) + pre_tend 578 160 579 ENDDO 161 580 ENDDO … … 168 587 DO j = nys, nyn 169 588 DO k = nzb_s_inner(j,i)+1, pch_index 170 tend(k,j,i) = tend(k,j,i) + & 171 ( canopy_heat_flux(k,j,i) - & 172 canopy_heat_flux(k-1,j,i) ) / & 173 dzw(k) 589 tend(k,j,i) = tend(k,j,i) + & 590 ( canopy_heat_flux(k,j,i) - & 591 canopy_heat_flux(k-1,j,i) ) / dzw(k) 174 592 ENDDO 175 593 ENDDO … … 182 600 DO j = nys, nyn 183 601 DO k = nzb_s_inner(j,i)+1, pch_index 184 tend(k,j,i) = tend(k,j,i) - & 185 sec(k,j,i) * lad_s(k,j,i) * & 186 SQRT( ( ( u(k,j,i) + & 187 u(k,j,i+1) ) & 188 / 2.0_wp )**2 + & 189 ( ( v(k,j,i) + & 190 v(k,j+1,i) ) & 191 / 2.0_wp )**2 + & 192 ( ( w(k-1,j,i) + & 193 w(k,j,i) ) & 194 / 2.0_wp )**2 ) * & 195 ( q(k,j,i) - sls(k,j,i) ) 602 tend(k,j,i) = tend(k,j,i) - & 603 lsec * & 604 lad_s(k,j,i) * & 605 SQRT( ( 0.5_wp * ( u(k,j,i) + & 606 u(k,j,i+1) ) & 607 )**2 + & 608 ( 0.5_wp * ( v(k,j,i) + & 609 v(k,j+1,i) ) & 610 )**2 + & 611 ( 0.5_wp * ( w(k-1,j,i) + & 612 w(k,j,i) ) & 613 )**2 & 614 ) * & 615 ( q(k,j,i) - lsc ) 196 616 ENDDO 197 617 ENDDO … … 204 624 DO j = nys, nyn 205 625 DO k = nzb_s_inner(j,i)+1, pch_index 206 tend(k,j,i) = tend(k,j,i) - & 207 2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * & 208 SQRT( ( ( u(k,j,i) + & 209 u(k,j,i+1) ) & 210 / 2.0_wp )**2 + & 211 ( ( v(k,j,i) + & 212 v(k,j+1,i) ) & 213 / 2.0_wp )**2 + & 214 ( ( w(k,j,i) + & 215 w(k+1,j,i) ) & 216 / 2.0_wp )**2 ) * & 217 e(k,j,i) 626 tend(k,j,i) = tend(k,j,i) - & 627 2.0_wp * cdc * & 628 lad_s(k,j,i) * & 629 SQRT( ( 0.5_wp * ( u(k,j,i) + & 630 u(k,j,i+1) ) & 631 )**2 + & 632 ( 0.5_wp * ( v(k,j,i) + & 633 v(k,j+1,i) ) & 634 )**2 + & 635 ( 0.5_wp * ( w(k,j,i) + & 636 w(k+1,j,i) ) & 637 )**2 & 638 ) * & 639 e(k,j,i) 218 640 ENDDO 219 641 ENDDO 220 642 ENDDO 221 643 644 222 645 CASE DEFAULT 223 646 … … 231 654 232 655 !------------------------------------------------------------------------------! 233 ! Call for grid point i,j 656 ! Description: 657 ! ------------ 658 !-- Calculation of the tendency terms, accounting for the effect of the plant 659 !-- canopy on momentum and scalar quantities. 660 !-- 661 !-- The canopy is located where the leaf area density lad_s(k,j,i) > 0.0 662 !-- (defined on scalar grid), as initialized in subroutine init_plant_canopy. 663 !-- The lad on the w-grid is vertically interpolated from the surrounding 664 !-- lad_s. The upper boundary of the canopy is defined on the w-grid at 665 !-- k = pch_index. Here, the lad is zero. 666 !-- 667 !-- The canopy drag must be limited (previously accounted for by calculation of 668 !-- a limiting canopy timestep for the determination of the maximum LES timestep 669 !-- in subroutine timestep), since it is physically impossible that the canopy 670 !-- drag alone can locally change the sign of a velocity component. This 671 !-- limitation is realized by calculating preliminary tendencies and velocities. 672 !-- It is subsequently checked if the preliminary new velocity has a different 673 !-- sign than the current velocity. If so, the tendency is limited in a way that 674 !-- the velocity can at maximum be reduced to zero by the canopy drag. 675 !-- 676 !-- 677 !-- Call for grid point i,j 234 678 !------------------------------------------------------------------------------! 235 679 SUBROUTINE plant_canopy_model_ij( i, j, component ) 236 680 237 USE arrays_3d, &238 ONLY: canopy_heat_flux, cdc, dzw, e, lad_s, lad_u, lad_v, lad_w, &239 q, sec, sls, tend, u, v, w240 681 241 682 USE control_parameters, & 242 ONLY: pch_index, message_string 243 244 USE indices, & 245 ONLY: nxl, nxlu, nxr, nys, nysv, nyn, nzb_s_inner, nzb_u_inner, & 246 nzb_v_inner, nzb_w_inner 683 ONLY: dt_3d, message_string 247 684 248 685 USE kinds … … 250 687 IMPLICIT NONE 251 688 252 INTEGER(iwp) :: component !: 253 INTEGER(iwp) :: i !: 254 INTEGER(iwp) :: j !: 255 INTEGER(iwp) :: k !: 256 257 ! 258 !-- Compute drag for the three velocity components 689 INTEGER(iwp) :: component !: prognostic variable (u,v,w,pt,q,e) 690 INTEGER(iwp) :: i !: running index 691 INTEGER(iwp) :: j !: running index 692 INTEGER(iwp) :: k !: running index 693 694 REAL(wp) :: ddt_3d !: inverse of the LES timestep (dt_3d) 695 REAL(wp) :: lad_local !: local lad value 696 REAL(wp) :: pre_tend !: preliminary tendency 697 REAL(wp) :: pre_u !: preliminary u-value 698 REAL(wp) :: pre_v !: preliminary v-value 699 REAL(wp) :: pre_w !: preliminary w-value 700 701 702 ddt_3d = 1.0_wp / dt_3d 703 704 ! 705 !-- Compute drag for the three velocity components and the SGS-TKE 259 706 SELECT CASE ( component ) 260 707 261 708 ! 262 709 !-- u-component 263 CASE ( 1 ) 264 DO k = nzb_u_inner(j,i)+1, pch_index 265 tend(k,j,i) = tend(k,j,i) - & 266 cdc(k,j,i) * lad_u(k,j,i) * & 267 SQRT( u(k,j,i)**2 + & 268 ( ( v(k,j,i-1) + & 269 v(k,j,i) + & 270 v(k,j+1,i) + & 271 v(k,j+1,i-1) ) & 272 / 4.0_wp )**2 + & 273 ( ( w(k-1,j,i-1) + & 274 w(k-1,j,i) + & 275 w(k,j,i-1) + & 276 w(k,j,i) ) & 277 / 4.0_wp )**2 ) * & 278 u(k,j,i) 279 ENDDO 710 CASE ( 1 ) 711 DO k = nzb_u_inner(j,i)+1, pch_index 712 713 ! 714 !-- In order to create sharp boundaries of the plant canopy, 715 !-- the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i), 716 !-- rather than being interpolated from the surrounding lad_s, 717 !-- because this would yield smaller lad at the canopy boundaries 718 !-- than inside of the canopy. 719 !-- For the same reason, the lad at the rightmost(i+1)canopy 720 !-- boundary on the u-grid equals lad_s(k,j,i). 721 lad_local = lad_s(k,j,i) 722 IF ( lad_local == 0.0_wp .AND. & 723 lad_s(k,j,i-1) > 0.0_wp ) THEN 724 lad_local = lad_s(k,j,i-1) 725 ENDIF 726 727 pre_tend = 0.0_wp 728 pre_u = 0.0_wp 729 ! 730 !-- Calculate preliminary value (pre_tend) of the tendency 731 pre_tend = - cdc * & 732 lad_local * & 733 SQRT( u(k,j,i)**2 + & 734 ( 0.25_wp * ( v(k,j,i-1) + & 735 v(k,j,i) + & 736 v(k,j+1,i) + & 737 v(k,j+1,i-1) ) & 738 )**2 + & 739 ( 0.25_wp * ( w(k-1,j,i-1) + & 740 w(k-1,j,i) + & 741 w(k,j,i-1) + & 742 w(k,j,i) ) & 743 )**2 & 744 ) * & 745 u(k,j,i) 746 747 ! 748 !-- Calculate preliminary new velocity, based on pre_tend 749 pre_u = u(k,j,i) + dt_3d * pre_tend 750 ! 751 !-- Compare sign of old velocity and new preliminary velocity, 752 !-- and in case the signs are different, limit the tendency 753 IF ( SIGN(pre_u,u(k,j,i)) /= pre_u ) THEN 754 pre_tend = - u(k,j,i) * ddt_3d 755 ELSE 756 pre_tend = pre_tend 757 ENDIF 758 ! 759 !-- Calculate final tendency 760 tend(k,j,i) = tend(k,j,i) + pre_tend 761 ENDDO 762 280 763 281 764 ! 282 765 !-- v-component 283 CASE ( 2 ) 284 DO k = nzb_v_inner(j,i)+1, pch_index 285 tend(k,j,i) = tend(k,j,i) - & 286 cdc(k,j,i) * lad_v(k,j,i) * & 287 SQRT( ( ( u(k,j-1,i) + & 288 u(k,j-1,i+1) + & 289 u(k,j,i) + & 290 u(k,j,i+1) ) & 291 / 4.0_wp )**2 + & 292 v(k,j,i)**2 + & 293 ( ( w(k-1,j-1,i) + & 294 w(k-1,j,i) + & 295 w(k,j-1,i) + & 296 w(k,j,i) ) & 297 / 4.0_wp )**2 ) * & 298 v(k,j,i) 299 ENDDO 766 CASE ( 2 ) 767 DO k = nzb_v_inner(j,i)+1, pch_index 768 769 ! 770 !-- In order to create sharp boundaries of the plant canopy, 771 !-- the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i), 772 !-- rather than being interpolated from the surrounding lad_s, 773 !-- because this would yield smaller lad at the canopy boundaries 774 !-- than inside of the canopy. 775 !-- For the same reason, the lad at the northmost(j+1)canopy 776 !-- boundary on the v-grid equals lad_s(k,j,i). 777 lad_local = lad_s(k,j,i) 778 IF ( lad_local == 0.0_wp .AND. & 779 lad_s(k,j-1,i) > 0.0_wp ) THEN 780 lad_local = lad_s(k,j-1,i) 781 ENDIF 782 783 pre_tend = 0.0_wp 784 pre_v = 0.0_wp 785 ! 786 !-- Calculate preliminary value (pre_tend) of the tendency 787 pre_tend = - cdc * & 788 lad_local * & 789 SQRT( ( 0.25_wp * ( u(k,j-1,i) + & 790 u(k,j-1,i+1) + & 791 u(k,j,i) + & 792 u(k,j,i+1) ) & 793 )**2 + & 794 v(k,j,i)**2 + & 795 ( 0.25_wp * ( w(k-1,j-1,i) + & 796 w(k-1,j,i) + & 797 w(k,j-1,i) + & 798 w(k,j,i) ) & 799 )**2 & 800 ) * & 801 v(k,j,i) 802 803 ! 804 !-- Calculate preliminary new velocity, based on pre_tend 805 pre_v = v(k,j,i) + dt_3d * pre_tend 806 ! 807 !-- Compare sign of old velocity and new preliminary velocity, 808 !-- and in case the signs are different, limit the tendency 809 IF ( SIGN(pre_v,v(k,j,i)) /= pre_v ) THEN 810 pre_tend = - v(k,j,i) * ddt_3d 811 ELSE 812 pre_tend = pre_tend 813 ENDIF 814 ! 815 !-- Calculate final tendency 816 tend(k,j,i) = tend(k,j,i) + pre_tend 817 ENDDO 818 300 819 301 820 ! 302 821 !-- w-component 303 CASE ( 3 ) 304 DO k = nzb_w_inner(j,i)+1, pch_index 305 tend(k,j,i) = tend(k,j,i) - & 306 cdc(k,j,i) * lad_w(k,j,i) * & 307 SQRT( ( ( u(k,j,i) + & 308 u(k,j,i+1) + & 309 u(k+1,j,i) + & 310 u(k+1,j,i+1) ) & 311 / 4.0_wp )**2 + & 312 ( ( v(k,j,i) + & 313 v(k,j+1,i) + & 314 v(k+1,j,i) + & 315 v(k+1,j+1,i) ) & 316 / 4.0_wp )**2 + & 317 w(k,j,i)**2 ) * & 318 w(k,j,i) 319 320 ENDDO 822 CASE ( 3 ) 823 DO k = nzb_w_inner(j,i)+1, pch_index-1 824 825 pre_tend = 0.0_wp 826 pre_w = 0.0_wp 827 ! 828 !-- Calculate preliminary value (pre_tend) of the tendency 829 pre_tend = - cdc * & 830 (0.5_wp * & 831 ( lad_s(k+1,j,i) + lad_s(k,j,i) )) * & 832 SQRT( ( 0.25_wp * ( u(k,j,i) + & 833 u(k,j,i+1) + & 834 u(k+1,j,i) + & 835 u(k+1,j,i+1) ) & 836 )**2 + & 837 ( 0.25_wp * ( v(k,j,i) + & 838 v(k,j+1,i) + & 839 v(k+1,j,i) + & 840 v(k+1,j+1,i) ) & 841 )**2 + & 842 w(k,j,i)**2 & 843 ) * & 844 w(k,j,i) 845 ! 846 !-- Calculate preliminary new velocity, based on pre_tend 847 pre_w = w(k,j,i) + dt_3d * pre_tend 848 ! 849 !-- Compare sign of old velocity and new preliminary velocity, 850 !-- and in case the signs are different, limit the tendency 851 IF ( SIGN(pre_w,w(k,j,i)) /= pre_w ) THEN 852 pre_tend = - w(k,j,i) * ddt_3d 853 ELSE 854 pre_tend = pre_tend 855 ENDIF 856 ! 857 !-- Calculate final tendency 858 tend(k,j,i) = tend(k,j,i) + pre_tend 859 ENDDO 321 860 322 861 ! … … 324 863 CASE ( 4 ) 325 864 DO k = nzb_s_inner(j,i)+1, pch_index 326 tend(k,j,i) = tend(k,j,i) + & 327 ( canopy_heat_flux(k,j,i) - & 328 canopy_heat_flux(k-1,j,i) ) / & 329 dzw(k) 865 tend(k,j,i) = tend(k,j,i) + & 866 ( canopy_heat_flux(k,j,i) - & 867 canopy_heat_flux(k-1,j,i) ) / dzw(k) 330 868 ENDDO 331 869 … … 335 873 CASE ( 5 ) 336 874 DO k = nzb_s_inner(j,i)+1, pch_index 337 tend(k,j,i) = tend(k,j,i) - & 338 sec(k,j,i) * lad_s(k,j,i) * & 339 SQRT( ( ( u(k,j,i) + & 340 u(k,j,i+1) ) & 341 / 2.0_wp )**2 + & 342 ( ( v(k,j,i) + & 343 v(k,j+1,i) ) & 344 / 2.0_wp )**2 + & 345 ( ( w(k-1,j,i) + & 346 w(k,j,i) ) & 347 / 2.0_wp )**2 ) * & 348 ( q(k,j,i) - sls(k,j,i) ) 875 tend(k,j,i) = tend(k,j,i) - & 876 lsec * & 877 lad_s(k,j,i) * & 878 SQRT( ( 0.5_wp * ( u(k,j,i) + & 879 u(k,j,i+1) ) & 880 )**2 + & 881 ( 0.5_wp * ( v(k,j,i) + & 882 v(k,j+1,i) ) & 883 )**2 + & 884 ( 0.5_wp * ( w(k-1,j,i) + & 885 w(k,j,i) ) & 886 )**2 & 887 ) * & 888 ( q(k,j,i) - lsc ) 349 889 ENDDO 350 890 351 891 ! 352 892 !-- sgs-tke 353 CASE ( 6 ) 354 DO k = nzb_s_inner(j,i)+1, pch_index 355 tend(k,j,i) = tend(k,j,i) - & 356 2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * & 357 SQRT( ( ( u(k,j,i) + & 358 u(k,j,i+1) ) & 359 / 2.0_wp )**2 + & 360 ( ( v(k,j,i) + & 361 v(k,j+1,i) ) & 362 / 2.0_wp )**2 + & 363 ( ( w(k,j,i) + & 364 w(k+1,j,i) ) & 365 / 2.0_wp )**2 ) * & 366 e(k,j,i) 367 ENDDO 893 CASE ( 6 ) 894 DO k = nzb_s_inner(j,i)+1, pch_index 895 tend(k,j,i) = tend(k,j,i) - & 896 2.0_wp * cdc * & 897 lad_s(k,j,i) * & 898 SQRT( ( 0.5_wp * ( u(k,j,i) + & 899 u(k,j,i+1) ) & 900 )**2 + & 901 ( 0.5_wp * ( v(k,j,i) + & 902 v(k,j+1,i) ) & 903 )**2 + & 904 ( 0.5_wp * ( w(k,j,i) + & 905 w(k+1,j,i) ) & 906 )**2 & 907 ) * & 908 e(k,j,i) 909 ENDDO 368 910 369 911 CASE DEFAULT -
palm/trunk/SOURCE/prognostic_equations.f90
r1410 r1484 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Changes due to new module structure of the plant canopy model: 23 ! parameters cthf and plant_canopy moved to module plant_canopy_model_mod. 24 ! Removed double-listing of use_upstream_for_tke in ONLY-list of module 25 ! control_parameters 23 26 ! 24 27 ! Former revisions: … … 163 166 USE control_parameters, & 164 167 ONLY: call_microphysics_at_all_substeps, cloud_physics, & 165 constant_diffusion, cthf, dp_external,&168 constant_diffusion, dp_external, & 166 169 dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d, humidity, & 167 170 icloud_scheme, inflow_l, intermediate_timestep_count, & 168 171 intermediate_timestep_count_max, large_scale_forcing, & 169 172 large_scale_subsidence, neutral, nudging, ocean, outflow_l, & 170 outflow_s, passive_scalar, p lant_canopy, precipitation,&173 outflow_s, passive_scalar, precipitation, & 171 174 prho_reference, prho_reference, prho_reference, pt_reference, & 172 175 pt_reference, pt_reference, radiation, scalar_advec, & 173 176 scalar_advec, simulated_time, sloping_surface, timestep_scheme, & 174 177 tsc, use_subsidence_tendencies, use_upstream_for_tke, & 175 use_upstream_for_tke, use_upstream_for_tke, wall_heatflux,&178 wall_heatflux, & 176 179 wall_nrflux, wall_qflux, wall_qflux, wall_qflux, wall_qrflux, & 177 180 wall_salinityflux, ws_scheme_mom, ws_scheme_sca … … 257 260 258 261 USE plant_canopy_model_mod, & 259 ONLY: plant_canopy_model262 ONLY: cthf, plant_canopy, plant_canopy_model 260 263 261 264 USE production_e_mod, & -
palm/trunk/SOURCE/read_var_list.f90
r1323 r1484 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Changes in the course of the canopy-model modularization: 23 ! parameters alpha_lad, beta_lad, lai_beta added, 24 ! module plant_canopy_model_mod added, 25 ! drag_coefficient, leaf_surface_concentration and scalar_exchange_coefficient 26 ! renamed to canopy_drag_coeff, leaf_surface_conc and leaf_scalar_exch_coeff 23 27 ! 24 28 ! Former revisions: … … 117 121 118 122 USE arrays_3d, & 119 ONLY: inflow_damping_factor, lad, mean_inflow_profiles, pt_init,&123 ONLY: inflow_damping_factor, mean_inflow_profiles, pt_init, & 120 124 q_init, ref_state, sa_init, u_init, ug, v_init, vg 121 125 … … 141 145 142 146 USE pegrid 147 148 USE plant_canopy_model_mod, & 149 ONLY: alpha_lad, beta_lad, canopy_drag_coeff, canopy_mode, cthf, lad, & 150 lad_surface, lad_vertical_gradient, lad_vertical_gradient_level,& 151 lad_vertical_gradient_level_ind, & 152 lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, pch_index, & 153 plant_canopy 143 154 144 155 USE statistics, & … … 242 253 CASE ( 'advected_distance_y' ) 243 254 READ ( 13 ) advected_distance_y 255 CASE ( 'alpha_lad' ) 256 READ ( 13 ) alpha_lad 244 257 CASE ( 'alpha_surface' ) 245 258 READ ( 13 ) alpha_surface … … 282 295 CASE ( 'bc_uv_t' ) 283 296 READ ( 13 ) bc_uv_t 297 CASE ( 'beta_lad' ) 298 READ ( 13 ) beta_lad 284 299 CASE ( 'bottom_salinityflux' ) 285 300 READ ( 13 ) bottom_salinityflux … … 296 311 CASE ( 'call_psolver_at_all_substeps' ) 297 312 READ ( 13 ) call_psolver_at_all_substeps 313 CASE ( 'canopy_drag_coeff' ) 314 READ ( 13 ) canopy_drag_coeff 298 315 CASE ( 'canopy_mode' ) 299 316 READ ( 13 ) canopy_mode … … 354 371 CASE ( 'dpdxy' ) 355 372 READ ( 13 ) dpdxy 356 CASE ( 'drag_coefficient' )357 READ ( 13 ) drag_coefficient358 373 CASE ( 'drizzle' ) 359 374 READ ( 13 ) drizzle … … 419 434 CASE ( 'lad_vertical_gradient_level_in' ) 420 435 READ ( 13 ) lad_vertical_gradient_level_ind 436 CASE ( 'lai_beta' ) 437 READ ( 13 ) lai_beta 421 438 CASE ( 'large_scale_forcing' ) 422 439 READ ( 13 ) large_scale_forcing 423 440 CASE ( 'large_scale_subsidence' ) 424 441 READ ( 13 ) large_scale_subsidence 425 CASE ( 'leaf_surface_concentration' ) 426 READ ( 13 ) leaf_surface_concentration 442 CASE ( 'leaf_scalar_exch_coeff' ) 443 READ ( 13 ) leaf_scalar_exch_coeff 444 CASE ( 'leaf_surface_conc' ) 445 READ ( 13 ) leaf_surface_conc 427 446 CASE ( 'limiter_sedimentation' ) 428 447 READ ( 13 ) limiter_sedimentation … … 558 577 CASE ( 'scalar_advec' ) 559 578 READ ( 13 ) scalar_advec 560 CASE ( 'scalar_exchange_coefficient' )561 READ ( 13 ) scalar_exchange_coefficient562 579 CASE ( 'simulated_time' ) 563 580 READ ( 13 ) simulated_time … … 731 748 732 749 USE arrays_3d, & 733 ONLY: inflow_damping_factor, lad, mean_inflow_profiles, pt_init,&750 ONLY: inflow_damping_factor, mean_inflow_profiles, pt_init, & 734 751 q_init, ref_state, sa_init, u_init, ug, v_init, vg 735 752 -
palm/trunk/SOURCE/timestep.f90
r1343 r1484 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Changes due to new module structure of the plant canopy model: 23 ! calculations and parameters related to the plant canopy model removed 24 ! (the limitation of the canopy drag, i.e. that the canopy drag itself should 25 ! not change the sign of the velocity components, is now assured for in the 26 ! calculation of the canopy tendency terms in subroutine plant_canopy_model) 23 27 ! 24 28 ! Former revisions: … … 76 80 77 81 USE arrays_3d, & 78 ONLY: cdc, dzu, dzw, kh, km, lad_u, lad_v, lad_w, u, v, w82 ONLY: dzu, dzw, kh, km, u, v, w 79 83 80 84 USE cloud_parameters, & … … 83 87 USE control_parameters, & 84 88 ONLY: cfl_factor, coupling_mode, dt_3d, dt_fixed, dt_max, & 85 galilei_transformation, old_dt, plant_canopy, message_string,&89 galilei_transformation, old_dt, message_string, & 86 90 stop_dt, terminate_coupled, terminate_coupled_remote, & 87 91 timestep_reason, u_gtrans, use_ug_for_galilei_tr, v_gtrans … … 115 119 REAL(wp) :: dt_diff !: 116 120 REAL(wp) :: dt_diff_l !: 117 REAL(wp) :: dt_plant_canopy !:118 REAL(wp) :: dt_plant_canopy_l !:119 REAL(wp) :: dt_plant_canopy_u !:120 REAL(wp) :: dt_plant_canopy_v !:121 REAL(wp) :: dt_plant_canopy_w !:122 121 REAL(wp) :: dt_u !: 123 122 REAL(wp) :: dt_u_l !: … … 413 412 414 413 ! 415 !-- Additional timestep criterion with plant canopies:416 !-- it is not allowed to extract more than the available momentum417 IF ( plant_canopy ) THEN418 419 dt_plant_canopy_l = 0.0_wp420 DO i = nxl, nxr421 DO j = nys, nyn422 DO k = nzb+1, nzt423 dt_plant_canopy_u = cdc(k,j,i) * lad_u(k,j,i) * &424 SQRT( u(k,j,i)**2 + &425 ( ( v(k,j,i-1) + &426 v(k,j,i) + &427 v(k,j+1,i) + &428 v(k,j+1,i-1) ) &429 / 4.0_wp )**2 + &430 ( ( w(k-1,j,i-1) + &431 w(k-1,j,i) + &432 w(k,j,i-1) + &433 w(k,j,i) ) &434 / 4.0_wp )**2 )435 IF ( dt_plant_canopy_u > dt_plant_canopy_l ) THEN436 dt_plant_canopy_l = dt_plant_canopy_u437 ENDIF438 dt_plant_canopy_v = cdc(k,j,i) * lad_v(k,j,i) * &439 SQRT( ( ( u(k,j-1,i) + &440 u(k,j-1,i+1) + &441 u(k,j,i) + &442 u(k,j,i+1) ) &443 / 4.0_wp )**2 + &444 v(k,j,i)**2 + &445 ( ( w(k-1,j-1,i) + &446 w(k-1,j,i) + &447 w(k,j-1,i) + &448 w(k,j,i) ) &449 / 4.0_wp )**2 )450 IF ( dt_plant_canopy_v > dt_plant_canopy_l ) THEN451 dt_plant_canopy_l = dt_plant_canopy_v452 ENDIF453 dt_plant_canopy_w = cdc(k,j,i) * lad_w(k,j,i) * &454 SQRT( ( ( u(k,j,i) + &455 u(k,j,i+1) + &456 u(k+1,j,i) + &457 u(k+1,j,i+1) ) &458 / 4.0_wp )**2 + &459 ( ( v(k,j,i) + &460 v(k,j+1,i) + &461 v(k+1,j,i) + &462 v(k+1,j+1,i) ) &463 / 4.0_wp )**2 + &464 w(k,j,i)**2 )465 IF ( dt_plant_canopy_w > dt_plant_canopy_l ) THEN466 dt_plant_canopy_l = dt_plant_canopy_w467 ENDIF468 ENDDO469 ENDDO470 ENDDO471 472 IF ( dt_plant_canopy_l > 0.0_wp ) THEN473 !474 !-- Invert dt_plant_canopy_l and apply a security timestep factor 0.1475 dt_plant_canopy_l = 0.1_wp / dt_plant_canopy_l476 ELSE477 !478 !-- In case of inhomogeneous plant canopy, some processors may have no479 !-- canopy at all. Then use dt_max as dummy instead.480 dt_plant_canopy_l = dt_max481 ENDIF482 483 !484 !-- Determine the global minumum485 #if defined( __parallel )486 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )487 CALL MPI_ALLREDUCE( dt_plant_canopy_l, dt_plant_canopy, 1, MPI_REAL, &488 MPI_MIN, comm2d, ierr )489 #else490 dt_plant_canopy = dt_plant_canopy_l491 #endif492 493 ELSE494 !495 !-- Use dt_diff as dummy value to avoid extra IF branches further below496 dt_plant_canopy = dt_diff497 498 ENDIF499 500 !501 414 !-- The time step is the minimum of the 3-4 components and the diffusion time 502 415 !-- step minus a reduction (cfl_factor) to be on the safe side. 503 416 !-- The time step must not exceed the maximum allowed value. 504 dt_3d = cfl_factor * MIN( dt_diff, dt_ plant_canopy, dt_u, dt_v, dt_w, &417 dt_3d = cfl_factor * MIN( dt_diff, dt_u, dt_v, dt_w, & 505 418 dt_precipitation ) 506 419 dt_3d = MIN( dt_3d, dt_max ) … … 508 421 ! 509 422 !-- Remember the restricting time step criterion for later output. 510 IF ( MIN( dt_u, dt_v, dt_w ) < MIN( dt_diff, dt_plant_canopy )) THEN423 IF ( MIN( dt_u, dt_v, dt_w ) < dt_diff ) THEN 511 424 timestep_reason = 'A' 512 ELSEIF ( dt_plant_canopy < dt_diff ) THEN513 timestep_reason = 'C'514 425 ELSE 515 426 timestep_reason = 'D' … … 528 439 '&dt_w = ', dt_w, ' s', & 529 440 '&dt_diff = ', dt_diff, ' s', & 530 '&dt_plant_canopy = ', dt_plant_canopy, ' s', &531 441 '&u_max = ', u_max, ' m/s k=', u_max_ijk(1), & 532 442 ' j=', u_max_ijk(2), ' i=', u_max_ijk(3), & -
palm/trunk/SOURCE/user_init_plant_canopy.f90
r1354 r1484 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Changes in the course of the canopy-model modularization: 23 ! module plant_canopy_model_mod added, 24 ! definition of array cdc (canopy drag coefficient) removed, since it is now 25 ! defined purely as a single constant value (see module plant_canopy_model) 23 26 ! 24 27 ! Former revisions: … … 57 60 58 61 USE kinds 62 63 USE plant_canopy_model_mod 59 64 60 65 USE user … … 62 67 IMPLICIT NONE 63 68 64 INTEGER(iwp) :: i !: 65 INTEGER(iwp) :: j !: 69 INTEGER(iwp) :: i !: running index 70 INTEGER(iwp) :: j !: running index 66 71 67 72 ! … … 69 74 70 75 ! 71 !-- Set the 3D-array s lad_s and cdcfor user defined canopies76 !-- Set the 3D-array lad_s for user defined canopies 72 77 SELECT CASE ( TRIM( canopy_mode ) ) 73 78 … … 78 83 CASE ( 'user_defined_canopy_1' ) 79 84 ! 80 !-- Here the user can define his own topography. 81 !-- The following lines contain an example in that the 82 !-- plant canopy extends only over the second half of the 83 !-- model domain along x 85 !-- Here the user can define his own forest topography. 86 !-- The following lines contain an example, where the plant canopy extends 87 !-- only over the second half of the model domain along x. 88 !-- Attention: DO-loops have to include the ghost points (nxlg-nxrg, 89 !-- nysg-nyng), because no exchange of ghost point information is intended, 90 !-- in order to minimize communication between CPUs 84 91 ! DO i = nxlg, nxrg 85 92 ! IF ( i >= INT(nx+1/2) ) THEN 86 93 ! DO j = nysg, nyng 87 94 ! lad_s(:,j,i) = lad(:) 88 ! cdc(:,j,i) = drag_coefficient89 95 ! ENDDO 90 96 ! ELSE 91 97 ! lad_s(:,:,i) = 0.0_wp 92 ! cdc(:,:,i) = 0.0_wp93 98 ! ENDIF 94 99 ! ENDDO 100 ! 95 101 !-- After definition, please 96 102 !-- remove the following three lines! -
palm/trunk/SOURCE/write_var_list.f90
r1353 r1484 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Changes in the course of the canopy-model modularization: 23 ! parameters alpha_lad, beta_lad, lai_beta added, 24 ! module plant_canopy_model_mod added, 25 ! drag_coefficient, leaf_surface_concentration and scalar_exchange_coefficient 26 ! renamed to canopy_drag_coeff, leaf_surface_conc and leaf_scalar_exch_coeff 23 27 ! 24 28 ! Former revisions: … … 104 108 105 109 USE arrays_3d, & 106 ONLY: inflow_damping_factor, lad, mean_inflow_profiles, pt_init,&110 ONLY: inflow_damping_factor, mean_inflow_profiles, pt_init, & 107 111 q_init, ref_state, sa_init, u_init, ug, v_init, vg 108 112 … … 128 132 129 133 USE pegrid 134 135 USE plant_canopy_model_mod, & 136 ONLY: alpha_lad, beta_lad, canopy_drag_coeff, canopy_mode, cthf, lad, & 137 lad_surface, lad_vertical_gradient, lad_vertical_gradient_level,& 138 lad_vertical_gradient_level_ind, & 139 lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, pch_index, & 140 plant_canopy 130 141 131 142 USE statistics, & … … 163 174 WRITE ( 14 ) 'advected_distance_y ' 164 175 WRITE ( 14 ) advected_distance_y 176 WRITE ( 14 ) 'alpha_lad ' 177 WRITE ( 14 ) alpha_lad 165 178 WRITE ( 14 ) 'alpha_surface ' 166 179 WRITE ( 14 ) alpha_surface … … 203 216 WRITE ( 14 ) 'bc_uv_t ' 204 217 WRITE ( 14 ) bc_uv_t 218 WRITE ( 14 ) 'beta_lad ' 219 WRITE ( 14 ) beta_lad 205 220 WRITE ( 14 ) 'bottom_salinityflux ' 206 221 WRITE ( 14 ) bottom_salinityflux … … 217 232 WRITE ( 14 ) 'call_psolver_at_all_substeps ' 218 233 WRITE ( 14 ) call_psolver_at_all_substeps 234 WRITE ( 14 ) 'canopy_drag_coeff ' 235 WRITE ( 14 ) canopy_drag_coeff 219 236 WRITE ( 14 ) 'canopy_mode ' 220 237 WRITE ( 14 ) canopy_mode … … 275 292 WRITE ( 14 ) 'dpdxy ' 276 293 WRITE ( 14 ) dpdxy 277 WRITE ( 14 ) 'drag_coefficient '278 WRITE ( 14 ) drag_coefficient279 294 WRITE ( 14 ) 'drizzle ' 280 295 WRITE ( 14 ) drizzle … … 339 354 WRITE ( 14 ) 'lad_vertical_gradient_level_in' 340 355 WRITE ( 14 ) lad_vertical_gradient_level_ind 356 WRITE ( 14 ) 'lai_beta ' 357 WRITE ( 14 ) lai_beta 341 358 WRITE ( 14 ) 'large_scale_forcing ' 342 359 WRITE ( 14 ) large_scale_forcing 343 360 WRITE ( 14 ) 'large_scale_subsidence ' 344 361 WRITE ( 14 ) large_scale_subsidence 345 WRITE ( 14 ) 'leaf_surface_concentration ' 346 WRITE ( 14 ) leaf_surface_concentration 362 WRITE ( 14 ) 'leaf_scalar_exch_coeff ' 363 WRITE ( 14 ) leaf_scalar_exch_coeff 364 WRITE ( 14 ) 'leaf_surface_conc ' 365 WRITE ( 14 ) leaf_surface_conc 347 366 WRITE ( 14 ) 'limiter_sedimentation ' 348 367 WRITE ( 14 ) limiter_sedimentation … … 475 494 WRITE ( 14 ) 'scalar_advec ' 476 495 WRITE ( 14 ) scalar_advec 477 WRITE ( 14 ) 'scalar_exchange_coefficient '478 WRITE ( 14 ) scalar_exchange_coefficient479 496 WRITE ( 14 ) 'simulated_time ' 480 497 WRITE ( 14 ) simulated_time
Note: See TracChangeset
for help on using the changeset viewer.
