source: palm/trunk/SOURCE/boundary_conds.f90 @ 4275

Last change on this file since 4275 was 4272, checked in by schwenkel, 4 years ago

further modularization of boundary conditions: moving boundary conditions to their modules

  • Property svn:keywords set to Id
File size: 30.4 KB
Line 
1!> @file boundary_conds.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
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.
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-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: boundary_conds.f90 4272 2019-10-23 15:18:57Z schwenkel $
27! Further modularization of boundary conditions: moved boundary conditions to
28! respective modules
29!
30! 4268 2019-10-17 11:29:38Z schwenkel
31! Removing bulk cloud variables to respective module
32!
33! 4182 2019-08-22 15:20:23Z scharf
34! Corrected "Former revisions" section
35!
36! 4102 2019-07-17 16:00:03Z suehring
37! - Revise setting for boundary conditions at horizontal walls, use the offset
38!   index that belongs to the data structure instead of pre-calculating
39!   the offset index for each facing.
40! - Set boundary conditions for bulk-cloud quantities also at downward facing
41!   surfaces
42!
43! 4087 2019-07-11 11:35:23Z gronemeier
44! Add comment
45!
46! 4086 2019-07-11 05:55:44Z gronemeier
47! Bugfix: use constant-flux layer condition for e in all rans modes
48!
49! 3879 2019-04-08 20:25:23Z knoop
50! Bugfix, do not set boundary conditions for potential temperature in neutral
51! runs.
52!
53! 3655 2019-01-07 16:51:22Z knoop
54! OpenACC port for SPEC
55!
56! Revision 1.1  1997/09/12 06:21:34  raasch
57! Initial revision
58!
59!
60! Description:
61! ------------
62!> Boundary conditions for the prognostic quantities.
63!> One additional bottom boundary condition is applied for the TKE (=(u*)**2)
64!> in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
65!> handled in routine exchange_horiz. Pressure boundary conditions are
66!> explicitly set in routines pres, poisfft, poismg and sor.
67!------------------------------------------------------------------------------!
68 SUBROUTINE boundary_conds
69 
70
71    USE arrays_3d,                                                             &
72        ONLY:  c_u, c_u_m, c_u_m_l, c_v, c_v_m, c_v_m_l, c_w, c_w_m, c_w_m_l,  &
73               dzu, pt, pt_init, pt_p, q,                                      &
74               q_p, s, s_p, u, u_init, u_m_l, u_m_n,                 &
75               u_m_r, u_m_s, u_p, v, v_init, v_m_l, v_m_n, v_m_r, v_m_s, v_p,  &
76               w, w_p, w_m_l, w_m_n, w_m_r, w_m_s
77
78    USE control_parameters,                                                    &
79        ONLY:  bc_dirichlet_l,                                  &
80               bc_dirichlet_s, bc_radiation_l, bc_radiation_n, bc_radiation_r, &
81               bc_radiation_s, bc_pt_t_val, bc_q_t_val, bc_s_t_val,            &
82               child_domain, coupling_mode, dt_3d,                             &
83               humidity, ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b,        &
84               ibc_s_t, ibc_uv_b, ibc_uv_t, intermediate_timestep_count,       &
85               nesting_offline, neutral, nudging, passive_scalar,  &
86               tsc, use_cmax
87
88    USE grid_variables,                                                        &
89        ONLY:  ddx, ddy, dx, dy
90
91    USE indices,                                                               &
92        ONLY:  nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, nzt
93
94    USE kinds
95
96    USE pegrid
97
98    USE pmc_interface,                                                         &
99        ONLY : nesting_mode
100 
101    USE surface_mod,                                                           &
102        ONLY :  bc_h
103
104    USE turbulence_closure_mod,                                                &
105        ONLY:  tcm_boundary_conds
106
107    IMPLICIT NONE
108
109    INTEGER(iwp) ::  i  !< grid index x direction
110    INTEGER(iwp) ::  j  !< grid index y direction
111    INTEGER(iwp) ::  k  !< grid index z direction
112    INTEGER(iwp) ::  l  !< running index boundary type, for up- and downward-facing walls
113    INTEGER(iwp) ::  m  !< running index surface elements
114
115    REAL(wp)    ::  c_max !< maximum phase velocity allowed by CFL criterion, used for outflow boundary condition
116    REAL(wp)    ::  denom !< horizontal gradient of velocity component normal to the outflow boundary
117
118
119!
120!-- Bottom boundary
121    IF ( ibc_uv_b == 1 )  THEN
122       u_p(nzb,:,:) = u_p(nzb+1,:,:)
123       v_p(nzb,:,:) = v_p(nzb+1,:,:)
124    ENDIF
125!
126!-- Set zero vertical velocity at topography top (l=0), or bottom (l=1) in case
127!-- of downward-facing surfaces.
128    DO  l = 0, 1
129       !$OMP PARALLEL DO PRIVATE( i, j, k )
130       !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
131       !$ACC PRESENT(bc_h, w_p)
132       DO  m = 1, bc_h(l)%ns
133          i = bc_h(l)%i(m)
134          j = bc_h(l)%j(m)
135          k = bc_h(l)%k(m)
136          w_p(k+bc_h(l)%koff,j,i) = 0.0_wp
137       ENDDO
138    ENDDO
139
140!
141!-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings.
142    IF ( ibc_uv_t == 0 )  THEN
143        !$ACC KERNELS PRESENT(u_p, v_p, u_init, v_init)
144        u_p(nzt+1,:,:) = u_init(nzt+1)
145        v_p(nzt+1,:,:) = v_init(nzt+1)
146        !$ACC END KERNELS
147    ELSEIF ( ibc_uv_t == 1 )  THEN
148        u_p(nzt+1,:,:) = u_p(nzt,:,:)
149        v_p(nzt+1,:,:) = v_p(nzt,:,:)
150    ENDIF
151
152!
153!-- Vertical nesting: Vertical velocity not zero at the top of the fine grid
154    IF (  .NOT.  child_domain  .AND.  .NOT.  nesting_offline  .AND.            &
155                 TRIM(coupling_mode) /= 'vnested_fine' )  THEN
156       !$ACC KERNELS PRESENT(w_p)
157       w_p(nzt:nzt+1,:,:) = 0.0_wp  !< nzt is not a prognostic level (but cf. pres)
158       !$ACC END KERNELS
159    ENDIF
160
161!
162!-- Temperature at bottom and top boundary.
163!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by
164!-- the sea surface temperature of the coupled ocean model.
165!-- Dirichlet
166    IF ( .NOT. neutral )  THEN
167       IF ( ibc_pt_b == 0 )  THEN
168          DO  l = 0, 1
169             !$OMP PARALLEL DO PRIVATE( i, j, k )
170             DO  m = 1, bc_h(l)%ns
171                i = bc_h(l)%i(m)
172                j = bc_h(l)%j(m)
173                k = bc_h(l)%k(m)
174                pt_p(k+bc_h(l)%koff,j,i) = pt(k+bc_h(l)%koff,j,i)
175             ENDDO
176          ENDDO
177!     
178!--    Neumann, zero-gradient
179       ELSEIF ( ibc_pt_b == 1 )  THEN
180          DO  l = 0, 1
181             !$OMP PARALLEL DO PRIVATE( i, j, k )
182             !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
183             !$ACC PRESENT(bc_h, pt_p)
184             DO  m = 1, bc_h(l)%ns
185                i = bc_h(l)%i(m)
186                j = bc_h(l)%j(m)
187                k = bc_h(l)%k(m)
188                pt_p(k+bc_h(l)%koff,j,i) = pt_p(k,j,i)
189             ENDDO
190          ENDDO
191       ENDIF
192       
193!     
194!--    Temperature at top boundary
195       IF ( ibc_pt_t == 0 )  THEN
196           pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
197!     
198!--        In case of nudging adjust top boundary to pt which is
199!--        read in from NUDGING-DATA
200           IF ( nudging )  THEN
201              pt_p(nzt+1,:,:) = pt_init(nzt+1)
202           ENDIF
203       ELSEIF ( ibc_pt_t == 1 )  THEN
204           pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
205       ELSEIF ( ibc_pt_t == 2 )  THEN
206           !$ACC KERNELS PRESENT(pt_p, dzu)
207           pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
208           !$ACC END KERNELS
209       ENDIF
210    ENDIF
211!
212!-- Boundary conditions for total water content,
213!-- bottom and top boundary (see also temperature)
214    IF ( humidity )  THEN
215!
216!--    Surface conditions for constant_humidity_flux
217!--    Run loop over all non-natural and natural walls. Note, in wall-datatype
218!--    the k coordinate belongs to the atmospheric grid point, therefore, set
219!--    q_p at k-1
220       IF ( ibc_q_b == 0 ) THEN
221
222          DO  l = 0, 1
223             !$OMP PARALLEL DO PRIVATE( i, j, k )
224             DO  m = 1, bc_h(l)%ns
225                i = bc_h(l)%i(m)
226                j = bc_h(l)%j(m)
227                k = bc_h(l)%k(m)
228                q_p(k+bc_h(l)%koff,j,i) = q(k+bc_h(l)%koff,j,i)
229             ENDDO
230          ENDDO
231
232       ELSE
233
234          DO  l = 0, 1
235             !$OMP PARALLEL DO PRIVATE( i, j, k )
236             DO  m = 1, bc_h(l)%ns
237                i = bc_h(l)%i(m)
238                j = bc_h(l)%j(m)
239                k = bc_h(l)%k(m)
240                q_p(k+bc_h(l)%koff,j,i) = q_p(k,j,i)
241             ENDDO
242          ENDDO
243       ENDIF
244!
245!--    Top boundary
246       IF ( ibc_q_t == 0 ) THEN
247          q_p(nzt+1,:,:) = q(nzt+1,:,:)
248       ELSEIF ( ibc_q_t == 1 ) THEN
249          q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
250       ENDIF
251    ENDIF
252!
253!-- Boundary conditions for scalar,
254!-- bottom and top boundary (see also temperature)
255    IF ( passive_scalar )  THEN
256!
257!--    Surface conditions for constant_humidity_flux
258!--    Run loop over all non-natural and natural walls. Note, in wall-datatype
259!--    the k coordinate belongs to the atmospheric grid point, therefore, set
260!--    s_p at k-1
261       IF ( ibc_s_b == 0 ) THEN
262
263          DO  l = 0, 1
264             !$OMP PARALLEL DO PRIVATE( i, j, k )
265             DO  m = 1, bc_h(l)%ns
266                i = bc_h(l)%i(m)
267                j = bc_h(l)%j(m)
268                k = bc_h(l)%k(m)
269                s_p(k+bc_h(l)%koff,j,i) = s(k+bc_h(l)%koff,j,i)
270             ENDDO
271          ENDDO
272
273       ELSE
274
275          DO  l = 0, 1
276             !$OMP PARALLEL DO PRIVATE( i, j, k )
277             DO  m = 1, bc_h(l)%ns
278                i = bc_h(l)%i(m)
279                j = bc_h(l)%j(m)
280                k = bc_h(l)%k(m)
281                s_p(k+bc_h(l)%koff,j,i) = s_p(k,j,i)
282             ENDDO
283          ENDDO
284       ENDIF
285!
286!--    Top boundary condition
287       IF ( ibc_s_t == 0 )  THEN
288          s_p(nzt+1,:,:) = s(nzt+1,:,:)
289       ELSEIF ( ibc_s_t == 1 )  THEN
290          s_p(nzt+1,:,:) = s_p(nzt,:,:)
291       ELSEIF ( ibc_s_t == 2 )  THEN
292          s_p(nzt+1,:,:) = s_p(nzt,:,:) + bc_s_t_val * dzu(nzt+1)
293       ENDIF
294
295    ENDIF 
296!
297!-- Set boundary conditions for subgrid TKE and dissipation (RANS mode)
298    CALL tcm_boundary_conds
299!
300!-- In case of inflow or nest boundary at the south boundary the boundary for v
301!-- is at nys and in case of inflow or nest boundary at the left boundary the
302!-- boundary for u is at nxl. Since in prognostic_equations (cache optimized
303!-- version) these levels are handled as a prognostic level, boundary values
304!-- have to be restored here.
305    IF ( bc_dirichlet_s )  THEN
306       v_p(:,nys,:) = v_p(:,nys-1,:)
307    ELSEIF ( bc_dirichlet_l ) THEN
308       u_p(:,:,nxl) = u_p(:,:,nxl-1)
309    ENDIF
310
311!
312!-- The same restoration for u at i=nxl and v at j=nys as above must be made
313!-- in case of nest boundaries. This must not be done in case of vertical nesting
314!-- mode as in that case the lateral boundaries are actually cyclic.
315!-- Lateral oundary conditions for TKE and dissipation are set
316!-- in tcm_boundary_conds.
317    IF ( nesting_mode /= 'vertical'  .OR.  nesting_offline )  THEN
318       IF ( bc_dirichlet_s )  THEN
319          v_p(:,nys,:) = v_p(:,nys-1,:)
320       ENDIF
321       IF ( bc_dirichlet_l )  THEN
322          u_p(:,:,nxl) = u_p(:,:,nxl-1)
323       ENDIF
324    ENDIF
325
326!
327!-- Lateral boundary conditions for scalar quantities at the outflow.
328!-- Lateral oundary conditions for TKE and dissipation are set
329!-- in tcm_boundary_conds.
330    IF ( bc_radiation_s )  THEN
331       pt_p(:,nys-1,:)     = pt_p(:,nys,:)
332       IF ( humidity )  THEN
333          q_p(:,nys-1,:) = q_p(:,nys,:)
334       ENDIF
335       IF ( passive_scalar )  s_p(:,nys-1,:) = s_p(:,nys,:)
336    ELSEIF ( bc_radiation_n )  THEN
337       pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
338       IF ( humidity )  THEN
339          q_p(:,nyn+1,:) = q_p(:,nyn,:)
340       ENDIF
341       IF ( passive_scalar )  s_p(:,nyn+1,:) = s_p(:,nyn,:)
342    ELSEIF ( bc_radiation_l )  THEN
343       pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
344       IF ( humidity )  THEN
345          q_p(:,:,nxl-1) = q_p(:,:,nxl)
346       ENDIF
347       IF ( passive_scalar )  s_p(:,:,nxl-1) = s_p(:,:,nxl)
348    ELSEIF ( bc_radiation_r )  THEN
349       pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
350       IF ( humidity )  THEN
351          q_p(:,:,nxr+1) = q_p(:,:,nxr)
352       ENDIF
353       IF ( passive_scalar )  s_p(:,:,nxr+1) = s_p(:,:,nxr)
354    ENDIF
355
356!
357!-- Radiation boundary conditions for the velocities at the respective outflow.
358!-- The phase velocity is either assumed to the maximum phase velocity that
359!-- ensures numerical stability (CFL-condition) or calculated after
360!-- Orlanski(1976) and averaged along the outflow boundary.
361    IF ( bc_radiation_s )  THEN
362
363       IF ( use_cmax )  THEN
364          u_p(:,-1,:) = u(:,0,:)
365          v_p(:,0,:)  = v(:,1,:)
366          w_p(:,-1,:) = w(:,0,:)
367       ELSEIF ( .NOT. use_cmax )  THEN
368
369          c_max = dy / dt_3d
370
371          c_u_m_l = 0.0_wp 
372          c_v_m_l = 0.0_wp
373          c_w_m_l = 0.0_wp
374
375          c_u_m = 0.0_wp 
376          c_v_m = 0.0_wp
377          c_w_m = 0.0_wp
378
379!
380!--       Calculate the phase speeds for u, v, and w, first local and then
381!--       average along the outflow boundary.
382          DO  k = nzb+1, nzt+1
383             DO  i = nxl, nxr
384
385                denom = u_m_s(k,0,i) - u_m_s(k,1,i)
386
387                IF ( denom /= 0.0_wp )  THEN
388                   c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / ( denom * tsc(2) )
389                   IF ( c_u(k,i) < 0.0_wp )  THEN
390                      c_u(k,i) = 0.0_wp
391                   ELSEIF ( c_u(k,i) > c_max )  THEN
392                      c_u(k,i) = c_max
393                   ENDIF
394                ELSE
395                   c_u(k,i) = c_max
396                ENDIF
397
398                denom = v_m_s(k,1,i) - v_m_s(k,2,i)
399
400                IF ( denom /= 0.0_wp )  THEN
401                   c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / ( denom * tsc(2) )
402                   IF ( c_v(k,i) < 0.0_wp )  THEN
403                      c_v(k,i) = 0.0_wp
404                   ELSEIF ( c_v(k,i) > c_max )  THEN
405                      c_v(k,i) = c_max
406                   ENDIF
407                ELSE
408                   c_v(k,i) = c_max
409                ENDIF
410
411                denom = w_m_s(k,0,i) - w_m_s(k,1,i)
412
413                IF ( denom /= 0.0_wp )  THEN
414                   c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) )
415                   IF ( c_w(k,i) < 0.0_wp )  THEN
416                      c_w(k,i) = 0.0_wp
417                   ELSEIF ( c_w(k,i) > c_max )  THEN
418                      c_w(k,i) = c_max
419                   ENDIF
420                ELSE
421                   c_w(k,i) = c_max
422                ENDIF
423
424                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
425                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
426                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
427
428             ENDDO
429          ENDDO
430
431#if defined( __parallel )
432          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
433          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
434                              MPI_SUM, comm1dx, ierr )   
435          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
436          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
437                              MPI_SUM, comm1dx, ierr ) 
438          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
439          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
440                              MPI_SUM, comm1dx, ierr ) 
441#else
442          c_u_m = c_u_m_l
443          c_v_m = c_v_m_l
444          c_w_m = c_w_m_l
445#endif
446
447          c_u_m = c_u_m / (nx+1)
448          c_v_m = c_v_m / (nx+1)
449          c_w_m = c_w_m / (nx+1)
450
451!
452!--       Save old timelevels for the next timestep
453          IF ( intermediate_timestep_count == 1 )  THEN
454             u_m_s(:,:,:) = u(:,0:1,:)
455             v_m_s(:,:,:) = v(:,1:2,:)
456             w_m_s(:,:,:) = w(:,0:1,:)
457          ENDIF
458
459!
460!--       Calculate the new velocities
461          DO  k = nzb+1, nzt+1
462             DO  i = nxlg, nxrg
463                u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) *          &
464                                       ( u(k,-1,i) - u(k,0,i) ) * ddy
465
466                v_p(k,0,i)  = v(k,0,i)  - dt_3d * tsc(2) * c_v_m(k) *          &
467                                       ( v(k,0,i) - v(k,1,i) ) * ddy
468
469                w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) *          &
470                                       ( w(k,-1,i) - w(k,0,i) ) * ddy
471             ENDDO
472          ENDDO
473
474!
475!--       Bottom boundary at the outflow
476          IF ( ibc_uv_b == 0 )  THEN
477             u_p(nzb,-1,:) = 0.0_wp 
478             v_p(nzb,0,:)  = 0.0_wp 
479          ELSE
480             u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
481             v_p(nzb,0,:)  =  v_p(nzb+1,0,:)
482          ENDIF
483          w_p(nzb,-1,:) = 0.0_wp
484
485!
486!--       Top boundary at the outflow
487          IF ( ibc_uv_t == 0 )  THEN
488             u_p(nzt+1,-1,:) = u_init(nzt+1)
489             v_p(nzt+1,0,:)  = v_init(nzt+1)
490          ELSE
491             u_p(nzt+1,-1,:) = u_p(nzt,-1,:)
492             v_p(nzt+1,0,:)  = v_p(nzt,0,:)
493          ENDIF
494          w_p(nzt:nzt+1,-1,:) = 0.0_wp
495
496       ENDIF
497
498    ENDIF
499
500    IF ( bc_radiation_n )  THEN
501
502       IF ( use_cmax )  THEN
503          u_p(:,ny+1,:) = u(:,ny,:)
504          v_p(:,ny+1,:) = v(:,ny,:)
505          w_p(:,ny+1,:) = w(:,ny,:)
506       ELSEIF ( .NOT. use_cmax )  THEN
507
508          c_max = dy / dt_3d
509
510          c_u_m_l = 0.0_wp 
511          c_v_m_l = 0.0_wp
512          c_w_m_l = 0.0_wp
513
514          c_u_m = 0.0_wp 
515          c_v_m = 0.0_wp
516          c_w_m = 0.0_wp
517
518!
519!--       Calculate the phase speeds for u, v, and w, first local and then
520!--       average along the outflow boundary.
521          DO  k = nzb+1, nzt+1
522             DO  i = nxl, nxr
523
524                denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
525
526                IF ( denom /= 0.0_wp )  THEN
527                   c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) )
528                   IF ( c_u(k,i) < 0.0_wp )  THEN
529                      c_u(k,i) = 0.0_wp
530                   ELSEIF ( c_u(k,i) > c_max )  THEN
531                      c_u(k,i) = c_max
532                   ENDIF
533                ELSE
534                   c_u(k,i) = c_max
535                ENDIF
536
537                denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
538
539                IF ( denom /= 0.0_wp )  THEN
540                   c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) )
541                   IF ( c_v(k,i) < 0.0_wp )  THEN
542                      c_v(k,i) = 0.0_wp
543                   ELSEIF ( c_v(k,i) > c_max )  THEN
544                      c_v(k,i) = c_max
545                   ENDIF
546                ELSE
547                   c_v(k,i) = c_max
548                ENDIF
549
550                denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
551
552                IF ( denom /= 0.0_wp )  THEN
553                   c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) )
554                   IF ( c_w(k,i) < 0.0_wp )  THEN
555                      c_w(k,i) = 0.0_wp
556                   ELSEIF ( c_w(k,i) > c_max )  THEN
557                      c_w(k,i) = c_max
558                   ENDIF
559                ELSE
560                   c_w(k,i) = c_max
561                ENDIF
562
563                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
564                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
565                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
566
567             ENDDO
568          ENDDO
569
570#if defined( __parallel )   
571          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
572          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
573                              MPI_SUM, comm1dx, ierr )   
574          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
575          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
576                              MPI_SUM, comm1dx, ierr ) 
577          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
578          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
579                              MPI_SUM, comm1dx, ierr ) 
580#else
581          c_u_m = c_u_m_l
582          c_v_m = c_v_m_l
583          c_w_m = c_w_m_l
584#endif
585
586          c_u_m = c_u_m / (nx+1)
587          c_v_m = c_v_m / (nx+1)
588          c_w_m = c_w_m / (nx+1)
589
590!
591!--       Save old timelevels for the next timestep
592          IF ( intermediate_timestep_count == 1 )  THEN
593                u_m_n(:,:,:) = u(:,ny-1:ny,:)
594                v_m_n(:,:,:) = v(:,ny-1:ny,:)
595                w_m_n(:,:,:) = w(:,ny-1:ny,:)
596          ENDIF
597
598!
599!--       Calculate the new velocities
600          DO  k = nzb+1, nzt+1
601             DO  i = nxlg, nxrg
602                u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) *      &
603                                       ( u(k,ny+1,i) - u(k,ny,i) ) * ddy
604
605                v_p(k,ny+1,i) = v(k,ny+1,i)  - dt_3d * tsc(2) * c_v_m(k) *     &
606                                       ( v(k,ny+1,i) - v(k,ny,i) ) * ddy
607
608                w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) *      &
609                                       ( w(k,ny+1,i) - w(k,ny,i) ) * ddy
610             ENDDO
611          ENDDO
612
613!
614!--       Bottom boundary at the outflow
615          IF ( ibc_uv_b == 0 )  THEN
616             u_p(nzb,ny+1,:) = 0.0_wp
617             v_p(nzb,ny+1,:) = 0.0_wp   
618          ELSE
619             u_p(nzb,ny+1,:) =  u_p(nzb+1,ny+1,:)
620             v_p(nzb,ny+1,:) =  v_p(nzb+1,ny+1,:)
621          ENDIF
622          w_p(nzb,ny+1,:) = 0.0_wp
623
624!
625!--       Top boundary at the outflow
626          IF ( ibc_uv_t == 0 )  THEN
627             u_p(nzt+1,ny+1,:) = u_init(nzt+1)
628             v_p(nzt+1,ny+1,:) = v_init(nzt+1)
629          ELSE
630             u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
631             v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
632          ENDIF
633          w_p(nzt:nzt+1,ny+1,:) = 0.0_wp
634
635       ENDIF
636
637    ENDIF
638
639    IF ( bc_radiation_l )  THEN
640
641       IF ( use_cmax )  THEN
642          u_p(:,:,0)  = u(:,:,1)
643          v_p(:,:,-1) = v(:,:,0)
644          w_p(:,:,-1) = w(:,:,0)
645       ELSEIF ( .NOT. use_cmax )  THEN
646
647          c_max = dx / dt_3d
648
649          c_u_m_l = 0.0_wp 
650          c_v_m_l = 0.0_wp
651          c_w_m_l = 0.0_wp
652
653          c_u_m = 0.0_wp 
654          c_v_m = 0.0_wp
655          c_w_m = 0.0_wp
656
657!
658!--       Calculate the phase speeds for u, v, and w, first local and then
659!--       average along the outflow boundary.
660          DO  k = nzb+1, nzt+1
661             DO  j = nys, nyn
662
663                denom = u_m_l(k,j,1) - u_m_l(k,j,2)
664
665                IF ( denom /= 0.0_wp )  THEN
666                   c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) )
667                   IF ( c_u(k,j) < 0.0_wp )  THEN
668                      c_u(k,j) = 0.0_wp
669                   ELSEIF ( c_u(k,j) > c_max )  THEN
670                      c_u(k,j) = c_max
671                   ENDIF
672                ELSE
673                   c_u(k,j) = c_max
674                ENDIF
675
676                denom = v_m_l(k,j,0) - v_m_l(k,j,1)
677
678                IF ( denom /= 0.0_wp )  THEN
679                   c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) )
680                   IF ( c_v(k,j) < 0.0_wp )  THEN
681                      c_v(k,j) = 0.0_wp
682                   ELSEIF ( c_v(k,j) > c_max )  THEN
683                      c_v(k,j) = c_max
684                   ENDIF
685                ELSE
686                   c_v(k,j) = c_max
687                ENDIF
688
689                denom = w_m_l(k,j,0) - w_m_l(k,j,1)
690
691                IF ( denom /= 0.0_wp )  THEN
692                   c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) )
693                   IF ( c_w(k,j) < 0.0_wp )  THEN
694                      c_w(k,j) = 0.0_wp
695                   ELSEIF ( c_w(k,j) > c_max )  THEN
696                      c_w(k,j) = c_max
697                   ENDIF
698                ELSE
699                   c_w(k,j) = c_max
700                ENDIF
701
702                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
703                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
704                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
705
706             ENDDO
707          ENDDO
708
709#if defined( __parallel )   
710          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
711          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
712                              MPI_SUM, comm1dy, ierr )   
713          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
714          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
715                              MPI_SUM, comm1dy, ierr ) 
716          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
717          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
718                              MPI_SUM, comm1dy, ierr ) 
719#else
720          c_u_m = c_u_m_l
721          c_v_m = c_v_m_l
722          c_w_m = c_w_m_l
723#endif
724
725          c_u_m = c_u_m / (ny+1)
726          c_v_m = c_v_m / (ny+1)
727          c_w_m = c_w_m / (ny+1)
728
729!
730!--       Save old timelevels for the next timestep
731          IF ( intermediate_timestep_count == 1 )  THEN
732                u_m_l(:,:,:) = u(:,:,1:2)
733                v_m_l(:,:,:) = v(:,:,0:1)
734                w_m_l(:,:,:) = w(:,:,0:1)
735          ENDIF
736
737!
738!--       Calculate the new velocities
739          DO  k = nzb+1, nzt+1
740             DO  j = nysg, nyng
741                u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) *            &
742                                       ( u(k,j,0) - u(k,j,1) ) * ddx
743
744                v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) *          &
745                                       ( v(k,j,-1) - v(k,j,0) ) * ddx
746
747                w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) *          &
748                                       ( w(k,j,-1) - w(k,j,0) ) * ddx
749             ENDDO
750          ENDDO
751
752!
753!--       Bottom boundary at the outflow
754          IF ( ibc_uv_b == 0 )  THEN
755             u_p(nzb,:,0)  = 0.0_wp 
756             v_p(nzb,:,-1) = 0.0_wp
757          ELSE
758             u_p(nzb,:,0)  =  u_p(nzb+1,:,0)
759             v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
760          ENDIF
761          w_p(nzb,:,-1) = 0.0_wp
762
763!
764!--       Top boundary at the outflow
765          IF ( ibc_uv_t == 0 )  THEN
766             u_p(nzt+1,:,0)  = u_init(nzt+1)
767             v_p(nzt+1,:,-1) = v_init(nzt+1)
768          ELSE
769             u_p(nzt+1,:,0)  = u_p(nzt,:,0)
770             v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
771          ENDIF
772          w_p(nzt:nzt+1,:,-1) = 0.0_wp
773
774       ENDIF
775
776    ENDIF
777
778    IF ( bc_radiation_r )  THEN
779
780       IF ( use_cmax )  THEN
781          u_p(:,:,nx+1) = u(:,:,nx)
782          v_p(:,:,nx+1) = v(:,:,nx)
783          w_p(:,:,nx+1) = w(:,:,nx)
784       ELSEIF ( .NOT. use_cmax )  THEN
785
786          c_max = dx / dt_3d
787
788          c_u_m_l = 0.0_wp 
789          c_v_m_l = 0.0_wp
790          c_w_m_l = 0.0_wp
791
792          c_u_m = 0.0_wp 
793          c_v_m = 0.0_wp
794          c_w_m = 0.0_wp
795
796!
797!--       Calculate the phase speeds for u, v, and w, first local and then
798!--       average along the outflow boundary.
799          DO  k = nzb+1, nzt+1
800             DO  j = nys, nyn
801
802                denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
803
804                IF ( denom /= 0.0_wp )  THEN
805                   c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) )
806                   IF ( c_u(k,j) < 0.0_wp )  THEN
807                      c_u(k,j) = 0.0_wp
808                   ELSEIF ( c_u(k,j) > c_max )  THEN
809                      c_u(k,j) = c_max
810                   ENDIF
811                ELSE
812                   c_u(k,j) = c_max
813                ENDIF
814
815                denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
816
817                IF ( denom /= 0.0_wp )  THEN
818                   c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) )
819                   IF ( c_v(k,j) < 0.0_wp )  THEN
820                      c_v(k,j) = 0.0_wp
821                   ELSEIF ( c_v(k,j) > c_max )  THEN
822                      c_v(k,j) = c_max
823                   ENDIF
824                ELSE
825                   c_v(k,j) = c_max
826                ENDIF
827
828                denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
829
830                IF ( denom /= 0.0_wp )  THEN
831                   c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) )
832                   IF ( c_w(k,j) < 0.0_wp )  THEN
833                      c_w(k,j) = 0.0_wp
834                   ELSEIF ( c_w(k,j) > c_max )  THEN
835                      c_w(k,j) = c_max
836                   ENDIF
837                ELSE
838                   c_w(k,j) = c_max
839                ENDIF
840
841                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
842                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
843                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
844
845             ENDDO
846          ENDDO
847
848#if defined( __parallel )
849          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
850          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
851                              MPI_SUM, comm1dy, ierr )
852          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
853          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
854                              MPI_SUM, comm1dy, ierr )
855          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
856          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
857                              MPI_SUM, comm1dy, ierr )
858#else
859          c_u_m = c_u_m_l
860          c_v_m = c_v_m_l
861          c_w_m = c_w_m_l
862#endif
863
864          c_u_m = c_u_m / (ny+1)
865          c_v_m = c_v_m / (ny+1)
866          c_w_m = c_w_m / (ny+1)
867
868!
869!--       Save old timelevels for the next timestep
870          IF ( intermediate_timestep_count == 1 )  THEN
871                u_m_r(:,:,:) = u(:,:,nx-1:nx)
872                v_m_r(:,:,:) = v(:,:,nx-1:nx)
873                w_m_r(:,:,:) = w(:,:,nx-1:nx)
874          ENDIF
875
876!
877!--       Calculate the new velocities
878          DO  k = nzb+1, nzt+1
879             DO  j = nysg, nyng
880                u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) *      &
881                                       ( u(k,j,nx+1) - u(k,j,nx) ) * ddx
882
883                v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) *      &
884                                       ( v(k,j,nx+1) - v(k,j,nx) ) * ddx
885
886                w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) *      &
887                                       ( w(k,j,nx+1) - w(k,j,nx) ) * ddx
888             ENDDO
889          ENDDO
890
891!
892!--       Bottom boundary at the outflow
893          IF ( ibc_uv_b == 0 )  THEN
894             u_p(nzb,:,nx+1) = 0.0_wp
895             v_p(nzb,:,nx+1) = 0.0_wp 
896          ELSE
897             u_p(nzb,:,nx+1) =  u_p(nzb+1,:,nx+1)
898             v_p(nzb,:,nx+1) =  v_p(nzb+1,:,nx+1)
899          ENDIF
900          w_p(nzb,:,nx+1) = 0.0_wp
901
902!
903!--       Top boundary at the outflow
904          IF ( ibc_uv_t == 0 )  THEN
905             u_p(nzt+1,:,nx+1) = u_init(nzt+1)
906             v_p(nzt+1,:,nx+1) = v_init(nzt+1)
907          ELSE
908             u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
909             v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
910          ENDIF
911          w_p(nzt:nzt+1,:,nx+1) = 0.0_wp
912
913       ENDIF
914
915    ENDIF
916
917 END SUBROUTINE boundary_conds
Note: See TracBrowser for help on using the repository browser.