[1682] | 1 | !> @file boundary_conds.f90 |
---|
[1036] | 2 | !--------------------------------------------------------------------------------! |
---|
| 3 | ! This file is part of PALM. |
---|
| 4 | ! |
---|
| 5 | ! PALM is free software: you can redistribute it and/or modify it under the terms |
---|
| 6 | ! of the GNU General Public License as published by the Free Software Foundation, |
---|
| 7 | ! either version 3 of the License, or (at your option) any later version. |
---|
| 8 | ! |
---|
| 9 | ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY |
---|
| 10 | ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR |
---|
| 11 | ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. |
---|
| 12 | ! |
---|
| 13 | ! You should have received a copy of the GNU General Public License along with |
---|
| 14 | ! PALM. If not, see <http://www.gnu.org/licenses/>. |
---|
| 15 | ! |
---|
[1818] | 16 | ! Copyright 1997-2016 Leibniz Universitaet Hannover |
---|
[1036] | 17 | !--------------------------------------------------------------------------------! |
---|
| 18 | ! |
---|
[484] | 19 | ! Current revisions: |
---|
[1] | 20 | ! ----------------- |
---|
[1718] | 21 | ! |
---|
[1823] | 22 | ! |
---|
[1321] | 23 | ! Former revisions: |
---|
| 24 | ! ----------------- |
---|
| 25 | ! $Id: boundary_conds.f90 1823 2016-04-07 08:57:52Z raasch $ |
---|
| 26 | ! |
---|
[1823] | 27 | ! 1822 2016-04-07 07:49:42Z hoffmann |
---|
| 28 | ! icloud_scheme removed. microphyisics_seifert added. |
---|
| 29 | ! |
---|
[1765] | 30 | ! 1764 2016-02-28 12:45:19Z raasch |
---|
| 31 | ! index bug for u_p at left outflow removed |
---|
| 32 | ! |
---|
[1763] | 33 | ! 1762 2016-02-25 12:31:13Z hellstea |
---|
| 34 | ! Introduction of nested domain feature |
---|
| 35 | ! |
---|
[1744] | 36 | ! 1742 2016-01-13 09:50:06Z raasch |
---|
| 37 | ! bugfix for outflow Neumann boundary conditions at bottom and top |
---|
| 38 | ! |
---|
[1718] | 39 | ! 1717 2015-11-11 15:09:47Z raasch |
---|
| 40 | ! Bugfix: index error in outflow conditions for left boundary |
---|
| 41 | ! |
---|
[1683] | 42 | ! 1682 2015-10-07 23:56:08Z knoop |
---|
| 43 | ! Code annotations made doxygen readable |
---|
| 44 | ! |
---|
[1717] | 45 | ! 1410 2014-05-23 12:16:18Z suehring |
---|
[1463] | 46 | ! Bugfix: set dirichlet boundary condition for passive_scalar at model domain |
---|
| 47 | ! top |
---|
| 48 | ! |
---|
[1410] | 49 | ! 1399 2014-05-07 11:16:25Z heinze |
---|
| 50 | ! Bugfix: set inflow boundary conditions also if no humidity or passive_scalar |
---|
| 51 | ! is used. |
---|
| 52 | ! |
---|
[1399] | 53 | ! 1398 2014-05-07 11:15:00Z heinze |
---|
| 54 | ! Dirichlet-condition at the top for u and v changed to u_init and v_init also |
---|
| 55 | ! for large_scale_forcing |
---|
| 56 | ! |
---|
[1381] | 57 | ! 1380 2014-04-28 12:40:45Z heinze |
---|
| 58 | ! Adjust Dirichlet-condition at the top for pt in case of nudging |
---|
| 59 | ! |
---|
[1362] | 60 | ! 1361 2014-04-16 15:17:48Z hoffmann |
---|
| 61 | ! Bottom and top boundary conditions of rain water content (qr) and |
---|
| 62 | ! rain drop concentration (nr) changed to Dirichlet |
---|
| 63 | ! |
---|
[1354] | 64 | ! 1353 2014-04-08 15:21:23Z heinze |
---|
| 65 | ! REAL constants provided with KIND-attribute |
---|
| 66 | ! |
---|
[1321] | 67 | ! 1320 2014-03-20 08:40:49Z raasch |
---|
[1320] | 68 | ! ONLY-attribute added to USE-statements, |
---|
| 69 | ! kind-parameters added to all INTEGER and REAL declaration statements, |
---|
| 70 | ! kinds are defined in new module kinds, |
---|
| 71 | ! revision history before 2012 removed, |
---|
| 72 | ! comment fields (!:) to be used for variable explanations added to |
---|
| 73 | ! all variable declaration statements |
---|
[1160] | 74 | ! |
---|
[1258] | 75 | ! 1257 2013-11-08 15:18:40Z raasch |
---|
| 76 | ! loop independent clauses added |
---|
| 77 | ! |
---|
[1242] | 78 | ! 1241 2013-10-30 11:36:58Z heinze |
---|
| 79 | ! Adjust ug and vg at each timestep in case of large_scale_forcing |
---|
| 80 | ! |
---|
[1160] | 81 | ! 1159 2013-05-21 11:58:22Z fricke |
---|
[1159] | 82 | ! Bugfix: Neumann boundary conditions for the velocity components at the |
---|
| 83 | ! outflow are in fact radiation boundary conditions using the maximum phase |
---|
| 84 | ! velocity that ensures numerical stability (CFL-condition). |
---|
| 85 | ! Hence, logical operator use_cmax is now used instead of bc_lr_dirneu/_neudir. |
---|
| 86 | ! Bugfix: In case of use_cmax at the outflow, u, v, w are replaced by |
---|
| 87 | ! u_p, v_p, w_p |
---|
[1116] | 88 | ! |
---|
| 89 | ! 1115 2013-03-26 18:16:16Z hoffmann |
---|
| 90 | ! boundary conditions of two-moment cloud scheme are restricted to Neumann- |
---|
| 91 | ! boundary-conditions |
---|
| 92 | ! |
---|
[1114] | 93 | ! 1113 2013-03-10 02:48:14Z raasch |
---|
| 94 | ! GPU-porting |
---|
| 95 | ! dummy argument "range" removed |
---|
| 96 | ! Bugfix: wrong index in loops of radiation boundary condition |
---|
[1113] | 97 | ! |
---|
[1054] | 98 | ! 1053 2012-11-13 17:11:03Z hoffmann |
---|
| 99 | ! boundary conditions for the two new prognostic equations (nr, qr) of the |
---|
| 100 | ! two-moment cloud scheme |
---|
| 101 | ! |
---|
[1037] | 102 | ! 1036 2012-10-22 13:43:42Z raasch |
---|
| 103 | ! code put under GPL (PALM 3.9) |
---|
| 104 | ! |
---|
[997] | 105 | ! 996 2012-09-07 10:41:47Z raasch |
---|
| 106 | ! little reformatting |
---|
| 107 | ! |
---|
[979] | 108 | ! 978 2012-08-09 08:28:32Z fricke |
---|
| 109 | ! Neumann boudnary conditions are added at the inflow boundary for the SGS-TKE. |
---|
| 110 | ! Outflow boundary conditions for the velocity components can be set to Neumann |
---|
| 111 | ! conditions or to radiation conditions with a horizontal averaged phase |
---|
| 112 | ! velocity. |
---|
| 113 | ! |
---|
[876] | 114 | ! 875 2012-04-02 15:35:15Z gryschka |
---|
| 115 | ! Bugfix in case of dirichlet inflow bc at the right or north boundary |
---|
| 116 | ! |
---|
[1] | 117 | ! Revision 1.1 1997/09/12 06:21:34 raasch |
---|
| 118 | ! Initial revision |
---|
| 119 | ! |
---|
| 120 | ! |
---|
| 121 | ! Description: |
---|
| 122 | ! ------------ |
---|
[1682] | 123 | !> Boundary conditions for the prognostic quantities. |
---|
| 124 | !> One additional bottom boundary condition is applied for the TKE (=(u*)**2) |
---|
| 125 | !> in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly |
---|
| 126 | !> handled in routine exchange_horiz. Pressure boundary conditions are |
---|
| 127 | !> explicitly set in routines pres, poisfft, poismg and sor. |
---|
[1] | 128 | !------------------------------------------------------------------------------! |
---|
[1682] | 129 | SUBROUTINE boundary_conds |
---|
| 130 | |
---|
[1] | 131 | |
---|
[1320] | 132 | USE arrays_3d, & |
---|
| 133 | 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, & |
---|
| 134 | dzu, e_p, nr_p, pt, pt_p, q, q_p, qr_p, sa, sa_p, & |
---|
| 135 | u, ug, u_init, u_m_l, u_m_n, u_m_r, u_m_s, u_p, & |
---|
| 136 | v, vg, v_init, v_m_l, v_m_n, v_m_r, v_m_s, v_p, & |
---|
[1380] | 137 | w, w_p, w_m_l, w_m_n, w_m_r, w_m_s,& |
---|
| 138 | pt_init |
---|
[1320] | 139 | |
---|
| 140 | USE control_parameters, & |
---|
| 141 | ONLY: bc_pt_t_val, bc_q_t_val, constant_diffusion, & |
---|
| 142 | cloud_physics, dt_3d, humidity, & |
---|
[1462] | 143 | ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_sa_t, ibc_uv_b, & |
---|
[1822] | 144 | ibc_uv_t, inflow_l, inflow_n, inflow_r, inflow_s, & |
---|
| 145 | intermediate_timestep_count, large_scale_forcing, & |
---|
| 146 | microphysics_seifert, nest_domain, nest_bound_l, nest_bound_s, & |
---|
| 147 | nudging, ocean, outflow_l, outflow_n, outflow_r, outflow_s, & |
---|
| 148 | passive_scalar, tsc, use_cmax |
---|
[1320] | 149 | |
---|
| 150 | USE grid_variables, & |
---|
| 151 | ONLY: ddx, ddy, dx, dy |
---|
| 152 | |
---|
| 153 | USE indices, & |
---|
| 154 | ONLY: nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, & |
---|
| 155 | nzb, nzb_s_inner, nzb_w_inner, nzt |
---|
| 156 | |
---|
| 157 | USE kinds |
---|
| 158 | |
---|
[1] | 159 | USE pegrid |
---|
| 160 | |
---|
[1320] | 161 | |
---|
[1] | 162 | IMPLICIT NONE |
---|
| 163 | |
---|
[1682] | 164 | INTEGER(iwp) :: i !< |
---|
| 165 | INTEGER(iwp) :: j !< |
---|
| 166 | INTEGER(iwp) :: k !< |
---|
[1] | 167 | |
---|
[1682] | 168 | REAL(wp) :: c_max !< |
---|
| 169 | REAL(wp) :: denom !< |
---|
[1] | 170 | |
---|
[73] | 171 | |
---|
[1] | 172 | ! |
---|
[1113] | 173 | !-- Bottom boundary |
---|
| 174 | IF ( ibc_uv_b == 1 ) THEN |
---|
| 175 | !$acc kernels present( u_p, v_p ) |
---|
| 176 | u_p(nzb,:,:) = u_p(nzb+1,:,:) |
---|
| 177 | v_p(nzb,:,:) = v_p(nzb+1,:,:) |
---|
| 178 | !$acc end kernels |
---|
| 179 | ENDIF |
---|
| 180 | |
---|
| 181 | !$acc kernels present( nzb_w_inner, w_p ) |
---|
| 182 | DO i = nxlg, nxrg |
---|
| 183 | DO j = nysg, nyng |
---|
[1353] | 184 | w_p(nzb_w_inner(j,i),j,i) = 0.0_wp |
---|
[1113] | 185 | ENDDO |
---|
| 186 | ENDDO |
---|
| 187 | !$acc end kernels |
---|
| 188 | |
---|
| 189 | ! |
---|
[1762] | 190 | !-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings. |
---|
[1113] | 191 | IF ( ibc_uv_t == 0 ) THEN |
---|
| 192 | !$acc kernels present( u_init, u_p, v_init, v_p ) |
---|
| 193 | u_p(nzt+1,:,:) = u_init(nzt+1) |
---|
| 194 | v_p(nzt+1,:,:) = v_init(nzt+1) |
---|
| 195 | !$acc end kernels |
---|
[1762] | 196 | ELSEIF ( ibc_uv_t == 1 ) THEN |
---|
[1113] | 197 | !$acc kernels present( u_p, v_p ) |
---|
| 198 | u_p(nzt+1,:,:) = u_p(nzt,:,:) |
---|
| 199 | v_p(nzt+1,:,:) = v_p(nzt,:,:) |
---|
| 200 | !$acc end kernels |
---|
| 201 | ENDIF |
---|
| 202 | |
---|
[1762] | 203 | IF ( .NOT. nest_domain ) THEN |
---|
| 204 | !$acc kernels present( w_p ) |
---|
| 205 | w_p(nzt:nzt+1,:,:) = 0.0_wp ! nzt is not a prognostic level (but cf. pres) |
---|
| 206 | !$acc end kernels |
---|
| 207 | ENDIF |
---|
| 208 | |
---|
[1113] | 209 | ! |
---|
| 210 | !-- Temperature at bottom boundary. |
---|
| 211 | !-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by |
---|
| 212 | !-- the sea surface temperature of the coupled ocean model. |
---|
| 213 | IF ( ibc_pt_b == 0 ) THEN |
---|
| 214 | !$acc kernels present( nzb_s_inner, pt, pt_p ) |
---|
[1257] | 215 | !$acc loop independent |
---|
[667] | 216 | DO i = nxlg, nxrg |
---|
[1257] | 217 | !$acc loop independent |
---|
[667] | 218 | DO j = nysg, nyng |
---|
[1113] | 219 | pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i) |
---|
[1] | 220 | ENDDO |
---|
| 221 | ENDDO |
---|
[1113] | 222 | !$acc end kernels |
---|
| 223 | ELSEIF ( ibc_pt_b == 1 ) THEN |
---|
| 224 | !$acc kernels present( nzb_s_inner, pt_p ) |
---|
[1257] | 225 | !$acc loop independent |
---|
[1113] | 226 | DO i = nxlg, nxrg |
---|
[1257] | 227 | !$acc loop independent |
---|
[1113] | 228 | DO j = nysg, nyng |
---|
| 229 | pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i) |
---|
| 230 | ENDDO |
---|
| 231 | ENDDO |
---|
| 232 | !$acc end kernels |
---|
| 233 | ENDIF |
---|
[1] | 234 | |
---|
| 235 | ! |
---|
[1113] | 236 | !-- Temperature at top boundary |
---|
| 237 | IF ( ibc_pt_t == 0 ) THEN |
---|
| 238 | !$acc kernels present( pt, pt_p ) |
---|
| 239 | pt_p(nzt+1,:,:) = pt(nzt+1,:,:) |
---|
[1380] | 240 | ! |
---|
| 241 | !-- In case of nudging adjust top boundary to pt which is |
---|
| 242 | !-- read in from NUDGING-DATA |
---|
| 243 | IF ( nudging ) THEN |
---|
| 244 | pt_p(nzt+1,:,:) = pt_init(nzt+1) |
---|
| 245 | ENDIF |
---|
[1113] | 246 | !$acc end kernels |
---|
| 247 | ELSEIF ( ibc_pt_t == 1 ) THEN |
---|
| 248 | !$acc kernels present( pt_p ) |
---|
| 249 | pt_p(nzt+1,:,:) = pt_p(nzt,:,:) |
---|
| 250 | !$acc end kernels |
---|
| 251 | ELSEIF ( ibc_pt_t == 2 ) THEN |
---|
| 252 | !$acc kernels present( dzu, pt_p ) |
---|
| 253 | pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1) |
---|
| 254 | !$acc end kernels |
---|
| 255 | ENDIF |
---|
[1] | 256 | |
---|
| 257 | ! |
---|
[1113] | 258 | !-- Boundary conditions for TKE |
---|
| 259 | !-- Generally Neumann conditions with de/dz=0 are assumed |
---|
| 260 | IF ( .NOT. constant_diffusion ) THEN |
---|
| 261 | !$acc kernels present( e_p, nzb_s_inner ) |
---|
[1257] | 262 | !$acc loop independent |
---|
[1113] | 263 | DO i = nxlg, nxrg |
---|
[1257] | 264 | !$acc loop independent |
---|
[1113] | 265 | DO j = nysg, nyng |
---|
| 266 | e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i) |
---|
[73] | 267 | ENDDO |
---|
[1113] | 268 | ENDDO |
---|
[1762] | 269 | IF ( .NOT. nest_domain ) THEN |
---|
| 270 | e_p(nzt+1,:,:) = e_p(nzt,:,:) |
---|
| 271 | ENDIF |
---|
[1113] | 272 | !$acc end kernels |
---|
| 273 | ENDIF |
---|
| 274 | |
---|
| 275 | ! |
---|
| 276 | !-- Boundary conditions for salinity |
---|
| 277 | IF ( ocean ) THEN |
---|
| 278 | ! |
---|
| 279 | !-- Bottom boundary: Neumann condition because salinity flux is always |
---|
| 280 | !-- given |
---|
| 281 | DO i = nxlg, nxrg |
---|
| 282 | DO j = nysg, nyng |
---|
| 283 | sa_p(nzb_s_inner(j,i),j,i) = sa_p(nzb_s_inner(j,i)+1,j,i) |
---|
[1] | 284 | ENDDO |
---|
[1113] | 285 | ENDDO |
---|
[1] | 286 | |
---|
| 287 | ! |
---|
[1113] | 288 | !-- Top boundary: Dirichlet or Neumann |
---|
| 289 | IF ( ibc_sa_t == 0 ) THEN |
---|
| 290 | sa_p(nzt+1,:,:) = sa(nzt+1,:,:) |
---|
| 291 | ELSEIF ( ibc_sa_t == 1 ) THEN |
---|
| 292 | sa_p(nzt+1,:,:) = sa_p(nzt,:,:) |
---|
[1] | 293 | ENDIF |
---|
| 294 | |
---|
[1113] | 295 | ENDIF |
---|
| 296 | |
---|
[1] | 297 | ! |
---|
[1113] | 298 | !-- Boundary conditions for total water content or scalar, |
---|
| 299 | !-- bottom and top boundary (see also temperature) |
---|
| 300 | IF ( humidity .OR. passive_scalar ) THEN |
---|
| 301 | ! |
---|
| 302 | !-- Surface conditions for constant_humidity_flux |
---|
| 303 | IF ( ibc_q_b == 0 ) THEN |
---|
[667] | 304 | DO i = nxlg, nxrg |
---|
| 305 | DO j = nysg, nyng |
---|
[1113] | 306 | q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i) |
---|
[1] | 307 | ENDDO |
---|
| 308 | ENDDO |
---|
[1113] | 309 | ELSE |
---|
[667] | 310 | DO i = nxlg, nxrg |
---|
| 311 | DO j = nysg, nyng |
---|
[1113] | 312 | q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i) |
---|
[95] | 313 | ENDDO |
---|
| 314 | ENDDO |
---|
[1113] | 315 | ENDIF |
---|
[95] | 316 | ! |
---|
[1113] | 317 | !-- Top boundary |
---|
[1462] | 318 | IF ( ibc_q_t == 0 ) THEN |
---|
| 319 | q_p(nzt+1,:,:) = q(nzt+1,:,:) |
---|
| 320 | ELSEIF ( ibc_q_t == 1 ) THEN |
---|
| 321 | q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1) |
---|
| 322 | ENDIF |
---|
[95] | 323 | |
---|
[1822] | 324 | IF ( cloud_physics .AND. microphysics_seifert ) THEN |
---|
[1113] | 325 | ! |
---|
[1361] | 326 | !-- Surface conditions rain water (Dirichlet) |
---|
[1115] | 327 | DO i = nxlg, nxrg |
---|
| 328 | DO j = nysg, nyng |
---|
[1361] | 329 | qr_p(nzb_s_inner(j,i),j,i) = 0.0_wp |
---|
| 330 | nr_p(nzb_s_inner(j,i),j,i) = 0.0_wp |
---|
[73] | 331 | ENDDO |
---|
[1115] | 332 | ENDDO |
---|
[1] | 333 | ! |
---|
[1361] | 334 | !-- Top boundary condition for rain water (Dirichlet) |
---|
| 335 | qr_p(nzt+1,:,:) = 0.0_wp |
---|
| 336 | nr_p(nzt+1,:,:) = 0.0_wp |
---|
[1115] | 337 | |
---|
[1] | 338 | ENDIF |
---|
[1409] | 339 | ENDIF |
---|
[1] | 340 | ! |
---|
[1762] | 341 | !-- In case of inflow or nest boundary at the south boundary the boundary for v |
---|
| 342 | !-- is at nys and in case of inflow or nest boundary at the left boundary the |
---|
| 343 | !-- boundary for u is at nxl. Since in prognostic_equations (cache optimized |
---|
| 344 | !-- version) these levels are handled as a prognostic level, boundary values |
---|
| 345 | !-- have to be restored here. |
---|
[1409] | 346 | !-- For the SGS-TKE, Neumann boundary conditions are used at the inflow. |
---|
| 347 | IF ( inflow_s ) THEN |
---|
| 348 | v_p(:,nys,:) = v_p(:,nys-1,:) |
---|
| 349 | IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:) |
---|
| 350 | ELSEIF ( inflow_n ) THEN |
---|
| 351 | IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:) |
---|
| 352 | ELSEIF ( inflow_l ) THEN |
---|
| 353 | u_p(:,:,nxl) = u_p(:,:,nxl-1) |
---|
| 354 | IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl) |
---|
| 355 | ELSEIF ( inflow_r ) THEN |
---|
| 356 | IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr) |
---|
| 357 | ENDIF |
---|
[1] | 358 | |
---|
| 359 | ! |
---|
[1762] | 360 | !-- The same restoration for u at i=nxl and v at j=nys as above must be made |
---|
| 361 | !-- in case of nest boundaries. Note however, that the above ELSEIF-structure is |
---|
| 362 | !-- not appropriate here as there may be more than one nest boundary on a |
---|
| 363 | !-- PE-domain. Furthermore Neumann conditions for SGS-TKE are not required here. |
---|
| 364 | IF ( nest_bound_s ) THEN |
---|
| 365 | v_p(:,nys,:) = v_p(:,nys-1,:) |
---|
| 366 | ENDIF |
---|
| 367 | IF ( nest_bound_l ) THEN |
---|
| 368 | u_p(:,:,nxl) = u_p(:,:,nxl-1) |
---|
| 369 | ENDIF |
---|
| 370 | |
---|
| 371 | ! |
---|
[1409] | 372 | !-- Lateral boundary conditions for scalar quantities at the outflow |
---|
| 373 | IF ( outflow_s ) THEN |
---|
| 374 | pt_p(:,nys-1,:) = pt_p(:,nys,:) |
---|
| 375 | IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:) |
---|
| 376 | IF ( humidity .OR. passive_scalar ) THEN |
---|
| 377 | q_p(:,nys-1,:) = q_p(:,nys,:) |
---|
[1822] | 378 | IF ( cloud_physics .AND. microphysics_seifert ) THEN |
---|
[1409] | 379 | qr_p(:,nys-1,:) = qr_p(:,nys,:) |
---|
| 380 | nr_p(:,nys-1,:) = nr_p(:,nys,:) |
---|
[1053] | 381 | ENDIF |
---|
[1409] | 382 | ENDIF |
---|
| 383 | ELSEIF ( outflow_n ) THEN |
---|
| 384 | pt_p(:,nyn+1,:) = pt_p(:,nyn,:) |
---|
| 385 | IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:) |
---|
| 386 | IF ( humidity .OR. passive_scalar ) THEN |
---|
| 387 | q_p(:,nyn+1,:) = q_p(:,nyn,:) |
---|
[1822] | 388 | IF ( cloud_physics .AND. microphysics_seifert ) THEN |
---|
[1409] | 389 | qr_p(:,nyn+1,:) = qr_p(:,nyn,:) |
---|
| 390 | nr_p(:,nyn+1,:) = nr_p(:,nyn,:) |
---|
[1053] | 391 | ENDIF |
---|
[1409] | 392 | ENDIF |
---|
| 393 | ELSEIF ( outflow_l ) THEN |
---|
| 394 | pt_p(:,:,nxl-1) = pt_p(:,:,nxl) |
---|
| 395 | IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl) |
---|
| 396 | IF ( humidity .OR. passive_scalar ) THEN |
---|
| 397 | q_p(:,:,nxl-1) = q_p(:,:,nxl) |
---|
[1822] | 398 | IF ( cloud_physics .AND. microphysics_seifert ) THEN |
---|
[1409] | 399 | qr_p(:,:,nxl-1) = qr_p(:,:,nxl) |
---|
| 400 | nr_p(:,:,nxl-1) = nr_p(:,:,nxl) |
---|
[1053] | 401 | ENDIF |
---|
[1409] | 402 | ENDIF |
---|
| 403 | ELSEIF ( outflow_r ) THEN |
---|
| 404 | pt_p(:,:,nxr+1) = pt_p(:,:,nxr) |
---|
| 405 | IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr) |
---|
| 406 | IF ( humidity .OR. passive_scalar ) THEN |
---|
| 407 | q_p(:,:,nxr+1) = q_p(:,:,nxr) |
---|
[1822] | 408 | IF ( cloud_physics .AND. microphysics_seifert ) THEN |
---|
[1409] | 409 | qr_p(:,:,nxr+1) = qr_p(:,:,nxr) |
---|
| 410 | nr_p(:,:,nxr+1) = nr_p(:,:,nxr) |
---|
[1053] | 411 | ENDIF |
---|
[1] | 412 | ENDIF |
---|
| 413 | ENDIF |
---|
| 414 | |
---|
| 415 | ! |
---|
[1159] | 416 | !-- Radiation boundary conditions for the velocities at the respective outflow. |
---|
| 417 | !-- The phase velocity is either assumed to the maximum phase velocity that |
---|
| 418 | !-- ensures numerical stability (CFL-condition) or calculated after |
---|
| 419 | !-- Orlanski(1976) and averaged along the outflow boundary. |
---|
[106] | 420 | IF ( outflow_s ) THEN |
---|
[75] | 421 | |
---|
[1159] | 422 | IF ( use_cmax ) THEN |
---|
| 423 | u_p(:,-1,:) = u(:,0,:) |
---|
| 424 | v_p(:,0,:) = v(:,1,:) |
---|
| 425 | w_p(:,-1,:) = w(:,0,:) |
---|
| 426 | ELSEIF ( .NOT. use_cmax ) THEN |
---|
[75] | 427 | |
---|
[978] | 428 | c_max = dy / dt_3d |
---|
[75] | 429 | |
---|
[1353] | 430 | c_u_m_l = 0.0_wp |
---|
| 431 | c_v_m_l = 0.0_wp |
---|
| 432 | c_w_m_l = 0.0_wp |
---|
[978] | 433 | |
---|
[1353] | 434 | c_u_m = 0.0_wp |
---|
| 435 | c_v_m = 0.0_wp |
---|
| 436 | c_w_m = 0.0_wp |
---|
[978] | 437 | |
---|
[75] | 438 | ! |
---|
[996] | 439 | !-- Calculate the phase speeds for u, v, and w, first local and then |
---|
| 440 | !-- average along the outflow boundary. |
---|
| 441 | DO k = nzb+1, nzt+1 |
---|
| 442 | DO i = nxl, nxr |
---|
[75] | 443 | |
---|
[106] | 444 | denom = u_m_s(k,0,i) - u_m_s(k,1,i) |
---|
| 445 | |
---|
[1353] | 446 | IF ( denom /= 0.0_wp ) THEN |
---|
[996] | 447 | c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / ( denom * tsc(2) ) |
---|
[1353] | 448 | IF ( c_u(k,i) < 0.0_wp ) THEN |
---|
| 449 | c_u(k,i) = 0.0_wp |
---|
[106] | 450 | ELSEIF ( c_u(k,i) > c_max ) THEN |
---|
| 451 | c_u(k,i) = c_max |
---|
| 452 | ENDIF |
---|
| 453 | ELSE |
---|
| 454 | c_u(k,i) = c_max |
---|
[75] | 455 | ENDIF |
---|
| 456 | |
---|
[106] | 457 | denom = v_m_s(k,1,i) - v_m_s(k,2,i) |
---|
| 458 | |
---|
[1353] | 459 | IF ( denom /= 0.0_wp ) THEN |
---|
[996] | 460 | c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / ( denom * tsc(2) ) |
---|
[1353] | 461 | IF ( c_v(k,i) < 0.0_wp ) THEN |
---|
| 462 | c_v(k,i) = 0.0_wp |
---|
[106] | 463 | ELSEIF ( c_v(k,i) > c_max ) THEN |
---|
| 464 | c_v(k,i) = c_max |
---|
| 465 | ENDIF |
---|
| 466 | ELSE |
---|
| 467 | c_v(k,i) = c_max |
---|
[75] | 468 | ENDIF |
---|
| 469 | |
---|
[106] | 470 | denom = w_m_s(k,0,i) - w_m_s(k,1,i) |
---|
[75] | 471 | |
---|
[1353] | 472 | IF ( denom /= 0.0_wp ) THEN |
---|
[996] | 473 | c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) ) |
---|
[1353] | 474 | IF ( c_w(k,i) < 0.0_wp ) THEN |
---|
| 475 | c_w(k,i) = 0.0_wp |
---|
[106] | 476 | ELSEIF ( c_w(k,i) > c_max ) THEN |
---|
| 477 | c_w(k,i) = c_max |
---|
| 478 | ENDIF |
---|
| 479 | ELSE |
---|
| 480 | c_w(k,i) = c_max |
---|
[75] | 481 | ENDIF |
---|
[106] | 482 | |
---|
[978] | 483 | c_u_m_l(k) = c_u_m_l(k) + c_u(k,i) |
---|
| 484 | c_v_m_l(k) = c_v_m_l(k) + c_v(k,i) |
---|
| 485 | c_w_m_l(k) = c_w_m_l(k) + c_w(k,i) |
---|
[106] | 486 | |
---|
[978] | 487 | ENDDO |
---|
| 488 | ENDDO |
---|
[75] | 489 | |
---|
[978] | 490 | #if defined( __parallel ) |
---|
| 491 | IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) |
---|
| 492 | CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, & |
---|
| 493 | MPI_SUM, comm1dx, ierr ) |
---|
| 494 | IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) |
---|
| 495 | CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, & |
---|
| 496 | MPI_SUM, comm1dx, ierr ) |
---|
| 497 | IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) |
---|
| 498 | CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, & |
---|
| 499 | MPI_SUM, comm1dx, ierr ) |
---|
| 500 | #else |
---|
| 501 | c_u_m = c_u_m_l |
---|
| 502 | c_v_m = c_v_m_l |
---|
| 503 | c_w_m = c_w_m_l |
---|
| 504 | #endif |
---|
| 505 | |
---|
| 506 | c_u_m = c_u_m / (nx+1) |
---|
| 507 | c_v_m = c_v_m / (nx+1) |
---|
| 508 | c_w_m = c_w_m / (nx+1) |
---|
| 509 | |
---|
[75] | 510 | ! |
---|
[978] | 511 | !-- Save old timelevels for the next timestep |
---|
| 512 | IF ( intermediate_timestep_count == 1 ) THEN |
---|
| 513 | u_m_s(:,:,:) = u(:,0:1,:) |
---|
| 514 | v_m_s(:,:,:) = v(:,1:2,:) |
---|
| 515 | w_m_s(:,:,:) = w(:,0:1,:) |
---|
| 516 | ENDIF |
---|
| 517 | |
---|
| 518 | ! |
---|
| 519 | !-- Calculate the new velocities |
---|
[996] | 520 | DO k = nzb+1, nzt+1 |
---|
| 521 | DO i = nxlg, nxrg |
---|
[978] | 522 | u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) * & |
---|
[75] | 523 | ( u(k,-1,i) - u(k,0,i) ) * ddy |
---|
| 524 | |
---|
[978] | 525 | v_p(k,0,i) = v(k,0,i) - dt_3d * tsc(2) * c_v_m(k) * & |
---|
[106] | 526 | ( v(k,0,i) - v(k,1,i) ) * ddy |
---|
[75] | 527 | |
---|
[978] | 528 | w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) * & |
---|
[75] | 529 | ( w(k,-1,i) - w(k,0,i) ) * ddy |
---|
[978] | 530 | ENDDO |
---|
[75] | 531 | ENDDO |
---|
| 532 | |
---|
| 533 | ! |
---|
[978] | 534 | !-- Bottom boundary at the outflow |
---|
| 535 | IF ( ibc_uv_b == 0 ) THEN |
---|
[1353] | 536 | u_p(nzb,-1,:) = 0.0_wp |
---|
| 537 | v_p(nzb,0,:) = 0.0_wp |
---|
[978] | 538 | ELSE |
---|
| 539 | u_p(nzb,-1,:) = u_p(nzb+1,-1,:) |
---|
| 540 | v_p(nzb,0,:) = v_p(nzb+1,0,:) |
---|
| 541 | ENDIF |
---|
[1353] | 542 | w_p(nzb,-1,:) = 0.0_wp |
---|
[73] | 543 | |
---|
[75] | 544 | ! |
---|
[978] | 545 | !-- Top boundary at the outflow |
---|
| 546 | IF ( ibc_uv_t == 0 ) THEN |
---|
| 547 | u_p(nzt+1,-1,:) = u_init(nzt+1) |
---|
| 548 | v_p(nzt+1,0,:) = v_init(nzt+1) |
---|
| 549 | ELSE |
---|
[1742] | 550 | u_p(nzt+1,-1,:) = u_p(nzt,-1,:) |
---|
| 551 | v_p(nzt+1,0,:) = v_p(nzt,0,:) |
---|
[978] | 552 | ENDIF |
---|
[1353] | 553 | w_p(nzt:nzt+1,-1,:) = 0.0_wp |
---|
[978] | 554 | |
---|
[75] | 555 | ENDIF |
---|
[73] | 556 | |
---|
[75] | 557 | ENDIF |
---|
[73] | 558 | |
---|
[106] | 559 | IF ( outflow_n ) THEN |
---|
[73] | 560 | |
---|
[1159] | 561 | IF ( use_cmax ) THEN |
---|
| 562 | u_p(:,ny+1,:) = u(:,ny,:) |
---|
| 563 | v_p(:,ny+1,:) = v(:,ny,:) |
---|
| 564 | w_p(:,ny+1,:) = w(:,ny,:) |
---|
| 565 | ELSEIF ( .NOT. use_cmax ) THEN |
---|
[75] | 566 | |
---|
[978] | 567 | c_max = dy / dt_3d |
---|
[75] | 568 | |
---|
[1353] | 569 | c_u_m_l = 0.0_wp |
---|
| 570 | c_v_m_l = 0.0_wp |
---|
| 571 | c_w_m_l = 0.0_wp |
---|
[978] | 572 | |
---|
[1353] | 573 | c_u_m = 0.0_wp |
---|
| 574 | c_v_m = 0.0_wp |
---|
| 575 | c_w_m = 0.0_wp |
---|
[978] | 576 | |
---|
[1] | 577 | ! |
---|
[996] | 578 | !-- Calculate the phase speeds for u, v, and w, first local and then |
---|
| 579 | !-- average along the outflow boundary. |
---|
| 580 | DO k = nzb+1, nzt+1 |
---|
| 581 | DO i = nxl, nxr |
---|
[73] | 582 | |
---|
[106] | 583 | denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i) |
---|
| 584 | |
---|
[1353] | 585 | IF ( denom /= 0.0_wp ) THEN |
---|
[996] | 586 | c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) ) |
---|
[1353] | 587 | IF ( c_u(k,i) < 0.0_wp ) THEN |
---|
| 588 | c_u(k,i) = 0.0_wp |
---|
[106] | 589 | ELSEIF ( c_u(k,i) > c_max ) THEN |
---|
| 590 | c_u(k,i) = c_max |
---|
| 591 | ENDIF |
---|
| 592 | ELSE |
---|
| 593 | c_u(k,i) = c_max |
---|
[73] | 594 | ENDIF |
---|
| 595 | |
---|
[106] | 596 | denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i) |
---|
[73] | 597 | |
---|
[1353] | 598 | IF ( denom /= 0.0_wp ) THEN |
---|
[996] | 599 | c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) ) |
---|
[1353] | 600 | IF ( c_v(k,i) < 0.0_wp ) THEN |
---|
| 601 | c_v(k,i) = 0.0_wp |
---|
[106] | 602 | ELSEIF ( c_v(k,i) > c_max ) THEN |
---|
| 603 | c_v(k,i) = c_max |
---|
| 604 | ENDIF |
---|
| 605 | ELSE |
---|
| 606 | c_v(k,i) = c_max |
---|
[73] | 607 | ENDIF |
---|
| 608 | |
---|
[106] | 609 | denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i) |
---|
[73] | 610 | |
---|
[1353] | 611 | IF ( denom /= 0.0_wp ) THEN |
---|
[996] | 612 | c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) ) |
---|
[1353] | 613 | IF ( c_w(k,i) < 0.0_wp ) THEN |
---|
| 614 | c_w(k,i) = 0.0_wp |
---|
[106] | 615 | ELSEIF ( c_w(k,i) > c_max ) THEN |
---|
| 616 | c_w(k,i) = c_max |
---|
| 617 | ENDIF |
---|
| 618 | ELSE |
---|
| 619 | c_w(k,i) = c_max |
---|
[73] | 620 | ENDIF |
---|
[106] | 621 | |
---|
[978] | 622 | c_u_m_l(k) = c_u_m_l(k) + c_u(k,i) |
---|
| 623 | c_v_m_l(k) = c_v_m_l(k) + c_v(k,i) |
---|
| 624 | c_w_m_l(k) = c_w_m_l(k) + c_w(k,i) |
---|
[106] | 625 | |
---|
[978] | 626 | ENDDO |
---|
| 627 | ENDDO |
---|
[73] | 628 | |
---|
[978] | 629 | #if defined( __parallel ) |
---|
| 630 | IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) |
---|
| 631 | CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, & |
---|
| 632 | MPI_SUM, comm1dx, ierr ) |
---|
| 633 | IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) |
---|
| 634 | CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, & |
---|
| 635 | MPI_SUM, comm1dx, ierr ) |
---|
| 636 | IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) |
---|
| 637 | CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, & |
---|
| 638 | MPI_SUM, comm1dx, ierr ) |
---|
| 639 | #else |
---|
| 640 | c_u_m = c_u_m_l |
---|
| 641 | c_v_m = c_v_m_l |
---|
| 642 | c_w_m = c_w_m_l |
---|
| 643 | #endif |
---|
| 644 | |
---|
| 645 | c_u_m = c_u_m / (nx+1) |
---|
| 646 | c_v_m = c_v_m / (nx+1) |
---|
| 647 | c_w_m = c_w_m / (nx+1) |
---|
| 648 | |
---|
[73] | 649 | ! |
---|
[978] | 650 | !-- Save old timelevels for the next timestep |
---|
| 651 | IF ( intermediate_timestep_count == 1 ) THEN |
---|
| 652 | u_m_n(:,:,:) = u(:,ny-1:ny,:) |
---|
| 653 | v_m_n(:,:,:) = v(:,ny-1:ny,:) |
---|
| 654 | w_m_n(:,:,:) = w(:,ny-1:ny,:) |
---|
| 655 | ENDIF |
---|
[73] | 656 | |
---|
[978] | 657 | ! |
---|
| 658 | !-- Calculate the new velocities |
---|
[996] | 659 | DO k = nzb+1, nzt+1 |
---|
| 660 | DO i = nxlg, nxrg |
---|
[978] | 661 | u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) * & |
---|
| 662 | ( u(k,ny+1,i) - u(k,ny,i) ) * ddy |
---|
[73] | 663 | |
---|
[978] | 664 | v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v_m(k) * & |
---|
| 665 | ( v(k,ny+1,i) - v(k,ny,i) ) * ddy |
---|
[73] | 666 | |
---|
[978] | 667 | w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) * & |
---|
| 668 | ( w(k,ny+1,i) - w(k,ny,i) ) * ddy |
---|
| 669 | ENDDO |
---|
[1] | 670 | ENDDO |
---|
| 671 | |
---|
| 672 | ! |
---|
[978] | 673 | !-- Bottom boundary at the outflow |
---|
| 674 | IF ( ibc_uv_b == 0 ) THEN |
---|
[1353] | 675 | u_p(nzb,ny+1,:) = 0.0_wp |
---|
| 676 | v_p(nzb,ny+1,:) = 0.0_wp |
---|
[978] | 677 | ELSE |
---|
| 678 | u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:) |
---|
| 679 | v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:) |
---|
| 680 | ENDIF |
---|
[1353] | 681 | w_p(nzb,ny+1,:) = 0.0_wp |
---|
[73] | 682 | |
---|
| 683 | ! |
---|
[978] | 684 | !-- Top boundary at the outflow |
---|
| 685 | IF ( ibc_uv_t == 0 ) THEN |
---|
| 686 | u_p(nzt+1,ny+1,:) = u_init(nzt+1) |
---|
| 687 | v_p(nzt+1,ny+1,:) = v_init(nzt+1) |
---|
| 688 | ELSE |
---|
| 689 | u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:) |
---|
| 690 | v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:) |
---|
| 691 | ENDIF |
---|
[1353] | 692 | w_p(nzt:nzt+1,ny+1,:) = 0.0_wp |
---|
[978] | 693 | |
---|
[1] | 694 | ENDIF |
---|
| 695 | |
---|
[75] | 696 | ENDIF |
---|
| 697 | |
---|
[106] | 698 | IF ( outflow_l ) THEN |
---|
[75] | 699 | |
---|
[1159] | 700 | IF ( use_cmax ) THEN |
---|
[1717] | 701 | u_p(:,:,0) = u(:,:,1) |
---|
| 702 | v_p(:,:,-1) = v(:,:,0) |
---|
[1159] | 703 | w_p(:,:,-1) = w(:,:,0) |
---|
| 704 | ELSEIF ( .NOT. use_cmax ) THEN |
---|
[75] | 705 | |
---|
[978] | 706 | c_max = dx / dt_3d |
---|
[75] | 707 | |
---|
[1353] | 708 | c_u_m_l = 0.0_wp |
---|
| 709 | c_v_m_l = 0.0_wp |
---|
| 710 | c_w_m_l = 0.0_wp |
---|
[978] | 711 | |
---|
[1353] | 712 | c_u_m = 0.0_wp |
---|
| 713 | c_v_m = 0.0_wp |
---|
| 714 | c_w_m = 0.0_wp |
---|
[978] | 715 | |
---|
[1] | 716 | ! |
---|
[996] | 717 | !-- Calculate the phase speeds for u, v, and w, first local and then |
---|
| 718 | !-- average along the outflow boundary. |
---|
| 719 | DO k = nzb+1, nzt+1 |
---|
| 720 | DO j = nys, nyn |
---|
[75] | 721 | |
---|
[106] | 722 | denom = u_m_l(k,j,1) - u_m_l(k,j,2) |
---|
| 723 | |
---|
[1353] | 724 | IF ( denom /= 0.0_wp ) THEN |
---|
[996] | 725 | c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) ) |
---|
[1353] | 726 | IF ( c_u(k,j) < 0.0_wp ) THEN |
---|
| 727 | c_u(k,j) = 0.0_wp |
---|
[107] | 728 | ELSEIF ( c_u(k,j) > c_max ) THEN |
---|
| 729 | c_u(k,j) = c_max |
---|
[106] | 730 | ENDIF |
---|
| 731 | ELSE |
---|
[107] | 732 | c_u(k,j) = c_max |
---|
[75] | 733 | ENDIF |
---|
| 734 | |
---|
[106] | 735 | denom = v_m_l(k,j,0) - v_m_l(k,j,1) |
---|
[75] | 736 | |
---|
[1353] | 737 | IF ( denom /= 0.0_wp ) THEN |
---|
[996] | 738 | c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) ) |
---|
[1353] | 739 | IF ( c_v(k,j) < 0.0_wp ) THEN |
---|
| 740 | c_v(k,j) = 0.0_wp |
---|
[106] | 741 | ELSEIF ( c_v(k,j) > c_max ) THEN |
---|
| 742 | c_v(k,j) = c_max |
---|
| 743 | ENDIF |
---|
| 744 | ELSE |
---|
| 745 | c_v(k,j) = c_max |
---|
[75] | 746 | ENDIF |
---|
| 747 | |
---|
[106] | 748 | denom = w_m_l(k,j,0) - w_m_l(k,j,1) |
---|
[75] | 749 | |
---|
[1353] | 750 | IF ( denom /= 0.0_wp ) THEN |
---|
[996] | 751 | c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) ) |
---|
[1353] | 752 | IF ( c_w(k,j) < 0.0_wp ) THEN |
---|
| 753 | c_w(k,j) = 0.0_wp |
---|
[106] | 754 | ELSEIF ( c_w(k,j) > c_max ) THEN |
---|
| 755 | c_w(k,j) = c_max |
---|
| 756 | ENDIF |
---|
| 757 | ELSE |
---|
| 758 | c_w(k,j) = c_max |
---|
[75] | 759 | ENDIF |
---|
[106] | 760 | |
---|
[978] | 761 | c_u_m_l(k) = c_u_m_l(k) + c_u(k,j) |
---|
| 762 | c_v_m_l(k) = c_v_m_l(k) + c_v(k,j) |
---|
| 763 | c_w_m_l(k) = c_w_m_l(k) + c_w(k,j) |
---|
[106] | 764 | |
---|
[978] | 765 | ENDDO |
---|
| 766 | ENDDO |
---|
[75] | 767 | |
---|
[978] | 768 | #if defined( __parallel ) |
---|
| 769 | IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) |
---|
| 770 | CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, & |
---|
| 771 | MPI_SUM, comm1dy, ierr ) |
---|
| 772 | IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) |
---|
| 773 | CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, & |
---|
| 774 | MPI_SUM, comm1dy, ierr ) |
---|
| 775 | IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) |
---|
| 776 | CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, & |
---|
| 777 | MPI_SUM, comm1dy, ierr ) |
---|
| 778 | #else |
---|
| 779 | c_u_m = c_u_m_l |
---|
| 780 | c_v_m = c_v_m_l |
---|
| 781 | c_w_m = c_w_m_l |
---|
| 782 | #endif |
---|
| 783 | |
---|
| 784 | c_u_m = c_u_m / (ny+1) |
---|
| 785 | c_v_m = c_v_m / (ny+1) |
---|
| 786 | c_w_m = c_w_m / (ny+1) |
---|
| 787 | |
---|
[73] | 788 | ! |
---|
[978] | 789 | !-- Save old timelevels for the next timestep |
---|
| 790 | IF ( intermediate_timestep_count == 1 ) THEN |
---|
| 791 | u_m_l(:,:,:) = u(:,:,1:2) |
---|
| 792 | v_m_l(:,:,:) = v(:,:,0:1) |
---|
| 793 | w_m_l(:,:,:) = w(:,:,0:1) |
---|
| 794 | ENDIF |
---|
| 795 | |
---|
| 796 | ! |
---|
| 797 | !-- Calculate the new velocities |
---|
[996] | 798 | DO k = nzb+1, nzt+1 |
---|
[1113] | 799 | DO j = nysg, nyng |
---|
[978] | 800 | u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) * & |
---|
[106] | 801 | ( u(k,j,0) - u(k,j,1) ) * ddx |
---|
[75] | 802 | |
---|
[978] | 803 | v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) * & |
---|
[75] | 804 | ( v(k,j,-1) - v(k,j,0) ) * ddx |
---|
| 805 | |
---|
[978] | 806 | w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) * & |
---|
[75] | 807 | ( w(k,j,-1) - w(k,j,0) ) * ddx |
---|
[978] | 808 | ENDDO |
---|
[75] | 809 | ENDDO |
---|
| 810 | |
---|
| 811 | ! |
---|
[978] | 812 | !-- Bottom boundary at the outflow |
---|
| 813 | IF ( ibc_uv_b == 0 ) THEN |
---|
[1353] | 814 | u_p(nzb,:,0) = 0.0_wp |
---|
| 815 | v_p(nzb,:,-1) = 0.0_wp |
---|
[978] | 816 | ELSE |
---|
| 817 | u_p(nzb,:,0) = u_p(nzb+1,:,0) |
---|
| 818 | v_p(nzb,:,-1) = v_p(nzb+1,:,-1) |
---|
| 819 | ENDIF |
---|
[1353] | 820 | w_p(nzb,:,-1) = 0.0_wp |
---|
[1] | 821 | |
---|
[75] | 822 | ! |
---|
[978] | 823 | !-- Top boundary at the outflow |
---|
| 824 | IF ( ibc_uv_t == 0 ) THEN |
---|
[1764] | 825 | u_p(nzt+1,:,0) = u_init(nzt+1) |
---|
[978] | 826 | v_p(nzt+1,:,-1) = v_init(nzt+1) |
---|
| 827 | ELSE |
---|
[1764] | 828 | u_p(nzt+1,:,0) = u_p(nzt,:,0) |
---|
[978] | 829 | v_p(nzt+1,:,-1) = v_p(nzt,:,-1) |
---|
| 830 | ENDIF |
---|
[1353] | 831 | w_p(nzt:nzt+1,:,-1) = 0.0_wp |
---|
[978] | 832 | |
---|
[75] | 833 | ENDIF |
---|
[73] | 834 | |
---|
[75] | 835 | ENDIF |
---|
[73] | 836 | |
---|
[106] | 837 | IF ( outflow_r ) THEN |
---|
[73] | 838 | |
---|
[1159] | 839 | IF ( use_cmax ) THEN |
---|
| 840 | u_p(:,:,nx+1) = u(:,:,nx) |
---|
| 841 | v_p(:,:,nx+1) = v(:,:,nx) |
---|
| 842 | w_p(:,:,nx+1) = w(:,:,nx) |
---|
| 843 | ELSEIF ( .NOT. use_cmax ) THEN |
---|
[75] | 844 | |
---|
[978] | 845 | c_max = dx / dt_3d |
---|
[75] | 846 | |
---|
[1353] | 847 | c_u_m_l = 0.0_wp |
---|
| 848 | c_v_m_l = 0.0_wp |
---|
| 849 | c_w_m_l = 0.0_wp |
---|
[978] | 850 | |
---|
[1353] | 851 | c_u_m = 0.0_wp |
---|
| 852 | c_v_m = 0.0_wp |
---|
| 853 | c_w_m = 0.0_wp |
---|
[978] | 854 | |
---|
[1] | 855 | ! |
---|
[996] | 856 | !-- Calculate the phase speeds for u, v, and w, first local and then |
---|
| 857 | !-- average along the outflow boundary. |
---|
| 858 | DO k = nzb+1, nzt+1 |
---|
| 859 | DO j = nys, nyn |
---|
[73] | 860 | |
---|
[106] | 861 | denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1) |
---|
| 862 | |
---|
[1353] | 863 | IF ( denom /= 0.0_wp ) THEN |
---|
[996] | 864 | c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) ) |
---|
[1353] | 865 | IF ( c_u(k,j) < 0.0_wp ) THEN |
---|
| 866 | c_u(k,j) = 0.0_wp |
---|
[106] | 867 | ELSEIF ( c_u(k,j) > c_max ) THEN |
---|
| 868 | c_u(k,j) = c_max |
---|
| 869 | ENDIF |
---|
| 870 | ELSE |
---|
| 871 | c_u(k,j) = c_max |
---|
[73] | 872 | ENDIF |
---|
| 873 | |
---|
[106] | 874 | denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1) |
---|
[73] | 875 | |
---|
[1353] | 876 | IF ( denom /= 0.0_wp ) THEN |
---|
[996] | 877 | c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) ) |
---|
[1353] | 878 | IF ( c_v(k,j) < 0.0_wp ) THEN |
---|
| 879 | c_v(k,j) = 0.0_wp |
---|
[106] | 880 | ELSEIF ( c_v(k,j) > c_max ) THEN |
---|
| 881 | c_v(k,j) = c_max |
---|
| 882 | ENDIF |
---|
| 883 | ELSE |
---|
| 884 | c_v(k,j) = c_max |
---|
[73] | 885 | ENDIF |
---|
| 886 | |
---|
[106] | 887 | denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1) |
---|
[73] | 888 | |
---|
[1353] | 889 | IF ( denom /= 0.0_wp ) THEN |
---|
[996] | 890 | c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) ) |
---|
[1353] | 891 | IF ( c_w(k,j) < 0.0_wp ) THEN |
---|
| 892 | c_w(k,j) = 0.0_wp |
---|
[106] | 893 | ELSEIF ( c_w(k,j) > c_max ) THEN |
---|
| 894 | c_w(k,j) = c_max |
---|
| 895 | ENDIF |
---|
| 896 | ELSE |
---|
| 897 | c_w(k,j) = c_max |
---|
[73] | 898 | ENDIF |
---|
[106] | 899 | |
---|
[978] | 900 | c_u_m_l(k) = c_u_m_l(k) + c_u(k,j) |
---|
| 901 | c_v_m_l(k) = c_v_m_l(k) + c_v(k,j) |
---|
| 902 | c_w_m_l(k) = c_w_m_l(k) + c_w(k,j) |
---|
[106] | 903 | |
---|
[978] | 904 | ENDDO |
---|
| 905 | ENDDO |
---|
[73] | 906 | |
---|
[978] | 907 | #if defined( __parallel ) |
---|
| 908 | IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) |
---|
| 909 | CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, & |
---|
| 910 | MPI_SUM, comm1dy, ierr ) |
---|
| 911 | IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) |
---|
| 912 | CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, & |
---|
| 913 | MPI_SUM, comm1dy, ierr ) |
---|
| 914 | IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) |
---|
| 915 | CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, & |
---|
| 916 | MPI_SUM, comm1dy, ierr ) |
---|
| 917 | #else |
---|
| 918 | c_u_m = c_u_m_l |
---|
| 919 | c_v_m = c_v_m_l |
---|
| 920 | c_w_m = c_w_m_l |
---|
| 921 | #endif |
---|
| 922 | |
---|
| 923 | c_u_m = c_u_m / (ny+1) |
---|
| 924 | c_v_m = c_v_m / (ny+1) |
---|
| 925 | c_w_m = c_w_m / (ny+1) |
---|
| 926 | |
---|
[73] | 927 | ! |
---|
[978] | 928 | !-- Save old timelevels for the next timestep |
---|
| 929 | IF ( intermediate_timestep_count == 1 ) THEN |
---|
| 930 | u_m_r(:,:,:) = u(:,:,nx-1:nx) |
---|
| 931 | v_m_r(:,:,:) = v(:,:,nx-1:nx) |
---|
| 932 | w_m_r(:,:,:) = w(:,:,nx-1:nx) |
---|
| 933 | ENDIF |
---|
[73] | 934 | |
---|
[978] | 935 | ! |
---|
| 936 | !-- Calculate the new velocities |
---|
[996] | 937 | DO k = nzb+1, nzt+1 |
---|
[1113] | 938 | DO j = nysg, nyng |
---|
[978] | 939 | u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) * & |
---|
| 940 | ( u(k,j,nx+1) - u(k,j,nx) ) * ddx |
---|
[73] | 941 | |
---|
[978] | 942 | v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) * & |
---|
| 943 | ( v(k,j,nx+1) - v(k,j,nx) ) * ddx |
---|
[73] | 944 | |
---|
[978] | 945 | w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) * & |
---|
| 946 | ( w(k,j,nx+1) - w(k,j,nx) ) * ddx |
---|
| 947 | ENDDO |
---|
[73] | 948 | ENDDO |
---|
| 949 | |
---|
| 950 | ! |
---|
[978] | 951 | !-- Bottom boundary at the outflow |
---|
| 952 | IF ( ibc_uv_b == 0 ) THEN |
---|
[1353] | 953 | u_p(nzb,:,nx+1) = 0.0_wp |
---|
| 954 | v_p(nzb,:,nx+1) = 0.0_wp |
---|
[978] | 955 | ELSE |
---|
| 956 | u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1) |
---|
| 957 | v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1) |
---|
| 958 | ENDIF |
---|
[1353] | 959 | w_p(nzb,:,nx+1) = 0.0_wp |
---|
[73] | 960 | |
---|
| 961 | ! |
---|
[978] | 962 | !-- Top boundary at the outflow |
---|
| 963 | IF ( ibc_uv_t == 0 ) THEN |
---|
| 964 | u_p(nzt+1,:,nx+1) = u_init(nzt+1) |
---|
| 965 | v_p(nzt+1,:,nx+1) = v_init(nzt+1) |
---|
| 966 | ELSE |
---|
| 967 | u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1) |
---|
| 968 | v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1) |
---|
| 969 | ENDIF |
---|
[1742] | 970 | w_p(nzt:nzt+1,:,nx+1) = 0.0_wp |
---|
[978] | 971 | |
---|
[1] | 972 | ENDIF |
---|
| 973 | |
---|
| 974 | ENDIF |
---|
| 975 | |
---|
| 976 | END SUBROUTINE boundary_conds |
---|