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

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

last commit documented

  • Property svn:keywords set to Id
File size: 15.4 KB
RevLine 
[138]1 MODULE plant_canopy_model_mod
2
[1036]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!
[1310]17! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]18!--------------------------------------------------------------------------------!
19!
[257]20! Current revisions:
[138]21! -----------------
[1341]22!
23!
[1321]24! Former revisions:
25! -----------------
26! $Id: plant_canopy_model.f90 1341 2014-03-25 19:48:09Z raasch $
27!
[1341]28! 1340 2014-03-25 19:45:13Z kanani
29! REAL constants defined as wp-kind
30!
[1321]31! 1320 2014-03-20 08:40:49Z raasch
[1320]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
[153]39!
[1037]40! 1036 2012-10-22 13:43:42Z raasch
41! code put under GPL (PALM 3.9)
42!
[139]43! 138 2007-11-28 10:03:58Z letzel
44! Initial revision
45!
[138]46! Description:
47! ------------
[153]48! Evaluation of sinks and sources of momentum, heat and scalar concentration
49! due to canopy elements
[138]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
[1320]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
[138]71
[1320]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
[138]81       IMPLICIT NONE
82
[1320]83       INTEGER(iwp) ::  component  !:
84       INTEGER(iwp) ::  i          !:
85       INTEGER(iwp) ::  j          !:
86       INTEGER(iwp) ::  k          !:
[138]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)      +  &
[153]104                                              v(k,j+1,i-1) )     &
[1340]105                                            / 4.0_wp )**2     +  &
[138]106                                          ( ( w(k-1,j,i-1)    +  &
107                                              w(k-1,j,i)      +  &
108                                              w(k,j,i-1)      +  &
109                                              w(k,j,i) )         &
[1340]110                                            / 4.0_wp )**2 )   *  &
[138]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) )       &
[1340]128                                            / 4.0_wp )**2     +  &
[138]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) )         &
[1340]134                                            / 4.0_wp )**2 )   *  &
[138]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) )     &
[1340]152                                            / 4.0_wp )**2     +  &
[138]153                                          ( ( v(k,j,i)        +  &
154                                              v(k,j+1,i)      +  &
155                                              v(k+1,j,i)      +  &
156                                              v(k+1,j+1,i) )     &
[1340]157                                            / 4.0_wp )**2     +  &
[138]158                                              w(k,j,i)**2 )   *  &
159                                    w(k,j,i)
160                   ENDDO
161                ENDDO
162             ENDDO
163
164!
[153]165!--       potential temperature
[138]166          CASE ( 4 )
167             DO  i = nxl, nxr
168                DO  j = nys, nyn
169                   DO  k = nzb_s_inner(j,i)+1, pch_index
[1320]170                      tend(k,j,i) = tend(k,j,i) +                   &
[153]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
[138]184                      tend(k,j,i) = tend(k,j,i) -                     &
[153]185                                    sec(k,j,i) * lad_s(k,j,i) *       &
186                                    SQRT( ( ( u(k,j,i)        +       &
187                                              u(k,j,i+1) )            &
[1340]188                                            / 2.0_wp )**2     +       &
[153]189                                          ( ( v(k,j,i)        +       &
190                                              v(k,j+1,i) )            &
[1340]191                                            / 2.0_wp )**2     +       &
[153]192                                          ( ( w(k-1,j,i)      +       & 
193                                              w(k,j,i) )              &
[1340]194                                            / 2.0_wp )**2 )   *       &
[153]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) -                     &
[1340]207                                    2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * &
[138]208                                    SQRT( ( ( u(k,j,i)              + &
209                                              u(k,j,i+1) )            &
[1340]210                                            / 2.0_wp )**2           + &
[138]211                                          ( ( v(k,j,i)              + &
212                                              v(k,j+1,i) )            &
[1340]213                                            / 2.0_wp )**2           + &
[138]214                                          ( ( w(k,j,i)              + &
215                                              w(k+1,j,i) )            &
[1340]216                                            / 2.0_wp )**2 )         * &
[138]217                                    e(k,j,i)
218                   ENDDO
219                ENDDO
220             ENDDO 
221                       
222          CASE DEFAULT
223
[257]224             WRITE( message_string, * ) 'wrong component: ', component
225             CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 ) 
[138]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
[1320]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
[138]240
[1320]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
[138]250       IMPLICIT NONE
251
[1320]252       INTEGER(iwp) ::  component  !:
253       INTEGER(iwp) ::  i          !:
254       INTEGER(iwp) ::  j          !:
255       INTEGER(iwp) ::  k          !:
[138]256
257!
[1320]258!--    Compute drag for the three velocity components
[142]259       SELECT CASE ( component )
[138]260
261!
[142]262!--       u-component
263       CASE ( 1 )
264          DO  k = nzb_u_inner(j,i)+1, pch_index
[1320]265             tend(k,j,i) = tend(k,j,i) -                     &
[142]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)  +        &
[153]271                                        v(k,j+1,i-1) )       &
[1340]272                                      / 4.0_wp )**2 +        &
[142]273                                    ( ( w(k-1,j,i-1) +       &
274                                        w(k-1,j,i)   +       &
275                                        w(k,j,i-1)   +       &
276                                        w(k,j,i) )           &
[1340]277                                      / 4.0_wp )**2 ) *      &
[142]278                              u(k,j,i)
279          ENDDO
[138]280
281!
[142]282!--       v-component
283       CASE ( 2 )
284          DO  k = nzb_v_inner(j,i)+1, pch_index
[1320]285             tend(k,j,i) = tend(k,j,i) -                     &
[142]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) )         &
[1340]291                                      / 4.0_wp )**2  +       &
[142]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) )           &
[1340]297                                      / 4.0_wp )**2 ) *      &
[142]298                              v(k,j,i)
299          ENDDO
[138]300
301!
[142]302!--       w-component
303       CASE ( 3 )
304          DO  k = nzb_w_inner(j,i)+1, pch_index
[1320]305             tend(k,j,i) = tend(k,j,i) -                     &
[142]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) )       &
[1340]311                                      / 4.0_wp )**2 +        &
[142]312                                    ( ( v(k,j,i)    +        &
313                                        v(k,j+1,i)  +        &
314                                        v(k+1,j,i)  +        &
315                                        v(k+1,j+1,i) )       &
[1340]316                                      / 4.0_wp )**2 +        &
[142]317                                        w(k,j,i)**2 ) *      &
318                              w(k,j,i)
[138]319   
[142]320          ENDDO
[138]321
322!
[153]323!--       potential temperature
324          CASE ( 4 )
325             DO  k = nzb_s_inner(j,i)+1, pch_index
[1320]326                tend(k,j,i) = tend(k,j,i) +                   &
[153]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) )            &
[1340]341                                      / 2.0_wp )**2     +       &
[153]342                                    ( ( v(k,j,i)        +       &
343                                        v(k,j+1,i) )            &
[1340]344                                      / 2.0_wp )**2     +       &
[153]345                                    ( ( w(k-1,j,i)      +       &
346                                        w(k,j,i) )              &
[1340]347                                      / 2.0_wp )**2 )   *       &
[153]348                              ( q(k,j,i) - sls(k,j,i) )
349             ENDDO   
350
351!
[142]352!--       sgs-tke
[153]353       CASE ( 6 )
[142]354          DO  k = nzb_s_inner(j,i)+1, pch_index   
[1320]355             tend(k,j,i) = tend(k,j,i) -                        &
[1340]356                              2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * &
[142]357                              SQRT( ( ( u(k,j,i)           +    &
358                                        u(k,j,i+1) )            &
[1340]359                                      / 2.0_wp )**2        +    & 
[142]360                                    ( ( v(k,j,i)           +    &
361                                        v(k,j+1,i) )            &
[1340]362                                      / 2.0_wp )**2        +    &
[142]363                                    ( ( w(k,j,i)           +    &
364                                        w(k+1,j,i) )            &
[1340]365                                      / 2.0_wp )**2 )      *    &
[142]366                              e(k,j,i)
367          ENDDO
[138]368
[142]369       CASE DEFAULT
[138]370
[257]371          WRITE( message_string, * ) 'wrong component: ', component
372          CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 ) 
[138]373
[142]374       END SELECT
[138]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.