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

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

FORTRAN an OpenMP errors removed
misplaced cpp-directive fixed
small E- and F-FORMAT changes to avoid informative compiler messages about insufficient field width

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