source: palm/trunk/SOURCE/lpm_init_sgs_tke.f90 @ 3436

Last change on this file since 3436 was 3241, checked in by raasch, 6 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

  • Property svn:keywords set to Id
File size: 11.6 KB
RevLine 
[1682]1!> @file lpm_init_sgs_tke.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[2000]5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
[1036]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!
[2718]17! Copyright 1997-2018 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[849]20! Current revisions:
21! ------------------
[1360]22!
[2233]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: lpm_init_sgs_tke.f90 3241 2018-09-12 15:02:00Z gronemeier $
[3241]27! unused variables removed
28!
29! 2718 2018-01-02 08:49:38Z maronga
[2716]30! Corrected "Former revisions" section
31!
32! 2696 2017-12-14 17:12:51Z kanani
33! Change in file header (GPL part)
[1321]34!
[2716]35! 2233 2017-05-30 18:08:54Z suehring
36!
[2233]37! 2232 2017-05-30 17:47:52Z suehring
38! Adjustments according to new topography realization
39!
[2001]40! 2000 2016-08-20 18:09:15Z knoop
41! Forced header and separation lines into 80 columns
42!
[1930]43! 1929 2016-06-09 16:25:25Z suehring
44! sgs_wfu_par, sgs_wfv_par and sgs_wfw_par are replaced by sgs_wf_par
45!
[1823]46! 1822 2016-04-07 07:49:42Z hoffmann
47! Unused variables removed.
48!
[1683]49! 1682 2015-10-07 23:56:08Z knoop
50! Code annotations made doxygen readable
51!
[1360]52! 1359 2014-04-11 17:15:14Z hoffmann
53! New particle structure integrated.
54! Kind definition added to all floating point numbers.
55!
[1321]56! 1320 2014-03-20 08:40:49Z raasch
[1320]57! ONLY-attribute added to USE-statements,
58! kind-parameters added to all INTEGER and REAL declaration statements,
59! kinds are defined in new module kinds,
60! comment fields (!:) to be used for variable explanations added to
61! all variable declaration statements
[849]62!
[1037]63! 1036 2012-10-22 13:43:42Z raasch
64! code put under GPL (PALM 3.9)
65!
[850]66! 849 2012-03-15 10:35:09Z raasch
67! initial revision (former part of advec_particles)
[849]68!
[850]69!
[849]70! Description:
71! ------------
[1682]72!> Calculates quantities required for considering the SGS velocity fluctuations
73!> in the particle transport by a stochastic approach. The respective
74!> quantities are: SGS-TKE gradients and horizontally averaged profiles of the
75!> SGS TKE and the resolved-scale velocity variances.
[849]76!------------------------------------------------------------------------------!
[1682]77 SUBROUTINE lpm_init_sgs_tke
78 
[849]79
[1320]80    USE arrays_3d,                                                             &
81        ONLY:  de_dx, de_dy, de_dz, diss, e, u, v, w, zu
82
83    USE grid_variables,                                                        &
84        ONLY:  ddx, ddy
85
86    USE indices,                                                               &
[3241]87        ONLY:  nbgp, ngp_2dh_outer, nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
[1320]88
89    USE kinds
90
91    USE particle_attributes,                                                   &
[1929]92        ONLY:  sgs_wf_part
[1320]93
[849]94    USE pegrid
95
[1320]96    USE statistics,                                                            &
97        ONLY:  flow_statistics_called, hom, sums, sums_l
98
[2232]99    USE surface_mod,                                                           &
100        ONLY:  bc_h
101
[849]102    IMPLICIT NONE
103
[2232]104    INTEGER(iwp) ::  i      !< index variable along x
105    INTEGER(iwp) ::  j      !< index variable along y
106    INTEGER(iwp) ::  k      !< index variable along z
107    INTEGER(iwp) ::  m      !< running index for the surface elements
[849]108
[2232]109    REAL(wp) ::  flag1      !< flag to mask topography
[849]110
111!
112!-- TKE gradient along x and y
113    DO  i = nxl, nxr
114       DO  j = nys, nyn
115          DO  k = nzb, nzt+1
116
[2232]117             IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 )  .AND.               &
118                        BTEST( wall_flags_0(k,j,i), 0   )  .AND.               &
119                        BTEST( wall_flags_0(k,j,i+1), 0 ) )                    &
[849]120             THEN
[1929]121                de_dx(k,j,i) = 2.0_wp * sgs_wf_part *                          &
[1359]122                               ( e(k,j,i+1) - e(k,j,i) ) * ddx
[2232]123             ELSEIF ( BTEST( wall_flags_0(k,j,i-1), 0 )  .AND.                 &
124                      BTEST( wall_flags_0(k,j,i), 0   )  .AND.                 &
125                .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) )                      &
[849]126             THEN
[1929]127                de_dx(k,j,i) = 2.0_wp * sgs_wf_part *                          &
[1359]128                               ( e(k,j,i) - e(k,j,i-1) ) * ddx
[2232]129             ELSEIF ( .NOT. BTEST( wall_flags_0(k,j,i), 22   )  .AND.          &
130                      .NOT. BTEST( wall_flags_0(k,j,i+1), 22 ) )               &   
[849]131             THEN
[1359]132                de_dx(k,j,i) = 0.0_wp
[2232]133             ELSEIF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 22 )  .AND.          &
134                      .NOT. BTEST( wall_flags_0(k,j,i), 22   ) )               &
[849]135             THEN
[1359]136                de_dx(k,j,i) = 0.0_wp
[849]137             ELSE
[1929]138                de_dx(k,j,i) = sgs_wf_part * ( e(k,j,i+1) - e(k,j,i-1) ) * ddx
[849]139             ENDIF
140
[2232]141             IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 )  .AND.               &
142                        BTEST( wall_flags_0(k,j,i), 0   )  .AND.               &
143                        BTEST( wall_flags_0(k,j+1,i), 0 ) )                    &
[849]144             THEN
[1929]145                de_dy(k,j,i) = 2.0_wp * sgs_wf_part *                          &
[1359]146                               ( e(k,j+1,i) - e(k,j,i) ) * ddy
[2232]147             ELSEIF ( BTEST( wall_flags_0(k,j-1,i), 0 )  .AND.                 &
148                      BTEST( wall_flags_0(k,j,i), 0   )  .AND.                 &
149                .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) )                      &
[849]150             THEN
[1929]151                de_dy(k,j,i) = 2.0_wp * sgs_wf_part *                          &
[1359]152                               ( e(k,j,i) - e(k,j-1,i) ) * ddy
[2232]153             ELSEIF ( .NOT. BTEST( wall_flags_0(k,j,i), 22   )  .AND.          &
154                      .NOT. BTEST( wall_flags_0(k,j+1,i), 22 ) )               &   
[849]155             THEN
[1359]156                de_dy(k,j,i) = 0.0_wp
[2232]157             ELSEIF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 22 )  .AND.          &
158                      .NOT. BTEST( wall_flags_0(k,j,i), 22   ) )               &
[849]159             THEN
[1359]160                de_dy(k,j,i) = 0.0_wp
[849]161             ELSE
[1929]162                de_dy(k,j,i) = sgs_wf_part * ( e(k,j+1,i) - e(k,j-1,i) ) * ddy
[849]163             ENDIF
164
165          ENDDO
166       ENDDO
167    ENDDO
168
169!
[2232]170!-- TKE gradient along z at topograhy and  including bottom and top boundary conditions
[849]171    DO  i = nxl, nxr
172       DO  j = nys, nyn
[2232]173          DO  k = nzb+1, nzt-1
174!
175!--          Flag to mask topography
176             flag1 = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0  ) )
[849]177
[2232]178             de_dz(k,j,i) = 2.0_wp * sgs_wf_part *                             &
179                           ( e(k+1,j,i) - e(k-1,j,i) ) / ( zu(k+1) - zu(k-1) ) &
180                                                 * flag1 
[849]181          ENDDO
[2232]182!
183!--       upward-facing surfaces
184          DO  m = bc_h(0)%start_index(j,i), bc_h(0)%end_index(j,i)
185             k            = bc_h(0)%k(m)
186             de_dz(k,j,i) = 2.0_wp * sgs_wf_part *                             &
187                           ( e(k+1,j,i) - e(k,j,i)   ) / ( zu(k+1) - zu(k) )
188          ENDDO
189!
190!--       downward-facing surfaces
191          DO  m = bc_h(1)%start_index(j,i), bc_h(1)%end_index(j,i)
192             k            = bc_h(1)%k(m)
193             de_dz(k,j,i) = 2.0_wp * sgs_wf_part *                             &
194                           ( e(k,j,i) - e(k-1,j,i)   ) / ( zu(k) - zu(k-1) )
195          ENDDO
[849]196
[2232]197          de_dz(nzb,j,i)   = 0.0_wp
[1359]198          de_dz(nzt,j,i)   = 0.0_wp
199          de_dz(nzt+1,j,i) = 0.0_wp
[849]200       ENDDO
201    ENDDO
202!
[2232]203!-- Ghost point exchange
[849]204    CALL exchange_horiz( de_dx, nbgp )
205    CALL exchange_horiz( de_dy, nbgp )
206    CALL exchange_horiz( de_dz, nbgp )
207    CALL exchange_horiz( diss, nbgp  )
208
209
210!
211!-- Calculate the horizontally averaged profiles of SGS TKE and resolved
212!-- velocity variances (they may have been already calculated in routine
213!-- flow_statistics).
214    IF ( .NOT. flow_statistics_called )  THEN
215
216!
217!--    First calculate horizontally averaged profiles of the horizontal
218!--    velocities.
[1359]219       sums_l(:,1,0) = 0.0_wp
220       sums_l(:,2,0) = 0.0_wp
[849]221
222       DO  i = nxl, nxr
223          DO  j =  nys, nyn
[2232]224             DO  k = nzb, nzt+1
225!
[3241]226!--             Flag indicating vicinity of wall
[2232]227                flag1 = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 24 ) )
228
229                sums_l(k,1,0)  = sums_l(k,1,0)  + u(k,j,i) * flag1
230                sums_l(k,2,0)  = sums_l(k,2,0)  + v(k,j,i) * flag1
[849]231             ENDDO
232          ENDDO
233       ENDDO
234
235#if defined( __parallel )
236!
237!--    Compute total sum from local sums
238       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
239       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, &
240                           MPI_REAL, MPI_SUM, comm2d, ierr )
241       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
242       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, &
243                              MPI_REAL, MPI_SUM, comm2d, ierr )
244#else
245       sums(:,1) = sums_l(:,1,0)
246       sums(:,2) = sums_l(:,2,0)
247#endif
248
249!
250!--    Final values are obtained by division by the total number of grid
251!--    points used for the summation.
252       hom(:,1,1,0) = sums(:,1) / ngp_2dh_outer(:,0)   ! u
253       hom(:,1,2,0) = sums(:,2) / ngp_2dh_outer(:,0)   ! v
254
255!
256!--    Now calculate the profiles of SGS TKE and the resolved-scale
257!--    velocity variances
[1359]258       sums_l(:,8,0)  = 0.0_wp
259       sums_l(:,30,0) = 0.0_wp
260       sums_l(:,31,0) = 0.0_wp
261       sums_l(:,32,0) = 0.0_wp
[849]262       DO  i = nxl, nxr
263          DO  j = nys, nyn
[2232]264             DO  k = nzb, nzt+1
265!
[3241]266!--             Flag indicating vicinity of wall
[2232]267                flag1 = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 24 ) )
268
269                sums_l(k,8,0)  = sums_l(k,8,0)  + e(k,j,i)                       * flag1
270                sums_l(k,30,0) = sums_l(k,30,0) + ( u(k,j,i) - hom(k,1,1,0) )**2 * flag1
271                sums_l(k,31,0) = sums_l(k,31,0) + ( v(k,j,i) - hom(k,1,2,0) )**2 * flag1
272                sums_l(k,32,0) = sums_l(k,32,0) + w(k,j,i)**2                    * flag1
[849]273             ENDDO
274          ENDDO
275       ENDDO
276
277#if defined( __parallel )
278!
279!--    Compute total sum from local sums
280       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
281       CALL MPI_ALLREDUCE( sums_l(nzb,8,0), sums(nzb,8), nzt+2-nzb, &
282                           MPI_REAL, MPI_SUM, comm2d, ierr )
283       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
284       CALL MPI_ALLREDUCE( sums_l(nzb,30,0), sums(nzb,30), nzt+2-nzb, &
285                           MPI_REAL, MPI_SUM, comm2d, ierr )
286       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
287       CALL MPI_ALLREDUCE( sums_l(nzb,31,0), sums(nzb,31), nzt+2-nzb, &
288                           MPI_REAL, MPI_SUM, comm2d, ierr )
289       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
290       CALL MPI_ALLREDUCE( sums_l(nzb,32,0), sums(nzb,32), nzt+2-nzb, &
291                           MPI_REAL, MPI_SUM, comm2d, ierr )
292
293#else
294       sums(:,8)  = sums_l(:,8,0)
295       sums(:,30) = sums_l(:,30,0)
296       sums(:,31) = sums_l(:,31,0)
297       sums(:,32) = sums_l(:,32,0)
298#endif
299
300!
301!--    Final values are obtained by division by the total number of grid
302!--    points used for the summation.
303       hom(:,1,8,0)  = sums(:,8)  / ngp_2dh_outer(:,0)   ! e
304       hom(:,1,30,0) = sums(:,30) / ngp_2dh_outer(:,0)   ! u*2
305       hom(:,1,31,0) = sums(:,31) / ngp_2dh_outer(:,0)   ! v*2
306       hom(:,1,32,0) = sums(:,32) / ngp_2dh_outer(:,0)   ! w*2
307
308    ENDIF
309
310 END SUBROUTINE lpm_init_sgs_tke
Note: See TracBrowser for help on using the repository browser.