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

Last change on this file since 3014 was 3014, checked in by maronga, 3 years ago

series of bugfixes

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