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

Last change on this file since 1757 was 1757, checked in by maronga, 6 years ago

some changes in land surface model, radiation model, nudging and some minor updates

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