[1] | 1 | SUBROUTINE boundary_conds( range ) |
---|
| 2 | |
---|
| 3 | !------------------------------------------------------------------------------! |
---|
| 4 | ! Actual revisions: |
---|
| 5 | ! ----------------- |
---|
[102] | 6 | ! Boundary conditions for temperature adjusted for coupled runs |
---|
[1] | 7 | ! |
---|
| 8 | ! Former revisions: |
---|
| 9 | ! ----------------- |
---|
[3] | 10 | ! $Id: boundary_conds.f90 102 2007-07-27 09:09:17Z raasch $ |
---|
[39] | 11 | ! |
---|
[98] | 12 | ! 95 2007-06-02 16:48:38Z raasch |
---|
| 13 | ! Boundary conditions for salinity added |
---|
| 14 | ! |
---|
[77] | 15 | ! 75 2007-03-22 09:54:05Z raasch |
---|
| 16 | ! The "main" part sets conditions for time level t+dt instead of level t, |
---|
| 17 | ! outflow boundary conditions changed from Neumann to radiation condition, |
---|
| 18 | ! uxrp, vynp eliminated, moisture renamed humidity |
---|
| 19 | ! |
---|
[39] | 20 | ! 19 2007-02-23 04:53:48Z raasch |
---|
| 21 | ! Boundary conditions for e(nzt), pt(nzt), and q(nzt) removed because these |
---|
| 22 | ! gridpoints are now calculated by the prognostic equation, |
---|
| 23 | ! Dirichlet and zero gradient condition for pt established at top boundary |
---|
| 24 | ! |
---|
[3] | 25 | ! RCS Log replace by Id keyword, revision history cleaned up |
---|
| 26 | ! |
---|
[1] | 27 | ! Revision 1.15 2006/02/23 09:54:55 raasch |
---|
| 28 | ! Surface boundary conditions in case of topography: nzb replaced by |
---|
| 29 | ! 2d-k-index-arrays (nzb_w_inner, etc.). Conditions for u and v remain |
---|
| 30 | ! unchanged (still using nzb) because a non-flat topography must use a |
---|
| 31 | ! Prandtl-layer, which don't requires explicit setting of the surface values. |
---|
| 32 | ! |
---|
| 33 | ! Revision 1.1 1997/09/12 06:21:34 raasch |
---|
| 34 | ! Initial revision |
---|
| 35 | ! |
---|
| 36 | ! |
---|
| 37 | ! Description: |
---|
| 38 | ! ------------ |
---|
| 39 | ! Boundary conditions for the prognostic quantities (range='main'). |
---|
| 40 | ! In case of non-cyclic lateral boundaries the conditions for velocities at |
---|
| 41 | ! the outflow are set after the pressure solver has been called (range= |
---|
| 42 | ! 'outflow_uvw'). |
---|
| 43 | ! One additional bottom boundary condition is applied for the TKE (=(u*)**2) |
---|
| 44 | ! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly |
---|
| 45 | ! handled in routine exchange_horiz. Pressure boundary conditions are |
---|
| 46 | ! explicitly set in routines pres, poisfft, poismg and sor. |
---|
| 47 | !------------------------------------------------------------------------------! |
---|
| 48 | |
---|
| 49 | USE arrays_3d |
---|
| 50 | USE control_parameters |
---|
| 51 | USE grid_variables |
---|
| 52 | USE indices |
---|
| 53 | USE pegrid |
---|
| 54 | |
---|
| 55 | IMPLICIT NONE |
---|
| 56 | |
---|
| 57 | CHARACTER (LEN=*) :: range |
---|
| 58 | |
---|
| 59 | INTEGER :: i, j, k |
---|
| 60 | |
---|
[73] | 61 | REAL :: c_max, c_u, c_v, c_w, denom |
---|
[1] | 62 | |
---|
[73] | 63 | |
---|
[1] | 64 | IF ( range == 'main') THEN |
---|
| 65 | ! |
---|
| 66 | !-- Bottom boundary |
---|
| 67 | IF ( ibc_uv_b == 0 ) THEN |
---|
[73] | 68 | ! |
---|
| 69 | !-- Satisfying the Dirichlet condition with an extra layer below the |
---|
| 70 | !-- surface where the u and v component change their sign |
---|
| 71 | u_p(nzb,:,:) = -u_p(nzb+1,:,:) |
---|
| 72 | v_p(nzb,:,:) = -v_p(nzb+1,:,:) |
---|
| 73 | ELSE |
---|
| 74 | u_p(nzb,:,:) = u_p(nzb+1,:,:) |
---|
| 75 | v_p(nzb,:,:) = v_p(nzb+1,:,:) |
---|
[1] | 76 | ENDIF |
---|
| 77 | DO i = nxl-1, nxr+1 |
---|
| 78 | DO j = nys-1, nyn+1 |
---|
[73] | 79 | w_p(nzb_w_inner(j,i),j,i) = 0.0 |
---|
[1] | 80 | ENDDO |
---|
| 81 | ENDDO |
---|
| 82 | |
---|
| 83 | ! |
---|
| 84 | !-- Top boundary |
---|
| 85 | IF ( ibc_uv_t == 0 ) THEN |
---|
[73] | 86 | u_p(nzt+1,:,:) = ug(nzt+1) |
---|
| 87 | v_p(nzt+1,:,:) = vg(nzt+1) |
---|
[1] | 88 | ELSE |
---|
[73] | 89 | u_p(nzt+1,:,:) = u_p(nzt,:,:) |
---|
| 90 | v_p(nzt+1,:,:) = v_p(nzt,:,:) |
---|
[1] | 91 | ENDIF |
---|
[73] | 92 | w_p(nzt:nzt+1,:,:) = 0.0 ! nzt is not a prognostic level (but cf. pres) |
---|
[1] | 93 | |
---|
| 94 | ! |
---|
[102] | 95 | !-- Temperature at bottom boundary. |
---|
| 96 | !-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by |
---|
| 97 | !-- the sea surface temperature of the coupled ocean model. |
---|
[1] | 98 | IF ( ibc_pt_b == 0 ) THEN |
---|
[73] | 99 | DO i = nxl-1, nxr+1 |
---|
| 100 | DO j = nys-1, nyn+1 |
---|
| 101 | pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i) |
---|
[1] | 102 | ENDDO |
---|
[73] | 103 | ENDDO |
---|
[102] | 104 | ELSEIF ( ibc_pt_b == 1 ) THEN |
---|
[1] | 105 | DO i = nxl-1, nxr+1 |
---|
| 106 | DO j = nys-1, nyn+1 |
---|
[73] | 107 | pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i) |
---|
[1] | 108 | ENDDO |
---|
| 109 | ENDDO |
---|
| 110 | ENDIF |
---|
| 111 | |
---|
| 112 | ! |
---|
| 113 | !-- Temperature at top boundary |
---|
[19] | 114 | IF ( ibc_pt_t == 0 ) THEN |
---|
[73] | 115 | pt_p(nzt+1,:,:) = pt(nzt+1,:,:) |
---|
[19] | 116 | ELSEIF ( ibc_pt_t == 1 ) THEN |
---|
[73] | 117 | pt_p(nzt+1,:,:) = pt_p(nzt,:,:) |
---|
[19] | 118 | ELSEIF ( ibc_pt_t == 2 ) THEN |
---|
[73] | 119 | pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1) |
---|
[1] | 120 | ENDIF |
---|
| 121 | |
---|
| 122 | ! |
---|
| 123 | !-- Boundary conditions for TKE |
---|
| 124 | !-- Generally Neumann conditions with de/dz=0 are assumed |
---|
| 125 | IF ( .NOT. constant_diffusion ) THEN |
---|
| 126 | DO i = nxl-1, nxr+1 |
---|
| 127 | DO j = nys-1, nyn+1 |
---|
[73] | 128 | e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i) |
---|
[1] | 129 | ENDDO |
---|
| 130 | ENDDO |
---|
[73] | 131 | e_p(nzt+1,:,:) = e_p(nzt,:,:) |
---|
[1] | 132 | ENDIF |
---|
| 133 | |
---|
| 134 | ! |
---|
[95] | 135 | !-- Boundary conditions for salinity |
---|
| 136 | IF ( ocean ) THEN |
---|
| 137 | ! |
---|
| 138 | !-- Bottom boundary: Neumann condition because salinity flux is always |
---|
| 139 | !-- given |
---|
| 140 | DO i = nxl-1, nxr+1 |
---|
| 141 | DO j = nys-1, nyn+1 |
---|
| 142 | sa_p(nzb_s_inner(j,i),j,i) = sa_p(nzb_s_inner(j,i)+1,j,i) |
---|
| 143 | ENDDO |
---|
| 144 | ENDDO |
---|
| 145 | |
---|
| 146 | ! |
---|
| 147 | !-- Top boundary: Dirichlet or Neumann |
---|
| 148 | IF ( ibc_sa_t == 0 ) THEN |
---|
| 149 | sa_p(nzt+1,:,:) = sa(nzt+1,:,:) |
---|
| 150 | ELSEIF ( ibc_sa_t == 1 ) THEN |
---|
| 151 | sa_p(nzt+1,:,:) = sa_p(nzt,:,:) |
---|
| 152 | ENDIF |
---|
| 153 | |
---|
| 154 | ENDIF |
---|
| 155 | |
---|
| 156 | ! |
---|
[1] | 157 | !-- Boundary conditions for total water content or scalar, |
---|
[95] | 158 | !-- bottom and top boundary (see also temperature) |
---|
[75] | 159 | IF ( humidity .OR. passive_scalar ) THEN |
---|
[1] | 160 | ! |
---|
[75] | 161 | !-- Surface conditions for constant_humidity_flux |
---|
[1] | 162 | IF ( ibc_q_b == 0 ) THEN |
---|
[73] | 163 | DO i = nxl-1, nxr+1 |
---|
| 164 | DO j = nys-1, nyn+1 |
---|
| 165 | q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i) |
---|
[1] | 166 | ENDDO |
---|
[73] | 167 | ENDDO |
---|
[1] | 168 | ELSE |
---|
| 169 | DO i = nxl-1, nxr+1 |
---|
| 170 | DO j = nys-1, nyn+1 |
---|
[73] | 171 | q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i) |
---|
[1] | 172 | ENDDO |
---|
| 173 | ENDDO |
---|
| 174 | ENDIF |
---|
| 175 | ! |
---|
| 176 | !-- Top boundary |
---|
[73] | 177 | q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1) |
---|
[1] | 178 | ENDIF |
---|
| 179 | |
---|
| 180 | ! |
---|
| 181 | !-- Lateral boundary conditions at the inflow. Quasi Neumann conditions |
---|
| 182 | !-- are needed for the wall normal velocity in order to ensure zero |
---|
| 183 | !-- divergence. Dirichlet conditions are used for all other quantities. |
---|
| 184 | IF ( inflow_s ) THEN |
---|
[73] | 185 | v_p(:,nys,:) = v_p(:,nys-1,:) |
---|
[1] | 186 | ELSEIF ( inflow_n ) THEN |
---|
[75] | 187 | v_p(:,nyn,:) = v_p(:,nyn+1,:) |
---|
[1] | 188 | ELSEIF ( inflow_l ) THEN |
---|
[73] | 189 | u_p(:,:,nxl) = u_p(:,:,nxl-1) |
---|
[1] | 190 | ELSEIF ( inflow_r ) THEN |
---|
[75] | 191 | u_p(:,:,nxr) = u_p(:,:,nxr+1) |
---|
[1] | 192 | ENDIF |
---|
| 193 | |
---|
| 194 | ! |
---|
| 195 | !-- Lateral boundary conditions for scalar quantities at the outflow |
---|
| 196 | IF ( outflow_s ) THEN |
---|
[73] | 197 | pt_p(:,nys-1,:) = pt_p(:,nys,:) |
---|
| 198 | IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:) |
---|
[75] | 199 | IF ( humidity .OR. passive_scalar ) q_p(:,nys-1,:) = q_p(:,nys,:) |
---|
[1] | 200 | ELSEIF ( outflow_n ) THEN |
---|
[73] | 201 | pt_p(:,nyn+1,:) = pt_p(:,nyn,:) |
---|
| 202 | IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:) |
---|
[75] | 203 | IF ( humidity .OR. passive_scalar ) q_p(:,nyn+1,:) = q_p(:,nyn,:) |
---|
[1] | 204 | ELSEIF ( outflow_l ) THEN |
---|
[73] | 205 | pt_p(:,:,nxl-1) = pt_p(:,:,nxl) |
---|
| 206 | IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl) |
---|
[75] | 207 | IF ( humidity .OR. passive_scalar ) q_p(:,:,nxl-1) = q_p(:,:,nxl) |
---|
[1] | 208 | ELSEIF ( outflow_r ) THEN |
---|
[73] | 209 | pt_p(:,:,nxr+1) = pt_p(:,:,nxr) |
---|
| 210 | IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr) |
---|
[75] | 211 | IF ( humidity .OR. passive_scalar ) q_p(:,:,nxr+1) = q_p(:,:,nxr) |
---|
[1] | 212 | ENDIF |
---|
| 213 | |
---|
| 214 | ENDIF |
---|
| 215 | |
---|
| 216 | ! |
---|
[75] | 217 | !-- Radiation boundary condition for the velocities at the respective outflow |
---|
| 218 | IF ( outflow_s .AND. & |
---|
| 219 | intermediate_timestep_count == intermediate_timestep_count_max ) & |
---|
| 220 | THEN |
---|
| 221 | |
---|
| 222 | c_max = dy / dt_3d |
---|
| 223 | |
---|
| 224 | DO i = nxl-1, nxr+1 |
---|
| 225 | DO k = nzb+1, nzt+1 |
---|
| 226 | |
---|
| 227 | ! |
---|
| 228 | !-- First calculate the phase speeds for u,v, and w |
---|
| 229 | denom = u_m_s(k,0,i) - u_m_s(k,1,i) |
---|
| 230 | |
---|
| 231 | IF ( denom /= 0.0 ) THEN |
---|
| 232 | c_u = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / denom |
---|
| 233 | IF ( c_u < 0.0 ) THEN |
---|
| 234 | c_u = 0.0 |
---|
| 235 | ELSEIF ( c_u > c_max ) THEN |
---|
| 236 | c_u = c_max |
---|
| 237 | ENDIF |
---|
| 238 | ELSE |
---|
| 239 | c_u = c_max |
---|
| 240 | ENDIF |
---|
| 241 | denom = v_m_s(k,0,i) - v_m_s(k,1,i) |
---|
| 242 | |
---|
| 243 | IF ( denom /= 0.0 ) THEN |
---|
| 244 | c_v = -c_max * ( v(k,0,i) - v_m_s(k,0,i) ) / denom |
---|
| 245 | IF ( c_v < 0.0 ) THEN |
---|
| 246 | c_v = 0.0 |
---|
| 247 | ELSEIF ( c_v > c_max ) THEN |
---|
| 248 | c_v = c_max |
---|
| 249 | ENDIF |
---|
| 250 | ELSE |
---|
| 251 | c_v = c_max |
---|
| 252 | ENDIF |
---|
| 253 | |
---|
| 254 | denom = w_m_s(k,0,i) - w_m_s(k,1,i) |
---|
| 255 | |
---|
| 256 | IF ( denom /= 0.0 ) THEN |
---|
| 257 | c_w = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / denom |
---|
| 258 | IF ( c_w < 0.0 ) THEN |
---|
| 259 | c_w = 0.0 |
---|
| 260 | ELSEIF ( c_w > c_max ) THEN |
---|
| 261 | c_w = c_max |
---|
| 262 | ENDIF |
---|
| 263 | ELSE |
---|
| 264 | c_w = c_max |
---|
| 265 | ENDIF |
---|
| 266 | |
---|
| 267 | ! |
---|
| 268 | !-- Calculate the new velocities |
---|
| 269 | u_p(k,-1,i) = u(k,-1,i) + dt_3d * c_u * & |
---|
| 270 | ( u(k,-1,i) - u(k,0,i) ) * ddy |
---|
| 271 | |
---|
| 272 | v_p(k,-1,i) = v(k,-1,i) + dt_3d * c_v * & |
---|
| 273 | ( v(k,-1,i) - v_m_s(k,0,i) ) * ddy |
---|
| 274 | |
---|
| 275 | w_p(k,-1,i) = w(k,-1,i) + dt_3d * c_w * & |
---|
| 276 | ( w(k,-1,i) - w(k,0,i) ) * ddy |
---|
| 277 | |
---|
| 278 | ! |
---|
| 279 | !-- Save old timelevels for the next timestep |
---|
| 280 | u_m_s(k,:,i) = u(k,-1:1,i) |
---|
| 281 | v_m_s(k,:,i) = v(k,-1:1,i) |
---|
| 282 | w_m_s(k,:,i) = w(k,-1:1,i) |
---|
| 283 | |
---|
| 284 | ENDDO |
---|
| 285 | ENDDO |
---|
| 286 | |
---|
| 287 | ! |
---|
| 288 | !-- Bottom boundary at the outflow |
---|
| 289 | IF ( ibc_uv_b == 0 ) THEN |
---|
| 290 | u_p(nzb,-1,:) = -u_p(nzb+1,-1,:) |
---|
| 291 | v_p(nzb,-1,:) = -v_p(nzb+1,-1,:) |
---|
| 292 | ELSE |
---|
| 293 | u_p(nzb,-1,:) = u_p(nzb+1,-1,:) |
---|
| 294 | v_p(nzb,-1,:) = v_p(nzb+1,-1,:) |
---|
[73] | 295 | ENDIF |
---|
[75] | 296 | w_p(nzb,ny+1,:) = 0.0 |
---|
[73] | 297 | |
---|
[75] | 298 | ! |
---|
| 299 | !-- Top boundary at the outflow |
---|
| 300 | IF ( ibc_uv_t == 0 ) THEN |
---|
| 301 | u_p(nzt+1,-1,:) = ug(nzt+1) |
---|
| 302 | v_p(nzt+1,-1,:) = vg(nzt+1) |
---|
| 303 | ELSE |
---|
| 304 | u_p(nzt+1,-1,:) = u(nzt,-1,:) |
---|
| 305 | v_p(nzt+1,-1,:) = v(nzt,-1,:) |
---|
| 306 | ENDIF |
---|
| 307 | w_p(nzt:nzt+1,-1,:) = 0.0 |
---|
[73] | 308 | |
---|
[75] | 309 | ENDIF |
---|
[73] | 310 | |
---|
[75] | 311 | IF ( outflow_n .AND. & |
---|
| 312 | intermediate_timestep_count == intermediate_timestep_count_max ) & |
---|
| 313 | THEN |
---|
[73] | 314 | |
---|
[75] | 315 | c_max = dy / dt_3d |
---|
| 316 | |
---|
| 317 | DO i = nxl-1, nxr+1 |
---|
| 318 | DO k = nzb+1, nzt+1 |
---|
| 319 | |
---|
[1] | 320 | ! |
---|
[75] | 321 | !-- First calculate the phase speeds for u,v, and w |
---|
| 322 | denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i) |
---|
[73] | 323 | |
---|
[75] | 324 | IF ( denom /= 0.0 ) THEN |
---|
| 325 | c_u = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / denom |
---|
| 326 | IF ( c_u < 0.0 ) THEN |
---|
| 327 | c_u = 0.0 |
---|
| 328 | ELSEIF ( c_u > c_max ) THEN |
---|
[73] | 329 | c_u = c_max |
---|
| 330 | ENDIF |
---|
[75] | 331 | ELSE |
---|
| 332 | c_u = c_max |
---|
| 333 | ENDIF |
---|
[73] | 334 | |
---|
[75] | 335 | denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i) |
---|
[73] | 336 | |
---|
[75] | 337 | IF ( denom /= 0.0 ) THEN |
---|
| 338 | c_v = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / denom |
---|
| 339 | IF ( c_v < 0.0 ) THEN |
---|
| 340 | c_v = 0.0 |
---|
| 341 | ELSEIF ( c_v > c_max ) THEN |
---|
[73] | 342 | c_v = c_max |
---|
| 343 | ENDIF |
---|
[75] | 344 | ELSE |
---|
| 345 | c_v = c_max |
---|
| 346 | ENDIF |
---|
[73] | 347 | |
---|
[75] | 348 | denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i) |
---|
[73] | 349 | |
---|
[75] | 350 | IF ( denom /= 0.0 ) THEN |
---|
| 351 | c_w = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / denom |
---|
| 352 | IF ( c_w < 0.0 ) THEN |
---|
| 353 | c_w = 0.0 |
---|
| 354 | ELSEIF ( c_w > c_max ) THEN |
---|
[73] | 355 | c_w = c_max |
---|
| 356 | ENDIF |
---|
[75] | 357 | ELSE |
---|
| 358 | c_w = c_max |
---|
| 359 | ENDIF |
---|
[73] | 360 | |
---|
| 361 | ! |
---|
[75] | 362 | !-- Calculate the new velocities |
---|
| 363 | u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * c_u * & |
---|
| 364 | ( u(k,ny+1,i) - u(k,ny,i) ) * ddy |
---|
[73] | 365 | |
---|
[75] | 366 | v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * c_v * & |
---|
| 367 | ( v(k,ny+1,i) - v(k,ny,i) ) * ddy |
---|
[73] | 368 | |
---|
[75] | 369 | w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * c_w * & |
---|
| 370 | ( w(k,ny+1,i) - w(k,ny,i) ) * ddy |
---|
[73] | 371 | |
---|
| 372 | ! |
---|
[75] | 373 | !-- Swap timelevels for the next timestep |
---|
| 374 | u_m_n(k,:,i) = u(k,ny-1:ny+1,i) |
---|
| 375 | v_m_n(k,:,i) = v(k,ny-1:ny+1,i) |
---|
| 376 | w_m_n(k,:,i) = w(k,ny-1:ny+1,i) |
---|
[73] | 377 | |
---|
[1] | 378 | ENDDO |
---|
[75] | 379 | ENDDO |
---|
[1] | 380 | |
---|
| 381 | ! |
---|
[75] | 382 | !-- Bottom boundary at the outflow |
---|
| 383 | IF ( ibc_uv_b == 0 ) THEN |
---|
| 384 | u_p(nzb,ny+1,:) = -u_p(nzb+1,ny+1,:) |
---|
| 385 | v_p(nzb,ny+1,:) = -v_p(nzb+1,ny+1,:) |
---|
| 386 | ELSE |
---|
| 387 | u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:) |
---|
| 388 | v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:) |
---|
| 389 | ENDIF |
---|
| 390 | w_p(nzb,ny+1,:) = 0.0 |
---|
[73] | 391 | |
---|
| 392 | ! |
---|
[75] | 393 | !-- Top boundary at the outflow |
---|
| 394 | IF ( ibc_uv_t == 0 ) THEN |
---|
| 395 | u_p(nzt+1,ny+1,:) = ug(nzt+1) |
---|
| 396 | v_p(nzt+1,ny+1,:) = vg(nzt+1) |
---|
| 397 | ELSE |
---|
| 398 | u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:) |
---|
| 399 | v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:) |
---|
[1] | 400 | ENDIF |
---|
[75] | 401 | w_p(nzt:nzt+1,ny+1,:) = 0.0 |
---|
[1] | 402 | |
---|
[75] | 403 | ENDIF |
---|
| 404 | |
---|
| 405 | IF ( outflow_l .AND. & |
---|
| 406 | intermediate_timestep_count == intermediate_timestep_count_max ) & |
---|
| 407 | THEN |
---|
| 408 | |
---|
| 409 | c_max = dx / dt_3d |
---|
| 410 | |
---|
| 411 | DO j = nys-1, nyn+1 |
---|
| 412 | DO k = nzb+1, nzt+1 |
---|
| 413 | |
---|
[1] | 414 | ! |
---|
[75] | 415 | !-- First calculate the phase speeds for u,v, and w |
---|
| 416 | denom = u_m_l(k,j,0) - u_m_l(k,j,1) |
---|
| 417 | |
---|
| 418 | IF ( denom /= 0.0 ) THEN |
---|
| 419 | c_u = -c_max * ( u(k,j,0) - u_m_r(k,j,0) ) / denom |
---|
| 420 | IF ( c_u > 0.0 ) THEN |
---|
| 421 | c_u = 0.0 |
---|
| 422 | ELSEIF ( c_u < -c_max ) THEN |
---|
| 423 | c_u = -c_max |
---|
| 424 | ENDIF |
---|
| 425 | ELSE |
---|
| 426 | c_u = -c_max |
---|
| 427 | ENDIF |
---|
| 428 | |
---|
| 429 | denom = v_m_l(k,j,0) - v_m_l(k,j,1) |
---|
| 430 | |
---|
| 431 | IF ( denom /= 0.0 ) THEN |
---|
| 432 | c_v = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / denom |
---|
| 433 | IF ( c_v < 0.0 ) THEN |
---|
| 434 | c_v = 0.0 |
---|
| 435 | ELSEIF ( c_v > c_max ) THEN |
---|
| 436 | c_v = c_max |
---|
| 437 | ENDIF |
---|
| 438 | ELSE |
---|
| 439 | c_v = c_max |
---|
| 440 | ENDIF |
---|
| 441 | |
---|
| 442 | denom = w_m_l(k,j,0) - w_m_l(k,j,1) |
---|
| 443 | |
---|
| 444 | IF ( denom /= 0.0 ) THEN |
---|
| 445 | c_w = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / denom |
---|
| 446 | IF ( c_w < 0.0 ) THEN |
---|
| 447 | c_w = 0.0 |
---|
| 448 | ELSEIF ( c_w > c_max ) THEN |
---|
| 449 | c_w = c_max |
---|
| 450 | ENDIF |
---|
| 451 | ELSE |
---|
| 452 | c_w = c_max |
---|
| 453 | ENDIF |
---|
| 454 | |
---|
[73] | 455 | ! |
---|
[75] | 456 | !-- Calculate the new velocities |
---|
| 457 | u_p(k,j,-1) = u(k,j,-1) + dt_3d * c_u * & |
---|
| 458 | ( u(k,j,-1) - u(k,j,0) ) * ddx |
---|
| 459 | |
---|
| 460 | v_p(k,j,-1) = v(k,j,-1) + dt_3d * c_v * & |
---|
| 461 | ( v(k,j,-1) - v(k,j,0) ) * ddx |
---|
| 462 | |
---|
| 463 | w_p(k,j,-1) = w(k,j,-1) + dt_3d * c_w * & |
---|
| 464 | ( w(k,j,-1) - w(k,j,0) ) * ddx |
---|
| 465 | |
---|
| 466 | ! |
---|
| 467 | !-- Swap timelevels for the next timestep |
---|
| 468 | u_m_l(k,j,:) = u(k,j,-1:1) |
---|
| 469 | v_m_l(k,j,:) = v(k,j,-1:1) |
---|
| 470 | w_m_l(k,j,:) = w(k,j,-1:1) |
---|
| 471 | |
---|
| 472 | ENDDO |
---|
| 473 | ENDDO |
---|
| 474 | |
---|
| 475 | ! |
---|
| 476 | !-- Bottom boundary at the outflow |
---|
| 477 | IF ( ibc_uv_b == 0 ) THEN |
---|
| 478 | u_p(nzb,:,-1) = -u_p(nzb+1,:,-1) |
---|
| 479 | v_p(nzb,:,-1) = -v_p(nzb+1,:,-1) |
---|
| 480 | ELSE |
---|
| 481 | u_p(nzb,:,-1) = u_p(nzb+1,:,-1) |
---|
| 482 | v_p(nzb,:,-1) = v_p(nzb+1,:,-1) |
---|
[1] | 483 | ENDIF |
---|
[75] | 484 | w_p(nzb,:,-1) = 0.0 |
---|
[1] | 485 | |
---|
[75] | 486 | ! |
---|
| 487 | !-- Top boundary at the outflow |
---|
| 488 | IF ( ibc_uv_t == 0 ) THEN |
---|
| 489 | u_p(nzt+1,:,-1) = ug(nzt+1) |
---|
| 490 | v_p(nzt+1,:,-1) = vg(nzt+1) |
---|
| 491 | ELSE |
---|
| 492 | u_p(nzt+1,:,-1) = u_p(nzt,:,-1) |
---|
| 493 | v_p(nzt+1,:,-1) = v_p(nzt,:,-1) |
---|
| 494 | ENDIF |
---|
| 495 | w_p(nzt:nzt+1,:,-1) = 0.0 |
---|
[73] | 496 | |
---|
[75] | 497 | ENDIF |
---|
[73] | 498 | |
---|
[75] | 499 | IF ( outflow_r .AND. & |
---|
| 500 | intermediate_timestep_count == intermediate_timestep_count_max ) & |
---|
| 501 | THEN |
---|
[73] | 502 | |
---|
[75] | 503 | c_max = dx / dt_3d |
---|
| 504 | |
---|
| 505 | DO j = nys-1, nyn+1 |
---|
| 506 | DO k = nzb+1, nzt+1 |
---|
| 507 | |
---|
[1] | 508 | ! |
---|
[75] | 509 | !-- First calculate the phase speeds for u,v, and w |
---|
| 510 | denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1) |
---|
[73] | 511 | |
---|
[75] | 512 | IF ( denom /= 0.0 ) THEN |
---|
| 513 | c_u = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / denom |
---|
| 514 | IF ( c_u < 0.0 ) THEN |
---|
| 515 | c_u = 0.0 |
---|
| 516 | ELSEIF ( c_u > c_max ) THEN |
---|
[73] | 517 | c_u = c_max |
---|
| 518 | ENDIF |
---|
[75] | 519 | ELSE |
---|
| 520 | c_u = c_max |
---|
| 521 | ENDIF |
---|
[73] | 522 | |
---|
[75] | 523 | denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1) |
---|
[73] | 524 | |
---|
[75] | 525 | IF ( denom /= 0.0 ) THEN |
---|
| 526 | c_v = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / denom |
---|
| 527 | IF ( c_v < 0.0 ) THEN |
---|
| 528 | c_v = 0.0 |
---|
| 529 | ELSEIF ( c_v > c_max ) THEN |
---|
[73] | 530 | c_v = c_max |
---|
| 531 | ENDIF |
---|
[75] | 532 | ELSE |
---|
| 533 | c_v = c_max |
---|
| 534 | ENDIF |
---|
[73] | 535 | |
---|
[75] | 536 | denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1) |
---|
[73] | 537 | |
---|
[75] | 538 | IF ( denom /= 0.0 ) THEN |
---|
| 539 | c_w = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / denom |
---|
| 540 | IF ( c_w < 0.0 ) THEN |
---|
| 541 | c_w = 0.0 |
---|
| 542 | ELSEIF ( c_w > c_max ) THEN |
---|
[73] | 543 | c_w = c_max |
---|
| 544 | ENDIF |
---|
[75] | 545 | ELSE |
---|
| 546 | c_w = c_max |
---|
| 547 | ENDIF |
---|
[73] | 548 | |
---|
| 549 | ! |
---|
[75] | 550 | !-- Calculate the new velocities |
---|
| 551 | u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * c_u * & |
---|
| 552 | ( u(k,j,nx+1) - u(k,j,nx) ) * ddx |
---|
[73] | 553 | |
---|
[75] | 554 | v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * c_v * & |
---|
| 555 | ( v(k,j,nx+1) - v(k,j,nx) ) * ddx |
---|
[73] | 556 | |
---|
[75] | 557 | w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * c_w * & |
---|
| 558 | ( w(k,j,nx+1) - w(k,j,nx) ) * ddx |
---|
[73] | 559 | |
---|
| 560 | ! |
---|
[75] | 561 | !-- Swap timelevels for the next timestep |
---|
| 562 | u_m_r(k,j,:) = u(k,j,nx-1:nx+1) |
---|
| 563 | v_m_r(k,j,:) = v(k,j,nx-1:nx+1) |
---|
| 564 | w_m_r(k,j,:) = w(k,j,nx-1:nx+1) |
---|
[73] | 565 | |
---|
| 566 | ENDDO |
---|
[75] | 567 | ENDDO |
---|
[73] | 568 | |
---|
| 569 | ! |
---|
[75] | 570 | !-- Bottom boundary at the outflow |
---|
| 571 | IF ( ibc_uv_b == 0 ) THEN |
---|
| 572 | u_p(nzb,:,nx+1) = -u_p(nzb+1,:,nx+1) |
---|
| 573 | v_p(nzb,:,nx+1) = -v_p(nzb+1,:,nx+1) |
---|
| 574 | ELSE |
---|
| 575 | u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1) |
---|
| 576 | v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1) |
---|
| 577 | ENDIF |
---|
| 578 | w_p(nzb,:,nx+1) = 0.0 |
---|
[73] | 579 | |
---|
| 580 | ! |
---|
[75] | 581 | !-- Top boundary at the outflow |
---|
| 582 | IF ( ibc_uv_t == 0 ) THEN |
---|
| 583 | u_p(nzt+1,:,nx+1) = ug(nzt+1) |
---|
| 584 | v_p(nzt+1,:,nx+1) = vg(nzt+1) |
---|
| 585 | ELSE |
---|
| 586 | u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1) |
---|
| 587 | v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1) |
---|
[1] | 588 | ENDIF |
---|
[75] | 589 | w(nzt:nzt+1,:,nx+1) = 0.0 |
---|
[1] | 590 | |
---|
| 591 | ENDIF |
---|
| 592 | |
---|
| 593 | |
---|
| 594 | END SUBROUTINE boundary_conds |
---|