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

Last change on this file since 3246 was 3246, checked in by sward, 3 years ago

Added error handling for wrong input parameters

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