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

Last change on this file since 1924 was 1857, checked in by maronga, 8 years ago

last commit documented

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