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

Last change on this file since 3291 was 3274, checked in by knoop, 6 years ago

Modularization of all bulk cloud physics code components

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