Ignore:
Timestamp:
Sep 27, 2017 10:36:13 AM (7 years ago)
Author:
maronga
Message:

lsm now allows for moist soil and roots below paved surfaces

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/land_surface_model_mod.f90

    r2476 r2504  
    2525! -----------------
    2626! $Id$
     27! Support roots and water under pavement. Added several pavement types.
     28!
     29! 2476 2017-09-18 07:54:32Z maronga
    2730! Bugfix for last commit
    2831!
     
    558561!
    559562!-- 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
    569580                                                                 /)                                                                 
    570                                                                  
     581   
     582
     583
     584
     585                                                             
    571586!
    572587!-- Water classes
     
    587602!-- r_canopy_min,     lai,   c_veg,     g_d         z0,         z0h, lambda_s_s, lambda_s_u, f_sw_in,  c_surface, albedo_type, emissivity
    588603    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, 18.0_wp, 0.94_wp, & !  1
     604          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
    590605        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
    591606        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
     
    649664!-- TO BE FILLED
    650665!-- 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 /) )                             
    660683
    661684!
     
    24262449                surf_lsm_h%lambda_h_def(m)    = pavement_heat_conduct                                     
    24272450                surf_lsm_h%rho_c_def(m)       = pavement_heat_capacity                                               
    2428                
     2451
    24292452             ELSEIF ( surf_lsm_h%vegetation_surface(m) )  THEN
    24302453
     
    26922715             j   = surf_lsm_h%j(m)
    26932716             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
    26962733          ENDDO
    26972734          DO  l = 0, 3
     
    27012738
    27022739                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
    27042745                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
    27052755             ENDDO
    27062756           ENDDO
     
    31613211
    31623212!
    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
    31883249
    31893250             ENDDO
     
    32013262!
    32023263!--          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
    32063266!
    32073267!--             Calculate soil diffusivity (lambda_w) at the _layer level
     
    32093269!--             ECMWF-IFS Eq. 8.81
    32103270                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
    32173289                ENDDO
    32183290
     
    32213293!--             In case of a closed bottom (= water content is conserved),
    32223294!--             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
    32253300                   surf%gamma_w(nzt_soil,m) = 0.0_wp
    32263301                ELSE
     
    32943369                                       + dt_3d * ( tsc(2) * tend(:)            &
    32953370                                       + tsc(3) * surf_tm_soil_m%var_2d(:,m) )   
    3296    
    32973371!
    32983372!--             Account for dry soils (find a better solution here!)
Note: See TracChangeset for help on using the changeset viewer.