Changeset 4583 for palm/trunk/SOURCE/diffusion_v.f90
- Timestamp:
- Jun 29, 2020 12:36:47 PM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/diffusion_v.f90 (modified) (28 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/diffusion_v.f90
r4360 r4583 1 1 !> @file diffusion_v.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 v-component 46 !------------------------------------------------------------------------------ !48 !--------------------------------------------------------------------------------------------------! 47 49 MODULE diffusion_v_mod 48 50 49 51 50 52 PRIVATE … … 59 61 60 62 61 !------------------------------------------------------------------------------ !63 !--------------------------------------------------------------------------------------------------! 62 64 ! Description: 63 65 ! ------------ 64 66 !> Call for all grid points 65 !------------------------------------------------------------------------------ !67 !--------------------------------------------------------------------------------------------------! 66 68 SUBROUTINE diffusion_v 67 69 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, & 76 77 ONLY: ddx, ddy, ddy2 77 78 USE indices, &78 79 USE indices, & 79 80 ONLY: nxl, nxr, nyn, nysv, nzb, nzt, wall_flags_total_0 80 81 81 82 USE kinds 82 83 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 86 86 87 87 IMPLICIT NONE … … 100 100 REAL(wp) :: kmzm !< diffusion coefficient on bottom of the gridbox - interpolated onto yv-zw grid 101 101 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 106 106 107 107 !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) & … … 122 122 123 123 ! 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 ) ) 128 128 mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 2 ) ) 129 129 mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 2 ) ) … … 133 133 kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 134 134 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 * flag135 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 149 149 150 150 ENDDO 151 151 152 152 ! 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. 157 156 !-- Default-type surfaces 158 157 DO l = 2, 3 … … 161 160 DO m = surf_s, surf_e 162 161 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 166 164 ENDDO 167 165 ! … … 172 170 DO m = surf_s, surf_e 173 171 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 177 174 ENDDO 178 175 ! … … 183 180 DO m = surf_s, surf_e 184 181 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 188 184 ENDDO 189 185 ! 190 !-- Compute vertical diffusion. In case of simulating a surface layer, 191 !-- respective grid diffusive fluxes are masked (flag 10) within this192 !-- loop, and added further below, else, simple gradient approachis193 !-- 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. 194 190 DO k = nzb+1, nzt 195 191 ! 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 ) ) 208 199 ! 209 200 !-- 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 223 212 ENDDO 224 213 225 214 ! 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. 233 220 IF ( use_surface_fluxes ) THEN 234 221 ! … … 239 226 k = surf_def_h(0)%k(m) 240 227 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) 244 230 ENDDO 245 231 ! … … 250 236 k = surf_def_h(1)%k(m) 251 237 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) 255 240 ENDDO 256 241 ! … … 261 246 k = surf_lsm_h%k(m) 262 247 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) 266 250 267 251 ENDDO … … 273 257 k = surf_usm_h%k(m) 274 258 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) 278 261 279 262 ENDDO … … 288 271 k = surf_def_h(2)%k(m) 289 272 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) 292 275 ENDDO 293 276 ENDIF … … 299 282 300 283 301 !------------------------------------------------------------------------------ !284 !--------------------------------------------------------------------------------------------------! 302 285 ! Description: 303 286 ! ------------ 304 287 !> Call for grid point i,j 305 !------------------------------------------------------------------------------ !288 !--------------------------------------------------------------------------------------------------! 306 289 SUBROUTINE diffusion_v_ij( i, j ) 307 290 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, & 316 298 ONLY: ddx, ddy, ddy2 317 318 USE indices, &299 300 USE indices, & 319 301 ONLY: nzb, nzt, wall_flags_total_0 320 302 321 303 USE kinds 322 304 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 326 307 327 308 IMPLICIT NONE … … 341 322 REAL(wp) :: kmzm !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid 342 323 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 346 327 REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface 347 328 … … 350 331 DO k = nzb+1, nzt 351 332 ! 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 ) ) 356 337 mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 2 ) ) 357 338 mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 2 ) ) … … 361 342 kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 362 343 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) ) & 376 357 ) * ddy2 * flag 377 358 ENDDO 378 359 379 360 ! 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. 384 364 !-- Default-type surfaces 385 365 DO l = 2, 3 … … 389 369 k = surf_def_v(l)%k(m) 390 370 tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddx 391 ENDDO 371 ENDDO 392 372 ENDDO 393 373 ! … … 399 379 k = surf_lsm_v(l)%k(m) 400 380 tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddx 401 ENDDO 381 ENDDO 402 382 ENDDO 403 383 ! … … 409 389 k = surf_usm_v(l)%k(m) 410 390 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. 418 397 DO k = nzb+1, nzt 419 398 ! 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 ) ) 432 406 ! 433 407 !-- Interpolate eddy diffusivities on staggered gridpoints … … 435 409 kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 436 410 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. 455 427 IF ( use_surface_fluxes ) THEN 456 428 ! … … 461 433 k = surf_def_h(0)%k(m) 462 434 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) 466 436 ENDDO 467 437 ! … … 472 442 k = surf_def_h(1)%k(m) 473 443 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) 477 445 ENDDO 478 446 ! … … 483 451 k = surf_lsm_h%k(m) 484 452 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) 488 454 489 455 ENDDO … … 495 461 k = surf_usm_h%k(m) 496 462 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) 500 464 501 465 ENDDO … … 510 474 k = surf_def_h(2)%k(m) 511 475 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) 514 477 ENDDO 515 478 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
