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

Last change on this file since 3046 was 3046, checked in by Giersch, 3 years ago

Remaining error messages revised, comments extended

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