source: palm/trunk/SOURCE/plant_canopy_model.f90 @ 1266

Last change on this file since 1266 was 1037, checked in by raasch, 12 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 14.2 KB
Line 
1 MODULE plant_canopy_model_mod
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later 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-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: plant_canopy_model.f90 1037 2012-10-22 14:10:22Z heinze $
27!
28! 1036 2012-10-22 13:43:42Z raasch
29! code put under GPL (PALM 3.9)
30!
31! 257 2009-03-11 15:17:42Z heinze
32! Output of messages replaced by message handling routine.
33! Bugfix: remove IF statement in plant_canopy_model_ij
34!
35! 153 2008-03-19 09:41:30Z steinfeld
36! heat sources within the forest canopy are added, which represent the
37! rate of heat input into the air from the forest leaves, evaluation of sinks
38! and sources for scalar concentration due to canopy elements
39!
40! 138 2007-11-28 10:03:58Z letzel
41! Initial revision
42!
43! Description:
44! ------------
45! Evaluation of sinks and sources of momentum, heat and scalar concentration
46! due to canopy elements
47!------------------------------------------------------------------------------!
48
49    PRIVATE
50    PUBLIC plant_canopy_model
51
52    INTERFACE plant_canopy_model
53       MODULE PROCEDURE plant_canopy_model
54       MODULE PROCEDURE plant_canopy_model_ij
55    END INTERFACE plant_canopy_model
56
57 CONTAINS
58
59
60!------------------------------------------------------------------------------!
61! Call for all grid points
62!------------------------------------------------------------------------------!
63    SUBROUTINE plant_canopy_model( component )
64
65       USE arrays_3d
66       USE control_parameters
67       USE indices
68       USE pegrid
69
70       IMPLICIT NONE
71
72       INTEGER ::  component, i, j, k
73 
74!
75!--    Compute drag for the three velocity components and the SGS-TKE
76       SELECT CASE ( component )
77
78!
79!--       u-component
80          CASE ( 1 )
81             DO  i = nxlu, nxr
82                DO  j = nys, nyn
83                   DO  k = nzb_u_inner(j,i)+1, pch_index
84                      tend(k,j,i) = tend(k,j,i) -                &
85                                    cdc(k,j,i) * lad_u(k,j,i) *  &
86                                    SQRT(     u(k,j,i)**2     +  &
87                                          ( ( v(k,j,i-1)      +  &
88                                              v(k,j,i)        +  &
89                                              v(k,j+1,i)      +  &
90                                              v(k,j+1,i-1) )     &
91                                            / 4.0 )**2        +  &
92                                          ( ( w(k-1,j,i-1)    +  &
93                                              w(k-1,j,i)      +  &
94                                              w(k,j,i-1)      +  &
95                                              w(k,j,i) )         &
96                                            / 4.0 )**2 )      *  &
97                                    u(k,j,i)
98                   ENDDO
99                ENDDO
100             ENDDO
101
102!
103!--       v-component
104          CASE ( 2 )
105             DO  i = nxl, nxr
106                DO  j = nysv, nyn
107                   DO  k = nzb_v_inner(j,i)+1, pch_index
108                      tend(k,j,i) = tend(k,j,i) -                &
109                                    cdc(k,j,i) * lad_v(k,j,i) *  &
110                                    SQRT( ( ( u(k,j-1,i)      +  &
111                                              u(k,j-1,i+1)    +  &
112                                              u(k,j,i)        +  &
113                                              u(k,j,i+1) )       &
114                                            / 4.0 )**2        +  &
115                                              v(k,j,i)**2     +  &
116                                          ( ( w(k-1,j-1,i)    +  &
117                                              w(k-1,j,i)      +  &
118                                              w(k,j-1,i)      +  &
119                                              w(k,j,i) )         &
120                                            / 4.0 )**2 )      *  &
121                                    v(k,j,i) 
122                   ENDDO
123                ENDDO
124             ENDDO
125
126!
127!--       w-component
128          CASE ( 3 )
129             DO  i = nxl, nxr
130                DO  j = nys, nyn
131                   DO  k = nzb_w_inner(j,i)+1, pch_index
132                      tend(k,j,i) = tend(k,j,i) -                & 
133                                    cdc(k,j,i) * lad_w(k,j,i) *  &
134                                    SQRT( ( ( u(k,j,i)        +  &
135                                              u(k,j,i+1)      +  &
136                                              u(k+1,j,i)      +  &
137                                              u(k+1,j,i+1) )     &
138                                            / 4.0 )**2        +  &
139                                          ( ( v(k,j,i)        +  &
140                                              v(k,j+1,i)      +  &
141                                              v(k+1,j,i)      +  &
142                                              v(k+1,j+1,i) )     &
143                                            / 4.0 )**2        +  &
144                                              w(k,j,i)**2 )   *  &
145                                    w(k,j,i)
146                   ENDDO
147                ENDDO
148             ENDDO
149
150!
151!--       potential temperature
152          CASE ( 4 )
153             DO  i = nxl, nxr
154                DO  j = nys, nyn
155                   DO  k = nzb_s_inner(j,i)+1, pch_index
156                      tend(k,j,i) = tend(k,j,i) +                     &
157                                    ( canopy_heat_flux(k,j,i) -     &
158                                      canopy_heat_flux(k-1,j,i) ) / &
159                                      dzw(k)
160                   ENDDO
161                ENDDO
162             ENDDO
163
164!
165!--       scalar concentration
166          CASE ( 5 )
167             DO  i = nxl, nxr
168                DO  j = nys, nyn
169                   DO  k = nzb_s_inner(j,i)+1, pch_index
170                      tend(k,j,i) = tend(k,j,i) -                     &
171                                    sec(k,j,i) * lad_s(k,j,i) *       &
172                                    SQRT( ( ( u(k,j,i)        +       &
173                                              u(k,j,i+1) )            &
174                                            / 2.0 )**2        +       &
175                                          ( ( v(k,j,i)        +       &
176                                              v(k,j+1,i) )            &
177                                            / 2.0 )**2        +       &
178                                          ( ( w(k-1,j,i)      +       & 
179                                              w(k,j,i) )              &
180                                            / 2.0 )**2 )      *       &
181                                    ( q(k,j,i) - sls(k,j,i) )
182                   ENDDO
183                ENDDO
184             ENDDO
185
186!
187!--       sgs-tke
188          CASE ( 6 )
189             DO  i = nxl, nxr
190                DO  j = nys, nyn
191                   DO  k = nzb_s_inner(j,i)+1, pch_index
192                      tend(k,j,i) = tend(k,j,i) -                     &
193                                    2.0 * cdc(k,j,i) * lad_s(k,j,i) * &
194                                    SQRT( ( ( u(k,j,i)              + &
195                                              u(k,j,i+1) )            &
196                                            / 2.0 )**2              + &
197                                          ( ( v(k,j,i)              + &
198                                              v(k,j+1,i) )            &
199                                            / 2.0 )**2              + &
200                                          ( ( w(k,j,i)              + &
201                                              w(k+1,j,i) )            &
202                                            / 2.0 )**2 )            * &
203                                    e(k,j,i)
204                   ENDDO
205                ENDDO
206             ENDDO 
207                       
208          CASE DEFAULT
209
210             WRITE( message_string, * ) 'wrong component: ', component
211             CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 ) 
212
213       END SELECT
214
215    END SUBROUTINE plant_canopy_model
216
217
218!------------------------------------------------------------------------------!
219! Call for grid point i,j
220!------------------------------------------------------------------------------!
221    SUBROUTINE plant_canopy_model_ij( i, j, component )
222
223       USE arrays_3d
224       USE control_parameters
225       USE indices
226       USE pegrid
227
228       IMPLICIT NONE
229
230       INTEGER ::  component, i, j, k
231
232!
233!--    Compute drag for the three velocity components
234       SELECT CASE ( component )
235
236!
237!--       u-component
238       CASE ( 1 )
239          DO  k = nzb_u_inner(j,i)+1, pch_index
240             tend(k,j,i) = tend(k,j,i) -                  &
241                              cdc(k,j,i) * lad_u(k,j,i) *    &   
242                              SQRT(     u(k,j,i)**2 +        &
243                                    ( ( v(k,j,i-1)  +        &
244                                        v(k,j,i)    +        &
245                                        v(k,j+1,i)  +        &
246                                        v(k,j+1,i-1) )       &
247                                      / 4.0 )**2    +        &
248                                    ( ( w(k-1,j,i-1) +       &
249                                        w(k-1,j,i)   +       &
250                                        w(k,j,i-1)   +       &
251                                        w(k,j,i) )           &
252                                      / 4.0 )**2 ) *         &
253                              u(k,j,i)
254          ENDDO
255
256!
257!--       v-component
258       CASE ( 2 )
259          DO  k = nzb_v_inner(j,i)+1, pch_index
260             tend(k,j,i) = tend(k,j,i) -                  &
261                              cdc(k,j,i) * lad_v(k,j,i) *    &
262                              SQRT( ( ( u(k,j-1,i)   +       &
263                                        u(k,j-1,i+1) +       &
264                                        u(k,j,i)     +       &
265                                        u(k,j,i+1) )         &
266                                      / 4.0 )**2     +       &
267                                        v(k,j,i)**2  +       &
268                                    ( ( w(k-1,j-1,i) +       &
269                                        w(k-1,j,i)   +       &
270                                        w(k,j-1,i)   +       &
271                                        w(k,j,i) )           &
272                                      / 4.0 )**2 ) *         &
273                              v(k,j,i)
274          ENDDO
275
276!
277!--       w-component
278       CASE ( 3 )
279          DO  k = nzb_w_inner(j,i)+1, pch_index
280             tend(k,j,i) = tend(k,j,i) -                  &
281                              cdc(k,j,i) * lad_w(k,j,i) *    & 
282                              SQRT( ( ( u(k,j,i)    +        & 
283                                        u(k,j,i+1)  +        &
284                                        u(k+1,j,i)  +        &
285                                        u(k+1,j,i+1) )       &
286                                      / 4.0 )**2    +        &
287                                    ( ( v(k,j,i)    +        &
288                                        v(k,j+1,i)  +        &
289                                        v(k+1,j,i)  +        &
290                                        v(k+1,j+1,i) )       &
291                                      / 4.0 )**2    +        &
292                                        w(k,j,i)**2 ) *      &
293                              w(k,j,i)
294   
295          ENDDO
296
297!
298!--       potential temperature
299          CASE ( 4 )
300             DO  k = nzb_s_inner(j,i)+1, pch_index
301                tend(k,j,i) = tend(k,j,i) +                     &
302                              ( canopy_heat_flux(k,j,i) -     &
303                                canopy_heat_flux(k-1,j,i) ) / &
304                                dzw(k)
305             ENDDO
306
307
308!
309!--       scalar concentration
310          CASE ( 5 )
311             DO  k = nzb_s_inner(j,i)+1, pch_index
312                tend(k,j,i) = tend(k,j,i) -                     &
313                              sec(k,j,i) * lad_s(k,j,i) *       &
314                              SQRT( ( ( u(k,j,i)        +       &
315                                        u(k,j,i+1) )            &
316                                      / 2.0 )**2        +       &
317                                    ( ( v(k,j,i)        +       &
318                                        v(k,j+1,i) )            &
319                                      / 2.0 )**2        +       &
320                                    ( ( w(k-1,j,i)      +       &
321                                        w(k,j,i) )              &
322                                      / 2.0 )**2 )      *       &
323                              ( q(k,j,i) - sls(k,j,i) )
324             ENDDO   
325
326!
327!--       sgs-tke
328       CASE ( 6 )
329          DO  k = nzb_s_inner(j,i)+1, pch_index   
330             tend(k,j,i) = tend(k,j,i) -                     &
331                              2.0 * cdc(k,j,i) * lad_s(k,j,i) * &
332                              SQRT( ( ( u(k,j,i)           +    &
333                                        u(k,j,i+1) )            &
334                                      / 2.0 )**2           +    & 
335                                    ( ( v(k,j,i)           +    &
336                                        v(k,j+1,i) )            &
337                                      / 2.0 )**2           +    &
338                                    ( ( w(k,j,i)           +    &
339                                        w(k+1,j,i) )            &
340                                      / 2.0 )**2 )         *    &
341                              e(k,j,i)
342          ENDDO
343
344       CASE DEFAULT
345
346          WRITE( message_string, * ) 'wrong component: ', component
347          CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 ) 
348
349       END SELECT
350
351    END SUBROUTINE plant_canopy_model_ij
352
353 END MODULE plant_canopy_model_mod
Note: See TracBrowser for help on using the repository browser.