source: palm/trunk/SOURCE/land_surface_model.f90 @ 1788

Last change on this file since 1788 was 1788, checked in by maronga, 6 years ago

added support for water and paved surfaced in land surface model / minor changes

  • Property svn:keywords set to Id
File size: 85.6 KB
Line 
1!> @file land_surface_model.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2015 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21! Bugfix: calculate lambda_surface based on temperature gradient between skin
22! layer and soil layer instead of Obukhov length
23! Changed: moved calculation of surface specific humidity to energy balance solver
24! New: water surfaces are available by using a fixed sea surface temperature.
25! The roughness lengths are calculated dynamically using the Charnock
26! parameterization. This involves the new roughness length for moisture z0q.
27! New: modified solution of the energy balance solver and soil model for
28! paved surfaces (i.e. asphalt concrete).
29! Syntax layout improved.
30! Changed: parameter dewfall removed.
31!
32! Former revisions:
33! -----------------
34! $Id: land_surface_model.f90 1788 2016-03-10 11:01:04Z maronga $
35!
36! 1783 2016-03-06 18:36:17Z raasch
37! netcdf variables moved to netcdf module
38!
39! 1757 2016-02-22 15:49:32Z maronga
40! Bugfix: set tm_soil_m to zero after allocation. Added parameter
41! unscheduled_radiation_calls to control calls of the radiation model based on
42! the skin temperature change during one time step (preliminary version). Set
43! qsws_soil_eb to zero at model start (previously set to qsws_eb). Removed MAX
44! function as it cannot be vectorized.
45!
46! 1709 2015-11-04 14:47:01Z maronga
47! Renamed pt_1 and qv_1 to pt1 and qv1.
48! Bugfix: set initial values for t_surface_p in case of restart runs
49! Bugfix: zero resistance caused crash when using radiation_scheme = 'clear-sky'
50! Bugfix: calculation of rad_net when using radiation_scheme = 'clear-sky'
51! Added todo action
52!
53! 1697 2015-10-28 17:14:10Z raasch
54! bugfix: misplaced cpp-directive
55!
56! 1695 2015-10-27 10:03:11Z maronga
57! Bugfix: REAL constants provided with KIND-attribute in call of
58! Replaced rif with ol
59!
60! 1691 2015-10-26 16:17:44Z maronga
61! Added skip_time_do_lsm to allow for spin-ups without LSM. Various bugfixes:
62! Soil temperatures are now defined at the edges of the layers, calculation of
63! shb_eb corrected, prognostic equation for skin temperature corrected. Surface
64! fluxes are now directly transfered to atmosphere
65!
66! 1682 2015-10-07 23:56:08Z knoop
67! Code annotations made doxygen readable
68!
69! 1590 2015-05-08 13:56:27Z maronga
70! Bugfix: definition of character strings requires same length for all elements
71!
72! 1585 2015-04-30 07:05:52Z maronga
73! Modifications for RRTMG. Changed tables to PARAMETER type.
74!
75! 1571 2015-03-12 16:12:49Z maronga
76! Removed upper-case variable names. Corrected distribution of precipitation to
77! the liquid water reservoir and the bare soil fractions.
78!
79! 1555 2015-03-04 17:44:27Z maronga
80! Added output of r_a and r_s
81!
82! 1553 2015-03-03 17:33:54Z maronga
83! Improved better treatment of roughness lengths. Added default soil temperature
84! profile
85!
86! 1551 2015-03-03 14:18:16Z maronga
87! Flux calculation is now done in prandtl_fluxes. Added support for data output.
88! Vertical indices have been replaced. Restart runs are now possible. Some
89! variables have beem renamed. Bugfix in the prognostic equation for the surface
90! temperature. Introduced z0_eb and z0h_eb, which overwrite the setting of
91! roughness_length and z0_factor. Added Clapp & Hornberger parametrization for
92! the hydraulic conductivity. Bugfix for root fraction and extraction
93! calculation
94!
95! intrinsic function MAX and MIN
96!
97! 1500 2014-12-03 17:42:41Z maronga
98! Corrected calculation of aerodynamic resistance (r_a).
99! Precipitation is now added to liquid water reservoir using LE_liq.
100! Added support for dry runs.
101!
102! 1496 2014-12-02 17:25:50Z maronga
103! Initial revision
104!
105!
106! Description:
107! ------------
108!> Land surface model, consisting of a solver for the energy balance at the
109!> surface and a four layer soil scheme. The scheme is similar to the TESSEL
110!> scheme implemented in the ECMWF IFS model, with modifications according to
111!> H-TESSEL. The implementation is based on the formulation implemented in the
112!> DALES and UCLA-LES models.
113!>
114!> @todo Consider partial absorption of the net shortwave radiation by the
115!>       skin layer.
116!> @todo Improve surface water parameterization
117!> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0,
118!>       nzt_soil=3)).
119!> @todo Implement surface runoff model (required when performing long-term LES
120!>       with considerable precipitation.
121!> @todo Fix crashes with radiation_scheme == 'constant'
122!>
123!> @note No time step criterion is required as long as the soil layers do not
124!>       become too thin.
125!------------------------------------------------------------------------------!
126 MODULE land_surface_model_mod
127 
128    USE arrays_3d,                                                             &
129        ONLY:  hyp, ol, pt, pt_p, q, q_p, ql, qsws, shf, ts, us, vpt, z0, z0h, &
130               z0q
131
132    USE cloud_parameters,                                                      &
133        ONLY:  cp, hyrho, l_d_cp, l_d_r, l_v, prr, pt_d_t, rho_l, r_d, r_v
134
135    USE control_parameters,                                                    &
136        ONLY:  cloud_physics, dt_3d, humidity, intermediate_timestep_count,    &
137               initializing_actions, intermediate_timestep_count_max,          &
138               max_masks, precipitation, pt_surface, rho_surface,              &
139               roughness_length, surface_pressure, timestep_scheme, tsc,       &
140               z0h_factor, time_since_reference_point
141
142    USE indices,                                                               &
143        ONLY:  nbgp, nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner
144
145    USE kinds
146
147    USE pegrid
148
149    USE radiation_model_mod,                                                   &
150        ONLY:  force_radiation_call, radiation_scheme, rad_net, rad_sw_in,     &
151               rad_lw_out, rad_lw_out_change_0, sigma_sb,                      &
152               unscheduled_radiation_calls
153       
154#if defined ( __rrtmg )
155    USE radiation_model_mod,                                                   &
156        ONLY:  rrtm_idrv
157#endif
158
159    USE statistics,                                                            &
160        ONLY:  hom, statistic_regions
161
162    IMPLICIT NONE
163
164!
165!-- LSM model constants
166    INTEGER(iwp), PARAMETER :: nzb_soil = 0, & !< bottom of the soil model (to be switched)
167                               nzt_soil = 3, & !< top of the soil model (to be switched)
168                               nzs = 4         !< number of soil layers (fixed for now)
169
170    REAL(wp), PARAMETER ::                     &
171              b_ch               = 6.04_wp,    & ! Clapp & Hornberger exponent
172              lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil   
173              lambda_h_sm        = 3.44_wp,    & ! heat conductivity of the soil matrix
174              lambda_h_water     = 0.57_wp,    & ! heat conductivity of water
175              psi_sat            = -0.388_wp,  & ! soil matrix potential at saturation
176              rho_c_soil         = 2.19E6_wp,  & ! volumetric heat capacity of soil
177              rho_c_water        = 4.20E6_wp,  & ! volumetric heat capacity of water
178              m_max_depth        = 0.0002_wp     ! Maximum capacity of the water reservoir (m)
179
180
181!
182!-- LSM variables
183    INTEGER(iwp) :: veg_type  = 2, & !< NAMELIST veg_type_2d
184                    soil_type = 3    !< NAMELIST soil_type_2d
185
186    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  soil_type_2d, &  !< soil type, 0: user-defined, 1-7: generic (see list)
187                                                  veg_type_2d      !< vegetation type, 0: user-defined, 1-19: generic (see list)
188
189    LOGICAL, DIMENSION(:,:), ALLOCATABLE :: water_surface,     & !< flag parameter for water surfaces (classes 14+15)
190                                            pave_surface,      & !< flag parameter for pavements (asphalt etc.) (class 20)
191                                            building_surface     !< flag parameter indicating that the surface element is covered by buildings (no LSM actions, not implemented yet)
192
193    LOGICAL :: conserve_water_content = .TRUE.,  & !< open or closed bottom surface for the soil model
194               force_radiation_call_l = .FALSE., & !< flag parameter for unscheduled radiation model calls
195               land_surface = .FALSE.              !< flag parameter indicating wheather the lsm is used
196
197!   value 9999999.9_wp -> generic available or user-defined value must be set
198!   otherwise -> no generic variable and user setting is optional
199    REAL(wp) :: alpha_vangenuchten = 9999999.9_wp,      & !< NAMELIST alpha_vg
200                canopy_resistance_coefficient = 9999999.9_wp, & !< NAMELIST g_d
201                c_surface   = 20000.0_wp,               & !< Surface (skin) heat capacity
202                drho_l_lv,                              & !< (rho_l * l_v)**-1
203                exn,                                    & !< value of the Exner function
204                e_s = 0.0_wp,                           & !< saturation water vapour pressure
205                field_capacity = 9999999.9_wp,          & !< NAMELIST m_fc
206                f_shortwave_incoming = 9999999.9_wp,    & !< NAMELIST f_sw_in
207                hydraulic_conductivity = 9999999.9_wp,  & !< NAMELIST gamma_w_sat
208                ke = 0.0_wp,                            & !< Kersten number
209                lambda_h_sat = 0.0_wp,                  & !< heat conductivity for saturated soil
210                lambda_surface_stable = 9999999.9_wp,   & !< NAMELIST lambda_surface_s
211                lambda_surface_unstable = 9999999.9_wp, & !< NAMELIST lambda_surface_u
212                leaf_area_index = 9999999.9_wp,         & !< NAMELIST lai
213                l_vangenuchten = 9999999.9_wp,          & !< NAMELIST l_vg
214                min_canopy_resistance = 9999999.9_wp,   & !< NAMELIST r_canopy_min
215                min_soil_resistance = 50.0_wp,          & !< NAMELIST r_soil_min
216                m_total = 0.0_wp,                       & !< weighted total water content of the soil (m3/m3)
217                n_vangenuchten = 9999999.9_wp,          & !< NAMELIST n_vg
218                pave_depth = 9999999.9_wp,              & !< depth of the pavement
219                pave_heat_capacity = 1.94E6_wp,         & !< volumetric heat capacity of pavement (e.g. roads)
220                pave_heat_conductivity = 1.00_wp,       & !< heat conductivity for pavements (e.g. roads)
221                q_s = 0.0_wp,                           & !< saturation specific humidity
222                residual_moisture = 9999999.9_wp,       & !< NAMELIST m_res
223                rho_cp,                                 & !< rho_surface * cp
224                rho_lv,                                 & !< rho * l_v
225                rd_d_rv,                                & !< r_d / r_v
226                saturation_moisture = 9999999.9_wp,     & !< NAMELIST m_sat
227                skip_time_do_lsm = 0.0_wp,              & !< LSM is not called before this time
228                vegetation_coverage = 9999999.9_wp,     & !< NAMELIST c_veg
229                wilting_point = 9999999.9_wp,           & !< NAMELIST m_wilt
230                z0_eb  = 9999999.9_wp,                  & !< NAMELIST z0 (lsm_par)
231                z0h_eb = 9999999.9_wp,                  & !< NAMELIST z0h (lsm_par)
232                z0q_eb = 9999999.9_wp                     !< NAMELIST z0q (lsm_par)
233
234    REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: &
235              ddz_soil_stag,                  & !< 1/dz_soil_stag
236              dz_soil_stag,                   & !< soil grid spacing (center-center)
237              root_extr = 0.0_wp,             & !< root extraction
238              root_fraction = (/9999999.9_wp, 9999999.9_wp,    &
239                                9999999.9_wp, 9999999.9_wp /), & !< distribution of root surface area to the individual soil layers
240              zs = (/0.07_wp, 0.28_wp, 1.00_wp,  2.89_wp/),    & !< soil layer depths (m)
241              soil_moisture = 0.0_wp          !< soil moisture content (m3/m3)
242
243    REAL(wp), DIMENSION(nzb_soil:nzt_soil+1) ::   &
244              soil_temperature = (/290.0_wp, 287.0_wp, 285.0_wp,  283.0_wp,    & !< soil temperature (K)
245                                   283.0_wp /),                                &                                   
246              ddz_soil,                                                        & !< 1/dz_soil
247              dz_soil                                                            !< soil grid spacing (edge-edge)
248
249#if defined( __nopointer )
250    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface,   & !< surface temperature (K)
251                                                     t_surface_p, & !< progn. surface temperature (K)
252                                                     m_liq_eb,    & !< liquid water reservoir (m)
253                                                     m_liq_eb_av, & !< liquid water reservoir (m)
254                                                     m_liq_eb_p     !< progn. liquid water reservoir (m)
255#else
256    REAL(wp), DIMENSION(:,:), POINTER :: t_surface,      &
257                                         t_surface_p,    & 
258                                         m_liq_eb,       & 
259                                         m_liq_eb_p
260
261    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface_1, t_surface_2, &
262                                                     m_liq_eb_av,              &
263                                                     m_liq_eb_1, m_liq_eb_2
264#endif
265
266!
267!-- Temporal tendencies for time stepping           
268    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_surface_m,  & !< surface temperature tendency (K)
269                                             tm_liq_eb_m      !< liquid water reservoir tendency (m)
270
271!
272!-- Energy balance variables               
273    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
274              alpha_vg,         & !< coef. of Van Genuchten
275              c_liq,            & !< liquid water coverage (of vegetated area)
276              c_liq_av,         & !< average of c_liq
277              c_soil_av,        & !< average of c_soil
278              c_veg,            & !< vegetation coverage
279              c_veg_av,         & !< average of c_veg
280              f_sw_in,          & !< fraction of absorbed shortwave radiation by the surface layer (not implemented yet)
281              ghf_eb,           & !< ground heat flux
282              ghf_eb_av,        & !< average of ghf_eb
283              gamma_w_sat,      & !< hydraulic conductivity at saturation
284              g_d,              & !< coefficient for dependence of r_canopy on water vapour pressure deficit
285              lai,              & !< leaf area index
286              lai_av,           & !< average of lai
287              lambda_surface_s, & !< coupling between surface and soil (depends on vegetation type)
288              lambda_surface_u, & !< coupling between surface and soil (depends on vegetation type)
289              l_vg,             & !< coef. of Van Genuchten
290              m_fc,             & !< soil moisture at field capacity (m3/m3)
291              m_res,            & !< residual soil moisture
292              m_sat,            & !< saturation soil moisture (m3/m3)
293              m_wilt,           & !< soil moisture at permanent wilting point (m3/m3)
294              n_vg,             & !< coef. Van Genuchten 
295              qsws_eb,          & !< surface flux of latent heat (total)
296              qsws_eb_av,       & !< average of qsws_eb
297              qsws_liq_eb,      & !< surface flux of latent heat (liquid water portion)
298              qsws_liq_eb_av,   & !< average of qsws_liq_eb
299              qsws_soil_eb,     & !< surface flux of latent heat (soil portion)
300              qsws_soil_eb_av,  & !< average of qsws_soil_eb
301              qsws_veg_eb,      & !< surface flux of latent heat (vegetation portion)
302              qsws_veg_eb_av,   & !< average of qsws_veg_eb
303              rad_net_l,        & !< local copy of rad_net (net radiation at surface)
304              r_a,              & !< aerodynamic resistance
305              r_a_av,           & !< average of r_a
306              r_canopy,         & !< canopy resistance
307              r_soil,           & !< soil resistance
308              r_soil_min,       & !< minimum soil resistance
309              r_s,              & !< total surface resistance (combination of r_soil and r_canopy)
310              r_s_av,           & !< average of r_s
311              r_canopy_min,     & !< minimum canopy (stomatal) resistance
312              shf_eb,           & !< surface flux of sensible heat
313              shf_eb_av           !< average of shf_eb
314
315
316    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
317              lambda_h, &   !< heat conductivity of soil (W/m/K)                           
318              lambda_w, &   !< hydraulic diffusivity of soil (?)
319              gamma_w,  &   !< hydraulic conductivity of soil (W/m/K)
320              rho_c_total   !< volumetric heat capacity of the actual soil matrix (?)
321
322#if defined( __nopointer )
323    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
324              t_soil,    & !< Soil temperature (K)
325              t_soil_av, & !< Average of t_soil
326              t_soil_p,  & !< Prog. soil temperature (K)
327              m_soil,    & !< Soil moisture (m3/m3)
328              m_soil_av, & !< Average of m_soil
329              m_soil_p     !< Prog. soil moisture (m3/m3)
330#else
331    REAL(wp), DIMENSION(:,:,:), POINTER ::                                     &
332              t_soil, t_soil_p, &
333              m_soil, m_soil_p   
334
335    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
336              t_soil_av, t_soil_1, t_soil_2,                                   &
337              m_soil_av, m_soil_1, m_soil_2
338#endif
339
340
341    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
342              tt_soil_m, & !< t_soil storage array
343              tm_soil_m, & !< m_soil storage array
344              root_fr      !< root fraction (sum=1)
345
346
347!
348!-- Predefined Land surface classes (veg_type)
349    CHARACTER(26), DIMENSION(0:20), PARAMETER :: veg_type_name = (/ &
350                                   'user defined              ',    & ! 0
351                                   'crops, mixed farming      ',    & !  1
352                                   'short grass               ',    & !  2
353                                   'evergreen needleleaf trees',    & !  3
354                                   'deciduous needleleaf trees',    & !  4
355                                   'evergreen broadleaf trees ',    & !  5
356                                   'deciduous broadleaf trees ',    & !  6
357                                   'tall grass                ',    & !  7
358                                   'desert                    ',    & !  8
359                                   'tundra                    ',    & !  9
360                                   'irrigated crops           ',    & ! 10
361                                   'semidesert                ',    & ! 11
362                                   'ice caps and glaciers     ',    & ! 12
363                                   'bogs and marshes          ',    & ! 13
364                                   'inland water              ',    & ! 14
365                                   'ocean                     ',    & ! 15
366                                   'evergreen shrubs          ',    & ! 16
367                                   'deciduous shrubs          ',    & ! 17
368                                   'mixed forest/woodland     ',    & ! 18
369                                   'interrupted forest        ',    & ! 19
370                                   'pavements/roads           '     & ! 20
371                                                                 /)
372
373!
374!-- Soil model classes (soil_type)
375    CHARACTER(12), DIMENSION(0:7), PARAMETER :: soil_type_name = (/ &
376                                   'user defined',                  & ! 0
377                                   'coarse      ',                  & ! 1
378                                   'medium      ',                  & ! 2
379                                   'medium-fine ',                  & ! 3
380                                   'fine        ',                  & ! 4
381                                   'very fine   ',                  & ! 5
382                                   'organic     ',                  & ! 6
383                                   'loamy (CH)  '                   & ! 7
384                                                                 /)
385!
386!-- Land surface parameters according to the respective classes (veg_type)
387
388!
389!-- Land surface parameters I
390!--                          r_canopy_min,     lai,   c_veg,     g_d
391    REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: veg_pars = RESHAPE( (/ &
392                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & !  1
393                                 110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp,  & !  2
394                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  3
395                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  4
396                                 175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  5
397                                 240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp,  & !  6
398                                 100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp,  & !  7
399                                 250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp,  & !  8
400                                  80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp,  & !  9
401                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & ! 10
402                                 150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp,  & ! 11
403                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 12
404                                 240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp,  & ! 13
405                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 14
406                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 15
407                                 225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp,  & ! 16
408                                 225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp,  & ! 17
409                                 250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & ! 18
410                                 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp,  & ! 19
411                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp   & ! 20
412                                 /), (/ 4, 20 /) )
413
414!
415!-- Land surface parameters II          z0,         z0h,         z0q
416    REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: roughness_par = RESHAPE( (/ & 
417                                   0.25_wp,  0.25E-2_wp,  0.25E-2_wp,       & !  1
418                                   0.20_wp,  0.20E-2_wp,  0.20E-2_wp,       & !  2
419                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  3
420                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  4
421                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  5
422                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  6
423                                   0.47_wp,  0.47E-2_wp,  0.47E-2_wp,       & !  7
424                                  0.013_wp, 0.013E-2_wp, 0.013E-2_wp,       & !  8
425                                  0.034_wp, 0.034E-2_wp, 0.034E-2_wp,       & !  9
426                                    0.5_wp,  0.50E-2_wp,  0.50E-2_wp,       & ! 10
427                                   0.17_wp,  0.17E-2_wp,  0.17E-2_wp,       & ! 11
428                                 1.3E-3_wp,   1.3E-4_wp,   1.3E-4_wp,       & ! 12
429                                   0.83_wp,  0.83E-2_wp,  0.83E-2_wp,       & ! 13
430                                   0.00_wp,     0.00_wp,     0.00_wp,       & ! 14
431                                   0.00_wp,     0.00_wp,     0.00_wp,       & ! 15
432                                   0.10_wp,  0.10E-2_wp,  0.10E-2_wp,       & ! 16
433                                   0.25_wp,  0.25E-2_wp,  0.25E-2_wp,       & ! 17
434                                   2.00_wp,  2.00E-2_wp,  2.00E-2_wp,       & ! 18
435                                   1.10_wp,  1.10E-2_wp,  1.10E-2_wp,       & ! 19
436                                 1.0E-4_wp,   1.0E-5_wp,   1.0E-5_wp        & ! 20
437                                 /), (/ 3, 20 /) )
438
439!
440!-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in
441    REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: surface_pars = RESHAPE( (/ &
442                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  1
443                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  2
444                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  3
445                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  4
446                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  5
447                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  6
448                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  7
449                                      15.0_wp,       15.0_wp, 0.00_wp,     & !  8
450                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  9
451                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 10
452                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 11
453                                      58.0_wp,       58.0_wp, 0.00_wp,     & ! 12
454                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 13
455                                    1.0E10_wp,     1.0E10_wp, 0.00_wp,     & ! 14
456                                    1.0E10_wp,     1.0E10_wp, 0.00_wp,     & ! 15
457                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 16
458                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 17
459                                      20.0_wp,       15.0_wp, 0.03_wp,     & ! 18
460                                      20.0_wp,       15.0_wp, 0.03_wp,     & ! 19
461                                       0.0_wp,        0.0_wp, 0.00_wp      & ! 20
462                                      /), (/ 3, 20 /) )
463
464!
465!-- Root distribution (sum = 1)  level 1, level 2, level 3, level 4,
466    REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: root_distribution = RESHAPE( (/ &
467                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & !  1
468                                 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp,            & !  2
469                                 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp,            & !  3
470                                 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp,            & !  4
471                                 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp,            & !  5
472                                 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp,            & !  6
473                                 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp,            & !  7
474                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  8
475                                 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp,            & !  9
476                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & ! 10
477                                 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp,            & ! 11
478                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 12
479                                 0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp,            & ! 13
480                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 14
481                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 15
482                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 16
483                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 17
484                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 18
485                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 19
486                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp             & ! 20
487                                 /), (/ 4, 20 /) )
488
489!
490!-- Soil parameters according to the following porosity classes (soil_type)
491
492!
493!-- Soil parameters I           alpha_vg,      l_vg,    n_vg, gamma_w_sat
494    REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: soil_pars = RESHAPE( (/     &
495                                 3.83_wp,  1.250_wp, 1.38_wp,  6.94E-6_wp, & ! 1
496                                 3.14_wp, -2.342_wp, 1.28_wp,  1.16E-6_wp, & ! 2
497                                 0.83_wp, -0.588_wp, 1.25_wp,  0.26E-6_wp, & ! 3
498                                 3.67_wp, -1.977_wp, 1.10_wp,  2.87E-6_wp, & ! 4
499                                 2.65_wp,  2.500_wp, 1.10_wp,  1.74E-6_wp, & ! 5
500                                 1.30_wp,  0.400_wp, 1.20_wp,  0.93E-6_wp, & ! 6
501                                 0.00_wp,  0.00_wp,  0.00_wp,  0.57E-6_wp  & ! 7
502                                 /), (/ 4, 7 /) )
503
504!
505!-- Soil parameters II              m_sat,     m_fc,   m_wilt,    m_res 
506    REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: m_soil_pars = RESHAPE( (/ &
507                                 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & ! 1
508                                 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & ! 2
509                                 0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp, & ! 3
510                                 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, & ! 4
511                                 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, & ! 5
512                                 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp, & ! 6
513                                 0.472_wp, 0.323_wp, 0.171_wp, 0.000_wp  & ! 7
514                                 /), (/ 4, 7 /) )
515
516
517    SAVE
518
519
520    PRIVATE
521
522
523!
524!-- Public parameters, constants and initial values
525    PUBLIC alpha_vangenuchten, c_surface, canopy_resistance_coefficient,       &
526           conserve_water_content, field_capacity,                             &
527           f_shortwave_incoming, hydraulic_conductivity, init_lsm,             &
528           init_lsm_arrays, lambda_surface_stable, lambda_surface_unstable,    &
529           land_surface, leaf_area_index, lsm_energy_balance, lsm_soil_model,  &
530           lsm_swap_timelevel, l_vangenuchten, min_canopy_resistance,          &
531           min_soil_resistance, n_vangenuchten, pave_heat_capacity,            &
532           pave_depth, pave_heat_conductivity, residual_moisture, rho_cp,      &
533           rho_lv, root_fraction, saturation_moisture, skip_time_do_lsm,       &
534           soil_moisture, soil_temperature, soil_type, soil_type_name,         &
535           vegetation_coverage, veg_type, veg_type_name, wilting_point, z0_eb, &
536           z0h_eb, z0q_eb
537
538!
539!-- Public grid variables
540    PUBLIC nzb_soil, nzs, nzt_soil, zs
541
542!
543!-- Public 2D output variables
544    PUBLIC c_liq, c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb, ghf_eb_av,     &
545           lai, lai_av, qsws_eb, qsws_eb_av, qsws_liq_eb, qsws_liq_eb_av,      &
546           qsws_soil_eb, qsws_soil_eb_av, qsws_veg_eb, qsws_veg_eb_av,         &
547           r_a, r_a_av, r_s, r_s_av, shf_eb, shf_eb_av
548
549!
550!-- Public prognostic variables
551    PUBLIC m_liq_eb, m_liq_eb_av, m_soil, m_soil_av, t_soil, t_soil_av
552
553    INTERFACE init_lsm
554       MODULE PROCEDURE init_lsm
555    END INTERFACE init_lsm
556
557    INTERFACE lsm_energy_balance
558       MODULE PROCEDURE lsm_energy_balance
559    END INTERFACE lsm_energy_balance
560
561    INTERFACE lsm_soil_model
562       MODULE PROCEDURE lsm_soil_model
563    END INTERFACE lsm_soil_model
564
565    INTERFACE lsm_swap_timelevel
566       MODULE PROCEDURE lsm_swap_timelevel
567    END INTERFACE lsm_swap_timelevel
568
569 CONTAINS
570
571
572!------------------------------------------------------------------------------!
573! Description:
574! ------------
575!> Allocate land surface model arrays and define pointers
576!------------------------------------------------------------------------------!
577    SUBROUTINE init_lsm_arrays
578   
579
580       IMPLICIT NONE
581
582!
583!--    Allocate surface and soil temperature / humidity
584#if defined( __nopointer )
585       ALLOCATE ( m_liq_eb(nysg:nyng,nxlg:nxrg) )
586       ALLOCATE ( m_liq_eb_p(nysg:nyng,nxlg:nxrg) )
587       ALLOCATE ( m_soil(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
588       ALLOCATE ( m_soil_p(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
589       ALLOCATE ( t_surface(nysg:nyng,nxlg:nxrg) )
590       ALLOCATE ( t_surface_p(nysg:nyng,nxlg:nxrg) )
591       ALLOCATE ( t_soil(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
592       ALLOCATE ( t_soil_p(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
593#else
594       ALLOCATE ( m_liq_eb_1(nysg:nyng,nxlg:nxrg) )
595       ALLOCATE ( m_liq_eb_2(nysg:nyng,nxlg:nxrg) )
596       ALLOCATE ( m_soil_1(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
597       ALLOCATE ( m_soil_2(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
598       ALLOCATE ( t_surface_1(nysg:nyng,nxlg:nxrg) )
599       ALLOCATE ( t_surface_2(nysg:nyng,nxlg:nxrg) )
600       ALLOCATE ( t_soil_1(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
601       ALLOCATE ( t_soil_2(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
602#endif
603
604!
605!--    Allocate intermediate timestep arrays
606       ALLOCATE ( tm_liq_eb_m(nysg:nyng,nxlg:nxrg) )
607       ALLOCATE ( tm_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
608       ALLOCATE ( tt_surface_m(nysg:nyng,nxlg:nxrg) )
609       ALLOCATE ( tt_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
610
611!
612!--    Allocate 2D vegetation model arrays
613       ALLOCATE ( alpha_vg(nysg:nyng,nxlg:nxrg) )
614       ALLOCATE ( building_surface(nysg:nyng,nxlg:nxrg) )
615       ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) )
616       ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) )
617       ALLOCATE ( f_sw_in(nysg:nyng,nxlg:nxrg) )
618       ALLOCATE ( ghf_eb(nysg:nyng,nxlg:nxrg) )
619       ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) )
620       ALLOCATE ( g_d(nysg:nyng,nxlg:nxrg) )
621       ALLOCATE ( lai(nysg:nyng,nxlg:nxrg) )
622       ALLOCATE ( l_vg(nysg:nyng,nxlg:nxrg) )
623       ALLOCATE ( lambda_surface_u(nysg:nyng,nxlg:nxrg) )
624       ALLOCATE ( lambda_surface_s(nysg:nyng,nxlg:nxrg) )
625       ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) )
626       ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) )
627       ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) )
628       ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) )
629       ALLOCATE ( n_vg(nysg:nyng,nxlg:nxrg) )
630       ALLOCATE ( pave_surface(nysg:nyng,nxlg:nxrg) )
631       ALLOCATE ( qsws_eb(nysg:nyng,nxlg:nxrg) )
632       ALLOCATE ( qsws_soil_eb(nysg:nyng,nxlg:nxrg) )
633       ALLOCATE ( qsws_liq_eb(nysg:nyng,nxlg:nxrg) )
634       ALLOCATE ( qsws_veg_eb(nysg:nyng,nxlg:nxrg) )
635       ALLOCATE ( rad_net_l(nysg:nyng,nxlg:nxrg) )
636       ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) )
637       ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) )
638       ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) )
639       ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) )
640       ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) )
641       ALLOCATE ( r_canopy_min(nysg:nyng,nxlg:nxrg) )
642       ALLOCATE ( shf_eb(nysg:nyng,nxlg:nxrg) )
643       ALLOCATE ( soil_type_2d(nysg:nyng,nxlg:nxrg) )
644       ALLOCATE ( veg_type_2d(nysg:nyng,nxlg:nxrg) )
645       ALLOCATE ( water_surface(nysg:nyng,nxlg:nxrg) )
646
647#if ! defined( __nopointer )
648!
649!--    Initial assignment of the pointers
650       t_soil    => t_soil_1;    t_soil_p    => t_soil_2
651       t_surface => t_surface_1; t_surface_p => t_surface_2
652       m_soil    => m_soil_1;    m_soil_p    => m_soil_2
653       m_liq_eb  => m_liq_eb_1;  m_liq_eb_p  => m_liq_eb_2
654#endif
655
656
657    END SUBROUTINE init_lsm_arrays
658
659!------------------------------------------------------------------------------!
660! Description:
661! ------------
662!> Initialization of the land surface model
663!------------------------------------------------------------------------------!
664    SUBROUTINE init_lsm
665   
666
667       IMPLICIT NONE
668
669       INTEGER(iwp) ::  i !< running index
670       INTEGER(iwp) ::  j !< running index
671       INTEGER(iwp) ::  k !< running index
672
673       REAL(wp) :: pt1   !< potential temperature at first grid level
674
675
676!
677!--    Calculate Exner function
678       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
679
680
681!
682!--    If no cloud physics is used, rho_surface has not been calculated before
683       IF (  .NOT.  cloud_physics )  THEN
684          rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn )
685       ENDIF
686
687!
688!--    Calculate frequently used parameters
689       rho_cp    = cp * rho_surface
690       rd_d_rv   = r_d / r_v
691       rho_lv    = rho_surface * l_v
692       drho_l_lv = 1.0_wp / (rho_l * l_v)
693
694!
695!--    Set inital values for prognostic quantities
696       tt_surface_m = 0.0_wp
697       tt_soil_m    = 0.0_wp
698       tm_soil_m    = 0.0_wp
699       tm_liq_eb_m  = 0.0_wp
700       c_liq        = 0.0_wp
701
702       ghf_eb = 0.0_wp
703       shf_eb = rho_cp * shf
704
705       IF ( humidity )  THEN
706          qsws_eb = rho_l * l_v * qsws
707       ELSE
708          qsws_eb = 0.0_wp
709       ENDIF
710
711       qsws_liq_eb  = 0.0_wp
712       qsws_soil_eb = 0.0_wp
713       qsws_veg_eb  = 0.0_wp
714
715       r_a        = 50.0_wp
716       r_s        = 50.0_wp
717       r_canopy   = 0.0_wp
718       r_soil     = 0.0_wp
719
720!
721!--    Allocate 3D soil model arrays
722       ALLOCATE ( root_fr(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
723       ALLOCATE ( lambda_h(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
724       ALLOCATE ( rho_c_total(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
725
726       lambda_h = 0.0_wp
727!
728!--    If required, allocate humidity-related variables for the soil model
729       IF ( humidity )  THEN
730          ALLOCATE ( lambda_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
731          ALLOCATE ( gamma_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )   
732
733          lambda_w = 0.0_wp 
734       ENDIF
735
736!
737!--    Calculate grid spacings. Temperature and moisture are defined at
738!--    the edges of the soil layers (_stag), whereas gradients/fluxes are defined
739!--    at the centers
740       dz_soil(nzb_soil) = zs(nzb_soil)
741
742       DO  k = nzb_soil+1, nzt_soil
743          dz_soil(k) = zs(k) - zs(k-1)
744       ENDDO
745       dz_soil(nzt_soil+1) = dz_soil(nzt_soil)
746
747       DO  k = nzb_soil, nzt_soil-1
748          dz_soil_stag(k) = 0.5_wp * (dz_soil(k+1) + dz_soil(k))
749       ENDDO
750       dz_soil_stag(nzt_soil) = dz_soil(nzt_soil)
751
752       ddz_soil      = 1.0_wp / dz_soil
753       ddz_soil_stag = 1.0_wp / dz_soil_stag
754
755!
756!--    Initialize standard soil types. It is possible to overwrite each
757!--    parameter by setting the respecticy NAMELIST variable to a
758!--    value /= 9999999.9.
759       IF ( soil_type /= 0 )  THEN 
760 
761          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
762             alpha_vangenuchten = soil_pars(0,soil_type)
763          ENDIF
764
765          IF ( l_vangenuchten == 9999999.9_wp )  THEN
766             l_vangenuchten = soil_pars(1,soil_type)
767          ENDIF
768
769          IF ( n_vangenuchten == 9999999.9_wp )  THEN
770             n_vangenuchten = soil_pars(2,soil_type)           
771          ENDIF
772
773          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
774             hydraulic_conductivity = soil_pars(3,soil_type)           
775          ENDIF
776
777          IF ( saturation_moisture == 9999999.9_wp )  THEN
778             saturation_moisture = m_soil_pars(0,soil_type)           
779          ENDIF
780
781          IF ( field_capacity == 9999999.9_wp )  THEN
782             field_capacity = m_soil_pars(1,soil_type)           
783          ENDIF
784
785          IF ( wilting_point == 9999999.9_wp )  THEN
786             wilting_point = m_soil_pars(2,soil_type)           
787          ENDIF
788
789          IF ( residual_moisture == 9999999.9_wp )  THEN
790             residual_moisture = m_soil_pars(3,soil_type)       
791          ENDIF
792
793       ENDIF   
794
795!
796!--    Map values to the respective 2D arrays
797       alpha_vg      = alpha_vangenuchten
798       l_vg          = l_vangenuchten
799       n_vg          = n_vangenuchten
800       gamma_w_sat   = hydraulic_conductivity
801       m_sat         = saturation_moisture
802       m_fc          = field_capacity
803       m_wilt        = wilting_point
804       m_res         = residual_moisture
805       r_soil_min    = min_soil_resistance
806
807!
808!--    Initial run actions
809       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
810
811          t_soil    = 0.0_wp
812          m_liq_eb  = 0.0_wp
813          m_soil    = 0.0_wp
814
815!
816!--       Map user settings of T and q for each soil layer
817!--       (make sure that the soil moisture does not drop below the permanent
818!--       wilting point) -> problems with devision by zero)
819          DO  k = nzb_soil, nzt_soil
820             t_soil(k,:,:)    = soil_temperature(k)
821             m_soil(k,:,:)    = MAX(soil_moisture(k),m_wilt(:,:))
822             soil_moisture(k) = MAX(soil_moisture(k),wilting_point)
823          ENDDO
824          t_soil(nzt_soil+1,:,:) = soil_temperature(nzt_soil+1)
825
826!
827!--       Calculate surface temperature
828          t_surface   = pt_surface * exn
829
830!
831!--       Set artifical values for ts and us so that r_a has its initial value
832!--       for the first time step
833          DO  i = nxlg, nxrg
834             DO  j = nysg, nyng
835                k = nzb_s_inner(j,i)
836
837                IF ( cloud_physics )  THEN
838                   pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
839                ELSE
840                   pt1 = pt(k+1,j,i)
841                ENDIF
842
843!
844!--             Assure that r_a cannot be zero at model start
845                IF ( pt1 == pt(k,j,i) )  pt1 = pt1 + 1.0E-10_wp
846
847                us(j,i)  = 0.1_wp
848                ts(j,i)  = (pt1 - pt(k,j,i)) / r_a(j,i)
849                shf(j,i) = - us(j,i) * ts(j,i)
850             ENDDO
851          ENDDO
852
853!
854!--    Actions for restart runs
855       ELSE
856
857          DO  i = nxlg, nxrg
858             DO  j = nysg, nyng
859                k = nzb_s_inner(j,i)               
860                t_surface(j,i) = pt(k,j,i) * exn
861             ENDDO
862          ENDDO
863
864       ENDIF
865
866       DO  k = nzb_soil, nzt_soil
867          root_fr(k,:,:) = root_fraction(k)
868       ENDDO
869
870       IF ( veg_type /= 0 )  THEN
871          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
872             min_canopy_resistance = veg_pars(0,veg_type)
873          ENDIF
874          IF ( leaf_area_index == 9999999.9_wp )  THEN
875             leaf_area_index = veg_pars(1,veg_type)         
876          ENDIF
877          IF ( vegetation_coverage == 9999999.9_wp )  THEN
878             vegetation_coverage = veg_pars(2,veg_type)     
879          ENDIF
880          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
881              canopy_resistance_coefficient= veg_pars(3,veg_type)     
882          ENDIF
883          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
884             lambda_surface_stable = surface_pars(0,veg_type)         
885          ENDIF
886          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
887             lambda_surface_unstable = surface_pars(1,veg_type)       
888          ENDIF
889          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
890             f_shortwave_incoming = surface_pars(2,veg_type)       
891          ENDIF
892          IF ( z0_eb == 9999999.9_wp )  THEN
893             roughness_length = roughness_par(0,veg_type) 
894             z0_eb            = roughness_par(0,veg_type) 
895          ENDIF
896          IF ( z0h_eb == 9999999.9_wp )  THEN
897             z0h_eb = roughness_par(1,veg_type)
898          ENDIF
899          IF ( z0q_eb == 9999999.9_wp )  THEN
900             z0q_eb = roughness_par(2,veg_type)
901          ENDIF
902          z0h_factor = z0h_eb / ( z0_eb + 1.0E-20_wp )
903
904          IF ( ANY( root_fraction == 9999999.9_wp ) )  THEN
905             DO  k = nzb_soil, nzt_soil
906                root_fr(k,:,:) = root_distribution(k,veg_type)
907                root_fraction(k) = root_distribution(k,veg_type)
908             ENDDO
909          ENDIF
910
911       ELSE
912
913          IF ( z0_eb == 9999999.9_wp )  THEN
914             z0_eb = roughness_length
915          ENDIF
916          IF ( z0h_eb == 9999999.9_wp )  THEN
917             z0h_eb = z0_eb * z0h_factor
918          ENDIF
919          IF ( z0q_eb == 9999999.9_wp )  THEN
920             z0q_eb = z0_eb * z0h_factor
921          ENDIF
922
923       ENDIF
924
925!
926!--    For surfaces covered with pavement, set depth of the pavement (with dry
927!--    soil below). The depth must be greater than the first soil layer depth
928       IF ( veg_type == 20 )  THEN
929          IF ( pave_depth == 9999999.9_wp )  THEN
930             pave_depth = zs(nzb_soil) 
931          ELSE
932             pave_depth = MAX( zs(nzb_soil), pave_depth )
933          ENDIF
934       ENDIF
935
936!
937!--    Map vegetation and soil types to 2D array to allow for heterogeneous
938!--    surfaces via user interface see below
939       veg_type_2d = veg_type
940       soil_type_2d = soil_type
941
942!
943!--    Map vegetation parameters to the respective 2D arrays
944       r_canopy_min         = min_canopy_resistance
945       lai                  = leaf_area_index
946       c_veg                = vegetation_coverage
947       g_d                  = canopy_resistance_coefficient
948       lambda_surface_s     = lambda_surface_stable
949       lambda_surface_u     = lambda_surface_unstable
950       f_sw_in              = f_shortwave_incoming
951       z0                   = z0_eb
952       z0h                  = z0h_eb
953       z0q                  = z0q_eb
954
955!
956!--    Possibly do user-defined actions (e.g. define heterogeneous land surface)
957       CALL user_init_land_surface
958
959!
960!--    Set flag parameter if vegetation type was set to a water surface. Also
961!--    set temperature to a constant value in all "soil" layers.
962       DO  i = nxlg, nxrg
963          DO  j = nysg, nyng
964             IF ( veg_type_2d(j,i) == 14  .OR.  veg_type_2d(j,i) == 15 )  THEN
965                water_surface(j,i) = .TRUE.
966             ELSEIF ( veg_type_2d(j,i) == 20 )  THEN
967                pave_surface(j,i) = .TRUE.
968                m_soil(:,j,i) = 0.0_wp
969             ENDIF
970
971          ENDDO
972       ENDDO
973
974!
975!--    Calculate new roughness lengths (for water surfaces only)
976       CALL calc_z0_water_surface
977
978       t_soil_p    = t_soil
979       m_soil_p    = m_soil
980       m_liq_eb_p  = m_liq_eb
981       t_surface_p = t_surface
982
983
984
985!--    Store initial profiles of t_soil and m_soil (assuming they are
986!--    horizontally homogeneous on this PE)
987       hom(nzb_soil:nzt_soil,1,90,:)  = SPREAD( t_soil(nzb_soil:nzt_soil,      &
988                                                nysg,nxlg), 2,                 &
989                                                statistic_regions+1 )
990       hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil(nzb_soil:nzt_soil,      &
991                                                nysg,nxlg), 2,                 &
992                                                statistic_regions+1 )
993
994    END SUBROUTINE init_lsm
995
996
997
998!------------------------------------------------------------------------------!
999! Description:
1000! ------------
1001!> Solver for the energy balance at the surface.
1002!------------------------------------------------------------------------------!
1003    SUBROUTINE lsm_energy_balance
1004
1005
1006       IMPLICIT NONE
1007
1008       INTEGER(iwp) ::  i         !< running index
1009       INTEGER(iwp) ::  j         !< running index
1010       INTEGER(iwp) ::  k, ks     !< running index
1011
1012       REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface
1013                   f1,          & !< resistance correction term 1
1014                   f2,          & !< resistance correction term 2
1015                   f3,          & !< resistance correction term 3
1016                   m_min,       & !< minimum soil moisture
1017                   e,           & !< water vapour pressure
1018                   e_s,         & !< water vapour saturation pressure
1019                   e_s_dt,      & !< derivate of e_s with respect to T
1020                   tend,        & !< tendency
1021                   dq_s_dt,     & !< derivate of q_s with respect to T
1022                   coef_1,      & !< coef. for prognostic equation
1023                   coef_2,      & !< coef. for prognostic equation
1024                   f_qsws,      & !< factor for qsws_eb
1025                   f_qsws_veg,  & !< factor for qsws_veg_eb
1026                   f_qsws_soil, & !< factor for qsws_soil_eb
1027                   f_qsws_liq,  & !< factor for qsws_liq_eb
1028                   f_shf,       & !< factor for shf_eb
1029                   lambda_surface, & !< Current value of lambda_surface
1030                   m_liq_eb_max,   & !< maxmimum value of the liq. water reservoir
1031                   pt1,         & !< potential temperature at first grid level
1032                   qv1            !< specific humidity at first grid level
1033
1034!
1035!--    Calculate the exner function for the current time step
1036       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
1037
1038       DO  i = nxlg, nxrg
1039          DO  j = nysg, nyng
1040             k = nzb_s_inner(j,i)
1041
1042!
1043!--          Set lambda_surface according to stratification between skin layer and soil
1044             IF (  .NOT.  pave_surface(j,i) )  THEN
1045
1046                c_surface_tmp = c_surface
1047
1048                IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i))  THEN
1049                   lambda_surface = lambda_surface_s(j,i)
1050                ELSE
1051                   lambda_surface = lambda_surface_u(j,i)
1052                ENDIF
1053             ELSE
1054
1055                c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp
1056                lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil)
1057
1058             ENDIF
1059
1060!
1061!--          First step: calculate aerodyamic resistance. As pt, us, ts
1062!--          are not available for the prognostic time step, data from the last
1063!--          time step is used here. Note that this formulation is the
1064!--          equivalent to the ECMWF formulation using drag coefficients
1065             IF ( cloud_physics )  THEN
1066                pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
1067                qv1 = q(k+1,j,i) - ql(k+1,j,i)
1068             ELSE
1069                pt1 = pt(k+1,j,i)
1070                qv1 = q(k+1,j,i)
1071             ENDIF
1072
1073             r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp)
1074
1075!
1076!--          Make sure that the resistance does not drop to zero
1077             IF ( ABS(r_a(j,i)) < 1.0E-10_wp )  r_a(j,i) = 1.0E-10_wp
1078
1079!
1080!--          Second step: calculate canopy resistance r_canopy
1081!--          f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
1082 
1083!--          f1: correction for incoming shortwave radiation (stomata close at
1084!--          night)
1085             IF ( radiation_scheme /= 'constant' )  THEN
1086                f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /  &
1087                              (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)          &
1088                               + 1.0_wp)) )
1089             ELSE
1090                f1 = 1.0_wp
1091             ENDIF
1092
1093
1094!
1095!--          f2: correction for soil moisture availability to plants (the
1096!--          integrated soil moisture must thus be considered here)
1097!--          f2 = 0 for very dry soils
1098             m_total = 0.0_wp
1099             DO  ks = nzb_soil, nzt_soil
1100                 m_total = m_total + root_fr(ks,j,i)                           &
1101                           * MAX(m_soil(ks,j,i),m_wilt(j,i))
1102             ENDDO
1103
1104             IF ( m_total > m_wilt(j,i)  .AND.  m_total < m_fc(j,i) )  THEN
1105                f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
1106             ELSEIF ( m_total >= m_fc(j,i) )  THEN
1107                f2 = 1.0_wp
1108             ELSE
1109                f2 = 1.0E-20_wp
1110             ENDIF
1111
1112!
1113!--          Calculate water vapour pressure at saturation
1114             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)     &
1115                           - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) )
1116
1117!
1118!--          f3: correction for vapour pressure deficit
1119             IF ( g_d(j,i) /= 0.0_wp )  THEN
1120!
1121!--             Calculate vapour pressure
1122                e  = qv1 * surface_pressure / 0.622_wp
1123                f3 = EXP ( -g_d(j,i) * (e_s - e) )
1124             ELSE
1125                f3 = 1.0_wp
1126             ENDIF
1127
1128!
1129!--          Calculate canopy resistance. In case that c_veg is 0 (bare soils),
1130!--          this calculation is obsolete, as r_canopy is not used below.
1131!--          To do: check for very dry soil -> r_canopy goes to infinity
1132             r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3      &
1133                                             + 1.0E-20_wp)
1134
1135!
1136!--          Third step: calculate bare soil resistance r_soil. The Clapp &
1137!--          Hornberger parametrization does not consider c_veg.
1138             IF ( soil_type_2d(j,i) /= 7 )  THEN
1139                m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *     &
1140                        m_res(j,i)
1141             ELSE
1142                m_min = m_wilt(j,i)
1143             ENDIF
1144
1145             f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min )
1146             f2 = MAX(f2,1.0E-20_wp)
1147             f2 = MIN(f2,1.0_wp)
1148
1149             r_soil(j,i) = r_soil_min(j,i) / f2
1150
1151!
1152!--          Calculate the maximum possible liquid water amount on plants and
1153!--          bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is
1154!--          assumed, while paved surfaces might hold up 1 mm of water. The
1155!--          liquid water fraction for paved surfaces is calculated after
1156!--          Noilhan & Planton (1989), while the ECMWF formulation is used for
1157!--          vegetated surfaces and bare soils.
1158             IF ( pave_surface(j,i) )  THEN
1159                m_liq_eb_max = m_max_depth * 5.0_wp
1160                c_liq(j,i) = MIN( 1.0_wp, (m_liq_eb(j,i) / m_liq_eb_max)**0.67 )
1161             ELSE
1162                m_liq_eb_max = m_max_depth * ( c_veg(j,i) * lai(j,i)           &
1163                               + (1.0_wp - c_veg(j,i)) )
1164                c_liq(j,i) = MIN( 1.0_wp, m_liq_eb(j,i) / m_liq_eb_max )
1165             ENDIF
1166
1167!
1168!--          Calculate saturation specific humidity
1169             q_s = 0.622_wp * e_s / surface_pressure
1170
1171!
1172!--          In case of dewfall, set evapotranspiration to zero
1173!--          All super-saturated water is then removed from the air
1174             IF ( humidity  .AND.  q_s <= qv1 )  THEN
1175                r_canopy(j,i) = 0.0_wp
1176                r_soil(j,i)   = 0.0_wp
1177             ENDIF
1178
1179!
1180!--          Calculate coefficients for the total evapotranspiration
1181!--          In case of water surface, set vegetation and soil fluxes to zero.
1182!--          For pavements, only evaporation of liquid water is possible.
1183             IF ( water_surface(j,i) )  THEN
1184                f_qsws_veg  = 0.0_wp
1185                f_qsws_soil = 0.0_wp
1186                f_qsws_liq  = rho_lv / r_a(j,i)
1187             ELSEIF ( pave_surface (j,i) )  THEN
1188                f_qsws_veg  = 0.0_wp
1189                f_qsws_soil = 0.0_wp
1190                f_qsws_liq  = rho_lv * c_liq(j,i) / r_a(j,i)
1191             ELSE
1192                f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/     &
1193                              (r_a(j,i) + r_canopy(j,i))
1194                f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +     &
1195                                                          r_soil(j,i))
1196                f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
1197             ENDIF
1198!
1199!--          If soil moisture is below wilting point, plants do no longer
1200!--          transpirate.
1201!              IF ( m_soil(k,j,i) < m_wilt(j,i) )  THEN
1202!                 f_qsws_veg = 0.0_wp
1203!              ENDIF
1204
1205             f_shf  = rho_cp / r_a(j,i)
1206             f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
1207
1208!
1209!--          Calculate derivative of q_s for Taylor series expansion
1210             e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) -        &
1211                              17.269_wp*(t_surface(j,i) - 273.16_wp)           &
1212                              / (t_surface(j,i) - 35.86_wp)**2 )
1213
1214             dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
1215
1216!
1217!--          Add LW up so that it can be removed in prognostic equation
1218             rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i)
1219
1220!
1221!--          Calculate new skin temperature
1222             IF ( humidity )  THEN
1223#if defined ( __rrtmg )
1224!
1225!--             Numerator of the prognostic equation
1226                coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
1227                         * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
1228                         + f_shf * pt1 + f_qsws * ( qv1 - q_s                  &
1229                         + dq_s_dt * t_surface(j,i) ) + lambda_surface         &
1230                         * t_soil(nzb_soil,j,i)
1231
1232!
1233!--             Denominator of the prognostic equation
1234                coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt           &
1235                         + lambda_surface + f_shf / exn
1236#else
1237
1238!
1239!--             Numerator of the prognostic equation
1240                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
1241                         * t_surface(j,i) ** 4                                 &
1242                         + f_shf * pt1 + f_qsws * ( qv1                        &
1243                         - q_s + dq_s_dt * t_surface(j,i) )                    &
1244                         + lambda_surface * t_soil(nzb_soil,j,i)
1245
1246!
1247!--             Denominator of the prognostic equation
1248                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws      &
1249                         * dq_s_dt + lambda_surface + f_shf / exn
1250 
1251#endif
1252             ELSE
1253
1254#if defined ( __rrtmg )
1255!
1256!--             Numerator of the prognostic equation
1257                coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
1258                         * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
1259                         + f_shf * pt1  + lambda_surface                       &
1260                         * t_soil(nzb_soil,j,i)
1261
1262!
1263!--             Denominator of the prognostic equation
1264                coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
1265#else
1266
1267!
1268!--             Numerator of the prognostic equation
1269                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
1270                          * t_surface(j,i) ** 4 + f_shf * pt1                  &
1271                          + lambda_surface * t_soil(nzb_soil,j,i)
1272
1273!
1274!--             Denominator of the prognostic equation
1275                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3               &
1276                         + lambda_surface + f_shf / exn
1277#endif
1278             ENDIF
1279
1280             tend = 0.0_wp
1281
1282!
1283!--          Implicit solution when the surface layer has no heat capacity,
1284!--          otherwise use RK3 scheme.
1285             t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *    &
1286                                t_surface(j,i) ) / ( c_surface_tmp + coef_2    &
1287                                * dt_3d * tsc(2) ) 
1288
1289!
1290!--          Add RK3 term
1291             IF ( c_surface_tmp /= 0.0_wp )  THEN
1292
1293                t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)           &
1294                                   * tt_surface_m(j,i)
1295
1296!
1297!--             Calculate true tendency
1298                tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3)     &
1299                       * tt_surface_m(j,i)) / (dt_3d  * tsc(2))
1300!
1301!--             Calculate t_surface tendencies for the next Runge-Kutta step
1302                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1303                   IF ( intermediate_timestep_count == 1 )  THEN
1304                      tt_surface_m(j,i) = tend
1305                   ELSEIF ( intermediate_timestep_count <                      &
1306                            intermediate_timestep_count_max )  THEN
1307                      tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp        &
1308                                          * tt_surface_m(j,i)
1309                   ENDIF
1310                ENDIF
1311             ENDIF
1312
1313!
1314!--          In case of fast changes in the skin temperature, it is possible to
1315!--          update the radiative fluxes independently from the prescribed
1316!--          radiation call frequency. This effectively prevents oscillations,
1317!--          especially when setting skip_time_do_radiation /= 0. The threshold
1318!--          value of 0.2 used here is just a first guess. This method should be
1319!--          revised in the future as tests have shown that the threshold is
1320!--          often reached, when no oscillations would occur (causes immense
1321!--          computing time for the radiation code).
1322             IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp  .AND.     &
1323                  unscheduled_radiation_calls )  THEN
1324                force_radiation_call_l = .TRUE.
1325             ENDIF
1326
1327             pt(k,j,i) = t_surface_p(j,i) / exn
1328
1329!
1330!--          Calculate fluxes
1331#if defined ( __rrtmg )
1332             rad_net_l(j,i)   = rad_net_l(j,i) + rad_lw_out_change_0(j,i)      &
1333                                * t_surface(j,i) - rad_lw_out(nzb,j,i)         &
1334                                - rad_lw_out_change_0(j,i) * t_surface_p(j,i)
1335
1336             IF ( rrtm_idrv == 1 )  THEN
1337                rad_net(j,i) = rad_net_l(j,i)
1338                rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i)                      &
1339                                      + rad_lw_out_change_0(j,i)               &
1340                                      * ( t_surface_p(j,i) - t_surface(j,i) )
1341             ENDIF
1342#else
1343             rad_net_l(j,i)   = rad_net_l(j,i) + 3.0_wp * sigma_sb             &
1344                                * t_surface(j,i)**4 - 4.0_wp * sigma_sb        &
1345                                * t_surface(j,i)**3 * t_surface_p(j,i)
1346#endif
1347
1348             ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)               &
1349                              - t_soil(nzb_soil,j,i))
1350
1351             shf_eb(j,i)    = - f_shf * ( pt1 - pt(k,j,i) )
1352
1353             shf(j,i) = shf_eb(j,i) / rho_cp
1354
1355             IF ( humidity )  THEN
1356                qsws_eb(j,i)  = - f_qsws    * ( qv1 - q_s + dq_s_dt            &
1357                                * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
1358
1359                qsws(j,i) = qsws_eb(j,i) / rho_lv
1360
1361                qsws_veg_eb(j,i)  = - f_qsws_veg  * ( qv1 - q_s                &
1362                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1363                                    * t_surface_p(j,i) )
1364
1365                qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s                &
1366                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1367                                    * t_surface_p(j,i) )
1368
1369                qsws_liq_eb(j,i)  = - f_qsws_liq  * ( qv1 - q_s                &
1370                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1371                                    * t_surface_p(j,i) )
1372             ENDIF
1373
1374!
1375!--          Calculate the true surface resistance
1376             IF ( qsws_eb(j,i) == 0.0_wp )  THEN
1377                r_s(j,i) = 1.0E10_wp
1378             ELSE
1379                r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt                    &
1380                           * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )     &
1381                           / qsws_eb(j,i) - r_a(j,i)
1382             ENDIF
1383
1384!
1385!--          Calculate change in liquid water reservoir due to dew fall or
1386!--          evaporation of liquid water
1387             IF ( humidity )  THEN
1388!
1389!--             If precipitation is activated, add rain water to qsws_liq_eb
1390!--             and qsws_soil_eb according the the vegetation coverage.
1391!--             precipitation_rate is given in mm.
1392                IF ( precipitation )  THEN
1393
1394!
1395!--                Add precipitation to liquid water reservoir, if possible.
1396!--                Otherwise, add the water to soil. In case of
1397!--                pavements, the exceeding water amount is implicitely removed
1398!--                as runoff as qsws_soil_eb is then not used in the soil model
1399                   IF ( m_liq_eb(j,i) /= m_liq_eb_max )  THEN
1400                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i)                      &
1401                                       + c_veg(j,i) * prr(k,j,i) * hyrho(k)    &
1402                                       * 0.001_wp * rho_l * l_v
1403                   ELSE
1404                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
1405                                        + c_veg(j,i) * prr(k,j,i) * hyrho(k)   &
1406                                        * 0.001_wp * rho_l * l_v
1407                   ENDIF
1408
1409!--                Add precipitation to bare soil according to the bare soil
1410!--                coverage.
1411                   qsws_soil_eb(j,i) = qsws_soil_eb(j,i) * (1.0_wp             &
1412                                       - c_veg(j,i)) * prr(k,j,i) * hyrho(k)   &
1413                                       * 0.001_wp * rho_l * l_v
1414                ENDIF
1415
1416!
1417!--             If the air is saturated, check the reservoir water level
1418                IF ( qsws_eb(j,i) < 0.0_wp )  THEN
1419
1420!
1421!--                Check if reservoir is full (avoid values > m_liq_eb_max)
1422!--                In that case, qsws_liq_eb goes to qsws_soil_eb. In this
1423!--                case qsws_veg_eb is zero anyway (because c_liq = 1),       
1424!--                so that tend is zero and no further check is needed
1425                   IF ( m_liq_eb(j,i) == m_liq_eb_max )  THEN
1426                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
1427                                           + qsws_liq_eb(j,i)
1428                      qsws_liq_eb(j,i)  = 0.0_wp
1429                   ENDIF
1430
1431!
1432!--                In case qsws_veg_eb becomes negative (unphysical behavior),
1433!--                let the water enter the liquid water reservoir as dew on the
1434!--                plant
1435                   IF ( qsws_veg_eb(j,i) < 0.0_wp )  THEN
1436                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i)
1437                      qsws_veg_eb(j,i) = 0.0_wp
1438                   ENDIF
1439                ENDIF                   
1440 
1441                tend = - qsws_liq_eb(j,i) * drho_l_lv
1442
1443                m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend      &
1444                                                   + tsc(3) * tm_liq_eb_m(j,i) )
1445
1446!
1447!--             Check if reservoir is overfull -> reduce to maximum
1448!--             (conservation of water is violated here)
1449                m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max)
1450
1451!
1452!--             Check if reservoir is empty (avoid values < 0.0)
1453!--             (conservation of water is violated here)
1454                m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp)
1455
1456
1457!
1458!--             Calculate m_liq_eb tendencies for the next Runge-Kutta step
1459                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1460                   IF ( intermediate_timestep_count == 1 )  THEN
1461                      tm_liq_eb_m(j,i) = tend
1462                   ELSEIF ( intermediate_timestep_count <                      &
1463                            intermediate_timestep_count_max )  THEN
1464                      tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp         &
1465                                      * tm_liq_eb_m(j,i)
1466                   ENDIF
1467                ENDIF
1468
1469             ENDIF
1470
1471          ENDDO
1472       ENDDO
1473
1474!
1475!--    Make a logical OR for all processes. Force radiation call if at
1476!--    least one processor reached the threshold change in skin temperature
1477       IF ( unscheduled_radiation_calls  .AND.  intermediate_timestep_count    &
1478            == intermediate_timestep_count_max-1 )  THEN
1479#if defined( __parallel )
1480          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1481          CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,    &
1482                              1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
1483#else
1484          force_radiation_call = force_radiation_call_l
1485#endif
1486          force_radiation_call_l = .FALSE.
1487       ENDIF
1488
1489!
1490!--    Calculate surface specific humidity
1491       IF ( humidity )  THEN
1492          CALL calc_q_surface
1493       ENDIF
1494
1495!
1496!--    Calculate new roughness lengths (for water surfaces only)
1497       CALL calc_z0_water_surface
1498
1499
1500    END SUBROUTINE lsm_energy_balance
1501
1502
1503!------------------------------------------------------------------------------!
1504! Description:
1505! ------------
1506!> Soil model as part of the land surface model. The model predicts soil
1507!> temperature and water content.
1508!------------------------------------------------------------------------------!
1509    SUBROUTINE lsm_soil_model
1510
1511
1512       IMPLICIT NONE
1513
1514       INTEGER(iwp) ::  i   !< running index
1515       INTEGER(iwp) ::  j   !< running index
1516       INTEGER(iwp) ::  k   !< running index
1517
1518       REAL(wp)     :: h_vg !< Van Genuchten coef. h
1519
1520       REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp,  & !< temp. gamma
1521                                                 lambda_temp, & !< temp. lambda
1522                                                 tend           !< tendency
1523
1524       DO  i = nxlg, nxrg
1525          DO  j = nysg, nyng
1526
1527             IF ( pave_surface(j,i) )  THEN
1528                rho_c_total(nzb_soil,j,i) = pave_heat_capacity
1529                lambda_temp(nzb_soil)     = pave_heat_conductivity
1530             ENDIF
1531
1532             IF (  .NOT.  water_surface(j,i) )  THEN
1533                DO  k = nzb_soil, nzt_soil
1534
1535
1536                   IF ( pave_surface(j,i)  .AND.  zs(k) <= pave_depth )  THEN
1537                   
1538                      rho_c_total(k,j,i) = pave_heat_capacity
1539                      lambda_temp(k)     = pave_heat_conductivity   
1540
1541                   ELSE           
1542!
1543!--                   Calculate volumetric heat capacity of the soil, taking
1544!--                   into account water content
1545                      rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i)) &
1546                                           + rho_c_water * m_soil(k,j,i))
1547
1548!
1549!--                   Calculate soil heat conductivity at the center of the soil
1550!--                   layers
1551                      lambda_h_sat = lambda_h_sm ** (1.0_wp - m_sat(j,i)) *    &
1552                                     lambda_h_water ** m_soil(k,j,i)
1553
1554                      ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i)))
1555
1556                      lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) +    &
1557                                       lambda_h_dry
1558                   ENDIF
1559
1560                ENDDO
1561
1562!
1563!--             Calculate soil heat conductivity (lambda_h) at the _stag level
1564!--             using linear interpolation. For pavement surface, the
1565!--             true pavement depth is considered
1566                DO  k = nzb_soil, nzt_soil-1
1567                   IF ( pave_surface(j,i)  .AND.  zs(k)   < pave_depth         &
1568                                           .AND.  zs(k+1) > pave_depth )  THEN
1569                      lambda_h(k,j,i) = ( pave_depth - zs(k) ) / dz_soil(k+1)  &
1570                                        * lambda_temp(k)                       &
1571                                        + ( 1.0_wp - ( pave_depth - zs(k) )    &
1572                                        / dz_soil(k+1) ) * lambda_temp(k+1)
1573                   ELSE
1574                      lambda_h(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )  &
1575                                        * 0.5_wp
1576                   ENDIF
1577                ENDDO
1578                lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil)
1579
1580
1581
1582
1583!
1584!--             Prognostic equation for soil temperature t_soil
1585                tend(:) = 0.0_wp
1586
1587                tend(nzb_soil) = (1.0_wp/rho_c_total(nzb_soil,j,i)) *          &
1588                          ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i)  &
1589                            - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil+1)    &
1590                            + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil)
1591
1592                DO  k = nzb_soil+1, nzt_soil
1593                   tend(k) = (1.0_wp/rho_c_total(k,j,i))                       &
1594                             * (   lambda_h(k,j,i)                             &
1595                                 * ( t_soil(k+1,j,i) - t_soil(k,j,i) )         &
1596                                 * ddz_soil(k+1)                               &
1597                                 - lambda_h(k-1,j,i)                           &
1598                                 * ( t_soil(k,j,i) - t_soil(k-1,j,i) )         &
1599                                 * ddz_soil(k)                                 &
1600                               ) * ddz_soil_stag(k)
1601
1602                ENDDO
1603
1604                t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)&
1605                                                  + dt_3d * ( tsc(2)           &
1606                                                  * tend(nzb_soil:nzt_soil)    & 
1607                                                  + tsc(3)                     &
1608                                                  * tt_soil_m(:,j,i) )   
1609
1610!
1611!--             Calculate t_soil tendencies for the next Runge-Kutta step
1612                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1613                   IF ( intermediate_timestep_count == 1 )  THEN
1614                      DO  k = nzb_soil, nzt_soil
1615                         tt_soil_m(k,j,i) = tend(k)
1616                      ENDDO
1617                   ELSEIF ( intermediate_timestep_count <                      &
1618                            intermediate_timestep_count_max )  THEN
1619                      DO  k = nzb_soil, nzt_soil
1620                         tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp   &
1621                                         * tt_soil_m(k,j,i)
1622                      ENDDO
1623                   ENDIF
1624                ENDIF
1625
1626
1627                DO  k = nzb_soil, nzt_soil
1628
1629!
1630!--                Calculate soil diffusivity at the center of the soil layers
1631                   lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat       &
1632                                     / m_sat(j,i) ) * ( MAX( m_soil(k,j,i),    &
1633                                     m_wilt(j,i) ) / m_sat(j,i) )**(           &
1634                                     b_ch + 2.0_wp )
1635
1636!
1637!--                Parametrization of Van Genuchten
1638                   IF ( soil_type /= 7 )  THEN
1639!
1640!--                   Calculate the hydraulic conductivity after Van Genuchten
1641!--                   (1980)
1642                      h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -    &
1643                                 MAX( m_soil(k,j,i), m_wilt(j,i) ) ) )**(      &
1644                                 n_vg(j,i) / (n_vg(j,i) - 1.0_wp ) ) - 1.0_wp  &
1645                             )**( 1.0_wp / n_vg(j,i) ) / alpha_vg(j,i)
1646
1647
1648                      gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp +         &
1649                                      ( alpha_vg(j,i) * h_vg )**n_vg(j,i))**(  &
1650                                      1.0_wp - 1.0_wp / n_vg(j,i) ) - (        &
1651                                      alpha_vg(j,i) * h_vg )**( n_vg(j,i)      &
1652                                      - 1.0_wp) )**2 )                         &
1653                                      / ( ( 1.0_wp + ( alpha_vg(j,i) * h_vg    &
1654                                      )**n_vg(j,i) )**( ( 1.0_wp  - 1.0_wp     &
1655                                      / n_vg(j,i) ) *( l_vg(j,i) + 2.0_wp) ) )
1656
1657!
1658!--                Parametrization of Clapp & Hornberger
1659                   ELSE
1660                      gamma_temp(k) = gamma_w_sat(j,i) * ( m_soil(k,j,i)       &
1661                                      / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp)
1662                   ENDIF
1663
1664                ENDDO
1665
1666!
1667!--             Prognostic equation for soil moisture content. Only performed,
1668!--             when humidity is enabled in the atmosphere and the surface type
1669!--             is not pavement (implies dry soil below).
1670                IF ( humidity  .AND.  .NOT.  pave_surface(j,i) )  THEN
1671!
1672!--                Calculate soil diffusivity (lambda_w) at the _stag level
1673!--                using linear interpolation. To do: replace this with
1674!--                ECMWF-IFS Eq. 8.81
1675                   DO  k = nzb_soil, nzt_soil-1
1676                     
1677                      lambda_w(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )  &
1678                                        * 0.5_wp
1679                      gamma_w(k,j,i)  = ( gamma_temp(k+1) + gamma_temp(k) )    &
1680                                        * 0.5_wp
1681
1682                   ENDDO
1683
1684!
1685!
1686!--                In case of a closed bottom (= water content is conserved), set
1687!--                hydraulic conductivity to zero to that no water will be lost
1688!--                in the bottom layer.
1689                   IF ( conserve_water_content )  THEN
1690                      gamma_w(nzt_soil,j,i) = 0.0_wp
1691                   ELSE
1692                      gamma_w(nzt_soil,j,i) = gamma_temp(nzt_soil)
1693                   ENDIF     
1694
1695!--                The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v))
1696!--                ensures the mass conservation for water. The transpiration of
1697!--                plants equals the cumulative withdrawals by the roots in the
1698!--                soil. The scheme takes into account the availability of water
1699!--                in the soil layers as well as the root fraction in the
1700!--                respective layer. Layer with moisture below wilting point will
1701!--                not contribute, which reflects the preference of plants to
1702!--                take water from moister layers.
1703
1704!
1705!--                Calculate the root extraction (ECMWF 7.69, the sum of root_extr
1706!--                = 1). The energy balance solver guarantees a positive
1707!--                transpiration, so that there is no need for an additional check.
1708                   DO  k = nzb_soil, nzt_soil
1709                       IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
1710                          m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i)
1711                       ENDIF
1712                   ENDDO 
1713
1714                   IF ( m_total > 0.0_wp )  THEN
1715                      DO  k = nzb_soil, nzt_soil
1716                         IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
1717                            root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total
1718                         ELSE
1719                            root_extr(k) = 0.0_wp
1720                         ENDIF
1721                      ENDDO
1722                   ENDIF
1723
1724!
1725!--                Prognostic equation for soil water content m_soil.
1726                   tend(:) = 0.0_wp
1727
1728                   tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * (               &
1729                            m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) )    &
1730                            * ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - ( &
1731                               root_extr(nzb_soil) * qsws_veg_eb(j,i)          &
1732                               + qsws_soil_eb(j,i) ) * drho_l_lv )             &
1733                               * ddz_soil_stag(nzb_soil)
1734
1735                   DO  k = nzb_soil+1, nzt_soil-1
1736                      tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)          &
1737                                - m_soil(k,j,i) ) * ddz_soil(k+1)              &
1738                                - gamma_w(k,j,i)                               &
1739                                - lambda_w(k-1,j,i) * (m_soil(k,j,i) -         &
1740                                m_soil(k-1,j,i)) * ddz_soil(k)                 &
1741                                + gamma_w(k-1,j,i) - (root_extr(k)             &
1742                                * qsws_veg_eb(j,i) * drho_l_lv)                &
1743                                ) * ddz_soil_stag(k)
1744
1745                   ENDDO
1746                   tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i)                  &
1747                                           - lambda_w(nzt_soil-1,j,i)          &
1748                                           * (m_soil(nzt_soil,j,i)             &
1749                                           - m_soil(nzt_soil-1,j,i))           &
1750                                           * ddz_soil(nzt_soil)                &
1751                                           + gamma_w(nzt_soil-1,j,i) - (       &
1752                                             root_extr(nzt_soil)               &
1753                                           * qsws_veg_eb(j,i) * drho_l_lv  )   &
1754                                     ) * ddz_soil_stag(nzt_soil)             
1755
1756                   m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)&
1757                                                   + dt_3d * ( tsc(2) * tend(:)   &
1758                                                   + tsc(3) * tm_soil_m(:,j,i) )   
1759   
1760!
1761!--                Account for dry soils (find a better solution here!)
1762                   DO  k = nzb_soil, nzt_soil
1763                      IF ( m_soil_p(k,j,i) < 0.0_wp )  m_soil_p(k,j,i) = 0.0_wp
1764                   ENDDO
1765
1766!
1767!--                Calculate m_soil tendencies for the next Runge-Kutta step
1768                   IF ( timestep_scheme(1:5) == 'runge' )  THEN
1769                      IF ( intermediate_timestep_count == 1 )  THEN
1770                         DO  k = nzb_soil, nzt_soil
1771                            tm_soil_m(k,j,i) = tend(k)
1772                         ENDDO
1773                      ELSEIF ( intermediate_timestep_count <                   &
1774                               intermediate_timestep_count_max )  THEN
1775                         DO  k = nzb_soil, nzt_soil
1776                            tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp&
1777                                     * tm_soil_m(k,j,i)
1778                         ENDDO
1779                      ENDIF
1780                   ENDIF
1781
1782                ENDIF
1783
1784             ENDIF
1785
1786          ENDDO
1787       ENDDO
1788
1789    END SUBROUTINE lsm_soil_model
1790
1791
1792!------------------------------------------------------------------------------!
1793! Description:
1794! ------------
1795!> Calculation of roughness length for open water (lakes, ocean). The
1796!> parameterization follows Charnock (1955). Two different implementations
1797!> are available: as in ECMWF-IFS (Beljaars 1994) or as in FLake (Subin et al.
1798!> 2012)
1799!------------------------------------------------------------------------------!
1800    SUBROUTINE calc_z0_water_surface
1801
1802       USE control_parameters,                                                 &
1803           ONLY: g, kappa, molecular_viscosity
1804
1805       IMPLICIT NONE
1806
1807       INTEGER :: i  !< running index
1808       INTEGER :: j  !< running index
1809
1810       REAL(wp), PARAMETER :: alpha_ch  = 0.018_wp !< Charnock constant (0.01-0.11). Use 0.01 for FLake and 0.018 for ECMWF
1811!       REAL(wp), PARAMETER :: pr_number = 0.71_wp !< molecular Prandtl number in the Charnock parameterization (differs from prandtl_number)
1812!       REAL(wp), PARAMETER :: sc_number = 0.66_wp !< molecular Schmidt number in the Charnock parameterization
1813!       REAL(wp) :: re_0 !< near-surface roughness Reynolds number
1814
1815
1816       DO  i = nxlg, nxrg   
1817          DO  j = nysg, nyng
1818             IF ( water_surface(j,i) )  THEN
1819
1820!
1821!--             Disabled: FLake parameterization. Ideally, the Charnock
1822!--             coefficient should depend on the water depth and the fetch
1823!--             length
1824!                re_0 = z0(j,i) * us(j,i) / molecular_viscosity
1825!       
1826!                z0(j,i) = MAX( 0.1_wp * molecular_viscosity / us(j,i),            &
1827!                              alpha_ch * us(j,i) / g )
1828!
1829!                z0h(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 3.2_wp ) )
1830!                z0q(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 4.2_wp ) )
1831
1832!
1833!--              Set minimum roughness length for u* > 0.2
1834!                IF ( us(j,i) > 0.2_wp )  THEN
1835!                   z0h(j,i) = MAX( 1.0E-5_wp, z0h(j,i) )
1836!                   z0q(j,i) = MAX( 1.0E-5_wp, z0q(j,i) )
1837!                ENDIF
1838
1839!
1840!--             ECMWF IFS model parameterization after Beljaars (1994). At low
1841!--             wind speed, the sea surface becomes aerodynamically smooth and
1842!--             the roughness scales with the viscosity. At high wind speed, the
1843!--             Charnock relation is used.
1844                z0(j,i) =   ( 0.11_wp * molecular_viscosity / us(j,i) )        &
1845                          + ( alpha_ch * us(j,i)**2 / g )
1846
1847                z0h(j,i) = 0.40_wp * molecular_viscosity / us(j,i)
1848                z0q(j,i) = 0.62_wp * molecular_viscosity / us(j,i)
1849
1850             ENDIF
1851          ENDDO
1852       ENDDO
1853
1854    END SUBROUTINE calc_z0_water_surface
1855
1856
1857!------------------------------------------------------------------------------!
1858! Description:
1859! ------------
1860!> Calculation of specific humidity of the skin layer (surface). It is assumend
1861!> that the skin is always saturated.
1862!------------------------------------------------------------------------------!
1863    SUBROUTINE calc_q_surface
1864
1865       IMPLICIT NONE
1866
1867       INTEGER :: i              !< running index
1868       INTEGER :: j              !< running index
1869       INTEGER :: k              !< running index
1870
1871       REAL(wp) :: resistance    !< aerodynamic and soil resistance term
1872
1873       DO  i = nxlg, nxrg   
1874          DO  j = nysg, nyng
1875             k = nzb_s_inner(j,i)
1876
1877!
1878!--          Calculate water vapour pressure at saturation
1879             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface_p(j,i)   &
1880                   - 273.16_wp ) / ( t_surface_p(j,i) - 35.86_wp ) )
1881
1882!
1883!--          Calculate specific humidity at saturation
1884             q_s = 0.622_wp * e_s / surface_pressure
1885
1886             resistance = r_a(j,i) / (r_a(j,i) + r_s(j,i))
1887
1888!
1889!--          Calculate specific humidity at surface
1890             IF ( cloud_physics )  THEN
1891                q(k,j,i) = resistance * q_s + (1.0_wp - resistance)            &
1892                             * ( q(k+1,j,i) - ql(k+1,j,i) )
1893             ELSE
1894                q(k,j,i) = resistance * q_s + (1.0_wp - resistance)            &
1895                             * q(k+1,j,i)
1896             ENDIF
1897
1898!
1899!--          Update virtual potential temperature
1900             vpt(k,j,i) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) )
1901
1902          ENDDO
1903       ENDDO
1904
1905    END SUBROUTINE calc_q_surface
1906
1907
1908!------------------------------------------------------------------------------!
1909! Description:
1910! ------------
1911!> Swapping of timelevels
1912!------------------------------------------------------------------------------!
1913    SUBROUTINE lsm_swap_timelevel ( mod_count )
1914
1915       IMPLICIT NONE
1916
1917       INTEGER, INTENT(IN) :: mod_count
1918
1919#if defined( __nopointer )
1920
1921       t_surface    = t_surface_p
1922       t_soil       = t_soil_p
1923       IF ( humidity )  THEN
1924          m_soil    = m_soil_p
1925          m_liq_eb  = m_liq_eb_p
1926       ENDIF
1927
1928#else
1929   
1930       SELECT CASE ( mod_count )
1931
1932          CASE ( 0 )
1933
1934             t_surface  => t_surface_1; t_surface_p  => t_surface_2
1935             t_soil     => t_soil_1;    t_soil_p     => t_soil_2
1936             IF ( humidity )  THEN
1937                m_soil    => m_soil_1;   m_soil_p    => m_soil_2
1938                m_liq_eb  => m_liq_eb_1; m_liq_eb_p  => m_liq_eb_2
1939             ENDIF
1940
1941
1942          CASE ( 1 )
1943
1944             t_surface  => t_surface_2; t_surface_p  => t_surface_1
1945             t_soil     => t_soil_2;    t_soil_p     => t_soil_1
1946             IF ( humidity )  THEN
1947                m_soil    => m_soil_2;   m_soil_p    => m_soil_1
1948                m_liq_eb  => m_liq_eb_2; m_liq_eb_p  => m_liq_eb_1
1949             ENDIF
1950
1951       END SELECT
1952#endif
1953
1954    END SUBROUTINE lsm_swap_timelevel
1955
1956
1957 END MODULE land_surface_model_mod
Note: See TracBrowser for help on using the repository browser.