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

Last change on this file since 3414 was 3337, checked in by kanani, 5 years ago

reintegrate branch resler to trunk

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