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_w.f90

    r2119 r2232  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Adjustments to new topography and surface concept
    2323!
    2424! Former revisions:
     
    9797 
    9898
    99     USE wall_fluxes_mod,                                                       &
    100         ONLY :  wall_fluxes
    101 
    10299    PRIVATE
    103100    PUBLIC diffusion_w
     
    125122           
    126123       USE grid_variables,                                                     &     
    127            ONLY :  ddx, ddy, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y
     124           ONLY :  ddx, ddy
    128125           
    129126       USE indices,                                                            &           
    130            ONLY :  nxl, nxr, nyn, nys, nzb, nzb_w_inner, nzb_w_outer, nzt
     127           ONLY :  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
    131128           
    132129       USE kinds
    133130
     131       USE surface_mod,                                                        &
     132           ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
     133
    134134       IMPLICIT NONE
    135135
    136        INTEGER(iwp) ::  i     !<
    137        INTEGER(iwp) ::  j     !<
    138        INTEGER(iwp) ::  k     !<
     136       INTEGER(iwp) ::  i             !< running index x direction
     137       INTEGER(iwp) ::  j             !< running index y direction
     138       INTEGER(iwp) ::  k             !< running index z direction
     139       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
     140       INTEGER(iwp) ::  m             !< running index surface elements
     141       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
     142       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
    139143       
    140        REAL(wp) ::  kmxm  !<
    141        REAL(wp) ::  kmxp  !<
    142        REAL(wp) ::  kmym  !<
    143        REAL(wp) ::  kmyp  !<
    144 
    145        REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus  !<
    146        REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsvs  !<
    147 
    148 
    149 !
    150 !--    First calculate horizontal momentum flux w'u' and/or w'v' at vertical
    151 !--    walls, if neccessary
    152        IF ( topography /= 'flat' )  THEN
    153           CALL wall_fluxes( wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, nzb_w_inner,             &
    154                             nzb_w_outer, wall_w_x )
    155           CALL wall_fluxes( wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, nzb_w_inner,             &
    156                             nzb_w_outer, wall_w_y )
    157        ENDIF
     144       REAL(wp) ::  flag              !< flag to mask topography grid points
     145       REAL(wp) ::  kmxm              !<
     146       REAL(wp) ::  kmxp              !<
     147       REAL(wp) ::  kmym              !<
     148       REAL(wp) ::  kmyp              !<
     149       REAL(wp) ::  mask_west         !< flag to mask vertical wall west of the grid point
     150       REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point
     151       REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point
     152       REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point
     153
     154
    158155
    159156       DO  i = nxl, nxr
    160157          DO  j = nys, nyn
    161              DO  k = nzb_w_outer(j,i)+1, nzt-1
     158             DO  k = nzb+1, nzt-1
     159!
     160!--             Predetermine flag to mask topography and wall-bounded grid points.
     161                flag       = MERGE( 1.0_wp, 0.0_wp,                            &
     162                                    BTEST( wall_flags_0(k,j,i),   3 ) )
     163                mask_east  = MERGE( 1.0_wp, 0.0_wp,                            &
     164                                    BTEST( wall_flags_0(k,j,i+1), 3 ) )
     165                mask_west  = MERGE( 1.0_wp, 0.0_wp,                            &
     166                                    BTEST( wall_flags_0(k,j,i-1), 3 ) )
     167                mask_south = MERGE( 1.0_wp, 0.0_wp,                            &
     168                                    BTEST( wall_flags_0(k,j-1,i), 3 ) )
     169                mask_north = MERGE( 1.0_wp, 0.0_wp,                            &
     170                                    BTEST( wall_flags_0(k,j+1,i), 3 ) )
    162171!
    163172!--             Interpolate eddy diffusivities on staggered gridpoints
    164                 kmxp = 0.25_wp *                                               &
    165                        ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
    166                 kmxm = 0.25_wp *                                               &
    167                        ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
    168                 kmyp = 0.25_wp *                                               &
    169                        ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
    170                 kmym = 0.25_wp *                                               &
    171                        ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
     173                kmxp = 0.25_wp * ( km(k,j,i)   +   km(k,j,i+1) +               &
     174                                   km(k+1,j,i) + km(k+1,j,i+1) )
     175                kmxm = 0.25_wp * ( km(k,j,i)   + km(k,j,i-1)   +               &
     176                                   km(k+1,j,i) + km(k+1,j,i-1) )
     177                kmyp = 0.25_wp * ( km(k,j,i)   + km(k+1,j,i)   +               &
     178                                   km(k,j+1,i) + km(k+1,j+1,i) )
     179                kmym = 0.25_wp * ( km(k,j,i)   + km(k+1,j,i)   +               &
     180                                   km(k,j-1,i) + km(k+1,j-1,i) )
    172181
    173182                tend(k,j,i) = tend(k,j,i)                                      &
    174                       & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
    175                       &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
    176                       &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
    177                       &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
    178                       &   ) * ddx                                              &
    179                       & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
    180                       &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
    181                       &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
    182                       &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
    183                       &   ) * ddy                                              &
    184                       & + 2.0_wp * (                                           &
    185                       &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
    186                       &               * rho_air(k+1)                           &
    187                       & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
    188                       &               * rho_air(k)                             &
    189                       &            ) * ddzu(k+1) * drho_air_zw(k)
    190              ENDDO
    191 
    192 !
    193 !--          Wall functions at all vertical walls, where necessary
    194              IF ( wall_w_x(j,i) /= 0.0_wp  .OR.  wall_w_y(j,i) /= 0.0_wp )  THEN
    195 
    196                 DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
    197 !
    198 !--                Interpolate eddy diffusivities on staggered gridpoints
    199                    kmxp = 0.25_wp *                                            &
    200                           ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
    201                    kmxm = 0.25_wp *                                            &
    202                           ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
    203                    kmyp = 0.25_wp *                                            &
    204                           ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
    205                    kmym = 0.25_wp *                                            &
    206                           ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    207 
    208                    tend(k,j,i) = tend(k,j,i)                                   &
    209                                  + (   fwxp(j,i) * (                           &
    210                             kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
    211                           + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
    212                                                    )                           &
    213                                      - fwxm(j,i) * (                           &
    214                             kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
    215                           + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
    216                                                    )                           &
    217                                      + wall_w_x(j,i) * wsus(k,j,i)             &
    218                                    ) * ddx                                     &
    219                                  + (   fwyp(j,i) * (                           &
    220                             kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
    221                           + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
    222                                                    )                           &
    223                                      - fwym(j,i) * (                           &
    224                             kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
    225                           + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
    226                                                    )                           &
    227                                      + wall_w_y(j,i) * wsvs(k,j,i)             &
    228                                    ) * ddy                                     &
    229                                  + 2.0_wp * (                                  &
    230                            km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
    231                                        * rho_air(k+1)                          &
    232                          - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
    233                                        * rho_air(k)                            &
    234                                             ) * ddzu(k+1) * drho_air_zw(k)
    235                 ENDDO
    236              ENDIF
     183                       + ( mask_east *  kmxp * (                               &
     184                                   ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
     185                                 + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
     186                                               )                               &
     187                         - mask_west * kmxm *  (                               &
     188                                   ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
     189                                 + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
     190                                               )                               &
     191                         ) * ddx                                 * flag        &
     192                       + ( mask_north * kmyp * (                               &
     193                                   ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
     194                                 + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
     195                                               )                               &
     196                         - mask_south * kmym * (                               &
     197                                   ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
     198                                 + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
     199                                               )                               &
     200                         ) * ddy                                 * flag        &
     201                       + 2.0_wp * (                                            &
     202                         km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) )   * ddzw(k+1) &
     203                                     * rho_air(k+1)                            &
     204                       - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)   &
     205                                     * rho_air(k)                              &
     206                                  ) * ddzu(k+1) * drho_air_zw(k) * flag
     207             ENDDO
     208
     209!
     210!--          Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1)
     211!--          surfaces. Note, in the the flat case, loops won't be entered as
     212!--          start_index > end_index. Furtermore, note, no vertical natural surfaces
     213!--          so far.           
     214!--          Default-type surfaces
     215             DO  l = 0, 1
     216                surf_s = surf_def_v(l)%start_index(j,i)
     217                surf_e = surf_def_v(l)%end_index(j,i)
     218                DO  m = surf_s, surf_e
     219                   k           = surf_def_v(l)%k(m)
     220                   tend(k,j,i) = tend(k,j,i) +                                 &
     221                                     surf_def_v(l)%mom_flux_w(m) * ddy
     222                ENDDO   
     223             ENDDO
     224!
     225!--          Natural-type surfaces
     226             DO  l = 0, 1
     227                surf_s = surf_lsm_v(l)%start_index(j,i)
     228                surf_e = surf_lsm_v(l)%end_index(j,i)
     229                DO  m = surf_s, surf_e
     230                   k           = surf_lsm_v(l)%k(m)
     231                   tend(k,j,i) = tend(k,j,i) +                                 &
     232                                     surf_lsm_v(l)%mom_flux_w(m) * ddy
     233                ENDDO   
     234             ENDDO
     235!
     236!--          Urban-type surfaces
     237             DO  l = 0, 1
     238                surf_s = surf_usm_v(l)%start_index(j,i)
     239                surf_e = surf_usm_v(l)%end_index(j,i)
     240                DO  m = surf_s, surf_e
     241                   k           = surf_usm_v(l)%k(m)
     242                   tend(k,j,i) = tend(k,j,i) +                                 &
     243                                     surf_usm_v(l)%mom_flux_w(m) * ddy
     244                ENDDO   
     245             ENDDO
     246!
     247!--          Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3)
     248!--          surface.
     249!--          Default-type surfaces
     250             DO  l = 2, 3
     251                surf_s = surf_def_v(l)%start_index(j,i)
     252                surf_e = surf_def_v(l)%end_index(j,i)
     253                DO  m = surf_s, surf_e
     254                   k           = surf_def_v(l)%k(m)
     255                   tend(k,j,i) = tend(k,j,i) +                                 &
     256                                     surf_def_v(l)%mom_flux_w(m) * ddx
     257                ENDDO   
     258             ENDDO
     259!
     260!--          Natural-type surfaces
     261             DO  l = 2, 3
     262                surf_s = surf_lsm_v(l)%start_index(j,i)
     263                surf_e = surf_lsm_v(l)%end_index(j,i)
     264                DO  m = surf_s, surf_e
     265                   k           = surf_lsm_v(l)%k(m)
     266                   tend(k,j,i) = tend(k,j,i) +                                 &
     267                                     surf_lsm_v(l)%mom_flux_w(m) * ddx
     268                ENDDO   
     269             ENDDO
     270!
     271!--          Urban-type surfaces
     272             DO  l = 2, 3
     273                surf_s = surf_usm_v(l)%start_index(j,i)
     274                surf_e = surf_usm_v(l)%end_index(j,i)
     275                DO  m = surf_s, surf_e
     276                   k           = surf_usm_v(l)%k(m)
     277                   tend(k,j,i) = tend(k,j,i) +                                 &
     278                                     surf_usm_v(l)%mom_flux_w(m) * ddx
     279                ENDDO   
     280             ENDDO
    237281
    238282          ENDDO
     
    256300           
    257301       USE grid_variables,                                                     &     
    258            ONLY :  ddx, ddy, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y
     302           ONLY :  ddx, ddy
    259303           
    260304       USE indices,                                                            &           
    261            ONLY :  nxl, nxr, nyn, nys, nzb, nzb_w_inner, nzb_w_outer, nzt
     305           ONLY :  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
    262306           
    263307       USE kinds
    264308
     309       USE surface_mod,                                                        &
     310           ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
     311
    265312       IMPLICIT NONE
    266313
    267        INTEGER(iwp) ::  i     !<
    268        INTEGER(iwp) ::  j     !<
    269        INTEGER(iwp) ::  k     !<
     314
     315       INTEGER(iwp) ::  i             !< running index x direction
     316       INTEGER(iwp) ::  j             !< running index y direction
     317       INTEGER(iwp) ::  k             !< running index z direction
     318       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
     319       INTEGER(iwp) ::  m             !< running index surface elements
     320       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
     321       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
    270322       
    271        REAL(wp) ::  kmxm  !<
    272        REAL(wp) ::  kmxp  !<
    273        REAL(wp) ::  kmym  !<
    274        REAL(wp) ::  kmyp  !<
    275 
    276        REAL(wp), DIMENSION(nzb:nzt+1) ::  wsus
    277        REAL(wp), DIMENSION(nzb:nzt+1) ::  wsvs
    278 
    279 
    280        DO  k = nzb_w_outer(j,i)+1, nzt-1
     323       REAL(wp) ::  flag              !< flag to mask topography grid points
     324       REAL(wp) ::  kmxm              !<
     325       REAL(wp) ::  kmxp              !<
     326       REAL(wp) ::  kmym              !<
     327       REAL(wp) ::  kmyp              !<
     328       REAL(wp) ::  mask_west         !< flag to mask vertical wall west of the grid point
     329       REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point
     330       REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point
     331       REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point
     332
     333
     334       DO  k = nzb+1, nzt-1
     335!
     336!--       Predetermine flag to mask topography and wall-bounded grid points.
     337          flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i),   3 ) )
     338          mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 3 ) )
     339          mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 3 ) )
     340          mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 3 ) )
     341          mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 3 ) )
    281342!
    282343!--       Interpolate eddy diffusivities on staggered gridpoints
     
    287348
    288349          tend(k,j,i) = tend(k,j,i)                                            &
    289                       & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
    290                       &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
    291                       &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
    292                       &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
    293                       &   ) * ddx                                              &
    294                       & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
    295                       &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
    296                       &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
    297                       &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
    298                       &   ) * ddy                                              &
    299                       & + 2.0_wp * (                                           &
    300                       &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
    301                       &               * rho_air(k+1)                           &
    302                       & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
    303                       &               * rho_air(k)                             &
    304                       &            ) * ddzu(k+1) * drho_air_zw(k)
    305        ENDDO
    306 
    307 !
    308 !--    Wall functions at all vertical walls, where necessary
    309        IF ( wall_w_x(j,i) /= 0.0_wp  .OR.  wall_w_y(j,i) /= 0.0_wp )  THEN
    310 
    311 !
    312 !--       Calculate the horizontal momentum fluxes w'u' and/or w'v'
    313           IF ( wall_w_x(j,i) /= 0.0_wp )  THEN
    314              CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i),     &
    315                                wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
    316           ELSE
    317              wsus = 0.0_wp
    318           ENDIF
    319 
    320           IF ( wall_w_y(j,i) /= 0.0_wp )  THEN
    321              CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i),     &
    322                                wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
    323           ELSE
    324              wsvs = 0.0_wp
    325           ENDIF
    326 
    327           DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
    328 !
    329 !--          Interpolate eddy diffusivities on staggered gridpoints
    330              kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
    331              kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
    332              kmyp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
    333              kmym = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    334 
    335              tend(k,j,i) = tend(k,j,i)                                         &
    336                                  + (   fwxp(j,i) * (                           &
    337                             kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
    338                           + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
    339                                                    )                           &
    340                                      - fwxm(j,i) * (                           &
    341                             kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
    342                           + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
    343                                                    )                           &
    344                                      + wall_w_x(j,i) * wsus(k)                 &
    345                                    ) * ddx                                     &
    346                                  + (   fwyp(j,i) * (                           &
    347                             kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
    348                           + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
    349                                                    )                           &
    350                                      - fwym(j,i) * (                           &
    351                             kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
    352                           + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
    353                                                    )                           &
    354                                      + wall_w_y(j,i) * wsvs(k)                 &
    355                                    ) * ddy                                     &
    356                                  + 2.0_wp * (                                  &
    357                            km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
    358                                        * rho_air(k+1)                          &
    359                          - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
    360                                        * rho_air(k)                            &
    361                                             ) * ddzu(k+1) * drho_air_zw(k)
    362           ENDDO
    363        ENDIF
     350                       + ( mask_east *  kmxp * (                               &
     351                                   ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
     352                                 + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
     353                                               )                               &
     354                         - mask_west * kmxm *  (                               &
     355                                   ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
     356                                 + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
     357                                               )                               &
     358                         ) * ddx                                 * flag        &
     359                       + ( mask_north * kmyp * (                               &
     360                                   ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
     361                                 + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
     362                                               )                               &
     363                         - mask_south * kmym * (                               &
     364                                   ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
     365                                 + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
     366                                               )                               &
     367                         ) * ddy                                 * flag        &
     368                       + 2.0_wp * (                                            &
     369                         km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)   &
     370                                     * rho_air(k+1)                            &
     371                       - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)   &
     372                                     * rho_air(k)                              &
     373                                  ) * ddzu(k+1) * drho_air_zw(k) * flag
     374       ENDDO
     375!
     376!--    Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1)
     377!--    surfaces. Note, in the the flat case, loops won't be entered as
     378!--    start_index > end_index. Furtermore, note, no vertical natural surfaces
     379!--    so far.           
     380!--    Default-type surfaces
     381       DO  l = 0, 1
     382          surf_s = surf_def_v(l)%start_index(j,i)
     383          surf_e = surf_def_v(l)%end_index(j,i)
     384          DO  m = surf_s, surf_e
     385             k           = surf_def_v(l)%k(m)
     386             tend(k,j,i) = tend(k,j,i) +                                       &
     387                                     surf_def_v(l)%mom_flux_w(m) * ddy
     388          ENDDO   
     389       ENDDO
     390!
     391!--    Natural-type surfaces
     392       DO  l = 0, 1
     393          surf_s = surf_lsm_v(l)%start_index(j,i)
     394          surf_e = surf_lsm_v(l)%end_index(j,i)
     395          DO  m = surf_s, surf_e
     396             k           = surf_lsm_v(l)%k(m)
     397             tend(k,j,i) = tend(k,j,i) +                                       &
     398                                     surf_lsm_v(l)%mom_flux_w(m) * ddy
     399          ENDDO   
     400       ENDDO
     401!
     402!--    Urban-type surfaces
     403       DO  l = 0, 1
     404          surf_s = surf_usm_v(l)%start_index(j,i)
     405          surf_e = surf_usm_v(l)%end_index(j,i)
     406          DO  m = surf_s, surf_e
     407             k           = surf_usm_v(l)%k(m)
     408             tend(k,j,i) = tend(k,j,i) +                                       &
     409                                     surf_usm_v(l)%mom_flux_w(m) * ddy
     410          ENDDO   
     411       ENDDO
     412!
     413!--    Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3)
     414!--    surfaces.
     415!--    Default-type surfaces
     416       DO  l = 2, 3
     417          surf_s = surf_def_v(l)%start_index(j,i)
     418          surf_e = surf_def_v(l)%end_index(j,i)
     419          DO  m = surf_s, surf_e
     420             k           = surf_def_v(l)%k(m)
     421             tend(k,j,i) = tend(k,j,i) +                                       &
     422                                     surf_def_v(l)%mom_flux_w(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) +                                       &
     433                                     surf_lsm_v(l)%mom_flux_w(m) * ddx
     434          ENDDO   
     435       ENDDO
     436!
     437!--    Urban-type surfaces
     438       DO  l = 2, 3
     439          surf_s = surf_usm_v(l)%start_index(j,i)
     440          surf_e = surf_usm_v(l)%end_index(j,i)
     441          DO  m = surf_s, surf_e
     442             k           = surf_usm_v(l)%k(m)
     443             tend(k,j,i) = tend(k,j,i) +                                       &
     444                                     surf_usm_v(l)%mom_flux_w(m) * ddx
     445          ENDDO   
     446       ENDDO
     447
    364448
    365449    END SUBROUTINE diffusion_w_ij
Note: See TracChangeset for help on using the changeset viewer.