source: palm/trunk/SOURCE/plant_canopy_model.f90 @ 1721

Last change on this file since 1721 was 1721, checked in by raasch, 6 years ago

bugfixes: shf is reduced in areas covered with canopy only, canopy is set on top of topography

  • Property svn:keywords set to Id
File size: 42.2 KB
Line 
1!> @file plant_canopy_model.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2014 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21! bugfixes: shf is reduced in areas covered with canopy only,
22!           canopy is set on top of topography
23!
24! Former revisions:
25! -----------------
26! $Id: plant_canopy_model.f90 1721 2015-11-16 12:56:48Z raasch $
27!
28! 1682 2015-10-07 23:56:08Z knoop
29! Code annotations made doxygen readable
30!
31! 1484 2014-10-21 10:53:05Z kanani
32! Changes due to new module structure of the plant canopy model:
33!   module plant_canopy_model_mod now contains a subroutine for the
34!   initialization of the canopy model (init_plant_canopy),
35!   limitation of the canopy drag (previously accounted for by calculation of
36!   a limiting canopy timestep for the determination of the maximum LES timestep
37!   in subroutine timestep) is now realized by the calculation of pre-tendencies
38!   and preliminary velocities in subroutine plant_canopy_model,
39!   some redundant MPI communication removed in subroutine init_plant_canopy
40!   (was previously in init_3d_model),
41!   unnecessary 3d-arrays lad_u, lad_v, lad_w removed - lad information on the
42!   respective grid is now provided only by lad_s (e.g. in the calculation of
43!   the tendency terms or of cum_lai_hf),
44!   drag_coefficient, lai, leaf_surface_concentration,
45!   scalar_exchange_coefficient, sec and sls renamed to canopy_drag_coeff,
46!   cum_lai_hf, leaf_surface_conc, leaf_scalar_exch_coeff, lsec and lsc,
47!   respectively,
48!   unnecessary 3d-arrays cdc, lsc and lsec now defined as single-value constants,
49!   USE-statements and ONLY-lists modified accordingly
50!
51! 1340 2014-03-25 19:45:13Z kanani
52! REAL constants defined as wp-kind
53!
54! 1320 2014-03-20 08:40:49Z raasch
55! ONLY-attribute added to USE-statements,
56! kind-parameters added to all INTEGER and REAL declaration statements,
57! kinds are defined in new module kinds,
58! old module precision_kind is removed,
59! revision history before 2012 removed,
60! comment fields (!:) to be used for variable explanations added to
61! all variable declaration statements
62!
63! 1036 2012-10-22 13:43:42Z raasch
64! code put under GPL (PALM 3.9)
65!
66! 138 2007-11-28 10:03:58Z letzel
67! Initial revision
68!
69! Description:
70! ------------
71!> 1) Initialization of the canopy model, e.g. construction of leaf area density
72!> profile (subroutine init_plant_canopy).
73!> 2) Calculation of sinks and sources of momentum, heat and scalar concentration
74!> due to canopy elements (subroutine plant_canopy_model).
75!------------------------------------------------------------------------------!
76 MODULE plant_canopy_model_mod
77 
78    USE arrays_3d,                                                             &
79        ONLY:  dzu, dzw, e, q, shf, tend, u, v, w, zu, zw
80
81    USE indices,                                                               &
82        ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
83               nz, nzb, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt
84
85    USE kinds
86
87
88    IMPLICIT NONE
89
90
91    CHARACTER (LEN=20)   ::  canopy_mode = 'block' !< canopy coverage
92
93    INTEGER(iwp) ::  pch_index = 0                 !< plant canopy height/top index
94    INTEGER(iwp) ::                                                            &
95       lad_vertical_gradient_level_ind(10) = -9999 !< lad-profile levels (index)
96
97    LOGICAL ::  calc_beta_lad_profile = .FALSE. !< switch for calc. of lad from beta func.
98    LOGICAL ::  plant_canopy = .FALSE.          !< switch for use of canopy model
99
100    REAL(wp) ::  alpha_lad = 9999999.9_wp   !< coefficient for lad calculation
101    REAL(wp) ::  beta_lad = 9999999.9_wp    !< coefficient for lad calculation
102    REAL(wp) ::  canopy_drag_coeff = 0.0_wp !< canopy drag coefficient (parameter)
103    REAL(wp) ::  cdc = 0.0_wp               !< canopy drag coeff. (abbreviation used in equations)
104    REAL(wp) ::  cthf = 0.0_wp              !< canopy top heat flux
105    REAL(wp) ::  dt_plant_canopy = 0.0_wp   !< timestep account. for canopy drag
106    REAL(wp) ::  lad_surface = 0.0_wp       !< lad surface value
107    REAL(wp) ::  lai_beta = 0.0_wp          !< leaf area index (lai) for lad calc.
108    REAL(wp) ::                                                                &
109       leaf_scalar_exch_coeff = 0.0_wp      !< canopy scalar exchange coeff.
110    REAL(wp) ::                                                                &
111       leaf_surface_conc = 0.0_wp           !< leaf surface concentration
112    REAL(wp) ::  lsec = 0.0_wp              !< leaf scalar exchange coeff.
113    REAL(wp) ::  lsc = 0.0_wp               !< leaf surface concentration
114
115    REAL(wp) ::                                                                &
116       lad_vertical_gradient(10) = 0.0_wp              !< lad gradient
117    REAL(wp) ::                                                                &
118       lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m)
119
120    REAL(wp), DIMENSION(:), ALLOCATABLE ::  lad            !< leaf area density
121    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pre_lad        !< preliminary lad
122   
123    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
124       canopy_heat_flux                                    !< canopy heat flux
125    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  cum_lai_hf !< cumulative lai for heatflux calc.
126    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  lad_s      !< lad on scalar-grid
127
128
129    SAVE
130
131
132    PRIVATE
133    PUBLIC alpha_lad, beta_lad, calc_beta_lad_profile, canopy_drag_coeff,      &
134           canopy_mode, cdc, cthf, dt_plant_canopy, init_plant_canopy, lad,    &
135           lad_s, lad_surface, lad_vertical_gradient,                          &
136           lad_vertical_gradient_level, lad_vertical_gradient_level_ind,       &
137           lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, lsc, lsec,     &
138           pch_index, plant_canopy, plant_canopy_model
139
140
141    INTERFACE init_plant_canopy
142       MODULE PROCEDURE init_plant_canopy
143    END INTERFACE init_plant_canopy
144
145    INTERFACE plant_canopy_model
146       MODULE PROCEDURE plant_canopy_model
147       MODULE PROCEDURE plant_canopy_model_ij
148    END INTERFACE plant_canopy_model
149
150
151
152
153 CONTAINS
154
155
156!------------------------------------------------------------------------------!
157! Description:
158! ------------
159!> Initialization of the plant canopy model
160!------------------------------------------------------------------------------!
161    SUBROUTINE init_plant_canopy
162   
163
164       USE control_parameters,                                                 &
165           ONLY: dz, ocean, passive_scalar
166
167
168       IMPLICIT NONE
169
170       INTEGER(iwp) ::  i !< running index
171       INTEGER(iwp) ::  j !< running index
172       INTEGER(iwp) ::  k !< running index
173
174       REAL(wp) ::  int_bpdf      !< vertical integral for lad-profile construction
175       REAL(wp) ::  dzh           !< vertical grid spacing in units of canopy height
176       REAL(wp) ::  gradient      !< gradient for lad-profile construction
177       REAL(wp) ::  canopy_height !< canopy height for lad-profile construction
178
179!
180!--    Allocate one-dimensional arrays for the computation of the
181!--    leaf area density (lad) profile
182       ALLOCATE( lad(0:nz+1), pre_lad(0:nz+1) )
183       lad = 0.0_wp
184       pre_lad = 0.0_wp
185
186!
187!--    Compute the profile of leaf area density used in the plant
188!--    canopy model. The profile can either be constructed from
189!--    prescribed vertical gradients of the leaf area density or by
190!--    using a beta probability density function (see e.g. Markkanen et al.,
191!--    2003: Boundary-Layer Meteorology, 106, 437-459)
192       IF (  .NOT.  calc_beta_lad_profile )  THEN   
193
194!
195!--       Use vertical gradients for lad-profile construction   
196          i = 1
197          gradient = 0.0_wp
198
199          IF (  .NOT.  ocean )  THEN
200
201             lad(0) = lad_surface
202             lad_vertical_gradient_level_ind(1) = 0
203 
204             DO k = 1, pch_index
205                IF ( i < 11 )  THEN
206                   IF ( lad_vertical_gradient_level(i) < zu(k)  .AND.          &
207                        lad_vertical_gradient_level(i) >= 0.0_wp )  THEN
208                      gradient = lad_vertical_gradient(i)
209                      lad_vertical_gradient_level_ind(i) = k - 1
210                      i = i + 1
211                   ENDIF
212                ENDIF
213                IF ( gradient /= 0.0_wp )  THEN
214                   IF ( k /= 1 )  THEN
215                      lad(k) = lad(k-1) + dzu(k) * gradient
216                   ELSE
217                      lad(k) = lad_surface + dzu(k) * gradient
218                   ENDIF
219                ELSE
220                   lad(k) = lad(k-1)
221                ENDIF
222             ENDDO
223
224          ENDIF
225
226!
227!--       In case of no given leaf area density gradients, choose a vanishing
228!--       gradient. This information is used for the HEADER and the RUN_CONTROL
229!--       file.
230          IF ( lad_vertical_gradient_level(1) == -9999999.9_wp )  THEN
231             lad_vertical_gradient_level(1) = 0.0_wp
232          ENDIF
233
234       ELSE
235
236!
237!--       Use beta function for lad-profile construction
238          int_bpdf = 0.0_wp
239          canopy_height = pch_index * dz
240
241          DO k = nzb, pch_index
242             int_bpdf = int_bpdf +                                             &
243                           ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) * &
244                           ( ( 1.0_wp - ( zw(k) / canopy_height ) )**( beta_lad-1.0_wp ) ) *   &
245                           ( ( zw(k+1)-zw(k) ) / canopy_height ) )
246          ENDDO
247
248!
249!--       Preliminary lad profile (defined on w-grid)
250          DO k = nzb, pch_index
251             pre_lad(k) =  lai_beta *                                              &
252                           (   ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) *          &
253                               ( ( 1.0_wp - ( zw(k) / canopy_height ) )**( beta_lad-1.0_wp ) ) /   &
254                               int_bpdf                                            &
255                           ) / canopy_height
256          ENDDO
257
258!
259!--       Final lad profile (defined on scalar-grid level, since most prognostic
260!--       quantities are defined there, hence, less interpolation is required
261!--       when calculating the canopy tendencies)
262          lad(0) = pre_lad(0)
263          DO k = nzb+1, pch_index
264             lad(k) = 0.5 * ( pre_lad(k-1) + pre_lad(k) )
265          ENDDO         
266
267       ENDIF
268
269!
270!--    Allocate 3D-array for the leaf area density (lad_s). In case of a
271!--    prescribed canopy-top heat flux (cthf), allocate 3D-arrays for
272!--    the cumulative leaf area index (cum_lai_hf) and the canopy heat flux.
273       ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
274
275       IF ( cthf /= 0.0_wp )  THEN
276          ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                 &
277                    canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
278       ENDIF
279
280!
281!--    Initialize canopy parameters cdc (canopy drag coefficient),
282!--    lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration)
283!--    with the prescribed values
284       cdc = canopy_drag_coeff
285       lsec = leaf_scalar_exch_coeff
286       lsc = leaf_surface_conc
287
288!
289!--    Initialization of the canopy coverage in the model domain:
290!--    Setting the parameter canopy_mode = 'block' initializes a canopy, which
291!--    fully covers the domain surface
292       SELECT CASE ( TRIM( canopy_mode ) )
293
294          CASE( 'block' )
295
296             DO  i = nxlg, nxrg
297                DO  j = nysg, nyng
298                   lad_s(:,j,i) = lad(:)
299                ENDDO
300             ENDDO
301
302          CASE DEFAULT
303
304!
305!--       The DEFAULT case is reached either if the parameter
306!--       canopy mode contains a wrong character string or if the
307!--       user has coded a special case in the user interface.
308!--       There, the subroutine user_init_plant_canopy checks
309!--       which of these two conditions applies.
310          CALL user_init_plant_canopy
311 
312       END SELECT
313
314!
315!--    Initialization of the canopy heat source distribution
316       IF ( cthf /= 0.0_wp )  THEN
317!
318!--       Piecewise calculation of the leaf area index by vertical
319!--       integration of the leaf area density
320          cum_lai_hf(:,:,:) = 0.0_wp
321          DO  i = nxlg, nxrg
322             DO  j = nysg, nyng
323                DO  k = pch_index-1, 0, -1
324                   IF ( k == pch_index-1 )  THEN
325                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
326                         ( 0.5_wp * lad_s(k+1,j,i) *                           &
327                           ( zw(k+1) - zu(k+1) ) )  +                          &
328                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
329                                                 lad_s(k,j,i) ) +              &
330                                      lad_s(k+1,j,i) ) *                       &
331                           ( zu(k+1) - zw(k) ) ) 
332                   ELSE
333                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
334                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) +              &
335                                                 lad_s(k+1,j,i) ) +            &
336                                      lad_s(k+1,j,i) ) *                       &
337                           ( zw(k+1) - zu(k+1) ) )  +                          &
338                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
339                                                 lad_s(k,j,i) ) +              &
340                                      lad_s(k+1,j,i) ) *                       &
341                           ( zu(k+1) - zw(k) ) )
342                   ENDIF
343                ENDDO
344             ENDDO
345          ENDDO
346
347!
348!--       Calculation of the upward kinematic vertical heat flux within the
349!--       canopy
350          DO  i = nxlg, nxrg
351             DO  j = nysg, nyng
352                DO  k = 0, pch_index
353                   canopy_heat_flux(k,j,i) = cthf *                            &
354                                             exp( -0.6_wp * cum_lai_hf(k,j,i) )
355                ENDDO
356             ENDDO
357          ENDDO
358
359!
360!--       In areas covered with canopy, the surface heat flux is set to
361!--       the surface value of the above calculated in-canopy heat flux
362!--       distribution
363          DO  i = nxlg,nxrg
364             DO  j = nysg, nyng
365                IF ( canopy_heat_flux(0,j,i) /= cthf )  THEN
366                   shf(j,i) = canopy_heat_flux(0,j,i)
367                ENDIF
368             ENDDO
369          ENDDO
370
371       ENDIF
372
373
374
375    END SUBROUTINE init_plant_canopy
376
377
378
379!------------------------------------------------------------------------------!
380! Description:
381! ------------
382!> Calculation of the tendency terms, accounting for the effect of the plant
383!> canopy on momentum and scalar quantities.
384!>
385!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
386!> (defined on scalar grid), as initialized in subroutine init_plant_canopy.
387!> The lad on the w-grid is vertically interpolated from the surrounding
388!> lad_s. The upper boundary of the canopy is defined on the w-grid at
389!> k = pch_index. Here, the lad is zero.
390!>
391!> The canopy drag must be limited (previously accounted for by calculation of
392!> a limiting canopy timestep for the determination of the maximum LES timestep
393!> in subroutine timestep), since it is physically impossible that the canopy
394!> drag alone can locally change the sign of a velocity component. This
395!> limitation is realized by calculating preliminary tendencies and velocities.
396!> It is subsequently checked if the preliminary new velocity has a different
397!> sign than the current velocity. If so, the tendency is limited in a way that
398!> the velocity can at maximum be reduced to zero by the canopy drag.
399!>
400!>
401!> Call for all grid points
402!------------------------------------------------------------------------------!
403    SUBROUTINE plant_canopy_model( component )
404
405
406       USE control_parameters,                                                 &
407           ONLY:  dt_3d, message_string
408
409       USE kinds
410
411       IMPLICIT NONE
412
413       INTEGER(iwp) ::  component !< prognostic variable (u,v,w,pt,q,e)
414       INTEGER(iwp) ::  i         !< running index
415       INTEGER(iwp) ::  j         !< running index
416       INTEGER(iwp) ::  k         !< running index
417       INTEGER(iwp) ::  kk        !< running index for flat lad arrays
418
419       REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
420       REAL(wp) ::  lad_local !< local lad value
421       REAL(wp) ::  pre_tend  !< preliminary tendency
422       REAL(wp) ::  pre_u     !< preliminary u-value
423       REAL(wp) ::  pre_v     !< preliminary v-value
424       REAL(wp) ::  pre_w     !< preliminary w-value
425
426
427       ddt_3d = 1.0_wp / dt_3d
428 
429!
430!--    Compute drag for the three velocity components and the SGS-TKE:
431       SELECT CASE ( component )
432
433!
434!--       u-component
435          CASE ( 1 )
436             DO  i = nxlu, nxr
437                DO  j = nys, nyn
438                   DO  k = nzb_u_inner(j,i)+1, nzb_u_inner(j,i)+pch_index
439
440                      kk = k - nzb_u_inner(j,i)  !- lad arrays are defined flat
441!
442!--                   In order to create sharp boundaries of the plant canopy,
443!--                   the lad on the u-grid at index (k,j,i) is equal to
444!--                   lad_s(k,j,i), rather than being interpolated from the
445!--                   surrounding lad_s, because this would yield smaller lad
446!--                   at the canopy boundaries than inside of the canopy.
447!--                   For the same reason, the lad at the rightmost(i+1)canopy
448!--                   boundary on the u-grid equals lad_s(k,j,i).
449                      lad_local = lad_s(kk,j,i)
450                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp )&
451                      THEN
452                         lad_local = lad_s(kk,j,i-1)
453                      ENDIF
454
455                      pre_tend = 0.0_wp
456                      pre_u = 0.0_wp
457!
458!--                   Calculate preliminary value (pre_tend) of the tendency
459                      pre_tend = - cdc *                                       &
460                                   lad_local *                                 &
461                                   SQRT( u(k,j,i)**2 +                         &
462                                         ( 0.25_wp * ( v(k,j,i-1) +            &
463                                                       v(k,j,i)   +            &
464                                                       v(k,j+1,i) +            &
465                                                       v(k,j+1,i-1) )          &
466                                         )**2 +                                &
467                                         ( 0.25_wp * ( w(k-1,j,i-1) +          &
468                                                       w(k-1,j,i)   +          &
469                                                       w(k,j,i-1)   +          &
470                                                       w(k,j,i) )              &
471                                         )**2                                  &
472                                       ) *                                     &
473                                   u(k,j,i)
474
475!
476!--                   Calculate preliminary new velocity, based on pre_tend
477                      pre_u = u(k,j,i) + dt_3d * pre_tend
478!
479!--                   Compare sign of old velocity and new preliminary velocity,
480!--                   and in case the signs are different, limit the tendency
481                      IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
482                         pre_tend = - u(k,j,i) * ddt_3d
483                      ELSE
484                         pre_tend = pre_tend
485                      ENDIF
486!
487!--                   Calculate final tendency
488                      tend(k,j,i) = tend(k,j,i) + pre_tend
489
490                   ENDDO
491                ENDDO
492             ENDDO
493
494!
495!--       v-component
496          CASE ( 2 )
497             DO  i = nxl, nxr
498                DO  j = nysv, nyn
499                   DO  k = nzb_v_inner(j,i)+1, nzb_v_inner(j,i)+pch_index
500
501                      kk = k - nzb_v_inner(j,i)  !- lad arrays are defined flat
502!
503!--                   In order to create sharp boundaries of the plant canopy,
504!--                   the lad on the v-grid at index (k,j,i) is equal to
505!--                   lad_s(k,j,i), rather than being interpolated from the
506!--                   surrounding lad_s, because this would yield smaller lad
507!--                   at the canopy boundaries than inside of the canopy.
508!--                   For the same reason, the lad at the northmost(j+1) canopy
509!--                   boundary on the v-grid equals lad_s(k,j,i).
510                      lad_local = lad_s(kk,j,i)
511                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp )&
512                      THEN
513                         lad_local = lad_s(kk,j-1,i)
514                      ENDIF
515
516                      pre_tend = 0.0_wp
517                      pre_v = 0.0_wp
518!
519!--                   Calculate preliminary value (pre_tend) of the tendency
520                      pre_tend = - cdc *                                       &
521                                   lad_local *                                 &
522                                   SQRT( ( 0.25_wp * ( u(k,j-1,i)   +          &
523                                                       u(k,j-1,i+1) +          &
524                                                       u(k,j,i)     +          &
525                                                       u(k,j,i+1) )            &
526                                         )**2 +                                &
527                                         v(k,j,i)**2 +                         &
528                                         ( 0.25_wp * ( w(k-1,j-1,i) +          &
529                                                       w(k-1,j,i)   +          &
530                                                       w(k,j-1,i)   +          &
531                                                       w(k,j,i) )              &
532                                         )**2                                  &
533                                       ) *                                     &
534                                   v(k,j,i)
535
536!
537!--                   Calculate preliminary new velocity, based on pre_tend
538                      pre_v = v(k,j,i) + dt_3d * pre_tend
539!
540!--                   Compare sign of old velocity and new preliminary velocity,
541!--                   and in case the signs are different, limit the tendency
542                      IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
543                         pre_tend = - v(k,j,i) * ddt_3d
544                      ELSE
545                         pre_tend = pre_tend
546                      ENDIF
547!
548!--                   Calculate final tendency
549                      tend(k,j,i) = tend(k,j,i) + pre_tend
550
551                   ENDDO
552                ENDDO
553             ENDDO
554
555!
556!--       w-component
557          CASE ( 3 )
558             DO  i = nxl, nxr
559                DO  j = nys, nyn
560                   DO  k = nzb_w_inner(j,i)+1, nzb_w_inner(j,i)+pch_index-1
561
562                      kk = k - nzb_w_inner(j,i)  !- lad arrays are defined flat
563
564                      pre_tend = 0.0_wp
565                      pre_w = 0.0_wp
566!
567!--                   Calculate preliminary value (pre_tend) of the tendency
568                      pre_tend = - cdc *                                       &
569                                   (0.5_wp *                                   &
570                                      ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) *   &
571                                   SQRT( ( 0.25_wp * ( u(k,j,i)   +            &
572                                                       u(k,j,i+1) +            &
573                                                       u(k+1,j,i) +            &
574                                                       u(k+1,j,i+1) )          &
575                                         )**2 +                                &
576                                         ( 0.25_wp * ( v(k,j,i)   +            &
577                                                       v(k,j+1,i) +            &
578                                                       v(k+1,j,i) +            &
579                                                       v(k+1,j+1,i) )          &
580                                         )**2 +                                &
581                                         w(k,j,i)**2                           &
582                                       ) *                                     &
583                                   w(k,j,i)
584!
585!--                   Calculate preliminary new velocity, based on pre_tend
586                      pre_w = w(k,j,i) + dt_3d * pre_tend
587!
588!--                   Compare sign of old velocity and new preliminary velocity,
589!--                   and in case the signs are different, limit the tendency
590                      IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
591                         pre_tend = - w(k,j,i) * ddt_3d
592                      ELSE
593                         pre_tend = pre_tend
594                      ENDIF
595!
596!--                   Calculate final tendency
597                      tend(k,j,i) = tend(k,j,i) + pre_tend
598
599                   ENDDO
600                ENDDO
601             ENDDO
602
603!
604!--       potential temperature
605          CASE ( 4 )
606             DO  i = nxl, nxr
607                DO  j = nys, nyn
608                   DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
609                      kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
610                      tend(k,j,i) = tend(k,j,i) +                              &
611                                       ( canopy_heat_flux(kk,j,i) -            &
612                                         canopy_heat_flux(kk-1,j,i) ) / dzw(k)
613                   ENDDO
614                ENDDO
615             ENDDO
616
617!
618!--       scalar concentration
619          CASE ( 5 )
620             DO  i = nxl, nxr
621                DO  j = nys, nyn
622                   DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
623                      kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
624                      tend(k,j,i) = tend(k,j,i) -                              &
625                                       lsec *                                  &
626                                       lad_s(kk,j,i) *                         &
627                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
628                                                          u(k,j,i+1) )         &
629                                             )**2 +                            &
630                                             ( 0.5_wp * ( v(k,j,i) +           &
631                                                          v(k,j+1,i) )         &
632                                             )**2 +                            &
633                                             ( 0.5_wp * ( w(k-1,j,i) +         & 
634                                                          w(k,j,i) )           &
635                                             )**2                              &
636                                           ) *                                 &
637                                       ( q(k,j,i) - lsc )
638                   ENDDO
639                ENDDO
640             ENDDO
641
642!
643!--       sgs-tke
644          CASE ( 6 )
645             DO  i = nxl, nxr
646                DO  j = nys, nyn
647                   DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
648                      kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
649                      tend(k,j,i) = tend(k,j,i) -                              &
650                                       2.0_wp * cdc *                          &
651                                       lad_s(kk,j,i) *                         &
652                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
653                                                          u(k,j,i+1) )         &
654                                             )**2 +                            &
655                                             ( 0.5_wp * ( v(k,j,i) +           &
656                                                          v(k,j+1,i) )         &
657                                             )**2 +                            &
658                                             ( 0.5_wp * ( w(k,j,i) +           &
659                                                          w(k+1,j,i) )         &
660                                             )**2                              &
661                                           ) *                                 &
662                                       e(k,j,i)
663                   ENDDO
664                ENDDO
665             ENDDO
666
667
668          CASE DEFAULT
669
670             WRITE( message_string, * ) 'wrong component: ', component
671             CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 ) 
672
673       END SELECT
674
675    END SUBROUTINE plant_canopy_model
676
677
678!------------------------------------------------------------------------------!
679! Description:
680! ------------
681!> Calculation of the tendency terms, accounting for the effect of the plant
682!> canopy on momentum and scalar quantities.
683!>
684!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
685!> (defined on scalar grid), as initialized in subroutine init_plant_canopy.
686!> The lad on the w-grid is vertically interpolated from the surrounding
687!> lad_s. The upper boundary of the canopy is defined on the w-grid at
688!> k = pch_index. Here, the lad is zero.
689!>
690!> The canopy drag must be limited (previously accounted for by calculation of
691!> a limiting canopy timestep for the determination of the maximum LES timestep
692!> in subroutine timestep), since it is physically impossible that the canopy
693!> drag alone can locally change the sign of a velocity component. This
694!> limitation is realized by calculating preliminary tendencies and velocities.
695!> It is subsequently checked if the preliminary new velocity has a different
696!> sign than the current velocity. If so, the tendency is limited in a way that
697!> the velocity can at maximum be reduced to zero by the canopy drag.
698!>
699!>
700!> Call for grid point i,j
701!------------------------------------------------------------------------------!
702    SUBROUTINE plant_canopy_model_ij( i, j, component )
703
704
705       USE control_parameters,                                                 &
706           ONLY:  dt_3d, message_string
707
708       USE kinds
709
710       IMPLICIT NONE
711
712       INTEGER(iwp) ::  component !< prognostic variable (u,v,w,pt,q,e)
713       INTEGER(iwp) ::  i         !< running index
714       INTEGER(iwp) ::  j         !< running index
715       INTEGER(iwp) ::  k         !< running index
716       INTEGER(iwp) ::  kk        !< running index for flat lad arrays
717
718       REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
719       REAL(wp) ::  lad_local !< local lad value
720       REAL(wp) ::  pre_tend  !< preliminary tendency
721       REAL(wp) ::  pre_u     !< preliminary u-value
722       REAL(wp) ::  pre_v     !< preliminary v-value
723       REAL(wp) ::  pre_w     !< preliminary w-value
724
725
726       ddt_3d = 1.0_wp / dt_3d
727
728!
729!--    Compute drag for the three velocity components and the SGS-TKE
730       SELECT CASE ( component )
731
732!
733!--       u-component
734          CASE ( 1 )
735             DO  k = nzb_u_inner(j,i)+1, nzb_u_inner(j,i)+pch_index
736
737                kk = k - nzb_u_inner(j,i)  !- lad arrays are defined flat
738!
739!--             In order to create sharp boundaries of the plant canopy,
740!--             the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
741!--             rather than being interpolated from the surrounding lad_s,
742!--             because this would yield smaller lad at the canopy boundaries
743!--             than inside of the canopy.
744!--             For the same reason, the lad at the rightmost(i+1)canopy
745!--             boundary on the u-grid equals lad_s(k,j,i).
746                lad_local = lad_s(kk,j,i)
747                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp )  THEN
748                   lad_local = lad_s(kk,j,i-1)
749                ENDIF
750
751                pre_tend = 0.0_wp
752                pre_u = 0.0_wp
753!
754!--             Calculate preliminary value (pre_tend) of the tendency
755                pre_tend = - cdc *                                             &
756                             lad_local *                                       &   
757                             SQRT( u(k,j,i)**2 +                               &
758                                   ( 0.25_wp * ( v(k,j,i-1)  +                 &
759                                                 v(k,j,i)    +                 &
760                                                 v(k,j+1,i)  +                 &
761                                                 v(k,j+1,i-1) )                &
762                                   )**2 +                                      &
763                                   ( 0.25_wp * ( w(k-1,j,i-1) +                &
764                                                 w(k-1,j,i)   +                &
765                                                 w(k,j,i-1)   +                &
766                                                 w(k,j,i) )                    &
767                                   )**2                                        &
768                                 ) *                                           &
769                             u(k,j,i)
770
771!
772!--             Calculate preliminary new velocity, based on pre_tend
773                pre_u = u(k,j,i) + dt_3d * pre_tend
774!
775!--             Compare sign of old velocity and new preliminary velocity,
776!--             and in case the signs are different, limit the tendency
777                IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
778                   pre_tend = - u(k,j,i) * ddt_3d
779                ELSE
780                   pre_tend = pre_tend
781                ENDIF
782!
783!--             Calculate final tendency
784                tend(k,j,i) = tend(k,j,i) + pre_tend
785             ENDDO
786
787
788!
789!--       v-component
790          CASE ( 2 )
791             DO  k = nzb_v_inner(j,i)+1, nzb_v_inner(j,i)+pch_index
792
793                kk = k - nzb_v_inner(j,i)  !- lad arrays are defined flat
794!
795!--             In order to create sharp boundaries of the plant canopy,
796!--             the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
797!--             rather than being interpolated from the surrounding lad_s,
798!--             because this would yield smaller lad at the canopy boundaries
799!--             than inside of the canopy.
800!--             For the same reason, the lad at the northmost(j+1)canopy
801!--             boundary on the v-grid equals lad_s(k,j,i).
802                lad_local = lad_s(kk,j,i)
803                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp )  THEN
804                   lad_local = lad_s(kk,j-1,i)
805                ENDIF
806
807                pre_tend = 0.0_wp
808                pre_v = 0.0_wp
809!
810!--             Calculate preliminary value (pre_tend) of the tendency
811                pre_tend = - cdc *                                             &
812                             lad_local *                                       &
813                             SQRT( ( 0.25_wp * ( u(k,j-1,i)   +                &
814                                                 u(k,j-1,i+1) +                &
815                                                 u(k,j,i)     +                &
816                                                 u(k,j,i+1) )                  &
817                                   )**2 +                                      &
818                                   v(k,j,i)**2 +                               &
819                                   ( 0.25_wp * ( w(k-1,j-1,i) +                &
820                                                 w(k-1,j,i)   +                &
821                                                 w(k,j-1,i)   +                &
822                                                 w(k,j,i) )                    &
823                                   )**2                                        &
824                                 ) *                                           &
825                             v(k,j,i)
826
827!
828!--             Calculate preliminary new velocity, based on pre_tend
829                pre_v = v(k,j,i) + dt_3d * pre_tend
830!
831!--             Compare sign of old velocity and new preliminary velocity,
832!--             and in case the signs are different, limit the tendency
833                IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
834                   pre_tend = - v(k,j,i) * ddt_3d
835                ELSE
836                   pre_tend = pre_tend
837                ENDIF
838!
839!--             Calculate final tendency
840                tend(k,j,i) = tend(k,j,i) + pre_tend
841             ENDDO
842
843
844!
845!--       w-component
846          CASE ( 3 )
847             DO  k = nzb_w_inner(j,i)+1, nzb_w_inner(j,i)+pch_index-1
848
849                kk = k - nzb_w_inner(j,i)  !- lad arrays are defined flat
850
851                pre_tend = 0.0_wp
852                pre_w = 0.0_wp
853!
854!--             Calculate preliminary value (pre_tend) of the tendency
855                pre_tend = - cdc *                                             &
856                             (0.5_wp *                                         &
857                                ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) *         &
858                             SQRT( ( 0.25_wp * ( u(k,j,i)    +                 & 
859                                                 u(k,j,i+1)  +                 &
860                                                 u(k+1,j,i)  +                 &
861                                                 u(k+1,j,i+1) )                &
862                                   )**2 +                                      &
863                                   ( 0.25_wp * ( v(k,j,i)    +                 &
864                                                 v(k,j+1,i)  +                 &
865                                                 v(k+1,j,i)  +                 &
866                                                 v(k+1,j+1,i) )                &
867                                   )**2 +                                      &
868                                   w(k,j,i)**2                                 &
869                                 ) *                                           &
870                             w(k,j,i)
871!
872!--             Calculate preliminary new velocity, based on pre_tend
873                pre_w = w(k,j,i) + dt_3d * pre_tend
874!
875!--             Compare sign of old velocity and new preliminary velocity,
876!--             and in case the signs are different, limit the tendency
877                IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
878                   pre_tend = - w(k,j,i) * ddt_3d
879                ELSE
880                   pre_tend = pre_tend
881                ENDIF
882!
883!--             Calculate final tendency
884                tend(k,j,i) = tend(k,j,i) + pre_tend
885             ENDDO
886
887!
888!--       potential temperature
889          CASE ( 4 )
890             DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
891                kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
892                tend(k,j,i) = tend(k,j,i) +                                    &
893                                 ( canopy_heat_flux(kk,j,i) -                  &
894                                   canopy_heat_flux(kk-1,j,i) ) / dzw(k)
895             ENDDO
896
897
898!
899!--       scalar concentration
900          CASE ( 5 )
901             DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
902                kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
903                tend(k,j,i) = tend(k,j,i) -                                    &
904                                 lsec *                                        &
905                                 lad_s(kk,j,i) *                               &
906                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
907                                                    u(k,j,i+1) )               &
908                                       )**2  +                                 &
909                                       ( 0.5_wp * ( v(k,j,i) +                 &
910                                                    v(k,j+1,i) )               &
911                                       )**2 +                                  &
912                                       ( 0.5_wp * ( w(k-1,j,i) +               &
913                                                    w(k,j,i) )                 &
914                                       )**2                                    &
915                                     ) *                                       &
916                                 ( q(k,j,i) - lsc )
917             ENDDO   
918
919!
920!--       sgs-tke
921          CASE ( 6 )
922             DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
923                kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
924                tend(k,j,i) = tend(k,j,i) -                                    &
925                                 2.0_wp * cdc *                                &
926                                 lad_s(kk,j,i) *                               &
927                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
928                                                    u(k,j,i+1) )               &
929                                       )**2 +                                  & 
930                                       ( 0.5_wp * ( v(k,j,i) +                 &
931                                                    v(k,j+1,i) )               &
932                                       )**2 +                                  &
933                                       ( 0.5_wp * ( w(k,j,i) +                 &
934                                                    w(k+1,j,i) )               &
935                                       )**2                                    &
936                                     ) *                                       &
937                                 e(k,j,i)
938             ENDDO
939
940       CASE DEFAULT
941
942          WRITE( message_string, * ) 'wrong component: ', component
943          CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 ) 
944
945       END SELECT
946
947    END SUBROUTINE plant_canopy_model_ij
948
949 END MODULE plant_canopy_model_mod
Note: See TracBrowser for help on using the repository browser.