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

Last change on this file since 1784 was 1784, checked in by raasch, 6 years ago

last commit documented

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