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

Last change on this file since 3022 was 3022, checked in by suehring, 3 years ago

Revise recent bugfix in nested runs at left and south boundary; bugfix in advection of u in case of OpenMP parallelization; bugfix in plant transpiration

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