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

Last change on this file since 4180 was 4180, checked in by scharf, 2 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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