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

Last change on this file since 3768 was 3761, checked in by raasch, 2 years ago

unused variables removed, OpenACC directives re-formatted, statements added to avoid compiler warnings

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