Changeset 73
- Timestamp:
- Mar 20, 2007 8:33:14 AM (18 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 14 edited
-
CURRENT_MODIFICATIONS (modified) (4 diffs)
-
boundary_conds.f90 (modified) (11 diffs)
-
calc_precipitation.f90 (modified) (1 diff)
-
check_parameters.f90 (modified) (4 diffs)
-
exchange_horiz_2d.f90 (modified) (2 diffs)
-
init_3d_model.f90 (modified) (4 diffs)
-
init_pegrid.f90 (modified) (6 diffs)
-
modules.f90 (modified) (2 diffs)
-
pres.f90 (modified) (1 diff)
-
production_e.f90 (modified) (1 diff)
-
prognostic_equations.f90 (modified) (4 diffs)
-
read_3d_binary.f90 (modified) (4 diffs)
-
time_integration.f90 (modified) (1 diff)
-
write_3d_binary.f90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/CURRENT_MODIFICATIONS
r72 r73 33 33 Changed: 34 34 ------- 35 General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d 36 35 37 +age_m in particle_type 36 38 … … 44 46 Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc. 45 47 48 q is not allowed to become negative (prognostic_equations). 49 46 50 __vtk directives removed from main program. 47 51 48 52 The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV. 49 53 50 check_parameters, data_output_dvrp, data_output_ptseries, data_output_ts, flow_statistics, header, init_particles, init_3d_model, modules, palm, package_parin, parin, time_integration54 check_parameters, data_output_dvrp, data_output_ptseries, data_output_ts, exchange_horiz_2d, flow_statistics, header, init_particles, init_pegrid, init_3d_model, modules, palm, package_parin, parin, prognostic_equations, read_3d_binary, time_integration, write_3d_binary 51 55 52 56 … … 58 62 Bugfix in sample for reading user defined data from restart file (user_init) 59 63 64 Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model) 65 60 66 Check for possible negative humidities in the initial humidity profile. 61 67 … … 64 70 65 71 Makefile 66 check_parameters, init_1d_model, user_interface72 check_parameters, init_1d_model, init_3d_model, user_interface -
palm/trunk/SOURCE/boundary_conds.f90
r39 r73 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! The "main" part sets conditions for time level t+dt insteat of level t, 7 ! outflow boundary conditions changed from Neumann to radiation condition 7 8 ! 8 9 ! Former revisions: … … 51 52 INTEGER :: i, j, k 52 53 54 REAL :: c_max, c_u, c_v, c_w, denom 55 53 56 54 57 IF ( range == 'main') THEN … … 56 59 !-- Bottom boundary 57 60 IF ( ibc_uv_b == 0 ) THEN 58 u(nzb,:,:) = -u(nzb+1,:,:) ! satisfying the Dirichlet condition with 59 v(nzb,:,:) = -v(nzb+1,:,:) ! an extra layer below the surface where 60 ELSE ! the u and v component change their sign 61 u(nzb,:,:) = u(nzb+1,:,:) 62 v(nzb,:,:) = v(nzb+1,:,:) 61 ! 62 !-- Satisfying the Dirichlet condition with an extra layer below the 63 !-- surface where the u and v component change their sign 64 u_p(nzb,:,:) = -u_p(nzb+1,:,:) 65 v_p(nzb,:,:) = -v_p(nzb+1,:,:) 66 ELSE 67 u_p(nzb,:,:) = u_p(nzb+1,:,:) 68 v_p(nzb,:,:) = v_p(nzb+1,:,:) 63 69 ENDIF 64 70 DO i = nxl-1, nxr+1 65 71 DO j = nys-1, nyn+1 66 w (nzb_w_inner(j,i),j,i) = 0.072 w_p(nzb_w_inner(j,i),j,i) = 0.0 67 73 ENDDO 68 74 ENDDO … … 71 77 !-- Top boundary 72 78 IF ( ibc_uv_t == 0 ) THEN 73 u (nzt+1,:,:) = ug(nzt+1)74 v (nzt+1,:,:) = vg(nzt+1)79 u_p(nzt+1,:,:) = ug(nzt+1) 80 v_p(nzt+1,:,:) = vg(nzt+1) 75 81 ELSE 76 u (nzt+1,:,:) = u(nzt,:,:)77 v (nzt+1,:,:) = v(nzt,:,:)78 ENDIF 79 w (nzt:nzt+1,:,:) = 0.0 ! nzt is not a prognostic level (but cf. pres)82 u_p(nzt+1,:,:) = u_p(nzt,:,:) 83 v_p(nzt+1,:,:) = v_p(nzt,:,:) 84 ENDIF 85 w_p(nzt:nzt+1,:,:) = 0.0 ! nzt is not a prognostic level (but cf. pres) 80 86 81 87 ! 82 88 !-- Temperature at bottom boundary 83 89 IF ( ibc_pt_b == 0 ) THEN 84 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 85 DO i = nxl-1, nxr+1 86 DO j = nys-1, nyn+1 87 pt(nzb_s_inner(j,i),j,i) = pt_m(nzb_s_inner(j,i),j,i) 88 ENDDO 89 ENDDO 90 ELSE 91 DO i = nxl-1, nxr+1 92 DO j = nys-1, nyn+1 93 pt(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i),j,i) 94 ! pt_m is not used for Runge-Kutta 95 ENDDO 96 ENDDO 97 ENDIF 90 DO i = nxl-1, nxr+1 91 DO j = nys-1, nyn+1 92 pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i) 93 ENDDO 94 ENDDO 98 95 ELSE 99 96 DO i = nxl-1, nxr+1 100 97 DO j = nys-1, nyn+1 101 pt (nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i)+1,j,i)98 pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i) 102 99 ENDDO 103 100 ENDDO … … 107 104 !-- Temperature at top boundary 108 105 IF ( ibc_pt_t == 0 ) THEN 109 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 110 pt(nzt+1,:,:) = pt_m(nzt+1,:,:) 111 ELSE 112 pt(nzt+1,:,:) = pt_p(nzt+1,:,:) ! pt_m not used for Runge-Kutta 113 ENDIF 106 pt_p(nzt+1,:,:) = pt(nzt+1,:,:) 114 107 ELSEIF ( ibc_pt_t == 1 ) THEN 115 pt (nzt+1,:,:) = pt(nzt,:,:)108 pt_p(nzt+1,:,:) = pt_p(nzt,:,:) 116 109 ELSEIF ( ibc_pt_t == 2 ) THEN 117 pt (nzt+1,:,:) = pt(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)110 pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1) 118 111 ENDIF 119 112 … … 124 117 DO i = nxl-1, nxr+1 125 118 DO j = nys-1, nyn+1 126 e (nzb_s_inner(j,i),j,i) = e(nzb_s_inner(j,i)+1,j,i)127 ENDDO 128 ENDDO 129 e (nzt+1,:,:) = e(nzt,:,:)119 e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i) 120 ENDDO 121 ENDDO 122 e_p(nzt+1,:,:) = e_p(nzt,:,:) 130 123 ENDIF 131 124 … … 137 130 !-- Surface conditions for constant_moisture_flux 138 131 IF ( ibc_q_b == 0 ) THEN 139 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 140 DO i = nxl-1, nxr+1 141 DO j = nys-1, nyn+1 142 q(nzb_s_inner(j,i),j,i) = q_m(nzb_s_inner(j,i),j,i) 143 ENDDO 132 DO i = nxl-1, nxr+1 133 DO j = nys-1, nyn+1 134 q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i) 144 135 ENDDO 145 ELSE 146 DO i = nxl-1, nxr+1 147 DO j = nys-1, nyn+1 148 q(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i),j,i) 149 ENDDO ! q_m is not used for Runge-Kutta 150 ENDDO 151 ENDIF 136 ENDDO 152 137 ELSE 153 138 DO i = nxl-1, nxr+1 154 139 DO j = nys-1, nyn+1 155 q (nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i)+1,j,i)140 q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i) 156 141 ENDDO 157 142 ENDDO … … 159 144 ! 160 145 !-- Top boundary 161 q (nzt+1,:,:) = q(nzt,:,:) + bc_q_t_val * dzu(nzt+1)146 q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1) 162 147 ENDIF 163 148 … … 167 152 !-- divergence. Dirichlet conditions are used for all other quantities. 168 153 IF ( inflow_s ) THEN 169 v (:,nys,:) = v(:,nys-1,:)154 v_p(:,nys,:) = v_p(:,nys-1,:) 170 155 ELSEIF ( inflow_n ) THEN 171 v (:,nyn+vynp,:) = v(:,nyn+vynp+1,:)156 v_p(:,nyn+vynp,:) = v_p(:,nyn+vynp+1,:) 172 157 ELSEIF ( inflow_l ) THEN 173 u (:,:,nxl) = u(:,:,nxl-1)158 u_p(:,:,nxl) = u_p(:,:,nxl-1) 174 159 ELSEIF ( inflow_r ) THEN 175 u (:,:,nxr+uxrp) = u(:,:,nxr+uxrp+1)160 u_p(:,:,nxr+uxrp) = u_p(:,:,nxr+uxrp+1) 176 161 ENDIF 177 162 … … 179 164 !-- Lateral boundary conditions for scalar quantities at the outflow 180 165 IF ( outflow_s ) THEN 181 pt (:,nys-1,:) = pt(:,nys,:)182 IF ( .NOT. constant_diffusion ) e (:,nys-1,:) = e(:,nys,:)183 IF ( moisture .OR. passive_scalar ) q (:,nys-1,:) = q(:,nys,:)166 pt_p(:,nys-1,:) = pt_p(:,nys,:) 167 IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:) 168 IF ( moisture .OR. passive_scalar ) q_p(:,nys-1,:) = q_p(:,nys,:) 184 169 ELSEIF ( outflow_n ) THEN 185 pt (:,nyn+1,:) = pt(:,nyn,:)186 IF ( .NOT. constant_diffusion ) e (:,nyn+1,:) = e(:,nyn,:)187 IF ( moisture .OR. passive_scalar ) q (:,nyn+1,:) = q(:,nyn,:)170 pt_p(:,nyn+1,:) = pt_p(:,nyn,:) 171 IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:) 172 IF ( moisture .OR. passive_scalar ) q_p(:,nyn+1,:) = q_p(:,nyn,:) 188 173 ELSEIF ( outflow_l ) THEN 189 pt (:,:,nxl-1) = pt(:,:,nxl)190 IF ( .NOT. constant_diffusion ) e (:,:,nxl-1) = e(:,:,nxl)191 IF ( moisture .OR. passive_scalar ) q (:,:,nxl-1) = q(:,:,nxl)174 pt_p(:,:,nxl-1) = pt_p(:,:,nxl) 175 IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl) 176 IF ( moisture .OR. passive_scalar ) q_p(:,:,nxl-1) = q_p(:,:,nxl) 192 177 ELSEIF ( outflow_r ) THEN 193 pt (:,:,nxr+1) = pt(:,:,nxr)194 IF ( .NOT. constant_diffusion ) e (:,:,nxr+1) = e(:,:,nxr)195 IF ( moisture .OR. passive_scalar ) q (:,:,nxr+1) = q(:,:,nxr)178 pt_p(:,:,nxr+1) = pt_p(:,:,nxr) 179 IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr) 180 IF ( moisture .OR. passive_scalar ) q_p(:,:,nxr+1) = q_p(:,:,nxr) 196 181 ENDIF 197 182 … … 200 185 IF ( range == 'outflow_uvw' ) THEN 201 186 ! 202 !-- Horizontal boundary conditions for the velocities at the outflow. 203 !-- A Neumann condition is used for the wall normal velocity. The vertical 204 !-- velocity is assumed as zero and a horizontal average along the wall is 205 !-- used for the wall parallel horizontal velocity. The combination of all 206 !-- three conditions ensures that the velocity field is free of divergence 207 !-- at the outflow (uvmean_outflow_l is calculated in pres). 187 !-- Radiation boundary condition for the velocities at the respective outflow 208 188 IF ( outflow_s ) THEN 209 v(:,nys-1,:) = v(:,nys,:) 210 w(:,nys-1,:) = 0.0 211 ! 212 !-- Compute the mean horizontal wind parallel to and within the outflow 213 !-- wall and use this as boundary condition for u 214 #if defined( __parallel ) 215 CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, & 216 MPI_REAL, MPI_SUM, comm1dx, ierr ) 217 uvmean_outflow = uvmean_outflow / ( nx + 1.0 ) 218 #else 219 uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 ) 220 #endif 221 DO k = nzb, nzt+1 222 u(k,nys-1,:) = uvmean_outflow(k) 223 ENDDO 224 ENDIF 225 226 IF ( outflow_n ) THEN 227 v(:,nyn+vynp+1,:) = v(:,nyn+vynp,:) 228 w(:,nyn+1,:) = 0.0 229 ! 230 !-- Compute the mean horizontal wind parallel to and within the outflow 231 !-- wall and use this as boundary condition for u 232 #if defined( __parallel ) 233 CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, & 234 MPI_REAL, MPI_SUM, comm1dx, ierr ) 235 uvmean_outflow = uvmean_outflow / ( nx + 1.0 ) 236 #else 237 uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 ) 238 #endif 239 DO k = nzb, nzt+1 240 u(k,nyn+1,:) = uvmean_outflow(k) 241 ENDDO 189 ! v(:,nys-1,:) = v(:,nys,:) 190 ! w(:,nys-1,:) = 0.0 191 !! 192 !!-- Compute the mean horizontal wind parallel to and within the outflow 193 !!-- wall and use this as boundary condition for u 194 !#if defined( __parallel ) 195 ! CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, & 196 ! MPI_REAL, MPI_SUM, comm1dx, ierr ) 197 ! uvmean_outflow = uvmean_outflow / ( nx + 1.0 ) 198 !#else 199 ! uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 ) 200 !#endif 201 ! DO k = nzb, nzt+1 202 ! u(k,nys-1,:) = uvmean_outflow(k) 203 ! ENDDO 204 ENDIF 205 206 IF ( outflow_n .AND. & 207 intermediate_timestep_count == intermediate_timestep_count_max ) & 208 THEN 209 210 c_max = dy / dt_3d 211 212 DO i = nxl-1, nxr+1 213 DO k = nzb+1, nzt+1 214 215 ! 216 !-- First calculate the phase speeds for u,v, and w 217 denom = u_m_n(k,ny,i,-2) - u_m_n(k,ny-1,i,-2) 218 219 IF ( denom /= 0.0 ) THEN 220 c_u = -c_max * ( u_m_n(k,ny,i,-1)-u_m_n(k,ny,i,-2) ) / denom 221 IF ( c_u < 0.0 ) THEN 222 c_u = 0.0 223 ELSEIF ( c_u > c_max ) THEN 224 c_u = c_max 225 ENDIF 226 ELSE 227 c_u = c_max 228 ENDIF 229 230 denom = v_m_n(k,ny,i,-2) - v_m_n(k,ny-1,i,-2) 231 232 IF ( denom /= 0.0 ) THEN 233 c_v = -c_max * ( v_m_n(k,ny,i,-1)-v_m_n(k,ny,i,-2) ) / denom 234 IF ( c_v < 0.0 ) THEN 235 c_v = 0.0 236 ELSEIF ( c_v > c_max ) THEN 237 c_v = c_max 238 ENDIF 239 ELSE 240 c_v = c_max 241 ENDIF 242 243 denom = w_m_n(k,ny,i,-2) - w_m_n(k,ny-1,i,-2) 244 245 IF ( denom /= 0.0 ) THEN 246 c_w = -c_max * ( w_m_n(k,ny,i,-1)-w_m_n(k,ny,i,-2) ) / denom 247 IF ( c_w < 0.0 ) THEN 248 c_w = 0.0 249 ELSEIF ( c_w > c_max ) THEN 250 c_w = c_max 251 ENDIF 252 ELSE 253 c_w = c_max 254 ENDIF 255 256 ! 257 !-- Calculate the new velocities 258 u(k,ny+1,i) = u_m_n(k,ny+1,i,-1) - dt_3d * c_u * & 259 ( u_m_n(k,ny+1,i,-1) - u_m_n(k,ny,i,-1) ) * ddy 260 261 v(k,ny+1,i) = v_m_n(k,ny+1,i,-1) - dt_3d * c_v * & 262 ( v_m_n(k,ny+1,i,-1) - v_m_n(k,ny,i,-1) ) * ddy 263 264 w(k,ny+1,i) = w_m_n(k,ny+1,i,-1) - dt_3d * c_w * & 265 ( w_m_n(k,ny+1,i,-1) - w_m_n(k,ny,i,-1) ) * ddy 266 267 ! 268 !-- Swap timelevels for the next timestep 269 u_m_n(k,:,i,-2) = u_m_n(k,:,i,-1) 270 u_m_n(k,:,i,-1) = u(k,ny-1:ny+1,i) 271 v_m_n(k,:,i,-2) = v_m_n(k,:,i,-1) 272 v_m_n(k,:,i,-1) = v(k,ny-1:ny+1,i) 273 w_m_n(k,:,i,-2) = w_m_n(k,:,i,-1) 274 w_m_n(k,:,i,-1) = w(k,ny-1:ny+1,i) 275 276 ENDDO 277 ENDDO 278 279 ! 280 !-- Bottom boundary at the outflow 281 IF ( ibc_uv_b == 0 ) THEN 282 u(nzb,ny+1,:) = -u(nzb+1,ny+1,:) 283 v(nzb,ny+1,:) = -v(nzb+1,ny+1,:) 284 ELSE 285 u(nzb,ny+1,:) = u(nzb+1,ny+1,:) 286 v(nzb,ny+1,:) = v(nzb+1,ny+1,:) 287 ENDIF 288 w(nzb,ny+1,:) = 0.0 289 290 ! 291 !-- Top boundary at the outflow 292 IF ( ibc_uv_t == 0 ) THEN 293 u(nzt+1,ny+1,:) = ug(nzt+1) 294 v(nzt+1,ny+1,:) = vg(nzt+1) 295 ELSE 296 u(nzt+1,ny+1,:) = u(nzt,nyn+1,:) 297 v(nzt+1,ny+1,:) = v(nzt,nyn+1,:) 298 ENDIF 299 w(nzt:nzt+1,ny+1,:) = 0.0 300 242 301 ENDIF 243 302 244 303 IF ( outflow_l ) THEN 245 u(:,:,nxl-1) = u(:,:,nxl)246 w(:,:,nxl-1) = 0.0304 ! u(:,:,nxl-1) = u(:,:,nxl) 305 ! w(:,:,nxl-1) = 0.0 247 306 ! 248 307 !-- Compute the mean horizontal wind parallel to and within the outflow 249 308 !-- wall and use this as boundary condition for v 250 #if defined( __parallel ) 251 CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, & 252 MPI_REAL, MPI_SUM, comm1dy, ierr ) 253 uvmean_outflow = uvmean_outflow / ( ny + 1.0 ) 254 #else 255 uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 ) 256 #endif 257 DO k = nzb, nzt+1 258 v(k,:,nxl-1) = uvmean_outflow(k) 259 ENDDO 260 261 ENDIF 262 263 IF ( outflow_r ) THEN 264 u(:,:,nxr+uxrp+1) = u(:,:,nxr+uxrp) 265 w(:,:,nxr+1) = 0.0 266 ! 267 !-- Compute the mean horizontal wind parallel to and within the outflow 268 !-- wall and use this as boundary condition for v 269 #if defined( __parallel ) 270 CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, & 271 MPI_REAL, MPI_SUM, comm1dy, ierr ) 272 uvmean_outflow = uvmean_outflow / ( ny + 1.0 ) 273 #else 274 uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 ) 275 #endif 276 DO k = nzb, nzt+1 277 v(k,:,nxr+1) = uvmean_outflow(k) 278 ENDDO 309 !#if defined( __parallel ) 310 ! CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, & 311 ! MPI_REAL, MPI_SUM, comm1dy, ierr ) 312 ! uvmean_outflow = uvmean_outflow / ( ny + 1.0 ) 313 !#else 314 ! uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 ) 315 !#endif 316 ! DO k = nzb, nzt+1 317 ! v(k,:,nxl-1) = uvmean_outflow(k) 318 ! ENDDO 319 ! 320 ENDIF 321 322 IF ( outflow_r .AND. & 323 intermediate_timestep_count == intermediate_timestep_count_max ) & 324 THEN 325 326 c_max = dx / dt_3d 327 328 DO j = nys-1, nyn+1 329 DO k = nzb+1, nzt+1 330 331 ! 332 !-- First calculate the phase speeds for u,v, and w 333 denom = u_m_r(k,j,nx,-2) - u_m_r(k,j,nx-1,-2) 334 335 IF ( denom /= 0.0 ) THEN 336 c_u = -c_max * ( u_m_r(k,j,nx,-1)-u_m_r(k,j,nx,-2) ) / denom 337 IF ( c_u < 0.0 ) THEN 338 c_u = 0.0 339 ELSEIF ( c_u > c_max ) THEN 340 c_u = c_max 341 ENDIF 342 ELSE 343 c_u = c_max 344 ENDIF 345 346 denom = v_m_r(k,j,nx,-2) - v_m_r(k,j,nx-1,-2) 347 348 IF ( denom /= 0.0 ) THEN 349 c_v = -c_max * ( v_m_r(k,j,nx,-1)-v_m_r(k,j,nx,-2) ) / denom 350 IF ( c_v < 0.0 ) THEN 351 c_v = 0.0 352 ELSEIF ( c_v > c_max ) THEN 353 c_v = c_max 354 ENDIF 355 ELSE 356 c_v = c_max 357 ENDIF 358 359 denom = w_m_r(k,j,nx,-2) - w_m_r(k,j,nx-1,-2) 360 361 IF ( denom /= 0.0 ) THEN 362 c_w = -c_max * ( w_m_r(k,j,nx,-1)-w_m_n(k,j,nx,-2) ) / denom 363 IF ( c_w < 0.0 ) THEN 364 c_w = 0.0 365 ELSEIF ( c_w > c_max ) THEN 366 c_w = c_max 367 ENDIF 368 ELSE 369 c_w = c_max 370 ENDIF 371 372 ! 373 !-- Calculate the new velocities 374 u(k,j,nx+1) = u_m_r(k,j,nx+1,-1) - dt_3d * c_u * & 375 ( u_m_r(k,j,nx+1,-1) - u_m_r(k,j,nx,-1) ) * ddx 376 377 v(k,j,nx+1) = v_m_r(k,j,nx+1,-1) - dt_3d * c_v * & 378 ( v_m_r(k,j,nx+1,-1) - v_m_r(k,j,nx,-1) ) * ddx 379 380 w(k,j,nx+1) = w_m_r(k,j,nx+1,-1) - dt_3d * c_w * & 381 ( w_m_r(k,j,nx+1,-1) - w_m_r(k,j,nx,-1) ) * ddx 382 383 ! 384 !-- Swap timelevels for the next timestep 385 u_m_r(k,j,:,-2) = u_m_r(k,j,:,-1) 386 u_m_r(k,j,:,-1) = u(k,j,nx-1:nx+1) 387 v_m_r(k,j,:,-2) = v_m_r(k,j,:,-1) 388 v_m_r(k,j,:,-1) = v(k,j,nx-1:nx+1) 389 w_m_r(k,j,:,-2) = w_m_r(k,j,:,-1) 390 w_m_r(k,j,:,-1) = w(k,j,nx-1:nx+1) 391 392 ENDDO 393 ENDDO 394 395 ! 396 !-- Bottom boundary at the outflow 397 IF ( ibc_uv_b == 0 ) THEN 398 u(nzb,ny+1,:) = -u(nzb+1,ny+1,:) 399 v(nzb,ny+1,:) = -v(nzb+1,ny+1,:) 400 ELSE 401 u(nzb,ny+1,:) = u(nzb+1,ny+1,:) 402 v(nzb,ny+1,:) = v(nzb+1,ny+1,:) 403 ENDIF 404 w(nzb,ny+1,:) = 0.0 405 406 ! 407 !-- Top boundary at the outflow 408 IF ( ibc_uv_t == 0 ) THEN 409 u(nzt+1,ny+1,:) = ug(nzt+1) 410 v(nzt+1,ny+1,:) = vg(nzt+1) 411 ELSE 412 u(nzt+1,ny+1,:) = u(nzt,nyn+1,:) 413 v(nzt+1,ny+1,:) = v(nzt,nyn+1,:) 414 ENDIF 415 w(nzt:nzt+1,ny+1,:) = 0.0 416 279 417 ENDIF 280 418 -
palm/trunk/SOURCE/calc_precipitation.f90
r72 r73 71 71 tend(k,j,i) = tend(k,j,i) - dqdt_precip 72 72 ! 73 !-- Precipitation rate in (kg * 0.001) / m**2 / s (because 1kg 74 !-- gives 1 mm) 73 !-- Precipitation rate in kg / m**2 / s (= mm/s) 75 74 precipitation_rate(j,i) = precipitation_rate(j,i) + & 76 dqdt_precip * dzw(k) * 0.00175 dqdt_precip * dzw(k) 77 76 78 77 ENDDO 79 78 ! 80 !-- Sum up the precipitation amount (unit kg * 0.001 / m**2)79 !-- Sum up the precipitation amount, unit kg / m**2 (= mm) 81 80 IF ( intermediate_timestep_count == & 82 81 intermediate_timestep_count_max .AND. & -
palm/trunk/SOURCE/check_parameters.f90
r72 r73 8 8 ! and pt_reference are checked, 9 9 ! output of precipitation amount/rate and roughnes length + check 10 ! possible negative humidities are avoided in initial profile 10 ! possible negative humidities are avoided in initial profile, 11 ! dirichlet/neumann changed to dirichlet/radiation, etc. 11 12 ! 12 13 ! Former revisions: … … 705 706 !-- Check boundary conditions and set internal variables: 706 707 !-- Lateral boundary conditions 707 IF ( bc_lr /= 'cyclic' .AND. bc_lr /= 'dirichlet/ neumann' .AND. &708 bc_lr /= ' neumann/dirichlet' ) THEN708 IF ( bc_lr /= 'cyclic' .AND. bc_lr /= 'dirichlet/radiation' .AND. & 709 bc_lr /= 'radiation/dirichlet' ) THEN 709 710 IF ( myid == 0 ) THEN 710 711 PRINT*, '+++ check_parameters:' … … 713 714 CALL local_stop 714 715 ENDIF 715 IF ( bc_ns /= 'cyclic' .AND. bc_ns /= 'dirichlet/ neumann' .AND. &716 bc_ns /= ' neumann/dirichlet' ) THEN716 IF ( bc_ns /= 'cyclic' .AND. bc_ns /= 'dirichlet/radiation' .AND. & 717 bc_ns /= 'radiation/dirichlet' ) THEN 717 718 IF ( myid == 0 ) THEN 718 719 PRINT*, '+++ check_parameters:' … … 2284 2285 ENDIF 2285 2286 2286 IF ( outflow_l) THEN2287 IF ( bc_lr == 'radiation/dirichlet' ) THEN 2287 2288 dist_nxr = nx - inflow_disturbance_begin 2288 2289 dist_nxl(1) = nx - inflow_disturbance_end 2289 ELSEIF ( outflow_r) THEN2290 ELSEIF ( bc_lr == 'dirichlet/radiation' ) THEN 2290 2291 dist_nxl = inflow_disturbance_begin 2291 2292 dist_nxr(1) = inflow_disturbance_end 2292 ELSEIF ( outflow_s ) THEN 2293 ENDIF 2294 IF ( bc_ns == 'dirichlet/radiation' ) THEN 2293 2295 dist_nyn = ny - inflow_disturbance_begin 2294 2296 dist_nys(1) = ny - inflow_disturbance_end 2295 ELSEIF ( outflow_n) THEN2297 ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN 2296 2298 dist_nys = inflow_disturbance_begin 2297 2299 dist_nyn(1) = inflow_disturbance_end -
palm/trunk/SOURCE/exchange_horiz_2d.f90
r4 r73 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary 7 ! conditions 7 8 ! 8 9 ! Former revisions: … … 98 99 #endif 99 100 101 ! 102 !-- Neumann-conditions at inflow/outflow in case of non-cyclic boundary 103 !-- conditions 104 IF ( inflow_l .OR. outflow_l ) ar(:,nxl-1) = ar(:,nxl) 105 IF ( inflow_r .OR. outflow_r ) ar(:,nxr+1) = ar(:,nxr) 106 IF ( inflow_s .OR. outflow_s ) ar(nys-1,:) = ar(nys,:) 107 IF ( inflow_n .OR. outflow_n ) ar(nyn+1,:) = ar(nyn,:) 108 100 109 CALL cpu_log( log_point_s(13), 'exchange_horiz_2d', 'stop' ) 101 110 -
palm/trunk/SOURCE/init_3d_model.f90
r72 r73 7 7 ! Actual revisions: 8 8 ! ----------------- 9 ! Arrays for radiation boundary conditions are allocated (u_m_l, u_m_r, etc.), 10 ! bugfix for cases with the outflow damping layer extending over more than one 11 ! subdomain, 9 12 ! New initializing action "by_user" calls user_init_3d_model, 10 13 ! precipitation_amount/rate, ts_value are allocated, +module netcdf_control, … … 206 209 207 210 ! 211 !-- Arrays to store velocity data from t-dt needed for radiation boundary 212 !-- conditions 213 IF ( outflow_l ) THEN 214 ALLOCATE( u_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1,-2:-1), & 215 v_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1,-2:-1), & 216 w_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1,-2:-1) ) 217 ENDIF 218 IF ( outflow_r ) THEN 219 ALLOCATE( u_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1,-2:-1), & 220 v_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1,-2:-1), & 221 w_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1,-2:-1) ) 222 ENDIF 223 IF ( outflow_s ) THEN 224 ALLOCATE( u_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1,-2:-1), & 225 v_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1,-2:-1), & 226 w_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1,-2:-1) ) 227 ENDIF 228 IF ( outflow_n ) THEN 229 ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1,-2:-1), & 230 v_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1,-2:-1), & 231 w_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1,-2:-1) ) 232 ENDIF 233 234 ! 208 235 !-- Initial assignment of the pointers 209 236 IF ( timestep_scheme(1:5) /= 'runge' ) THEN … … 726 753 ENDIF 727 754 755 ! 756 !-- Initialize old timelevels needed for radiation boundary conditions 757 IF ( outflow_l ) THEN 758 u_m_l(:,:,:,-2) = u(:,:,-1:1) 759 v_m_l(:,:,:,-2) = v(:,:,-1:1) 760 w_m_l(:,:,:,-2) = w(:,:,-1:1) 761 u_m_l(:,:,:,-1) = u(:,:,-1:1) 762 v_m_l(:,:,:,-1) = v(:,:,-1:1) 763 w_m_l(:,:,:,-1) = w(:,:,-1:1) 764 ENDIF 765 IF ( outflow_r ) THEN 766 u_m_r(:,:,:,-2) = u(:,:,nx-1:nx+1) 767 v_m_r(:,:,:,-2) = v(:,:,nx-1:nx+1) 768 w_m_r(:,:,:,-2) = w(:,:,nx-1:nx+1) 769 u_m_r(:,:,:,-1) = u(:,:,nx-1:nx+1) 770 v_m_r(:,:,:,-1) = v(:,:,nx-1:nx+1) 771 w_m_r(:,:,:,-1) = w(:,:,nx-1:nx+1) 772 ENDIF 773 IF ( outflow_s ) THEN 774 u_m_s(:,:,:,-2) = u(:,-1:1,:) 775 v_m_s(:,:,:,-2) = v(:,-1:1,:) 776 w_m_s(:,:,:,-2) = w(:,-1:1,:) 777 u_m_s(:,:,:,-1) = u(:,-1:1,:) 778 v_m_s(:,:,:,-1) = v(:,-1:1,:) 779 w_m_s(:,:,:,-1) = w(:,-1:1,:) 780 ENDIF 781 IF ( outflow_n ) THEN 782 u_m_n(:,:,:,-2) = u(:,ny-1:ny+1,:) 783 v_m_n(:,:,:,-2) = v(:,ny-1:ny+1,:) 784 w_m_n(:,:,:,-2) = w(:,ny-1:ny+1,:) 785 u_m_n(:,:,:,-1) = u(:,ny-1:ny+1,:) 786 v_m_n(:,:,:,-1) = v(:,ny-1:ny+1,:) 787 w_m_n(:,:,:,-1) = w(:,ny-1:ny+1,:) 788 ENDIF 789 728 790 ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' ) & 729 791 THEN … … 788 850 !-- non-cyclic lateral boundaries. A linear increase is assumed over the first 789 851 !-- half of the width of the damping layer 790 IF ( bc_lr /= 'cyclic' ) THEN852 IF ( bc_lr == 'dirichlet/radiation' ) THEN 791 853 792 854 DO i = nxl-1, nxr+1 793 794 IF ( outflow_r ) THEN 795 796 IF ( i >= nx - outflow_damping_width ) THEN 797 km_damp_x(i) = km_damp_max * MIN( 1.0, & 798 ( i - ( nx - outflow_damping_width ) ) / & 799 REAL( outflow_damping_width/2 ) & 800 ) 801 ELSE 802 km_damp_x(i) = 0.0 803 ENDIF 804 805 ELSEIF ( outflow_l ) THEN 806 807 IF ( i <= outflow_damping_width ) THEN 808 km_damp_x(i) = km_damp_max * MIN( 1.0, & 809 ( outflow_damping_width - i ) / & 810 REAL( outflow_damping_width/2 ) & 811 ) 812 ELSE 813 km_damp_x(i) = 0.0 814 ENDIF 815 816 ENDIF 817 855 IF ( i >= nx - outflow_damping_width ) THEN 856 km_damp_x(i) = km_damp_max * MIN( 1.0, & 857 ( i - ( nx - outflow_damping_width ) ) / & 858 REAL( outflow_damping_width/2 ) & 859 ) 860 ELSE 861 km_damp_x(i) = 0.0 862 ENDIF 818 863 ENDDO 819 ENDIF 820 821 IF ( bc_ns /= 'cyclic' ) THEN 864 865 ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN 866 867 DO i = nxl-1, nxr+1 868 IF ( i <= outflow_damping_width ) THEN 869 km_damp_x(i) = km_damp_max * MIN( 1.0, & 870 ( outflow_damping_width - i ) / & 871 REAL( outflow_damping_width/2 ) & 872 ) 873 ELSE 874 km_damp_x(i) = 0.0 875 ENDIF 876 ENDDO 877 878 ENDIF 879 880 IF ( bc_ns == 'dirichlet/radiation' ) THEN 822 881 823 882 DO j = nys-1, nyn+1 824 825 IF ( outflow_n ) THEN 826 827 IF ( j >= ny - outflow_damping_width ) THEN 828 km_damp_y(j) = km_damp_max * MIN( 1.0, & 829 ( j - ( ny - outflow_damping_width ) ) / & 830 REAL( outflow_damping_width/2 ) & 831 ) 832 ELSE 833 km_damp_y(j) = 0.0 834 ENDIF 835 836 ELSEIF ( outflow_s ) THEN 837 838 IF ( j <= outflow_damping_width ) THEN 839 km_damp_y(j) = km_damp_max * MIN( 1.0, & 840 ( outflow_damping_width - j ) / & 841 REAL( outflow_damping_width/2 ) & 842 ) 843 ELSE 844 km_damp_y(j) = 0.0 845 ENDIF 846 847 ENDIF 848 883 IF ( j >= ny - outflow_damping_width ) THEN 884 km_damp_y(j) = km_damp_max * MIN( 1.0, & 885 ( j - ( ny - outflow_damping_width ) ) / & 886 REAL( outflow_damping_width/2 ) & 887 ) 888 ELSE 889 km_damp_y(j) = 0.0 890 ENDIF 849 891 ENDDO 892 893 ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN 894 895 DO j = nys-1, nyn+1 896 IF ( j <= outflow_damping_width ) THEN 897 km_damp_y(j) = km_damp_max * MIN( 1.0, & 898 ( outflow_damping_width - j ) / & 899 REAL( outflow_damping_width/2 ) & 900 ) 901 ELSE 902 km_damp_y(j) = 0.0 903 ENDIF 904 ENDDO 905 850 906 ENDIF 851 907 -
palm/trunk/SOURCE/init_pegrid.f90
r4 r73 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! uxrp, vynp eliminated, 7 ! dirichlet/neumann changed to dirichlet/radiation, etc. 7 8 ! 8 9 ! Former revisions: … … 495 496 ! 496 497 !-- For non-cyclic boundaries extend array u (v) by one gridpoint 497 IF ( bc_lr /= 'cyclic' ) uxrp = 1498 IF ( bc_ns /= 'cyclic' ) vynp = 1498 ! IF ( bc_lr /= 'cyclic' ) uxrp = 1 499 ! IF ( bc_ns /= 'cyclic' ) vynp = 1 499 500 500 501 ! … … 743 744 !-- by one gridpoint on the left/rightmost (northest/southest) processor 744 745 IF ( pleft == MPI_PROC_NULL ) THEN 745 IF ( bc_lr == 'dirichlet/ neumann' ) THEN746 IF ( bc_lr == 'dirichlet/radiation' ) THEN 746 747 inflow_l = .TRUE. 747 ELSEIF ( bc_lr == ' neumann/dirichlet' ) THEN748 ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN 748 749 outflow_l = .TRUE. 749 750 ENDIF … … 751 752 752 753 IF ( pright == MPI_PROC_NULL ) THEN 753 IF ( bc_lr == 'dirichlet/ neumann' ) THEN754 IF ( bc_lr == 'dirichlet/radiation' ) THEN 754 755 outflow_r = .TRUE. 755 ELSEIF ( bc_lr == ' neumann/dirichlet' ) THEN756 ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN 756 757 inflow_r = .TRUE. 757 758 ENDIF 758 uxrp = 1759 ! uxrp = 1 759 760 ENDIF 760 761 761 762 IF ( psouth == MPI_PROC_NULL ) THEN 762 IF ( bc_ns == 'dirichlet/ neumann' ) THEN763 IF ( bc_ns == 'dirichlet/radiation' ) THEN 763 764 outflow_s = .TRUE. 764 ELSEIF ( bc_ns == ' neumann/dirichlet' ) THEN765 ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN 765 766 inflow_s = .TRUE. 766 767 ENDIF … … 768 769 769 770 IF ( pnorth == MPI_PROC_NULL ) THEN 770 IF ( bc_ns == 'dirichlet/ neumann' ) THEN771 IF ( bc_ns == 'dirichlet/radiation' ) THEN 771 772 inflow_n = .TRUE. 772 ELSEIF ( bc_ns == ' neumann/dirichlet' ) THEN773 ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN 773 774 outflow_n = .TRUE. 774 775 ENDIF 775 vynp = 1776 ! vynp = 1 776 777 ENDIF 777 778 … … 787 788 ENDIF 788 789 #else 789 IF ( bc_lr == 'dirichlet/ neumann' ) THEN790 IF ( bc_lr == 'dirichlet/radiation' ) THEN 790 791 inflow_l = .TRUE. 791 792 outflow_r = .TRUE. 792 uxrp = 1793 ELSEIF ( bc_lr == ' neumann/dirichlet' ) THEN793 ! uxrp = 1 794 ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN 794 795 outflow_l = .TRUE. 795 796 inflow_r = .TRUE. 796 797 ENDIF 797 798 798 IF ( bc_ns == 'dirichlet/ neumann' ) THEN799 IF ( bc_ns == 'dirichlet/radiation' ) THEN 799 800 inflow_n = .TRUE. 800 801 outflow_s = .TRUE. 801 ELSEIF ( bc_ns == ' neumann/dirichlet' ) THEN802 ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN 802 803 outflow_n = .TRUE. 803 804 inflow_s = .TRUE. 804 vynp = 1805 ! vynp = 1 805 806 ENDIF 806 807 #endif -
palm/trunk/SOURCE/modules.f90
r72 r73 6 6 ! ----------------- 7 7 ! +arrays precipitation_amount, precipitation_rate, precipitation_rate_av, 8 ! rif_wall, z0_av 8 ! rif_wall, z0_av, +arrays u_m_l, u_m_r, etc. for radiation boundary conditions, 9 9 ! +loop_optimization, netcdf_64bit_3d, zu_s_inner, zw_w_inner, id_var_zusi_*, 10 10 ! id_var_zwwi_*, ts_value, u_nzb_p1_for_vfc, v_nzb_p1_for_vfc, pt_reference, … … 111 111 vpt, vpt_m, w, w_m, w_p 112 112 113 REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: rif_wall 113 REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: rif_wall, u_m_l, u_m_n, u_m_r, & 114 u_m_s, v_m_l, v_m_n, v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s 114 115 115 116 -
palm/trunk/SOURCE/pres.f90
r4 r73 396 396 !-- of non-cyclic boundary conditions). The respective mean velocity 397 397 !-- is calculated from this in routine boundary_conds. 398 IF ( outflow_l .AND. i == nxl ) THEN399 !$OMP CRITICAL400 DO k = nzb, nzt+1401 uvmean_outflow_l(k) = uvmean_outflow_l(k) + v(k,j,nxl)402 ENDDO403 !$OMP END CRITICAL404 ELSEIF ( outflow_r .AND. i == nxr ) THEN405 !$OMP CRITICAL406 DO k = nzb, nzt+1407 uvmean_outflow_l(k) = uvmean_outflow_l(k) + v(k,j,nxr)408 ENDDO409 !$OMP END CRITICAL410 ELSEIF ( outflow_s .AND. j == nys ) THEN411 !$OMP CRITICAL412 DO k = nzb, nzt+1413 uvmean_outflow_l(k) = uvmean_outflow_l(k) + u(k,nys,i)414 ENDDO415 !$OMP END CRITICAL416 ELSEIF ( outflow_n .AND. j == nyn ) THEN417 !$OMP CRITICAL418 DO k = nzb, nzt+1419 uvmean_outflow_l(k) = uvmean_outflow_l(k) + u(k,nyn,i)420 ENDDO421 !$OMP END CRITICAL422 ENDIF398 ! IF ( outflow_l .AND. i == nxl ) THEN 399 ! !$OMP CRITICAL 400 ! DO k = nzb, nzt+1 401 ! uvmean_outflow_l(k) = uvmean_outflow_l(k) + v(k,j,nxl) 402 ! ENDDO 403 ! !$OMP END CRITICAL 404 ! ELSEIF ( outflow_r .AND. i == nxr ) THEN 405 ! !$OMP CRITICAL 406 ! DO k = nzb, nzt+1 407 ! uvmean_outflow_l(k) = uvmean_outflow_l(k) + v(k,j,nxr) 408 ! ENDDO 409 ! !$OMP END CRITICAL 410 ! ELSEIF ( outflow_s .AND. j == nys ) THEN 411 ! !$OMP CRITICAL 412 ! DO k = nzb, nzt+1 413 ! uvmean_outflow_l(k) = uvmean_outflow_l(k) + u(k,nys,i) 414 ! ENDDO 415 ! !$OMP END CRITICAL 416 ! ELSEIF ( outflow_n .AND. j == nyn ) THEN 417 ! !$OMP CRITICAL 418 ! DO k = nzb, nzt+1 419 ! uvmean_outflow_l(k) = uvmean_outflow_l(k) + u(k,nyn,i) 420 ! ENDDO 421 ! !$OMP END CRITICAL 422 ! ENDIF 423 423 424 424 ! -
palm/trunk/SOURCE/production_e.f90
r72 r73 910 910 !-- of the eddy diffusivity by km = u* * kappa * zp / phi_m. 911 911 !$OMP PARALLEL DO PRIVATE( ku, kv ) 912 DO i = nxl, nxr +1913 DO j = nys, nyn +1912 DO i = nxl, nxr 913 DO j = nys, nyn 914 914 915 915 ku = nzb_u_inner(j,i)+1 -
palm/trunk/SOURCE/prognostic_equations.f90
r70 r73 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! checking for negative q and limiting for positive values, 6 7 ! z0 removed from arguments in calls of diffusion_u/v/w, 7 8 ! subroutine names changed to .._noopt, .._cache, and .._vector … … 481 482 ) - & 482 483 tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) ) 484 IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i) 483 485 ENDDO 484 486 … … 915 917 ) - & 916 918 tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) ) 919 IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i) 917 920 ENDDO 918 921 … … 1419 1422 ) - & 1420 1423 tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) ) 1424 IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i) 1421 1425 ENDDO 1422 1426 ENDDO -
palm/trunk/SOURCE/read_3d_binary.f90
r72 r73 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! +precipitation_amount, precipitation_rate_av, rif_wall, z0_av 6 ! +precipitation_amount, precipitation_rate_av, rif_wall, u_m_l, u_m_r, etc., 7 ! z0_av 7 8 ! 8 9 ! Former revisions: … … 302 303 CASE ( 'u_m' ) 303 304 READ ( 13 ) u_m 305 CASE ( 'u_m_l' ) 306 READ ( 13 ) u_m_l 307 CASE ( 'u_m_n' ) 308 READ ( 13 ) u_m_n 309 CASE ( 'u_m_r' ) 310 READ ( 13 ) u_m_r 311 CASE ( 'u_m_s' ) 312 READ ( 13 ) u_m_s 304 313 CASE ( 'us' ) 305 314 READ ( 13 ) us … … 322 331 CASE ( 'v_m' ) 323 332 READ (13 ) v_m 333 CASE ( 'v_m_l' ) 334 READ ( 13 ) v_m_l 335 CASE ( 'v_m_n' ) 336 READ ( 13 ) v_m_n 337 CASE ( 'v_m_r' ) 338 READ ( 13 ) v_m_r 339 CASE ( 'v_m_s' ) 340 READ ( 13 ) v_m_s 324 341 CASE ( 'vpt' ) 325 342 READ ( 13 ) vpt … … 340 357 CASE ( 'w_m' ) 341 358 READ ( 13 ) w_m 359 CASE ( 'w_m_l' ) 360 READ ( 13 ) w_m_l 361 CASE ( 'w_m_n' ) 362 READ ( 13 ) w_m_n 363 CASE ( 'w_m_r' ) 364 READ ( 13 ) w_m_r 365 CASE ( 'w_m_s' ) 366 READ ( 13 ) w_m_s 342 367 CASE ( 'z0' ) 343 368 READ ( 13 ) z0 -
palm/trunk/SOURCE/time_integration.f90
r63 r73 158 158 159 159 ! 160 !-- Swap the time levels in preparation for the next time step.161 CALL swap_timelevel162 163 !164 160 !-- Boundary conditions for the prognostic quantities (except of the 165 161 !-- velocities at the outflow in case of a non-cyclic lateral wall) 166 162 CALL boundary_conds( 'main' ) 163 164 ! 165 !-- Swap the time levels in preparation for the next time step. 166 CALL swap_timelevel 167 167 168 168 ! -
palm/trunk/SOURCE/write_3d_binary.f90
r72 r73 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! +precipitation_amount, precipitation_rate_av, rif_wall, z0_av 6 ! +precipitation_amount, precipitation_rate_av, rif_wall, u_m_l, u_m_r, etc., 7 ! z0_av 7 8 ! 8 9 ! Former revisions: … … 169 170 ENDIF 170 171 WRITE ( 14 ) 'u_m '; WRITE ( 14 ) u_m 172 IF ( ALLOCATED( u_m_l ) ) THEN 173 WRITE ( 14 ) 'u_m_l '; WRITE ( 14 ) u_m_l 174 ENDIF 175 IF ( ALLOCATED( u_m_n ) ) THEN 176 WRITE ( 14 ) 'u_m_n '; WRITE ( 14 ) u_m_n 177 ENDIF 178 IF ( ALLOCATED( u_m_r ) ) THEN 179 WRITE ( 14 ) 'u_m_r '; WRITE ( 14 ) u_m_r 180 ENDIF 181 IF ( ALLOCATED( u_m_s ) ) THEN 182 WRITE ( 14 ) 'u_m_s '; WRITE ( 14 ) u_m_s 183 ENDIF 171 184 WRITE ( 14 ) 'us '; WRITE ( 14 ) us 172 185 WRITE ( 14 ) 'usws '; WRITE ( 14 ) usws … … 184 197 ENDIF 185 198 WRITE ( 14 ) 'v_m '; WRITE ( 14 ) v_m 199 IF ( ALLOCATED( v_m_l ) ) THEN 200 WRITE ( 14 ) 'v_m_l '; WRITE ( 14 ) v_m_l 201 ENDIF 202 IF ( ALLOCATED( v_m_n ) ) THEN 203 WRITE ( 14 ) 'v_m_n '; WRITE ( 14 ) v_m_n 204 ENDIF 205 IF ( ALLOCATED( v_m_r ) ) THEN 206 WRITE ( 14 ) 'v_m_r '; WRITE ( 14 ) v_m_r 207 ENDIF 208 IF ( ALLOCATED( v_m_s ) ) THEN 209 WRITE ( 14 ) 'v_m_s '; WRITE ( 14 ) v_m_s 210 ENDIF 186 211 IF ( moisture ) THEN 187 212 WRITE ( 14 ) 'vpt '; WRITE ( 14 ) vpt … … 202 227 ENDIF 203 228 WRITE ( 14 ) 'w_m '; WRITE ( 14 ) w_m 229 IF ( ALLOCATED( w_m_l ) ) THEN 230 WRITE ( 14 ) 'w_m_l '; WRITE ( 14 ) w_m_l 231 ENDIF 232 IF ( ALLOCATED( w_m_n ) ) THEN 233 WRITE ( 14 ) 'w_m_n '; WRITE ( 14 ) w_m_n 234 ENDIF 235 IF ( ALLOCATED( w_m_r ) ) THEN 236 WRITE ( 14 ) 'w_m_r '; WRITE ( 14 ) w_m_r 237 ENDIF 238 IF ( ALLOCATED( w_m_s ) ) THEN 239 WRITE ( 14 ) 'w_m_s '; WRITE ( 14 ) w_m_s 240 ENDIF 204 241 WRITE ( 14 ) 'z0 '; WRITE ( 14 ) z0 205 242 IF ( ALLOCATED( z0_av ) ) THEN
Note: See TracChangeset
for help on using the changeset viewer.
