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

Last change on this file since 1682 was 1682, checked in by knoop, 6 years ago

Code annotations made doxygen readable

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