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

Last change on this file since 2210 was 2210, checked in by kanani, 4 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 63.2 KB
Line 
1!> @file plant_canopy_model_mod.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
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
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-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: plant_canopy_model_mod.f90 2210 2017-04-19 09:35:50Z kanani $
27!
28! 2209 2017-04-19 09:34:46Z kanani
29! Added 3d output of leaf area density (pcm_lad) and canopy
30! heat rate (pcm_heatrate)
31!
32! 2024 2016-10-12 16:42:37Z kanani
33! Added missing lad_s initialization
34!
35! 2011 2016-09-19 17:29:57Z kanani
36! Renamed canopy_heat_flux to pc_heating_rate, since the original meaning/
37! calculation of the quantity has changed, related to the urban surface model
38! and similar future applications.
39!
40! 2007 2016-08-24 15:47:17Z kanani
41! Added SUBROUTINE pcm_read_plant_canopy_3d for reading 3d plant canopy data
42! from file (new case canopy_mode=read_from_file_3d) in the course of
43! introduction of urban surface model,
44! introduced variable ext_coef,
45! resorted SUBROUTINEs to alphabetical order
46!
47! 2000 2016-08-20 18:09:15Z knoop
48! Forced header and separation lines into 80 columns
49!
50! 1960 2016-07-12 16:34:24Z suehring
51! Separate humidity and passive scalar
52!
53! 1953 2016-06-21 09:28:42Z suehring
54! Bugfix, lad_s and lad must be public
55!
56! 1826 2016-04-07 12:01:39Z maronga
57! Further modularization
58!
59! 1721 2015-11-16 12:56:48Z raasch
60! bugfixes: shf is reduced in areas covered with canopy only,
61!           canopy is set on top of topography
62!
63! 1682 2015-10-07 23:56:08Z knoop
64! Code annotations made doxygen readable
65!
66! 1484 2014-10-21 10:53:05Z kanani
67! Changes due to new module structure of the plant canopy model:
68!   module plant_canopy_model_mod now contains a subroutine for the
69!   initialization of the canopy model (pcm_init),
70!   limitation of the canopy drag (previously accounted for by calculation of
71!   a limiting canopy timestep for the determination of the maximum LES timestep
72!   in subroutine timestep) is now realized by the calculation of pre-tendencies
73!   and preliminary velocities in subroutine pcm_tendency,
74!   some redundant MPI communication removed in subroutine pcm_init
75!   (was previously in init_3d_model),
76!   unnecessary 3d-arrays lad_u, lad_v, lad_w removed - lad information on the
77!   respective grid is now provided only by lad_s (e.g. in the calculation of
78!   the tendency terms or of cum_lai_hf),
79!   drag_coefficient, lai, leaf_surface_concentration,
80!   scalar_exchange_coefficient, sec and sls renamed to canopy_drag_coeff,
81!   cum_lai_hf, leaf_surface_conc, leaf_scalar_exch_coeff, lsec and lsc,
82!   respectively,
83!   unnecessary 3d-arrays cdc, lsc and lsec now defined as single-value constants,
84!   USE-statements and ONLY-lists modified accordingly
85!
86! 1340 2014-03-25 19:45:13Z kanani
87! REAL constants defined as wp-kind
88!
89! 1320 2014-03-20 08:40:49Z raasch
90! ONLY-attribute added to USE-statements,
91! kind-parameters added to all INTEGER and REAL declaration statements,
92! kinds are defined in new module kinds,
93! old module precision_kind is removed,
94! revision history before 2012 removed,
95! comment fields (!:) to be used for variable explanations added to
96! all variable declaration statements
97!
98! 1036 2012-10-22 13:43:42Z raasch
99! code put under GPL (PALM 3.9)
100!
101! 138 2007-11-28 10:03:58Z letzel
102! Initial revision
103!
104! Description:
105! ------------
106!> 1) Initialization of the canopy model, e.g. construction of leaf area density
107!> profile (subroutine pcm_init).
108!> 2) Calculation of sinks and sources of momentum, heat and scalar concentration
109!> due to canopy elements (subroutine pcm_tendency).
110!------------------------------------------------------------------------------!
111 MODULE plant_canopy_model_mod
112 
113    USE arrays_3d,                                                             &
114        ONLY:  dzu, dzw, e, q, s, shf, tend, u, v, w, zu, zw
115
116    USE indices,                                                               &
117        ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
118               nz, nzb, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt
119
120    USE kinds
121
122
123    IMPLICIT NONE
124
125
126    CHARACTER (LEN=20)   ::  canopy_mode = 'block' !< canopy coverage
127
128    INTEGER(iwp) ::  pch_index = 0                 !< plant canopy height/top index
129    INTEGER(iwp) ::                                                            &
130       lad_vertical_gradient_level_ind(10) = -9999 !< lad-profile levels (index)
131
132    LOGICAL ::  calc_beta_lad_profile = .FALSE. !< switch for calc. of lad from beta func.
133    LOGICAL ::  plant_canopy = .FALSE.          !< switch for use of canopy model
134
135    REAL(wp) ::  alpha_lad = 9999999.9_wp   !< coefficient for lad calculation
136    REAL(wp) ::  beta_lad = 9999999.9_wp    !< coefficient for lad calculation
137    REAL(wp) ::  canopy_drag_coeff = 0.0_wp !< canopy drag coefficient (parameter)
138    REAL(wp) ::  cdc = 0.0_wp               !< canopy drag coeff. (abbreviation used in equations)
139    REAL(wp) ::  cthf = 0.0_wp              !< canopy top heat flux
140    REAL(wp) ::  dt_plant_canopy = 0.0_wp   !< timestep account. for canopy drag
141    REAL(wp) ::  ext_coef = 0.6_wp          !< extinction coefficient
142    REAL(wp) ::  lad_surface = 0.0_wp       !< lad surface value
143    REAL(wp) ::  lai_beta = 0.0_wp          !< leaf area index (lai) for lad calc.
144    REAL(wp) ::                                                                &
145       leaf_scalar_exch_coeff = 0.0_wp      !< canopy scalar exchange coeff.
146    REAL(wp) ::                                                                &
147       leaf_surface_conc = 0.0_wp           !< leaf surface concentration
148    REAL(wp) ::  lsec = 0.0_wp              !< leaf scalar exchange coeff.
149    REAL(wp) ::  lsc = 0.0_wp               !< leaf surface concentration
150
151    REAL(wp) ::                                                                &
152       lad_vertical_gradient(10) = 0.0_wp              !< lad gradient
153    REAL(wp) ::                                                                &
154       lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m)
155
156    REAL(wp), DIMENSION(:), ALLOCATABLE ::  lad            !< leaf area density
157    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pre_lad        !< preliminary lad
158   
159    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
160       pc_heating_rate                                    !< plant canopy heating rate
161    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  cum_lai_hf !< cumulative lai for heatflux calc.
162    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  lad_s      !< lad on scalar-grid
163
164
165    SAVE
166
167
168    PRIVATE
169 
170!
171!-- Public functions
172    PUBLIC pcm_check_data_output, pcm_check_parameters, pcm_data_output_3d,    &
173           pcm_define_netcdf_grid, pcm_header, pcm_init, pcm_parin, pcm_tendency
174
175!
176!-- Public variables and constants
177    PUBLIC pc_heating_rate, canopy_mode, cthf, dt_plant_canopy, lad, lad_s,   &
178           pch_index, plant_canopy
179           
180
181    INTERFACE pcm_check_data_output
182       MODULE PROCEDURE pcm_check_data_output
183    END INTERFACE pcm_check_data_output
184   
185    INTERFACE pcm_check_parameters
186       MODULE PROCEDURE pcm_check_parameters
187    END INTERFACE pcm_check_parameters
188
189    INTERFACE pcm_data_output_3d
190       MODULE PROCEDURE pcm_data_output_3d
191    END INTERFACE pcm_data_output_3d
192
193    INTERFACE pcm_define_netcdf_grid
194       MODULE PROCEDURE pcm_define_netcdf_grid
195    END INTERFACE pcm_define_netcdf_grid
196   
197     INTERFACE pcm_header
198       MODULE PROCEDURE pcm_header
199    END INTERFACE pcm_header       
200   
201    INTERFACE pcm_init
202       MODULE PROCEDURE pcm_init
203    END INTERFACE pcm_init
204
205    INTERFACE pcm_parin
206       MODULE PROCEDURE pcm_parin
207    END INTERFACE pcm_parin
208
209    INTERFACE pcm_read_plant_canopy_3d
210       MODULE PROCEDURE pcm_read_plant_canopy_3d
211    END INTERFACE pcm_read_plant_canopy_3d
212   
213    INTERFACE pcm_tendency
214       MODULE PROCEDURE pcm_tendency
215       MODULE PROCEDURE pcm_tendency_ij
216    END INTERFACE pcm_tendency
217
218
219 CONTAINS
220
221
222!------------------------------------------------------------------------------!
223! Description:
224! ------------
225!> Check data output for plant canopy model
226!------------------------------------------------------------------------------!
227 SUBROUTINE pcm_check_data_output( var, unit )
228 
229 
230    USE control_parameters,                                                 &
231        ONLY:  data_output, message_string
232
233    IMPLICIT NONE
234
235    CHARACTER (LEN=*) ::  unit  !<
236    CHARACTER (LEN=*) ::  var   !<
237
238
239    SELECT CASE ( TRIM( var ) )
240
241       CASE ( 'pcm_heatrate' )
242          unit = 'K s-1'
243   
244       CASE ( 'pcm_lad' )
245          unit = 'm2 m-3'
246
247
248       CASE DEFAULT
249          unit = 'illegal'
250
251    END SELECT
252
253
254 END SUBROUTINE pcm_check_data_output
255 
256 
257!------------------------------------------------------------------------------!
258! Description:
259! ------------
260!> Check parameters routine for plant canopy model
261!------------------------------------------------------------------------------!
262    SUBROUTINE pcm_check_parameters
263
264       USE control_parameters,                                                 &
265           ONLY: cloud_physics, message_string, microphysics_seifert
266                 
267   
268       IMPLICIT NONE
269
270   
271       IF ( canopy_drag_coeff == 0.0_wp )  THEN
272          message_string = 'plant_canopy = .TRUE. requires a non-zero drag '// &
273                           'coefficient & given value is canopy_drag_coeff = 0.0'
274          CALL message( 'check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
275       ENDIF
276   
277       IF ( ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad == 9999999.9_wp )  .OR.&
278              beta_lad /= 9999999.9_wp   .AND.  alpha_lad == 9999999.9_wp )  THEN
279          message_string = 'using the beta function for the construction ' //  &
280                           'of the leaf area density profile requires '    //  &
281                           'both alpha_lad and beta_lad to be /= 9999999.9'
282          CALL message( 'check_parameters', 'PA0118', 1, 2, 0, 6, 0 )
283       ENDIF
284   
285       IF ( calc_beta_lad_profile  .AND.  lai_beta == 0.0_wp )  THEN
286          message_string = 'using the beta function for the construction ' //  &
287                           'of the leaf area density profile requires '    //  &
288                           'a non-zero lai_beta, but given value is '      //  &
289                           'lai_beta = 0.0'
290          CALL message( 'check_parameters', 'PA0119', 1, 2, 0, 6, 0 )
291       ENDIF
292
293       IF ( calc_beta_lad_profile  .AND.  lad_surface /= 0.0_wp )  THEN
294          message_string = 'simultaneous setting of alpha_lad /= 9999999.9' // &
295                           'and lad_surface /= 0.0 is not possible, '       // &
296                           'use either vertical gradients or the beta '     // &
297                           'function for the construction of the leaf area '// &
298                           'density profile'
299          CALL message( 'check_parameters', 'PA0120', 1, 2, 0, 6, 0 )
300       ENDIF
301
302       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
303          message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' // &
304                          ' seifert_beheng'
305          CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
306       ENDIF
307
308 
309    END SUBROUTINE pcm_check_parameters
310 
311
312!------------------------------------------------------------------------------!
313!
314! Description:
315! ------------
316!> Subroutine defining 3D output variables
317!------------------------------------------------------------------------------!
318 SUBROUTINE pcm_data_output_3d( av, variable, found, local_pf )
319 
320    USE control_parameters,                                                    &
321        ONLY :  nz_do3d
322 
323    USE indices
324
325    USE kinds
326
327
328    IMPLICIT NONE
329
330    CHARACTER (LEN=*) ::  variable !<
331
332    INTEGER(iwp) ::  av    !<
333    INTEGER(iwp) ::  i     !<
334    INTEGER(iwp) ::  j     !<
335    INTEGER(iwp) ::  k     !<
336
337    LOGICAL      ::  found !<
338
339    REAL(sp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nz_do3d) ::  local_pf !<
340
341
342    found = .TRUE.
343
344
345    SELECT CASE ( TRIM( variable ) )
346
347      CASE ( 'pcm_heatrate' )
348         IF ( av == 0 )  THEN
349            DO  i = nxlg, nxrg
350               DO  j = nysg, nyng
351                  DO  k = nzb_s_inner(j,i), nz_do3d
352                     local_pf(i,j,k) = pc_heating_rate(k,j,i)
353                  ENDDO
354               ENDDO
355            ENDDO
356         ENDIF
357   
358   
359      CASE ( 'pcm_lad' )
360
361         IF ( av == 0 )  THEN
362            DO  i = nxlg, nxrg
363               DO  j = nysg, nyng
364                  DO  k = nzb_s_inner(j,i), nz_do3d
365                     local_pf(i,j,k) = lad_s(k,j,i)
366                  ENDDO
367               ENDDO
368            ENDDO
369         ENDIF
370                 
371         
372       CASE DEFAULT
373          found = .FALSE.
374
375    END SELECT
376
377
378 END SUBROUTINE pcm_data_output_3d
379         
380!------------------------------------------------------------------------------!
381!
382! Description:
383! ------------
384!> Subroutine defining appropriate grid for netcdf variables.
385!> It is called from subroutine netcdf.
386!------------------------------------------------------------------------------!
387 SUBROUTINE pcm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
388   
389     IMPLICIT NONE
390
391     CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
392     LOGICAL, INTENT(OUT)           ::  found       !<
393     CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
394     CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
395     CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
396
397     found  = .TRUE.
398
399!
400!--  Check for the grid
401     SELECT CASE ( TRIM( var ) )
402
403        CASE ( 'pcm_heatrate', 'pcm_lad' )
404           grid_x = 'x'
405           grid_y = 'y'
406           grid_z = 'zu'
407
408        CASE DEFAULT
409           found  = .FALSE.
410           grid_x = 'none'
411           grid_y = 'none'
412           grid_z = 'none'
413     END SELECT
414
415 END SUBROUTINE pcm_define_netcdf_grid
416 
417 
418!------------------------------------------------------------------------------!
419! Description:
420! ------------
421!> Header output for plant canopy model
422!------------------------------------------------------------------------------!
423    SUBROUTINE pcm_header ( io )
424
425       USE control_parameters,                                                 &
426           ONLY: dz, passive_scalar
427
428
429       IMPLICIT NONE
430 
431       CHARACTER (LEN=10) ::  coor_chr            !<
432
433       CHARACTER (LEN=86) ::  coordinates         !<
434       CHARACTER (LEN=86) ::  gradients           !<
435       CHARACTER (LEN=86) ::  leaf_area_density   !<
436       CHARACTER (LEN=86) ::  slices              !<
437 
438       INTEGER(iwp) :: i                !<
439       INTEGER(iwp),  INTENT(IN) ::  io !< Unit of the output file
440       INTEGER(iwp) :: k                !<       
441   
442       REAL(wp) ::  canopy_height       !< canopy height (in m)
443       
444       canopy_height = pch_index * dz
445
446       WRITE ( io, 1 )  canopy_mode, canopy_height, pch_index,                 &
447                          canopy_drag_coeff
448       IF ( passive_scalar )  THEN
449          WRITE ( io, 2 )  leaf_scalar_exch_coeff,                             &
450                             leaf_surface_conc
451       ENDIF
452
453!
454!--    Heat flux at the top of vegetation
455       WRITE ( io, 3 )  cthf
456
457!
458!--    Leaf area density profile, calculated either from given vertical
459!--    gradients or from beta probability density function.
460       IF (  .NOT.  calc_beta_lad_profile )  THEN
461
462!--       Building output strings, starting with surface value
463          WRITE ( leaf_area_density, '(F7.4)' )  lad_surface
464          gradients = '------'
465          slices = '     0'
466          coordinates = '   0.0'
467          i = 1
468          DO  WHILE ( i < 11  .AND.  lad_vertical_gradient_level_ind(i)        &
469                      /= -9999 )
470
471             WRITE (coor_chr,'(F7.2)')  lad(lad_vertical_gradient_level_ind(i))
472             leaf_area_density = TRIM( leaf_area_density ) // ' ' //           &
473                                 TRIM( coor_chr )
474 
475             WRITE (coor_chr,'(F7.2)')  lad_vertical_gradient(i)
476             gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
477
478             WRITE (coor_chr,'(I7)')  lad_vertical_gradient_level_ind(i)
479             slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
480
481             WRITE (coor_chr,'(F7.1)')  lad_vertical_gradient_level(i)
482             coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
483
484             i = i + 1
485          ENDDO
486
487          WRITE ( io, 4 )  TRIM( coordinates ), TRIM( leaf_area_density ),     &
488                             TRIM( gradients ), TRIM( slices )
489
490       ELSE
491       
492          WRITE ( leaf_area_density, '(F7.4)' )  lad_surface
493          coordinates = '   0.0'
494         
495          DO  k = 1, pch_index
496
497             WRITE (coor_chr,'(F7.2)')  lad(k)
498             leaf_area_density = TRIM( leaf_area_density ) // ' ' //           &
499                                 TRIM( coor_chr )
500 
501             WRITE (coor_chr,'(F7.1)')  zu(k)
502             coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
503
504          ENDDO       
505
506          WRITE ( io, 5 ) TRIM( coordinates ), TRIM( leaf_area_density ),      &
507                          alpha_lad, beta_lad, lai_beta
508
509       ENDIF 
510
5111 FORMAT (//' Vegetation canopy (drag) model:'/                                &
512              ' ------------------------------'//                              &
513              ' Canopy mode: ', A /                                            &
514              ' Canopy height: ',F6.2,'m (',I4,' grid points)' /               &
515              ' Leaf drag coefficient: ',F6.2 /)
5162 FORMAT (/ ' Scalar exchange coefficient: ',F6.2 /                            &
517              ' Scalar concentration at leaf surfaces in kg/m**3: ',F6.2 /)
5183 FORMAT (' Predefined constant heatflux at the top of the vegetation: ',F6.2, &
519          ' K m/s')
5204 FORMAT (/ ' Characteristic levels of the leaf area density:'//               &
521              ' Height:              ',A,'  m'/                                &
522              ' Leaf area density:   ',A,'  m**2/m**3'/                        &
523              ' Gradient:            ',A,'  m**2/m**4'/                        &
524              ' Gridpoint:           ',A)
5255 FORMAT (//' Characteristic levels of the leaf area density and coefficients:'&
526          //  ' Height:              ',A,'  m'/                                &
527              ' Leaf area density:   ',A,'  m**2/m**3'/                        &
528              ' Coefficient alpha: ',F6.2 /                                    &
529              ' Coefficient beta: ',F6.2 /                                     &
530              ' Leaf area index: ',F6.2,'  m**2/m**2' /)   
531       
532    END SUBROUTINE pcm_header
533 
534 
535!------------------------------------------------------------------------------!
536! Description:
537! ------------
538!> Initialization of the plant canopy model
539!------------------------------------------------------------------------------!
540    SUBROUTINE pcm_init
541   
542
543       USE control_parameters,                                                 &
544           ONLY: coupling_char, dz, humidity, io_blocks, io_group,             &
545                 message_string, ocean, passive_scalar
546
547
548       IMPLICIT NONE
549
550       CHARACTER(10) :: pct
551       
552       INTEGER(iwp) ::  i   !< running index
553       INTEGER(iwp) ::  ii  !< index       
554       INTEGER(iwp) ::  j   !< running index
555       INTEGER(iwp) ::  k   !< running index
556
557       REAL(wp) ::  int_bpdf        !< vertical integral for lad-profile construction
558       REAL(wp) ::  dzh             !< vertical grid spacing in units of canopy height
559       REAL(wp) ::  gradient        !< gradient for lad-profile construction
560       REAL(wp) ::  canopy_height   !< canopy height for lad-profile construction
561       REAL(wp) ::  pcv(nzb:nzt+1)  !<
562       
563!
564!--    Allocate one-dimensional arrays for the computation of the
565!--    leaf area density (lad) profile
566       ALLOCATE( lad(0:nz+1), pre_lad(0:nz+1) )
567       lad = 0.0_wp
568       pre_lad = 0.0_wp
569
570!
571!--    Set flag that indicates that the lad-profile shall be calculated by using
572!--    a beta probability density function
573       IF ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad /= 9999999.9_wp )  THEN
574          calc_beta_lad_profile = .TRUE.
575       ENDIF
576       
577       
578!
579!--    Compute the profile of leaf area density used in the plant
580!--    canopy model. The profile can either be constructed from
581!--    prescribed vertical gradients of the leaf area density or by
582!--    using a beta probability density function (see e.g. Markkanen et al.,
583!--    2003: Boundary-Layer Meteorology, 106, 437-459)
584       IF (  .NOT.  calc_beta_lad_profile )  THEN   
585
586!
587!--       Use vertical gradients for lad-profile construction   
588          i = 1
589          gradient = 0.0_wp
590
591          IF (  .NOT.  ocean )  THEN
592
593             lad(0) = lad_surface
594             lad_vertical_gradient_level_ind(1) = 0
595 
596             DO k = 1, pch_index
597                IF ( i < 11 )  THEN
598                   IF ( lad_vertical_gradient_level(i) < zu(k)  .AND.          &
599                        lad_vertical_gradient_level(i) >= 0.0_wp )  THEN
600                      gradient = lad_vertical_gradient(i)
601                      lad_vertical_gradient_level_ind(i) = k - 1
602                      i = i + 1
603                   ENDIF
604                ENDIF
605                IF ( gradient /= 0.0_wp )  THEN
606                   IF ( k /= 1 )  THEN
607                      lad(k) = lad(k-1) + dzu(k) * gradient
608                   ELSE
609                      lad(k) = lad_surface + dzu(k) * gradient
610                   ENDIF
611                ELSE
612                   lad(k) = lad(k-1)
613                ENDIF
614             ENDDO
615
616          ENDIF
617
618!
619!--       In case of no given leaf area density gradients, choose a vanishing
620!--       gradient. This information is used for the HEADER and the RUN_CONTROL
621!--       file.
622          IF ( lad_vertical_gradient_level(1) == -9999999.9_wp )  THEN
623             lad_vertical_gradient_level(1) = 0.0_wp
624          ENDIF
625
626       ELSE
627
628!
629!--       Use beta function for lad-profile construction
630          int_bpdf = 0.0_wp
631          canopy_height = pch_index * dz
632
633          DO k = nzb, pch_index
634             int_bpdf = int_bpdf +                                             &
635                      ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) *  &
636                      ( ( 1.0_wp - ( zw(k) / canopy_height ) )**(              &
637                          beta_lad-1.0_wp ) )                                  &
638                      * ( ( zw(k+1)-zw(k) ) / canopy_height ) )
639          ENDDO
640
641!
642!--       Preliminary lad profile (defined on w-grid)
643          DO k = nzb, pch_index
644             pre_lad(k) =  lai_beta *                                          &
645                        ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) )  &
646                        * ( ( 1.0_wp - ( zw(k) / canopy_height ) )**(          &
647                              beta_lad-1.0_wp ) ) / int_bpdf                   &
648                        ) / canopy_height
649          ENDDO
650
651!
652!--       Final lad profile (defined on scalar-grid level, since most prognostic
653!--       quantities are defined there, hence, less interpolation is required
654!--       when calculating the canopy tendencies)
655          lad(0) = pre_lad(0)
656          DO k = nzb+1, pch_index
657             lad(k) = 0.5 * ( pre_lad(k-1) + pre_lad(k) )
658          ENDDO         
659
660       ENDIF
661
662!
663!--    Allocate 3D-array for the leaf area density (lad_s). In case of a
664!--    prescribed canopy-top heat flux (cthf), allocate 3D-arrays for
665!--    the cumulative leaf area index (cum_lai_hf) and the canopy heat flux.
666       ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
667
668       IF ( cthf /= 0.0_wp )  THEN
669          ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                 &
670                    pc_heating_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
671       ENDIF
672
673!
674!--    Initialize canopy parameters cdc (canopy drag coefficient),
675!--    lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration)
676!--    with the prescribed values
677       cdc = canopy_drag_coeff
678       lsec = leaf_scalar_exch_coeff
679       lsc = leaf_surface_conc
680
681!
682!--    Initialization of the canopy coverage in the model domain:
683!--    Setting the parameter canopy_mode = 'block' initializes a canopy, which
684!--    fully covers the domain surface
685       SELECT CASE ( TRIM( canopy_mode ) )
686
687          CASE( 'block' )
688
689             DO  i = nxlg, nxrg
690                DO  j = nysg, nyng
691                   lad_s(:,j,i) = lad(:)
692                ENDDO
693             ENDDO
694
695          CASE ( 'read_from_file_3d' )
696!
697!--          Initialize canopy parameters cdc (canopy drag coefficient),
698!--          lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration)
699!--          from file which contains complete 3D data (separate vertical profiles for
700!--          each location).
701             CALL pcm_read_plant_canopy_3d
702
703          CASE DEFAULT
704!
705!--          The DEFAULT case is reached either if the parameter
706!--          canopy mode contains a wrong character string or if the
707!--          user has coded a special case in the user interface.
708!--          There, the subroutine user_init_plant_canopy checks
709!--          which of these two conditions applies.
710             CALL user_init_plant_canopy
711 
712       END SELECT
713
714!
715!--    Initialization of the canopy heat source distribution due to heating
716!--    of the canopy layers by incoming solar radiation, in case that a non-zero
717!--    value is set for the canopy top heat flux (cthf), which equals the
718!--    available net radiation at canopy top.
719!--    The heat source distribution is calculated by a decaying exponential
720!--    function of the downward cumulative leaf area index (cum_lai_hf),
721!--    assuming that the foliage inside the plant canopy is heated by solar
722!--    radiation penetrating the canopy layers according to the distribution
723!--    of net radiation as suggested by Brown & Covey (1966; Agric. Meteorol. 3,
724!--    73–96). This approach has been applied e.g. by Shaw & Schumann (1992;
725!--    Bound.-Layer Meteorol. 61, 47–64)
726       IF ( cthf /= 0.0_wp )  THEN
727!
728!--       Piecewise calculation of the cumulative leaf area index by vertical
729!--       integration of the leaf area density
730          cum_lai_hf(:,:,:) = 0.0_wp
731          DO  i = nxlg, nxrg
732             DO  j = nysg, nyng
733                DO  k = pch_index-1, 0, -1
734                   IF ( k == pch_index-1 )  THEN
735                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
736                         ( 0.5_wp * lad_s(k+1,j,i) *                           &
737                           ( zw(k+1) - zu(k+1) ) )  +                          &
738                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
739                                                 lad_s(k,j,i) ) +              &
740                                      lad_s(k+1,j,i) ) *                       &
741                           ( zu(k+1) - zw(k) ) ) 
742                   ELSE
743                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
744                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) +              &
745                                                 lad_s(k+1,j,i) ) +            &
746                                      lad_s(k+1,j,i) ) *                       &
747                           ( zw(k+1) - zu(k+1) ) )  +                          &
748                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
749                                                 lad_s(k,j,i) ) +              &
750                                      lad_s(k+1,j,i) ) *                       &
751                           ( zu(k+1) - zw(k) ) )
752                   ENDIF
753                ENDDO
754             ENDDO
755          ENDDO
756
757!
758!--       Calculation of the heating rate (K/s) within the different layers of
759!--       the plant canopy
760          DO  i = nxlg, nxrg
761             DO  j = nysg, nyng
762!
763!--             Calculation only necessary in areas covered with canopy
764                IF ( cum_lai_hf(0,j,i) /= 0.0_wp )  THEN
765!--             
766!--                In areas with canopy the surface value of the canopy heat
767!--                flux distribution overrides the surface heat flux (shf)
768                   shf(j,i) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
769!
770!--                Within the different canopy layers the plant-canopy heating
771!--                rate (pc_heating_rate) is calculated as the vertical
772!--                divergence of the canopy heat fluxes at the top and bottom
773!--                of the respective layer
774                   DO  k = 1, pch_index
775                      pc_heating_rate(k,j,i) = cthf *                         &
776                                                ( exp(-ext_coef*cum_lai_hf(k,j,i)) -  &
777                                                  exp(-ext_coef*cum_lai_hf(k-1,j,i)) ) / dzw(k)
778                   ENDDO
779                ENDIF
780             ENDDO
781          ENDDO
782
783       ENDIF
784
785
786
787    END SUBROUTINE pcm_init
788
789
790!------------------------------------------------------------------------------!
791! Description:
792! ------------
793!> Parin for &canopy_par for plant canopy model
794!------------------------------------------------------------------------------!
795    SUBROUTINE pcm_parin
796
797
798       IMPLICIT NONE
799
800       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
801       
802       NAMELIST /canopy_par/      alpha_lad, beta_lad, canopy_drag_coeff,      &
803                                  canopy_mode, cthf,                           &
804                                  lad_surface,                                 & 
805                                  lad_vertical_gradient,                       &
806                                  lad_vertical_gradient_level,                 &
807                                  lai_beta,                                    &
808                                  leaf_scalar_exch_coeff,                      &
809                                  leaf_surface_conc, pch_index
810       
811       line = ' '
812       
813!
814!--    Try to find radiation model package
815       REWIND ( 11 )
816       line = ' '
817       DO   WHILE ( INDEX( line, '&canopy_par' ) == 0 )
818          READ ( 11, '(A)', END=10 )  line
819       ENDDO
820       BACKSPACE ( 11 )
821
822!
823!--    Read user-defined namelist
824       READ ( 11, canopy_par )
825
826!
827!--    Set flag that indicates that the radiation model is switched on
828       plant_canopy = .TRUE.
829
830 10    CONTINUE
831       
832
833    END SUBROUTINE pcm_parin
834
835
836
837!------------------------------------------------------------------------------!
838! Description:
839! ------------
840!
841!> Loads 3D plant canopy data from file. File format is as follows:
842!>
843!> num_levels
844!> dtype,x,y,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
845!> dtype,x,y,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
846!> dtype,x,y,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
847!> ...
848!>
849!> i.e. first line determines number of levels and further lines represent plant
850!> canopy data, one line per column and variable. In each data line,
851!> dtype represents variable to be set:
852!>
853!> dtype=1: leaf area density (lad_s)
854!> dtype=2: canopy drag coefficient (cdc)
855!> dtype=3: leaf scalar exchange coefficient (lsec)
856!> dtype=4: leaf surface concentration (lsc)
857!>
858!> Zeros are added automatically above num_levels until top of domain.  Any
859!> non-specified (x,y) columns have zero values as default.
860!------------------------------------------------------------------------------!
861    SUBROUTINE pcm_read_plant_canopy_3d
862        USE control_parameters, &
863            ONLY: passive_scalar, message_string
864        IMPLICIT NONE
865
866        INTEGER(iwp)                            :: i, j, dtype, nzp, nzpltop, nzpl, kk
867        REAL(wp), DIMENSION(:), ALLOCATABLE     :: col
868
869        lad_s = 0.0_wp
870        OPEN(152, file='PLANT_CANOPY_DATA_3D', access='SEQUENTIAL', &
871                action='READ', status='OLD', form='FORMATTED', err=515)
872        READ(152, *, err=516, end=517) nzp   !< read first line = number of vertical layers
873        ALLOCATE(col(1:nzp))
874        nzpltop = MIN(nzt+1, nzb+nzp-1)
875        nzpl = nzpltop - nzb + 1    !< no. of layers to assign
876
877        DO
878            READ(152, *, err=516, end=517) dtype, i, j, col(:)
879            IF ( i < nxlg .or. i > nxrg .or. j < nysg .or. j > nyng ) CYCLE
880
881            SELECT CASE (dtype)
882              CASE( 1 ) !< leaf area density
883                !-- only lad_s has flat z-coordinate, others have regular
884                kk = nzb_s_inner(j, i)
885                lad_s(nzb:nzpltop-kk, j, i) = col(1+kk:nzpl)
886!               CASE( 2 ) !< canopy drag coefficient
887!                 cdc(nzb:nzpltop, j, i) = col(1:nzpl)
888!               CASE( 3 ) !< leaf scalar exchange coefficient
889!                 lsec(nzb:nzpltop, j, i) = col(1:nzpl)
890!               CASE( 4 ) !< leaf surface concentration
891!                 lsc(nzb:nzpltop, j, i) = col(1:nzpl)
892              CASE DEFAULT
893                write(message_string, '(a,i2,a)')   &
894                    'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"'
895                CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 )
896            END SELECT
897        ENDDO
898
899515     message_string = 'error opening file PLANT_CANOPY_DATA_3D'
900        CALL message( 'pcm_read_plant_canopy_3d', 'PA0531', 1, 2, 0, 6, 0 )
901
902516     message_string = 'error reading file PLANT_CANOPY_DATA_3D'
903        CALL message( 'pcm_read_plant_canopy_3d', 'PA0532', 1, 2, 0, 6, 0 )
904
905517     CLOSE(152)
906        DEALLOCATE(col)
907       
908    END SUBROUTINE pcm_read_plant_canopy_3d
909   
910   
911
912!------------------------------------------------------------------------------!
913! Description:
914! ------------
915!> Calculation of the tendency terms, accounting for the effect of the plant
916!> canopy on momentum and scalar quantities.
917!>
918!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
919!> (defined on scalar grid), as initialized in subroutine pcm_init.
920!> The lad on the w-grid is vertically interpolated from the surrounding
921!> lad_s. The upper boundary of the canopy is defined on the w-grid at
922!> k = pch_index. Here, the lad is zero.
923!>
924!> The canopy drag must be limited (previously accounted for by calculation of
925!> a limiting canopy timestep for the determination of the maximum LES timestep
926!> in subroutine timestep), since it is physically impossible that the canopy
927!> drag alone can locally change the sign of a velocity component. This
928!> limitation is realized by calculating preliminary tendencies and velocities.
929!> It is subsequently checked if the preliminary new velocity has a different
930!> sign than the current velocity. If so, the tendency is limited in a way that
931!> the velocity can at maximum be reduced to zero by the canopy drag.
932!>
933!>
934!> Call for all grid points
935!------------------------------------------------------------------------------!
936    SUBROUTINE pcm_tendency( component )
937
938
939       USE control_parameters,                                                 &
940           ONLY:  dt_3d, message_string
941
942       USE kinds
943
944       IMPLICIT NONE
945
946       INTEGER(iwp) ::  component !< prognostic variable (u,v,w,pt,q,e)
947       INTEGER(iwp) ::  i         !< running index
948       INTEGER(iwp) ::  j         !< running index
949       INTEGER(iwp) ::  k         !< running index
950       INTEGER(iwp) ::  kk        !< running index for flat lad arrays
951
952       REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
953       REAL(wp) ::  lad_local !< local lad value
954       REAL(wp) ::  pre_tend  !< preliminary tendency
955       REAL(wp) ::  pre_u     !< preliminary u-value
956       REAL(wp) ::  pre_v     !< preliminary v-value
957       REAL(wp) ::  pre_w     !< preliminary w-value
958
959
960       ddt_3d = 1.0_wp / dt_3d
961 
962!
963!--    Compute drag for the three velocity components and the SGS-TKE:
964       SELECT CASE ( component )
965
966!
967!--       u-component
968          CASE ( 1 )
969             DO  i = nxlu, nxr
970                DO  j = nys, nyn
971                   DO  k = nzb_u_inner(j,i)+1, nzb_u_inner(j,i)+pch_index
972
973                      kk = k - nzb_u_inner(j,i)  !- lad arrays are defined flat
974!
975!--                   In order to create sharp boundaries of the plant canopy,
976!--                   the lad on the u-grid at index (k,j,i) is equal to
977!--                   lad_s(k,j,i), rather than being interpolated from the
978!--                   surrounding lad_s, because this would yield smaller lad
979!--                   at the canopy boundaries than inside of the canopy.
980!--                   For the same reason, the lad at the rightmost(i+1)canopy
981!--                   boundary on the u-grid equals lad_s(k,j,i).
982                      lad_local = lad_s(kk,j,i)
983                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp )&
984                      THEN
985                         lad_local = lad_s(kk,j,i-1)
986                      ENDIF
987
988                      pre_tend = 0.0_wp
989                      pre_u = 0.0_wp
990!
991!--                   Calculate preliminary value (pre_tend) of the tendency
992                      pre_tend = - cdc *                                       &
993                                   lad_local *                                 &
994                                   SQRT( u(k,j,i)**2 +                         &
995                                         ( 0.25_wp * ( v(k,j,i-1) +            &
996                                                       v(k,j,i)   +            &
997                                                       v(k,j+1,i) +            &
998                                                       v(k,j+1,i-1) )          &
999                                         )**2 +                                &
1000                                         ( 0.25_wp * ( w(k-1,j,i-1) +          &
1001                                                       w(k-1,j,i)   +          &
1002                                                       w(k,j,i-1)   +          &
1003                                                       w(k,j,i) )              &
1004                                         )**2                                  &
1005                                       ) *                                     &
1006                                   u(k,j,i)
1007
1008!
1009!--                   Calculate preliminary new velocity, based on pre_tend
1010                      pre_u = u(k,j,i) + dt_3d * pre_tend
1011!
1012!--                   Compare sign of old velocity and new preliminary velocity,
1013!--                   and in case the signs are different, limit the tendency
1014                      IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
1015                         pre_tend = - u(k,j,i) * ddt_3d
1016                      ELSE
1017                         pre_tend = pre_tend
1018                      ENDIF
1019!
1020!--                   Calculate final tendency
1021                      tend(k,j,i) = tend(k,j,i) + pre_tend
1022
1023                   ENDDO
1024                ENDDO
1025             ENDDO
1026
1027!
1028!--       v-component
1029          CASE ( 2 )
1030             DO  i = nxl, nxr
1031                DO  j = nysv, nyn
1032                   DO  k = nzb_v_inner(j,i)+1, nzb_v_inner(j,i)+pch_index
1033
1034                      kk = k - nzb_v_inner(j,i)  !- lad arrays are defined flat
1035!
1036!--                   In order to create sharp boundaries of the plant canopy,
1037!--                   the lad on the v-grid at index (k,j,i) is equal to
1038!--                   lad_s(k,j,i), rather than being interpolated from the
1039!--                   surrounding lad_s, because this would yield smaller lad
1040!--                   at the canopy boundaries than inside of the canopy.
1041!--                   For the same reason, the lad at the northmost(j+1) canopy
1042!--                   boundary on the v-grid equals lad_s(k,j,i).
1043                      lad_local = lad_s(kk,j,i)
1044                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp )&
1045                      THEN
1046                         lad_local = lad_s(kk,j-1,i)
1047                      ENDIF
1048
1049                      pre_tend = 0.0_wp
1050                      pre_v = 0.0_wp
1051!
1052!--                   Calculate preliminary value (pre_tend) of the tendency
1053                      pre_tend = - cdc *                                       &
1054                                   lad_local *                                 &
1055                                   SQRT( ( 0.25_wp * ( u(k,j-1,i)   +          &
1056                                                       u(k,j-1,i+1) +          &
1057                                                       u(k,j,i)     +          &
1058                                                       u(k,j,i+1) )            &
1059                                         )**2 +                                &
1060                                         v(k,j,i)**2 +                         &
1061                                         ( 0.25_wp * ( w(k-1,j-1,i) +          &
1062                                                       w(k-1,j,i)   +          &
1063                                                       w(k,j-1,i)   +          &
1064                                                       w(k,j,i) )              &
1065                                         )**2                                  &
1066                                       ) *                                     &
1067                                   v(k,j,i)
1068
1069!
1070!--                   Calculate preliminary new velocity, based on pre_tend
1071                      pre_v = v(k,j,i) + dt_3d * pre_tend
1072!
1073!--                   Compare sign of old velocity and new preliminary velocity,
1074!--                   and in case the signs are different, limit the tendency
1075                      IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
1076                         pre_tend = - v(k,j,i) * ddt_3d
1077                      ELSE
1078                         pre_tend = pre_tend
1079                      ENDIF
1080!
1081!--                   Calculate final tendency
1082                      tend(k,j,i) = tend(k,j,i) + pre_tend
1083
1084                   ENDDO
1085                ENDDO
1086             ENDDO
1087
1088!
1089!--       w-component
1090          CASE ( 3 )
1091             DO  i = nxl, nxr
1092                DO  j = nys, nyn
1093                   DO  k = nzb_w_inner(j,i)+1, nzb_w_inner(j,i)+pch_index-1
1094
1095                      kk = k - nzb_w_inner(j,i)  !- lad arrays are defined flat
1096
1097                      pre_tend = 0.0_wp
1098                      pre_w = 0.0_wp
1099!
1100!--                   Calculate preliminary value (pre_tend) of the tendency
1101                      pre_tend = - cdc *                                       &
1102                                   (0.5_wp *                                   &
1103                                      ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) *   &
1104                                   SQRT( ( 0.25_wp * ( u(k,j,i)   +            &
1105                                                       u(k,j,i+1) +            &
1106                                                       u(k+1,j,i) +            &
1107                                                       u(k+1,j,i+1) )          &
1108                                         )**2 +                                &
1109                                         ( 0.25_wp * ( v(k,j,i)   +            &
1110                                                       v(k,j+1,i) +            &
1111                                                       v(k+1,j,i) +            &
1112                                                       v(k+1,j+1,i) )          &
1113                                         )**2 +                                &
1114                                         w(k,j,i)**2                           &
1115                                       ) *                                     &
1116                                   w(k,j,i)
1117!
1118!--                   Calculate preliminary new velocity, based on pre_tend
1119                      pre_w = w(k,j,i) + dt_3d * pre_tend
1120!
1121!--                   Compare sign of old velocity and new preliminary velocity,
1122!--                   and in case the signs are different, limit the tendency
1123                      IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
1124                         pre_tend = - w(k,j,i) * ddt_3d
1125                      ELSE
1126                         pre_tend = pre_tend
1127                      ENDIF
1128!
1129!--                   Calculate final tendency
1130                      tend(k,j,i) = tend(k,j,i) + pre_tend
1131
1132                   ENDDO
1133                ENDDO
1134             ENDDO
1135
1136!
1137!--       potential temperature
1138          CASE ( 4 )
1139             DO  i = nxl, nxr
1140                DO  j = nys, nyn
1141                   DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
1142                      kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
1143                      tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
1144                   ENDDO
1145                ENDDO
1146             ENDDO
1147
1148!
1149!--       humidity
1150          CASE ( 5 )
1151             DO  i = nxl, nxr
1152                DO  j = nys, nyn
1153                   DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
1154                      kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
1155                      tend(k,j,i) = tend(k,j,i) -                              &
1156                                       lsec *                                  &
1157                                       lad_s(kk,j,i) *                         &
1158                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
1159                                                          u(k,j,i+1) )         &
1160                                             )**2 +                            &
1161                                             ( 0.5_wp * ( v(k,j,i) +           &
1162                                                          v(k,j+1,i) )         &
1163                                             )**2 +                            &
1164                                             ( 0.5_wp * ( w(k-1,j,i) +         & 
1165                                                          w(k,j,i) )           &
1166                                             )**2                              &
1167                                           ) *                                 &
1168                                       ( q(k,j,i) - lsc )
1169                   ENDDO
1170                ENDDO
1171             ENDDO
1172
1173!
1174!--       sgs-tke
1175          CASE ( 6 )
1176             DO  i = nxl, nxr
1177                DO  j = nys, nyn
1178                   DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
1179                      kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
1180                      tend(k,j,i) = tend(k,j,i) -                              &
1181                                       2.0_wp * cdc *                          &
1182                                       lad_s(kk,j,i) *                         &
1183                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
1184                                                          u(k,j,i+1) )         &
1185                                             )**2 +                            &
1186                                             ( 0.5_wp * ( v(k,j,i) +           &
1187                                                          v(k,j+1,i) )         &
1188                                             )**2 +                            &
1189                                             ( 0.5_wp * ( w(k,j,i) +           &
1190                                                          w(k+1,j,i) )         &
1191                                             )**2                              &
1192                                           ) *                                 &
1193                                       e(k,j,i)
1194                   ENDDO
1195                ENDDO
1196             ENDDO
1197!
1198!--       scalar concentration
1199          CASE ( 7 )
1200             DO  i = nxl, nxr
1201                DO  j = nys, nyn
1202                   DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
1203                      kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
1204                      tend(k,j,i) = tend(k,j,i) -                              &
1205                                       lsec *                                  &
1206                                       lad_s(kk,j,i) *                         &
1207                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
1208                                                          u(k,j,i+1) )         &
1209                                             )**2 +                            &
1210                                             ( 0.5_wp * ( v(k,j,i) +           &
1211                                                          v(k,j+1,i) )         &
1212                                             )**2 +                            &
1213                                             ( 0.5_wp * ( w(k-1,j,i) +         & 
1214                                                          w(k,j,i) )           &
1215                                             )**2                              &
1216                                           ) *                                 &
1217                                       ( s(k,j,i) - lsc )
1218                   ENDDO
1219                ENDDO
1220             ENDDO   
1221
1222
1223
1224          CASE DEFAULT
1225
1226             WRITE( message_string, * ) 'wrong component: ', component
1227             CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 ) 
1228
1229       END SELECT
1230
1231    END SUBROUTINE pcm_tendency
1232
1233
1234!------------------------------------------------------------------------------!
1235! Description:
1236! ------------
1237!> Calculation of the tendency terms, accounting for the effect of the plant
1238!> canopy on momentum and scalar quantities.
1239!>
1240!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
1241!> (defined on scalar grid), as initialized in subroutine pcm_init.
1242!> The lad on the w-grid is vertically interpolated from the surrounding
1243!> lad_s. The upper boundary of the canopy is defined on the w-grid at
1244!> k = pch_index. Here, the lad is zero.
1245!>
1246!> The canopy drag must be limited (previously accounted for by calculation of
1247!> a limiting canopy timestep for the determination of the maximum LES timestep
1248!> in subroutine timestep), since it is physically impossible that the canopy
1249!> drag alone can locally change the sign of a velocity component. This
1250!> limitation is realized by calculating preliminary tendencies and velocities.
1251!> It is subsequently checked if the preliminary new velocity has a different
1252!> sign than the current velocity. If so, the tendency is limited in a way that
1253!> the velocity can at maximum be reduced to zero by the canopy drag.
1254!>
1255!>
1256!> Call for grid point i,j
1257!------------------------------------------------------------------------------!
1258    SUBROUTINE pcm_tendency_ij( i, j, component )
1259
1260
1261       USE control_parameters,                                                 &
1262           ONLY:  dt_3d, message_string
1263
1264       USE kinds
1265
1266       IMPLICIT NONE
1267
1268       INTEGER(iwp) ::  component !< prognostic variable (u,v,w,pt,q,e)
1269       INTEGER(iwp) ::  i         !< running index
1270       INTEGER(iwp) ::  j         !< running index
1271       INTEGER(iwp) ::  k         !< running index
1272       INTEGER(iwp) ::  kk        !< running index for flat lad arrays
1273
1274       REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
1275       REAL(wp) ::  lad_local !< local lad value
1276       REAL(wp) ::  pre_tend  !< preliminary tendency
1277       REAL(wp) ::  pre_u     !< preliminary u-value
1278       REAL(wp) ::  pre_v     !< preliminary v-value
1279       REAL(wp) ::  pre_w     !< preliminary w-value
1280
1281
1282       ddt_3d = 1.0_wp / dt_3d
1283
1284!
1285!--    Compute drag for the three velocity components and the SGS-TKE
1286       SELECT CASE ( component )
1287
1288!
1289!--       u-component
1290          CASE ( 1 )
1291             DO  k = nzb_u_inner(j,i)+1, nzb_u_inner(j,i)+pch_index
1292
1293                kk = k - nzb_u_inner(j,i)  !- lad arrays are defined flat
1294!
1295!--             In order to create sharp boundaries of the plant canopy,
1296!--             the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
1297!--             rather than being interpolated from the surrounding lad_s,
1298!--             because this would yield smaller lad at the canopy boundaries
1299!--             than inside of the canopy.
1300!--             For the same reason, the lad at the rightmost(i+1)canopy
1301!--             boundary on the u-grid equals lad_s(k,j,i).
1302                lad_local = lad_s(kk,j,i)
1303                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp )  THEN
1304                   lad_local = lad_s(kk,j,i-1)
1305                ENDIF
1306
1307                pre_tend = 0.0_wp
1308                pre_u = 0.0_wp
1309!
1310!--             Calculate preliminary value (pre_tend) of the tendency
1311                pre_tend = - cdc *                                             &
1312                             lad_local *                                       &   
1313                             SQRT( u(k,j,i)**2 +                               &
1314                                   ( 0.25_wp * ( v(k,j,i-1)  +                 &
1315                                                 v(k,j,i)    +                 &
1316                                                 v(k,j+1,i)  +                 &
1317                                                 v(k,j+1,i-1) )                &
1318                                   )**2 +                                      &
1319                                   ( 0.25_wp * ( w(k-1,j,i-1) +                &
1320                                                 w(k-1,j,i)   +                &
1321                                                 w(k,j,i-1)   +                &
1322                                                 w(k,j,i) )                    &
1323                                   )**2                                        &
1324                                 ) *                                           &
1325                             u(k,j,i)
1326
1327!
1328!--             Calculate preliminary new velocity, based on pre_tend
1329                pre_u = u(k,j,i) + dt_3d * pre_tend
1330!
1331!--             Compare sign of old velocity and new preliminary velocity,
1332!--             and in case the signs are different, limit the tendency
1333                IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
1334                   pre_tend = - u(k,j,i) * ddt_3d
1335                ELSE
1336                   pre_tend = pre_tend
1337                ENDIF
1338!
1339!--             Calculate final tendency
1340                tend(k,j,i) = tend(k,j,i) + pre_tend
1341             ENDDO
1342
1343
1344!
1345!--       v-component
1346          CASE ( 2 )
1347             DO  k = nzb_v_inner(j,i)+1, nzb_v_inner(j,i)+pch_index
1348
1349                kk = k - nzb_v_inner(j,i)  !- lad arrays are defined flat
1350!
1351!--             In order to create sharp boundaries of the plant canopy,
1352!--             the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
1353!--             rather than being interpolated from the surrounding lad_s,
1354!--             because this would yield smaller lad at the canopy boundaries
1355!--             than inside of the canopy.
1356!--             For the same reason, the lad at the northmost(j+1)canopy
1357!--             boundary on the v-grid equals lad_s(k,j,i).
1358                lad_local = lad_s(kk,j,i)
1359                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp )  THEN
1360                   lad_local = lad_s(kk,j-1,i)
1361                ENDIF
1362
1363                pre_tend = 0.0_wp
1364                pre_v = 0.0_wp
1365!
1366!--             Calculate preliminary value (pre_tend) of the tendency
1367                pre_tend = - cdc *                                             &
1368                             lad_local *                                       &
1369                             SQRT( ( 0.25_wp * ( u(k,j-1,i)   +                &
1370                                                 u(k,j-1,i+1) +                &
1371                                                 u(k,j,i)     +                &
1372                                                 u(k,j,i+1) )                  &
1373                                   )**2 +                                      &
1374                                   v(k,j,i)**2 +                               &
1375                                   ( 0.25_wp * ( w(k-1,j-1,i) +                &
1376                                                 w(k-1,j,i)   +                &
1377                                                 w(k,j-1,i)   +                &
1378                                                 w(k,j,i) )                    &
1379                                   )**2                                        &
1380                                 ) *                                           &
1381                             v(k,j,i)
1382
1383!
1384!--             Calculate preliminary new velocity, based on pre_tend
1385                pre_v = v(k,j,i) + dt_3d * pre_tend
1386!
1387!--             Compare sign of old velocity and new preliminary velocity,
1388!--             and in case the signs are different, limit the tendency
1389                IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
1390                   pre_tend = - v(k,j,i) * ddt_3d
1391                ELSE
1392                   pre_tend = pre_tend
1393                ENDIF
1394!
1395!--             Calculate final tendency
1396                tend(k,j,i) = tend(k,j,i) + pre_tend
1397             ENDDO
1398
1399
1400!
1401!--       w-component
1402          CASE ( 3 )
1403             DO  k = nzb_w_inner(j,i)+1, nzb_w_inner(j,i)+pch_index-1
1404
1405                kk = k - nzb_w_inner(j,i)  !- lad arrays are defined flat
1406
1407                pre_tend = 0.0_wp
1408                pre_w = 0.0_wp
1409!
1410!--             Calculate preliminary value (pre_tend) of the tendency
1411                pre_tend = - cdc *                                             &
1412                             (0.5_wp *                                         &
1413                                ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) *         &
1414                             SQRT( ( 0.25_wp * ( u(k,j,i)    +                 & 
1415                                                 u(k,j,i+1)  +                 &
1416                                                 u(k+1,j,i)  +                 &
1417                                                 u(k+1,j,i+1) )                &
1418                                   )**2 +                                      &
1419                                   ( 0.25_wp * ( v(k,j,i)    +                 &
1420                                                 v(k,j+1,i)  +                 &
1421                                                 v(k+1,j,i)  +                 &
1422                                                 v(k+1,j+1,i) )                &
1423                                   )**2 +                                      &
1424                                   w(k,j,i)**2                                 &
1425                                 ) *                                           &
1426                             w(k,j,i)
1427!
1428!--             Calculate preliminary new velocity, based on pre_tend
1429                pre_w = w(k,j,i) + dt_3d * pre_tend
1430!
1431!--             Compare sign of old velocity and new preliminary velocity,
1432!--             and in case the signs are different, limit the tendency
1433                IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
1434                   pre_tend = - w(k,j,i) * ddt_3d
1435                ELSE
1436                   pre_tend = pre_tend
1437                ENDIF
1438!
1439!--             Calculate final tendency
1440                tend(k,j,i) = tend(k,j,i) + pre_tend
1441             ENDDO
1442
1443!
1444!--       potential temperature
1445          CASE ( 4 )
1446             DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
1447                kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
1448                tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
1449             ENDDO
1450
1451
1452!
1453!--       humidity
1454          CASE ( 5 )
1455             DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
1456                kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
1457                tend(k,j,i) = tend(k,j,i) -                                    &
1458                                 lsec *                                        &
1459                                 lad_s(kk,j,i) *                               &
1460                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
1461                                                    u(k,j,i+1) )               &
1462                                       )**2  +                                 &
1463                                       ( 0.5_wp * ( v(k,j,i) +                 &
1464                                                    v(k,j+1,i) )               &
1465                                       )**2 +                                  &
1466                                       ( 0.5_wp * ( w(k-1,j,i) +               &
1467                                                    w(k,j,i) )                 &
1468                                       )**2                                    &
1469                                     ) *                                       &
1470                                 ( q(k,j,i) - lsc )
1471             ENDDO   
1472
1473!
1474!--       sgs-tke
1475          CASE ( 6 )
1476             DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
1477                kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
1478                tend(k,j,i) = tend(k,j,i) -                                    &
1479                                 2.0_wp * cdc *                                &
1480                                 lad_s(kk,j,i) *                               &
1481                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
1482                                                    u(k,j,i+1) )               &
1483                                       )**2 +                                  & 
1484                                       ( 0.5_wp * ( v(k,j,i) +                 &
1485                                                    v(k,j+1,i) )               &
1486                                       )**2 +                                  &
1487                                       ( 0.5_wp * ( w(k,j,i) +                 &
1488                                                    w(k+1,j,i) )               &
1489                                       )**2                                    &
1490                                     ) *                                       &
1491                                 e(k,j,i)
1492             ENDDO
1493             
1494!
1495!--       scalar concentration
1496          CASE ( 7 )
1497             DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
1498                kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
1499                tend(k,j,i) = tend(k,j,i) -                                    &
1500                                 lsec *                                        &
1501                                 lad_s(kk,j,i) *                               &
1502                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
1503                                                    u(k,j,i+1) )               &
1504                                       )**2  +                                 &
1505                                       ( 0.5_wp * ( v(k,j,i) +                 &
1506                                                    v(k,j+1,i) )               &
1507                                       )**2 +                                  &
1508                                       ( 0.5_wp * ( w(k-1,j,i) +               &
1509                                                    w(k,j,i) )                 &
1510                                       )**2                                    &
1511                                     ) *                                       &
1512                                 ( s(k,j,i) - lsc )
1513             ENDDO               
1514
1515       CASE DEFAULT
1516
1517          WRITE( message_string, * ) 'wrong component: ', component
1518          CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 ) 
1519
1520       END SELECT
1521
1522    END SUBROUTINE pcm_tendency_ij
1523
1524
1525
1526 END MODULE plant_canopy_model_mod
Note: See TracBrowser for help on using the repository browser.