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

Last change on this file since 1723 was 1710, checked in by maronga, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 75.5 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 1710 2015-11-04 14:48:36Z boeske $
26!
27! 1709 2015-11-04 14:47:01Z maronga
28! Renamed pt_1 and qv_1 to pt1 and qv1.
29! Bugfix: set initial values for t_surface_p in case of restart runs
30! Bugfix: zero resistance caused crash when using radiation_scheme = 'clear-sky'
31! Bugfix: calculation of rad_net when using radiation_scheme = 'clear-sky'
32! Added todo action
33!
34! 1697 2015-10-28 17:14:10Z raasch
35! bugfix: misplaced cpp-directive
36!
37! 1695 2015-10-27 10:03:11Z maronga
38! Bugfix: REAL constants provided with KIND-attribute in call of
39! Replaced rif with ol
40!
41! 1691 2015-10-26 16:17:44Z maronga
42! Added skip_time_do_lsm to allow for spin-ups without LSM. Various bugfixes:
43! Soil temperatures are now defined at the edges of the layers, calculation of
44! shb_eb corrected, prognostic equation for skin temperature corrected. Surface
45! fluxes are now directly transfered to atmosphere
46!
47! 1682 2015-10-07 23:56:08Z knoop
48! Code annotations made doxygen readable
49!
50! 1590 2015-05-08 13:56:27Z maronga
51! Bugfix: definition of character strings requires same length for all elements
52!
53! 1585 2015-04-30 07:05:52Z maronga
54! Modifications for RRTMG. Changed tables to PARAMETER type.
55!
56! 1571 2015-03-12 16:12:49Z maronga
57! Removed upper-case variable names. Corrected distribution of precipitation to
58! the liquid water reservoir and the bare soil fractions.
59!
60! 1555 2015-03-04 17:44:27Z maronga
61! Added output of r_a and r_s
62!
63! 1553 2015-03-03 17:33:54Z maronga
64! Improved better treatment of roughness lengths. Added default soil temperature
65! profile
66!
67! 1551 2015-03-03 14:18:16Z maronga
68! Flux calculation is now done in prandtl_fluxes. Added support for data output.
69! Vertical indices have been replaced. Restart runs are now possible. Some
70! variables have beem renamed. Bugfix in the prognostic equation for the surface
71! temperature. Introduced z0_eb and z0h_eb, which overwrite the setting of
72! roughness_length and z0_factor. Added Clapp & Hornberger parametrization for
73! the hydraulic conductivity. Bugfix for root fraction and extraction
74! calculation
75!
76! intrinsic function MAX and MIN
77!
78! 1500 2014-12-03 17:42:41Z maronga
79! Corrected calculation of aerodynamic resistance (r_a).
80! Precipitation is now added to liquid water reservoir using LE_liq.
81! Added support for dry runs.
82!
83! 1496 2014-12-02 17:25:50Z maronga
84! Initial revision
85!
86!
87! Description:
88! ------------
89!> Land surface model, consisting of a solver for the energy balance at the
90!> surface and a four layer soil scheme. The scheme is similar to the TESSEL
91!> scheme implemented in the ECMWF IFS model, with modifications according to
92!> H-TESSEL. The implementation is based on the formulation implemented in the
93!> DALES and UCLA-LES models.
94!>
95!> @todo Consider partial absorption of the net shortwave radiation by the
96!>       skin layer.
97!> @todo Allow for water surfaces
98!> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0,
99!>       nzt_soil=3)).
100!> @todo Implement surface runoff model (required when performing long-term LES
101!>       with considerable precipitation.
102!> @todo Fix crashes with radiation_scheme == 'constant'
103!>
104!> @note No time step criterion is required as long as the soil layers do not
105!>       become too thin.
106!------------------------------------------------------------------------------!
107 MODULE land_surface_model_mod
108 
109    USE arrays_3d,                                                             &
110        ONLY:  hyp, ol, pt, pt_p, q, q_p, ql, qsws, shf, ts, us, vpt, z0, z0h
111
112    USE cloud_parameters,                                                      &
113        ONLY:  cp, hyrho, l_d_cp, l_d_r, l_v, prr, pt_d_t, rho_l, r_d, r_v
114
115    USE control_parameters,                                                    &
116        ONLY:  cloud_physics, dt_3d, humidity, intermediate_timestep_count,    &
117               initializing_actions, intermediate_timestep_count_max,          &
118               max_masks, precipitation, pt_surface, rho_surface,              &
119               roughness_length, surface_pressure, timestep_scheme, tsc,       &
120               z0h_factor, time_since_reference_point
121
122    USE indices,                                                               &
123        ONLY:  nbgp, nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner 
124
125    USE kinds
126
127    USE netcdf_control,                                                        &
128        ONLY:  dots_label, dots_num, dots_unit
129
130    USE pegrid
131
132    USE radiation_model_mod,                                                   &
133        ONLY:  force_radiation_call, radiation_scheme, rad_net, rad_sw_in,     &
134               rad_lw_out, rad_lw_out_change_0, sigma_sb
135       
136#if defined ( __rrtmg )
137    USE radiation_model_mod,                                                   &
138        ONLY:  rrtm_idrv
139#endif
140
141    USE statistics,                                                            &
142        ONLY:  hom, statistic_regions
143
144    IMPLICIT NONE
145
146!
147!-- LSM model constants
148    INTEGER(iwp), PARAMETER :: nzb_soil = 0, & !< bottom of the soil model (to be switched)
149                               nzt_soil = 3, & !< top of the soil model (to be switched)
150                               nzs = 4         !< number of soil layers (fixed for now)
151
152    INTEGER(iwp) :: dots_soil = 0  !< starting index for timeseries output
153
154    INTEGER(iwp), DIMENSION(0:1) :: id_dim_zs_xy, id_dim_zs_xz, id_dim_zs_yz,  &
155                                    id_dim_zs_3d, id_var_zs_xy,                & 
156                                    id_var_zs_xz, id_var_zs_yz, id_var_zs_3d
157                                   
158    INTEGER(iwp), DIMENSION(1:max_masks,0:1) :: id_dim_zs_mask, id_var_zs_mask
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 and NetCDF variables
514    PUBLIC dots_soil, id_dim_zs_xy, id_dim_zs_xz, id_dim_zs_yz,                &
515           id_dim_zs_3d, id_dim_zs_mask, id_var_zs_xy, id_var_zs_xz,           &
516           id_var_zs_yz, id_var_zs_3d, id_var_zs_mask, nzb_soil, nzs, nzt_soil,&
517           zs
518
519!
520!-- Public 2D output variables
521    PUBLIC c_liq, c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb, ghf_eb_av,     &
522           lai, lai_av, qsws_eb, qsws_eb_av, qsws_liq_eb, qsws_liq_eb_av,      &
523           qsws_soil_eb, qsws_soil_eb_av, qsws_veg_eb, qsws_veg_eb_av,         &
524           r_a, r_a_av, r_s, r_s_av, shf_eb, shf_eb_av
525
526!
527!-- Public prognostic variables
528    PUBLIC m_liq_eb, m_liq_eb_av, m_soil, m_soil_av, t_soil, t_soil_av
529
530    INTERFACE init_lsm
531       MODULE PROCEDURE init_lsm
532    END INTERFACE init_lsm
533
534    INTERFACE lsm_energy_balance
535       MODULE PROCEDURE lsm_energy_balance
536    END INTERFACE lsm_energy_balance
537
538    INTERFACE lsm_soil_model
539       MODULE PROCEDURE lsm_soil_model
540    END INTERFACE lsm_soil_model
541
542    INTERFACE lsm_swap_timelevel
543       MODULE PROCEDURE lsm_swap_timelevel
544    END INTERFACE lsm_swap_timelevel
545
546 CONTAINS
547
548
549!------------------------------------------------------------------------------!
550! Description:
551! ------------
552!> Allocate land surface model arrays and define pointers
553!------------------------------------------------------------------------------!
554    SUBROUTINE init_lsm_arrays
555   
556
557       IMPLICIT NONE
558
559!
560!--    Allocate surface and soil temperature / humidity
561#if defined( __nopointer )
562       ALLOCATE ( m_liq_eb(nysg:nyng,nxlg:nxrg) )
563       ALLOCATE ( m_liq_eb_p(nysg:nyng,nxlg:nxrg) )
564       ALLOCATE ( m_soil(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
565       ALLOCATE ( m_soil_p(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
566       ALLOCATE ( t_surface(nysg:nyng,nxlg:nxrg) )
567       ALLOCATE ( t_surface_p(nysg:nyng,nxlg:nxrg) )
568       ALLOCATE ( t_soil(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
569       ALLOCATE ( t_soil_p(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
570#else
571       ALLOCATE ( m_liq_eb_1(nysg:nyng,nxlg:nxrg) )
572       ALLOCATE ( m_liq_eb_2(nysg:nyng,nxlg:nxrg) )
573       ALLOCATE ( m_soil_1(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
574       ALLOCATE ( m_soil_2(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
575       ALLOCATE ( t_surface_1(nysg:nyng,nxlg:nxrg) )
576       ALLOCATE ( t_surface_2(nysg:nyng,nxlg:nxrg) )
577       ALLOCATE ( t_soil_1(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
578       ALLOCATE ( t_soil_2(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
579#endif
580
581!
582!--    Allocate intermediate timestep arrays
583       ALLOCATE ( tm_liq_eb_m(nysg:nyng,nxlg:nxrg) )
584       ALLOCATE ( tm_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
585       ALLOCATE ( tt_surface_m(nysg:nyng,nxlg:nxrg) )
586       ALLOCATE ( tt_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
587
588!
589!--    Allocate 2D vegetation model arrays
590       ALLOCATE ( alpha_vg(nysg:nyng,nxlg:nxrg) )
591       ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) )
592       ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) )
593       ALLOCATE ( f_sw_in(nysg:nyng,nxlg:nxrg) )
594       ALLOCATE ( ghf_eb(nysg:nyng,nxlg:nxrg) )
595       ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) )
596       ALLOCATE ( g_d(nysg:nyng,nxlg:nxrg) )
597       ALLOCATE ( lai(nysg:nyng,nxlg:nxrg) )
598       ALLOCATE ( l_vg(nysg:nyng,nxlg:nxrg) )
599       ALLOCATE ( lambda_surface_u(nysg:nyng,nxlg:nxrg) )
600       ALLOCATE ( lambda_surface_s(nysg:nyng,nxlg:nxrg) )
601       ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) )
602       ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) )
603       ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) )
604       ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) )
605       ALLOCATE ( n_vg(nysg:nyng,nxlg:nxrg) )
606       ALLOCATE ( qsws_eb(nysg:nyng,nxlg:nxrg) )
607       ALLOCATE ( qsws_soil_eb(nysg:nyng,nxlg:nxrg) )
608       ALLOCATE ( qsws_liq_eb(nysg:nyng,nxlg:nxrg) )
609       ALLOCATE ( qsws_veg_eb(nysg:nyng,nxlg:nxrg) )
610       ALLOCATE ( rad_net_l(nysg:nyng,nxlg:nxrg) )
611       ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) )
612       ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) )
613       ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) )
614       ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) )
615       ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) )
616       ALLOCATE ( r_canopy_min(nysg:nyng,nxlg:nxrg) )
617       ALLOCATE ( shf_eb(nysg:nyng,nxlg:nxrg) )
618
619#if ! defined( __nopointer )
620!
621!--    Initial assignment of the pointers
622       t_soil    => t_soil_1;    t_soil_p    => t_soil_2
623       t_surface => t_surface_1; t_surface_p => t_surface_2
624       m_soil    => m_soil_1;    m_soil_p    => m_soil_2
625       m_liq_eb  => m_liq_eb_1;  m_liq_eb_p  => m_liq_eb_2
626#endif
627
628
629    END SUBROUTINE init_lsm_arrays
630
631!------------------------------------------------------------------------------!
632! Description:
633! ------------
634!> Initialization of the land surface model
635!------------------------------------------------------------------------------!
636    SUBROUTINE init_lsm
637   
638
639       IMPLICIT NONE
640
641       INTEGER(iwp) ::  i !< running index
642       INTEGER(iwp) ::  j !< running index
643       INTEGER(iwp) ::  k !< running index
644
645       REAL(wp) :: pt1   !< potential temperature at first grid level
646
647
648!
649!--    Calculate Exner function
650       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
651
652
653!
654!--    If no cloud physics is used, rho_surface has not been calculated before
655       IF ( .NOT. cloud_physics )  THEN
656          rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn )
657       ENDIF
658
659!
660!--    Calculate frequently used parameters
661       rho_cp    = cp * rho_surface
662       rd_d_rv   = r_d / r_v
663       rho_lv    = rho_surface * l_v
664       drho_l_lv = 1.0_wp / (rho_l * l_v)
665
666!
667!--    Set inital values for prognostic quantities
668       tt_surface_m = 0.0_wp
669       tt_soil_m    = 0.0_wp
670       tm_liq_eb_m  = 0.0_wp
671       c_liq        = 0.0_wp
672
673       ghf_eb = 0.0_wp
674       shf_eb = rho_cp * shf
675
676       IF ( humidity )  THEN
677          qsws_eb = rho_l * l_v * qsws
678       ELSE
679          qsws_eb = 0.0_wp
680       ENDIF
681
682       qsws_liq_eb  = 0.0_wp
683       qsws_soil_eb = qsws_eb
684       qsws_veg_eb  = 0.0_wp
685
686       r_a        = 50.0_wp
687       r_s        = 50.0_wp
688       r_canopy   = 0.0_wp
689       r_soil     = 0.0_wp
690
691!
692!--    Allocate 3D soil model arrays
693       ALLOCATE ( root_fr(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
694       ALLOCATE ( lambda_h(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
695       ALLOCATE ( rho_c_total(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
696
697       lambda_h = 0.0_wp
698!
699!--    If required, allocate humidity-related variables for the soil model
700       IF ( humidity )  THEN
701          ALLOCATE ( lambda_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
702          ALLOCATE ( gamma_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )   
703
704          lambda_w = 0.0_wp 
705       ENDIF
706
707!
708!--    Calculate grid spacings. Temperature and moisture are defined at
709!--    the edges of the soil layers (_stag), whereas gradients/fluxes are defined
710!--    at the centers
711       dz_soil(nzb_soil) = zs(nzb_soil)
712
713       DO  k = nzb_soil+1, nzt_soil
714          dz_soil(k) = zs(k) - zs(k-1)
715       ENDDO
716       dz_soil(nzt_soil+1) = dz_soil(nzt_soil)
717
718       DO  k = nzb_soil, nzt_soil-1
719          dz_soil_stag(k) = 0.5_wp * (dz_soil(k+1) + dz_soil(k))
720       ENDDO
721       dz_soil_stag(nzt_soil) = dz_soil(nzt_soil)
722
723       ddz_soil      = 1.0_wp / dz_soil
724       ddz_soil_stag = 1.0_wp / dz_soil_stag
725
726!
727!--    Initialize standard soil types. It is possible to overwrite each
728!--    parameter by setting the respecticy NAMELIST variable to a
729!--    value /= 9999999.9.
730       IF ( soil_type /= 0 )  THEN 
731 
732          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
733             alpha_vangenuchten = soil_pars(0,soil_type)
734          ENDIF
735
736          IF ( l_vangenuchten == 9999999.9_wp )  THEN
737             l_vangenuchten = soil_pars(1,soil_type)
738          ENDIF
739
740          IF ( n_vangenuchten == 9999999.9_wp )  THEN
741             n_vangenuchten = soil_pars(2,soil_type)           
742          ENDIF
743
744          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
745             hydraulic_conductivity = soil_pars(3,soil_type)           
746          ENDIF
747
748          IF ( saturation_moisture == 9999999.9_wp )  THEN
749             saturation_moisture = m_soil_pars(0,soil_type)           
750          ENDIF
751
752          IF ( field_capacity == 9999999.9_wp )  THEN
753             field_capacity = m_soil_pars(1,soil_type)           
754          ENDIF
755
756          IF ( wilting_point == 9999999.9_wp )  THEN
757             wilting_point = m_soil_pars(2,soil_type)           
758          ENDIF
759
760          IF ( residual_moisture == 9999999.9_wp )  THEN
761             residual_moisture = m_soil_pars(3,soil_type)       
762          ENDIF
763
764       ENDIF   
765
766       alpha_vg      = alpha_vangenuchten
767       l_vg          = l_vangenuchten
768       n_vg          = n_vangenuchten 
769       gamma_w_sat   = hydraulic_conductivity
770       m_sat         = saturation_moisture
771       m_fc          = field_capacity
772       m_wilt        = wilting_point
773       m_res         = residual_moisture
774       r_soil_min    = min_soil_resistance
775
776!
777!--    Initial run actions
778       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
779
780          t_soil    = 0.0_wp
781          m_liq_eb  = 0.0_wp
782          m_soil    = 0.0_wp
783
784!
785!--       Map user settings of T and q for each soil layer
786!--       (make sure that the soil moisture does not drop below the permanent
787!--       wilting point) -> problems with devision by zero)
788          DO  k = nzb_soil, nzt_soil
789             t_soil(k,:,:)    = soil_temperature(k)
790             m_soil(k,:,:)    = MAX(soil_moisture(k),m_wilt(:,:))
791             soil_moisture(k) = MAX(soil_moisture(k),wilting_point)
792          ENDDO
793          t_soil(nzt_soil+1,:,:) = soil_temperature(nzt_soil+1)
794
795!
796!--       Calculate surface temperature
797          t_surface   = pt_surface * exn
798
799!
800!--       Set artifical values for ts and us so that r_a has its initial value for
801!--       the first time step
802          DO  i = nxlg, nxrg
803             DO  j = nysg, nyng
804                k = nzb_s_inner(j,i)
805
806                IF ( cloud_physics )  THEN
807                   pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
808                ELSE
809                   pt1 = pt(k+1,j,i)
810                ENDIF
811
812!
813!--             Assure that r_a cannot be zero at model start
814                IF ( pt1 == pt(k,j,i) )  pt1 = pt1 + 1.0E-10_wp
815
816                us(j,i)  = 0.1_wp
817                ts(j,i)  = (pt1 - pt(k,j,i)) / r_a(j,i)
818                shf(j,i) = - us(j,i) * ts(j,i)
819             ENDDO
820          ENDDO
821
822!
823!--    Actions for restart runs
824       ELSE
825
826          DO  i = nxlg, nxrg
827             DO  j = nysg, nyng
828                k = nzb_s_inner(j,i)               
829                t_surface(j,i) = pt(k,j,i) * exn
830             ENDDO
831          ENDDO
832
833       ENDIF
834
835       DO  k = nzb_soil, nzt_soil
836          root_fr(k,:,:) = root_fraction(k)
837       ENDDO
838
839       IF ( veg_type /= 0 )  THEN
840          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
841             min_canopy_resistance = veg_pars(0,veg_type)
842          ENDIF
843          IF ( leaf_area_index == 9999999.9_wp )  THEN
844             leaf_area_index = veg_pars(1,veg_type)         
845          ENDIF
846          IF ( vegetation_coverage == 9999999.9_wp )  THEN
847             vegetation_coverage = veg_pars(2,veg_type)     
848          ENDIF
849          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
850              canopy_resistance_coefficient= veg_pars(3,veg_type)     
851          ENDIF
852          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
853             lambda_surface_stable = surface_pars(0,veg_type)         
854          ENDIF
855          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
856             lambda_surface_unstable = surface_pars(1,veg_type)       
857          ENDIF
858          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
859             f_shortwave_incoming = surface_pars(2,veg_type)       
860          ENDIF
861          IF ( z0_eb == 9999999.9_wp )  THEN
862             roughness_length = roughness_par(0,veg_type) 
863             z0_eb            = roughness_par(0,veg_type) 
864          ENDIF
865          IF ( z0h_eb == 9999999.9_wp )  THEN
866             z0h_eb = roughness_par(1,veg_type)
867          ENDIF
868          z0h_factor = z0h_eb / z0_eb
869
870          IF ( ANY( root_fraction == 9999999.9_wp ) )  THEN
871             DO  k = nzb_soil, nzt_soil
872                root_fr(k,:,:) = root_distribution(k,veg_type)
873                root_fraction(k) = root_distribution(k,veg_type)
874             ENDDO
875          ENDIF
876
877       ELSE
878
879          IF ( z0_eb == 9999999.9_wp )  THEN
880             z0_eb = roughness_length
881          ENDIF
882          IF ( z0h_eb == 9999999.9_wp )  THEN
883             z0h_eb = z0_eb * z0h_factor
884          ENDIF
885
886       ENDIF
887
888!
889!--    Initialize vegetation
890       r_canopy_min         = min_canopy_resistance
891       lai                  = leaf_area_index
892       c_veg                = vegetation_coverage
893       g_d                  = canopy_resistance_coefficient
894       lambda_surface_s     = lambda_surface_stable
895       lambda_surface_u     = lambda_surface_unstable
896       f_sw_in              = f_shortwave_incoming
897       z0                   = z0_eb
898       z0h                  = z0h_eb
899
900!
901!--    Possibly do user-defined actions (e.g. define heterogeneous land surface)
902       CALL user_init_land_surface
903
904
905       t_soil_p    = t_soil
906       m_soil_p    = m_soil
907       m_liq_eb_p  = m_liq_eb
908       t_surface_p = t_surface
909
910
911
912!--    Store initial profiles of t_soil and m_soil (assuming they are
913!--    horizontally homogeneous on this PE)
914       hom(nzb_soil:nzt_soil,1,90,:)  = SPREAD( t_soil(nzb_soil:nzt_soil,      &
915                                                nysg,nxlg), 2,                 &
916                                                statistic_regions+1 )
917       hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil(nzb_soil:nzt_soil,      &
918                                                nysg,nxlg), 2,                 &
919                                                statistic_regions+1 )
920
921!
922!--    Add timeseries for land surface model
923       dots_soil = dots_num + 1
924       dots_num  = dots_num + 8
925
926       dots_label(dots_soil) = "ghf_eb"
927       dots_label(dots_soil+1) = "shf_eb"
928       dots_label(dots_soil+2) = "qsws_eb"
929       dots_label(dots_soil+3) = "qsws_liq_eb"
930       dots_label(dots_soil+4) = "qsws_soil_eb"
931       dots_label(dots_soil+5) = "qsws_veg_eb"
932       dots_label(dots_soil+6) = "r_a"
933       dots_label(dots_soil+7) = "r_s"
934
935       dots_unit(dots_soil:dots_soil+5) = "W/m2"
936       dots_unit(dots_soil+6:dots_soil+7) = "s/m"
937
938
939    END SUBROUTINE init_lsm
940
941
942
943!------------------------------------------------------------------------------!
944! Description:
945! ------------
946!> Solver for the energy balance at the surface.
947!------------------------------------------------------------------------------!
948    SUBROUTINE lsm_energy_balance
949
950
951       IMPLICIT NONE
952
953       INTEGER(iwp) ::  i         !< running index
954       INTEGER(iwp) ::  j         !< running index
955       INTEGER(iwp) ::  k, ks     !< running index
956
957       REAL(wp) :: f1,          & !< resistance correction term 1
958                   f2,          & !< resistance correction term 2
959                   f3,          & !< resistance correction term 3
960                   m_min,       & !< minimum soil moisture
961                   e,           & !< water vapour pressure
962                   e_s,         & !< water vapour saturation pressure
963                   e_s_dt,      & !< derivate of e_s with respect to T
964                   tend,        & !< tendency
965                   dq_s_dt,     & !< derivate of q_s with respect to T
966                   coef_1,      & !< coef. for prognostic equation
967                   coef_2,      & !< coef. for prognostic equation
968                   f_qsws,      & !< factor for qsws_eb
969                   f_qsws_veg,  & !< factor for qsws_veg_eb
970                   f_qsws_soil, & !< factor for qsws_soil_eb
971                   f_qsws_liq,  & !< factor for qsws_liq_eb
972                   f_shf,       & !< factor for shf_eb
973                   lambda_surface, & !< Current value of lambda_surface
974                   m_liq_eb_max,   & !< maxmimum value of the liq. water reservoir
975                   pt1,         & !< potential temperature at first grid level
976                   qv1            !< specific humidity at first grid level
977
978!
979!--    Calculate the exner function for the current time step
980       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
981
982       DO  i = nxlg, nxrg
983          DO  j = nysg, nyng
984             k = nzb_s_inner(j,i)
985
986!
987!--          Set lambda_surface according to stratification
988             IF ( ol(j,i) >= 0.0_wp )  THEN
989                lambda_surface = lambda_surface_s(j,i)
990             ELSE
991                lambda_surface = lambda_surface_u(j,i)
992             ENDIF
993
994!
995!--          First step: calculate aerodyamic resistance. As pt, us, ts
996!--          are not available for the prognostic time step, data from the last
997!--          time step is used here. Note that this formulation is the
998!--          equivalent to the ECMWF formulation using drag coefficients
999             IF ( cloud_physics )  THEN
1000                pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
1001                qv1 = q(k+1,j,i) - ql(k+1,j,i)
1002             ELSE
1003                pt1 = pt(k+1,j,i)
1004                qv1 = q(k+1,j,i)
1005             ENDIF
1006
1007             r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp)
1008
1009!
1010!--          Make sure that the resistance does not drop to zero
1011             IF ( ABS(r_a(j,i)) < 1.0E-10_wp )  r_a(j,i) = 1.0E-10_wp
1012
1013!
1014!--          Second step: calculate canopy resistance r_canopy
1015!--          f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
1016 
1017!--          f1: correction for incoming shortwave radiation (stomata close at
1018!--          night)
1019             IF ( radiation_scheme /= 'constant' )  THEN
1020                f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /  &
1021                              (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)          &
1022                               + 1.0_wp)) )
1023             ELSE
1024                f1 = 1.0_wp
1025             ENDIF
1026
1027
1028!
1029!--          f2: correction for soil moisture availability to plants (the
1030!--          integrated soil moisture must thus be considered here)
1031!--          f2 = 0 for very dry soils
1032             m_total = 0.0_wp
1033             DO  ks = nzb_soil, nzt_soil
1034                 m_total = m_total + root_fr(ks,j,i)                           &
1035                           * MAX(m_soil(ks,j,i),m_wilt(j,i))
1036             ENDDO 
1037
1038             IF ( m_total > m_wilt(j,i) .AND. m_total < m_fc(j,i) )  THEN
1039                f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
1040             ELSEIF ( m_total >= m_fc(j,i) )  THEN
1041                f2 = 1.0_wp
1042             ELSE
1043                f2 = 1.0E-20_wp
1044             ENDIF
1045
1046!
1047!--          Calculate water vapour pressure at saturation
1048             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)     &
1049                           - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) )
1050
1051!
1052!--          f3: correction for vapour pressure deficit
1053             IF ( g_d(j,i) /= 0.0_wp )  THEN
1054!
1055!--             Calculate vapour pressure
1056                e  = qv1 * surface_pressure / 0.622_wp
1057                f3 = EXP ( -g_d(j,i) * (e_s - e) )
1058             ELSE
1059                f3 = 1.0_wp
1060             ENDIF
1061
1062!
1063!--          Calculate canopy resistance. In case that c_veg is 0 (bare soils),
1064!--          this calculation is obsolete, as r_canopy is not used below.
1065!--          To do: check for very dry soil -> r_canopy goes to infinity
1066             r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3      &
1067                                             + 1.0E-20_wp)
1068
1069!
1070!--          Third step: calculate bare soil resistance r_soil. The Clapp &
1071!--          Hornberger parametrization does not consider c_veg.
1072             IF ( soil_type /= 7 )  THEN
1073                m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *     &
1074                        m_res(j,i)
1075             ELSE
1076                m_min = m_wilt(j,i)
1077             ENDIF
1078
1079             f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min )
1080             f2 = MAX(f2,1.0E-20_wp)
1081             f2 = MIN(f2,1.0_wp)
1082
1083             r_soil(j,i) = r_soil_min(j,i) / f2
1084
1085!
1086!--          Calculate fraction of liquid water reservoir
1087             m_liq_eb_max = m_max_depth * lai(j,i)
1088             c_liq(j,i) = MIN(1.0_wp, m_liq_eb(j,i) / (m_liq_eb_max+1.0E-20_wp))
1089             
1090
1091!
1092!--          Calculate saturation specific humidity
1093             q_s = 0.622_wp * e_s / surface_pressure
1094
1095!
1096!--          In case of dewfall, set evapotranspiration to zero
1097!--          All super-saturated water is then removed from the air
1098             IF ( humidity .AND. q_s <= qv1 )  THEN
1099                r_canopy(j,i) = 0.0_wp
1100                r_soil(j,i)   = 0.0_wp
1101             ENDIF
1102
1103!
1104!--          Calculate coefficients for the total evapotranspiration
1105             f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/        &
1106                           (r_a(j,i) + r_canopy(j,i))
1107
1108             f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +        &
1109                                                          r_soil(j,i))
1110             f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
1111
1112
1113!
1114!--          If soil moisture is below wilting point, plants do no longer
1115!--          transpirate.
1116!              IF ( m_soil(k,j,i) < m_wilt(j,i) )  THEN
1117!                 f_qsws_veg = 0.0_wp
1118!              ENDIF
1119
1120!
1121!--          If dewfall is deactivated, vegetation, soil and liquid water
1122!--          reservoir are not allowed to take up water from the super-saturated
1123!--          air.
1124             IF ( humidity .AND. q_s <= qv1 .AND. .NOT. dewfall )  THEN
1125                f_qsws_veg  = 0.0_wp
1126                f_qsws_soil = 0.0_wp
1127                f_qsws_liq  = 0.0_wp
1128             ENDIF
1129
1130             f_shf  = rho_cp / r_a(j,i)
1131             f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
1132
1133!
1134!--          Calculate derivative of q_s for Taylor series expansion
1135             e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) -        &
1136                              17.269_wp*(t_surface(j,i) - 273.16_wp)           &
1137                              / (t_surface(j,i) - 35.86_wp)**2 )
1138
1139             dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
1140
1141!
1142!--          Add LW up so that it can be removed in prognostic equation
1143             rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i)
1144
1145             IF ( humidity )  THEN
1146#if defined ( __rrtmg )
1147!
1148!--             Numerator of the prognostic equation
1149                coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
1150                         * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
1151                         + f_shf * pt1 + f_qsws * ( qv1 - q_s                  &
1152                         + dq_s_dt * t_surface(j,i) ) + lambda_surface         &
1153                         * t_soil(nzb_soil,j,i)
1154
1155!
1156!--             Denominator of the prognostic equation
1157                coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt           &
1158                         + lambda_surface + f_shf / exn
1159#else
1160
1161!
1162!--             Numerator of the prognostic equation
1163                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
1164                         * t_surface(j,i) ** 4                                 &
1165                         + f_shf * pt1 + f_qsws * ( qv1                        &
1166                         - q_s + dq_s_dt * t_surface(j,i) )                    &
1167                         + lambda_surface * t_soil(nzb_soil,j,i)
1168
1169!
1170!--             Denominator of the prognostic equation
1171                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws      &
1172                         * dq_s_dt + lambda_surface + f_shf / exn
1173 
1174#endif
1175             ELSE
1176
1177#if defined ( __rrtmg )
1178!
1179!--             Numerator of the prognostic equation
1180                coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
1181                         * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
1182                         + f_shf * pt1  + lambda_surface                       &
1183                         * t_soil(nzb_soil,j,i)
1184
1185!
1186!--             Denominator of the prognostic equation
1187                coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
1188#else
1189
1190!
1191!--             Numerator of the prognostic equation
1192                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
1193                          * t_surface(j,i) ** 4 + f_shf * pt1                  &
1194                          + lambda_surface * t_soil(nzb_soil,j,i)
1195
1196!
1197!--             Denominator of the prognostic equation
1198                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3               &
1199                         + lambda_surface + f_shf / exn
1200#endif
1201             ENDIF
1202
1203             tend = 0.0_wp
1204
1205!
1206!--          Implicit solution when the surface layer has no heat capacity,
1207!--          otherwise use RK3 scheme.
1208             t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface *        &
1209                                t_surface(j,i) ) / ( c_surface + coef_2 * dt_3d&
1210                                * tsc(2) ) 
1211!
1212!--          Add RK3 term
1213             IF ( c_surface /= 0.0_wp )  THEN
1214
1215                t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)           &
1216                                   * tt_surface_m(j,i)
1217
1218!
1219!--             Calculate true tendency
1220                tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3)     &
1221                       * tt_surface_m(j,i)) / (dt_3d  * tsc(2))
1222!
1223!--             Calculate t_surface tendencies for the next Runge-Kutta step
1224                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1225                   IF ( intermediate_timestep_count == 1 )  THEN
1226                      tt_surface_m(j,i) = tend
1227                   ELSEIF ( intermediate_timestep_count <                      &
1228                            intermediate_timestep_count_max )  THEN
1229                      tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp        &
1230                                          * tt_surface_m(j,i)
1231                   ENDIF
1232                ENDIF
1233
1234             ENDIF
1235
1236!
1237!--          In case of fast changes in the skin temperature, it is required to
1238!--          update the radiative fluxes in order to keep the solution stable
1239             IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp )  THEN
1240                force_radiation_call_l = .TRUE.
1241             ENDIF
1242
1243             pt(k,j,i) = t_surface_p(j,i) / exn
1244
1245!
1246!--          Calculate fluxes
1247#if defined ( __rrtmg )
1248             rad_net_l(j,i)   = rad_net_l(j,i) + rad_lw_out_change_0(j,i)      &
1249                                * t_surface(j,i) - rad_lw_out(nzb,j,i)         &
1250                                - rad_lw_out_change_0(j,i) * t_surface_p(j,i)
1251
1252             IF ( rrtm_idrv == 1 )  THEN
1253                rad_net(j,i) = rad_net_l(j,i)
1254                rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i)                      &
1255                                      + rad_lw_out_change_0(j,i)               &
1256                                      * ( t_surface_p(j,i) - t_surface(j,i) )
1257             ENDIF
1258#else
1259             rad_net_l(j,i)   = rad_net_l(j,i) + 3.0_wp * sigma_sb             &
1260                                * t_surface(j,i)**4 - 4.0_wp * sigma_sb        &
1261                                * t_surface(j,i)**3 * t_surface_p(j,i)
1262#endif
1263
1264
1265             ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)               &
1266                              - t_soil(nzb_soil,j,i))
1267
1268             shf_eb(j,i)    = - f_shf * ( pt1 - pt(k,j,i) )
1269
1270             shf(j,i) = shf_eb(j,i) / rho_cp
1271
1272             IF ( humidity )  THEN
1273                qsws_eb(j,i)  = - f_qsws    * ( qv1 - q_s + dq_s_dt            &
1274                                * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
1275
1276                qsws(j,i) = qsws_eb(j,i) / rho_lv
1277
1278                qsws_veg_eb(j,i)  = - f_qsws_veg  * ( qv1 - q_s                &
1279                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1280                                    * t_surface_p(j,i) )
1281
1282                qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s                &
1283                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1284                                    * t_surface_p(j,i) )
1285
1286                qsws_liq_eb(j,i)  = - f_qsws_liq  * ( qv1 - q_s                &
1287                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1288                                    * t_surface_p(j,i) )
1289             ENDIF
1290!
1291!--          Calculate the true surface resistance
1292             IF ( qsws_eb(j,i) == 0.0_wp )  THEN
1293                r_s(j,i) = 1.0E10_wp
1294             ELSE
1295                r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt                   &
1296                           * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )     &
1297                           / qsws_eb(j,i) - r_a(j,i)
1298             ENDIF
1299
1300!
1301!--          Calculate change in liquid water reservoir due to dew fall or
1302!--          evaporation of liquid water
1303             IF ( humidity )  THEN
1304!
1305!--             If precipitation is activated, add rain water to qsws_liq_eb
1306!--             and qsws_soil_eb according the the vegetation coverage.
1307!--             precipitation_rate is given in mm.
1308                IF ( precipitation )  THEN
1309
1310!
1311!--                Add precipitation to liquid water reservoir, if possible.
1312!--                Otherwise, add the water to bare soil fraction.
1313                   IF ( m_liq_eb(j,i) /= m_liq_eb_max )  THEN
1314                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i)                      &
1315                                       + c_veg(j,i) * prr(k,j,i) * hyrho(k)    &
1316                                       * 0.001_wp * rho_l * l_v
1317                   ELSE
1318                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
1319                                        + c_veg(j,i) * prr(k,j,i) * hyrho(k)   &
1320                                        * 0.001_wp * rho_l * l_v
1321                   ENDIF
1322
1323!--                Add precipitation to bare soil according to the bare soil
1324!--                coverage.
1325                   qsws_soil_eb(j,i) = qsws_soil_eb(j,i) * (1.0_wp             &
1326                                       - c_veg(j,i)) * prr(k,j,i) * hyrho(k)   &
1327                                       * 0.001_wp * rho_l * l_v
1328                ENDIF
1329
1330!
1331!--             If the air is saturated, check the reservoir water level
1332                IF ( qsws_eb(j,i) < 0.0_wp )  THEN
1333
1334!
1335!--                Check if reservoir is full (avoid values > m_liq_eb_max)
1336!--                In that case, qsws_liq_eb goes to qsws_soil_eb. In this
1337!--                case qsws_veg_eb is zero anyway (because c_liq = 1),       
1338!--                so that tend is zero and no further check is needed
1339                   IF ( m_liq_eb(j,i) == m_liq_eb_max )  THEN
1340                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
1341                                           + qsws_liq_eb(j,i)
1342                      qsws_liq_eb(j,i)  = 0.0_wp
1343                   ENDIF
1344
1345!
1346!--                In case qsws_veg_eb becomes negative (unphysical behavior),
1347!--                let the water enter the liquid water reservoir as dew on the
1348!--                plant
1349                   IF ( qsws_veg_eb(j,i) < 0.0_wp )  THEN
1350                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i)
1351                      qsws_veg_eb(j,i) = 0.0_wp
1352                   ENDIF
1353                ENDIF                   
1354 
1355                tend = - qsws_liq_eb(j,i) * drho_l_lv
1356
1357                m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend      &
1358                                                   + tsc(3) * tm_liq_eb_m(j,i) )
1359
1360!
1361!--             Check if reservoir is overfull -> reduce to maximum
1362!--             (conservation of water is violated here)
1363                m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max)
1364
1365!
1366!--             Check if reservoir is empty (avoid values < 0.0)
1367!--             (conservation of water is violated here)
1368                m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp)
1369
1370
1371!
1372!--             Calculate m_liq_eb tendencies for the next Runge-Kutta step
1373                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1374                   IF ( intermediate_timestep_count == 1 )  THEN
1375                      tm_liq_eb_m(j,i) = tend
1376                   ELSEIF ( intermediate_timestep_count <                      &
1377                            intermediate_timestep_count_max )  THEN
1378                      tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp         &
1379                                      * tm_liq_eb_m(j,i)
1380                   ENDIF
1381                ENDIF
1382
1383             ENDIF
1384
1385          ENDDO
1386       ENDDO
1387
1388!
1389!--    Make a logical OR for all processes. Force radiation call if at
1390!--    least one processor
1391       IF ( intermediate_timestep_count == intermediate_timestep_count_max-1 ) &
1392       THEN
1393#if defined( __parallel )
1394          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1395          CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,    &
1396                              1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
1397#else
1398          force_radiation_call = force_radiation_call_l
1399#endif
1400          force_radiation_call_l = .FALSE.
1401       ENDIF
1402
1403
1404    END SUBROUTINE lsm_energy_balance
1405
1406
1407!------------------------------------------------------------------------------!
1408! Description:
1409! ------------
1410!> Soil model as part of the land surface model. The model predicts soil
1411!> temperature and water content.
1412!------------------------------------------------------------------------------!
1413    SUBROUTINE lsm_soil_model
1414
1415
1416       IMPLICIT NONE
1417
1418       INTEGER(iwp) ::  i   !< running index
1419       INTEGER(iwp) ::  j   !< running index
1420       INTEGER(iwp) ::  k   !< running index
1421
1422       REAL(wp)     :: h_vg !< Van Genuchten coef. h
1423
1424       REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp,  & !< temp. gamma
1425                                                 lambda_temp, & !< temp. lambda
1426                                                 tend           !< tendency
1427
1428       DO  i = nxlg, nxrg
1429          DO  j = nysg, nyng
1430             DO  k = nzb_soil, nzt_soil
1431!
1432!--             Calculate volumetric heat capacity of the soil, taking into
1433!--             account water content
1434                rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i))       &
1435                                     + rho_c_water * m_soil(k,j,i))
1436
1437!
1438!--             Calculate soil heat conductivity at the center of the soil
1439!--             layers
1440                lambda_h_sat = lambda_h_sm ** (1.0_wp - m_sat(j,i)) *          &
1441                               lambda_h_water ** m_soil(k,j,i)
1442
1443                ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i)))
1444
1445                lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) +          &
1446                                 lambda_h_dry
1447
1448             ENDDO
1449
1450!
1451!--          Calculate soil heat conductivity (lambda_h) at the _stag level
1452!--          using linear interpolation
1453             DO  k = nzb_soil, nzt_soil-1
1454                lambda_h(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) ) * 0.5_wp
1455             ENDDO
1456             lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil)
1457
1458!
1459!--          Prognostic equation for soil temperature t_soil
1460             tend(:) = 0.0_wp
1461             tend(nzb_soil) = (1.0_wp/rho_c_total(nzb_soil,j,i)) *             &
1462                       ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i)     &
1463                         - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil+1)       &
1464                         + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil)
1465
1466
1467
1468
1469             DO  k = nzb_soil+1, nzt_soil
1470                tend(k) = (1.0_wp/rho_c_total(k,j,i))                          &
1471                          * (   lambda_h(k,j,i)                                &
1472                              * ( t_soil(k+1,j,i) - t_soil(k,j,i) )            &
1473                              * ddz_soil(k+1)                                  &
1474                              - lambda_h(k-1,j,i)                              &
1475                              * ( t_soil(k,j,i) - t_soil(k-1,j,i) )            &
1476                              * ddz_soil(k)                                    &
1477                            ) * ddz_soil_stag(k)
1478             ENDDO
1479
1480             t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)   &
1481                                             + dt_3d * ( tsc(2)                &
1482                                             * tend(:) + tsc(3)                &
1483                                             * tt_soil_m(:,j,i) )   
1484
1485!
1486!--          Calculate t_soil tendencies for the next Runge-Kutta step
1487             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1488                IF ( intermediate_timestep_count == 1 )  THEN
1489                   DO  k = nzb_soil, nzt_soil
1490                      tt_soil_m(k,j,i) = tend(k)
1491                   ENDDO
1492                ELSEIF ( intermediate_timestep_count <                         &
1493                         intermediate_timestep_count_max )  THEN
1494                   DO  k = nzb_soil, nzt_soil
1495                      tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp      &
1496                                         * tt_soil_m(k,j,i)
1497                   ENDDO
1498                ENDIF
1499             ENDIF
1500
1501
1502             DO  k = nzb_soil, nzt_soil
1503
1504!
1505!--             Calculate soil diffusivity at the center of the soil layers
1506                lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat          &
1507                                  / m_sat(j,i) ) * ( MAX(m_soil(k,j,i),        &
1508                                  m_wilt(j,i)) / m_sat(j,i) )**(b_ch + 2.0_wp)
1509
1510!
1511!--             Parametrization of Van Genuchten
1512                IF ( soil_type /= 7 )  THEN
1513!
1514!--                Calculate the hydraulic conductivity after Van Genuchten
1515!--                (1980)
1516                   h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -       &
1517                              MAX(m_soil(k,j,i),m_wilt(j,i)) ) )**(n_vg(j,i)   &
1518                              / (n_vg(j,i)-1.0_wp)) - 1.0_wp                   &
1519                          )**(1.0_wp/n_vg(j,i)) / alpha_vg(j,i)
1520
1521
1522                   gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp +            &
1523                                   (alpha_vg(j,i)*h_vg)**n_vg(j,i))**(1.0_wp   &
1524                                   -1.0_wp/n_vg(j,i)) - (alpha_vg(j,i)*h_vg    &
1525                                   )**(n_vg(j,i)-1.0_wp))**2 )                 &
1526                                   / ( (1.0_wp + (alpha_vg(j,i)*h_vg           &
1527                                   )**n_vg(j,i))**((1.0_wp - 1.0_wp/n_vg(j,i)) &
1528                                   *(l_vg(j,i) + 2.0_wp)) )
1529
1530!
1531!--             Parametrization of Clapp & Hornberger
1532                ELSE
1533                   gamma_temp(k) = gamma_w_sat(j,i) * (m_soil(k,j,i)           &
1534                                   / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp)
1535                ENDIF
1536
1537             ENDDO
1538
1539
1540             IF ( humidity )  THEN
1541!
1542!--             Calculate soil diffusivity (lambda_w) at the _stag level
1543!--             using linear interpolation. To do: replace this with
1544!--             ECMWF-IFS Eq. 8.81
1545                DO  k = nzb_soil, nzt_soil-1
1546                     
1547                   lambda_w(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )     &
1548                                     * 0.5_wp
1549                   gamma_w(k,j,i)  = ( gamma_temp(k+1) + gamma_temp(k) )       &
1550                                     * 0.5_wp
1551
1552                ENDDO
1553
1554!
1555!
1556!--             In case of a closed bottom (= water content is conserved), set
1557!--             hydraulic conductivity to zero to that no water will be lost
1558!--             in the bottom layer.
1559                IF ( conserve_water_content )  THEN
1560                   gamma_w(nzt_soil,j,i) = 0.0_wp
1561                ELSE
1562                   gamma_w(nzt_soil,j,i) = gamma_temp(nzt_soil)
1563                ENDIF     
1564
1565!--             The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v))
1566!--             ensures the mass conservation for water. The transpiration of
1567!--             plants equals the cumulative withdrawals by the roots in the
1568!--             soil. The scheme takes into account the availability of water
1569!--             in the soil layers as well as the root fraction in the
1570!--             respective layer. Layer with moisture below wilting point will
1571!--             not contribute, which reflects the preference of plants to
1572!--             take water from moister layers.
1573
1574!
1575!--             Calculate the root extraction (ECMWF 7.69, the sum of root_extr
1576!--             = 1). The energy balance solver guarantees a positive
1577!--             transpiration, so that there is no need for an additional check.
1578                DO  k = nzb_soil, nzt_soil
1579                    IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
1580                       m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i)
1581                    ENDIF
1582                ENDDO 
1583
1584                IF ( m_total > 0.0_wp )  THEN
1585                   DO  k = nzb_soil, nzt_soil
1586                      IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
1587                         root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total
1588                      ELSE
1589                         root_extr(k) = 0.0_wp
1590                      ENDIF
1591                   ENDDO
1592                ENDIF
1593
1594!
1595!--             Prognostic equation for soil water content m_soil
1596                tend(:) = 0.0_wp
1597                tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * (                  &
1598                            m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) )    &
1599                            * ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - ( &
1600                            root_extr(nzb_soil) * qsws_veg_eb(j,i)             &
1601                            + qsws_soil_eb(j,i) ) * drho_l_lv )                &
1602                            * ddz_soil_stag(nzb_soil)
1603
1604                DO  k = nzb_soil+1, nzt_soil-1
1605                   tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)             &
1606                             - m_soil(k,j,i) ) * ddz_soil(k+1) - gamma_w(k,j,i)&
1607                             - lambda_w(k-1,j,i) * (m_soil(k,j,i) -            &
1608                             m_soil(k-1,j,i)) * ddz_soil(k)                    &
1609                             + gamma_w(k-1,j,i) - (root_extr(k)                &
1610                             * qsws_veg_eb(j,i) * drho_l_lv)                   &
1611                             ) * ddz_soil_stag(k)
1612
1613                ENDDO
1614                tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i)                     &
1615                                        - lambda_w(nzt_soil-1,j,i)             &
1616                                        * (m_soil(nzt_soil,j,i)                &
1617                                        - m_soil(nzt_soil-1,j,i))              &
1618                                        * ddz_soil(nzt_soil)                   &
1619                                        + gamma_w(nzt_soil-1,j,i) - (          &
1620                                          root_extr(nzt_soil)                  &
1621                                        * qsws_veg_eb(j,i) * drho_l_lv  )      &
1622                                      ) * ddz_soil_stag(nzt_soil)             
1623
1624                m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)&
1625                                                + dt_3d * ( tsc(2) * tend(:)   &
1626                                                + tsc(3) * tm_soil_m(:,j,i) )   
1627
1628!
1629!--             Account for dry soils (find a better solution here!)
1630                m_soil_p(:,j,i) = MAX(m_soil_p(:,j,i),0.0_wp)
1631
1632!
1633!--             Calculate m_soil tendencies for the next Runge-Kutta step
1634                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1635                   IF ( intermediate_timestep_count == 1 )  THEN
1636                      DO  k = nzb_soil, nzt_soil
1637                         tm_soil_m(k,j,i) = tend(k)
1638                      ENDDO
1639                   ELSEIF ( intermediate_timestep_count <                      &
1640                            intermediate_timestep_count_max )  THEN
1641                      DO  k = nzb_soil, nzt_soil
1642                         tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp   &
1643                                     * tm_soil_m(k,j,i)
1644                      ENDDO
1645                   ENDIF
1646                ENDIF
1647
1648             ENDIF
1649
1650          ENDDO
1651       ENDDO
1652
1653!
1654!--    Calculate surface specific humidity
1655       IF ( humidity )  THEN
1656          CALL calc_q_surface
1657       ENDIF
1658
1659
1660    END SUBROUTINE lsm_soil_model
1661
1662
1663!------------------------------------------------------------------------------!
1664! Description:
1665! ------------
1666!> Calculation of specific humidity of the surface layer (surface)
1667!------------------------------------------------------------------------------!
1668    SUBROUTINE calc_q_surface
1669
1670       IMPLICIT NONE
1671
1672       INTEGER :: i              !< running index
1673       INTEGER :: j              !< running index
1674       INTEGER :: k              !< running index
1675       REAL(wp) :: resistance    !< aerodynamic and soil resistance term
1676
1677       DO  i = nxlg, nxrg   
1678          DO  j = nysg, nyng
1679             k = nzb_s_inner(j,i)
1680
1681!
1682!--          Calculate water vapour pressure at saturation
1683             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface_p(j,i)   &
1684                   - 273.16_wp ) / ( t_surface_p(j,i) - 35.86_wp ) )
1685
1686!
1687!--          Calculate specific humidity at saturation
1688             q_s = 0.622_wp * e_s / surface_pressure
1689
1690             resistance = r_a(j,i) / (r_a(j,i) + r_s(j,i))
1691
1692!
1693!--          Calculate specific humidity at surface
1694             IF ( cloud_physics )  THEN
1695                q(k,j,i) = resistance * q_s + (1.0_wp - resistance)            &
1696                             * ( q(k+1,j,i) - ql(k+1,j,i) )
1697             ELSE
1698                q(k,j,i) = resistance * q_s + (1.0_wp - resistance)            &
1699                             * q(k+1,j,i)
1700             ENDIF
1701
1702!
1703!--          Update virtual potential temperature
1704             vpt(k,j,i) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) )
1705
1706          ENDDO
1707       ENDDO
1708
1709    END SUBROUTINE calc_q_surface
1710
1711!------------------------------------------------------------------------------!
1712! Description:
1713! ------------
1714!> Swapping of timelevels
1715!------------------------------------------------------------------------------!
1716    SUBROUTINE lsm_swap_timelevel ( mod_count )
1717
1718       IMPLICIT NONE
1719
1720       INTEGER, INTENT(IN) :: mod_count
1721
1722#if defined( __nopointer )
1723
1724       t_surface    = t_surface_p
1725       t_soil       = t_soil_p
1726       IF ( humidity )  THEN
1727          m_soil    = m_soil_p
1728          m_liq_eb  = m_liq_eb_p
1729       ENDIF
1730
1731#else
1732   
1733       SELECT CASE ( mod_count )
1734
1735          CASE ( 0 )
1736
1737             t_surface  => t_surface_1; t_surface_p  => t_surface_2
1738             t_soil     => t_soil_1;    t_soil_p     => t_soil_2
1739             IF ( humidity )  THEN
1740                m_soil    => m_soil_1;   m_soil_p    => m_soil_2
1741                m_liq_eb  => m_liq_eb_1; m_liq_eb_p  => m_liq_eb_2
1742             ENDIF
1743
1744
1745          CASE ( 1 )
1746
1747             t_surface  => t_surface_2; t_surface_p  => t_surface_1
1748             t_soil     => t_soil_2;    t_soil_p     => t_soil_1
1749             IF ( humidity )  THEN
1750                m_soil    => m_soil_2;   m_soil_p    => m_soil_1
1751                m_liq_eb  => m_liq_eb_2; m_liq_eb_p  => m_liq_eb_1
1752             ENDIF
1753
1754       END SELECT
1755#endif
1756
1757    END SUBROUTINE lsm_swap_timelevel
1758
1759
1760 END MODULE land_surface_model_mod
Note: See TracBrowser for help on using the repository browser.