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

Last change on this file since 2750 was 2746, checked in by suehring, 7 years ago

Read information from statitic driver for resolved vegetation independently from land- or urban-surface model

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