Changeset 978


Ignore:
Timestamp:
Aug 9, 2012 8:28:32 AM (12 years ago)
Author:
fricke
Message:

merge fricke branch back into trunk

Location:
palm/trunk
Files:
27 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk

  • palm/trunk/EXAMPLES/dvr

  • palm/trunk/SOURCE

  • palm/trunk/SOURCE/average_3d_data.f90

    r772 r978  
    44! Current revisions:
    55! -----------------
     6! +z0h_av
    67!
    78! Former revisions:
     
    284285             ENDDO
    285286
     287          CASE ( 'z0h*' )
     288             DO  i = nxlg, nxrg
     289                DO  j = nysg, nyng
     290                   z0h_av(j,i) = z0h_av(j,i) / REAL( average_count_3d )
     291                ENDDO
     292             ENDDO
     293
    286294          CASE DEFAULT
    287295!
  • palm/trunk/SOURCE/boundary_conds.f90

    r876 r978  
    44! Current revisions:
    55! -----------------
    6 !
     6! Neumann boudnary conditions are added at the inflow boundary for the SGS-TKE.
     7! Outflow boundary conditions for the velocity components can be set to Neumann
     8! conditions or to radiation conditions with a horizontal averaged phase
     9! velocity.
    710!
    811
     
    190193          q_p(nzt+1,:,:) = q_p(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
    191194
    192 
    193195       ENDIF
    194196
     
    198200!--    Since in prognostic_equations (cache optimized version) these levels are
    199201!--    handled as a prognostic level, boundary values have to be restored here.
     202!--    For the SGS-TKE, Neumann boundary conditions are used at the inflow.
    200203       IF ( inflow_s )  THEN
    201204          v_p(:,nys,:) = v_p(:,nys-1,:)
     205          IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
     206       ELSEIF ( inflow_n )  THEN
     207          IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
    202208       ELSEIF ( inflow_l ) THEN
    203209          u_p(:,:,nxl) = u_p(:,:,nxl-1)
     210          IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
     211       ELSEIF ( inflow_r )  THEN
     212          IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
    204213       ENDIF
    205214
     
    227236
    228237!
    229 !-- Radiation boundary condition for the velocities at the respective outflow
     238!-- Neumann or Radiation boundary condition for the velocities at the
     239!-- respective outflow
    230240    IF ( outflow_s )  THEN
    231241
    232        c_max = dy / dt_3d
    233 
    234        DO i = nxlg, nxrg
     242       IF ( bc_ns_dirneu )  THEN
     243          u(:,-1,:) = u(:,0,:)
     244          v(:,0,:)  = v(:,1,:)
     245          w(:,-1,:) = w(:,0,:)         
     246       ELSEIF ( bc_ns_dirrad )  THEN
     247
     248          c_max = dy / dt_3d
     249
     250          c_u_m_l = 0.0
     251          c_v_m_l = 0.0
     252          c_w_m_l = 0.0
     253
     254          c_u_m = 0.0
     255          c_v_m = 0.0
     256          c_w_m = 0.0
     257
     258!
     259!--       Calculate the phase speeds for u,v, and w, first local and then
     260!--       average parallel along the outflow boundary. 
    235261          DO k = nzb+1, nzt+1
    236 
    237 !
    238 !--          Calculate the phase speeds for u,v, and w. In case of using a
    239 !--          Runge-Kutta scheme, do this for the first substep only and then
    240 !--          reuse this values for the further substeps.
    241              IF ( intermediate_timestep_count == 1 )  THEN
     262             DO i = nxl, nxr
    242263
    243264                denom = u_m_s(k,0,i) - u_m_s(k,1,i)
    244265
    245266                IF ( denom /= 0.0 )  THEN
    246                    c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / denom
     267                   c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) )             &
     268                              / ( denom * tsc(2) )
    247269                   IF ( c_u(k,i) < 0.0 )  THEN
    248270                      c_u(k,i) = 0.0
     
    257279
    258280                IF ( denom /= 0.0 )  THEN
    259                    c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / denom
     281                   c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) )             &
     282                              / ( denom * tsc(2) )
    260283                   IF ( c_v(k,i) < 0.0 )  THEN
    261284                      c_v(k,i) = 0.0
     
    270293
    271294                IF ( denom /= 0.0 )  THEN
    272                    c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / denom
     295                   c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) )             &
     296                              / ( denom * tsc(2) )
    273297                   IF ( c_w(k,i) < 0.0 )  THEN
    274298                      c_w(k,i) = 0.0
     
    280304                ENDIF
    281305
    282 !
    283 !--             Save old timelevels for the next timestep
    284                 u_m_s(k,:,i) = u(k,0:1,i)
    285                 v_m_s(k,:,i) = v(k,1:2,i)
    286                 w_m_s(k,:,i) = w(k,0:1,i)
    287 
    288              ENDIF
    289 
    290 !
    291 !--          Calculate the new velocities
    292              u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u(k,i) * &
     306                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
     307                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
     308                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
     309
     310             ENDDO
     311          ENDDO
     312
     313#if defined( __parallel )   
     314          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
     315          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
     316                              MPI_SUM, comm1dx, ierr )   
     317          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
     318          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
     319                              MPI_SUM, comm1dx, ierr ) 
     320          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
     321          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
     322                              MPI_SUM, comm1dx, ierr ) 
     323#else
     324          c_u_m = c_u_m_l
     325          c_v_m = c_v_m_l
     326          c_w_m = c_w_m_l
     327#endif
     328
     329          c_u_m = c_u_m / (nx+1)
     330          c_v_m = c_v_m / (nx+1)
     331          c_w_m = c_w_m / (nx+1)
     332
     333!
     334!--       Save old timelevels for the next timestep
     335          IF ( intermediate_timestep_count == 1 )  THEN
     336             u_m_s(:,:,:) = u(:,0:1,:)
     337             v_m_s(:,:,:) = v(:,1:2,:)
     338             w_m_s(:,:,:) = w(:,0:1,:)
     339          ENDIF
     340
     341!
     342!--       Calculate the new velocities
     343          DO k = nzb+1, nzt+1
     344             DO i = nxlg, nxrg
     345                u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) *          &
    293346                                       ( u(k,-1,i) - u(k,0,i) ) * ddy
    294347
    295              v_p(k,0,i)  = v(k,0,i)  - dt_3d * tsc(2) * c_v(k,i) * &
     348                v_p(k,0,i)  = v(k,0,i)  - dt_3d * tsc(2) * c_v_m(k) *          &
    296349                                       ( v(k,0,i) - v(k,1,i) ) * ddy
    297350
    298              w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w(k,i) * &
     351                w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) *          &
    299352                                       ( w(k,-1,i) - w(k,0,i) ) * ddy
    300 
    301           ENDDO
    302        ENDDO
    303 
    304 !
    305 !--    Bottom boundary at the outflow
    306        IF ( ibc_uv_b == 0 )  THEN
    307           u_p(nzb,-1,:) = 0.0
    308           v_p(nzb,0,:)  = 0.0 
    309        ELSE                   
    310           u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
    311           v_p(nzb,0,:)  =  v_p(nzb+1,0,:)
    312        ENDIF
    313        w_p(nzb,-1,:) = 0.0
    314 
    315 !
    316 !--    Top boundary at the outflow
    317        IF ( ibc_uv_t == 0 )  THEN
    318           u_p(nzt+1,-1,:) = u_init(nzt+1)
    319           v_p(nzt+1,0,:)  = v_init(nzt+1)
    320        ELSE
    321           u_p(nzt+1,-1,:) = u(nzt,-1,:)
    322           v_p(nzt+1,0,:)  = v(nzt,0,:)
    323        ENDIF
    324        w_p(nzt:nzt+1,-1,:) = 0.0
     353             ENDDO
     354          ENDDO
     355
     356!
     357!--       Bottom boundary at the outflow
     358          IF ( ibc_uv_b == 0 )  THEN
     359             u_p(nzb,-1,:) = 0.0
     360             v_p(nzb,0,:)  = 0.0 
     361          ELSE                   
     362             u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
     363             v_p(nzb,0,:)  =  v_p(nzb+1,0,:)
     364          ENDIF
     365          w_p(nzb,-1,:) = 0.0
     366
     367!
     368!--       Top boundary at the outflow
     369          IF ( ibc_uv_t == 0 )  THEN
     370             u_p(nzt+1,-1,:) = u_init(nzt+1)
     371             v_p(nzt+1,0,:)  = v_init(nzt+1)
     372          ELSE
     373             u_p(nzt+1,-1,:) = u(nzt,-1,:)
     374             v_p(nzt+1,0,:)  = v(nzt,0,:)
     375          ENDIF
     376          w_p(nzt:nzt+1,-1,:) = 0.0
     377
     378       ENDIF
    325379
    326380    ENDIF
     
    328382    IF ( outflow_n )  THEN
    329383
    330        c_max = dy / dt_3d
    331 
    332        DO i = nxlg, nxrg
     384       IF ( bc_ns_neudir )  THEN
     385          u(:,ny+1,:) = u(:,ny,:)
     386          v(:,ny+1,:) = v(:,ny,:)
     387          w(:,ny+1,:) = w(:,ny,:)         
     388       ELSEIF ( bc_ns_dirrad )  THEN
     389
     390          c_max = dy / dt_3d
     391
     392          c_u_m_l = 0.0
     393          c_v_m_l = 0.0
     394          c_w_m_l = 0.0
     395
     396          c_u_m = 0.0
     397          c_v_m = 0.0
     398          c_w_m = 0.0
     399
     400!
     401!--       Calculate the phase speeds for u,v, and w, first local and then
     402!--       average parallel along the outflow boundary. 
    333403          DO k = nzb+1, nzt+1
    334 
    335 !
    336 !--          Calculate the phase speeds for u,v, and w. In case of using a
    337 !--          Runge-Kutta scheme, do this for the first substep only and then
    338 !--          reuse this values for the further substeps.
    339              IF ( intermediate_timestep_count == 1 )  THEN
     404             DO i = nxl, nxr
    340405
    341406                denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
    342407
    343408                IF ( denom /= 0.0 )  THEN
    344                    c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / denom
     409                   c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) )           &
     410                              / ( denom * tsc(2) )
    345411                   IF ( c_u(k,i) < 0.0 )  THEN
    346412                      c_u(k,i) = 0.0
     
    355421
    356422                IF ( denom /= 0.0 )  THEN
    357                    c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / denom
     423                   c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) )           &
     424                              / ( denom * tsc(2) )
    358425                   IF ( c_v(k,i) < 0.0 )  THEN
    359426                      c_v(k,i) = 0.0
     
    368435
    369436                IF ( denom /= 0.0 )  THEN
    370                    c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / denom
     437                   c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) )           &
     438                              / ( denom * tsc(2) )
    371439                   IF ( c_w(k,i) < 0.0 )  THEN
    372440                      c_w(k,i) = 0.0
     
    378446                ENDIF
    379447
    380 !
    381 !--             Swap timelevels for the next timestep
    382                 u_m_n(k,:,i) = u(k,ny-1:ny,i)
    383                 v_m_n(k,:,i) = v(k,ny-1:ny,i)
    384                 w_m_n(k,:,i) = w(k,ny-1:ny,i)
    385 
    386              ENDIF
    387 
    388 !
    389 !--          Calculate the new velocities
    390              u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u(k,i) * &
    391                                            ( u(k,ny+1,i) - u(k,ny,i) ) * ddy
    392 
    393              v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v(k,i) * &
    394                                            ( v(k,ny+1,i) - v(k,ny,i) ) * ddy
    395 
    396              w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w(k,i) * &
    397                                            ( w(k,ny+1,i) - w(k,ny,i) ) * ddy
    398 
    399           ENDDO
    400        ENDDO
    401 
    402 !
    403 !--    Bottom boundary at the outflow
    404        IF ( ibc_uv_b == 0 )  THEN
    405           u_p(nzb,ny+1,:) = 0.0
    406           v_p(nzb,ny+1,:) = 0.0   
    407        ELSE                   
    408           u_p(nzb,ny+1,:) =  u_p(nzb+1,ny+1,:)
    409           v_p(nzb,ny+1,:) =  v_p(nzb+1,ny+1,:)
    410        ENDIF
    411        w_p(nzb,ny+1,:) = 0.0
    412 
    413 !
    414 !--    Top boundary at the outflow
    415        IF ( ibc_uv_t == 0 )  THEN
    416           u_p(nzt+1,ny+1,:) = u_init(nzt+1)
    417           v_p(nzt+1,ny+1,:) = v_init(nzt+1)
    418        ELSE
    419           u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
    420           v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
    421        ENDIF
    422        w_p(nzt:nzt+1,ny+1,:) = 0.0
     448                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
     449                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
     450                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
     451
     452             ENDDO
     453          ENDDO
     454
     455#if defined( __parallel )   
     456          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
     457          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
     458                              MPI_SUM, comm1dx, ierr )   
     459          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
     460          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
     461                              MPI_SUM, comm1dx, ierr ) 
     462          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
     463          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
     464                              MPI_SUM, comm1dx, ierr ) 
     465#else
     466          c_u_m = c_u_m_l
     467          c_v_m = c_v_m_l
     468          c_w_m = c_w_m_l
     469#endif
     470
     471          c_u_m = c_u_m / (nx+1)
     472          c_v_m = c_v_m / (nx+1)
     473          c_w_m = c_w_m / (nx+1)
     474
     475!
     476!--       Save old timelevels for the next timestep
     477          IF ( intermediate_timestep_count == 1 )  THEN
     478                u_m_n(:,:,:) = u(:,ny-1:ny,:)
     479                v_m_n(:,:,:) = v(:,ny-1:ny,:)
     480                w_m_n(:,:,:) = w(:,ny-1:ny,:)
     481          ENDIF
     482
     483!
     484!--       Calculate the new velocities
     485          DO k = nzb+1, nzt+1
     486             DO i = nxlg, nxrg
     487                u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) *      &
     488                                       ( u(k,ny+1,i) - u(k,ny,i) ) * ddy
     489
     490                v_p(k,ny+1,i) = v(k,ny+1,i)  - dt_3d * tsc(2) * c_v_m(k) *     &
     491                                       ( v(k,ny+1,i) - v(k,ny,i) ) * ddy
     492
     493                w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) *      &
     494                                       ( w(k,ny+1,i) - w(k,ny,i) ) * ddy
     495             ENDDO
     496          ENDDO
     497
     498!
     499!--       Bottom boundary at the outflow
     500          IF ( ibc_uv_b == 0 )  THEN
     501             u_p(nzb,ny+1,:) = 0.0
     502             v_p(nzb,ny+1,:) = 0.0   
     503          ELSE                   
     504             u_p(nzb,ny+1,:) =  u_p(nzb+1,ny+1,:)
     505             v_p(nzb,ny+1,:) =  v_p(nzb+1,ny+1,:)
     506          ENDIF
     507          w_p(nzb,ny+1,:) = 0.0
     508
     509!
     510!--       Top boundary at the outflow
     511          IF ( ibc_uv_t == 0 )  THEN
     512             u_p(nzt+1,ny+1,:) = u_init(nzt+1)
     513             v_p(nzt+1,ny+1,:) = v_init(nzt+1)
     514          ELSE
     515             u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
     516             v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
     517          ENDIF
     518          w_p(nzt:nzt+1,ny+1,:) = 0.0
     519
     520       ENDIF
    423521
    424522    ENDIF
     
    426524    IF ( outflow_l )  THEN
    427525
    428        c_max = dx / dt_3d
    429 
    430        DO j = nysg, nyng
     526       IF ( bc_lr_neudir )  THEN
     527          u(:,:,-1) = u(:,:,0)
     528          v(:,:,0)  = v(:,:,1)
     529          w(:,:,-1) = w(:,:,0)         
     530       ELSEIF ( bc_ns_dirrad )  THEN
     531
     532          c_max = dx / dt_3d
     533
     534          c_u_m_l = 0.0
     535          c_v_m_l = 0.0
     536          c_w_m_l = 0.0
     537
     538          c_u_m = 0.0
     539          c_v_m = 0.0
     540          c_w_m = 0.0
     541
     542!
     543!--       Calculate the phase speeds for u,v, and w, first local and then
     544!--       average parallel along the outflow boundary. 
    431545          DO k = nzb+1, nzt+1
    432 
    433 !
    434 !--          Calculate the phase speeds for u,v, and w. In case of using a
    435 !--          Runge-Kutta scheme, do this for the first substep only and then
    436 !--          reuse this values for the further substeps.
    437              IF ( intermediate_timestep_count == 1 )  THEN
     546             DO j = nys, nyn
    438547
    439548                denom = u_m_l(k,j,1) - u_m_l(k,j,2)
    440549
    441550                IF ( denom /= 0.0 )  THEN
    442                    c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / denom
     551                   c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) )             &
     552                              / ( denom * tsc(2) )
    443553                   IF ( c_u(k,j) < 0.0 )  THEN
    444554                      c_u(k,j) = 0.0
     
    453563
    454564                IF ( denom /= 0.0 )  THEN
    455                    c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / denom
     565                   c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) )             &
     566                              / ( denom * tsc(2) )
    456567                   IF ( c_v(k,j) < 0.0 )  THEN
    457568                      c_v(k,j) = 0.0
     
    466577
    467578                IF ( denom /= 0.0 )  THEN
    468                    c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / denom
     579                   c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) )             &
     580                              / ( denom * tsc(2) )
    469581                   IF ( c_w(k,j) < 0.0 )  THEN
    470582                      c_w(k,j) = 0.0
     
    476588                ENDIF
    477589
    478 !
    479 !--             Swap timelevels for the next timestep
    480                 u_m_l(k,j,:) = u(k,j,1:2)
    481                 v_m_l(k,j,:) = v(k,j,0:1)
    482                 w_m_l(k,j,:) = w(k,j,0:1)
    483 
    484              ENDIF
    485 
    486 !
    487 !--          Calculate the new velocities
    488              u_p(k,j,0)  = u(k,j,0)  - dt_3d * tsc(2) * c_u(k,j) * &
     590                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
     591                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
     592                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
     593
     594             ENDDO
     595          ENDDO
     596
     597#if defined( __parallel )   
     598          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
     599          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
     600                              MPI_SUM, comm1dy, ierr )   
     601          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
     602          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
     603                              MPI_SUM, comm1dy, ierr ) 
     604          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
     605          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
     606                              MPI_SUM, comm1dy, ierr ) 
     607#else
     608          c_u_m = c_u_m_l
     609          c_v_m = c_v_m_l
     610          c_w_m = c_w_m_l
     611#endif
     612
     613          c_u_m = c_u_m / (ny+1)
     614          c_v_m = c_v_m / (ny+1)
     615          c_w_m = c_w_m / (ny+1)
     616
     617!
     618!--       Save old timelevels for the next timestep
     619          IF ( intermediate_timestep_count == 1 )  THEN
     620                u_m_l(:,:,:) = u(:,:,1:2)
     621                v_m_l(:,:,:) = v(:,:,0:1)
     622                w_m_l(:,:,:) = w(:,:,0:1)
     623          ENDIF
     624
     625!
     626!--       Calculate the new velocities
     627          DO k = nzb+1, nzt+1
     628             DO i = nxlg, nxrg
     629                u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) *            &
    489630                                       ( u(k,j,0) - u(k,j,1) ) * ddx
    490631
    491              v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v(k,j) * &
     632                v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) *          &
    492633                                       ( v(k,j,-1) - v(k,j,0) ) * ddx
    493634
    494              w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w(k,j) * &
     635                w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) *          &
    495636                                       ( w(k,j,-1) - w(k,j,0) ) * ddx
    496 
    497           ENDDO
    498        ENDDO
    499 
    500 !
    501 !--    Bottom boundary at the outflow
    502        IF ( ibc_uv_b == 0 )  THEN
    503           u_p(nzb,:,0)  = 0.0
    504           v_p(nzb,:,-1) = 0.0
    505        ELSE                   
    506           u_p(nzb,:,0)  =  u_p(nzb+1,:,0)
    507           v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
    508        ENDIF
    509        w_p(nzb,:,-1) = 0.0
    510 
    511 !
    512 !--    Top boundary at the outflow
    513        IF ( ibc_uv_t == 0 )  THEN
    514           u_p(nzt+1,:,-1) = u_init(nzt+1)
    515           v_p(nzt+1,:,-1) = v_init(nzt+1)
    516        ELSE
    517           u_p(nzt+1,:,-1) = u_p(nzt,:,-1)
    518           v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
    519        ENDIF
    520        w_p(nzt:nzt+1,:,-1) = 0.0
     637             ENDDO
     638          ENDDO
     639
     640!
     641!--       Bottom boundary at the outflow
     642          IF ( ibc_uv_b == 0 )  THEN
     643             u_p(nzb,:,0)  = 0.0
     644             v_p(nzb,:,-1) = 0.0
     645          ELSE                   
     646             u_p(nzb,:,0)  =  u_p(nzb+1,:,0)
     647             v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
     648          ENDIF
     649          w_p(nzb,:,-1) = 0.0
     650
     651!
     652!--       Top boundary at the outflow
     653          IF ( ibc_uv_t == 0 )  THEN
     654             u_p(nzt+1,:,-1) = u_init(nzt+1)
     655             v_p(nzt+1,:,-1) = v_init(nzt+1)
     656          ELSE
     657             u_p(nzt+1,:,-1) = u_p(nzt,:,-1)
     658             v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
     659          ENDIF
     660          w_p(nzt:nzt+1,:,-1) = 0.0
     661
     662       ENDIF
    521663
    522664    ENDIF
     
    524666    IF ( outflow_r )  THEN
    525667
    526        c_max = dx / dt_3d
    527 
    528        DO j = nysg, nyng
     668       IF ( bc_lr_dirneu )  THEN
     669          u(:,:,nx+1) = u(:,:,nx)
     670          v(:,:,nx+1) = v(:,:,nx)
     671          w(:,:,nx+1) = w(:,:,nx)         
     672       ELSEIF ( bc_ns_dirrad )  THEN
     673
     674          c_max = dx / dt_3d
     675
     676          c_u_m_l = 0.0
     677          c_v_m_l = 0.0
     678          c_w_m_l = 0.0
     679
     680          c_u_m = 0.0
     681          c_v_m = 0.0
     682          c_w_m = 0.0
     683
     684!
     685!--       Calculate the phase speeds for u,v, and w, first local and then
     686!--       average parallel along the outflow boundary. 
    529687          DO k = nzb+1, nzt+1
    530 
    531 !
    532 !--          Calculate the phase speeds for u,v, and w. In case of using a
    533 !--          Runge-Kutta scheme, do this for the first substep only and then
    534 !--          reuse this values for the further substeps.
    535              IF ( intermediate_timestep_count == 1 )  THEN
     688             DO j = nys, nyn
    536689
    537690                denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
    538691
    539692                IF ( denom /= 0.0 )  THEN
    540                    c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / denom
     693                   c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) )           &
     694                              / ( denom * tsc(2) )
    541695                   IF ( c_u(k,j) < 0.0 )  THEN
    542696                      c_u(k,j) = 0.0
     
    551705
    552706                IF ( denom /= 0.0 )  THEN
    553                    c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / denom
     707                   c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) )           &
     708                              / ( denom * tsc(2) )
    554709                   IF ( c_v(k,j) < 0.0 )  THEN
    555710                      c_v(k,j) = 0.0
     
    564719
    565720                IF ( denom /= 0.0 )  THEN
    566                    c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / denom
     721                   c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) )           &
     722                              / ( denom * tsc(2) )
    567723                   IF ( c_w(k,j) < 0.0 )  THEN
    568724                      c_w(k,j) = 0.0
     
    574730                ENDIF
    575731
    576 !
    577 !--             Swap timelevels for the next timestep
    578                 u_m_r(k,j,:) = u(k,j,nx-1:nx)
    579                 v_m_r(k,j,:) = v(k,j,nx-1:nx)
    580                 w_m_r(k,j,:) = w(k,j,nx-1:nx)
    581 
    582              ENDIF
    583 
    584 !
    585 !--          Calculate the new velocities
    586              u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u(k,j) * &
    587                                            ( u(k,j,nx+1) - u(k,j,nx) ) * ddx
    588 
    589              v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v(k,j) * &
    590                                            ( v(k,j,nx+1) - v(k,j,nx) ) * ddx
    591 
    592              w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w(k,j) * &
    593                                            ( w(k,j,nx+1) - w(k,j,nx) ) * ddx
    594 
    595           ENDDO
    596        ENDDO
    597 
    598 !
    599 !--    Bottom boundary at the outflow
    600        IF ( ibc_uv_b == 0 )  THEN
    601           u_p(nzb,:,nx+1) = 0.0
    602           v_p(nzb,:,nx+1) = 0.0
    603        ELSE                   
    604           u_p(nzb,:,nx+1) =  u_p(nzb+1,:,nx+1)
    605           v_p(nzb,:,nx+1) =  v_p(nzb+1,:,nx+1)
    606        ENDIF
    607        w_p(nzb,:,nx+1) = 0.0
    608 
    609 !
    610 !--    Top boundary at the outflow
    611        IF ( ibc_uv_t == 0 )  THEN
    612           u_p(nzt+1,:,nx+1) = u_init(nzt+1)
    613           v_p(nzt+1,:,nx+1) = v_init(nzt+1)
    614        ELSE
    615           u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
    616           v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
    617        ENDIF
    618        w(nzt:nzt+1,:,nx+1) = 0.0
     732                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
     733                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
     734                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
     735
     736             ENDDO
     737          ENDDO
     738
     739#if defined( __parallel )   
     740          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
     741          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
     742                              MPI_SUM, comm1dy, ierr )   
     743          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
     744          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
     745                              MPI_SUM, comm1dy, ierr ) 
     746          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
     747          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
     748                              MPI_SUM, comm1dy, ierr ) 
     749#else
     750          c_u_m = c_u_m_l
     751          c_v_m = c_v_m_l
     752          c_w_m = c_w_m_l
     753#endif
     754
     755          c_u_m = c_u_m / (ny+1)
     756          c_v_m = c_v_m / (ny+1)
     757          c_w_m = c_w_m / (ny+1)
     758
     759!
     760!--       Save old timelevels for the next timestep
     761          IF ( intermediate_timestep_count == 1 )  THEN
     762                u_m_r(:,:,:) = u(:,:,nx-1:nx)
     763                v_m_r(:,:,:) = v(:,:,nx-1:nx)
     764                w_m_r(:,:,:) = w(:,:,nx-1:nx)
     765          ENDIF
     766
     767!
     768!--       Calculate the new velocities
     769          DO k = nzb+1, nzt+1
     770             DO i = nxlg, nxrg
     771                u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) *      &
     772                                       ( u(k,j,nx+1) - u(k,j,nx) ) * ddx
     773
     774                v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) *      &
     775                                       ( v(k,j,nx+1) - v(k,j,nx) ) * ddx
     776
     777                w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) *      &
     778                                       ( w(k,j,nx+1) - w(k,j,nx) ) * ddx
     779             ENDDO
     780          ENDDO
     781
     782!
     783!--       Bottom boundary at the outflow
     784          IF ( ibc_uv_b == 0 )  THEN
     785             u_p(nzb,:,nx+1) = 0.0
     786             v_p(nzb,:,nx+1) = 0.0
     787          ELSE                   
     788             u_p(nzb,:,nx+1) =  u_p(nzb+1,:,nx+1)
     789             v_p(nzb,:,nx+1) =  v_p(nzb+1,:,nx+1)
     790          ENDIF
     791          w_p(nzb,:,nx+1) = 0.0
     792
     793!
     794!--       Top boundary at the outflow
     795          IF ( ibc_uv_t == 0 )  THEN
     796             u_p(nzt+1,:,nx+1) = u_init(nzt+1)
     797             v_p(nzt+1,:,nx+1) = v_init(nzt+1)
     798          ELSE
     799             u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
     800             v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
     801          ENDIF
     802          w(nzt:nzt+1,:,nx+1) = 0.0
     803
     804       ENDIF
    619805
    620806    ENDIF
    621807
    622  
    623808 END SUBROUTINE boundary_conds
  • palm/trunk/SOURCE/check_parameters.f90

    r965 r978  
    44! Current revisions:
    55! -----------------
    6 !
     6! setting of bc_lr/ns_dirneu/neudir
     7! outflow damping layer removed
     8! check for z0h*
     9! check for pt_damping_width
    710!
    811! Former revisions:
     
    13791382!-- Lateral boundary conditions
    13801383    IF ( bc_lr /= 'cyclic'  .AND.  bc_lr /= 'dirichlet/radiation'  .AND. &
    1381          bc_lr /= 'radiation/dirichlet' )  THEN
     1384         bc_lr /= 'radiation/dirichlet' .AND. bc_lr /= 'dirichlet/neumann' &
     1385         .AND. bc_lr /= 'neumann/dirichlet' )  THEN
    13821386       message_string = 'unknown boundary condition: bc_lr = "' // &
    13831387                        TRIM( bc_lr ) // '"'
     
    13851389    ENDIF
    13861390    IF ( bc_ns /= 'cyclic'  .AND.  bc_ns /= 'dirichlet/radiation'  .AND. &
    1387          bc_ns /= 'radiation/dirichlet' )  THEN
     1391         bc_ns /= 'radiation/dirichlet' .AND. bc_ns /= 'dirichlet/neumann' &
     1392         .AND. bc_ns /= 'neumann/dirichlet' )  THEN
    13881393       message_string = 'unknown boundary condition: bc_ns = "' // &
    13891394                        TRIM( bc_ns ) // '"'
     
    13961401    IF ( bc_lr == 'dirichlet/radiation' )  bc_lr_dirrad = .TRUE.
    13971402    IF ( bc_lr == 'radiation/dirichlet' )  bc_lr_raddir = .TRUE.
     1403    IF ( bc_lr == 'dirichlet/neumann' )    bc_lr_dirneu = .TRUE.
     1404    IF ( bc_lr == 'neumann/dirichlet' )    bc_lr_neudir = .TRUE.
    13981405    IF ( bc_ns /= 'cyclic' )               bc_ns_cyc    = .FALSE.
    13991406    IF ( bc_ns == 'dirichlet/radiation' )  bc_ns_dirrad = .TRUE.
    14001407    IF ( bc_ns == 'radiation/dirichlet' )  bc_ns_raddir = .TRUE.
     1408    IF ( bc_ns == 'dirichlet/neumann' )    bc_ns_dirneu = .TRUE.
     1409    IF ( bc_ns == 'neumann/dirichlet' )    bc_ns_neudir = .TRUE.
    14011410
    14021411!
     
    26592668             unit = 'psu'
    26602669
    2661           CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'qsws*', 'shf*', 'z0*' )
     2670          CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'qsws*', 'shf*', 'z0*', 'z0h*' )
    26622671             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
    26632672                message_string = 'illegal value for data_output: "' // &
     
    27002709             IF ( TRIM( var ) == 'u*'     )  unit = 'm/s'
    27012710             IF ( TRIM( var ) == 'z0*'    )  unit = 'm'
     2711             IF ( TRIM( var ) == 'z0h*'    )  unit = 'm'
    27022712
    27032713
     
    29642974
    29652975!
    2966 !-- In case of non-cyclic lateral boundaries, set the default maximum value
    2967 !-- for the horizontal diffusivity used within the outflow damping layer,
    2968 !-- and check/set the width of the damping layer
     2976!-- In case of non-cyclic lateral boundaries and a damping layer for the
     2977!-- potential temperature, check the width of the damping layer
    29692978    IF ( bc_lr /= 'cyclic' ) THEN
    2970        IF ( km_damp_max == -1.0 )  THEN
    2971           km_damp_max = 0.5 * dx
    2972        ENDIF
    2973        IF ( outflow_damping_width == -1.0 )  THEN
    2974           outflow_damping_width = MIN( 20, nx/2 )
    2975        ENDIF
    2976        IF ( outflow_damping_width <= 0  .OR.  outflow_damping_width > nx )  THEN
    2977           message_string = 'outflow_damping width out of range'
     2979       IF ( pt_damping_width < 0.0 .OR. pt_damping_width > REAL( nx * dx ) )  THEN
     2980          message_string = 'pt_damping_width out of range'
    29782981          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
    29792982       ENDIF
     
    29812984
    29822985    IF ( bc_ns /= 'cyclic' )  THEN
    2983        IF ( km_damp_max == -1.0 )  THEN
    2984           km_damp_max = 0.5 * dy
    2985        ENDIF
    2986        IF ( outflow_damping_width == -1.0 )  THEN
    2987           outflow_damping_width = MIN( 20, ny/2 )
    2988        ENDIF
    2989        IF ( outflow_damping_width <= 0  .OR.  outflow_damping_width > ny )  THEN
    2990           message_string = 'outflow_damping width out of range'
     2986       IF ( pt_damping_width < 0.0 .OR. pt_damping_width > REAL( ny * dy ) )  THEN
     2987          message_string = 'pt_damping_width out of range'
    29912988          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
    29922989       ENDIF
     
    31083105    ENDIF
    31093106
    3110     IF ( bc_lr == 'radiation/dirichlet' )  THEN
     3107    IF ( bc_lr == 'radiation/dirichlet' .OR. bc_lr == 'neumann/dirichlet' )  THEN
    31113108       dist_nxr    = nx - inflow_disturbance_begin
    31123109       dist_nxl(1) = nx - inflow_disturbance_end
    3113     ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
     3110    ELSEIF ( bc_lr == 'dirichlet/radiation' .OR. bc_lr == 'dirichlet/neumann' )  THEN
    31143111       dist_nxl    = inflow_disturbance_begin
    31153112       dist_nxr(1) = inflow_disturbance_end
    31163113    ENDIF
    3117     IF ( bc_ns == 'dirichlet/radiation' )  THEN
     3114    IF ( bc_ns == 'dirichlet/radiation' .OR. bc_ns == 'dirichlet/neumann' )  THEN
    31183115       dist_nyn    = ny - inflow_disturbance_begin
    31193116       dist_nys(1) = ny - inflow_disturbance_end
    3120     ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
     3117    ELSEIF ( bc_ns == 'radiation/dirichlet' .OR. bc_ns == 'neumann/dirichlet' )  THEN
    31213118       dist_nys    = inflow_disturbance_begin
    31223119       dist_nyn(1) = inflow_disturbance_end
     
    31263123!-- A turbulent inflow requires Dirichlet conditions at the respective inflow
    31273124!-- boundary (so far, a turbulent inflow is realized from the left side only)
    3128     IF ( turbulent_inflow  .AND.  bc_lr /= 'dirichlet/radiation' )  THEN
     3125    IF ( turbulent_inflow  .AND.  bc_lr /= 'dirichlet/radiation' .AND.  bc_lr /= 'dirichlet/neumann' )  THEN
    31293126       message_string = 'turbulent_inflow = .T. requires a Dirichlet ' // &
    31303127                        'condition at the inflow boundary'
  • palm/trunk/SOURCE/data_output_2d.f90

    r791 r978  
    44! Current revisions:
    55! -----------------
    6 !
     6! +z0h
    77!
    88! Former revisions:
     
    594594                      DO  j = nysg, nyng
    595595                         local_pf(i,j,nzb+1) =  z0_av(j,i)
     596                      ENDDO
     597                   ENDDO
     598                ENDIF
     599                resorted = .TRUE.
     600                two_d = .TRUE.
     601                level_z(nzb+1) = zu(nzb+1)
     602
     603             CASE ( 'z0h*_xy' )        ! 2d-array
     604                IF ( av == 0 ) THEN
     605                   DO  i = nxlg, nxrg
     606                      DO  j = nysg, nyng
     607                         local_pf(i,j,nzb+1) =  z0h(j,i)
     608                      ENDDO
     609                   ENDDO
     610                ELSE
     611                   DO  i = nxlg, nxrg
     612                      DO  j = nysg, nyng
     613                         local_pf(i,j,nzb+1) =  z0h_av(j,i)
    596614                      ENDDO
    597615                   ENDDO
  • palm/trunk/SOURCE/diffusion_u.f90

    r668 r978  
    44! Current revisions:
    55! -----------------
     6! outflow damping layer removed
     7! kmym_x/_y and kmyp_x/_y change to kmym and kmyp
    68!
    79! Former revisions:
     
    6365! Call for all grid points
    6466!------------------------------------------------------------------------------!
    65     SUBROUTINE diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, uswst, &
    66                             v, w )
     67    SUBROUTINE diffusion_u( ddzu, ddzw, km, tend, u, usws, uswst, v, w )
    6768
    6869       USE control_parameters
     
    7374
    7475       INTEGER ::  i, j, k
    75        REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
    76        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nysg:nyng)
     76       REAL    ::  kmym, kmyp, kmzm, kmzp
     77       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
    7778       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    7879       REAL, DIMENSION(:,:),   POINTER ::  usws, uswst
     
    9596!
    9697!--             Interpolate eddy diffusivities on staggered gridpoints
    97                 kmyp_x = 0.25 * &
    98                          ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
    99                 kmym_x = 0.25 * &
    100                          ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
    101                 kmyp_y = kmyp_x
    102                 kmym_y = kmym_x
    103 !
    104 !--             Increase diffusion at the outflow boundary in case of
    105 !--             non-cyclic lateral boundaries. Damping is only needed for
    106 !--             velocity components parallel to the outflow boundary in
    107 !--             the direction normal to the outflow boundary.
    108                 IF ( .NOT. bc_ns_cyc )  THEN
    109                    kmyp_y = MAX( kmyp_y, km_damp_y(j) )
    110                    kmym_y = MAX( kmym_y, km_damp_y(j) )
    111                 ENDIF
     98                kmyp = 0.25 * &
     99                       ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
     100                kmym = 0.25 * &
     101                       ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
    112102
    113103                tend(k,j,i) = tend(k,j,i)                                    &
     
    116106                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
    117107                      &         ) * ddx2                                     &
    118                       & + ( kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy     &
    119                       &   + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx     &
    120                       &   - kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
    121                       &   - kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
     108                      & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy       &
     109                      &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx       &
     110                      &   - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy           &
     111                      &   - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx           &
    122112                      &   ) * ddy
    123113             ENDDO
     
    128118
    129119                DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
    130                    kmyp_x = 0.25 * &
    131                             ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
    132                    kmym_x = 0.25 * &
    133                             ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
    134                    kmyp_y = kmyp_x
    135                    kmym_y = kmym_x
    136 !
    137 !--                Increase diffusion at the outflow boundary in case of
    138 !--                non-cyclic lateral boundaries. Damping is only needed for
    139 !--                velocity components parallel to the outflow boundary in
    140 !--                the direction normal to the outflow boundary.
    141                    IF ( .NOT. bc_ns_cyc )  THEN
    142                       kmyp_y = MAX( kmyp_y, km_damp_y(j) )
    143                       kmym_y = MAX( kmym_y, km_damp_y(j) )
    144                    ENDIF
     120                   kmyp = 0.25 * &
     121                          ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
     122                   kmym = 0.25 * &
     123                          ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
    145124
    146125                   tend(k,j,i) = tend(k,j,i)                                   &
     
    150129                                         ) * ddx2                              &
    151130                                 + (   fyp(j,i) * (                            &
    152                                   kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy &
    153                                 + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
     131                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy  &
     132                                + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx  &
    154133                                                  )                            &
    155134                                     - fym(j,i) * (                            &
    156                                   kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
    157                                 + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
     135                                  kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
     136                                + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
    158137                                                  )                            &
    159138                                     + wall_u(j,i) * usvs(k,j,i)               &
     
    239218! Call for grid point i,j
    240219!------------------------------------------------------------------------------!
    241     SUBROUTINE diffusion_u_ij( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
     220    SUBROUTINE diffusion_u_ij( i, j, ddzu, ddzw, km, tend, u, usws, &
    242221                               uswst, v, w )
    243222
     
    249228
    250229       INTEGER ::  i, j, k
    251        REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
    252        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nysg:nyng)
     230       REAL    ::  kmym, kmyp, kmzm, kmzp
     231       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
    253232       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    254233       REAL, DIMENSION(nzb:nzt+1)      ::  usvs
     
    261240!
    262241!--       Interpolate eddy diffusivities on staggered gridpoints
    263           kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
    264           kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
    265           kmyp_y = kmyp_x
    266           kmym_y = kmym_x
    267 
    268 !
    269 !--       Increase diffusion at the outflow boundary in case of non-cyclic
    270 !--       lateral boundaries. Damping is only needed for velocity components
    271 !--       parallel to the outflow boundary in the direction normal to the
    272 !--       outflow boundary.
    273           IF ( .NOT. bc_ns_cyc )  THEN
    274              kmyp_y = MAX( kmyp_y, km_damp_y(j) )
    275              kmym_y = MAX( kmym_y, km_damp_y(j) )
    276           ENDIF
     242          kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
     243          kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
    277244
    278245          tend(k,j,i) = tend(k,j,i)                                          &
     
    281248                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
    282249                      &         ) * ddx2                                     &
    283                       & + ( kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy     &
    284                       &   + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx     &
    285                       &   - kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
    286                       &   - kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
     250                      & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy       &
     251                      &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx       &
     252                      &   - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy           &
     253                      &   - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx           &
    287254                      &   ) * ddy
    288255       ENDDO
     
    298265
    299266          DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
    300              kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
    301              kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
    302              kmyp_y = kmyp_x
    303              kmym_y = kmym_x
    304 !
    305 !--          Increase diffusion at the outflow boundary in case of
    306 !--          non-cyclic lateral boundaries. Damping is only needed for
    307 !--          velocity components parallel to the outflow boundary in
    308 !--          the direction normal to the outflow boundary.
    309              IF ( .NOT. bc_ns_cyc )  THEN
    310                 kmyp_y = MAX( kmyp_y, km_damp_y(j) )
    311                 kmym_y = MAX( kmym_y, km_damp_y(j) )
    312              ENDIF
     267             kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
     268             kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
    313269
    314270             tend(k,j,i) = tend(k,j,i)                                         &
     
    318274                                         ) * ddx2                              &
    319275                                 + (   fyp(j,i) * (                            &
    320                                   kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy &
    321                                 + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
     276                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy  &
     277                                + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx  &
    322278                                                  )                            &
    323279                                     - fym(j,i) * (                            &
    324                                   kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
    325                                 + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
     280                                  kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
     281                                + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
    326282                                                  )                            &
    327283                                     + wall_u(j,i) * usvs(k)                   &
  • palm/trunk/SOURCE/diffusion_v.f90

    r668 r978  
    44! Current revisions:
    55! -----------------
     6! outflow damping layer removed
     7! kmxm_x/_y and kmxp_x/_y change to kmxm and kmxp
    68!
    79! Former revisions:
     
    6163! Call for all grid points
    6264!------------------------------------------------------------------------------!
    63     SUBROUTINE diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, &
    64                             vswst, w )
     65    SUBROUTINE diffusion_v( ddzu, ddzw, km, tend, u, v, vsws, vswst, w )
    6566
    6667       USE control_parameters
     
    7172
    7273       INTEGER ::  i, j, k
    73        REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp
    74        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg)
     74       REAL    ::  kmxm, kmxp, kmzm, kmzp
     75       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
    7576       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    7677       REAL, DIMENSION(:,:),   POINTER ::  vsws, vswst
     
    9394!
    9495!--             Interpolate eddy diffusivities on staggered gridpoints
    95                 kmxp_x = 0.25 * &
    96                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
    97                 kmxm_x = 0.25 * &
    98                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    99                 kmxp_y = kmxp_x
    100                 kmxm_y = kmxm_x
    101 !
    102 !--             Increase diffusion at the outflow boundary in case of
    103 !--             non-cyclic lateral boundaries. Damping is only needed for
    104 !--             velocity components parallel to the outflow boundary in
    105 !--             the direction normal to the outflow boundary.
    106                 IF ( .NOT. bc_lr_cyc )  THEN
    107                    kmxp_x = MAX( kmxp_x, km_damp_x(i) )
    108                    kmxm_x = MAX( kmxm_x, km_damp_x(i) )
    109                 ENDIF
     96                kmxp = 0.25 * &
     97                       ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
     98                kmxm = 0.25 * &
     99                       ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    110100
    111101                tend(k,j,i) = tend(k,j,i)                                    &
    112                       & + ( kmxp_x * ( v(k,j,i+1) - v(k,j,i)     ) * ddx     &
    113                       &   + kmxp_y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy     &
    114                       &   - kmxm_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
    115                       &   - kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
     102                      & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx       &
     103                      &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy       &
     104                      &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx           &
     105                      &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy           &
    116106                      &   ) * ddx                                            &
    117107                      & + 2.0 * (                                            &
     
    126116
    127117                DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
    128                    kmxp_x = 0.25 * &
    129                             ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
    130                    kmxm_x = 0.25 * &
    131                             ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    132                    kmxp_y = kmxp_x
    133                    kmxm_y = kmxm_x
    134 !
    135 !--                Increase diffusion at the outflow boundary in case of
    136 !--                non-cyclic lateral boundaries. Damping is only needed for
    137 !--                velocity components parallel to the outflow boundary in
    138 !--                the direction normal to the outflow boundary.
    139                    IF ( .NOT. bc_lr_cyc )  THEN
    140                       kmxp_x = MAX( kmxp_x, km_damp_x(i) )
    141                       kmxm_x = MAX( kmxm_x, km_damp_x(i) )
    142                    ENDIF
    143 
     118                   kmxp = 0.25 * &
     119                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
     120                   kmxm = 0.25 * &
     121                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
     122                   
    144123                   tend(k,j,i) = tend(k,j,i)                                   &
    145124                                 + 2.0 * (                                     &
     
    148127                                         ) * ddy2                              &
    149128                                 + (   fxp(j,i) * (                            &
    150                                   kmxp_x * ( v(k,j,i+1) - v(k,j,i)     ) * ddx &
    151                                 + kmxp_y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
     129                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx  &
     130                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy  &
    152131                                                  )                            &
    153132                                     - fxm(j,i) * (                            &
    154                                   kmxm_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
    155                                 + kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
     133                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
     134                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
    156135                                                  )                            &
    157136                                     + wall_v(j,i) * vsus(k,j,i)               &
     
    237216! Call for grid point i,j
    238217!------------------------------------------------------------------------------!
    239     SUBROUTINE diffusion_v_ij( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
     218    SUBROUTINE diffusion_v_ij( i, j, ddzu, ddzw, km, tend, u, v, &
    240219                               vsws, vswst, w )
    241220
     
    247226
    248227       INTEGER ::  i, j, k
    249        REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp
    250        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg)
     228       REAL    ::  kmxm, kmxp, kmzm, kmzp
     229       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
    251230       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    252231       REAL, DIMENSION(nzb:nzt+1)      ::  vsus
     
    259238!
    260239!--       Interpolate eddy diffusivities on staggered gridpoints
    261           kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
    262           kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    263           kmxp_y = kmxp_x
    264           kmxm_y = kmxm_x
    265 !
    266 !--       Increase diffusion at the outflow boundary in case of non-cyclic
    267 !--       lateral boundaries. Damping is only needed for velocity components
    268 !--       parallel to the outflow boundary in the direction normal to the
    269 !--       outflow boundary.
    270           IF ( .NOT. bc_lr_cyc )  THEN
    271              kmxp_x = MAX( kmxp_x, km_damp_x(i) )
    272              kmxm_x = MAX( kmxm_x, km_damp_x(i) )
    273           ENDIF
     240          kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
     241          kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    274242
    275243          tend(k,j,i) = tend(k,j,i)                                          &
    276                       & + ( kmxp_x * ( v(k,j,i+1) - v(k,j,i)     ) * ddx     &
    277                       &   + kmxp_y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy     &
    278                       &   - kmxm_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
    279                       &   - kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
     244                      & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx       &
     245                      &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy       &
     246                      &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx           &
     247                      &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy           &
    280248                      &   ) * ddx                                            &
    281249                      & + 2.0 * (                                            &
     
    295263
    296264          DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
    297              kmxp_x = 0.25 * &
    298                       ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
    299              kmxm_x = 0.25 * &
    300                       ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    301              kmxp_y = kmxp_x
    302              kmxm_y = kmxm_x
    303 !
    304 !--          Increase diffusion at the outflow boundary in case of
    305 !--          non-cyclic lateral boundaries. Damping is only needed for
    306 !--          velocity components parallel to the outflow boundary in
    307 !--          the direction normal to the outflow boundary.
    308              IF ( .NOT. bc_lr_cyc )  THEN
    309                 kmxp_x = MAX( kmxp_x, km_damp_x(i) )
    310                 kmxm_x = MAX( kmxm_x, km_damp_x(i) )
    311              ENDIF
     265             kmxp = 0.25 * &
     266                    ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
     267             kmxm = 0.25 * &
     268                    ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    312269
    313270             tend(k,j,i) = tend(k,j,i)                                         &
     
    317274                                         ) * ddy2                              &
    318275                                 + (   fxp(j,i) * (                            &
    319                                   kmxp_x * ( v(k,j,i+1) - v(k,j,i)     ) * ddx &
    320                                 + kmxp_y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
     276                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx  &
     277                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy  &
    321278                                                  )                            &
    322279                                     - fxm(j,i) * (                            &
    323                                   kmxm_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
    324                                 + kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
     280                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
     281                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
    325282                                                  )                            &
    326283                                     + wall_v(j,i) * vsus(k)                   &
  • palm/trunk/SOURCE/diffusion_w.f90

    r668 r978  
    44! Current revisions:
    55! -----------------
     6! outflow damping layer removed
     7! kmxm_x/_z and kmxp_x/_z change to kmxm and kmxp
     8! kmym_y/_z and kmyp_y/_z change to kmym and kmyp
    69!
    710! Former revisions:
     
    5457! Call for all grid points
    5558!------------------------------------------------------------------------------!
    56     SUBROUTINE diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, &
    57                             w )
     59    SUBROUTINE diffusion_w( ddzu, ddzw, km, tend, u, v, w )
    5860
    5961       USE control_parameters
     
    6466
    6567       INTEGER ::  i, j, k
    66        REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
    67                    kmyp_z
    68        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg),        &
    69                    km_damp_y(nysg:nyng)
     68       REAL    ::  kmxm, kmxp, kmym, kmyp
     69       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
    7070       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    7171       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
     
    8888!
    8989!--             Interpolate eddy diffusivities on staggered gridpoints
    90                 kmxp_x = 0.25 * &
    91                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
    92                 kmxm_x = 0.25 * &
    93                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
    94                 kmxp_z = kmxp_x
    95                 kmxm_z = kmxm_x
    96                 kmyp_y = 0.25 * &
    97                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
    98                 kmym_y = 0.25 * &
    99                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    100                 kmyp_z = kmyp_y
    101                 kmym_z = kmym_y
    102 !
    103 !--             Increase diffusion at the outflow boundary in case of
    104 !--             non-cyclic lateral boundaries. Damping is only needed for
    105 !--             velocity components parallel to the outflow boundary in
    106 !--             the direction normal to the outflow boundary.
    107                 IF ( .NOT. bc_lr_cyc )  THEN
    108                    kmxp_x = MAX( kmxp_x, km_damp_x(i) )
    109                    kmxm_x = MAX( kmxm_x, km_damp_x(i) )
    110                 ENDIF
    111                 IF ( .NOT. bc_ns_cyc )  THEN
    112                    kmyp_y = MAX( kmyp_y, km_damp_y(j) )
    113                    kmym_y = MAX( kmym_y, km_damp_y(j) )
    114                 ENDIF
     90                kmxp = 0.25 * &
     91                       ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
     92                kmxm = 0.25 * &
     93                       ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
     94                kmyp = 0.25 * &
     95                       ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
     96                kmym = 0.25 * &
     97                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    11598
    11699                tend(k,j,i) = tend(k,j,i)                                      &
    117                       & + ( kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
    118                       &   + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
    119                       &   - kmxm_x * ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
    120                       &   - kmxm_z * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
     100                      & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
     101                      &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)  &
     102                      &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
     103                      &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
    121104                      &   ) * ddx                                              &
    122                       & + ( kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
    123                       &   + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
    124                       &   - kmym_y * ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
    125                       &   - kmym_z * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
     105                      & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
     106                      &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)  &
     107                      &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
     108                      &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
    126109                      &   ) * ddy                                              &
    127110                      & + 2.0 * (                                              &
     
    138121!
    139122!--                Interpolate eddy diffusivities on staggered gridpoints
    140                    kmxp_x = 0.25 * &
    141                             ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
    142                    kmxm_x = 0.25 * &
    143                             ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
    144                    kmxp_z = kmxp_x
    145                    kmxm_z = kmxm_x
    146                    kmyp_y = 0.25 * &
    147                             ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
    148                    kmym_y = 0.25 * &
    149                             ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    150                    kmyp_z = kmyp_y
    151                    kmym_z = kmym_y
    152 !
    153 !--                Increase diffusion at the outflow boundary in case of
    154 !--                non-cyclic lateral boundaries. Damping is only needed for
    155 !--                velocity components parallel to the outflow boundary in
    156 !--                the direction normal to the outflow boundary.
    157                    IF ( .NOT. bc_lr_cyc )  THEN
    158                       kmxp_x = MAX( kmxp_x, km_damp_x(i) )
    159                       kmxm_x = MAX( kmxm_x, km_damp_x(i) )
    160                    ENDIF
    161                    IF ( .NOT. bc_ns_cyc )  THEN
    162                       kmyp_y = MAX( kmyp_y, km_damp_y(j) )
    163                       kmym_y = MAX( kmym_y, km_damp_y(j) )
    164                    ENDIF
     123                   kmxp = 0.25 * &
     124                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
     125                   kmxm = 0.25 * &
     126                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
     127                   kmyp = 0.25 * &
     128                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
     129                   kmym = 0.25 * &
     130                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    165131
    166132                   tend(k,j,i) = tend(k,j,i)                                   &
    167133                                 + (   fwxp(j,i) * (                           &
    168                             kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
    169                           + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
     134                            kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
     135                          + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)  &
    170136                                                   )                           &
    171137                                     - fwxm(j,i) * (                           &
    172                             kmxm_x * ( w(k,j,i)     - w(k,j,i-1) ) * ddx       &
    173                           + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
     138                            kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
     139                          + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)  &
    174140                                                   )                           &
    175141                                     + wall_w_x(j,i) * wsus(k,j,i)             &
    176142                                   ) * ddx                                     &
    177143                                 + (   fwyp(j,i) * (                           &
    178                             kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
    179                           + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
     144                            kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
     145                          + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)  &
    180146                                                   )                           &
    181147                                     - fwym(j,i) * (                           &
    182                             kmym_y * ( w(k,j,i)     - w(k,j-1,i) ) * ddy       &
    183                           + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
     148                            kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
     149                          + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)  &
    184150                                                   )                           &
    185151                                     + wall_w_y(j,i) * wsvs(k,j,i)             &
     
    201167! Call for grid point i,j
    202168!------------------------------------------------------------------------------!
    203     SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
    204                                tend, u, v, w )
     169    SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, tend, u, v, w )
    205170
    206171       USE control_parameters
     
    211176
    212177       INTEGER ::  i, j, k
    213        REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
    214                    kmyp_z
    215        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg),        &
    216                    km_damp_y(nysg:nyng)
     178       REAL    ::  kmxm, kmxp, kmym, kmyp
     179       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
    217180       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    218181       REAL, DIMENSION(nzb:nzt+1)      ::  wsus, wsvs
     
    223186!
    224187!--       Interpolate eddy diffusivities on staggered gridpoints
    225           kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
    226           kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
    227           kmxp_z = kmxp_x
    228           kmxm_z = kmxm_x
    229           kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
    230           kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    231           kmyp_z = kmyp_y
    232           kmym_z = kmym_y
    233 !
    234 !--       Increase diffusion at the outflow boundary in case of non-cyclic
    235 !--       lateral boundaries. Damping is only needed for velocity components
    236 !--       parallel to the outflow boundary in the direction normal to the
    237 !--       outflow boundary.
    238           IF ( .NOT. bc_lr_cyc )  THEN
    239              kmxp_x = MAX( kmxp_x, km_damp_x(i) )
    240              kmxm_x = MAX( kmxm_x, km_damp_x(i) )
    241           ENDIF
    242           IF ( .NOT. bc_ns_cyc )  THEN
    243              kmyp_y = MAX( kmyp_y, km_damp_y(j) )
    244              kmym_y = MAX( kmym_y, km_damp_y(j) )
    245           ENDIF
     188          kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
     189          kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
     190          kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
     191          kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    246192
    247193          tend(k,j,i) = tend(k,j,i)                                            &
    248                       & + ( kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
    249                       &   + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
    250                       &   - kmxm_x * ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
    251                       &   - kmxm_z * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
     194                      & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
     195                      &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)  &
     196                      &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
     197                      &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
    252198                      &   ) * ddx                                              &
    253                       & + ( kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
    254                       &   + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
    255                       &   - kmym_y * ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
    256                       &   - kmym_z * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
     199                      & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
     200                      &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)  &
     201                      &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
     202                      &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
    257203                      &   ) * ddy                                              &
    258204                      & + 2.0 * (                                              &
     
    285231!
    286232!--          Interpolate eddy diffusivities on staggered gridpoints
    287              kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
    288              kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
    289              kmxp_z = kmxp_x
    290              kmxm_z = kmxm_x
    291              kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
    292              kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    293              kmyp_z = kmyp_y
    294              kmym_z = kmym_y
    295 !
    296 !--          Increase diffusion at the outflow boundary in case of
    297 !--          non-cyclic lateral boundaries. Damping is only needed for
    298 !--          velocity components parallel to the outflow boundary in
    299 !--          the direction normal to the outflow boundary.
    300              IF ( .NOT. bc_lr_cyc )  THEN
    301                 kmxp_x = MAX( kmxp_x, km_damp_x(i) )
    302                 kmxm_x = MAX( kmxm_x, km_damp_x(i) )
    303              ENDIF
    304              IF ( .NOT. bc_ns_cyc )  THEN
    305                 kmyp_y = MAX( kmyp_y, km_damp_y(j) )
    306                 kmym_y = MAX( kmym_y, km_damp_y(j) )
    307              ENDIF
     233             kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
     234             kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
     235             kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
     236             kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    308237
    309238             tend(k,j,i) = tend(k,j,i)                                         &
    310239                                 + (   fwxp(j,i) * (                           &
    311                             kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
    312                           + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
     240                            kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
     241                          + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)  &
    313242                                                   )                           &
    314243                                     - fwxm(j,i) * (                           &
    315                             kmxm_x * ( w(k,j,i)     - w(k,j,i-1) ) * ddx       &
    316                           + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
     244                            kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
     245                          + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)  &
    317246                                                   )                           &
    318247                                     + wall_w_x(j,i) * wsus(k)                 &
    319248                                   ) * ddx                                     &
    320249                                 + (   fwyp(j,i) * (                           &
    321                             kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
    322                           + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
     250                            kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
     251                          + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)  &
    323252                                                   )                           &
    324253                                     - fwym(j,i) * (                           &
    325                             kmym_y * ( w(k,j,i)     - w(k,j-1,i) ) * ddy       &
    326                           + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
     254                            kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
     255                          + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)  &
    327256                                                   )                           &
    328257                                     + wall_w_y(j,i) * wsvs(k)                 &
  • palm/trunk/SOURCE/header.f90

    r965 r978  
    44! Current revisions:
    55! -----------------
    6 !
     6! -km_damp_max, outflow_damping_width
     7! +pt_damping_factor, pt_damping_width
     8! +z0h
    79!
    810! Former revisions:
     
    704706
    705707    IF ( prandtl_layer )  THEN
    706        WRITE ( io, 305 )  (zu(1)-zu(0)), roughness_length, kappa, &
     708       WRITE ( io, 305 )  (zu(1)-zu(0)), roughness_length, &
     709                          z0h_factor*roughness_length, kappa, &
    707710                          rif_min, rif_max
    708711       IF ( .NOT. constant_heatflux )  WRITE ( io, 308 )
     
    721724    WRITE ( io, 317 )  bc_lr, bc_ns
    722725    IF ( .NOT. bc_lr_cyc  .OR.  .NOT. bc_ns_cyc )  THEN
    723        WRITE ( io, 318 )  outflow_damping_width, km_damp_max
     726       WRITE ( io, 318 )  pt_damping_width, pt_damping_factor       
    724727       IF ( turbulent_inflow )  THEN
    725728          WRITE ( io, 319 )  recycling_width, recycling_plane, &
     
    17791782305 FORMAT (//'    Prandtl-Layer between bottom surface and first ', &
    17801783               'computational u,v-level:'// &
    1781              '       zp = ',F6.2,' m   z0 = ',F6.4,' m   kappa = ',F4.2/ &
     1784             '       zp = ',F6.2,' m   z0 = ',F6.4,' m   z0h = ',F7.5,&
     1785             ' m   kappa = ',F4.2/ &
    17821786             '       Rif value range:   ',F6.2,' <= rif <=',F6.2)
    17831787306 FORMAT ('       Predefined constant heatflux:   ',F9.6,' K m/s')
     
    17971801            '       left/right:  ',A/    &
    17981802            '       north/south: ',A)
    1799 318 FORMAT (/'       outflow damping layer width: ',I3,' gridpoints with km_', &
    1800                     'max =',F5.1,' m**2/s')
     1803318 FORMAT (/'       pt damping layer width = ',F7.2,' m, pt ', &
     1804                    'damping factor = ',F6.4)
    18011805319 FORMAT ('       turbulence recycling at inflow switched on'/ &
    18021806            '       width of recycling domain: ',F7.1,' m   grid index: ',I4/ &
  • palm/trunk/SOURCE/init_1d_model.f90

    r668 r978  
    44! Current revisions:
    55! -----------------
    6 !
     6! roughness length for scalar quantities z0h1d added
    77!
    88! Former revisions:
     
    149149    vsws1d = 0.0; vsws1d_m = 0.0
    150150    z01d = roughness_length
     151    z0h1d = z0h_factor * z01d
    151152    IF ( humidity .OR. passive_scalar )  qs1d = 0.0
    152153
     
    445446!--                Stable stratification
    446447                   ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /      &
    447                           ( LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * &
    448                                           ( zu(nzb+1) - z01d ) / zu(nzb+1) &
     448                          ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * &
     449                                          ( zu(nzb+1) - z0h1d ) / zu(nzb+1) &
    449450                          )
    450451                ELSE
     
    452453!--                Unstable stratification
    453454                   a = SQRT( 1.0 - 16.0 * rif1d(nzb+1) )
    454                    b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z01d )
     455                   b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z0h1d )
    455456!
    456457!--                In the borderline case the formula for stable stratification
     
    459460                   IF ( a == 0.0  .OR.  b == 0.0 )  THEN
    460461                      ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /      &
    461                              ( LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * &
    462                                              ( zu(nzb+1) - z01d ) / zu(nzb+1) &
     462                             ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * &
     463                                             ( zu(nzb+1) - z0h1d ) / zu(nzb+1) &
    463464                             )
    464465                   ELSE
     
    578579!--                Stable stratification
    579580                   qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /        &
    580                           ( LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * &
    581                                           ( zu(nzb+1) - z01d ) / zu(nzb+1) &
     581                          ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * &
     582                                          ( zu(nzb+1) - z0h1d ) / zu(nzb+1) &
    582583                          )
    583584                ELSE
     
    585586!--                Unstable stratification
    586587                   a = SQRT( 1.0 - 16.0 * rif1d(nzb+1) )
    587                    b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z01d )
     588                   b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z0h1d )
    588589!
    589590!--                In the borderline case the formula for stable stratification
     
    592593                   IF ( a == 1.0  .OR.  b == 1.0 )  THEN
    593594                      qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /        &
    594                              ( LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * &
    595                                              ( zu(nzb+1) - z01d ) / zu(nzb+1) &
     595                             ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * &
     596                                             ( zu(nzb+1) - z0h1d ) / zu(nzb+1) &
    596597                             )
    597598                   ELSE
  • palm/trunk/SOURCE/init_3d_model.f90

    r850 r978  
    77! Current revisions:
    88! ------------------
    9 !
     9! outflow damping layer removed
     10! roughness length for scalar quantites z0h added
     11! damping zone for the potential temperatur in case of non-cyclic lateral
     12! boundaries added
     13! initialization of ptdf_x, ptdf_y
     14! initialization of c_u_m, c_u_m_l, c_v_m, c_v_m_l, c_w_m, c_w_m_l
    1015!
    1116! Former revisions:
     
    164169    USE control_parameters
    165170    USE cpulog
     171    USE grid_variables
    166172    USE indices
    167173    USE interfaces
     
    211217              sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions),                &
    212218              ts_value(dots_max,0:statistic_regions) )
    213     ALLOCATE( km_damp_x(nxlg:nxrg), km_damp_y(nysg:nyng) )
     219    ALLOCATE( ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng) )
    214220
    215221    ALLOCATE( rif_1(nysg:nyng,nxlg:nxrg), shf_1(nysg:nyng,nxlg:nxrg), &
     
    218224              uswst_1(nysg:nyng,nxlg:nxrg),                               &
    219225              vsws_1(nysg:nyng,nxlg:nxrg),                                &
    220               vswst_1(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg) )
     226              vswst_1(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg),      &
     227              z0h(nysg:nyng,nxlg:nxrg) )
    221228
    222229    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     
    413420                 c_w(nzb:nzt+1,nxlg:nxrg) )
    414421    ENDIF
     422    IF ( outflow_l  .OR.  outflow_r .OR. outflow_s .OR. outflow_n )  THEN
     423       ALLOCATE( c_u_m_l(nzb:nzt+1), c_v_m_l(nzb:nzt+1), c_w_m_l(nzb:nzt+1) )                   
     424       ALLOCATE( c_u_m(nzb:nzt+1), c_v_m(nzb:nzt+1), c_w_m(nzb:nzt+1) )
     425    ENDIF
     426
    415427
    416428!
     
    847859
    848860          z0 = roughness_length
     861          z0h = z0h_factor * z0
    849862
    850863          IF ( .NOT. constant_heatflux )  THEN
     
    15331546
    15341547!
    1535 !-- Initialize diffusivities used within the outflow damping layer in case of
    1536 !-- non-cyclic lateral boundaries. A linear increase is assumed over the first
    1537 !-- half of the width of the damping layer
    1538     IF ( bc_lr_dirrad )  THEN
    1539 
    1540        DO  i = nxlg, nxrg
    1541           IF ( i >= nx - outflow_damping_width )  THEN
    1542              km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
    1543                             ( i - ( nx - outflow_damping_width ) ) /   &
    1544                             REAL( outflow_damping_width/2 )            &
    1545                                              )
    1546           ELSE
    1547              km_damp_x(i) = 0.0
     1548!-- Initialize damping zone for the potential temperature in case of
     1549!-- non-cyclic lateral boundaries. The damping zone has the maximum value
     1550!-- at the inflow boundary and decreases to zero at pt_damping_width.
     1551    ptdf_x = 0.0
     1552    ptdf_y = 0.0
     1553    IF ( bc_lr_dirrad .OR. bc_lr_dirneu )  THEN
     1554       DO i = nxl, nxr
     1555          IF ( ( i * dx ) < pt_damping_width )  THEN
     1556             ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5 *        &
     1557                         REAL( pt_damping_width - i * dx ) / (        &
     1558                         REAL( pt_damping_width )            ) ) )**2
    15481559          ENDIF
    15491560       ENDDO
    1550 
    1551     ELSEIF ( bc_lr_raddir )  THEN
    1552 
    1553        DO  i = nxlg, nxrg
    1554           IF ( i <= outflow_damping_width )  THEN
    1555              km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
    1556                             ( outflow_damping_width - i ) /            &
    1557                             REAL( outflow_damping_width/2 )            &
    1558                                              )
    1559           ELSE
    1560              km_damp_x(i) = 0.0
     1561    ELSEIF ( bc_lr_raddir .OR. bc_lr_neudir )  THEN
     1562       DO i = nxl, nxr
     1563          IF ( ( i * dx ) > ( nx * dx - pt_damping_width ) )  THEN
     1564             ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5 *                 &
     1565                         REAL( ( i - nx ) * dx + pt_damping_width ) / (        &
     1566                         REAL( pt_damping_width )                     ) ) )**2       
     1567          ENDIF
     1568       ENDDO
     1569    ELSEIF ( bc_ns_dirrad .OR. bc_ns_dirneu )  THEN
     1570       DO j = nys, nyn
     1571          IF ( ( j * dy ) > ( ny * dy - pt_damping_width ) )  THEN
     1572             ptdf_y(j) = pt_damping_factor * ( SIN( pi * 0.5 *                 &
     1573                         REAL( ( j - ny ) * dy + pt_damping_width ) / (        &
     1574                         REAL( pt_damping_width )                     ) ) )**2       
     1575          ENDIF
     1576       ENDDO
     1577    ELSEIF ( bc_ns_raddir .OR. bc_ns_neudir )  THEN
     1578       DO j = nys, nyn
     1579          IF ( ( j * dy ) < pt_damping_width )  THEN
     1580             ptdf_y(j) = pt_damping_factor * ( SIN( pi * 0.5 *        &
     1581                         REAL( pt_damping_width - j * dy ) / (        &
     1582                         REAL( pt_damping_width )            ) ) )**2       
    15611583          ENDIF
    15621584       ENDDO
    1563 
    1564     ENDIF
    1565 
    1566     IF ( bc_ns_dirrad )  THEN
    1567 
    1568        DO  j = nysg, nyng
    1569           IF ( j >= ny - outflow_damping_width )  THEN
    1570              km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
    1571                             ( j - ( ny - outflow_damping_width ) ) /   &
    1572                             REAL( outflow_damping_width/2 )            &
    1573                                              )
    1574           ELSE
    1575              km_damp_y(j) = 0.0
    1576           ENDIF
    1577        ENDDO
    1578 
    1579     ELSEIF ( bc_ns_raddir )  THEN
    1580 
    1581        DO  j = nysg, nyng
    1582           IF ( j <= outflow_damping_width )  THEN
    1583              km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
    1584                             ( outflow_damping_width - j ) /            &
    1585                             REAL( outflow_damping_width/2 )            &
    1586                                              )
    1587           ELSE
    1588              km_damp_y(j) = 0.0
    1589           ENDIF
    1590        ENDDO
    1591 
    15921585    ENDIF
    15931586
  • palm/trunk/SOURCE/init_grid.f90

    r928 r978  
    44! Current revisions:
    55! -----------------
    6 !
     6! Bugfix: nzb_max is set to nzt at non-cyclic lateral boundaries
     7! Bugfix: Set wall_flags_0 for inflow boundary
    78!
    89! Former revisions:
     
    608609!-- Determine the maximum level of topography. Furthermore it is used for
    609610!-- steering the degradation of order of the applied advection scheme.
     611!-- In case of non-cyclic lateral boundaries, the order of the advection
     612!-- scheme is reduced at the lateral boundaries up to nzt.
    610613    nzb_max = MAXVAL( nzb_local )
     614    IF ( inflow_l .OR. outflow_l .OR. inflow_r .OR. outflow_r .OR.    &
     615         inflow_n .OR. outflow_n .OR. inflow_s .OR. outflow_s )  THEN
     616         nzb_max = nzt
     617    ENDIF
     618
    611619!
    612620!-- Consistency checks and index array initialization are only required for
     
    10791087!--             scalar - x-direction
    10801088!--             WS1 (0), WS3 (1), WS5 (2)
    1081                 IF ( k <= nzb_s_inner(j,i+1) .OR. ( outflow_l .AND. i == nxl ) &
    1082                      .OR. ( outflow_r .AND. i == nxr ) )  THEN
     1089                IF ( k <= nzb_s_inner(j,i+1) .OR. ( ( inflow_l .OR. outflow_l )&
     1090                     .AND. i == nxl ) .OR. ( ( inflow_r .OR. outflow_r )       &
     1091                     .AND. i == nxr ) )  THEN
    10831092                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 0 )
    10841093                ELSEIF ( k <= nzb_s_inner(j,i+2) .OR. k <= nzb_s_inner(j,i-1)  &
    1085                          .OR. ( outflow_r .AND. i == nxr-1 )                   &
    1086                          .OR. ( outflow_l .AND. i == nxlu  ) )  THEN
     1094                         .OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr-1 ) &
     1095                         .OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu  ) &
     1096                       )  THEN
    10871097                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 1 )
    10881098                ELSE
     
    10921102!--             scalar - y-direction
    10931103!--             WS1 (3), WS3 (4), WS5 (5)
    1094                 IF ( k <= nzb_s_inner(j+1,i) .OR. ( outflow_s .AND. j == nys ) &
    1095                      .OR. ( outflow_n .AND. j == nyn ) )  THEN
     1104                IF ( k <= nzb_s_inner(j+1,i) .OR. ( ( inflow_s .OR. outflow_s )&
     1105                     .AND. j == nys ) .OR. ( ( inflow_n .OR. outflow_n )       &
     1106                     .AND. j == nyn ) )  THEN
    10961107                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 3 )
    10971108!--             WS3
    10981109                ELSEIF ( k <= nzb_s_inner(j+2,i) .OR. k <= nzb_s_inner(j-1,i)  &
    1099                         .OR. ( outflow_s .AND. j == nysv )                     &
    1100                         .OR. ( outflow_n .AND. j == nyn-1  ) )  THEN
     1110                         .OR. ( ( inflow_s .OR. outflow_s ) .AND. j == nysv  ) &
     1111                         .OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn-1 ) &
     1112                       )  THEN
    11011113                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 4 )
    11021114!--             WS5
     
    11351147!--             WS1 (9), WS3 (10), WS5 (11)
    11361148                IF ( k <= nzb_u_inner(j,i+1)                                  &
    1137                      .OR. ( outflow_r .AND. i == nxr ) )  THEN
     1149                     .OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu )     &
     1150                     .OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr  )     &
     1151                   )  THEN
    11381152                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 9 )
    11391153                ELSEIF ( k <= nzb_u_inner(j,i+2) .OR. k <= nzb_u_inner(j,i-1) &
    1140                          .OR. ( outflow_r .AND. i == nxr-1 )                  &
    1141                          .OR. ( outflow_l .AND. i == nxlu  ) )  THEN
     1154                         .OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr-1 )&
     1155                         .OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu+1)&
     1156                       )  THEN
    11421157                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 10 )
    11431158                ELSE
    11441159                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 11 )
    11451160                ENDIF
     1161
    11461162!
    11471163!--             u component - y-direction
    11481164!--             WS1 (12), WS3 (13), WS5 (14)
    1149                 IF ( k <= nzb_u_inner(j+1,i) .OR. ( outflow_s .AND. j == nys ) &
    1150                     .OR. ( outflow_n .AND. j == nyn ) )  THEN
     1165                IF ( k <= nzb_u_inner(j+1,i) .OR. ( ( inflow_s .OR. outflow_s )&
     1166                     .AND. j == nys ) .OR. ( ( inflow_n .OR. outflow_n )       &
     1167                     .AND. j == nyn ) )  THEN
    11511168                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 )
    11521169                ELSEIF ( k <= nzb_u_inner(j+2,i) .OR. k <= nzb_u_inner(j-1,i)  &
    1153                     .OR. ( outflow_s .AND. j == nysv )                         &
    1154                     .OR. ( outflow_n .AND. j == nyn-1  ) )  THEN
     1170                         .OR. ( ( inflow_s .OR. outflow_s ) .AND. j == nysv  ) &
     1171                         .OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn-1 ) &
     1172                       )  THEN
    11551173                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 13 )
    11561174                ELSE
     
    11811199!--             v component - x-direction
    11821200!--             WS1 (18), WS3 (19), WS5 (20)
    1183                 IF ( k <= nzb_v_inner(j,i+1) .OR. ( outflow_l .AND. i == nxl ) &
    1184                     .OR. ( outflow_r .AND. i == nxr ) )  THEN
     1201                IF ( k <= nzb_v_inner(j,i+1) .OR. ( ( inflow_l .OR. outflow_l )&
     1202                     .AND. i == nxl ) .OR. (( inflow_r .OR. outflow_r )        &
     1203                     .AND. i == nxr ) )  THEN
    11851204                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 )
    11861205!--             WS3
    11871206                ELSEIF ( k <= nzb_v_inner(j,i+2) .OR. k <= nzb_v_inner(j,i-1)  &
    1188                     .OR. ( outflow_r .AND. i == nxr-1 )                        &
    1189                     .OR. ( outflow_l .AND. i == nxlu  ) )  THEN
     1207                         .OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr-1 ) &
     1208                         .OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu  ) &
     1209                       )  THEN
    11901210                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 19 )
    11911211                ELSE
     
    11961216!--             WS1 (21), WS3 (22), WS5 (23)
    11971217                IF ( k <= nzb_v_inner(j+1,i)                                   &
    1198                      .OR. ( outflow_n .AND. j == nyn ) )  THEN
     1218                     .OR. ( ( inflow_s .OR. outflow_s ) .AND. i == nysv )      &
     1219                     .OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn  )      &
     1220                   )  THEN
    11991221                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 )
    12001222                ELSEIF ( k <= nzb_v_inner(j+2,i) .OR. k <= nzb_v_inner(j-1,i)  &
    1201                          .OR. ( outflow_s .AND. j == nysv )                    &
    1202                          .OR. ( outflow_n .AND. j == nyn-1  ) )  THEN
     1223                         .OR. ( ( inflow_s .OR. outflow_s ) .AND. j == nysv+1 )&
     1224                         .OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn-1  )&
     1225                       )  THEN
    12031226                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 22 )
    12041227                ELSE
     
    12281251!--             w component - x-direction
    12291252!--             WS1 (27), WS3 (28), WS5 (29)
    1230                 IF ( k <= nzb_w_inner(j,i+1) .OR. ( outflow_l .AND. i == nxl ) &
    1231                     .OR. ( outflow_r .AND. i == nxr ) )  THEN
     1253                IF ( k <= nzb_w_inner(j,i+1) .OR. ( ( inflow_l .OR. outflow_l )&
     1254                     .AND. i == nxl ) .OR. ( ( inflow_r .OR. outflow_r )       &
     1255                     .AND. i == nxr ) )  THEN
    12321256                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 )
    12331257                ELSEIF ( k <= nzb_w_inner(j,i+2) .OR. k <= nzb_w_inner(j,i-1)  &
    1234                          .OR. ( outflow_r .AND. i == nxr-1 )                   &
    1235                          .OR. ( outflow_l .AND. i == nxlu  ) )  THEN
     1258                         .OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr-1 ) &
     1259                         .OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu  ) &
     1260                       )  THEN
    12361261                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 28 )
    12371262                ELSE
     
    12411266!--             w component - y-direction
    12421267!--             WS1 (30), WS3 (31), WS5 (32)
    1243                 IF ( k <= nzb_w_inner(j+1,i) .OR. ( outflow_s .AND. j == nys ) &
    1244                      .OR. ( outflow_n .AND. j == nyn ) )  THEN
     1268                IF ( k <= nzb_w_inner(j+1,i) .OR. ( ( inflow_s .OR. outflow_s )&
     1269                     .AND. j == nys ) .OR. ( ( inflow_n .OR. outflow_n )       &
     1270                     .AND. j == nyn ) )  THEN
    12451271                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 )
    12461272                ELSEIF ( k <= nzb_w_inner(j+2,i) .OR. k <= nzb_w_inner(j-1,i)  &
    1247                          .OR. ( outflow_s .AND. j == nysv )                    &
    1248                          .OR. ( outflow_n .AND. j == nyn-1  ) )  THEN
     1273                         .OR. ( ( inflow_s .OR. outflow_s ) .AND. j == nysv  ) &
     1274                         .OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn-1 ) &
     1275                       )  THEN
    12491276                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 31 )
    12501277                ELSE
  • palm/trunk/SOURCE/init_masks.f90

    r810 r978  
    44! Current revisions:
    55! -----------------
    6 !
     6! +z0h*
    77!
    88! Former revisions:
     
    257257                unit = 'psu'
    258258
    259              CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'z0*' )
     259             CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'z0*', 'z0h*' )
    260260                WRITE ( message_string, * ) 'illegal value for data_',&
    261261                     'output: "', TRIM( var ), '" is only allowed', &
  • palm/trunk/SOURCE/init_pegrid.f90

    r810 r978  
    44! Current revisions:
    55! -----------------
     6! dirichlet/neumann and neumann/dirichlet added
     7! nxlu and nysv are also calculated for inflow boundary
    68!
    79! ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
     
    11081110!-- horizontal boundary conditions.
    11091111    IF ( pleft == MPI_PROC_NULL )  THEN
    1110        IF ( bc_lr == 'dirichlet/radiation' )  THEN
     1112       IF ( bc_lr == 'dirichlet/radiation' .OR. bc_lr == 'dirichlet/neumann' )  THEN
    11111113          inflow_l  = .TRUE.
    1112        ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
     1114       ELSEIF ( bc_lr == 'radiation/dirichlet' .OR. bc_lr == 'neumann/dirichlet' )  THEN
    11131115          outflow_l = .TRUE.
    11141116       ENDIF
     
    11161118
    11171119    IF ( pright == MPI_PROC_NULL )  THEN
    1118        IF ( bc_lr == 'dirichlet/radiation' )  THEN
     1120       IF ( bc_lr == 'dirichlet/radiation' .OR. bc_lr == 'dirichlet/neumann' )  THEN
    11191121          outflow_r = .TRUE.
    1120        ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
     1122       ELSEIF ( bc_lr == 'radiation/dirichlet' .OR. bc_lr == 'neumann/dirichlet' )  THEN
    11211123          inflow_r  = .TRUE.
    11221124       ENDIF
     
    11241126
    11251127    IF ( psouth == MPI_PROC_NULL )  THEN
    1126        IF ( bc_ns == 'dirichlet/radiation' )  THEN
     1128       IF ( bc_ns == 'dirichlet/radiation' .OR. bc_ns == 'dirichlet/neumann' )  THEN
    11271129          outflow_s = .TRUE.
    1128        ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
     1130       ELSEIF ( bc_ns == 'radiation/dirichlet' .OR. bc_ns == 'neumann/dirichlet' )  THEN
    11291131          inflow_s  = .TRUE.
    11301132       ENDIF
     
    11321134
    11331135    IF ( pnorth == MPI_PROC_NULL )  THEN
    1134        IF ( bc_ns == 'dirichlet/radiation' )  THEN
     1136       IF ( bc_ns == 'dirichlet/radiation' .OR. bc_ns == 'dirichlet/neumann' )  THEN
    11351137          inflow_n  = .TRUE.
    1136        ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
     1138       ELSEIF ( bc_ns == 'radiation/dirichlet' .OR. bc_ns == 'neumann/dirichlet' )  THEN
    11371139          outflow_n = .TRUE.
    11381140       ENDIF
     
    11641166
    11651167#elif ! defined ( __parallel )
    1166     IF ( bc_lr == 'dirichlet/radiation' )  THEN
     1168    IF ( bc_lr == 'dirichlet/radiation' .OR. bc_lr == 'dirichlet/neumann' )  THEN
    11671169       inflow_l  = .TRUE.
    11681170       outflow_r = .TRUE.
    1169     ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
     1171    ELSEIF ( bc_lr == 'radiation/dirichlet' .OR. bc_lr == 'neumann/dirichlet' )  THEN
    11701172       outflow_l = .TRUE.
    11711173       inflow_r  = .TRUE.
    11721174    ENDIF
    11731175
    1174     IF ( bc_ns == 'dirichlet/radiation' )  THEN
     1176    IF ( bc_ns == 'dirichlet/radiation' .OR. bc_ns == 'dirichlet/neumann' )  THEN
    11751177       inflow_n  = .TRUE.
    11761178       outflow_s = .TRUE.
    1177     ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
     1179    ELSEIF ( bc_ns == 'radiation/dirichlet' .OR. bc_ns == 'neumann/dirichlet' )  THEN
    11781180       outflow_n = .TRUE.
    11791181       inflow_s  = .TRUE.
     
    11821184
    11831185!
    1184 !-- At the outflow, u or v, respectively, have to be calculated for one more
    1185 !-- grid point.
    1186     IF ( outflow_l )  THEN
     1186!-- At the inflow or outflow, u or v, respectively, have to be calculated for
     1187!-- one more grid point.
     1188    IF ( inflow_l .OR. outflow_l )  THEN
    11871189       nxlu = nxl + 1
    11881190    ELSE
    11891191       nxlu = nxl
    11901192    ENDIF
    1191     IF ( outflow_s )  THEN
     1193    IF ( inflow_s .OR. outflow_s )  THEN
    11921194       nysv = nys + 1
    11931195    ELSE
  • palm/trunk/SOURCE/modules.f90

    r965 r978  
    44! Current revisions:
    55! -----------------
    6 !
     6! +c_u_m, c_u_m_l, c_v_m, c_v_m_l, c_w_m, c_w_m_l,
     7! +bc_lr_dirneu, bc_lr_neudir, bc_ns_dirneu, bc_ns_neudir
     8! -km_damp_x, km_damp_y, km_damp_max, outflow_damping_width
     9! +z0h, z0h_av, z0h_factor, z0h1d
     10! +ptdf_x, ptdf_y, pt_damping_width, pt_damping_factor
    711!
    812! Former revisions:
     
    324328
    325329    REAL, DIMENSION(:), ALLOCATABLE ::                                         &
    326           ddzu, ddzu_pres, dd2zu, dzu, ddzw, dzw, hyp, inflow_damping_factor,  &
    327           km_damp_x, km_damp_y, lad, l_grid, pt_init, q_init, rdf, rdf_sc,     &
    328           sa_init, ug, u_init, u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, &
    329           w_subs, zu, zw
     330          c_u_m, c_u_m_l, c_v_m, c_v_m_l, c_w_m, c_w_m_l, ddzu, ddzu_pres,     &
     331          dd2zu, dzu, ddzw, dzw, hyp, inflow_damping_factor, lad, l_grid,      &
     332          ptdf_x, ptdf_y, pt_init, q_init, rdf, rdf_sc, sa_init, ug, u_init,  &
     333          u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, w_subs, zu, zw
    330334
    331335    REAL, DIMENSION(:,:), ALLOCATABLE ::                                       &
    332           c_u, c_v, c_w, dzu_mg, dzw_mg, f1_mg, f2_mg, f3_mg,                  &
    333           mean_inflow_profiles, pt_slope_ref, qs, qswst_remote, ts, us, z0,    &
    334           flux_s_pt, diss_s_pt, flux_s_e, diss_s_e, flux_s_q, diss_s_q,        &
    335           flux_s_sa, diss_s_sa, flux_s_u, flux_s_v, flux_s_w, diss_s_u,        &
    336           diss_s_v, diss_s_w, total_2d_o, total_2d_a
     336          c_u, c_v, c_w, diss_s_e, diss_s_pt, diss_s_q, diss_s_sa, diss_s_u,   &
     337          diss_s_v, diss_s_w, dzu_mg, dzw_mg, flux_s_e, flux_s_pt, flux_s_q,   &
     338          flux_s_sa, flux_s_u, flux_s_v, flux_s_w, f1_mg, f2_mg, f3_mg,        &
     339          mean_inflow_profiles, pt_slope_ref, qs, qswst_remote, total_2d_a,    &
     340          total_2d_o, ts, us, z0, z0h
    337341
    338342    REAL, DIMENSION(:,:), ALLOCATABLE, TARGET ::                               &
     
    347351
    348352    REAL, DIMENSION(:,:,:), ALLOCATABLE ::                                     &
    349           canopy_heat_flux, cdc, d, de_dx, de_dy, de_dz, diss, lad_s, lad_u,   &
    350           lad_v, lad_w, lai, l_wall, p_loc, sec, sls, tend, u_m_l, u_m_n,      &
    351           u_m_r, u_m_s, v_m_l, v_m_n, v_m_r, v_m_s, w_m_l, w_m_n, w_m_r,       &
    352           w_m_s, flux_l_pt, diss_l_pt, flux_l_e, diss_l_e, flux_l_q, diss_l_q, &
    353           flux_l_sa, diss_l_sa, flux_l_u, flux_l_v, flux_l_w, diss_l_u,        &
    354           diss_l_v, diss_l_w
     353          canopy_heat_flux, cdc, d, de_dx, de_dy, de_dz, diss, diss_l_e,       &
     354          diss_l_pt, diss_l_q, diss_l_sa, diss_l_u, diss_l_v, diss_l_w,        &
     355          flux_l_e, flux_l_pt, flux_l_q, flux_l_sa, flux_l_u, flux_l_v,        &
     356          flux_l_w, lad_s, lad_u, lad_v, lad_w, lai, l_wall, p_loc, sec, sls, &
     357          tend, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n, v_m_r, v_m_s, w_m_l, &
     358          w_m_n, w_m_r, w_m_s
    355359
    356360    REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                             &
     
    388392
    389393    REAL, DIMENSION(:,:), ALLOCATABLE ::  lwp_av, precipitation_rate_av,       &
    390                                           qsws_av, shf_av,ts_av, us_av, z0_av
     394                                          qsws_av, shf_av,ts_av, us_av, z0_av, &
     395                                          z0h_av
    391396
    392397    REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
     
    549554                netcdf_data_format = 2, ngsrb = 2, nsor = 20, &
    550555                nsor_ini = 100, n_sor, normalizing_region = 0, &
    551                 nz_do3d = -9999, outflow_damping_width = -1, &
    552                 pch_index = 0, prt_time_count = 0, recycling_plane, runnr = 0, &
     556                nz_do3d = -9999, pch_index = 0, prt_time_count = 0, &
     557                recycling_plane, runnr = 0, &
    553558                skip_do_avs = 0, subdomain_size, terminate_coupled = 0, &
    554559                terminate_coupled_remote = 0, timestep_count = 0
     
    582587
    583588    LOGICAL ::  adjust_mixing_length = .FALSE., avs_output = .FALSE., &
    584                 bc_lr_cyc =.TRUE., bc_lr_dirrad = .FALSE., &
     589                bc_lr_cyc =.TRUE., bc_lr_dirneu = .FALSE., &
     590                bc_lr_dirrad = .FALSE., bc_lr_neudir = .FALSE., &
    585591                bc_lr_raddir = .FALSE., bc_ns_cyc = .TRUE., &
    586                 bc_ns_dirrad = .FALSE., bc_ns_raddir = .FALSE., &
     592                bc_ns_dirneu = .FALSE., bc_ns_dirrad = .FALSE., &
     593                bc_ns_neudir = .FALSE., bc_ns_raddir = .FALSE., &
    587594                call_psolver_at_all_substeps = .TRUE., &
    588595                cloud_droplets = .FALSE., cloud_physics = .FALSE., &
     
    659666             f = 0.0, fs = 0.0, g = 9.81, inflow_damping_height = 9999999.9, &
    660667             inflow_damping_width = 9999999.9, kappa = 0.4, km_constant = -1.0,&
    661              km_damp_max = -1.0, lad_surface = 0.0,  &
     668             lad_surface = 0.0,  &
    662669             leaf_surface_concentration = 0.0, long_filter_factor = 0.0, &
    663670             mask_scale_x = 1.0, mask_scale_y = 1.0, mask_scale_z = 1.0, &
     
    670677             phi = 55.0, prandtl_number = 1.0, &
    671678             precipitation_amount_interval = 9999999.9, prho_reference, &
     679             pt_damping_factor = 0.0, pt_damping_width = 0.0, &
    672680             pt_reference = 9999999.9, pt_slope_offset = 0.0, &
    673681             pt_surface = 300.0, pt_surface_initial_change = 0.0, &
     
    700708             ups_limit_v = 0.0, ups_limit_w = 0.0, vg_surface = 0.0, &
    701709             v_bulk = 0.0, v_gtrans = 0.0, wall_adjustment_factor = 1.8, &
    702              z_max_do2d = -1.0
     710             z_max_do2d = -1.0, z0h_factor = 1.0
    703711
    704712    REAL ::  do2d_xy_last_time(0:1) = -1.0, do2d_xz_last_time(0:1) = -1.0, &
     
    10571065                qs1d, simulated_time_1d = 0.0, time_pr_1d = 0.0, &
    10581066                time_run_control_1d = 0.0, ts1d, us1d, usws1d, usws1d_m, &
    1059                 vsws1d, vsws1d_m, z01d
     1067                vsws1d, vsws1d_m, z01d, z0h1d
    10601068
    10611069
  • palm/trunk/SOURCE/parin.f90

    r965 r978  
    44! Current revisions:
    55! -----------------
    6 !
     6! -km_damp_max, outflow_damping_width
     7! +pt_damping_factor, pt_damping_width
     8! +z0h_factor
    79!
    810! Former revisions:
     
    168170             inflow_damping_height, inflow_damping_width, &
    169171             inflow_disturbance_begin, inflow_disturbance_end, &
    170              initializing_actions, km_constant, km_damp_max, lad_surface, &
     172             initializing_actions, km_constant, lad_surface, &
    171173             lad_vertical_gradient, lad_vertical_gradient_level, &
    172174             leaf_surface_concentration, long_filter_factor, &
     
    175177             netcdf_precision, neutral, ngsrb, nsor, &
    176178             nsor_ini, nx, ny, nz, ocean, omega, omega_sor, &
    177              outflow_damping_width, overshoot_limit_e, overshoot_limit_pt, &
     179             overshoot_limit_e, overshoot_limit_pt, &
    178180             overshoot_limit_u, overshoot_limit_v, overshoot_limit_w, &
    179181             passive_scalar, pch_index, phi, plant_canopy, prandtl_layer, &
    180              prandtl_number, precipitation, psolver, pt_reference, pt_surface, &
     182             prandtl_number, precipitation, psolver, pt_damping_factor, &
     183             pt_damping_width, pt_reference, pt_surface, &
    181184             pt_surface_initial_change, pt_vertical_gradient, &
    182185             pt_vertical_gradient_level, q_surface, q_surface_initial_change, &
     
    200203             uv_heights, u_bulk, u_profile, vg_surface, vg_vertical_gradient, &
    201204             vg_vertical_gradient_level, v_bulk, v_profile, wall_adjustment, &
    202              wall_heatflux, wall_humidityflux, wall_scalarflux
     205             wall_heatflux, wall_humidityflux, wall_scalarflux, z0h_factor
    203206
    204207
  • palm/trunk/SOURCE/poismg.f90

    r881 r978  
    77! Current revisions:
    88! -----------------
    9 !
     9! bc_lr/ns_dirneu/neudir added
     10!
    1011!
    1112! Former revisions:
     
    12291230!--       outflow conditions have to be used on all PEs after the switch,
    12301231!--       because then they have the total domain.
    1231           IF ( bc_lr_dirrad )  THEN
     1232          IF ( bc_lr_dirrad .OR. bc_lr_dirneu )  THEN
    12321233             inflow_l  = .TRUE.
    12331234             inflow_r  = .FALSE.
    12341235             outflow_l = .FALSE.
    12351236             outflow_r = .TRUE.
    1236           ELSEIF ( bc_lr_raddir )  THEN
     1237          ELSEIF ( bc_lr_raddir .OR. bc_lr_neudir )  THEN
    12371238             inflow_l  = .FALSE.
    12381239             inflow_r  = .TRUE.
     
    12411242          ENDIF
    12421243
    1243           IF ( bc_ns_dirrad )  THEN
     1244          IF ( bc_ns_dirrad .OR. bc_ns_dirneu )  THEN
    12441245             inflow_n  = .TRUE.
    12451246             inflow_s  = .FALSE.
    12461247             outflow_n = .FALSE.
    12471248             outflow_s = .TRUE.
    1248           ELSEIF ( bc_ns_raddir )  THEN
     1249          ELSEIF ( bc_ns_raddir .OR. bc_ns_neudir )  THEN
    12491250             inflow_n  = .FALSE.
    12501251             inflow_s  = .TRUE.
     
    13161317
    13171318          IF ( pleft == MPI_PROC_NULL )  THEN
    1318              IF ( bc_lr_dirrad )  THEN
     1319             IF ( bc_lr_dirrad .OR. bc_lr_dirneu )  THEN
    13191320                inflow_l  = .TRUE.
    1320              ELSEIF ( bc_lr_raddir )  THEN
     1321             ELSEIF ( bc_lr_raddir .OR. bc_lr_neudir )  THEN
    13211322                outflow_l = .TRUE.
    13221323             ENDIF
     
    13241325
    13251326          IF ( pright == MPI_PROC_NULL )  THEN
    1326              IF ( bc_lr_dirrad )  THEN
     1327             IF ( bc_lr_dirrad .OR. bc_lr_dirneu )  THEN
    13271328                outflow_r = .TRUE.
    1328              ELSEIF ( bc_lr_raddir )  THEN
     1329             ELSEIF ( bc_lr_raddir .OR. bc_lr_neudir )  THEN
    13291330                inflow_r  = .TRUE.
    13301331             ENDIF
     
    13321333
    13331334          IF ( psouth == MPI_PROC_NULL )  THEN
    1334              IF ( bc_ns_dirrad )  THEN
     1335             IF ( bc_ns_dirrad .OR. bc_ns_dirneu )  THEN
    13351336                outflow_s = .TRUE.
    1336              ELSEIF ( bc_ns_raddir )  THEN
     1337             ELSEIF ( bc_ns_raddir .OR. bc_ns_neudir )  THEN
    13371338                inflow_s  = .TRUE.
    13381339             ENDIF
     
    13401341
    13411342          IF ( pnorth == MPI_PROC_NULL )  THEN
    1342              IF ( bc_ns_dirrad )  THEN
     1343             IF ( bc_ns_dirrad .OR. bc_ns_dirneu )  THEN
    13431344                inflow_n  = .TRUE.
    1344              ELSEIF ( bc_ns_raddir )  THEN
     1345             ELSEIF ( bc_ns_raddir .OR. bc_ns_neudir )  THEN
    13451346                outflow_n = .TRUE.
    13461347             ENDIF
  • palm/trunk/SOURCE/prandtl_fluxes.f90

    r760 r978  
    44! Current revisions:
    55! -----------------
    6 !
     6! roughness length for scalar quantities z0h added
    77!
    88! Former revisions:
     
    9696!
    9797!--             Stable stratification
    98                 ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (          &
    99                                   LOG( z_p / z0(j,i) ) +                   &
    100                                   5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
     98                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (           &
     99                                  LOG( z_p / z0h(j,i) ) +                   &
     100                                  5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
    101101                                                                )
    102102             ELSE
     
    104104!--             Unstable stratification
    105105                a = SQRT( 1.0 - 16.0 * rif(j,i) )
    106                 b = SQRT( 1.0 - 16.0 * rif(j,i) * z0(j,i) / z_p )
    107 
    108                 ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) /  (         &
    109                           LOG( z_p / z0(j,i) ) -                           &
     106                b = SQRT( 1.0 - 16.0 * rif(j,i) * z0h(j,i) / z_p )
     107
     108                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) /  (          &
     109                          LOG( z_p / z0h(j,i) ) -                           &
    110110                          2.0 * LOG( ( 1.0 + a ) / ( 1.0 + b ) ) )
    111111             ENDIF
     
    308308!
    309309!--                Stable stratification
    310                    qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (         &
    311                                   LOG( z_p / z0(j,i) ) +                   &
    312                                   5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
     310                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (          &
     311                                  LOG( z_p / z0h(j,i) ) +                   &
     312                                  5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
    313313                                                                 )
    314314                ELSE
     
    316316!--                Unstable stratification
    317317                   a = SQRT( 1.0 - 16.0 * rif(j,i) )
    318                    b = SQRT( 1.0 - 16.0 * rif(j,i) * z0(j,i) / z_p )
     318                   b = SQRT( 1.0 - 16.0 * rif(j,i) * z0h(j,i) / z_p )
    319319 
    320                    qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) /   (   &
    321                              LOG( z_p / z0(j,i) ) -                    &
     320                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) /   (        &
     321                             LOG( z_p / z0h(j,i) ) -                        &
    322322                              2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
    323323                ENDIF
  • palm/trunk/SOURCE/prognostic_equations.f90

    r941 r978  
    44! Current revisions:
    55! -----------------
    6 !
     6! km_damp_x and km_damp_y removed in calls of diffusion_u and diffusion_v
     7! add ptdf_x, ptdf_y for damping the potential temperature at the inflow
     8! boundary in case of non-cyclic lateral boundaries
     9! Bugfix: first thread index changes for WS-scheme at the inflow
    710!
    811! Former revisions:
     
    201204          ENDIF
    202205          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
    203              CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, u_m,  &
    204                                usws_m, uswst_m, v_m, w_m )
     206             CALL diffusion_u( i, j, ddzu, ddzw, km_m, tend, u_m, usws_m, &
     207                               uswst_m, v_m, w_m )
    205208          ELSE
    206              CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
    207                                uswst, v, w )
     209             CALL diffusion_u( i, j, ddzu, ddzw, km, tend, u, usws, uswst, &
     210                               v, w )
    208211          ENDIF
    209212          CALL coriolis( i, j, 1 )
     
    289292          ENDIF
    290293          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
    291              CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, u_m, &
    292                                v_m, vsws_m, vswst_m, w_m )
     294             CALL diffusion_v( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, &
     295                               vsws_m, vswst_m, w_m )
    293296          ELSE
    294              CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
    295                                vsws, vswst, w )
     297             CALL diffusion_v( i, j, ddzu, ddzw, km, tend, u, v, vsws, &
     298                               vswst, w )
    296299          ENDIF
    297300          CALL coriolis( i, j, 2 )
     
    373376          ENDIF
    374377          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
    375              CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, km_damp_y,  &
    376                                tend, u_m, v_m, w_m )
     378             CALL diffusion_w( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, w_m )
    377379          ELSE
    378              CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y,    &
    379                                tend, u, v, w )
     380             CALL diffusion_w( i, j, ddzu, ddzw, km, tend, u, v, w )
    380381          ENDIF
    381382          CALL coriolis( i, j, 3 )
     
    524525             DO  k = nzb_s_inner(j,i)+1, nzt
    525526                pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
    526                               dt_3d * (                                           &
     527                              dt_3d * (                                        &
    527528                                       sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
    528                                       ) -                                         &
    529                               tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) )
     529                                      ) - &
     530                              tsc(5) * ( pt(k,j,i) - pt_init(k) ) * &
     531                              ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) ) 
    530532             ENDDO
    531533
     
    983985             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
    984986                IF ( ws_scheme_mom )  THEN
    985                    IF ( outflow_l .AND. i_omp_start == nxl )  THEN
     987                   IF ( ( inflow_l .OR. outflow_l ) .AND. i_omp_start == nxl )  THEN
    986988!                     CALL local_diss( i, j, u)    ! dissipation control
    987989                      CALL advec_u_ws( i, j, i_omp_start + 1, tn )
     
    997999             IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
    9981000             THEN
    999                 CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, &
    1000                                   u_m, usws_m, uswst_m, v_m, w_m )
    1001              ELSE
    1002                 CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, &
    1003                                   usws, uswst, v, w )
     1001                CALL diffusion_u( i, j, ddzu, ddzw, km_m, tend,  u_m, &
     1002                                  usws_m, uswst_m, v_m, w_m )
     1003             ELSE
     1004                CALL diffusion_u( i, j, ddzu, ddzw, km, tend, u, usws, &
     1005                                  uswst, v, w )
    10041006             ENDIF
    10051007             CALL coriolis( i, j, 1 )
     
    10661068             IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
    10671069             THEN
    1068                 CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend,    &
    1069                                   u_m, v_m, vsws_m, vswst_m, w_m )
    1070              ELSE
    1071                 CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
    1072                                   vsws, vswst, w )
     1070                CALL diffusion_v( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, &
     1071                                  vsws_m, vswst_m, w_m )
     1072             ELSE
     1073                CALL diffusion_v( i, j, ddzu, ddzw, km, tend, u, v, vsws, &
     1074                                  vswst, w )
    10731075             ENDIF
    10741076             CALL coriolis( i, j, 2 )
     
    11301132          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
    11311133          THEN
    1132              CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x,          &
    1133                                km_damp_y, tend, u_m, v_m, w_m )
     1134             CALL diffusion_w( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, &
     1135                               w_m )
    11341136          ELSE
    1135              CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
    1136                                tend, u, v, w )
     1137             CALL diffusion_w( i, j, ddzu, ddzw, km, tend, u, v, w )
    11371138          ENDIF
    11381139          CALL coriolis( i, j, 3 )
     
    12401241                              dt_3d * (                                        &
    12411242                                  tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
    1242                                       ) -                                      &
    1243                               tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) )
     1243                                      ) - &
     1244                              tsc(5) * ( pt(k,j,i) - pt_init(k) ) * &
     1245                              ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) ) 
    12441246             ENDDO
    12451247
     
    15401542    ENDIF
    15411543    IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
    1542        CALL diffusion_u( ddzu, ddzw, km_m, km_damp_y, tend, u_m, usws_m, &
    1543                          uswst_m, v_m, w_m )
     1544       CALL diffusion_u( ddzu, ddzw, km_m, tend, u_m, usws_m, uswst_m, &
     1545                         v_m, w_m )
    15441546    ELSE
    1545        CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, uswst, v, w )
     1547       CALL diffusion_u( ddzu, ddzw, km, tend, u, usws, uswst, v, w )
    15461548    ENDIF
    15471549    CALL coriolis( 1 )
     
    16341636    ENDIF
    16351637    IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
    1636        CALL diffusion_v( ddzu, ddzw, km_m, km_damp_x, tend, u_m, v_m, vsws_m, &
    1637                          vswst_m, w_m )
     1638       CALL diffusion_v( ddzu, ddzw, km_m, tend, u_m, v_m, vsws_m, vswst_m, &
     1639                         w_m )
    16381640    ELSE
    1639        CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, vswst, w )
     1641       CALL diffusion_v( ddzu, ddzw, km, tend, u, v, vsws, vswst, w )
    16401642    ENDIF
    16411643    CALL coriolis( 2 )
     
    17251727    ENDIF
    17261728    IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
    1727        CALL diffusion_w( ddzu, ddzw, km_m, km_damp_x, km_damp_y, tend, u_m, &
    1728                          v_m, w_m )
     1729       CALL diffusion_w( ddzu, ddzw, km_m, tend, u_m, v_m, w_m )
    17291730    ELSE
    1730        CALL diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, w )
     1731       CALL diffusion_w( ddzu, ddzw, km, tend, u, v, w )
    17311732    ENDIF
    17321733    CALL coriolis( 3 )
     
    18771878          DO  j = nys, nyn
    18781879             DO  k = nzb_s_inner(j,i)+1, nzt
    1879                 pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) +       &
    1880                               dt_3d * (                                           &
    1881                                        sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
    1882                                       ) -                                         &
    1883                               tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) )
     1880                pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) +    &
     1881                              dt_3d * (                                        &
     1882                                       sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i)&
     1883                                      ) - &
     1884                              tsc(5) * ( pt(k,j,i) - pt_init(k) ) * &
     1885                              ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) ) 
    18841886             ENDDO
    18851887          ENDDO
  • palm/trunk/SOURCE/read_3d_binary.f90

    r777 r978  
    44! Current revisions:
    55! -----------------
     6! +z0h, z0h_av
    67!
    78! Former revisions:
     
    963964                       tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    964965
     966                CASE ( 'z0h' )
     967                   IF ( k == 1 )  READ ( 13 )  tmp_2d
     968                   z0h(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
     969                     tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     970
     971                CASE ( 'z0h_av' )
     972                   IF ( .NOT. ALLOCATED( z0h_av ) )  THEN
     973                      ALLOCATE( z0h_av(nysg:nyng,nxlg:nxrg) )
     974                   ENDIF
     975                   IF ( k == 1 )  READ ( 13 )  tmp_2d
     976                   z0h_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
     977                       tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     978
    965979                CASE DEFAULT
    966980                   WRITE( message_string, * ) 'unknown field named "', &
  • palm/trunk/SOURCE/read_var_list.f90

    r941 r978  
    44! Current revisions:
    55! ------------------
    6 !
     6! -km_damp_max, outflow_damping_width
     7! +pt_damping_factor, pt_damping_width
     8! +z0h_factor
    79!
    810! Former revisions:
     
    385387          CASE ( 'km_constant' )
    386388             READ ( 13 )  km_constant
    387           CASE ( 'km_damp_max' )
    388              READ ( 13 )  km_damp_max
    389389          CASE ( 'lad' )
    390390             READ ( 13 )  lad
     
    446446          CASE ( 'omega_sor' )
    447447             READ ( 13 )  omega_sor
    448           CASE ( 'outflow_damping_width' )
    449              READ ( 13 )  outflow_damping_width
    450448          CASE ( 'output_for_t0' )
    451449             READ (13)    output_for_t0
     
    476474          CASE ( 'psolver' )
    477475             READ ( 13 )  psolver
     476          CASE ( 'pt_damping_factor' )
     477             READ ( 13 )  pt_damping_factor
     478          CASE ( 'pt_damping_width' )
     479             READ ( 13 )  pt_damping_width
    478480          CASE ( 'pt_init' )
    479481             READ ( 13 )  pt_init
     
    684686          CASE ( 'w_max_ijk' )
    685687             READ ( 13 )  w_max_ijk
     688          CASE ( 'z0h_factor' )
     689             READ ( 13 )  z0h_factor
    686690
    687691          CASE DEFAULT
  • palm/trunk/SOURCE/sum_up_3d_data.f90

    r791 r978  
    44! Current revisions:
    55! -----------------
    6 !
     6! +z0h*
    77!
    88! Former revisions:
     
    227227                ENDIF
    228228                z0_av = 0.0
     229
     230             CASE ( 'z0h*' )
     231                IF ( .NOT. ALLOCATED( z0h_av ) )  THEN
     232                   ALLOCATE( z0h_av(nysg:nyng,nxlg:nxrg) )
     233                ENDIF
     234                z0h_av = 0.0
    229235
    230236             CASE DEFAULT
     
    494500             ENDDO
    495501
     502          CASE ( 'z0h*' )
     503             DO  i = nxlg, nxrg
     504                DO  j = nysg, nyng
     505                   z0h_av(j,i) = z0h_av(j,i) + z0h(j,i)
     506                ENDDO
     507             ENDDO
     508
    496509          CASE DEFAULT
    497510!
  • palm/trunk/SOURCE/timestep.f90

    r867 r978  
    44! Current revisions:
    55! ------------------
    6 !
     6! restriction of the outflow damping layer in the diffusion criterion removed
    77!
    88! Former revisions:
     
    185185       dt_diff = dt_diff_l
    186186#endif
    187 
    188 !
    189 !--    In case of non-cyclic lateral boundaries, the diffusion time step
    190 !--    may be further restricted by the lateral damping layer (damping only
    191 !--    along x and y)
    192        IF ( .NOT. bc_lr_cyc )  THEN
    193           dt_diff = MIN( dt_diff, 0.125 * dx2 / ( km_damp_max + 1E-20 ) )
    194        ELSEIF ( .NOT. bc_ns_cyc )  THEN
    195           dt_diff = MIN( dt_diff, 0.125 * dy2 / ( km_damp_max + 1E-20 ) )
    196        ENDIF
    197187
    198188!
  • palm/trunk/SOURCE/write_3d_binary.f90

    r791 r978  
    44! Current revisions:
    55! -----------------
    6 !
     6! +z0h, z0h_av
    77!
    88! Former revisions:
     
    298298       WRITE ( 14 )  'z0_av               ';  WRITE ( 14 )  z0_av
    299299    ENDIF
     300    WRITE ( 14 )  'z0h                 ';  WRITE ( 14 )  z0h
     301    IF ( ALLOCATED( z0h_av ) )  THEN
     302       WRITE ( 14 )  'z0h_av              ';  WRITE ( 14 )  z0h_av
     303    ENDIF
    300304
    301305!
  • palm/trunk/SOURCE/write_var_list.f90

    r941 r978  
    44! Current revisions:
    55! -----------------
    6 !
     6! -km_damp_max, outflow_damping_width
     7! +pt_damping_factor, pt_damping_width
     8! +z0h_factor
    79!
    810! Former revisions:
     
    300302    WRITE ( 14 )  'km_constant                   '
    301303    WRITE ( 14 )  km_constant
    302     WRITE ( 14 )  'km_damp_max                   '
    303     WRITE ( 14 )  km_damp_max
    304304    WRITE ( 14 )  'lad                           '
    305305    WRITE ( 14 )  lad
     
    358358    WRITE ( 14 )  'omega_sor                     '
    359359    WRITE ( 14 )  omega_sor
    360     WRITE ( 14 )  'outflow_damping_width         '
    361     WRITE ( 14 )  outflow_damping_width
    362360    WRITE ( 14 )  'output_for_t0                 '
    363361    WRITE ( 14 )  output_for_t0
     
    388386    WRITE ( 14 )  'psolver                       '
    389387    WRITE ( 14 )  psolver
     388    WRITE ( 14 )  'pt_damping_factor             '
     389    WRITE ( 14 )  pt_damping_factor
     390    WRITE ( 14 )  'pt_damping_width              '
     391    WRITE ( 14 )  pt_damping_width
    390392    WRITE ( 14 )  'pt_init                       '
    391393    WRITE ( 14 )  pt_init
     
    596598    WRITE ( 14 )  'w_max_ijk                     '
    597599    WRITE ( 14 )  w_max_ijk
     600    WRITE ( 14 )  'z0h_factor                    '
     601    WRITE ( 14 )  z0h_factor
    598602
    599603!
Note: See TracChangeset for help on using the changeset viewer.