Changeset 56


Ignore:
Timestamp:
Mar 8, 2007 1:57:07 PM (15 years ago)
Author:
raasch
Message:

further checkin of preliminary changes

Location:
palm/trunk/SOURCE
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile

    r54 r56  
    153153diffusion_e.o: modules.o
    154154diffusion_s.o: modules.o
    155 diffusion_u.o: modules.o
    156 diffusion_v.o: modules.o
    157 diffusion_w.o: modules.o
     155diffusion_u.o: modules.o wall_fluxes.o
     156diffusion_v.o: modules.o wall_fluxes.o
     157diffusion_w.o: modules.o wall_fluxes.o
    158158diffusivities.o: modules.o
    159159disturb_field.o: modules.o random_function.o
     
    193193pres.o: modules.o poisfft.o poisfft_hybrid.o
    194194print_1d.o: modules.o
    195 production_e.o: modules.o
     195production_e.o: modules.o wall_fluxes.o
    196196prognostic_equations.o: modules.o advec_s_pw.o advec_s_up.o advec_u_pw.o \
    197197        advec_u_up.o advec_v_pw.o advec_v_up.o advec_w_pw.o advec_w_up.o  \
  • palm/trunk/SOURCE/diffusion_u.f90

    r53 r56  
    3434!------------------------------------------------------------------------------!
    3535
     36    USE wall_fluxes_mod
     37
    3638    PRIVATE
    3739    PUBLIC diffusion_u
     
    6163       REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
    6264       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
    63        REAL, DIMENSION(nzb:nzt+1)      ::  usvs
    6465       REAL, DIMENSION(:,:),   POINTER ::  usws
    6566       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
     67       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr+uxrp) ::  usvs
     68
     69!
     70!--    First calculate horizontal momentum flux u'v' at vertical walls,
     71!--    if neccessary
     72       IF ( topography /= 'flat' )  THEN
     73          CALL wall_fluxes( usvs, 1.0, 0.0, 0.0, 0.0, uxrp, 0, nzb_u_inner, &
     74                            nzb_u_outer, wall_u )
     75       ENDIF
    6676
    6777       DO  i = nxl, nxr+uxrp
     
    103113!--          Wall functions at the north and south walls, respectively
    104114             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 )
    110115
    111116                DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
     
    139144                                + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
    140145                                                  )                            &
    141                                      + wall_u(j,i) * usvs(k)                   &
     146                                     + wall_u(j,i) * usvs(k,j,i)               &
    142147                                   ) * ddy
    143148                ENDDO
  • palm/trunk/SOURCE/diffusion_v.f90

    r53 r56  
    3232!------------------------------------------------------------------------------!
    3333
     34    USE wall_fluxes_mod
     35
    3436    PRIVATE
    3537    PUBLIC diffusion_v
     
    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)
    61        REAL, DIMENSION(nzb:nzt+1)      ::  vsus
    6263       REAL, DIMENSION(:,:),   POINTER ::  vsws
    6364       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
     65       REAL, DIMENSION(nzb:nzt+1,nys:nyn+vynp,nxl:nxr) ::  vsus
     66
     67!
     68!--    First calculate horizontal momentum flux v'u' at vertical walls,
     69!--    if neccessary
     70       IF ( topography /= 'flat' )  THEN
     71          CALL wall_fluxes( vsus, 0.0, 1.0, 0.0, 0.0, 0, vynp, nzb_v_inner, &
     72                            nzb_v_outer, wall_v )
     73       ENDIF
    6474
    6575       DO  i = nxl, nxr
     
    101111!--          Wall functions at the left and right walls, respectively
    102112             IF ( wall_v(j,i) /= 0.0 )  THEN
    103 
    104 !
    105 !--             Calculate the horizontal momentum flux v'u'
    106                 CALL wall_fluxes( i, j, nzb_v_inner(j,i)+1, nzb_v_outer(j,i), &
    107                                   vsus, 0.0, 1.0, 0.0, 0.0 )
    108113
    109114                DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
     
    137142                                + kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
    138143                                                  )                            &
    139                                      + wall_v(j,i) * vsus(k)                   &
     144                                     + wall_v(j,i) * vsus(k,j,i)               &
    140145                                   ) * ddx
    141146                ENDDO
  • palm/trunk/SOURCE/diffusion_w.f90

    r53 r56  
    2828! Diffusion term of the w-component
    2929!------------------------------------------------------------------------------!
     30
     31    USE wall_fluxes_mod
    3032
    3133    PRIVATE
     
    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)
    61        REAL, DIMENSION(nzb:nzt+1)      ::  wsus, wsvs
    6263       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
    63 
     64       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus, wsvs
     65
     66
     67!
     68!--    First calculate horizontal momentum flux w'u' and/or w'v' at vertical
     69!--    walls, if neccessary
     70       IF ( topography /= 'flat' )  THEN
     71          CALL wall_fluxes( wsus, 0.0, 0.0, 0.0, 1.0, 0, 0, nzb_w_inner, &
     72                            nzb_w_outer, wall_w_x )
     73          CALL wall_fluxes( wsvs, 0.0, 0.0, 1.0, 0.0, 0, 0, nzb_w_inner, &
     74                            nzb_w_outer, wall_w_y )
     75       ENDIF
    6476
    6577       DO  i = nxl, nxr
     
    115127             IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
    116128
    117 !
    118 !--             Calculate the horizontal momentum fluxes w'u' and/or w'v'
    119                 IF ( wall_w_x(j,i) /= 0.0 )  THEN
    120                    CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1,              &
    121                                      nzb_w_outer(j,i), wsus, 0.0, 0.0, 0.0, &
    122                                      1.0 )
    123                 ELSE
    124                    wsus = 0.0
    125                 ENDIF
    126 
    127                 IF ( wall_w_y(j,i) /= 0.0 )  THEN
    128                    CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1,              &
    129                                      nzb_w_outer(j,i), wsvs, 0.0, 0.0, 1.0, &
    130                                      0.0 )
    131                 ELSE
    132                    wsvs = 0.0
    133                 ENDIF
    134 
    135129                DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
    136130!
     
    171165                          + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
    172166                                                   )                           &
    173                                      + wall_w_x(j,i) * wsus(k)                 &
     167                                     + wall_w_x(j,i) * wsus(k,j,i)             &
    174168                                   ) * ddx                                     &
    175169                                 + (   fwyp(j,i) * (                           &
     
    181175                          + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
    182176                                                   )                           &
    183                                      + wall_w_y(j,i) * wsvs(k)                 &
     177                                     + wall_w_y(j,i) * wsvs(k,j,i)             &
    184178                                   ) * ddy                                     &
    185179                                 + 2.0 * (                                     &
  • palm/trunk/SOURCE/production_e.f90

    r55 r56  
    44! Actual revisions:
    55! -----------------
    6 ! Wall functions now include diabatic conditions, call of routine wall_fluxes
     6! Wall functions now include diabatic conditions, call of routine wall_fluxes_e
    77!
    88! Former revisions:
     
    3131!------------------------------------------------------------------------------!
    3232
     33    USE wall_fluxes_mod
     34
    3335    PRIVATE
    3436    PUBLIC production_e, production_e_init
    35    
     37
    3638    LOGICAL, SAVE ::  first_call = .TRUE.
    3739
     
    6971                   k1, k2, theta, temp
    7072
    71        REAL, DIMENSION(nzb:nzt+1) ::  usvs, vsus, wsus, wsvs
    72 
     73!       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs, vsus, wsus, wsvs
     74       REAL, DIMENSION(nzb:nzt+1) ::   usvs, vsus, wsus, wsvs
     75
     76!
     77!--    First calculate horizontal momentum flux u'v', w'v', v'u', w'u' at
     78!--    vertical walls, if neccessary
     79!--    So far, results are slightly different from the ij-Version.
     80!--    Therefore, ij-Version is called further below within the ij-loops.
     81!       IF ( topography /= 'flat' )  THEN
     82!          CALL wall_fluxes_e( usvs, 1.0, 0.0, 0.0, 0.0, wall_e_y )
     83!          CALL wall_fluxes_e( wsvs, 0.0, 0.0, 1.0, 0.0, wall_e_y )
     84!          CALL wall_fluxes_e( vsus, 0.0, 1.0, 0.0, 0.0, wall_e_x )
     85!          CALL wall_fluxes_e( wsus, 0.0, 0.0, 0.0, 1.0, wall_e_x )
     86!       ENDIF
    7387
    7488!
     
    131145                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
    132146                                          usvs, 1.0, 0.0, 0.0, 0.0 )
    133 
    134147                      dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
     148!                      dudy = wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
    135149                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
    136150                                          wsvs, 0.0, 0.0, 1.0, 0.0 )
    137151                      dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
     152!                      dwdy = wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
    138153                   ELSE
    139154                      dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     
    147162                                          vsus, 0.0, 1.0, 0.0, 0.0 )
    148163                      dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
     164!                      dvdx = wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
    149165                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
    150166                                          wsus, 0.0, 0.0, 0.0, 1.0 )
    151167                      dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
     168!                      dwdx = wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
    152169                   ELSE
    153170                      dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     
    185202                      IF ( wall_e_y(j,i) /= 0.0 )  THEN
    186203                         dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
     204!                         dudy = wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
    187205                         dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
     206!                         dwdy = wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
    188207                      ELSE
    189208                         dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     
    195214                      IF ( wall_e_x(j,i) /= 0.0 )  THEN
    196215                         dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
     216!                         dvdx = wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
    197217                         dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
     218!                         dwdx = wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
    198219                      ELSE
    199220                         dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     
    478499                   k1, k2, theta, temp
    479500
    480        REAL, DIMENSION(nzb:nzt+1) ::  usvs, vsus, wsus,wsvs
     501       REAL, DIMENSION(nzb:nzt+1) ::  usvs, vsus, wsus, wsvs
    481502
    482503!
  • palm/trunk/SOURCE/wall_fluxes.f90

    r55 r56  
    1  SUBROUTINE wall_fluxes( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
     1 MODULE wall_fluxes_mod
    22!------------------------------------------------------------------------------!
    33! Actual revisions:
     
    1515! similarity.
    1616! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
    17 !------------------------------------------------------------------------------!
    18 
    19     USE arrays_3d
    20     USE control_parameters
    21     USE grid_variables
    22     USE indices
    23     USE statistics
    24     USE user
    25 
    26     IMPLICIT NONE
    27 
    28     INTEGER ::  i, j, k, nzb_w, nzt_w, wall_index
    29     REAL    ::  a, b, c1, c2, h1, h2, zp
    30 
    31     REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
    32 
    33     REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
    34 
    35 
    36     zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
    37     wall_flux  = 0.0
    38     wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
    39 
    40 !
    41 !-- All subsequent variables are computed for the respective location where
    42 !-- the relevant variable is defined
    43     DO  k = nzb_w, nzt_w
    44 
    45 !
    46 !--    (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
    47        rifs  = rif_wall(k,j,i,wall_index)
    48 
    49        u_i   = a * u(k,j,i) + c1 * 0.25 * &
    50                ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
    51 
    52        v_i   = b * v(k,j,i) + c2 * 0.25 * &
    53                ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
    54 
    55        ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                             &
    56                    a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
    57                  + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
    58                                                )
    59        pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + b * pt(k,j-1,i)  &
    60                        + ( c1 + c2 ) * pt(k+1,j,i) )
    61 
    62        pts   = pt_i - hom(k,1,4,0)
    63        wspts = ws * pts
    64 
    65 !
    66 !--    (2) Compute wall-parallel absolute velocity vel_total
    67        vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 )
    68 
    69 !
    70 !--    (3) Compute wall friction velocity us_wall
    71        IF ( rifs >= 0.0 )  THEN
    72 
    73 !
    74 !--       Stable stratification (and neutral)
    75           us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +               &
    76                                           5.0 * rifs * ( zp - z0(j,i) ) / zp  &
    77                                         )
    78        ELSE
    79 
    80 !
    81 !--       Unstable stratification
    82           h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
    83           h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
    84 
    85 !
    86 !--       If a borderline case occurs, the formula for stable stratification
    87 !--       must be used anyway, or else a zero division would occur in the
    88 !--       argument of the logarithm.
    89           IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
     17! The all-gridpoint version of wall_fluxes_e is not used so far, because
     18! it gives slightly different results from the ij-version for some unknown
     19! reason.
     20!------------------------------------------------------------------------------!
     21    PRIVATE
     22    PUBLIC wall_fluxes, wall_fluxes_e
     23   
     24    INTERFACE wall_fluxes
     25       MODULE PROCEDURE wall_fluxes
     26       MODULE PROCEDURE wall_fluxes_ij
     27    END INTERFACE wall_fluxes
     28   
     29    INTERFACE wall_fluxes_e
     30       MODULE PROCEDURE wall_fluxes_e
     31       MODULE PROCEDURE wall_fluxes_e_ij
     32    END INTERFACE wall_fluxes_e
     33 
     34 CONTAINS
     35
     36!------------------------------------------------------------------------------!
     37! Call for all grid points
     38!------------------------------------------------------------------------------!
     39    SUBROUTINE wall_fluxes( wall_flux, a, b, c1, c2, ixp, jyp, nzb_uvw_inner, &
     40                            nzb_uvw_outer, wall )
     41
     42       USE arrays_3d
     43       USE control_parameters
     44       USE grid_variables
     45       USE indices
     46       USE statistics
     47       USE user
     48
     49       IMPLICIT NONE
     50
     51       INTEGER ::  i, ixp, j, jyp, k, wall_index
     52
     53       INTEGER, DIMENSION(nys-1:nyn+1,nxl-1:nxr+1) ::  nzb_uvw_inner, &
     54                                                       nzb_uvw_outer
     55       REAL ::  a, b, c1, c2, h1, h2, zp
     56       REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
     57
     58       REAL, DIMENSION(nys-1:nyn+1,nxl-1:nxr+1)           ::  wall
     59       REAL, DIMENSION(nzb:nzt+1,nys:nyn+jyp,nxl:nxr+ixp) ::  wall_flux
     60
     61
     62       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
     63       wall_flux  = 0.0
     64       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
     65
     66       DO  i = nxl, nxr+ixp
     67          DO  j = nys, nyn+jyp
     68
     69             IF ( wall(j,i) /= 0.0 )  THEN
     70!
     71!--             All subsequent variables are computed for the respective
     72!--             location where the relevant variable is defined
     73                DO  k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
     74
     75!
     76!--                (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
     77                   rifs  = rif_wall(k,j,i,wall_index)
     78
     79                   u_i   = a * u(k,j,i) + c1 * 0.25 * &
     80                           ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
     81
     82                   v_i   = b * v(k,j,i) + c2 * 0.25 * &
     83                           ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
     84
     85                   ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                   &
     86                     a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
     87                   + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
     88                                                           )
     89                   pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + &
     90                                   b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
     91
     92                   pts   = pt_i - hom(k,1,4,0)
     93                   wspts = ws * pts
     94
     95!
     96!--                (2) Compute wall-parallel absolute velocity vel_total
     97                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
     98
     99!
     100!--                (3) Compute wall friction velocity us_wall
     101                   IF ( rifs >= 0.0 )  THEN
     102
     103!
     104!--                   Stable stratification (and neutral)
     105                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
     106                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
     107                                                    )
     108                   ELSE
     109
     110!
     111!--                   Unstable stratification
     112                      h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     113                      h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ))
     114
     115!
     116!--                   If a borderline case occurs, the formula for stable
     117!--                   stratification must be used anyway, or else a zero
     118!--                   division would occur in the argument of the logarithm.
     119                      IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
     120                         us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
     121                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
     122                                                       )
     123                      ELSE
     124                         us_wall = kappa * vel_total / (                       &
     125                            LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + &
     126                            2.0 * ( ATAN( h2 ) - ATAN( h1 ) )                  &
     127                                                       )
     128                      ENDIF
     129
     130                   ENDIF
     131
     132!
     133!--                (4) Compute zp/L (corresponds to neutral Richardson flux
     134!--                    number rifs)
     135                   rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * &
     136                                                        ( us_wall**3 + 1E-30 ) )
     137
     138!
     139!--                Limit the value range of the Richardson numbers.
     140!--                This is necessary for very small velocities (u,w --> 0),
     141!--                because the absolute value of rif can then become very
     142!--                large, which in consequence would result in very large
     143!--                shear stresses and very small momentum fluxes (both are
     144!--                generally unrealistic).
     145                   IF ( rifs < rif_min )  rifs = rif_min
     146                   IF ( rifs > rif_max )  rifs = rif_max
     147
     148!
     149!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
     150                   IF ( rifs >= 0.0 )  THEN
     151
     152!
     153!--                   Stable stratification (and neutral)
     154                      wall_flux(k,j,i) = kappa *                               &
     155                              ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
     156                              (  LOG( zp / z0(j,i) ) +                         &
     157                                 5.0 * rifs * ( zp - z0(j,i) ) / zp            &
     158                              )
     159                   ELSE
     160
     161!
     162!--                   Unstable stratification
     163                      h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     164                      h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ))
     165
     166!
     167!--                   If a borderline case occurs, the formula for stable
     168!--                   stratification must be used anyway, or else a zero
     169!--                   division would occur in the argument of the logarithm.
     170                      IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
     171                         wall_flux(k,j,i) = kappa *                            &
     172                              ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
     173                              ( LOG( zp / z0(j,i) ) +                          &
     174                                5.0 * rifs * ( zp - z0(j,i) ) / zp             &
     175                              )
     176                      ELSE
     177                         wall_flux(k,j,i) = kappa *                            &
     178                             ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) /  &
     179                             ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) )&
     180                               + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )             &
     181                             )
     182                      ENDIF
     183
     184                   ENDIF
     185                   wall_flux(k,j,i) = -wall_flux(k,j,i) * ABS(wall_flux(k,j,i))
     186
     187!
     188!--                store rifs for next time step
     189                   rif_wall(k,j,i,wall_index) = rifs
     190
     191                ENDDO
     192
     193             ENDIF
     194
     195          ENDDO
     196       ENDDO
     197
     198    END SUBROUTINE wall_fluxes
     199
     200
     201
     202!------------------------------------------------------------------------------!
     203! Call for all grid point i,j
     204!------------------------------------------------------------------------------!
     205    SUBROUTINE wall_fluxes_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
     206
     207       USE arrays_3d
     208       USE control_parameters
     209       USE grid_variables
     210       USE indices
     211       USE statistics
     212       USE user
     213
     214       IMPLICIT NONE
     215
     216       INTEGER ::  i, j, k, nzb_w, nzt_w, wall_index
     217       REAL    ::  a, b, c1, c2, h1, h2, zp
     218
     219       REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
     220
     221       REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
     222
     223
     224       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
     225       wall_flux  = 0.0
     226       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
     227
     228!
     229!--    All subsequent variables are computed for the respective location where
     230!--    the relevant variable is defined
     231       DO  k = nzb_w, nzt_w
     232
     233!
     234!--       (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
     235          rifs  = rif_wall(k,j,i,wall_index)
     236
     237          u_i   = a * u(k,j,i) + c1 * 0.25 * &
     238                  ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
     239
     240          v_i   = b * v(k,j,i) + c2 * 0.25 * &
     241                  ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
     242
     243          ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                            &
     244                     a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
     245                   + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
     246                                                  )
     247          pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + b * pt(k,j-1,i)  &
     248                          + ( c1 + c2 ) * pt(k+1,j,i) )
     249
     250          pts   = pt_i - hom(k,1,4,0)
     251          wspts = ws * pts
     252
     253!
     254!--       (2) Compute wall-parallel absolute velocity vel_total
     255          vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 )
     256
     257!
     258!--       (3) Compute wall friction velocity us_wall
     259          IF ( rifs >= 0.0 )  THEN
     260
     261!
     262!--          Stable stratification (and neutral)
    90263             us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             &
    91264                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
    92265                                           )
    93266          ELSE
    94              us_wall = kappa * vel_total / (                                   &
     267
     268!
     269!--          Unstable stratification
     270             h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     271             h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
     272
     273!
     274!--          If a borderline case occurs, the formula for stable stratification
     275!--          must be used anyway, or else a zero division would occur in the
     276!--          argument of the logarithm.
     277             IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
     278                us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +          &
     279                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
     280                                              )
     281             ELSE
     282                us_wall = kappa * vel_total / (                                &
    95283                            LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + &
    96284                            2.0 * ( ATAN( h2 ) - ATAN( h1 ) )                  &
    97                                            )
     285                                              )
     286             ENDIF
     287
    98288          ENDIF
    99289
    100        ENDIF
    101 
    102 !
    103 !--    (4) Compute zp/L (corresponds to neutral Richardson flux number
    104 !--        rifs)
    105        rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * ( us_wall**3 + 1E-30 ) )
    106 
    107 !
    108 !--    Limit the value range of the Richardson numbers.
    109 !--    This is necessary for very small velocities (u,w --> 0), because
    110 !--    the absolute value of rif can then become very large, which in
    111 !--    consequence would result in very large shear stresses and very
    112 !--    small momentum fluxes (both are generally unrealistic).
    113        IF ( rifs < rif_min )  rifs = rif_min
    114        IF ( rifs > rif_max )  rifs = rif_max
    115 
    116 !
    117 !--    (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
    118        IF ( rifs >= 0.0 )  THEN
    119 
    120 !
    121 !--       Stable stratification (and neutral)
    122           wall_flux(k) = kappa *                                          &
    123                          ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
    124                          (  LOG( zp / z0(j,i) ) +                         &
    125                             5.0 * rifs * ( zp - z0(j,i) ) / zp            &
    126                          )
    127        ELSE
    128 
    129 !
    130 !--       Unstable stratification
    131           h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
    132           h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
    133 
    134 !
    135 !--       If a borderline case occurs, the formula for stable stratification
    136 !--       must be used anyway, or else a zero division would occur in the
    137 !--       argument of the logarithm.
    138           IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
     290!
     291!--       (4) Compute zp/L (corresponds to neutral Richardson flux number
     292!--           rifs)
     293          rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * (us_wall**3 + 1E-30) )
     294
     295!
     296!--       Limit the value range of the Richardson numbers.
     297!--       This is necessary for very small velocities (u,w --> 0), because
     298!--       the absolute value of rif can then become very large, which in
     299!--       consequence would result in very large shear stresses and very
     300!--       small momentum fluxes (both are generally unrealistic).
     301          IF ( rifs < rif_min )  rifs = rif_min
     302          IF ( rifs > rif_max )  rifs = rif_max
     303
     304!
     305!--       (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
     306          IF ( rifs >= 0.0 )  THEN
     307
     308!
     309!--          Stable stratification (and neutral)
    139310             wall_flux(k) = kappa *                                          &
    140311                            ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
    141                             ( LOG( zp / z0(j,i) ) +                          &
    142                               5.0 * rifs * ( zp - z0(j,i) ) / zp             &
     312                            (  LOG( zp / z0(j,i) ) +                         &
     313                               5.0 * rifs * ( zp - z0(j,i) ) / zp            &
    143314                            )
    144315          ELSE
    145              wall_flux(k) = kappa *                                            &
    146                             ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) /   &
    147                             ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) &
    148                               + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )              &
    149                             )
     316
     317!
     318!--          Unstable stratification
     319             h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     320             h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
     321
     322!
     323!--          If a borderline case occurs, the formula for stable stratification
     324!--          must be used anyway, or else a zero division would occur in the
     325!--          argument of the logarithm.
     326             IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
     327                wall_flux(k) = kappa *                                         &
     328                              ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
     329                              ( LOG( zp / z0(j,i) ) +                          &
     330                                5.0 * rifs * ( zp - z0(j,i) ) / zp             &
     331                              )
     332             ELSE
     333                wall_flux(k) = kappa *                                         &
     334                             ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) /  &
     335                             ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) )&
     336                               + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )             &
     337                             )
     338             ENDIF
     339
    150340          ENDIF
    151 
    152        ENDIF
    153        wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) )
    154 
    155 !
    156 !--    store rifs for next time step
    157        rif_wall(k,j,i,wall_index) = rifs
    158 
    159     ENDDO
    160 
    161  END SUBROUTINE wall_fluxes
    162 
    163 
    164 
    165  SUBROUTINE wall_fluxes_e( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
    166 !------------------------------------------------------------------------------!
    167 ! Actual revisions:
    168 ! -----------------
    169 !
    170 !
    171 ! Former revisions:
    172 ! -----------------
    173 ! Initial version (2007/03/07)
    174 !
     341          wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) )
     342
     343!
     344!--       store rifs for next time step
     345          rif_wall(k,j,i,wall_index) = rifs
     346
     347       ENDDO
     348
     349    END SUBROUTINE wall_fluxes_ij
     350
     351
     352
     353!------------------------------------------------------------------------------!
     354! Call for all grid points
     355!------------------------------------------------------------------------------!
     356    SUBROUTINE wall_fluxes_e( wall_flux, a, b, c1, c2, wall )
     357
     358!------------------------------------------------------------------------------!
    175359! Description:
    176360! ------------
     
    180364!------------------------------------------------------------------------------!
    181365
    182     USE arrays_3d
    183     USE control_parameters
    184     USE grid_variables
    185     USE indices
    186     USE statistics
    187     USE user
    188 
    189     IMPLICIT NONE
    190 
    191     INTEGER ::  i, j, k, kk, nzb_w, nzt_w, wall_index
    192     REAL    ::  a, b, c1, c2, h1, h2, vel_zp, zp
    193 
    194     REAL ::  rifs
    195 
    196     REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
    197 
    198 
    199     zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
    200     wall_flux  = 0.0
    201     wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
    202 
    203 !
    204 !-- All subsequent variables are computed for the respective location where
    205 !-- the relevant variable is defined
    206     DO  k = nzb_w, nzt_w
    207 
    208 !
    209 !--    (1) Compute rifs
    210        IF ( k == nzb_w )  THEN
    211           kk = nzb_w
    212        ELSE
    213           kk = k-1
    214        ENDIF
    215        rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                         &
    216                        a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
    217                        c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
    218                      )
    219 
    220 !
    221 !--    Skip (2) to (4) of wall_fluxes, because here rifs is already available
    222 !--    from (1)
    223 
    224 !
    225 !--    (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
    226        vel_zp = 0.5 * (       a * ( u(k,j,i) + u(k,j,i+1) ) +  &
    227                               b * ( v(k,j,i) + v(k,j+1,i) ) +  &
    228                         (c1+c2) * ( w(k,j,i) + w(k-1,j,i) )    &
    229                       )
    230 
    231        IF ( rifs >= 0.0 )  THEN
    232 
    233 !
    234 !--       Stable stratification (and neutral)
    235           wall_flux(k) = kappa *  vel_zp / &
    236                          ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
    237        ELSE
    238 
    239 !
    240 !--       Unstable stratification
    241           h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
    242           h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
    243 
    244 !
    245 !--       If a borderline case occurs, the formula for stable stratification
    246 !--       must be used anyway, or else a zero division would occur in the
    247 !--       argument of the logarithm.
    248           IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
    249              wall_flux(k) = kappa * vel_zp /                                 &
    250                             ( LOG( zp / z0(j,i) ) +                          &
    251                               5.0 * rifs * ( zp - z0(j,i) ) / zp             &
    252                             )
    253           ELSE
    254              wall_flux(k) = kappa * vel_zp /                                   &
     366       USE arrays_3d
     367       USE control_parameters
     368       USE grid_variables
     369       USE indices
     370       USE statistics
     371       USE user
     372
     373       IMPLICIT NONE
     374
     375       INTEGER ::  i, j, k, kk, wall_index
     376       REAL    ::  a, b, c1, c2, h1, h2, vel_zp, zp
     377
     378       REAL ::  rifs
     379
     380       REAL, DIMENSION(nys-1:nyn+1,nxl-1:nxr+1)   ::  wall
     381       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
     382
     383
     384       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
     385       wall_flux  = 0.0
     386       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
     387
     388       DO  i = nxl, nxr
     389          DO  j = nys, nyn
     390
     391             IF ( wall(j,i) /= 0.0 )  THEN
     392!
     393!--             All subsequent variables are computed for the respective
     394!--             location where the relevant variable is defined
     395                DO  k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2
     396
     397!
     398!--                (1) Compute rifs
     399                   IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
     400                      kk = nzb_diff_s_inner(j,i)-1
     401                   ELSE
     402                      kk = k-1
     403                   ENDIF
     404                   rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                &
     405                          a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
     406                          c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
     407                                 )
     408
     409!
     410!--                Skip (2) to (4) of wall_fluxes, because here rifs is
     411!--                already available from (1)
     412
     413!
     414!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
     415                   vel_zp = 0.5 * (       a * ( u(k,j,i) + u(k,j,i+1) ) +  &
     416                                          b * ( v(k,j,i) + v(k,j+1,i) ) +  &
     417                                    (c1+c2) * ( w(k,j,i) + w(k-1,j,i) )    &
     418                                  )
     419
     420                   IF ( rifs >= 0.0 )  THEN
     421
     422!
     423!--                   Stable stratification (and neutral)
     424                      wall_flux(k,j,i) = kappa *  vel_zp / &
     425                          ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
     426                   ELSE
     427
     428!
     429!--                   Unstable stratification
     430                      h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     431                      h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ))
     432
     433!
     434!--                   If a borderline case occurs, the formula for stable
     435!--                   stratification must be used anyway, or else a zero
     436!--                   division would occur in the argument of the logarithm.
     437                      IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
     438                         wall_flux(k,j,i) = kappa * vel_zp /                 &
     439                                        ( LOG( zp / z0(j,i) ) +              &
     440                                          5.0 * rifs * ( zp - z0(j,i) ) / zp &
     441                                        )
     442                      ELSE
     443                         wall_flux(k,j,i) = kappa * vel_zp /                   &
    255444                            ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) &
    256445                              + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )              &
    257446                            )
     447                      ENDIF
     448
     449                   ENDIF
     450                   wall_flux(k,j,i) = wall_flux(k,j,i) * ABS( wall_flux(k,j,i) )
     451
     452!
     453!--                Store rifs for next time step
     454                   rif_wall(k,j,i,wall_index) = rifs
     455
     456                ENDDO
     457
     458             ENDIF
     459
     460          ENDDO
     461       ENDDO
     462
     463    END SUBROUTINE wall_fluxes_e
     464
     465
     466
     467!------------------------------------------------------------------------------!
     468! Call for grid point i,j
     469!------------------------------------------------------------------------------!
     470    SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
     471
     472       USE arrays_3d
     473       USE control_parameters
     474       USE grid_variables
     475       USE indices
     476       USE statistics
     477       USE user
     478
     479       IMPLICIT NONE
     480
     481       INTEGER ::  i, j, k, kk, nzb_w, nzt_w, wall_index
     482       REAL    ::  a, b, c1, c2, h1, h2, vel_zp, zp
     483
     484       REAL ::  rifs
     485
     486       REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
     487
     488
     489       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
     490       wall_flux  = 0.0
     491       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
     492
     493!
     494!--    All subsequent variables are computed for the respective location where
     495!--    the relevant variable is defined
     496       DO  k = nzb_w, nzt_w
     497
     498!
     499!--       (1) Compute rifs
     500          IF ( k == nzb_w )  THEN
     501             kk = nzb_w
     502          ELSE
     503             kk = k-1
    258504          ENDIF
    259 
    260        ENDIF
    261        wall_flux(k) = wall_flux(k) * ABS( wall_flux(k) )
    262 
    263 !
    264 !--    store rifs for next time step
    265        rif_wall(k,j,i,wall_index) = rifs
    266 
    267     ENDDO
    268 
    269  END SUBROUTINE wall_fluxes_e
     505          rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                         &
     506                          a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
     507                          c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
     508                        )
     509
     510!
     511!--       Skip (2) to (4) of wall_fluxes, because here rifs is already available
     512!--       from (1)
     513
     514!
     515!--       (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
     516          vel_zp = 0.5 * (       a * ( u(k,j,i) + u(k,j,i+1) ) +  &
     517                                 b * ( v(k,j,i) + v(k,j+1,i) ) +  &
     518                           (c1+c2) * ( w(k,j,i) + w(k-1,j,i) )    &
     519                         )
     520
     521          IF ( rifs >= 0.0 )  THEN
     522
     523!
     524!--          Stable stratification (and neutral)
     525             wall_flux(k) = kappa *  vel_zp / &
     526                          ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
     527          ELSE
     528
     529!
     530!--          Unstable stratification
     531             h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     532             h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
     533
     534!
     535!--          If a borderline case occurs, the formula for stable stratification
     536!--          must be used anyway, or else a zero division would occur in the
     537!--          argument of the logarithm.
     538             IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
     539                wall_flux(k) = kappa * vel_zp /                     &
     540                               ( LOG( zp / z0(j,i) ) +              &
     541                                 5.0 * rifs * ( zp - z0(j,i) ) / zp &
     542                               )
     543             ELSE
     544                wall_flux(k) = kappa * vel_zp /                                &
     545                            ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) &
     546                              + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )              &
     547                            )
     548             ENDIF
     549
     550          ENDIF
     551          wall_flux(k) = wall_flux(k) * ABS( wall_flux(k) )
     552
     553!
     554!--       Store rifs for next time step
     555          rif_wall(k,j,i,wall_index) = rifs
     556
     557       ENDDO
     558
     559    END SUBROUTINE wall_fluxes_e_ij
     560
     561 END MODULE wall_fluxes_mod
Note: See TracChangeset for help on using the changeset viewer.