source: palm/tags/release-3.4a/SOURCE/plant_canopy_model.f90 @ 710

Last change on this file since 710 was 139, checked in by raasch, 16 years ago

New:
---

Plant canopy model of Watanabe (2004,BLM 112,307-341) added.
It can be switched on by the inipar parameter plant_canopy.
The inipar parameter canopy_mode can be used to prescribe a
plant canopy type. The default case is a homogeneous plant
canopy. Heterogeneous distributions of the leaf area
density and the canopy drag coefficient can be defined in the
new routine user_init_plant_canopy (user_interface).
The inipar parameters lad_surface, lad_vertical_gradient and
lad_vertical_gradient_level can be used in order to
prescribe the vertical profile of leaf area density. The
inipar parameter drag_coefficient determines the canopy
drag coefficient.
Finally, the inipar parameter pch_index determines the
index of the upper boundary of the plant canopy.

Allow new case bc_uv_t = 'dirichlet_0' for channel flow.

For unknown variables (CASE DEFAULT) call new subroutine user_data_output_dvrp

Pressure boundary conditions for vertical walls added to the multigrid solver.
They are applied using new wall flag arrays (wall_flags_..) which are defined
for each grid level. New argument gls added to routine user_init_grid
(user_interface).

Frequence of sorting particles can be controlled with new particles_par
parameter dt_sort_particles. Sorting is moved from the SGS timestep loop in
advec_particles after the end of this loop.

advec_particles, check_parameters, data_output_dvrp, header, init_3d_model, init_grid, init_particles, init_pegrid, modules, package_parin, parin, plant_canopy_model, read_var_list, read_3d_binary, user_interface, write_var_list, write_3d_binary

Changed:


Redefine initial nzb_local as the actual total size of topography (later the
extent of topography in nzb_local is reduced by 1dx at the E topography walls
and by 1dy at the N topography walls to form the basis for nzb_s_inner);
for consistency redefine 'single_building' case.

Vertical profiles now based on nzb_s_inner; they are divided by
ngp_2dh_s_inner (scalars, procucts of scalars) and ngp_2dh (staggered velocity
components and their products, procucts of scalars and velocity components),
respectively.

Allow two instead of one digit to specify isosurface and slicer variables.

Status of 3D-volume NetCDF data file only depends on switch netcdf_64bit_3d (check_open)

prognostic_equations include the respective wall_*flux in the parameter list of
calls of diffusion_s. Same as before, only the values of wall_heatflux(0:4)
can be assigned. At present, wall_humidityflux, wall_qflux, wall_salinityflux,
and wall_scalarflux are kept zero. diffusion_s uses the respective wall_*flux
instead of wall_heatflux. This update serves two purposes:

  • it avoids errors in calculations with humidity/scalar/salinity and prescribed

non-zero wall_heatflux,

  • it prepares PALM for a possible assignment of wall fluxes of

humidity/scalar/salinity in a future release.

buoyancy, check_open, data_output_dvrp, diffusion_s, diffusivities, flow_statistics, header, init_3d_model, init_dvrp, init_grid, modules, prognostic_equations

Errors:


Bugfix: summation of sums_l_l in diffusivities.

Several bugfixes in the ocean part: Initial density rho is calculated
(init_ocean). Error in initializing u_init and v_init removed
(check_parameters). Calculation of density flux now starts from
nzb+1 (production_e).

Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
numbers along y, small bugfixes in the SGS part (advec_particles)

Bugfix: model_string needed a default value (combine_plot_fields)

Bugfix: wavenumber calculation for even nx in routines maketri (poisfft)

Bugfix: assignment of fluxes at walls

Bugfix: absolute value of f must be used when calculating the Blackadar mixing length (init_1d_model)

advec_particles, check_parameters, combine_plot_fields, diffusion_s, diffusivities, init_ocean, init_1d_model, poisfft, production_e

  • Property svn:keywords set to Id
File size: 10.5 KB
Line 
1 MODULE plant_canopy_model_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: plant_canopy_model.f90 139 2007-11-29 09:37:41Z raasch $
11!
12! 138 2007-11-28 10:03:58Z letzel
13! Initial revision
14!
15! Description:
16! ------------
17! Evaluation of the drag due to vegetation
18!------------------------------------------------------------------------------!
19
20    PRIVATE
21    PUBLIC plant_canopy_model
22
23    INTERFACE plant_canopy_model
24       MODULE PROCEDURE plant_canopy_model
25       MODULE PROCEDURE plant_canopy_model_ij
26    END INTERFACE plant_canopy_model
27
28 CONTAINS
29
30
31!------------------------------------------------------------------------------!
32! Call for all grid points
33!------------------------------------------------------------------------------!
34    SUBROUTINE plant_canopy_model( component )
35
36       USE arrays_3d
37       USE control_parameters
38       USE indices
39       USE pegrid
40
41       IMPLICIT NONE
42
43       INTEGER ::  component, i, j, k
44 
45!
46!--    Compute drag for the three velocity components and the SGS-TKE
47       SELECT CASE ( component )
48
49!
50!--       u-component
51          CASE ( 1 )
52             DO  i = nxlu, nxr
53                DO  j = nys, nyn
54                   DO  k = nzb_u_inner(j,i)+1, pch_index
55                      tend(k,j,i) = tend(k,j,i) -                &
56                                    cdc(k,j,i) * lad_u(k,j,i) *  &
57                                    SQRT(     u(k,j,i)**2     +  &
58                                          ( ( v(k,j,i-1)      +  &
59                                              v(k,j,i)        +  &
60                                              v(k,j+1,i)      +  &
61                                              v(k,j+1,i+1) )     &
62                                            / 4.0 )**2        +  &
63                                          ( ( w(k-1,j,i-1)    +  &
64                                              w(k-1,j,i)      +  &
65                                              w(k,j,i-1)      +  &
66                                              w(k,j,i) )         &
67                                            / 4.0 )**2 )      *  &
68                                    u(k,j,i)
69                   ENDDO
70                ENDDO
71             ENDDO
72
73!
74!--       v-component
75          CASE ( 2 )
76             DO  i = nxl, nxr
77                DO  j = nysv, nyn
78                   DO  k = nzb_v_inner(j,i)+1, pch_index
79                      tend(k,j,i) = tend(k,j,i) -                &
80                                    cdc(k,j,i) * lad_v(k,j,i) *  &
81                                    SQRT( ( ( u(k,j-1,i)      +  &
82                                              u(k,j-1,i+1)    +  &
83                                              u(k,j,i)        +  &
84                                              u(k,j,i+1) )       &
85                                            / 4.0 )**2        +  &
86                                              v(k,j,i)**2     +  &
87                                          ( ( w(k-1,j-1,i)    +  &
88                                              w(k-1,j,i)      +  &
89                                              w(k,j-1,i)      +  &
90                                              w(k,j,i) )         &
91                                            / 4.0 )**2 )      *  &
92                                    v(k,j,i) 
93                   ENDDO
94                ENDDO
95             ENDDO
96
97!
98!--       w-component
99          CASE ( 3 )
100             DO  i = nxl, nxr
101                DO  j = nys, nyn
102                   DO  k = nzb_w_inner(j,i)+1, pch_index
103                      tend(k,j,i) = tend(k,j,i) -                & 
104                                    cdc(k,j,i) * lad_w(k,j,i) *  &
105                                    SQRT( ( ( u(k,j,i)        +  &
106                                              u(k,j,i+1)      +  &
107                                              u(k+1,j,i)      +  &
108                                              u(k+1,j,i+1) )     &
109                                            / 4.0 )**2        +  &
110                                          ( ( v(k,j,i)        +  &
111                                              v(k,j+1,i)      +  &
112                                              v(k+1,j,i)      +  &
113                                              v(k+1,j+1,i) )     &
114                                            / 4.0 )**2        +  &
115                                              w(k,j,i)**2 )   *  &
116                                    w(k,j,i)
117                   ENDDO
118                ENDDO
119             ENDDO
120
121!
122!--       sgs-tke
123          CASE ( 4 )
124             DO  i = nxl, nxr
125                DO  j = nys, nyn
126                   DO  k = nzb_s_inner(j,i)+1, pch_index
127                      tend(k,j,i) = tend(k,j,i) -                     &
128                                    2.0 * cdc(k,j,i) * lad_s(k,j,i) * &
129                                    SQRT( ( ( u(k,j,i)              + &
130                                              u(k,j,i+1) )            &
131                                            / 2.0 )**2              + &
132                                          ( ( v(k,j,i)              + &
133                                              v(k,j+1,i) )            &
134                                            / 2.0 )**2              + &
135                                          ( ( w(k,j,i)              + &
136                                              w(k+1,j,i) )            &
137                                            / 2.0 )**2 )            * &
138                                    e(k,j,i)
139                   ENDDO
140                ENDDO
141             ENDDO 
142                       
143          CASE DEFAULT
144
145             IF ( myid == 0 )  PRINT*,'+++ pcm:  wrong component: ', &
146                                      component
147             CALL local_stop
148
149       END SELECT
150
151    END SUBROUTINE plant_canopy_model
152
153
154!------------------------------------------------------------------------------!
155! Call for grid point i,j
156!------------------------------------------------------------------------------!
157    SUBROUTINE plant_canopy_model_ij( i, j, component )
158
159       USE arrays_3d
160       USE control_parameters
161       USE indices
162       USE pegrid
163
164       IMPLICIT NONE
165
166       INTEGER ::  component, i, j, k
167
168       IF ( drag_coefficient /= 0.0 ) THEN
169
170!
171!--       Compute drag for the three velocity components
172          SELECT CASE ( component )
173
174!
175!--          u-component
176             CASE ( 1 )
177                DO  k = nzb_u_inner(j,i)+1, pch_index
178                   tend(k,j,i) = tend(k,j,i) -                  &
179                                 cdc(k,j,i) * lad_u(k,j,i) *    &   
180                                 SQRT(     u(k,j,i)**2 +        &
181                                       ( ( v(k,j,i-1)  +        &
182                                           v(k,j,i)    +        &
183                                           v(k,j+1,i)  +        &
184                                           v(k,j+1,i+1) )       &
185                                         / 4.0 )**2    +        &
186                                       ( ( w(k-1,j,i-1) +       &
187                                           w(k-1,j,i)   +       &
188                                           w(k,j,i-1)   +       &
189                                           w(k,j,i) )           &
190                                         / 4.0 )**2 ) *         &
191                                 u(k,j,i)
192                ENDDO
193
194!
195!--          v-component
196             CASE ( 2 )
197                DO  k = nzb_v_inner(j,i)+1, pch_index
198                   tend(k,j,i) = tend(k,j,i) -                  &
199                                 cdc(k,j,i) * lad_v(k,j,i) *    &
200                                 SQRT( ( ( u(k,j-1,i)   +       &
201                                           u(k,j-1,i+1) +       &
202                                           u(k,j,i)     +       &
203                                           u(k,j,i+1) )         &
204                                         / 4.0 )**2     +       &
205                                           v(k,j,i)**2  +       &
206                                       ( ( w(k-1,j-1,i) +       &
207                                           w(k-1,j,i)   +       &
208                                           w(k,j-1,i)   +       &
209                                           w(k,j,i) )           &
210                                         / 4.0 )**2 ) *         &
211                                 v(k,j,i)
212                ENDDO
213
214!
215!--          w-component
216             CASE ( 3 )
217                DO  k = nzb_w_inner(j,i)+1, pch_index
218                   tend(k,j,i) = tend(k,j,i) -                  &
219                                 cdc(k,j,i) * lad_w(k,j,i) *    & 
220                                 SQRT( ( ( u(k,j,i)    +        & 
221                                           u(k,j,i+1)  +        &
222                                           u(k+1,j,i)  +        &
223                                           u(k+1,j,i+1) )       &
224                                         / 4.0 )**2    +        &
225                                       ( ( v(k,j,i)    +        &
226                                           v(k,j+1,i)  +        &
227                                           v(k+1,j,i)  +        &
228                                           v(k+1,j+1,i) )       &
229                                         / 4.0 )**2    +        &
230                                           w(k,j,i)**2 ) *      &
231                                 w(k,j,i)
232   
233                ENDDO
234
235!
236!--          sgs-tke
237             CASE ( 4 )
238                DO  k = nzb_s_inner(j,i)+1, pch_index   
239                   tend(k,j,i) = tend(k,j,i) -                     &
240                                 2.0 * cdc(k,j,i) * lad_s(k,j,i) * &
241                                 SQRT( ( ( u(k,j,i)           +    &
242                                           u(k,j,i+1) )            &
243                                         / 2.0 )**2           +    & 
244                                       ( ( v(k,j,i)           +    &
245                                           v(k,j+1,i) )            &
246                                         / 2.0 )**2           +    &
247                                       ( ( w(k,j,i)           +    &
248                                           w(k+1,j,i) )            &
249                                         / 2.0 )**2 )         *    &
250                                 e(k,j,i)
251
252                ENDDO
253
254             CASE DEFAULT
255
256                IF ( myid == 0 )  PRINT*,'+++ pcm:  wrong component: ', &
257                                         component
258                CALL local_stop
259
260          END SELECT
261
262       ENDIF
263
264    END SUBROUTINE plant_canopy_model_ij
265
266 END MODULE plant_canopy_model_mod
Note: See TracBrowser for help on using the repository browser.