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

Last change on this file since 2766 was 2766, checked in by kanani, 4 years ago

Removal of chem directive, plus minor changes

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