Changeset 4583 for palm/trunk/SOURCE/diffusion_u.f90
- Timestamp:
- Jun 29, 2020 12:36:47 PM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/diffusion_u.f90 (modified) (26 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/diffusion_u.f90
r4360 r4583 1 1 !> @file diffusion_u.f90 2 !------------------------------------------------------------------------------ !2 !--------------------------------------------------------------------------------------------------! 3 3 ! This file is part of the PALM model system. 4 4 ! 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/>. 16 15 ! 17 16 ! Copyright 1997-2020 Leibniz Universitaet Hannover 18 !------------------------------------------------------------------------------ !17 !--------------------------------------------------------------------------------------------------! 19 18 ! 20 19 ! Current revisions: … … 25 24 ! ----------------- 26 25 ! $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 ! 30 32 ! 4329 2019-12-10 15:46:36Z motisi 31 33 ! Renamed wall_flags_0 to wall_flags_static_0 32 ! 34 ! 33 35 ! 4182 2019-08-22 15:20:23Z scharf 34 36 ! Corrected "Former revisions" section 35 ! 37 ! 36 38 ! 3655 2019-01-07 16:51:22Z knoop 37 39 ! OpenACC port for SPEC … … 44 46 ! ------------ 45 47 !> Diffusion term of the u-component 46 !> @todo additional damping (needed for non-cyclic bc) causes bad vectorization 47 ! > and slows down thespeed on NEC about 5-10%48 !------------------------------------------------------------------------------ !48 !> @todo additional damping (needed for non-cyclic bc) causes bad vectorization and slows down the 49 ! speed on NEC about 5-10% 50 !--------------------------------------------------------------------------------------------------! 49 51 MODULE diffusion_u_mod 50 52 51 53 52 54 PRIVATE … … 61 63 62 64 63 !------------------------------------------------------------------------------ !65 !--------------------------------------------------------------------------------------------------! 64 66 ! Description: 65 67 ! ------------ 66 68 !> Call for all grid points 67 !------------------------------------------------------------------------------ !69 !--------------------------------------------------------------------------------------------------! 68 70 SUBROUTINE diffusion_u 69 71 70 USE arrays_3d, & 71 ONLY: ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw 72 73 USE control_parameters, & 74 ONLY: constant_top_momentumflux, use_surface_fluxes, & 75 use_top_fluxes 76 77 USE grid_variables, & 72 USE arrays_3d, & 73 ONLY: ddzu, ddzw, drho_air, km, rho_air_zw, tend, u, v, w 74 75 USE control_parameters, & 76 ONLY: constant_top_momentumflux, use_surface_fluxes, use_top_fluxes 77 78 USE grid_variables, & 78 79 ONLY: ddx, ddx2, ddy 79 80 USE indices, &80 81 USE indices, & 81 82 ONLY: nxlu, nxr, nyn, nys, nzb, nzt, wall_flags_total_0 82 83 83 84 USE kinds 84 85 85 USE surface_mod, & 86 ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & 87 surf_usm_v 86 USE surface_mod, & 87 ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v 88 88 89 89 IMPLICIT NONE … … 102 102 REAL(wp) :: kmzm !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid 103 103 REAL(wp) :: kmzp !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid 104 REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface 105 REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point 106 REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point 107 REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface 104 REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface 105 REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point 106 REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point 107 REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface 108 108 109 109 … … 125 125 DO k = nzb+1, nzt 126 126 ! 127 !-- Predetermine flag to mask topography and wall-bounded grid points. 128 !-- It is sufficient to masked only north- and south-facing surfaces, which 129 !-- need special treatment for the u-component.130 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 127 !-- Predetermine flag to mask topography and wall-bounded grid points. 128 !-- It is sufficient to masked only north- and south-facing surfaces, which need special 129 !-- treatment for the u-component. 130 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 131 131 mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 1 ) ) 132 132 mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 1 ) ) 133 133 ! 134 134 !-- Interpolate eddy diffusivities on staggered gridpoints 135 kmyp = 0.25_wp * & 136 ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) ) 137 kmym = 0.25_wp * & 138 ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) ) 139 140 tend(k,j,i) = tend(k,j,i) & 141 + 2.0_wp * ( & 142 km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) & 143 - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) & 144 ) * ddx2 * flag & 145 + ( mask_north * ( & 146 kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy & 147 + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx & 148 ) & 149 - mask_south * ( & 150 kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 151 + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx & 152 ) & 153 ) * ddy * flag 135 kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) ) 136 kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) ) 137 138 tend(k,j,i) = tend(k,j,i) & 139 + 2.0_wp * ( & 140 km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) & 141 - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) & 142 ) * ddx2 * flag & 143 + ( mask_north * ( & 144 kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy & 145 + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx & 146 ) & 147 - mask_south * ( & 148 kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 149 + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx & 150 ) & 151 ) * ddy * flag 154 152 ENDDO 155 153 ! 156 !-- Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1) 157 !-- surfaces. Note, in the the flat case, loops won't be entered as 158 !-- start_index > end_index. Furtermore, note, no vertical natural surfaces 159 !-- so far. 154 !-- Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1) surfaces. 155 !-- Note, in the the flat case, loops won't be entered as start_index > end_index. 156 !-- Furtermore, note, no vertical natural surfaces so far. 160 157 !-- Default-type surfaces 161 158 DO l = 0, 1 … … 164 161 DO m = surf_s, surf_e 165 162 k = surf_def_v(l)%k(m) 166 tend(k,j,i) = tend(k,j,i) + & 167 surf_def_v(l)%mom_flux_uv(m) * ddy 168 ENDDO 163 tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddy 164 ENDDO 169 165 ENDDO 170 166 ! … … 175 171 DO m = surf_s, surf_e 176 172 k = surf_lsm_v(l)%k(m) 177 tend(k,j,i) = tend(k,j,i) + & 178 surf_lsm_v(l)%mom_flux_uv(m) * ddy 179 ENDDO 173 tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddy 174 ENDDO 180 175 ENDDO 181 176 ! … … 186 181 DO m = surf_s, surf_e 187 182 k = surf_usm_v(l)%k(m) 188 tend(k,j,i) = tend(k,j,i) + & 189 surf_usm_v(l)%mom_flux_uv(m) * ddy 190 ENDDO 183 tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddy 184 ENDDO 191 185 ENDDO 192 186 193 187 ! 194 !-- Compute vertical diffusion. In case of simulating a surface layer, 195 !-- respective grid diffusive fluxes are masked (flag 8) within this196 !-- loop, and added further below, else, simple gradient approachis197 !-- applied. Model top is also mask if top-momentum flux is given.188 !-- Compute vertical diffusion. In case of simulating a surface layer, respective grid 189 !-- diffusive fluxes are masked (flag 8) within this loop, and added further below, else, 190 !-- simple gradient approach is applied. Model top is also mask if top-momentum flux is 191 !-- given. 198 192 DO k = nzb+1, nzt 199 193 ! 200 !-- Determine flags to mask topography below and above. Flag 1 is 201 !-- used to mask topography in general, and flag 8 implies 202 !-- information about use_surface_fluxes. Flag 9 is used to control 203 !-- momentum flux at model top. 204 mask_bottom = MERGE( 1.0_wp, 0.0_wp, & 205 BTEST( wall_flags_total_0(k-1,j,i), 8 ) ) 206 mask_top = MERGE( 1.0_wp, 0.0_wp, & 207 BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) * & 208 MERGE( 1.0_wp, 0.0_wp, & 209 BTEST( wall_flags_total_0(k+1,j,i), 9 ) ) 210 flag = MERGE( 1.0_wp, 0.0_wp, & 211 BTEST( wall_flags_total_0(k,j,i), 1 ) ) 194 !-- Determine flags to mask topography below and above. Flag 1 is used to mask 195 !-- topography in general, and flag 8 implies information about use_surface_fluxes. 196 !-- Flag 9 is used to control momentum flux at model top. 197 mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) ) 198 mask_top = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) * & 199 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) ) 200 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 212 201 ! 213 202 !-- Interpolate eddy diffusivities on staggered gridpoints 214 kmzp = 0.25_wp * & 215 ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) ) 216 kmzm = 0.25_wp * & 217 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) ) 218 219 tend(k,j,i) = tend(k,j,i) & 220 + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 221 + ( w(k,j,i) - w(k,j,i-1) ) * ddx & 222 ) * rho_air_zw(k) * mask_top & 223 - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) & 224 + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx & 225 ) * rho_air_zw(k-1) * mask_bottom & 226 ) * ddzw(k) * drho_air(k) * flag 203 kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) ) 204 kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) ) 205 206 tend(k,j,i) = tend(k,j,i) & 207 + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 208 + ( w(k,j,i) - w(k,j,i-1) ) * ddx & 209 ) * rho_air_zw(k) * mask_top & 210 - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) & 211 + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx & 212 ) * rho_air_zw(k-1) * mask_bottom & 213 ) * ddzw(k) * drho_air(k) * flag 227 214 ENDDO 228 215 229 216 ! 230 !-- Vertical diffusion at the first grid point above the surface, 231 !-- if the momentum flux at the bottom is given by the Prandtl law or 232 !-- if it is prescribed by the user. 233 !-- Difference quotient of the momentum flux is not formed over half 234 !-- of the grid spacing (2.0*ddzw(k)) any more, since the comparison 235 !-- with other (LES) models showed that the values of the momentum 236 !-- flux becomes too large in this case. 237 !-- The term containing w(k-1,..) (see above equation) is removed here 238 !-- because the vertical velocity is assumed to be zero at the surface. 217 !-- Vertical diffusion at the first grid point above the surface, if the momentum flux at 218 !-- the bottom is given by the Prandtl law or if it is prescribed by the user. 219 !-- Difference quotient of the momentum flux is not formed over half of the grid spacing 220 !-- (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the 221 !-- values of the momentum flux becomes too large in this case. 222 !-- The term containing w(k-1,..) (see above equation) is removed here because the vertical 223 !-- velocity is assumed to be zero at the surface. 239 224 IF ( use_surface_fluxes ) THEN 240 225 ! … … 246 231 k = surf_def_h(0)%k(m) 247 232 248 tend(k,j,i) = tend(k,j,i) & 249 + ( - ( - surf_def_h(0)%usws(m) ) & 250 ) * ddzw(k) * drho_air(k) 233 tend(k,j,i) = tend(k,j,i) & 234 + ( - ( - surf_def_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k) 251 235 ENDDO 252 236 ! … … 258 242 k = surf_def_h(1)%k(m) 259 243 260 tend(k,j,i) = tend(k,j,i) & 261 + ( - surf_def_h(1)%usws(m) & 262 ) * ddzw(k) * drho_air(k) 244 tend(k,j,i) = tend(k,j,i) & 245 + ( - surf_def_h(1)%usws(m) ) * ddzw(k) * drho_air(k) 263 246 ENDDO 264 247 ! … … 270 253 k = surf_lsm_h%k(m) 271 254 272 tend(k,j,i) = tend(k,j,i) & 273 + ( - ( - surf_lsm_h%usws(m) ) & 274 ) * ddzw(k) * drho_air(k) 255 tend(k,j,i) = tend(k,j,i) & 256 + ( - ( - surf_lsm_h%usws(m) ) ) * ddzw(k) * drho_air(k) 275 257 ENDDO 276 258 ! … … 282 264 k = surf_usm_h%k(m) 283 265 284 tend(k,j,i) = tend(k,j,i) & 285 + ( - ( - surf_usm_h%usws(m) ) & 286 ) * ddzw(k) * drho_air(k) 266 tend(k,j,i) = tend(k,j,i) & 267 + ( - ( - surf_usm_h%usws(m) ) ) * ddzw(k) * drho_air(k) 287 268 ENDDO 288 269 … … 297 278 k = surf_def_h(2)%k(m) 298 279 299 tend(k,j,i) = tend(k,j,i) &300 + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)280 tend(k,j,i) = tend(k,j,i) & 281 + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k) 301 282 ENDDO 302 283 ENDIF … … 308 289 309 290 310 !------------------------------------------------------------------------------ !291 !--------------------------------------------------------------------------------------------------! 311 292 ! Description: 312 293 ! ------------ 313 294 !> Call for grid point i,j 314 !------------------------------------------------------------------------------ !295 !--------------------------------------------------------------------------------------------------! 315 296 SUBROUTINE diffusion_u_ij( i, j ) 316 297 317 USE arrays_3d, & 318 ONLY: ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw 319 320 USE control_parameters, & 321 ONLY: constant_top_momentumflux, use_surface_fluxes, & 322 use_top_fluxes 323 324 USE grid_variables, & 298 USE arrays_3d, & 299 ONLY: ddzu, ddzw, drho_air, km, tend, u, v, w, rho_air_zw 300 301 USE control_parameters, & 302 ONLY: constant_top_momentumflux, use_surface_fluxes, use_top_fluxes 303 304 USE grid_variables, & 325 305 ONLY: ddx, ddx2, ddy 326 327 USE indices, &306 307 USE indices, & 328 308 ONLY: nzb, nzt, wall_flags_total_0 329 309 330 310 USE kinds 331 311 332 USE surface_mod, & 333 ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & 334 surf_usm_v 312 USE surface_mod, & 313 ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v 335 314 336 315 IMPLICIT NONE … … 349 328 REAL(wp) :: kmzm !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid 350 329 REAL(wp) :: kmzp !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid 351 REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface 352 REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point 353 REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point 354 REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface 355 ! 330 REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface 331 REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point 332 REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point 333 REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface 334 ! 356 335 !-- Compute horizontal diffusion 357 336 DO k = nzb+1, nzt 358 337 ! 359 !-- Predetermine flag to mask topography and wall-bounded grid points. 360 !-- It is sufficient to masked only north- and south-facing surfaces, which 361 !-- need special treatment for the u-component.362 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 338 !-- Predetermine flag to mask topography and wall-bounded grid points. 339 !-- It is sufficient to masked only north- and south-facing surfaces, which need special 340 !-- treatment for the u-component. 341 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 363 342 mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 1 ) ) 364 343 mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 1 ) ) … … 368 347 kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) ) 369 348 370 tend(k,j,i) = tend(k,j,i) & 371 + 2.0_wp * ( & 372 km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) & 373 - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) & 374 ) * ddx2 * flag & 375 + ( & 376 mask_north * kmyp * ( ( u(k,j+1,i) - u(k,j,i) ) * ddy & 377 + ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx & 378 ) & 379 - mask_south * kmym * ( ( u(k,j,i) - u(k,j-1,i) ) * ddy & 380 + ( v(k,j,i) - v(k,j,i-1) ) * ddx & 381 ) & 382 ) * ddy * flag 383 ENDDO 384 385 ! 386 !-- Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1) 387 !-- surfaces. Note, in the the flat case, loops won't be entered as 388 !-- start_index > end_index. Furtermore, note, no vertical natural surfaces 389 !-- so far. 349 tend(k,j,i) = tend(k,j,i) & 350 + 2.0_wp * ( & 351 km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) & 352 - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) & 353 ) * ddx2 * flag & 354 + ( & 355 mask_north * kmyp * ( ( u(k,j+1,i) - u(k,j,i) ) * ddy & 356 + ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx & 357 ) & 358 - mask_south * kmym * ( ( u(k,j,i) - u(k,j-1,i) ) * ddy & 359 + ( v(k,j,i) - v(k,j,i-1) ) * ddx & 360 ) & 361 ) * ddy * flag 362 ENDDO 363 364 ! 365 !-- Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1) surfaces. Note, in 366 !-- the the flat case, loops won't be entered as start_index > end_index. Furtermore, note, no 367 !-- vertical natural surfaces so far. 390 368 !-- Default-type surfaces 391 369 DO l = 0, 1 … … 395 373 k = surf_def_v(l)%k(m) 396 374 tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddy 397 ENDDO 375 ENDDO 398 376 ENDDO 399 377 ! … … 405 383 k = surf_lsm_v(l)%k(m) 406 384 tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddy 407 ENDDO 385 ENDDO 408 386 ENDDO 409 387 ! … … 415 393 k = surf_usm_v(l)%k(m) 416 394 tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddy 417 ENDDO 418 ENDDO 419 ! 420 !-- Compute vertical diffusion. In case of simulating a surface layer, 421 !-- respective grid diffusive fluxes are masked (flag 8) within this 422 !-- loop, and added further below, else, simple gradient approach is 423 !-- applied. Model top is also mask if top-momentum flux is given. 395 ENDDO 396 ENDDO 397 ! 398 !-- Compute vertical diffusion. In case of simulating a surface layer, respective grid diffusive 399 !-- fluxes are masked (flag 8) within this loop, and added further below, else, simple gradient 400 !-- approach is applied. Model top is also mask if top-momentum flux is given. 424 401 DO k = nzb+1, nzt 425 402 ! 426 !-- Determine flags to mask topography below and above. Flag 1 is 427 !-- used to mask topography in general, and flag 8 implies 428 !-- information about use_surface_fluxes. Flag 9 is used to control 429 !-- momentum flux at model top. 430 mask_bottom = MERGE( 1.0_wp, 0.0_wp, & 431 BTEST( wall_flags_total_0(k-1,j,i), 8 ) ) 432 mask_top = MERGE( 1.0_wp, 0.0_wp, & 433 BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) * & 434 MERGE( 1.0_wp, 0.0_wp, & 435 BTEST( wall_flags_total_0(k+1,j,i), 9 ) ) 436 flag = MERGE( 1.0_wp, 0.0_wp, & 437 BTEST( wall_flags_total_0(k,j,i), 1 ) ) 403 !-- Determine flags to mask topography below and above. Flag 1 is used to mask topography in 404 !-- general, and flag 8 implies information about use_surface_fluxes. Flag 9 is used to 405 !-- control momentum flux at model top. 406 mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) ) 407 mask_top = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) * & 408 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) ) 409 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 438 410 ! 439 411 !-- Interpolate eddy diffusivities on staggered gridpoints … … 441 413 kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) ) 442 414 443 tend(k,j,i) = tend(k,j,i) &444 + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &445 + ( w(k,j,i) - w(k,j,i-1) ) * ddx &446 ) * rho_air_zw(k) * mask_top &447 - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &448 + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &449 ) * rho_air_zw(k-1) * mask_bottom &415 tend(k,j,i) = tend(k,j,i) & 416 + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 417 + ( w(k,j,i) - w(k,j,i-1) ) * ddx & 418 ) * rho_air_zw(k) * mask_top & 419 - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) & 420 + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx & 421 ) * rho_air_zw(k-1) * mask_bottom & 450 422 ) * ddzw(k) * drho_air(k) * flag 451 423 ENDDO 452 424 453 425 ! 454 !-- Vertical diffusion at the first surface grid points, if the 455 !-- momentum flux at the bottom is given by the Prandtl law or if it is 456 !-- prescribed by the user. 457 !-- Difference quotient of the momentum flux is not formed over half of 458 !-- the grid spacing (2.0*ddzw(k)) any more, since the comparison with 459 !-- other (LES) models showed that the values of the momentum flux becomes 460 !-- too large in this case. 426 !-- Vertical diffusion at the first surface grid points, if the momentum flux at the bottom is 427 !-- given by the Prandtl law or if it is prescribed by the user. 428 !-- Difference quotient of the momentum flux is not formed over half of the grid spacing 429 !-- (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the values 430 !-- of the momentum flux becomes too large in this case. 461 431 IF ( use_surface_fluxes ) THEN 462 432 ! … … 468 438 k = surf_def_h(0)%k(m) 469 439 470 tend(k,j,i) = tend(k,j,i) & 471 + ( - ( - surf_def_h(0)%usws(m) ) & 472 ) * ddzw(k) * drho_air(k) 440 tend(k,j,i) = tend(k,j,i) + ( - ( - surf_def_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k) 473 441 ENDDO 474 442 ! … … 480 448 k = surf_def_h(1)%k(m) 481 449 482 tend(k,j,i) = tend(k,j,i) & 483 + ( - surf_def_h(1)%usws(m) & 484 ) * ddzw(k) * drho_air(k) 450 tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(1)%usws(m) ) * ddzw(k) * drho_air(k) 485 451 ENDDO 486 452 ! … … 492 458 k = surf_lsm_h%k(m) 493 459 494 tend(k,j,i) = tend(k,j,i) & 495 + ( - ( - surf_lsm_h%usws(m) ) & 496 ) * ddzw(k) * drho_air(k) 460 tend(k,j,i) = tend(k,j,i) + ( - ( - surf_lsm_h%usws(m) ) ) * ddzw(k) * drho_air(k) 497 461 ENDDO 498 462 ! … … 504 468 k = surf_usm_h%k(m) 505 469 506 tend(k,j,i) = tend(k,j,i) & 507 + ( - ( - surf_usm_h%usws(m) ) & 508 ) * ddzw(k) * drho_air(k) 470 tend(k,j,i) = tend(k,j,i) + ( - ( - surf_usm_h%usws(m) ) ) * ddzw(k) * drho_air(k) 509 471 ENDDO 510 472 … … 519 481 k = surf_def_h(2)%k(m) 520 482 521 tend(k,j,i) = tend(k,j,i) & 522 + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k) 483 tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k) 523 484 ENDDO 524 485 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
