| 1 | SUBROUTINE boundary_conds( range ) |
|---|
| 2 | |
|---|
| 3 | !------------------------------------------------------------------------------! |
|---|
| 4 | ! Actual revisions: |
|---|
| 5 | ! ----------------- |
|---|
| 6 | ! |
|---|
| 7 | ! |
|---|
| 8 | ! Former revisions: |
|---|
| 9 | ! ----------------- |
|---|
| 10 | ! $Id: boundary_conds.f90 4 2007-02-13 11:33:16Z raasch $ |
|---|
| 11 | ! RCS Log replace by Id keyword, revision history cleaned up |
|---|
| 12 | ! |
|---|
| 13 | ! Revision 1.15 2006/02/23 09:54:55 raasch |
|---|
| 14 | ! Surface boundary conditions in case of topography: nzb replaced by |
|---|
| 15 | ! 2d-k-index-arrays (nzb_w_inner, etc.). Conditions for u and v remain |
|---|
| 16 | ! unchanged (still using nzb) because a non-flat topography must use a |
|---|
| 17 | ! Prandtl-layer, which don't requires explicit setting of the surface values. |
|---|
| 18 | ! |
|---|
| 19 | ! Revision 1.1 1997/09/12 06:21:34 raasch |
|---|
| 20 | ! Initial revision |
|---|
| 21 | ! |
|---|
| 22 | ! |
|---|
| 23 | ! Description: |
|---|
| 24 | ! ------------ |
|---|
| 25 | ! Boundary conditions for the prognostic quantities (range='main'). |
|---|
| 26 | ! In case of non-cyclic lateral boundaries the conditions for velocities at |
|---|
| 27 | ! the outflow are set after the pressure solver has been called (range= |
|---|
| 28 | ! 'outflow_uvw'). |
|---|
| 29 | ! One additional bottom boundary condition is applied for the TKE (=(u*)**2) |
|---|
| 30 | ! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly |
|---|
| 31 | ! handled in routine exchange_horiz. Pressure boundary conditions are |
|---|
| 32 | ! explicitly set in routines pres, poisfft, poismg and sor. |
|---|
| 33 | !------------------------------------------------------------------------------! |
|---|
| 34 | |
|---|
| 35 | USE arrays_3d |
|---|
| 36 | USE control_parameters |
|---|
| 37 | USE grid_variables |
|---|
| 38 | USE indices |
|---|
| 39 | USE pegrid |
|---|
| 40 | |
|---|
| 41 | IMPLICIT NONE |
|---|
| 42 | |
|---|
| 43 | CHARACTER (LEN=*) :: range |
|---|
| 44 | |
|---|
| 45 | INTEGER :: i, j, k |
|---|
| 46 | |
|---|
| 47 | |
|---|
| 48 | IF ( range == 'main') THEN |
|---|
| 49 | ! |
|---|
| 50 | !-- Bottom boundary |
|---|
| 51 | IF ( ibc_uv_b == 0 ) THEN |
|---|
| 52 | u(nzb,:,:) = -u(nzb+1,:,:) ! satisfying the Dirichlet condition with |
|---|
| 53 | v(nzb,:,:) = -v(nzb+1,:,:) ! an extra layer below the surface where |
|---|
| 54 | ELSE ! the u and v component change their sign |
|---|
| 55 | u(nzb,:,:) = u(nzb+1,:,:) |
|---|
| 56 | v(nzb,:,:) = v(nzb+1,:,:) |
|---|
| 57 | ENDIF |
|---|
| 58 | DO i = nxl-1, nxr+1 |
|---|
| 59 | DO j = nys-1, nyn+1 |
|---|
| 60 | w(nzb_w_inner(j,i),j,i) = 0.0 |
|---|
| 61 | ENDDO |
|---|
| 62 | ENDDO |
|---|
| 63 | |
|---|
| 64 | ! |
|---|
| 65 | !-- Top boundary |
|---|
| 66 | IF ( ibc_uv_t == 0 ) THEN |
|---|
| 67 | u(nzt+1,:,:) = ug(nzt+1) |
|---|
| 68 | v(nzt+1,:,:) = vg(nzt+1) |
|---|
| 69 | ELSE |
|---|
| 70 | u(nzt+1,:,:) = u(nzt,:,:) |
|---|
| 71 | v(nzt+1,:,:) = v(nzt,:,:) |
|---|
| 72 | ENDIF |
|---|
| 73 | w(nzt:nzt+1,:,:) = 0.0 ! nzt is not a prognostic level (but cf. pres) |
|---|
| 74 | |
|---|
| 75 | ! |
|---|
| 76 | !-- Temperature at bottom boundary |
|---|
| 77 | IF ( ibc_pt_b == 0 ) THEN |
|---|
| 78 | IF ( timestep_scheme(1:5) /= 'runge' ) THEN |
|---|
| 79 | DO i = nxl-1, nxr+1 |
|---|
| 80 | DO j = nys-1, nyn+1 |
|---|
| 81 | pt(nzb_s_inner(j,i),j,i) = pt_m(nzb_s_inner(j,i),j,i) |
|---|
| 82 | ENDDO |
|---|
| 83 | ENDDO |
|---|
| 84 | ELSE |
|---|
| 85 | DO i = nxl-1, nxr+1 |
|---|
| 86 | DO j = nys-1, nyn+1 |
|---|
| 87 | pt(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i),j,i) |
|---|
| 88 | ! pt_m is not used for Runge-Kutta |
|---|
| 89 | ENDDO |
|---|
| 90 | ENDDO |
|---|
| 91 | ENDIF |
|---|
| 92 | ELSE |
|---|
| 93 | DO i = nxl-1, nxr+1 |
|---|
| 94 | DO j = nys-1, nyn+1 |
|---|
| 95 | pt(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i)+1,j,i) |
|---|
| 96 | ENDDO |
|---|
| 97 | ENDDO |
|---|
| 98 | ENDIF |
|---|
| 99 | |
|---|
| 100 | ! |
|---|
| 101 | !-- Temperature at top boundary |
|---|
| 102 | IF ( ibc_pt_t == 1 ) THEN |
|---|
| 103 | pt(nzt,:,:) = pt(nzt-1,:,:) + bc_pt_t_val * dzu(nzt) |
|---|
| 104 | pt(nzt+1,:,:) = pt(nzt,:,:) + bc_pt_t_val * dzu(nzt+1) |
|---|
| 105 | ENDIF |
|---|
| 106 | |
|---|
| 107 | ! |
|---|
| 108 | !-- Boundary conditions for TKE |
|---|
| 109 | !-- Generally Neumann conditions with de/dz=0 are assumed |
|---|
| 110 | IF ( .NOT. constant_diffusion ) THEN |
|---|
| 111 | DO i = nxl-1, nxr+1 |
|---|
| 112 | DO j = nys-1, nyn+1 |
|---|
| 113 | e(nzb_s_inner(j,i),j,i) = e(nzb_s_inner(j,i)+1,j,i) |
|---|
| 114 | ENDDO |
|---|
| 115 | ENDDO |
|---|
| 116 | e(nzt,:,:) = e(nzt-1,:,:) |
|---|
| 117 | e(nzt+1,:,:) = e(nzt,:,:) |
|---|
| 118 | ENDIF |
|---|
| 119 | |
|---|
| 120 | ! |
|---|
| 121 | !-- Boundary conditions for total water content or scalar, |
|---|
| 122 | !-- bottom and surface boundary (see also temperature) |
|---|
| 123 | IF ( moisture .OR. passive_scalar ) THEN |
|---|
| 124 | ! |
|---|
| 125 | !-- Surface conditions for constant_moisture_flux |
|---|
| 126 | IF ( ibc_q_b == 0 ) THEN |
|---|
| 127 | IF ( timestep_scheme(1:5) /= 'runge' ) THEN |
|---|
| 128 | DO i = nxl-1, nxr+1 |
|---|
| 129 | DO j = nys-1, nyn+1 |
|---|
| 130 | q(nzb_s_inner(j,i),j,i) = q_m(nzb_s_inner(j,i),j,i) |
|---|
| 131 | ENDDO |
|---|
| 132 | ENDDO |
|---|
| 133 | ELSE |
|---|
| 134 | DO i = nxl-1, nxr+1 |
|---|
| 135 | DO j = nys-1, nyn+1 |
|---|
| 136 | q(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i),j,i) |
|---|
| 137 | ENDDO ! q_m is not used for Runge-Kutta |
|---|
| 138 | ENDDO |
|---|
| 139 | ENDIF |
|---|
| 140 | ELSE |
|---|
| 141 | DO i = nxl-1, nxr+1 |
|---|
| 142 | DO j = nys-1, nyn+1 |
|---|
| 143 | q(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i)+1,j,i) |
|---|
| 144 | ENDDO |
|---|
| 145 | ENDDO |
|---|
| 146 | ENDIF |
|---|
| 147 | ! |
|---|
| 148 | !-- Top boundary |
|---|
| 149 | q(nzt,:,:) = q(nzt-1,:,:) + bc_q_t_val * dzu(nzt) |
|---|
| 150 | q(nzt+1,:,:) = q(nzt,:,:) + bc_q_t_val * dzu(nzt+1) |
|---|
| 151 | ENDIF |
|---|
| 152 | |
|---|
| 153 | ! |
|---|
| 154 | !-- Lateral boundary conditions at the inflow. Quasi Neumann conditions |
|---|
| 155 | !-- are needed for the wall normal velocity in order to ensure zero |
|---|
| 156 | !-- divergence. Dirichlet conditions are used for all other quantities. |
|---|
| 157 | IF ( inflow_s ) THEN |
|---|
| 158 | v(:,nys,:) = v(:,nys-1,:) |
|---|
| 159 | ELSEIF ( inflow_n ) THEN |
|---|
| 160 | v(:,nyn+vynp,:) = v(:,nyn+vynp+1,:) |
|---|
| 161 | ELSEIF ( inflow_l ) THEN |
|---|
| 162 | u(:,:,nxl) = u(:,:,nxl-1) |
|---|
| 163 | ELSEIF ( inflow_r ) THEN |
|---|
| 164 | u(:,:,nxr+uxrp) = u(:,:,nxr+uxrp+1) |
|---|
| 165 | ENDIF |
|---|
| 166 | |
|---|
| 167 | ! |
|---|
| 168 | !-- Lateral boundary conditions for scalar quantities at the outflow |
|---|
| 169 | IF ( outflow_s ) THEN |
|---|
| 170 | pt(:,nys-1,:) = pt(:,nys,:) |
|---|
| 171 | IF ( .NOT. constant_diffusion ) e(:,nys-1,:) = e(:,nys,:) |
|---|
| 172 | IF ( moisture .OR. passive_scalar ) q(:,nys-1,:) = q(:,nys,:) |
|---|
| 173 | ELSEIF ( outflow_n ) THEN |
|---|
| 174 | pt(:,nyn+1,:) = pt(:,nyn,:) |
|---|
| 175 | IF ( .NOT. constant_diffusion ) e(:,nyn+1,:) = e(:,nyn,:) |
|---|
| 176 | IF ( moisture .OR. passive_scalar ) q(:,nyn+1,:) = q(:,nyn,:) |
|---|
| 177 | ELSEIF ( outflow_l ) THEN |
|---|
| 178 | pt(:,:,nxl-1) = pt(:,:,nxl) |
|---|
| 179 | IF ( .NOT. constant_diffusion ) e(:,:,nxl-1) = e(:,:,nxl) |
|---|
| 180 | IF ( moisture .OR. passive_scalar ) q(:,:,nxl-1) = q(:,:,nxl) |
|---|
| 181 | ELSEIF ( outflow_r ) THEN |
|---|
| 182 | pt(:,:,nxr+1) = pt(:,:,nxr) |
|---|
| 183 | IF ( .NOT. constant_diffusion ) e(:,:,nxr+1) = e(:,:,nxr) |
|---|
| 184 | IF ( moisture .OR. passive_scalar ) q(:,:,nxr+1) = q(:,:,nxr) |
|---|
| 185 | ENDIF |
|---|
| 186 | |
|---|
| 187 | ENDIF |
|---|
| 188 | |
|---|
| 189 | IF ( range == 'outflow_uvw' ) THEN |
|---|
| 190 | ! |
|---|
| 191 | !-- Horizontal boundary conditions for the velocities at the outflow. |
|---|
| 192 | !-- A Neumann condition is used for the wall normal velocity. The vertical |
|---|
| 193 | !-- velocity is assumed as zero and a horizontal average along the wall is |
|---|
| 194 | !-- used for the wall parallel horizontal velocity. The combination of all |
|---|
| 195 | !-- three conditions ensures that the velocity field is free of divergence |
|---|
| 196 | !-- at the outflow (uvmean_outflow_l is calculated in pres). |
|---|
| 197 | IF ( outflow_s ) THEN |
|---|
| 198 | v(:,nys-1,:) = v(:,nys,:) |
|---|
| 199 | w(:,nys-1,:) = 0.0 |
|---|
| 200 | ! |
|---|
| 201 | !-- Compute the mean horizontal wind parallel to and within the outflow |
|---|
| 202 | !-- wall and use this as boundary condition for u |
|---|
| 203 | #if defined( __parallel ) |
|---|
| 204 | CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, & |
|---|
| 205 | MPI_REAL, MPI_SUM, comm1dx, ierr ) |
|---|
| 206 | uvmean_outflow = uvmean_outflow / ( nx + 1.0 ) |
|---|
| 207 | #else |
|---|
| 208 | uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 ) |
|---|
| 209 | #endif |
|---|
| 210 | DO k = nzb, nzt+1 |
|---|
| 211 | u(k,nys-1,:) = uvmean_outflow(k) |
|---|
| 212 | ENDDO |
|---|
| 213 | ENDIF |
|---|
| 214 | |
|---|
| 215 | IF ( outflow_n ) THEN |
|---|
| 216 | v(:,nyn+vynp+1,:) = v(:,nyn+vynp,:) |
|---|
| 217 | w(:,nyn+1,:) = 0.0 |
|---|
| 218 | ! |
|---|
| 219 | !-- Compute the mean horizontal wind parallel to and within the outflow |
|---|
| 220 | !-- wall and use this as boundary condition for u |
|---|
| 221 | #if defined( __parallel ) |
|---|
| 222 | CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, & |
|---|
| 223 | MPI_REAL, MPI_SUM, comm1dx, ierr ) |
|---|
| 224 | uvmean_outflow = uvmean_outflow / ( nx + 1.0 ) |
|---|
| 225 | #else |
|---|
| 226 | uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 ) |
|---|
| 227 | #endif |
|---|
| 228 | DO k = nzb, nzt+1 |
|---|
| 229 | u(k,nyn+1,:) = uvmean_outflow(k) |
|---|
| 230 | ENDDO |
|---|
| 231 | ENDIF |
|---|
| 232 | |
|---|
| 233 | IF ( outflow_l ) THEN |
|---|
| 234 | u(:,:,nxl-1) = u(:,:,nxl) |
|---|
| 235 | w(:,:,nxl-1) = 0.0 |
|---|
| 236 | ! |
|---|
| 237 | !-- Compute the mean horizontal wind parallel to and within the outflow |
|---|
| 238 | !-- wall and use this as boundary condition for v |
|---|
| 239 | #if defined( __parallel ) |
|---|
| 240 | CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, & |
|---|
| 241 | MPI_REAL, MPI_SUM, comm1dy, ierr ) |
|---|
| 242 | uvmean_outflow = uvmean_outflow / ( ny + 1.0 ) |
|---|
| 243 | #else |
|---|
| 244 | uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 ) |
|---|
| 245 | #endif |
|---|
| 246 | DO k = nzb, nzt+1 |
|---|
| 247 | v(k,:,nxl-1) = uvmean_outflow(k) |
|---|
| 248 | ENDDO |
|---|
| 249 | |
|---|
| 250 | ENDIF |
|---|
| 251 | |
|---|
| 252 | IF ( outflow_r ) THEN |
|---|
| 253 | u(:,:,nxr+uxrp+1) = u(:,:,nxr+uxrp) |
|---|
| 254 | w(:,:,nxr+1) = 0.0 |
|---|
| 255 | ! |
|---|
| 256 | !-- Compute the mean horizontal wind parallel to and within the outflow |
|---|
| 257 | !-- wall and use this as boundary condition for v |
|---|
| 258 | #if defined( __parallel ) |
|---|
| 259 | CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, & |
|---|
| 260 | MPI_REAL, MPI_SUM, comm1dy, ierr ) |
|---|
| 261 | uvmean_outflow = uvmean_outflow / ( ny + 1.0 ) |
|---|
| 262 | #else |
|---|
| 263 | uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 ) |
|---|
| 264 | #endif |
|---|
| 265 | DO k = nzb, nzt+1 |
|---|
| 266 | v(k,:,nxr+1) = uvmean_outflow(k) |
|---|
| 267 | ENDDO |
|---|
| 268 | ENDIF |
|---|
| 269 | |
|---|
| 270 | ENDIF |
|---|
| 271 | |
|---|
| 272 | |
|---|
| 273 | END SUBROUTINE boundary_conds |
|---|
