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

preliminary version, several changes to be explained later

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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 ) +           &
Note: See TracChangeset for help on using the changeset viewer.