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

Last change on this file since 4309 was 4309, checked in by suehring, 18 months ago

Synthetic turbulence generator: Computation of velocity seeds optimized. This implies that random numbers are computed now using the parallel random number generator. Random number are now only computed and normalized locally, while distributed over all mpi ranks afterwards, instead of computing random numbers on a global array. urther, the number of calls for the time-consuming velocity-seed generation is reduced - now the left and right, as well as the north and south boundary share the same velocity-seed matrices.

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