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

Last change on this file since 4218 was 4182, checked in by scharf, 5 years ago
  • corrected "Former revisions" section
  • minor formatting in "Former revisions" section
  • added "Author" section
  • Property svn:keywords set to Id
File size: 35.0 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 4182 2019-08-22 15:20:23Z knoop $
27! Corrected "Former revisions" section
28!
29! 4102 2019-07-17 16:00:03Z suehring
30! - Revise setting for boundary conditions at horizontal walls, use the offset
31!   index that belongs to the data structure instead of pre-calculating
32!   the offset index for each facing.
33! - Set boundary conditions for bulk-cloud quantities also at downward facing
34!   surfaces
35!
36! 4087 2019-07-11 11:35:23Z gronemeier
37! Add comment
38!
39! 4086 2019-07-11 05:55:44Z gronemeier
40! Bugfix: use constant-flux layer condition for e in all rans modes
41!
42! 3879 2019-04-08 20:25:23Z knoop
43! Bugfix, do not set boundary conditions for potential temperature in neutral
44! runs.
45!
46! 3655 2019-01-07 16:51:22Z knoop
47! OpenACC port for SPEC
48!
49! Revision 1.1  1997/09/12 06:21:34  raasch
50! Initial revision
51!
52!
53! Description:
54! ------------
55!> Boundary conditions for the prognostic quantities.
56!> One additional bottom boundary condition is applied for the TKE (=(u*)**2)
57!> in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
58!> handled in routine exchange_horiz. Pressure boundary conditions are
59!> explicitly set in routines pres, poisfft, poismg and sor.
60!------------------------------------------------------------------------------!
61 SUBROUTINE boundary_conds
62 
63
64    USE arrays_3d,                                                             &
65        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,  &
66               dzu, nc_p, nr_p, pt, pt_init, pt_p, q,                          &
67               q_p, qc_p, qr_p, s, s_p, sa, sa_p, u, u_init, u_m_l, u_m_n,     &
68               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,  &
69               w, w_p, w_m_l, w_m_n, w_m_r, w_m_s
70
71    USE bulk_cloud_model_mod,                                                  &
72        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
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
282       IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
283!             
284!--       Surface conditions cloud water (Dirichlet)
285!--       Run loop over all non-natural and natural walls. Note, in wall-datatype
286!--       the k coordinate belongs to the atmospheric grid point, therefore, set
287!--       qr_p and nr_p at upward (k-1) and downward-facing (k+1) walls
288          DO  l = 0, 1
289          !$OMP PARALLEL DO PRIVATE( i, j, k )
290             DO  m = 1, bc_h(l)%ns
291                i = bc_h(l)%i(m)           
292                j = bc_h(l)%j(m)
293                k = bc_h(l)%k(m)
294                qc_p(k+bc_h(l)%koff,j,i) = 0.0_wp
295                nc_p(k+bc_h(l)%koff,j,i) = 0.0_wp
296             ENDDO
297          ENDDO
298!
299!--       Top boundary condition for cloud water (Dirichlet)
300          qc_p(nzt+1,:,:) = 0.0_wp
301          nc_p(nzt+1,:,:) = 0.0_wp
302           
303       ENDIF
304
305       IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
306!             
307!--       Surface conditions rain water (Dirichlet)
308!--       Run loop over all non-natural and natural walls. Note, in wall-datatype
309!--       the k coordinate belongs to the atmospheric grid point, therefore, set
310!--       qr_p and nr_p at upward (k-1) and downward-facing (k+1) walls
311          DO  l = 0, 1
312          !$OMP PARALLEL DO PRIVATE( i, j, k )
313             DO  m = 1, bc_h(l)%ns
314                i = bc_h(l)%i(m)           
315                j = bc_h(l)%j(m)
316                k = bc_h(l)%k(m)
317                qr_p(k+bc_h(l)%koff,j,i) = 0.0_wp
318                nr_p(k+bc_h(l)%koff,j,i) = 0.0_wp
319             ENDDO
320          ENDDO
321!
322!--       Top boundary condition for rain water (Dirichlet)
323          qr_p(nzt+1,:,:) = 0.0_wp
324          nr_p(nzt+1,:,:) = 0.0_wp
325           
326       ENDIF
327    ENDIF
328!
329!-- Boundary conditions for scalar,
330!-- bottom and top boundary (see also temperature)
331    IF ( passive_scalar )  THEN
332!
333!--    Surface conditions for constant_humidity_flux
334!--    Run loop over all non-natural and natural walls. Note, in wall-datatype
335!--    the k coordinate belongs to the atmospheric grid point, therefore, set
336!--    s_p at k-1
337       IF ( ibc_s_b == 0 ) THEN
338         
339          DO  l = 0, 1
340             !$OMP PARALLEL DO PRIVATE( i, j, k )
341             DO  m = 1, bc_h(l)%ns
342                i = bc_h(l)%i(m)           
343                j = bc_h(l)%j(m)
344                k = bc_h(l)%k(m)
345                s_p(k+bc_h(l)%koff,j,i) = s(k+bc_h(l)%koff,j,i)
346             ENDDO
347          ENDDO
348         
349       ELSE
350         
351          DO  l = 0, 1
352             !$OMP PARALLEL DO PRIVATE( i, j, k )
353             DO  m = 1, bc_h(l)%ns
354                i = bc_h(l)%i(m)           
355                j = bc_h(l)%j(m)
356                k = bc_h(l)%k(m)
357                s_p(k+bc_h(l)%koff,j,i) = s_p(k,j,i)
358             ENDDO
359          ENDDO
360       ENDIF
361!
362!--    Top boundary condition
363       IF ( ibc_s_t == 0 )  THEN
364          s_p(nzt+1,:,:) = s(nzt+1,:,:)
365       ELSEIF ( ibc_s_t == 1 )  THEN
366          s_p(nzt+1,:,:) = s_p(nzt,:,:)
367       ELSEIF ( ibc_s_t == 2 )  THEN
368          s_p(nzt+1,:,:) = s_p(nzt,:,:) + bc_s_t_val * dzu(nzt+1)
369       ENDIF
370
371    ENDIF 
372!
373!-- Set boundary conditions for subgrid TKE and dissipation (RANS mode)
374    CALL tcm_boundary_conds
375!
376!-- Top/bottom boundary conditions for chemical species
377    IF ( air_chemistry )  CALL chem_boundary_conds( 'set_bc_bottomtop' )
378!
379!-- In case of inflow or nest boundary at the south boundary the boundary for v
380!-- is at nys and in case of inflow or nest boundary at the left boundary the
381!-- boundary for u is at nxl. Since in prognostic_equations (cache optimized
382!-- version) these levels are handled as a prognostic level, boundary values
383!-- have to be restored here.
384    IF ( bc_dirichlet_s )  THEN
385       v_p(:,nys,:) = v_p(:,nys-1,:)
386    ELSEIF ( bc_dirichlet_l ) THEN
387       u_p(:,:,nxl) = u_p(:,:,nxl-1)
388    ENDIF
389
390!
391!-- The same restoration for u at i=nxl and v at j=nys as above must be made
392!-- in case of nest boundaries. This must not be done in case of vertical nesting
393!-- mode as in that case the lateral boundaries are actually cyclic.
394!-- Lateral oundary conditions for TKE and dissipation are set
395!-- in tcm_boundary_conds.
396    IF ( nesting_mode /= 'vertical'  .OR.  nesting_offline )  THEN
397       IF ( bc_dirichlet_s )  THEN
398          v_p(:,nys,:) = v_p(:,nys-1,:)
399       ENDIF
400       IF ( bc_dirichlet_l )  THEN
401          u_p(:,:,nxl) = u_p(:,:,nxl-1)
402       ENDIF
403    ENDIF
404
405!
406!-- Lateral boundary conditions for scalar quantities at the outflow.
407!-- Lateral oundary conditions for TKE and dissipation are set
408!-- in tcm_boundary_conds.
409    IF ( bc_radiation_s )  THEN
410       pt_p(:,nys-1,:)     = pt_p(:,nys,:)
411       IF ( humidity )  THEN
412          q_p(:,nys-1,:) = q_p(:,nys,:)
413          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
414             qc_p(:,nys-1,:) = qc_p(:,nys,:)
415             nc_p(:,nys-1,:) = nc_p(:,nys,:)
416          ENDIF
417          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
418             qr_p(:,nys-1,:) = qr_p(:,nys,:)
419             nr_p(:,nys-1,:) = nr_p(:,nys,:)
420          ENDIF
421       ENDIF
422       IF ( passive_scalar )  s_p(:,nys-1,:) = s_p(:,nys,:)
423    ELSEIF ( bc_radiation_n )  THEN
424       pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
425       IF ( humidity )  THEN
426          q_p(:,nyn+1,:) = q_p(:,nyn,:)
427          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
428             qc_p(:,nyn+1,:) = qc_p(:,nyn,:)
429             nc_p(:,nyn+1,:) = nc_p(:,nyn,:)
430          ENDIF
431          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
432             qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
433             nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
434          ENDIF
435       ENDIF
436       IF ( passive_scalar )  s_p(:,nyn+1,:) = s_p(:,nyn,:)
437    ELSEIF ( bc_radiation_l )  THEN
438       pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
439       IF ( humidity )  THEN
440          q_p(:,:,nxl-1) = q_p(:,:,nxl)
441          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
442             qc_p(:,:,nxl-1) = qc_p(:,:,nxl)
443             nc_p(:,:,nxl-1) = nc_p(:,:,nxl)
444          ENDIF
445          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
446             qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
447             nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
448          ENDIF
449       ENDIF
450       IF ( passive_scalar )  s_p(:,:,nxl-1) = s_p(:,:,nxl)
451    ELSEIF ( bc_radiation_r )  THEN
452       pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
453       IF ( humidity )  THEN
454          q_p(:,:,nxr+1) = q_p(:,:,nxr)
455          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
456             qc_p(:,:,nxr+1) = qc_p(:,:,nxr)
457             nc_p(:,:,nxr+1) = nc_p(:,:,nxr)
458          ENDIF
459          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
460             qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
461             nr_p(:,:,nxr+1) = nr_p(:,:,nxr)
462          ENDIF
463       ENDIF
464       IF ( passive_scalar )  s_p(:,:,nxr+1) = s_p(:,:,nxr)
465    ENDIF
466
467!
468!-- Lateral boundary conditions for chemical species
469    IF ( air_chemistry )  CALL chem_boundary_conds( 'set_bc_lateral' )   
470
471!
472!-- Radiation boundary conditions for the velocities at the respective outflow.
473!-- The phase velocity is either assumed to the maximum phase velocity that
474!-- ensures numerical stability (CFL-condition) or calculated after
475!-- Orlanski(1976) and averaged along the outflow boundary.
476    IF ( bc_radiation_s )  THEN
477
478       IF ( use_cmax )  THEN
479          u_p(:,-1,:) = u(:,0,:)
480          v_p(:,0,:)  = v(:,1,:)
481          w_p(:,-1,:) = w(:,0,:)         
482       ELSEIF ( .NOT. use_cmax )  THEN
483
484          c_max = dy / dt_3d
485
486          c_u_m_l = 0.0_wp 
487          c_v_m_l = 0.0_wp
488          c_w_m_l = 0.0_wp
489
490          c_u_m = 0.0_wp 
491          c_v_m = 0.0_wp
492          c_w_m = 0.0_wp
493
494!
495!--       Calculate the phase speeds for u, v, and w, first local and then
496!--       average along the outflow boundary.
497          DO  k = nzb+1, nzt+1
498             DO  i = nxl, nxr
499
500                denom = u_m_s(k,0,i) - u_m_s(k,1,i)
501
502                IF ( denom /= 0.0_wp )  THEN
503                   c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / ( denom * tsc(2) )
504                   IF ( c_u(k,i) < 0.0_wp )  THEN
505                      c_u(k,i) = 0.0_wp
506                   ELSEIF ( c_u(k,i) > c_max )  THEN
507                      c_u(k,i) = c_max
508                   ENDIF
509                ELSE
510                   c_u(k,i) = c_max
511                ENDIF
512
513                denom = v_m_s(k,1,i) - v_m_s(k,2,i)
514
515                IF ( denom /= 0.0_wp )  THEN
516                   c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / ( denom * tsc(2) )
517                   IF ( c_v(k,i) < 0.0_wp )  THEN
518                      c_v(k,i) = 0.0_wp
519                   ELSEIF ( c_v(k,i) > c_max )  THEN
520                      c_v(k,i) = c_max
521                   ENDIF
522                ELSE
523                   c_v(k,i) = c_max
524                ENDIF
525
526                denom = w_m_s(k,0,i) - w_m_s(k,1,i)
527
528                IF ( denom /= 0.0_wp )  THEN
529                   c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) )
530                   IF ( c_w(k,i) < 0.0_wp )  THEN
531                      c_w(k,i) = 0.0_wp
532                   ELSEIF ( c_w(k,i) > c_max )  THEN
533                      c_w(k,i) = c_max
534                   ENDIF
535                ELSE
536                   c_w(k,i) = c_max
537                ENDIF
538
539                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
540                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
541                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
542
543             ENDDO
544          ENDDO
545
546#if defined( __parallel )   
547          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
548          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
549                              MPI_SUM, comm1dx, ierr )   
550          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
551          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
552                              MPI_SUM, comm1dx, ierr ) 
553          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
554          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
555                              MPI_SUM, comm1dx, ierr ) 
556#else
557          c_u_m = c_u_m_l
558          c_v_m = c_v_m_l
559          c_w_m = c_w_m_l
560#endif
561
562          c_u_m = c_u_m / (nx+1)
563          c_v_m = c_v_m / (nx+1)
564          c_w_m = c_w_m / (nx+1)
565
566!
567!--       Save old timelevels for the next timestep
568          IF ( intermediate_timestep_count == 1 )  THEN
569             u_m_s(:,:,:) = u(:,0:1,:)
570             v_m_s(:,:,:) = v(:,1:2,:)
571             w_m_s(:,:,:) = w(:,0:1,:)
572          ENDIF
573
574!
575!--       Calculate the new velocities
576          DO  k = nzb+1, nzt+1
577             DO  i = nxlg, nxrg
578                u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) *          &
579                                       ( u(k,-1,i) - u(k,0,i) ) * ddy
580
581                v_p(k,0,i)  = v(k,0,i)  - dt_3d * tsc(2) * c_v_m(k) *          &
582                                       ( v(k,0,i) - v(k,1,i) ) * ddy
583
584                w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) *          &
585                                       ( w(k,-1,i) - w(k,0,i) ) * ddy
586             ENDDO
587          ENDDO
588
589!
590!--       Bottom boundary at the outflow
591          IF ( ibc_uv_b == 0 )  THEN
592             u_p(nzb,-1,:) = 0.0_wp 
593             v_p(nzb,0,:)  = 0.0_wp 
594          ELSE                   
595             u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
596             v_p(nzb,0,:)  =  v_p(nzb+1,0,:)
597          ENDIF
598          w_p(nzb,-1,:) = 0.0_wp
599
600!
601!--       Top boundary at the outflow
602          IF ( ibc_uv_t == 0 )  THEN
603             u_p(nzt+1,-1,:) = u_init(nzt+1)
604             v_p(nzt+1,0,:)  = v_init(nzt+1)
605          ELSE
606             u_p(nzt+1,-1,:) = u_p(nzt,-1,:)
607             v_p(nzt+1,0,:)  = v_p(nzt,0,:)
608          ENDIF
609          w_p(nzt:nzt+1,-1,:) = 0.0_wp
610
611       ENDIF
612
613    ENDIF
614
615    IF ( bc_radiation_n )  THEN
616
617       IF ( use_cmax )  THEN
618          u_p(:,ny+1,:) = u(:,ny,:)
619          v_p(:,ny+1,:) = v(:,ny,:)
620          w_p(:,ny+1,:) = w(:,ny,:)         
621       ELSEIF ( .NOT. use_cmax )  THEN
622
623          c_max = dy / dt_3d
624
625          c_u_m_l = 0.0_wp 
626          c_v_m_l = 0.0_wp
627          c_w_m_l = 0.0_wp
628
629          c_u_m = 0.0_wp 
630          c_v_m = 0.0_wp
631          c_w_m = 0.0_wp
632
633!
634!--       Calculate the phase speeds for u, v, and w, first local and then
635!--       average along the outflow boundary.
636          DO  k = nzb+1, nzt+1
637             DO  i = nxl, nxr
638
639                denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
640
641                IF ( denom /= 0.0_wp )  THEN
642                   c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) )
643                   IF ( c_u(k,i) < 0.0_wp )  THEN
644                      c_u(k,i) = 0.0_wp
645                   ELSEIF ( c_u(k,i) > c_max )  THEN
646                      c_u(k,i) = c_max
647                   ENDIF
648                ELSE
649                   c_u(k,i) = c_max
650                ENDIF
651
652                denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
653
654                IF ( denom /= 0.0_wp )  THEN
655                   c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) )
656                   IF ( c_v(k,i) < 0.0_wp )  THEN
657                      c_v(k,i) = 0.0_wp
658                   ELSEIF ( c_v(k,i) > c_max )  THEN
659                      c_v(k,i) = c_max
660                   ENDIF
661                ELSE
662                   c_v(k,i) = c_max
663                ENDIF
664
665                denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
666
667                IF ( denom /= 0.0_wp )  THEN
668                   c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) )
669                   IF ( c_w(k,i) < 0.0_wp )  THEN
670                      c_w(k,i) = 0.0_wp
671                   ELSEIF ( c_w(k,i) > c_max )  THEN
672                      c_w(k,i) = c_max
673                   ENDIF
674                ELSE
675                   c_w(k,i) = c_max
676                ENDIF
677
678                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
679                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
680                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
681
682             ENDDO
683          ENDDO
684
685#if defined( __parallel )   
686          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
687          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
688                              MPI_SUM, comm1dx, ierr )   
689          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
690          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
691                              MPI_SUM, comm1dx, ierr ) 
692          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
693          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
694                              MPI_SUM, comm1dx, ierr ) 
695#else
696          c_u_m = c_u_m_l
697          c_v_m = c_v_m_l
698          c_w_m = c_w_m_l
699#endif
700
701          c_u_m = c_u_m / (nx+1)
702          c_v_m = c_v_m / (nx+1)
703          c_w_m = c_w_m / (nx+1)
704
705!
706!--       Save old timelevels for the next timestep
707          IF ( intermediate_timestep_count == 1 )  THEN
708                u_m_n(:,:,:) = u(:,ny-1:ny,:)
709                v_m_n(:,:,:) = v(:,ny-1:ny,:)
710                w_m_n(:,:,:) = w(:,ny-1:ny,:)
711          ENDIF
712
713!
714!--       Calculate the new velocities
715          DO  k = nzb+1, nzt+1
716             DO  i = nxlg, nxrg
717                u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) *      &
718                                       ( u(k,ny+1,i) - u(k,ny,i) ) * ddy
719
720                v_p(k,ny+1,i) = v(k,ny+1,i)  - dt_3d * tsc(2) * c_v_m(k) *     &
721                                       ( v(k,ny+1,i) - v(k,ny,i) ) * ddy
722
723                w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) *      &
724                                       ( w(k,ny+1,i) - w(k,ny,i) ) * ddy
725             ENDDO
726          ENDDO
727
728!
729!--       Bottom boundary at the outflow
730          IF ( ibc_uv_b == 0 )  THEN
731             u_p(nzb,ny+1,:) = 0.0_wp
732             v_p(nzb,ny+1,:) = 0.0_wp   
733          ELSE                   
734             u_p(nzb,ny+1,:) =  u_p(nzb+1,ny+1,:)
735             v_p(nzb,ny+1,:) =  v_p(nzb+1,ny+1,:)
736          ENDIF
737          w_p(nzb,ny+1,:) = 0.0_wp
738
739!
740!--       Top boundary at the outflow
741          IF ( ibc_uv_t == 0 )  THEN
742             u_p(nzt+1,ny+1,:) = u_init(nzt+1)
743             v_p(nzt+1,ny+1,:) = v_init(nzt+1)
744          ELSE
745             u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
746             v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
747          ENDIF
748          w_p(nzt:nzt+1,ny+1,:) = 0.0_wp
749
750       ENDIF
751
752    ENDIF
753
754    IF ( bc_radiation_l )  THEN
755
756       IF ( use_cmax )  THEN
757          u_p(:,:,0)  = u(:,:,1)
758          v_p(:,:,-1) = v(:,:,0)
759          w_p(:,:,-1) = w(:,:,0)         
760       ELSEIF ( .NOT. use_cmax )  THEN
761
762          c_max = dx / dt_3d
763
764          c_u_m_l = 0.0_wp 
765          c_v_m_l = 0.0_wp
766          c_w_m_l = 0.0_wp
767
768          c_u_m = 0.0_wp 
769          c_v_m = 0.0_wp
770          c_w_m = 0.0_wp
771
772!
773!--       Calculate the phase speeds for u, v, and w, first local and then
774!--       average along the outflow boundary.
775          DO  k = nzb+1, nzt+1
776             DO  j = nys, nyn
777
778                denom = u_m_l(k,j,1) - u_m_l(k,j,2)
779
780                IF ( denom /= 0.0_wp )  THEN
781                   c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) )
782                   IF ( c_u(k,j) < 0.0_wp )  THEN
783                      c_u(k,j) = 0.0_wp
784                   ELSEIF ( c_u(k,j) > c_max )  THEN
785                      c_u(k,j) = c_max
786                   ENDIF
787                ELSE
788                   c_u(k,j) = c_max
789                ENDIF
790
791                denom = v_m_l(k,j,0) - v_m_l(k,j,1)
792
793                IF ( denom /= 0.0_wp )  THEN
794                   c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) )
795                   IF ( c_v(k,j) < 0.0_wp )  THEN
796                      c_v(k,j) = 0.0_wp
797                   ELSEIF ( c_v(k,j) > c_max )  THEN
798                      c_v(k,j) = c_max
799                   ENDIF
800                ELSE
801                   c_v(k,j) = c_max
802                ENDIF
803
804                denom = w_m_l(k,j,0) - w_m_l(k,j,1)
805
806                IF ( denom /= 0.0_wp )  THEN
807                   c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) )
808                   IF ( c_w(k,j) < 0.0_wp )  THEN
809                      c_w(k,j) = 0.0_wp
810                   ELSEIF ( c_w(k,j) > c_max )  THEN
811                      c_w(k,j) = c_max
812                   ENDIF
813                ELSE
814                   c_w(k,j) = c_max
815                ENDIF
816
817                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
818                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
819                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
820
821             ENDDO
822          ENDDO
823
824#if defined( __parallel )   
825          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
826          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
827                              MPI_SUM, comm1dy, ierr )   
828          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
829          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
830                              MPI_SUM, comm1dy, ierr ) 
831          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
832          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
833                              MPI_SUM, comm1dy, ierr ) 
834#else
835          c_u_m = c_u_m_l
836          c_v_m = c_v_m_l
837          c_w_m = c_w_m_l
838#endif
839
840          c_u_m = c_u_m / (ny+1)
841          c_v_m = c_v_m / (ny+1)
842          c_w_m = c_w_m / (ny+1)
843
844!
845!--       Save old timelevels for the next timestep
846          IF ( intermediate_timestep_count == 1 )  THEN
847                u_m_l(:,:,:) = u(:,:,1:2)
848                v_m_l(:,:,:) = v(:,:,0:1)
849                w_m_l(:,:,:) = w(:,:,0:1)
850          ENDIF
851
852!
853!--       Calculate the new velocities
854          DO  k = nzb+1, nzt+1
855             DO  j = nysg, nyng
856                u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) *            &
857                                       ( u(k,j,0) - u(k,j,1) ) * ddx
858
859                v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) *          &
860                                       ( v(k,j,-1) - v(k,j,0) ) * ddx
861
862                w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) *          &
863                                       ( w(k,j,-1) - w(k,j,0) ) * ddx
864             ENDDO
865          ENDDO
866
867!
868!--       Bottom boundary at the outflow
869          IF ( ibc_uv_b == 0 )  THEN
870             u_p(nzb,:,0)  = 0.0_wp 
871             v_p(nzb,:,-1) = 0.0_wp
872          ELSE                   
873             u_p(nzb,:,0)  =  u_p(nzb+1,:,0)
874             v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
875          ENDIF
876          w_p(nzb,:,-1) = 0.0_wp
877
878!
879!--       Top boundary at the outflow
880          IF ( ibc_uv_t == 0 )  THEN
881             u_p(nzt+1,:,0)  = u_init(nzt+1)
882             v_p(nzt+1,:,-1) = v_init(nzt+1)
883          ELSE
884             u_p(nzt+1,:,0)  = u_p(nzt,:,0)
885             v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
886          ENDIF
887          w_p(nzt:nzt+1,:,-1) = 0.0_wp
888
889       ENDIF
890
891    ENDIF
892
893    IF ( bc_radiation_r )  THEN
894
895       IF ( use_cmax )  THEN
896          u_p(:,:,nx+1) = u(:,:,nx)
897          v_p(:,:,nx+1) = v(:,:,nx)
898          w_p(:,:,nx+1) = w(:,:,nx)         
899       ELSEIF ( .NOT. use_cmax )  THEN
900
901          c_max = dx / dt_3d
902
903          c_u_m_l = 0.0_wp 
904          c_v_m_l = 0.0_wp
905          c_w_m_l = 0.0_wp
906
907          c_u_m = 0.0_wp 
908          c_v_m = 0.0_wp
909          c_w_m = 0.0_wp
910
911!
912!--       Calculate the phase speeds for u, v, and w, first local and then
913!--       average along the outflow boundary.
914          DO  k = nzb+1, nzt+1
915             DO  j = nys, nyn
916
917                denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
918
919                IF ( denom /= 0.0_wp )  THEN
920                   c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) )
921                   IF ( c_u(k,j) < 0.0_wp )  THEN
922                      c_u(k,j) = 0.0_wp
923                   ELSEIF ( c_u(k,j) > c_max )  THEN
924                      c_u(k,j) = c_max
925                   ENDIF
926                ELSE
927                   c_u(k,j) = c_max
928                ENDIF
929
930                denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
931
932                IF ( denom /= 0.0_wp )  THEN
933                   c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) )
934                   IF ( c_v(k,j) < 0.0_wp )  THEN
935                      c_v(k,j) = 0.0_wp
936                   ELSEIF ( c_v(k,j) > c_max )  THEN
937                      c_v(k,j) = c_max
938                   ENDIF
939                ELSE
940                   c_v(k,j) = c_max
941                ENDIF
942
943                denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
944
945                IF ( denom /= 0.0_wp )  THEN
946                   c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) )
947                   IF ( c_w(k,j) < 0.0_wp )  THEN
948                      c_w(k,j) = 0.0_wp
949                   ELSEIF ( c_w(k,j) > c_max )  THEN
950                      c_w(k,j) = c_max
951                   ENDIF
952                ELSE
953                   c_w(k,j) = c_max
954                ENDIF
955
956                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
957                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
958                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
959
960             ENDDO
961          ENDDO
962
963#if defined( __parallel )   
964          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
965          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
966                              MPI_SUM, comm1dy, ierr )   
967          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
968          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
969                              MPI_SUM, comm1dy, ierr ) 
970          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
971          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
972                              MPI_SUM, comm1dy, ierr ) 
973#else
974          c_u_m = c_u_m_l
975          c_v_m = c_v_m_l
976          c_w_m = c_w_m_l
977#endif
978
979          c_u_m = c_u_m / (ny+1)
980          c_v_m = c_v_m / (ny+1)
981          c_w_m = c_w_m / (ny+1)
982
983!
984!--       Save old timelevels for the next timestep
985          IF ( intermediate_timestep_count == 1 )  THEN
986                u_m_r(:,:,:) = u(:,:,nx-1:nx)
987                v_m_r(:,:,:) = v(:,:,nx-1:nx)
988                w_m_r(:,:,:) = w(:,:,nx-1:nx)
989          ENDIF
990
991!
992!--       Calculate the new velocities
993          DO  k = nzb+1, nzt+1
994             DO  j = nysg, nyng
995                u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) *      &
996                                       ( u(k,j,nx+1) - u(k,j,nx) ) * ddx
997
998                v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) *      &
999                                       ( v(k,j,nx+1) - v(k,j,nx) ) * ddx
1000
1001                w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) *      &
1002                                       ( w(k,j,nx+1) - w(k,j,nx) ) * ddx
1003             ENDDO
1004          ENDDO
1005
1006!
1007!--       Bottom boundary at the outflow
1008          IF ( ibc_uv_b == 0 )  THEN
1009             u_p(nzb,:,nx+1) = 0.0_wp
1010             v_p(nzb,:,nx+1) = 0.0_wp 
1011          ELSE                   
1012             u_p(nzb,:,nx+1) =  u_p(nzb+1,:,nx+1)
1013             v_p(nzb,:,nx+1) =  v_p(nzb+1,:,nx+1)
1014          ENDIF
1015          w_p(nzb,:,nx+1) = 0.0_wp
1016
1017!
1018!--       Top boundary at the outflow
1019          IF ( ibc_uv_t == 0 )  THEN
1020             u_p(nzt+1,:,nx+1) = u_init(nzt+1)
1021             v_p(nzt+1,:,nx+1) = v_init(nzt+1)
1022          ELSE
1023             u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
1024             v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
1025          ENDIF
1026          w_p(nzt:nzt+1,:,nx+1) = 0.0_wp
1027
1028       ENDIF
1029
1030    ENDIF
1031
1032    IF ( salsa )  THEN
1033       CALL salsa_boundary_conds
1034    ENDIF
1035
1036 END SUBROUTINE boundary_conds
Note: See TracBrowser for help on using the repository browser.