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

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

several bugfixes related to the new surface layer routine and land-surface-radiation interaction

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