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

Last change on this file since 3885 was 3885, checked in by kanani, 2 years ago

restructure/add location/debug messages

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