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

Last change on this file since 2925 was 2920, checked in by kanani, 6 years ago

Optimize SVF calculation, clean-up, bugfixes

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