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

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

Minor format changes

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