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

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

Introducing module interface for boundary conditions and move module specific boundary conditions into their modules

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