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

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

Remove write statement

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