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

Last change on this file since 1766 was 1758, checked in by maronga, 8 years ago

last commit documented

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