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

Last change on this file since 2769 was 2768, checked in by kanani, 7 years ago

Added parameter check, reduced line length, some formatting

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