MODULE land_surface_model_mod !--------------------------------------------------------------------------------! ! This file is part of PALM. ! ! PALM is free software: you can redistribute it and/or modify it under the terms ! of the GNU General Public License as published by the Free Software Foundation, ! either version 3 of the License, or (at your option) any later version. ! ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License along with ! PALM. If not, see . ! ! Copyright 1997-2014 Leibniz Universitaet Hannover !--------------------------------------------------------------------------------! ! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: land_surface_model.f90 1514 2014-12-19 09:14:55Z maronga $ ! Bugfix: REAL constants provided with KIND-attribute in call of ! intrinsic function MAX and MIN ! ! 1500 2014-12-03 17:42:41Z maronga ! Corrected calculation of aerodynamic resistance (r_a). ! Precipitation is now added to liquid water reservoir using LE_liq. ! Added support for dry runs. ! ! 1496 2014-12-02 17:25:50Z maronga ! Initial revision ! ! ! Description: ! ------------ ! Land surface model, consisting of a solver for the energy balance at the ! surface and a four layer soil scheme. The scheme is similar to the TESSEL ! scheme implemented in the ECMWF IFS model, with modifications according to ! H-TESSEL. The implementation is based on the formulation implemented in the ! DALES model. ! ! To do list: ! ----------- ! - Add support for binary I/O support ! - Add support for lsm data output ! - Check for time step criterion ! - Check use with RK-2 and Euler time-stepping ! - Adaption for use with cloud physics (liquid water potential temperature) ! - Check reaction of plants at wilting point and at atmospheric saturation ! - Consider partial absorption of the net shortwave radiation by the skin layer ! - Allow for water surfaces, check performance for bare soils !------------------------------------------------------------------------------! USE arrays_3d, & ONLY: pt, pt_p, q, q_p, qsws, rif, shf, ts, us, z0, z0h USE cloud_parameters, & ONLY: cp, l_d_r, l_v, precipitation_rate, rho_l, r_d, r_v USE control_parameters, & ONLY: dt_3d, humidity, intermediate_timestep_count, & intermediate_timestep_count_max, precipitation, pt_surface, & rho_surface, surface_pressure, timestep_scheme, tsc USE indices, & ONLY: nxlg, nxrg, nyng, nysg, nzb_s_inner USE kinds USE radiation_model_mod, & ONLY: Rn, SW_in, sigma_SB IMPLICIT NONE ! !-- LSM model constants INTEGER(iwp), PARAMETER :: soil_layers = 4 !: number of soil layers (fixed for now) REAL(wp), PARAMETER :: & b_CH = 6.04_wp, & ! Clapp & Hornberger exponent lambda_h_dry = 0.19_wp, & ! heat conductivity for dry soil lambda_h_sm = 3.44_wp, & ! heat conductivity of the soil matrix lambda_h_water = 0.57_wp, & ! heat conductivity of water psi_sat = -0.388_wp, & ! soil matrix potential at saturation rhoC_soil = 2.19E6_wp, & ! volumetric heat capacity of soil rhoC_water = 4.20E6_wp, & ! volumetric heat capacity of water m_max_depth = 0.0002_wp ! Maximum capacity of the water reservoir (m) ! !-- LSM variables INTEGER(iwp) :: veg_type = 2, & !: vegetation type, 0: user-defined, 1-19: generic (see list) soil_type = 3 !: soil type, 0: user-defined, 1-6: generic (see list) LOGICAL :: conserve_water_content = .TRUE., & !: open or closed bottom surface for the soil model land_surface = .FALSE. !: flag parameter indicating wheather the lsm is used ! value 9999999.9_wp -> generic available or user-defined value must be set ! otherwise -> no generic variable and user setting is optional REAL(wp) :: alpha_VanGenuchten = 0.0_wp, & !: NAMELIST alpha_VG canopy_resistance_coefficient = 0.0_wp, & !: NAMELIST gD C_skin = 20000.0_wp, & !: Skin heat capacity drho_l_lv, & !: (rho_l * l_v)**-1 exn, & !: value of the Exner function e_s = 0.0_wp, & !: saturation water vapour pressure field_capacity = 0.0_wp, & !: NAMELIST m_fc f_shortwave_incoming = 9999999.9_wp, & !: NAMELIST f_SW_in hydraulic_conductivity = 0.0_wp, & !: NAMELIST gamma_w_sat Ke = 0.0_wp, & !: Kersten number lambda_skin_stable = 9999999.9_wp, & !: NAMELIST lambda_skin_s lambda_skin_unstable = 9999999.9_wp, & !: NAMELIST lambda_skin_u leaf_area_index = 9999999.9_wp, & !: NAMELIST LAI l_VanGenuchten = 0.0_WP, & !: NAMELIST l_VG min_canopy_resistance = 110.0_wp, & !: NAMELIST r_s_min m_total = 0.0_wp, & !: weighed total water content of the soil (m3/m3) n_VanGenuchten = 0.0_WP, & !: NAMELIST n_VG q_s = 0.0_wp, & !: saturation specific humidity residual_moisture = 0.0_wp, & !: NAMELIST m_res rho_cp, & !: rho_surface * cp rho_lv, & !: rho * l_v rd_d_rv, & !: r_d / r_v saturation_moisture = 0.0_wp, & !: NAMELIST m_sat vegetation_coverage = 9999999.9_wp, & !: NAMELIST c_veg wilting_point = 0.0_wp !: NAMELIST m_wilt REAL(wp), DIMENSION(0:soil_layers-1) :: & ddz_soil, & !: 1/dz_soil ddz_soil_stag, & !: 1/dz_soil_stag dz_soil, & !: soil grid spacing (center-center) dz_soil_stag, & !: soil grid spacing (edge-edge) root_extr = 0.0_wp, & !: root extraction root_fraction = (/0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp/), & !: distribution of root surface area to the individual soil layers soil_level = (/0.07_wp, 0.28_wp, 1.00_wp, 2.89_wp/), & !: soil layer depths (m) soil_moisture = 0.0_wp !: soil moisture content (m3/m3) REAL(wp), DIMENSION(0:soil_layers) :: & soil_temperature = 9999999.9_wp !: soil temperature (K) #if defined( __nopointer ) REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: T_0, & !: skin temperature (K) T_0_p, & !: progn. skin temperature (K) m_liq, & !: liquid water reservoir (m) m_liq_p !: progn. liquid water reservoir (m) #else REAL(wp), DIMENSION(:,:), POINTER :: T_0, & T_0_p, & m_liq, & m_liq_p REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: T_0_1, T_0_2, & m_liq_1, m_liq_2 #endif ! !-- Temporal tendencies for time stepping REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tT_0_m, & !: skin temperature tendency (K) tm_liq_m !: liquid water reservoir tendency (m) ! !-- Energy balance variables REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & alpha_VG, & !: coef. of Van Genuchten c_liq, & !: liquid water coverage (of vegetated area) c_veg, & !: vegetation coverage f_SW_in, & !: ? G, & !: surface soil heat flux H, & !: surface flux of sensible heat gamma_w_sat, & !: hydraulic conductivity at saturation gD, & !: coefficient for dependence of r_canopy on water vapour pressure deficit LAI, & !: leaf area index LE, & !: surface flux of latent heat (total) LE_veg, & !: surface flux of latent heat (vegetation portion) LE_soil, & !: surface flux of latent heat (soil portion) LE_liq, & !: surface flux of latent heat (liquid water portion) lambda_h_sat, & !: heat conductivity for dry soil lambda_skin_s, & !: coupling between skin and soil (depends on vegetation type) lambda_skin_u, & !: coupling between skin and soil (depends on vegetation type) l_VG, & !: coef. of Van Genuchten m_fc, & !: soil moisture at field capacity (m3/m3) m_res, & !: residual soil moisture m_sat, & !: saturation soil moisture (m3/m3) m_wilt, & !: soil moisture at permanent wilting point (m3/m3) n_VG, & !: coef. Van Genuchten r_a, & !: aerodynamic resistance r_canopy, & !: canopy resistance r_soil, & !: soil resitance r_soil_min, & !: minimum soil resistance r_s, & !: total surface resistance (combination of r_soil and r_canopy) r_s_min !: minimum canopy resistance REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & lambda_h, & !: heat conductivity of soil (?) lambda_w, & !: hydraulic diffusivity of soil (?) gamma_w, & !: hydraulic conductivity of soil (?) rhoC_total !: volumetric heat capacity of the actual soil matrix (?) #if defined( __nopointer ) REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: & T_soil, & !: Soil temperature (K) T_soil_p, & !: Prog. soil temperature (K) m_soil, & !: Soil moisture (m3/m3) m_soil_p !: Prog. soil moisture (m3/m3) #else REAL(wp), DIMENSION(:,:,:), POINTER :: & T_soil, T_soil_p, & m_soil, m_soil_p REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: & T_soil_1, T_soil_2, & m_soil_1, m_soil_2 #endif REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & tT_soil_m, & !: T_soil storage array tm_soil_m, & !: m_soil storage array root_fr !: root fraction (sum=1) ! !-- Land surface parameters according to the following classes (veg_type) !-- (0 user defined) !-- 1 crops, mixed farming !-- 2 short grass !-- 3 evergreen needleleaf trees !-- 4 deciduous needleleaf trees !-- 5 evergreen broadleaf trees !-- 6 deciduous broadleaf trees !-- 7 tall grass !-- 8 desert !-- 9 tundra !-- 10 irrigated crops !-- 11 semidesert !-- 12 ice caps and glaciers !-- 13 bogs and marshes !-- 14 inland water !-- 15 ocean !-- 16 evergreen shrubs !-- 17 deciduous shrubs !-- 18 mixed forest/woodland !-- 19 interrupted forest ! !-- Land surface parameters I r_s_min, LAI, c_veg, gD REAL(wp), DIMENSION(0:3,1:19) :: veg_pars = RESHAPE( (/ & 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & ! 1 110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp, & ! 2 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 3 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 4 175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 5 240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp, & ! 6 100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp, & ! 7 250.0_wp, 0.50_wp, 0.00_wp, 0.00_wp, & ! 8 80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp, & ! 9 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & ! 10 150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp, & ! 11 0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 12 240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp, & ! 13 0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 14 0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 15 225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp, & ! 16 225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp, & ! 17 250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 18 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp & ! 19 /), (/ 4, 19 /) ) ! !-- Land surface parameters II z0, z0h REAL(wp), DIMENSION(0:1,1:19) :: roughness_par = RESHAPE( (/ & 0.25_wp, 0.25E-2_wp, & ! 1 0.20_wp, 0.20E-2_wp, & ! 2 2.00_wp, 2.00_wp, & ! 3 2.00_wp, 2.00_wp, & ! 4 2.00_wp, 2.00_wp, & ! 5 2.00_wp, 2.00_wp, & ! 6 0.47_wp, 0.47E-2_wp, & ! 7 0.013_wp, 0.013E-2_wp, & ! 8 0.034_wp, 0.034E-2_wp, & ! 9 0.5_wp, 0.50E-2_wp, & ! 10 0.17_wp, 0.17E-2_wp, & ! 11 1.3E-3_wp, 1.3E-4_wp, & ! 12 0.83_wp, 0.83E-2_wp, & ! 13 0.00_wp, 0.00E-2_wp, & ! 14 0.00_wp, 0.00E-2_wp, & ! 15 0.10_wp, 0.10E-2_wp, & ! 16 0.25_wp, 0.25E-2_wp, & ! 17 2.00_wp, 2.00E-2_wp, & ! 18 1.10_wp, 1.10E-2_wp & ! 19 /), (/ 2, 19 /) ) ! !-- Land surface parameters III lambda_skin_s, lambda_skin_u, f_SW_in REAL(wp), DIMENSION(0:2,1:19) :: skin_pars = RESHAPE( (/ & 10.0_wp, 10.0_wp, 0.05_wp, & ! 1 10.0_wp, 10.0_wp, 0.05_wp, & ! 2 20.0_wp, 15.0_wp, 0.03_wp, & ! 3 20.0_wp, 15.0_wp, 0.03_wp, & ! 4 20.0_wp, 15.0_wp, 0.03_wp, & ! 5 20.0_wp, 15.0_wp, 0.03_wp, & ! 6 10.0_wp, 10.0_wp, 0.05_wp, & ! 7 15.0_wp, 15.0_wp, 0.00_wp, & ! 8 10.0_wp, 10.0_wp, 0.05_wp, & ! 9 10.0_wp, 10.0_wp, 0.05_wp, & ! 10 10.0_wp, 10.0_wp, 0.05_wp, & ! 11 58.0_wp, 58.0_wp, 0.00_wp, & ! 12 10.0_wp, 10.0_wp, 0.05_wp, & ! 13 1.0E20_wp, 1.0E20_wp, 0.00_wp, & ! 14 1.0E20_wp, 1.0E20_wp, 0.00_wp, & ! 15 10.0_wp, 10.0_wp, 0.05_wp, & ! 16 10.0_wp, 10.0_wp, 0.05_wp, & ! 17 20.0_wp, 15.0_wp, 0.03_wp, & ! 18 20.0_wp, 15.0_wp, 0.03_wp & ! 19 /), (/ 3, 19 /) ) ! !-- Root distribution (sum = 1) level 1, level 2, level 3, level 4, REAL(wp), DIMENSION(0:3,1:19) :: root_distribution = RESHAPE( (/ & 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 1 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp, & ! 2 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp, & ! 3 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp, & ! 4 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp, & ! 5 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp, & ! 6 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp, & ! 7 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 8 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp, & ! 9 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 10 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp, & ! 11 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 12 0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp, & ! 13 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 14 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 15 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 16 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 17 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp, & ! 18 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp & ! 19 /), (/ 4, 19 /) ) ! !-- Soil parameters according to the following porosity classes (soil_type) !-- (0 user defined) !-- 1 coarse !-- 2 medium !-- 3 medium-fine !-- 4 fine !-- 5 very fine !-- 6 organic ! !-- Soil parameters I alpha_VG, l_VG, n_VG, gamma_w_sat REAL(wp), DIMENSION(0:3,1:6) :: soil_pars = RESHAPE( (/ & 3.83_wp, 1.250_wp, 1.38_wp, 6.94E-6_wp, & ! 1 3.14_wp, -2.342_wp, 1.28_wp, 1.16E-6_wp, & ! 2 0.83_wp, -0.588_wp, 1.25_wp, 0.26E-6_wp, & ! 3 3.67_wp, -1.977_wp, 1.10_wp, 2.87E-6_wp, & ! 4 2.65_wp, 2.500_wp, 1.10_wp, 1.74E-6_wp, & ! 5 1.30_wp, 0.400_wp, 1.20_wp, 0.93E-6_wp & ! 6 /), (/ 4, 6 /) ) ! !-- Soil parameters II m_sat, m_fc, m_wilt, m_res REAL(wp), DIMENSION(0:3,1:6) :: m_soil_pars = RESHAPE( (/ & 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & ! 1 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & ! 2 0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp, & ! 3 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, & ! 4 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, & ! 5 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp & ! 6 /), (/ 4, 6 /) ) SAVE PRIVATE PUBLIC alpha_VanGenuchten, C_skin, canopy_resistance_coefficient, & conserve_water_content, field_capacity, f_shortwave_incoming, & hydraulic_conductivity, init_lsm, lambda_skin_stable, & lambda_skin_unstable, land_surface, leaf_area_index, & lsm_energy_balance, lsm_soil_model, l_VanGenuchten, & min_canopy_resistance, n_VanGenuchten, residual_moisture, & root_fraction, saturation_moisture, soil_level, soil_moisture, & soil_temperature, soil_type, vegetation_coverage, veg_type, & wilting_point #if defined( __nopointer ) PUBLIC m_liq, m_liq_p, m_soil, m_soil_p, T_0, T_0_p, T_soil, T_soil_p #else PUBLIC m_liq, m_liq_1, m_liq_2, m_liq_p, m_soil, m_soil_1, m_soil_2, & m_soil_p, T_0, T_0_1, T_0_2, T_0_p, T_soil, T_soil_1, T_soil_2, & T_soil_p #endif INTERFACE init_lsm MODULE PROCEDURE init_lsm END INTERFACE init_lsm INTERFACE lsm_energy_balance MODULE PROCEDURE lsm_energy_balance END INTERFACE lsm_energy_balance INTERFACE lsm_soil_model MODULE PROCEDURE lsm_soil_model END INTERFACE lsm_soil_model CONTAINS !------------------------------------------------------------------------------! ! Description: ! ------------ !-- Initialization of the land surface model !------------------------------------------------------------------------------! SUBROUTINE init_lsm IMPLICIT NONE INTEGER(iwp) :: i !: running index INTEGER(iwp) :: j !: running index INTEGER(iwp) :: k !: running index ! !-- Calculate frequently used parameters rho_cp = cp * rho_surface rd_d_rv = r_d / r_v rho_lv = rho_surface * l_v drho_l_lv = 1.0 / (rho_l * l_v) ! !-- Allocate skin and soil temperature / humidity #if defined( __nopointer ) ALLOCATE ( T_0(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( T_0_p(nysg:nyng,nxlg:nxrg) ) #else ALLOCATE ( T_0_1(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( T_0_2(nysg:nyng,nxlg:nxrg) ) #endif ALLOCATE ( tT_0_m(nysg:nyng,nxlg:nxrg) ) #if defined( __nopointer ) ALLOCATE ( T_soil(0:soil_layers,nysg:nyng,nxlg:nxrg) ) ALLOCATE ( T_soil_p(0:soil_layers,nysg:nyng,nxlg:nxrg) ) #else ALLOCATE ( T_soil_1(0:soil_layers,nysg:nyng,nxlg:nxrg) ) ALLOCATE ( T_soil_2(0:soil_layers,nysg:nyng,nxlg:nxrg) ) #endif ALLOCATE ( tT_soil_m(0:soil_layers-1,nysg:nyng,nxlg:nxrg) ) #if defined( __nopointer ) ALLOCATE ( m_liq(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( m_liq_p(nysg:nyng,nxlg:nxrg) ) #else ALLOCATE ( m_liq_1(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( m_liq_2(nysg:nyng,nxlg:nxrg) ) #endif ALLOCATE ( tm_liq_m(nysg:nyng,nxlg:nxrg) ) #if defined( __nopointer ) ALLOCATE ( m_soil(0:soil_layers-1,nysg:nyng,nxlg:nxrg) ) ALLOCATE ( m_soil_p(0:soil_layers-1,nysg:nyng,nxlg:nxrg) ) #else ALLOCATE ( m_soil_1(0:soil_layers-1,nysg:nyng,nxlg:nxrg) ) ALLOCATE ( m_soil_2(0:soil_layers-1,nysg:nyng,nxlg:nxrg) ) #endif ALLOCATE ( tm_soil_m(0:soil_layers-1,nysg:nyng,nxlg:nxrg) ) #if ! defined( __nopointer ) ! !-- Initial assignment of the pointers T_soil => T_soil_1; T_soil_p => T_soil_2 T_0 => T_0_1; T_0_p => T_0_2 m_soil => m_soil_1; m_soil_p => m_soil_2 m_liq => m_liq_1; m_liq_p => m_liq_2 #endif T_0 = 0.0_wp T_0_p = 0.0_wp tT_0_m = 0.0_wp T_soil = 0.0_wp T_soil_p = 0.0_wp tT_soil_m = 0.0_wp m_liq = 0.0_wp m_liq_p = 0.0_wp tm_liq_m = 0.0_wp m_soil = 0.0_wp m_soil_p = 0.0_wp tm_soil_m = 0.0_wp ! !-- Allocate 2D vegetation model arrays ALLOCATE ( alpha_VG(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( f_SW_in(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( G(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( H(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( gD(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( LAI(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( LE(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( LE_veg(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( LE_soil(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( LE_liq(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( l_VG(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( lambda_h_sat(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( lambda_skin_u(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( lambda_skin_s(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( n_VG(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) ) ALLOCATE ( r_s_min(nysg:nyng,nxlg:nxrg) ) ! !-- Set initial and default values c_liq = 0.0_wp c_veg = 0.0_wp f_SW_in = 0.05_wp gD = 0.0_wp LAI = 0.0_wp lambda_skin_u = 10.0_wp lambda_skin_s = 10.0_wp G = 0.0_wp H = rho_cp * shf IF ( humidity ) THEN LE = rho_l * l_v * qsws ELSE LE = 0.0_wp ENDIF LE_veg = 0.0_wp LE_soil = LE LE_liq = 0.0_wp r_a = 50.0_wp r_canopy = 0.0_wp r_soil = 0.0_wp r_soil_min = 50.0_wp r_s = 110.0_wp r_s_min = min_canopy_resistance ! !-- Allocate 3D soil model arrays ALLOCATE ( root_fr(0:soil_layers-1,nysg:nyng,nxlg:nxrg) ) ALLOCATE ( lambda_h(0:soil_layers-1,nysg:nyng,nxlg:nxrg) ) ALLOCATE ( rhoC_total(0:soil_layers-1,nysg:nyng,nxlg:nxrg) ) lambda_h = 0.0_wp ! !-- If required, allocate humidity-related variables for the soil model IF ( humidity ) THEN ALLOCATE ( lambda_w(0:soil_layers-1,nysg:nyng,nxlg:nxrg) ) ALLOCATE ( gamma_w(0:soil_layers-1,nysg:nyng,nxlg:nxrg) ) lambda_w = 0.0_wp ENDIF ! !-- Calculate grid spacings. Temperature and moisture are defined at !-- the center of the soil layers, whereas gradients/fluxes are defined !-- at the edges (_stag) dz_soil_stag(0) = soil_level(0) DO k = 1, soil_layers-1 dz_soil_stag(k) = soil_level(k) - soil_level(k-1) ENDDO DO k = 0, soil_layers-2 dz_soil(k) = 0.5 * (dz_soil_stag(k+1) + dz_soil_stag(k)) ENDDO dz_soil(soil_layers-1) = dz_soil_stag(soil_layers-1) ddz_soil = 1.0 / dz_soil ddz_soil_stag = 1.0 / dz_soil_stag ! !-- Initialize soil IF ( soil_type .NE. 0 ) THEN alpha_VG = soil_pars(0,soil_type) l_VG = soil_pars(1,soil_type) n_VG = soil_pars(2,soil_type) gamma_w_sat = soil_pars(3,soil_type) m_sat = m_soil_pars(0,soil_type) m_fc = m_soil_pars(1,soil_type) m_wilt = m_soil_pars(2,soil_type) m_res = m_soil_pars(3,soil_type) ELSE alpha_VG = alpha_VanGenuchten l_VG = l_VanGenuchten n_VG = n_VanGenuchten gamma_w_sat = hydraulic_conductivity m_sat = saturation_moisture m_fc = field_capacity m_wilt = wilting_point m_res = residual_moisture ENDIF ! !-- Map user settings of T and q for each soil layer !-- (make sure that the soil moisture does not drop below the permanent !-- wilting point) -> problems with devision by zero) DO k = 0, soil_layers-1 T_soil(k,:,:) = soil_temperature(k) m_soil(k,:,:) = MAX(soil_moisture(k),m_wilt(:,:)) ENDDO T_soil(soil_layers,:,:) = soil_temperature(soil_layers) exn = ( surface_pressure / 1000.0_wp )**0.286_wp T_0 = pt_surface * exn T_soil_p = T_soil m_soil_p = m_soil ! !-- Calculate saturation soil heat conductivity lambda_h_sat(:,:) = lambda_h_sm ** (1.0_wp - m_sat(:,:)) * & lambda_h_water ** m_sat(:,:) ! !-- Initialize vegetation IF ( veg_type .NE. 0 ) THEN r_s_min = veg_pars(0,veg_type) LAI = veg_pars(1,veg_type) c_veg = veg_pars(2,veg_type) gD = veg_pars(3,veg_type) lambda_skin_s = skin_pars(0,veg_type) lambda_skin_u = skin_pars(1,veg_type) f_SW_in = skin_pars(2,veg_type) z0 = roughness_par(0,veg_type) z0h = roughness_par(1,veg_type) DO k = 0, soil_layers-1 root_fr(k,:,:) = root_distribution(k,veg_type) ENDDO ELSE DO k = 0, soil_layers-1 root_fr(k,:,:) = root_fraction(k) ENDDO ENDIF ! !-- Possibly do user-defined actions (e.g. define heterogeneous land surface) CALL user_init_land_surface ! !-- Set artifical values for ts and us so that r_a has its initial value for !-- the first time step DO i = nxlg, nxrg DO j = nysg, nyng k = nzb_s_inner(j,i) ! !-- Assure that r_a cannot be zero at model start IF ( pt(k+1,j,i) == pt(k,j,i) ) pt(k+1,j,i) = pt(k+1,j,i) + 1.0E-10_wp us(j,i) = 0.1_wp ts(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / r_a(j,i) shf(j,i) = - us(j,i) * ts(j,i) ENDDO ENDDO ! !-- Calculate humidity at the surface IF ( humidity ) THEN CALL calc_q0 ENDIF RETURN END SUBROUTINE init_lsm !------------------------------------------------------------------------------! ! Description: ! ------------ ! !------------------------------------------------------------------------------! SUBROUTINE lsm_energy_balance IMPLICIT NONE INTEGER(iwp) :: i !: running index INTEGER(iwp) :: j !: running index INTEGER(iwp) :: k, ks !: running index REAL(wp) :: f1, & !: resistance correction term 1 f2, & !: resistance correction term 2 f3, & !: resistance correction term 3 m_min, & !: minimum soil moisture T_1, & !: actual temperature at first grid point e, & !: water vapour pressure e_s, & !: water vapour saturation pressure e_s_dT, & !: derivate of e_s with respect to T tend, & !: tendency dq_s_dT, & !: derivate of q_s with respect to T coef_1, & !: coef. for prognostic equation coef_2, & !: coef. for prognostic equation f_LE, & !: factor for LE f_LE_veg, & !: factor for LE_veg f_LE_soil, & !: factor for LE_soil f_LE_liq, & !: factor for LE_liq f_H, & !: factor for H lambda_skin, & !: Current value of lambda_skin m_liq_max !: maxmimum value of the liquid water reservoir ! !-- Calculate the exner function for the current time step exn = ( surface_pressure / 1000.0_wp )**0.286_wp DO i = nxlg, nxrg DO j = nysg, nyng k = nzb_s_inner(j,i) ! !-- Set lambda_skin according to stratification IF ( rif(j,i) >= 0.0_wp ) THEN lambda_skin = lambda_skin_s(j,i) ELSE lambda_skin = lambda_skin_u(j,i) ENDIF ! !-- First step: calculate aerodyamic resistance. As pt, us, ts !-- are not available for the prognostic time step, data from the last !-- time step is used here. Note that this formulation is the !-- equivalent to the ECMWF formulation using drag coefficients r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20) ! !-- Second step: calculate canopy resistance r_canopy !-- f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation !-- f1: correction for incoming shortwave radiation f1 = MIN(1.0_wp, ( 0.004_wp * SW_in(j,i) + 0.05_wp ) / & (0.81_wp * (0.004_wp * SW_in(j,i) + 1.0_wp) ) ) ! !-- f2: correction for soil moisture f2=0 for very dry soil m_total = 0.0_wp DO ks = 0, soil_layers-1 m_total = m_total + root_fr(ks,j,i) * m_soil(ks,j,i) ENDDO IF ( m_total .GT. m_wilt(j,i) .AND. m_total .LE. m_fc(j,i) ) THEN f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) ) ELSE f2 = 1.0E-20_wp ENDIF ! !-- Calculate water vapour pressure at saturation !-- (T_0 should be replaced by liquid water temp?!) e_s = 0.01 * 610.78_wp * EXP( 17.269_wp * ( T_0(j,i) - 273.16_wp )& / ( T_0(j,i) - 35.86_wp ) ) ! !-- f3: correction for vapour pressure deficit IF ( gD(j,i) .NE. 0.0_wp ) THEN ! !-- Calculate vapour pressure e = q_p(k+1,j,i) * surface_pressure / 0.622 f3 = EXP ( -gD(j,i) * (e_s - e) ) ELSE f3 = 1.0_wp ENDIF ! !-- To do: check for very dry soil -> r_canopy goes to infinity r_canopy(j,i) = r_s_min(j,i) / (LAI(j,i) * f1 * f2 * f3 + 1.0E-20) ! !-- Third step: calculate bare soil resistance r_soil m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) * & m_res(j,i) f2 = ( m_soil(0,j,i) - m_min ) / ( m_fc(j,i) - m_min ) f2 = MAX(f2,1.0E-20_wp) r_soil(j,i) = r_soil_min(j,i) / f2 ! !-- Calculate fraction of liquid water reservoir m_liq_max = m_max_depth * LAI(j,i) c_liq(j,i) = MIN(1.0_wp, m_liq(j,i)/m_liq_max) q_s = 0.622_wp * e_s / surface_pressure ! !-- In case of dew fall, set resistances to zero. !-- To do: what does that physically reasoning behind this? IF ( humidity ) THEN IF ( q_s .LE. q_p(k+1,j,i) ) THEN r_canopy(j,i) = 0.0_wp r_soil(j,i) = 0.0_wp ENDIF ENDIF ! !-- Calculate coefficients for the total evapotranspiration f_LE_veg = rho_lv * c_veg(j,i) * (1.0 - c_liq(j,i)) / (r_a(j,i) & + r_canopy(j,i)) f_LE_soil = rho_lv * (1.0 - c_veg(j,i)) / (r_a(j,i) + r_soil(j,i)) f_LE_liq = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i) ! !-- If soil moisture is below wilting point, plants do no longer !-- transpirate. IF ( m_soil(k,j,i) .LT. m_wilt(j,i) ) THEN f_LE_veg = 0.0 ENDIF f_H = rho_cp / r_a(j,i) f_LE = f_LE_veg + f_LE_soil + f_LE_liq ! !-- Calculate derivative of q_s for Taylor series expansion e_s_dT = e_s * ( 17.269_wp / (T_0(j,i) - 35.86_wp) - & 17.269_wp*(T_0(j,i) - 273.16_wp) / (T_0(j,i) & - 35.86_wp)**2 ) dq_s_dT = 0.622_wp * e_s_dT / surface_pressure T_1 = pt_p(k+1,j,i) * exn ! !-- Add LW up so that it can be removed in prognostic equation Rn(j,i) = Rn(j,i) + sigma_SB * T_0(j,i) ** 4 IF ( humidity ) THEN ! !-- Numerator of the prognostic equation coef_1 = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i) ** 4 + f_H & / exn * T_1 + f_LE * ( q_p(k+1,j,i) - q_s + dq_s_dT & * T_0(j,i) ) + lambda_skin * T_soil(0,j,i) ! !-- Denominator of the prognostic equation coef_2 = 4.0_wp * sigma_SB * T_0(j,i) ** 3 + f_LE * dq_s_dT & + lambda_skin + f_H / exn ELSE ! !-- Numerator of the prognostic equation coef_1 = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i) ** 4 + f_H / & exn * T_1 + lambda_skin * T_soil(0,j,i) ! !-- Denominator of the prognostic equation coef_2 = 4.0_wp * sigma_SB * T_0(j,i) ** 3 & + lambda_skin + f_H / exn ENDIF tend = 0.0_wp ! !-- Implicit solution when the skin layer has no heat capacity, !-- otherwise use RK3 scheme. T_0_p(j,i) = ( coef_1 * dt_3d * tsc(2) + C_skin * T_0(j,i) ) / & ( C_skin + coef_2 * dt_3d * tsc(2) ) ! !-- Add RK3 term T_0_p(j,i) = T_0_p(j,i) + dt_3d * tsc(3) * tT_soil_m(0,j,i) ! !-- Calculate true tendency tend = (T_0_p(j,i) - T_0(j,i) - tsc(3) * tT_0_m(j,i)) / (dt_3d & * tsc(2)) ! !-- Calculate T_0 tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN tT_0_m(j,i) = tend ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN tT_0_m(j,i) = -9.5625_wp * tend + 5.3125_wp * tT_0_m(j,i) ENDIF ENDIF pt_p(k,j,i) = T_0_p(j,i) / exn ! !-- Calculate fluxes Rn(j,i) = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i)**4 & - 4.0_wp * sigma_SB * T_0(j,i)**3 * T_0_p(j,i) G(j,i) = lambda_skin * (T_0_p(j,i) - T_soil(0,j,i)) H(j,i) = - f_H * ( pt_p(k+1,j,i) - pt_p(k,j,i) ) IF ( humidity ) THEN LE(j,i) = - f_LE * ( q_p(k+1,j,i) - q_s + dq_s_dT & * T_0(j,i) - dq_s_dT * T_0_p(j,i) ) LE_veg(j,i) = - f_LE_veg * ( q_p(k+1,j,i) - q_s + dq_s_dT & * T_0(j,i) - dq_s_dT * T_0_p(j,i) ) LE_soil(j,i) = - f_LE_soil * ( q_p(k+1,j,i) - q_s + dq_s_dT & * T_0(j,i) - dq_s_dT * T_0_p(j,i) ) LE_liq(j,i) = - f_LE_liq * ( q_p(k+1,j,i) - q_s + dq_s_dT & * T_0(j,i) - dq_s_dT * T_0_p(j,i) ) ENDIF ! IF ( i == 1 .AND. j == 1 ) THEN ! PRINT*, "Rn", Rn(j,i) ! PRINT*, "H", H(j,i) ! PRINT*, "LE", LE(j,i) ! PRINT*, "LE_liq", LE_liq(j,i) ! PRINT*, "LE_veg", LE_veg(j,i) ! PRINT*, "LE_soil", LE_soil(j,i) ! PRINT*, "G", G(j,i) ! ENDIF ! !-- Calculate the true surface resistance IF ( LE(j,i) .EQ. 0.0 ) THEN r_s(j,i) = 1.0E10 ELSE r_s(j,i) = - rho_lv * ( q_p(k+1,j,i) - q_s + dq_s_dT * T_0(j,i)& - dq_s_dT * T_0_p(j,i) ) / LE(j,i) - r_a(j,i) ENDIF ! !-- Calculate fluxes in the atmosphere shf(j,i) = H(j,i) / rho_cp ! !-- Calculate change in liquid water reservoir due to dew fall or !-- evaporation of liquid water IF ( humidity ) THEN ! !-- If precipitation is activated, add rain water to LE_liq. !-- precipitation_rate is given in mm. IF ( precipitation ) THEN LE_liq(j,i) = LE_liq(j,i) + precipitation_rate(j,i) & * 0.001_wp * rho_l * l_v ENDIF ! !-- If the air is saturated, check the reservoir water level IF ( q_s .LE. q_p(k+1,j,i)) THEN ! !-- Check if reservoir is full (avoid values > m_liq_max) !-- In that case, LE_liq goes to LE_soil. In this case !-- LE_veg is zero anyway (because c_liq = 1), so that tend is !-- zero and no further check is needed IF ( m_liq(j,i) .EQ. m_liq_max ) THEN LE_soil(j,i) = LE_soil(j,i) + LE_liq(j,i) LE_liq(j,i) = 0.0_wp ENDIF ! !-- In case LE_veg becomes negative (unphysical behavior), let !-- the water enter the liquid water reservoir as dew on the !-- plant IF ( LE_veg(j,i) .LT. 0.0_wp ) THEN LE_liq(j,i) = LE_liq(j,i) + LE_veg(j,i) LE_veg(j,i) = 0.0_wp ENDIF ENDIF tend = - LE_liq(j,i) * drho_l_lv m_liq_p(j,i) = m_liq(j,i) + dt_3d * ( tsc(2) * tend & + tsc(3) * tm_liq_m(j,i) ) ! !-- Check if reservoir is overfull -> reduce to maximum !-- (conservation of water is violated here) m_liq_p(j,i) = MIN(m_liq_p(j,i),m_liq_max) ! !-- Check if reservoir is empty (avoid values < 0.0) !-- (conservation of water is violated here) m_liq_p(j,i) = MAX(m_liq_p(j,i),0.0_wp) ! !-- Calculate m_liq tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN tm_liq_m(j,i) = tend ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN tm_liq_m(j,i) = -9.5625_wp * tend + 5.3125_wp & * tm_liq_m(j,i) ENDIF ENDIF ! !-- Calculate moisture flux in the atmosphere qsws(j,i) = LE(j,i) / rho_lv ENDIF ENDDO ENDDO END SUBROUTINE lsm_energy_balance !------------------------------------------------------------------------------! ! Description: ! ------------ ! !------------------------------------------------------------------------------! SUBROUTINE lsm_soil_model IMPLICIT NONE INTEGER(iwp) :: i !: running index INTEGER(iwp) :: j !: running index INTEGER(iwp) :: k !: running index REAL(wp) :: h_VG !: Van Genuchten coef. h REAL(wp), DIMENSION(0:soil_layers-1) :: gamma_temp, & !: temp. gamma lambda_temp, & !: temp. lambda tend !: tendency DO i = nxlg, nxrg DO j = nysg, nyng DO k = 0, soil_layers-1 ! !-- Calculate volumetric heat capacity of the soil, taking into !-- account water content rhoC_total(k,j,i) = (rhoC_soil * (1.0 - m_sat(j,i)) & + rhoC_water * m_soil(k,j,i)) ! !-- Calculate soil heat conductivity at the center of the soil !-- layers Ke = 1.0 + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i))) lambda_temp(k) = Ke * (lambda_h_sat(j,i) + lambda_h_dry) + & lambda_h_dry ENDDO ! !-- Calculate soil heat conductivity (lambda_h) at the _stag level !-- using linear interpolation DO k = 0, soil_layers-2 lambda_h(k,j,i) = lambda_temp(k) + & ( lambda_temp(k+1) - lambda_temp(k) ) & * 0.5 * dz_soil_stag(k) * ddz_soil(k+1) ENDDO lambda_h(soil_layers-1,j,i) = lambda_temp(soil_layers-1) ! !-- Prognostic equation for soil temperature T_soil tend(:) = 0.0_wp tend(0) = (1.0/rhoC_total(0,j,i)) * & ( lambda_h(0,j,i) * ( T_soil(1,j,i) - T_soil(0,j,i) ) & * ddz_soil(0) + G(j,i) ) * ddz_soil_stag(0) DO k = 1, soil_layers-1 tend(k) = (1.0/rhoC_total(k,j,i)) & * ( lambda_h(k,j,i) & * ( T_soil(k+1,j,i) - T_soil(k,j,i) ) & * ddz_soil(k) & - lambda_h(k-1,j,i) & * ( T_soil(k,j,i) - T_soil(k-1,j,i) ) & * ddz_soil(k-1) & ) * ddz_soil_stag(k) ENDDO T_soil_p(0:soil_layers-1,j,i) = T_soil(0:soil_layers-1,j,i) & + dt_3d * ( tsc(2) & * tend(:) + tsc(3) & * tT_soil_m(:,j,i) ) ! !-- Calculate T_soil tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO k = 0, soil_layers-1 tT_soil_m(k,j,i) = tend(k) ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO k = 0, soil_layers-1 tT_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp & * tT_soil_m(k,j,i) ENDDO ENDIF ENDIF DO k = 0, soil_layers-1 ! !-- Calculate soil diffusivity at the center of the soil layers lambda_temp(k) = (- b_CH * gamma_w_sat(j,i) * psi_sat & / m_sat(j,i) ) * ( MAX(m_soil(k,j,i), & m_wilt(j,i)) / m_sat(j,i) )**(b_CH + 2.0_wp) ! !-- Calculate the hydraulic conductivity after Van Genuchten (1980) h_VG = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) - & MAX(m_soil(k,j,i),m_wilt(j,i)) ) )**(n_VG(j,i) & / (n_VG(j,i)-1.0_wp)) - 1.0_wp & )**(1.0_wp/n_VG(j,i)) / alpha_VG(j,i) gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp + & (alpha_VG(j,i)*h_VG)**n_VG(j,i))**(1.0_wp & -1.0_wp/n_VG(j,i)) - (alpha_VG(j,i)*h_VG & )**(n_VG(j,i)-1.0_wp))**2 ) & / ( (1.0_wp + (alpha_VG(j,i)*h_VG)**n_VG(j,i) & )**((1.0_wp - 1.0_wp/n_VG(j,i))*(l_VG(j,i) & + 2.0)) ) ENDDO IF ( humidity ) THEN ! !-- Calculate soil diffusivity (lambda_w) at the _stag level !-- using linear interpolation DO k = 0, soil_layers-2 lambda_w(k,j,i) = lambda_temp(k) + & ( lambda_temp(k+1) - lambda_temp(k) ) & * 0.5 * dz_soil_stag(k) * ddz_soil(k+1) gamma_w(k,j,i) = gamma_temp(k) + & ( gamma_temp(k+1) - gamma_temp(k) ) & * 0.5 * dz_soil_stag(k) * ddz_soil(k+1) ENDDO ! ! !-- In case of a closed bottom (= water content is conserved), set !-- hydraulic conductivity to zero to that no water will be lost !-- in the bottom layer. IF ( conserve_water_content ) THEN gamma_w(soil_layers-1,j,i) = 0.0_wp ELSE gamma_w(soil_layers-1,j,i) = lambda_temp(soil_layers-1) ENDIF !-- The root extraction (= root_extr * LE_veg / (rho_l * l_v)) !-- ensures the mass conservation for water. The transpiration of !-- plants equals the cumulative withdrawals by the roots in the !-- soil. The scheme takes into account the availability of water !-- in the soil layers as well as the root fraction in the !-- respective layer ! !-- Calculate the root extraction (ECMWF 7.69, with some !-- modifications) m_total = 0.0_wp DO k = 0, soil_layers-1 m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i) * & dz_soil_stag(k) ENDDO ! !-- For conservation of mass, the sum of root_extr must be 1 DO k = 0, soil_layers-1 root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) & * dz_soil_stag(k) / m_total ENDDO ! !-- Prognostic equation for soil water content m_soil tend(:) = 0.0_wp tend(0) = ( lambda_w(0,j,i) * ( m_soil(1,j,i) - m_soil(0,j,i) )& * ddz_soil(0) - gamma_w(0,j,i) - ( root_extr(0) & * LE_veg(j,i) + LE_soil(j,i) ) * drho_l_lv & ) * ddz_soil_stag(0) DO k = 1, soil_layers-2 tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i) & - m_soil(k,j,i) ) * ddz_soil(k) - gamma_w(k,j,i)& - lambda_w(k-1,j,i) * (m_soil(k,j,i) - & m_soil(k-1,j,i)) * ddz_soil(k-1) & + gamma_w(k-1,j,i) - (root_extr(k) * LE_veg(j,i)& * drho_l_lv) & ) * ddz_soil_stag(k) ENDDO tend(soil_layers-1) = ( - gamma_w(soil_layers-1,j,i) & - lambda_w(soil_layers-2,j,i) & * (m_soil(soil_layers-1,j,i) & - m_soil(soil_layers-2,j,i)) & * ddz_soil(soil_layers-2) & + gamma_w(soil_layers-2,j,i) - ( & root_extr(soil_layers-1) & * LE_veg(j,i) * drho_l_lv ) & ) * ddz_soil_stag(soil_layers-1) m_soil_p(0:soil_layers-1,j,i) = m_soil(0:soil_layers-1,j,i) & + dt_3d * ( tsc(2) * tend(:) & + tsc(3) * tm_soil_m(:,j,i) ) ! !-- Account for dry soils (find a better solution here!) m_soil_p(:,j,i) = MAX(m_soil_p(:,j,i),0.0_wp) ! !-- Calculate m_soil tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO k = 0, soil_layers-1 tm_soil_m(k,j,i) = tend(k) ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO k = 0, soil_layers-1 tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp & * tm_soil_m(k,j,i) ENDDO ENDIF ENDIF ENDIF ENDDO ENDDO ! !-- Calculate surface specific humidity IF ( humidity ) THEN CALL calc_q0 ENDIF END SUBROUTINE lsm_soil_model !------------------------------------------------------------------------------! ! Description: ! ------------ ! !------------------------------------------------------------------------------! SUBROUTINE calc_q0 IMPLICIT NONE INTEGER :: i !: running index INTEGER :: j !: running index INTEGER :: k !: running index REAL(wp) :: resistance !: aerodynamic and soil resistance term DO i = nxlg, nxrg DO j = nysg, nyng k = nzb_s_inner(j,i) ! !-- Temporary solution as long as T_0 is prescribed pt_p(k,j,i) = T_0(j,i) / exn ! !-- Calculate water vapour pressure at saturation e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( T_0(j,i) - & 273.16_wp ) / ( T_0(j,i) - & 35.86_wp ) ) ! !-- Calculate specific humidity at saturation q_s = 0.622_wp * e_s / surface_pressure resistance = r_a(j,i) / (r_a(j,i) + r_s(j,i)) ! !-- Calculate specific humidity at surface q_p(k,j,i) = resistance * q_s + (1.0_wp - resistance) & * q_p(k+1,j,i) ENDDO ENDDO END SUBROUTINE calc_q0 END MODULE land_surface_model_mod