Changeset 4583 for palm/trunk/SOURCE/diffusion_w.f90
- Timestamp:
- Jun 29, 2020 12:36:47 PM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/diffusion_w.f90 (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/diffusion_w.f90
r4360 r4583 1 1 !> @file diffusion_w.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 w-component 46 !------------------------------------------------------------------------------ !48 !--------------------------------------------------------------------------------------------------! 47 49 MODULE diffusion_w_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_w 67 69 68 USE arrays_3d, &69 ONLY : ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air70 71 USE grid_variables, &70 USE arrays_3d, & 71 ONLY : ddzu, ddzw, drho_air_zw, km, rho_air, tend, u, v, w 72 73 USE grid_variables, & 72 74 ONLY : ddx, ddy 73 74 USE indices, &75 76 USE indices, & 75 77 ONLY : nxl, nxr, nyn, nys, nzb, nzt, wall_flags_total_0 76 78 77 79 USE kinds 78 80 79 USE surface_mod, &81 USE surface_mod, & 80 82 ONLY : surf_def_v, surf_lsm_v, surf_usm_v 81 83 … … 89 91 INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint 90 92 INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint 91 93 92 94 REAL(wp) :: flag !< flag to mask topography grid points 93 95 REAL(wp) :: kmxm !< diffusion coefficient on leftward side of the w-gridbox - interpolated onto xu-y grid … … 95 97 REAL(wp) :: kmym !< diffusion coefficient on southward side of the w-gridbox - interpolated onto x-yv grid 96 98 REAL(wp) :: kmyp !< diffusion coefficient on northward side of the w-gridbox - interpolated onto x-yv grid 99 REAL(wp) :: mask_east !< flag to mask vertical wall east of the grid point 100 REAL(wp) :: mask_north !< flag to mask vertical wall north of the grid point 101 REAL(wp) :: mask_south !< flag to mask vertical wall south of the grid point 97 102 REAL(wp) :: mask_west !< flag to mask vertical wall west of the grid point 98 REAL(wp) :: mask_east !< flag to mask vertical wall east of the grid point99 REAL(wp) :: mask_south !< flag to mask vertical wall south of the grid point100 REAL(wp) :: mask_north !< flag to mask vertical wall north of the grid point101 103 102 104 … … 116 118 DO k = nzb+1, nzt-1 117 119 ! 118 !-- Predetermine flag to mask topography and wall-bounded grid points. 119 flag = MERGE( 1.0_wp, 0.0_wp, & 120 BTEST( wall_flags_total_0(k,j,i), 3 ) ) 121 mask_east = MERGE( 1.0_wp, 0.0_wp, & 122 BTEST( wall_flags_total_0(k,j,i+1), 3 ) ) 123 mask_west = MERGE( 1.0_wp, 0.0_wp, & 124 BTEST( wall_flags_total_0(k,j,i-1), 3 ) ) 125 mask_south = MERGE( 1.0_wp, 0.0_wp, & 126 BTEST( wall_flags_total_0(k,j-1,i), 3 ) ) 127 mask_north = MERGE( 1.0_wp, 0.0_wp, & 128 BTEST( wall_flags_total_0(k,j+1,i), 3 ) ) 120 !-- Predetermine flag to mask topography and wall-bounded grid points. 121 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) ) 122 mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 3 ) ) 123 mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 3 ) ) 124 mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 3 ) ) 125 mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 3 ) ) 129 126 ! 130 127 !-- Interpolate eddy diffusivities on staggered gridpoints 131 kmxp = 0.25_wp * ( km(k,j,i) + km(k,j,i+1) + &128 kmxp = 0.25_wp * ( km(k,j,i) + km(k,j,i+1) + & 132 129 km(k+1,j,i) + km(k+1,j,i+1) ) 133 kmxm = 0.25_wp * ( km(k,j,i) + km(k,j,i-1) + &130 kmxm = 0.25_wp * ( km(k,j,i) + km(k,j,i-1) + & 134 131 km(k+1,j,i) + km(k+1,j,i-1) ) 135 kmyp = 0.25_wp * ( km(k,j,i) + km(k+1,j,i) + &132 kmyp = 0.25_wp * ( km(k,j,i) + km(k+1,j,i) + & 136 133 km(k,j+1,i) + km(k+1,j+1,i) ) 137 kmym = 0.25_wp * ( km(k,j,i) + km(k+1,j,i) + &134 kmym = 0.25_wp * ( km(k,j,i) + km(k+1,j,i) + & 138 135 km(k,j-1,i) + km(k+1,j-1,i) ) 139 136 140 tend(k,j,i) = tend(k,j,i) & 141 + ( mask_east * kmxp * ( & 142 ( w(k,j,i+1) - w(k,j,i) ) * ddx & 143 + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) & 144 ) & 145 - mask_west * kmxm * ( & 146 ( w(k,j,i) - w(k,j,i-1) ) * ddx & 147 + ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 148 ) & 149 ) * ddx * flag & 150 + ( mask_north * kmyp * ( & 151 ( w(k,j+1,i) - w(k,j,i) ) * ddy & 152 + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) & 153 ) & 154 - mask_south * kmym * ( & 155 ( w(k,j,i) - w(k,j-1,i) ) * ddy & 156 + ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 157 ) & 158 ) * ddy * flag & 159 + 2.0_wp * ( & 160 km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 161 * rho_air(k+1) & 162 - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 163 * rho_air(k) & 164 ) * ddzu(k+1) * drho_air_zw(k) * flag 165 ENDDO 166 167 ! 168 !-- Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1) 169 !-- surfaces. Note, in the the flat case, loops won't be entered as 170 !-- start_index > end_index. Furtermore, note, no vertical natural surfaces 171 !-- so far. 137 tend(k,j,i) = tend(k,j,i) & 138 + ( mask_east * kmxp * ( & 139 ( w(k,j,i+1) - w(k,j,i) ) * ddx & 140 + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) & 141 ) & 142 - mask_west * kmxm * ( & 143 ( w(k,j,i) - w(k,j,i-1) ) * ddx & 144 + ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 145 ) & 146 ) * ddx * flag & 147 + ( mask_north * kmyp * ( & 148 ( w(k,j+1,i) - w(k,j,i) ) * ddy & 149 + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) & 150 ) & 151 - mask_south * kmym * ( & 152 ( w(k,j,i) - w(k,j-1,i) ) * ddy & 153 + ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 154 ) & 155 ) * ddy * flag & 156 + 2.0_wp * ( & 157 km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 158 * rho_air(k+1) & 159 - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 160 * rho_air(k) & 161 ) * ddzu(k+1) * drho_air_zw(k) * flag 162 ENDDO 163 164 ! 165 !-- Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1) surfaces. 166 !-- Note, in the the flat case, loops won't be entered as start_index > end_index. 167 !-- Furtermore, note, no vertical natural surfaces so far. 172 168 !-- Default-type surfaces 173 169 DO l = 0, 1 … … 176 172 DO m = surf_s, surf_e 177 173 k = surf_def_v(l)%k(m) 178 tend(k,j,i) = tend(k,j,i) + & 179 surf_def_v(l)%mom_flux_w(m) * ddy 180 ENDDO 174 tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddy 175 ENDDO 181 176 ENDDO 182 177 ! … … 187 182 DO m = surf_s, surf_e 188 183 k = surf_lsm_v(l)%k(m) 189 tend(k,j,i) = tend(k,j,i) + & 190 surf_lsm_v(l)%mom_flux_w(m) * ddy 191 ENDDO 184 tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddy 185 ENDDO 192 186 ENDDO 193 187 ! … … 198 192 DO m = surf_s, surf_e 199 193 k = surf_usm_v(l)%k(m) 200 tend(k,j,i) = tend(k,j,i) + & 201 surf_usm_v(l)%mom_flux_w(m) * ddy 202 ENDDO 203 ENDDO 204 ! 205 !-- Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3) 206 !-- surface. 194 tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddy 195 ENDDO 196 ENDDO 197 ! 198 !-- Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3) surface. 207 199 !-- Default-type surfaces 208 200 DO l = 2, 3 … … 211 203 DO m = surf_s, surf_e 212 204 k = surf_def_v(l)%k(m) 213 tend(k,j,i) = tend(k,j,i) + & 214 surf_def_v(l)%mom_flux_w(m) * ddx 215 ENDDO 205 tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddx 206 ENDDO 216 207 ENDDO 217 208 ! … … 222 213 DO m = surf_s, surf_e 223 214 k = surf_lsm_v(l)%k(m) 224 tend(k,j,i) = tend(k,j,i) + & 225 surf_lsm_v(l)%mom_flux_w(m) * ddx 226 ENDDO 215 tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddx 216 ENDDO 227 217 ENDDO 228 218 ! … … 233 223 DO m = surf_s, surf_e 234 224 k = surf_usm_v(l)%k(m) 235 tend(k,j,i) = tend(k,j,i) + & 236 surf_usm_v(l)%mom_flux_w(m) * ddx 237 ENDDO 225 tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddx 226 ENDDO 238 227 ENDDO 239 228 … … 244 233 245 234 246 !------------------------------------------------------------------------------ !235 !--------------------------------------------------------------------------------------------------! 247 236 ! Description: 248 237 ! ------------ 249 238 !> Call for grid point i,j 250 !------------------------------------------------------------------------------ !239 !--------------------------------------------------------------------------------------------------! 251 240 SUBROUTINE diffusion_w_ij( i, j ) 252 241 253 USE arrays_3d, &254 ONLY : ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air255 256 USE grid_variables, &242 USE arrays_3d, & 243 ONLY : ddzu, ddzw, drho_air_zw, km, rho_air, tend, u, v, w 244 245 USE grid_variables, & 257 246 ONLY : ddx, ddy 258 259 USE indices, &247 248 USE indices, & 260 249 ONLY : nzb, nzt, wall_flags_total_0 261 250 262 251 USE kinds 263 252 264 USE surface_mod, &253 USE surface_mod, & 265 254 ONLY : surf_def_v, surf_lsm_v, surf_usm_v 266 255 … … 275 264 INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint 276 265 INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint 277 266 278 267 REAL(wp) :: flag !< flag to mask topography grid points 279 268 REAL(wp) :: kmxm !< diffusion coefficient on leftward side of the w-gridbox - interpolated onto xu-y grid … … 281 270 REAL(wp) :: kmym !< diffusion coefficient on southward side of the w-gridbox - interpolated onto x-yv grid 282 271 REAL(wp) :: kmyp !< diffusion coefficient on northward side of the w-gridbox - interpolated onto x-yv grid 272 REAL(wp) :: mask_east !< flag to mask vertical wall east of the grid point 273 REAL(wp) :: mask_north !< flag to mask vertical wall north of the grid point 274 REAL(wp) :: mask_south !< flag to mask vertical wall south of the grid point 283 275 REAL(wp) :: mask_west !< flag to mask vertical wall west of the grid point 284 REAL(wp) :: mask_east !< flag to mask vertical wall east of the grid point285 REAL(wp) :: mask_south !< flag to mask vertical wall south of the grid point286 REAL(wp) :: mask_north !< flag to mask vertical wall north of the grid point287 276 288 277 289 278 DO k = nzb+1, nzt-1 290 279 ! 291 !-- Predetermine flag to mask topography and wall-bounded grid points. 292 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) ) 280 !-- Predetermine flag to mask topography and wall-bounded grid points. 281 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) ) 293 282 mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 3 ) ) 294 283 mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 3 ) ) … … 302 291 kmym = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 303 292 304 tend(k,j,i) = tend(k,j,i) & 305 + ( mask_east * kmxp * ( & 306 ( w(k,j,i+1) - w(k,j,i) ) * ddx & 307 + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) & 308 ) & 309 - mask_west * kmxm * ( & 310 ( w(k,j,i) - w(k,j,i-1) ) * ddx & 311 + ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 312 ) & 313 ) * ddx * flag & 314 + ( mask_north * kmyp * ( & 315 ( w(k,j+1,i) - w(k,j,i) ) * ddy & 316 + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) & 317 ) & 318 - mask_south * kmym * ( & 319 ( w(k,j,i) - w(k,j-1,i) ) * ddy & 320 + ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 321 ) & 322 ) * ddy * flag & 323 + 2.0_wp * ( & 324 km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 325 * rho_air(k+1) & 326 - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 327 * rho_air(k) & 328 ) * ddzu(k+1) * drho_air_zw(k) * flag 329 ENDDO 330 ! 331 !-- Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1) 332 !-- surfaces. Note, in the the flat case, loops won't be entered as 333 !-- start_index > end_index. Furtermore, note, no vertical natural surfaces 334 !-- so far. 293 tend(k,j,i) = tend(k,j,i) & 294 + ( mask_east * kmxp * ( & 295 ( w(k,j,i+1) - w(k,j,i) ) * ddx & 296 + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) & 297 ) & 298 - mask_west * kmxm * ( & 299 ( w(k,j,i) - w(k,j,i-1) ) * ddx & 300 + ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 301 ) & 302 ) * ddx * flag & 303 + ( mask_north * kmyp * ( & 304 ( w(k,j+1,i) - w(k,j,i) ) * ddy & 305 + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) & 306 ) & 307 - mask_south * kmym * ( & 308 ( w(k,j,i) - w(k,j-1,i) ) * ddy & 309 + ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 310 ) & 311 ) * ddy * flag & 312 + 2.0_wp * ( & 313 km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 314 * rho_air(k+1) & 315 - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 316 * rho_air(k) & 317 ) * ddzu(k+1) * drho_air_zw(k) * flag 318 ENDDO 319 ! 320 !-- Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1) surfaces. Note, in 321 !-- the the flat case, loops won't be entered as start_index > end_index. Furtermore, note, no 322 !-- vertical natural surfaces so far. 335 323 !-- Default-type surfaces 336 324 DO l = 0, 1 … … 339 327 DO m = surf_s, surf_e 340 328 k = surf_def_v(l)%k(m) 341 tend(k,j,i) = tend(k,j,i) + & 342 surf_def_v(l)%mom_flux_w(m) * ddy 343 ENDDO 329 tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddy 330 ENDDO 344 331 ENDDO 345 332 ! … … 350 337 DO m = surf_s, surf_e 351 338 k = surf_lsm_v(l)%k(m) 352 tend(k,j,i) = tend(k,j,i) + & 353 surf_lsm_v(l)%mom_flux_w(m) * ddy 354 ENDDO 339 tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddy 340 ENDDO 355 341 ENDDO 356 342 ! … … 361 347 DO m = surf_s, surf_e 362 348 k = surf_usm_v(l)%k(m) 363 tend(k,j,i) = tend(k,j,i) + & 364 surf_usm_v(l)%mom_flux_w(m) * ddy 365 ENDDO 366 ENDDO 367 ! 368 !-- Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3) 369 !-- surfaces. 349 tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddy 350 ENDDO 351 ENDDO 352 ! 353 !-- Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3) surfaces. 370 354 !-- Default-type surfaces 371 355 DO l = 2, 3 … … 374 358 DO m = surf_s, surf_e 375 359 k = surf_def_v(l)%k(m) 376 tend(k,j,i) = tend(k,j,i) + & 377 surf_def_v(l)%mom_flux_w(m) * ddx 378 ENDDO 360 tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddx 361 ENDDO 379 362 ENDDO 380 363 ! … … 385 368 DO m = surf_s, surf_e 386 369 k = surf_lsm_v(l)%k(m) 387 tend(k,j,i) = tend(k,j,i) + & 388 surf_lsm_v(l)%mom_flux_w(m) * ddx 389 ENDDO 370 tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddx 371 ENDDO 390 372 ENDDO 391 373 ! … … 396 378 DO m = surf_s, surf_e 397 379 k = surf_usm_v(l)%k(m) 398 tend(k,j,i) = tend(k,j,i) + & 399 surf_usm_v(l)%mom_flux_w(m) * ddx 400 ENDDO 380 tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddx 381 ENDDO 401 382 ENDDO 402 383
Note: See TracChangeset
for help on using the changeset viewer.
