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

Last change on this file since 3524 was 3524, checked in by raasch, 3 years ago

unused variables removed, missing working precision added, missing preprocessor directives added, bugfix concerning allocation of t_surf_wall_v in nopointer case, declaration statements rearranged to avoid compile time errors, mpi_abort arguments replaced to avoid compile errors

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