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

Last change on this file since 3703 was 3685, checked in by knoop, 5 years ago

Some interface calls moved to module_interface + cleanup

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