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

Last change on this file since 153 was 153, checked in by steinfeld, 14 years ago

Update for the plant canopy model

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