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

Last change on this file since 331 was 257, checked in by heinze, 16 years ago

Output of messages replaced by message handling routine

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