Ignore:
Timestamp:
Jun 29, 2020 12:36:47 PM (15 months ago)
Author:
raasch
Message:

files re-formatted to follow the PALM coding standard

File:
1 edited

Legend:

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

    r4360 r4583  
    11!> @file diffusion_v.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
    9 !
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
    13 !
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    1918!
    2019! Current revisions:
     
    2524! -----------------
    2625! $Id$
    27 ! Introduction of wall_flags_total_0, which currently sets bits based on static
    28 ! topography information used in wall_flags_static_0
    29 !
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4360 2020-01-07 11:25:50Z suehring
     29! Introduction of wall_flags_total_0, which currently sets bits based on static topography
     30! information used in wall_flags_static_0
     31!
    3032! 4329 2019-12-10 15:46:36Z motisi
    3133! Renamed wall_flags_0 to wall_flags_static_0
    32 ! 
     34!
    3335! 4182 2019-08-22 15:20:23Z scharf
    3436! Corrected "Former revisions" section
    35 ! 
     37!
    3638! 3655 2019-01-07 16:51:22Z knoop
    3739! OpenACC port for SPEC
     
    4446! ------------
    4547!> Diffusion term of the v-component
    46 !------------------------------------------------------------------------------!
     48!--------------------------------------------------------------------------------------------------!
    4749 MODULE diffusion_v_mod
    48  
     50
    4951
    5052    PRIVATE
     
    5961
    6062
    61 !------------------------------------------------------------------------------!
     63!--------------------------------------------------------------------------------------------------!
    6264! Description:
    6365! ------------
    6466!> Call for all grid points
    65 !------------------------------------------------------------------------------!
     67!--------------------------------------------------------------------------------------------------!
    6668    SUBROUTINE diffusion_v
    6769
    68        USE arrays_3d,                                                          &
    69            ONLY:  ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw
    70        
    71        USE control_parameters,                                                 &
    72            ONLY:  constant_top_momentumflux, use_surface_fluxes,               &
    73                   use_top_fluxes
    74        
    75        USE grid_variables,                                                     &
     70       USE arrays_3d,                                                                              &
     71           ONLY:  ddzu, ddzw, drho_air, km, rho_air_zw, tend, u, v, w
     72
     73       USE control_parameters,                                                                     &
     74           ONLY:  constant_top_momentumflux, use_surface_fluxes,  use_top_fluxes
     75
     76       USE grid_variables,                                                                         &
    7677           ONLY:  ddx, ddy, ddy2
    77        
    78        USE indices,                                                            &
     78
     79       USE indices,                                                                                &
    7980           ONLY:  nxl, nxr, nyn, nysv, nzb, nzt, wall_flags_total_0
    80        
     81
    8182       USE kinds
    8283
    83        USE surface_mod,                                                        &
    84            ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
    85                    surf_usm_v
     84       USE surface_mod,                                                                            &
     85           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
    8686
    8787       IMPLICIT NONE
     
    100100       REAL(wp)     ::  kmzm          !< diffusion coefficient on bottom of the gridbox - interpolated onto yv-zw grid
    101101       REAL(wp)     ::  kmzp          !< diffusion coefficient on top of the gridbox - interpolated onto yv-zw grid
    102        REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface 
    103        REAL(wp)     ::  mask_east     !< flag to mask vertical surface south of the grid point 
    104        REAL(wp)     ::  mask_west     !< flag to mask vertical surface north of the grid point 
    105        REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface     
     102       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface
     103       REAL(wp)     ::  mask_east     !< flag to mask vertical surface south of the grid point
     104       REAL(wp)     ::  mask_west     !< flag to mask vertical surface north of the grid point
     105       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
    106106
    107107       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
     
    122122
    123123!
    124 !--             Predetermine flag to mask topography and wall-bounded grid points. 
    125 !--             It is sufficient to masked only east- and west-facing surfaces, which
    126 !--             need special treatment for the v-component.
    127                 flag      = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   2 ) ) 
     124!--             Predetermine flag to mask topography and wall-bounded grid points.
     125!--             It is sufficient to masked only east- and west-facing surfaces, which need special
     126!--             treatment for the v-component.
     127                flag      = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   2 ) )
    128128                mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 2 ) )
    129129                mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 2 ) )
     
    133133                kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    134134
    135                 tend(k,j,i) = tend(k,j,i) +    (                             &
    136                           mask_east * kmxp * (                               &
    137                                  ( v(k,j,i+1) - v(k,j,i)     ) * ddx         &
    138                                + ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy         &
    139                                              )                               &
    140                         - mask_west * kmxm * (                               &
    141                                  ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
    142                                + ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
    143                                              )                               &
    144                                                ) * ddx  * flag               &
    145                                     + 2.0_wp * (                             &
    146                                   km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i)   )  &
    147                                 - km(k,j-1,i) * ( v(k,j,i)   - v(k,j-1,i) )  &
    148                                                ) * ddy2 * flag
     135                tend(k,j,i) = tend(k,j,i) +          (                                             &
     136                                mask_east * kmxp * (                                               &
     137                                       ( v(k,j,i+1) - v(k,j,i)     ) * ddx                         &
     138                                     + ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy                         &
     139                                                   )                                               &
     140                              - mask_west * kmxm * (                                               &
     141                                       ( v(k,j,i) - v(k,j,i-1) ) * ddx                             &
     142                                     + ( u(k,j,i) - u(k,j-1,i) ) * ddy                             &
     143                                                   )                                               &
     144                                                     ) * ddx  * flag                               &
     145                                          + 2.0_wp * (                                             &
     146                                        km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i)   )                  &
     147                                      - km(k,j-1,i) * ( v(k,j,i)   - v(k,j-1,i) )                  &
     148                                                     ) * ddy2 * flag
    149149
    150150             ENDDO
    151151
    152152!
    153 !--          Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3)
    154 !--          surfaces. Note, in the the flat case, loops won't be entered as
    155 !--          start_index > end_index. Furtermore, note, no vertical natural surfaces
    156 !--          so far.           
     153!--          Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3) surfaces. Note,
     154!--          in the the flat case, loops won't be entered as start_index > end_index. Furtermore,
     155!--          note, no vertical natural surfaces so far.
    157156!--          Default-type surfaces
    158157             DO  l = 2, 3
     
    161160                DO  m = surf_s, surf_e
    162161                   k           = surf_def_v(l)%k(m)
    163                    tend(k,j,i) = tend(k,j,i) +                                 &
    164                                     surf_def_v(l)%mom_flux_uv(m) * ddx
    165                 ENDDO   
     162                   tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddx
     163                ENDDO
    166164             ENDDO
    167165!
     
    172170                DO  m = surf_s, surf_e
    173171                   k           = surf_lsm_v(l)%k(m)
    174                    tend(k,j,i) = tend(k,j,i) +                                 &
    175                                     surf_lsm_v(l)%mom_flux_uv(m) * ddx
    176                 ENDDO   
     172                   tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddx
     173                ENDDO
    177174             ENDDO
    178175!
     
    183180                DO  m = surf_s, surf_e
    184181                   k           = surf_usm_v(l)%k(m)
    185                    tend(k,j,i) = tend(k,j,i) +                                 &
    186                                     surf_usm_v(l)%mom_flux_uv(m) * ddx
    187                 ENDDO   
     182                   tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddx
     183                ENDDO
    188184             ENDDO
    189185!
    190 !--          Compute vertical diffusion. In case of simulating a surface layer,
    191 !--          respective grid diffusive fluxes are masked (flag 10) within this
    192 !--          loop, and added further below, else, simple gradient approach is
    193 !--          applied. Model top is also mask if top-momentum flux is given.
     186!--          Compute vertical diffusion. In case of simulating a surface layer, respective grid
     187!--          diffusive fluxes are masked (flag 10) within this loop, and added further below, else,
     188!--          simple gradient approach is applied. Model top is also mask if top-momentum flux is
     189!--          given.
    194190             DO  k = nzb+1, nzt
    195191!
    196 !--             Determine flags to mask topography below and above. Flag 2 is
    197 !--             used to mask topography in general, while flag 8 implies also
    198 !--             information about use_surface_fluxes. Flag 9 is used to control
    199 !--             momentum flux at model top. 
    200                 mask_bottom = MERGE( 1.0_wp, 0.0_wp,                           &
    201                                  BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
    202                 mask_top    = MERGE( 1.0_wp, 0.0_wp,                           &
    203                                  BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *   &
    204                               MERGE( 1.0_wp, 0.0_wp,                           &
    205                                  BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
    206                 flag        = MERGE( 1.0_wp, 0.0_wp,                           &
    207                                  BTEST( wall_flags_total_0(k,j,i), 2 ) )
     192!--             Determine flags to mask topography below and above. Flag 2 is used to mask
     193!--             topography in general, while flag 8 implies also information about
     194!--             use_surface_fluxes. Flag 9 is used to control momentum flux at model top.
     195                mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
     196                mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *   &
     197                              MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
     198                flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
    208199!
    209200!--             Interpolate eddy diffusivities on staggered gridpoints
    210                 kmzp = 0.25_wp * &
    211                        ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    212                 kmzm = 0.25_wp * &
    213                        ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
    214 
    215                 tend(k,j,i) = tend(k,j,i)                                      &
    216                       & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
    217                       &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
    218                       &            ) * rho_air_zw(k)   * mask_top              &
    219                       &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
    220                       &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
    221                       &            ) * rho_air_zw(k-1) * mask_bottom           &
    222                       &   ) * ddzw(k) * drho_air(k) * flag
     201                kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
     202                kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
     203
     204                tend(k,j,i) = tend(k,j,i)                                                          &
     205                              & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)                 &
     206                              &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy                       &
     207                              &            ) * rho_air_zw(k)   * mask_top                          &
     208                              &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)               &
     209                              &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy                   &
     210                              &            ) * rho_air_zw(k-1) * mask_bottom                       &
     211                              &   ) * ddzw(k) * drho_air(k) * flag
    223212             ENDDO
    224213
    225214!
    226 !--          Vertical diffusion at the first grid point above the surface,
    227 !--          if the momentum flux at the bottom is given by the Prandtl law
    228 !--          or if it is prescribed by the user.
    229 !--          Difference quotient of the momentum flux is not formed over
    230 !--          half of the grid spacing (2.0*ddzw(k)) any more, since the
    231 !--          comparison with other (LES) models showed that the values of
    232 !--          the momentum flux becomes too large in this case.
     215!--          Vertical diffusion at the first grid point above the surface, if the momentum flux at
     216!--          the bottom is given by the Prandtl law or if it is prescribed by the user.
     217!--          Difference quotient of the momentum flux is not formed over half of the grid spacing
     218!--          (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the
     219!--          values of the momentum flux becomes too large in this case.
    233220             IF ( use_surface_fluxes )  THEN
    234221!
     
    239226                   k   = surf_def_h(0)%k(m)
    240227
    241                    tend(k,j,i) = tend(k,j,i)                                   &
    242                         + ( - ( - surf_def_h(0)%vsws(m) )                      &
    243                           ) * ddzw(k) * drho_air(k)
     228                   tend(k,j,i) = tend(k,j,i)                                                       &
     229                                 + ( - ( - surf_def_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
    244230                ENDDO
    245231!
     
    250236                   k   = surf_def_h(1)%k(m)
    251237
    252                    tend(k,j,i) = tend(k,j,i)                                   &
    253                         + ( - surf_def_h(1)%vsws(m)                            &
    254                           ) * ddzw(k) * drho_air(k)
     238                   tend(k,j,i) = tend(k,j,i)                                                       &
     239                                 + ( - surf_def_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
    255240                ENDDO
    256241!
     
    261246                   k   = surf_lsm_h%k(m)
    262247
    263                    tend(k,j,i) = tend(k,j,i)                                   &
    264                         + ( - ( - surf_lsm_h%vsws(m) )                         &
    265                           ) * ddzw(k) * drho_air(k)
     248                   tend(k,j,i) = tend(k,j,i)                                                       &
     249                                 + ( - ( - surf_lsm_h%vsws(m) ) ) * ddzw(k) * drho_air(k)
    266250
    267251                ENDDO
     
    273257                   k   = surf_usm_h%k(m)
    274258
    275                    tend(k,j,i) = tend(k,j,i)                                   &
    276                         + ( - ( - surf_usm_h%vsws(m) )                         &
    277                           ) * ddzw(k) * drho_air(k)
     259                   tend(k,j,i) = tend(k,j,i)                                                       &
     260                                 + ( - ( - surf_usm_h%vsws(m) ) ) * ddzw(k) * drho_air(k)
    278261
    279262                ENDDO
     
    288271                   k   = surf_def_h(2)%k(m)
    289272
    290                    tend(k,j,i) = tend(k,j,i)                                   &
    291                            + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
     273                   tend(k,j,i) = tend(k,j,i)                                                       &
     274                                 + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
    292275                ENDDO
    293276             ENDIF
     
    299282
    300283
    301 !------------------------------------------------------------------------------!
     284!--------------------------------------------------------------------------------------------------!
    302285! Description:
    303286! ------------
    304287!> Call for grid point i,j
    305 !------------------------------------------------------------------------------!
     288!--------------------------------------------------------------------------------------------------!
    306289    SUBROUTINE diffusion_v_ij( i, j )
    307290
    308        USE arrays_3d,                                                          &
    309            ONLY:  ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw
    310        
    311        USE control_parameters,                                                 &
    312            ONLY:  constant_top_momentumflux, use_surface_fluxes,               &
    313                   use_top_fluxes
    314        
    315        USE grid_variables,                                                     &
     291       USE arrays_3d,                                                                              &
     292           ONLY:  ddzu, ddzw, drho_air, km, tend, u, v, w, rho_air_zw
     293
     294       USE control_parameters,                                                                     &
     295           ONLY:  constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
     296
     297       USE grid_variables,                                                                         &
    316298           ONLY:  ddx, ddy, ddy2
    317        
    318        USE indices,                                                            &
     299
     300       USE indices,                                                                                &
    319301           ONLY:  nzb, nzt, wall_flags_total_0
    320        
     302
    321303       USE kinds
    322304
    323        USE surface_mod,                                                        &
    324            ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
    325                    surf_usm_v
     305       USE surface_mod,                                                                            &
     306           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
    326307
    327308       IMPLICIT NONE
     
    341322       REAL(wp)     ::  kmzm          !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid
    342323       REAL(wp)     ::  kmzp          !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid
    343        REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface 
    344        REAL(wp)     ::  mask_east     !< flag to mask vertical surface south of the grid point 
    345        REAL(wp)     ::  mask_west     !< flag to mask vertical surface north of the grid point 
     324       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface
     325       REAL(wp)     ::  mask_east     !< flag to mask vertical surface south of the grid point
     326       REAL(wp)     ::  mask_west     !< flag to mask vertical surface north of the grid point
    346327       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
    347328
     
    350331       DO  k = nzb+1, nzt
    351332!
    352 !--       Predetermine flag to mask topography and wall-bounded grid points. 
    353 !--       It is sufficient to masked only east- and west-facing surfaces, which
    354 !--       need special treatment for the v-component.
    355           flag      = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   2 ) ) 
     333!--       Predetermine flag to mask topography and wall-bounded grid points.
     334!--       It is sufficient to masked only east- and west-facing surfaces, which need special
     335!--       treatment for the v-component.
     336          flag      = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   2 ) )
    356337          mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 2 ) )
    357338          mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 2 ) )
     
    361342          kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    362343
    363           tend(k,j,i) = tend(k,j,i) +          (                             &
    364                           mask_east * kmxp * (                               &
    365                                  ( v(k,j,i+1) - v(k,j,i)     ) * ddx         &
    366                                + ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy         &
    367                                              )                               &
    368                         - mask_west * kmxm * (                               &
    369                                  ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
    370                                + ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
    371                                              )                               &
    372                                                ) * ddx  * flag               &
    373                                     + 2.0_wp * (                             &
    374                                   km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i)   )  &
    375                                 - km(k,j-1,i) * ( v(k,j,i)   - v(k,j-1,i) )  &
     344          tend(k,j,i) = tend(k,j,i) +          (                                                   &
     345                          mask_east * kmxp * (                                                     &
     346                                 ( v(k,j,i+1) - v(k,j,i)     ) * ddx                               &
     347                               + ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy                               &
     348                                             )                                                     &
     349                        - mask_west * kmxm * (                                                     &
     350                                 ( v(k,j,i) - v(k,j,i-1) ) * ddx                                   &
     351                               + ( u(k,j,i) - u(k,j-1,i) ) * ddy                                   &
     352                                             )                                                     &
     353                                               ) * ddx  * flag                                     &
     354                                    + 2.0_wp * (                                                   &
     355                                  km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i)   )                        &
     356                                - km(k,j-1,i) * ( v(k,j,i)   - v(k,j-1,i) )                        &
    376357                                               ) * ddy2 * flag
    377358       ENDDO
    378359
    379360!
    380 !--    Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3)
    381 !--    surfaces. Note, in the the flat case, loops won't be entered as
    382 !--    start_index > end_index. Furtermore, note, no vertical natural surfaces
    383 !--    so far.           
     361!--    Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3) surfaces. Note, in the
     362!--    the flat case, loops won't be entered as start_index > end_index. Furtermore, note, no
     363!--    vertical natural surfaces so far.
    384364!--    Default-type surfaces
    385365       DO  l = 2, 3
     
    389369             k           = surf_def_v(l)%k(m)
    390370             tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddx
    391           ENDDO   
     371          ENDDO
    392372       ENDDO
    393373!
     
    399379             k           = surf_lsm_v(l)%k(m)
    400380             tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddx
    401           ENDDO   
     381          ENDDO
    402382       ENDDO
    403383!
     
    409389             k           = surf_usm_v(l)%k(m)
    410390             tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddx
    411           ENDDO   
    412        ENDDO
    413 !
    414 !--    Compute vertical diffusion. In case of simulating a surface layer,
    415 !--    respective grid diffusive fluxes are masked (flag 8) within this
    416 !--    loop, and added further below, else, simple gradient approach is
    417 !--    applied. Model top is also mask if top-momentum flux is given.
     391          ENDDO
     392       ENDDO
     393!
     394!--    Compute vertical diffusion. In case of simulating a surface layer, respective grid diffusive
     395!--    fluxes are masked (flag 8) within this loop, and added further below, else, simple gradient
     396!--    approach is applied. Model top is also mask if top-momentum flux is given.
    418397       DO  k = nzb+1, nzt
    419398!
    420 !--       Determine flags to mask topography below and above. Flag 2 is
    421 !--       used to mask topography in general, while flag 10 implies also
    422 !--       information about use_surface_fluxes. Flag 9 is used to control
    423 !--       momentum flux at model top. 
    424           mask_bottom = MERGE( 1.0_wp, 0.0_wp,                                 &
    425                                BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
    426           mask_top    = MERGE( 1.0_wp, 0.0_wp,                                 &
    427                                BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *     &
    428                         MERGE( 1.0_wp, 0.0_wp,                                 &
    429                                BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
    430           flag        = MERGE( 1.0_wp, 0.0_wp,                                 &
    431                                BTEST( wall_flags_total_0(k,j,i), 2 ) )
     399!--       Determine flags to mask topography below and above. Flag 2 is used to mask topography in
     400!--       general, while flag 10 implies also information about use_surface_fluxes. Flag 9 is used
     401!--       to control momentum flux at model top.
     402          mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
     403          mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *         &
     404                        MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
     405          flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
    432406!
    433407!--       Interpolate eddy diffusivities on staggered gridpoints
     
    435409          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
    436410
    437           tend(k,j,i) = tend(k,j,i)                                            &
    438                       & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
    439                       &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
    440                       &            ) * rho_air_zw(k)   * mask_top              &
    441                       &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
    442                       &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
    443                       &            ) * rho_air_zw(k-1) * mask_bottom           &
    444                       &   ) * ddzw(k) * drho_air(k) * flag
    445        ENDDO
    446 
    447 !
    448 !--    Vertical diffusion at the first grid point above the surface, if the
    449 !--    momentum flux at the bottom is given by the Prandtl law or if it is
    450 !--    prescribed by the user.
    451 !--    Difference quotient of the momentum flux is not formed over half of
    452 !--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
    453 !--    other (LES) models showed that the values of the momentum flux becomes
    454 !--    too large in this case.
     411          tend(k,j,i) = tend(k,j,i)                                                                &
     412                        & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)                       &
     413                        &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy                             &
     414                        &            ) * rho_air_zw(k)   * mask_top                                &
     415                        &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)                     &
     416                        &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy                         &
     417                        &            ) * rho_air_zw(k-1) * mask_bottom                             &
     418                        &   ) * ddzw(k) * drho_air(k) * flag
     419       ENDDO
     420
     421!
     422!--    Vertical diffusion at the first grid point above the surface, if the momentum flux at the
     423!--    bottom is given by the Prandtl law or if it is prescribed by the user.
     424!--    Difference quotient of the momentum flux is not formed over half of the grid spacing
     425!--    (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the values
     426!--    of the momentum flux becomes too large in this case.
    455427       IF ( use_surface_fluxes )  THEN
    456428!
     
    461433             k   = surf_def_h(0)%k(m)
    462434
    463              tend(k,j,i) = tend(k,j,i)                                         &
    464                         + ( - ( - surf_def_h(0)%vsws(m) )                      &
    465                           ) * ddzw(k) * drho_air(k)
     435             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_def_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
    466436          ENDDO
    467437!
     
    472442             k   = surf_def_h(1)%k(m)
    473443
    474              tend(k,j,i) = tend(k,j,i)                                         &
    475                         + ( - surf_def_h(1)%vsws(m)                            &
    476                           ) * ddzw(k) * drho_air(k)
     444             tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
    477445          ENDDO
    478446!
     
    483451             k   = surf_lsm_h%k(m)
    484452
    485              tend(k,j,i) = tend(k,j,i)                                         &
    486                         + ( - ( - surf_lsm_h%vsws(m) )                         &
    487                           ) * ddzw(k) * drho_air(k)
     453             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_lsm_h%vsws(m) ) ) * ddzw(k) * drho_air(k)
    488454
    489455          ENDDO
     
    495461             k   = surf_usm_h%k(m)
    496462
    497              tend(k,j,i) = tend(k,j,i)                                         &
    498                         + ( - ( - surf_usm_h%vsws(m) )                         &
    499                           ) * ddzw(k) * drho_air(k)
     463             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_usm_h%vsws(m) ) ) * ddzw(k) * drho_air(k)
    500464
    501465          ENDDO
     
    510474             k   = surf_def_h(2)%k(m)
    511475
    512              tend(k,j,i) = tend(k,j,i)                                         &
    513                            + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
     476             tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
    514477          ENDDO
    515478       ENDIF
Note: See TracChangeset for help on using the changeset viewer.