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

Last change on this file since 3655 was 3655, checked in by knoop, 3 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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