source: palm/trunk/SOURCE/plant_canopy_model_mod.f90 @ 3498

Last change on this file since 3498 was 3498, checked in by gronemeier, 4 years ago

Bugfix: print error message by processor which encounters the error

  • Property svn:keywords set to Id
File size: 88.2 KB
Line 
1!> @file plant_canopy_model_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2018 Leibniz Universitaet Hannover
18! Copyright 2018 Institute of Computer Science of the
19!                Czech Academy of Sciences, Prague
20!------------------------------------------------------------------------------!
21!
22! Current revisions:
23! ------------------
24!
25!
26! Former revisions:
27! -----------------
28! $Id: plant_canopy_model_mod.f90 3498 2018-11-07 10:53:03Z gronemeier $
29! corrected revisions section
30!
31! 3497 2018-11-06 17:15:18Z suehring
32! Remove write statement.
33!
34! 3467 2018-10-30 19:05:21Z suehring
35! Add calculation of transpiration for resolved plant canopy (trees, shrubs)
36! (V. Fuka, MFF UK Prague, J.Resler, ICS AS, Prague)
37!
38! Fix reading plant canopy over buildings
39!
40! 3337 2018-10-12 15:17:09Z kanani
41! Fix reading plant canopy over buildings
42!
43! 3294 2018-10-01 02:37:10Z raasch
44! ocean renamed ocean_mode
45!
46! 3274 2018-09-24 15:42:55Z knoop
47! Modularization of all bulk cloud physics code components
48!
49! 3248 2018-09-14 09:42:06Z sward
50! Minor formating changes
51!
52! 3246 2018-09-13 15:14:50Z sward
53! Added error handling for input namelist via parin_fail_message
54!
55! 3241 2018-09-12 15:02:00Z raasch
56! unused variables removed
57!
58! 3065 2018-06-12 07:03:02Z Giersch
59! dz was replaced by the help of zw to allow for vertical stretching
60!
61! 3049 2018-05-29 13:52:36Z Giersch
62! Error messages revised
63!
64! 3045 2018-05-28 07:55:41Z Giersch
65! Error message revised
66!
67! 3022 2018-05-18 11:12:35Z suehring
68! Bugfix in allocation of transpiration rate
69!
70! 3014 2018-05-09 08:42:38Z maronga
71! Bugfix: nzb_do and nzt_do were not used for 3d data output
72! Added pc_transpiration_rate
73!
74! 2977 2018-04-17 10:27:57Z kanani
75! Implement changes from branch radiation (r2948-2971) with minor modifications,
76! plus some formatting.
77! (moh.hefny):
78! Add plant canopy type to account for changes in LAD (based on the changes
79! done by Resler & Pavel) and correct the error message to PALM Standard.
80!
81! 2932 2018-03-26 09:39:22Z maronga
82! renamed canopy_par to plant_canopy_parameters
83!
84! 2920 2018-03-22 11:22:01Z kanani
85! Move usm_lad_rma and prototype_lad to radiation_model (moh.hefny)
86!
87! 2892 2018-03-14 15:06:29Z suehring
88! Bugfix, read separate ASCII LAD files for parent and child model.
89!
90! 2770 2018-01-25 15:10:09Z kanani
91! Correction of parameter check
92!
93! 2768 2018-01-24 15:38:29Z kanani
94! Added check for output quantity pcm_heatrate, some formatting
95!
96! 2766 2018-01-22 17:17:47Z kanani
97! Increased LEN of canopy mode to 30
98!
99! 2746 2018-01-15 12:06:04Z suehring
100! Move flag plant canopy to modules
101!
102! 2718 2018-01-02 08:49:38Z maronga
103! Corrected "Former revisions" section
104!
105! 2701 2017-12-15 15:40:50Z suehring
106! Changes from last commit documented
107!
108! 2698 2017-12-14 18:46:24Z suehring
109! Bugfix in get_topography_top_index
110!
111! 2696 2017-12-14 17:12:51Z kanani
112! Change in file header (GPL part)
113! Bugfix for vertical loop index pch_index in case of Netcdf input
114! Introduce 2D index array incorporate canopy top index
115! Check if canopy on top of topography do not exceed vertical dimension
116! Add check for canopy_mode in case of Netcdf input.
117! Enable _FillValue output for 3d quantities
118! Bugfix in reading of PIDS leaf area density (MS)
119!
120! 2669 2017-12-06 16:03:27Z raasch
121! coupling_char removed
122!
123! 2512 2017-10-04 08:26:59Z raasch
124! upper bounds of 3d output changed from nx+1,ny+1 to nx,ny
125! no output of ghost layer data
126!
127! 2318 2017-07-20 17:27:44Z suehring
128! Get topography top index via Function call
129!
130! 2317 2017-07-20 17:27:19Z suehring
131! Changed error messages
132!
133! 2233 2017-05-30 18:08:54Z suehring
134!
135! 2232 2017-05-30 17:47:52Z suehring
136! Adjustments to new topography concept
137!
138! 2213 2017-04-24 15:10:35Z kanani
139! Bugfix: exchange of ghost points in array pc_heating_rate needed for output
140! of pcm_heatrate, onetime ghost point exchange of lad_s after initialization.
141! Formatting and clean-up of subroutine pcm_read_plant_canopy_3d,
142! minor re-organization of canopy-heating initialization.
143!
144! 2209 2017-04-19 09:34:46Z kanani
145! Added 3d output of leaf area density (pcm_lad) and canopy
146! heat rate (pcm_heatrate)
147!
148! 2024 2016-10-12 16:42:37Z kanani
149! Added missing lad_s initialization
150!
151! 2011 2016-09-19 17:29:57Z kanani
152! Renamed canopy_heat_flux to pc_heating_rate, since the original meaning/
153! calculation of the quantity has changed, related to the urban surface model
154! and similar future applications.
155!
156! 2007 2016-08-24 15:47:17Z kanani
157! Added SUBROUTINE pcm_read_plant_canopy_3d for reading 3d plant canopy data
158! from file (new case canopy_mode=read_from_file_3d) in the course of
159! introduction of urban surface model,
160! introduced variable ext_coef,
161! resorted SUBROUTINEs to alphabetical order
162!
163! 2000 2016-08-20 18:09:15Z knoop
164! Forced header and separation lines into 80 columns
165!
166! 1960 2016-07-12 16:34:24Z suehring
167! Separate humidity and passive scalar
168!
169! 1953 2016-06-21 09:28:42Z suehring
170! Bugfix, lad_s and lad must be public
171!
172! 1826 2016-04-07 12:01:39Z maronga
173! Further modularization
174!
175! 1721 2015-11-16 12:56:48Z raasch
176! bugfixes: shf is reduced in areas covered with canopy only,
177!           canopy is set on top of topography
178!
179! 1682 2015-10-07 23:56:08Z knoop
180! Code annotations made doxygen readable
181!
182! 1484 2014-10-21 10:53:05Z kanani
183! Changes due to new module structure of the plant canopy model:
184!   module plant_canopy_model_mod now contains a subroutine for the
185!   initialization of the canopy model (pcm_init),
186!   limitation of the canopy drag (previously accounted for by calculation of
187!   a limiting canopy timestep for the determination of the maximum LES timestep
188!   in subroutine timestep) is now realized by the calculation of pre-tendencies
189!   and preliminary velocities in subroutine pcm_tendency,
190!   some redundant MPI communication removed in subroutine pcm_init
191!   (was previously in init_3d_model),
192!   unnecessary 3d-arrays lad_u, lad_v, lad_w removed - lad information on the
193!   respective grid is now provided only by lad_s (e.g. in the calculation of
194!   the tendency terms or of cum_lai_hf),
195!   drag_coefficient, lai, leaf_surface_concentration,
196!   scalar_exchange_coefficient, sec and sls renamed to canopy_drag_coeff,
197!   cum_lai_hf, leaf_surface_conc, leaf_scalar_exch_coeff, lsec and lsc,
198!   respectively,
199!   unnecessary 3d-arrays cdc, lsc and lsec now defined as single-value constants,
200!   USE-statements and ONLY-lists modified accordingly
201!
202! 1340 2014-03-25 19:45:13Z kanani
203! REAL constants defined as wp-kind
204!
205! 1320 2014-03-20 08:40:49Z raasch
206! ONLY-attribute added to USE-statements,
207! kind-parameters added to all INTEGER and REAL declaration statements,
208! kinds are defined in new module kinds,
209! old module precision_kind is removed,
210! revision history before 2012 removed,
211! comment fields (!:) to be used for variable explanations added to
212! all variable declaration statements
213!
214! 1036 2012-10-22 13:43:42Z raasch
215! code put under GPL (PALM 3.9)
216!
217! 138 2007-11-28 10:03:58Z letzel
218! Initial revision
219!
220! Description:
221! ------------
222!> 1) Initialization of the canopy model, e.g. construction of leaf area density
223!> profile (subroutine pcm_init).
224!> 2) Calculation of sinks and sources of momentum, heat and scalar concentration
225!> due to canopy elements (subroutine pcm_tendency).
226!------------------------------------------------------------------------------!
227 MODULE plant_canopy_model_mod
228 
229    USE arrays_3d,                                                             &
230        ONLY:  dzu, dzw, e, exner, hyp, pt, q, s, tend, u, v, w, zu, zw
231
232    USE basic_constants_and_equations_mod,                                     &
233        ONLY:  c_p, degc_to_k, l_v, lv_d_cp, r_d, rd_d_rv
234
235    USE control_parameters,                                                   &
236        ONLY: humidity
237
238    USE indices,                                                               &
239        ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
240               nz, nzb, nzt
241
242    USE kinds
243
244    USE pegrid
245
246    USE surface_mod,                                                           &
247        ONLY:  get_topography_top_index_ji
248
249
250    IMPLICIT NONE
251
252
253    CHARACTER (LEN=30)   ::  canopy_mode = 'block'                 !< canopy coverage
254    LOGICAL              ::  plant_canopy_transpiration = .FALSE.  !< flag to switch calculation of transpiration and corresponding latent heat
255                                                                   !< for resolved plant canopy inside radiation model
256                                                                   !< (calls subroutine pcm_calc_transpiration_rate from module plant_canopy_mod)
257
258    INTEGER(iwp) ::  pch_index = 0                                 !< plant canopy height/top index
259    INTEGER(iwp) ::  lad_vertical_gradient_level_ind(10) = -9999   !< lad-profile levels (index)
260
261    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  pch_index_ji     !< local plant canopy top
262
263    LOGICAL ::  calc_beta_lad_profile = .FALSE.   !< switch for calc. of lad from beta func.
264
265    REAL(wp) ::  alpha_lad = 9999999.9_wp         !< coefficient for lad calculation
266    REAL(wp) ::  beta_lad = 9999999.9_wp          !< coefficient for lad calculation
267    REAL(wp) ::  canopy_drag_coeff = 0.0_wp       !< canopy drag coefficient (parameter)
268    REAL(wp) ::  cdc = 0.0_wp                     !< canopy drag coeff. (abbreviation used in equations)
269    REAL(wp) ::  cthf = 0.0_wp                    !< canopy top heat flux
270    REAL(wp) ::  dt_plant_canopy = 0.0_wp         !< timestep account. for canopy drag
271    REAL(wp) ::  ext_coef = 0.6_wp                !< extinction coefficient
272    REAL(wp) ::  lad_surface = 0.0_wp             !< lad surface value
273    REAL(wp) ::  lai_beta = 0.0_wp                !< leaf area index (lai) for lad calc.
274    REAL(wp) ::  leaf_scalar_exch_coeff = 0.0_wp  !< canopy scalar exchange coeff.
275    REAL(wp) ::  leaf_surface_conc = 0.0_wp       !< leaf surface concentration
276    REAL(wp) ::  lsc = 0.0_wp                     !< leaf surface concentration
277    REAL(wp) ::  lsec = 0.0_wp                    !< leaf scalar exchange coeff.
278
279    REAL(wp) ::  lad_vertical_gradient(10) = 0.0_wp              !< lad gradient
280    REAL(wp) ::  lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m)
281
282    REAL(wp) ::  lad_type_coef(0:10) = 1.0_wp   !< multiplicative coeficients for particular types
283                                                !< of plant canopy (e.g. deciduous tree during winter)
284
285    REAL(wp), DIMENSION(:), ALLOCATABLE ::  lad            !< leaf area density
286    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pre_lad        !< preliminary lad
287   
288    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  cum_lai_hf             !< cumulative lai for heatflux calc.
289    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  lad_s                  !< lad on scalar-grid
290    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pc_heating_rate        !< plant canopy heating rate
291    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pc_transpiration_rate  !< plant canopy transpiration rate
292    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pc_latent_rate         !< plant canopy latent heating rate
293
294    SAVE
295
296
297    PRIVATE
298 
299!
300!-- Public functions
301    PUBLIC pcm_calc_transpiration_rate, pcm_check_data_output,                &
302           pcm_check_parameters, pcm_data_output_3d, pcm_define_netcdf_grid,  &
303           pcm_header, pcm_init, pcm_parin, pcm_tendency
304
305!
306!-- Public variables and constants
307    PUBLIC cdc, pc_heating_rate, pc_transpiration_rate, pc_latent_rate,       &
308           canopy_mode, cthf, dt_plant_canopy, lad, lad_s, pch_index,         &
309           plant_canopy_transpiration
310
311    INTERFACE pcm_calc_transpiration_rate
312       MODULE PROCEDURE pcm_calc_transpiration_rate
313    END INTERFACE pcm_calc_transpiration_rate
314
315    INTERFACE pcm_check_data_output
316       MODULE PROCEDURE pcm_check_data_output
317    END INTERFACE pcm_check_data_output
318   
319    INTERFACE pcm_check_parameters
320       MODULE PROCEDURE pcm_check_parameters
321    END INTERFACE pcm_check_parameters
322
323    INTERFACE pcm_data_output_3d
324       MODULE PROCEDURE pcm_data_output_3d
325    END INTERFACE pcm_data_output_3d
326
327    INTERFACE pcm_define_netcdf_grid
328       MODULE PROCEDURE pcm_define_netcdf_grid
329    END INTERFACE pcm_define_netcdf_grid
330   
331     INTERFACE pcm_header
332       MODULE PROCEDURE pcm_header
333    END INTERFACE pcm_header       
334   
335    INTERFACE pcm_init
336       MODULE PROCEDURE pcm_init
337    END INTERFACE pcm_init
338
339    INTERFACE pcm_parin
340       MODULE PROCEDURE pcm_parin
341    END INTERFACE pcm_parin
342
343    INTERFACE pcm_read_plant_canopy_3d
344       MODULE PROCEDURE pcm_read_plant_canopy_3d
345    END INTERFACE pcm_read_plant_canopy_3d
346   
347    INTERFACE pcm_tendency
348       MODULE PROCEDURE pcm_tendency
349       MODULE PROCEDURE pcm_tendency_ij
350    END INTERFACE pcm_tendency
351
352
353 CONTAINS
354
355
356
357!------------------------------------------------------------------------------!
358! Description:
359! ------------
360!> Calculation of the plant canopy transpiration rate based on the Jarvis-Stewart
361!> with parametrizations described in Daudet et al. (1999; Agricult. and Forest
362!> Meteorol. 97) and Ngao, Adam and Saudreau (2017;  Agricult. and Forest Meteorol
363!> 237-238). Model functions f1-f4 were adapted from Stewart (1998; Agric.
364!> and Forest. Meteorol. 43) instead, because they are valid for broader intervals
365!> of values. Funcion f4 used in form present in van Wijk et al. (1998;
366!> Tree Physiology 20).
367!>
368!> This subroutine is called from subroutine radiation_interaction
369!> after the calculation of radiation in plant canopy boxes.
370!> (arrays pcbinsw and pcbinlw).
371!>
372!------------------------------------------------------------------------------!
373 SUBROUTINE pcm_calc_transpiration_rate(i, j, k, kk, pcbsw, pcblw, pcbtr, pcblh)
374
375     USE control_parameters,                                                   &
376        ONLY: dz
377
378     USE grid_variables,                                                       &
379        ONLY:  dx, dy
380
381     IMPLICIT NONE
382!--  input parameters
383     INTEGER(iwp), INTENT(IN)          :: i, j, k, kk        !< indices of the pc gridbox
384     REAL(wp), INTENT(IN)              :: pcbsw              !< sw radiation in gridbox (W)
385     REAL(wp), INTENT(IN)              :: pcblw              !< lw radiation in gridbox (W)
386     REAL(wp), INTENT(OUT)             :: pcbtr              !< transpiration rate dq/dt (kg/kg/s)
387     REAL(wp), INTENT(OUT)             :: pcblh              !< latent heat from transpiration dT/dt (K/s)
388
389!--  variables and parameters for calculation of transpiration rate
390     REAL(wp)                          :: sat_press, sat_press_d, temp, v_lad
391     REAL(wp)                          :: d_fact, g_b, g_s, wind_speed, evapor_rate
392     REAL(wp)                          :: f1, f2, f3, f4, vpd, rswc, e_eq, e_imp, rad
393     REAL(wp), PARAMETER               :: gama_psychr = 66   !< psychrometric constant (Pa/K)
394     REAL(wp), PARAMETER               :: g_s_max = 0.01     !< maximum stomatal conductivity (m/s)
395     REAL(wp), PARAMETER               :: m_soil = 0.4_wp    !< soil water content (needs to adjust or take from LSM)
396     REAL(wp), PARAMETER               :: m_wilt = 0.01_wp   !< wilting point soil water content (needs to adjust or take from LSM)
397     REAL(wp), PARAMETER               :: m_sat = 0.51_wp    !< saturation soil water content (needs to adjust or take from LSM)
398     REAL(wp), PARAMETER               :: t2_min = 0.0_wp    !< minimal temperature for calculation of f2
399     REAL(wp), PARAMETER               :: t2_max = 40.0_wp   !< maximal temperature for calculation of f2
400
401
402!--  Temperature (deg C)
403     temp = pt(k,j,i) * exner(k) - degc_to_k
404!--  Coefficient for conversion of radiation to grid to radiation to unit leaves surface
405     v_lad = 1.0_wp / ( MAX(lad_s(kk,j,i), 1.0e-10) * dx * dy * dz(1) )
406!--  Magnus formula for the saturation pressure (see Ngao, Adam and Saudreau (2017) eq. 1)
407!--  There are updated formulas available, kept consistent with the rest of the parametrization
408     sat_press = 610.8_wp * exp(17.27_wp * temp/(temp + 237.3_wp))
409!--  Saturation pressure derivative (derivative of the above)
410     sat_press_d = sat_press * 17.27_wp * 237.3_wp / (temp + 237.3_wp)**2
411!--  Wind speed
412     wind_speed = SQRT( ((u(k,j,i) + u(k,j,i+1))/2.0_wp )**2 +           &
413                        ((v(k,j,i) + v(k,j,i+1))/2.0_wp )**2 +           &
414                        ((w(k,j,i) + w(k,j,i+1))/2.0_wp )**2 )
415!--  Aerodynamic conductivity (Daudet et al. (1999) eq. 14
416     g_b = 0.01_wp * wind_speed + 0.0071_wp
417!--  Radiation flux per leaf surface unit
418     rad = pcbsw * v_lad
419!--  First function for calculation of stomatal conductivity (radiation dependency)
420!--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 17
421     f1 = rad * (1000._wp+42.1_wp) / 1000._wp / (rad+42.1_wp)
422!--  Second function for calculation of stomatal conductivity (temperature dependency)
423!--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 21
424     f2 = MAX(t2_min, (temp-t2_min) * (t2_max-temp)**((t2_max-16.9_wp)/(16.9_wp-t2_min)) / &
425              ((16.9_wp-t2_min) * (t2_max-16.9_wp)**((t2_max-16.9_wp)/(16.9_wp-t2_min))) )
426!--  Water pressure deficit
427!--  Ngao, Adam and Saudreau (2017) eq. 6 but with water vapour partial pressure
428     vpd = max( sat_press - q(k,j,i) * hyp(k) / rd_d_rv, 0._wp )
429!--  Third function for calculation of stomatal conductivity (water pressure deficit dependency)
430!--  Ngao, Adam and Saudreau (2017) Table 1, limited from below according to Stewart (1988)
431!--  The coefficients of the linear dependence should better correspond to broad-leaved trees
432!--  than the coefficients from Stewart (1988) which correspond to conifer trees.
433     vpd = MIN(MAX(vpd,770.0_wp),3820.0_wp)
434     f3 = -2e-4_wp * vpd + 1.154_wp
435!--  Fourth function for calculation of stomatal conductivity (soil moisture dependency)
436!--  Residual soil water content
437!--  van Wijk et al. (1998; Tree Physiology 20) eq. 7
438!--  TODO - over LSM surface might be calculated from LSM parameters
439     rswc = ( m_sat - m_soil ) / ( m_sat - m_wilt )
440!--  van Wijk et al. (1998; Tree Physiology 20) eq. 5-6 (it is a reformulation of eq. 22-23 of Stewart(1988))
441     f4 = MAX(0._wp, MIN(1.0_wp - 0.041_wp * EXP(3.2_wp * rswc), 1.0_wp - 0.041_wp))
442!--  Stomatal conductivity
443!--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 12
444!--  (notation according to Ngao, Adam and Saudreau (2017) and others)
445     g_s = g_s_max * f1 * f2 * f3 * f4 + 1.0e-10_wp
446!--  Decoupling factor
447!--  Daudet et al. (1999) eq. 6
448     d_fact = (sat_press_d / gama_psychr + 2._wp ) /                        &
449              (sat_press_d / gama_psychr + 2._wp + 2 * g_b / g_s )
450!--  Equilibrium evaporation rate
451!--  Daudet et al. (1999) eq. 4
452     e_eq = (pcbsw + pcblw) * v_lad * sat_press_d /                         &
453                 gama_psychr /( sat_press_d / gama_psychr + 2.0_wp ) / l_v
454!--  Imposed evaporation rate
455!--  Daudet et al. (1999) eq. 5
456     e_imp = r_d * pt(k,j,i) * exner(k) / hyp(k) * c_p * g_s * vpd / gama_psychr / l_v
457!--  Evaporation rate
458!--  Daudet et al. (1999) eq. 3
459!--  (evaporation rate is limited to non-negative values)
460     evapor_rate = MAX(d_fact * e_eq + ( 1.0_wp - d_fact ) * e_imp, 0.0_wp)
461!--  Conversion of evaporation rate to q tendency in gridbox
462!--  dq/dt = E * LAD * V_g / (rho_air * V_g)
463     pcbtr = evapor_rate * r_d * pt(k,j,i) * exner(k) * lad_s(kk,j,i) / hyp(k)  !-- = dq/dt
464!--  latent heat from evaporation
465     pcblh = pcbtr * lv_d_cp  !-- = - dT/dt
466
467 END SUBROUTINE pcm_calc_transpiration_rate
468
469
470!------------------------------------------------------------------------------!
471! Description:
472! ------------
473!> Check data output for plant canopy model
474!------------------------------------------------------------------------------!
475 SUBROUTINE pcm_check_data_output( var, unit )
476 
477 
478    USE control_parameters,                                                 &
479        ONLY:  message_string, urban_surface
480
481    IMPLICIT NONE
482
483    CHARACTER (LEN=*) ::  unit  !<
484    CHARACTER (LEN=*) ::  var   !<
485
486
487    SELECT CASE ( TRIM( var ) )
488
489       CASE ( 'pcm_heatrate' )
490          IF ( cthf == 0.0_wp  .AND. .NOT.  urban_surface )  THEN
491             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
492                              'res setting of parameter cthf /= 0.0'
493             CALL message( 'pcm_check_data_output', 'PA1000', 1, 2, 0, 6, 0 )
494          ENDIF
495          unit = 'K s-1'
496   
497       CASE ( 'pcm_transpirationrate' )
498          unit = 'kg kg-1 s-1'
499
500       CASE ( 'pcm_latentrate' )
501          unit = 'K s-1'
502
503       CASE ( 'pcm_bowenratio' )
504          unit = 'K s-1'
505
506       CASE ( 'pcm_lad' )
507          unit = 'm2 m-3'
508
509
510       CASE DEFAULT
511          unit = 'illegal'
512
513    END SELECT
514
515
516 END SUBROUTINE pcm_check_data_output
517 
518 
519!------------------------------------------------------------------------------!
520! Description:
521! ------------
522!> Check parameters routine for plant canopy model
523!------------------------------------------------------------------------------!
524    SUBROUTINE pcm_check_parameters
525
526       USE control_parameters,                                                 &
527           ONLY: coupling_char, message_string
528
529       USE bulk_cloud_model_mod,                                               &
530           ONLY: bulk_cloud_model, microphysics_seifert
531
532       USE netcdf_data_input_mod,                                              &
533           ONLY:  input_file_static, input_pids_static
534                 
535   
536       IMPLICIT NONE
537
538   
539       IF ( canopy_drag_coeff == 0.0_wp )  THEN
540          message_string = 'plant_canopy = .TRUE. requires a non-zero drag '// &
541                           'coefficient & given value is canopy_drag_coeff = 0.0'
542          CALL message( 'pcm_check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
543       ENDIF
544   
545       IF ( ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad == 9999999.9_wp ) .OR.&
546              beta_lad /= 9999999.9_wp   .AND.  alpha_lad == 9999999.9_wp )  THEN
547          message_string = 'using the beta function for the construction ' //  &
548                           'of the leaf area density profile requires '    //  &
549                           'both alpha_lad and beta_lad to be /= 9999999.9'
550          CALL message( 'pcm_check_parameters', 'PA0118', 1, 2, 0, 6, 0 )
551       ENDIF
552   
553       IF ( calc_beta_lad_profile  .AND.  lai_beta == 0.0_wp )  THEN
554          message_string = 'using the beta function for the construction ' //  &
555                           'of the leaf area density profile requires '    //  &
556                           'a non-zero lai_beta, but given value is '      //  &
557                           'lai_beta = 0.0'
558          CALL message( 'pcm_check_parameters', 'PA0119', 1, 2, 0, 6, 0 )
559       ENDIF
560
561       IF ( calc_beta_lad_profile  .AND.  lad_surface /= 0.0_wp )  THEN
562          message_string = 'simultaneous setting of alpha_lad /= 9999999.9 '// &
563                           'combined with beta_lad /= 9999999.9 '           // &
564                           'and lad_surface /= 0.0 is not possible, '       // &
565                           'use either vertical gradients or the beta '     // &
566                           'function for the construction of the leaf area '// &
567                           'density profile'
568          CALL message( 'pcm_check_parameters', 'PA0120', 1, 2, 0, 6, 0 )
569       ENDIF
570
571       IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
572          message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' // &
573                          ' seifert_beheng'
574          CALL message( 'pcm_check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
575       ENDIF
576!
577!--    If dynamic input file is used, canopy need to be read from file
578       IF ( input_pids_static  .AND.                                           &
579            TRIM( canopy_mode ) /= 'read_from_file_3d' )  THEN
580          message_string = 'Usage of dynamic input file ' //                   &
581                           TRIM( input_file_static ) //                        &
582                           TRIM( coupling_char ) // ' requires ' //            &
583                           'canopy_mode = read_from_file_3d'
584          CALL message( 'pcm_check_parameters', 'PA0999', 1, 2, 0, 6, 0 )
585       ENDIF
586
587 
588    END SUBROUTINE pcm_check_parameters 
589 
590
591!------------------------------------------------------------------------------!
592!
593! Description:
594! ------------
595!> Subroutine defining 3D output variables
596!------------------------------------------------------------------------------!
597 SUBROUTINE pcm_data_output_3d( av, variable, found, local_pf, fill_value,     &
598                                nzb_do, nzt_do )
599 
600    USE indices
601
602    USE kinds
603
604
605    IMPLICIT NONE
606
607    CHARACTER (LEN=*) ::  variable !<
608
609    INTEGER(iwp) ::  av     !<
610    INTEGER(iwp) ::  i      !<
611    INTEGER(iwp) ::  j      !<
612    INTEGER(iwp) ::  k      !<
613    INTEGER(iwp) ::  k_topo !< topography top index
614    INTEGER(iwp) ::  nzb_do !< lower limit of the data output (usually 0)
615    INTEGER(iwp) ::  nzt_do !< vertical upper limit of the data output (usually nz_do3d)
616
617    LOGICAL      ::  found !<
618
619    REAL(wp)     ::  fill_value
620    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
621
622
623    found = .TRUE.
624
625    local_pf = REAL( fill_value, KIND = 4 )
626
627    SELECT CASE ( TRIM( variable ) )
628
629      CASE ( 'pcm_heatrate' )
630         IF ( av == 0 )  THEN
631            DO  i = nxl, nxr
632               DO  j = nys, nyn
633                  IF ( pch_index_ji(j,i) /= 0 )  THEN
634                     k_topo = get_topography_top_index_ji( j, i, 's' )
635                     DO  k = k_topo, k_topo + pch_index_ji(j,i)
636                        local_pf(i,j,k) = pc_heating_rate(k-k_topo,j,i)
637                     ENDDO
638                  ENDIF
639               ENDDO
640            ENDDO
641         ENDIF
642   
643       CASE ( 'pcm_transpirationrate' )
644         IF ( av == 0 )  THEN
645            DO  i = nxl, nxr
646               DO  j = nys, nyn
647                  IF ( pch_index_ji(j,i) /= 0 )  THEN
648                     k_topo = get_topography_top_index_ji( j, i, 's' )
649                     DO  k = k_topo, k_topo + pch_index_ji(j,i)
650                        local_pf(i,j,k) = pc_transpiration_rate(k-k_topo,j,i)
651                     ENDDO
652                  ENDIF
653               ENDDO
654            ENDDO
655         ENDIF
656
657       CASE ( 'pcm_latentrate' )
658         IF ( av == 0 )  THEN
659            DO  i = nxl, nxr
660               DO  j = nys, nyn
661                  IF ( pch_index_ji(j,i) /= 0 )  THEN
662                     k_topo = get_topography_top_index_ji( j, i, 's' )
663                     DO  k = k_topo, k_topo + pch_index_ji(j,i)
664                        local_pf(i,j,k) = pc_latent_rate(k-k_topo,j,i)
665                     ENDDO
666                  ENDIF
667               ENDDO
668            ENDDO
669         ENDIF
670
671       CASE ( 'pcm_bowenratio' )
672         IF ( av == 0 )  THEN
673            DO  i = nxl, nxr
674               DO  j = nys, nyn
675                  IF ( pch_index_ji(j,i) /= 0 )  THEN
676                     k_topo = get_topography_top_index_ji( j, i, 's' )
677                     DO  k = k_topo, k_topo + pch_index_ji(j,i)
678                        IF ( pc_latent_rate(k-k_topo,j,i) /= 0._wp ) THEN
679                           local_pf(i,j,k) = pc_heating_rate(k-k_topo,j,i) / &
680                                             pc_latent_rate(k-k_topo,j,i)
681                        ENDIF
682                     ENDDO
683                  ENDIF
684               ENDDO
685            ENDDO
686         ENDIF
687
688      CASE ( 'pcm_lad' )
689         IF ( av == 0 )  THEN
690            DO  i = nxl, nxr
691               DO  j = nys, nyn
692                  IF ( pch_index_ji(j,i) /= 0 )  THEN
693                     k_topo = get_topography_top_index_ji( j, i, 's' )
694                     DO  k = k_topo, k_topo + pch_index_ji(j,i)
695                        local_pf(i,j,k) = lad_s(k-k_topo,j,i)
696                     ENDDO
697                  ENDIF
698               ENDDO
699            ENDDO
700         ENDIF
701                 
702         
703       CASE DEFAULT
704          found = .FALSE.
705
706    END SELECT
707
708
709 END SUBROUTINE pcm_data_output_3d 
710         
711!------------------------------------------------------------------------------!
712!
713! Description:
714! ------------
715!> Subroutine defining appropriate grid for netcdf variables.
716!> It is called from subroutine netcdf.
717!------------------------------------------------------------------------------!
718 SUBROUTINE pcm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
719   
720     IMPLICIT NONE
721
722     CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
723     LOGICAL, INTENT(OUT)           ::  found       !<
724     CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
725     CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
726     CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
727
728     found  = .TRUE.
729
730!
731!--  Check for the grid
732     SELECT CASE ( TRIM( var ) )
733
734        CASE ( 'pcm_heatrate', 'pcm_lad', 'pcm_transpirationrate', 'pcm_latentrate', 'pcm_bowenratio')
735           grid_x = 'x'
736           grid_y = 'y'
737           grid_z = 'zu'
738
739        CASE DEFAULT
740           found  = .FALSE.
741           grid_x = 'none'
742           grid_y = 'none'
743           grid_z = 'none'
744     END SELECT
745
746 END SUBROUTINE pcm_define_netcdf_grid
747 
748 
749!------------------------------------------------------------------------------!
750! Description:
751! ------------
752!> Header output for plant canopy model
753!------------------------------------------------------------------------------!
754    SUBROUTINE pcm_header ( io )
755
756       USE control_parameters,                                                 &
757           ONLY: passive_scalar
758
759
760       IMPLICIT NONE
761 
762       CHARACTER (LEN=10) ::  coor_chr            !<
763
764       CHARACTER (LEN=86) ::  coordinates         !<
765       CHARACTER (LEN=86) ::  gradients           !<
766       CHARACTER (LEN=86) ::  leaf_area_density   !<
767       CHARACTER (LEN=86) ::  slices              !<
768 
769       INTEGER(iwp) :: i                !<
770       INTEGER(iwp),  INTENT(IN) ::  io !< Unit of the output file
771       INTEGER(iwp) :: k                !<       
772   
773       REAL(wp) ::  canopy_height       !< canopy height (in m)
774       
775       canopy_height = zw(pch_index)
776
777       WRITE ( io, 1 )  canopy_mode, canopy_height, pch_index,                 &
778                          canopy_drag_coeff
779       IF ( passive_scalar )  THEN
780          WRITE ( io, 2 )  leaf_scalar_exch_coeff,                             &
781                             leaf_surface_conc
782       ENDIF
783
784!
785!--    Heat flux at the top of vegetation
786       WRITE ( io, 3 )  cthf
787
788!
789!--    Leaf area density profile, calculated either from given vertical
790!--    gradients or from beta probability density function.
791       IF (  .NOT.  calc_beta_lad_profile )  THEN
792
793!--       Building output strings, starting with surface value
794          WRITE ( leaf_area_density, '(F7.4)' )  lad_surface
795          gradients = '------'
796          slices = '     0'
797          coordinates = '   0.0'
798          i = 1
799          DO  WHILE ( i < 11  .AND.  lad_vertical_gradient_level_ind(i)        &
800                      /= -9999 )
801
802             WRITE (coor_chr,'(F7.2)')  lad(lad_vertical_gradient_level_ind(i))
803             leaf_area_density = TRIM( leaf_area_density ) // ' ' //           &
804                                 TRIM( coor_chr )
805 
806             WRITE (coor_chr,'(F7.2)')  lad_vertical_gradient(i)
807             gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
808
809             WRITE (coor_chr,'(I7)')  lad_vertical_gradient_level_ind(i)
810             slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
811
812             WRITE (coor_chr,'(F7.1)')  lad_vertical_gradient_level(i)
813             coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
814
815             i = i + 1
816          ENDDO
817
818          WRITE ( io, 4 )  TRIM( coordinates ), TRIM( leaf_area_density ),     &
819                             TRIM( gradients ), TRIM( slices )
820
821       ELSE
822       
823          WRITE ( leaf_area_density, '(F7.4)' )  lad_surface
824          coordinates = '   0.0'
825         
826          DO  k = 1, pch_index
827
828             WRITE (coor_chr,'(F7.2)')  lad(k)
829             leaf_area_density = TRIM( leaf_area_density ) // ' ' //           &
830                                 TRIM( coor_chr )
831 
832             WRITE (coor_chr,'(F7.1)')  zu(k)
833             coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
834
835          ENDDO       
836
837          WRITE ( io, 5 ) TRIM( coordinates ), TRIM( leaf_area_density ),      &
838                          alpha_lad, beta_lad, lai_beta
839
840       ENDIF 
841
8421 FORMAT (//' Vegetation canopy (drag) model:'/                                &
843              ' ------------------------------'//                              &
844              ' Canopy mode: ', A /                                            &
845              ' Canopy height: ',F6.2,'m (',I4,' grid points)' /               &
846              ' Leaf drag coefficient: ',F6.2 /)
8472 FORMAT (/ ' Scalar exchange coefficient: ',F6.2 /                            &
848              ' Scalar concentration at leaf surfaces in kg/m**3: ',F6.2 /)
8493 FORMAT (' Predefined constant heatflux at the top of the vegetation: ',F6.2, &
850          ' K m/s')
8514 FORMAT (/ ' Characteristic levels of the leaf area density:'//               &
852              ' Height:              ',A,'  m'/                                &
853              ' Leaf area density:   ',A,'  m**2/m**3'/                        &
854              ' Gradient:            ',A,'  m**2/m**4'/                        &
855              ' Gridpoint:           ',A)
8565 FORMAT (//' Characteristic levels of the leaf area density and coefficients:'&
857          //  ' Height:              ',A,'  m'/                                &
858              ' Leaf area density:   ',A,'  m**2/m**3'/                        &
859              ' Coefficient alpha: ',F6.2 /                                    &
860              ' Coefficient beta: ',F6.2 /                                     &
861              ' Leaf area index: ',F6.2,'  m**2/m**2' /)   
862       
863    END SUBROUTINE pcm_header
864 
865 
866!------------------------------------------------------------------------------!
867! Description:
868! ------------
869!> Initialization of the plant canopy model
870!------------------------------------------------------------------------------!
871    SUBROUTINE pcm_init
872   
873
874       USE control_parameters,                                                 &
875           ONLY: message_string, ocean_mode, urban_surface
876
877       USE netcdf_data_input_mod,                                              &
878           ONLY:  leaf_area_density_f
879
880       USE surface_mod,                                                        &
881           ONLY: surf_def_h, surf_lsm_h, surf_usm_h
882
883       IMPLICIT NONE
884
885       INTEGER(iwp) ::  i   !< running index
886       INTEGER(iwp) ::  j   !< running index
887       INTEGER(iwp) ::  k   !< running index
888       INTEGER(iwp) ::  m   !< running index
889       INTEGER(iwp) ::  pch_index_l
890
891       REAL(wp) ::  int_bpdf        !< vertical integral for lad-profile construction
892       REAL(wp) ::  gradient        !< gradient for lad-profile construction
893       REAL(wp) ::  canopy_height   !< canopy height for lad-profile construction
894
895!
896!--    Allocate one-dimensional arrays for the computation of the
897!--    leaf area density (lad) profile
898       ALLOCATE( lad(0:nz+1), pre_lad(0:nz+1) )
899       lad = 0.0_wp
900       pre_lad = 0.0_wp
901
902!
903!--    Set flag that indicates that the lad-profile shall be calculated by using
904!--    a beta probability density function
905       IF ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad /= 9999999.9_wp )  THEN
906          calc_beta_lad_profile = .TRUE.
907       ENDIF
908       
909       
910!
911!--    Compute the profile of leaf area density used in the plant
912!--    canopy model. The profile can either be constructed from
913!--    prescribed vertical gradients of the leaf area density or by
914!--    using a beta probability density function (see e.g. Markkanen et al.,
915!--    2003: Boundary-Layer Meteorology, 106, 437-459)
916       IF (  .NOT.  calc_beta_lad_profile )  THEN   
917
918!
919!--       Use vertical gradients for lad-profile construction   
920          i = 1
921          gradient = 0.0_wp
922
923          IF (  .NOT.  ocean_mode )  THEN
924
925             lad(0) = lad_surface
926             lad_vertical_gradient_level_ind(1) = 0
927 
928             DO k = 1, pch_index
929                IF ( i < 11 )  THEN
930                   IF ( lad_vertical_gradient_level(i) < zu(k)  .AND.          &
931                        lad_vertical_gradient_level(i) >= 0.0_wp )  THEN
932                      gradient = lad_vertical_gradient(i)
933                      lad_vertical_gradient_level_ind(i) = k - 1
934                      i = i + 1
935                   ENDIF
936                ENDIF
937                IF ( gradient /= 0.0_wp )  THEN
938                   IF ( k /= 1 )  THEN
939                      lad(k) = lad(k-1) + dzu(k) * gradient
940                   ELSE
941                      lad(k) = lad_surface + dzu(k) * gradient
942                   ENDIF
943                ELSE
944                   lad(k) = lad(k-1)
945                ENDIF
946             ENDDO
947
948          ENDIF
949
950!
951!--       In case of no given leaf area density gradients, choose a vanishing
952!--       gradient. This information is used for the HEADER and the RUN_CONTROL
953!--       file.
954          IF ( lad_vertical_gradient_level(1) == -9999999.9_wp )  THEN
955             lad_vertical_gradient_level(1) = 0.0_wp
956          ENDIF
957
958       ELSE
959
960!
961!--       Use beta function for lad-profile construction
962          int_bpdf = 0.0_wp
963          canopy_height = zw(pch_index)
964
965          DO k = 0, pch_index
966             int_bpdf = int_bpdf +                                             &
967                      ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) *  &
968                      ( ( 1.0_wp - ( zw(k) / canopy_height ) )**(              &
969                          beta_lad-1.0_wp ) )                                  &
970                      * ( ( zw(k+1)-zw(k) ) / canopy_height ) )
971          ENDDO
972
973!
974!--       Preliminary lad profile (defined on w-grid)
975          DO k = 0, pch_index
976             pre_lad(k) =  lai_beta *                                          &
977                        ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) )  &
978                        * ( ( 1.0_wp - ( zw(k) / canopy_height ) )**(          &
979                              beta_lad-1.0_wp ) ) / int_bpdf                   &
980                        ) / canopy_height
981          ENDDO
982
983!
984!--       Final lad profile (defined on scalar-grid level, since most prognostic
985!--       quantities are defined there, hence, less interpolation is required
986!--       when calculating the canopy tendencies)
987          lad(0) = pre_lad(0)
988          DO k = 1, pch_index
989             lad(k) = 0.5 * ( pre_lad(k-1) + pre_lad(k) )
990          ENDDO         
991
992       ENDIF
993
994!
995!--    Allocate 3D-array for the leaf area density (lad_s).
996       ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
997
998!
999!--    Initialize canopy parameters cdc (canopy drag coefficient),
1000!--    lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration)
1001!--    with the prescribed values
1002       cdc = canopy_drag_coeff
1003       lsec = leaf_scalar_exch_coeff
1004       lsc = leaf_surface_conc
1005
1006!
1007!--    Initialization of the canopy coverage in the model domain:
1008!--    Setting the parameter canopy_mode = 'block' initializes a canopy, which
1009!--    fully covers the domain surface
1010       SELECT CASE ( TRIM( canopy_mode ) )
1011
1012          CASE( 'block' )
1013
1014             DO  i = nxlg, nxrg
1015                DO  j = nysg, nyng
1016                   lad_s(:,j,i) = lad(:)
1017                ENDDO
1018             ENDDO
1019
1020          CASE ( 'read_from_file_3d' )
1021!
1022!--          Initialize LAD with data from file. If LAD is given in NetCDF file,
1023!--          use these values, else take LAD profiles from ASCII file.
1024!--          Please note, in NetCDF file LAD is only given up to the maximum
1025!--          canopy top, indicated by leaf_area_density_f%nz. 
1026             lad_s = 0.0_wp
1027             IF ( leaf_area_density_f%from_file )  THEN
1028!
1029!--             Set also pch_index, used to be the upper bound of the vertical
1030!--             loops. Therefore, use the global top of the canopy layer.
1031                pch_index = leaf_area_density_f%nz - 1
1032
1033                DO  i = nxl, nxr
1034                   DO  j = nys, nyn
1035                      DO  k = 0, leaf_area_density_f%nz - 1
1036                      IF ( leaf_area_density_f%var(k,j,i) /=                   &
1037                           leaf_area_density_f%fill )                          &
1038                         lad_s(k,j,i) = leaf_area_density_f%var(k,j,i)
1039                      ENDDO
1040                   ENDDO
1041                ENDDO
1042                CALL exchange_horiz( lad_s, nbgp )
1043!
1044!            ASCII file
1045!--          Initialize canopy parameters cdc (canopy drag coefficient),
1046!--          lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration)
1047!--          from file which contains complete 3D data (separate vertical profiles for
1048!--          each location).
1049             ELSE
1050                CALL pcm_read_plant_canopy_3d
1051             ENDIF
1052
1053          CASE DEFAULT
1054!
1055!--          The DEFAULT case is reached either if the parameter
1056!--          canopy mode contains a wrong character string or if the
1057!--          user has coded a special case in the user interface.
1058!--          There, the subroutine user_init_plant_canopy checks
1059!--          which of these two conditions applies.
1060             CALL user_init_plant_canopy
1061 
1062       END SELECT
1063!
1064!--    Initialize 2D index array indicating canopy top index.
1065       ALLOCATE( pch_index_ji(nysg:nyng,nxlg:nxrg) )
1066       pch_index_ji = 0
1067
1068       DO  i = nxl, nxr
1069          DO  j = nys, nyn
1070             DO  k = 0, pch_index
1071                IF ( lad_s(k,j,i) /= 0 )  pch_index_ji(j,i) = k
1072             ENDDO
1073!
1074!--          Check whether topography and local vegetation on top exceed
1075!--          height of the model domain.
1076             k = get_topography_top_index_ji( j, i, 's' )
1077             IF ( k + pch_index_ji(j,i) >= nzt + 1 )  THEN
1078                message_string =  'Local vegetation height on top of ' //      &
1079                                  'topography exceeds height of model domain.'
1080                CALL message( 'pcm_init', 'PA0999', 2, 2, 0, 6, 0 )
1081             ENDIF
1082
1083          ENDDO
1084       ENDDO
1085
1086       CALL exchange_horiz_2d_int( pch_index_ji, nys, nyn, nxl, nxr, nbgp )
1087!
1088!--    Calculate global pch_index value (index of top of plant canopy from ground)
1089       pch_index = MAXVAL( pch_index_ji )
1090!
1091!--    Exchange pch_index from all processors
1092#if defined( __parallel )
1093       CALL MPI_ALLREDUCE( MPI_IN_PLACE, pch_index, 1, MPI_INTEGER,            &
1094                           MPI_MAX, comm2d, ierr)
1095#endif
1096
1097!--    Allocation of arrays pc_heating_rate, pc_transpiration_rate and pc_latent_rate
1098       ALLOCATE( pc_heating_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
1099       IF ( humidity )  THEN
1100          ALLOCATE( pc_transpiration_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
1101          pc_transpiration_rate = 0.0_wp
1102          ALLOCATE( pc_latent_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
1103          pc_latent_rate = 0.0_wp
1104       ENDIF
1105
1106!
1107!--    Initialization of the canopy heat source distribution due to heating
1108!--    of the canopy layers by incoming solar radiation, in case that a non-zero
1109!--    value is set for the canopy top heat flux (cthf), which equals the
1110!--    available net radiation at canopy top.
1111!--    The heat source distribution is calculated by a decaying exponential
1112!--    function of the downward cumulative leaf area index (cum_lai_hf),
1113!--    assuming that the foliage inside the plant canopy is heated by solar
1114!--    radiation penetrating the canopy layers according to the distribution
1115!--    of net radiation as suggested by Brown & Covey (1966; Agric. Meteorol. 3,
1116!--    73–96). This approach has been applied e.g. by Shaw & Schumann (1992;
1117!--    Bound.-Layer Meteorol. 61, 47–64).
1118!--    When using the radiation_interactions, canopy heating (pc_heating_rate)
1119!--    and plant canopy transpiration (pc_transpiration_rate, pc_latent_rate)
1120!--    are calculated in the RTM after the calculation of radiation.
1121!--    We cannot use variable radiation_interactions here to determine the situation
1122!--    as it is assigned in init_3d_model after the call of pcm_init.
1123       IF ( cthf /= 0.0_wp )  THEN
1124
1125          ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1126!
1127!--       Piecewise calculation of the cumulative leaf area index by vertical
1128!--       integration of the leaf area density
1129          cum_lai_hf(:,:,:) = 0.0_wp
1130          DO  i = nxlg, nxrg
1131             DO  j = nysg, nyng
1132                DO  k = pch_index_ji(j,i)-1, 0, -1
1133                   IF ( k == pch_index_ji(j,i)-1 )  THEN
1134                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
1135                         ( 0.5_wp * lad_s(k+1,j,i) *                           &
1136                           ( zw(k+1) - zu(k+1) ) )  +                          &
1137                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
1138                                                 lad_s(k,j,i) ) +              &
1139                                      lad_s(k+1,j,i) ) *                       &
1140                           ( zu(k+1) - zw(k) ) ) 
1141                   ELSE
1142                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
1143                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) +              &
1144                                                 lad_s(k+1,j,i) ) +            &
1145                                      lad_s(k+1,j,i) ) *                       &
1146                           ( zw(k+1) - zu(k+1) ) )  +                          &
1147                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
1148                                                 lad_s(k,j,i) ) +              &
1149                                      lad_s(k+1,j,i) ) *                       &
1150                           ( zu(k+1) - zw(k) ) )
1151                   ENDIF
1152                ENDDO
1153             ENDDO
1154          ENDDO
1155
1156!           
1157!--       In areas with canopy the surface value of the canopy heat
1158!--       flux distribution overrides the surface heat flux (shf)
1159!--       Start with default surface type
1160          DO  m = 1, surf_def_h(0)%ns
1161             k = surf_def_h(0)%k(m)
1162             IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                &
1163                surf_def_h(0)%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
1164          ENDDO
1165!
1166!--       Natural surfaces
1167          DO  m = 1, surf_lsm_h%ns
1168             k = surf_lsm_h%k(m)
1169             IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                &
1170                surf_lsm_h%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
1171          ENDDO
1172!
1173!--       Urban surfaces
1174          DO  m = 1, surf_usm_h%ns
1175             k = surf_usm_h%k(m)
1176             IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                &
1177                surf_usm_h%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
1178          ENDDO
1179!
1180!
1181!--       Calculation of the heating rate (K/s) within the different layers of
1182!--       the plant canopy. Calculation is only necessary in areas covered with
1183!--       canopy.
1184!--       Within the different canopy layers the plant-canopy heating
1185!--       rate (pc_heating_rate) is calculated as the vertical
1186!--       divergence of the canopy heat fluxes at the top and bottom
1187!--       of the respective layer
1188          DO  i = nxlg, nxrg
1189             DO  j = nysg, nyng
1190                DO  k = 1, pch_index_ji(j,i)
1191                   IF ( cum_lai_hf(0,j,i) /= 0.0_wp )  THEN
1192                      pc_heating_rate(k,j,i) = cthf *                          &
1193                                ( exp(-ext_coef*cum_lai_hf(k,j,i)) -           &
1194                                  exp(-ext_coef*cum_lai_hf(k-1,j,i) ) ) / dzw(k)
1195                   ENDIF
1196                ENDDO
1197             ENDDO
1198          ENDDO
1199
1200       ENDIF
1201
1202
1203    END SUBROUTINE pcm_init
1204
1205
1206!------------------------------------------------------------------------------!
1207! Description:
1208! ------------
1209!> Parin for &plant_canopy_parameters for plant canopy model
1210!------------------------------------------------------------------------------!
1211    SUBROUTINE pcm_parin
1212
1213       USE control_parameters,                                                 &
1214           ONLY:  message_string, plant_canopy
1215
1216       IMPLICIT NONE
1217
1218       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
1219       
1220       NAMELIST /plant_canopy_parameters/                                      &
1221                                  alpha_lad, beta_lad, canopy_drag_coeff,      &
1222                                  canopy_mode, cthf,                           &
1223                                  lad_surface, lad_type_coef,                  & 
1224                                  lad_vertical_gradient,                       &
1225                                  lad_vertical_gradient_level,                 &
1226                                  lai_beta,                                    &
1227                                  leaf_scalar_exch_coeff,                      &
1228                                  leaf_surface_conc, pch_index,                &
1229                                  plant_canopy_transpiration
1230
1231       NAMELIST /canopy_par/      alpha_lad, beta_lad, canopy_drag_coeff,      &
1232                                  canopy_mode, cthf,                           &
1233                                  lad_surface, lad_type_coef,                  & 
1234                                  lad_vertical_gradient,                       &
1235                                  lad_vertical_gradient_level,                 &
1236                                  lai_beta,                                    &
1237                                  leaf_scalar_exch_coeff,                      &
1238                                  leaf_surface_conc, pch_index,                &
1239                                  plant_canopy_transpiration
1240
1241       line = ' '
1242
1243!
1244!--    Try to find radiation model package
1245       REWIND ( 11 )
1246       line = ' '
1247       DO WHILE ( INDEX( line, '&plant_canopy_parameters' ) == 0 )
1248          READ ( 11, '(A)', END=12 )  line
1249       ENDDO
1250       BACKSPACE ( 11 )
1251
1252!
1253!--    Read user-defined namelist
1254       READ ( 11, plant_canopy_parameters, ERR = 10 )
1255
1256!
1257!--    Set flag that indicates that the radiation model is switched on
1258       plant_canopy = .TRUE.
1259
1260       GOTO 14
1261
1262 10    BACKSPACE( 11 )
1263       READ( 11 , '(A)') line
1264       CALL parin_fail_message( 'plant_canopy_parameters', line )
1265!
1266!--    Try to find old namelist
1267 12    REWIND ( 11 )
1268       line = ' '
1269       DO WHILE ( INDEX( line, '&canopy_par' ) == 0 )
1270          READ ( 11, '(A)', END=14 )  line
1271       ENDDO
1272       BACKSPACE ( 11 )
1273
1274!
1275!--    Read user-defined namelist
1276       READ ( 11, canopy_par, ERR = 13, END = 14 )
1277
1278       message_string = 'namelist canopy_par is deprecated and will be ' // &
1279                     'removed in near future. Please use namelist ' //      &
1280                     'plant_canopy_parameters instead'
1281       CALL message( 'pcm_parin', 'PA0487', 0, 1, 0, 6, 0 )
1282
1283!
1284!--    Set flag that indicates that the radiation model is switched on
1285       plant_canopy = .TRUE.
1286
1287       GOTO 14
1288
1289 13    BACKSPACE( 11 )
1290       READ( 11 , '(A)') line
1291       CALL parin_fail_message( 'canopy_par', line )
1292
1293 14    CONTINUE
1294
1295
1296    END SUBROUTINE pcm_parin
1297
1298
1299
1300!------------------------------------------------------------------------------!
1301! Description:
1302! ------------
1303!
1304!> Loads 3D plant canopy data from file. File format is as follows:
1305!>
1306!> num_levels
1307!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
1308!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
1309!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
1310!> ...
1311!>
1312!> i.e. first line determines number of levels and further lines represent plant
1313!> canopy data, one line per column and variable. In each data line,
1314!> dtype represents variable to be set:
1315!>
1316!> dtype=1: leaf area density (lad_s)
1317!> dtype=2....n: some additional plant canopy input data quantity
1318!>
1319!> Zeros are added automatically above num_levels until top of domain.  Any
1320!> non-specified (x,y) columns have zero values as default.
1321!------------------------------------------------------------------------------!
1322    SUBROUTINE pcm_read_plant_canopy_3d
1323   
1324       USE control_parameters,                                                 &
1325           ONLY:  coupling_char, message_string
1326
1327       USE indices,                                                            &
1328           ONLY:  nbgp
1329           
1330       IMPLICIT NONE
1331
1332       INTEGER(iwp)                        ::  dtype     !< type of input data (1=lad)
1333       INTEGER(iwp)                        ::  pctype    !< type of plant canopy (deciduous,non-deciduous,...)
1334       INTEGER(iwp)                        ::  i, j      !< running index
1335       INTEGER(iwp)                        ::  nzp       !< number of vertical layers of plant canopy
1336       INTEGER(iwp)                        ::  nzpltop   !<
1337       INTEGER(iwp)                        ::  nzpl      !<
1338       INTEGER(iwp)                        ::  kk        !<
1339       
1340       REAL(wp), DIMENSION(:), ALLOCATABLE ::  col   !< vertical column of input data
1341
1342!
1343!--    Initialize lad_s array
1344       lad_s = 0.0_wp
1345       
1346!
1347!--    Open and read plant canopy input data
1348       OPEN(152, FILE='PLANT_CANOPY_DATA_3D' // TRIM( coupling_char ),         &
1349                 ACCESS='SEQUENTIAL', ACTION='READ', STATUS='OLD',             &
1350                 FORM='FORMATTED', ERR=515)
1351       READ(152, *, ERR=516, END=517)  nzp   !< read first line = number of vertical layers
1352       nzpltop = MIN(nzt+1, nzb+nzp-1)
1353       nzpl = nzpltop - nzb + 1    !< no. of layers to assign
1354       ALLOCATE( col(0:nzp-1) )
1355
1356       DO
1357          READ(152, *, ERR=516, END=517) dtype, i, j, pctype, col(:)
1358          IF ( i < nxlg  .OR.  i > nxrg  .OR.  j < nysg  .OR.  j > nyng )  CYCLE
1359         
1360          SELECT CASE (dtype)
1361             CASE( 1 )   !< leaf area density
1362!
1363!--             This is just the pure canopy layer assumed to be grounded to
1364!--             a flat domain surface. At locations where plant canopy sits
1365!--             on top of any kind of topography, the vertical plant column
1366!--             must be "lifted", which is done in SUBROUTINE pcm_tendency.           
1367                IF ( pctype < 0  .OR.  pctype > 10 )  THEN   !< incorrect plant canopy type
1368                   WRITE( message_string, * ) 'Incorrect type of plant canopy. '   //  &
1369                                              'Allowed values 0 <= pctype <= 10, ' //  &
1370                                              'but pctype is ', pctype
1371                   CALL message( 'pcm_read_plant_canopy_3d', 'PA0349', 1, 2, 0, 6, 0 )
1372                ENDIF
1373                kk = get_topography_top_index_ji( j, i, 's' )
1374                lad_s(nzb:nzpltop-kk, j, i) = col(kk:nzpl-1)*lad_type_coef(pctype)
1375             CASE DEFAULT
1376                WRITE(message_string, '(a,i2,a)')   &
1377                     'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"'
1378                CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 )
1379          END SELECT
1380       ENDDO
1381
1382515    message_string = 'error opening file PLANT_CANOPY_DATA_3D'
1383       CALL message( 'pcm_read_plant_canopy_3d', 'PA0531', 1, 2, 0, 6, 0 )
1384
1385516    message_string = 'error reading file PLANT_CANOPY_DATA_3D'
1386       CALL message( 'pcm_read_plant_canopy_3d', 'PA0532', 1, 2, 0, 6, 0 )
1387
1388517    CLOSE(152)
1389       DEALLOCATE( col )
1390       
1391       CALL exchange_horiz( lad_s, nbgp )
1392       
1393    END SUBROUTINE pcm_read_plant_canopy_3d
1394   
1395   
1396
1397!------------------------------------------------------------------------------!
1398! Description:
1399! ------------
1400!> Calculation of the tendency terms, accounting for the effect of the plant
1401!> canopy on momentum and scalar quantities.
1402!>
1403!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
1404!> (defined on scalar grid), as initialized in subroutine pcm_init.
1405!> The lad on the w-grid is vertically interpolated from the surrounding
1406!> lad_s. The upper boundary of the canopy is defined on the w-grid at
1407!> k = pch_index. Here, the lad is zero.
1408!>
1409!> The canopy drag must be limited (previously accounted for by calculation of
1410!> a limiting canopy timestep for the determination of the maximum LES timestep
1411!> in subroutine timestep), since it is physically impossible that the canopy
1412!> drag alone can locally change the sign of a velocity component. This
1413!> limitation is realized by calculating preliminary tendencies and velocities.
1414!> It is subsequently checked if the preliminary new velocity has a different
1415!> sign than the current velocity. If so, the tendency is limited in a way that
1416!> the velocity can at maximum be reduced to zero by the canopy drag.
1417!>
1418!>
1419!> Call for all grid points
1420!------------------------------------------------------------------------------!
1421    SUBROUTINE pcm_tendency( component )
1422
1423
1424       USE control_parameters,                                                 &
1425           ONLY:  dt_3d, message_string
1426
1427       USE kinds
1428
1429       IMPLICIT NONE
1430
1431       INTEGER(iwp) ::  component !< prognostic variable (u,v,w,pt,q,e)
1432       INTEGER(iwp) ::  i         !< running index
1433       INTEGER(iwp) ::  j         !< running index
1434       INTEGER(iwp) ::  k         !< running index
1435       INTEGER(iwp) ::  k_wall    !< vertical index of topography top
1436       INTEGER(iwp) ::  kk        !< running index for flat lad arrays
1437
1438       REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
1439       REAL(wp) ::  lad_local !< local lad value
1440       REAL(wp) ::  pre_tend  !< preliminary tendency
1441       REAL(wp) ::  pre_u     !< preliminary u-value
1442       REAL(wp) ::  pre_v     !< preliminary v-value
1443       REAL(wp) ::  pre_w     !< preliminary w-value
1444
1445
1446       ddt_3d = 1.0_wp / dt_3d
1447 
1448!
1449!--    Compute drag for the three velocity components and the SGS-TKE:
1450       SELECT CASE ( component )
1451
1452!
1453!--       u-component
1454          CASE ( 1 )
1455             DO  i = nxlu, nxr
1456                DO  j = nys, nyn
1457!
1458!--                Determine topography-top index on u-grid
1459                   k_wall = get_topography_top_index_ji( j, i, 'u' )
1460                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
1461
1462                      kk = k - k_wall   !- lad arrays are defined flat
1463!
1464!--                   In order to create sharp boundaries of the plant canopy,
1465!--                   the lad on the u-grid at index (k,j,i) is equal to
1466!--                   lad_s(k,j,i), rather than being interpolated from the
1467!--                   surrounding lad_s, because this would yield smaller lad
1468!--                   at the canopy boundaries than inside of the canopy.
1469!--                   For the same reason, the lad at the rightmost(i+1)canopy
1470!--                   boundary on the u-grid equals lad_s(k,j,i).
1471                      lad_local = lad_s(kk,j,i)
1472                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp )&
1473                      THEN
1474                         lad_local = lad_s(kk,j,i-1)
1475                      ENDIF
1476
1477                      pre_tend = 0.0_wp
1478                      pre_u = 0.0_wp
1479!
1480!--                   Calculate preliminary value (pre_tend) of the tendency
1481                      pre_tend = - cdc *                                       &
1482                                   lad_local *                                 &
1483                                   SQRT( u(k,j,i)**2 +                         &
1484                                         ( 0.25_wp * ( v(k,j,i-1) +            &
1485                                                       v(k,j,i)   +            &
1486                                                       v(k,j+1,i) +            &
1487                                                       v(k,j+1,i-1) )          &
1488                                         )**2 +                                &
1489                                         ( 0.25_wp * ( w(k-1,j,i-1) +          &
1490                                                       w(k-1,j,i)   +          &
1491                                                       w(k,j,i-1)   +          &
1492                                                       w(k,j,i) )              &
1493                                         )**2                                  &
1494                                       ) *                                     &
1495                                   u(k,j,i)
1496
1497!
1498!--                   Calculate preliminary new velocity, based on pre_tend
1499                      pre_u = u(k,j,i) + dt_3d * pre_tend
1500!
1501!--                   Compare sign of old velocity and new preliminary velocity,
1502!--                   and in case the signs are different, limit the tendency
1503                      IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
1504                         pre_tend = - u(k,j,i) * ddt_3d
1505                      ELSE
1506                         pre_tend = pre_tend
1507                      ENDIF
1508!
1509!--                   Calculate final tendency
1510                      tend(k,j,i) = tend(k,j,i) + pre_tend
1511
1512                   ENDDO
1513                ENDDO
1514             ENDDO
1515
1516!
1517!--       v-component
1518          CASE ( 2 )
1519             DO  i = nxl, nxr
1520                DO  j = nysv, nyn
1521!
1522!--                Determine topography-top index on v-grid
1523                   k_wall = get_topography_top_index_ji( j, i, 'v' )
1524
1525                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
1526
1527                      kk = k - k_wall   !- lad arrays are defined flat
1528!
1529!--                   In order to create sharp boundaries of the plant canopy,
1530!--                   the lad on the v-grid at index (k,j,i) is equal to
1531!--                   lad_s(k,j,i), rather than being interpolated from the
1532!--                   surrounding lad_s, because this would yield smaller lad
1533!--                   at the canopy boundaries than inside of the canopy.
1534!--                   For the same reason, the lad at the northmost(j+1) canopy
1535!--                   boundary on the v-grid equals lad_s(k,j,i).
1536                      lad_local = lad_s(kk,j,i)
1537                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp )&
1538                      THEN
1539                         lad_local = lad_s(kk,j-1,i)
1540                      ENDIF
1541
1542                      pre_tend = 0.0_wp
1543                      pre_v = 0.0_wp
1544!
1545!--                   Calculate preliminary value (pre_tend) of the tendency
1546                      pre_tend = - cdc *                                       &
1547                                   lad_local *                                 &
1548                                   SQRT( ( 0.25_wp * ( u(k,j-1,i)   +          &
1549                                                       u(k,j-1,i+1) +          &
1550                                                       u(k,j,i)     +          &
1551                                                       u(k,j,i+1) )            &
1552                                         )**2 +                                &
1553                                         v(k,j,i)**2 +                         &
1554                                         ( 0.25_wp * ( w(k-1,j-1,i) +          &
1555                                                       w(k-1,j,i)   +          &
1556                                                       w(k,j-1,i)   +          &
1557                                                       w(k,j,i) )              &
1558                                         )**2                                  &
1559                                       ) *                                     &
1560                                   v(k,j,i)
1561
1562!
1563!--                   Calculate preliminary new velocity, based on pre_tend
1564                      pre_v = v(k,j,i) + dt_3d * pre_tend
1565!
1566!--                   Compare sign of old velocity and new preliminary velocity,
1567!--                   and in case the signs are different, limit the tendency
1568                      IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
1569                         pre_tend = - v(k,j,i) * ddt_3d
1570                      ELSE
1571                         pre_tend = pre_tend
1572                      ENDIF
1573!
1574!--                   Calculate final tendency
1575                      tend(k,j,i) = tend(k,j,i) + pre_tend
1576
1577                   ENDDO
1578                ENDDO
1579             ENDDO
1580
1581!
1582!--       w-component
1583          CASE ( 3 )
1584             DO  i = nxl, nxr
1585                DO  j = nys, nyn
1586!
1587!--                Determine topography-top index on w-grid
1588                   k_wall = get_topography_top_index_ji( j, i, 'w' )
1589
1590                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i) - 1
1591
1592                      kk = k - k_wall   !- lad arrays are defined flat
1593
1594                      pre_tend = 0.0_wp
1595                      pre_w = 0.0_wp
1596!
1597!--                   Calculate preliminary value (pre_tend) of the tendency
1598                      pre_tend = - cdc *                                       &
1599                                   (0.5_wp *                                   &
1600                                      ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) *   &
1601                                   SQRT( ( 0.25_wp * ( u(k,j,i)   +            &
1602                                                       u(k,j,i+1) +            &
1603                                                       u(k+1,j,i) +            &
1604                                                       u(k+1,j,i+1) )          &
1605                                         )**2 +                                &
1606                                         ( 0.25_wp * ( v(k,j,i)   +            &
1607                                                       v(k,j+1,i) +            &
1608                                                       v(k+1,j,i) +            &
1609                                                       v(k+1,j+1,i) )          &
1610                                         )**2 +                                &
1611                                         w(k,j,i)**2                           &
1612                                       ) *                                     &
1613                                   w(k,j,i)
1614!
1615!--                   Calculate preliminary new velocity, based on pre_tend
1616                      pre_w = w(k,j,i) + dt_3d * pre_tend
1617!
1618!--                   Compare sign of old velocity and new preliminary velocity,
1619!--                   and in case the signs are different, limit the tendency
1620                      IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
1621                         pre_tend = - w(k,j,i) * ddt_3d
1622                      ELSE
1623                         pre_tend = pre_tend
1624                      ENDIF
1625!
1626!--                   Calculate final tendency
1627                      tend(k,j,i) = tend(k,j,i) + pre_tend
1628
1629                   ENDDO
1630                ENDDO
1631             ENDDO
1632
1633!
1634!--       potential temperature
1635          CASE ( 4 )
1636             IF ( humidity ) THEN
1637                DO  i = nxl, nxr
1638                   DO  j = nys, nyn
1639!--                   Determine topography-top index on scalar-grid
1640                      k_wall = get_topography_top_index_ji( j, i, 's' )
1641                      DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
1642                         kk = k - k_wall   !- lad arrays are defined flat
1643                         tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) - pc_latent_rate(kk,j,i)
1644                      ENDDO
1645                   ENDDO
1646                ENDDO
1647             ELSE
1648                DO  i = nxl, nxr
1649                   DO  j = nys, nyn
1650!--                   Determine topography-top index on scalar-grid
1651                      k_wall = get_topography_top_index_ji( j, i, 's' )
1652                      DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
1653                         kk = k - k_wall   !- lad arrays are defined flat
1654                         tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
1655                      ENDDO
1656                   ENDDO
1657                ENDDO
1658             ENDIF
1659
1660!
1661!--       humidity
1662          CASE ( 5 )
1663             DO  i = nxl, nxr
1664                DO  j = nys, nyn
1665!
1666!--                Determine topography-top index on scalar-grid
1667                   k_wall = get_topography_top_index_ji( j, i, 's' )
1668
1669                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
1670
1671                      kk = k - k_wall   !- lad arrays are defined flat
1672
1673                      IF ( .NOT. plant_canopy_transpiration ) THEN
1674                         ! pc_transpiration_rate is calculated in radiation model
1675                         ! in case of plant_canopy_transpiration = .T.
1676                         ! to include also the dependecy to the radiation
1677                         ! in the plant canopy box
1678                         pc_transpiration_rate(kk,j,i) =  - lsec                  &
1679                                          * lad_s(kk,j,i) *                       &
1680                                          SQRT( ( 0.5_wp * ( u(k,j,i) +           &
1681                                                             u(k,j,i+1) )         &
1682                                                )**2 +                            &
1683                                                ( 0.5_wp * ( v(k,j,i) +           &
1684                                                             v(k,j+1,i) )         &
1685                                                )**2 +                            &
1686                                                ( 0.5_wp * ( w(k-1,j,i) +         &
1687                                                             w(k,j,i) )           &
1688                                                )**2                              &
1689                                              ) *                                 &
1690                                          ( q(k,j,i) - lsc )
1691                      ENDIF
1692
1693                      tend(k,j,i) = tend(k,j,i) + pc_transpiration_rate(kk,j,i)
1694                   ENDDO
1695                ENDDO
1696             ENDDO
1697
1698!
1699!--       sgs-tke
1700          CASE ( 6 )
1701             DO  i = nxl, nxr
1702                DO  j = nys, nyn
1703!
1704!--                Determine topography-top index on scalar-grid
1705                   k_wall = get_topography_top_index_ji( j, i, 's' )
1706
1707                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
1708
1709                      kk = k - k_wall   !- lad arrays are defined flat
1710                      tend(k,j,i) = tend(k,j,i) -                              &
1711                                       2.0_wp * cdc *                          &
1712                                       lad_s(kk,j,i) *                         &
1713                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
1714                                                          u(k,j,i+1) )         &
1715                                             )**2 +                            &
1716                                             ( 0.5_wp * ( v(k,j,i) +           &
1717                                                          v(k,j+1,i) )         &
1718                                             )**2 +                            &
1719                                             ( 0.5_wp * ( w(k,j,i) +           &
1720                                                          w(k+1,j,i) )         &
1721                                             )**2                              &
1722                                           ) *                                 &
1723                                       e(k,j,i)
1724                   ENDDO
1725                ENDDO
1726             ENDDO 
1727!
1728!--       scalar concentration
1729          CASE ( 7 )
1730             DO  i = nxl, nxr
1731                DO  j = nys, nyn
1732!
1733!--                Determine topography-top index on scalar-grid
1734                   k_wall = get_topography_top_index_ji( j, i, 's' )
1735
1736                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
1737
1738                      kk = k - k_wall   !- lad arrays are defined flat
1739                      tend(k,j,i) = tend(k,j,i) -                              &
1740                                       lsec *                                  &
1741                                       lad_s(kk,j,i) *                         &
1742                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
1743                                                          u(k,j,i+1) )         &
1744                                             )**2 +                            &
1745                                             ( 0.5_wp * ( v(k,j,i) +           &
1746                                                          v(k,j+1,i) )         &
1747                                             )**2 +                            &
1748                                             ( 0.5_wp * ( w(k-1,j,i) +         & 
1749                                                          w(k,j,i) )           &
1750                                             )**2                              &
1751                                           ) *                                 &
1752                                       ( s(k,j,i) - lsc )
1753                   ENDDO
1754                ENDDO
1755             ENDDO   
1756
1757
1758
1759          CASE DEFAULT
1760
1761             WRITE( message_string, * ) 'wrong component: ', component
1762             CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 ) 
1763
1764       END SELECT
1765
1766    END SUBROUTINE pcm_tendency
1767
1768
1769!------------------------------------------------------------------------------!
1770! Description:
1771! ------------
1772!> Calculation of the tendency terms, accounting for the effect of the plant
1773!> canopy on momentum and scalar quantities.
1774!>
1775!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
1776!> (defined on scalar grid), as initialized in subroutine pcm_init.
1777!> The lad on the w-grid is vertically interpolated from the surrounding
1778!> lad_s. The upper boundary of the canopy is defined on the w-grid at
1779!> k = pch_index. Here, the lad is zero.
1780!>
1781!> The canopy drag must be limited (previously accounted for by calculation of
1782!> a limiting canopy timestep for the determination of the maximum LES timestep
1783!> in subroutine timestep), since it is physically impossible that the canopy
1784!> drag alone can locally change the sign of a velocity component. This
1785!> limitation is realized by calculating preliminary tendencies and velocities.
1786!> It is subsequently checked if the preliminary new velocity has a different
1787!> sign than the current velocity. If so, the tendency is limited in a way that
1788!> the velocity can at maximum be reduced to zero by the canopy drag.
1789!>
1790!>
1791!> Call for grid point i,j
1792!------------------------------------------------------------------------------!
1793    SUBROUTINE pcm_tendency_ij( i, j, component )
1794
1795
1796       USE control_parameters,                                                 &
1797           ONLY:  dt_3d, message_string
1798
1799       USE kinds
1800
1801       IMPLICIT NONE
1802
1803       INTEGER(iwp) ::  component !< prognostic variable (u,v,w,pt,q,e)
1804       INTEGER(iwp) ::  i         !< running index
1805       INTEGER(iwp) ::  j         !< running index
1806       INTEGER(iwp) ::  k         !< running index
1807       INTEGER(iwp) ::  k_wall    !< vertical index of topography top
1808       INTEGER(iwp) ::  kk        !< running index for flat lad arrays
1809
1810       REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
1811       REAL(wp) ::  lad_local !< local lad value
1812       REAL(wp) ::  pre_tend  !< preliminary tendency
1813       REAL(wp) ::  pre_u     !< preliminary u-value
1814       REAL(wp) ::  pre_v     !< preliminary v-value
1815       REAL(wp) ::  pre_w     !< preliminary w-value
1816
1817
1818       ddt_3d = 1.0_wp / dt_3d
1819!
1820!--    Compute drag for the three velocity components and the SGS-TKE
1821       SELECT CASE ( component )
1822
1823!
1824!--       u-component
1825          CASE ( 1 )
1826!
1827!--          Determine topography-top index on u-grid
1828             k_wall = get_topography_top_index_ji( j, i, 'u' )
1829             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
1830
1831                kk = k - k_wall  !- lad arrays are defined flat
1832
1833!
1834!--             In order to create sharp boundaries of the plant canopy,
1835!--             the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
1836!--             rather than being interpolated from the surrounding lad_s,
1837!--             because this would yield smaller lad at the canopy boundaries
1838!--             than inside of the canopy.
1839!--             For the same reason, the lad at the rightmost(i+1)canopy
1840!--             boundary on the u-grid equals lad_s(k,j,i).
1841                lad_local = lad_s(kk,j,i)
1842                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp )  THEN
1843                   lad_local = lad_s(kk,j,i-1)
1844                ENDIF
1845
1846                pre_tend = 0.0_wp
1847                pre_u = 0.0_wp
1848!
1849!--             Calculate preliminary value (pre_tend) of the tendency
1850                pre_tend = - cdc *                                             &
1851                             lad_local *                                       &   
1852                             SQRT( u(k,j,i)**2 +                               &
1853                                   ( 0.25_wp * ( v(k,j,i-1)  +                 &
1854                                                 v(k,j,i)    +                 &
1855                                                 v(k,j+1,i)  +                 &
1856                                                 v(k,j+1,i-1) )                &
1857                                   )**2 +                                      &
1858                                   ( 0.25_wp * ( w(k-1,j,i-1) +                &
1859                                                 w(k-1,j,i)   +                &
1860                                                 w(k,j,i-1)   +                &
1861                                                 w(k,j,i) )                    &
1862                                   )**2                                        &
1863                                 ) *                                           &
1864                             u(k,j,i)
1865
1866!
1867!--             Calculate preliminary new velocity, based on pre_tend
1868                pre_u = u(k,j,i) + dt_3d * pre_tend
1869!
1870!--             Compare sign of old velocity and new preliminary velocity,
1871!--             and in case the signs are different, limit the tendency
1872                IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
1873                   pre_tend = - u(k,j,i) * ddt_3d
1874                ELSE
1875                   pre_tend = pre_tend
1876                ENDIF
1877!
1878!--             Calculate final tendency
1879                tend(k,j,i) = tend(k,j,i) + pre_tend
1880             ENDDO
1881
1882
1883!
1884!--       v-component
1885          CASE ( 2 )
1886!
1887!--          Determine topography-top index on v-grid
1888             k_wall = get_topography_top_index_ji( j, i, 'v' )
1889
1890             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
1891
1892                kk = k - k_wall  !- lad arrays are defined flat
1893!
1894!--             In order to create sharp boundaries of the plant canopy,
1895!--             the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
1896!--             rather than being interpolated from the surrounding lad_s,
1897!--             because this would yield smaller lad at the canopy boundaries
1898!--             than inside of the canopy.
1899!--             For the same reason, the lad at the northmost(j+1)canopy
1900!--             boundary on the v-grid equals lad_s(k,j,i).
1901                lad_local = lad_s(kk,j,i)
1902                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp )  THEN
1903                   lad_local = lad_s(kk,j-1,i)
1904                ENDIF
1905
1906                pre_tend = 0.0_wp
1907                pre_v = 0.0_wp
1908!
1909!--             Calculate preliminary value (pre_tend) of the tendency
1910                pre_tend = - cdc *                                             &
1911                             lad_local *                                       &
1912                             SQRT( ( 0.25_wp * ( u(k,j-1,i)   +                &
1913                                                 u(k,j-1,i+1) +                &
1914                                                 u(k,j,i)     +                &
1915                                                 u(k,j,i+1) )                  &
1916                                   )**2 +                                      &
1917                                   v(k,j,i)**2 +                               &
1918                                   ( 0.25_wp * ( w(k-1,j-1,i) +                &
1919                                                 w(k-1,j,i)   +                &
1920                                                 w(k,j-1,i)   +                &
1921                                                 w(k,j,i) )                    &
1922                                   )**2                                        &
1923                                 ) *                                           &
1924                             v(k,j,i)
1925
1926!
1927!--             Calculate preliminary new velocity, based on pre_tend
1928                pre_v = v(k,j,i) + dt_3d * pre_tend
1929!
1930!--             Compare sign of old velocity and new preliminary velocity,
1931!--             and in case the signs are different, limit the tendency
1932                IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
1933                   pre_tend = - v(k,j,i) * ddt_3d
1934                ELSE
1935                   pre_tend = pre_tend
1936                ENDIF
1937!
1938!--             Calculate final tendency
1939                tend(k,j,i) = tend(k,j,i) + pre_tend
1940             ENDDO
1941
1942
1943!
1944!--       w-component
1945          CASE ( 3 )
1946!
1947!--          Determine topography-top index on w-grid
1948             k_wall = get_topography_top_index_ji( j, i, 'w' )
1949
1950             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i) - 1
1951
1952                kk = k - k_wall  !- lad arrays are defined flat
1953
1954                pre_tend = 0.0_wp
1955                pre_w = 0.0_wp
1956!
1957!--             Calculate preliminary value (pre_tend) of the tendency
1958                pre_tend = - cdc *                                             &
1959                             (0.5_wp *                                         &
1960                                ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) *         &
1961                             SQRT( ( 0.25_wp * ( u(k,j,i)    +                 & 
1962                                                 u(k,j,i+1)  +                 &
1963                                                 u(k+1,j,i)  +                 &
1964                                                 u(k+1,j,i+1) )                &
1965                                   )**2 +                                      &
1966                                   ( 0.25_wp * ( v(k,j,i)    +                 &
1967                                                 v(k,j+1,i)  +                 &
1968                                                 v(k+1,j,i)  +                 &
1969                                                 v(k+1,j+1,i) )                &
1970                                   )**2 +                                      &
1971                                   w(k,j,i)**2                                 &
1972                                 ) *                                           &
1973                             w(k,j,i)
1974!
1975!--             Calculate preliminary new velocity, based on pre_tend
1976                pre_w = w(k,j,i) + dt_3d * pre_tend
1977!
1978!--             Compare sign of old velocity and new preliminary velocity,
1979!--             and in case the signs are different, limit the tendency
1980                IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
1981                   pre_tend = - w(k,j,i) * ddt_3d
1982                ELSE
1983                   pre_tend = pre_tend
1984                ENDIF
1985!
1986!--             Calculate final tendency
1987                tend(k,j,i) = tend(k,j,i) + pre_tend
1988             ENDDO
1989
1990!
1991!--       potential temperature
1992          CASE ( 4 )
1993!
1994!--          Determine topography-top index on scalar grid
1995             k_wall = get_topography_top_index_ji( j, i, 's' )
1996
1997             IF ( humidity ) THEN
1998                DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
1999                   kk = k - k_wall  !- lad arrays are defined flat
2000                   tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) - pc_latent_rate(kk,j,i)
2001                ENDDO
2002             ELSE
2003                DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
2004                   kk = k - k_wall  !- lad arrays are defined flat
2005                   tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
2006                ENDDO
2007             ENDIF
2008
2009!
2010!--       humidity
2011          CASE ( 5 )
2012!
2013!--          Determine topography-top index on scalar grid
2014             k_wall = get_topography_top_index_ji( j, i, 's' )
2015
2016             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
2017                kk = k - k_wall  !- lad arrays are defined flat
2018                IF ( .NOT. plant_canopy_transpiration ) THEN
2019                   ! pc_transpiration_rate is calculated in radiation model
2020                   ! in case of plant_canopy_transpiration = .T.
2021                   ! to include also the dependecy to the radiation
2022                   ! in the plant canopy box
2023                   pc_transpiration_rate(kk,j,i) = - lsec                         &
2024                                    * lad_s(kk,j,i) *                             &
2025                                    SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
2026                                                       u(k,j,i+1) )               &
2027                                          )**2  +                                 &
2028                                          ( 0.5_wp * ( v(k,j,i) +                 &
2029                                                       v(k,j+1,i) )               &
2030                                          )**2 +                                  &
2031                                          ( 0.5_wp * ( w(k-1,j,i) +               &
2032                                                       w(k,j,i) )                 &
2033                                          )**2                                    &
2034                                        ) *                                       &
2035                                    ( q(k,j,i) - lsc )
2036                ENDIF
2037
2038                tend(k,j,i) = tend(k,j,i) + pc_transpiration_rate(kk,j,i)
2039
2040             ENDDO   
2041
2042!
2043!--       sgs-tke
2044          CASE ( 6 )
2045!
2046!--          Determine topography-top index on scalar grid
2047             k_wall = get_topography_top_index_ji( j, i, 's' )
2048
2049             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
2050
2051                kk = k - k_wall
2052                tend(k,j,i) = tend(k,j,i) -                                    &
2053                                 2.0_wp * cdc *                                &
2054                                 lad_s(kk,j,i) *                               &
2055                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
2056                                                    u(k,j,i+1) )               &
2057                                       )**2 +                                  & 
2058                                       ( 0.5_wp * ( v(k,j,i) +                 &
2059                                                    v(k,j+1,i) )               &
2060                                       )**2 +                                  &
2061                                       ( 0.5_wp * ( w(k,j,i) +                 &
2062                                                    w(k+1,j,i) )               &
2063                                       )**2                                    &
2064                                     ) *                                       &
2065                                 e(k,j,i)
2066             ENDDO
2067             
2068!
2069!--       scalar concentration
2070          CASE ( 7 )
2071!
2072!--          Determine topography-top index on scalar grid
2073             k_wall = get_topography_top_index_ji( j, i, 's' )
2074
2075             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
2076
2077                kk = k - k_wall
2078                tend(k,j,i) = tend(k,j,i) -                                    &
2079                                 lsec *                                        &
2080                                 lad_s(kk,j,i) *                               &
2081                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
2082                                                    u(k,j,i+1) )               &
2083                                       )**2  +                                 &
2084                                       ( 0.5_wp * ( v(k,j,i) +                 &
2085                                                    v(k,j+1,i) )               &
2086                                       )**2 +                                  &
2087                                       ( 0.5_wp * ( w(k-1,j,i) +               &
2088                                                    w(k,j,i) )                 &
2089                                       )**2                                    &
2090                                     ) *                                       &
2091                                 ( s(k,j,i) - lsc )
2092             ENDDO               
2093
2094       CASE DEFAULT
2095
2096          WRITE( message_string, * ) 'wrong component: ', component
2097          CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 ) 
2098
2099       END SELECT
2100
2101    END SUBROUTINE pcm_tendency_ij
2102
2103
2104
2105 END MODULE plant_canopy_model_mod
Note: See TracBrowser for help on using the repository browser.