Changeset 3449 for palm/trunk/SOURCE/plant_canopy_model_mod.f90
- Timestamp:
- Oct 29, 2018 7:36:56 PM (5 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/plant_canopy_model_mod.f90 (modified) (21 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/plant_canopy_model_mod.f90
r3337 r3449 16 16 ! 17 17 ! Copyright 1997-2018 Leibniz Universitaet Hannover 18 ! Copyright 2018 Institute of Computer Science of the 19 ! Czech Academy of Sciences, Prague 18 20 !------------------------------------------------------------------------------! 19 21 ! … … 25 27 ! ----------------- 26 28 ! $Id$ 29 ! Add calculation of transpiration for resolved plant canopy (trees, shrubs) 30 ! (V. Fuka, MFF UK Prague, J.Resler, ICS AS, Prague) 31 ! 32 ! Fix reading plant canopy over buildings 33 ! 34 ! 3337 2018-10-12 15:17:09Z kanani 27 35 ! Fix reading plant canopy over buildings 28 36 ! … … 214 222 215 223 USE arrays_3d, & 216 ONLY: dzu, dzw, e, q, s, tend, u, v, w, zu, zw 224 ONLY: dzu, dzw, e, exner, hyp, pt, q, s, tend, u, v, w, zu, zw 225 226 USE basic_constants_and_equations_mod, & 227 ONLY: c_p, degc_to_k, l_v, lv_d_cp, r_d, rd_d_rv 228 229 USE control_parameters, & 230 ONLY: humidity 217 231 218 232 USE indices, & … … 222 236 USE kinds 223 237 238 USE pegrid 239 224 240 USE surface_mod, & 225 241 ONLY: get_topography_top_index_ji … … 229 245 230 246 231 CHARACTER (LEN=30) :: canopy_mode = 'block' !< canopy coverage 232 233 INTEGER(iwp) :: pch_index = 0 !< plant canopy height/top index 234 INTEGER(iwp) :: lad_vertical_gradient_level_ind(10) = -9999 !< lad-profile levels (index) 235 236 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pch_index_ji !< local plant canopy top 237 238 LOGICAL :: calc_beta_lad_profile = .FALSE. !< switch for calc. of lad from beta func. 247 CHARACTER (LEN=30) :: canopy_mode = 'block' !< canopy coverage 248 LOGICAL :: plant_canopy_transpiration = .FALSE. !< flag to switch calculation of transpiration and corresponding latent heat 249 !< for resolved plant canopy inside radiation model 250 !< (calls subroutine pcm_calc_transpiration_rate from module plant_canopy_mod) 251 252 INTEGER(iwp) :: pch_index = 0 !< plant canopy height/top index 253 INTEGER(iwp) :: lad_vertical_gradient_level_ind(10) = -9999 !< lad-profile levels (index) 254 255 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pch_index_ji !< local plant canopy top 256 257 LOGICAL :: calc_beta_lad_profile = .FALSE. !< switch for calc. of lad from beta func. 239 258 240 259 REAL(wp) :: alpha_lad = 9999999.9_wp !< coefficient for lad calculation … … 261 280 REAL(wp), DIMENSION(:), ALLOCATABLE :: pre_lad !< preliminary lad 262 281 263 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: cum_lai_hf !< cumulative lai for heatflux calc.264 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: lad_s !< lad on scalar-grid265 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pc_heating_rate !< plant canopy heating rate282 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: cum_lai_hf !< cumulative lai for heatflux calc. 283 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: lad_s !< lad on scalar-grid 284 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pc_heating_rate !< plant canopy heating rate 266 285 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pc_transpiration_rate !< plant canopy transpiration rate 286 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pc_latent_rate !< plant canopy latent heating rate 267 287 268 288 SAVE … … 273 293 ! 274 294 !-- Public functions 275 PUBLIC pcm_check_data_output, pcm_check_parameters, pcm_data_output_3d, & 276 pcm_define_netcdf_grid, pcm_header, pcm_init, pcm_parin, pcm_tendency 295 PUBLIC pcm_calc_transpiration_rate, pcm_check_data_output, & 296 pcm_check_parameters, pcm_data_output_3d, pcm_define_netcdf_grid, & 297 pcm_header, pcm_init, pcm_parin, pcm_tendency 277 298 278 299 ! 279 300 !-- Public variables and constants 280 PUBLIC pc_heating_rate, pc_transpiration_rate, canopy_mode, cthf, dt_plant_canopy, lad, lad_s, & 281 pch_index 282 301 PUBLIC pc_heating_rate, pc_transpiration_rate, pc_latent_rate, & 302 canopy_mode, cthf, dt_plant_canopy, lad, lad_s, pch_index, & 303 plant_canopy_transpiration 304 305 INTERFACE pcm_calc_transpiration_rate 306 MODULE PROCEDURE pcm_calc_transpiration_rate 307 END INTERFACE pcm_calc_transpiration_rate 283 308 284 309 INTERFACE pcm_check_data_output … … 323 348 324 349 350 351 !------------------------------------------------------------------------------! 352 ! Description: 353 ! ------------ 354 !> Calculation of the plant canopy transpiration rate based on the Jarvis-Stewart 355 !> with parametrizations described in Daudet et al. (1999; Agricult. and Forest 356 !> Meteorol. 97) and Ngao, Adam and Saudreau (2017; Agricult. and Forest Meteorol 357 !> 237-238). Model functions f1-f4 were adapted from Stewart (1998; Agric. 358 !> and Forest. Meteorol. 43) instead, because they are valid for broader intervals 359 !> of values. Funcion f4 used in form present in van Wijk et al. (1998; 360 !> Tree Physiology 20). 361 !> 362 !> This subroutine is called from subroutine radiation_interaction 363 !> after the calculation of radiation in plant canopy boxes. 364 !> (arrays pcbinsw and pcbinlw). 365 !> 366 !------------------------------------------------------------------------------! 367 SUBROUTINE pcm_calc_transpiration_rate(i, j, k, kk, pcbsw, pcblw, pcbtr, pcblh) 368 369 USE control_parameters, & 370 ONLY: dz 371 372 USE grid_variables, & 373 ONLY: dx, dy 374 375 IMPLICIT NONE 376 !-- input parameters 377 INTEGER(iwp), INTENT(IN) :: i, j, k, kk !< indices of the pc gridbox 378 REAL(wp), INTENT(IN) :: pcbsw !< sw radiation in gridbox (W) 379 REAL(wp), INTENT(IN) :: pcblw !< lw radiation in gridbox (W) 380 REAL(wp), INTENT(OUT) :: pcbtr !< transpiration rate dq/dt (kg/kg/s) 381 REAL(wp), INTENT(OUT) :: pcblh !< latent heat from transpiration dT/dt (K/s) 382 383 !-- variables and parameters for calculation of transpiration rate 384 REAL(wp) :: sat_press, sat_press_d, temp, v_lad 385 REAL(wp) :: d_fact, g_b, g_s, wind_speed, evapor_rate 386 REAL(wp) :: f1, f2, f3, f4, vpd, rswc, e_eq, e_imp, rad 387 REAL(wp), PARAMETER :: gama_psychr = 66 !< psychrometric constant (Pa/K) 388 REAL(wp), PARAMETER :: g_s_max = 0.01 !< maximum stomatal conductivity (m/s) 389 REAL(wp), PARAMETER :: m_soil = 0.4_wp !< soil water content (needs to adjust or take from LSM) 390 REAL(wp), PARAMETER :: m_wilt = 0.01_wp !< wilting point soil water content (needs to adjust or take from LSM) 391 REAL(wp), PARAMETER :: m_sat = 0.51_wp !< saturation soil water content (needs to adjust or take from LSM) 392 REAL(wp), PARAMETER :: t2_min = 0.0_wp !< minimal temperature for calculation of f2 393 REAL(wp), PARAMETER :: t2_max = 40.0_wp !< maximal temperature for calculation of f2 394 395 396 !-- Temperature (deg C) 397 temp = pt(k,j,i) * exner(k) - degc_to_k 398 !-- Coefficient for conversion of radiation to grid to radiation to unit leaves surface 399 v_lad = 1.0_wp / ( MAX(lad_s(kk,j,i), 1.0e-10) * dx * dy * dz(1) ) 400 !-- Magnus formula for the saturation pressure (see Ngao, Adam and Saudreau (2017) eq. 1) 401 !-- There are updated formulas available, kept consistent with the rest of the parametrization 402 sat_press = 610.8_wp * exp(17.27_wp * temp/(temp + 237.3_wp)) 403 !-- Saturation pressure derivative (derivative of the above) 404 sat_press_d = sat_press * 17.27_wp * 237.3_wp / (temp + 237.3_wp)**2 405 !-- Wind speed 406 wind_speed = SQRT( ((u(k,j,i) + u(k,j,i+1))/2.0_wp )**2 + & 407 ((v(k,j,i) + v(k,j,i+1))/2.0_wp )**2 + & 408 ((w(k,j,i) + w(k,j,i+1))/2.0_wp )**2 ) 409 !-- Aerodynamic conductivity (Daudet et al. (1999) eq. 14 410 g_b = 0.01_wp * wind_speed + 0.0071_wp 411 !-- Radiation flux per leaf surface unit 412 rad = pcbsw * v_lad 413 !-- First function for calculation of stomatal conductivity (radiation dependency) 414 !-- Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 17 415 f1 = rad * (1000._wp+42.1_wp) / 1000._wp / (rad+42.1_wp) 416 !-- Second function for calculation of stomatal conductivity (temperature dependency) 417 !-- Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 21 418 f2 = MAX(t2_min, (temp-t2_min) * (t2_max-temp)**((t2_max-16.9_wp)/(16.9_wp-t2_min)) / & 419 ((16.9_wp-t2_min) * (t2_max-16.9_wp)**((t2_max-16.9_wp)/(16.9_wp-t2_min))) ) 420 !-- Water pressure deficit 421 !-- Ngao, Adam and Saudreau (2017) eq. 6 but with water vapour partial pressure 422 vpd = max( sat_press - q(k,j,i) * hyp(k) / rd_d_rv, 0._wp ) 423 !-- Third function for calculation of stomatal conductivity (water pressure deficit dependency) 424 !-- Ngao, Adam and Saudreau (2017) Table 1, limited from below according to Stewart (1988) 425 !-- The coefficients of the linear dependence should better correspond to broad-leaved trees 426 !-- than the coefficients from Stewart (1988) which correspond to conifer trees. 427 vpd = MIN(MAX(vpd,770.0_wp),3820.0_wp) 428 f3 = -2e-4_wp * vpd + 1.154_wp 429 !-- Fourth function for calculation of stomatal conductivity (soil moisture dependency) 430 !-- Residual soil water content 431 !-- van Wijk et al. (1998; Tree Physiology 20) eq. 7 432 !-- TODO - over LSM surface might be calculated from LSM parameters 433 rswc = ( m_sat - m_soil ) / ( m_sat - m_wilt ) 434 !-- van Wijk et al. (1998; Tree Physiology 20) eq. 5-6 (it is a reformulation of eq. 22-23 of Stewart(1988)) 435 f4 = MAX(0._wp, MIN(1.0_wp - 0.041_wp * EXP(3.2_wp * rswc), 1.0_wp - 0.041_wp)) 436 !-- Stomatal conductivity 437 !-- Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 12 438 !-- (notation according to Ngao, Adam and Saudreau (2017) and others) 439 g_s = g_s_max * f1 * f2 * f3 * f4 + 1.0e-10_wp 440 !-- Decoupling factor 441 !-- Daudet et al. (1999) eq. 6 442 d_fact = (sat_press_d / gama_psychr + 2._wp ) / & 443 (sat_press_d / gama_psychr + 2._wp + 2 * g_b / g_s ) 444 !-- Equilibrium evaporation rate 445 !-- Daudet et al. (1999) eq. 4 446 e_eq = (pcbsw + pcblw) * v_lad * sat_press_d / & 447 gama_psychr /( sat_press_d / gama_psychr + 2.0_wp ) / l_v 448 !-- Imposed evaporation rate 449 !-- Daudet et al. (1999) eq. 5 450 e_imp = r_d * pt(k,j,i) * exner(k) / hyp(k) * c_p * g_s * vpd / gama_psychr / l_v 451 !-- Evaporation rate 452 !-- Daudet et al. (1999) eq. 3 453 !-- (evaporation rate is limited to non-negative values) 454 evapor_rate = MAX(d_fact * e_eq + ( 1.0_wp - d_fact ) * e_imp, 0.0_wp) 455 !-- Conversion of evaporation rate to q tendency in gridbox 456 !-- dq/dt = E * LAD * V_g / (rho_air * V_g) 457 pcbtr = evapor_rate * r_d * pt(k,j,i) * exner(k) * lad_s(kk,j,i) / hyp(k) !-- = dq/dt 458 !-- latent heat from evaporation 459 pcblh = pcbtr * lv_d_cp !-- = - dT/dt 460 461 END SUBROUTINE pcm_calc_transpiration_rate 462 463 325 464 !------------------------------------------------------------------------------! 326 465 ! Description: … … 352 491 CASE ( 'pcm_transpirationrate' ) 353 492 unit = 'kg kg-1 s-1' 493 494 CASE ( 'pcm_latentrate' ) 495 unit = 'K s-1' 496 497 CASE ( 'pcm_bowenratio' ) 498 unit = 'K s-1' 354 499 355 500 CASE ( 'pcm_lad' ) … … 503 648 ENDDO 504 649 ENDIF 505 650 651 CASE ( 'pcm_latentrate' ) 652 IF ( av == 0 ) THEN 653 DO i = nxl, nxr 654 DO j = nys, nyn 655 IF ( pch_index_ji(j,i) /= 0 ) THEN 656 k_topo = get_topography_top_index_ji( j, i, 's' ) 657 DO k = k_topo, k_topo + pch_index_ji(j,i) 658 local_pf(i,j,k) = pc_latent_rate(k-k_topo,j,i) 659 ENDDO 660 ENDIF 661 ENDDO 662 ENDDO 663 ENDIF 664 665 CASE ( 'pcm_bowenratio' ) 666 IF ( av == 0 ) THEN 667 DO i = nxl, nxr 668 DO j = nys, nyn 669 IF ( pch_index_ji(j,i) /= 0 ) THEN 670 k_topo = get_topography_top_index_ji( j, i, 's' ) 671 DO k = k_topo, k_topo + pch_index_ji(j,i) 672 IF ( pc_latent_rate(k-k_topo,j,i) /= 0._wp ) THEN 673 local_pf(i,j,k) = pc_heating_rate(k-k_topo,j,i) / & 674 pc_latent_rate(k-k_topo,j,i) 675 ENDIF 676 ENDDO 677 ENDIF 678 ENDDO 679 ENDDO 680 ENDIF 681 506 682 CASE ( 'pcm_lad' ) 507 683 IF ( av == 0 ) THEN … … 550 726 SELECT CASE ( TRIM( var ) ) 551 727 552 CASE ( 'pcm_heatrate', 'pcm_lad', 'pcm_transpirationrate' )728 CASE ( 'pcm_heatrate', 'pcm_lad', 'pcm_transpirationrate', 'pcm_latentrate', 'pcm_bowenratio') 553 729 grid_x = 'x' 554 730 grid_y = 'y' … … 691 867 692 868 USE control_parameters, & 693 ONLY: humidity,message_string, ocean_mode, urban_surface869 ONLY: message_string, ocean_mode, urban_surface 694 870 695 871 USE netcdf_data_input_mod, & … … 705 881 INTEGER(iwp) :: k !< running index 706 882 INTEGER(iwp) :: m !< running index 883 INTEGER(iwp) :: pch_index_l 707 884 708 885 REAL(wp) :: int_bpdf !< vertical integral for lad-profile construction … … 903 1080 CALL exchange_horiz_2d_int( pch_index_ji, nys, nyn, nxl, nxr, nbgp ) 904 1081 1082 !-- Calculate global pch_index value (index of top of plant canopy from ground) 1083 pch_index = MAXVAL(pch_index_ji) 1084 !-- Exchange pch_index from all processors 1085 #if defined( __parallel ) 1086 CALL MPI_Allreduce(MPI_IN_PLACE, pch_index, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr) 1087 IF ( ierr /= 0 ) THEN 1088 WRITE (9,*) 'Error MPI_Allreduce pch_index:', ierr 1089 FLUSH(9) 1090 ENDIF 1091 #endif 1092 1093 !-- Allocation of arrays pc_heating_rate, pc_transpiration_rate and pc_latent_rate 1094 ALLOCATE( pc_heating_rate(0:pch_index,nysg:nyng,nxlg:nxrg) ) 1095 IF ( humidity ) THEN 1096 ALLOCATE( pc_transpiration_rate(0:pch_index,nysg:nyng,nxlg:nxrg) ) 1097 pc_transpiration_rate = 0.0_wp 1098 ALLOCATE( pc_latent_rate(0:pch_index,nysg:nyng,nxlg:nxrg) ) 1099 pc_latent_rate = 0.0_wp 1100 ENDIF 1101 905 1102 ! 906 1103 !-- Initialization of the canopy heat source distribution due to heating … … 915 1112 !-- 73â96). This approach has been applied e.g. by Shaw & Schumann (1992; 916 1113 !-- Bound.-Layer Meteorol. 61, 47â64). 917 !-- When using the urban surface model (USM), canopy heating (pc_heating_rate) 918 !-- by radiation is calculated in the USM. 919 IF ( cthf /= 0.0_wp .AND. .NOT. urban_surface ) THEN 920 921 ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 922 pc_heating_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1114 !-- When using the radiation_interactions, canopy heating (pc_heating_rate) 1115 !-- and plant canopy transpiration (pc_transpiration_rate, pc_latent_rate) 1116 !-- are calculated in the RTM after the calculation of radiation. 1117 !-- We cannot use variable radiation_interactions here to determine the situation 1118 !-- as it is assigned in init_3d_model after the call of pcm_init. 1119 IF ( cthf /= 0.0_wp ) THEN 1120 1121 ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 923 1122 ! 924 1123 !-- Piecewise calculation of the cumulative leaf area index by vertical … … 996 1195 997 1196 ENDIF 998 !999 !-- Allocate transpiration rate1000 IF ( humidity ) &1001 ALLOCATE( pc_transpiration_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )1002 1003 1197 1004 1198 … … 1028 1222 lai_beta, & 1029 1223 leaf_scalar_exch_coeff, & 1030 leaf_surface_conc, pch_index 1224 leaf_surface_conc, pch_index, & 1225 plant_canopy_transpiration 1031 1226 1032 1227 NAMELIST /canopy_par/ alpha_lad, beta_lad, canopy_drag_coeff, & … … 1037 1232 lai_beta, & 1038 1233 leaf_scalar_exch_coeff, & 1039 leaf_surface_conc, pch_index 1234 leaf_surface_conc, pch_index, & 1235 plant_canopy_transpiration 1040 1236 1041 1237 line = ' ' … … 1434 1630 !-- potential temperature 1435 1631 CASE ( 4 ) 1436 DO i = nxl, nxr 1437 DO j = nys, nyn 1438 ! 1439 !-- Determine topography-top index on scalar-grid 1440 k_wall = get_topography_top_index_ji( j, i, 's' ) 1441 1442 DO k = k_wall+1, k_wall + pch_index_ji(j,i) 1443 1444 kk = k - k_wall !- lad arrays are defined flat 1445 tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) 1632 IF ( humidity ) THEN 1633 DO i = nxl, nxr 1634 DO j = nys, nyn 1635 !-- Determine topography-top index on scalar-grid 1636 k_wall = get_topography_top_index_ji( j, i, 's' ) 1637 DO k = k_wall+1, k_wall + pch_index_ji(j,i) 1638 kk = k - k_wall !- lad arrays are defined flat 1639 tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) - pc_latent_rate(kk,j,i) 1640 ENDDO 1446 1641 ENDDO 1447 1642 ENDDO 1448 ENDDO 1643 ELSE 1644 DO i = nxl, nxr 1645 DO j = nys, nyn 1646 !-- Determine topography-top index on scalar-grid 1647 k_wall = get_topography_top_index_ji( j, i, 's' ) 1648 DO k = k_wall+1, k_wall + pch_index_ji(j,i) 1649 kk = k - k_wall !- lad arrays are defined flat 1650 tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) 1651 ENDDO 1652 ENDDO 1653 ENDDO 1654 ENDIF 1449 1655 1450 1656 ! … … 1460 1666 1461 1667 kk = k - k_wall !- lad arrays are defined flat 1462 pc_transpiration_rate(kk,j,i) = - lsec & 1463 * lad_s(kk,j,i) * & 1464 SQRT( ( 0.5_wp * ( u(k,j,i) + & 1465 u(k,j,i+1) ) & 1466 )**2 + & 1467 ( 0.5_wp * ( v(k,j,i) + & 1468 v(k,j+1,i) ) & 1469 )**2 + & 1470 ( 0.5_wp * ( w(k-1,j,i) + & 1471 w(k,j,i) ) & 1472 )**2 & 1473 ) * & 1474 ( q(k,j,i) - lsc ) 1668 1669 IF ( .NOT. plant_canopy_transpiration ) THEN 1670 ! pc_transpiration_rate is calculated in radiation model 1671 ! in case of plant_canopy_transpiration = .T. 1672 ! to include also the dependecy to the radiation 1673 ! in the plant canopy box 1674 pc_transpiration_rate(kk,j,i) = - lsec & 1675 * lad_s(kk,j,i) * & 1676 SQRT( ( 0.5_wp * ( u(k,j,i) + & 1677 u(k,j,i+1) ) & 1678 )**2 + & 1679 ( 0.5_wp * ( v(k,j,i) + & 1680 v(k,j+1,i) ) & 1681 )**2 + & 1682 ( 0.5_wp * ( w(k-1,j,i) + & 1683 w(k,j,i) ) & 1684 )**2 & 1685 ) * & 1686 ( q(k,j,i) - lsc ) 1687 ENDIF 1475 1688 1476 1689 tend(k,j,i) = tend(k,j,i) + pc_transpiration_rate(kk,j,i) … … 1778 1991 k_wall = get_topography_top_index_ji( j, i, 's' ) 1779 1992 1993 IF ( humidity ) THEN 1994 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) 1995 kk = k - k_wall !- lad arrays are defined flat 1996 tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) - pc_latent_rate(kk,j,i) 1997 ENDDO 1998 ELSE 1999 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) 2000 kk = k - k_wall !- lad arrays are defined flat 2001 tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) 2002 ENDDO 2003 ENDIF 2004 2005 ! 2006 !-- humidity 2007 CASE ( 5 ) 2008 ! 2009 !-- Determine topography-top index on scalar grid 2010 k_wall = get_topography_top_index_ji( j, i, 's' ) 2011 1780 2012 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) 1781 2013 kk = k - k_wall !- lad arrays are defined flat 1782 tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) 1783 ENDDO 1784 1785 1786 ! 1787 !-- humidity 1788 CASE ( 5 ) 1789 ! 1790 !-- Determine topography-top index on scalar grid 1791 k_wall = get_topography_top_index_ji( j, i, 's' ) 1792 1793 DO k = k_wall + 1, k_wall + pch_index_ji(j,i) 1794 kk = k - k_wall !- lad arrays are defined flat 1795 1796 pc_transpiration_rate(kk,j,i) = - lsec & 1797 * lad_s(kk,j,i) * & 1798 SQRT( ( 0.5_wp * ( u(k,j,i) + & 1799 u(k,j,i+1) ) & 1800 )**2 + & 1801 ( 0.5_wp * ( v(k,j,i) + & 1802 v(k,j+1,i) ) & 1803 )**2 + & 1804 ( 0.5_wp * ( w(k-1,j,i) + & 1805 w(k,j,i) ) & 1806 )**2 & 1807 ) * & 1808 ( q(k,j,i) - lsc ) 2014 IF ( .NOT. plant_canopy_transpiration ) THEN 2015 ! pc_transpiration_rate is calculated in radiation model 2016 ! in case of plant_canopy_transpiration = .T. 2017 ! to include also the dependecy to the radiation 2018 ! in the plant canopy box 2019 pc_transpiration_rate(kk,j,i) = - lsec & 2020 * lad_s(kk,j,i) * & 2021 SQRT( ( 0.5_wp * ( u(k,j,i) + & 2022 u(k,j,i+1) ) & 2023 )**2 + & 2024 ( 0.5_wp * ( v(k,j,i) + & 2025 v(k,j+1,i) ) & 2026 )**2 + & 2027 ( 0.5_wp * ( w(k-1,j,i) + & 2028 w(k,j,i) ) & 2029 )**2 & 2030 ) * & 2031 ( q(k,j,i) - lsc ) 2032 ENDIF 1809 2033 1810 2034 tend(k,j,i) = tend(k,j,i) + pc_transpiration_rate(kk,j,i)
Note: See TracChangeset
for help on using the changeset viewer.
