Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Changeset 53


Ignore:
Timestamp:
Mar 7, 2007 12:33:47 PM (18 years ago)
Author:
raasch
Message:

preliminary version, several changes to be explained later

Location:
palm/trunk/SOURCE
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/CURRENT_MODIFICATIONS

    r51 r53  
    22---
    33
    4 Wall functions for vertical walls now include diabatic conditions. New subroutine wall_fluxes. New 4D-array rif_wall.
     4Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.
    55
    66new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

  • palm/trunk/SOURCE/diffusion_u.f90

    r51 r53  
    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
    77!
    88! Former revisions:

  • palm/trunk/SOURCE/diffusion_v.f90

    r51 r53  
    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
    77!
    88! Former revisions:

  • palm/trunk/SOURCE/diffusion_w.f90

    r51 r53  
    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
    77!
    88! Former revisions:

  • palm/trunk/SOURCE/production_e.f90

    r39 r53  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Wall functions now include diabatic conditions, call of routine wall_fluxes
    77!
    88! Former revisions:
    … …  
    6767
    6868       REAL    ::  def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
    69                    k1, k2, theta, temp, usvs, vsus, wsus, wsvs
     69                   k1, k2, theta, temp
     70
     71       REAL, DIMENSION(nzb:nzt+1) ::  usvs, vsus, wsus, wsvs
    7072
    7173
    … …  
    119121                   k = nzb_diff_s_inner(j,i) - 1
    120122                   dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
     123                   dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     124                                  u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
     125                   dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
     126                   dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     127                                  v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
     128                   dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
     129
    121130                   IF ( wall_e_y(j,i) /= 0.0 )  THEN
    122                       usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) &
    123                              / LOG( 0.5 * dy / z0(j,i) )
    124                       usvs = usvs * ABS( usvs )
    125                       dudy = wall_e_y(j,i) * usvs / km(k,j,i)
     131                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
     132                                          usvs, 1.0, 0.0, 0.0, 0.0 )
     133
     134                      dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
     135                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
     136                                          wsvs, 0.0, 0.0, 1.0, 0.0 )
     137                      dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
    126138                   ELSE
    127139                      dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    128140                                      u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     141                      dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     142                                      w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    129143                   ENDIF
    130                    dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    131                                   u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
    132144
    133145                   IF ( wall_e_x(j,i) /= 0.0 )  THEN
    134                       vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) &
    135                              / LOG( 0.5 * dx / z0(j,i))
    136                       vsus = vsus * ABS( vsus )
    137                       dvdx = wall_e_x(j,i) * vsus / km(k,j,i)
     146                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
     147                                          vsus, 0.0, 1.0, 0.0, 0.0 )
     148                      dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
     149                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
     150                                          wsus, 0.0, 0.0, 0.0, 1.0 )
     151                      dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
    138152                   ELSE
    139153                      dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    140154                                      v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    141                    ENDIF
    142                    dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
    143                    dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    144                                   v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
    145 
    146                    IF ( wall_e_x(j,i) /= 0.0 )  THEN
    147                       wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) &
    148                              / LOG( 0.5 * dx / z0(j,i))
    149                       wsus = wsus * ABS( wsus )
    150                       dwdx = wall_e_x(j,i) * wsus / km(k,j,i)
    151                    ELSE
    152155                      dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    153156                                      w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    154157                   ENDIF
    155                    IF ( wall_e_y(j,i) /= 0.0 ) THEN
    156                       wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) &
    157                              / LOG( 0.5 * dy / z0(j,i))
    158                       wsvs = wsvs * ABS( wsvs )
    159                       dwdy = wall_e_y(j,i) * wsvs / km(k,j,i)
    160                    ELSE
    161                       dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    162                                       w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    163                    ENDIF
    164                    dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
    165158
    166159                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
    … …  
    172165                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
    173166
    174                 ENDIF
    175 
    176              ENDDO
    177 
    178 !
    179 !--          3 - Wird nur ausgefuehrt, wenn mindestens ein Niveau
    180 !--          zwischen 2 und 4 liegt, d.h. ab einer Gebaeudemindest-
    181 !--          hoehe von 2 dz. 'Nur Wand: Wall functions'
    182              DO  j = nys, nyn
    183 
    184                 IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
    185                 THEN
     167
     168!
     169!--                3 - Wird nur ausgefuehrt, wenn mindestens ein Niveau
     170!--                zwischen 2 und 4 liegt, d.h. ab einer Gebaeudemindest-
     171!--                hoehe von 2 dz. Die Wall fluxes sind in diesem Fall
     172!--                schon unter (2) berechnet worden. 'Nur Wand: Wall functions'
    186173
    187174                   DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
    188175
    189176                      dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
     177                      dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     178                                     u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     179                      dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     180                      dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     181                                     v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     182                      dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
     183
    190184                      IF ( wall_e_y(j,i) /= 0.0 )  THEN
    191                          usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) &
    192                                 / LOG( 0.5 * dy / z0(j,i) )
    193                          usvs = usvs * ABS( usvs )
    194                          dudy = wall_e_y(j,i) * usvs / km(k,j,i)
     185                         dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
     186                         dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
    195187                      ELSE
    196188                         dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    197189                                         u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     190                         dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     191                                         w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    198192                      ENDIF
    199                       dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    200                                      u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    201193
    202194                      IF ( wall_e_x(j,i) /= 0.0 )  THEN
    203                          vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) &
    204                                 / LOG( 0.5 * dx / z0(j,i))
    205                          vsus = vsus * ABS( vsus )
    206                          dvdx = wall_e_x(j,i) * vsus / km(k,j,i)
     195                         dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
     196                         dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
    207197                      ELSE
    208198                         dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    209199                                         v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    210                       ENDIF
    211                       dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    212                       dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    213                                      v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    214 
    215                       IF ( wall_e_x(j,i) /= 0.0 )  THEN
    216                          wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) &
    217                                 / LOG( 0.5 * dx / z0(j,i))
    218                          wsus = wsus * ABS( wsus )
    219                          dwdx = wall_e_x(j,i) * wsus / km(k,j,i)
    220                       ELSE
    221200                         dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    222201                                         w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    223202                      ENDIF
    224                       IF ( wall_e_y(j,i) /= 0.0 )  THEN
    225                          wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) &
    226                                 / LOG( 0.5 * dy / z0(j,i))
    227                          wsvs = wsvs * ABS( wsvs )
    228                          dwdy = wall_e_y(j,i) * wsvs / km(k,j,i)
    229                       ELSE
    230                          dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    231                                          w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    232                       ENDIF
    233                       dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
    234203
    235204                      def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
    … …  
    423392                DO  j = nys, nyn
    424393
    425                    k = nzb_diff_s_inner(j,i)
     394                   k = nzb_diff_s_inner(j,i)-1
    426395
    427396                   IF ( .NOT. cloud_physics )  THEN
    … …  
    505474
    506475       REAL    ::  def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
    507                    k1, k2, theta, temp, usvs, vsus, wsus,wsvs
     476                   k1, k2, theta, temp
     477
     478       REAL, DIMENSION(nzb:nzt+1) ::  usvs, vsus, wsus,wsvs
    508479
    509480!
    … …  
    549520
    550521             dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
     522             dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     523                            u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
     524             dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     525             dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     526                            v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
     527             dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
     528
    551529             IF ( wall_e_y(j,i) /= 0.0 )  THEN
    552                 usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) &
    553                        / LOG( 0.5 * dy / z0(j,i) )
    554                 usvs = usvs * ABS( usvs )
    555                 dudy = wall_e_y(j,i) * usvs / km(k,j,i)
     530                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
     531                                    usvs, 1.0, 0.0, 0.0, 0.0 )
     532                dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
     533                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
     534                                    wsvs, 0.0, 0.0, 1.0, 0.0 )
     535                dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
    556536             ELSE
    557537                dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    558538                                u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     539                dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     540                                w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    559541             ENDIF
    560              dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    561                             u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
    562542
    563543             IF ( wall_e_x(j,i) /= 0.0 )  THEN
    564                 vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) &
    565                        / LOG( 0.5 * dx / z0(j,i))
    566                 vsus = vsus * ABS( vsus )
    567                 dvdx = wall_e_x(j,i) * vsus / km(k,j,i)
     544                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
     545                                    vsus, 0.0, 1.0, 0.0, 0.0 )
     546                dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
     547                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
     548                                    wsus, 0.0, 0.0, 0.0, 1.0 )
     549                dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
    568550             ELSE
    569551                dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    570552                                v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    571              ENDIF
    572              dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    573              dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    574                             v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
    575 
    576              IF ( wall_e_x(j,i) /= 0.0 )  THEN
    577                 wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) &
    578                        / LOG( 0.5 * dx / z0(j,i))
    579                 wsus = wsus * ABS( wsus )
    580                 dwdx = wall_e_x(j,i) * wsus / km(k,j,i)
    581              ELSE
    582553                dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    583554                                w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    584555             ENDIF
    585              IF ( wall_e_y(j,i) /= 0.0 )  THEN
    586                 wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) &
    587                        / LOG( 0.5 * dy / z0(j,i))
    588                 wsvs = wsvs * ABS( wsvs )
    589                 dwdy = wall_e_y(j,i) * wsvs / km(k,j,i)
    590              ELSE
    591                 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    592                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    593              ENDIF
    594              dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
    595556
    596557             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
    … …  
    605566!--          3 - Wird nur ausgefuehrt, wenn mindestens ein Niveau
    606567!--          zwischen 2 und 4 liegt, d.h. ab einer Gebaeudemindest-
    607 !--          hoehe von 2 dz. 'Nur Wand: Wall functions'
     568!--          hoehe von 2 dz. Di Wall fluxes sind in diesem Fall bereits vorab
     569!--          unter (2) berechnet worden. 'Nur Wand: Wall functions'
    608570             DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
    609571
    610572                dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
     573                dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     574                               u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     575                dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     576                dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     577                               v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     578                dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
     579
    611580                IF ( wall_e_y(j,i) /= 0.0 )  THEN
    612                    usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) &
    613                           / LOG( 0.5 * dy / z0(j,i) )
    614                    usvs = usvs * ABS( usvs )
    615                    dudy = wall_e_y(j,i) * usvs / km(k,j,i)
     581                   dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
     582                   dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
    616583                ELSE
    617584                   dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    618585                                   u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    619                 ENDIF
    620                 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    621                                u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     586                   dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     587                                   w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     588                ENDIF
    622589
    623590                IF ( wall_e_x(j,i) /= 0.0 )  THEN
    624                    vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) &
    625                           / LOG( 0.5 * dx / z0(j,i))
    626                    vsus = vsus * ABS( vsus )
    627                    dvdx = wall_e_x(j,i) * vsus / km(k,j,i)
     591                   dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
     592                   dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
    628593                ELSE
    629594                   dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    630595                                   v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    631                 ENDIF
    632                 dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    633                 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    634                                v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    635 
    636                 IF ( wall_e_x(j,i) /= 0.0 )  THEN
    637                    wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) &
    638                           / LOG( 0.5 * dx / z0(j,i))
    639                    wsus = wsus * ABS( wsus )
    640                    dwdx = wall_e_x(j,i) * wsus / km(k,j,i)
    641                 ELSE
    642596                   dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    643597                                   w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    644598                ENDIF
    645                 IF ( wall_e_y(j,i) /= 0.0 )  THEN
    646                    wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) &
    647                           / LOG( 0.5 * dy / z0(j,i))
    648                    wsvs = wsvs * ABS( wsvs )
    649                    dwdy = wall_e_y(j,i) * wsvs / km(k,j,i)
    650                 ELSE
    651                    dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    652                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    653                 ENDIF
    654                 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
    655599
    656600                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &

  • palm/trunk/SOURCE/wall_fluxes.f90

    r52 r53  
    2626    IMPLICIT NONE
    2727
    28     INTEGER ::  i, j, k, nzb_w, nzt_w
    29     REAL    ::  a, b, c1, c2, h1, h2, delta_p
     28    INTEGER ::  i, j, k, nzb_w, nzt_w, wall_index
     29    REAL    ::  a, b, c1, c2, h1, h2, zp
    3030
    3131    REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
    … …  
    3434
    3535
    36     delta_p   = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
    37     wall_flux = 0.0
     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 )
    3839
    3940!
    … …  
    4142!-- the relevant variable is defined
    4243    DO  k = nzb_w, nzt_w
     44
    4345!
    4446!--    (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
    45        rifs    = rif_wall(k,j,i,NINT(a+2*b+3*c1+4*c2))
    46        u_i     = a * u(k,j,i) +  &
    47             c1 * 0.25 * ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
    48        v_i     = b * v(k,j,i) +  &
    49             c2 * 0.25 * ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
    50        ws      = ( c1 + c2 ) * w(k,j,i) +  &
    51             a * 0.25 * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) + &
    52             b * 0.25 * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) )
    53        pt_i    = 0.5 * ( pt(k,j,i) +  &
    54             a *  pt(k,j,i-1) + b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
    55        pts     = pt_i - hom(k,1,4,0)
    56        wspts   = ws * pts
     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
    5765!
    5866!--    (2) Compute wall-parallel absolute velocity vel_total
    5967       vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 )
     68
    6069!
    6170!--    (3) Compute wall friction velocity us_wall
    6271       IF ( rifs >= 0.0 )  THEN
     72
    6373!
    6474!--       Stable stratification (and neutral)
    65           us_wall = kappa * vel_total / (                         &
    66                     LOG( delta_p / z0(j,i) ) +                    &
    67                     5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p  &
     75          us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +               &
     76                                          5.0 * rifs * ( zp - z0(j,i) ) / zp  &
    6877                                        )
    6978       ELSE
     79
    7080!
    7181!--       Unstable stratification
    7282          h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
    73           h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / delta_p * z0(j,i) ) )
     83          h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
     84
    7485!
    7586!--       If a borderline case occurs, the formula for stable stratification
    … …  
    7788!--       argument of the logarithm.
    7889          IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
    79              us_wall = kappa * vel_total / (                           &
    80                        LOG( delta_p / z0(j,i) ) +                         &
    81                        5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p &
    82                                                    )
     90             us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             &
     91                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
     92                                           )
    8393          ELSE
    84              us_wall = kappa * vel_total / (                           &
     94             us_wall = kappa * vel_total / (                                   &
    8595                            LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + &
    8696                            2.0 * ( ATAN( h2 ) - ATAN( h1 ) )                  &
    87                                                    )
     97                                           )
    8898          ENDIF
    89        ENDIF
    90 !
    91 !--    (4) Compute delta_p/L (corresponds to neutral Richardson flux number
     99
     100       ENDIF
     101
     102!
     103!--    (4) Compute zp/L (corresponds to neutral Richardson flux number
    92104!--        rifs)
    93        rifs = -1.0 * delta_p * kappa * g * wspts / &
    94                        ( pt_i * ( us_wall**3 + 1E-30 ) )
     105       rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * ( us_wall**3 + 1E-30 ) )
     106
    95107!
    96108!--    Limit the value range of the Richardson numbers.
    … …  
    101113       IF ( rifs < rif_min )  rifs = rif_min
    102114       IF ( rifs > rif_max )  rifs = rif_max
     115
    103116!
    104117!--    (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
    105118       IF ( rifs >= 0.0 )  THEN
     119
    106120!
    107121!--       Stable stratification (and neutral)
    108           wall_flux(k) = kappa *  &
    109                ( a * u(k,j,i) + b * v(k,j,i) + (c1 + c2 ) * w(k,j,i) ) / (  &
    110                LOG( delta_p / z0(j,i) ) +                                &
    111                5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p        &
    112                                                                          )
    113        ELSE
     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
    114129!
    115130!--       Unstable stratification
    116131          h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
    117           h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / delta_p * z0(j,i) ) )
     132          h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
     133
    118134!
    119135!--       If a borderline case occurs, the formula for stable stratification
    … …  
    121137!--       argument of the logarithm.
    122138          IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
    123              wall_flux(k) = kappa *  &
    124                   ( a * u(k,j,i) + b * v(k,j,i) + (c1 + c2 ) * w(k,j,i) ) / ( &
    125                   LOG( delta_p / z0(j,i) ) +                                &
    126                   5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p        &
    127                                                                             )
     139             wall_flux(k) = kappa *                                          &
     140                            ( 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             &
     143                            )
    128144          ELSE
    129              wall_flux(k) = kappa *  &
    130                   ( a * u(k,j,i) + b * v(k,j,i) + (c1 + c2 ) * w(k,j,i) ) / (  &
    131                   LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) +          &
    132                   2.0 * ( ATAN( h2 ) - ATAN( h1 ) )                            &
    133                                                                             )
     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                            )
    134150          ENDIF
     151
    135152       ENDIF
    136153       wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) )
    … …  
    138155!
    139156!--    store rifs for next time step
    140        rif_wall(k,j,i,NINT(a+2*b+3*c1+4*c2)) = rifs
     157       rif_wall(k,j,i,wall_index) = rifs
    141158
    142159    ENDDO
     160
    143161 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!
     175! Description:
     176! ------------
     177! Calculates momentum fluxes at vertical walls for routine production_e
     178! assuming Monin-Obukhov similarity.
     179! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
     180!------------------------------------------------------------------------------!
     181
     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, 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       IF ( rifs >= 0.0 )  THEN
     227
     228!
     229!--       Stable stratification (and neutral)
     230          wall_flux(k) = kappa *                                          &
     231                         ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
     232                         (  LOG( zp / z0(j,i) ) +                         &
     233                            5.0 * rifs * ( zp - z0(j,i) ) / zp            &
     234                         )
     235       ELSE
     236
     237!
     238!--       Unstable stratification
     239          h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     240          h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
     241
     242!
     243!--       If a borderline case occurs, the formula for stable stratification
     244!--       must be used anyway, or else a zero division would occur in the
     245!--       argument of the logarithm.
     246          IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
     247             wall_flux(k) = kappa *                                          &
     248                            ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
     249                            ( LOG( zp / z0(j,i) ) +                          &
     250                              5.0 * rifs * ( zp - z0(j,i) ) / zp             &
     251                            )
     252          ELSE
     253             wall_flux(k) = kappa *                                            &
     254                            ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) /   &
     255                            ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) &
     256                              + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )              &
     257                            )
     258          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
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |