Changeset 2504 for palm/trunk/SOURCE/land_surface_model_mod.f90
- Timestamp:
- Sep 27, 2017 10:36:13 AM (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/land_surface_model_mod.f90
r2476 r2504 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Support roots and water under pavement. Added several pavement types. 28 ! 29 ! 2476 2017-09-18 07:54:32Z maronga 27 30 ! Bugfix for last commit 28 31 ! … … 558 561 ! 559 562 !-- Pavement classes 560 CHARACTER(20), DIMENSION(0:7), PARAMETER :: pavement_type_name = (/ & 561 'user defined ', & ! 0 562 'asphalt ', & ! 1 563 'concrete ', & ! 2 564 'asphalt/concrete mix', & ! 3 565 'brick pavers ', & ! 4 566 'cobblestone pavers ', & ! 5 567 'sett pavers ', & ! 6 568 'gravel pavers ' & ! 7 563 CHARACTER(29), DIMENSION(0:15), PARAMETER :: pavement_type_name = (/ & 564 'user defined ', & ! 0 565 'asphalt/concrete mix ', & ! 1 566 'asphalt (asphalt concrete) ', & ! 2 567 'concrete (Portland concrete) ', & ! 3 568 'sett ', & ! 4 569 'paving stones ', & ! 5 570 'cobblestone ', & ! 6 571 'metal ', & ! 7 572 'wood ', & ! 8 573 'gravel ', & ! 9 574 'fine gravel ', & ! 10 575 'pebblestone ', & ! 11 576 'woodchips ', & ! 12 577 'tartan (sports) ', & ! 13 578 'artifical turf (sports) ', & ! 14 579 'clay (sports) ' & ! 15 569 580 /) 570 581 582 583 584 585 571 586 ! 572 587 !-- Water classes … … 587 602 !-- r_canopy_min, lai, c_veg, g_d z0, z0h, lambda_s_s, lambda_s_u, f_sw_in, c_surface, albedo_type, emissivity 588 603 REAL(wp), DIMENSION(0:11,1:18), PARAMETER :: vegetation_pars = RESHAPE( (/ & 589 0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, 0.005_wp, 0.5E-4_wp, 0.0_wp, 0.0_wp, 0.00_wp, 0.00_wp, 1 8.0_wp, 0.94_wp, & ! 1604 0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, 0.005_wp, 0.5E-4_wp, 0.0_wp, 0.0_wp, 0.00_wp, 0.00_wp, 17.0_wp, 0.94_wp, & ! 1 590 605 180.0_wp, 3.00_wp, 1.00_wp, 0.00_wp, 0.10_wp, 0.001_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 2.0_wp, 0.95_wp, & ! 2 591 606 110.0_wp, 2.00_wp, 1.00_wp, 0.00_wp, 0.03_wp, 0.3E-4_wp, 10.0_wp, 10.0_wp, 0.05_wp, 0.00_wp, 2.0_wp, 0.95_wp, & ! 3 … … 649 664 !-- TO BE FILLED 650 665 !-- Pavement parameters depth, z0, z0h, lambda_h, rho_c, albedo_type, emissivity 651 REAL(wp), DIMENSION(0:6,1:7), PARAMETER :: pavement_pars = RESHAPE( (/ & 652 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 17.0_wp, 0.97_wp, & ! 1 653 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 18.0_wp, 0.94_wp, & ! 2 654 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 19.0_wp, 0.98_wp, & ! 3 655 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 20.0_wp, 0.93_wp, & ! 4 656 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 21.0_wp, 0.97_wp, & ! 5 657 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 22.0_wp, 0.97_wp, & ! 6 658 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 23.0_wp, 0.97_wp & ! 7 659 /), (/ 7, 7 /) ) 666 REAL(wp), DIMENSION(0:6,1:15), PARAMETER :: pavement_pars = RESHAPE( (/ & 667 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 18.0_wp, 0.97_wp, & ! 1 668 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 19.0_wp, 0.94_wp, & ! 2 669 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 20.0_wp, 0.98_wp, & ! 3 670 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 21.0_wp, 0.93_wp, & ! 4 671 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 22.0_wp, 0.97_wp, & ! 5 672 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 23.0_wp, 0.97_wp, & ! 6 673 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 24.0_wp, 0.97_wp, & ! 7 674 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 25.0_wp, 0.94_wp, & ! 8 675 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 26.0_wp, 0.98_wp, & ! 9 676 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 27.0_wp, 0.93_wp, & ! 10 677 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 28.0_wp, 0.97_wp, & ! 11 678 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 29.0_wp, 0.97_wp, & ! 12 679 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 30.0_wp, 0.97_wp, & ! 13 680 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 31.0_wp, 0.94_wp, & ! 14 681 0.050_wp, 1.0E-4_wp, 1.0E-5_wp, 1.00_wp, 1.94E6_wp, 32.0_wp, 0.98_wp & ! 15 682 /), (/ 7, 15 /) ) 660 683 661 684 ! … … 2426 2449 surf_lsm_h%lambda_h_def(m) = pavement_heat_conduct 2427 2450 surf_lsm_h%rho_c_def(m) = pavement_heat_capacity 2428 2451 2429 2452 ELSEIF ( surf_lsm_h%vegetation_surface(m) ) THEN 2430 2453 … … 2692 2715 j = surf_lsm_h%j(m) 2693 2716 DO k = nzb_soil, nzt_soil 2694 surf_lsm_h%root_fr(k,m) = root_fraction(k) 2695 ENDDO 2717 IF ( surf_lsm_h%pavement_surface(m) .AND. zs(k) <= pavement_depth ) THEN 2718 surf_lsm_h%root_fr(k,m) = 0.0_wp 2719 ELSE 2720 surf_lsm_h%root_fr(k,m) = root_fraction(k) 2721 ENDIF 2722 ENDDO 2723 2724 ! 2725 !-- Normalize so that the sum = 1. Only relevant when the root distribution was 2726 !-- set to zero due to pavement at some layers. 2727 IF ( SUM( surf_lsm_h%root_fr(:,m) ) > 0.0_wp ) THEN 2728 DO k = nzb_soil, nzt_soil 2729 surf_lsm_h%root_fr(k,m) = surf_lsm_h%root_fr(k,m) / SUM( surf_lsm_h%root_fr(:,m) ) 2730 ENDDO 2731 ENDIF 2732 2696 2733 ENDDO 2697 2734 DO l = 0, 3 … … 2701 2738 2702 2739 DO k = nzb_soil, nzt_soil 2703 surf_lsm_v(l)%root_fr(k,m) = root_fraction(k) 2740 IF ( surf_lsm_v(l)%pavement_surface(m) .AND. zs(k) <= pavement_depth ) THEN 2741 surf_lsm_v(l)%root_fr(k,m) = 0.0_wp 2742 ELSE 2743 surf_lsm_v(l)%root_fr(k,m) = root_fraction(k) 2744 ENDIF 2704 2745 ENDDO 2746 ! 2747 !-- Normalize so that the sum = 1. Only relevant when the root distribution was 2748 !-- sset to zero due to pavement at some layers. 2749 IF ( SUM( surf_lsm_v(l)%root_fr(:,m) ) > 0.0_wp ) THEN 2750 DO k = nzb_soil, nzt_soil 2751 surf_lsm_v(l)%root_fr(k,m) = surf_lsm_v(l)%root_fr(k,m) / SUM( surf_lsm_v(l)%root_fr(:,m) ) 2752 ENDDO 2753 ENDIF 2754 2705 2755 ENDDO 2706 2756 ENDDO … … 3161 3211 3162 3212 ! 3163 !-- Calculate soil diffusivity at the center of the soil layers 3164 lambda_temp(k) = (- b_ch * surf%gamma_w_sat(k,m) * psi_sat & 3165 / surf%m_sat(k,m) ) * ( & 3166 MAX( surf_m_soil%var_2d(k,m), & 3167 surf%m_wilt(k,m) ) / surf%m_sat(k,m) )**( & 3168 b_ch + 2.0_wp ) 3169 3170 ! 3171 !-- Parametrization of Van Genuchten 3172 !-- Calculate the hydraulic conductivity after Van Genuchten (1980) 3173 h_vg = ( ( ( surf%m_res(k,m) - surf%m_sat(k,m) ) / & 3174 ( surf%m_res(k,m) - & 3175 MAX( surf_m_soil%var_2d(k,m), surf%m_wilt(k,m) ) & 3176 ) & 3177 )**( surf%n_vg(k,m) / ( surf%n_vg(k,m) - 1.0_wp ) & 3178 ) - 1.0_wp & 3179 )**( 1.0_wp / surf%n_vg(k,m) ) / surf%alpha_vg(k,m) 3180 3181 gamma_temp(k) = surf%gamma_w_sat(k,m) * ( ( ( 1.0_wp + & 3182 ( surf%alpha_vg(k,m) * h_vg )**surf%n_vg(k,m) & 3183 )**( 1.0_wp - 1.0_wp / surf%n_vg(k,m) ) & 3184 - ( surf%alpha_vg(k,m) * h_vg )**( surf%n_vg(k,m) & 3185 - 1.0_wp) )**2 ) / ( ( 1.0_wp + ( surf%alpha_vg(k,m) & 3186 * h_vg )**surf%n_vg(k,m) )**( ( 1.0_wp - 1.0_wp & 3187 / surf%n_vg(k,m) ) * ( surf%l_vg(k,m) + 2.0_wp) ) ) 3213 !-- In order to prevent water tranport through paved surfaces, 3214 !-- conductivity and diffusivity are set to zero 3215 IF ( surf%pavement_surface(m) .AND. zs(k) <= pavement_depth ) THEN 3216 3217 lambda_temp(k) = 0.0_wp 3218 gamma_temp(k) = 0.0_wp 3219 3220 ELSE 3221 3222 ! 3223 !-- Calculate soil diffusivity at the center of the soil layers 3224 lambda_temp(k) = (- b_ch * surf%gamma_w_sat(k,m) * psi_sat & 3225 / surf%m_sat(k,m) ) * ( & 3226 MAX( surf_m_soil%var_2d(k,m), & 3227 surf%m_wilt(k,m) ) / surf%m_sat(k,m) )**( & 3228 b_ch + 2.0_wp ) 3229 3230 ! 3231 !-- Parametrization of Van Genuchten 3232 !-- Calculate the hydraulic conductivity after Van Genuchten (1980) 3233 h_vg = ( ( ( surf%m_res(k,m) - surf%m_sat(k,m) ) / & 3234 ( surf%m_res(k,m) - & 3235 MAX(surf_m_soil%var_2d(k,m), surf%m_wilt(k,m)) & 3236 ) & 3237 )**( surf%n_vg(k,m) / ( surf%n_vg(k,m) - 1.0_wp ) & 3238 ) - 1.0_wp & 3239 )**( 1.0_wp / surf%n_vg(k,m) ) / surf%alpha_vg(k,m) 3240 3241 gamma_temp(k) = surf%gamma_w_sat(k,m) * ( ( ( 1.0_wp + & 3242 ( surf%alpha_vg(k,m) * h_vg )**surf%n_vg(k,m) & 3243 )**( 1.0_wp - 1.0_wp / surf%n_vg(k,m) ) & 3244 - ( surf%alpha_vg(k,m) * h_vg )**( surf%n_vg(k,m) & 3245 - 1.0_wp) )**2 ) / ( ( 1.0_wp + ( surf%alpha_vg(k,m) & 3246 * h_vg )**surf%n_vg(k,m) )**(( 1.0_wp - 1.0_wp & 3247 / surf%n_vg(k,m) ) * ( surf%l_vg(k,m) + 2.0_wp) )) 3248 ENDIF 3188 3249 3189 3250 ENDDO … … 3201 3262 ! 3202 3263 !-- Prognostic equation for soil moisture content. Only performed, 3203 !-- when humidity is enabled in the atmosphere and the surface type 3204 !-- is not pavement (implies dry soil below). 3205 IF ( humidity .AND. .NOT. surf%pavement_surface(m) ) THEN 3264 !-- when humidity is enabled in the atmosphere. 3265 IF ( humidity ) THEN 3206 3266 ! 3207 3267 !-- Calculate soil diffusivity (lambda_w) at the _layer level … … 3209 3269 !-- ECMWF-IFS Eq. 8.81 3210 3270 DO k = nzb_soil, nzt_soil-1 3211 3212 surf%lambda_w(k,m) = ( lambda_temp(k+1) + & 3213 lambda_temp(k) ) * 0.5_wp 3214 surf%gamma_w(k,m) = ( gamma_temp(k+1) + & 3215 gamma_temp(k) ) * 0.5_wp 3216 3271 IF ( surf%pavement_surface(m) .AND. zs(k) < pavement_depth& 3272 .AND. zs(k+1) > pavement_depth ) & 3273 THEN 3274 surf%lambda_w(k,m) = ( pavement_depth - zs(k) ) & 3275 * ddz_soil_center(k+1) * lambda_temp(k) & 3276 + ( zs(k+1) - pavement_depth ) & 3277 * ddz_soil_center(k+1) * lambda_temp(k+1) 3278 3279 surf%gamma_w(k,m) = ( pavement_depth - zs(k) ) & 3280 * ddz_soil_center(k+1) * gamma_temp(k) & 3281 + ( zs(k+1) - pavement_depth ) & 3282 * ddz_soil_center(k+1) * gamma_temp(k+1) 3283 ELSE 3284 surf%lambda_w(k,m) = ( lambda_temp(k+1) + & 3285 lambda_temp(k) ) * 0.5_wp 3286 surf%gamma_w(k,m) = ( gamma_temp(k+1) + & 3287 gamma_temp(k) ) * 0.5_wp 3288 ENDIF 3217 3289 ENDDO 3218 3290 … … 3221 3293 !-- In case of a closed bottom (= water content is conserved), 3222 3294 !-- set hydraulic conductivity to zero to that no water will be 3223 !-- lost in the bottom layer. 3224 IF ( conserve_water_content ) THEN 3295 !-- lost in the bottom layer. As gamma_w is always a positive value, 3296 !-- it cannot be set to zero in case of purely dry soil since this 3297 !-- would cause accumulation of (non-existing) water in the lowest 3298 !-- soil layer 3299 IF ( conserve_water_content .AND. surf_m_soil%var_2d(nzt_soil,m) /= 0.0_wp ) THEN 3225 3300 surf%gamma_w(nzt_soil,m) = 0.0_wp 3226 3301 ELSE … … 3294 3369 + dt_3d * ( tsc(2) * tend(:) & 3295 3370 + tsc(3) * surf_tm_soil_m%var_2d(:,m) ) 3296 3297 3371 ! 3298 3372 !-- Account for dry soils (find a better solution here!)
Note: See TracChangeset
for help on using the changeset viewer.