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

Last change on this file since 1698 was 1698, checked in by raasch, 8 years ago

last commit documented
example_cbl_rc updated

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