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

Last change on this file since 3302 was 3294, checked in by raasch, 6 years ago

modularization of the ocean code

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