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

Last change on this file since 4793 was 4770, checked in by suehring, 4 years ago

Consider basal-area density as an additional sink for momentum in the prognostic equations for momentum and SGS TKE

  • Property svn:keywords set to Id
File size: 123.6 KB
RevLine 
[1826]1!> @file plant_canopy_model_mod.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[2000]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.
[1036]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!
[4360]17! Copyright 1997-2020 Leibniz Universitaet Hannover
[4392]18! Copyright 2017-2020 Institute of Computer Science of the
[3885]19!                     Czech Academy of Sciences, Prague
[2000]20!------------------------------------------------------------------------------!
[1036]21!
[257]22! Current revisions:
[2977]23! ------------------
[4671]24!
25!
[2214]26! Former revisions:
27! -----------------
[4495]28! $Id: plant_canopy_model_mod.f90 4770 2020-11-03 14:04:53Z maronga $
[4770]29! Consider basal-area density as an additional sink for momentum in the
30! prognostic equations for momentum and SGS TKE
31!
32! 4768 2020-11-02 19:11:23Z suehring
[4768]33! Enable 3D data output also with 64-bit precision
34!
35! 4671 2020-09-09 20:27:58Z pavelkrc
[4671]36! Implementation of downward facing USM and LSM surfaces
37!
38! 4535 2020-05-15 12:07:23Z raasch
[4535]39! bugfix for restart data format query
40!
41! 4525 2020-05-10 17:05:07Z raasch
[4525]42! bugfix for reading/writing pcm_...rate_av with MPI-IO
43!
44! 4517 2020-05-03 14:29:30Z raasch
[4517]45! added restart with MPI-IO for reading local arrays
46!
47! 4515 2020-04-30 16:37:18Z suehring
[4515]48! Rename error number again since this was recently given in -r 4511
49!
50! 4514 2020-04-30 16:29:59Z suehring
[4514]51! - Bugfix in output of pcm_heatrate_av in a restart run. In order to fix this,
52!   pch_index is now output for a restart run. Therefore, define global restart
53!   routines.
54! - Error message number renamed and check for PA0505 revised in order to also
55!   consider natural surfaces with plant-canopy.
56!
57! 4495 2020-04-13 20:11:20Z raasch
[4495]58! restart data handling with MPI-IO added
[4448]59!
[4495]60! 4457 2020-03-11 14:20:43Z raasch
61!
[4457]62! use statement for exchange horiz added
[4448]63! (salim) removed the error message PA0672 to consider PC 3d data via ascii file
64!
65! 4392 2020-01-31 16:14:57Z pavelkrc (resler)
66! Make pcm_heatrate_av, pcm_latentrate_av public to allow calculation
[4392]67! of averaged Bowen ratio in the user procedure
68!
69! 4381 2020-01-20 13:51:46Z suehring
[4381]70! Give error message 313 only once
71!
72! 4363 2020-01-07 18:11:28Z suehring
[4363]73! Fix for last commit
74!
75! 4362 2020-01-07 17:15:02Z suehring
[4362]76! Input of plant canopy variables from static driver moved to plant-canopy
77! model
78!
79! 4361 2020-01-07 12:22:38Z suehring
[4361]80! - Remove unused arrays in pmc_rrd_local
81! - Remove one exchange of ghost points
82!
83! 4360 2020-01-07 11:25:50Z suehring
[4360]84! - Bugfix, read restart data for time-averaged pcm output quantities
85! - Output of plant-canopy quantities will fill values
86!
87! 4356 2019-12-20 17:09:33Z suehring
[4356]88! Correct single message call, local check must be given by the respective
89! mpi rank.
90!
91! 4346 2019-12-18 11:55:56Z motisi
[4346]92! Introduction of wall_flags_total_0, which currently sets bits based on static
93! topography information used in wall_flags_static_0
94!
95! 4342 2019-12-16 13:49:14Z Giersch
[4342]96! Use statements moved to module level, ocean dependency removed, redundant
97! variables removed
98!
99! 4341 2019-12-16 10:43:49Z motisi
[4341]100! - Unification of variable names: pc_-variables now pcm_-variables
101!   (pc_latent_rate, pc_heating_rate, pc_transpiration_rate)
102! - Removal of pcm_bowenratio output
103! - Renamed canopy-mode 'block' to 'homogeneous'
104! - Renamed value 'read_from_file_3d' to 'read_from_file'
105! - Removal of confusing comment lines
106! - Replacement of k_wall by topo_top_ind
107! - Removal of Else-Statement in tendency-calculation
108!
109! 4335 2019-12-12 16:39:05Z suehring
[4335]110! Fix for LAD at building edges also implemented in vector branch.
111!
112! 4331 2019-12-10 18:25:02Z suehring
[4331]113! Typo corrected
114!
115! 4329 2019-12-10 15:46:36Z motisi
[4329]116! Renamed wall_flags_0 to wall_flags_static_0
117!
118! 4314 2019-11-29 10:29:20Z suehring
[4314]119! - Bugfix, plant canopy was still considered at building edges on for the u-
120!   and v-component.
121! - Relax restriction of LAD on building tops. LAD is only omitted at
122!   locations where building grid points emerged artificially by the
123!   topography filtering.
124!
125! 4309 2019-11-26 18:49:59Z suehring
[4309]126! Typo
127!
128! 4302 2019-11-22 13:15:56Z suehring
[4302]129! Omit tall canopy mapped on top of buildings
130!
131! 4279 2019-10-29 08:48:17Z scharf
[4279]132! unused variables removed
133!
134! 4258 2019-10-07 13:29:08Z scharf
[4278]135! changed check for static driver and fixed bugs in initialization and header
136!
137! 4258 2019-10-07 13:29:08Z suehring
[4258]138! Check if any LAD is prescribed when plant-canopy model is applied.
139!
140! 4226 2019-09-10 17:03:24Z suehring
[4226]141! Bugfix, missing initialization of heating rate
142!
143! 4221 2019-09-09 08:50:35Z suehring
[4220]144! Further bugfix in 3d data output for plant canopy
145!
146! 4216 2019-09-04 09:09:03Z suehring
[4216]147! Bugfixes in 3d data output
148!
149! 4205 2019-08-30 13:25:00Z suehring
[4205]150! Missing working precision + bugfix in calculation of wind speed
151!
152! 4188 2019-08-26 14:15:47Z suehring
[4188]153! Minor adjustment in error number
154!
155! 4187 2019-08-26 12:43:15Z suehring
[4187]156! Give specific error numbers instead of PA0999
157!
158! 4182 2019-08-22 15:20:23Z scharf
[4182]159! Corrected "Former revisions" section
160!
161! 4168 2019-08-16 13:50:17Z suehring
[4168]162! Replace function get_topography_top_index by topo_top_ind
163!
164! 4127 2019-07-30 14:47:10Z suehring
[4127]165! Output of 3D plant canopy variables changed. It is now relative to the local
166! terrain rather than located at the acutal vertical level in the model. This
167! way, the vertical dimension of the output can be significantly reduced.
168! (merge from branch resler)
169!
170! 3885 2019-04-11 11:29:34Z kanani
[3885]171! Changes related to global restructuring of location messages and introduction
172! of additional debug messages
173!
174! 3864 2019-04-05 09:01:56Z monakurppa
[3761]175! unsed variables removed
176!
177! 3745 2019-02-15 18:57:56Z suehring
[3745]178! Bugfix in transpiration, floating invalid when temperature
179! becomes > 40 degrees
180!
181! 3744 2019-02-15 18:38:58Z suehring
[3685]182! Some interface calls moved to module_interface + cleanup
183!
184! 3655 2019-01-07 16:51:22Z knoop
[3614]185! unused variables removed
[3498]186!
[4182]187! 138 2007-11-28 10:03:58Z letzel
188! Initial revision
189!
[138]190! Description:
191! ------------
[1682]192!> 1) Initialization of the canopy model, e.g. construction of leaf area density
[1826]193!> profile (subroutine pcm_init).
[1682]194!> 2) Calculation of sinks and sources of momentum, heat and scalar concentration
[1826]195!> due to canopy elements (subroutine pcm_tendency).
[3744]196!
197! @todo - precalculate constant terms in pcm_calc_transpiration_rate
[4216]198! @todo - unify variable names (pcm_, pc_, ...)
[4514]199! @todo - get rid-off dependency on radiation model
[138]200!------------------------------------------------------------------------------!
[1682]201 MODULE plant_canopy_model_mod
[4360]202
[1484]203    USE arrays_3d,                                                             &
[3449]204        ONLY:  dzu, dzw, e, exner, hyp, pt, q, s, tend, u, v, w, zu, zw
[138]205
[3449]206    USE basic_constants_and_equations_mod,                                     &
207        ONLY:  c_p, degc_to_k, l_v, lv_d_cp, r_d, rd_d_rv
[4360]208
[4342]209    USE bulk_cloud_model_mod,                                                  &
210        ONLY: bulk_cloud_model, microphysics_seifert
[3449]211
[3885]212    USE control_parameters,                                                    &
[4360]213        ONLY: average_count_3d,                                                &
214              coupling_char,                                                   &
215              debug_output,                                                    &
216              dt_3d,                                                           &
217              dz,                                                              &
218              humidity,                                                        &
[4514]219              land_surface,                                                    &
[4360]220              length,                                                          &
221              message_string,                                                  &
222              ocean_mode,                                                      &
223              passive_scalar,                                                  &
224              plant_canopy,                                                    &
[4495]225              restart_data_format_output,                                      &
[4360]226              restart_string,                                                  &
227              urban_surface
228
[4342]229    USE grid_variables,                                                        &
230        ONLY:  dx, dy
[3449]231
[1484]232    USE indices,                                                               &
233        ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
[4346]234               nz, nzb, nzt, topo_top_ind, wall_flags_total_0
[1484]235
236    USE kinds
[4360]237
[4342]238    USE netcdf_data_input_mod,                                                 &
[4362]239        ONLY:  input_pids_static,                                              &
240               char_fill,                                                      &
241               check_existence,                                                &
242               close_input_file,                                               &
243               get_attribute,                                                  &
244               get_dimension_length,                                           &
245               get_variable,                                                   &
246               inquire_num_variables,                                          &
247               inquire_variable_names,                                         &
248               input_file_static,                                              &
249               num_var_pids,                                                   &
250               open_read_file,                                                 &
251               pids_id,                                                        &
252               real_3d,                                                        &
253               vars_pids
[1484]254
[3449]255    USE pegrid
[4360]256
[4495]257    USE restart_data_mpi_io_mod,                                                                   &
[4517]258        ONLY:  rd_mpi_io_check_array,                                                              &
259               rrd_mpi_io,                                                                         &
[4514]260               wrd_mpi_io
[4495]261
[4342]262    USE surface_mod,                                                           &
263        ONLY: surf_def_h, surf_lsm_h, surf_usm_h
[3449]264
[1484]265
266    IMPLICIT NONE
267
[4341]268    CHARACTER (LEN=30)   ::  canopy_mode = 'homogeneous'           !< canopy coverage
[3449]269    LOGICAL              ::  plant_canopy_transpiration = .FALSE.  !< flag to switch calculation of transpiration and corresponding latent heat
270                                                                   !< for resolved plant canopy inside radiation model
271                                                                   !< (calls subroutine pcm_calc_transpiration_rate from module plant_canopy_mod)
[1484]272
[3449]273    INTEGER(iwp) ::  pch_index = 0                                 !< plant canopy height/top index
274    INTEGER(iwp) ::  lad_vertical_gradient_level_ind(10) = -9999   !< lad-profile levels (index)
[1484]275
[3449]276    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  pch_index_ji     !< local plant canopy top
[2696]277
[3449]278    LOGICAL ::  calc_beta_lad_profile = .FALSE.   !< switch for calc. of lad from beta func.
[1484]279
[2696]280    REAL(wp) ::  alpha_lad = 9999999.9_wp         !< coefficient for lad calculation
281    REAL(wp) ::  beta_lad = 9999999.9_wp          !< coefficient for lad calculation
282    REAL(wp) ::  canopy_drag_coeff = 0.0_wp       !< canopy drag coefficient (parameter)
283    REAL(wp) ::  cthf = 0.0_wp                    !< canopy top heat flux
284    REAL(wp) ::  dt_plant_canopy = 0.0_wp         !< timestep account. for canopy drag
285    REAL(wp) ::  ext_coef = 0.6_wp                !< extinction coefficient
286    REAL(wp) ::  lad_surface = 0.0_wp             !< lad surface value
287    REAL(wp) ::  lai_beta = 0.0_wp                !< leaf area index (lai) for lad calc.
288    REAL(wp) ::  leaf_scalar_exch_coeff = 0.0_wp  !< canopy scalar exchange coeff.
289    REAL(wp) ::  leaf_surface_conc = 0.0_wp       !< leaf surface concentration
[1484]290
[2696]291    REAL(wp) ::  lad_vertical_gradient(10) = 0.0_wp              !< lad gradient
292    REAL(wp) ::  lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m)
[1484]293
[2977]294    REAL(wp) ::  lad_type_coef(0:10) = 1.0_wp   !< multiplicative coeficients for particular types
295                                                !< of plant canopy (e.g. deciduous tree during winter)
296
[1682]297    REAL(wp), DIMENSION(:), ALLOCATABLE ::  lad            !< leaf area density
298    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pre_lad        !< preliminary lad
[4770]299
300    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  bad_s                    !< basal-area density on scalar-grid
[4127]301    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  cum_lai_hf               !< cumulative lai for heatflux calc.
302    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  lad_s                    !< lad on scalar-grid
[4341]303    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_heating_rate         !< plant canopy heating rate
304    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_transpiration_rate   !< plant canopy transpiration rate
305    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_latent_rate          !< plant canopy latent heating rate
[1484]306
[4127]307    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_heatrate_av          !< array for averaging plant canopy sensible heating rate
308    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_latentrate_av        !< array for averaging plant canopy latent heating rate
309    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_transpirationrate_av !< array for averaging plant canopy transpiration rate
310
[4362]311    TYPE(real_3d) ::  basal_area_density_f    !< input variable for basal area density - resolved vegetation
312    TYPE(real_3d) ::  leaf_area_density_f     !< input variable for leaf area density - resolved vegetation
313    TYPE(real_3d) ::  root_area_density_lad_f !< input variable for root area density - resolved vegetation
314
[1484]315    SAVE
316
[138]317    PRIVATE
[4342]318   
[1826]319!
320!-- Public functions
[4360]321    PUBLIC pcm_calc_transpiration_rate,                                       &
322           pcm_check_data_output,                                             &
323           pcm_check_parameters,                                              &
324           pcm_3d_data_averaging,                                             &
325           pcm_data_output_3d,                                                &
326           pcm_define_netcdf_grid,                                            &
327           pcm_header,                                                        &
328           pcm_init,                                                          &
329           pcm_parin,                                                         &
[4514]330           pcm_rrd_global,                                                    &
[4360]331           pcm_rrd_local,                                                     &
332           pcm_tendency,                                                      &
[4514]333           pcm_wrd_global,                                                    &
[4360]334           pcm_wrd_local
[138]335
[1826]336!
337!-- Public variables and constants
[4342]338    PUBLIC canopy_drag_coeff, pcm_heating_rate, pcm_transpiration_rate,       &
339           pcm_latent_rate, canopy_mode, cthf, dt_plant_canopy, lad, lad_s,   &
[4392]340           pch_index, plant_canopy_transpiration,                             &
341           pcm_heatrate_av, pcm_latentrate_av
[1484]342
[3449]343    INTERFACE pcm_calc_transpiration_rate
344       MODULE PROCEDURE pcm_calc_transpiration_rate
345    END INTERFACE pcm_calc_transpiration_rate
346
[2209]347    INTERFACE pcm_check_data_output
348       MODULE PROCEDURE pcm_check_data_output
349    END INTERFACE pcm_check_data_output
350   
[1826]351    INTERFACE pcm_check_parameters
352       MODULE PROCEDURE pcm_check_parameters
[2209]353    END INTERFACE pcm_check_parameters
354
[4127]355    INTERFACE pcm_3d_data_averaging
356       MODULE PROCEDURE pcm_3d_data_averaging
357    END INTERFACE pcm_3d_data_averaging
358
[2209]359    INTERFACE pcm_data_output_3d
360       MODULE PROCEDURE pcm_data_output_3d
361    END INTERFACE pcm_data_output_3d
362
363    INTERFACE pcm_define_netcdf_grid
364       MODULE PROCEDURE pcm_define_netcdf_grid
365    END INTERFACE pcm_define_netcdf_grid
[1826]366   
367     INTERFACE pcm_header
368       MODULE PROCEDURE pcm_header
369    END INTERFACE pcm_header       
370   
371    INTERFACE pcm_init
372       MODULE PROCEDURE pcm_init
373    END INTERFACE pcm_init
[138]374
[1826]375    INTERFACE pcm_parin
376       MODULE PROCEDURE pcm_parin
[2007]377    END INTERFACE pcm_parin
378
379    INTERFACE pcm_read_plant_canopy_3d
380       MODULE PROCEDURE pcm_read_plant_canopy_3d
381    END INTERFACE pcm_read_plant_canopy_3d
[4360]382
383    INTERFACE pcm_rrd_local
[4517]384       MODULE PROCEDURE pcm_rrd_local_ftn
385       MODULE PROCEDURE pcm_rrd_local_mpi
[4360]386    END INTERFACE pcm_rrd_local
387
[4514]388    INTERFACE pcm_rrd_global
389       MODULE PROCEDURE pcm_rrd_global_ftn
390       MODULE PROCEDURE pcm_rrd_global_mpi
391    END INTERFACE pcm_rrd_global
392
[1826]393    INTERFACE pcm_tendency
394       MODULE PROCEDURE pcm_tendency
395       MODULE PROCEDURE pcm_tendency_ij
396    END INTERFACE pcm_tendency
[1484]397
[4360]398    INTERFACE pcm_wrd_local
399       MODULE PROCEDURE pcm_wrd_local
400    END INTERFACE pcm_wrd_local
[1484]401
[4514]402    INTERFACE pcm_wrd_global
403       MODULE PROCEDURE pcm_wrd_global
404    END INTERFACE pcm_wrd_global
[4360]405
[4514]406
[138]407 CONTAINS
[4342]408 
409 
[2209]410!------------------------------------------------------------------------------!
411! Description:
412! ------------
[3449]413!> Calculation of the plant canopy transpiration rate based on the Jarvis-Stewart
414!> with parametrizations described in Daudet et al. (1999; Agricult. and Forest
415!> Meteorol. 97) and Ngao, Adam and Saudreau (2017;  Agricult. and Forest Meteorol
416!> 237-238). Model functions f1-f4 were adapted from Stewart (1998; Agric.
417!> and Forest. Meteorol. 43) instead, because they are valid for broader intervals
418!> of values. Funcion f4 used in form present in van Wijk et al. (1998;
419!> Tree Physiology 20).
420!>
421!> This subroutine is called from subroutine radiation_interaction
422!> after the calculation of radiation in plant canopy boxes.
423!> (arrays pcbinsw and pcbinlw).
424!>
425!------------------------------------------------------------------------------!
426 SUBROUTINE pcm_calc_transpiration_rate(i, j, k, kk, pcbsw, pcblw, pcbtr, pcblh)
427
[4342]428!
[3449]429!--  input parameters
[4205]430     INTEGER(iwp), INTENT(IN) ::  i, j, k, kk        !< indices of the pc gridbox
431     REAL(wp), INTENT(IN)     ::  pcbsw              !< sw radiation in gridbox (W)
432     REAL(wp), INTENT(IN)     ::  pcblw              !< lw radiation in gridbox (W)
433     REAL(wp), INTENT(OUT)    ::  pcbtr              !< transpiration rate dq/dt (kg/kg/s)
434     REAL(wp), INTENT(OUT)    ::  pcblh              !< latent heat from transpiration dT/dt (K/s)
[3449]435
436!--  variables and parameters for calculation of transpiration rate
[4205]437     REAL(wp)             ::  sat_press, sat_press_d, temp, v_lad
438     REAL(wp)             ::  d_fact, g_b, g_s, wind_speed, evapor_rate
439     REAL(wp)             ::  f1, f2, f3, f4, vpd, rswc, e_eq, e_imp, rad
440     REAL(wp), PARAMETER  ::  gama_psychr = 66.0_wp !< psychrometric constant (Pa/K)
441     REAL(wp), PARAMETER  ::  g_s_max = 0.01        !< maximum stomatal conductivity (m/s)
442     REAL(wp), PARAMETER  ::  m_soil = 0.4_wp       !< soil water content (needs to adjust or take from LSM)
443     REAL(wp), PARAMETER  ::  m_wilt = 0.01_wp      !< wilting point soil water content (needs to adjust or take from LSM)
444     REAL(wp), PARAMETER  ::  m_sat = 0.51_wp       !< saturation soil water content (needs to adjust or take from LSM)
445     REAL(wp), PARAMETER  ::  t2_min = 0.0_wp       !< minimal temperature for calculation of f2
446     REAL(wp), PARAMETER  ::  t2_max = 40.0_wp      !< maximal temperature for calculation of f2
[3449]447
448
449!--  Temperature (deg C)
450     temp = pt(k,j,i) * exner(k) - degc_to_k
451!--  Coefficient for conversion of radiation to grid to radiation to unit leaves surface
[4205]452     v_lad = 1.0_wp / ( MAX( lad_s(kk,j,i), 1.0E-10_wp ) * dx * dy * dz(1) )
[3449]453!--  Magnus formula for the saturation pressure (see Ngao, Adam and Saudreau (2017) eq. 1)
454!--  There are updated formulas available, kept consistent with the rest of the parametrization
455     sat_press = 610.8_wp * exp(17.27_wp * temp/(temp + 237.3_wp))
456!--  Saturation pressure derivative (derivative of the above)
457     sat_press_d = sat_press * 17.27_wp * 237.3_wp / (temp + 237.3_wp)**2
458!--  Wind speed
[3744]459     wind_speed = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 +            &
[4205]460                        ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 +            &
461                        ( 0.5_wp * ( w(k,j,i) + w(k-1,j,i) ) )**2 )
[3449]462!--  Aerodynamic conductivity (Daudet et al. (1999) eq. 14
463     g_b = 0.01_wp * wind_speed + 0.0071_wp
464!--  Radiation flux per leaf surface unit
465     rad = pcbsw * v_lad
466!--  First function for calculation of stomatal conductivity (radiation dependency)
467!--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 17
[4205]468     f1 = rad * (1000.0_wp+42.1_wp) / 1000.0_wp / (rad+42.1_wp)
[3449]469!--  Second function for calculation of stomatal conductivity (temperature dependency)
470!--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 21
[3744]471     f2 = MAX(t2_min, (temp-t2_min) * MAX(0.0_wp,t2_max-temp)**((t2_max-16.9_wp)/(16.9_wp-t2_min)) / &
[3449]472              ((16.9_wp-t2_min) * (t2_max-16.9_wp)**((t2_max-16.9_wp)/(16.9_wp-t2_min))) )
473!--  Water pressure deficit
474!--  Ngao, Adam and Saudreau (2017) eq. 6 but with water vapour partial pressure
475     vpd = max( sat_press - q(k,j,i) * hyp(k) / rd_d_rv, 0._wp )
476!--  Third function for calculation of stomatal conductivity (water pressure deficit dependency)
477!--  Ngao, Adam and Saudreau (2017) Table 1, limited from below according to Stewart (1988)
478!--  The coefficients of the linear dependence should better correspond to broad-leaved trees
479!--  than the coefficients from Stewart (1988) which correspond to conifer trees.
480     vpd = MIN(MAX(vpd,770.0_wp),3820.0_wp)
[4205]481     f3 = -2E-4_wp * vpd + 1.154_wp
[3449]482!--  Fourth function for calculation of stomatal conductivity (soil moisture dependency)
483!--  Residual soil water content
484!--  van Wijk et al. (1998; Tree Physiology 20) eq. 7
485!--  TODO - over LSM surface might be calculated from LSM parameters
486     rswc = ( m_sat - m_soil ) / ( m_sat - m_wilt )
487!--  van Wijk et al. (1998; Tree Physiology 20) eq. 5-6 (it is a reformulation of eq. 22-23 of Stewart(1988))
[4205]488     f4 = MAX(0.0_wp, MIN(1.0_wp - 0.041_wp * EXP(3.2_wp * rswc), 1.0_wp - 0.041_wp))
[3449]489!--  Stomatal conductivity
490!--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 12
491!--  (notation according to Ngao, Adam and Saudreau (2017) and others)
[4205]492     g_s = g_s_max * f1 * f2 * f3 * f4 + 1.0E-10_wp
[3449]493!--  Decoupling factor
494!--  Daudet et al. (1999) eq. 6
[4205]495     d_fact = (sat_press_d / gama_psychr + 2.0_wp ) /                        &
496              (sat_press_d / gama_psychr + 2.0_wp + 2.0_wp * g_b / g_s )
[3449]497!--  Equilibrium evaporation rate
498!--  Daudet et al. (1999) eq. 4
499     e_eq = (pcbsw + pcblw) * v_lad * sat_press_d /                         &
500                 gama_psychr /( sat_press_d / gama_psychr + 2.0_wp ) / l_v
501!--  Imposed evaporation rate
502!--  Daudet et al. (1999) eq. 5
503     e_imp = r_d * pt(k,j,i) * exner(k) / hyp(k) * c_p * g_s * vpd / gama_psychr / l_v
504!--  Evaporation rate
505!--  Daudet et al. (1999) eq. 3
506!--  (evaporation rate is limited to non-negative values)
507     evapor_rate = MAX(d_fact * e_eq + ( 1.0_wp - d_fact ) * e_imp, 0.0_wp)
508!--  Conversion of evaporation rate to q tendency in gridbox
509!--  dq/dt = E * LAD * V_g / (rho_air * V_g)
510     pcbtr = evapor_rate * r_d * pt(k,j,i) * exner(k) * lad_s(kk,j,i) / hyp(k)  !-- = dq/dt
511!--  latent heat from evaporation
512     pcblh = pcbtr * lv_d_cp  !-- = - dT/dt
513
514 END SUBROUTINE pcm_calc_transpiration_rate
515
516
517!------------------------------------------------------------------------------!
518! Description:
519! ------------
[2209]520!> Check data output for plant canopy model
521!------------------------------------------------------------------------------!
522 SUBROUTINE pcm_check_data_output( var, unit )
[4279]523
[2209]524    CHARACTER (LEN=*) ::  unit  !<
525    CHARACTER (LEN=*) ::  var   !<
526
527
528    SELECT CASE ( TRIM( var ) )
529
530       CASE ( 'pcm_heatrate' )
[4514]531!
532!--       Output of heatrate can be only done if it is explicitely set by cthf,
533!--       or parametrized by absorption of radiation. The latter, however, is
534!--       only available if radiation_interactions are on. Note, these are
535!--       enabled if land-surface or urban-surface is switched-on. Using
536!--       radiation_interactions_on directly is not possible since it belongs
537!--       to the radition_model, which in turn depends on the plant-canopy model,
538!--       creating circular dependencies.
539          IF ( cthf == 0.0_wp  .AND. (  .NOT.  urban_surface  .AND.             &
540                                        .NOT.  land_surface ) )  THEN
541             message_string = 'output of "' // TRIM( var ) // '" requi' //      &
[2768]542                              'res setting of parameter cthf /= 0.0'
[4515]543             CALL message( 'pcm_check_data_output', 'PA0718', 1, 2, 0, 6, 0 )
[2768]544          ENDIF
[2209]545          unit = 'K s-1'
546   
[3014]547       CASE ( 'pcm_transpirationrate' )
548          unit = 'kg kg-1 s-1'
549
[3449]550       CASE ( 'pcm_latentrate' )
551          unit = 'K s-1'
552
[4770]553       CASE ( 'pcm_bad', 'pcm_lad' )
[2209]554          unit = 'm2 m-3'
555
556
557       CASE DEFAULT
558          unit = 'illegal'
559
560    END SELECT
561
562
563 END SUBROUTINE pcm_check_data_output
564 
565 
[1826]566!------------------------------------------------------------------------------!
567! Description:
568! ------------
569!> Check parameters routine for plant canopy model
570!------------------------------------------------------------------------------!
571    SUBROUTINE pcm_check_parameters
[4342]572           
573       IF ( ocean_mode )  THEN
574          message_string = 'plant_canopy = .TRUE. is not allowed in the '//    &
575                           'ocean'
576          CALL message( 'pcm_check_parameters', 'PA0696', 1, 2, 0, 6, 0 )
577       ENDIF
578           
[1826]579       IF ( canopy_drag_coeff == 0.0_wp )  THEN
580          message_string = 'plant_canopy = .TRUE. requires a non-zero drag '// &
[3046]581                           'coefficient & given value is canopy_drag_coeff = 0.0'
[2768]582          CALL message( 'pcm_check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
[1826]583       ENDIF
[4279]584
[3045]585       IF ( ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad == 9999999.9_wp ) .OR.&
[1826]586              beta_lad /= 9999999.9_wp   .AND.  alpha_lad == 9999999.9_wp )  THEN
587          message_string = 'using the beta function for the construction ' //  &
588                           'of the leaf area density profile requires '    //  &
589                           'both alpha_lad and beta_lad to be /= 9999999.9'
[2768]590          CALL message( 'pcm_check_parameters', 'PA0118', 1, 2, 0, 6, 0 )
[1826]591       ENDIF
[4279]592
[1826]593       IF ( calc_beta_lad_profile  .AND.  lai_beta == 0.0_wp )  THEN
594          message_string = 'using the beta function for the construction ' //  &
595                           'of the leaf area density profile requires '    //  &
596                           'a non-zero lai_beta, but given value is '      //  &
597                           'lai_beta = 0.0'
[2768]598          CALL message( 'pcm_check_parameters', 'PA0119', 1, 2, 0, 6, 0 )
[1826]599       ENDIF
600
601       IF ( calc_beta_lad_profile  .AND.  lad_surface /= 0.0_wp )  THEN
[2274]602          message_string = 'simultaneous setting of alpha_lad /= 9999999.9 '// &
603                           'combined with beta_lad /= 9999999.9 '           // &
[1826]604                           'and lad_surface /= 0.0 is not possible, '       // &
605                           'use either vertical gradients or the beta '     // &
606                           'function for the construction of the leaf area '// &
607                           'density profile'
[2768]608          CALL message( 'pcm_check_parameters', 'PA0120', 1, 2, 0, 6, 0 )
[1826]609       ENDIF
610
[3274]611       IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
[1826]612          message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' // &
613                          ' seifert_beheng'
[2768]614          CALL message( 'pcm_check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
[1826]615       ENDIF
616
617    END SUBROUTINE pcm_check_parameters 
618 
619
[138]620!------------------------------------------------------------------------------!
[2209]621!
[1484]622! Description:
623! ------------
[4127]624!> Subroutine for averaging 3D data
[2209]625!------------------------------------------------------------------------------!
[4216]626 SUBROUTINE pcm_3d_data_averaging( mode, variable )
[4127]627
628    CHARACTER (LEN=*) ::  mode    !<
629    CHARACTER (LEN=*) :: variable !<
630
631    INTEGER(iwp) ::  i            !<
632    INTEGER(iwp) ::  j            !<
633    INTEGER(iwp) ::  k            !<
634
635
636    IF ( mode == 'allocate' )  THEN
637
638       SELECT CASE ( TRIM( variable ) )
639
640          CASE ( 'pcm_heatrate' )
641             IF ( .NOT. ALLOCATED( pcm_heatrate_av ) )  THEN
642                ALLOCATE( pcm_heatrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
643             ENDIF
644             pcm_heatrate_av = 0.0_wp
645
646
647          CASE ( 'pcm_latentrate' )
648             IF ( .NOT. ALLOCATED( pcm_latentrate_av ) )  THEN
649                ALLOCATE( pcm_latentrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
650             ENDIF
651             pcm_latentrate_av = 0.0_wp
652
653
654          CASE ( 'pcm_transpirationrate' )
655             IF ( .NOT. ALLOCATED( pcm_transpirationrate_av ) )  THEN
656                ALLOCATE( pcm_transpirationrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
657             ENDIF
658             pcm_transpirationrate_av = 0.0_wp
659
660          CASE DEFAULT
661             CONTINUE
662
663       END SELECT
664
665    ELSEIF ( mode == 'sum' )  THEN
666
667       SELECT CASE ( TRIM( variable ) )
668
669          CASE ( 'pcm_heatrate' )
670             IF ( ALLOCATED( pcm_heatrate_av ) ) THEN
671                DO  i = nxl, nxr
672                   DO  j = nys, nyn
673                      IF ( pch_index_ji(j,i) /= 0 )  THEN
674                         DO  k = 0, pch_index_ji(j,i)
[4341]675                            pcm_heatrate_av(k,j,i) = pcm_heatrate_av(k,j,i) + pcm_heating_rate(k,j,i)
[4127]676                         ENDDO
677                      ENDIF
678                   ENDDO
679                ENDDO
680             ENDIF
681
682
683          CASE ( 'pcm_latentrate' )
684             IF ( ALLOCATED( pcm_latentrate_av ) ) THEN
685                DO  i = nxl, nxr
686                   DO  j = nys, nyn
687                      IF ( pch_index_ji(j,i) /= 0 )  THEN
688                         DO  k = 0, pch_index_ji(j,i)
[4341]689                            pcm_latentrate_av(k,j,i) = pcm_latentrate_av(k,j,i) + pcm_latent_rate(k,j,i)
[4127]690                         ENDDO
691                      ENDIF
692                   ENDDO
693                ENDDO
694             ENDIF
695
696
697          CASE ( 'pcm_transpirationrate' )
698             IF ( ALLOCATED( pcm_transpirationrate_av ) ) THEN
699                DO  i = nxl, nxr
700                   DO  j = nys, nyn
701                      IF ( pch_index_ji(j,i) /= 0 )  THEN
702                         DO  k = 0, pch_index_ji(j,i)
[4341]703                            pcm_transpirationrate_av(k,j,i) = pcm_transpirationrate_av(k,j,i) + pcm_transpiration_rate(k,j,i)
[4127]704                         ENDDO
705                      ENDIF
706                   ENDDO
707                ENDDO
708             ENDIF
709
710          CASE DEFAULT
711             CONTINUE
712
713       END SELECT
714
715    ELSEIF ( mode == 'average' )  THEN
716
717       SELECT CASE ( TRIM( variable ) )
718
719          CASE ( 'pcm_heatrate' )
720             IF ( ALLOCATED( pcm_heatrate_av ) ) THEN
721                DO  i = nxlg, nxrg
722                   DO  j = nysg, nyng
723                      IF ( pch_index_ji(j,i) /= 0 )  THEN
724                         DO  k = 0, pch_index_ji(j,i)
725                            pcm_heatrate_av(k,j,i) = pcm_heatrate_av(k,j,i)                 &
726                                                     / REAL( average_count_3d, KIND=wp )
727                         ENDDO
728                      ENDIF
729                   ENDDO
730                ENDDO
731             ENDIF
732
733
734          CASE ( 'pcm_latentrate' )
735             IF ( ALLOCATED( pcm_latentrate_av ) ) THEN
736                DO  i = nxlg, nxrg
737                   DO  j = nysg, nyng
738                      IF ( pch_index_ji(j,i) /= 0 )  THEN
739                         DO  k = 0, pch_index_ji(j,i)
740                            pcm_latentrate_av(k,j,i) = pcm_latentrate_av(k,j,i)              &
741                                                       / REAL( average_count_3d, KIND=wp )
742                         ENDDO
743                      ENDIF
744                   ENDDO
745                ENDDO
746             ENDIF
747
748
749          CASE ( 'pcm_transpirationrate' )
750             IF ( ALLOCATED( pcm_transpirationrate_av ) ) THEN
751                DO  i = nxlg, nxrg
752                   DO  j = nysg, nyng
753                      IF ( pch_index_ji(j,i) /= 0 )  THEN
754                         DO  k = 0, pch_index_ji(j,i)
755                            pcm_transpirationrate_av(k,j,i) = pcm_transpirationrate_av(k,j,i)  &
756                                                              / REAL( average_count_3d, KIND=wp )
757                         ENDDO
758                      ENDIF
759                   ENDDO
760                ENDDO
761             ENDIF
762
763       END SELECT
764
765    ENDIF
766
[4216]767 END SUBROUTINE pcm_3d_data_averaging
[4127]768
769!------------------------------------------------------------------------------!
770!
771! Description:
772! ------------
773!> Subroutine defining 3D output variables.
774!> Note, 3D plant-canopy output has it's own vertical output dimension, meaning
775!> that 3D output is relative to the model surface now rather than at the actual
776!> grid point where the plant canopy is located.
777!------------------------------------------------------------------------------!
[3014]778 SUBROUTINE pcm_data_output_3d( av, variable, found, local_pf, fill_value,     &
779                                nzb_do, nzt_do )
[2209]780
[4216]781    CHARACTER (LEN=*) ::  variable !< treated variable
[2209]782
[4216]783    INTEGER(iwp) ::  av     !< flag indicating instantaneous or averaged data output
784    INTEGER(iwp) ::  i      !< grid index x-direction
785    INTEGER(iwp) ::  j      !< grid index y-direction
786    INTEGER(iwp) ::  k      !< grid index z-direction
[3014]787    INTEGER(iwp) ::  nzb_do !< lower limit of the data output (usually 0)
788    INTEGER(iwp) ::  nzt_do !< vertical upper limit of the data output (usually nz_do3d)
[2209]789
[4216]790    LOGICAL      ::  found  !< flag indicating if variable is found
[2209]791
[4216]792    REAL(wp)     ::  fill_value !< fill value
[4768]793    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< data output array
[2209]794
795
796    found = .TRUE.
797
[2696]798    local_pf = REAL( fill_value, KIND = 4 )
[2209]799
800    SELECT CASE ( TRIM( variable ) )
[4216]801!
802!--    Note, to save memory arrays for heating are allocated from 0:pch_index.
803!--    Thus, output must be relative to these array indices. Further, check
804!--    whether the output is within the vertical output range,
[4220]805!--    i.e. nzb_do:nzt_do, which is necessary as local_pf is only allocated
806!--    for this index space. Note, plant-canopy output has a separate
807!--    vertical output coordinate zlad, so that output is mapped down to the
808!--    surface.
[4127]809       CASE ( 'pcm_heatrate' )
810          IF ( av == 0 )  THEN
811             DO  i = nxl, nxr
812                DO  j = nys, nyn
[4360]813                   DO  k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
[4341]814                      local_pf(i,j,k) = pcm_heating_rate(k,j,i)
[4216]815                   ENDDO
[4127]816                ENDDO
817             ENDDO
818          ELSE
819             DO  i = nxl, nxr
820                DO  j = nys, nyn
[4360]821                   DO  k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
[4220]822                      local_pf(i,j,k) = pcm_heatrate_av(k,j,i)
[4127]823                   ENDDO
824                ENDDO
825             ENDDO
826          ENDIF
[3449]827
828       CASE ( 'pcm_latentrate' )
[4127]829          IF ( av == 0 )  THEN
830             DO  i = nxl, nxr
831                DO  j = nys, nyn
[4360]832                   DO  k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
[4341]833                      local_pf(i,j,k) = pcm_latent_rate(k,j,i)
[4216]834                   ENDDO
[4127]835                ENDDO
836             ENDDO
837          ELSE
838             DO  i = nxl, nxr
839                DO  j = nys, nyn
[4360]840                   DO  k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
[4220]841                      local_pf(i,j,k) = pcm_latentrate_av(k,j,i)
[4127]842                   ENDDO
843                ENDDO
844             ENDDO
845          ENDIF
[3449]846
[4127]847       CASE ( 'pcm_transpirationrate' )
848          IF ( av == 0 )  THEN
849             DO  i = nxl, nxr
850                DO  j = nys, nyn
[4360]851                   DO  k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
[4341]852                      local_pf(i,j,k) = pcm_transpiration_rate(k,j,i)
[4216]853                   ENDDO
[4127]854                ENDDO
855             ENDDO
856          ELSE
857             DO  i = nxl, nxr
858                DO  j = nys, nyn
[4360]859                   DO  k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
[4220]860                      local_pf(i,j,k) = pcm_transpirationrate_av(k,j,i)
[4127]861                   ENDDO
862                ENDDO
863             ENDDO
864          ENDIF
865
866       CASE ( 'pcm_lad' )
867          IF ( av == 0 )  THEN
868             DO  i = nxl, nxr
869                DO  j = nys, nyn
[4360]870                   DO  k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
[4220]871                      local_pf(i,j,k) = lad_s(k,j,i)
[4216]872                   ENDDO
[4127]873                ENDDO
874             ENDDO
875          ENDIF
876
[4770]877       CASE ( 'pcm_bad' )
878          IF ( av == 0 )  THEN
879             DO  i = nxl, nxr
880                DO  j = nys, nyn
881                   DO  k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
882                      local_pf(i,j,k) = bad_s(k,j,i)
883                   ENDDO
884                ENDDO
885             ENDDO
886          ENDIF
887
[2209]888       CASE DEFAULT
889          found = .FALSE.
890
891    END SELECT
892
893 END SUBROUTINE pcm_data_output_3d 
894         
895!------------------------------------------------------------------------------!
896!
897! Description:
898! ------------
899!> Subroutine defining appropriate grid for netcdf variables.
900!> It is called from subroutine netcdf.
901!------------------------------------------------------------------------------!
902 SUBROUTINE pcm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
903
904     CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
905     LOGICAL, INTENT(OUT)           ::  found       !<
906     CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
907     CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
908     CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
909
910     found  = .TRUE.
911
912!
[4342]913!--  Check for the grid. zpc is zu(nzb:nzb+pch_index)
[2209]914     SELECT CASE ( TRIM( var ) )
915
[4770]916        CASE ( 'pcm_heatrate', 'pcm_bad', 'pcm_lad', 'pcm_transpirationrate', 'pcm_latentrate')
[2209]917           grid_x = 'x'
918           grid_y = 'y'
[4127]919           grid_z = 'zpc'
[2209]920
921        CASE DEFAULT
922           found  = .FALSE.
923           grid_x = 'none'
924           grid_y = 'none'
925           grid_z = 'none'
926     END SELECT
927
928 END SUBROUTINE pcm_define_netcdf_grid
929 
930 
931!------------------------------------------------------------------------------!
932! Description:
933! ------------
[1826]934!> Header output for plant canopy model
935!------------------------------------------------------------------------------!
[4362]936 SUBROUTINE pcm_header ( io )
[1826]937 
[4362]938    CHARACTER (LEN=10) ::  coor_chr            !<
[1826]939
[4362]940    CHARACTER (LEN=86) ::  coordinates         !<
941    CHARACTER (LEN=86) ::  gradients           !<
942    CHARACTER (LEN=86) ::  leaf_area_density   !<
943    CHARACTER (LEN=86) ::  slices              !<
[1826]944 
[4362]945    INTEGER(iwp) :: i                !<
946    INTEGER(iwp),  INTENT(IN) ::  io !< Unit of the output file
947    INTEGER(iwp) :: k                !<       
948 
949    REAL(wp) ::  canopy_height       !< canopy height (in m)
950   
951    canopy_height = zw(pch_index)
[1826]952
[4362]953    WRITE ( io, 1 )  canopy_mode, canopy_height, pch_index,                    &
954                       canopy_drag_coeff                                       
955    IF ( passive_scalar )  THEN                                               
956       WRITE ( io, 2 )  leaf_scalar_exch_coeff,                                &
957                          leaf_surface_conc
958    ENDIF
[1826]959
960!
[4362]961!   Heat flux at the top of vegetation
962    WRITE ( io, 3 )  cthf
[1826]963
964!
[4362]965!   Leaf area density profile, calculated either from given vertical
966!   gradients or from beta probability density function.
967    IF (  .NOT.  calc_beta_lad_profile )  THEN
[1826]968
[4362]969!      Building output strings, starting with surface value
970       WRITE ( leaf_area_density, '(F7.4)' )  lad_surface
971       gradients = '------'
972       slices = '     0'
973       coordinates = '   0.0'
974       DO i = 1, UBOUND(lad_vertical_gradient_level_ind, DIM=1)
975          IF  ( lad_vertical_gradient_level_ind(i) /= -9999 ) THEN
[1826]976
[4362]977             WRITE (coor_chr,'(F7.2)') lad(lad_vertical_gradient_level_ind(i))
978             leaf_area_density = TRIM( leaf_area_density ) // ' ' // TRIM( coor_chr )
[1826]979
[4362]980             WRITE (coor_chr,'(F7.2)') lad_vertical_gradient(i)
981             gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
[1826]982
[4362]983             WRITE (coor_chr,'(I7)') lad_vertical_gradient_level_ind(i)
984             slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
[1826]985
[4362]986             WRITE (coor_chr,'(F7.1)') lad_vertical_gradient_level(i)
987             coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
988          ELSE
989             EXIT
990          ENDIF
991       ENDDO
[1826]992
[4362]993       WRITE ( io, 4 )  TRIM( coordinates ), TRIM( leaf_area_density ),        &
994                          TRIM( gradients ), TRIM( slices )
[1826]995
[4362]996    ELSE
997   
998       WRITE ( leaf_area_density, '(F7.4)' )  lad_surface
999       coordinates = '   0.0'
[1826]1000       
[4362]1001       DO  k = 1, pch_index
[1826]1002
[4362]1003          WRITE (coor_chr,'(F7.2)')  lad(k)
1004          leaf_area_density = TRIM( leaf_area_density ) // ' ' //              &
1005                              TRIM( coor_chr )                                 
1006                                                                               
1007          WRITE (coor_chr,'(F7.1)')  zu(k)                                     
1008          coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )       
1009                                                                               
1010       ENDDO                                                                   
1011                                                                               
1012       WRITE ( io, 5 ) TRIM( coordinates ), TRIM( leaf_area_density ),         &
1013                       alpha_lad, beta_lad, lai_beta
[1826]1014
[4362]1015    ENDIF 
[1826]1016
[4362]10171   FORMAT (//' Vegetation canopy (drag) model:'/                              &
[1826]1018              ' ------------------------------'//                              &
1019              ' Canopy mode: ', A /                                            &
1020              ' Canopy height: ',F6.2,'m (',I4,' grid points)' /               &
1021              ' Leaf drag coefficient: ',F6.2 /)
[4362]10222   FORMAT (/ ' Scalar exchange coefficient: ',F6.2 /                          &
[1826]1023              ' Scalar concentration at leaf surfaces in kg/m**3: ',F6.2 /)
[4362]10243   FORMAT (' Predefined constant heatflux at the top of the vegetation: ',    &
1025             F6.2, ' K m/s')
10264   FORMAT (/ ' Characteristic levels of the leaf area density:'//             &
[1826]1027              ' Height:              ',A,'  m'/                                &
1028              ' Leaf area density:   ',A,'  m**2/m**3'/                        &
1029              ' Gradient:            ',A,'  m**2/m**4'/                        &
1030              ' Gridpoint:           ',A)
[4362]10315   FORMAT (//' Characteristic levels of the leaf area density and coefficients:'&
[1826]1032          //  ' Height:              ',A,'  m'/                                &
1033              ' Leaf area density:   ',A,'  m**2/m**3'/                        &
1034              ' Coefficient alpha: ',F6.2 /                                    &
1035              ' Coefficient beta: ',F6.2 /                                     &
1036              ' Leaf area index: ',F6.2,'  m**2/m**2' /)   
1037       
1038    END SUBROUTINE pcm_header
1039 
1040 
1041!------------------------------------------------------------------------------!
1042! Description:
1043! ------------
[1682]1044!> Initialization of the plant canopy model
[138]1045!------------------------------------------------------------------------------!
[1826]1046    SUBROUTINE pcm_init
[1484]1047
[4457]1048       USE exchange_horiz_mod,                                                    &
1049           ONLY:  exchange_horiz
1050
[2007]1051       INTEGER(iwp) ::  i   !< running index
1052       INTEGER(iwp) ::  j   !< running index
1053       INTEGER(iwp) ::  k   !< running index
[2232]1054       INTEGER(iwp) ::  m   !< running index
[1484]1055
[4381]1056       LOGICAL      ::  lad_on_top = .FALSE.  !< dummy flag to indicate that LAD is defined on a building roof
[4770]1057       LOGICAL      ::  bad_on_top = .FALSE.  !< dummy flag to indicate that BAD is defined on a building roof
[4381]1058
[4258]1059       REAL(wp) ::  canopy_height   !< canopy height for lad-profile construction
[2007]1060       REAL(wp) ::  gradient        !< gradient for lad-profile construction
[4258]1061       REAL(wp) ::  int_bpdf        !< vertical integral for lad-profile construction     
1062       REAL(wp) ::  lad_max         !< maximum LAD value in the model domain, used to perform a check
[3241]1063
[3885]1064       IF ( debug_output )  CALL debug_message( 'pcm_init', 'start' )
[1484]1065!
1066!--    Allocate one-dimensional arrays for the computation of the
1067!--    leaf area density (lad) profile
1068       ALLOCATE( lad(0:nz+1), pre_lad(0:nz+1) )
1069       lad = 0.0_wp
1070       pre_lad = 0.0_wp
1071
1072!
[1826]1073!--    Set flag that indicates that the lad-profile shall be calculated by using
1074!--    a beta probability density function
1075       IF ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad /= 9999999.9_wp )  THEN
1076          calc_beta_lad_profile = .TRUE.
1077       ENDIF
1078       
1079       
1080!
[1484]1081!--    Compute the profile of leaf area density used in the plant
1082!--    canopy model. The profile can either be constructed from
1083!--    prescribed vertical gradients of the leaf area density or by
1084!--    using a beta probability density function (see e.g. Markkanen et al.,
1085!--    2003: Boundary-Layer Meteorology, 106, 437-459)
1086       IF (  .NOT.  calc_beta_lad_profile )  THEN   
1087
1088!
1089!--       Use vertical gradients for lad-profile construction   
1090          i = 1
1091          gradient = 0.0_wp
1092
[4342]1093          lad(0) = lad_surface
1094          lad_vertical_gradient_level_ind(1) = 0
[1484]1095
[4342]1096          DO k = 1, pch_index
1097             IF ( i < 11 )  THEN
1098                IF ( lad_vertical_gradient_level(i) < zu(k)  .AND.          &
1099                     lad_vertical_gradient_level(i) >= 0.0_wp )  THEN
1100                   gradient = lad_vertical_gradient(i)
1101                   lad_vertical_gradient_level_ind(i) = k - 1
1102                   i = i + 1
[1484]1103                ENDIF
[4342]1104             ENDIF
1105             IF ( gradient /= 0.0_wp )  THEN
1106                IF ( k /= 1 )  THEN
1107                   lad(k) = lad(k-1) + dzu(k) * gradient
[1484]1108                ELSE
[4342]1109                   lad(k) = lad_surface + dzu(k) * gradient
[1484]1110                ENDIF
[4342]1111             ELSE
1112                lad(k) = lad(k-1)
1113             ENDIF
1114          ENDDO
[1484]1115
1116!
1117!--       In case of no given leaf area density gradients, choose a vanishing
1118!--       gradient. This information is used for the HEADER and the RUN_CONTROL
1119!--       file.
1120          IF ( lad_vertical_gradient_level(1) == -9999999.9_wp )  THEN
1121             lad_vertical_gradient_level(1) = 0.0_wp
1122          ENDIF
1123
1124       ELSE
1125
1126!
1127!--       Use beta function for lad-profile construction
1128          int_bpdf = 0.0_wp
[3065]1129          canopy_height = zw(pch_index)
[1484]1130
[2232]1131          DO k = 0, pch_index
[1484]1132             int_bpdf = int_bpdf +                                             &
[1826]1133                      ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) *  &
1134                      ( ( 1.0_wp - ( zw(k) / canopy_height ) )**(              &
1135                          beta_lad-1.0_wp ) )                                  &
1136                      * ( ( zw(k+1)-zw(k) ) / canopy_height ) )
[1484]1137          ENDDO
1138
1139!
1140!--       Preliminary lad profile (defined on w-grid)
[2232]1141          DO k = 0, pch_index
[1826]1142             pre_lad(k) =  lai_beta *                                          &
1143                        ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) )  &
1144                        * ( ( 1.0_wp - ( zw(k) / canopy_height ) )**(          &
1145                              beta_lad-1.0_wp ) ) / int_bpdf                   &
1146                        ) / canopy_height
[1484]1147          ENDDO
1148
1149!
1150!--       Final lad profile (defined on scalar-grid level, since most prognostic
1151!--       quantities are defined there, hence, less interpolation is required
1152!--       when calculating the canopy tendencies)
1153          lad(0) = pre_lad(0)
[2232]1154          DO k = 1, pch_index
[1484]1155             lad(k) = 0.5 * ( pre_lad(k-1) + pre_lad(k) )
[4302]1156          ENDDO
[1484]1157
1158       ENDIF
1159
1160!
[4770]1161!--    Allocate 3D-array for the leaf-area density (lad_s) as well as
1162!--    for basal-area densitiy (bad_s). Note, by default bad_s is zero.
[1484]1163       ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4770]1164       ALLOCATE( bad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1165       bad_s = 0.0_wp
[1484]1166!
1167!--    Initialization of the canopy coverage in the model domain:
[4341]1168!--    Setting the parameter canopy_mode = 'homogeneous' initializes a canopy, which
[1484]1169!--    fully covers the domain surface
1170       SELECT CASE ( TRIM( canopy_mode ) )
1171
[4341]1172          CASE( 'homogeneous' )
[1484]1173
1174             DO  i = nxlg, nxrg
1175                DO  j = nysg, nyng
1176                   lad_s(:,j,i) = lad(:)
1177                ENDDO
1178             ENDDO
1179
[4341]1180          CASE ( 'read_from_file' )
[2007]1181!
[4362]1182!--          Read plant canopy
1183             IF ( input_pids_static )  THEN
1184!
1185!--             Open the static input file
1186#if defined( __netcdf )
1187                CALL open_read_file( TRIM( input_file_static ) //              &
1188                                     TRIM( coupling_char ),                    &
1189                                     pids_id )
1190
1191                CALL inquire_num_variables( pids_id, num_var_pids )
1192!
1193!--             Allocate memory to store variable names and read them
1194                ALLOCATE( vars_pids(1:num_var_pids) )
1195                CALL inquire_variable_names( pids_id, vars_pids )
1196!
1197!--             Read leaf area density - resolved vegetation
1198                IF ( check_existence( vars_pids, 'lad' ) )  THEN
1199                   leaf_area_density_f%from_file = .TRUE.
1200                   CALL get_attribute( pids_id, char_fill,                     &
1201                                       leaf_area_density_f%fill,               &
1202                                       .FALSE., 'lad' )
1203!
1204!--                Inquire number of vertical vegetation layer
1205                   CALL get_dimension_length( pids_id,                         &
1206                                              leaf_area_density_f%nz,          &
1207                                              'zlad' )
1208!
1209!--                Allocate variable for leaf-area density
1210                   ALLOCATE( leaf_area_density_f%var                           &
1211                                                (0:leaf_area_density_f%nz-1,   &
1212                                                 nys:nyn,nxl:nxr) )
1213
1214                   CALL get_variable( pids_id, 'lad', leaf_area_density_f%var, &
1215                                      nxl, nxr, nys, nyn,                      &
1216                                      0, leaf_area_density_f%nz-1 )
1217
1218                ELSE
1219                   leaf_area_density_f%from_file = .FALSE.
1220                ENDIF
1221!
1222!--             Read basal area density - resolved vegetation
1223                IF ( check_existence( vars_pids, 'bad' ) )  THEN
1224                   basal_area_density_f%from_file = .TRUE.
1225                   CALL get_attribute( pids_id, char_fill,                     &
1226                                       basal_area_density_f%fill,              &
1227                                       .FALSE., 'bad' )
1228!
1229!--                Inquire number of vertical vegetation layer
1230                   CALL get_dimension_length( pids_id,                         &
1231                                              basal_area_density_f%nz,         & 
1232                                              'zlad' )
1233!
1234!--                Allocate variable
1235                   ALLOCATE( basal_area_density_f%var                          &
1236                                              (0:basal_area_density_f%nz-1,    &
1237                                               nys:nyn,nxl:nxr) )
1238
1239                   CALL get_variable( pids_id, 'bad', basal_area_density_f%var,&
1240                                      nxl, nxr, nys, nyn,                      &
1241                                      0,  basal_area_density_f%nz-1 )
1242                ELSE
1243                   basal_area_density_f%from_file = .FALSE.
1244                ENDIF
1245!
1246!--             Read root area density - resolved vegetation
1247                IF ( check_existence( vars_pids, 'root_area_dens_r' ) )  THEN
1248                   root_area_density_lad_f%from_file = .TRUE.
1249                   CALL get_attribute( pids_id, char_fill,                     &
1250                                       root_area_density_lad_f%fill,           &
1251                                       .FALSE., 'root_area_dens_r' )
1252!
1253!--                Inquire number of vertical soil layers
1254                   CALL get_dimension_length( pids_id,                         &
1255                                              root_area_density_lad_f%nz,      &
1256                                              'zsoil' )
1257!
1258!--                Allocate variable
1259                   ALLOCATE( root_area_density_lad_f%var                       &
1260                                               (0:root_area_density_lad_f%nz-1,&
1261                                                nys:nyn,nxl:nxr) )
1262
1263                   CALL get_variable( pids_id, 'root_area_dens_r',             &
1264                                      root_area_density_lad_f%var,             &
1265                                      nxl, nxr, nys, nyn,                      &
1266                                      0,  root_area_density_lad_f%nz-1 )
1267                ELSE
1268                   root_area_density_lad_f%from_file = .FALSE.
1269                ENDIF
1270
1271                DEALLOCATE( vars_pids )
1272!
[4363]1273!--             Finally, close the input file and deallocate temporary array
1274                CALL close_input_file( pids_id )
[4362]1275#endif
[4363]1276             ENDIF
1277
[4362]1278!
[2696]1279!--          Initialize LAD with data from file. If LAD is given in NetCDF file,
1280!--          use these values, else take LAD profiles from ASCII file.
1281!--          Please note, in NetCDF file LAD is only given up to the maximum
1282!--          canopy top, indicated by leaf_area_density_f%nz. 
1283             lad_s = 0.0_wp
1284             IF ( leaf_area_density_f%from_file )  THEN
1285!
1286!--             Set also pch_index, used to be the upper bound of the vertical
1287!--             loops. Therefore, use the global top of the canopy layer.
1288                pch_index = leaf_area_density_f%nz - 1
1289
1290                DO  i = nxl, nxr
1291                   DO  j = nys, nyn
1292                      DO  k = 0, leaf_area_density_f%nz - 1
[3864]1293                         IF ( leaf_area_density_f%var(k,j,i) /=                &
1294                              leaf_area_density_f%fill )                       &
[2696]1295                         lad_s(k,j,i) = leaf_area_density_f%var(k,j,i)
1296                      ENDDO
[4302]1297!
1298!--                   Check if resolved vegetation is mapped onto buildings.
[4314]1299!--                   In general, this is allowed and also meaningful, e.g.
1300!--                   when trees carry across roofs. However, due to the
1301!--                   topography filtering, new building grid points can emerge
1302!--                   at location where also plant canopy is defined. As a
1303!--                   result, plant canopy is mapped on top of roofs, with
1304!--                   siginficant impact on the downstream flow field and the
1305!--                   nearby surface radiation. In order to avoid that
1306!--                   plant canopy is mistakenly mapped onto building roofs,
1307!--                   check for building grid points (bit 6) that emerge from
1308!--                   the filtering (bit 4) and set LAD to zero at these
1309!--                   artificially  created building grid points. This case,
1310!--                   an informative message is given.
[4302]1311                      IF ( ANY( lad_s(:,j,i) /= 0.0_wp )  .AND.                &
[4346]1312                           ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) ) .AND.  &
1313                           ANY( BTEST( wall_flags_total_0(:,j,i), 4 ) ) )  THEN
[4302]1314                         lad_s(:,j,i) = 0.0_wp
[4381]1315                         lad_on_top   = .TRUE.
1316                      ENDIF
1317                   ENDDO
1318                ENDDO
1319#if defined( __parallel )
1320               CALL MPI_ALLREDUCE( MPI_IN_PLACE, lad_on_top, 1, MPI_LOGICAL,  &
1321                                   MPI_LOR, comm2d, ierr)
1322#endif
1323                IF ( lad_on_top )  THEN
[4770]1324                   WRITE( message_string, * )                                                      &
1325                                        'Resolved plant-canopy is defined on top of an ' //        &
1326                                        'artificially created building grid point(s) '//           &
1327                                        'the filtering) - LAD/BAD profile is omitted at this / ' //&
[4381]1328                                        'these grid point(s).'
1329                   CALL message( 'pcm_init', 'PA0313', 0, 0, 0, 6, 0 )
1330                ENDIF
[2696]1331                CALL exchange_horiz( lad_s, nbgp )
1332!
1333!            ASCII file
[4342]1334!--          Initialize canopy parameters canopy_drag_coeff,
1335!--          leaf_scalar_exch_coeff, leaf_surface_conc
[2007]1336!--          from file which contains complete 3D data (separate vertical profiles for
1337!--          each location).
[2696]1338             ELSE
1339                CALL pcm_read_plant_canopy_3d
1340             ENDIF
[4770]1341!
1342!--          Initialize LAD with data from file. If LAD is given in NetCDF file,
1343!--          use these values, else take LAD profiles from ASCII file.
1344!--          Please note, in NetCDF file LAD is only given up to the maximum
1345!--          canopy top, indicated by basal_area_density_f%nz. 
1346             bad_s = 0.0_wp
1347             IF ( basal_area_density_f%from_file )  THEN
1348!
1349!--             Set also pch_index, used to be the upper bound of the vertical
1350!--             loops. Therefore, use the global top of the canopy layer.
1351                pch_index = basal_area_density_f%nz - 1
[2007]1352
[4770]1353                DO  i = nxl, nxr
1354                   DO  j = nys, nyn
1355                      DO  k = 0, basal_area_density_f%nz - 1
1356                         IF ( basal_area_density_f%var(k,j,i) /= basal_area_density_f%fill )       &
1357                            bad_s(k,j,i) = basal_area_density_f%var(k,j,i)
1358                      ENDDO
1359!
1360!--                   Check if resolved vegetation is mapped onto buildings.
1361!--                   Please see comment for leaf_area density
1362                      IF ( ANY( bad_s(:,j,i) /= 0.0_wp )  .AND.                &
1363                           ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) ) .AND.  &
1364                           ANY( BTEST( wall_flags_total_0(:,j,i), 4 ) ) )  THEN
1365                         bad_s(:,j,i) = 0.0_wp
1366                         bad_on_top   = .TRUE.
1367                      ENDIF
1368                   ENDDO
1369                ENDDO
1370#if defined( __parallel )
1371               CALL MPI_ALLREDUCE( MPI_IN_PLACE, bad_on_top, 1, MPI_LOGICAL,  &
1372                                   MPI_LOR, comm2d, ierr)
1373#endif
1374                IF ( bad_on_top )  THEN
1375                   WRITE( message_string, * )                                                      &
1376                                        'Resolved plant-canopy is defined on top of an ' //        &
1377                                        'artificially created building grid point(s) '//           &
1378                                        'the filtering) - LAD/BAD profile is omitted at this / ' //&
1379                                        'these grid point(s).'
1380                   CALL message( 'pcm_init', 'PA0313', 0, 0, 0, 6, 0 )
1381                ENDIF
1382                CALL exchange_horiz( bad_s, nbgp )
1383             ENDIF
1384
[1484]1385          CASE DEFAULT
1386!
[2007]1387!--          The DEFAULT case is reached either if the parameter
1388!--          canopy mode contains a wrong character string or if the
1389!--          user has coded a special case in the user interface.
1390!--          There, the subroutine user_init_plant_canopy checks
1391!--          which of these two conditions applies.
1392             CALL user_init_plant_canopy
[1484]1393 
1394       END SELECT
[2696]1395!
[4258]1396!--    Check that at least one grid point has an LAD /= 0, else this may
1397!--    cause errors in the radiation model.
1398       lad_max = MAXVAL( lad_s )
1399#if defined( __parallel )
1400       CALL MPI_ALLREDUCE( MPI_IN_PLACE, lad_max, 1, MPI_REAL, MPI_MAX,        &
1401                           comm2d, ierr)
1402#endif
1403       IF ( lad_max <= 0.0_wp )  THEN
1404          message_string = 'Plant-canopy model is switched-on but no ' //      &
1405                           'plant canopy is present in the model domain.'
1406          CALL message( 'pcm_init', 'PA0685', 1, 2, 0, 6, 0 )
1407       ENDIF
1408   
1409!
[2696]1410!--    Initialize 2D index array indicating canopy top index.
1411       ALLOCATE( pch_index_ji(nysg:nyng,nxlg:nxrg) )
1412       pch_index_ji = 0
[4187]1413       
[4361]1414       DO  i = nxlg, nxrg
1415          DO  j = nysg, nyng
[2696]1416             DO  k = 0, pch_index
[4770]1417                IF ( lad_s(k,j,i) /= 0.0_wp  .OR.  bad_s(k,j,i) /= 0.0_wp )  pch_index_ji(j,i) = k
[2696]1418             ENDDO
[1484]1419!
[2696]1420!--          Check whether topography and local vegetation on top exceed
1421!--          height of the model domain.
[4356]1422             IF ( topo_top_ind(j,i,0) + pch_index_ji(j,i) >= nzt + 1 )  THEN
[2696]1423                message_string =  'Local vegetation height on top of ' //      &
1424                                  'topography exceeds height of model domain.'
[4356]1425                CALL message( 'pcm_init', 'PA0674', 2, 2, myid, 6, 0 )
[2696]1426             ENDIF
1427
1428          ENDDO
1429       ENDDO
[3497]1430!
[3449]1431!--    Calculate global pch_index value (index of top of plant canopy from ground)
[3497]1432       pch_index = MAXVAL( pch_index_ji )
1433!
[3449]1434!--    Exchange pch_index from all processors
1435#if defined( __parallel )
[3497]1436       CALL MPI_ALLREDUCE( MPI_IN_PLACE, pch_index, 1, MPI_INTEGER,            &
1437                           MPI_MAX, comm2d, ierr)
[3449]1438#endif
[4360]1439!
[4341]1440!--    Allocation of arrays pcm_heating_rate, pcm_transpiration_rate and pcm_latent_rate
1441       ALLOCATE( pcm_heating_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
1442       pcm_heating_rate = 0.0_wp
[4221]1443       
[3449]1444       IF ( humidity )  THEN
[4341]1445          ALLOCATE( pcm_transpiration_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
1446          pcm_transpiration_rate = 0.0_wp
1447          ALLOCATE( pcm_latent_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
1448          pcm_latent_rate = 0.0_wp
[3449]1449       ENDIF
[2696]1450!
[2011]1451!--    Initialization of the canopy heat source distribution due to heating
1452!--    of the canopy layers by incoming solar radiation, in case that a non-zero
1453!--    value is set for the canopy top heat flux (cthf), which equals the
1454!--    available net radiation at canopy top.
1455!--    The heat source distribution is calculated by a decaying exponential
1456!--    function of the downward cumulative leaf area index (cum_lai_hf),
1457!--    assuming that the foliage inside the plant canopy is heated by solar
1458!--    radiation penetrating the canopy layers according to the distribution
1459!--    of net radiation as suggested by Brown & Covey (1966; Agric. Meteorol. 3,
1460!--    73–96). This approach has been applied e.g. by Shaw & Schumann (1992;
[2213]1461!--    Bound.-Layer Meteorol. 61, 47–64).
[4341]1462!--    When using the radiation_interactions, canopy heating (pcm_heating_rate)
1463!--    and plant canopy transpiration (pcm_transpiration_rate, pcm_latent_rate)
[3449]1464!--    are calculated in the RTM after the calculation of radiation.
1465       IF ( cthf /= 0.0_wp )  THEN
[2213]1466
[3449]1467          ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1484]1468!
[2011]1469!--       Piecewise calculation of the cumulative leaf area index by vertical
[1484]1470!--       integration of the leaf area density
1471          cum_lai_hf(:,:,:) = 0.0_wp
1472          DO  i = nxlg, nxrg
1473             DO  j = nysg, nyng
[2696]1474                DO  k = pch_index_ji(j,i)-1, 0, -1
1475                   IF ( k == pch_index_ji(j,i)-1 )  THEN
[1484]1476                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
1477                         ( 0.5_wp * lad_s(k+1,j,i) *                           &
1478                           ( zw(k+1) - zu(k+1) ) )  +                          &
1479                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
1480                                                 lad_s(k,j,i) ) +              &
1481                                      lad_s(k+1,j,i) ) *                       &
1482                           ( zu(k+1) - zw(k) ) ) 
1483                   ELSE
1484                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
1485                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) +              &
1486                                                 lad_s(k+1,j,i) ) +            &
1487                                      lad_s(k+1,j,i) ) *                       &
1488                           ( zw(k+1) - zu(k+1) ) )  +                          &
1489                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
1490                                                 lad_s(k,j,i) ) +              &
1491                                      lad_s(k+1,j,i) ) *                       &
1492                           ( zu(k+1) - zw(k) ) )
1493                   ENDIF
1494                ENDDO
1495             ENDDO
1496          ENDDO
1497
[4278]1498!
[2232]1499!--       In areas with canopy the surface value of the canopy heat
1500!--       flux distribution overrides the surface heat flux (shf)
1501!--       Start with default surface type
1502          DO  m = 1, surf_def_h(0)%ns
[4278]1503             i = surf_def_h(0)%i(m)
1504             j = surf_def_h(0)%j(m)
[2232]1505             IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                &
1506                surf_def_h(0)%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
1507          ENDDO
[1484]1508!
[2232]1509!--       Natural surfaces
[4671]1510          DO  m = 1, surf_lsm_h(0)%ns
1511             i = surf_lsm_h(0)%i(m)
1512             j = surf_lsm_h(0)%j(m)
[2232]1513             IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                &
[4671]1514                surf_lsm_h(0)%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
[2232]1515          ENDDO
1516!
1517!--       Urban surfaces
[4671]1518          DO  m = 1, surf_usm_h(0)%ns
1519             i = surf_usm_h(0)%i(m)
1520             j = surf_usm_h(0)%j(m)
[2232]1521             IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                &
[4671]1522                surf_usm_h(0)%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
[2232]1523          ENDDO
1524!
1525!
[2011]1526!--       Calculation of the heating rate (K/s) within the different layers of
[2232]1527!--       the plant canopy. Calculation is only necessary in areas covered with
1528!--       canopy.
1529!--       Within the different canopy layers the plant-canopy heating
[4341]1530!--       rate (pcm_heating_rate) is calculated as the vertical
[2232]1531!--       divergence of the canopy heat fluxes at the top and bottom
1532!--       of the respective layer
[1484]1533          DO  i = nxlg, nxrg
1534             DO  j = nysg, nyng
[2696]1535                DO  k = 1, pch_index_ji(j,i)
[2232]1536                   IF ( cum_lai_hf(0,j,i) /= 0.0_wp )  THEN
[4341]1537                      pcm_heating_rate(k,j,i) = cthf *                          &
[3022]1538                                ( exp(-ext_coef*cum_lai_hf(k,j,i)) -           &
[2232]1539                                  exp(-ext_coef*cum_lai_hf(k-1,j,i) ) ) / dzw(k)
1540                   ENDIF
1541                ENDDO
[1721]1542             ENDDO
1543          ENDDO
[1484]1544
1545       ENDIF
1546
[3885]1547       IF ( debug_output )  CALL debug_message( 'pcm_init', 'end' )
[1484]1548
[1826]1549    END SUBROUTINE pcm_init
[1484]1550
1551
[2007]1552!------------------------------------------------------------------------------!
1553! Description:
1554! ------------
[2932]1555!> Parin for &plant_canopy_parameters for plant canopy model
[2007]1556!------------------------------------------------------------------------------!
1557    SUBROUTINE pcm_parin
[1484]1558
[2007]1559       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
1560       
[2932]1561       NAMELIST /plant_canopy_parameters/                                      &
1562                                  alpha_lad, beta_lad, canopy_drag_coeff,      &
1563                                  canopy_mode, cthf,                           &
[2977]1564                                  lad_surface, lad_type_coef,                  & 
[2932]1565                                  lad_vertical_gradient,                       &
1566                                  lad_vertical_gradient_level,                 &
1567                                  lai_beta,                                    &
1568                                  leaf_scalar_exch_coeff,                      &
[3449]1569                                  leaf_surface_conc, pch_index,                &
1570                                  plant_canopy_transpiration
[2932]1571
[2007]1572       NAMELIST /canopy_par/      alpha_lad, beta_lad, canopy_drag_coeff,      &
1573                                  canopy_mode, cthf,                           &
[2977]1574                                  lad_surface, lad_type_coef,                  & 
[2007]1575                                  lad_vertical_gradient,                       &
1576                                  lad_vertical_gradient_level,                 &
1577                                  lai_beta,                                    &
1578                                  leaf_scalar_exch_coeff,                      &
[3449]1579                                  leaf_surface_conc, pch_index,                &
1580                                  plant_canopy_transpiration
[3246]1581
[2007]1582       line = ' '
[3246]1583
[2007]1584!
[4258]1585!--    Try to find plant-canopy model package
[2007]1586       REWIND ( 11 )
1587       line = ' '
[3248]1588       DO WHILE ( INDEX( line, '&plant_canopy_parameters' ) == 0 )
[3246]1589          READ ( 11, '(A)', END=12 )  line
[2007]1590       ENDDO
1591       BACKSPACE ( 11 )
1592
1593!
1594!--    Read user-defined namelist
[3246]1595       READ ( 11, plant_canopy_parameters, ERR = 10 )
[2932]1596
1597!
[4258]1598!--    Set flag that indicates that the plant-canopy model is switched on
[2932]1599       plant_canopy = .TRUE.
[3246]1600
1601       GOTO 14
1602
1603 10    BACKSPACE( 11 )
[3248]1604       READ( 11 , '(A)') line
1605       CALL parin_fail_message( 'plant_canopy_parameters', line )
[2932]1606!
1607!--    Try to find old namelist
[3246]1608 12    REWIND ( 11 )
[2932]1609       line = ' '
[3248]1610       DO WHILE ( INDEX( line, '&canopy_par' ) == 0 )
[3246]1611          READ ( 11, '(A)', END=14 )  line
[2932]1612       ENDDO
1613       BACKSPACE ( 11 )
1614
1615!
1616!--    Read user-defined namelist
[3246]1617       READ ( 11, canopy_par, ERR = 13, END = 14 )
[2007]1618
[2932]1619       message_string = 'namelist canopy_par is deprecated and will be ' // &
[3046]1620                     'removed in near future. Please use namelist ' //      &
[2932]1621                     'plant_canopy_parameters instead'
1622       CALL message( 'pcm_parin', 'PA0487', 0, 1, 0, 6, 0 )
[3246]1623
[2007]1624!
[4258]1625!--    Set flag that indicates that the plant-canopy model is switched on
[2007]1626       plant_canopy = .TRUE.
1627
[3246]1628       GOTO 14
[2007]1629
[3246]1630 13    BACKSPACE( 11 )
[3248]1631       READ( 11 , '(A)') line
1632       CALL parin_fail_message( 'canopy_par', line )
[3246]1633
1634 14    CONTINUE
1635
[2007]1636    END SUBROUTINE pcm_parin
1637
1638
[1484]1639!------------------------------------------------------------------------------!
1640! Description:
1641! ------------
[2007]1642!
1643!> Loads 3D plant canopy data from file. File format is as follows:
1644!>
1645!> num_levels
[2977]1646!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
1647!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
1648!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
[2007]1649!> ...
1650!>
1651!> i.e. first line determines number of levels and further lines represent plant
1652!> canopy data, one line per column and variable. In each data line,
1653!> dtype represents variable to be set:
1654!>
1655!> dtype=1: leaf area density (lad_s)
[2213]1656!> dtype=2....n: some additional plant canopy input data quantity
[2007]1657!>
1658!> Zeros are added automatically above num_levels until top of domain.  Any
1659!> non-specified (x,y) columns have zero values as default.
1660!------------------------------------------------------------------------------!
1661    SUBROUTINE pcm_read_plant_canopy_3d
1662
[4457]1663       USE exchange_horiz_mod,                                                    &
1664           ONLY:  exchange_horiz
1665
[2213]1666       INTEGER(iwp)                        ::  dtype     !< type of input data (1=lad)
[2977]1667       INTEGER(iwp)                        ::  pctype    !< type of plant canopy (deciduous,non-deciduous,...)
[2213]1668       INTEGER(iwp)                        ::  i, j      !< running index
1669       INTEGER(iwp)                        ::  nzp       !< number of vertical layers of plant canopy
[3337]1670       INTEGER(iwp)                        ::  nzpltop   !<
1671       INTEGER(iwp)                        ::  nzpl      !<
1672       INTEGER(iwp)                        ::  kk        !<
[2213]1673       
1674       REAL(wp), DIMENSION(:), ALLOCATABLE ::  col   !< vertical column of input data
[2007]1675
[2213]1676!
1677!--    Initialize lad_s array
1678       lad_s = 0.0_wp
1679       
1680!
1681!--    Open and read plant canopy input data
[2977]1682       OPEN(152, FILE='PLANT_CANOPY_DATA_3D' // TRIM( coupling_char ),         &
1683                 ACCESS='SEQUENTIAL', ACTION='READ', STATUS='OLD',             &
1684                 FORM='FORMATTED', ERR=515)
1685       READ(152, *, ERR=516, END=517)  nzp   !< read first line = number of vertical layers
[3337]1686       nzpltop = MIN(nzt+1, nzb+nzp-1)
1687       nzpl = nzpltop - nzb + 1    !< no. of layers to assign
[2977]1688       ALLOCATE( col(0:nzp-1) )
[2007]1689
[2213]1690       DO
[2977]1691          READ(152, *, ERR=516, END=517) dtype, i, j, pctype, col(:)
1692          IF ( i < nxlg  .OR.  i > nxrg  .OR.  j < nysg  .OR.  j > nyng )  CYCLE
1693         
1694          SELECT CASE (dtype)
1695             CASE( 1 )   !< leaf area density
[2213]1696!
[2977]1697!--             This is just the pure canopy layer assumed to be grounded to
1698!--             a flat domain surface. At locations where plant canopy sits
1699!--             on top of any kind of topography, the vertical plant column
1700!--             must be "lifted", which is done in SUBROUTINE pcm_tendency.           
1701                IF ( pctype < 0  .OR.  pctype > 10 )  THEN   !< incorrect plant canopy type
1702                   WRITE( message_string, * ) 'Incorrect type of plant canopy. '   //  &
1703                                              'Allowed values 0 <= pctype <= 10, ' //  &
1704                                              'but pctype is ', pctype
1705                   CALL message( 'pcm_read_plant_canopy_3d', 'PA0349', 1, 2, 0, 6, 0 )
1706                ENDIF
[4168]1707                kk = topo_top_ind(j,i,0)
[3337]1708                lad_s(nzb:nzpltop-kk, j, i) = col(kk:nzpl-1)*lad_type_coef(pctype)
[2977]1709             CASE DEFAULT
1710                WRITE(message_string, '(a,i2,a)')   &
1711                     'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"'
1712                CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 )
1713          END SELECT
[2213]1714       ENDDO
[2007]1715
[2213]1716515    message_string = 'error opening file PLANT_CANOPY_DATA_3D'
1717       CALL message( 'pcm_read_plant_canopy_3d', 'PA0531', 1, 2, 0, 6, 0 )
[2007]1718
[2213]1719516    message_string = 'error reading file PLANT_CANOPY_DATA_3D'
1720       CALL message( 'pcm_read_plant_canopy_3d', 'PA0532', 1, 2, 0, 6, 0 )
1721
1722517    CLOSE(152)
[2977]1723       DEALLOCATE( col )
[2213]1724       
1725       CALL exchange_horiz( lad_s, nbgp )
[2007]1726       
1727    END SUBROUTINE pcm_read_plant_canopy_3d
[4360]1728
[2007]1729!------------------------------------------------------------------------------!
1730! Description:
1731! ------------
[4514]1732!> Read module-specific global restart data (Fortran binary format).
1733!------------------------------------------------------------------------------!
1734    SUBROUTINE pcm_rrd_global_ftn( found )
1735
1736       LOGICAL, INTENT(OUT)  ::  found
1737
1738       found = .TRUE.
1739
1740       SELECT CASE ( restart_string(1:length) )
1741
1742          CASE ( 'pch_index' )
1743             READ ( 13 )  pch_index
1744
1745          CASE DEFAULT
1746
1747             found = .FALSE.
1748
1749       END SELECT
1750
1751    END SUBROUTINE pcm_rrd_global_ftn
1752
1753!------------------------------------------------------------------------------!
1754! Description:
1755! ------------
1756!> Read module-specific global restart data (MPI-IO).
1757!------------------------------------------------------------------------------!
1758    SUBROUTINE pcm_rrd_global_mpi
1759
1760       CALL rrd_mpi_io( 'pch_index', pch_index )
1761
1762    END SUBROUTINE pcm_rrd_global_mpi
1763
[4517]1764
[4514]1765!------------------------------------------------------------------------------!
1766! Description:
1767! ------------
[4517]1768!> Read module-specific local restart data arrays (Fortran binary format).
[4360]1769!------------------------------------------------------------------------------!
[4517]1770    SUBROUTINE pcm_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,          &
1771                                  nxr_on_file, nynf, nync, nyn_on_file, nysf,      &
1772                                  nysc, nys_on_file, found )
[4360]1773
1774       INTEGER(iwp) ::  k               !<
1775       INTEGER(iwp) ::  nxlc            !<
1776       INTEGER(iwp) ::  nxlf            !<
1777       INTEGER(iwp) ::  nxl_on_file     !<
1778       INTEGER(iwp) ::  nxrc            !<
1779       INTEGER(iwp) ::  nxrf            !<
1780       INTEGER(iwp) ::  nxr_on_file     !<
1781       INTEGER(iwp) ::  nync            !<
1782       INTEGER(iwp) ::  nynf            !<
1783       INTEGER(iwp) ::  nyn_on_file     !<
1784       INTEGER(iwp) ::  nysc            !<
1785       INTEGER(iwp) ::  nysf            !<
1786       INTEGER(iwp) ::  nys_on_file     !<
1787
1788       LOGICAL, INTENT(OUT)  :: found
1789
1790       REAL(wp), DIMENSION(0:pch_index,                                        &
1791                           nys_on_file-nbgp:nyn_on_file+nbgp,                  &
1792                           nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d2  !< temporary 3D array for entire vertical
1793                                                                          !< extension of canopy layer
1794       found = .TRUE.
1795
1796       SELECT CASE ( restart_string(1:length) )
1797
1798          CASE ( 'pcm_heatrate_av' )
1799             IF ( .NOT. ALLOCATED( pcm_heatrate_av ) )  THEN
1800                ALLOCATE( pcm_heatrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
1801                pcm_heatrate_av = 0.0_wp
1802             ENDIF 
1803             IF ( k == 1 )  READ ( 13 )  tmp_3d2
1804             pcm_heatrate_av(0:pch_index,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1805                        tmp_3d2(0:pch_index,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1806
1807          CASE ( 'pcm_latentrate_av' )
1808             IF ( .NOT. ALLOCATED( pcm_latentrate_av ) )  THEN
1809                ALLOCATE( pcm_latentrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
1810                pcm_latentrate_av = 0.0_wp
1811             ENDIF 
1812             IF ( k == 1 )  READ ( 13 )  tmp_3d2
1813             pcm_latentrate_av(0:pch_index,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1814                        tmp_3d2(0:pch_index,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1815
1816          CASE ( 'pcm_transpirationrate_av' )
1817             IF ( .NOT. ALLOCATED( pcm_transpirationrate_av ) )  THEN
1818                ALLOCATE( pcm_transpirationrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
1819                pcm_transpirationrate_av = 0.0_wp
1820             ENDIF 
1821             IF ( k == 1 )  READ ( 13 )  tmp_3d2
1822             pcm_transpirationrate_av(0:pch_index,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1823                        tmp_3d2(0:pch_index,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1824
1825          CASE DEFAULT
1826
1827             found = .FALSE.
1828
1829       END SELECT
1830
[4517]1831    END SUBROUTINE pcm_rrd_local_ftn
[4360]1832
[4517]1833
[4360]1834!------------------------------------------------------------------------------!
1835! Description:
1836! ------------
[4517]1837!> Read module-specific local restart data arrays (MPI-IO).
1838!------------------------------------------------------------------------------!
1839    SUBROUTINE pcm_rrd_local_mpi
1840
1841       IMPLICIT NONE
1842
1843       LOGICAL ::  array_found  !<
1844
[4525]1845       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  tmp_3d  !< temporary array to store pcm data with
1846                                                           !< non-standard vertical index bounds
[4517]1847
[4525]1848!
1849!--    Plant canopy arrays have non standard reduced vertical index bounds. They are stored with
1850!--    full vertical bounds (bzb:nzt+1) in the restart file and must be re-stored after reading.
[4517]1851       CALL rd_mpi_io_check_array( 'pcm_heatrate_av' , found = array_found )
1852       IF ( array_found )  THEN
[4525]1853          IF ( .NOT. ALLOCATED( pcm_heatrate_av ) )  THEN
1854             ALLOCATE( pcm_heatrate_av(nzb:pch_index,nysg:nyng,nxlg:nxrg) )
1855          ENDIF
1856          ALLOCATE( tmp_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1857          CALL rrd_mpi_io( 'pcm_heatrate_av', tmp_3d )
1858          pcm_heatrate_av = tmp_3d(nzb:pch_index,:,:)
1859          DEALLOCATE( tmp_3d )
[4517]1860       ENDIF
1861
1862       CALL rd_mpi_io_check_array( 'pcm_latentrate_av' , found = array_found )
1863       IF ( array_found )  THEN
[4525]1864          IF ( .NOT. ALLOCATED( pcm_latentrate_av ) )  THEN
1865             ALLOCATE( pcm_latentrate_av(nzb:pch_index,nysg:nyng,nxlg:nxrg) )
1866          ENDIF
1867          ALLOCATE( tmp_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1868          CALL rrd_mpi_io( 'pcm_latentrate_av', tmp_3d )
1869          pcm_latentrate_av = tmp_3d(nzb:pch_index,:,:)
1870          DEALLOCATE( tmp_3d )
[4517]1871       ENDIF
1872
1873       CALL rd_mpi_io_check_array( 'pcm_transpirationrate_av' , found = array_found )
1874       IF ( array_found )  THEN
[4525]1875          IF ( .NOT. ALLOCATED( pcm_transpirationrate_av ) )  THEN
1876             ALLOCATE( pcm_transpirationrate_av(nzb:pch_index,nysg:nyng,nxlg:nxrg) )
1877          ENDIF
1878          ALLOCATE( tmp_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1879          CALL rrd_mpi_io( 'pcm_transpirationrate_av', tmp_3d )
1880          pcm_transpirationrate_av = tmp_3d(nzb:pch_index,:,:)
1881          DEALLOCATE( tmp_3d )
[4517]1882       ENDIF
1883
1884    END SUBROUTINE pcm_rrd_local_mpi
1885
1886
1887!------------------------------------------------------------------------------!
1888! Description:
1889! ------------
[1682]1890!> Calculation of the tendency terms, accounting for the effect of the plant
1891!> canopy on momentum and scalar quantities.
1892!>
1893!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
[1826]1894!> (defined on scalar grid), as initialized in subroutine pcm_init.
[1682]1895!> The lad on the w-grid is vertically interpolated from the surrounding
1896!> lad_s. The upper boundary of the canopy is defined on the w-grid at
1897!> k = pch_index. Here, the lad is zero.
1898!>
1899!> The canopy drag must be limited (previously accounted for by calculation of
1900!> a limiting canopy timestep for the determination of the maximum LES timestep
1901!> in subroutine timestep), since it is physically impossible that the canopy
1902!> drag alone can locally change the sign of a velocity component. This
1903!> limitation is realized by calculating preliminary tendencies and velocities.
1904!> It is subsequently checked if the preliminary new velocity has a different
1905!> sign than the current velocity. If so, the tendency is limited in a way that
1906!> the velocity can at maximum be reduced to zero by the canopy drag.
1907!>
1908!>
1909!> Call for all grid points
[1484]1910!------------------------------------------------------------------------------!
[1826]1911    SUBROUTINE pcm_tendency( component )
[138]1912
[1682]1913       INTEGER(iwp) ::  component !< prognostic variable (u,v,w,pt,q,e)
1914       INTEGER(iwp) ::  i         !< running index
1915       INTEGER(iwp) ::  j         !< running index
1916       INTEGER(iwp) ::  k         !< running index
[1721]1917       INTEGER(iwp) ::  kk        !< running index for flat lad arrays
[1484]1918
[4335]1919       LOGICAL ::  building_edge_e !< control flag indicating an eastward-facing building edge
1920       LOGICAL ::  building_edge_n !< control flag indicating a north-facing building edge
1921       LOGICAL ::  building_edge_s !< control flag indicating a south-facing building edge
1922       LOGICAL ::  building_edge_w !< control flag indicating a westward-facing building edge
1923
[4770]1924       REAL(wp) ::  bad_local !< local bad value
[1682]1925       REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
1926       REAL(wp) ::  lad_local !< local lad value
1927       REAL(wp) ::  pre_tend  !< preliminary tendency
1928       REAL(wp) ::  pre_u     !< preliminary u-value
1929       REAL(wp) ::  pre_v     !< preliminary v-value
1930       REAL(wp) ::  pre_w     !< preliminary w-value
[1484]1931
1932
1933       ddt_3d = 1.0_wp / dt_3d
[138]1934 
1935!
[1484]1936!--    Compute drag for the three velocity components and the SGS-TKE:
[138]1937       SELECT CASE ( component )
1938
1939!
1940!--       u-component
1941          CASE ( 1 )
1942             DO  i = nxlu, nxr
1943                DO  j = nys, nyn
[2232]1944!
[4335]1945!--                Set control flags indicating east- and westward-orientated
1946!--                building edges. Note, building_egde_w is set from the perspective
1947!--                of the potential rooftop grid point, while building_edge_e is
1948!--                is set from the perspective of the non-building grid point.
[4346]1949                   building_edge_w = ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )&
1950                        .AND.  .NOT. ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )
1951                   building_edge_e = ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )&
1952                        .AND.  .NOT. ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )
[4335]1953!
[2232]1954!--                Determine topography-top index on u-grid
[4341]1955                   DO  k = topo_top_ind(j,i,1)+1, topo_top_ind(j,i,1) + pch_index_ji(j,i)
[1484]1956
[4341]1957                      kk = k - topo_top_ind(j,i,1)   !- lad arrays are defined flat
[1484]1958!
1959!--                   In order to create sharp boundaries of the plant canopy,
[4335]1960!--                   the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
1961!--                   rather than being interpolated from the surrounding lad_s,
1962!--                   because this would yield smaller lad at the canopy boundaries
1963!--                   than inside of the canopy.
[1484]1964!--                   For the same reason, the lad at the rightmost(i+1)canopy
[4335]1965!--                   boundary on the u-grid equals lad_s(k,j,i), which is considered
1966!--                   in the next if-statement. Note, at left-sided building edges
1967!--                   this is not applied, here the LAD is equals the LAD at grid
1968!--                   point (k,j,i), in order to avoid that LAD is mistakenly mapped
[4770]1969!--                   on top of a roof where (usually) is no LAD is defined. The same is
1970!--                   also valid for bad_s.
[1721]1971                      lad_local = lad_s(kk,j,i)
[4770]1972                      IF ( lad_local == 0.0_wp  .AND.  lad_s(kk,j,i-1) > 0.0_wp                    &
[4335]1973                           .AND.  .NOT. building_edge_w )  lad_local = lad_s(kk,j,i-1)
[4770]1974
1975                      bad_local = bad_s(kk,j,i)
1976                      IF ( bad_local == 0.0_wp  .AND.  bad_s(kk,j,i-1) > 0.0_wp                    &
1977                           .AND.  .NOT. building_edge_w )  bad_local = bad_s(kk,j,i-1)
[4335]1978!
1979!--                   In order to avoid that LAD is mistakenly considered at right-
1980!--                   sided building edges (here the topography-top index for the
1981!--                   u-component at index j,i is still on the building while the
1982!--                   topography top for the scalar isn't), LAD is taken from grid
[4770]1983!--                   point (j,i-1). The same is also valid for bad_s.
1984                      IF ( lad_local > 0.0_wp  .AND.  lad_s(kk,j,i-1) == 0.0_wp                    &
[4335]1985                           .AND.  building_edge_e )  lad_local = lad_s(kk,j,i-1)
[4770]1986                      IF ( bad_local > 0.0_wp  .AND.  bad_s(kk,j,i-1) == 0.0_wp                    &
1987                           .AND.  building_edge_e )  bad_local = bad_s(kk,j,i-1)
[1484]1988
1989                      pre_tend = 0.0_wp
1990                      pre_u = 0.0_wp
1991!
1992!--                   Calculate preliminary value (pre_tend) of the tendency
[4342]1993                      pre_tend = - canopy_drag_coeff *                         &
[4770]1994                                   ( lad_local + bad_local ) *                 &
[1484]1995                                   SQRT( u(k,j,i)**2 +                         &
1996                                         ( 0.25_wp * ( v(k,j,i-1) +            &
1997                                                       v(k,j,i)   +            &
1998                                                       v(k,j+1,i) +            &
1999                                                       v(k,j+1,i-1) )          &
2000                                         )**2 +                                &
2001                                         ( 0.25_wp * ( w(k-1,j,i-1) +          &
2002                                                       w(k-1,j,i)   +          &
2003                                                       w(k,j,i-1)   +          &
2004                                                       w(k,j,i) )              &
2005                                         )**2                                  &
2006                                       ) *                                     &
2007                                   u(k,j,i)
2008
2009!
2010!--                   Calculate preliminary new velocity, based on pre_tend
2011                      pre_u = u(k,j,i) + dt_3d * pre_tend
2012!
2013!--                   Compare sign of old velocity and new preliminary velocity,
2014!--                   and in case the signs are different, limit the tendency
2015                      IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
2016                         pre_tend = - u(k,j,i) * ddt_3d
2017                      ENDIF
2018!
2019!--                   Calculate final tendency
2020                      tend(k,j,i) = tend(k,j,i) + pre_tend
2021
[138]2022                   ENDDO
2023                ENDDO
2024             ENDDO
2025
2026!
2027!--       v-component
2028          CASE ( 2 )
2029             DO  i = nxl, nxr
2030                DO  j = nysv, nyn
[2232]2031!
[4335]2032!--                Set control flags indicating north- and southward-orientated
2033!--                building edges. Note, building_egde_s is set from the perspective
2034!--                of the potential rooftop grid point, while building_edge_n is
2035!--                is set from the perspective of the non-building grid point.
[4346]2036                   building_edge_s = ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )&
2037                        .AND.  .NOT. ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )
2038                   building_edge_n = ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )&
2039                        .AND.  .NOT. ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )
[4335]2040!
[2232]2041!--                Determine topography-top index on v-grid
[4341]2042                   DO  k = topo_top_ind(j,i,2)+1, topo_top_ind(j,i,2) + pch_index_ji(j,i)
[2317]2043
[4341]2044                      kk = k - topo_top_ind(j,i,2)   !- lad arrays are defined flat
[1484]2045!
2046!--                   In order to create sharp boundaries of the plant canopy,
[4335]2047!--                   the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
2048!--                   rather than being interpolated from the surrounding lad_s,
2049!--                   because this would yield smaller lad at the canopy boundaries
2050!--                   than inside of the canopy.
2051!--                   For the same reason, the lad at the northmost(j+1)canopy
2052!--                   boundary on the v-grid equals lad_s(k,j,i), which is considered
2053!--                   in the next if-statement. Note, at left-sided building edges
2054!--                   this is not applied, here the LAD is equals the LAD at grid
2055!--                   point (k,j,i), in order to avoid that LAD is mistakenly mapped
2056!--                   on top of a roof where (usually) is no LAD is defined.
[4770]2057!--                   The same is also valid for bad_s.
[1721]2058                      lad_local = lad_s(kk,j,i)
[4770]2059                      IF ( lad_local == 0.0_wp  .AND.  lad_s(kk,j-1,i) > 0.0_wp                    &
[4335]2060                       .AND.  .NOT. building_edge_s )  lad_local = lad_s(kk,j-1,i)
[4770]2061
2062                      bad_local = bad_s(kk,j,i)
2063                      IF ( bad_local == 0.0_wp  .AND.  bad_s(kk,j-1,i) > 0.0_wp                    &
2064                       .AND.  .NOT. building_edge_s )  bad_local = bad_s(kk,j-1,i)
[4335]2065!
2066!--                   In order to avoid that LAD is mistakenly considered at right-
2067!--                   sided building edges (here the topography-top index for the
2068!--                   u-component at index j,i is still on the building while the
2069!--                   topography top for the scalar isn't), LAD is taken from grid
[4770]2070!--                   point (j,i-1). The same is also valid for bad_s.
2071                      IF ( lad_local > 0.0_wp  .AND.  lad_s(kk,j-1,i) == 0.0_wp                    &
[4335]2072                       .AND.  building_edge_n )  lad_local = lad_s(kk,j-1,i)
[1484]2073
[4770]2074                      IF ( bad_local > 0.0_wp  .AND.  bad_s(kk,j-1,i) == 0.0_wp                    &
2075                       .AND.  building_edge_n )  bad_local = bad_s(kk,j-1,i)
2076
[1484]2077                      pre_tend = 0.0_wp
2078                      pre_v = 0.0_wp
2079!
2080!--                   Calculate preliminary value (pre_tend) of the tendency
[4342]2081                      pre_tend = - canopy_drag_coeff *                         &
[4770]2082                                   ( lad_local + bad_local ) *                 &
[1484]2083                                   SQRT( ( 0.25_wp * ( u(k,j-1,i)   +          &
2084                                                       u(k,j-1,i+1) +          &
2085                                                       u(k,j,i)     +          &
2086                                                       u(k,j,i+1) )            &
2087                                         )**2 +                                &
2088                                         v(k,j,i)**2 +                         &
2089                                         ( 0.25_wp * ( w(k-1,j-1,i) +          &
2090                                                       w(k-1,j,i)   +          &
2091                                                       w(k,j-1,i)   +          &
2092                                                       w(k,j,i) )              &
2093                                         )**2                                  &
2094                                       ) *                                     &
2095                                   v(k,j,i)
2096
2097!
2098!--                   Calculate preliminary new velocity, based on pre_tend
2099                      pre_v = v(k,j,i) + dt_3d * pre_tend
2100!
2101!--                   Compare sign of old velocity and new preliminary velocity,
2102!--                   and in case the signs are different, limit the tendency
2103                      IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
2104                         pre_tend = - v(k,j,i) * ddt_3d
2105                      ELSE
2106                         pre_tend = pre_tend
2107                      ENDIF
2108!
2109!--                   Calculate final tendency
2110                      tend(k,j,i) = tend(k,j,i) + pre_tend
2111
[138]2112                   ENDDO
2113                ENDDO
2114             ENDDO
2115
2116!
2117!--       w-component
2118          CASE ( 3 )
2119             DO  i = nxl, nxr
2120                DO  j = nys, nyn
[2232]2121!
2122!--                Determine topography-top index on w-grid
[4341]2123                   DO  k = topo_top_ind(j,i,3)+1, topo_top_ind(j,i,3) + pch_index_ji(j,i) - 1
[2317]2124
[4341]2125                      kk = k - topo_top_ind(j,i,3)   !- lad arrays are defined flat
[1484]2126
2127                      pre_tend = 0.0_wp
2128                      pre_w = 0.0_wp
2129!
2130!--                   Calculate preliminary value (pre_tend) of the tendency
[4342]2131                      pre_tend = - canopy_drag_coeff *                         &
[4770]2132                                   ( 0.5_wp *                                  &
2133                                      ( lad_s(kk+1,j,i) + lad_s(kk,j,i) ) +    &
2134                                     0.5_wp *                                  &
2135                                      ( bad_s(kk+1,j,i) + bad_s(kk,j,i) ) ) *  &
[1484]2136                                   SQRT( ( 0.25_wp * ( u(k,j,i)   +            &
2137                                                       u(k,j,i+1) +            &
2138                                                       u(k+1,j,i) +            &
2139                                                       u(k+1,j,i+1) )          &
2140                                         )**2 +                                &
2141                                         ( 0.25_wp * ( v(k,j,i)   +            &
2142                                                       v(k,j+1,i) +            &
2143                                                       v(k+1,j,i) +            &
2144                                                       v(k+1,j+1,i) )          &
2145                                         )**2 +                                &
2146                                         w(k,j,i)**2                           &
2147                                       ) *                                     &
2148                                   w(k,j,i)
2149!
2150!--                   Calculate preliminary new velocity, based on pre_tend
2151                      pre_w = w(k,j,i) + dt_3d * pre_tend
2152!
2153!--                   Compare sign of old velocity and new preliminary velocity,
2154!--                   and in case the signs are different, limit the tendency
2155                      IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
2156                         pre_tend = - w(k,j,i) * ddt_3d
2157                      ELSE
2158                         pre_tend = pre_tend
2159                      ENDIF
2160!
2161!--                   Calculate final tendency
2162                      tend(k,j,i) = tend(k,j,i) + pre_tend
2163
[138]2164                   ENDDO
2165                ENDDO
2166             ENDDO
2167
2168!
[153]2169!--       potential temperature
[138]2170          CASE ( 4 )
[3449]2171             IF ( humidity ) THEN
2172                DO  i = nxl, nxr
2173                   DO  j = nys, nyn
2174!--                   Determine topography-top index on scalar-grid
[4341]2175                      DO  k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2176                         kk = k - topo_top_ind(j,i,0)   !- lad arrays are defined flat
2177                         tend(k,j,i) = tend(k,j,i) + pcm_heating_rate(kk,j,i) - pcm_latent_rate(kk,j,i)
[3449]2178                      ENDDO
[153]2179                   ENDDO
2180                ENDDO
[3449]2181             ELSE
2182                DO  i = nxl, nxr
2183                   DO  j = nys, nyn
2184!--                   Determine topography-top index on scalar-grid
[4341]2185                      DO  k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2186                         kk = k - topo_top_ind(j,i,0)   !- lad arrays are defined flat
2187                         tend(k,j,i) = tend(k,j,i) + pcm_heating_rate(kk,j,i)
[3449]2188                      ENDDO
2189                   ENDDO
2190                ENDDO
2191             ENDIF
[153]2192
2193!
[1960]2194!--       humidity
[153]2195          CASE ( 5 )
2196             DO  i = nxl, nxr
2197                DO  j = nys, nyn
[2232]2198!
2199!--                Determine topography-top index on scalar-grid
[4341]2200                   DO  k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
[2317]2201
[4341]2202                      kk = k - topo_top_ind(j,i,0)   !- lad arrays are defined flat
[2232]2203
[3449]2204                      IF ( .NOT. plant_canopy_transpiration ) THEN
[4341]2205                         ! pcm_transpiration_rate is calculated in radiation model
[3449]2206                         ! in case of plant_canopy_transpiration = .T.
2207                         ! to include also the dependecy to the radiation
2208                         ! in the plant canopy box
[4342]2209                         pcm_transpiration_rate(kk,j,i) =  - leaf_scalar_exch_coeff &
2210                                          * lad_s(kk,j,i) *                         &
2211                                          SQRT( ( 0.5_wp * ( u(k,j,i) +             &
2212                                                             u(k,j,i+1) )           &
2213                                                )**2 +                              &
2214                                                ( 0.5_wp * ( v(k,j,i) +             &
2215                                                             v(k,j+1,i) )           &
2216                                                )**2 +                              &
2217                                                ( 0.5_wp * ( w(k-1,j,i) +           &
2218                                                             w(k,j,i) )             &
2219                                                )**2                                &
2220                                              ) *                                   &
2221                                          ( q(k,j,i) - leaf_surface_conc )
[3449]2222                      ENDIF
2223
[4341]2224                      tend(k,j,i) = tend(k,j,i) + pcm_transpiration_rate(kk,j,i)
[153]2225                   ENDDO
2226                ENDDO
2227             ENDDO
2228
2229!
2230!--       sgs-tke
2231          CASE ( 6 )
2232             DO  i = nxl, nxr
2233                DO  j = nys, nyn
[2232]2234!
2235!--                Determine topography-top index on scalar-grid
[4341]2236                   DO  k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
[2317]2237
[4341]2238                      kk = k - topo_top_ind(j,i,0)   !- lad arrays are defined flat
[1484]2239                      tend(k,j,i) = tend(k,j,i) -                              &
[4342]2240                                       2.0_wp * canopy_drag_coeff *            &
[4770]2241                                       ( lad_s(kk,j,i) + bad_s(kk,j,i) ) *     &
[1484]2242                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
2243                                                          u(k,j,i+1) )         &
2244                                             )**2 +                            &
2245                                             ( 0.5_wp * ( v(k,j,i) +           &
2246                                                          v(k,j+1,i) )         &
2247                                             )**2 +                            &
2248                                             ( 0.5_wp * ( w(k,j,i) +           &
2249                                                          w(k+1,j,i) )         &
2250                                             )**2                              &
2251                                           ) *                                 &
2252                                       e(k,j,i)
[138]2253                   ENDDO
2254                ENDDO
2255             ENDDO 
[1960]2256!
2257!--       scalar concentration
2258          CASE ( 7 )
2259             DO  i = nxl, nxr
2260                DO  j = nys, nyn
[2232]2261!
2262!--                Determine topography-top index on scalar-grid
[4341]2263                   DO  k = topo_top_ind(j,i,0)+1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
[2317]2264
[4341]2265                      kk = k - topo_top_ind(j,i,0)   !- lad arrays are defined flat
[1960]2266                      tend(k,j,i) = tend(k,j,i) -                              &
[4342]2267                                       leaf_scalar_exch_coeff *                &
[1960]2268                                       lad_s(kk,j,i) *                         &
2269                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
2270                                                          u(k,j,i+1) )         &
2271                                             )**2 +                            &
2272                                             ( 0.5_wp * ( v(k,j,i) +           &
2273                                                          v(k,j+1,i) )         &
2274                                             )**2 +                            &
2275                                             ( 0.5_wp * ( w(k-1,j,i) +         & 
2276                                                          w(k,j,i) )           &
2277                                             )**2                              &
2278                                           ) *                                 &
[4342]2279                                       ( s(k,j,i) - leaf_surface_conc )
[1960]2280                   ENDDO
2281                ENDDO
2282             ENDDO   
[1484]2283
2284
[1960]2285
[138]2286          CASE DEFAULT
2287
[257]2288             WRITE( message_string, * ) 'wrong component: ', component
[1826]2289             CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 ) 
[138]2290
2291       END SELECT
2292
[1826]2293    END SUBROUTINE pcm_tendency
[138]2294
2295
2296!------------------------------------------------------------------------------!
[1484]2297! Description:
2298! ------------
[1682]2299!> Calculation of the tendency terms, accounting for the effect of the plant
2300!> canopy on momentum and scalar quantities.
2301!>
2302!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
[1826]2303!> (defined on scalar grid), as initialized in subroutine pcm_init.
[1682]2304!> The lad on the w-grid is vertically interpolated from the surrounding
2305!> lad_s. The upper boundary of the canopy is defined on the w-grid at
2306!> k = pch_index. Here, the lad is zero.
2307!>
2308!> The canopy drag must be limited (previously accounted for by calculation of
2309!> a limiting canopy timestep for the determination of the maximum LES timestep
2310!> in subroutine timestep), since it is physically impossible that the canopy
2311!> drag alone can locally change the sign of a velocity component. This
2312!> limitation is realized by calculating preliminary tendencies and velocities.
2313!> It is subsequently checked if the preliminary new velocity has a different
2314!> sign than the current velocity. If so, the tendency is limited in a way that
2315!> the velocity can at maximum be reduced to zero by the canopy drag.
2316!>
2317!>
2318!> Call for grid point i,j
[138]2319!------------------------------------------------------------------------------!
[1826]2320    SUBROUTINE pcm_tendency_ij( i, j, component )
[138]2321
[1682]2322       INTEGER(iwp) ::  component !< prognostic variable (u,v,w,pt,q,e)
2323       INTEGER(iwp) ::  i         !< running index
2324       INTEGER(iwp) ::  j         !< running index
2325       INTEGER(iwp) ::  k         !< running index
[1721]2326       INTEGER(iwp) ::  kk        !< running index for flat lad arrays
[138]2327
[4314]2328       LOGICAL ::  building_edge_e !< control flag indicating an eastward-facing building edge
2329       LOGICAL ::  building_edge_n !< control flag indicating a north-facing building edge
2330       LOGICAL ::  building_edge_s !< control flag indicating a south-facing building edge
2331       LOGICAL ::  building_edge_w !< control flag indicating a westward-facing building edge
2332
[4770]2333       REAL(wp) ::  bad_local !< local lad value
[1682]2334       REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
2335       REAL(wp) ::  lad_local !< local lad value
2336       REAL(wp) ::  pre_tend  !< preliminary tendency
2337       REAL(wp) ::  pre_u     !< preliminary u-value
2338       REAL(wp) ::  pre_v     !< preliminary v-value
2339       REAL(wp) ::  pre_w     !< preliminary w-value
[1484]2340
2341
2342       ddt_3d = 1.0_wp / dt_3d
[138]2343!
[1484]2344!--    Compute drag for the three velocity components and the SGS-TKE
[142]2345       SELECT CASE ( component )
[138]2346
2347!
[142]2348!--       u-component
[1484]2349          CASE ( 1 )
[2232]2350!
[4314]2351!--          Set control flags indicating east- and westward-orientated
2352!--          building edges. Note, building_egde_w is set from the perspective
2353!--          of the potential rooftop grid point, while building_edge_e is
2354!--          is set from the perspective of the non-building grid point.
[4346]2355             building_edge_w = ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )  .AND. &
2356                         .NOT. ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )
2357             building_edge_e = ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )  .AND. &
2358                         .NOT. ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )
[4314]2359!
[2232]2360!--          Determine topography-top index on u-grid
[4341]2361             DO  k = topo_top_ind(j,i,1) + 1, topo_top_ind(j,i,1) + pch_index_ji(j,i)
[2317]2362
[4341]2363                kk = k - topo_top_ind(j,i,1)  !- lad arrays are defined flat
[138]2364
2365!
[1484]2366!--             In order to create sharp boundaries of the plant canopy,
2367!--             the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
2368!--             rather than being interpolated from the surrounding lad_s,
2369!--             because this would yield smaller lad at the canopy boundaries
2370!--             than inside of the canopy.
2371!--             For the same reason, the lad at the rightmost(i+1)canopy
[4314]2372!--             boundary on the u-grid equals lad_s(k,j,i), which is considered
2373!--             in the next if-statement. Note, at left-sided building edges
2374!--             this is not applied, here the LAD is equals the LAD at grid
2375!--             point (k,j,i), in order to avoid that LAD is mistakenly mapped
2376!--             on top of a roof where (usually) is no LAD is defined.
[4770]2377!--             The same is also valid for bad_s.
[1721]2378                lad_local = lad_s(kk,j,i)
[4770]2379                IF ( lad_local == 0.0_wp  .AND.  lad_s(kk,j,i-1) > 0.0_wp  .AND.                   &
[4314]2380                     .NOT. building_edge_w )  lad_local = lad_s(kk,j,i-1)
[4770]2381
2382                bad_local = bad_s(kk,j,i)
2383                IF ( bad_local == 0.0_wp  .AND.  bad_s(kk,j,i-1) > 0.0_wp  .AND.                   &
2384                     .NOT. building_edge_w )  bad_local = bad_s(kk,j,i-1)
[4314]2385!
2386!--             In order to avoid that LAD is mistakenly considered at right-
2387!--             sided building edges (here the topography-top index for the
2388!--             u-component at index j,i is still on the building while the
2389!--             topography top for the scalar isn't), LAD is taken from grid
[4770]2390!--             point (j,i-1). The same is also valid for bad_s.
2391                IF ( lad_local > 0.0_wp  .AND.  lad_s(kk,j,i-1) == 0.0_wp  .AND.                   &
[4314]2392                     building_edge_e )  lad_local = lad_s(kk,j,i-1)
[4770]2393                IF ( bad_local > 0.0_wp  .AND.  bad_s(kk,j,i-1) == 0.0_wp  .AND.                   &
2394                     building_edge_e )  bad_local = bad_s(kk,j,i-1)
[1484]2395
2396                pre_tend = 0.0_wp
2397                pre_u = 0.0_wp
2398!
2399!--             Calculate preliminary value (pre_tend) of the tendency
[4342]2400                pre_tend = - canopy_drag_coeff *                               &
[4770]2401                             ( lad_local + bad_local ) *                       &
[1484]2402                             SQRT( u(k,j,i)**2 +                               &
2403                                   ( 0.25_wp * ( v(k,j,i-1)  +                 &
2404                                                 v(k,j,i)    +                 &
2405                                                 v(k,j+1,i)  +                 &
2406                                                 v(k,j+1,i-1) )                &
2407                                   )**2 +                                      &
2408                                   ( 0.25_wp * ( w(k-1,j,i-1) +                &
2409                                                 w(k-1,j,i)   +                &
2410                                                 w(k,j,i-1)   +                &
2411                                                 w(k,j,i) )                    &
2412                                   )**2                                        &
2413                                 ) *                                           &
2414                             u(k,j,i)
2415
2416!
2417!--             Calculate preliminary new velocity, based on pre_tend
2418                pre_u = u(k,j,i) + dt_3d * pre_tend
2419!
2420!--             Compare sign of old velocity and new preliminary velocity,
2421!--             and in case the signs are different, limit the tendency
2422                IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
2423                   pre_tend = - u(k,j,i) * ddt_3d
2424                ELSE
2425                   pre_tend = pre_tend
2426                ENDIF
2427!
2428!--             Calculate final tendency
2429                tend(k,j,i) = tend(k,j,i) + pre_tend
2430             ENDDO
2431
2432
2433!
[142]2434!--       v-component
[1484]2435          CASE ( 2 )
[2232]2436!
[4314]2437!--          Set control flags indicating north- and southward-orientated
2438!--          building edges. Note, building_egde_s is set from the perspective
2439!--          of the potential rooftop grid point, while building_edge_n is
2440!--          is set from the perspective of the non-building grid point.
[4346]2441             building_edge_s = ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )  .AND. &
2442                         .NOT. ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )
2443             building_edge_n = ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )  .AND. &
2444                         .NOT. ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )
[4314]2445!
[2232]2446!--          Determine topography-top index on v-grid
[4341]2447             DO  k = topo_top_ind(j,i,2) + 1, topo_top_ind(j,i,2) + pch_index_ji(j,i)
[138]2448
[4341]2449                kk = k - topo_top_ind(j,i,2)  !- lad arrays are defined flat
[138]2450!
[1484]2451!--             In order to create sharp boundaries of the plant canopy,
2452!--             the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
2453!--             rather than being interpolated from the surrounding lad_s,
2454!--             because this would yield smaller lad at the canopy boundaries
2455!--             than inside of the canopy.
2456!--             For the same reason, the lad at the northmost(j+1)canopy
[4314]2457!--             boundary on the v-grid equals lad_s(k,j,i), which is considered
2458!--             in the next if-statement. Note, at left-sided building edges
2459!--             this is not applied, here the LAD is equals the LAD at grid
2460!--             point (k,j,i), in order to avoid that LAD is mistakenly mapped
2461!--             on top of a roof where (usually) is no LAD is defined.
[4770]2462!--             The same is also valid for bad_s.
[1721]2463                lad_local = lad_s(kk,j,i)
[4770]2464                IF ( lad_local == 0.0_wp  .AND.  lad_s(kk,j-1,i) > 0.0_wp  .AND.                   &
[4314]2465                     .NOT. building_edge_s )  lad_local = lad_s(kk,j-1,i)
[4770]2466
2467                bad_local = bad_s(kk,j,i)
2468                IF ( bad_local == 0.0_wp  .AND.  bad_s(kk,j-1,i) > 0.0_wp  .AND.                   &
2469                     .NOT. building_edge_s )  bad_local = bad_s(kk,j-1,i)
[4314]2470!
2471!--             In order to avoid that LAD is mistakenly considered at right-
2472!--             sided building edges (here the topography-top index for the
2473!--             u-component at index j,i is still on the building while the
2474!--             topography top for the scalar isn't), LAD is taken from grid
[4770]2475!--             point (j,i-1). The same is also valid for bad_s.
2476                IF ( lad_local > 0.0_wp  .AND.  lad_s(kk,j-1,i) == 0.0_wp  .AND.                   &
[4314]2477                     building_edge_n )  lad_local = lad_s(kk,j-1,i)
[4770]2478                IF ( bad_local > 0.0_wp  .AND.  bad_s(kk,j-1,i) == 0.0_wp  .AND.                   &
2479                     building_edge_n )  bad_local = bad_s(kk,j-1,i)
[1484]2480
2481                pre_tend = 0.0_wp
2482                pre_v = 0.0_wp
2483!
2484!--             Calculate preliminary value (pre_tend) of the tendency
[4342]2485                pre_tend = - canopy_drag_coeff *                               &
[4770]2486                             ( lad_local + bad_local ) *                       &
[1484]2487                             SQRT( ( 0.25_wp * ( u(k,j-1,i)   +                &
2488                                                 u(k,j-1,i+1) +                &
2489                                                 u(k,j,i)     +                &
2490                                                 u(k,j,i+1) )                  &
2491                                   )**2 +                                      &
2492                                   v(k,j,i)**2 +                               &
2493                                   ( 0.25_wp * ( w(k-1,j-1,i) +                &
2494                                                 w(k-1,j,i)   +                &
2495                                                 w(k,j-1,i)   +                &
2496                                                 w(k,j,i) )                    &
2497                                   )**2                                        &
2498                                 ) *                                           &
2499                             v(k,j,i)
2500
2501!
2502!--             Calculate preliminary new velocity, based on pre_tend
2503                pre_v = v(k,j,i) + dt_3d * pre_tend
2504!
2505!--             Compare sign of old velocity and new preliminary velocity,
2506!--             and in case the signs are different, limit the tendency
2507                IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
2508                   pre_tend = - v(k,j,i) * ddt_3d
2509                ELSE
2510                   pre_tend = pre_tend
2511                ENDIF
2512!
2513!--             Calculate final tendency
2514                tend(k,j,i) = tend(k,j,i) + pre_tend
2515             ENDDO
2516
2517
2518!
[142]2519!--       w-component
[1484]2520          CASE ( 3 )
[2232]2521!
2522!--          Determine topography-top index on w-grid
[4341]2523             DO  k = topo_top_ind(j,i,3) + 1, topo_top_ind(j,i,3) + pch_index_ji(j,i) - 1
[2317]2524
[4341]2525                kk = k - topo_top_ind(j,i,3)  !- lad arrays are defined flat
[138]2526
[1484]2527                pre_tend = 0.0_wp
2528                pre_w = 0.0_wp
[138]2529!
[1484]2530!--             Calculate preliminary value (pre_tend) of the tendency
[4342]2531                pre_tend = - canopy_drag_coeff *                               &
[4770]2532                             ( 0.5_wp *                                        &
2533                                ( lad_s(kk+1,j,i) + lad_s(kk,j,i) ) +          &
2534                               0.5_wp *                                        &
2535                                ( bad_s(kk+1,j,i) + bad_s(kk,j,i) ) ) *        &
2536                             SQRT( ( 0.25_wp * ( u(k,j,i)    +                 &
[1484]2537                                                 u(k,j,i+1)  +                 &
2538                                                 u(k+1,j,i)  +                 &
2539                                                 u(k+1,j,i+1) )                &
2540                                   )**2 +                                      &
2541                                   ( 0.25_wp * ( v(k,j,i)    +                 &
2542                                                 v(k,j+1,i)  +                 &
2543                                                 v(k+1,j,i)  +                 &
2544                                                 v(k+1,j+1,i) )                &
2545                                   )**2 +                                      &
2546                                   w(k,j,i)**2                                 &
2547                                 ) *                                           &
2548                             w(k,j,i)
2549!
2550!--             Calculate preliminary new velocity, based on pre_tend
2551                pre_w = w(k,j,i) + dt_3d * pre_tend
2552!
2553!--             Compare sign of old velocity and new preliminary velocity,
2554!--             and in case the signs are different, limit the tendency
2555                IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
2556                   pre_tend = - w(k,j,i) * ddt_3d
2557                ELSE
2558                   pre_tend = pre_tend
2559                ENDIF
2560!
2561!--             Calculate final tendency
2562                tend(k,j,i) = tend(k,j,i) + pre_tend
2563             ENDDO
2564
2565!
[153]2566!--       potential temperature
2567          CASE ( 4 )
[2232]2568!
2569!--          Determine topography-top index on scalar grid
[3449]2570             IF ( humidity ) THEN
[4341]2571                DO  k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2572                   kk = k - topo_top_ind(j,i,0)  !- lad arrays are defined flat
2573                   tend(k,j,i) = tend(k,j,i) + pcm_heating_rate(kk,j,i) -       &
2574                                 pcm_latent_rate(kk,j,i)
[3449]2575                ENDDO
2576             ELSE
[4341]2577                DO  k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2578                   kk = k - topo_top_ind(j,i,0)  !- lad arrays are defined flat
2579                   tend(k,j,i) = tend(k,j,i) + pcm_heating_rate(kk,j,i)
[3449]2580                ENDDO
2581             ENDIF
[153]2582
2583!
[1960]2584!--       humidity
[153]2585          CASE ( 5 )
[2232]2586!
2587!--          Determine topography-top index on scalar grid
[4341]2588             DO  k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
2589                kk = k - topo_top_ind(j,i,0)  !- lad arrays are defined flat
[3449]2590                IF ( .NOT. plant_canopy_transpiration ) THEN
[4341]2591                   ! pcm_transpiration_rate is calculated in radiation model
[3449]2592                   ! in case of plant_canopy_transpiration = .T.
2593                   ! to include also the dependecy to the radiation
2594                   ! in the plant canopy box
[4342]2595                   pcm_transpiration_rate(kk,j,i) = - leaf_scalar_exch_coeff   &
[3582]2596                                    * lad_s(kk,j,i) *                          &
2597                                    SQRT( ( 0.5_wp * ( u(k,j,i) +              &
2598                                                       u(k,j,i+1) )            &
2599                                          )**2  +                              &
2600                                          ( 0.5_wp * ( v(k,j,i) +              &
2601                                                       v(k,j+1,i) )            &
2602                                          )**2 +                               &
2603                                          ( 0.5_wp * ( w(k-1,j,i) +            &
2604                                                       w(k,j,i) )              &
2605                                          )**2                                 &
2606                                        ) *                                    &
[4342]2607                                    ( q(k,j,i) - leaf_surface_conc )
[3449]2608                ENDIF
[2232]2609
[4341]2610                tend(k,j,i) = tend(k,j,i) + pcm_transpiration_rate(kk,j,i)
[3014]2611
[153]2612             ENDDO   
2613
2614!
[142]2615!--       sgs-tke
[1484]2616          CASE ( 6 )
[2232]2617!
2618!--          Determine topography-top index on scalar grid
[4341]2619             DO  k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
[2317]2620
[4341]2621                kk = k - topo_top_ind(j,i,0)
[1484]2622                tend(k,j,i) = tend(k,j,i) -                                    &
[4342]2623                                 2.0_wp * canopy_drag_coeff *                  &
[4770]2624                                 ( lad_s(kk,j,i) + bad_s(kk,j,i) ) *           &
[1484]2625                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
2626                                                    u(k,j,i+1) )               &
2627                                       )**2 +                                  & 
2628                                       ( 0.5_wp * ( v(k,j,i) +                 &
2629                                                    v(k,j+1,i) )               &
2630                                       )**2 +                                  &
2631                                       ( 0.5_wp * ( w(k,j,i) +                 &
2632                                                    w(k+1,j,i) )               &
2633                                       )**2                                    &
2634                                     ) *                                       &
2635                                 e(k,j,i)
2636             ENDDO
[1960]2637!
2638!--       scalar concentration
2639          CASE ( 7 )
[2232]2640!
2641!--          Determine topography-top index on scalar grid
[4341]2642             DO  k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
[2317]2643
[4341]2644                kk = k - topo_top_ind(j,i,0)
[1960]2645                tend(k,j,i) = tend(k,j,i) -                                    &
[4342]2646                                 leaf_scalar_exch_coeff *                      &
[1960]2647                                 lad_s(kk,j,i) *                               &
2648                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
2649                                                    u(k,j,i+1) )               &
2650                                       )**2  +                                 &
2651                                       ( 0.5_wp * ( v(k,j,i) +                 &
2652                                                    v(k,j+1,i) )               &
2653                                       )**2 +                                  &
2654                                       ( 0.5_wp * ( w(k-1,j,i) +               &
2655                                                    w(k,j,i) )                 &
2656                                       )**2                                    &
2657                                     ) *                                       &
[4342]2658                                 ( s(k,j,i) - leaf_surface_conc )
[1960]2659             ENDDO               
[138]2660
[142]2661       CASE DEFAULT
[138]2662
[257]2663          WRITE( message_string, * ) 'wrong component: ', component
[1826]2664          CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 ) 
[138]2665
[142]2666       END SELECT
[138]2667
[1826]2668    END SUBROUTINE pcm_tendency_ij
[138]2669
[4360]2670!------------------------------------------------------------------------------!
2671! Description:
2672! ------------
[4514]2673!> Subroutine writes global restart data
2674!------------------------------------------------------------------------------!
2675    SUBROUTINE pcm_wrd_global
2676
2677       IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
2678
2679          CALL wrd_write_string( 'pch_index' )
2680          WRITE ( 14 )  pch_index
2681
[4535]2682       ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
[4514]2683
2684          CALL wrd_mpi_io( 'pch_index', pch_index )
2685
2686       ENDIF
2687
2688    END SUBROUTINE pcm_wrd_global
2689
2690!------------------------------------------------------------------------------!
2691! Description:
2692! ------------
[4360]2693!> Subroutine writes local (subdomain) restart data
2694!------------------------------------------------------------------------------!
2695    SUBROUTINE pcm_wrd_local
[2007]2696
[4525]2697       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  tmp_3d  !< temporary array to store pcm data with
2698                                                           !< non-standard vertical index bounds
2699
[4495]2700       IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
[4360]2701
[4495]2702          IF ( ALLOCATED( pcm_heatrate_av ) )  THEN
2703             CALL wrd_write_string( 'pcm_heatrate_av' )
2704             WRITE ( 14 )  pcm_heatrate_av
2705          ENDIF
[4360]2706
[4495]2707          IF ( ALLOCATED( pcm_latentrate_av ) )  THEN
2708             CALL wrd_write_string( 'pcm_latentrate_av' )
2709             WRITE ( 14 )  pcm_latentrate_av
2710          ENDIF
2711
2712          IF ( ALLOCATED( pcm_transpirationrate_av ) )  THEN
2713             CALL wrd_write_string( 'pcm_transpirationrate_av' )
2714             WRITE ( 14 )  pcm_transpirationrate_av
2715          ENDIF
2716
[4535]2717       ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
[4495]2718
[4525]2719!
2720!--       Plant canopy arrays have non standard reduced vertical index bounds. They are stored with
2721!--       full vertical bounds (bzb:nzt+1) in the restart file and must be re-stored before writing.
2722          IF ( ALLOCATED( pcm_heatrate_av ) )  THEN
2723             ALLOCATE( tmp_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2724             tmp_3d(nzb:pch_index,:,:)     = pcm_heatrate_av
2725             tmp_3d(pch_index+1:nzt+1,:,:) = 0.0_wp
2726             CALL wrd_mpi_io( 'pcm_heatrate_av', tmp_3d )
2727             DEALLOCATE( tmp_3d )
2728          ENDIF
[4495]2729          IF ( ALLOCATED( pcm_latentrate_av ) )  THEN
[4525]2730             ALLOCATE( tmp_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2731             tmp_3d(nzb:pch_index,:,:)     = pcm_latentrate_av
2732             tmp_3d(pch_index+1:nzt+1,:,:) = 0.0_wp
2733             CALL wrd_mpi_io( 'pcm_latentrate_av', tmp_3d )
2734             DEALLOCATE( tmp_3d )
[4495]2735          ENDIF
2736          IF ( ALLOCATED( pcm_transpirationrate_av ) )  THEN
[4525]2737             ALLOCATE( tmp_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2738             tmp_3d(nzb:pch_index,:,:)     = pcm_transpirationrate_av
2739             tmp_3d(pch_index+1:nzt+1,:,:) = 0.0_wp
2740             CALL wrd_mpi_io( 'pcm_transpirationrate_av', tmp_3d )
2741             DEALLOCATE( tmp_3d )
[4495]2742          ENDIF
2743
[4360]2744       ENDIF
2745
2746    END SUBROUTINE pcm_wrd_local
2747
2748
[138]2749 END MODULE plant_canopy_model_mod
Note: See TracBrowser for help on using the repository browser.