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

Last change on this file since 3753 was 3745, checked in by suehring, 6 years ago

document last commit

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