Changeset 2232 for palm/trunk/SOURCE/diffusion_w.f90
- Timestamp:
- May 30, 2017 5:47:52 PM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/diffusion_w.f90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/diffusion_w.f90
r2119 r2232 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Adjustments to new topography and surface concept 23 23 ! 24 24 ! Former revisions: … … 97 97 98 98 99 USE wall_fluxes_mod, &100 ONLY : wall_fluxes101 102 99 PRIVATE 103 100 PUBLIC diffusion_w … … 125 122 126 123 USE grid_variables, & 127 ONLY : ddx, ddy , fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y124 ONLY : ddx, ddy 128 125 129 126 USE indices, & 130 ONLY : nxl, nxr, nyn, nys, nzb, nz b_w_inner, nzb_w_outer, nzt127 ONLY : nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0 131 128 132 129 USE kinds 133 130 131 USE surface_mod, & 132 ONLY : surf_def_v, surf_lsm_v, surf_usm_v 133 134 134 IMPLICIT NONE 135 135 136 INTEGER(iwp) :: i !< 137 INTEGER(iwp) :: j !< 138 INTEGER(iwp) :: k !< 136 INTEGER(iwp) :: i !< running index x direction 137 INTEGER(iwp) :: j !< running index y direction 138 INTEGER(iwp) :: k !< running index z direction 139 INTEGER(iwp) :: l !< running index of surface type, south- or north-facing wall 140 INTEGER(iwp) :: m !< running index surface elements 141 INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint 142 INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint 139 143 140 REAL(wp) :: kmxm !< 141 REAL(wp) :: kmxp !< 142 REAL(wp) :: kmym !< 143 REAL(wp) :: kmyp !< 144 145 REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wsus !< 146 REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wsvs !< 147 148 149 ! 150 !-- First calculate horizontal momentum flux w'u' and/or w'v' at vertical 151 !-- walls, if neccessary 152 IF ( topography /= 'flat' ) THEN 153 CALL wall_fluxes( wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, nzb_w_inner, & 154 nzb_w_outer, wall_w_x ) 155 CALL wall_fluxes( wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, nzb_w_inner, & 156 nzb_w_outer, wall_w_y ) 157 ENDIF 144 REAL(wp) :: flag !< flag to mask topography grid points 145 REAL(wp) :: kmxm !< 146 REAL(wp) :: kmxp !< 147 REAL(wp) :: kmym !< 148 REAL(wp) :: kmyp !< 149 REAL(wp) :: mask_west !< flag to mask vertical wall west of the grid point 150 REAL(wp) :: mask_east !< flag to mask vertical wall east of the grid point 151 REAL(wp) :: mask_south !< flag to mask vertical wall south of the grid point 152 REAL(wp) :: mask_north !< flag to mask vertical wall north of the grid point 153 154 158 155 159 156 DO i = nxl, nxr 160 157 DO j = nys, nyn 161 DO k = nzb_w_outer(j,i)+1, nzt-1 158 DO k = nzb+1, nzt-1 159 ! 160 !-- Predetermine flag to mask topography and wall-bounded grid points. 161 flag = MERGE( 1.0_wp, 0.0_wp, & 162 BTEST( wall_flags_0(k,j,i), 3 ) ) 163 mask_east = MERGE( 1.0_wp, 0.0_wp, & 164 BTEST( wall_flags_0(k,j,i+1), 3 ) ) 165 mask_west = MERGE( 1.0_wp, 0.0_wp, & 166 BTEST( wall_flags_0(k,j,i-1), 3 ) ) 167 mask_south = MERGE( 1.0_wp, 0.0_wp, & 168 BTEST( wall_flags_0(k,j-1,i), 3 ) ) 169 mask_north = MERGE( 1.0_wp, 0.0_wp, & 170 BTEST( wall_flags_0(k,j+1,i), 3 ) ) 162 171 ! 163 172 !-- Interpolate eddy diffusivities on staggered gridpoints 164 kmxp = 0.25_wp * 165 ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )166 kmxm = 0.25_wp * 167 ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )168 kmyp = 0.25_wp * 169 ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )170 kmym = 0.25_wp * 171 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )173 kmxp = 0.25_wp * ( km(k,j,i) + km(k,j,i+1) + & 174 km(k+1,j,i) + km(k+1,j,i+1) ) 175 kmxm = 0.25_wp * ( km(k,j,i) + km(k,j,i-1) + & 176 km(k+1,j,i) + km(k+1,j,i-1) ) 177 kmyp = 0.25_wp * ( km(k,j,i) + km(k+1,j,i) + & 178 km(k,j+1,i) + km(k+1,j+1,i) ) 179 kmym = 0.25_wp * ( km(k,j,i) + km(k+1,j,i) + & 180 km(k,j-1,i) + km(k+1,j-1,i) ) 172 181 173 182 tend(k,j,i) = tend(k,j,i) & 174 & + ( kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx & 175 & + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) & 176 & - kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx & 177 & - kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 178 & ) * ddx & 179 & + ( kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy & 180 & + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) & 181 & - kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy & 182 & - kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 183 & ) * ddy & 184 & + 2.0_wp * ( & 185 & km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 186 & * rho_air(k+1) & 187 & - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 188 & * rho_air(k) & 189 & ) * ddzu(k+1) * drho_air_zw(k) 190 ENDDO 191 192 ! 193 !-- Wall functions at all vertical walls, where necessary 194 IF ( wall_w_x(j,i) /= 0.0_wp .OR. wall_w_y(j,i) /= 0.0_wp ) THEN 195 196 DO k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i) 197 ! 198 !-- Interpolate eddy diffusivities on staggered gridpoints 199 kmxp = 0.25_wp * & 200 ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) ) 201 kmxm = 0.25_wp * & 202 ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) ) 203 kmyp = 0.25_wp * & 204 ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) ) 205 kmym = 0.25_wp * & 206 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 207 208 tend(k,j,i) = tend(k,j,i) & 209 + ( fwxp(j,i) * ( & 210 kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx & 211 + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) & 212 ) & 213 - fwxm(j,i) * ( & 214 kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx & 215 + kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 216 ) & 217 + wall_w_x(j,i) * wsus(k,j,i) & 218 ) * ddx & 219 + ( fwyp(j,i) * ( & 220 kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy & 221 + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) & 222 ) & 223 - fwym(j,i) * ( & 224 kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy & 225 + kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 226 ) & 227 + wall_w_y(j,i) * wsvs(k,j,i) & 228 ) * ddy & 229 + 2.0_wp * ( & 230 km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 231 * rho_air(k+1) & 232 - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 233 * rho_air(k) & 234 ) * ddzu(k+1) * drho_air_zw(k) 235 ENDDO 236 ENDIF 183 + ( mask_east * kmxp * ( & 184 ( w(k,j,i+1) - w(k,j,i) ) * ddx & 185 + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) & 186 ) & 187 - mask_west * kmxm * ( & 188 ( w(k,j,i) - w(k,j,i-1) ) * ddx & 189 + ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 190 ) & 191 ) * ddx * flag & 192 + ( mask_north * kmyp * ( & 193 ( w(k,j+1,i) - w(k,j,i) ) * ddy & 194 + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) & 195 ) & 196 - mask_south * kmym * ( & 197 ( w(k,j,i) - w(k,j-1,i) ) * ddy & 198 + ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 199 ) & 200 ) * ddy * flag & 201 + 2.0_wp * ( & 202 km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 203 * rho_air(k+1) & 204 - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 205 * rho_air(k) & 206 ) * ddzu(k+1) * drho_air_zw(k) * flag 207 ENDDO 208 209 ! 210 !-- Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1) 211 !-- surfaces. Note, in the the flat case, loops won't be entered as 212 !-- start_index > end_index. Furtermore, note, no vertical natural surfaces 213 !-- so far. 214 !-- Default-type surfaces 215 DO l = 0, 1 216 surf_s = surf_def_v(l)%start_index(j,i) 217 surf_e = surf_def_v(l)%end_index(j,i) 218 DO m = surf_s, surf_e 219 k = surf_def_v(l)%k(m) 220 tend(k,j,i) = tend(k,j,i) + & 221 surf_def_v(l)%mom_flux_w(m) * ddy 222 ENDDO 223 ENDDO 224 ! 225 !-- Natural-type surfaces 226 DO l = 0, 1 227 surf_s = surf_lsm_v(l)%start_index(j,i) 228 surf_e = surf_lsm_v(l)%end_index(j,i) 229 DO m = surf_s, surf_e 230 k = surf_lsm_v(l)%k(m) 231 tend(k,j,i) = tend(k,j,i) + & 232 surf_lsm_v(l)%mom_flux_w(m) * ddy 233 ENDDO 234 ENDDO 235 ! 236 !-- Urban-type surfaces 237 DO l = 0, 1 238 surf_s = surf_usm_v(l)%start_index(j,i) 239 surf_e = surf_usm_v(l)%end_index(j,i) 240 DO m = surf_s, surf_e 241 k = surf_usm_v(l)%k(m) 242 tend(k,j,i) = tend(k,j,i) + & 243 surf_usm_v(l)%mom_flux_w(m) * ddy 244 ENDDO 245 ENDDO 246 ! 247 !-- Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3) 248 !-- surface. 249 !-- Default-type surfaces 250 DO l = 2, 3 251 surf_s = surf_def_v(l)%start_index(j,i) 252 surf_e = surf_def_v(l)%end_index(j,i) 253 DO m = surf_s, surf_e 254 k = surf_def_v(l)%k(m) 255 tend(k,j,i) = tend(k,j,i) + & 256 surf_def_v(l)%mom_flux_w(m) * ddx 257 ENDDO 258 ENDDO 259 ! 260 !-- Natural-type surfaces 261 DO l = 2, 3 262 surf_s = surf_lsm_v(l)%start_index(j,i) 263 surf_e = surf_lsm_v(l)%end_index(j,i) 264 DO m = surf_s, surf_e 265 k = surf_lsm_v(l)%k(m) 266 tend(k,j,i) = tend(k,j,i) + & 267 surf_lsm_v(l)%mom_flux_w(m) * ddx 268 ENDDO 269 ENDDO 270 ! 271 !-- Urban-type surfaces 272 DO l = 2, 3 273 surf_s = surf_usm_v(l)%start_index(j,i) 274 surf_e = surf_usm_v(l)%end_index(j,i) 275 DO m = surf_s, surf_e 276 k = surf_usm_v(l)%k(m) 277 tend(k,j,i) = tend(k,j,i) + & 278 surf_usm_v(l)%mom_flux_w(m) * ddx 279 ENDDO 280 ENDDO 237 281 238 282 ENDDO … … 256 300 257 301 USE grid_variables, & 258 ONLY : ddx, ddy , fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y302 ONLY : ddx, ddy 259 303 260 304 USE indices, & 261 ONLY : nxl, nxr, nyn, nys, nzb, nz b_w_inner, nzb_w_outer, nzt305 ONLY : nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0 262 306 263 307 USE kinds 264 308 309 USE surface_mod, & 310 ONLY : surf_def_v, surf_lsm_v, surf_usm_v 311 265 312 IMPLICIT NONE 266 313 267 INTEGER(iwp) :: i !< 268 INTEGER(iwp) :: j !< 269 INTEGER(iwp) :: k !< 314 315 INTEGER(iwp) :: i !< running index x direction 316 INTEGER(iwp) :: j !< running index y direction 317 INTEGER(iwp) :: k !< running index z direction 318 INTEGER(iwp) :: l !< running index of surface type, south- or north-facing wall 319 INTEGER(iwp) :: m !< running index surface elements 320 INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint 321 INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint 270 322 271 REAL(wp) :: kmxm !< 272 REAL(wp) :: kmxp !< 273 REAL(wp) :: kmym !< 274 REAL(wp) :: kmyp !< 275 276 REAL(wp), DIMENSION(nzb:nzt+1) :: wsus 277 REAL(wp), DIMENSION(nzb:nzt+1) :: wsvs 278 279 280 DO k = nzb_w_outer(j,i)+1, nzt-1 323 REAL(wp) :: flag !< flag to mask topography grid points 324 REAL(wp) :: kmxm !< 325 REAL(wp) :: kmxp !< 326 REAL(wp) :: kmym !< 327 REAL(wp) :: kmyp !< 328 REAL(wp) :: mask_west !< flag to mask vertical wall west of the grid point 329 REAL(wp) :: mask_east !< flag to mask vertical wall east of the grid point 330 REAL(wp) :: mask_south !< flag to mask vertical wall south of the grid point 331 REAL(wp) :: mask_north !< flag to mask vertical wall north of the grid point 332 333 334 DO k = nzb+1, nzt-1 335 ! 336 !-- Predetermine flag to mask topography and wall-bounded grid points. 337 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 3 ) ) 338 mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 3 ) ) 339 mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 3 ) ) 340 mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 3 ) ) 341 mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 3 ) ) 281 342 ! 282 343 !-- Interpolate eddy diffusivities on staggered gridpoints … … 287 348 288 349 tend(k,j,i) = tend(k,j,i) & 289 & + ( kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx & 290 & + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) & 291 & - kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx & 292 & - kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 293 & ) * ddx & 294 & + ( kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy & 295 & + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) & 296 & - kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy & 297 & - kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 298 & ) * ddy & 299 & + 2.0_wp * ( & 300 & km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 301 & * rho_air(k+1) & 302 & - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 303 & * rho_air(k) & 304 & ) * ddzu(k+1) * drho_air_zw(k) 305 ENDDO 306 307 ! 308 !-- Wall functions at all vertical walls, where necessary 309 IF ( wall_w_x(j,i) /= 0.0_wp .OR. wall_w_y(j,i) /= 0.0_wp ) THEN 310 311 ! 312 !-- Calculate the horizontal momentum fluxes w'u' and/or w'v' 313 IF ( wall_w_x(j,i) /= 0.0_wp ) THEN 314 CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i), & 315 wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp ) 316 ELSE 317 wsus = 0.0_wp 318 ENDIF 319 320 IF ( wall_w_y(j,i) /= 0.0_wp ) THEN 321 CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i), & 322 wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp ) 323 ELSE 324 wsvs = 0.0_wp 325 ENDIF 326 327 DO k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i) 328 ! 329 !-- Interpolate eddy diffusivities on staggered gridpoints 330 kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) ) 331 kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) ) 332 kmyp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) ) 333 kmym = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 334 335 tend(k,j,i) = tend(k,j,i) & 336 + ( fwxp(j,i) * ( & 337 kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx & 338 + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) & 339 ) & 340 - fwxm(j,i) * ( & 341 kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx & 342 + kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 343 ) & 344 + wall_w_x(j,i) * wsus(k) & 345 ) * ddx & 346 + ( fwyp(j,i) * ( & 347 kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy & 348 + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) & 349 ) & 350 - fwym(j,i) * ( & 351 kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy & 352 + kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 353 ) & 354 + wall_w_y(j,i) * wsvs(k) & 355 ) * ddy & 356 + 2.0_wp * ( & 357 km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 358 * rho_air(k+1) & 359 - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 360 * rho_air(k) & 361 ) * ddzu(k+1) * drho_air_zw(k) 362 ENDDO 363 ENDIF 350 + ( mask_east * kmxp * ( & 351 ( w(k,j,i+1) - w(k,j,i) ) * ddx & 352 + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) & 353 ) & 354 - mask_west * kmxm * ( & 355 ( w(k,j,i) - w(k,j,i-1) ) * ddx & 356 + ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 357 ) & 358 ) * ddx * flag & 359 + ( mask_north * kmyp * ( & 360 ( w(k,j+1,i) - w(k,j,i) ) * ddy & 361 + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) & 362 ) & 363 - mask_south * kmym * ( & 364 ( w(k,j,i) - w(k,j-1,i) ) * ddy & 365 + ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 366 ) & 367 ) * ddy * flag & 368 + 2.0_wp * ( & 369 km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) & 370 * rho_air(k+1) & 371 - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) & 372 * rho_air(k) & 373 ) * ddzu(k+1) * drho_air_zw(k) * flag 374 ENDDO 375 ! 376 !-- Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1) 377 !-- surfaces. Note, in the the flat case, loops won't be entered as 378 !-- start_index > end_index. Furtermore, note, no vertical natural surfaces 379 !-- so far. 380 !-- Default-type surfaces 381 DO l = 0, 1 382 surf_s = surf_def_v(l)%start_index(j,i) 383 surf_e = surf_def_v(l)%end_index(j,i) 384 DO m = surf_s, surf_e 385 k = surf_def_v(l)%k(m) 386 tend(k,j,i) = tend(k,j,i) + & 387 surf_def_v(l)%mom_flux_w(m) * ddy 388 ENDDO 389 ENDDO 390 ! 391 !-- Natural-type surfaces 392 DO l = 0, 1 393 surf_s = surf_lsm_v(l)%start_index(j,i) 394 surf_e = surf_lsm_v(l)%end_index(j,i) 395 DO m = surf_s, surf_e 396 k = surf_lsm_v(l)%k(m) 397 tend(k,j,i) = tend(k,j,i) + & 398 surf_lsm_v(l)%mom_flux_w(m) * ddy 399 ENDDO 400 ENDDO 401 ! 402 !-- Urban-type surfaces 403 DO l = 0, 1 404 surf_s = surf_usm_v(l)%start_index(j,i) 405 surf_e = surf_usm_v(l)%end_index(j,i) 406 DO m = surf_s, surf_e 407 k = surf_usm_v(l)%k(m) 408 tend(k,j,i) = tend(k,j,i) + & 409 surf_usm_v(l)%mom_flux_w(m) * ddy 410 ENDDO 411 ENDDO 412 ! 413 !-- Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3) 414 !-- surfaces. 415 !-- Default-type surfaces 416 DO l = 2, 3 417 surf_s = surf_def_v(l)%start_index(j,i) 418 surf_e = surf_def_v(l)%end_index(j,i) 419 DO m = surf_s, surf_e 420 k = surf_def_v(l)%k(m) 421 tend(k,j,i) = tend(k,j,i) + & 422 surf_def_v(l)%mom_flux_w(m) * ddx 423 ENDDO 424 ENDDO 425 ! 426 !-- Natural-type surfaces 427 DO l = 2, 3 428 surf_s = surf_lsm_v(l)%start_index(j,i) 429 surf_e = surf_lsm_v(l)%end_index(j,i) 430 DO m = surf_s, surf_e 431 k = surf_lsm_v(l)%k(m) 432 tend(k,j,i) = tend(k,j,i) + & 433 surf_lsm_v(l)%mom_flux_w(m) * ddx 434 ENDDO 435 ENDDO 436 ! 437 !-- Urban-type surfaces 438 DO l = 2, 3 439 surf_s = surf_usm_v(l)%start_index(j,i) 440 surf_e = surf_usm_v(l)%end_index(j,i) 441 DO m = surf_s, surf_e 442 k = surf_usm_v(l)%k(m) 443 tend(k,j,i) = tend(k,j,i) + & 444 surf_usm_v(l)%mom_flux_w(m) * ddx 445 ENDDO 446 ENDDO 447 364 448 365 449 END SUBROUTINE diffusion_w_ij
Note: See TracChangeset
for help on using the changeset viewer.
