source: palm/tags/release-5.0/SOURCE/plant_canopy_model_mod.f90 @ 4095

Last change on this file since 4095 was 2701, checked in by suehring, 6 years ago

changes from last commit documented

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