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

Last change on this file since 4149 was 4127, checked in by suehring, 5 years ago

Merge with branch resler: biomet- output of bio_mrt added; plant_canopy - separate vertical dimension for 3D output (to save disk space); radiation - remove unused plant canopy variables; urban-surface model - do not add anthropogenic heat during wall spin-up

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