Ignore:
Timestamp:
Mar 7, 2007 8:38:00 AM (15 years ago)
Author:
raasch
Message:

preliminary version, several changes to be explained later

File:
1 edited

Legend:

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

    r39 r51  
    44! Actual revisions:
    55! -----------------
    6 !
     6! wall functions now include diabatic conditions, call of routine wall_fluxes
    77!
    88! Former revisions:
     
    3030! ------------
    3131! Diffusion term of the u-component
     32! To do: additional damping (needed for non-cyclic bc) causes bad vectorization
     33!        and slows down the speed on NEC about 5-10%
    3234!------------------------------------------------------------------------------!
    3335
     
    5557
    5658       INTEGER ::  i, j, k
    57        REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp, usvs
     59       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
    5860       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
    5961       REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
    6062       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
     63       REAL, DIMENSION(nzb:nzt+1)      ::  usvs
    6164       REAL, DIMENSION(:,:),   POINTER ::  usws
    6265       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
     
    100103!--          Wall functions at the north and south walls, respectively
    101104             IF ( wall_u(j,i) /= 0.0 )  THEN
     105
     106!
     107!--             Calculate the horizontal momentum flux u'v'
     108                CALL wall_fluxes( i, j, nzb_u_inner(j,i)+1, nzb_u_outer(j,i),  &
     109                                  usvs, 1.0, 0.0, 0.0, 0.0 )
     110
    102111                DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
    103                    usvs   = kappa * u(k,j,i) / LOG( 0.5 * dy / z0(j,i) )
    104                    usvs   = -usvs * ABS( usvs )
    105112                   kmyp_x = 0.25 * &
    106113                            ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
     
    132139                                + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
    133140                                                  )                            &
    134                                      + wall_u(j,i) * usvs                      &
     141                                     + wall_u(j,i) * usvs(k)                   &
    135142                                   ) * ddy
    136143                ENDDO
     
    204211
    205212       INTEGER ::  i, j, k
    206        REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp, usvs
     213       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
    207214       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
    208215       REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
    209216       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
     217       REAL, DIMENSION(nzb:nzt+1)      ::  usvs
    210218       REAL, DIMENSION(:,:),   POINTER ::  usws
    211219       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
     
    246254!--    Wall functions at the north and south walls, respectively
    247255       IF ( wall_u(j,i) .NE. 0.0 )  THEN
     256
     257!
     258!--       Calculate the horizontal momentum flux u'v'
     259          CALL wall_fluxes( i, j, nzb_u_inner(j,i)+1, nzb_u_outer(j,i),  &
     260                            usvs, 1.0, 0.0, 0.0, 0.0 )
     261
    248262          DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
    249              usvs   = kappa * u(k,j,i) / LOG( 0.5 * dy / z0(j,i) )
    250              usvs   = -usvs * ABS( usvs )
    251263             kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
    252264             kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
     
    276288                                + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
    277289                                                  )                            &
    278                                      + wall_u(j,i) * usvs                      &
     290                                     + wall_u(j,i) * usvs(k)                   &
    279291                                   ) * ddy
    280292          ENDDO
Note: See TracChangeset for help on using the changeset viewer.