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

Last change on this file since 3609 was 3589, checked in by suehring, 5 years ago

Remove erroneous UTF encoding; last commit documented

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