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

Last change on this file since 150 was 142, checked in by letzel, 16 years ago

Bugfix in plant_canopy_model: remove IF statement in plant_canopy_model_ij
Bugfix in flow_statistics: divide sums(k,8) (e) and sums(k,34) (e*) by
ngp_2dh_s_inner(k,sr) (like other scalars)

  • Property svn:keywords set to Id
File size: 10.2 KB
Line 
1 MODULE plant_canopy_model_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! Bugfix: remove IF statement in plant_canopy_model_ij
7!
8! Former revisions:
9! -----------------
10! $Id: plant_canopy_model.f90 142 2007-12-04 12:36:44Z 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!
169!--    Compute drag for the three velocity components
170       SELECT CASE ( component )
171
172!
173!--       u-component
174       CASE ( 1 )
175          DO  k = nzb_u_inner(j,i)+1, pch_index
176             tend(k,j,i) = tend(k,j,i) -                  &
177                              cdc(k,j,i) * lad_u(k,j,i) *    &   
178                              SQRT(     u(k,j,i)**2 +        &
179                                    ( ( v(k,j,i-1)  +        &
180                                        v(k,j,i)    +        &
181                                        v(k,j+1,i)  +        &
182                                        v(k,j+1,i+1) )       &
183                                      / 4.0 )**2    +        &
184                                    ( ( w(k-1,j,i-1) +       &
185                                        w(k-1,j,i)   +       &
186                                        w(k,j,i-1)   +       &
187                                        w(k,j,i) )           &
188                                      / 4.0 )**2 ) *         &
189                              u(k,j,i)
190          ENDDO
191
192!
193!--       v-component
194       CASE ( 2 )
195          DO  k = nzb_v_inner(j,i)+1, pch_index
196             tend(k,j,i) = tend(k,j,i) -                  &
197                              cdc(k,j,i) * lad_v(k,j,i) *    &
198                              SQRT( ( ( u(k,j-1,i)   +       &
199                                        u(k,j-1,i+1) +       &
200                                        u(k,j,i)     +       &
201                                        u(k,j,i+1) )         &
202                                      / 4.0 )**2     +       &
203                                        v(k,j,i)**2  +       &
204                                    ( ( w(k-1,j-1,i) +       &
205                                        w(k-1,j,i)   +       &
206                                        w(k,j-1,i)   +       &
207                                        w(k,j,i) )           &
208                                      / 4.0 )**2 ) *         &
209                              v(k,j,i)
210          ENDDO
211
212!
213!--       w-component
214       CASE ( 3 )
215          DO  k = nzb_w_inner(j,i)+1, pch_index
216             tend(k,j,i) = tend(k,j,i) -                  &
217                              cdc(k,j,i) * lad_w(k,j,i) *    & 
218                              SQRT( ( ( u(k,j,i)    +        & 
219                                        u(k,j,i+1)  +        &
220                                        u(k+1,j,i)  +        &
221                                        u(k+1,j,i+1) )       &
222                                      / 4.0 )**2    +        &
223                                    ( ( v(k,j,i)    +        &
224                                        v(k,j+1,i)  +        &
225                                        v(k+1,j,i)  +        &
226                                        v(k+1,j+1,i) )       &
227                                      / 4.0 )**2    +        &
228                                        w(k,j,i)**2 ) *      &
229                              w(k,j,i)
230   
231          ENDDO
232
233!
234!--       sgs-tke
235       CASE ( 4 )
236          DO  k = nzb_s_inner(j,i)+1, pch_index   
237             tend(k,j,i) = tend(k,j,i) -                     &
238                              2.0 * cdc(k,j,i) * lad_s(k,j,i) * &
239                              SQRT( ( ( u(k,j,i)           +    &
240                                        u(k,j,i+1) )            &
241                                      / 2.0 )**2           +    & 
242                                    ( ( v(k,j,i)           +    &
243                                        v(k,j+1,i) )            &
244                                      / 2.0 )**2           +    &
245                                    ( ( w(k,j,i)           +    &
246                                        w(k+1,j,i) )            &
247                                      / 2.0 )**2 )         *    &
248                              e(k,j,i)
249          ENDDO
250
251       CASE DEFAULT
252
253          IF ( myid == 0 )  PRINT*,'+++ pcm:  wrong component: ', &
254                                      component
255          CALL local_stop
256
257       END SELECT
258
259    END SUBROUTINE plant_canopy_model_ij
260
261 END MODULE plant_canopy_model_mod
Note: See TracBrowser for help on using the repository browser.