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

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

REAL constants defined as wp-kind

  • 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! -----------------
[1340]22! REAL constants defined as wp-kind
[1321]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
[1320]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
[153]36!
[1037]37! 1036 2012-10-22 13:43:42Z raasch
38! code put under GPL (PALM 3.9)
39!
[139]40! 138 2007-11-28 10:03:58Z letzel
41! Initial revision
42!
[138]43! Description:
44! ------------
[153]45! Evaluation of sinks and sources of momentum, heat and scalar concentration
46! due to canopy elements
[138]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
[1320]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
[138]68
[1320]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
[138]78       IMPLICIT NONE
79
[1320]80       INTEGER(iwp) ::  component  !:
81       INTEGER(iwp) ::  i          !:
82       INTEGER(iwp) ::  j          !:
83       INTEGER(iwp) ::  k          !:
[138]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)      +  &
[153]101                                              v(k,j+1,i-1) )     &
[1340]102                                            / 4.0_wp )**2     +  &
[138]103                                          ( ( w(k-1,j,i-1)    +  &
104                                              w(k-1,j,i)      +  &
105                                              w(k,j,i-1)      +  &
106                                              w(k,j,i) )         &
[1340]107                                            / 4.0_wp )**2 )   *  &
[138]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) )       &
[1340]125                                            / 4.0_wp )**2     +  &
[138]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) )         &
[1340]131                                            / 4.0_wp )**2 )   *  &
[138]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) )     &
[1340]149                                            / 4.0_wp )**2     +  &
[138]150                                          ( ( v(k,j,i)        +  &
151                                              v(k,j+1,i)      +  &
152                                              v(k+1,j,i)      +  &
153                                              v(k+1,j+1,i) )     &
[1340]154                                            / 4.0_wp )**2     +  &
[138]155                                              w(k,j,i)**2 )   *  &
156                                    w(k,j,i)
157                   ENDDO
158                ENDDO
159             ENDDO
160
161!
[153]162!--       potential temperature
[138]163          CASE ( 4 )
164             DO  i = nxl, nxr
165                DO  j = nys, nyn
166                   DO  k = nzb_s_inner(j,i)+1, pch_index
[1320]167                      tend(k,j,i) = tend(k,j,i) +                   &
[153]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
[138]181                      tend(k,j,i) = tend(k,j,i) -                     &
[153]182                                    sec(k,j,i) * lad_s(k,j,i) *       &
183                                    SQRT( ( ( u(k,j,i)        +       &
184                                              u(k,j,i+1) )            &
[1340]185                                            / 2.0_wp )**2     +       &
[153]186                                          ( ( v(k,j,i)        +       &
187                                              v(k,j+1,i) )            &
[1340]188                                            / 2.0_wp )**2     +       &
[153]189                                          ( ( w(k-1,j,i)      +       & 
190                                              w(k,j,i) )              &
[1340]191                                            / 2.0_wp )**2 )   *       &
[153]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) -                     &
[1340]204                                    2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * &
[138]205                                    SQRT( ( ( u(k,j,i)              + &
206                                              u(k,j,i+1) )            &
[1340]207                                            / 2.0_wp )**2           + &
[138]208                                          ( ( v(k,j,i)              + &
209                                              v(k,j+1,i) )            &
[1340]210                                            / 2.0_wp )**2           + &
[138]211                                          ( ( w(k,j,i)              + &
212                                              w(k+1,j,i) )            &
[1340]213                                            / 2.0_wp )**2 )         * &
[138]214                                    e(k,j,i)
215                   ENDDO
216                ENDDO
217             ENDDO 
218                       
219          CASE DEFAULT
220
[257]221             WRITE( message_string, * ) 'wrong component: ', component
222             CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 ) 
[138]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
[1320]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
[138]237
[1320]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
[138]247       IMPLICIT NONE
248
[1320]249       INTEGER(iwp) ::  component  !:
250       INTEGER(iwp) ::  i          !:
251       INTEGER(iwp) ::  j          !:
252       INTEGER(iwp) ::  k          !:
[138]253
254!
[1320]255!--    Compute drag for the three velocity components
[142]256       SELECT CASE ( component )
[138]257
258!
[142]259!--       u-component
260       CASE ( 1 )
261          DO  k = nzb_u_inner(j,i)+1, pch_index
[1320]262             tend(k,j,i) = tend(k,j,i) -                     &
[142]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)  +        &
[153]268                                        v(k,j+1,i-1) )       &
[1340]269                                      / 4.0_wp )**2 +        &
[142]270                                    ( ( w(k-1,j,i-1) +       &
271                                        w(k-1,j,i)   +       &
272                                        w(k,j,i-1)   +       &
273                                        w(k,j,i) )           &
[1340]274                                      / 4.0_wp )**2 ) *      &
[142]275                              u(k,j,i)
276          ENDDO
[138]277
278!
[142]279!--       v-component
280       CASE ( 2 )
281          DO  k = nzb_v_inner(j,i)+1, pch_index
[1320]282             tend(k,j,i) = tend(k,j,i) -                     &
[142]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) )         &
[1340]288                                      / 4.0_wp )**2  +       &
[142]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) )           &
[1340]294                                      / 4.0_wp )**2 ) *      &
[142]295                              v(k,j,i)
296          ENDDO
[138]297
298!
[142]299!--       w-component
300       CASE ( 3 )
301          DO  k = nzb_w_inner(j,i)+1, pch_index
[1320]302             tend(k,j,i) = tend(k,j,i) -                     &
[142]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) )       &
[1340]308                                      / 4.0_wp )**2 +        &
[142]309                                    ( ( v(k,j,i)    +        &
310                                        v(k,j+1,i)  +        &
311                                        v(k+1,j,i)  +        &
312                                        v(k+1,j+1,i) )       &
[1340]313                                      / 4.0_wp )**2 +        &
[142]314                                        w(k,j,i)**2 ) *      &
315                              w(k,j,i)
[138]316   
[142]317          ENDDO
[138]318
319!
[153]320!--       potential temperature
321          CASE ( 4 )
322             DO  k = nzb_s_inner(j,i)+1, pch_index
[1320]323                tend(k,j,i) = tend(k,j,i) +                   &
[153]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) )            &
[1340]338                                      / 2.0_wp )**2     +       &
[153]339                                    ( ( v(k,j,i)        +       &
340                                        v(k,j+1,i) )            &
[1340]341                                      / 2.0_wp )**2     +       &
[153]342                                    ( ( w(k-1,j,i)      +       &
343                                        w(k,j,i) )              &
[1340]344                                      / 2.0_wp )**2 )   *       &
[153]345                              ( q(k,j,i) - sls(k,j,i) )
346             ENDDO   
347
348!
[142]349!--       sgs-tke
[153]350       CASE ( 6 )
[142]351          DO  k = nzb_s_inner(j,i)+1, pch_index   
[1320]352             tend(k,j,i) = tend(k,j,i) -                        &
[1340]353                              2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * &
[142]354                              SQRT( ( ( u(k,j,i)           +    &
355                                        u(k,j,i+1) )            &
[1340]356                                      / 2.0_wp )**2        +    & 
[142]357                                    ( ( v(k,j,i)           +    &
358                                        v(k,j+1,i) )            &
[1340]359                                      / 2.0_wp )**2        +    &
[142]360                                    ( ( w(k,j,i)           +    &
361                                        w(k+1,j,i) )            &
[1340]362                                      / 2.0_wp )**2 )      *    &
[142]363                              e(k,j,i)
364          ENDDO
[138]365
[142]366       CASE DEFAULT
[138]367
[257]368          WRITE( message_string, * ) 'wrong component: ', component
369          CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 ) 
[138]370
[142]371       END SELECT
[138]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.