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

Last change on this file since 249 was 198, checked in by raasch, 16 years ago

file headers updated for the next release 3.5

  • 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!
7!
8! Bugfix: remove IF statement in plant_canopy_model_ij
9!
10! Former revisions:
11! -----------------
12! $Id: plant_canopy_model.f90 198 2008-09-17 08:55:28Z letzel $
13!
14! 153 2008-03-19 09:41:30Z steinfeld
15! heat sources within the forest canopy are added, which represent the
16! rate of heat input into the air from the forest leaves, evaluation of sinks
17! and sources for scalar concentration due to canopy elements
18!
19! 138 2007-11-28 10:03:58Z letzel
20! Initial revision
21!
22! Description:
23! ------------
24! Evaluation of sinks and sources of momentum, heat and scalar concentration
25! due to canopy elements
26!------------------------------------------------------------------------------!
27
28    PRIVATE
29    PUBLIC plant_canopy_model
30
31    INTERFACE plant_canopy_model
32       MODULE PROCEDURE plant_canopy_model
33       MODULE PROCEDURE plant_canopy_model_ij
34    END INTERFACE plant_canopy_model
35
36 CONTAINS
37
38
39!------------------------------------------------------------------------------!
40! Call for all grid points
41!------------------------------------------------------------------------------!
42    SUBROUTINE plant_canopy_model( component )
43
44       USE arrays_3d
45       USE control_parameters
46       USE indices
47       USE pegrid
48
49       IMPLICIT NONE
50
51       INTEGER ::  component, i, j, k
52 
53!
54!--    Compute drag for the three velocity components and the SGS-TKE
55       SELECT CASE ( component )
56
57!
58!--       u-component
59          CASE ( 1 )
60             DO  i = nxlu, nxr
61                DO  j = nys, nyn
62                   DO  k = nzb_u_inner(j,i)+1, pch_index
63                      tend(k,j,i) = tend(k,j,i) -                &
64                                    cdc(k,j,i) * lad_u(k,j,i) *  &
65                                    SQRT(     u(k,j,i)**2     +  &
66                                          ( ( v(k,j,i-1)      +  &
67                                              v(k,j,i)        +  &
68                                              v(k,j+1,i)      +  &
69                                              v(k,j+1,i-1) )     &
70                                            / 4.0 )**2        +  &
71                                          ( ( w(k-1,j,i-1)    +  &
72                                              w(k-1,j,i)      +  &
73                                              w(k,j,i-1)      +  &
74                                              w(k,j,i) )         &
75                                            / 4.0 )**2 )      *  &
76                                    u(k,j,i)
77                   ENDDO
78                ENDDO
79             ENDDO
80
81!
82!--       v-component
83          CASE ( 2 )
84             DO  i = nxl, nxr
85                DO  j = nysv, nyn
86                   DO  k = nzb_v_inner(j,i)+1, pch_index
87                      tend(k,j,i) = tend(k,j,i) -                &
88                                    cdc(k,j,i) * lad_v(k,j,i) *  &
89                                    SQRT( ( ( u(k,j-1,i)      +  &
90                                              u(k,j-1,i+1)    +  &
91                                              u(k,j,i)        +  &
92                                              u(k,j,i+1) )       &
93                                            / 4.0 )**2        +  &
94                                              v(k,j,i)**2     +  &
95                                          ( ( w(k-1,j-1,i)    +  &
96                                              w(k-1,j,i)      +  &
97                                              w(k,j-1,i)      +  &
98                                              w(k,j,i) )         &
99                                            / 4.0 )**2 )      *  &
100                                    v(k,j,i) 
101                   ENDDO
102                ENDDO
103             ENDDO
104
105!
106!--       w-component
107          CASE ( 3 )
108             DO  i = nxl, nxr
109                DO  j = nys, nyn
110                   DO  k = nzb_w_inner(j,i)+1, pch_index
111                      tend(k,j,i) = tend(k,j,i) -                & 
112                                    cdc(k,j,i) * lad_w(k,j,i) *  &
113                                    SQRT( ( ( u(k,j,i)        +  &
114                                              u(k,j,i+1)      +  &
115                                              u(k+1,j,i)      +  &
116                                              u(k+1,j,i+1) )     &
117                                            / 4.0 )**2        +  &
118                                          ( ( v(k,j,i)        +  &
119                                              v(k,j+1,i)      +  &
120                                              v(k+1,j,i)      +  &
121                                              v(k+1,j+1,i) )     &
122                                            / 4.0 )**2        +  &
123                                              w(k,j,i)**2 )   *  &
124                                    w(k,j,i)
125                   ENDDO
126                ENDDO
127             ENDDO
128
129!
130!--       potential temperature
131          CASE ( 4 )
132             DO  i = nxl, nxr
133                DO  j = nys, nyn
134                   DO  k = nzb_s_inner(j,i)+1, pch_index
135                      tend(k,j,i) = tend(k,j,i) +                     &
136                                    ( canopy_heat_flux(k,j,i) -     &
137                                      canopy_heat_flux(k-1,j,i) ) / &
138                                      dzw(k)
139                   ENDDO
140                ENDDO
141             ENDDO
142
143!
144!--       scalar concentration
145          CASE ( 5 )
146             DO  i = nxl, nxr
147                DO  j = nys, nyn
148                   DO  k = nzb_s_inner(j,i)+1, pch_index
149                      tend(k,j,i) = tend(k,j,i) -                     &
150                                    sec(k,j,i) * lad_s(k,j,i) *       &
151                                    SQRT( ( ( u(k,j,i)        +       &
152                                              u(k,j,i+1) )            &
153                                            / 2.0 )**2        +       &
154                                          ( ( v(k,j,i)        +       &
155                                              v(k,j+1,i) )            &
156                                            / 2.0 )**2        +       &
157                                          ( ( w(k-1,j,i)      +       & 
158                                              w(k,j,i) )              &
159                                            / 2.0 )**2 )      *       &
160                                    ( q(k,j,i) - sls(k,j,i) )
161                   ENDDO
162                ENDDO
163             ENDDO
164
165!
166!--       sgs-tke
167          CASE ( 6 )
168             DO  i = nxl, nxr
169                DO  j = nys, nyn
170                   DO  k = nzb_s_inner(j,i)+1, pch_index
171                      tend(k,j,i) = tend(k,j,i) -                     &
172                                    2.0 * cdc(k,j,i) * lad_s(k,j,i) * &
173                                    SQRT( ( ( u(k,j,i)              + &
174                                              u(k,j,i+1) )            &
175                                            / 2.0 )**2              + &
176                                          ( ( v(k,j,i)              + &
177                                              v(k,j+1,i) )            &
178                                            / 2.0 )**2              + &
179                                          ( ( w(k,j,i)              + &
180                                              w(k+1,j,i) )            &
181                                            / 2.0 )**2 )            * &
182                                    e(k,j,i)
183                   ENDDO
184                ENDDO
185             ENDDO 
186                       
187          CASE DEFAULT
188
189             IF ( myid == 0 )  PRINT*,'+++ pcm:  wrong component: ', &
190                                      component
191             CALL local_stop
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          IF ( myid == 0 )  PRINT*,'+++ pcm:  wrong component: ', &
327                                      component
328          CALL local_stop
329
330       END SELECT
331
332    END SUBROUTINE plant_canopy_model_ij
333
334 END MODULE plant_canopy_model_mod
Note: See TracBrowser for help on using the repository browser.