source: palm/tags/release-4.0/SOURCE/plant_canopy_model.f90 @ 3800

Last change on this file since 3800 was 1485, checked in by kanani, 9 years ago

last commit documented

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