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

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

various bugfixes and modifications of the atmosphere-land-surface-radiation interaction. Completely re-written routine to calculate surface fluxes (surface_layer_fluxes.f90) that replaces prandtl_fluxes. Minor formatting corrections and renamings

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