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

Last change on this file since 1444 was 1341, checked in by kanani, 11 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 15.4 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-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: plant_canopy_model.f90 1341 2014-03-25 19:48:09Z letzel $
27!
28! 1340 2014-03-25 19:45:13Z kanani
29! REAL constants defined as wp-kind
30!
31! 1320 2014-03-20 08:40:49Z raasch
32! ONLY-attribute added to USE-statements,
33! kind-parameters added to all INTEGER and REAL declaration statements,
34! kinds are defined in new module kinds,
35! old module precision_kind is removed,
36! revision history before 2012 removed,
37! comment fields (!:) to be used for variable explanations added to
38! all variable declaration statements
39!
40! 1036 2012-10-22 13:43:42Z raasch
41! code put under GPL (PALM 3.9)
42!
43! 138 2007-11-28 10:03:58Z letzel
44! Initial revision
45!
46! Description:
47! ------------
48! Evaluation of sinks and sources of momentum, heat and scalar concentration
49! due to canopy elements
50!------------------------------------------------------------------------------!
51
52    PRIVATE
53    PUBLIC plant_canopy_model
54
55    INTERFACE plant_canopy_model
56       MODULE PROCEDURE plant_canopy_model
57       MODULE PROCEDURE plant_canopy_model_ij
58    END INTERFACE plant_canopy_model
59
60 CONTAINS
61
62
63!------------------------------------------------------------------------------!
64! Call for all grid points
65!------------------------------------------------------------------------------!
66    SUBROUTINE plant_canopy_model( component )
67
68       USE arrays_3d,                                                          &
69           ONLY:  canopy_heat_flux, cdc, dzw, e, lad_s, lad_u, lad_v, lad_w,   &
70                  q, sec, sls, tend, u, v, w
71
72       USE control_parameters,                                                 &
73           ONLY:  pch_index, message_string
74
75       USE indices,                                                            &
76           ONLY:  nxl, nxlu, nxr, nys, nysv, nyn, nzb_s_inner, nzb_u_inner,    &
77                  nzb_v_inner, nzb_w_inner
78
79       USE kinds
80
81       IMPLICIT NONE
82
83       INTEGER(iwp) ::  component  !:
84       INTEGER(iwp) ::  i          !:
85       INTEGER(iwp) ::  j          !:
86       INTEGER(iwp) ::  k          !:
87 
88!
89!--    Compute drag for the three velocity components and the SGS-TKE
90       SELECT CASE ( component )
91
92!
93!--       u-component
94          CASE ( 1 )
95             DO  i = nxlu, nxr
96                DO  j = nys, nyn
97                   DO  k = nzb_u_inner(j,i)+1, pch_index
98                      tend(k,j,i) = tend(k,j,i) -                &
99                                    cdc(k,j,i) * lad_u(k,j,i) *  &
100                                    SQRT(     u(k,j,i)**2     +  &
101                                          ( ( v(k,j,i-1)      +  &
102                                              v(k,j,i)        +  &
103                                              v(k,j+1,i)      +  &
104                                              v(k,j+1,i-1) )     &
105                                            / 4.0_wp )**2     +  &
106                                          ( ( w(k-1,j,i-1)    +  &
107                                              w(k-1,j,i)      +  &
108                                              w(k,j,i-1)      +  &
109                                              w(k,j,i) )         &
110                                            / 4.0_wp )**2 )   *  &
111                                    u(k,j,i)
112                   ENDDO
113                ENDDO
114             ENDDO
115
116!
117!--       v-component
118          CASE ( 2 )
119             DO  i = nxl, nxr
120                DO  j = nysv, nyn
121                   DO  k = nzb_v_inner(j,i)+1, pch_index
122                      tend(k,j,i) = tend(k,j,i) -                &
123                                    cdc(k,j,i) * lad_v(k,j,i) *  &
124                                    SQRT( ( ( u(k,j-1,i)      +  &
125                                              u(k,j-1,i+1)    +  &
126                                              u(k,j,i)        +  &
127                                              u(k,j,i+1) )       &
128                                            / 4.0_wp )**2     +  &
129                                              v(k,j,i)**2     +  &
130                                          ( ( w(k-1,j-1,i)    +  &
131                                              w(k-1,j,i)      +  &
132                                              w(k,j-1,i)      +  &
133                                              w(k,j,i) )         &
134                                            / 4.0_wp )**2 )   *  &
135                                    v(k,j,i) 
136                   ENDDO
137                ENDDO
138             ENDDO
139
140!
141!--       w-component
142          CASE ( 3 )
143             DO  i = nxl, nxr
144                DO  j = nys, nyn
145                   DO  k = nzb_w_inner(j,i)+1, pch_index
146                      tend(k,j,i) = tend(k,j,i) -                & 
147                                    cdc(k,j,i) * lad_w(k,j,i) *  &
148                                    SQRT( ( ( u(k,j,i)        +  &
149                                              u(k,j,i+1)      +  &
150                                              u(k+1,j,i)      +  &
151                                              u(k+1,j,i+1) )     &
152                                            / 4.0_wp )**2     +  &
153                                          ( ( v(k,j,i)        +  &
154                                              v(k,j+1,i)      +  &
155                                              v(k+1,j,i)      +  &
156                                              v(k+1,j+1,i) )     &
157                                            / 4.0_wp )**2     +  &
158                                              w(k,j,i)**2 )   *  &
159                                    w(k,j,i)
160                   ENDDO
161                ENDDO
162             ENDDO
163
164!
165!--       potential temperature
166          CASE ( 4 )
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                                    ( canopy_heat_flux(k,j,i) -     &
172                                      canopy_heat_flux(k-1,j,i) ) / &
173                                      dzw(k)
174                   ENDDO
175                ENDDO
176             ENDDO
177
178!
179!--       scalar concentration
180          CASE ( 5 )
181             DO  i = nxl, nxr
182                DO  j = nys, nyn
183                   DO  k = nzb_s_inner(j,i)+1, pch_index
184                      tend(k,j,i) = tend(k,j,i) -                     &
185                                    sec(k,j,i) * lad_s(k,j,i) *       &
186                                    SQRT( ( ( u(k,j,i)        +       &
187                                              u(k,j,i+1) )            &
188                                            / 2.0_wp )**2     +       &
189                                          ( ( v(k,j,i)        +       &
190                                              v(k,j+1,i) )            &
191                                            / 2.0_wp )**2     +       &
192                                          ( ( w(k-1,j,i)      +       & 
193                                              w(k,j,i) )              &
194                                            / 2.0_wp )**2 )   *       &
195                                    ( q(k,j,i) - sls(k,j,i) )
196                   ENDDO
197                ENDDO
198             ENDDO
199
200!
201!--       sgs-tke
202          CASE ( 6 )
203             DO  i = nxl, nxr
204                DO  j = nys, nyn
205                   DO  k = nzb_s_inner(j,i)+1, pch_index
206                      tend(k,j,i) = tend(k,j,i) -                     &
207                                    2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * &
208                                    SQRT( ( ( u(k,j,i)              + &
209                                              u(k,j,i+1) )            &
210                                            / 2.0_wp )**2           + &
211                                          ( ( v(k,j,i)              + &
212                                              v(k,j+1,i) )            &
213                                            / 2.0_wp )**2           + &
214                                          ( ( w(k,j,i)              + &
215                                              w(k+1,j,i) )            &
216                                            / 2.0_wp )**2 )         * &
217                                    e(k,j,i)
218                   ENDDO
219                ENDDO
220             ENDDO 
221                       
222          CASE DEFAULT
223
224             WRITE( message_string, * ) 'wrong component: ', component
225             CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 ) 
226
227       END SELECT
228
229    END SUBROUTINE plant_canopy_model
230
231
232!------------------------------------------------------------------------------!
233! Call for grid point i,j
234!------------------------------------------------------------------------------!
235    SUBROUTINE plant_canopy_model_ij( i, j, component )
236
237       USE arrays_3d,                                                          &
238           ONLY:  canopy_heat_flux, cdc, dzw, e, lad_s, lad_u, lad_v, lad_w,   &
239                  q, sec, sls, tend, u, v, w
240
241       USE control_parameters,                                                 &
242           ONLY:  pch_index, message_string
243
244       USE indices,                                                            &
245           ONLY:  nxl, nxlu, nxr, nys, nysv, nyn, nzb_s_inner, nzb_u_inner,    &
246                  nzb_v_inner, nzb_w_inner
247
248       USE kinds
249
250       IMPLICIT NONE
251
252       INTEGER(iwp) ::  component  !:
253       INTEGER(iwp) ::  i          !:
254       INTEGER(iwp) ::  j          !:
255       INTEGER(iwp) ::  k          !:
256
257!
258!--    Compute drag for the three velocity components
259       SELECT CASE ( component )
260
261!
262!--       u-component
263       CASE ( 1 )
264          DO  k = nzb_u_inner(j,i)+1, pch_index
265             tend(k,j,i) = tend(k,j,i) -                     &
266                              cdc(k,j,i) * lad_u(k,j,i) *    &   
267                              SQRT(     u(k,j,i)**2 +        &
268                                    ( ( v(k,j,i-1)  +        &
269                                        v(k,j,i)    +        &
270                                        v(k,j+1,i)  +        &
271                                        v(k,j+1,i-1) )       &
272                                      / 4.0_wp )**2 +        &
273                                    ( ( w(k-1,j,i-1) +       &
274                                        w(k-1,j,i)   +       &
275                                        w(k,j,i-1)   +       &
276                                        w(k,j,i) )           &
277                                      / 4.0_wp )**2 ) *      &
278                              u(k,j,i)
279          ENDDO
280
281!
282!--       v-component
283       CASE ( 2 )
284          DO  k = nzb_v_inner(j,i)+1, pch_index
285             tend(k,j,i) = tend(k,j,i) -                     &
286                              cdc(k,j,i) * lad_v(k,j,i) *    &
287                              SQRT( ( ( u(k,j-1,i)   +       &
288                                        u(k,j-1,i+1) +       &
289                                        u(k,j,i)     +       &
290                                        u(k,j,i+1) )         &
291                                      / 4.0_wp )**2  +       &
292                                        v(k,j,i)**2  +       &
293                                    ( ( w(k-1,j-1,i) +       &
294                                        w(k-1,j,i)   +       &
295                                        w(k,j-1,i)   +       &
296                                        w(k,j,i) )           &
297                                      / 4.0_wp )**2 ) *      &
298                              v(k,j,i)
299          ENDDO
300
301!
302!--       w-component
303       CASE ( 3 )
304          DO  k = nzb_w_inner(j,i)+1, pch_index
305             tend(k,j,i) = tend(k,j,i) -                     &
306                              cdc(k,j,i) * lad_w(k,j,i) *    & 
307                              SQRT( ( ( u(k,j,i)    +        & 
308                                        u(k,j,i+1)  +        &
309                                        u(k+1,j,i)  +        &
310                                        u(k+1,j,i+1) )       &
311                                      / 4.0_wp )**2 +        &
312                                    ( ( v(k,j,i)    +        &
313                                        v(k,j+1,i)  +        &
314                                        v(k+1,j,i)  +        &
315                                        v(k+1,j+1,i) )       &
316                                      / 4.0_wp )**2 +        &
317                                        w(k,j,i)**2 ) *      &
318                              w(k,j,i)
319   
320          ENDDO
321
322!
323!--       potential temperature
324          CASE ( 4 )
325             DO  k = nzb_s_inner(j,i)+1, pch_index
326                tend(k,j,i) = tend(k,j,i) +                   &
327                              ( canopy_heat_flux(k,j,i) -     &
328                                canopy_heat_flux(k-1,j,i) ) / &
329                                dzw(k)
330             ENDDO
331
332
333!
334!--       scalar concentration
335          CASE ( 5 )
336             DO  k = nzb_s_inner(j,i)+1, pch_index
337                tend(k,j,i) = tend(k,j,i) -                     &
338                              sec(k,j,i) * lad_s(k,j,i) *       &
339                              SQRT( ( ( u(k,j,i)        +       &
340                                        u(k,j,i+1) )            &
341                                      / 2.0_wp )**2     +       &
342                                    ( ( v(k,j,i)        +       &
343                                        v(k,j+1,i) )            &
344                                      / 2.0_wp )**2     +       &
345                                    ( ( w(k-1,j,i)      +       &
346                                        w(k,j,i) )              &
347                                      / 2.0_wp )**2 )   *       &
348                              ( q(k,j,i) - sls(k,j,i) )
349             ENDDO   
350
351!
352!--       sgs-tke
353       CASE ( 6 )
354          DO  k = nzb_s_inner(j,i)+1, pch_index   
355             tend(k,j,i) = tend(k,j,i) -                        &
356                              2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * &
357                              SQRT( ( ( u(k,j,i)           +    &
358                                        u(k,j,i+1) )            &
359                                      / 2.0_wp )**2        +    & 
360                                    ( ( v(k,j,i)           +    &
361                                        v(k,j+1,i) )            &
362                                      / 2.0_wp )**2        +    &
363                                    ( ( w(k,j,i)           +    &
364                                        w(k+1,j,i) )            &
365                                      / 2.0_wp )**2 )      *    &
366                              e(k,j,i)
367          ENDDO
368
369       CASE DEFAULT
370
371          WRITE( message_string, * ) 'wrong component: ', component
372          CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 ) 
373
374       END SELECT
375
376    END SUBROUTINE plant_canopy_model_ij
377
378 END MODULE plant_canopy_model_mod
Note: See TracBrowser for help on using the repository browser.