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

    r2119 r2232  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Adjustments to new topography and surface concept
    2323!
    2424! Former revisions:
     
    110110!> Call for all grid points
    111111!------------------------------------------------------------------------------!
    112     SUBROUTINE diffusion_s( s, s_flux_b, s_flux_t, wall_s_flux )
     112    SUBROUTINE diffusion_s( s, s_flux_def_h_up,    s_flux_def_h_down,          &
     113                               s_flux_t,                                       &
     114                               s_flux_lsm_h_up,    s_flux_usm_h_up,            &
     115                               s_flux_def_v_north, s_flux_def_v_south,         &
     116                               s_flux_def_v_east,  s_flux_def_v_west,          &
     117                               s_flux_lsm_v_north, s_flux_lsm_v_south,         &
     118                               s_flux_lsm_v_east,  s_flux_lsm_v_west,          &
     119                               s_flux_usm_v_north, s_flux_usm_v_south,         &
     120                               s_flux_usm_v_east,  s_flux_usm_v_west )
    113121
    114122       USE arrays_3d,                                                          &
     
    119127       
    120128       USE grid_variables,                                                     &
    121            ONLY:  ddx2, ddy2, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y
     129           ONLY:  ddx2, ddy2
    122130       
    123131       USE indices,                                                            &
    124132           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb,             &
    125                   nzb_diff_s_inner, nzb_s_inner, nzb_s_outer, nzt, nzt_diff
     133                  nzt, wall_flags_0
    126134       
    127135       USE kinds
    128136
     137       USE surface_mod,                                                        &
     138           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
     139                   surf_usm_v
     140
    129141       IMPLICIT NONE
    130142
    131        INTEGER(iwp) ::  i                 !<
    132        INTEGER(iwp) ::  j                 !<
    133        INTEGER(iwp) ::  k                 !<
    134        REAL(wp)     ::  wall_s_flux(0:4)  !<
    135        REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_b, s_flux_t !<
     143       INTEGER(iwp) ::  i             !< running index x direction
     144       INTEGER(iwp) ::  j             !< running index y direction
     145       INTEGER(iwp) ::  k             !< running index z direction
     146       INTEGER(iwp) ::  m             !< running index surface elements
     147       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
     148       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
     149
     150       REAL(wp) ::  flag              !< flag to mask topography grid points
     151       REAL(wp) ::  mask_bottom       !< flag to mask vertical upward-facing surface     
     152       REAL(wp) ::  mask_east         !< flag to mask vertical surface east of the grid point
     153       REAL(wp) ::  mask_north        !< flag to mask vertical surface north of the grid point
     154       REAL(wp) ::  mask_south        !< flag to mask vertical surface south of the grid point
     155       REAL(wp) ::  mask_west         !< flag to mask vertical surface west of the grid point
     156       REAL(wp) ::  mask_top          !< flag to mask vertical downward-facing surface 
     157
     158       REAL(wp), DIMENSION(1:surf_def_v(0)%ns) ::  s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
     159       REAL(wp), DIMENSION(1:surf_def_v(1)%ns) ::  s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
     160       REAL(wp), DIMENSION(1:surf_def_v(2)%ns) ::  s_flux_def_v_east  !< flux at east-facing vertical default-type surfaces
     161       REAL(wp), DIMENSION(1:surf_def_v(3)%ns) ::  s_flux_def_v_west  !< flux at west-facing vertical default-type surfaces
     162       REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  s_flux_def_h_up    !< flux at horizontal upward-facing default-type surfaces
     163       REAL(wp), DIMENSION(1:surf_def_h(1)%ns) ::  s_flux_def_h_down  !< flux at horizontal donwward-facing default-type surfaces
     164       REAL(wp), DIMENSION(1:surf_lsm_h%ns)    ::  s_flux_lsm_h_up    !< flux at horizontal upward-facing natural-type surfaces
     165       REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) ::  s_flux_lsm_v_north !< flux at north-facing vertical natural-type surfaces
     166       REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) ::  s_flux_lsm_v_south !< flux at south-facing vertical natural-type surfaces
     167       REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) ::  s_flux_lsm_v_east  !< flux at east-facing vertical natural-type surfaces
     168       REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) ::  s_flux_lsm_v_west  !< flux at west-facing vertical natural-type surfaces
     169       REAL(wp), DIMENSION(1:surf_usm_h%ns)    ::  s_flux_usm_h_up    !< flux at horizontal upward-facing urban-type surfaces
     170       REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) ::  s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
     171       REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) ::  s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
     172       REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) ::  s_flux_usm_v_east  !< flux at east-facing vertical urban-type surfaces
     173       REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) ::  s_flux_usm_v_west  !< flux at west-facing vertical urban-type surfaces
     174       REAL(wp), DIMENSION(1:surf_def_h(2)%ns) ::  s_flux_t           !< flux at model top
    136175#if defined( __nopointer )
    137176       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s  !<
     
    144183!
    145184!--          Compute horizontal diffusion
    146              DO  k = nzb_s_outer(j,i)+1, nzt
     185             DO  k = nzb+1, nzt
     186!
     187!--             Predetermine flag to mask topography and wall-bounded grid points
     188                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     189!
     190!--             Predetermine flag to mask wall-bounded grid points, equivalent to
     191!--             former s_outer array
     192                mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 0 ) )
     193                mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 0 ) )
     194                mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 0 ) )
     195                mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 0 ) )
    147196
    148197                tend(k,j,i) = tend(k,j,i)                                      &
    149198                                          + 0.5_wp * (                         &
    150                         ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
    151                       - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
    152                                                      ) * ddx2                  &
     199                        mask_east  * ( kh(k,j,i) + kh(k,j,i+1) )               &
     200                                   * ( s(k,j,i+1) - s(k,j,i)   )               &
     201                      - mask_west  * ( kh(k,j,i) + kh(k,j,i-1) )               &
     202                                   * ( s(k,j,i)   - s(k,j,i-1) )               &
     203                                                     ) * ddx2 * flag           &
    153204                                          + 0.5_wp * (                         &
    154                         ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
    155                       - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
    156                                                      ) * ddy2
    157              ENDDO
    158 
    159 !
    160 !--          Apply prescribed horizontal wall heatflux where necessary
    161              IF ( ( wall_w_x(j,i) /= 0.0_wp ) .OR. ( wall_w_y(j,i) /= 0.0_wp ) &
    162                 )  THEN
    163                 DO  k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
    164 
     205                        mask_north * ( kh(k,j,i) + kh(k,j+1,i) )               &
     206                                   * ( s(k,j+1,i) - s(k,j,i)   )               &
     207                      - mask_south * ( kh(k,j,i) + kh(k,j-1,i) )               &
     208                                   * ( s(k,j,i)   - s(k,j-1,i) )               &
     209                                                     ) * ddy2 * flag
     210             ENDDO
     211
     212!
     213!--          Apply prescribed horizontal wall heatflux where necessary. First,
     214!--          determine start and end index for respective (j,i)-index. Please
     215!--          note, in the flat case following loop will not be entered, as
     216!--          surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces
     217!--          so far.
     218!--          First, for default-type surfaces
     219!--          North-facing vertical default-type surfaces
     220             surf_s = surf_def_v(0)%start_index(j,i)
     221             surf_e = surf_def_v(0)%end_index(j,i)
     222             DO  m = surf_s, surf_e
     223                k           = surf_def_v(0)%k(m)
     224                tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy2
     225             ENDDO
     226!
     227!--          South-facing vertical default-type surfaces
     228             surf_s = surf_def_v(1)%start_index(j,i)
     229             surf_e = surf_def_v(1)%end_index(j,i)
     230             DO  m = surf_s, surf_e
     231                k           = surf_def_v(1)%k(m)
     232                tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy2
     233             ENDDO
     234!
     235!--          East-facing vertical default-type surfaces
     236             surf_s = surf_def_v(2)%start_index(j,i)
     237             surf_e = surf_def_v(2)%end_index(j,i)
     238             DO  m = surf_s, surf_e
     239                k           = surf_def_v(2)%k(m)
     240                tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx2
     241             ENDDO
     242!
     243!--          West-facing vertical default-type surfaces
     244             surf_s = surf_def_v(3)%start_index(j,i)
     245             surf_e = surf_def_v(3)%end_index(j,i)
     246             DO  m = surf_s, surf_e
     247                k           = surf_def_v(3)%k(m)
     248                tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx2
     249             ENDDO
     250!
     251!--          Now, for natural-type surfaces.
     252!--          North-facing
     253             surf_s = surf_lsm_v(0)%start_index(j,i)
     254             surf_e = surf_lsm_v(0)%end_index(j,i)
     255             DO  m = surf_s, surf_e
     256                k           = surf_lsm_v(0)%k(m)
     257                tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy2
     258             ENDDO
     259!
     260!--          South-facing
     261             surf_s = surf_lsm_v(1)%start_index(j,i)
     262             surf_e = surf_lsm_v(1)%end_index(j,i)
     263             DO  m = surf_s, surf_e
     264                k           = surf_lsm_v(1)%k(m)
     265                tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy2
     266             ENDDO
     267!
     268!--          East-facing
     269             surf_s = surf_lsm_v(2)%start_index(j,i)
     270             surf_e = surf_lsm_v(2)%end_index(j,i)
     271             DO  m = surf_s, surf_e
     272                k           = surf_lsm_v(2)%k(m)
     273                tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx2
     274             ENDDO
     275!
     276!--          West-facing
     277             surf_s = surf_lsm_v(3)%start_index(j,i)
     278             surf_e = surf_lsm_v(3)%end_index(j,i)
     279             DO  m = surf_s, surf_e
     280                k           = surf_lsm_v(3)%k(m)
     281                tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx2
     282             ENDDO
     283!
     284!--          Now, for urban-type surfaces.
     285!--          North-facing
     286             surf_s = surf_usm_v(0)%start_index(j,i)
     287             surf_e = surf_usm_v(0)%end_index(j,i)
     288             DO  m = surf_s, surf_e
     289                k           = surf_usm_v(0)%k(m)
     290                tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy2
     291             ENDDO
     292!
     293!--          South-facing
     294             surf_s = surf_usm_v(1)%start_index(j,i)
     295             surf_e = surf_usm_v(1)%end_index(j,i)
     296             DO  m = surf_s, surf_e
     297                k           = surf_usm_v(1)%k(m)
     298                tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy2
     299             ENDDO
     300!
     301!--          East-facing
     302             surf_s = surf_usm_v(2)%start_index(j,i)
     303             surf_e = surf_usm_v(2)%end_index(j,i)
     304             DO  m = surf_s, surf_e
     305                k           = surf_usm_v(2)%k(m)
     306                tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx2
     307             ENDDO
     308!
     309!--          West-facing
     310             surf_s = surf_usm_v(3)%start_index(j,i)
     311             surf_e = surf_usm_v(3)%end_index(j,i)
     312             DO  m = surf_s, surf_e
     313                k           = surf_usm_v(3)%k(m)
     314                tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx2
     315             ENDDO
     316
     317!
     318!--          Compute vertical diffusion. In case that surface fluxes have been
     319!--          prescribed or computed at bottom and/or top, index k starts/ends at
     320!--          nzb+2 or nzt-1, respectively. Model top is also mask if top flux
     321!--          is given.
     322             DO  k = nzb+1, nzt
     323!
     324!--             Determine flags to mask topography below and above. Flag 0 is
     325!--             used to mask topography in general, and flag 8 implies
     326!--             information about use_surface_fluxes. Flag 9 is used to control
     327!--             flux at model top.
     328                mask_bottom = MERGE( 1.0_wp, 0.0_wp,                           &
     329                                     BTEST( wall_flags_0(k-1,j,i), 8 ) )
     330                mask_top    = MERGE( 1.0_wp, 0.0_wp,                           &
     331                                     BTEST( wall_flags_0(k+1,j,i), 8 ) ) *     &
     332                              MERGE( 1.0_wp, 0.0_wp,                           &
     333                                     BTEST( wall_flags_0(k+1,j,i), 9 ) )
     334                flag        = MERGE( 1.0_wp, 0.0_wp,                           &
     335                                     BTEST( wall_flags_0(k,j,i), 0 ) )
     336
     337                tend(k,j,i) = tend(k,j,i)                                      &
     338                                       + 0.5_wp * (                            &
     339                                      ( kh(k,j,i) + kh(k+1,j,i) ) *            &
     340                                          ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
     341                                                            * rho_air_zw(k)    &
     342                                                            * mask_top         &
     343                                    - ( kh(k,j,i) + kh(k-1,j,i) ) *            &
     344                                          ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
     345                                                            * rho_air_zw(k-1)  &
     346                                                            * mask_bottom      &
     347                                                  ) * ddzw(k) * drho_air(k)    &
     348                                                              * flag
     349             ENDDO
     350
     351!
     352!--          Vertical diffusion at horizontal walls.
     353             IF ( use_surface_fluxes )  THEN
     354!
     355!--             Default-type surfaces, upward-facing               
     356                surf_s = surf_def_h(0)%start_index(j,i)
     357                surf_e = surf_def_h(0)%end_index(j,i)
     358                DO  m = surf_s, surf_e
     359
     360                   k   = surf_def_h(0)%k(m)
     361                   tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m)              &
     362                                       * ddzw(k) * drho_air(k)
     363
     364                ENDDO
     365!
     366!--             Default-type surfaces, downward-facing               
     367                surf_s = surf_def_h(1)%start_index(j,i)
     368                surf_e = surf_def_h(1)%end_index(j,i)
     369                DO  m = surf_s, surf_e
     370
     371                   k   = surf_def_h(1)%k(m)
     372                   tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m)            &
     373                                       * ddzw(k) * drho_air(k)
     374
     375                ENDDO
     376!
     377!--             Natural-type surfaces, upward-facing 
     378                surf_s = surf_lsm_h%start_index(j,i)
     379                surf_e = surf_lsm_h%end_index(j,i)
     380                DO  m = surf_s, surf_e
     381
     382                   k   = surf_lsm_h%k(m)
     383                   tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m)              &
     384                                       * ddzw(k) * drho_air(k)
     385
     386                ENDDO
     387!
     388!--             Urban-type surfaces, upward-facing     
     389                surf_s = surf_usm_h%start_index(j,i)
     390                surf_e = surf_usm_h%end_index(j,i)
     391                DO  m = surf_s, surf_e
     392
     393                   k   = surf_usm_h%k(m)
     394                   tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m)              &
     395                                       * ddzw(k) * drho_air(k)
     396
     397                ENDDO
     398
     399             ENDIF
     400!
     401!--          Vertical diffusion at the last computational gridpoint along z-direction
     402             IF ( use_top_fluxes )  THEN
     403                surf_s = surf_def_h(2)%start_index(j,i)
     404                surf_e = surf_def_h(2)%end_index(j,i)
     405                DO  m = surf_s, surf_e
     406
     407                   k   = surf_def_h(2)%k(m)
    165408                   tend(k,j,i) = tend(k,j,i)                                   &
    166                                                 + ( fwxp(j,i) * 0.5_wp *       &
    167                         ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
    168                         + ( 1.0_wp - fwxp(j,i) ) * wall_s_flux(1)              &
    169                                                    -fwxm(j,i) * 0.5_wp *       &
    170                         ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
    171                         + ( 1.0_wp - fwxm(j,i) ) * wall_s_flux(2)              &
    172                                                   ) * ddx2                     &
    173                                                 + ( fwyp(j,i) * 0.5_wp *       &
    174                         ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
    175                         + ( 1.0_wp - fwyp(j,i) ) * wall_s_flux(3)              &
    176                                                    -fwym(j,i) * 0.5_wp *       &
    177                         ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
    178                         + ( 1.0_wp - fwym(j,i) ) * wall_s_flux(4)              &
    179                                                   ) * ddy2
     409                           + ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
    180410                ENDDO
    181411             ENDIF
    182412
    183 !
    184 !--          Compute vertical diffusion. In case that surface fluxes have been
    185 !--          prescribed or computed at bottom and/or top, index k starts/ends at
    186 !--          nzb+2 or nzt-1, respectively.
    187              DO  k = nzb_diff_s_inner(j,i), nzt_diff
    188 
    189                 tend(k,j,i) = tend(k,j,i)                                      &
    190                                        + 0.5_wp * (                            &
    191             ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
    192                                                             * rho_air_zw(k)    &
    193           - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
    194                                                             * rho_air_zw(k-1)  &
    195                                                   ) * ddzw(k) * drho_air(k)
    196              ENDDO
    197 
    198 !
    199 !--          Vertical diffusion at the first computational gridpoint along
    200 !--          z-direction
    201              IF ( use_surface_fluxes )  THEN
    202 
    203                 k = nzb_s_inner(j,i)+1
    204 
    205                 tend(k,j,i) = tend(k,j,i)                                      &
    206                                        + ( 0.5_wp * ( kh(k,j,i)+kh(k+1,j,i) )  &
    207                                                   * ( s(k+1,j,i)-s(k,j,i) )    &
    208                                                   * ddzu(k+1)                  &
    209                                                   * rho_air_zw(k)              &
    210                                            + s_flux_b(j,i)                     &
    211                                          ) * ddzw(k) * drho_air(k)
    212 
    213              ENDIF
    214 
    215 !
    216 !--          Vertical diffusion at the last computational gridpoint along
    217 !--          z-direction
    218              IF ( use_top_fluxes )  THEN
    219 
    220                 k = nzt
    221 
    222                 tend(k,j,i) = tend(k,j,i)                                      &
    223                                        + ( - s_flux_t(j,i)                     &
    224                                            - 0.5_wp * ( kh(k-1,j,i)+kh(k,j,i) )&
    225                                                     * ( s(k,j,i)-s(k-1,j,i) )  &
    226                                                     * ddzu(k)                  &
    227                                                     * rho_air_zw(k-1)          &
    228                                          ) * ddzw(k) * drho_air(k)
    229 
    230              ENDIF
    231 
    232413          ENDDO
    233414       ENDDO
    234415
    235416    END SUBROUTINE diffusion_s
    236 
    237417
    238418!------------------------------------------------------------------------------!
     
    241421!> Call for grid point i,j
    242422!------------------------------------------------------------------------------!
    243     SUBROUTINE diffusion_s_ij( i, j, s, s_flux_b, s_flux_t, wall_s_flux )
     423    SUBROUTINE diffusion_s_ij( i, j, s,                                        &
     424                               s_flux_def_h_up,    s_flux_def_h_down,          &
     425                               s_flux_t,                                       &
     426                               s_flux_lsm_h_up,    s_flux_usm_h_up,            &
     427                               s_flux_def_v_north, s_flux_def_v_south,         &
     428                               s_flux_def_v_east,  s_flux_def_v_west,          &
     429                               s_flux_lsm_v_north, s_flux_lsm_v_south,         &
     430                               s_flux_lsm_v_east,  s_flux_lsm_v_west,          &
     431                               s_flux_usm_v_north, s_flux_usm_v_south,         &
     432                               s_flux_usm_v_east,  s_flux_usm_v_west )       
    244433
    245434       USE arrays_3d,                                                          &
     
    250439       
    251440       USE grid_variables,                                                     &
    252            ONLY:  ddx2, ddy2, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y
     441           ONLY:  ddx2, ddy2
    253442       
    254443       USE indices,                                                            &
    255            ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_diff_s_inner, nzb_s_inner,  &
    256                   nzb_s_outer, nzt, nzt_diff
     444           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
    257445       
    258446       USE kinds
    259447
     448       USE surface_mod,                                                        &
     449           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
     450                   surf_usm_v
     451
    260452       IMPLICIT NONE
    261453
    262        INTEGER(iwp) ::  i                 !<
    263        INTEGER(iwp) ::  j                 !<
    264        INTEGER(iwp) ::  k                 !<
    265        REAL(wp)     ::  wall_s_flux(0:4)  !<
    266        REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_b  !<
    267        REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_t  !<
     454       INTEGER(iwp) ::  i             !< running index x direction
     455       INTEGER(iwp) ::  j             !< running index y direction
     456       INTEGER(iwp) ::  k             !< running index z direction
     457       INTEGER(iwp) ::  m             !< running index surface elements
     458       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
     459       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
     460
     461       REAL(wp) ::  flag              !< flag to mask topography grid points
     462       REAL(wp) ::  mask_bottom       !< flag to mask vertical upward-facing surface     
     463       REAL(wp) ::  mask_east         !< flag to mask vertical surface east of the grid point
     464       REAL(wp) ::  mask_north        !< flag to mask vertical surface north of the grid point
     465       REAL(wp) ::  mask_south        !< flag to mask vertical surface south of the grid point
     466       REAL(wp) ::  mask_west         !< flag to mask vertical surface west of the grid point
     467       REAL(wp) ::  mask_top          !< flag to mask vertical downward-facing surface 
     468
     469       REAL(wp), DIMENSION(1:surf_def_v(0)%ns) ::  s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
     470       REAL(wp), DIMENSION(1:surf_def_v(1)%ns) ::  s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
     471       REAL(wp), DIMENSION(1:surf_def_v(2)%ns) ::  s_flux_def_v_east  !< flux at east-facing vertical default-type surfaces
     472       REAL(wp), DIMENSION(1:surf_def_v(3)%ns) ::  s_flux_def_v_west  !< flux at west-facing vertical default-type surfaces
     473       REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  s_flux_def_h_up    !< flux at horizontal upward-facing default-type surfaces
     474       REAL(wp), DIMENSION(1:surf_def_h(1)%ns) ::  s_flux_def_h_down  !< flux at horizontal donwward-facing default-type surfaces
     475       REAL(wp), DIMENSION(1:surf_lsm_h%ns)    ::  s_flux_lsm_h_up    !< flux at horizontal upward-facing natural-type surfaces
     476       REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) ::  s_flux_lsm_v_north !< flux at north-facing vertical urban-type surfaces
     477       REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) ::  s_flux_lsm_v_south !< flux at south-facing vertical urban-type surfaces
     478       REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) ::  s_flux_lsm_v_east  !< flux at east-facing vertical urban-type surfaces
     479       REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) ::  s_flux_lsm_v_west  !< flux at west-facing vertical urban-type surfaces
     480       REAL(wp), DIMENSION(1:surf_usm_h%ns)    ::  s_flux_usm_h_up    !< flux at horizontal upward-facing urban-type surfaces
     481       REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) ::  s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
     482       REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) ::  s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
     483       REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) ::  s_flux_usm_v_east  !< flux at east-facing vertical urban-type surfaces
     484       REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) ::  s_flux_usm_v_west  !< flux at west-facing vertical urban-type surfaces
     485       REAL(wp), DIMENSION(1:surf_def_h(2)%ns) ::  s_flux_t           !< flux at model top
    268486#if defined( __nopointer )
    269487       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s !<
     
    274492!
    275493!--    Compute horizontal diffusion
    276        DO  k = nzb_s_outer(j,i)+1, nzt
     494       DO  k = nzb+1, nzt
     495!
     496!--       Predetermine flag to mask topography and wall-bounded grid points
     497          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     498!
     499!--       Predetermine flag to mask wall-bounded grid points, equivalent to
     500!--       former s_outer array
     501          mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 0 ) )
     502          mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 0 ) )
     503          mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 0 ) )
     504          mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 0 ) )
     505!
     506!--       Finally, determine flag to mask both topography itself as well
     507!--       as wall-bounded grid points, which will be treated further below
    277508
    278509          tend(k,j,i) = tend(k,j,i)                                            &
    279510                                          + 0.5_wp * (                         &
    280                         ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
    281                       - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
    282                                                      ) * ddx2                  &
     511                            mask_east  * ( kh(k,j,i) + kh(k,j,i+1) )           &
     512                                       * ( s(k,j,i+1) - s(k,j,i)   )           &
     513                          - mask_west  * ( kh(k,j,i) + kh(k,j,i-1) )           &
     514                                       * ( s(k,j,i)   - s(k,j,i-1) )           &
     515                                                     ) * ddx2 * flag           &
    283516                                          + 0.5_wp * (                         &
    284                         ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
    285                       - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
    286                                                      ) * ddy2
    287        ENDDO
    288 
    289 !
    290 !--    Apply prescribed horizontal wall heatflux where necessary
    291        IF ( ( wall_w_x(j,i) /= 0.0_wp ) .OR. ( wall_w_y(j,i) /= 0.0_wp ) )     &
    292        THEN
    293           DO  k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
    294 
    295              tend(k,j,i) = tend(k,j,i)                                         &
    296                                                 + ( fwxp(j,i) * 0.5_wp *       &
    297                         ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
    298                         + ( 1.0_wp - fwxp(j,i) ) * wall_s_flux(1)              &
    299                                                    -fwxm(j,i) * 0.5_wp *       &
    300                         ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
    301                         + ( 1.0_wp - fwxm(j,i) ) * wall_s_flux(2)              &
    302                                                   ) * ddx2                     &
    303                                                 + ( fwyp(j,i) * 0.5_wp *       &
    304                         ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
    305                         + ( 1.0_wp - fwyp(j,i) ) * wall_s_flux(3)              &
    306                                                    -fwym(j,i) * 0.5_wp *       &
    307                         ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
    308                         + ( 1.0_wp - fwym(j,i) ) * wall_s_flux(4)              &
    309                                                   ) * ddy2
     517                            mask_north * ( kh(k,j,i) + kh(k,j+1,i) )           &
     518                                       * ( s(k,j+1,i) - s(k,j,i)   )           &
     519                          - mask_south * ( kh(k,j,i) + kh(k,j-1,i) )           &
     520                                       * ( s(k,j,i)  - s(k,j-1,i)  )           &
     521                                                     ) * ddy2 * flag
     522       ENDDO
     523
     524!
     525!--    Apply prescribed horizontal wall heatflux where necessary. First,
     526!--    determine start and end index for respective (j,i)-index. Please
     527!--    note, in the flat case following loops will not be entered, as
     528!--    surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces
     529!--    so far.
     530!--    First, for default-type surfaces
     531!--    North-facing vertical default-type surfaces
     532       surf_s = surf_def_v(0)%start_index(j,i)
     533       surf_e = surf_def_v(0)%end_index(j,i)
     534       DO  m = surf_s, surf_e
     535          k           = surf_def_v(0)%k(m)
     536          tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy2
     537       ENDDO
     538!
     539!--    South-facing vertical default-type surfaces
     540       surf_s = surf_def_v(1)%start_index(j,i)
     541       surf_e = surf_def_v(1)%end_index(j,i)
     542       DO  m = surf_s, surf_e
     543          k           = surf_def_v(1)%k(m)
     544          tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy2
     545       ENDDO
     546!
     547!--    East-facing vertical default-type surfaces
     548       surf_s = surf_def_v(2)%start_index(j,i)
     549       surf_e = surf_def_v(2)%end_index(j,i)
     550       DO  m = surf_s, surf_e
     551          k           = surf_def_v(2)%k(m)
     552          tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx2
     553       ENDDO
     554!
     555!--    West-facing vertical default-type surfaces
     556       surf_s = surf_def_v(3)%start_index(j,i)
     557       surf_e = surf_def_v(3)%end_index(j,i)
     558       DO  m = surf_s, surf_e
     559          k           = surf_def_v(3)%k(m)
     560          tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx2
     561       ENDDO
     562!
     563!--    Now, for natural-type surfaces
     564!--    North-facing
     565       surf_s = surf_lsm_v(0)%start_index(j,i)
     566       surf_e = surf_lsm_v(0)%end_index(j,i)
     567       DO  m = surf_s, surf_e
     568          k           = surf_lsm_v(0)%k(m)
     569          tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy2
     570       ENDDO
     571!
     572!--    South-facing
     573       surf_s = surf_lsm_v(1)%start_index(j,i)
     574       surf_e = surf_lsm_v(1)%end_index(j,i)
     575       DO  m = surf_s, surf_e
     576          k           = surf_lsm_v(1)%k(m)
     577          tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy2
     578       ENDDO
     579!
     580!--    East-facing
     581       surf_s = surf_lsm_v(2)%start_index(j,i)
     582       surf_e = surf_lsm_v(2)%end_index(j,i)
     583       DO  m = surf_s, surf_e
     584          k           = surf_lsm_v(2)%k(m)
     585          tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx2
     586       ENDDO
     587!
     588!--    West-facing
     589       surf_s = surf_lsm_v(3)%start_index(j,i)
     590       surf_e = surf_lsm_v(3)%end_index(j,i)
     591       DO  m = surf_s, surf_e
     592          k           = surf_lsm_v(3)%k(m)
     593          tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx2
     594       ENDDO
     595!
     596!--    Now, for urban-type surfaces
     597!--    North-facing
     598       surf_s = surf_usm_v(0)%start_index(j,i)
     599       surf_e = surf_usm_v(0)%end_index(j,i)
     600       DO  m = surf_s, surf_e
     601          k           = surf_usm_v(0)%k(m)
     602          tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy2
     603       ENDDO
     604!
     605!--    South-facing
     606       surf_s = surf_usm_v(1)%start_index(j,i)
     607       surf_e = surf_usm_v(1)%end_index(j,i)
     608       DO  m = surf_s, surf_e
     609          k           = surf_usm_v(1)%k(m)
     610          tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy2
     611       ENDDO
     612!
     613!--    East-facing
     614       surf_s = surf_usm_v(2)%start_index(j,i)
     615       surf_e = surf_usm_v(2)%end_index(j,i)
     616       DO  m = surf_s, surf_e
     617          k           = surf_usm_v(2)%k(m)
     618          tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx2
     619       ENDDO
     620!
     621!--    West-facing
     622       surf_s = surf_usm_v(3)%start_index(j,i)
     623       surf_e = surf_usm_v(3)%end_index(j,i)
     624       DO  m = surf_s, surf_e
     625          k           = surf_usm_v(3)%k(m)
     626          tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx2
     627       ENDDO
     628
     629
     630!
     631!--    Compute vertical diffusion. In case that surface fluxes have been
     632!--    prescribed or computed at bottom and/or top, index k starts/ends at
     633!--    nzb+2 or nzt-1, respectively. Model top is also mask if top flux
     634!--    is given.
     635       DO  k = nzb+1, nzt
     636!
     637!--       Determine flags to mask topography below and above. Flag 0 is
     638!--       used to mask topography in general, and flag 8 implies
     639!--       information about use_surface_fluxes. Flag 9 is used to control
     640!--       flux at model top.   
     641          mask_bottom = MERGE( 1.0_wp, 0.0_wp,                                 &
     642                               BTEST( wall_flags_0(k-1,j,i), 8 ) )
     643          mask_top    = MERGE( 1.0_wp, 0.0_wp,                                 &
     644                               BTEST( wall_flags_0(k+1,j,i), 8 ) )  *          &
     645                        MERGE( 1.0_wp, 0.0_wp,                                 &
     646                               BTEST( wall_flags_0(k+1,j,i), 9 ) )
     647          flag        = MERGE( 1.0_wp, 0.0_wp,                                 &
     648                               BTEST( wall_flags_0(k,j,i), 0 ) )
     649
     650          tend(k,j,i) = tend(k,j,i)                                            &
     651                                       + 0.5_wp * (                            &
     652                                      ( kh(k,j,i) + kh(k+1,j,i) ) *            &
     653                                          ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
     654                                                            * rho_air_zw(k)    &
     655                                                            * mask_top         &
     656                                    - ( kh(k,j,i) + kh(k-1,j,i) ) *            &
     657                                          ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
     658                                                            * rho_air_zw(k-1)  &
     659                                                            * mask_bottom      &
     660                                                  ) * ddzw(k) * drho_air(k)    &
     661                                                              * flag
     662       ENDDO
     663
     664!
     665!--    Vertical diffusion at horizontal walls.
     666!--    TO DO: Adjust for downward facing walls and mask already in main loop
     667       IF ( use_surface_fluxes )  THEN
     668!
     669!--       Default-type surfaces, upward-facing
     670          surf_s = surf_def_h(0)%start_index(j,i)
     671          surf_e = surf_def_h(0)%end_index(j,i)
     672          DO  m = surf_s, surf_e
     673
     674             k   = surf_def_h(0)%k(m)
     675
     676             tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m)                    &
     677                                       * ddzw(k) * drho_air(k)
     678          ENDDO
     679!
     680!--       Default-type surfaces, downward-facing
     681          surf_s = surf_def_h(1)%start_index(j,i)
     682          surf_e = surf_def_h(1)%end_index(j,i)
     683          DO  m = surf_s, surf_e
     684
     685             k   = surf_def_h(1)%k(m)
     686
     687             tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m)                  &
     688                                       * ddzw(k) * drho_air(k)
     689          ENDDO
     690!
     691!--       Natural-type surfaces, upward-facing
     692          surf_s = surf_lsm_h%start_index(j,i)
     693          surf_e = surf_lsm_h%end_index(j,i)
     694          DO  m = surf_s, surf_e
     695             k   = surf_lsm_h%k(m)
     696
     697             tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m)                    &
     698                                       * ddzw(k) * drho_air(k)
     699          ENDDO
     700!
     701!--       Urban-type surfaces, upward-facing
     702          surf_s = surf_usm_h%start_index(j,i)
     703          surf_e = surf_usm_h%end_index(j,i)
     704          DO  m = surf_s, surf_e
     705             k   = surf_usm_h%k(m)
     706
     707             tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m)                    &
     708                                       * ddzw(k) * drho_air(k)
    310709          ENDDO
    311710       ENDIF
    312 
    313 !
    314 !--    Compute vertical diffusion. In case that surface fluxes have been
    315 !--    prescribed or computed at bottom and/or top, index k starts/ends at
    316 !--    nzb+2 or nzt-1, respectively.
    317        DO  k = nzb_diff_s_inner(j,i), nzt_diff
    318 
    319           tend(k,j,i) = tend(k,j,i)                                            &
    320                                        + 0.5_wp * (                            &
    321             ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
    322                                                             * rho_air_zw(k)    &
    323           - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
    324                                                             * rho_air_zw(k-1)  &
    325                                                   ) * ddzw(k) * drho_air(k)
    326        ENDDO
    327 
    328 !
    329 !--    Vertical diffusion at the first computational gridpoint along z-direction
    330        IF ( use_surface_fluxes )  THEN
    331 
    332           k = nzb_s_inner(j,i)+1
    333 
    334           tend(k,j,i) = tend(k,j,i) + ( 0.5_wp * ( kh(k,j,i)+kh(k+1,j,i) )     &
    335                                                * ( s(k+1,j,i)-s(k,j,i) )       &
    336                                                * ddzu(k+1)                     &
    337                                                * rho_air_zw(k)                 &
    338                                         + s_flux_b(j,i)                        &
    339                                       ) * ddzw(k) * drho_air(k)
    340 
    341        ENDIF
    342 
    343711!
    344712!--    Vertical diffusion at the last computational gridpoint along z-direction
    345713       IF ( use_top_fluxes )  THEN
    346 
    347           k = nzt
    348 
    349           tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(j,i)                        &
    350                                       - 0.5_wp * ( kh(k-1,j,i)+kh(k,j,i) )     &
    351                                                * ( s(k,j,i)-s(k-1,j,i) )       &
    352                                                * ddzu(k)                       &
    353                                                * rho_air_zw(k-1)               &
    354                                       ) * ddzw(k) * drho_air(k)
    355 
     714          surf_s = surf_def_h(2)%start_index(j,i)
     715          surf_e = surf_def_h(2)%end_index(j,i)
     716          DO  m = surf_s, surf_e
     717
     718             k   = surf_def_h(2)%k(m)
     719             tend(k,j,i) = tend(k,j,i)                                         &
     720                           + ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
     721          ENDDO
    356722       ENDIF
    357723
Note: See TracChangeset for help on using the changeset viewer.