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

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

small bugfix, formatting and new PCM output

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