Ignore:
Timestamp:
Aug 3, 2017 12:34:22 PM (7 years ago)
Author:
maronga
Message:

some improvements in land surface model

File:
1 edited

Legend:

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

    r2307 r2328  
    2525! -----------------
    2626! $Id$
     27! Revised skin layer concept.
     28! Bugfix for runs with pavement surface and humidity
     29! Revised some standard values in vegetation_pars
     30! Added emissivity and default albedo_type as variable to tables
     31! Changed default surface type to vegetation
     32! Revised input of soil layer configuration
     33!
     34! 2307 2017-07-07 11:32:10Z suehring
    2735! Bugfix, variable names corrected
    2836!
     
    257265
    258266    USE radiation_model_mod,                                                   &
    259         ONLY:  force_radiation_call, rad_net, rad_sw_in, rad_lw_out,           &
    260                rad_lw_out_change_0, radiation_scheme,                          &
     267        ONLY:  albedo, albedo_type, emissivity, force_radiation_call, rad_net,         &
     268               rad_sw_in, rad_lw_out, rad_lw_out_change_0, radiation_scheme,   &
    261269               unscheduled_radiation_calls
    262270       
     
    288296
    289297
    290     REAL(wp), DIMENSION(0:7), PARAMETER  :: zs_default =                       & ! default soil layer configuration
    291                                             (/ 0.005_wp, 0.02_wp, 0.04_wp,     &
    292                                                0.07_wp, 0.15_wp, 0.28_wp,      &
    293                                                1.00_wp, 2.89_wp /)
     298    REAL(wp), DIMENSION(0:7), PARAMETER  :: dz_soil_default =                  & ! default soil layer configuration
     299                                            (/ 0.01_wp, 0.02_wp, 0.04_wp,      &
     300                                               0.07_wp, 0.15_wp, 0.21_wp,      &
     301                                               0.72_wp, 1.89_wp/)
    294302
    295303
    296304!
    297305!-- LSM variables
    298     CHARACTER(10) :: surface_type = 'netcdf'  !< general classification. Allowed are:
    299                                              !< 'vegetation', 'pavement', ('building'),
    300                                              !< 'water', and 'netcdf'
     306    CHARACTER(10) :: surface_type = 'vegetation'  !< general classification. Allowed are:
     307                                                  !< 'vegetation', 'pavement', ('building'),
     308                                                  !< 'water', and 'netcdf'
    301309
    302310
     
    362370               
    363371               
    364     REAL(wp), DIMENSION(:), ALLOCATABLE  :: ddz_soil,         & !< 1/dz_soil
    365                                             ddz_soil_layer,   & !< 1/dz_soil_layer
    366                                             dz_soil,          & !< soil grid spacing (center-center)
    367                                             dz_soil_layer,    & !< soil grid spacing (edge-edge)
    368                                             root_extr           !< root extraction
     372    REAL(wp), DIMENSION(:), ALLOCATABLE  :: ddz_soil_center, & !< 1/dz_soil_center
     373                                            ddz_soil,        & !< 1/dz_soil
     374                                            dz_soil_center,  & !< soil grid spacing (center-center)
     375                                            zs,              & !< depth of the temperature/moisute levels
     376                                            root_extr          !< root extraction
    369377
    370378
     
    373381                                   soil_moisture = 0.0_wp,                     & !< NAMELIST soil moisture content (m3/m3)
    374382                                   soil_temperature = 300.0_wp,                & !< NAMELIST soil temperature (K) +1
    375                                    zs = 9999999.9_wp,                          & !< (NAMELIST) soil layer depths (center)
     383                                   dz_soil  = 9999999.9_wp,                    & !< (NAMELIST) soil layer depths (spacing)
    376384                                   zs_layer = 9999999.9_wp                       !< soil layer depths (edge)
    377385                                 
     
    536544!
    537545!-- Water classes
    538     CHARACTER(12), DIMENSION(0:4), PARAMETER :: water_type_name = (/ &
     546    CHARACTER(12), DIMENSION(0:5), PARAMETER :: water_type_name = (/ &
    539547                                   'user defined',                   & ! 0
    540548                                   'lake        ',                   & ! 1
    541549                                   'river       ',                   & ! 2
    542550                                   'ocean       ',                   & ! 3
    543                                    'pond        '                    & ! 4
     551                                   'pond        ',                   & ! 4
     552                                   'fountain    '                    & ! 5
    544553                                                                  /)                                                                 
    545554                                                                 
     
    549558!
    550559!-- Land surface parameters
    551 !-- r_canopy_min,     lai,   c_veg,     g_d         z0,         z0h, lambda_s_s, lambda_s_u, f_sw_in,  c_surface
    552     REAL(wp), DIMENSION(0:9,1:18), PARAMETER :: vegetation_pars = RESHAPE( (/ &
    553           0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, 0.5E-2_wp,   0.5E-4_wp,     0.0_wp,     0.0_wp, 0.00_wp, 0.00_wp, & !  1
    554         180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,   0.25_wp,  0.25E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & !  2
    555         110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp,   0.20_wp,  0.20E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & !  3
    556         500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & !  4
    557         500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & !  5
    558         175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & !  6
    559         240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & !  7
    560         100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp,   0.47_wp,  0.47E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & !  8
    561         250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp,  0.013_wp, 0.013E-2_wp,   196.0_wp,   196.0_wp, 0.00_wp, 0.00_wp, & !  9
    562          80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp,  0.034_wp, 0.034E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 10
    563         180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,    0.5_wp,  0.50E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 11
    564         150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp,   0.17_wp,  0.17E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 12
    565           0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, 1.3E-3_wp,   1.3E-4_wp,   812.0_wp,   812.0_wp, 0.00_wp, 0.00_wp, & ! 13
    566         240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp,   0.83_wp,  0.83E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 14
    567         225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp,   0.10_wp,  0.10E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 15
    568         225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp,   0.25_wp,  0.25E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 16
    569         250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & ! 17
    570         175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp,   1.10_wp,     1.10_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp  & ! 18
    571                                                                /), (/ 10, 18 /) )
     560!-- r_canopy_min,     lai,   c_veg,     g_d         z0,         z0h, lambda_s_s, lambda_s_u, f_sw_in,  c_surface, albedo_type, emissivity
     561    REAL(wp), DIMENSION(0:11,1:18), PARAMETER :: vegetation_pars = RESHAPE( (/ &
     562          0.0_wp, 0.00_wp, 1.00_wp, 0.00_wp,  0.005_wp,   0.5E-4_wp,     0.0_wp,    0.0_wp, 0.00_wp, 0.00_wp,  0.0_wp, 0.94_wp, & !  1
     563        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
     564        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
     565        500.0_wp, 5.00_wp, 1.00_wp, 0.03_wp,   2.00_wp,     2.00_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  5.0_wp, 0.97_wp, & !  4
     566        500.0_wp, 5.00_wp, 1.00_wp, 0.03_wp,   2.00_wp,     2.00_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  6.0_wp, 0.97_wp, & !  5
     567        175.0_wp, 5.00_wp, 1.00_wp, 0.03_wp,   2.00_wp,     2.00_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  8.0_wp, 0.97_wp, & !  6
     568        240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp,   2.00_wp,     2.00_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  9.0_wp, 0.97_wp, & !  7
     569        100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp,   0.47_wp,  0.47E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  8.0_wp, 0.97_wp, & !  8
     570        250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp,  0.013_wp, 0.013E-2_wp,    15.0_wp,   15.0_wp, 0.00_wp, 0.00_wp,  3.0_wp, 0.94_wp, & !  9
     571         80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp,  0.034_wp, 0.034E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp, 11.0_wp, 0.97_wp, & ! 10
     572        180.0_wp, 3.00_wp, 1.00_wp, 0.00_wp,    0.5_wp,  0.50E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp, 13.0_wp, 0.97_wp, & ! 11
     573        150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp,   0.17_wp,  0.17E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  2.0_wp, 0.97_wp, & ! 12
     574          0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, 1.3E-3_wp,   1.3E-4_wp,    58.0_wp,   58.0_wp, 0.00_wp, 0.00_wp, 11.0_wp, 0.97_wp, & ! 13
     575        240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp,   0.83_wp,  0.83E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  4.0_wp, 0.97_wp, & ! 14
     576        225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp,   0.10_wp,  0.10E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  4.0_wp, 0.97_wp, & ! 15
     577        225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp,   0.25_wp,  0.25E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  4.0_wp, 0.97_wp, & ! 16
     578        250.0_wp, 5.00_wp, 1.00_wp, 0.03_wp,   2.00_wp,     2.00_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  7.0_wp, 0.97_wp, & ! 17
     579        175.0_wp, 2.50_wp, 1.00_wp, 0.03_wp,   1.10_wp,     1.10_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  8.0_wp, 0.97_wp  & ! 18
     580                                                               /), (/ 12, 18 /) )
    572581
    573582                                   
     
    631640!
    632641!-- TO BE FILLED
    633 !-- Pavement parameters  depth,        z0,       z0h, lambda_h,      rho_c 
    634     REAL(wp), DIMENSION(0:4,1:7), PARAMETER :: pavement_pars = RESHAPE( (/ &
    635                       0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 1
    636                       0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 2
    637                       0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 3                                 
    638                       0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 4
    639                       0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 5
    640                       0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 6
    641                       0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp  & ! 7
    642                                                                  /), (/ 5, 7 /) )                             
     642!-- Pavement parameters  depth,        z0,       z0h, lambda_h,      rho_c, albedo_type, emissivity 
     643    REAL(wp), DIMENSION(0:5,1:7), PARAMETER :: pavement_pars = RESHAPE( (/ &
     644                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 17.0_wp, 0.97_wp, & ! 1
     645                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 18.0_wp, 0.94_wp,  & ! 2
     646                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 19.0_wp, 0.98_wp,  & ! 3                                 
     647                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 20.0_wp, 0.93_wp,  & ! 4
     648                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 21.0_wp, 0.97_wp,  & ! 5
     649                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 22.0_wp, 0.97_wp,  & ! 6
     650                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 23.0_wp, 0.97_wp  & ! 7
     651                                                                 /), (/ 6, 7 /) )                             
    643652
    644653!
    645654!-- TO BE FILLED
    646 !-- Water parameters                    temperature,     z0,      z0h
    647     REAL(wp), DIMENSION(0:2,1:4), PARAMETER :: water_pars = RESHAPE( (/ &
    648                                            283.0_wp, 0.01_wp, 0.001_wp, & ! 1
    649                                            283.0_wp, 0.01_wp, 0.001_wp, & ! 2
    650                                            283.0_wp, 0.01_wp, 0.001_wp, & ! 3
    651                                            283.0_wp, 10.01_wp, 0.001_wp & ! 4
    652                                                                 /), (/ 3, 4 /) )                                                                     
     655!-- Water parameters                    temperature,     z0,      z0h, albedo_type, emissivity,
     656    REAL(wp), DIMENSION(0:4,1:5), PARAMETER :: water_pars = RESHAPE( (/ &
     657                                           283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 1
     658                                           283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 2
     659                                           283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 3
     660                                           283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 4
     661                                           283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp  & ! 5
     662                                                                /), (/ 5, 5 /) )                                                                     
    653663                                                                                                                                     
    654664    SAVE
     
    10991109    ENDIF
    11001110   
    1101    
    1102      
    11031111    IF ( TRIM( surface_type ) == 'pavement' )  THEN
    11041112
     
    11441152    ENDIF
    11451153   
     1154!
     1155!-- Temporary message as long as NetCDF input is not available
     1156    IF ( TRIM( surface_type ) == 'netcdf' )  THEN
     1157             message_string = 'surface_type = netcdf '//                       &
     1158                              'is not supported at the moment'
     1159             CALL message( 'check_parameters', 'PA0465', 1, 2, 0, 6, 0 )
     1160    ENDIF
     1161
    11461162    IF ( soil_type == 0 )  THEN
    11471163
     
    12091225    nzb_soil = 0
    12101226    nzt_soil = -1
    1211     IF ( ALL( zs == 9999999.9_wp ) )  THEN
     1227    IF ( ALL( dz_soil == 9999999.9_wp ) )  THEN
    12121228       nzt_soil = 7
    1213        zs(nzb_soil:nzt_soil) = zs_default
     1229       dz_soil(nzb_soil:nzt_soil) = dz_soil_default
    12141230    ELSE
    12151231       DO k = 0, 19
    1216           IF ( zs(k) /= 9999999.9_wp )  THEN
     1232          IF ( dz_soil(k) /= 9999999.9_wp )  THEN
    12171233             nzt_soil = nzt_soil + 1
    12181234             IF ( root_fraction(k) == 9999999.9_wp )  THEN
    1219                 message_string = 'manual setting of zs '//                     &
     1235                message_string = 'manual setting of dz_soil '//                &
    12201236                                 'requires adequate setting of root_fraction'//&
    12211237                                 '/= 9999999.9 ' //                            &
     
    12921308                f_qsws_liq,  & !< factor for qsws_liq
    12931309                f_shf,       & !< factor for shf
     1310                lambda_soil, & !< Thermal conductivity of the uppermost soil layer (W/m2/K)
    12941311                lambda_surface, & !< Current value of lambda_surface (W/m2/K)
    12951312                m_liq_max,   & !< maxmimum value of the liq. water reservoir
     
    13701387!
    13711388!--    Define heat conductivity between surface and soil depending on surface
    1372 !--    type. For vegetation, a skin layer parameterization is used. For bare
    1373 !--    soil and pavements, no skin layer is applied. In these cases, the
    1374 !--    temperature is assumed to be constant between the surface and the first
    1375 !--    soil layer. The heat conductivity is then derived from the soil/
    1376 !--    pavement properties. For water surfaces, the conductivity is already set
    1377 !--    to 1E10.
     1389!--    type. For vegetation, a skin layer parameterization is used. The new
     1390!--    parameterization uses a combination of two conductivities: a constant
     1391!--    conductivity for the skin layer, and a conductivity according to the
     1392!--    uppermost soil layer. For bare soil and pavements, no skin layer is
     1393!--    applied. In these cases, the temperature is assumed to be constant
     1394!--    between the surface and the first soil layer. The heat conductivity is
     1395!--    then derived from the soil/pavement properties.
     1396!--    For water surfaces, the conductivity is already set to 1E10.
    13781397!--    Moreover, the heat capacity is set. For bare soil the heat capacity is
    1379 !--    the capacity of the uppermost soil layer
    1380        IF ( surf%lambda_surface_s(m) == 0.0_wp )  THEN
    1381        
     1398!--    the capacity of the uppermost soil layer, for pavement it is that of
     1399!--    the material involved.
     1400
     1401!
     1402!--    for vegetation type surfaces, the thermal conductivity of the soil is
     1403!--    needed
     1404       IF ( surf%vegetation_surface(m) )  THEN
     1405
    13821406          lambda_h_sat = lambda_h_sm**(1.0_wp - surf%m_sat(nzb_soil,m)) *      &
    13831407                         lambda_h_water ** surf_m_soil%var_2d(nzb_soil,m)
     
    13861410                                                     surf%m_sat(nzb_soil,m) ) )                   
    13871411                         
    1388           lambda_surface = (ke * (lambda_h_sat - lambda_h_dry) + lambda_h_dry )&
    1389                            * ddz_soil_layer(nzb_soil) * 2.0_wp
    1390          
    1391           c_surface_tmp = (rho_c_soil * ( 1.0_wp - surf%m_sat(nzb_soil,m) )    &
    1392                            + rho_c_water * surf_m_soil%var_2d(nzb_soil,m) )    &
    1393                            * dz_soil_layer(nzb_soil) * 0.5_wp               
    1394 
    1395        ELSEIF ( surf_t_surface%var_1d(m) >= surf_t_soil%var_2d(nzb_soil,m)) &
    1396        THEN
    1397 
     1412          lambda_soil = (ke * (lambda_h_sat - lambda_h_dry) + lambda_h_dry )   &
     1413                           * ddz_soil(nzb_soil) * 2.0_wp
     1414
     1415!
     1416!--       When bare soil is set without a thermal conductivity (no skin layer),
     1417!--       a heat capacity is that of the soil layer, otherwise it is a
     1418!--       combination of the conductivities from the skin and the soil layer
     1419          IF ( surf%lambda_surface_s(m) == 0.0_wp )  THEN
     1420            surf%c_surface(m) = (rho_c_soil * (1.0_wp - surf%m_sat(nzb_soil,m))&
     1421                              + rho_c_water * surf_m_soil%var_2d(nzb_soil,m) ) &
     1422                              * dz_soil(nzb_soil) * 0.5_wp   
     1423            lambda_surface = lambda_soil
     1424
     1425          ELSE IF ( surf_t_surface%var_1d(m) >= surf_t_soil%var_2d(nzb_soil,m))&
     1426          THEN
     1427             lambda_surface = surf%lambda_surface_s(m) * lambda_soil           &
     1428                              / ( surf%lambda_surface_s(m) + lambda_soil )
     1429          ELSE
     1430
     1431             lambda_surface = surf%lambda_surface_u(m) * lambda_soil           &
     1432                              / ( surf%lambda_surface_u(m) + lambda_soil )
     1433          ENDIF
     1434       ELSE
    13981435          lambda_surface = surf%lambda_surface_s(m)
    1399           c_surface_tmp = surf%c_surface(m)
    1400          
    1401        ELSE
    1402 
    1403           lambda_surface = surf%lambda_surface_u(m)
    1404           c_surface_tmp = surf%c_surface(m)
    1405          
    1406        ENDIF
     1436       ENDIF
     1437
     1438!
     1439!--    Set heat capacity of the skin/surface. It is ususally zero when a skin
     1440!--    layer is used, and non-zero otherwise.
     1441       c_surface_tmp = surf%c_surface(m)
    14071442
    14081443!
     
    16861721                          * surf_t_surface%var_1d(m) - dq_s_dt *               &
    16871722                            surf_t_surface_p%var_1d(m) ) /                     &
    1688                             surf%qsws(m) - surf%r_a(m)
     1723                            (surf%qsws(m) + 1.0E-20) - surf%r_a(m)
    16891724       ENDIF
    16901725
     
    20692104       ENDDO
    20702105
     2106
     2107
    20712108!
    20722109!--    Calculate grid spacings. Temperature and moisture are defined at
    20732110!--    the center of the soil layers, whereas gradients/fluxes are
    20742111!--    defined at the edges (_layer)
     2112       zs(nzb_soil) = 0.5_wp * dz_soil(nzb_soil)
     2113       zs_layer(nzb_soil) = dz_soil(nzb_soil)
     2114
     2115       DO  k = nzb_soil+1, nzt_soil
     2116          zs_layer(k) = zs_layer(k-1) + dz_soil(k)
     2117          zs(k) = (zs_layer(k) +  zs_layer(k-1)) * 0.5_wp
     2118
     2119       ENDDO
     2120
     2121       dz_soil(nzt_soil+1) = zs_layer(nzt_soil) + dz_soil(nzt_soil)
     2122       zs(nzt_soil+1) = zs_layer(nzt_soil) + 0.5_wp * dz_soil(nzt_soil)
     2123 
    20752124       DO  k = nzb_soil, nzt_soil-1
    2076           dz_soil(k) = zs(k+1) - zs(k)
     2125          dz_soil_center(k) = zs(k+1) - zs(k)
    20772126         
    2078           IF ( dz_soil(k) == 0.0_wp )  THEN
     2127          IF ( dz_soil_center(k) == 0.0_wp )  THEN
    20792128             message_string = 'invalid soil layer configuration found ' //     &
    2080                               '(dz_soil(k) = 0.0)'
     2129                              '(dz_soil_center(k) = 0.0)'
    20812130             CALL message( 'lsm_read_restart_data', 'PA0140', 1, 2, 0, 6, 0 )
    2082           ENDIF
    2083          
     2131          ENDIF 
    20842132       ENDDO
    2085 
    2086        dz_soil_layer(nzb_soil) = 2.0_wp * zs(nzb_soil)
    2087        zs_layer(nzb_soil)      = 2.0_wp * zs(nzb_soil)
    20882133 
    2089        DO  k = nzb_soil+1, nzt_soil
    2090           dz_soil_layer(k) = ( zs(k) - zs_layer(k-1) ) * 2.0_wp
    2091           zs_layer(k) = zs_layer(k-1) + dz_soil_layer(k)
    2092        ENDDO
    2093        dz_soil_layer(nzt_soil+1) = 0.5_wp * dz_soil_layer(nzt_soil)
    2094        dz_soil(nzt_soil) = zs_layer(k-1) + dz_soil_layer(k) - zs(nzt_soil)
     2134       dz_soil_center(nzt_soil) = zs_layer(k-1) + dz_soil(k) - zs(nzt_soil)
     2135
     2136       IF ( myid == 0 )  THEN
     2137          PRINT*, "zs:", zs
     2138          PRINT*, "zs_layer:", zs_layer
     2139          PRINT*, "dz_soil:", dz_soil
     2140          PRINT*, "dz_soil_center:", dz_soil_center
     2141       ENDIF
    20952142
    20962143       
    2097        ddz_soil       = 1.0_wp / dz_soil
    2098        ddz_soil_layer = 1.0_wp / dz_soil_layer
     2144       ddz_soil_center       = 1.0_wp / dz_soil_center
     2145       ddz_soil = 1.0_wp / dz_soil
     2146
     2147
     2148
     2149
     2150
     2151
     2152!        DO  k = nzb_soil, nzt_soil-1
     2153!           dz_soil_center(k) = zs(k+1) - zs(k)
     2154!           
     2155!           IF ( dz_soil_center(k) == 0.0_wp )  THEN
     2156!              message_string = 'invalid soil layer configuration found ' //     &
     2157!                               '(dz_soil_center(k) = 0.0)'
     2158!              CALL message( 'lsm_read_restart_data', 'PA0140', 1, 2, 0, 6, 0 )
     2159!           ENDIF
     2160!           
     2161!        ENDDO
     2162!
     2163!        dz_soil(nzb_soil) = 2.0_wp * zs(nzb_soil)
     2164!        zs_layer(nzb_soil)      = 2.0_wp * zs(nzb_soil)
     2165!   
     2166!        DO  k = nzb_soil+1, nzt_soil
     2167!           dz_soil(k) = ( zs(k) - zs_layer(k-1) ) * 2.0_wp
     2168!           zs_layer(k) = zs_layer(k-1) + dz_soil(k)
     2169!        ENDDO
     2170!        dz_soil(nzt_soil+1) = 0.5_wp * dz_soil(nzt_soil)
     2171!        dz_soil_center(nzt_soil) = zs_layer(k-1) + dz_soil(k) - zs(nzt_soil)
     2172!
     2173!       
     2174!        ddz_soil_center       = 1.0_wp / dz_soil_center
     2175!        ddz_soil = 1.0_wp / dz_soil
     2176
     2177
    20992178
    21002179
     
    22082287          ENDIF   
    22092288         
    2210 !
    2211 !--       The heat conductivity between canopy and first soil layer must
    2212 !--       depend on the depth of the soil layer. The value in the lookup
    2213 !--       table refers to the first soil temperature at zs = -0.005 m. Lambda
    2214 !--       is thus interpolated to the actual soil layer depth used.
    22152289          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
    2216              lambda_surface_stable = vegetation_pars(6,vegetation_type)        &
    2217                                      * 0.005_wp * ddz_soil_layer(nzb_soil)     &
    2218                                      * 2.0_wp   
     2290             lambda_surface_stable = vegetation_pars(6,vegetation_type)
    22192291          ENDIF
    22202292
    22212293          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
    2222              lambda_surface_unstable = vegetation_pars(7,vegetation_type)      &
    2223                                        * 0.005_wp  * ddz_soil_layer(nzb_soil)  &
    2224                                        * 2.0_wp             
     2294             lambda_surface_unstable = vegetation_pars(7,vegetation_type)           
    22252295          ENDIF
    22262296
     
    22322302             c_surface = vegetation_pars(9,vegetation_type)       
    22332303          ENDIF
    2234              
     2304 
     2305          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
     2306             albedo_type = INT(vegetation_pars(10,vegetation_type))       
     2307          ENDIF
     2308   
     2309          IF ( emissivity == 9999999.9_wp )  THEN
     2310             emissivity = vegetation_pars(11,vegetation_type)     
     2311          ENDIF
     2312         
    22352313       ENDIF
    22362314
     
    22422320
    22432321          IF ( z0_water == 9999999.9_wp )  THEN
    2244              z0_water = water_pars(1,water_type)       
     2322             z0_water = water_pars(1,water_type)
    22452323          ENDIF       
    22462324
    22472325          IF ( z0h_water == 9999999.9_wp )  THEN
    2248              z0h_water = water_pars(2,water_type)       
     2326             z0h_water = water_pars(2,water_type)   
    22492327          ENDIF 
     2328
     2329          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
     2330             albedo_type = INT(water_pars(3,water_type))       
     2331          ENDIF
     2332   
     2333          IF ( emissivity == 9999999.9_wp )  THEN
     2334             emissivity = water_pars(4,water_type)       
     2335          ENDIF
    22502336
    22512337       ENDIF
     
    22562342             pavement_depth = pavement_pars(0,pavement_type)
    22572343          ELSE
    2258              pavement_depth = MAX( zs(nzb_soil), pavement_depth )
     2344             pavement_depth = MAX( 0.5_wp * dz_soil(nzb_soil), pavement_depth )
    22592345          ENDIF
    22602346
     
    22752361          ENDIF   
    22762362   
     2363          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
     2364             albedo_type = INT(pavement_pars(5,pavement_type))       
     2365          ENDIF
     2366   
     2367          IF ( emissivity == 9999999.9_wp )  THEN
     2368             emissivity = pavement_pars(6,pavement_type)       
     2369          ENDIF
     2370
    22772371       ENDIF
    22782372!
     
    23402434                t_soil_h%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
    23412435                surf_lsm_h%lambda_surface_s(m) = pavement_heat_conduct         &
    2342                                                  * ddz_soil_layer(nzb_soil)    &
     2436                                                 * ddz_soil(nzb_soil)    &
    23432437                                                 * 2.0_wp
    23442438                surf_lsm_h%lambda_surface_u(m) = surf_lsm_h%lambda_surface_s(m)                                         
    23452439                surf_lsm_h%c_surface(m)        = pavement_heat_capacity        &
    2346                                                  * dz_soil_layer(nzb_soil)     &
     2440                                                 * dz_soil(nzb_soil)     &
    23472441                                                 * 0.25_wp
    23482442                surf_lsm_h%lambda_h_def(m)    = pavement_heat_conduct                                     
     
    24282522                   t_soil_v(l)%var_2d(nzt_soil+1,m)  = soil_temperature(nzt_soil+1)
    24292523                   surf_lsm_v(l)%lambda_surface_s(m) = pavement_heat_conduct      &
    2430                                                        * ddz_soil_layer(nzb_soil) &
     2524                                                       * ddz_soil(nzb_soil) &
    24312525                                                       * 2.0_wp
    24322526                   surf_lsm_v(l)%lambda_surface_u(m) = surf_lsm_v(l)%lambda_surface_s(m)
    24332527                   surf_lsm_v(l)%c_surface(m) = pavement_heat_capacity         &
    2434                                                 * dz_soil_layer(nzb_soil)      &
     2528                                                * dz_soil(nzb_soil)      &
    24352529                                                * 0.25_wp
    24362530                   surf_lsm_v(l)%lambda_h_def(m)   = pavement_heat_conduct                                     
     
    26422736!
    26432737!--    Allocate global 1D arrays
    2644        ALLOCATE ( ddz_soil(nzb_soil:nzt_soil) )
    2645        ALLOCATE ( ddz_soil_layer(nzb_soil:nzt_soil+1) )
    2646        ALLOCATE ( dz_soil(nzb_soil:nzt_soil) )
    2647        ALLOCATE ( dz_soil_layer(nzb_soil:nzt_soil+1) )
     2738       ALLOCATE ( ddz_soil_center(nzb_soil:nzt_soil) )
     2739       ALLOCATE ( ddz_soil(nzb_soil:nzt_soil+1) )
     2740       ALLOCATE ( dz_soil_center(nzb_soil:nzt_soil) )
    26482741       ALLOCATE ( root_extr(nzb_soil:nzt_soil) )
    2649  
     2742       ALLOCATE ( zs(nzb_soil:nzt_soil+1) )
     2743
    26502744       root_extr = 0.0_wp
    26512745       
     
    28292923                                  constant_roughness,                          &
    28302924                                  conserve_water_content,                      &
     2925                                  dz_soil,                                     &
    28312926                                  f_shortwave_incoming, field_capacity,        &
    28322927                                  aero_resist_kray, hydraulic_conductivity,    &
     
    28432938                                  vegetation_coverage, vegetation_type,        &
    28442939                                  water_temperature, water_type,               &
    2845                                   wilting_point, zs, z0_vegetation,            &
     2940                                  wilting_point, z0_vegetation,                &
    28462941                                  z0h_vegetation, z0q_vegetation, z0_water,    &
    28472942                                  z0h_water, z0q_water, z0_pavement,           &
     
    29793074                     .AND.  zs(k+1) > pavement_depth )                         &
    29803075                THEN
    2981                    surf%lambda_h(k,m) = ( pavement_depth - zs(k) ) * ddz_soil(k+1)   &
     3076                   surf%lambda_h(k,m) = ( pavement_depth - zs(k) ) * ddz_soil_center(k+1)   &
    29823077                                      * lambda_temp(k)                               &
    2983                                       + ( zs(k+1) - pavement_depth ) * ddz_soil(k+1) &
     3078                                      + ( zs(k+1) - pavement_depth ) * ddz_soil_center(k+1) &
    29843079                                      * lambda_temp(k+1)
    29853080                ELSE
     
    29973092             tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) *            &
    29983093                    ( surf%lambda_h(nzb_soil,m) * ( surf_t_soil%var_2d(nzb_soil+1,m) &
    2999                       - surf_t_soil%var_2d(nzb_soil,m) ) * ddz_soil(nzb_soil)        &
    3000                       + surf%ghf(m) ) * ddz_soil_layer(nzb_soil)
     3094                      - surf_t_soil%var_2d(nzb_soil,m) ) * ddz_soil_center(nzb_soil)        &
     3095                      + surf%ghf(m) ) * ddz_soil(nzb_soil)
    30013096
    30023097             DO  k = nzb_soil+1, nzt_soil
     
    30043099                          * (   surf%lambda_h(k,m)                                      &
    30053100                              * ( surf_t_soil%var_2d(k+1,m) - surf_t_soil%var_2d(k,m) ) &
    3006                               * ddz_soil(k)                                             &
     3101                              * ddz_soil_center(k)                                             &
    30073102                              - surf%lambda_h(k-1,m)                                    &
    30083103                              * ( surf_t_soil%var_2d(k,m) - surf_t_soil%var_2d(k-1,m) ) &
    3009                               * ddz_soil(k-1)                                           &
    3010                             ) * ddz_soil_layer(k)
     3104                              * ddz_soil_center(k-1)                                           &
     3105                            ) * ddz_soil(k)
    30113106
    30123107             ENDDO
     
    31433238                tend(nzb_soil) = ( surf%lambda_w(nzb_soil,m) *   (       &
    31443239                         surf_m_soil%var_2d(nzb_soil+1,m) - surf_m_soil%var_2d(nzb_soil,m) )           &
    3145                          * ddz_soil(nzb_soil) - surf%gamma_w(nzb_soil,m) - &
     3240                         * ddz_soil_center(nzb_soil) - surf%gamma_w(nzb_soil,m) - &
    31463241                                                     (                         &
    31473242                            root_extr(nzb_soil) * surf%qsws_veg(m)    &
    31483243                            + surf%qsws_soil(m) ) * drho_l_lv )       &
    3149                             * ddz_soil_layer(nzb_soil)
     3244                            * ddz_soil(nzb_soil)
    31503245
    31513246                DO  k = nzb_soil+1, nzt_soil-1
    31523247                   tend(k) = ( surf%lambda_w(k,m) * ( surf_m_soil%var_2d(k+1,m)      &
    3153                              - surf_m_soil%var_2d(k,m) ) * ddz_soil(k)                   &
     3248                             - surf_m_soil%var_2d(k,m) ) * ddz_soil_center(k)                   &
    31543249                             - surf%gamma_w(k,m)                         &
    31553250                             - surf%lambda_w(k-1,m) * ( surf_m_soil%var_2d(k,m) -    &
    3156                              surf_m_soil%var_2d(k-1,m)) * ddz_soil(k-1)                      &
     3251                             surf_m_soil%var_2d(k-1,m)) * ddz_soil_center(k-1)                      &
    31573252                             + surf%gamma_w(k-1,m) - (root_extr(k)       &
    31583253                             * surf%qsws_veg(m) * drho_l_lv)          &
    3159                              ) * ddz_soil_layer(k)
     3254                             ) * ddz_soil(k)
    31603255                ENDDO
    31613256                tend(nzt_soil) = ( - surf%gamma_w(nzt_soil,m)            &
     
    31633258                                     * ( surf_m_soil%var_2d(nzt_soil,m)                    &
    31643259                                     - surf_m_soil%var_2d(nzt_soil-1,m))                   &
    3165                                      * ddz_soil(nzt_soil-1)                      &
     3260                                     * ddz_soil_center(nzt_soil-1)                      &
    31663261                                     + surf%gamma_w(nzt_soil-1,m) - (    &
    31673262                                       root_extr(nzt_soil)                     &
    31683263                                     * surf%qsws_veg(m) * drho_l_lv ) &
    3169                                   ) * ddz_soil_layer(nzt_soil)             
     3264                                  ) * ddz_soil(nzt_soil)             
    31703265
    31713266                surf_m_soil_p%var_2d(nzb_soil:nzt_soil,m) = surf_m_soil%var_2d(nzb_soil:nzt_soil,m)    &
Note: See TracChangeset for help on using the changeset viewer.