Index: palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2339)
+++ palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2340)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Revised root_distribution tabel and implemented a pseudo-generic root fraction
+! calculation
+! 
+! 2333 2017-08-04 09:08:26Z maronga
 ! minor bugfixes
 ! 
@@ -306,4 +310,12 @@
                                                0.07_wp, 0.15_wp, 0.21_wp,      &
                                                0.72_wp, 1.89_wp/)
+
+    REAL(wp), DIMENSION(0:3), PARAMETER  :: dz_soil_ref =                      & ! reference four layer soil configuration used for estimating the root fractions
+                                            (/ 0.07_wp, 0.21_wp, 0.72_wp,      &
+                                               1.89_wp /)
+
+    REAL(wp), DIMENSION(0:3), PARAMETER  :: zs_ref =                           & ! reference four layer soil configuration used for estimating the root fractions
+                                            (/ 0.07_wp, 0.28_wp, 1.0_wp,       &
+                                               2.89_wp /)
 
 
@@ -589,44 +601,25 @@
 !
 !-- Root distribution for default soil layer configuration (sum = 1)
-!--                                level 1 - level 8
-    REAL(wp), DIMENSION(0:7,1:18), PARAMETER :: root_distribution = RESHAPE( (/&
-                                   1.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     &
-                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     & ! 1
-                                   0.035_wp, 0.069_wp, 0.069_wp, 0.108_wp,     &
-                                   0.195_wp, 0.214_wp, 0.284_wp, 0.026_wp,     & ! 2
-                                   0.050_wp, 0.100_wp, 0.100_wp, 0.136_wp,     &
-                                   0.181_wp, 0.192_wp, 0.215_wp, 0.026_wp,     & ! 3
-                                   0.038_wp, 0.075_wp, 0.075_wp, 0.111_wp,     &
-                                   0.185_wp, 0.203_wp, 0.273_wp, 0.040_wp,     & ! 4
-                                   0.038_wp, 0.075_wp, 0.075_wp, 0.110_wp,     &
-                                   0.180_wp, 0.199_wp, 0.277_wp, 0.046_wp,     & ! 5
-                                   0.035_wp, 0.069_wp, 0.069_wp, 0.105_wp,     &
-                                   0.180_wp, 0.201_wp, 0.295_wp, 0.046_wp,     & ! 6
-                                   0.035_wp, 0.072_wp, 0.072_wp, 0.105_wp,     &
-                                   0.161_wp, 0.180_wp, 0.282_wp, 0.093_wp,     & ! 7
-                                   0.040_wp, 0.077_wp, 0.077_wp, 0.112_wp,     &
-                                   0.176_wp, 0.192_wp, 0.266_wp, 0.060_wp,     & ! 8
-                                   0.142_wp, 0.286_wp, 0.286_wp, 0.286_wp,     &
-                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     & ! 9
-                                   0.068_wp, 0.134_wp, 0.134_wp, 0.177_wp,     &
-                                   0.214_wp, 0.203_wp, 0.070_wp, 0.000_wp,     & ! 10
-                                   0.035_wp, 0.068_wp, 0.068_wp, 0.108_wp,     &
-                                   0.195_wp, 0.215_wp, 0.285_wp, 0.026_wp,     & ! 11
-                                   0.025_wp, 0.048_wp, 0.048_wp, 0.078_wp,     &
-                                   0.147_wp, 0.175_wp, 0.353_wp, 0.126_wp,     & ! 12
-                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     &
-                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     & ! 13
-                                   0.036_wp, 0.072_wp, 0.072_wp, 0.103_wp,     &
-                                   0.163_wp, 0.180_wp, 0.273_wp, 0.074_wp,     & ! 14
-                                   0.032_wp, 0.066_wp, 0.066_wp, 0.100_wp,     &
-                                   0.172_wp, 0.192_wp, 0.299_wp, 0.073_wp,     & ! 15
-                                   0.032_wp, 0.066_wp, 0.066_wp, 0.100_wp,     &
-                                   0.172_wp, 0.192_wp, 0.299_wp, 0.073_wp,     & ! 16
-                                   0.028_wp, 0.055_wp, 0.055_wp, 0.087_wp,     &
-                                   0.166_wp, 0.195_wp, 0.348_wp, 0.066_wp,     & ! 17
-                                   0.028_wp, 0.055_wp, 0.055_wp, 0.087_wp,     &
-                                   0.166_wp, 0.195_wp, 0.348_wp, 0.066_wp      & ! 18
-                                 /), (/ 8, 18 /) )
-
+!--                                level 1 - level 4 according to zs_ref
+    REAL(wp), DIMENSION(0:3,1:18), PARAMETER :: root_distribution = RESHAPE( (/ &
+                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  1
+                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & !  2
+                                 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp,            & !  3
+                                 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp,            & !  4
+                                 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp,            & !  5
+                                 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp,            & !  6
+                                 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp,            & !  7
+                                 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp,            & !  8
+                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  9
+                                 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp,            & ! 10
+                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & ! 11
+                                 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp,            & ! 12
+                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 13
+                                 0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp,            & ! 14
+                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 15 
+                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 16 
+                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 17
+                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp             & ! 18
+                                 /), (/ 4, 18 /) )
 
 !
@@ -1238,11 +1231,4 @@
           IF ( dz_soil(k) /= 9999999.9_wp )  THEN
              nzt_soil = nzt_soil + 1
-             IF ( root_fraction(k) == 9999999.9_wp )  THEN
-                message_string = 'manual setting of dz_soil '//                &
-                                 'requires adequate setting of root_fraction'//&
-                                 '/= 9999999.9 ' //                            &
-                                 'and SUM(root_fraction) = 1'
-                CALL message( 'check_parameters', 'PA0452', 1, 2, 0, 6, 0 )
-             ENDIF
           ENDIF
        ENDDO    
@@ -2035,8 +2021,16 @@
        INTEGER(iwp) ::  j_off !< index offset of surface element, seen from atmospheric grid point 
        INTEGER(iwp) ::  k     !< running index
+       INTEGER(iwp) ::  kn    !< running index
+       INTEGER(iwp) ::  ko    !< running index
+       INTEGER(iwp) ::  kroot !< running index
+       INTEGER(iwp) ::  kzs   !< running index
        INTEGER(iwp) ::  l     !< running index surface facing
        INTEGER(iwp) ::  m     !< running index
+       INTEGER(iwp) ::  n_soil_layers_total     !< temperature variable, stores the total number of soil layers + 4 
+
 
        REAL(wp) :: pt1   !< potential temperature at first grid level
+
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  bound, bound_root_fr  !< temporary arrays for storing index bounds
 
 !
@@ -2613,15 +2607,67 @@
        ENDDO
 
-
-
+!
+!--    When no root distribution is given by the user, use look-up table to prescribe
+!--    the root fraction in the individual soil layers
        IF ( ALL( root_fraction == 9999999.9_wp ) )  THEN
+
+!
+!--       First, calculate the index bounds for integration
+          n_soil_layers_total = nzt_soil - nzb_soil + 6
+          ALLOCATE ( bound(0:n_soil_layers_total) )
+          ALLOCATE ( bound_root_fr(0:n_soil_layers_total) )
+
+          kn = 0
+          ko = 0
+          bound(0) = 0.0_wp
+          DO k = 1, n_soil_layers_total-1
+             IF ( zs_layer(kn) <= zs_ref(ko) )  THEN
+                bound(k) = zs_layer(kn)
+                bound_root_fr(k) = ko
+                kn = kn + 1
+                IF ( kn > nzt_soil+1 )  THEN
+                   kn = nzt_soil
+                ENDIF
+             ELSE
+                bound(k) = zs_ref(ko)
+                bound_root_fr(k) = ko
+                ko = ko + 1
+                IF ( ko > 3 )  THEN
+                   ko = 3
+                ENDIF
+             ENDIF
+
+          ENDDO
+
+!
+!--       Integrate over all soil layers based on the four-layer root fraction
+          kzs = 1
+          root_fraction = 0.0_wp
+          DO k = 0, n_soil_layers_total-2
+             kroot = bound_root_fr(k+1)
+             root_fraction(kzs-1) = root_fraction(kzs-1)                       &
+                                + root_distribution(kroot,vegetation_type)     &
+                                / dz_soil_ref(kroot) * ( bound(k+1) - bound(k) )
+
+             IF ( bound(k+1) == zs_layer(kzs-1) )  THEN
+                kzs = kzs+1
+             ENDIF
+          ENDDO
+
+
+!
+!--       Normalize so that the sum of all fractions equals one
+          root_fraction = root_fraction / SUM(root_fraction)
+
+          DEALLOCATE ( bound )
+          DEALLOCATE ( bound_root_fr )
+
+!
+!--       Map calculated root fractions
           DO  m = 1, surf_lsm_h%ns
              i   = surf_lsm_h%i(m)            
              j   = surf_lsm_h%j(m)
-
-
              DO  k = nzb_soil, nzt_soil 
-                surf_lsm_h%root_fr(k,m) = root_distribution(k,vegetation_type)
-                root_fraction(k)        = root_distribution(k,vegetation_type)
+                surf_lsm_h%root_fr(k,m) = root_fraction(k)
              ENDDO
           ENDDO
@@ -2632,11 +2678,9 @@
 
                 DO  k = nzb_soil, nzt_soil
-                   surf_lsm_v(l)%root_fr(k,m) = root_distribution(k,vegetation_type)
-                   root_fraction(k)           = root_distribution(k,vegetation_type)
+                   surf_lsm_v(l)%root_fr(k,m) = root_fraction(k)
                 ENDDO
              ENDDO
            ENDDO
        ENDIF
-
 !
 !--    Map vegetation parameters to the respective 2D arrays
