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

Last change on this file since 3049 was 3049, checked in by Giersch, 6 years ago

Revision history corrected

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