Changeset 1551 for palm/trunk/SOURCE/land_surface_model.f90
- Timestamp:
- Mar 3, 2015 2:18:16 PM (10 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/land_surface_model.f90 (modified) (43 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/land_surface_model.f90
r1514 r1551 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Flux calculation is now done in prandtl_fluxes. Added support for data output. 23 ! Vertical indices have been replaced. Restart runs are now possible. Some 24 ! variables have beem renamed. Bugfix in the prognostic equation for the surface 25 ! temperature. Introduced z0_eb and z0h_eb, which overwrite the setting of 26 ! roughness_length and z0_factor. Added Clapp & Hornberger parametrization for 27 ! the hydraulic conductivity. Bugfix for root fraction and extraction 28 ! calculation 23 29 ! 24 30 ! Former revisions: … … 43 49 ! scheme implemented in the ECMWF IFS model, with modifications according to 44 50 ! H-TESSEL. The implementation is based on the formulation implemented in the 45 ! DALES model.51 ! DALES and UCLA-LES models. 46 52 ! 47 53 ! To do list: 48 54 ! ----------- 49 ! - Add support for binary I/O support 50 ! - Add support for lsm data output 51 ! - Check for time step criterion 52 ! - Check use with RK-2 and Euler time-stepping 53 ! - Adaption for use with cloud physics (liquid water potential temperature) 54 ! - Check reaction of plants at wilting point and at atmospheric saturation 55 ! - Consider partial absorption of the net shortwave radiation by the skin layer 55 ! - Check dewfall parametrization for fog simulations 56 ! - Consider partial absorption of the net shortwave radiation by the surface layer 56 57 ! - Allow for water surfaces, check performance for bare soils 58 ! - Invert indices (running from -3 to 0. Currently: nzb_soil=0, nzt_soil=3)) 59 ! - Implement surface runoff model (required when performing long-term LES with 60 ! considerable precipitation 61 ! Notes: 62 ! ------ 63 ! - No time step criterion is required as long as the soil layers do not become 64 ! too thin 57 65 !------------------------------------------------------------------------------! 58 66 USE arrays_3d, & … … 60 68 61 69 USE cloud_parameters, & 62 ONLY: cp, l_d_r, l_v, precipitation_rate, rho_l, r_d, r_v70 ONLY: cp, l_d_r, l_v, precipitation_rate, rho_l, r_d, r_v 63 71 64 72 USE control_parameters, & 65 ONLY: dt_3d, humidity, intermediate_timestep_count, & 66 intermediate_timestep_count_max, precipitation, pt_surface, & 67 rho_surface, surface_pressure, timestep_scheme, tsc 73 ONLY: cloud_physics, dt_3d, humidity, intermediate_timestep_count, & 74 initializing_actions, intermediate_timestep_count_max, & 75 max_masks, precipitation, pt_surface, rho_surface, & 76 roughness_length, surface_pressure, timestep_scheme, tsc, & 77 z0h_factor 68 78 69 79 USE indices, & 70 ONLY: nxlg, nxrg, nyng, nysg, nzb_s_inner80 ONLY: nxlg, nxrg, nyng, nysg, nzb_s_inner 71 81 72 82 USE kinds 73 83 84 USE netcdf_control, & 85 ONLY: dots_label, dots_num, dots_unit 86 74 87 USE radiation_model_mod, & 75 ONLY: Rn, SW_in, sigma_SB 76 88 ONLY: irad_scheme, rad_net, rad_sw_in, sigma_SB 89 90 USE statistics, & 91 ONLY: hom, statistic_regions 77 92 78 93 IMPLICIT NONE … … 80 95 ! 81 96 !-- LSM model constants 82 INTEGER(iwp), PARAMETER :: soil_layers = 4 !: number of soil layers (fixed for now) 97 INTEGER(iwp), PARAMETER :: nzb_soil = 0, & !: bottom of the soil model (to be switched) 98 nzt_soil = 3, & !: top of the soil model (to be switched) 99 nzs = 4 !: number of soil layers (fixed for now) 100 101 INTEGER(iwp) :: dots_soil = 0 !: starting index for timeseries output 102 103 INTEGER(iwp), DIMENSION(0:1) :: id_dim_zs_xy, id_dim_zs_xz, id_dim_zs_yz, & 104 id_dim_zs_3d, id_var_zs_xy, & 105 id_var_zs_xz, id_var_zs_yz, id_var_zs_3d 106 107 INTEGER(iwp), DIMENSION(1:max_masks,0:1) :: id_dim_zs_mask, id_var_zs_mask 83 108 84 109 REAL(wp), PARAMETER :: & 85 b_ CH= 6.04_wp, & ! Clapp & Hornberger exponent86 lambda_h_dry = 0.19_wp, & ! heat conductivity for dry soil 110 b_ch = 6.04_wp, & ! Clapp & Hornberger exponent 111 lambda_h_dry = 0.19_wp, & ! heat conductivity for dry soil 87 112 lambda_h_sm = 3.44_wp, & ! heat conductivity of the soil matrix 88 113 lambda_h_water = 0.57_wp, & ! heat conductivity of water 89 114 psi_sat = -0.388_wp, & ! soil matrix potential at saturation 90 rho C_soil= 2.19E6_wp, & ! volumetric heat capacity of soil91 rho C_water= 4.20E6_wp, & ! volumetric heat capacity of water115 rho_c_soil = 2.19E6_wp, & ! volumetric heat capacity of soil 116 rho_c_water = 4.20E6_wp, & ! volumetric heat capacity of water 92 117 m_max_depth = 0.0002_wp ! Maximum capacity of the water reservoir (m) 93 118 … … 99 124 100 125 LOGICAL :: conserve_water_content = .TRUE., & !: open or closed bottom surface for the soil model 126 dewfall = .TRUE., & !: allow/inhibit dewfall 101 127 land_surface = .FALSE. !: flag parameter indicating wheather the lsm is used 102 128 103 129 ! value 9999999.9_wp -> generic available or user-defined value must be set 104 130 ! otherwise -> no generic variable and user setting is optional 105 REAL(wp) :: alpha_ VanGenuchten = 0.0_wp, & !: NAMELIST alpha_VG106 canopy_resistance_coefficient = 0.0_wp, & !: NAMELIST gD107 C_skin = 20000.0_wp, & !: Skinheat capacity131 REAL(wp) :: alpha_vangenuchten = 9999999.9_wp, & !: NAMELIST alpha_vg 132 canopy_resistance_coefficient = 9999999.9_wp, & !: NAMELIST g_d 133 c_surface = 20000.0_wp, & !: Surface (skin) heat capacity 108 134 drho_l_lv, & !: (rho_l * l_v)**-1 109 135 exn, & !: value of the Exner function 110 136 e_s = 0.0_wp, & !: saturation water vapour pressure 111 field_capacity = 0.0_wp, & !: NAMELIST m_fc 112 f_shortwave_incoming = 9999999.9_wp, & !: NAMELIST f_SW_in 113 hydraulic_conductivity = 0.0_wp, & !: NAMELIST gamma_w_sat 114 Ke = 0.0_wp, & !: Kersten number 115 lambda_skin_stable = 9999999.9_wp, & !: NAMELIST lambda_skin_s 116 lambda_skin_unstable = 9999999.9_wp, & !: NAMELIST lambda_skin_u 117 leaf_area_index = 9999999.9_wp, & !: NAMELIST LAI 118 l_VanGenuchten = 0.0_WP, & !: NAMELIST l_VG 119 min_canopy_resistance = 110.0_wp, & !: NAMELIST r_s_min 120 m_total = 0.0_wp, & !: weighed total water content of the soil (m3/m3) 121 n_VanGenuchten = 0.0_WP, & !: NAMELIST n_VG 137 field_capacity = 9999999.9_wp, & !: NAMELIST m_fc 138 f_shortwave_incoming = 9999999.9_wp, & !: NAMELIST f_sw_in 139 hydraulic_conductivity = 9999999.9_wp, & !: NAMELIST gamma_w_sat 140 ke = 0.0_wp, & !: Kersten number 141 lambda_surface_stable = 9999999.9_wp, & !: NAMELIST lambda_surface_s 142 lambda_surface_unstable = 9999999.9_wp, & !: NAMELIST lambda_surface_u 143 leaf_area_index = 9999999.9_wp, & !: NAMELIST lai 144 l_vangenuchten = 9999999.9_wp, & !: NAMELIST l_vg 145 min_canopy_resistance = 9999999.9_wp, & !: NAMELIST r_canopy_min 146 min_soil_resistance = 50.0_wp, & !: NAMELIST r_soil_min 147 m_total = 0.0_wp, & !: weighted total water content of the soil (m3/m3) 148 n_vangenuchten = 9999999.9_wp, & !: NAMELIST n_vg 122 149 q_s = 0.0_wp, & !: saturation specific humidity 123 residual_moisture = 0.0_wp,& !: NAMELIST m_res150 residual_moisture = 9999999.9_wp, & !: NAMELIST m_res 124 151 rho_cp, & !: rho_surface * cp 125 152 rho_lv, & !: rho * l_v 126 153 rd_d_rv, & !: r_d / r_v 127 saturation_moisture = 0.0_wp,& !: NAMELIST m_sat154 saturation_moisture = 9999999.9_wp, & !: NAMELIST m_sat 128 155 vegetation_coverage = 9999999.9_wp, & !: NAMELIST c_veg 129 wilting_point = 0.0_wp !: NAMELIST m_wilt 130 131 REAL(wp), DIMENSION(0:soil_layers-1) :: & 156 wilting_point = 9999999.9_wp, & !: NAMELIST m_wilt 157 z0_eb = 9999999.9_wp, & !: NAMELIST z0 (lsm_par) 158 z0h_eb = 9999999.9_wp !: NAMELIST z0h (lsm_par) 159 160 REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: & 132 161 ddz_soil, & !: 1/dz_soil 133 162 ddz_soil_stag, & !: 1/dz_soil_stag … … 135 164 dz_soil_stag, & !: soil grid spacing (edge-edge) 136 165 root_extr = 0.0_wp, & !: root extraction 137 root_fraction = (/0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp/), & !: distribution of root surface area to the individual soil layers 138 soil_level = (/0.07_wp, 0.28_wp, 1.00_wp, 2.89_wp/), & !: soil layer depths (m) 166 root_fraction = (/9999999.9_wp, 9999999.9_wp, & 167 9999999.9_wp, 9999999.9_wp /), & !: distribution of root surface area to the individual soil layers 168 zs = (/0.07_wp, 0.28_wp, 1.00_wp, 2.89_wp/), & !: soil layer depths (m) 139 169 soil_moisture = 0.0_wp !: soil moisture content (m3/m3) 140 170 141 REAL(wp), DIMENSION( 0:soil_layers) :: &171 REAL(wp), DIMENSION(nzb_soil:nzt_soil+1) :: & 142 172 soil_temperature = 9999999.9_wp !: soil temperature (K) 143 173 144 174 #if defined( __nopointer ) 145 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: T_0, & !: skin temperature (K) 146 T_0_p, & !: progn. skin temperature (K) 147 m_liq, & !: liquid water reservoir (m) 148 m_liq_p !: progn. liquid water reservoir (m) 175 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface, & !: surface temperature (K) 176 t_surface_p, & !: progn. surface temperature (K) 177 m_liq_eb, & !: liquid water reservoir (m) 178 m_liq_eb_av, & !: liquid water reservoir (m) 179 m_liq_eb_p !: progn. liquid water reservoir (m) 149 180 #else 150 REAL(wp), DIMENSION(:,:), POINTER :: T_0, & 151 T_0_p, & 152 m_liq, & 153 m_liq_p 154 155 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: T_0_1, T_0_2, & 156 m_liq_1, m_liq_2 181 REAL(wp), DIMENSION(:,:), POINTER :: t_surface, & 182 t_surface_p, & 183 m_liq_eb, & 184 m_liq_eb_p 185 186 REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface_1, t_surface_2, & 187 m_liq_eb_av, & 188 m_liq_eb_1, m_liq_eb_2 157 189 #endif 158 190 159 191 ! 160 192 !-- Temporal tendencies for time stepping 161 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t T_0_m, & !: skintemperature tendency (K)162 tm_liq_ m!: liquid water reservoir tendency (m)193 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_surface_m, & !: surface temperature tendency (K) 194 tm_liq_eb_m !: liquid water reservoir tendency (m) 163 195 164 196 ! 165 197 !-- Energy balance variables 166 198 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & 167 alpha_VG, & !: coef. of Van Genuchten 168 c_liq, & !: liquid water coverage (of vegetated area) 169 c_veg, & !: vegetation coverage 170 f_SW_in, & !: ? 171 G, & !: surface soil heat flux 172 H, & !: surface flux of sensible heat 173 gamma_w_sat, & !: hydraulic conductivity at saturation 174 gD, & !: coefficient for dependence of r_canopy on water vapour pressure deficit 175 LAI, & !: leaf area index 176 LE, & !: surface flux of latent heat (total) 177 LE_veg, & !: surface flux of latent heat (vegetation portion) 178 LE_soil, & !: surface flux of latent heat (soil portion) 179 LE_liq, & !: surface flux of latent heat (liquid water portion) 180 lambda_h_sat, & !: heat conductivity for dry soil 181 lambda_skin_s, & !: coupling between skin and soil (depends on vegetation type) 182 lambda_skin_u, & !: coupling between skin and soil (depends on vegetation type) 183 l_VG, & !: coef. of Van Genuchten 184 m_fc, & !: soil moisture at field capacity (m3/m3) 185 m_res, & !: residual soil moisture 186 m_sat, & !: saturation soil moisture (m3/m3) 187 m_wilt, & !: soil moisture at permanent wilting point (m3/m3) 188 n_VG, & !: coef. Van Genuchten 189 r_a, & !: aerodynamic resistance 190 r_canopy, & !: canopy resistance 191 r_soil, & !: soil resitance 192 r_soil_min, & !: minimum soil resistance 193 r_s, & !: total surface resistance (combination of r_soil and r_canopy) 194 r_s_min !: minimum canopy resistance 199 alpha_vg, & !: coef. of Van Genuchten 200 c_liq, & !: liquid water coverage (of vegetated area) 201 c_liq_av, & !: average of c_liq 202 c_soil_av, & !: average of c_soil 203 c_veg, & !: vegetation coverage 204 c_veg_av, & !: average of c_veg 205 f_sw_in, & !: fraction of absorbed shortwave radiation by the surface layer (not implemented yet) 206 ghf_eb, & !: ground heat flux 207 ghf_eb_av, & !: average of ghf_eb 208 gamma_w_sat, & !: hydraulic conductivity at saturation 209 g_d, & !: coefficient for dependence of r_canopy on water vapour pressure deficit 210 lai, & !: leaf area index 211 lai_av, & !: average of lai 212 lambda_h_sat, & !: heat conductivity for dry soil 213 lambda_surface_s, & !: coupling between surface and soil (depends on vegetation type) 214 lambda_surface_u, & !: coupling between surface and soil (depends on vegetation type) 215 l_vg, & !: coef. of Van Genuchten 216 m_fc, & !: soil moisture at field capacity (m3/m3) 217 m_res, & !: residual soil moisture 218 m_sat, & !: saturation soil moisture (m3/m3) 219 m_wilt, & !: soil moisture at permanent wilting point (m3/m3) 220 n_vg, & !: coef. Van Genuchten 221 qsws_eb, & !: surface flux of latent heat (total) 222 qsws_eb_av, & !: average of qsws_eb 223 qsws_liq_eb, & !: surface flux of latent heat (liquid water portion) 224 qsws_liq_eb_av, & !: average of qsws_liq_eb 225 qsws_soil_eb, & !: surface flux of latent heat (soil portion) 226 qsws_soil_eb_av, & !: average of qsws_soil_eb 227 qsws_veg_eb, & !: surface flux of latent heat (vegetation portion) 228 qsws_veg_eb_av, & !: average of qsws_veg_eb 229 r_a, & !: aerodynamic resistance 230 r_canopy, & !: canopy resistance 231 r_soil, & !: soil resitance 232 r_soil_min, & !: minimum soil resistance 233 r_s, & !: total surface resistance (combination of r_soil and r_canopy) 234 r_canopy_min, & !: minimum canopy (stomatal) resistance 235 shf_eb, & !: surface flux of sensible heat 236 shf_eb_av !: average of shf_eb 195 237 196 238 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & … … 198 240 lambda_w, & !: hydraulic diffusivity of soil (?) 199 241 gamma_w, & !: hydraulic conductivity of soil (?) 200 rho C_total!: volumetric heat capacity of the actual soil matrix (?)242 rho_c_total !: volumetric heat capacity of the actual soil matrix (?) 201 243 202 244 #if defined( __nopointer ) 203 245 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: & 204 T_soil, & !: Soil temperature (K) 205 T_soil_p, & !: Prog. soil temperature (K) 246 t_soil, & !: Soil temperature (K) 247 t_soil_av, & !: Average of t_soil 248 t_soil_p, & !: Prog. soil temperature (K) 206 249 m_soil, & !: Soil moisture (m3/m3) 250 m_soil_av, & !: Average of m_soil 207 251 m_soil_p !: Prog. soil moisture (m3/m3) 208 252 #else 209 253 REAL(wp), DIMENSION(:,:,:), POINTER :: & 210 T_soil, T_soil_p, &254 t_soil, t_soil_p, & 211 255 m_soil, m_soil_p 212 256 213 257 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: & 214 T_soil_1, T_soil_2,&215 m_soil_ 1, m_soil_2258 t_soil_av, t_soil_1, t_soil_2, & 259 m_soil_av, m_soil_1, m_soil_2 216 260 217 261 … … 220 264 221 265 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 222 t T_soil_m, & !: T_soil storage array266 tt_soil_m, & !: t_soil storage array 223 267 tm_soil_m, & !: m_soil storage array 224 268 root_fr !: root fraction (sum=1) 225 269 226 ! 227 !-- Land surface parameters according to the following classes (veg_type) 228 !-- (0 user defined) 229 !-- 1 crops, mixed farming 230 !-- 2 short grass 231 !-- 3 evergreen needleleaf trees 232 !-- 4 deciduous needleleaf trees 233 !-- 5 evergreen broadleaf trees 234 !-- 6 deciduous broadleaf trees 235 !-- 7 tall grass 236 !-- 8 desert 237 !-- 9 tundra 238 !-- 10 irrigated crops 239 !-- 11 semidesert 240 !-- 12 ice caps and glaciers 241 !-- 13 bogs and marshes 242 !-- 14 inland water 243 !-- 15 ocean 244 !-- 16 evergreen shrubs 245 !-- 17 deciduous shrubs 246 !-- 18 mixed forest/woodland 247 !-- 19 interrupted forest 248 249 ! 250 !-- Land surface parameters I r_s_min, LAI, c_veg, gD 270 271 ! 272 !-- Predefined Land surface classes (veg_type) 273 CHARACTER(26), DIMENSION(0:19) :: veg_type_name = (/ & 274 'user defined', & ! 0 275 'crops, mixed farming', & ! 1 276 'short grass', & ! 2 277 'evergreen needleleaf trees', & ! 3 278 'deciduous needleleaf trees', & ! 4 279 'evergreen broadleaf trees' , & ! 5 280 'deciduous broadleaf trees', & ! 6 281 'tall grass', & ! 7 282 'desert', & ! 8 283 'tundra', & ! 9 284 'irrigated crops', & ! 10 285 'semidesert', & ! 11 286 'ice caps and glaciers' , & ! 12 287 'bogs and marshes', & ! 13 288 'inland water', & ! 14 289 'ocean', & ! 15 290 'evergreen shrubs', & ! 16 291 'deciduous shrubs', & ! 17 292 'mixed forest/woodland', & ! 18 293 'interrupted forest' & ! 19 294 /) 295 296 ! 297 !-- Soil model classes (soil_type) 298 CHARACTER(12), DIMENSION(0:7) :: soil_type_name = (/ & 299 'user defined', & ! 0 300 'coarse', & ! 1 301 'medium', & ! 2 302 'medium-fine', & ! 3 303 'fine', & ! 4 304 'very fine' , & ! 5 305 'organic', & ! 6 306 'loamy (CH)' & ! 7 307 /) 308 ! 309 !-- Land surface parameters according to the respective classes (veg_type) 310 311 ! 312 !-- Land surface parameters I 313 !-- r_canopy_min, lai, c_veg, g_d 251 314 REAL(wp), DIMENSION(0:3,1:19) :: veg_pars = RESHAPE( (/ & 252 315 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & ! 1 … … 257 320 240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp, & ! 6 258 321 100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp, & ! 7 259 250.0_wp, 0. 50_wp, 0.00_wp, 0.00_wp, & ! 8322 250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp, & ! 8 260 323 80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp, & ! 9 261 324 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & ! 10 … … 296 359 297 360 ! 298 !-- Land surface parameters III lambda_s kin_s, lambda_skin_u, f_SW_in299 REAL(wp), DIMENSION(0:2,1:19) :: s kin_pars = RESHAPE( (/ &361 !-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in 362 REAL(wp), DIMENSION(0:2,1:19) :: surface_pars = RESHAPE( (/ & 300 363 10.0_wp, 10.0_wp, 0.05_wp, & ! 1 301 364 10.0_wp, 10.0_wp, 0.05_wp, & ! 2 … … 345 408 ! 346 409 !-- Soil parameters according to the following porosity classes (soil_type) 347 !-- (0 user defined) 348 !-- 1 coarse 349 !-- 2 medium 350 !-- 3 medium-fine 351 !-- 4 fine 352 !-- 5 very fine 353 !-- 6 organic 354 ! 355 !-- Soil parameters I alpha_VG, l_VG, n_VG, gamma_w_sat 356 REAL(wp), DIMENSION(0:3,1:6) :: soil_pars = RESHAPE( (/ & 410 411 ! 412 !-- Soil parameters I alpha_vg, l_vg, n_vg, gamma_w_sat 413 REAL(wp), DIMENSION(0:3,1:7) :: soil_pars = RESHAPE( (/ & 357 414 3.83_wp, 1.250_wp, 1.38_wp, 6.94E-6_wp, & ! 1 358 415 3.14_wp, -2.342_wp, 1.28_wp, 1.16E-6_wp, & ! 2 … … 360 417 3.67_wp, -1.977_wp, 1.10_wp, 2.87E-6_wp, & ! 4 361 418 2.65_wp, 2.500_wp, 1.10_wp, 1.74E-6_wp, & ! 5 362 1.30_wp, 0.400_wp, 1.20_wp, 0.93E-6_wp & ! 6 363 /), (/ 4, 6 /) ) 419 1.30_wp, 0.400_wp, 1.20_wp, 0.93E-6_wp, & ! 6 420 0.00_wp, 0.00_wp, 0.00_wp, 0.57E-6_wp & ! 7 421 /), (/ 4, 7 /) ) 364 422 365 423 ! 366 424 !-- Soil parameters II m_sat, m_fc, m_wilt, m_res 367 REAL(wp), DIMENSION(0:3,1: 6) :: m_soil_pars = RESHAPE( (/ &425 REAL(wp), DIMENSION(0:3,1:7) :: m_soil_pars = RESHAPE( (/ & 368 426 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & ! 1 369 427 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & ! 2 … … 371 429 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, & ! 4 372 430 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, & ! 5 373 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp & ! 6 374 /), (/ 4, 6 /) ) 431 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp, & ! 6 432 0.472_wp, 0.323_wp, 0.171_wp, 0.000_wp & ! 7 433 /), (/ 4, 7 /) ) 375 434 376 435 … … 381 440 382 441 383 PUBLIC alpha_VanGenuchten, C_skin, canopy_resistance_coefficient, & 384 conserve_water_content, field_capacity, f_shortwave_incoming, & 385 hydraulic_conductivity, init_lsm, lambda_skin_stable, & 386 lambda_skin_unstable, land_surface, leaf_area_index, & 387 lsm_energy_balance, lsm_soil_model, l_VanGenuchten, & 388 min_canopy_resistance, n_VanGenuchten, residual_moisture, & 389 root_fraction, saturation_moisture, soil_level, soil_moisture, & 390 soil_temperature, soil_type, vegetation_coverage, veg_type, & 391 wilting_point 392 393 #if defined( __nopointer ) 394 PUBLIC m_liq, m_liq_p, m_soil, m_soil_p, T_0, T_0_p, T_soil, T_soil_p 395 #else 396 PUBLIC m_liq, m_liq_1, m_liq_2, m_liq_p, m_soil, m_soil_1, m_soil_2, & 397 m_soil_p, T_0, T_0_1, T_0_2, T_0_p, T_soil, T_soil_1, T_soil_2, & 398 T_soil_p 399 #endif 442 ! 443 !-- Public parameters, constants and initial values 444 PUBLIC alpha_vangenuchten, c_surface, canopy_resistance_coefficient, & 445 conserve_water_content, dewfall, field_capacity, & 446 f_shortwave_incoming, hydraulic_conductivity, init_lsm, & 447 init_lsm_arrays, lambda_surface_stable, lambda_surface_unstable, & 448 land_surface, leaf_area_index, lsm_energy_balance, lsm_soil_model, & 449 lsm_swap_timelevel, l_vangenuchten, min_canopy_resistance, & 450 min_soil_resistance, n_vangenuchten, residual_moisture, rho_cp, & 451 rho_lv, root_fraction, saturation_moisture, soil_moisture, & 452 soil_temperature, soil_type, soil_type_name, vegetation_coverage, & 453 veg_type, veg_type_name, wilting_point, z0_eb, z0h_eb 454 455 ! 456 !-- Public grid and NetCDF variables 457 PUBLIC dots_soil, id_dim_zs_xy, id_dim_zs_xz, id_dim_zs_yz, & 458 id_dim_zs_3d, id_dim_zs_mask, id_var_zs_xy, id_var_zs_xz, & 459 id_var_zs_yz, id_var_zs_3d, id_var_zs_mask, nzb_soil, nzs, nzt_soil,& 460 zs 461 462 463 ! 464 !-- Public 2D output variables 465 PUBLIC c_liq, c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb, ghf_eb_av, & 466 lai, lai_av, qsws_eb, qsws_eb_av, qsws_liq_eb, qsws_liq_eb_av, & 467 qsws_soil_eb, qsws_soil_eb_av, qsws_veg_eb, qsws_veg_eb_av, & 468 shf_eb, shf_eb_av 469 470 471 ! 472 !-- Public prognostic variables 473 PUBLIC m_liq_eb, m_liq_eb_av, m_soil, m_soil_av, t_soil, t_soil_av 400 474 401 475 … … 412 486 END INTERFACE lsm_soil_model 413 487 488 INTERFACE lsm_swap_timelevel 489 MODULE PROCEDURE lsm_swap_timelevel 490 END INTERFACE lsm_swap_timelevel 414 491 415 492 CONTAINS 416 493 494 495 !------------------------------------------------------------------------------! 496 ! Description: 497 ! ------------ 498 !-- Allocate land surface model arrays and define pointers 499 !------------------------------------------------------------------------------! 500 SUBROUTINE init_lsm_arrays 501 502 503 IMPLICIT NONE 504 505 ! 506 !-- Allocate surface and soil temperature / humidity 507 #if defined( __nopointer ) 508 ALLOCATE ( m_liq_eb(nysg:nyng,nxlg:nxrg) ) 509 ALLOCATE ( m_liq_eb_p(nysg:nyng,nxlg:nxrg) ) 510 ALLOCATE ( m_soil(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 511 ALLOCATE ( m_soil_p(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 512 ALLOCATE ( t_surface(nysg:nyng,nxlg:nxrg) ) 513 ALLOCATE ( t_surface_p(nysg:nyng,nxlg:nxrg) ) 514 ALLOCATE ( t_soil(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) ) 515 ALLOCATE ( t_soil_p(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) ) 516 #else 517 ALLOCATE ( m_liq_eb_1(nysg:nyng,nxlg:nxrg) ) 518 ALLOCATE ( m_liq_eb_2(nysg:nyng,nxlg:nxrg) ) 519 ALLOCATE ( m_soil_1(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 520 ALLOCATE ( m_soil_2(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 521 ALLOCATE ( t_surface_1(nysg:nyng,nxlg:nxrg) ) 522 ALLOCATE ( t_surface_2(nysg:nyng,nxlg:nxrg) ) 523 ALLOCATE ( t_soil_1(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) ) 524 ALLOCATE ( t_soil_2(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) ) 525 #endif 526 527 ! 528 !-- Allocate intermediate timestep arrays 529 ALLOCATE ( tm_liq_eb_m(nysg:nyng,nxlg:nxrg) ) 530 ALLOCATE ( tm_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 531 ALLOCATE ( tt_surface_m(nysg:nyng,nxlg:nxrg) ) 532 ALLOCATE ( tt_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 533 534 ! 535 !-- Allocate 2D vegetation model arrays 536 ALLOCATE ( alpha_vg(nysg:nyng,nxlg:nxrg) ) 537 ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) ) 538 ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) ) 539 ALLOCATE ( f_sw_in(nysg:nyng,nxlg:nxrg) ) 540 ALLOCATE ( ghf_eb(nysg:nyng,nxlg:nxrg) ) 541 ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) ) 542 ALLOCATE ( g_d(nysg:nyng,nxlg:nxrg) ) 543 ALLOCATE ( lai(nysg:nyng,nxlg:nxrg) ) 544 ALLOCATE ( l_vg(nysg:nyng,nxlg:nxrg) ) 545 ALLOCATE ( lambda_h_sat(nysg:nyng,nxlg:nxrg) ) 546 ALLOCATE ( lambda_surface_u(nysg:nyng,nxlg:nxrg) ) 547 ALLOCATE ( lambda_surface_s(nysg:nyng,nxlg:nxrg) ) 548 ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) ) 549 ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) ) 550 ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) ) 551 ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) ) 552 ALLOCATE ( n_vg(nysg:nyng,nxlg:nxrg) ) 553 ALLOCATE ( qsws_eb(nysg:nyng,nxlg:nxrg) ) 554 ALLOCATE ( qsws_soil_eb(nysg:nyng,nxlg:nxrg) ) 555 ALLOCATE ( qsws_liq_eb(nysg:nyng,nxlg:nxrg) ) 556 ALLOCATE ( qsws_veg_eb(nysg:nyng,nxlg:nxrg) ) 557 ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) ) 558 ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) ) 559 ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) ) 560 ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) ) 561 ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) ) 562 ALLOCATE ( r_canopy_min(nysg:nyng,nxlg:nxrg) ) 563 ALLOCATE ( shf_eb(nysg:nyng,nxlg:nxrg) ) 564 565 #if ! defined( __nopointer ) 566 ! 567 !-- Initial assignment of the pointers 568 t_soil => t_soil_1; t_soil_p => t_soil_2 569 t_surface => t_surface_1; t_surface_p => t_surface_2 570 m_soil => m_soil_1; m_soil_p => m_soil_2 571 m_liq_eb => m_liq_eb_1; m_liq_eb_p => m_liq_eb_2 572 #endif 573 574 575 END SUBROUTINE init_lsm_arrays 417 576 418 577 !------------------------------------------------------------------------------! … … 432 591 433 592 ! 593 !-- Calculate Exner function 594 exn = ( surface_pressure / 1000.0_wp )**0.286_wp 595 596 597 ! 598 !-- If no cloud physics is used, rho_surface has not been calculated before 599 IF ( .NOT. cloud_physics ) THEN 600 rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn ) 601 ENDIF 602 603 ! 434 604 !-- Calculate frequently used parameters 435 605 rho_cp = cp * rho_surface 436 606 rd_d_rv = r_d / r_v 437 607 rho_lv = rho_surface * l_v 438 drho_l_lv = 1.0 / (rho_l * l_v) 439 440 ! 441 !-- Allocate skin and soil temperature / humidity 442 #if defined( __nopointer ) 443 ALLOCATE ( T_0(nysg:nyng,nxlg:nxrg) ) 444 ALLOCATE ( T_0_p(nysg:nyng,nxlg:nxrg) ) 445 #else 446 ALLOCATE ( T_0_1(nysg:nyng,nxlg:nxrg) ) 447 ALLOCATE ( T_0_2(nysg:nyng,nxlg:nxrg) ) 448 #endif 449 450 ALLOCATE ( tT_0_m(nysg:nyng,nxlg:nxrg) ) 451 452 #if defined( __nopointer ) 453 ALLOCATE ( T_soil(0:soil_layers,nysg:nyng,nxlg:nxrg) ) 454 ALLOCATE ( T_soil_p(0:soil_layers,nysg:nyng,nxlg:nxrg) ) 455 #else 456 ALLOCATE ( T_soil_1(0:soil_layers,nysg:nyng,nxlg:nxrg) ) 457 ALLOCATE ( T_soil_2(0:soil_layers,nysg:nyng,nxlg:nxrg) ) 458 #endif 459 460 ALLOCATE ( tT_soil_m(0:soil_layers-1,nysg:nyng,nxlg:nxrg) ) 461 462 #if defined( __nopointer ) 463 ALLOCATE ( m_liq(nysg:nyng,nxlg:nxrg) ) 464 ALLOCATE ( m_liq_p(nysg:nyng,nxlg:nxrg) ) 465 #else 466 ALLOCATE ( m_liq_1(nysg:nyng,nxlg:nxrg) ) 467 ALLOCATE ( m_liq_2(nysg:nyng,nxlg:nxrg) ) 468 #endif 469 470 ALLOCATE ( tm_liq_m(nysg:nyng,nxlg:nxrg) ) 471 472 #if defined( __nopointer ) 473 ALLOCATE ( m_soil(0:soil_layers-1,nysg:nyng,nxlg:nxrg) ) 474 ALLOCATE ( m_soil_p(0:soil_layers-1,nysg:nyng,nxlg:nxrg) ) 475 #else 476 ALLOCATE ( m_soil_1(0:soil_layers-1,nysg:nyng,nxlg:nxrg) ) 477 ALLOCATE ( m_soil_2(0:soil_layers-1,nysg:nyng,nxlg:nxrg) ) 478 #endif 479 480 ALLOCATE ( tm_soil_m(0:soil_layers-1,nysg:nyng,nxlg:nxrg) ) 481 482 483 #if ! defined( __nopointer ) 484 ! 485 !-- Initial assignment of the pointers 486 T_soil => T_soil_1; T_soil_p => T_soil_2 487 T_0 => T_0_1; T_0_p => T_0_2 488 m_soil => m_soil_1; m_soil_p => m_soil_2 489 m_liq => m_liq_1; m_liq_p => m_liq_2 490 #endif 491 492 T_0 = 0.0_wp 493 T_0_p = 0.0_wp 494 tT_0_m = 0.0_wp 495 496 T_soil = 0.0_wp 497 T_soil_p = 0.0_wp 498 tT_soil_m = 0.0_wp 499 500 m_liq = 0.0_wp 501 m_liq_p = 0.0_wp 502 tm_liq_m = 0.0_wp 503 504 m_soil = 0.0_wp 505 m_soil_p = 0.0_wp 506 tm_soil_m = 0.0_wp 507 508 ! 509 !-- Allocate 2D vegetation model arrays 510 ALLOCATE ( alpha_VG(nysg:nyng,nxlg:nxrg) ) 511 ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) ) 512 ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) ) 513 ALLOCATE ( f_SW_in(nysg:nyng,nxlg:nxrg) ) 514 ALLOCATE ( G(nysg:nyng,nxlg:nxrg) ) 515 ALLOCATE ( H(nysg:nyng,nxlg:nxrg) ) 516 ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) ) 517 ALLOCATE ( gD(nysg:nyng,nxlg:nxrg) ) 518 ALLOCATE ( LAI(nysg:nyng,nxlg:nxrg) ) 519 ALLOCATE ( LE(nysg:nyng,nxlg:nxrg) ) 520 ALLOCATE ( LE_veg(nysg:nyng,nxlg:nxrg) ) 521 ALLOCATE ( LE_soil(nysg:nyng,nxlg:nxrg) ) 522 ALLOCATE ( LE_liq(nysg:nyng,nxlg:nxrg) ) 523 ALLOCATE ( l_VG(nysg:nyng,nxlg:nxrg) ) 524 ALLOCATE ( lambda_h_sat(nysg:nyng,nxlg:nxrg) ) 525 ALLOCATE ( lambda_skin_u(nysg:nyng,nxlg:nxrg) ) 526 ALLOCATE ( lambda_skin_s(nysg:nyng,nxlg:nxrg) ) 527 ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) ) 528 ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) ) 529 ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) ) 530 ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) ) 531 ALLOCATE ( n_VG(nysg:nyng,nxlg:nxrg) ) 532 ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) ) 533 ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) ) 534 ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) ) 535 ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) ) 536 ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) ) 537 ALLOCATE ( r_s_min(nysg:nyng,nxlg:nxrg) ) 538 539 ! 540 !-- Set initial and default values 541 c_liq = 0.0_wp 542 c_veg = 0.0_wp 543 f_SW_in = 0.05_wp 544 gD = 0.0_wp 545 LAI = 0.0_wp 546 lambda_skin_u = 10.0_wp 547 lambda_skin_s = 10.0_wp 548 549 550 G = 0.0_wp 551 H = rho_cp * shf 608 drho_l_lv = 1.0_wp / (rho_l * l_v) 609 610 ! 611 !-- Set inital values for prognostic quantities 612 tt_surface_m = 0.0_wp 613 tt_soil_m = 0.0_wp 614 tm_liq_eb_m = 0.0_wp 615 c_liq = 0.0_wp 616 617 ghf_eb = 0.0_wp 618 shf_eb = rho_cp * shf 552 619 553 620 IF ( humidity ) THEN 554 LE= rho_l * l_v * qsws621 qsws_eb = rho_l * l_v * qsws 555 622 ELSE 556 LE= 0.0_wp623 qsws_eb = 0.0_wp 557 624 ENDIF 558 625 559 LE_veg= 0.0_wp560 LE_soil = LE561 LE_liq= 0.0_wp626 qsws_liq_eb = 0.0_wp 627 qsws_soil_eb = qsws_eb 628 qsws_veg_eb = 0.0_wp 562 629 563 630 r_a = 50.0_wp 631 r_s = 50.0_wp 564 632 r_canopy = 0.0_wp 565 633 r_soil = 0.0_wp 566 r_soil_min = 50.0_wp567 r_s = 110.0_wp568 r_s_min = min_canopy_resistance569 634 570 635 ! 571 636 !-- Allocate 3D soil model arrays 572 ALLOCATE ( root_fr( 0:soil_layers-1,nysg:nyng,nxlg:nxrg) )573 ALLOCATE ( lambda_h( 0:soil_layers-1,nysg:nyng,nxlg:nxrg) )574 ALLOCATE ( rho C_total(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )637 ALLOCATE ( root_fr(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 638 ALLOCATE ( lambda_h(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 639 ALLOCATE ( rho_c_total(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 575 640 576 641 lambda_h = 0.0_wp … … 578 643 !-- If required, allocate humidity-related variables for the soil model 579 644 IF ( humidity ) THEN 580 ALLOCATE ( lambda_w( 0:soil_layers-1,nysg:nyng,nxlg:nxrg) )581 ALLOCATE ( gamma_w( 0:soil_layers-1,nysg:nyng,nxlg:nxrg) )645 ALLOCATE ( lambda_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 646 ALLOCATE ( gamma_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ) 582 647 583 648 lambda_w = 0.0_wp … … 588 653 !-- the center of the soil layers, whereas gradients/fluxes are defined 589 654 !-- at the edges (_stag) 590 dz_soil_stag( 0) = soil_level(0)591 592 DO k = 1, soil_layers-1593 dz_soil_stag(k) = soil_level(k) - soil_level(k-1)655 dz_soil_stag(nzb_soil) = zs(nzb_soil) 656 657 DO k = nzb_soil+1, nzt_soil 658 dz_soil_stag(k) = zs(k) - zs(k-1) 594 659 ENDDO 595 660 596 DO k = 0, soil_layers-2661 DO k = nzb_soil, nzt_soil-1 597 662 dz_soil(k) = 0.5 * (dz_soil_stag(k+1) + dz_soil_stag(k)) 598 663 ENDDO 599 dz_soil(soil_layers-1) = dz_soil_stag(soil_layers-1) 600 601 ddz_soil = 1.0 / dz_soil 602 ddz_soil_stag = 1.0 / dz_soil_stag 603 ! 604 !-- Initialize soil 605 IF ( soil_type .NE. 0 ) THEN 606 alpha_VG = soil_pars(0,soil_type) 607 l_VG = soil_pars(1,soil_type) 608 n_VG = soil_pars(2,soil_type) 609 gamma_w_sat = soil_pars(3,soil_type) 610 m_sat = m_soil_pars(0,soil_type) 611 m_fc = m_soil_pars(1,soil_type) 612 m_wilt = m_soil_pars(2,soil_type) 613 m_res = m_soil_pars(3,soil_type) 664 dz_soil(nzt_soil) = dz_soil_stag(nzt_soil) 665 666 ddz_soil = 1.0_wp / dz_soil 667 ddz_soil_stag = 1.0_wp / dz_soil_stag 668 669 ! 670 !-- Initialize standard soil types. It is possible to overwrite each 671 !-- parameter by setting the respecticy NAMELIST variable to a 672 !-- value /= 9999999.9. 673 IF ( soil_type /= 0 ) THEN 674 675 IF ( alpha_vangenuchten == 9999999.9_wp ) THEN 676 alpha_vangenuchten = soil_pars(0,soil_type) 677 ENDIF 678 679 IF ( l_vangenuchten == 9999999.9_wp ) THEN 680 l_vangenuchten = soil_pars(1,soil_type) 681 ENDIF 682 683 IF ( n_vangenuchten == 9999999.9_wp ) THEN 684 n_vangenuchten = soil_pars(2,soil_type) 685 ENDIF 686 687 IF ( hydraulic_conductivity == 9999999.9_wp ) THEN 688 hydraulic_conductivity = soil_pars(3,soil_type) 689 ENDIF 690 691 IF ( saturation_moisture == 9999999.9_wp ) THEN 692 saturation_moisture = m_soil_pars(0,soil_type) 693 ENDIF 694 695 IF ( field_capacity == 9999999.9_wp ) THEN 696 field_capacity = m_soil_pars(1,soil_type) 697 ENDIF 698 699 IF ( wilting_point == 9999999.9_wp ) THEN 700 wilting_point = m_soil_pars(2,soil_type) 701 ENDIF 702 703 IF ( residual_moisture == 9999999.9_wp ) THEN 704 residual_moisture = m_soil_pars(3,soil_type) 705 ENDIF 706 707 ENDIF 708 709 alpha_vg = alpha_vangenuchten 710 l_vg = l_vangenuchten 711 n_vg = n_vangenuchten 712 gamma_w_sat = hydraulic_conductivity 713 m_sat = saturation_moisture 714 m_fc = field_capacity 715 m_wilt = wilting_point 716 m_res = residual_moisture 717 r_soil_min = min_soil_resistance 718 719 ! 720 !-- Initial run actions 721 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 722 723 t_soil = 0.0_wp 724 m_liq_eb = 0.0_wp 725 m_soil = 0.0_wp 726 727 ! 728 !-- Map user settings of T and q for each soil layer 729 !-- (make sure that the soil moisture does not drop below the permanent 730 !-- wilting point) -> problems with devision by zero) 731 DO k = nzb_soil, nzt_soil 732 t_soil(k,:,:) = soil_temperature(k) 733 m_soil(k,:,:) = MAX(soil_moisture(k),m_wilt(:,:)) 734 soil_moisture(k) = MAX(soil_moisture(k),wilting_point) 735 ENDDO 736 t_soil(nzt_soil+1,:,:) = soil_temperature(nzt_soil+1) 737 738 ! 739 !-- Calculate surface temperature 740 t_surface = pt_surface * exn 741 742 ! 743 !-- Set artifical values for ts and us so that r_a has its initial value for 744 !-- the first time step 745 DO i = nxlg, nxrg 746 DO j = nysg, nyng 747 k = nzb_s_inner(j,i) 748 749 ! 750 !-- Assure that r_a cannot be zero at model start 751 IF ( pt(k+1,j,i) == pt(k,j,i) ) pt(k+1,j,i) = pt(k+1,j,i) & 752 + 1.0E-10_wp 753 754 us(j,i) = 0.1_wp 755 ts(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / r_a(j,i) 756 shf(j,i) = - us(j,i) * ts(j,i) 757 ENDDO 758 ENDDO 759 760 ! 761 !-- Actions for restart runs 614 762 ELSE 615 alpha_VG = alpha_VanGenuchten 616 l_VG = l_VanGenuchten 617 n_VG = n_VanGenuchten 618 gamma_w_sat = hydraulic_conductivity 619 m_sat = saturation_moisture 620 m_fc = field_capacity 621 m_wilt = wilting_point 622 m_res = residual_moisture 623 ENDIF 624 625 ! 626 !-- Map user settings of T and q for each soil layer 627 !-- (make sure that the soil moisture does not drop below the permanent 628 !-- wilting point) -> problems with devision by zero) 629 DO k = 0, soil_layers-1 630 T_soil(k,:,:) = soil_temperature(k) 631 m_soil(k,:,:) = MAX(soil_moisture(k),m_wilt(:,:)) 632 ENDDO 633 T_soil(soil_layers,:,:) = soil_temperature(soil_layers) 634 635 636 exn = ( surface_pressure / 1000.0_wp )**0.286_wp 637 T_0 = pt_surface * exn 638 639 T_soil_p = T_soil 640 m_soil_p = m_soil 763 764 DO i = nxlg, nxrg 765 DO j = nysg, nyng 766 k = nzb_s_inner(j,i) 767 t_surface(j,i) = pt(k,j,i) * exn 768 ENDDO 769 ENDDO 770 771 ENDIF 641 772 642 773 ! … … 645 776 lambda_h_water ** m_sat(:,:) 646 777 778 779 780 781 DO k = nzb_soil, nzt_soil 782 root_fr(k,:,:) = root_fraction(k) 783 ENDDO 784 785 IF ( veg_type /= 0 ) THEN 786 IF ( min_canopy_resistance == 9999999.9_wp ) THEN 787 min_canopy_resistance = veg_pars(0,veg_type) 788 ENDIF 789 IF ( leaf_area_index == 9999999.9_wp ) THEN 790 leaf_area_index = veg_pars(1,veg_type) 791 ENDIF 792 IF ( vegetation_coverage == 9999999.9_wp ) THEN 793 vegetation_coverage = veg_pars(2,veg_type) 794 ENDIF 795 IF ( canopy_resistance_coefficient == 9999999.9_wp ) THEN 796 canopy_resistance_coefficient= veg_pars(3,veg_type) 797 ENDIF 798 IF ( lambda_surface_stable == 9999999.9_wp ) THEN 799 lambda_surface_stable = surface_pars(0,veg_type) 800 ENDIF 801 IF ( lambda_surface_unstable == 9999999.9_wp ) THEN 802 lambda_surface_unstable = surface_pars(1,veg_type) 803 ENDIF 804 IF ( f_shortwave_incoming == 9999999.9_wp ) THEN 805 f_shortwave_incoming = surface_pars(2,veg_type) 806 ENDIF 807 IF ( z0_eb == 9999999.9_wp ) THEN 808 roughness_length = roughness_par(0,veg_type) 809 z0_eb = roughness_par(0,veg_type) 810 ENDIF 811 IF ( z0h_eb == 9999999.9_wp ) THEN 812 z0h_eb = roughness_par(1,veg_type) 813 ENDIF 814 z0h_factor = z0h_eb / z0_eb 815 816 IF ( ANY( root_fraction == 9999999.9_wp ) ) THEN 817 DO k = nzb_soil, nzt_soil 818 root_fr(k,:,:) = root_distribution(k,veg_type) 819 root_fraction(k) = root_distribution(k,veg_type) 820 ENDDO 821 ENDIF 822 823 ENDIF 824 647 825 ! 648 826 !-- Initialize vegetation 649 IF ( veg_type .NE. 0 ) THEN 650 651 r_s_min = veg_pars(0,veg_type) 652 LAI = veg_pars(1,veg_type) 653 c_veg = veg_pars(2,veg_type) 654 gD = veg_pars(3,veg_type) 655 lambda_skin_s = skin_pars(0,veg_type) 656 lambda_skin_u = skin_pars(1,veg_type) 657 f_SW_in = skin_pars(2,veg_type) 658 z0 = roughness_par(0,veg_type) 659 z0h = roughness_par(1,veg_type) 660 661 662 DO k = 0, soil_layers-1 663 root_fr(k,:,:) = root_distribution(k,veg_type) 664 ENDDO 665 666 ELSE 667 668 DO k = 0, soil_layers-1 669 root_fr(k,:,:) = root_fraction(k) 670 ENDDO 671 672 ENDIF 827 r_canopy_min = min_canopy_resistance 828 lai = leaf_area_index 829 c_veg = vegetation_coverage 830 g_d = canopy_resistance_coefficient 831 lambda_surface_s = lambda_surface_stable 832 lambda_surface_u = lambda_surface_unstable 833 f_sw_in = f_shortwave_incoming 834 z0 = z0_eb 835 z0h = z0h_eb 673 836 674 837 ! … … 676 839 CALL user_init_land_surface 677 840 678 ! 679 !-- Set artifical values for ts and us so that r_a has its initial value for 680 !-- the first time step 681 DO i = nxlg, nxrg 682 DO j = nysg, nyng 683 k = nzb_s_inner(j,i) 684 ! 685 !-- Assure that r_a cannot be zero at model start 686 IF ( pt(k+1,j,i) == pt(k,j,i) ) pt(k+1,j,i) = pt(k+1,j,i) + 1.0E-10_wp 687 688 us(j,i) = 0.1_wp 689 ts(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / r_a(j,i) 690 shf(j,i) = - us(j,i) * ts(j,i) 691 ENDDO 692 ENDDO 841 842 t_soil_p = t_soil 843 m_soil_p = m_soil 844 m_liq_eb_p = m_liq_eb 845 846 !-- Store initial profiles of t_soil and m_soil (assuming they are 847 !-- horizontally homogeneous on this PE) 848 hom(nzb_soil:nzt_soil,1,90,:) = SPREAD( t_soil(nzb_soil:nzt_soil, & 849 nysg,nxlg), 2, & 850 statistic_regions+1 ) 851 hom(nzb_soil:nzt_soil,1,92,:) = SPREAD( m_soil(nzb_soil:nzt_soil, & 852 nysg,nxlg), 2, & 853 statistic_regions+1 ) 693 854 694 855 ! 695 856 !-- Calculate humidity at the surface 696 857 IF ( humidity ) THEN 697 CALL calc_q 0858 CALL calc_q_surface 698 859 ENDIF 860 861 862 863 ! 864 !-- Add timeseries for land surface model 865 dots_label(dots_num+1) = "ghf_eb" 866 dots_label(dots_num+2) = "shf_eb" 867 dots_label(dots_num+3) = "qsws_eb" 868 dots_label(dots_num+4) = "qsws_liq_eb" 869 dots_label(dots_num+5) = "qsws_soil_eb" 870 dots_label(dots_num+6) = "qsws_veg_eb" 871 dots_unit(dots_num+1:dots_num+6) = "W/m2" 872 873 dots_soil = dots_num + 1 874 dots_num = dots_num + 6 875 699 876 700 877 RETURN … … 707 884 ! Description: 708 885 ! ------------ 709 ! 886 ! Solver for the energy balance at the surface. 887 ! 888 ! Note: surface fluxes are calculated in the land surface model, but these are 889 ! not used in the atmospheric code. The fluxes are calculated afterwards in 890 ! prandtl_fluxes using the surface values of temperature and humidity as 891 ! provided by the land surface model. In this way, the fluxes in the land 892 ! surface model are not equal to the ones calculated in prandtl_fluxes 710 893 !------------------------------------------------------------------------------! 711 894 SUBROUTINE lsm_energy_balance … … 730 913 coef_1, & !: coef. for prognostic equation 731 914 coef_2, & !: coef. for prognostic equation 732 f_LE, & !: factor for LE 733 f_LE_veg, & !: factor for LE_veg 734 f_LE_soil, & !: factor for LE_soil 735 f_LE_liq, & !: factor for LE_liq 736 f_H, & !: factor for H 737 lambda_skin, & !: Current value of lambda_skin 738 m_liq_max !: maxmimum value of the liquid water reservoir 915 f_qsws, & !: factor for qsws_eb 916 f_qsws_veg, & !: factor for qsws_veg_eb 917 f_qsws_soil, & !: factor for qsws_soil_eb 918 f_qsws_liq, & !: factor for qsws_liq_eb 919 f_shf, & !: factor for shf_eb 920 lambda_surface, & !: Current value of lambda_surface 921 m_liq_eb_max !: maxmimum value of the liq. water reservoir 922 739 923 740 924 ! … … 748 932 749 933 ! 750 !-- Set lambda_s kinaccording to stratification934 !-- Set lambda_surface according to stratification 751 935 IF ( rif(j,i) >= 0.0_wp ) THEN 752 lambda_s kin = lambda_skin_s(j,i)936 lambda_surface = lambda_surface_s(j,i) 753 937 ELSE 754 lambda_s kin = lambda_skin_u(j,i)938 lambda_surface = lambda_surface_u(j,i) 755 939 ENDIF 756 940 … … 760 944 !-- time step is used here. Note that this formulation is the 761 945 !-- equivalent to the ECMWF formulation using drag coefficients 762 r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20) 946 r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / (ts(j,i) * us(j,i) + & 947 1.0E-20_wp) 763 948 764 949 ! … … 766 951 !-- f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation 767 952 768 !-- f1: correction for incoming shortwave radiation 769 f1 = MIN(1.0_wp, ( 0.004_wp * SW_in(j,i) + 0.05_wp ) / & 770 (0.81_wp * (0.004_wp * SW_in(j,i) + 1.0_wp) ) ) 771 772 ! 773 !-- f2: correction for soil moisture f2=0 for very dry soil 953 !-- f1: correction for incoming shortwave radiation (stomata close at 954 !-- night) 955 IF ( irad_scheme /= 0 ) THEN 956 f1 = MIN(1.0_wp, ( 0.004_wp * rad_sw_in(j,i) + 0.05_wp ) / & 957 (0.81_wp * (0.004_wp * rad_sw_in(j,i) + 1.0_wp) )) 958 ELSE 959 f1 = 1.0_wp 960 ENDIF 961 962 ! 963 !-- f2: correction for soil moisture availability to plants (the 964 !-- integrated soil moisture must thus be considered here) 965 !-- f2 = 0 for very dry soils 774 966 m_total = 0.0_wp 775 DO ks = 0, soil_layers-1 776 m_total = m_total + root_fr(ks,j,i) * m_soil(ks,j,i) 967 DO ks = nzb_soil, nzt_soil 968 m_total = m_total + root_fr(ks,j,i) & 969 * MAX(m_soil(ks,j,i),m_wilt(j,i)) 777 970 ENDDO 778 971 779 IF ( m_total .GT. m_wilt(j,i) .AND. m_total .LE. m_fc(j,i) ) THEN972 IF ( m_total .GT. m_wilt(j,i) .AND. m_total .LT. m_fc(j,i) ) THEN 780 973 f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) ) 974 ELSEIF ( m_total .GE. m_fc(j,i) ) THEN 975 f2 = 1.0_wp 781 976 ELSE 782 977 f2 = 1.0E-20_wp … … 785 980 ! 786 981 !-- Calculate water vapour pressure at saturation 787 !-- (T_0 should be replaced by liquid water temp?!) 788 e_s = 0.01 * 610.78_wp * EXP( 17.269_wp * ( T_0(j,i) - 273.16_wp )& 789 / ( T_0(j,i) - 35.86_wp ) ) 982 e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i) & 983 - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) ) 790 984 791 985 ! 792 986 !-- f3: correction for vapour pressure deficit 793 IF ( g D(j,i) .NE.0.0_wp ) THEN987 IF ( g_d(j,i) /= 0.0_wp ) THEN 794 988 ! 795 989 !-- Calculate vapour pressure 796 e = q_p(k+1,j,i) * surface_pressure / 0.622 797 f3 = EXP ( -g D(j,i) * (e_s - e) )990 e = q_p(k+1,j,i) * surface_pressure / 0.622_wp 991 f3 = EXP ( -g_d(j,i) * (e_s - e) ) 798 992 ELSE 799 993 f3 = 1.0_wp … … 801 995 802 996 ! 997 !-- Calculate canopy resistance. In case that c_veg is 0 (bare soils), 998 !-- this calculation is obsolete, as r_canopy is not used below. 803 999 !-- To do: check for very dry soil -> r_canopy goes to infinity 804 r_canopy(j,i) = r_s_min(j,i) / (LAI(j,i) * f1 * f2 * f3 + 1.0E-20) 805 806 ! 807 !-- Third step: calculate bare soil resistance r_soil 808 m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) * & 809 m_res(j,i) 810 811 f2 = ( m_soil(0,j,i) - m_min ) / ( m_fc(j,i) - m_min ) 1000 r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3 & 1001 + 1.0E-20_wp) 1002 1003 ! 1004 !-- Third step: calculate bare soil resistance r_soil. The Clapp & 1005 !-- Hornberger parametrization does not consider c_veg. 1006 IF ( soil_type /= 7 ) THEN 1007 m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) * & 1008 m_res(j,i) 1009 ELSE 1010 m_min = m_wilt(j,i) 1011 ENDIF 1012 1013 f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min ) 812 1014 f2 = MAX(f2,1.0E-20_wp) 1015 f2 = MIN(f2,1.0_wp) 813 1016 814 1017 r_soil(j,i) = r_soil_min(j,i) / f2 … … 816 1019 ! 817 1020 !-- Calculate fraction of liquid water reservoir 818 m_liq_max = m_max_depth * LAI(j,i) 819 c_liq(j,i) = MIN(1.0_wp, m_liq(j,i)/m_liq_max) 820 1021 m_liq_eb_max = m_max_depth * lai(j,i) 1022 c_liq(j,i) = MIN(1.0_wp, m_liq_eb(j,i) / (m_liq_eb_max+1.0E-20_wp)) 1023 1024 1025 ! 1026 !-- Calculate saturation specific humidity 821 1027 q_s = 0.622_wp * e_s / surface_pressure 822 1028 823 1029 ! 824 !-- In case of dew fall, set resistances to zero. 825 !-- To do: what does that physically reasoning behind this? 1030 !-- In case of dewfall, set evapotranspiration to zero 1031 !-- All super-saturated water is then removed from the air 1032 IF ( humidity .AND. dewfall .AND. q_s .LE. q_p(k+1,j,i) ) THEN 1033 r_canopy(j,i) = 0.0_wp 1034 r_soil(j,i) = 0.0_wp 1035 ENDIF 1036 1037 ! 1038 !-- Calculate coefficients for the total evapotranspiration 1039 f_qsws_veg = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/ & 1040 (r_a(j,i) + r_canopy(j,i)) 1041 1042 f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) + & 1043 r_soil(j,i)) 1044 f_qsws_liq = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i) 1045 1046 1047 ! 1048 !-- If soil moisture is below wilting point, plants do no longer 1049 !-- transpirate. 1050 ! IF ( m_soil(k,j,i) .LT. m_wilt(j,i) ) THEN 1051 ! f_qsws_veg = 0.0_wp 1052 ! ENDIF 1053 1054 ! 1055 !-- If dewfall is deactivated, vegetation, soil and liquid water 1056 !-- reservoir are not allowed to take up water from the super-saturated 1057 !-- air. 826 1058 IF ( humidity ) THEN 827 1059 IF ( q_s .LE. q_p(k+1,j,i) ) THEN 828 r_canopy(j,i) = 0.0_wp 829 r_soil(j,i) = 0.0_wp 1060 IF ( .NOT. dewfall ) THEN 1061 f_qsws_veg = 0.0_wp 1062 f_qsws_soil = 0.0_wp 1063 f_qsws_liq = 0.0_wp 1064 ENDIF 830 1065 ENDIF 831 1066 ENDIF 832 1067 833 834 ! 835 !-- Calculate coefficients for the total evapotranspiration 836 f_LE_veg = rho_lv * c_veg(j,i) * (1.0 - c_liq(j,i)) / (r_a(j,i) & 837 + r_canopy(j,i)) 838 f_LE_soil = rho_lv * (1.0 - c_veg(j,i)) / (r_a(j,i) + r_soil(j,i)) 839 f_LE_liq = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i) 840 841 842 ! 843 !-- If soil moisture is below wilting point, plants do no longer 844 !-- transpirate. 845 IF ( m_soil(k,j,i) .LT. m_wilt(j,i) ) THEN 846 f_LE_veg = 0.0 1068 f_shf = rho_cp / r_a(j,i) 1069 f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq 1070 1071 ! 1072 !-- Calculate derivative of q_s for Taylor series expansion 1073 e_s_dT = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) - & 1074 17.269_wp*(t_surface(j,i) - 273.16_wp) & 1075 / (t_surface(j,i) - 35.86_wp)**2 ) 1076 1077 dq_s_dT = 0.622_wp * e_s_dT / surface_pressure 1078 1079 T_1 = pt_p(k+1,j,i) * exn 1080 1081 ! 1082 !-- Add LW up so that it can be removed in prognostic equation 1083 rad_net(j,i) = rad_net(j,i) + sigma_SB * t_surface(j,i) ** 4 1084 1085 IF ( humidity ) THEN 1086 1087 ! 1088 !-- Numerator of the prognostic equation 1089 coef_1 = rad_net(j,i) + 3.0_wp * sigma_SB * t_surface(j,i) ** 4& 1090 + f_shf / exn * T_1 + f_qsws * ( q_p(k+1,j,i) - q_s & 1091 + dq_s_dT * t_surface(j,i) ) + lambda_surface & 1092 * t_soil(nzb_soil,j,i) 1093 1094 ! 1095 !-- Denominator of the prognostic equation 1096 coef_2 = 4.0_wp * sigma_SB * t_surface(j,i) ** 3 + f_qsws & 1097 * dq_s_dT + lambda_surface + f_shf / exn 1098 1099 ELSE 1100 1101 ! 1102 !-- Numerator of the prognostic equation 1103 coef_1 = rad_net(j,i) + 3.0_wp * sigma_SB * t_surface(j,i) ** 4& 1104 + f_shf / exn * T_1 + lambda_surface & 1105 * t_soil(nzb_soil,j,i) 1106 1107 ! 1108 !-- Denominator of the prognostic equation 1109 coef_2 = 4.0_wp * sigma_SB * t_surface(j,i) ** 3 & 1110 + lambda_surface + f_shf / exn 1111 847 1112 ENDIF 848 1113 849 f_H = rho_cp / r_a(j,i)850 f_LE = f_LE_veg + f_LE_soil + f_LE_liq851 852 !853 !-- Calculate derivative of q_s for Taylor series expansion854 e_s_dT = e_s * ( 17.269_wp / (T_0(j,i) - 35.86_wp) - &855 17.269_wp*(T_0(j,i) - 273.16_wp) / (T_0(j,i) &856 - 35.86_wp)**2 )857 858 dq_s_dT = 0.622_wp * e_s_dT / surface_pressure859 860 T_1 = pt_p(k+1,j,i) * exn861 862 !863 !-- Add LW up so that it can be removed in prognostic equation864 Rn(j,i) = Rn(j,i) + sigma_SB * T_0(j,i) ** 4865 866 IF ( humidity ) THEN867 868 !869 !-- Numerator of the prognostic equation870 coef_1 = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i) ** 4 + f_H &871 / exn * T_1 + f_LE * ( q_p(k+1,j,i) - q_s + dq_s_dT &872 * T_0(j,i) ) + lambda_skin * T_soil(0,j,i)873 874 !875 !-- Denominator of the prognostic equation876 coef_2 = 4.0_wp * sigma_SB * T_0(j,i) ** 3 + f_LE * dq_s_dT &877 + lambda_skin + f_H / exn878 879 ELSE880 881 !882 !-- Numerator of the prognostic equation883 coef_1 = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i) ** 4 + f_H / &884 exn * T_1 + lambda_skin * T_soil(0,j,i)885 886 !887 !-- Denominator of the prognostic equation888 coef_2 = 4.0_wp * sigma_SB * T_0(j,i) ** 3 &889 + lambda_skin + f_H / exn890 891 ENDIF892 893 1114 tend = 0.0_wp 894 1115 895 1116 ! 896 !-- Implicit solution when the s kinlayer has no heat capacity,1117 !-- Implicit solution when the surface layer has no heat capacity, 897 1118 !-- otherwise use RK3 scheme. 898 T_0_p(j,i) = ( coef_1 * dt_3d * tsc(2) + C_skin * T_0(j,i) ) /&899 ( C_skin + coef_2 * dt_3d * tsc(2) )900 1119 t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface * & 1120 t_surface(j,i) ) / ( c_surface + coef_2 * dt_3d& 1121 * tsc(2) ) 901 1122 ! 902 1123 !-- Add RK3 term 903 T_0_p(j,i) = T_0_p(j,i) + dt_3d * tsc(3) * tT_soil_m(0,j,i) 1124 t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3) & 1125 * tt_surface_m(j,i) 904 1126 905 1127 ! 906 1128 !-- Calculate true tendency 907 tend = (T_0_p(j,i) - T_0(j,i) - tsc(3) * tT_0_m(j,i)) / (dt_3d & 908 * tsc(2)) 909 910 ! 911 !-- Calculate T_0 tendencies for the next Runge-Kutta step 1129 tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3) & 1130 * tt_surface_m(j,i)) / (dt_3d * tsc(2)) 1131 ! 1132 !-- Calculate t_surface tendencies for the next Runge-Kutta step 912 1133 IF ( timestep_scheme(1:5) == 'runge' ) THEN 913 1134 IF ( intermediate_timestep_count == 1 ) THEN 914 t T_0_m(j,i) = tend1135 tt_surface_m(j,i) = tend 915 1136 ELSEIF ( intermediate_timestep_count < & 916 1137 intermediate_timestep_count_max ) THEN 917 tT_0_m(j,i) = -9.5625_wp * tend + 5.3125_wp * tT_0_m(j,i) 1138 tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp & 1139 * tt_surface_m(j,i) 918 1140 ENDIF 919 1141 ENDIF 920 1142 921 pt_p(k,j,i) = T_0_p(j,i) / exn1143 pt_p(k,j,i) = t_surface_p(j,i) / exn 922 1144 ! 923 1145 !-- Calculate fluxes 924 Rn(j,i) = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i)**4 & 925 - 4.0_wp * sigma_SB * T_0(j,i)**3 * T_0_p(j,i) 926 G(j,i) = lambda_skin * (T_0_p(j,i) - T_soil(0,j,i)) 927 H(j,i) = - f_H * ( pt_p(k+1,j,i) - pt_p(k,j,i) ) 1146 rad_net(j,i) = rad_net(j,i) + 3.0_wp * sigma_SB & 1147 * t_surface(j,i)**4 - 4.0_wp * sigma_SB & 1148 * t_surface(j,i)**3 * t_surface_p(j,i) 1149 ghf_eb(j,i) = lambda_surface * (t_surface_p(j,i) & 1150 - t_soil(nzb_soil,j,i)) 1151 shf_eb(j,i) = - f_shf * ( pt_p(k+1,j,i) - pt_p(k,j,i) ) 928 1152 929 1153 IF ( humidity ) THEN 930 LE(j,i) = - f_LE * ( q_p(k+1,j,i) - q_s + dq_s_dT & 931 * T_0(j,i) - dq_s_dT * T_0_p(j,i) ) 932 933 LE_veg(j,i) = - f_LE_veg * ( q_p(k+1,j,i) - q_s + dq_s_dT & 934 * T_0(j,i) - dq_s_dT * T_0_p(j,i) ) 935 LE_soil(j,i) = - f_LE_soil * ( q_p(k+1,j,i) - q_s + dq_s_dT & 936 * T_0(j,i) - dq_s_dT * T_0_p(j,i) ) 937 LE_liq(j,i) = - f_LE_liq * ( q_p(k+1,j,i) - q_s + dq_s_dT & 938 * T_0(j,i) - dq_s_dT * T_0_p(j,i) ) 1154 qsws_eb(j,i) = - f_qsws * ( q_p(k+1,j,i) - q_s + dq_s_dT & 1155 * t_surface(j,i) - dq_s_dT * t_surface_p(j,i) ) 1156 1157 qsws_veg_eb(j,i) = - f_qsws_veg * ( q_p(k+1,j,i) - q_s & 1158 + dq_s_dT * t_surface(j,i) - dq_s_dT & 1159 * t_surface_p(j,i) ) 1160 1161 qsws_soil_eb(j,i) = - f_qsws_soil * ( q_p(k+1,j,i) - q_s & 1162 + dq_s_dT * t_surface(j,i) - dq_s_dT & 1163 * t_surface_p(j,i) ) 1164 1165 qsws_liq_eb(j,i) = - f_qsws_liq * ( q_p(k+1,j,i) - q_s & 1166 + dq_s_dT * t_surface(j,i) - dq_s_dT & 1167 * t_surface_p(j,i) ) 939 1168 ENDIF 940 1169 941 ! IF ( i == 1 .AND. j == 1 ) THEN942 ! PRINT*, "Rn", Rn(j,i)943 ! PRINT*, "H", H(j,i)944 ! PRINT*, "LE", LE(j,i)945 ! PRINT*, "LE_liq", LE_liq(j,i)946 ! PRINT*, "LE_veg", LE_veg(j,i)947 ! PRINT*, "LE_soil", LE_soil(j,i)948 ! PRINT*, "G", G(j,i)949 ! ENDIF950 951 1170 ! 952 1171 !-- Calculate the true surface resistance 953 IF ( LE(j,i) .EQ. 0.0) THEN954 r_s(j,i) = 1.0E10 1172 IF ( qsws_eb(j,i) .EQ. 0.0_wp ) THEN 1173 r_s(j,i) = 1.0E10_wp 955 1174 ELSE 956 r_s(j,i) = - rho_lv * ( q_p(k+1,j,i) - q_s + dq_s_dT * T_0(j,i)& 957 - dq_s_dT * T_0_p(j,i) ) / LE(j,i) - r_a(j,i) 1175 r_s(j,i) = - rho_lv * ( q_p(k+1,j,i) - q_s + dq_s_dT & 1176 * t_surface(j,i) - dq_s_dT * t_surface_p(j,i) ) & 1177 / qsws_eb(j,i) - r_a(j,i) 958 1178 ENDIF 959 960 !961 !-- Calculate fluxes in the atmosphere962 shf(j,i) = H(j,i) / rho_cp963 1179 964 1180 ! … … 967 1183 IF ( humidity ) THEN 968 1184 ! 969 !-- If precipitation is activated, add rain water to LE_liq.1185 !-- If precipitation is activated, add rain water to qsws_liq_eb. 970 1186 !-- precipitation_rate is given in mm. 971 1187 IF ( precipitation ) THEN 972 LE_liq(j,i) = LE_liq(j,i) + precipitation_rate(j,i) & 973 * 0.001_wp * rho_l * l_v 1188 qsws_liq_eb(j,i) = qsws_liq_eb(j,i) & 1189 + precipitation_rate(j,i) * 0.001_wp & 1190 * rho_l * l_v 974 1191 ENDIF 975 1192 ! … … 977 1194 IF ( q_s .LE. q_p(k+1,j,i)) THEN 978 1195 ! 979 !-- Check if reservoir is full (avoid values > m_liq_max) 980 !-- In that case, LE_liq goes to LE_soil. In this case 981 !-- LE_veg is zero anyway (because c_liq = 1), so that tend is 982 !-- zero and no further check is needed 983 IF ( m_liq(j,i) .EQ. m_liq_max ) THEN 984 LE_soil(j,i) = LE_soil(j,i) + LE_liq(j,i) 985 LE_liq(j,i) = 0.0_wp 1196 !-- Check if reservoir is full (avoid values > m_liq_eb_max) 1197 !-- In that case, qsws_liq_eb goes to qsws_soil_eb. In this 1198 !-- case qsws_veg_eb is zero anyway (because c_liq = 1), 1199 !-- so that tend is zero and no further check is needed 1200 IF ( m_liq_eb(j,i) .EQ. m_liq_eb_max ) THEN 1201 qsws_soil_eb(j,i) = qsws_soil_eb(j,i) & 1202 + qsws_liq_eb(j,i) 1203 qsws_liq_eb(j,i) = 0.0_wp 986 1204 ENDIF 987 1205 988 1206 ! 989 !-- In case LE_veg becomes negative (unphysical behavior), let990 !-- the water enter the liquid water reservoir as dew on the1207 !-- In case qsws_veg_eb becomes negative (unphysical behavior), 1208 !-- let the water enter the liquid water reservoir as dew on the 991 1209 !-- plant 992 IF ( LE_veg(j,i) .LT. 0.0_wp ) THEN993 LE_liq(j,i) = LE_liq(j,i) + LE_veg(j,i)994 LE_veg(j,i) = 0.0_wp1210 IF ( qsws_veg_eb(j,i) .LT. 0.0_wp ) THEN 1211 qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i) 1212 qsws_veg_eb(j,i) = 0.0_wp 995 1213 ENDIF 996 1214 ENDIF 997 1215 998 tend = - LE_liq(j,i) * drho_l_lv999 1000 m_liq_ p(j,i) = m_liq(j,i) + dt_3d * ( tsc(2) * tend&1001 + tsc(3) * tm_liq_ m(j,i) )1216 tend = - qsws_liq_eb(j,i) * drho_l_lv 1217 1218 m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend & 1219 + tsc(3) * tm_liq_eb_m(j,i) ) 1002 1220 1003 1221 ! 1004 1222 !-- Check if reservoir is overfull -> reduce to maximum 1005 1223 !-- (conservation of water is violated here) 1006 m_liq_ p(j,i) = MIN(m_liq_p(j,i),m_liq_max)1224 m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max) 1007 1225 1008 1226 ! 1009 1227 !-- Check if reservoir is empty (avoid values < 0.0) 1010 1228 !-- (conservation of water is violated here) 1011 m_liq_ p(j,i) = MAX(m_liq_p(j,i),0.0_wp)1012 1013 1014 ! 1015 !-- Calculate m_liq tendencies for the next Runge-Kutta step1229 m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp) 1230 1231 1232 ! 1233 !-- Calculate m_liq_eb tendencies for the next Runge-Kutta step 1016 1234 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1017 1235 IF ( intermediate_timestep_count == 1 ) THEN 1018 tm_liq_ m(j,i) = tend1236 tm_liq_eb_m(j,i) = tend 1019 1237 ELSEIF ( intermediate_timestep_count < & 1020 1238 intermediate_timestep_count_max ) THEN 1021 tm_liq_ m(j,i) = -9.5625_wp * tend + 5.3125_wp&1022 * tm_liq_ m(j,i)1239 tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp & 1240 * tm_liq_eb_m(j,i) 1023 1241 ENDIF 1024 1242 ENDIF 1025 1243 1026 !1027 !-- Calculate moisture flux in the atmosphere1028 qsws(j,i) = LE(j,i) / rho_lv1029 1030 1244 ENDIF 1031 1245 1032 1246 ENDDO 1033 1247 ENDDO 1034 1035 1036 1248 1037 1249 END SUBROUTINE lsm_energy_balance … … 1042 1254 ! ------------ 1043 1255 ! 1256 ! Soil model as part of the land surface model. The model predicts soil 1257 ! temperature and water content. 1044 1258 !------------------------------------------------------------------------------! 1045 1259 SUBROUTINE lsm_soil_model … … 1052 1266 INTEGER(iwp) :: k !: running index 1053 1267 1054 REAL(wp) :: h_ VG!: Van Genuchten coef. h1055 1056 REAL(wp), DIMENSION( 0:soil_layers-1) :: gamma_temp, & !: temp. gamma1268 REAL(wp) :: h_vg !: Van Genuchten coef. h 1269 1270 REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp, & !: temp. gamma 1057 1271 lambda_temp, & !: temp. lambda 1058 1272 tend !: tendency … … 1060 1274 DO i = nxlg, nxrg 1061 1275 DO j = nysg, nyng 1062 DO k = 0, soil_layers-11276 DO k = nzb_soil, nzt_soil 1063 1277 ! 1064 1278 !-- Calculate volumetric heat capacity of the soil, taking into 1065 1279 !-- account water content 1066 rho C_total(k,j,i) = (rhoC_soil * (1.0 - m_sat(j,i))&1067 + rho C_water * m_soil(k,j,i))1280 rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i)) & 1281 + rho_c_water * m_soil(k,j,i)) 1068 1282 1069 1283 ! 1070 1284 !-- Calculate soil heat conductivity at the center of the soil 1071 1285 !-- layers 1072 Ke = 1.0+ LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i)))1073 lambda_temp(k) = Ke * (lambda_h_sat(j,i) + lambda_h_dry) + &1286 ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i))) 1287 lambda_temp(k) = ke * (lambda_h_sat(j,i) + lambda_h_dry) + & 1074 1288 lambda_h_dry 1075 1289 … … 1079 1293 !-- Calculate soil heat conductivity (lambda_h) at the _stag level 1080 1294 !-- using linear interpolation 1081 DO k = 0, soil_layers-21295 DO k = nzb_soil, nzt_soil-1 1082 1296 1083 1297 lambda_h(k,j,i) = lambda_temp(k) + & … … 1086 1300 1087 1301 ENDDO 1088 lambda_h( soil_layers-1,j,i) = lambda_temp(soil_layers-1)1089 1090 ! 1091 !-- Prognostic equation for soil temperature T_soil1302 lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil) 1303 1304 ! 1305 !-- Prognostic equation for soil temperature t_soil 1092 1306 tend(:) = 0.0_wp 1093 tend(0) = (1.0/rhoC_total(0,j,i)) * & 1094 ( lambda_h(0,j,i) * ( T_soil(1,j,i) - T_soil(0,j,i) ) & 1095 * ddz_soil(0) + G(j,i) ) * ddz_soil_stag(0) 1096 1097 DO k = 1, soil_layers-1 1098 tend(k) = (1.0/rhoC_total(k,j,i)) & 1307 tend(0) = (1.0_wp/rho_c_total(nzb_soil,j,i)) * & 1308 ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i) & 1309 - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil) & 1310 + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil) 1311 1312 DO k = 1, nzt_soil 1313 tend(k) = (1.0_wp/rho_c_total(k,j,i)) & 1099 1314 * ( lambda_h(k,j,i) & 1100 * ( T_soil(k+1,j,i) - T_soil(k,j,i) ) &1315 * ( t_soil(k+1,j,i) - t_soil(k,j,i) ) & 1101 1316 * ddz_soil(k) & 1102 1317 - lambda_h(k-1,j,i) & 1103 * ( T_soil(k,j,i) - T_soil(k-1,j,i) ) &1318 * ( t_soil(k,j,i) - t_soil(k-1,j,i) ) & 1104 1319 * ddz_soil(k-1) & 1105 1320 ) * ddz_soil_stag(k) 1106 1321 ENDDO 1107 1322 1108 T_soil_p(0:soil_layers-1,j,i) = T_soil(0:soil_layers-1,j,i)&1323 t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i) & 1109 1324 + dt_3d * ( tsc(2) & 1110 1325 * tend(:) + tsc(3) & 1111 * t T_soil_m(:,j,i) )1112 1113 ! 1114 !-- Calculate T_soil tendencies for the next Runge-Kutta step1326 * tt_soil_m(:,j,i) ) 1327 1328 ! 1329 !-- Calculate t_soil tendencies for the next Runge-Kutta step 1115 1330 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1116 1331 IF ( intermediate_timestep_count == 1 ) THEN 1117 DO k = 0, soil_layers-11118 t T_soil_m(k,j,i) = tend(k)1332 DO k = nzb_soil, nzt_soil 1333 tt_soil_m(k,j,i) = tend(k) 1119 1334 ENDDO 1120 1335 ELSEIF ( intermediate_timestep_count < & 1121 1336 intermediate_timestep_count_max ) THEN 1122 DO k = 0, soil_layers-11123 t T_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp &1124 * t T_soil_m(k,j,i)1337 DO k = nzb_soil, nzt_soil 1338 tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp & 1339 * tt_soil_m(k,j,i) 1125 1340 ENDDO 1126 1341 ENDIF … … 1128 1343 1129 1344 1130 DO k = 0, soil_layers-1 1345 DO k = nzb_soil, nzt_soil 1346 1131 1347 ! 1132 1348 !-- Calculate soil diffusivity at the center of the soil layers 1133 lambda_temp(k) = (- b_ CH* gamma_w_sat(j,i) * psi_sat &1349 lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat & 1134 1350 / m_sat(j,i) ) * ( MAX(m_soil(k,j,i), & 1135 m_wilt(j,i)) / m_sat(j,i) )**(b_CH + 2.0_wp) 1136 1137 ! 1138 !-- Calculate the hydraulic conductivity after Van Genuchten (1980) 1139 h_VG = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) - & 1140 MAX(m_soil(k,j,i),m_wilt(j,i)) ) )**(n_VG(j,i) & 1141 / (n_VG(j,i)-1.0_wp)) - 1.0_wp & 1142 )**(1.0_wp/n_VG(j,i)) / alpha_VG(j,i) 1143 1144 gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp + & 1145 (alpha_VG(j,i)*h_VG)**n_VG(j,i))**(1.0_wp & 1146 -1.0_wp/n_VG(j,i)) - (alpha_VG(j,i)*h_VG & 1147 )**(n_VG(j,i)-1.0_wp))**2 ) & 1148 / ( (1.0_wp + (alpha_VG(j,i)*h_VG)**n_VG(j,i) & 1149 )**((1.0_wp - 1.0_wp/n_VG(j,i))*(l_VG(j,i) & 1150 + 2.0)) ) 1351 m_wilt(j,i)) / m_sat(j,i) )**(b_ch + 2.0_wp) 1352 1353 ! 1354 !-- Parametrization of Van Genuchten 1355 IF ( soil_type /= 7 ) THEN 1356 ! 1357 !-- Calculate the hydraulic conductivity after Van Genuchten 1358 !-- (1980) 1359 h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) - & 1360 MAX(m_soil(k,j,i),m_wilt(j,i)) ) )**(n_vg(j,i) & 1361 / (n_vg(j,i)-1.0_wp)) - 1.0_wp & 1362 )**(1.0_wp/n_vg(j,i)) / alpha_vg(j,i) 1363 1364 gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp + & 1365 (alpha_vg(j,i)*h_vg)**n_vg(j,i))**(1.0_wp & 1366 -1.0_wp/n_vg(j,i)) - (alpha_vg(j,i)*h_vg & 1367 )**(n_vg(j,i)-1.0_wp))**2 ) & 1368 / ( (1.0_wp + (alpha_vg(j,i)*h_vg & 1369 )**n_vg(j,i))**((1.0_wp - 1.0_wp/n_vg(j,i)) & 1370 *(l_vg(j,i) + 2.0_wp)) ) 1371 1372 ! 1373 !-- Parametrization of Clapp & Hornberger 1374 ELSE 1375 gamma_temp(k) = gamma_w_sat(j,i) * (m_soil(k,j,i) & 1376 / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp) 1377 ENDIF 1151 1378 1152 1379 ENDDO … … 1156 1383 ! 1157 1384 !-- Calculate soil diffusivity (lambda_w) at the _stag level 1158 !-- using linear interpolation 1159 DO k = 0, soil_layers-2 1385 !-- using linear interpolation. To do: replace this with 1386 !-- ECMWF-IFS Eq. 8.81 1387 DO k = nzb_soil, nzt_soil-1 1160 1388 1161 1389 lambda_w(k,j,i) = lambda_temp(k) + & 1162 1390 ( lambda_temp(k+1) - lambda_temp(k) ) & 1163 * 0.5 * dz_soil_stag(k) * ddz_soil(k+1)1391 * 0.5_wp * dz_soil_stag(k) * ddz_soil(k+1) 1164 1392 gamma_w(k,j,i) = gamma_temp(k) + & 1165 1393 ( gamma_temp(k+1) - gamma_temp(k) ) & 1166 * 0.5 * dz_soil_stag(k) * ddz_soil(k+1)1394 * 0.5_wp * dz_soil_stag(k) * ddz_soil(k+1) 1167 1395 1168 1396 ENDDO … … 1174 1402 !-- in the bottom layer. 1175 1403 IF ( conserve_water_content ) THEN 1176 gamma_w( soil_layers-1,j,i) = 0.0_wp1404 gamma_w(nzt_soil,j,i) = 0.0_wp 1177 1405 ELSE 1178 gamma_w( soil_layers-1,j,i) = lambda_temp(soil_layers-1)1406 gamma_w(nzt_soil,j,i) = lambda_temp(nzt_soil) 1179 1407 ENDIF 1180 1408 1181 !-- The root extraction (= root_extr * LE_veg/ (rho_l * l_v))1409 !-- The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v)) 1182 1410 !-- ensures the mass conservation for water. The transpiration of 1183 1411 !-- plants equals the cumulative withdrawals by the roots in the 1184 1412 !-- soil. The scheme takes into account the availability of water 1185 1413 !-- in the soil layers as well as the root fraction in the 1186 !-- respective layer 1187 1188 ! 1189 !-- Calculate the root extraction (ECMWF 7.69, with some 1190 !-- modifications) 1414 !-- respective layer. Layer with moisture below wilting point will 1415 !-- not contribute, which reflects the preference of plants to 1416 !-- take water from moister layers. 1417 1418 ! 1419 !-- Calculate the root extraction (ECMWF 7.69, modified so that the 1420 !-- sum of root_extr = 1). The energy balance solver guarantees a 1421 !-- positive transpiration, so that there is no need for an 1422 !-- additional check. 1423 1191 1424 m_total = 0.0_wp 1192 DO k = 0, soil_layers-11193 m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i) * &1194 dz_soil_stag(k)1195 1425 DO k = nzb_soil, nzt_soil 1426 IF ( m_soil(k,j,i) .GT. m_wilt(j,i) ) THEN 1427 m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i) 1428 ENDIF 1196 1429 ENDDO 1197 1430 1198 ! 1199 !-- For conservation of mass, the sum of root_extr must be 1 1200 DO k = 0, soil_layers-1 1201 root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) & 1202 * dz_soil_stag(k) / m_total 1431 DO k = nzb_soil, nzt_soil 1432 IF ( m_soil(k,j,i) .GT. m_wilt(j,i) ) THEN 1433 root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total 1434 ELSE 1435 root_extr(k) = 0.0_wp 1436 ENDIF 1203 1437 ENDDO 1204 1205 1438 1206 1439 ! 1207 1440 !-- Prognostic equation for soil water content m_soil 1208 1441 tend(:) = 0.0_wp 1209 tend(0) = ( lambda_w(0,j,i) * ( m_soil(1,j,i) - m_soil(0,j,i) )& 1210 * ddz_soil(0) - gamma_w(0,j,i) - ( root_extr(0) & 1211 * LE_veg(j,i) + LE_soil(j,i) ) * drho_l_lv & 1212 ) * ddz_soil_stag(0) 1213 1214 DO k = 1, soil_layers-2 1442 tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * ( & 1443 m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) ) & 1444 * ddz_soil(nzb_soil) - gamma_w(nzb_soil,j,i) - ( & 1445 root_extr(nzb_soil) * qsws_veg_eb(j,i) & 1446 + qsws_soil_eb(j,i) ) * drho_l_lv ) & 1447 * ddz_soil_stag(nzb_soil) 1448 1449 DO k = nzb_soil+1, nzt_soil-1 1215 1450 tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i) & 1216 1451 - m_soil(k,j,i) ) * ddz_soil(k) - gamma_w(k,j,i)& 1217 1452 - lambda_w(k-1,j,i) * (m_soil(k,j,i) - & 1218 1453 m_soil(k-1,j,i)) * ddz_soil(k-1) & 1219 + gamma_w(k-1,j,i) - (root_extr(k) * LE_veg(j,i)&1220 * drho_l_lv)&1454 + gamma_w(k-1,j,i) - (root_extr(k) & 1455 * qsws_veg_eb(j,i) * drho_l_lv) & 1221 1456 ) * ddz_soil_stag(k) 1222 1457 1223 1458 ENDDO 1224 tend( soil_layers-1) = ( - gamma_w(soil_layers-1,j,i)&1225 - lambda_w( soil_layers-2,j,i)&1226 * (m_soil( soil_layers-1,j,i)&1227 - m_soil( soil_layers-2,j,i))&1228 * ddz_soil( soil_layers-2)&1229 + gamma_w( soil_layers-2,j,i) - (&1230 root_extr( soil_layers-1)&1231 * LE_veg(j,i) * drho_l_lv )&1232 ) * ddz_soil_stag( soil_layers-1)1233 1234 m_soil_p( 0:soil_layers-1,j,i) = m_soil(0:soil_layers-1,j,i)&1459 tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i) & 1460 - lambda_w(nzt_soil-1,j,i) & 1461 * (m_soil(nzt_soil,j,i) & 1462 - m_soil(nzt_soil-1,j,i)) & 1463 * ddz_soil(nzt_soil-1) & 1464 + gamma_w(nzt_soil-1,j,i) - ( & 1465 root_extr(nzt_soil) & 1466 * qsws_veg_eb(j,i) * drho_l_lv ) & 1467 ) * ddz_soil_stag(nzt_soil) 1468 1469 m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)& 1235 1470 + dt_3d * ( tsc(2) * tend(:) & 1236 1471 + tsc(3) * tm_soil_m(:,j,i) ) … … 1244 1479 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1245 1480 IF ( intermediate_timestep_count == 1 ) THEN 1246 DO k = 0, soil_layers-11481 DO k = nzb_soil, nzt_soil 1247 1482 tm_soil_m(k,j,i) = tend(k) 1248 1483 ENDDO 1249 1484 ELSEIF ( intermediate_timestep_count < & 1250 1485 intermediate_timestep_count_max ) THEN 1251 DO k = 0, soil_layers-11486 DO k = nzb_soil, nzt_soil 1252 1487 tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp & 1253 1488 * tm_soil_m(k,j,i) … … 1264 1499 !-- Calculate surface specific humidity 1265 1500 IF ( humidity ) THEN 1266 CALL calc_q 01501 CALL calc_q_surface 1267 1502 ENDIF 1268 1503 … … 1274 1509 ! Description: 1275 1510 ! ------------ 1276 ! 1511 ! Calculation of specific humidity of the surface layer (surface) 1277 1512 !------------------------------------------------------------------------------! 1278 SUBROUTINE calc_q 01513 SUBROUTINE calc_q_surface 1279 1514 1280 1515 IMPLICIT NONE … … 1288 1523 DO j = nysg, nyng 1289 1524 k = nzb_s_inner(j,i) 1290 ! 1291 !-- Temporary solution as long as T_0 is prescribed 1292 1293 pt_p(k,j,i) = T_0(j,i) / exn 1525 1294 1526 ! 1295 1527 !-- Calculate water vapour pressure at saturation 1296 e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( T_0(j,i) -&1297 273.16_wp ) / ( T_0(j,i) -&1298 35.86_wp ) )1528 e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i) & 1529 - 273.16_wp ) / ( t_surface(j,i) & 1530 - 35.86_wp ) ) 1299 1531 1300 1532 ! … … 1313 1545 ENDDO 1314 1546 1315 END SUBROUTINE calc_q0 1547 END SUBROUTINE calc_q_surface 1548 1549 !------------------------------------------------------------------------------! 1550 ! Description: 1551 ! ------------ 1552 ! Swapping of timelevels 1553 !------------------------------------------------------------------------------! 1554 SUBROUTINE lsm_swap_timelevel ( mod_count ) 1555 1556 IMPLICIT NONE 1557 1558 INTEGER, INTENT(IN) :: mod_count 1559 1560 #if defined( __nopointer ) 1561 1562 t_surface = t_surface_p 1563 t_soil = t_soil_p 1564 IF ( humidity ) THEN 1565 m_soil = m_soil_p 1566 m_liq_eb = m_liq_eb_p 1567 ENDIF 1568 1569 #else 1570 1571 SELECT CASE ( mod_count ) 1572 1573 CASE ( 0 ) 1574 1575 t_surface = t_surface_p 1576 t_soil = t_soil_p 1577 IF ( humidity ) THEN 1578 m_soil = m_soil_p 1579 m_liq_eb = m_liq_eb_p 1580 ENDIF 1581 1582 1583 CASE ( 1 ) 1584 1585 t_surface => t_surface_1; t_surface_p => t_surface_2 1586 t_soil => t_soil_1; t_soil_p => t_soil_2 1587 IF ( humidity ) THEN 1588 m_soil => m_soil_1; m_soil_p => m_soil_2 1589 m_liq_eb => m_liq_eb_1; m_liq_eb_p => m_liq_eb_2 1590 ENDIF 1591 1592 END SELECT 1593 #endif 1594 1595 END SUBROUTINE lsm_swap_timelevel 1316 1596 1317 1597
Note: See TracChangeset
for help on using the changeset viewer.
