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

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

last commit documented

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