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

Last change on this file since 1340 was 1340, checked in by kanani, 8 years ago

REAL constants defined as wp-kind

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