Ignore:
Timestamp:
Mar 22, 2007 9:54:05 AM (15 years ago)
Author:
raasch
Message:

preliminary update for changes concerning non-cyclic boundary conditions

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/boundary_conds.f90

    r73 r75  
    44! Actual revisions:
    55! -----------------
    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
     6! The "main" part sets conditions for time level t+dt instead of level t,
     7! outflow boundary conditions changed from Neumann to radiation condition,
     8! uxrp, vynp eliminated, moisture renamed humidity
    89!
    910! Former revisions:
     
    126127!--    Boundary conditions for total water content or scalar,
    127128!--    bottom and surface boundary (see also temperature)
    128        IF ( moisture  .OR.  passive_scalar )  THEN
    129 !
    130 !--       Surface conditions for constant_moisture_flux
     129       IF ( humidity  .OR.  passive_scalar )  THEN
     130!
     131!--       Surface conditions for constant_humidity_flux
    131132          IF ( ibc_q_b == 0 ) THEN
    132133             DO  i = nxl-1, nxr+1
     
    154155          v_p(:,nys,:) = v_p(:,nys-1,:)
    155156       ELSEIF ( inflow_n ) THEN
    156           v_p(:,nyn+vynp,:) = v_p(:,nyn+vynp+1,:)
     157          v_p(:,nyn,:) = v_p(:,nyn+1,:)
    157158       ELSEIF ( inflow_l ) THEN
    158159          u_p(:,:,nxl) = u_p(:,:,nxl-1)
    159160       ELSEIF ( inflow_r ) THEN
    160           u_p(:,:,nxr+uxrp) = u_p(:,:,nxr+uxrp+1)
     161          u_p(:,:,nxr) = u_p(:,:,nxr+1)
    161162       ENDIF
    162163
     
    166167          pt_p(:,nys-1,:)     = pt_p(:,nys,:)
    167168          IF ( .NOT. constant_diffusion     )  e_p(:,nys-1,:) = e_p(:,nys,:)
    168           IF ( moisture .OR. passive_scalar )  q_p(:,nys-1,:) = q_p(:,nys,:)
     169          IF ( humidity .OR. passive_scalar )  q_p(:,nys-1,:) = q_p(:,nys,:)
    169170       ELSEIF ( outflow_n )  THEN
    170171          pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
    171172          IF ( .NOT. constant_diffusion     )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
    172           IF ( moisture .OR. passive_scalar )  q_p(:,nyn+1,:) = q_p(:,nyn,:)
     173          IF ( humidity .OR. passive_scalar )  q_p(:,nyn+1,:) = q_p(:,nyn,:)
    173174       ELSEIF ( outflow_l )  THEN
    174175          pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
    175176          IF ( .NOT. constant_diffusion     )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
    176           IF ( moisture .OR. passive_scalar )  q_p(:,:,nxl-1) = q_p(:,:,nxl)     
     177          IF ( humidity .OR. passive_scalar )  q_p(:,:,nxl-1) = q_p(:,:,nxl)
    177178       ELSEIF ( outflow_r )  THEN
    178179          pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
    179180          IF ( .NOT. constant_diffusion     )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
    180           IF ( moisture .OR. passive_scalar )  q_p(:,:,nxr+1) = q_p(:,:,nxr)
     181          IF ( humidity .OR. passive_scalar )  q_p(:,:,nxr+1) = q_p(:,:,nxr)
    181182       ENDIF
    182183
    183184    ENDIF
    184185
    185     IF ( range == 'outflow_uvw' ) THEN
    186 !
    187 !--    Radiation boundary condition for the velocities at the respective outflow
    188        IF ( outflow_s ) THEN
    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
     186!
     187!-- Radiation boundary condition for the velocities at the respective outflow
     188    IF ( outflow_s  .AND.                                                 &
     189         intermediate_timestep_count == intermediate_timestep_count_max ) &
     190    THEN
     191
     192       c_max = dy / dt_3d
     193
     194       DO i = nxl-1, nxr+1
     195          DO k = nzb+1, nzt+1
     196
     197!
     198!--          First calculate the phase speeds for u,v, and w
     199             denom = u_m_s(k,0,i) - u_m_s(k,1,i)
     200
     201             IF ( denom /= 0.0 )  THEN
     202                c_u = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / denom
     203                IF ( c_u < 0.0 )  THEN
     204                   c_u = 0.0
     205                ELSEIF ( c_u > c_max )  THEN
    227206                   c_u = c_max
    228207                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
     208             ELSE
     209                c_u = c_max
     210             ENDIF
     211             denom = v_m_s(k,0,i) - v_m_s(k,1,i)
     212
     213             IF ( denom /= 0.0 )  THEN
     214                c_v = -c_max * ( v(k,0,i) - v_m_s(k,0,i) ) / denom
     215                IF ( c_v < 0.0 )  THEN
     216                   c_v = 0.0
     217                ELSEIF ( c_v > c_max )  THEN
    240218                   c_v = c_max
    241219                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
     220             ELSE
     221                c_v = c_max
     222             ENDIF
     223
     224             denom = w_m_s(k,0,i) - w_m_s(k,1,i)
     225
     226             IF ( denom /= 0.0 )  THEN
     227                c_w = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / denom
     228                IF ( c_w < 0.0 )  THEN
     229                   c_w = 0.0
     230                ELSEIF ( c_w > c_max )  THEN
    253231                   c_w = c_max
    254232                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 
    301        ENDIF
    302 
    303        IF ( outflow_l ) THEN
    304 !          u(:,:,nxl-1) = u(:,:,nxl)
    305 !          w(:,:,nxl-1) = 0.0
    306 !
    307 !--       Compute the mean horizontal wind parallel to and within the outflow
    308 !--       wall and use this as boundary condition for v
    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
     233             ELSE
     234                c_w = c_max
     235             ENDIF
     236
     237!
     238!--          Calculate the new velocities
     239             u_p(k,-1,i) = u(k,-1,i) + dt_3d * c_u * &
     240                                       ( u(k,-1,i) - u(k,0,i) ) * ddy
     241
     242             v_p(k,-1,i) = v(k,-1,i) + dt_3d * c_v * &
     243                                       ( v(k,-1,i) - v_m_s(k,0,i) ) * ddy
     244
     245             w_p(k,-1,i) = w(k,-1,i) + dt_3d * c_w * &
     246                                       ( w(k,-1,i) - w(k,0,i) ) * ddy
     247
     248!
     249!--          Save old timelevels for the next timestep
     250             u_m_s(k,:,i) = u(k,-1:1,i)
     251             v_m_s(k,:,i) = v(k,-1:1,i)
     252             w_m_s(k,:,i) = w(k,-1:1,i)
     253
     254          ENDDO
     255       ENDDO
     256
     257!
     258!--    Bottom boundary at the outflow
     259       IF ( ibc_uv_b == 0 )  THEN
     260          u_p(nzb,-1,:) = -u_p(nzb+1,-1,:)
     261          v_p(nzb,-1,:) = -v_p(nzb+1,-1,:) 
     262       ELSE                   
     263          u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
     264          v_p(nzb,-1,:) =  v_p(nzb+1,-1,:)
     265       ENDIF
     266       w_p(nzb,ny+1,:) = 0.0
     267
     268!
     269!--    Top boundary at the outflow
     270       IF ( ibc_uv_t == 0 )  THEN
     271          u_p(nzt+1,-1,:) = ug(nzt+1)
     272          v_p(nzt+1,-1,:) = vg(nzt+1)
     273       ELSE
     274          u_p(nzt+1,-1,:) = u(nzt,-1,:)
     275          v_p(nzt+1,-1,:) = v(nzt,-1,:)
     276       ENDIF
     277       w_p(nzt:nzt+1,-1,:) = 0.0
     278
     279    ENDIF
     280
     281    IF ( outflow_n  .AND.                                                 &
     282         intermediate_timestep_count == intermediate_timestep_count_max ) &
     283    THEN
     284
     285       c_max = dy / dt_3d
     286
     287       DO i = nxl-1, nxr+1
     288          DO k = nzb+1, nzt+1
     289
     290!
     291!--          First calculate the phase speeds for u,v, and w
     292             denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
     293
     294             IF ( denom /= 0.0 )  THEN
     295                c_u = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / denom
     296                IF ( c_u < 0.0 )  THEN
     297                   c_u = 0.0
     298                ELSEIF ( c_u > c_max )  THEN
    343299                   c_u = c_max
    344300                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
     301             ELSE
     302                c_u = c_max
     303             ENDIF
     304
     305             denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
     306
     307             IF ( denom /= 0.0 )  THEN
     308                c_v = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / denom
     309                IF ( c_v < 0.0 )  THEN
     310                   c_v = 0.0
     311                ELSEIF ( c_v > c_max )  THEN
    356312                   c_v = c_max
    357313                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
     314             ELSE
     315                c_v = c_max
     316             ENDIF
     317
     318             denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
     319
     320             IF ( denom /= 0.0 )  THEN
     321                c_w = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / denom
     322                IF ( c_w < 0.0 )  THEN
     323                   c_w = 0.0
     324                ELSEIF ( c_w > c_max )  THEN
    369325                   c_w = c_max
    370326                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 
    417        ENDIF
     327             ELSE
     328                c_w = c_max
     329             ENDIF
     330
     331!
     332!--          Calculate the new velocities
     333             u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * c_u * &
     334                                           ( u(k,ny+1,i) - u(k,ny,i) ) * ddy
     335
     336             v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * c_v * &
     337                                           ( v(k,ny+1,i) - v(k,ny,i) ) * ddy
     338
     339             w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * c_w * &
     340                                           ( w(k,ny+1,i) - w(k,ny,i) ) * ddy
     341
     342!
     343!--          Swap timelevels for the next timestep
     344             u_m_n(k,:,i) = u(k,ny-1:ny+1,i)
     345             v_m_n(k,:,i) = v(k,ny-1:ny+1,i)
     346             w_m_n(k,:,i) = w(k,ny-1:ny+1,i)
     347
     348          ENDDO
     349       ENDDO
     350
     351!
     352!--    Bottom boundary at the outflow
     353       IF ( ibc_uv_b == 0 )  THEN
     354          u_p(nzb,ny+1,:) = -u_p(nzb+1,ny+1,:)
     355          v_p(nzb,ny+1,:) = -v_p(nzb+1,ny+1,:) 
     356       ELSE                   
     357          u_p(nzb,ny+1,:) =  u_p(nzb+1,ny+1,:)
     358          v_p(nzb,ny+1,:) =  v_p(nzb+1,ny+1,:)
     359       ENDIF
     360       w_p(nzb,ny+1,:) = 0.0
     361
     362!
     363!--    Top boundary at the outflow
     364       IF ( ibc_uv_t == 0 )  THEN
     365          u_p(nzt+1,ny+1,:) = ug(nzt+1)
     366          v_p(nzt+1,ny+1,:) = vg(nzt+1)
     367       ELSE
     368          u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
     369          v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
     370       ENDIF
     371       w_p(nzt:nzt+1,ny+1,:) = 0.0
     372
     373    ENDIF
     374
     375    IF ( outflow_l  .AND.                                                 &
     376         intermediate_timestep_count == intermediate_timestep_count_max ) &
     377    THEN
     378
     379       c_max = dx / dt_3d
     380
     381       DO j = nys-1, nyn+1
     382          DO k = nzb+1, nzt+1
     383
     384!
     385!--          First calculate the phase speeds for u,v, and w
     386             denom = u_m_l(k,j,0) - u_m_l(k,j,1)
     387
     388             IF ( denom /= 0.0 )  THEN
     389                c_u = -c_max * ( u(k,j,0) - u_m_r(k,j,0) ) / denom
     390                IF ( c_u > 0.0 )  THEN
     391                   c_u = 0.0
     392                ELSEIF ( c_u < -c_max )  THEN
     393                   c_u = -c_max
     394                ENDIF
     395             ELSE
     396                c_u = -c_max
     397             ENDIF
     398
     399             denom = v_m_l(k,j,0) - v_m_l(k,j,1)
     400
     401             IF ( denom /= 0.0 )  THEN
     402                c_v = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / denom
     403                IF ( c_v < 0.0 )  THEN
     404                   c_v = 0.0
     405                ELSEIF ( c_v > c_max )  THEN
     406                   c_v = c_max
     407                ENDIF
     408             ELSE
     409                c_v = c_max
     410             ENDIF
     411
     412             denom = w_m_l(k,j,0) - w_m_l(k,j,1)
     413
     414             IF ( denom /= 0.0 )  THEN
     415                c_w = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / denom
     416                IF ( c_w < 0.0 )  THEN
     417                   c_w = 0.0
     418                ELSEIF ( c_w > c_max )  THEN
     419                   c_w = c_max
     420                ENDIF
     421             ELSE
     422                c_w = c_max
     423             ENDIF
     424
     425!
     426!--          Calculate the new velocities
     427             u_p(k,j,-1) = u(k,j,-1) + dt_3d * c_u * &
     428                                       ( u(k,j,-1) - u(k,j,0) ) * ddx
     429
     430             v_p(k,j,-1) = v(k,j,-1) + dt_3d * c_v * &
     431                                       ( v(k,j,-1) - v(k,j,0) ) * ddx
     432
     433             w_p(k,j,-1) = w(k,j,-1) + dt_3d * c_w * &
     434                                       ( w(k,j,-1) - w(k,j,0) ) * ddx
     435
     436!
     437!--          Swap timelevels for the next timestep
     438             u_m_l(k,j,:) = u(k,j,-1:1)
     439             v_m_l(k,j,:) = v(k,j,-1:1)
     440             w_m_l(k,j,:) = w(k,j,-1:1)
     441
     442          ENDDO
     443       ENDDO
     444
     445!
     446!--    Bottom boundary at the outflow
     447       IF ( ibc_uv_b == 0 )  THEN
     448          u_p(nzb,:,-1) = -u_p(nzb+1,:,-1)
     449          v_p(nzb,:,-1) = -v_p(nzb+1,:,-1) 
     450       ELSE                   
     451          u_p(nzb,:,-1) =  u_p(nzb+1,:,-1)
     452          v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
     453       ENDIF
     454       w_p(nzb,:,-1) = 0.0
     455
     456!
     457!--    Top boundary at the outflow
     458       IF ( ibc_uv_t == 0 )  THEN
     459          u_p(nzt+1,:,-1) = ug(nzt+1)
     460          v_p(nzt+1,:,-1) = vg(nzt+1)
     461       ELSE
     462          u_p(nzt+1,:,-1) = u_p(nzt,:,-1)
     463          v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
     464       ENDIF
     465       w_p(nzt:nzt+1,:,-1) = 0.0
     466
     467    ENDIF
     468
     469    IF ( outflow_r  .AND.                                                 &
     470         intermediate_timestep_count == intermediate_timestep_count_max ) &
     471    THEN
     472
     473       c_max = dx / dt_3d
     474
     475       DO j = nys-1, nyn+1
     476          DO k = nzb+1, nzt+1
     477
     478!
     479!--          First calculate the phase speeds for u,v, and w
     480             denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
     481
     482             IF ( denom /= 0.0 )  THEN
     483                c_u = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / denom
     484                IF ( c_u < 0.0 )  THEN
     485                   c_u = 0.0
     486                ELSEIF ( c_u > c_max )  THEN
     487                   c_u = c_max
     488                ENDIF
     489             ELSE
     490                c_u = c_max
     491             ENDIF
     492
     493             denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
     494
     495             IF ( denom /= 0.0 )  THEN
     496                c_v = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / denom
     497                IF ( c_v < 0.0 )  THEN
     498                   c_v = 0.0
     499                ELSEIF ( c_v > c_max )  THEN
     500                   c_v = c_max
     501                ENDIF
     502             ELSE
     503                c_v = c_max
     504             ENDIF
     505
     506             denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
     507
     508             IF ( denom /= 0.0 )  THEN
     509                c_w = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / denom
     510                IF ( c_w < 0.0 )  THEN
     511                   c_w = 0.0
     512                ELSEIF ( c_w > c_max )  THEN
     513                   c_w = c_max
     514                ENDIF
     515             ELSE
     516                c_w = c_max
     517             ENDIF
     518
     519!
     520!--          Calculate the new velocities
     521             u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * c_u * &
     522                                           ( u(k,j,nx+1) - u(k,j,nx) ) * ddx
     523
     524             v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * c_v * &
     525                                           ( v(k,j,nx+1) - v(k,j,nx) ) * ddx
     526
     527             w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * c_w * &
     528                                           ( w(k,j,nx+1) - w(k,j,nx) ) * ddx
     529
     530!
     531!--          Swap timelevels for the next timestep
     532             u_m_r(k,j,:) = u(k,j,nx-1:nx+1)
     533             v_m_r(k,j,:) = v(k,j,nx-1:nx+1)
     534             w_m_r(k,j,:) = w(k,j,nx-1:nx+1)
     535
     536          ENDDO
     537       ENDDO
     538
     539!
     540!--    Bottom boundary at the outflow
     541       IF ( ibc_uv_b == 0 )  THEN
     542          u_p(nzb,:,nx+1) = -u_p(nzb+1,:,nx+1)
     543          v_p(nzb,:,nx+1) = -v_p(nzb+1,:,nx+1) 
     544       ELSE                   
     545          u_p(nzb,:,nx+1) =  u_p(nzb+1,:,nx+1)
     546          v_p(nzb,:,nx+1) =  v_p(nzb+1,:,nx+1)
     547       ENDIF
     548       w_p(nzb,:,nx+1) = 0.0
     549
     550!
     551!--    Top boundary at the outflow
     552       IF ( ibc_uv_t == 0 )  THEN
     553          u_p(nzt+1,:,nx+1) = ug(nzt+1)
     554          v_p(nzt+1,:,nx+1) = vg(nzt+1)
     555       ELSE
     556          u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
     557          v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
     558       ENDIF
     559       w(nzt:nzt+1,:,nx+1) = 0.0
    418560
    419561    ENDIF
Note: See TracChangeset for help on using the changeset viewer.