Ignore:
Timestamp:
May 30, 2017 5:47:52 PM (4 years ago)
Author:
suehring
Message:

Adjustments according new topography and surface-modelling concept implemented

File:
1 edited

Legend:

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

    r2119 r2232  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Adjustments to new topography and surface concept
    2323!
    2424! Former revisions:
     
    9090 
    9191
    92     USE wall_fluxes_mod
    93 
    9492    PRIVATE
    9593    PUBLIC diffusion_v
     
    111109
    112110       USE arrays_3d,                                                          &
    113            ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w,                  &
    114                   drho_air, rho_air_zw
     111           ONLY:  ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw
    115112       
    116113       USE control_parameters,                                                 &
    117            ONLY:  constant_top_momentumflux, topography, use_surface_fluxes,   &
     114           ONLY:  constant_top_momentumflux, use_surface_fluxes,               &
    118115                  use_top_fluxes
    119116       
    120117       USE grid_variables,                                                     &
    121            ONLY:  ddx, ddy, ddy2, fxm, fxp, wall_v
     118           ONLY:  ddx, ddy, ddy2
    122119       
    123120       USE indices,                                                            &
    124            ONLY:  nxl, nxr, nyn, nys, nysv, nzb, nzb_diff_v, nzb_v_inner,      &
    125                   nzb_v_outer, nzt, nzt_diff
     121           ONLY:  nxl, nxr, nyn, nys, nysv, nzb, nzt, wall_flags_0
    126122       
    127123       USE kinds
    128124
     125       USE surface_mod,                                                        &
     126           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
     127                   surf_usm_v
     128
    129129       IMPLICIT NONE
    130130
    131        INTEGER(iwp) ::  i     !<
    132        INTEGER(iwp) ::  j     !<
    133        INTEGER(iwp) ::  k     !<
    134        REAL(wp)     ::  kmxm  !<
    135        REAL(wp)     ::  kmxp  !<
    136        REAL(wp)     ::  kmzm  !<
    137        REAL(wp)     ::  kmzp  !<
    138 
    139        REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus  !<
    140 
    141 !
    142 !--    First calculate horizontal momentum flux v'u' at vertical walls,
    143 !--    if neccessary
    144        IF ( topography /= 'flat' )  THEN
    145           CALL wall_fluxes( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, nzb_v_inner, &
    146                             nzb_v_outer, wall_v )
    147        ENDIF
     131       INTEGER(iwp) ::  i             !< running index x direction
     132       INTEGER(iwp) ::  j             !< running index y direction
     133       INTEGER(iwp) ::  k             !< running index z direction
     134       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
     135       INTEGER(iwp) ::  m             !< running index surface elements
     136       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
     137       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
     138
     139       REAL(wp)     ::  flag          !< flag to mask topography grid points
     140       REAL(wp)     ::  kmxm          !<
     141       REAL(wp)     ::  kmxp          !<
     142       REAL(wp)     ::  kmzm          !<
     143       REAL(wp)     ::  kmzp          !<
     144       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface 
     145       REAL(wp)     ::  mask_east     !< flag to mask vertical surface south of the grid point
     146       REAL(wp)     ::  mask_west     !< flag to mask vertical surface north of the grid point
     147       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface     
    148148
    149149       DO  i = nxl, nxr
     
    151151!
    152152!--          Compute horizontal diffusion
    153              DO  k = nzb_v_outer(j,i)+1, nzt
     153             DO  k = nzb+1, nzt
     154
     155!
     156!--             Predetermine flag to mask topography and wall-bounded grid points.
     157!--             It is sufficient to masked only east- and west-facing surfaces, which
     158!--             need special treatment for the v-component.
     159                flag      = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i),   2 ) )
     160                mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 2 ) )
     161                mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 2 ) )
    154162!
    155163!--             Interpolate eddy diffusivities on staggered gridpoints
    156                 kmxp = 0.25_wp * &
    157                        ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
    158                 kmxm = 0.25_wp * &
    159                        ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    160 
    161                 tend(k,j,i) = tend(k,j,i)                                      &
    162                       & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx         &
    163                       &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy         &
    164                       &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
    165                       &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
    166                       &   ) * ddx                                              &
    167                       & + 2.0_wp * (                                           &
    168                       &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )      &
    169                       &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )      &
    170                       &            ) * ddy2
     164                kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
     165                kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
     166
     167                tend(k,j,i) = tend(k,j,i) +    (                             &
     168                          mask_east * kmxp * (                               &
     169                                 ( v(k,j,i+1) - v(k,j,i)     ) * ddx         &
     170                               + ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy         &
     171                                             )                               &
     172                        - mask_west * kmxm * (                               &
     173                                 ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
     174                               + ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
     175                                             )                               &
     176                                               ) * ddx  * flag               &
     177                                    + 2.0_wp * (                             &
     178                                  km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i)   )  &
     179                                - km(k,j-1,i) * ( v(k,j,i)   - v(k,j-1,i) )  &
     180                                               ) * ddy2 * flag
     181
    171182             ENDDO
    172183
    173184!
    174 !--          Wall functions at the left and right walls, respectively
    175              IF ( wall_v(j,i) /= 0.0_wp )  THEN
    176 
    177                 DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
    178                    kmxp = 0.25_wp *                                            &
    179                           ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
    180                    kmxm = 0.25_wp *                                            &
    181                           ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    182                    
    183                    tend(k,j,i) = tend(k,j,i)                                   &
    184                                  + 2.0_wp * (                                  &
    185                                        km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
    186                                      - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
    187                                             ) * ddy2                           &
    188                                  + (   fxp(j,i) * (                            &
    189                                   kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
    190                                 + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
    191                                                   )                            &
    192                                      - fxm(j,i) * (                            &
    193                                   kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
    194                                 + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
    195                                                   )                            &
    196                                      + wall_v(j,i) * vsus(k,j,i)               &
    197                                    ) * ddx
    198                 ENDDO
    199              ENDIF
    200 
    201 !
    202 !--          Compute vertical diffusion. In case of simulating a Prandtl
    203 !--          layer, index k starts at nzb_v_inner+2.
    204              DO  k = nzb_diff_v(j,i), nzt_diff
     185!--          Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3)
     186!--          surfaces. Note, in the the flat case, loops won't be entered as
     187!--          start_index > end_index. Furtermore, note, no vertical natural surfaces
     188!--          so far.           
     189!--          Default-type surfaces
     190             DO  l = 2, 3
     191                surf_s = surf_def_v(l)%start_index(j,i)
     192                surf_e = surf_def_v(l)%end_index(j,i)
     193                DO  m = surf_s, surf_e
     194                   k           = surf_def_v(l)%k(m)
     195                   tend(k,j,i) = tend(k,j,i) +                                 &
     196                                    surf_def_v(l)%mom_flux_uv(m) * ddx
     197                ENDDO   
     198             ENDDO
     199!
     200!--          Natural-type surfaces
     201             DO  l = 2, 3
     202                surf_s = surf_lsm_v(l)%start_index(j,i)
     203                surf_e = surf_lsm_v(l)%end_index(j,i)
     204                DO  m = surf_s, surf_e
     205                   k           = surf_lsm_v(l)%k(m)
     206                   tend(k,j,i) = tend(k,j,i) +                                 &
     207                                    surf_lsm_v(l)%mom_flux_uv(m) * ddx
     208                ENDDO   
     209             ENDDO
     210!
     211!--          Urban-type surfaces
     212             DO  l = 2, 3
     213                surf_s = surf_usm_v(l)%start_index(j,i)
     214                surf_e = surf_usm_v(l)%end_index(j,i)
     215                DO  m = surf_s, surf_e
     216                   k           = surf_usm_v(l)%k(m)
     217                   tend(k,j,i) = tend(k,j,i) +                                 &
     218                                    surf_usm_v(l)%mom_flux_uv(m) * ddx
     219                ENDDO   
     220             ENDDO
     221!
     222!--          Compute vertical diffusion. In case of simulating a surface layer,
     223!--          respective grid diffusive fluxes are masked (flag 10) within this
     224!--          loop, and added further below, else, simple gradient approach is
     225!--          applied. Model top is also mask if top-momentum flux is given.
     226             DO  k = nzb+1, nzt
     227!
     228!--             Determine flags to mask topography below and above. Flag 2 is
     229!--             used to mask topography in general, while flag 8 implies also
     230!--             information about use_surface_fluxes. Flag 9 is used to control
     231!--             momentum flux at model top. 
     232                mask_bottom = MERGE( 1.0_wp, 0.0_wp,                           &
     233                                     BTEST( wall_flags_0(k-1,j,i), 8 ) )
     234                mask_top    = MERGE( 1.0_wp, 0.0_wp,                           &
     235                                     BTEST( wall_flags_0(k+1,j,i), 8 ) ) *     &
     236                              MERGE( 1.0_wp, 0.0_wp,                           &
     237                                     BTEST( wall_flags_0(k+1,j,i), 9 ) )
     238                flag        = MERGE( 1.0_wp, 0.0_wp,                           &
     239                                     BTEST( wall_flags_0(k,j,i), 2 ) )
    205240!
    206241!--             Interpolate eddy diffusivities on staggered gridpoints
     
    213248                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
    214249                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
    215                       &            ) * rho_air_zw(k)                           &
     250                      &            ) * rho_air_zw(k)   * mask_top              &
    216251                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
    217252                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
    218                       &            ) * rho_air_zw(k-1)                         &
    219                       &   ) * ddzw(k) * drho_air(k)
     253                      &            ) * rho_air_zw(k-1) * mask_bottom           &
     254                      &   ) * ddzw(k) * drho_air(k) * flag
    220255             ENDDO
    221256
     
    228263!--          comparison with other (LES) models showed that the values of
    229264!--          the momentum flux becomes too large in this case.
    230 !--          The term containing w(k-1,..) (see above equation) is removed here
    231 !--          because the vertical velocity is assumed to be zero at the surface.
    232265             IF ( use_surface_fluxes )  THEN
    233                 k = nzb_v_inner(j,i)+1
    234 !
    235 !--             Interpolate eddy diffusivities on staggered gridpoints
    236                 kmzp = 0.25_wp *                                               &
    237                        ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    238 
    239                 tend(k,j,i) = tend(k,j,i)                                      &
    240                       & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
    241                       &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
    242                       &            ) * rho_air_zw(k)                           &
    243                       &   - ( -vsws(j,i) )                                     &
    244                       &   ) * ddzw(k) * drho_air(k)
     266!
     267!--             Default-type surfaces, upward-facing
     268                surf_s = surf_def_h(0)%start_index(j,i)
     269                surf_e = surf_def_h(0)%end_index(j,i)
     270                DO  m = surf_s, surf_e
     271                   k   = surf_def_h(0)%k(m)
     272
     273                   tend(k,j,i) = tend(k,j,i)                                   &
     274                        + ( - ( - surf_def_h(0)%vsws(m) )                      &
     275                          ) * ddzw(k) * drho_air(k)
     276                ENDDO
     277!
     278!--             Default-type surfaces, dowward-facing
     279                surf_s = surf_def_h(1)%start_index(j,i)
     280                surf_e = surf_def_h(1)%end_index(j,i)
     281                DO  m = surf_s, surf_e
     282                   k   = surf_def_h(1)%k(m)
     283
     284                   tend(k,j,i) = tend(k,j,i)                                   &
     285                        + ( - surf_def_h(1)%vsws(m)                            &
     286                          ) * ddzw(k) * drho_air(k)
     287                ENDDO
     288!
     289!--             Natural-type surfaces, upward-facing
     290                surf_s = surf_lsm_h%start_index(j,i)
     291                surf_e = surf_lsm_h%end_index(j,i)
     292                DO  m = surf_s, surf_e
     293                   k   = surf_lsm_h%k(m)
     294
     295                   tend(k,j,i) = tend(k,j,i)                                   &
     296                        + ( - ( - surf_lsm_h%vsws(m) )                         &
     297                          ) * ddzw(k) * drho_air(k)
     298
     299                ENDDO
     300!
     301!--             Urban-type surfaces, upward-facing
     302                surf_s = surf_usm_h%start_index(j,i)
     303                surf_e = surf_usm_h%end_index(j,i)
     304                DO  m = surf_s, surf_e
     305                   k   = surf_usm_h%k(m)
     306
     307                   tend(k,j,i) = tend(k,j,i)                                   &
     308                        + ( - ( - surf_usm_h%vsws(m) )                         &
     309                          ) * ddzw(k) * drho_air(k)
     310
     311                ENDDO
    245312             ENDIF
    246 
    247 !
    248 !--          Vertical diffusion at the first gridpoint below the top boundary,
    249 !--          if the momentum flux at the top is prescribed by the user
    250              IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
    251                 k = nzt
    252 !
    253 !--             Interpolate eddy diffusivities on staggered gridpoints
    254                 kmzm = 0.25_wp *                                               &
    255                        ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
    256 
    257                 tend(k,j,i) = tend(k,j,i)                                      &
    258                       & + ( ( -vswst(j,i) )                                    &
    259                       &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
    260                       &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
    261                       &            ) * rho_air_zw(k-1)                         &
    262                       &   ) * ddzw(k) * drho_air(k)
     313!
     314!--          Add momentum flux at model top
     315             IF ( use_top_fluxes )  THEN
     316                surf_s = surf_def_h(2)%start_index(j,i)
     317                surf_e = surf_def_h(2)%end_index(j,i)
     318                DO  m = surf_s, surf_e
     319
     320                   k   = surf_def_h(2)%k(m)
     321
     322                   tend(k,j,i) = tend(k,j,i)                                   &
     323                           + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
     324                ENDDO
    263325             ENDIF
    264326
     
    277339
    278340       USE arrays_3d,                                                          &
    279            ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w,                  &
    280                   drho_air, rho_air_zw
     341           ONLY:  ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw
    281342       
    282343       USE control_parameters,                                                 &
    283            ONLY:  constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
     344           ONLY:  constant_top_momentumflux, use_surface_fluxes,               &
     345                  use_top_fluxes
    284346       
    285347       USE grid_variables,                                                     &
    286            ONLY:  ddx, ddy, ddy2, fxm, fxp, wall_v
     348           ONLY:  ddx, ddy, ddy2
    287349       
    288350       USE indices,                                                            &
    289            ONLY:  nzb, nzb_diff_v, nzb_v_inner, nzb_v_outer, nzt, nzt_diff
     351           ONLY:  nzb, nzt, wall_flags_0
    290352       
    291353       USE kinds
    292354
     355       USE surface_mod,                                                        &
     356           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
     357                   surf_usm_v
     358
    293359       IMPLICIT NONE
    294360
    295        INTEGER(iwp) ::  i     !<
    296        INTEGER(iwp) ::  j     !<
    297        INTEGER(iwp) ::  k     !<
    298        REAL(wp)     ::  kmxm  !<
    299        REAL(wp)     ::  kmxp  !<
    300        REAL(wp)     ::  kmzm  !<
    301        REAL(wp)     ::  kmzp  !<
    302 
    303        REAL(wp), DIMENSION(nzb:nzt+1) ::  vsus  !<
     361
     362       INTEGER(iwp) ::  i             !< running index x direction
     363       INTEGER(iwp) ::  j             !< running index y direction
     364       INTEGER(iwp) ::  k             !< running index z direction
     365       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
     366       INTEGER(iwp) ::  m             !< running index surface elements
     367       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
     368       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
     369
     370       REAL(wp)     ::  flag          !< flag to mask topography grid points
     371       REAL(wp)     ::  kmxm          !<
     372       REAL(wp)     ::  kmxp          !<
     373       REAL(wp)     ::  kmzm          !<
     374       REAL(wp)     ::  kmzp          !<
     375       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface 
     376       REAL(wp)     ::  mask_east     !< flag to mask vertical surface south of the grid point
     377       REAL(wp)     ::  mask_west     !< flag to mask vertical surface north of the grid point
     378       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
    304379
    305380!
    306381!--    Compute horizontal diffusion
    307        DO  k = nzb_v_outer(j,i)+1, nzt
     382       DO  k = nzb+1, nzt
     383!
     384!--       Predetermine flag to mask topography and wall-bounded grid points.
     385!--       It is sufficient to masked only east- and west-facing surfaces, which
     386!--       need special treatment for the v-component.
     387          flag      = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i),   2 ) )
     388          mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 2 ) )
     389          mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 2 ) )
    308390!
    309391!--       Interpolate eddy diffusivities on staggered gridpoints
     
    311393          kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    312394
    313           tend(k,j,i) = tend(k,j,i)                                            &
    314                       & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx         &
    315                       &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy         &
    316                       &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
    317                       &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
    318                       &   ) * ddx                                              &
    319                       & + 2.0_wp * (                                           &
    320                       &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )      &
    321                       &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )      &
    322                       &            ) * ddy2
     395          tend(k,j,i) = tend(k,j,i) +          (                             &
     396                          mask_east * kmxp * (                               &
     397                                 ( v(k,j,i+1) - v(k,j,i)     ) * ddx         &
     398                               + ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy         &
     399                                             )                               &
     400                        - mask_west * kmxm * (                               &
     401                                 ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
     402                               + ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
     403                                             )                               &
     404                                               ) * ddx  * flag               &
     405                                    + 2.0_wp * (                             &
     406                                  km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i)   )  &
     407                                - km(k,j-1,i) * ( v(k,j,i)   - v(k,j-1,i) )  &
     408                                               ) * ddy2 * flag
    323409       ENDDO
    324410
    325411!
    326 !--    Wall functions at the left and right walls, respectively
    327        IF ( wall_v(j,i) /= 0.0_wp )  THEN
    328 
    329 !
    330 !--       Calculate the horizontal momentum flux v'u'
    331           CALL wall_fluxes( i, j, nzb_v_inner(j,i)+1, nzb_v_outer(j,i),        &
    332                             vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp )
    333 
    334           DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
    335              kmxp = 0.25_wp *                                                  &
    336                     ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
    337              kmxm = 0.25_wp *                                                  &
    338                     ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    339 
    340              tend(k,j,i) = tend(k,j,i)                                         &
    341                                  + 2.0_wp * (                                  &
    342                                        km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
    343                                      - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
    344                                             ) * ddy2                           &
    345                                  + (   fxp(j,i) * (                            &
    346                                   kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
    347                                 + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
    348                                                   )                            &
    349                                      - fxm(j,i) * (                            &
    350                                   kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
    351                                 + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
    352                                                   )                            &
    353                                      + wall_v(j,i) * vsus(k)                   &
    354                                    ) * ddx
    355           ENDDO
    356        ENDIF
    357 
    358 !
    359 !--    Compute vertical diffusion. In case of simulating a Prandtl layer,
    360 !--    index k starts at nzb_v_inner+2.
    361        DO  k = nzb_diff_v(j,i), nzt_diff
     412!--    Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3)
     413!--    surfaces. Note, in the the flat case, loops won't be entered as
     414!--    start_index > end_index. Furtermore, note, no vertical natural surfaces
     415!--    so far.           
     416!--    Default-type surfaces
     417       DO  l = 2, 3
     418          surf_s = surf_def_v(l)%start_index(j,i)
     419          surf_e = surf_def_v(l)%end_index(j,i)
     420          DO  m = surf_s, surf_e
     421             k           = surf_def_v(l)%k(m)
     422             tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddx
     423          ENDDO   
     424       ENDDO
     425!
     426!--    Natural-type surfaces
     427       DO  l = 2, 3
     428          surf_s = surf_lsm_v(l)%start_index(j,i)
     429          surf_e = surf_lsm_v(l)%end_index(j,i)
     430          DO  m = surf_s, surf_e
     431             k           = surf_lsm_v(l)%k(m)
     432             tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddx
     433          ENDDO   
     434       ENDDO
     435!
     436!--    Urban-type surfaces
     437       DO  l = 2, 3
     438          surf_s = surf_usm_v(l)%start_index(j,i)
     439          surf_e = surf_usm_v(l)%end_index(j,i)
     440          DO  m = surf_s, surf_e
     441             k           = surf_usm_v(l)%k(m)
     442             tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddx
     443          ENDDO   
     444       ENDDO
     445!
     446!--    Compute vertical diffusion. In case of simulating a surface layer,
     447!--    respective grid diffusive fluxes are masked (flag 8) within this
     448!--    loop, and added further below, else, simple gradient approach is
     449!--    applied. Model top is also mask if top-momentum flux is given.
     450       DO  k = nzb+1, nzt
     451!
     452!--       Determine flags to mask topography below and above. Flag 2 is
     453!--       used to mask topography in general, while flag 10 implies also
     454!--       information about use_surface_fluxes. Flag 9 is used to control
     455!--       momentum flux at model top. 
     456          mask_bottom = MERGE( 1.0_wp, 0.0_wp,                                 &
     457                               BTEST( wall_flags_0(k-1,j,i), 8 ) )
     458          mask_top    = MERGE( 1.0_wp, 0.0_wp,                                 &
     459                               BTEST( wall_flags_0(k+1,j,i), 8 ) ) *           &
     460                        MERGE( 1.0_wp, 0.0_wp,                                 &
     461                               BTEST( wall_flags_0(k+1,j,i), 9 ) )
     462          flag        = MERGE( 1.0_wp, 0.0_wp,                                 &
     463                               BTEST( wall_flags_0(k,j,i), 2 ) )
    362464!
    363465!--       Interpolate eddy diffusivities on staggered gridpoints
     
    368470                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
    369471                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
    370                       &            ) * rho_air_zw(k)                           &
     472                      &            ) * rho_air_zw(k)   * mask_top              &
    371473                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
    372474                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
    373                       &            ) * rho_air_zw(k-1)                         &
    374                       &   ) * ddzw(k) * drho_air(k)
     475                      &            ) * rho_air_zw(k-1) * mask_bottom           &
     476                      &   ) * ddzw(k) * drho_air(k) * flag
    375477       ENDDO
    376478
     
    383485!--    other (LES) models showed that the values of the momentum flux becomes
    384486!--    too large in this case.
    385 !--    The term containing w(k-1,..) (see above equation) is removed here
    386 !--    because the vertical velocity is assumed to be zero at the surface.
    387487       IF ( use_surface_fluxes )  THEN
    388           k = nzb_v_inner(j,i)+1
    389 !
    390 !--       Interpolate eddy diffusivities on staggered gridpoints
    391           kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    392 
    393           tend(k,j,i) = tend(k,j,i)                                            &
    394                       & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
    395                       &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
    396                       &            ) * rho_air_zw(k)                           &
    397                       &   - ( -vsws(j,i) )                                     &
    398                       &   ) * ddzw(k) * drho_air(k)
     488!
     489!--       Default-type surfaces, upward-facing
     490          surf_s = surf_def_h(0)%start_index(j,i)
     491          surf_e = surf_def_h(0)%end_index(j,i)
     492          DO  m = surf_s, surf_e
     493             k   = surf_def_h(0)%k(m)
     494
     495             tend(k,j,i) = tend(k,j,i)                                         &
     496                        + ( - ( - surf_def_h(0)%vsws(m) )                      &
     497                          ) * ddzw(k) * drho_air(k)
     498          ENDDO
     499!
     500!--       Default-type surfaces, dowward-facing
     501          surf_s = surf_def_h(1)%start_index(j,i)
     502          surf_e = surf_def_h(1)%end_index(j,i)
     503          DO  m = surf_s, surf_e
     504             k   = surf_def_h(1)%k(m)
     505
     506             tend(k,j,i) = tend(k,j,i)                                         &
     507                        + ( - surf_def_h(1)%vsws(m)                            &
     508                          ) * ddzw(k) * drho_air(k)
     509          ENDDO
     510!
     511!--       Natural-type surfaces, upward-facing
     512          surf_s = surf_lsm_h%start_index(j,i)
     513          surf_e = surf_lsm_h%end_index(j,i)
     514          DO  m = surf_s, surf_e
     515             k   = surf_lsm_h%k(m)
     516
     517             tend(k,j,i) = tend(k,j,i)                                         &
     518                        + ( - ( - surf_lsm_h%vsws(m) )                         &
     519                          ) * ddzw(k) * drho_air(k)
     520
     521          ENDDO
     522!
     523!--       Urban-type surfaces, upward-facing
     524          surf_s = surf_usm_h%start_index(j,i)
     525          surf_e = surf_usm_h%end_index(j,i)
     526          DO  m = surf_s, surf_e
     527             k   = surf_usm_h%k(m)
     528
     529             tend(k,j,i) = tend(k,j,i)                                         &
     530                        + ( - ( - surf_usm_h%vsws(m) )                         &
     531                          ) * ddzw(k) * drho_air(k)
     532
     533          ENDDO
    399534       ENDIF
    400 
    401 !
    402 !--    Vertical diffusion at the first gridpoint below the top boundary,
    403 !--    if the momentum flux at the top is prescribed by the user
    404        IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
    405           k = nzt
    406 !
    407 !--       Interpolate eddy diffusivities on staggered gridpoints
    408           kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
    409 
    410           tend(k,j,i) = tend(k,j,i)                                            &
    411                       & + ( ( -vswst(j,i) )                                    &
    412                       &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
    413                       &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
    414                       &            ) * rho_air_zw(k-1)                         &
    415                       &   ) * ddzw(k) * drho_air(k)
     535!
     536!--    Add momentum flux at model top
     537       IF ( use_top_fluxes )  THEN
     538          surf_s = surf_def_h(2)%start_index(j,i)
     539          surf_e = surf_def_h(2)%end_index(j,i)
     540          DO  m = surf_s, surf_e
     541
     542             k   = surf_def_h(2)%k(m)
     543
     544             tend(k,j,i) = tend(k,j,i)                                         &
     545                           + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
     546          ENDDO
    416547       ENDIF
    417548
Note: See TracChangeset for help on using the changeset viewer.