Changeset 2328 for palm/trunk


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

some improvements in land surface model

Location:
palm/trunk/SOURCE
Files:
2 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)    &
  • palm/trunk/SOURCE/radiation_model_mod.f90

    r2318 r2328  
    2525! -----------------
    2626! $Id$
     27! Emissivity can now be set individually for each pixel.
     28! Albedo type can be inferred from land surface model.
     29! Added default albedo type for bare soil
     30!
     31! 2318 2017-07-20 17:27:44Z suehring
    2732! Get topography top index via Function call
    2833!
     
    211216!
    212217!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
    213     CHARACTER(37), DIMENSION(0:17), PARAMETER :: albedo_type_name = (/      &
     218    CHARACTER(37), DIMENSION(0:18), PARAMETER :: albedo_type_name = (/      &
    214219                                   'user defined                         ', & !  0
    215220                                   'ocean                                ', & !  1
     
    229234                                   'sea ice                              ', & ! 15
    230235                                   'snow                                 ', & ! 16
    231                                    'pavement/roads                       '  & ! 17
     236                                   'pavement/roads                       ', & ! 17
     237                                   'bare soil                            '  & ! 18
    232238                                                         /)
    233239
    234     INTEGER(iwp) :: albedo_type  = 5,    & !< Albedo surface type (default: short grassland)
    235                     day,                 & !< current day of the year
    236                     day_init     = 172,  & !< day of the year at model start (21/06)
    237                     dots_rad     = 0       !< starting index for timeseries output
     240    INTEGER(iwp) :: albedo_type  = 9999999, & !< Albedo surface type
     241                    day,                    & !< current day of the year
     242                    day_init     = 172,     & !< day of the year at model start (21/06)
     243                    dots_rad     = 0          !< starting index for timeseries output
    238244
    239245    LOGICAL ::  unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
     
    261267                decl_3,                          & !< declination coef. 3
    262268                dt_radiation = 0.0_wp,           & !< radiation model timestep
    263                 emissivity = 0.98_wp,            & !< NAMELIST surface emissivity
     269                emissivity = 9999999.9_wp,       & !< NAMELIST surface emissivity
    264270                lambda = 0.0_wp,                 & !< longitude in degrees
    265271                lon = 0.0_wp,                    & !< longitude in radians
     
    278284    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
    279285                alpha,                       & !< surface broadband albedo (used for clear-sky scheme)
     286                emis,                        & !< surface broadband emissitivity
    280287                rad_lw_out_change_0,         & !< change in LW out due to change in surface temperature
    281288                rad_net,                     & !< net radiation at the surface
     
    285292!-- Land surface albedos for solar zenith angle of 60° after Briegleb (1992)     
    286293!-- (shortwave, longwave, broadband):   sw,      lw,      bb,
    287     REAL(wp), DIMENSION(0:2,1:17), PARAMETER :: albedo_pars = RESHAPE( (/&
     294    REAL(wp), DIMENSION(0:2,1:18), PARAMETER :: albedo_pars = RESHAPE( (/&
    288295                                   0.06_wp, 0.06_wp, 0.06_wp,            & !  1
    289296                                   0.09_wp, 0.28_wp, 0.19_wp,            & !  2
     
    302309                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 15
    303310                                   0.95_wp, 0.70_wp, 0.82_wp,            & ! 16
    304                                    0.08_wp, 0.08_wp, 0.08_wp             & ! 17
    305                                  /), (/ 3, 17 /) )
     311                                   0.08_wp, 0.08_wp, 0.08_wp,            & ! 17
     312                                   0.17_wp, 0.17_wp, 0.17_wp             & ! 18
     313                                 /), (/ 3, 18 /) )
    306314
    307315    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
     
    510518!
    511519!-- Public variables and constants / NEEDS SORTING
    512     PUBLIC decl_1, decl_2, decl_3, dots_rad, dt_radiation, force_radiation_call,&
     520    PUBLIC albedo, albedo_type, decl_1, decl_2, decl_3, dots_rad, dt_radiation, emissivity, force_radiation_call,&
    513521           lat, lon, rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in,        &
    514522           rad_lw_in_av, rad_lw_out, rad_lw_out_av, rad_lw_out_change_0,       &
     
    888896
    889897!
     898!--    Allocate array for storing emissivity
     899       IF ( .NOT. ALLOCATED ( emis ) )  THEN
     900          ALLOCATE ( emis(nysg:nyng,nxlg:nxrg) )
     901          emis = emissivity
     902       ENDIF
     903
     904!
    890905!--    Allocate array for storing the surface net radiation
    891906       IF ( .NOT. ALLOCATED ( rad_net ) )  THEN
     
    956971!
    957972!--       Overwrite albedo if manually set in parameter file
    958           IF ( albedo_type /= 0 .AND. albedo == 9999999.9_wp )  THEN
     973          IF ( albedo_type /= 0  .AND.  albedo_type /= 9999999 .AND. albedo == 9999999.9_wp )  THEN
    959974             albedo = albedo_pars(2,albedo_type)
    960975          ENDIF
    961    
     976!
     977!--       Write albedo to 2d array alpha to allow surface heterogeneities   
    962978          alpha = albedo
    963979 
     
    11941210             rad_sw_in(0,j,i)  = solar_constant * sky_trans * zenith(0)
    11951211             rad_sw_out(0,j,i) = alpha(j,i) * rad_sw_in(0,j,i)
    1196              rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
     1212             rad_lw_out(0,j,i) = emis(j,i) * sigma_sb * (pt(k,j,i) * exn)**4
    11971213
    11981214             IF ( cloud_physics )  THEN
     
    12071223
    12081224
    1209              rad_lw_out_change_0(j,i) = 3.0_wp * sigma_sb * emissivity         &
     1225             rad_lw_out_change_0(j,i) = 3.0_wp * sigma_sb * emis(j,i)         &
    12101226                                        * (pt(k,j,i) * exn) ** 3
    12111227
     
    12541270             ENDIF
    12551271
    1256              rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
     1272             rad_lw_out(0,j,i) = emis(j,i) * sigma_sb * (pt(k,j,i) * exn)**4
    12571273
    12581274             rad_sw_in(0,j,i) = ( rad_net(j,i) - rad_lw_in(0,j,i)              &
     
    15331549             rrtm_tsfc = pt(nzb,j,i) * (surface_pressure / 1000.0_wp )**0.286_wp
    15341550
     1551!
     1552!--          Set surface emissivity
     1553             rrtm_emis = emis(j,i)
     1554
    15351555             IF ( lw_radiation )  THEN
    15361556               CALL rrtmg_lw( 1, nzt_rad      , rrtm_icld    , rrtm_idrv      ,&
     
    17231743                 rrtm_aldir(0,:,:) = aldif
    17241744                 rrtm_asdir(0,:,:) = asdif
     1745
     1746!
     1747!--        Bare soil
     1748           ELSEIF ( albedo_type == 18 )  THEN
     1749                 rrtm_aldir(0,:,:) = aldif
     1750                 rrtm_asdir(0,:,:) = asdif
     1751
    17251752!
    17261753!--        Land surfaces
Note: See TracChangeset for help on using the changeset viewer.