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

Last change on this file since 2636 was 2233, checked in by suehring, 7 years ago

last commit documented

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