source: palm/trunk/SOURCE/land_surface_model_mod.f90 @ 1834

Last change on this file since 1834 was 1827, checked in by maronga, 8 years ago

last commit documented

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