- Timestamp:
- Apr 29, 2020 2:19:18 PM (10 months ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/tridia_solver_mod.f90
r4360 r4510 1 1 !> @file tridia_solver_mod.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 !--------------------------------------------------------------------------------------------------! 18 ! 19 19 ! 20 20 ! Current revisions: 21 ! ----------------- -21 ! ----------------- 22 22 ! 23 23 ! … … 25 25 ! ----------------- 26 26 ! $Id$ 27 ! file re-formatted to follow the PALM coding standard 28 ! 29 ! 4360 2020-01-07 11:25:50Z suehring 27 30 ! Added missing OpenMP directives 28 ! 31 ! 29 32 ! 4182 2019-08-22 15:20:23Z scharf 30 33 ! Corrected "Former revisions" section 31 ! 34 ! 32 35 ! 3761 2019-02-25 15:31:42Z raasch 33 36 ! OpenACC modification to prevent compiler warning about unused variable 34 ! 37 ! 35 38 ! 3690 2019-01-22 22:56:42Z knoop 36 39 ! OpenACC port for SPEC … … 38 41 ! 1212 2013-08-15 08:46:27Z raasch 39 42 ! Initial revision. 40 ! Routines have been moved to seperate module from former file poisfft to here. 41 ! The tridiagonal matrix coefficients of array tri are calculated only once at42 ! the beginning, i.e. routine split iscalled within tridia_init.43 ! 44 ! 45 ! Description: 46 ! ------------ 47 !> solves the linear system of equations:43 ! Routines have been moved to seperate module from former file poisfft to here. The tridiagonal 44 ! matrix coefficients of array tri are calculated only once at the beginning, i.e. routine split is 45 ! called within tridia_init. 46 ! 47 ! 48 ! Description: 49 ! ------------ 50 !> Solves the linear system of equations: 48 51 !> 49 !> -(4 pi^2(i^2/(dx^2*nnx^2)+j^2/(dy^2*nny^2))+ 50 !> 1/(dzu(k)*dzw(k))+1/(dzu(k-1)*dzw(k)))*p(i,j,k)+ 52 !> -(4 pi^2(i^2/(dx^2*nnx^2)+j^2/(dy^2*nny^2))+ 1/(dzu(k)*dzw(k))+1/(dzu(k-1)*dzw(k)))*p(i,j,k)+ 51 53 !> 1/(dzu(k)*dzw(k))*p(i,j,k+1)+1/(dzu(k-1)*dzw(k))*p(i,j,k-1)=d(i,j,k) 52 54 !> 53 55 !> by using the Thomas algorithm 54 !------------------------------------------------------------------------------ !56 !--------------------------------------------------------------------------------------------------! 55 57 56 58 #define __acc_fft_device ( defined( _OPENACC ) && ( defined ( __cuda_fft ) ) ) 57 59 58 60 MODULE tridia_solver 59 60 61 USE basic_constants_and_equations_mod, &61 62 63 USE basic_constants_and_equations_mod, & 62 64 ONLY: pi 63 65 64 USE indices, & 65 ONLY: nx, ny, nz 66 USE indices, & 67 ONLY: nx, & 68 ny, & 69 nz 66 70 67 71 USE kinds 68 72 69 USE transpose_indices, & 70 ONLY: nxl_z, nyn_z, nxr_z, nys_z 71 72 IMPLICIT NONE 73 74 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddzuw !< 73 USE transpose_indices, & 74 ONLY: nxl_z, & 75 nyn_z, & 76 nxr_z, & 77 nys_z 78 79 IMPLICIT NONE 80 81 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddzuw !< 75 82 76 83 PRIVATE … … 89 96 90 97 91 !------------------------------------------------------------------------------ !98 !--------------------------------------------------------------------------------------------------! 92 99 ! Description: 93 100 ! ------------ 94 101 !> @todo Missing subroutine description. 95 !------------------------------------------------------------------------------! 96 SUBROUTINE tridia_init 97 98 USE arrays_3d, & 99 ONLY: ddzu_pres, ddzw, rho_air_zw 102 !--------------------------------------------------------------------------------------------------! 103 SUBROUTINE tridia_init 104 105 USE arrays_3d, & 106 ONLY: ddzu_pres, & 107 ddzw, & 108 rho_air_zw 100 109 101 110 #if defined( _OPENACC ) 102 USE arrays_3d, &111 USE arrays_3d, & 103 112 ONLY: tri 104 113 #endif 105 114 106 IMPLICIT NONE 107 108 INTEGER(iwp) :: k !< 109 110 ALLOCATE( ddzuw(0:nz-1,3) ) 111 112 DO k = 0, nz-1 113 ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k) 114 ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1) 115 ddzuw(k,3) = -1.0_wp * & 116 ( ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1) + & 117 ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k) ) 118 ENDDO 119 ! 120 !-- Calculate constant coefficients of the tridiagonal matrix 121 CALL maketri 122 CALL split 123 124 #if __acc_fft_device 125 !$ACC ENTER DATA & 126 !$ACC COPYIN(ddzuw(0:nz-1,1:3)) & 127 !$ACC COPYIN(tri(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,1:2)) 128 #endif 129 130 END SUBROUTINE tridia_init 131 132 133 !------------------------------------------------------------------------------! 134 ! Description: 135 ! ------------ 136 !> Computes the i- and j-dependent component of the matrix 137 !> Provide the constant coefficients of the tridiagonal matrix for solution 138 !> of the Poisson equation in Fourier space. 139 !> The coefficients are computed following the method of 140 !> Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan 141 !> Siano's original version by discretizing the Poisson equation, 142 !> before it is Fourier-transformed. 143 !------------------------------------------------------------------------------! 144 SUBROUTINE maketri 145 146 147 USE arrays_3d, & 148 ONLY: tric, rho_air 149 150 USE control_parameters, & 151 ONLY: ibc_p_b, ibc_p_t 152 153 USE grid_variables, & 154 ONLY: dx, dy 155 156 157 IMPLICIT NONE 158 159 INTEGER(iwp) :: i !< 160 INTEGER(iwp) :: j !< 161 INTEGER(iwp) :: k !< 162 INTEGER(iwp) :: nnxh !< 163 INTEGER(iwp) :: nnyh !< 164 165 REAL(wp) :: ll(nxl_z:nxr_z,nys_z:nyn_z) !< 166 167 168 nnxh = ( nx + 1 ) / 2 169 nnyh = ( ny + 1 ) / 2 170 171 DO j = nys_z, nyn_z 172 DO i = nxl_z, nxr_z 173 IF ( j >= 0 .AND. j <= nnyh ) THEN 174 IF ( i >= 0 .AND. i <= nnxh ) THEN 175 ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / & 176 REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + & 177 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / & 178 REAL( ny+1, KIND=wp ) ) ) / ( dy * dy ) 179 ELSE 180 ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / & 181 REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + & 182 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / & 183 REAL( ny+1, KIND=wp ) ) ) / ( dy * dy ) 184 ENDIF 185 ELSE 186 IF ( i >= 0 .AND. i <= nnxh ) THEN 187 ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / & 188 REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + & 189 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / & 190 REAL( ny+1, KIND=wp ) ) ) / ( dy * dy ) 191 ELSE 192 ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / & 193 REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + & 194 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / & 195 REAL( ny+1, KIND=wp ) ) ) / ( dy * dy ) 196 ENDIF 197 ENDIF 198 ENDDO 199 ENDDO 200 201 DO k = 0, nz-1 202 DO j = nys_z, nyn_z 203 DO i = nxl_z, nxr_z 204 tric(i,j,k) = ddzuw(k,3) - ll(i,j) * rho_air(k+1) 205 ENDDO 206 ENDDO 207 ENDDO 208 209 IF ( ibc_p_b == 1 ) THEN 210 DO j = nys_z, nyn_z 211 DO i = nxl_z, nxr_z 212 tric(i,j,0) = tric(i,j,0) + ddzuw(0,1) 213 ENDDO 214 ENDDO 215 ENDIF 216 IF ( ibc_p_t == 1 ) THEN 217 DO j = nys_z, nyn_z 218 DO i = nxl_z, nxr_z 219 tric(i,j,nz-1) = tric(i,j,nz-1) + ddzuw(nz-1,2) 220 ENDDO 221 ENDDO 222 ENDIF 223 224 END SUBROUTINE maketri 225 226 227 !------------------------------------------------------------------------------! 228 ! Description: 229 ! ------------ 230 !> Substitution (Forward and Backward) (Thomas algorithm) 231 !------------------------------------------------------------------------------! 232 SUBROUTINE tridia_substi( ar ) 233 234 235 USE arrays_3d, & 236 ONLY: tri 237 238 USE control_parameters, & 239 ONLY: ibc_p_b, ibc_p_t 240 241 IMPLICIT NONE 242 243 INTEGER(iwp) :: i !< 244 INTEGER(iwp) :: j !< 245 INTEGER(iwp) :: k !< 246 247 REAL(wp) :: ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !< 248 249 REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) :: ar1 !< 250 #if __acc_fft_device 251 !$ACC DECLARE CREATE(ar1) 252 #endif 253 254 !$OMP PARALLEL PRIVATE(i,j,k) 255 256 ! 257 !-- Forward substitution 258 #if __acc_fft_device 259 !$ACC PARALLEL PRESENT(ar, ar1, tri) PRIVATE(i,j,k) 260 #endif 261 DO k = 0, nz - 1 262 #if __acc_fft_device 263 !$ACC LOOP COLLAPSE(2) 264 #endif 265 !$OMP DO 266 DO j = nys_z, nyn_z 267 DO i = nxl_z, nxr_z 268 269 IF ( k == 0 ) THEN 270 ar1(i,j,k) = ar(i,j,k+1) 271 ELSE 272 ar1(i,j,k) = ar(i,j,k+1) - tri(i,j,k,2) * ar1(i,j,k-1) 273 ENDIF 274 275 ENDDO 276 ENDDO 277 ENDDO 278 #if __acc_fft_device 279 !$ACC END PARALLEL 280 #endif 281 282 ! 283 !-- Backward substitution 284 !-- Note, the 1.0E-20 in the denominator is due to avoid divisions 285 !-- by zero appearing if the pressure bc is set to neumann at the top of 286 !-- the model domain. 287 #if __acc_fft_device 288 !$ACC PARALLEL PRESENT(ar, ar1, ddzuw, tri) PRIVATE(i,j,k) 289 #endif 290 DO k = nz-1, 0, -1 291 #if __acc_fft_device 292 !$ACC LOOP COLLAPSE(2) 293 #endif 294 !$OMP DO 295 DO j = nys_z, nyn_z 296 DO i = nxl_z, nxr_z 297 298 IF ( k == nz-1 ) THEN 299 ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20_wp ) 300 ELSE 301 ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) & 302 / tri(i,j,k,1) 303 ENDIF 304 ENDDO 305 ENDDO 306 ENDDO 307 #if __acc_fft_device 308 !$ACC END PARALLEL 309 #endif 310 311 !$OMP END PARALLEL 312 313 ! 314 !-- Indices i=0, j=0 correspond to horizontally averaged pressure. 315 !-- The respective values of ar should be zero at all k-levels if 316 !-- acceleration of horizontally averaged vertical velocity is zero. 317 IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN 318 IF ( nys_z == 0 .AND. nxl_z == 0 ) THEN 319 #if __acc_fft_device 320 !$ACC PARALLEL LOOP PRESENT(ar) 321 #endif 322 DO k = 1, nz 323 ar(nxl_z,nys_z,k) = 0.0_wp 324 ENDDO 115 IMPLICIT NONE 116 117 INTEGER(iwp) :: k !< 118 119 ALLOCATE( ddzuw(0:nz-1,3) ) 120 121 DO k = 0, nz-1 122 ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k) 123 ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1) 124 ddzuw(k,3) = -1.0_wp * ( ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1) + & 125 ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k) ) 126 ENDDO 127 ! 128 !-- Calculate constant coefficients of the tridiagonal matrix 129 CALL maketri 130 CALL split 131 132 #if __acc_fft_device 133 !$ACC ENTER DATA & 134 !$ACC COPYIN(ddzuw(0:nz-1,1:3)) & 135 !$ACC COPYIN(tri(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,1:2)) 136 #endif 137 138 END SUBROUTINE tridia_init 139 140 141 !--------------------------------------------------------------------------------------------------! 142 ! Description: 143 ! ------------ 144 !> Computes the i- and j-dependent component of the matrix. 145 !> Provide the constant coefficients of the tridiagonal matrix for solution of the Poisson equation 146 !> in Fourier space. The coefficients are computed following the method of Schmidt et al. 147 !> (DFVLR-Mitteilung 84-15), which departs from Stephan Siano's original version by discretizing the 148 !> Poisson equation, before it is Fourier-transformed. 149 !--------------------------------------------------------------------------------------------------! 150 SUBROUTINE maketri 151 152 153 USE arrays_3d, & 154 ONLY: tric, & 155 rho_air 156 157 USE control_parameters, & 158 ONLY: ibc_p_b, & 159 ibc_p_t 160 161 USE grid_variables, & 162 ONLY: dx, & 163 dy 164 165 166 IMPLICIT NONE 167 168 INTEGER(iwp) :: i !< 169 INTEGER(iwp) :: j !< 170 INTEGER(iwp) :: k !< 171 INTEGER(iwp) :: nnxh !< 172 INTEGER(iwp) :: nnyh !< 173 174 REAL(wp) :: ll(nxl_z:nxr_z,nys_z:nyn_z) !< 175 176 177 nnxh = ( nx + 1 ) / 2 178 nnyh = ( ny + 1 ) / 2 179 180 DO j = nys_z, nyn_z 181 DO i = nxl_z, nxr_z 182 IF ( j >= 0 .AND. j <= nnyh ) THEN 183 IF ( i >= 0 .AND. i <= nnxh ) THEN 184 ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / & 185 REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + & 186 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / & 187 REAL( ny+1, KIND=wp ) ) ) / ( dy * dy ) 188 ELSE 189 ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / & 190 REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + & 191 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / & 192 REAL( ny+1, KIND=wp ) ) ) / ( dy * dy ) 193 ENDIF 194 ELSE 195 IF ( i >= 0 .AND. i <= nnxh ) THEN 196 ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / & 197 REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + & 198 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / & 199 REAL( ny+1, KIND=wp ) ) ) / ( dy * dy ) 200 ELSE 201 ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / & 202 REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + & 203 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / & 204 REAL( ny+1, KIND=wp ) ) ) / ( dy * dy ) 325 205 ENDIF 326 206 ENDIF 327 328 END SUBROUTINE tridia_substi 329 330 331 !------------------------------------------------------------------------------! 332 ! Description: 333 ! ------------ 334 !> Substitution (Forward and Backward) (Thomas algorithm) 335 !------------------------------------------------------------------------------! 336 SUBROUTINE tridia_substi_overlap( ar, jj ) 337 338 339 USE arrays_3d, & 340 ONLY: tri 341 342 USE control_parameters, & 343 ONLY: ibc_p_b, ibc_p_t 344 345 IMPLICIT NONE 346 347 INTEGER(iwp) :: i !< 348 INTEGER(iwp) :: j !< 349 INTEGER(iwp) :: jj !< 350 INTEGER(iwp) :: k !< 351 352 REAL(wp) :: ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !< 353 354 REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) :: ar1 !< 355 356 ! 357 !-- Forward substitution 358 DO k = 0, nz - 1 359 DO j = nys_z, nyn_z 360 DO i = nxl_z, nxr_z 361 362 IF ( k == 0 ) THEN 363 ar1(i,j,k) = ar(i,j,k+1) 364 ELSE 365 ar1(i,j,k) = ar(i,j,k+1) - tri(i,jj,k,2) * ar1(i,j,k-1) 366 ENDIF 367 368 ENDDO 369 ENDDO 370 ENDDO 371 372 ! 373 !-- Backward substitution 374 !-- Note, the 1.0E-20 in the denominator is due to avoid divisions 375 !-- by zero appearing if the pressure bc is set to neumann at the top of 376 !-- the model domain. 377 DO k = nz-1, 0, -1 378 DO j = nys_z, nyn_z 379 DO i = nxl_z, nxr_z 380 381 IF ( k == nz-1 ) THEN 382 ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20_wp ) 383 ELSE 384 ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) & 385 / tri(i,jj,k,1) 386 ENDIF 387 ENDDO 388 ENDDO 389 ENDDO 390 391 ! 392 !-- Indices i=0, j=0 correspond to horizontally averaged pressure. 393 !-- The respective values of ar should be zero at all k-levels if 394 !-- acceleration of horizontally averaged vertical velocity is zero. 395 IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN 396 IF ( nys_z == 0 .AND. nxl_z == 0 ) THEN 397 DO k = 1, nz 398 ar(nxl_z,nys_z,k) = 0.0_wp 399 ENDDO 207 ENDDO 208 ENDDO 209 210 DO k = 0, nz-1 211 DO j = nys_z, nyn_z 212 DO i = nxl_z, nxr_z 213 tric(i,j,k) = ddzuw(k,3) - ll(i,j) * rho_air(k+1) 214 ENDDO 215 ENDDO 216 ENDDO 217 218 IF ( ibc_p_b == 1 ) THEN 219 DO j = nys_z, nyn_z 220 DO i = nxl_z, nxr_z 221 tric(i,j,0) = tric(i,j,0) + ddzuw(0,1) 222 ENDDO 223 ENDDO 224 ENDIF 225 IF ( ibc_p_t == 1 ) THEN 226 DO j = nys_z, nyn_z 227 DO i = nxl_z, nxr_z 228 tric(i,j,nz-1) = tric(i,j,nz-1) + ddzuw(nz-1,2) 229 ENDDO 230 ENDDO 231 ENDIF 232 233 END SUBROUTINE maketri 234 235 236 !--------------------------------------------------------------------------------------------------! 237 ! Description: 238 ! ------------ 239 !> Substitution (Forward and Backward) (Thomas algorithm). 240 !--------------------------------------------------------------------------------------------------! 241 SUBROUTINE tridia_substi( ar ) 242 243 244 USE arrays_3d, & 245 ONLY: tri 246 247 USE control_parameters, & 248 ONLY: ibc_p_b, & 249 ibc_p_t 250 251 IMPLICIT NONE 252 253 INTEGER(iwp) :: i !< 254 INTEGER(iwp) :: j !< 255 INTEGER(iwp) :: k !< 256 257 REAL(wp) :: ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !< 258 259 REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) :: ar1 !< 260 #if __acc_fft_device 261 !$ACC DECLARE CREATE(ar1) 262 #endif 263 264 !$OMP PARALLEL PRIVATE(i,j,k) 265 266 ! 267 !-- Forward substitution 268 #if __acc_fft_device 269 !$ACC PARALLEL PRESENT(ar, ar1, tri) PRIVATE(i,j,k) 270 #endif 271 DO k = 0, nz - 1 272 #if __acc_fft_device 273 !$ACC LOOP COLLAPSE(2) 274 #endif 275 !$OMP DO 276 DO j = nys_z, nyn_z 277 DO i = nxl_z, nxr_z 278 279 IF ( k == 0 ) THEN 280 ar1(i,j,k) = ar(i,j,k+1) 281 ELSE 282 ar1(i,j,k) = ar(i,j,k+1) - tri(i,j,k,2) * ar1(i,j,k-1) 400 283 ENDIF 401 ENDIF 402 403 END SUBROUTINE tridia_substi_overlap 404 405 406 !------------------------------------------------------------------------------! 407 ! Description: 408 ! ------------ 409 !> Splitting of the tridiagonal matrix (Thomas algorithm) 410 !------------------------------------------------------------------------------! 411 SUBROUTINE split 412 413 414 USE arrays_3d, & 415 ONLY: tri, tric 416 417 IMPLICIT NONE 418 419 INTEGER(iwp) :: i !< 420 INTEGER(iwp) :: j !< 421 INTEGER(iwp) :: k !< 422 ! 423 !-- Splitting 424 DO j = nys_z, nyn_z 425 DO i = nxl_z, nxr_z 426 tri(i,j,0,1) = tric(i,j,0) 427 ENDDO 428 ENDDO 429 430 DO k = 1, nz-1 431 DO j = nys_z, nyn_z 432 DO i = nxl_z, nxr_z 433 tri(i,j,k,2) = ddzuw(k,1) / tri(i,j,k-1,1) 434 tri(i,j,k,1) = tric(i,j,k) - ddzuw(k-1,2) * tri(i,j,k,2) 435 ENDDO 436 ENDDO 437 ENDDO 438 439 END SUBROUTINE split 440 441 442 !------------------------------------------------------------------------------! 443 ! Description: 444 ! ------------ 445 !> Solves the linear system of equations for a 1d-decomposition along x (see 446 !> tridia) 284 285 ENDDO 286 ENDDO 287 ENDDO 288 #if __acc_fft_device 289 !$ACC END PARALLEL 290 #endif 291 292 ! 293 !-- Backward substitution 294 !-- Note, the 1.0E-20 in the denominator is due to avoid divisions by zero appearing if the 295 !-- pressure bc is set to neumann at the top of the model domain. 296 #if __acc_fft_device 297 !$ACC PARALLEL PRESENT(ar, ar1, ddzuw, tri) PRIVATE(i,j,k) 298 #endif 299 DO k = nz-1, 0, -1 300 #if __acc_fft_device 301 !$ACC LOOP COLLAPSE(2) 302 #endif 303 !$OMP DO 304 DO j = nys_z, nyn_z 305 DO i = nxl_z, nxr_z 306 307 IF ( k == nz-1 ) THEN 308 ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20_wp ) 309 ELSE 310 ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) / tri(i,j,k,1) 311 ENDIF 312 ENDDO 313 ENDDO 314 ENDDO 315 #if __acc_fft_device 316 !$ACC END PARALLEL 317 #endif 318 319 !$OMP END PARALLEL 320 321 ! 322 !-- Indices i=0, j=0 correspond to horizontally averaged pressure. The respective values of ar 323 !-- should be zero at all k-levels if acceleration of horizontally averaged vertical velocity 324 !-- is zero. 325 IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN 326 IF ( nys_z == 0 .AND. nxl_z == 0 ) THEN 327 #if __acc_fft_device 328 !$ACC PARALLEL LOOP PRESENT(ar) 329 #endif 330 DO k = 1, nz 331 ar(nxl_z,nys_z,k) = 0.0_wp 332 ENDDO 333 ENDIF 334 ENDIF 335 336 END SUBROUTINE tridia_substi 337 338 339 !--------------------------------------------------------------------------------------------------! 340 ! Description: 341 ! ------------ 342 !> Substitution (Forward and Backward) (Thomas algorithm). 343 !--------------------------------------------------------------------------------------------------! 344 SUBROUTINE tridia_substi_overlap( ar, jj ) 345 346 347 USE arrays_3d, & 348 ONLY: tri 349 350 USE control_parameters, & 351 ONLY: ibc_p_b, & 352 ibc_p_t 353 354 IMPLICIT NONE 355 356 INTEGER(iwp) :: i !< 357 INTEGER(iwp) :: j !< 358 INTEGER(iwp) :: jj !< 359 INTEGER(iwp) :: k !< 360 361 REAL(wp) :: ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !< 362 363 REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) :: ar1 !< 364 365 ! 366 !-- Forward substitution 367 DO k = 0, nz - 1 368 DO j = nys_z, nyn_z 369 DO i = nxl_z, nxr_z 370 371 IF ( k == 0 ) THEN 372 ar1(i,j,k) = ar(i,j,k+1) 373 ELSE 374 ar1(i,j,k) = ar(i,j,k+1) - tri(i,jj,k,2) * ar1(i,j,k-1) 375 ENDIF 376 377 ENDDO 378 ENDDO 379 ENDDO 380 381 ! 382 !-- Backward substitution 383 !-- Note, the 1.0E-20 in the denominator is due to avoid divisions by zero appearing if the 384 !-- pressure bc is set to neumann at the top of the model domain. 385 DO k = nz-1, 0, -1 386 DO j = nys_z, nyn_z 387 DO i = nxl_z, nxr_z 388 389 IF ( k == nz-1 ) THEN 390 ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20_wp ) 391 ELSE 392 ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) / tri(i,jj,k,1) 393 ENDIF 394 ENDDO 395 ENDDO 396 ENDDO 397 398 ! 399 !-- Indices i=0, j=0 correspond to horizontally averaged pressure. The respective values of ar 400 !-- should be zero at all k-levels if acceleration of horizontally averaged vertical velocity 401 !-- is zero. 402 IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN 403 IF ( nys_z == 0 .AND. nxl_z == 0 ) THEN 404 DO k = 1, nz 405 ar(nxl_z,nys_z,k) = 0.0_wp 406 ENDDO 407 ENDIF 408 ENDIF 409 410 END SUBROUTINE tridia_substi_overlap 411 412 413 !--------------------------------------------------------------------------------------------------! 414 ! Description: 415 ! ------------ 416 !> Splitting of the tridiagonal matrix (Thomas algorithm). 417 !--------------------------------------------------------------------------------------------------! 418 SUBROUTINE split 419 420 421 USE arrays_3d, & 422 ONLY: tri, & 423 tric 424 425 IMPLICIT NONE 426 427 INTEGER(iwp) :: i !< 428 INTEGER(iwp) :: j !< 429 INTEGER(iwp) :: k !< 430 ! 431 ! Splitting 432 DO j = nys_z, nyn_z 433 DO i = nxl_z, nxr_z 434 tri(i,j,0,1) = tric(i,j,0) 435 ENDDO 436 ENDDO 437 438 DO k = 1, nz-1 439 DO j = nys_z, nyn_z 440 DO i = nxl_z, nxr_z 441 tri(i,j,k,2) = ddzuw(k,1) / tri(i,j,k-1,1) 442 tri(i,j,k,1) = tric(i,j,k) - ddzuw(k-1,2) * tri(i,j,k,2) 443 ENDDO 444 ENDDO 445 ENDDO 446 447 END SUBROUTINE split 448 449 450 !--------------------------------------------------------------------------------------------------! 451 ! Description: 452 ! ------------ 453 !> Solves the linear system of equations for a 1d-decomposition along x (see tridia). 447 454 !> 448 !> @attention when using the intel compilers older than 12.0, array tri must 449 !> be passed as an argument to the contained subroutines. Otherwise 450 !> addres faults will occur. This feature can be activated with 451 !> cpp-switch __intel11 452 !> On NEC, tri should not be passed (except for routine substi_1dd) 453 !> because this causes very bad performance. 454 !------------------------------------------------------------------------------! 455 456 SUBROUTINE tridia_1dd( ddx2, ddy2, nx, ny, j, ar, tri_for_1d ) 457 458 459 USE arrays_3d, & 460 ONLY: ddzu_pres, ddzw, rho_air, rho_air_zw 461 462 USE control_parameters, & 463 ONLY: ibc_p_b, ibc_p_t 464 465 IMPLICIT NONE 466 467 INTEGER(iwp) :: i !< 468 INTEGER(iwp) :: j !< 469 INTEGER(iwp) :: k !< 470 INTEGER(iwp) :: nnyh !< 471 INTEGER(iwp) :: nx !< 472 INTEGER(iwp) :: ny !< 473 474 REAL(wp) :: ddx2 !< 475 REAL(wp) :: ddy2 !< 476 477 REAL(wp), DIMENSION(0:nx,1:nz) :: ar !< 478 REAL(wp), DIMENSION(5,0:nx,0:nz-1) :: tri_for_1d !< 479 480 481 nnyh = ( ny + 1 ) / 2 482 483 ! 484 !-- Define constant elements of the tridiagonal matrix. 485 !-- The compiler on SX6 does loop exchange. If 0:nx is a high power of 2, 486 !-- the exchanged loops create bank conflicts. The following directive 487 !-- prohibits loop exchange and the loops perform much better. 455 !> @attention When using intel compilers older than 12.0, array tri must be passed as an argument to 456 !> the contained subroutines. Otherwise address faults will occur. This feature can be 457 !> activated with cpp-switch __intel11. On NEC, tri should not be passed 458 !> (except for routine substi_1dd) because this causes very bad performance. 459 !--------------------------------------------------------------------------------------------------! 460 461 SUBROUTINE tridia_1dd( ddx2, ddy2, nx, ny, j, ar, tri_for_1d ) 462 463 464 USE arrays_3d, & 465 ONLY: ddzu_pres, & 466 ddzw, & 467 rho_air, & 468 rho_air_zw 469 470 USE control_parameters, & 471 ONLY: ibc_p_b, & 472 ibc_p_t 473 474 IMPLICIT NONE 475 476 INTEGER(iwp) :: i !< 477 INTEGER(iwp) :: j !< 478 INTEGER(iwp) :: k !< 479 INTEGER(iwp) :: nnyh !< 480 INTEGER(iwp) :: nx !< 481 INTEGER(iwp) :: ny !< 482 483 REAL(wp) :: ddx2 !< 484 REAL(wp) :: ddy2 !< 485 486 REAL(wp), DIMENSION(0:nx,1:nz) :: ar !< 487 REAL(wp), DIMENSION(5,0:nx,0:nz-1) :: tri_for_1d !< 488 489 490 nnyh = ( ny + 1 ) / 2 491 492 ! 493 !-- Define constant elements of the tridiagonal matrix. The compiler on SX6 does loop exchange. 494 !-- If 0:nx is a high power of 2, the exchanged loops create bank conflicts. The following directive 495 !-- prohibits loop exchange and the loops perform much better. 488 496 !CDIR NOLOOPCHG 489 DO k = 0, nz-1 490 DO i = 0,nx 491 tri_for_1d(2,i,k) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k) 492 tri_for_1d(3,i,k) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1) 493 ENDDO 494 ENDDO 495 496 IF ( j <= nnyh ) THEN 497 CALL maketri_1dd( j ) 497 DO k = 0, nz-1 498 DO i = 0,nx 499 tri_for_1d(2,i,k) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k) 500 tri_for_1d(3,i,k) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1) 501 ENDDO 502 ENDDO 503 504 IF ( j <= nnyh ) THEN 505 CALL maketri_1dd( j ) 506 ELSE 507 CALL maketri_1dd( ny+1-j ) 508 ENDIF 509 510 CALL split_1dd 511 CALL substi_1dd( ar, tri_for_1d ) 512 513 CONTAINS 514 515 516 !--------------------------------------------------------------------------------------------------! 517 ! Description: 518 ! ------------ 519 !> Computes the i- and j-dependent component of the matrix. 520 !--------------------------------------------------------------------------------------------------! 521 SUBROUTINE maketri_1dd( j ) 522 523 IMPLICIT NONE 524 525 INTEGER(iwp) :: i !< 526 INTEGER(iwp) :: j !< 527 INTEGER(iwp) :: k !< 528 INTEGER(iwp) :: nnxh !< 529 530 REAL(wp) :: a !< 531 REAL(wp) :: c !< 532 533 REAL(wp), DIMENSION(0:nx) :: l !< 534 535 536 nnxh = ( nx + 1 ) / 2 537 ! 538 !-- Provide the tridiagonal matrix for solution of the Poisson equation in Fourier space. 539 !-- The coefficients are computed following the method of Schmidt et al. (DFVLR-Mitteilung 84-15), 540 !-- which departs from Stephan Siano's original version by discretizing the Poisson equation, 541 !-- before it is Fourier-transformed. 542 DO i = 0, nx 543 IF ( i >= 0 .AND. i <= nnxh ) THEN 544 l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / & 545 REAL( nx+1, KIND=wp ) ) ) * ddx2 + & 546 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / & 547 REAL( ny+1, KIND=wp ) ) ) * ddy2 498 548 ELSE 499 CALL maketri_1dd( ny+1-j ) 549 l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / & 550 REAL( nx+1, KIND=wp ) ) ) * ddx2 + & 551 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / & 552 REAL( ny+1, KIND=wp ) ) ) * ddy2 500 553 ENDIF 501 502 CALL split_1dd 503 CALL substi_1dd( ar, tri_for_1d ) 504 505 CONTAINS 506 507 508 !------------------------------------------------------------------------------! 509 ! Description: 510 ! ------------ 511 !> computes the i- and j-dependent component of the matrix 512 !------------------------------------------------------------------------------! 513 SUBROUTINE maketri_1dd( j ) 514 515 IMPLICIT NONE 516 517 INTEGER(iwp) :: i !< 518 INTEGER(iwp) :: j !< 519 INTEGER(iwp) :: k !< 520 INTEGER(iwp) :: nnxh !< 521 522 REAL(wp) :: a !< 523 REAL(wp) :: c !< 524 525 REAL(wp), DIMENSION(0:nx) :: l !< 526 527 528 nnxh = ( nx + 1 ) / 2 529 ! 530 !-- Provide the tridiagonal matrix for solution of the Poisson equation in 531 !-- Fourier space. The coefficients are computed following the method of 532 !-- Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan 533 !-- Siano's original version by discretizing the Poisson equation, 534 !-- before it is Fourier-transformed 535 DO i = 0, nx 536 IF ( i >= 0 .AND. i <= nnxh ) THEN 537 l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / & 538 REAL( nx+1, KIND=wp ) ) ) * ddx2 + & 539 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / & 540 REAL( ny+1, KIND=wp ) ) ) * ddy2 541 ELSE 542 l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / & 543 REAL( nx+1, KIND=wp ) ) ) * ddx2 + & 544 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / & 545 REAL( ny+1, KIND=wp ) ) ) * ddy2 546 ENDIF 547 ENDDO 548 549 DO k = 0, nz-1 550 DO i = 0, nx 551 a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1) 552 c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k) 553 tri_for_1d(1,i,k) = a + c - l(i) * rho_air(k+1) 554 ENDDO 555 ENDDO 556 IF ( ibc_p_b == 1 ) THEN 557 DO i = 0, nx 558 tri_for_1d(1,i,0) = tri_for_1d(1,i,0) + tri_for_1d(2,i,0) 559 ENDDO 560 ENDIF 561 IF ( ibc_p_t == 1 ) THEN 562 DO i = 0, nx 563 tri_for_1d(1,i,nz-1) = tri_for_1d(1,i,nz-1) + tri_for_1d(3,i,nz-1) 564 ENDDO 565 ENDIF 566 567 END SUBROUTINE maketri_1dd 568 569 570 !------------------------------------------------------------------------------! 571 ! Description: 572 ! ------------ 573 !> Splitting of the tridiagonal matrix (Thomas algorithm) 574 !------------------------------------------------------------------------------! 575 SUBROUTINE split_1dd 576 577 IMPLICIT NONE 578 579 INTEGER(iwp) :: i !< 580 INTEGER(iwp) :: k !< 581 582 583 ! 584 !-- Splitting 585 DO i = 0, nx 586 tri_for_1d(4,i,0) = tri_for_1d(1,i,0) 587 ENDDO 588 DO k = 1, nz-1 589 DO i = 0, nx 590 tri_for_1d(5,i,k) = tri_for_1d(2,i,k) / tri_for_1d(4,i,k-1) 591 tri_for_1d(4,i,k) = tri_for_1d(1,i,k) - tri_for_1d(3,i,k-1) * tri_for_1d(5,i,k) 592 ENDDO 593 ENDDO 594 595 END SUBROUTINE split_1dd 596 597 598 !------------------------------------------------------------------------------! 599 ! Description: 600 ! ------------ 601 !> Substitution (Forward and Backward) (Thomas algorithm) 602 !------------------------------------------------------------------------------! 603 SUBROUTINE substi_1dd( ar, tri_for_1d ) 604 605 606 IMPLICIT NONE 607 608 INTEGER(iwp) :: i !< 609 INTEGER(iwp) :: k !< 610 611 REAL(wp), DIMENSION(0:nx,nz) :: ar !< 612 REAL(wp), DIMENSION(0:nx,0:nz-1) :: ar1 !< 613 REAL(wp), DIMENSION(5,0:nx,0:nz-1) :: tri_for_1d !< 614 615 ! 616 !-- Forward substitution 617 DO i = 0, nx 618 ar1(i,0) = ar(i,1) 619 ENDDO 620 DO k = 1, nz-1 621 DO i = 0, nx 622 ar1(i,k) = ar(i,k+1) - tri_for_1d(5,i,k) * ar1(i,k-1) 623 ENDDO 624 ENDDO 625 626 ! 627 !-- Backward substitution 628 !-- Note, the add of 1.0E-20 in the denominator is due to avoid divisions 629 !-- by zero appearing if the pressure bc is set to neumann at the top of 630 !-- the model domain. 631 DO i = 0, nx 632 ar(i,nz) = ar1(i,nz-1) / ( tri_for_1d(4,i,nz-1) + 1.0E-20_wp ) 633 ENDDO 634 DO k = nz-2, 0, -1 635 DO i = 0, nx 636 ar(i,k+1) = ( ar1(i,k) - tri_for_1d(3,i,k) * ar(i,k+2) ) & 637 / tri_for_1d(4,i,k) 638 ENDDO 639 ENDDO 640 641 ! 642 !-- Indices i=0, j=0 correspond to horizontally averaged pressure. 643 !-- The respective values of ar should be zero at all k-levels if 644 !-- acceleration of horizontally averaged vertical velocity is zero. 645 IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN 646 IF ( j == 0 ) THEN 647 DO k = 1, nz 648 ar(0,k) = 0.0_wp 649 ENDDO 650 ENDIF 651 ENDIF 652 653 END SUBROUTINE substi_1dd 654 655 END SUBROUTINE tridia_1dd 554 ENDDO 555 556 DO k = 0, nz-1 557 DO i = 0, nx 558 a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1) 559 c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k) 560 tri_for_1d(1,i,k) = a + c - l(i) * rho_air(k+1) 561 ENDDO 562 ENDDO 563 IF ( ibc_p_b == 1 ) THEN 564 DO i = 0, nx 565 tri_for_1d(1,i,0) = tri_for_1d(1,i,0) + tri_for_1d(2,i,0) 566 ENDDO 567 ENDIF 568 IF ( ibc_p_t == 1 ) THEN 569 DO i = 0, nx 570 tri_for_1d(1,i,nz-1) = tri_for_1d(1,i,nz-1) + tri_for_1d(3,i,nz-1) 571 ENDDO 572 ENDIF 573 574 END SUBROUTINE maketri_1dd 575 576 577 !--------------------------------------------------------------------------------------------------! 578 ! Description: 579 ! ------------ 580 !> Splitting of the tridiagonal matrix (Thomas algorithm). 581 !--------------------------------------------------------------------------------------------------! 582 SUBROUTINE split_1dd 583 584 IMPLICIT NONE 585 586 INTEGER(iwp) :: i !< 587 INTEGER(iwp) :: k !< 588 589 590 ! 591 !-- Splitting 592 DO i = 0, nx 593 tri_for_1d(4,i,0) = tri_for_1d(1,i,0) 594 ENDDO 595 DO k = 1, nz-1 596 DO i = 0, nx 597 tri_for_1d(5,i,k) = tri_for_1d(2,i,k) / tri_for_1d(4,i,k-1) 598 tri_for_1d(4,i,k) = tri_for_1d(1,i,k) - tri_for_1d(3,i,k-1) * tri_for_1d(5,i,k) 599 ENDDO 600 ENDDO 601 602 END SUBROUTINE split_1dd 603 604 605 !--------------------------------------------------------------------------------------------------! 606 ! Description: 607 ! ------------ 608 !> Substitution (Forward and Backward) (Thomas algorithm). 609 !--------------------------------------------------------------------------------------------------! 610 SUBROUTINE substi_1dd( ar, tri_for_1d ) 611 612 613 IMPLICIT NONE 614 615 INTEGER(iwp) :: i !< 616 INTEGER(iwp) :: k !< 617 618 REAL(wp), DIMENSION(0:nx,nz) :: ar !< 619 REAL(wp), DIMENSION(0:nx,0:nz-1) :: ar1 !< 620 REAL(wp), DIMENSION(5,0:nx,0:nz-1) :: tri_for_1d !< 621 622 ! 623 !-- Forward substitution 624 DO i = 0, nx 625 ar1(i,0) = ar(i,1) 626 ENDDO 627 DO k = 1, nz-1 628 DO i = 0, nx 629 ar1(i,k) = ar(i,k+1) - tri_for_1d(5,i,k) * ar1(i,k-1) 630 ENDDO 631 ENDDO 632 633 ! 634 !-- Backward substitution 635 !-- Note, the add of 1.0E-20 in the denominator is due to avoid divisions by zero appearing if the 636 !-- pressure bc is set to neumann at the top of the model domain. 637 DO i = 0, nx 638 ar(i,nz) = ar1(i,nz-1) / ( tri_for_1d(4,i,nz-1) + 1.0E-20_wp ) 639 ENDDO 640 DO k = nz-2, 0, -1 641 DO i = 0, nx 642 ar(i,k+1) = ( ar1(i,k) - tri_for_1d(3,i,k) * ar(i,k+2) ) / tri_for_1d(4,i,k) 643 ENDDO 644 ENDDO 645 646 ! 647 !-- Indices i=0, j=0 correspond to horizontally averaged pressure. The respective values of ar 648 !-- should be zero at all k-levels if acceleration of horizontally averaged vertical velocity is 649 !-- zero. 650 IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN 651 IF ( j == 0 ) THEN 652 DO k = 1, nz 653 ar(0,k) = 0.0_wp 654 ENDDO 655 ENDIF 656 ENDIF 657 658 END SUBROUTINE substi_1dd 659 660 END SUBROUTINE tridia_1dd 656 661 657 662 -
palm/trunk/SOURCE/turbulence_closure_mod.f90
r4495 r4510 1 1 !> @file turbulence_closure_mod.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 Software7 ! Foundation, either version 3 of the License, or (at your option) any later8 ! version.9 ! 10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR12 ! 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 with15 ! PALM. If not, see <http://www.gnu.org/licenses/>.16 ! 17 ! Copyright 2017-2020 Leibniz Universitaet Hannover18 ! --------------------------------------------------------------------------------!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/>. 15 ! 16 ! Copyright 1997-2020 Leibniz Universitaet Hannover 17 !--------------------------------------------------------------------------------------------------! 18 ! 19 19 ! 20 20 ! Current revisions: … … 25 25 ! ----------------- 26 26 ! $Id$ 27 ! file re-formatted to follow the PALM coding standard 28 ! 29 ! 4495 2020-04-13 20:11:20Z raasch 27 30 ! workaround for Intel14 compiler added 28 ! 31 ! 29 32 ! 4486 2020-04-02 20:45:12Z maronga 30 ! Bugfix: include topography in calculation of distance_to_wall (1.5-order-dai 31 ! closure) 32 ! 33 ! Bugfix: include topography in calculation of distance_to_wall (1.5-order-dai closure) 34 ! 33 35 ! 4481 2020-03-31 18:55:54Z maronga 34 ! - added new LES closure after Dai et al. (2020), which provides much better grid 35 ! convergence in stable boundary layer runs. The implementation is experimental36 ! at the moment and should be used with special care. The new SGS closure can be37 ! switched on viaturbulence_closure = '1.5-order-dai'38 ! - variable ml_wall_adjusted renamed to delta as it represents a grid size and 39 ! not a mixing length(see Equations 14 and 18 in Maronga et al. 2015, GMD)36 ! - added new LES closure after Dai et al. (2020), which provides much better grid convergence in 37 ! stable boundary layer runs. The implementation is experimental at the moment and should be used 38 ! with special care. The new SGS closure can be switched on via 39 ! turbulence_closure = '1.5-order-dai' 40 ! - variable ml_wall_adjusted renamed to delta as it represents a grid size and not a mixing length 41 ! (see Equations 14 and 18 in Maronga et al. 2015, GMD) 40 42 ! - nameing of turbulence closures revised: 41 43 ! 'Moeng_Wyngaard' to '1.5-order' … … 44 46 ! - LOGICAL steering variable renamed: 45 47 ! les_mw to les_default 46 ! 48 ! 47 49 ! 4473 2020-03-25 21:04:07Z gronemeier 48 50 ! - rename l-grid to gridsize-geometric-mean … … 54 56 ! l to ml 55 57 ! - adjust some comments 56 ! - corrected definition of wall-adjusted mixing length to include 57 ! gridsize-geometric-mean 58 ! - corrected definition of wall-adjusted mixing length to include gridsize-geometric-mean 58 59 ! - moved definition of wall_adjustment_factor to this module 59 60 ! … … 68 69 ! 69 70 ! 4346 2019-12-18 11:55:56Z motisi 70 ! Introduction of wall_flags_total_0, which currently sets bits based on static 71 ! topographyinformation used in wall_flags_static_071 ! Introduction of wall_flags_total_0, which currently sets bits based on static topography 72 ! information used in wall_flags_static_0 72 73 ! 73 74 ! 4329 2019-12-10 15:46:36Z motisi … … 81 82 ! 82 83 ! 4170 2019-08-19 17:12:31Z gronemeier 83 ! - add performance optimizations according to K. Ketelsen 84 ! t o diffusion_e and tcm_diffusivities_default84 ! - add performance optimizations according to K. Ketelsen to diffusion_e and 85 ! tcm_diffusivities_default 85 86 ! - bugfix in calculating l_wall for vertical walls 86 87 ! - bugfix in using l_wall in initialization (consider wall_adjustment_factor) … … 91 92 ! 92 93 ! 4110 2019-07-22 17:05:21Z suehring 93 ! pass integer flag array as well as boundary flags to WS scalar advection 94 ! routine 94 ! pass integer flag array as well as boundary flags to WS scalar advection routine 95 95 ! 96 96 ! 4109 2019-07-22 17:00:34Z suehring 97 97 ! - Modularize setting of boundary conditions for TKE and dissipation 98 98 ! - Neumann boundary condition for TKE at model top is set also in child domain 99 ! - Revise setting of Neumann boundary conditions at non-cyclic lateral 100 ! boundaries 101 ! - Bugfix, set Neumann boundary condition for TKE at vertical wall instead of 102 ! an implicit Dirichlet boundary condition which implied a sink of TKE 103 ! at vertical walls 99 ! - Revise setting of Neumann boundary conditions at non-cyclic lateral boundaries 100 ! - Bugfix, set Neumann boundary condition for TKE at vertical wall instead of an implicit Dirichlet 101 ! boundary condition which implied a sink of TKE at vertical walls 104 102 ! 105 103 ! 4048 2019-06-21 21:00:21Z knoop … … 113 111 ! 114 112 ! 3719 2019-02-06 13:10:18Z kanani 115 ! Changed log_point to log_point_s, otherwise this overlaps with 116 ! 'all progn.equations' cpumeasurement.113 ! Changed log_point to log_point_s, otherwise this overlaps with 'all progn.equations' cpu 114 ! measurement. 117 115 ! 118 116 ! 3684 2019-01-20 20:20:58Z knoop … … 136 134 !> @todo Check for random disturbances 137 135 !> @note <Enter notes on the module> 138 !----------------------------------------------------------------------------- !136 !--------------------------------------------------------------------------------------------------! 139 137 MODULE turbulence_closure_mod 140 138 141 139 142 USE arrays_3d, & 143 ONLY: diss, diss_1, diss_2, diss_3, diss_p, dzu, e, e_1, e_2, e_3, & 144 e_p, kh, km, mean_inflow_profiles, prho, pt, tdiss_m, & 145 te_m, tend, u, v, vpt, w 146 147 USE basic_constants_and_equations_mod, & 148 ONLY: g, kappa, lv_d_cp, lv_d_rd, rd_d_rv 149 150 USE control_parameters, & 151 ONLY: bc_dirichlet_l, & 152 bc_dirichlet_n, & 153 bc_dirichlet_r, & 154 bc_dirichlet_s, & 155 bc_radiation_l, & 156 bc_radiation_n, & 157 bc_radiation_r, & 158 bc_radiation_s, & 159 child_domain, & 160 constant_diffusion, dt_3d, e_init, humidity, & 161 initializing_actions, intermediate_timestep_count, & 162 intermediate_timestep_count_max, km_constant, & 163 les_dai, les_dynamic, les_default, & 164 ocean_mode, plant_canopy, prandtl_number, & 165 pt_reference, rans_mode, rans_tke_e, rans_tke_l, & 166 timestep_scheme, turbulence_closure, & 167 turbulent_inflow, use_upstream_for_tke, vpt_reference, & 168 ws_scheme_sca, current_timestep_number 169 170 USE advec_ws, & 140 USE arrays_3d, & 141 ONLY: diss, & 142 diss_1, & 143 diss_2, & 144 diss_3, & 145 diss_p, & 146 dzu, & 147 e, & 148 e_1, & 149 e_2, & 150 e_3, & 151 e_p, & 152 kh, & 153 km, & 154 mean_inflow_profiles, & 155 prho, & 156 pt, & 157 tdiss_m, & 158 te_m, & 159 tend, & 160 u, & 161 v, & 162 vpt, & 163 w 164 165 USE basic_constants_and_equations_mod, & 166 ONLY: g, & 167 kappa, & 168 lv_d_cp, & 169 lv_d_rd, & 170 rd_d_rv 171 172 USE control_parameters, & 173 ONLY: bc_dirichlet_l, & 174 bc_dirichlet_n, & 175 bc_dirichlet_r, & 176 bc_dirichlet_s, & 177 bc_radiation_l, & 178 bc_radiation_n, & 179 bc_radiation_r, & 180 bc_radiation_s, & 181 child_domain, & 182 constant_diffusion, & 183 current_timestep_number, & 184 dt_3d, & 185 e_init, & 186 humidity, & 187 initializing_actions, & 188 intermediate_timestep_count, & 189 intermediate_timestep_count_max, & 190 km_constant, & 191 les_dai, & 192 les_dynamic, & 193 les_default, & 194 ocean_mode, & 195 plant_canopy, & 196 prandtl_number, & 197 pt_reference, & 198 rans_mode, & 199 rans_tke_e, & 200 rans_tke_l, & 201 timestep_scheme, & 202 turbulence_closure, & 203 turbulent_inflow, & 204 use_upstream_for_tke, & 205 vpt_reference, & 206 ws_scheme_sca 207 208 209 USE advec_ws, & 171 210 ONLY: advec_s_ws 172 211 173 USE advec_s_bc_mod, &212 USE advec_s_bc_mod, & 174 213 ONLY: advec_s_bc 175 214 176 USE advec_s_pw_mod, &215 USE advec_s_pw_mod, & 177 216 ONLY: advec_s_pw 178 217 179 USE advec_s_up_mod, &218 USE advec_s_up_mod, & 180 219 ONLY: advec_s_up 181 220 182 USE cpulog, &221 USE cpulog, & 183 222 ONLY: cpu_log, log_point_s 184 223 185 USE indices, & 186 ONLY: advc_flags_s, & 187 nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, & 188 topo_top_ind, & 224 USE indices, & 225 ONLY: advc_flags_s, & 226 nbgp, & 227 nxl, & 228 nxlg, & 229 nxr, & 230 nxrg, & 231 nyn, & 232 nyng, & 233 nys, & 234 nysg, & 235 nzb, & 236 nzt, & 237 topo_top_ind, & 189 238 wall_flags_total_0 190 239 191 240 USE kinds 192 241 193 USE ocean_mod, &242 USE ocean_mod, & 194 243 ONLY: prho_reference 195 244 196 245 USE pegrid 197 246 198 USE plant_canopy_model_mod, &247 USE plant_canopy_model_mod, & 199 248 ONLY: pcm_tendency 200 249 201 USE statistics, & 202 ONLY: hom, hom_sum, statistic_regions 203 204 USE surface_mod, & 205 ONLY: bc_h, & 206 bc_v, & 207 surf_def_h, & 208 surf_def_v, & 209 surf_lsm_h, & 210 surf_lsm_v, & 211 surf_usm_h, & 250 USE statistics, & 251 ONLY: hom, & 252 hom_sum, & 253 statistic_regions 254 255 USE surface_mod, & 256 ONLY: bc_h, & 257 bc_v, & 258 surf_def_h, & 259 surf_def_v, & 260 surf_lsm_h, & 261 surf_lsm_v, & 262 surf_usm_h, & 212 263 surf_usm_v 213 264 … … 225 276 REAL(wp) :: wall_adjustment_factor = 1.8_wp !< adjustment factor for mixing length 226 277 227 REAL(wp), DIMENSION(0:4) :: rans_const_c = & !< model constants for RANS mode (namelist param)228 (/ 0.55_wp, 1.44_wp, 1.92_wp, 1.44_wp, 0.0_wp /) !> default values fit for standard-tke-e closure229 230 REAL(wp), DIMENSION(2) :: rans_const_sigma = & !< model constants for RANS mode, sigma values (namelist param)231 (/ 1.0_wp, 1.30_wp /) !> (sigma_e, sigma_diss)232 233 REAL(wp), ALLOCATABLE, DIMENSION(:) :: ml_blackadar 234 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: delta 235 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: distance_to_wall 278 REAL(wp), DIMENSION(0:4) :: rans_const_c = & !< model constants for RANS mode (namelist param) 279 (/ 0.55_wp, 1.44_wp, 1.92_wp, 1.44_wp, 0.0_wp /) !> default values fit for standard-tke-e closure 280 281 REAL(wp), DIMENSION(2) :: rans_const_sigma = & !< model constants for RANS mode, sigma values (namelist param) 282 (/ 1.0_wp, 1.30_wp /) !> (sigma_e, sigma_diss) 283 284 REAL(wp), ALLOCATABLE, DIMENSION(:) :: ml_blackadar !< mixing length according to Blackadar 285 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: delta !< grid size, possibly limited by wall adjustment factor 286 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: distance_to_wall !< distance to the surface/wall 236 287 237 288 ! 238 289 !-- Public variables 239 PUBLIC c_0, rans_const_c, rans_const_sigma 290 PUBLIC c_0, & 291 rans_const_c, & 292 rans_const_sigma 240 293 241 294 SAVE … … 244 297 ! 245 298 !-- Public subroutines 246 PUBLIC & 247 tcm_boundary_conds, & 248 tcm_check_parameters, & 249 tcm_check_data_output, & 250 tcm_define_netcdf_grid, & 251 tcm_init_arrays, & 252 tcm_init, & 253 tcm_actions, & 254 tcm_prognostic_equations, & 255 tcm_swap_timelevel, & 256 tcm_3d_data_averaging, & 257 tcm_data_output_2d, & 258 tcm_data_output_3d, & 259 tcm_diffusivities 299 PUBLIC tcm_actions, & 300 tcm_boundary_conds, & 301 tcm_check_parameters, & 302 tcm_check_data_output, & 303 tcm_data_output_2d, & 304 tcm_data_output_3d, & 305 tcm_define_netcdf_grid, & 306 tcm_diffusivities, & 307 tcm_init_arrays, & 308 tcm_init, & 309 tcm_prognostic_equations, & 310 tcm_swap_timelevel, & 311 tcm_3d_data_averaging 312 313 314 260 315 261 316 ! … … 342 397 CONTAINS 343 398 344 !------------------------------------------------------------------------------ !399 !--------------------------------------------------------------------------------------------------! 345 400 ! Description: 346 401 ! ------------ 347 402 !> Check parameters routine for turbulence closure module. 348 !------------------------------------------------------------------------------ !403 !--------------------------------------------------------------------------------------------------! 349 404 SUBROUTINE tcm_boundary_conds 350 405 351 USE pmc_interface, &406 USE pmc_interface, & 352 407 ONLY : rans_mode_parent 353 408 … … 364 419 ! 365 420 !-- In LES mode, Neumann conditions with de/x_i=0 are assumed at solid walls. 366 !-- Note, only TKE is prognostic in this case and dissipation is only 367 !-- a diagnostic quantity. 421 !-- Note, only TKE is prognostic in this case and dissipation is only a diagnostic quantity. 368 422 IF ( .NOT. rans_mode ) THEN 369 423 ! … … 384 438 DO l = 0, 3 385 439 ! 386 !-- Note concerning missing ACC directive for this loop: Even though 387 !-- the data structure bc_v is present, it may not contain any 388 !-- allocated arrays in the flat but also in a topography case, 389 !-- leading to a runtime error. Therefore, omit ACC directives 390 !-- for this loop, in contrast to the bc_h loop. 440 !-- Note concerning missing ACC directive for this loop: Even though the data structure 441 !-- bc_v is present, it may not contain any allocated arrays in the flat but also in a 442 !-- topography case, leading to a runtime error. Therefore, omit ACC directives for this 443 !-- loop, in contrast to the bc_h loop. 391 444 !$OMP PARALLEL DO PRIVATE( i, j, k ) 392 445 DO m = 1, bc_v(l)%ns … … 402 455 ! 403 456 !-- Use wall function within constant-flux layer 404 !-- Note, grid points listed in bc_h are not included in any calculations in RANS mode and 405 !-- aretherefore not set here.457 !-- Note, grid points listed in bc_h are not included in any calculations in RANS mode and are 458 !-- therefore not set here. 406 459 ! 407 460 !-- Upward-facing surfaces … … 459 512 ENDIF 460 513 ! 461 !-- Set Neumann boundary condition for TKE at model top. Do this also 462 !-- in case of a nested run. 514 !-- Set Neumann boundary condition for TKE at model top. Do this also in case of a nested run. 463 515 !$ACC KERNELS PRESENT(e_p) 464 516 e_p(nzt+1,:,:) = e_p(nzt,:,:) 465 517 !$ACC END KERNELS 466 518 ! 467 !-- Nesting case: if parent operates in RANS mode and child in LES mode, 468 !-- no TKE is transfered. This case, set Neumann conditions at lateral and 469 !-- top child boundaries. 470 !-- If not ( both either in RANS or in LES mode ), TKE boundary condition 471 !-- is treated in the nesting. 519 !-- Nesting case: if parent operates in RANS mode and child in LES mode, no TKE is transfered. 520 !-- This case, set Neumann conditions at lateral and top child boundaries. 521 !-- If not ( both either in RANS or in LES mode ), TKE boundary condition is treated in the 522 !-- nesting. 472 523 If ( child_domain ) THEN 473 524 IF ( rans_mode_parent .AND. .NOT. rans_mode ) THEN … … 482 533 ENDIF 483 534 ! 484 !-- At in- and outflow boundaries also set Neumann boundary conditions 485 !-- for the SGS-TKE. An exception is made for the child domain if 486 !-- both parent and child operate in RANS mode. This case no 487 !-- lateral Neumann boundary conditions will be set but Dirichlet 488 !-- conditions will be set in the nesting. 489 IF ( .NOT. child_domain .AND. .NOT. rans_mode_parent .AND. & 490 .NOT. rans_mode ) THEN 535 !-- At in- and outflow boundaries also set Neumann boundary conditions for the SGS-TKE. An 536 !-- exception is made for the child domain if both parent and child operate in RANS mode. This 537 !-- case no lateral Neumann boundary conditions will be set but Dirichlet conditions will be set 538 !-- in the nesting. 539 IF ( .NOT. child_domain .AND. .NOT. rans_mode_parent .AND. .NOT. rans_mode ) THEN 491 540 IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN 492 541 e_p(:,nys-1,:) = e_p(:,nys,:) … … 519 568 j = surf_def_h(0)%j(m) 520 569 k = surf_def_h(0)%k(m) 521 diss_p(k,j,i) = surf_def_h(0)%us(m)**3 & 522 / ( kappa * surf_def_h(0)%z_mo(m) ) 570 diss_p(k,j,i) = surf_def_h(0)%us(m)**3 / ( kappa * surf_def_h(0)%z_mo(m) ) 523 571 ENDDO 524 572 ! … … 528 576 j = surf_lsm_h%j(m) 529 577 k = surf_lsm_h%k(m) 530 diss_p(k,j,i) = surf_lsm_h%us(m)**3 & 531 / ( kappa * surf_lsm_h%z_mo(m) ) 578 diss_p(k,j,i) = surf_lsm_h%us(m)**3 / ( kappa * surf_lsm_h%z_mo(m) ) 532 579 ENDDO 533 580 ! … … 537 584 j = surf_usm_h%j(m) 538 585 k = surf_usm_h%k(m) 539 diss_p(k,j,i) = surf_usm_h%us(m)**3 & 540 / ( kappa * surf_usm_h%z_mo(m) ) 586 diss_p(k,j,i) = surf_usm_h%us(m)**3 / ( kappa * surf_usm_h%z_mo(m) ) 541 587 ENDDO 542 588 ! … … 549 595 j = surf_def_v(l)%j(m) 550 596 k = surf_def_v(l)%k(m) 551 diss_p(k,j,i) = surf_def_v(l)%us(m)**3 & 552 / ( kappa * surf_def_v(l)%z_mo(m) ) 597 diss_p(k,j,i) = surf_def_v(l)%us(m)**3 / ( kappa * surf_def_v(l)%z_mo(m) ) 553 598 ENDDO 554 599 ! … … 558 603 j = surf_lsm_v(l)%j(m) 559 604 k = surf_lsm_v(l)%k(m) 560 diss_p(k,j,i) = surf_lsm_v(l)%us(m)**3 & 561 / ( kappa * surf_lsm_v(l)%z_mo(m) ) 605 diss_p(k,j,i) = surf_lsm_v(l)%us(m)**3 / ( kappa * surf_lsm_v(l)%z_mo(m) ) 562 606 ENDDO 563 607 ! … … 567 611 j = surf_usm_v(l)%j(m) 568 612 k = surf_usm_v(l)%k(m) 569 diss_p(k,j,i) = surf_usm_v(l)%us(m)**3 & 570 / ( kappa * surf_usm_v(l)%z_mo(m) ) 613 diss_p(k,j,i) = surf_usm_v(l)%us(m)**3 / ( kappa * surf_usm_v(l)%z_mo(m) ) 571 614 ENDDO 572 615 ENDDO 573 616 ! 574 !-- Limit change of diss to be between -90% and +100%. Also, set an absolute 575 !-- minimum value 617 !-- Limit change of diss to be between -90% and +100%. Also, set an absolute minimum value 576 618 DO i = nxl, nxr 577 619 DO j = nys, nyn 578 620 DO k = nzb, nzt+1 579 diss_p(k,j,i) = MAX( MIN( diss_p(k,j,i), & 580 2.0_wp * diss(k,j,i) ), & 581 0.1_wp * diss(k,j,i), & 582 0.0001_wp ) 621 diss_p(k,j,i) = MAX( MIN( diss_p(k,j,i), 2.0_wp * diss(k,j,i) ), & 622 0.1_wp * diss(k,j,i), 0.0001_wp ) 583 623 ENDDO 584 624 ENDDO … … 591 631 END SUBROUTINE tcm_boundary_conds 592 632 593 !------------------------------------------------------------------------------ !633 !--------------------------------------------------------------------------------------------------! 594 634 ! Description: 595 635 ! ------------ 596 636 !> Check parameters routine for turbulence closure module. 597 !------------------------------------------------------------------------------ !637 !--------------------------------------------------------------------------------------------------! 598 638 SUBROUTINE tcm_check_parameters 599 639 600 USE control_parameters, & 601 ONLY: message_string, turbulent_inflow, turbulent_outflow 640 USE control_parameters, & 641 ONLY: message_string, & 642 turbulent_inflow, & 643 turbulent_outflow 602 644 603 645 IMPLICIT NONE … … 625 667 626 668 CASE DEFAULT 627 message_string = 'Unknown turbulence closure: ' // & 628 TRIM( turbulence_closure ) 669 message_string = 'Unknown turbulence closure: ' // TRIM( turbulence_closure ) 629 670 CALL message( 'tcm_check_parameters', 'PA0500', 1, 2, 0, 6, 0 ) 630 671 … … 641 682 c_1 = rans_const_c(1) 642 683 c_2 = rans_const_c(2) 643 c_3 = rans_const_c(3) !> @todo clarify how to switch between different models684 c_3 = rans_const_c(3) !> @todo Clarify how to switch between different models 644 685 c_4 = rans_const_c(4) 645 686 646 687 IF ( turbulent_inflow .OR. turbulent_outflow ) THEN 647 message_string = 'turbulent inflow/outflow is not yet '// & 648 'implemented for RANS mode' 688 message_string = 'turbulent inflow/outflow is not yet '// 'implemented for RANS mode' 649 689 CALL message( 'tcm_check_parameters', 'PA0501', 1, 2, 0, 6, 0 ) 650 690 ENDIF … … 653 693 ! 654 694 !-- LES mode 655 c_0 = 0.1_wp !according to Lilly (1967) and Deardorff (1980) 656 657 dsig_e = 1.0_wp !assure to use K_m to calculate TKE instead 658 !of K_e which is used in RANS mode 695 c_0 = 0.1_wp !According to Lilly (1967) and Deardorff (1980) 696 697 dsig_e = 1.0_wp !Assure to use K_m to calculate TKE instead of K_e which is used in RANS mode 659 698 660 699 ENDIF … … 662 701 END SUBROUTINE tcm_check_parameters 663 702 664 !------------------------------------------------------------------------------ !703 !--------------------------------------------------------------------------------------------------! 665 704 ! Description: 666 705 ! ------------ 667 706 !> Check data output. 668 !------------------------------------------------------------------------------ !707 !--------------------------------------------------------------------------------------------------! 669 708 SUBROUTINE tcm_check_data_output( var, unit ) 670 709 671 710 IMPLICIT NONE 672 711 673 CHARACTER (LEN=*) :: unit!< unit of output variable674 CHARACTER (LEN=*) :: var!< name of output variable712 CHARACTER(LEN=*) :: unit !< unit of output variable 713 CHARACTER(LEN=*) :: var !< name of output variable 675 714 676 715 … … 691 730 692 731 693 !------------------------------------------------------------------------------ !732 !--------------------------------------------------------------------------------------------------! 694 733 ! Description: 695 734 ! ------------ 696 !> Define appropriate grid for netcdf variables. 697 !> It is called out from subroutine netcdf. 698 !------------------------------------------------------------------------------! 735 !> Define appropriate grid for netcdf variables. It is called out from subroutine netcdf. 736 !--------------------------------------------------------------------------------------------------! 699 737 SUBROUTINE tcm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z ) 700 738 701 739 IMPLICIT NONE 702 740 703 CHARACTER (LEN=*), INTENT(OUT) :: grid_x!< x grid of output variable704 CHARACTER (LEN=*), INTENT(OUT) :: grid_y!< y grid of output variable705 CHARACTER (LEN=*), INTENT(OUT) :: grid_z!< z grid of output variable706 CHARACTER (LEN=*), INTENT(IN) :: var!< name of output variable741 CHARACTER(LEN=*), INTENT(OUT) :: grid_x !< x grid of output variable 742 CHARACTER(LEN=*), INTENT(OUT) :: grid_y !< y grid of output variable 743 CHARACTER(LEN=*), INTENT(OUT) :: grid_z !< z grid of output variable 744 CHARACTER(LEN=*), INTENT(IN) :: var !< name of output variable 707 745 708 746 LOGICAL, INTENT(OUT) :: found !< flag if output variable is found … … 740 778 741 779 742 !------------------------------------------------------------------------------ !780 !--------------------------------------------------------------------------------------------------! 743 781 ! Description: 744 782 ! ------------ 745 783 !> Average 3D data. 746 !------------------------------------------------------------------------------ !784 !--------------------------------------------------------------------------------------------------! 747 785 SUBROUTINE tcm_3d_data_averaging( mode, variable ) 748 786 749 787 750 USE averaging, & 751 ONLY: diss_av, kh_av, km_av 752 753 USE control_parameters, & 788 USE averaging, & 789 ONLY: diss_av, & 790 kh_av, & 791 km_av 792 793 USE control_parameters, & 754 794 ONLY: average_count_3d 755 795 756 796 IMPLICIT NONE 757 797 758 CHARACTER (LEN=*) :: mode!< flag defining mode 'allocate', 'sum' or 'average'759 CHARACTER (LEN=*) :: variable!< name of variable760 761 INTEGER(iwp) :: i 762 INTEGER(iwp) :: j 763 INTEGER(iwp) :: k 798 CHARACTER(LEN=*) :: mode !< flag defining mode 'allocate', 'sum' or 'average' 799 CHARACTER(LEN=*) :: variable !< name of variable 800 801 INTEGER(iwp) :: i !< loop index 802 INTEGER(iwp) :: j !< loop index 803 INTEGER(iwp) :: k !< loop index 764 804 765 805 IF ( mode == 'allocate' ) THEN … … 841 881 DO j = nysg, nyng 842 882 DO k = nzb, nzt+1 843 diss_av(k,j,i) = diss_av(k,j,i) & 844 / REAL( average_count_3d, KIND=wp ) 883 diss_av(k,j,i) = diss_av(k,j,i) / REAL( average_count_3d, KIND = wp ) 845 884 ENDDO 846 885 ENDDO … … 853 892 DO j = nysg, nyng 854 893 DO k = nzb, nzt+1 855 kh_av(k,j,i) = kh_av(k,j,i) & 856 / REAL( average_count_3d, KIND=wp ) 894 kh_av(k,j,i) = kh_av(k,j,i) / REAL( average_count_3d, KIND = wp ) 857 895 ENDDO 858 896 ENDDO … … 865 903 DO j = nysg, nyng 866 904 DO k = nzb, nzt+1 867 km_av(k,j,i) = km_av(k,j,i) & 868 / REAL( average_count_3d, KIND=wp ) 905 km_av(k,j,i) = km_av(k,j,i) / REAL( average_count_3d, KIND = wp ) 869 906 ENDDO 870 907 ENDDO … … 879 916 880 917 881 !------------------------------------------------------------------------------ !918 !--------------------------------------------------------------------------------------------------! 882 919 ! Description: 883 920 ! ------------ 884 921 !> Define 2D output variables. 885 !------------------------------------------------------------------------------! 886 SUBROUTINE tcm_data_output_2d( av, variable, found, grid, mode, local_pf, & 887 nzb_do, nzt_do ) 888 889 USE averaging, & 890 ONLY: diss_av, kh_av, km_av 922 !--------------------------------------------------------------------------------------------------! 923 SUBROUTINE tcm_data_output_2d( av, variable, found, grid, mode, local_pf, nzb_do, nzt_do ) 924 925 USE averaging, & 926 ONLY: diss_av, & 927 kh_av, & 928 km_av 891 929 892 930 IMPLICIT NONE 893 931 894 CHARACTER (LEN=*) :: grid!< name of vertical grid895 CHARACTER (LEN=*) :: mode!< either 'xy', 'xz' or 'yz'896 CHARACTER (LEN=*) :: variable!< name of variable897 898 INTEGER(iwp) :: av 899 INTEGER(iwp) :: flag_nr 900 INTEGER(iwp) :: i 901 INTEGER(iwp) :: j 902 INTEGER(iwp) :: k 903 INTEGER(iwp) :: nzb_do 904 INTEGER(iwp) :: nzt_do 932 CHARACTER(LEN=*) :: grid !< name of vertical grid 933 CHARACTER(LEN=*) :: mode !< either 'xy', 'xz' or 'yz' 934 CHARACTER(LEN=*) :: variable !< name of variable 935 936 INTEGER(iwp) :: av !< flag for (non-)average output 937 INTEGER(iwp) :: flag_nr !< number of masking flag 938 INTEGER(iwp) :: i !< loop index 939 INTEGER(iwp) :: j !< loop index 940 INTEGER(iwp) :: k !< loop index 941 INTEGER(iwp) :: nzb_do !< vertical output index (bottom) 942 INTEGER(iwp) :: nzt_do !< vertical output index (top) 905 943 906 944 LOGICAL :: found !< flag if output variable is found … … 909 947 REAL(wp) :: fill_value = -9999.0_wp !< value for the _FillValue attribute 910 948 911 REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local912 949 REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local 950 !< array to which output data is resorted to 913 951 914 952 REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable … … 926 964 to_be_resorted => diss 927 965 ELSE 928 IF ( .NOT. ALLOCATED( diss_av ) ) THEN966 IF ( .NOT. ALLOCATED( diss_av ) ) THEN 929 967 ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 930 968 diss_av = REAL( fill_value, KIND = wp ) … … 938 976 to_be_resorted => kh 939 977 ELSE 940 IF ( .NOT. ALLOCATED( kh_av ) ) THEN978 IF ( .NOT. ALLOCATED( kh_av ) ) THEN 941 979 ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 942 980 kh_av = REAL( fill_value, KIND = wp ) … … 950 988 to_be_resorted => km 951 989 ELSE 952 IF ( .NOT. ALLOCATED( km_av ) ) THEN990 IF ( .NOT. ALLOCATED( km_av ) ) THEN 953 991 ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 954 992 km_av = REAL( fill_value, KIND = wp ) … … 964 1002 END SELECT 965 1003 966 IF ( found .AND. .NOT.resorted ) THEN1004 IF ( found .AND. .NOT. resorted ) THEN 967 1005 DO i = nxl, nxr 968 1006 DO j = nys, nyn 969 1007 DO k = nzb_do, nzt_do 970 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), & 971 REAL( fill_value, KIND = wp ), & 1008 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, KIND = wp ), & 972 1009 BTEST( wall_flags_total_0(k,j,i), flag_nr ) ) 973 1010 ENDDO … … 979 1016 980 1017 981 !------------------------------------------------------------------------------ !1018 !--------------------------------------------------------------------------------------------------! 982 1019 ! Description: 983 1020 ! ------------ 984 1021 !> Define 3D output variables. 985 !------------------------------------------------------------------------------ !1022 !--------------------------------------------------------------------------------------------------! 986 1023 SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do ) 987 1024 988 1025 989 USE averaging, & 990 ONLY: diss_av, kh_av, km_av 1026 USE averaging, & 1027 ONLY: diss_av, & 1028 kh_av, & 1029 km_av 991 1030 992 1031 IMPLICIT NONE 993 1032 994 CHARACTER (LEN=*) :: variable!< name of variable995 996 INTEGER(iwp) :: av 997 INTEGER(iwp) :: flag_nr 998 INTEGER(iwp) :: i 999 INTEGER(iwp) :: j 1000 INTEGER(iwp) :: k 1001 INTEGER(iwp) :: nzb_do 1002 INTEGER(iwp) :: nzt_do 1033 CHARACTER(LEN=*) :: variable !< name of variable 1034 1035 INTEGER(iwp) :: av !< flag for (non-)average output 1036 INTEGER(iwp) :: flag_nr !< number of masking flag 1037 INTEGER(iwp) :: i !< loop index 1038 INTEGER(iwp) :: j !< loop index 1039 INTEGER(iwp) :: k !< loop index 1040 INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0) 1041 INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d) 1003 1042 1004 1043 LOGICAL :: found !< flag if output variable is found … … 1007 1046 REAL(wp) :: fill_value = -9999.0_wp !< value for the _FillValue attribute 1008 1047 1009 REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf 1048 REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local 1010 1049 !< array to which output data is resorted to 1011 1050 … … 1024 1063 to_be_resorted => diss 1025 1064 ELSE 1026 IF ( .NOT. ALLOCATED( diss_av ) ) THEN1065 IF ( .NOT. ALLOCATED( diss_av ) ) THEN 1027 1066 ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1028 1067 diss_av = REAL( fill_value, KIND = wp ) … … 1035 1074 to_be_resorted => kh 1036 1075 ELSE 1037 IF ( .NOT. ALLOCATED( kh_av ) ) THEN1076 IF ( .NOT. ALLOCATED( kh_av ) ) THEN 1038 1077 ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1039 1078 kh_av = REAL( fill_value, KIND = wp ) … … 1046 1085 to_be_resorted => km 1047 1086 ELSE 1048 IF ( .NOT. ALLOCATED( km_av ) ) THEN1087 IF ( .NOT. ALLOCATED( km_av ) ) THEN 1049 1088 ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1050 1089 km_av = REAL( fill_value, KIND = wp ) … … 1059 1098 1060 1099 1061 IF ( found .AND. .NOT.resorted ) THEN1100 IF ( found .AND. .NOT. resorted ) THEN 1062 1101 DO i = nxl, nxr 1063 1102 DO j = nys, nyn 1064 1103 DO k = nzb_do, nzt_do 1065 local_pf(i,j,k) = MERGE( & 1066 to_be_resorted(k,j,i), & 1067 REAL( fill_value, KIND = wp ), & 1068 BTEST( wall_flags_total_0(k,j,i), flag_nr ) ) 1104 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, KIND = wp ), & 1105 BTEST( wall_flags_total_0(k,j,i), flag_nr ) ) 1069 1106 ENDDO 1070 1107 ENDDO … … 1076 1113 1077 1114 1078 !------------------------------------------------------------------------------ !1115 !--------------------------------------------------------------------------------------------------! 1079 1116 ! Description: 1080 1117 ! ------------ 1081 1118 !> Allocate arrays and assign pointers. 1082 !------------------------------------------------------------------------------ !1119 !--------------------------------------------------------------------------------------------------! 1083 1120 SUBROUTINE tcm_init_arrays 1084 1121 1085 USE bulk_cloud_model_mod, &1122 USE bulk_cloud_model_mod, & 1086 1123 ONLY: collision_turbulence 1087 1124 1088 USE pmc_interface, &1125 USE pmc_interface, & 1089 1126 ONLY: nested_run 1090 1127 … … 1100 1137 ! 1101 1138 !-- Allocate arrays required for dissipation. 1102 !-- Please note, if it is a nested run, arrays need to be allocated even if 1103 !-- they do not necessarily need to be transferred, which is attributed to1104 !-- the design of the model coupler which allocatesmemory for each variable.1139 !-- Please note, if it is a nested run, arrays need to be allocated even if they do not necessarily 1140 !-- need to be transferred, which is attributed to the design of the model coupler which allocates 1141 !-- memory for each variable. 1105 1142 ALLOCATE( diss_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1106 1143 … … 1122 1159 1123 1160 1124 !------------------------------------------------------------------------------ !1161 !--------------------------------------------------------------------------------------------------! 1125 1162 ! Description: 1126 1163 ! ------------ 1127 1164 !> Initialization of turbulence closure module. 1128 !------------------------------------------------------------------------------ !1165 !--------------------------------------------------------------------------------------------------! 1129 1166 SUBROUTINE tcm_init 1130 1167 1131 USE control_parameters, & 1132 ONLY: bc_dirichlet_l, complex_terrain, topography 1133 1134 USE model_1d_mod, & 1135 ONLY: e1d, kh1d, km1d 1168 USE control_parameters, & 1169 ONLY: bc_dirichlet_l, & 1170 complex_terrain, & 1171 topography 1172 1173 USE model_1d_mod, & 1174 ONLY: e1d, & 1175 kh1d, & 1176 km1d 1136 1177 1137 1178 IMPLICIT NONE 1138 1179 1139 INTEGER(iwp) :: i !< loop index1140 INTEGER(iwp) :: j !< loop index1141 INTEGER(iwp) :: k !< loop index1142 INTEGER(iwp) :: nz_s_shift !< lower shift index for scalars1143 INTEGER(iwp) :: nz_s_shift_l !< local lower shift index in case of turbulent inflow1180 INTEGER(iwp) :: i !< loop index 1181 INTEGER(iwp) :: j !< loop index 1182 INTEGER(iwp) :: k !< loop index 1183 INTEGER(iwp) :: nz_s_shift !< lower shift index for scalars 1184 INTEGER(iwp) :: nz_s_shift_l !< local lower shift index in case of turbulent inflow 1144 1185 1145 1186 ! … … 1149 1190 ! 1150 1191 !-- Actions for initial runs 1151 IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &1192 IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. & 1152 1193 TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN 1153 1194 1154 1195 IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN 1155 1196 1156 IF ( .NOT. rans_tke_e ) THEN1197 IF ( .NOT. rans_tke_e ) THEN 1157 1198 ! 1158 1199 !-- Transfer initial profiles to the arrays of the 3D model … … 1173 1214 ELSE 1174 1215 ! 1175 !-- In case of TKE-e closure in RANS mode, do not use e, diss, and km 1176 !-- profiles from 1Dmodel. Instead, initialize with constant profiles1216 !-- In case of TKE-e closure in RANS mode, do not use e, diss, and km profiles from 1D 1217 !-- model. Instead, initialize with constant profiles 1177 1218 IF ( constant_diffusion ) THEN 1178 1219 km = km_constant … … 1183 1224 DO j = nysg, nyng 1184 1225 DO k = nzb+1, nzt 1185 km(k,j,i) = c_0 * SQRT( e_init ) * MIN( delta(k,j,i), & 1186 ml_blackadar(k) ) 1226 km(k,j,i) = c_0 * SQRT( e_init ) * MIN( delta(k,j,i), ml_blackadar(k) ) 1187 1227 ENDDO 1188 1228 ENDDO … … 1194 1234 ELSE 1195 1235 IF ( .NOT. ocean_mode ) THEN 1196 kh = 0.01_wp ! there must exist an initial diffusion, because 1197 km = 0.01_wp ! otherwise no TKE would be produced by the 1198 ! production terms, as long as not yet 1199 ! e = (u*/cm)**2 at k=nzb+1 1236 kh = 0.01_wp ! There must exist an initial diffusion, because otherwise no 1237 km = 0.01_wp ! TKE would be produced by the production terms, as long as not 1238 ! yet e = (u*/cm)**2 at k=nzb+1 1200 1239 ELSE 1201 1240 kh = 0.00001_wp … … 1217 1256 ENDIF 1218 1257 1219 ELSEIF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 .OR. &1258 ELSEIF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 .OR. & 1220 1259 INDEX( initializing_actions, 'inifor' ) /= 0 ) THEN 1221 1260 … … 1238 1277 ELSE 1239 1278 IF ( .NOT. ocean_mode ) THEN 1240 kh = 0.01_wp ! there must exist an initial diffusion, because 1241 km = 0.01_wp ! otherwise no TKE would be produced by the 1242 ! production terms, as long as not yet 1279 kh = 0.01_wp ! There must exist an initial diffusion, because otherwise no TKE 1280 km = 0.01_wp ! would be produced by the production terms, as long as not yet 1243 1281 ! e = (u*/cm)**2 at k=nzb+1 1244 1282 ELSE … … 1277 1315 ENDIF 1278 1316 1279 ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' .OR. & 1280 TRIM( initializing_actions ) == 'cyclic_fill' ) & 1281 THEN 1282 1283 ! 1284 !-- In case of complex terrain and cyclic fill method as initialization, 1285 !-- shift initial data in the vertical direction for each point in the 1286 !-- x-y-plane depending on local surface height 1287 IF ( complex_terrain .AND. & 1288 TRIM( initializing_actions ) == 'cyclic_fill' ) THEN 1317 ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' .OR. & 1318 TRIM( initializing_actions ) == 'cyclic_fill' ) THEN 1319 1320 ! 1321 !-- In case of complex terrain and cyclic fill method as initialization, shift initial data in 1322 !-- the vertical direction for each point in the x-y-plane depending on local surface height 1323 IF ( complex_terrain .AND. TRIM( initializing_actions ) == 'cyclic_fill' ) THEN 1289 1324 DO i = nxlg, nxrg 1290 1325 DO j = nysg, nyng … … 1311 1346 ! 1312 1347 !-- Initialization of the turbulence recycling method 1313 IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. & 1314 turbulent_inflow ) THEN 1348 IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. turbulent_inflow ) THEN 1315 1349 mean_inflow_profiles(:,5) = hom_sum(:,8,0) ! e 1316 1350 ! 1317 !-- In case of complex terrain, determine vertical displacement at inflow 1318 !-- boundary and adjustmean inflow profiles1351 !-- In case of complex terrain, determine vertical displacement at inflow boundary and adjust 1352 !-- mean inflow profiles 1319 1353 IF ( complex_terrain ) THEN 1320 IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. & 1321 nysg <= 0 .AND. nyng >= 0 ) THEN 1354 IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. nysg <= 0 .AND. nyng >= 0 ) THEN 1322 1355 nz_s_shift_l = topo_top_ind(0,0,0) 1323 1356 ELSE … … 1325 1358 ENDIF 1326 1359 #if defined( __parallel ) 1327 CALL MPI_ALLREDUCE(nz_s_shift_l, nz_s_shift, 1, MPI_INTEGER, & 1328 MPI_MAX, comm2d, ierr) 1360 CALL MPI_ALLREDUCE(nz_s_shift_l, nz_s_shift, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr) 1329 1361 #else 1330 1362 nz_s_shift = nz_s_shift_l 1331 1363 #endif 1332 mean_inflow_profiles(nz_s_shift:nzt+1,5) = & 1333 hom_sum(0:nzt+1-nz_s_shift,8,0) ! e 1364 mean_inflow_profiles(nz_s_shift:nzt+1,5) = hom_sum(0:nzt+1-nz_s_shift,8,0) ! e 1334 1365 ENDIF 1335 1366 ! 1336 !-- Use these mean profiles at the inflow (provided that Dirichlet 1337 !-- conditions are used) 1367 !-- Use these mean profiles at the inflow (provided that Dirichlet conditions are used) 1338 1368 IF ( bc_dirichlet_l ) THEN 1339 1369 DO j = nysg, nyng … … 1346 1376 ! 1347 1377 !-- Inside buildings set TKE back to zero 1348 IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. & 1349 topography /= 'flat' ) THEN 1378 IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. topography /= 'flat' ) THEN 1350 1379 ! 1351 1380 !-- Inside buildings set TKE back to zero. 1352 !-- Other scalars (km, kh,...) are ignored at present, 1353 !-- maybe revise later. 1381 !-- Other scalars (km, kh,...) are ignored at present, maybe revise later. 1354 1382 DO i = nxlg, nxrg 1355 1383 DO j = nysg, nyng 1356 1384 DO k = nzb, nzt 1357 e(k,j,i) = MERGE( e(k,j,i), 0.0_wp, & 1358 BTEST( wall_flags_total_0(k,j,i), 0 ) ) 1385 e(k,j,i) = MERGE( e(k,j,i), 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 1359 1386 ENDDO 1360 1387 ENDDO … … 1365 1392 DO j = nysg, nyng 1366 1393 DO k = nzb, nzt 1367 diss(k,j,i) = MERGE( diss(k,j,i), 0.0_wp,&1368 1394 diss(k,j,i) = MERGE( diss(k,j,i), 0.0_wp, & 1395 BTEST( wall_flags_total_0(k,j,i), 0 ) ) 1369 1396 ENDDO 1370 1397 ENDDO … … 1373 1400 ENDIF 1374 1401 ! 1375 !-- Initialize new time levels (only done in order to set boundary values 1376 !-- including ghost points) 1402 !-- Initialize new time levels (only done in order to set boundary values including ghost points) 1377 1403 e_p = e 1378 1404 ! 1379 !-- Allthough tendency arrays are set in prognostic_equations, they have 1380 !-- to be predefined here because there they are used (but multiplied with 0) 1381 !-- before they are set. 1405 !-- Although tendency arrays are set in prognostic_equations, they have to be predefined here 1406 !-- because there they are used (but multiplied with 0) before they are set. 1382 1407 te_m = 0.0_wp 1383 1408 … … 1392 1417 1393 1418 1394 !------------------------------------------------------------------------------ !1419 !--------------------------------------------------------------------------------------------------! 1395 1420 ! Description: 1396 1421 ! ------------ 1397 1422 !> Pre-computation of grid-dependent and near-wall mixing length. 1398 !> @todo consider walls in horizontal direction at a distance further than a 1399 !> single grid point(RANS mode)1400 !------------------------------------------------------------------------------ !1423 !> @todo consider walls in horizontal direction at a distance further than a single grid point 1424 !> (RANS mode) 1425 !--------------------------------------------------------------------------------------------------! 1401 1426 SUBROUTINE tcm_init_mixing_length 1402 1427 1403 USE arrays_3d, & 1404 ONLY: dzw, ug, vg, zu, zw 1405 1406 USE control_parameters, & 1407 ONLY: f, message_string, wall_adjustment 1408 1409 USE exchange_horiz_mod, & 1410 ONLY: exchange_horiz, exchange_horiz_int 1411 1412 USE grid_variables, & 1413 ONLY: dx, dy 1414 1415 USE indices, & 1416 ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, & 1417 nzt, wall_flags_total_0 1428 USE arrays_3d, & 1429 ONLY: dzw, & 1430 ug, & 1431 vg, & 1432 zu, & 1433 zw 1434 1435 USE control_parameters, & 1436 ONLY: f, & 1437 message_string, & 1438 wall_adjustment 1439 1440 USE exchange_horiz_mod, & 1441 ONLY: exchange_horiz, & 1442 exchange_horiz_int 1443 1444 USE grid_variables, & 1445 ONLY: dx, & 1446 dy 1447 1448 USE indices, & 1449 ONLY: nbgp, & 1450 nx, & 1451 nxl, & 1452 nxlg, & 1453 nxr, & 1454 nxrg, & 1455 ny, & 1456 nyn, & 1457 nyng, & 1458 nys, & 1459 nysg, & 1460 nzb, & 1461 nzt, & 1462 wall_flags_total_0 1418 1463 1419 1464 USE kinds … … 1422 1467 IMPLICIT NONE 1423 1468 1424 INTEGER(iwp) :: dist_dx !< found distance devided by dx1425 INTEGER(iwp) :: i !< index variable along x1426 INTEGER(iwp) :: ii !< index variable along x1427 INTEGER(iwp) :: j !< index variable along y1428 INTEGER(iwp) :: k !< index variable along z1429 INTEGER(iwp) :: k_max_topo !< index of maximum topography height1430 INTEGER(iwp) :: kk !< index variable along z1431 INTEGER(iwp) :: rad_i !< search radius in grid points along x1432 INTEGER(iwp) :: rad_i_l !< possible search radius to the left1433 INTEGER(iwp) :: rad_i_r !< possible search radius to the right1434 INTEGER(iwp) :: rad_j !< search radius in grid points along y1435 INTEGER(iwp) :: rad_j_n !< possible search radius to north1436 INTEGER(iwp) :: rad_j_s !< possible search radius to south1437 INTEGER(iwp) :: rad_k !< search radius in grid points along z1438 INTEGER(iwp) :: rad_k_b !< search radius in grid points along negative z1439 INTEGER(iwp) :: rad_k_t !< search radius in grid points along positive z1440 1441 INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE :: vic_yz!< contains a quarter of a single yz-slice of vicinity1442 1443 INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity!< contains topography information of the vicinity of (i/j/k)1444 1445 REAL(wp) :: distance_up !< distance of grid box center to its boundary in upper direction1446 REAL(wp) :: distance_down !< distance of grid box center to its boundary in lower direction1447 REAL(wp) :: distance_ns !< distance of grid box center to its boundary in y direction1448 REAL(wp) :: distance_lr !< distance of grid box center to its boundary in x direction1449 REAL(wp) :: distance_edge_yz_down !< distance of grid box center to its boundary along yz diagonal (down)1450 REAL(wp) :: distance_edge_yz_up !< distance of grid box center to its boundary along yz diagonal (up)1451 REAL(wp) :: distance_edge_xz_down !< distance of grid box center to its boundary along xz diagonal1452 REAL(wp) :: distance_edge_xz_up !< distance of grid box center to its boundary along xz diagonal (up)1453 REAL(wp) :: distance_edge_xy !< distance of grid box center to its boundary along xy diagonal1454 REAL(wp) :: distance_corners_down !< distance of grid box center to its upper corners1455 REAL(wp) :: distance_corners_up !< distance of grid box center to its lower corners1456 REAL(wp) :: radius !< search radius in meter1469 INTEGER(iwp) :: dist_dx !< found distance devided by dx 1470 INTEGER(iwp) :: i !< index variable along x 1471 INTEGER(iwp) :: ii !< index variable along x 1472 INTEGER(iwp) :: j !< index variable along y 1473 INTEGER(iwp) :: k !< index variable along z 1474 INTEGER(iwp) :: k_max_topo !< index of maximum topography height 1475 INTEGER(iwp) :: kk !< index variable along z 1476 INTEGER(iwp) :: rad_i !< search radius in grid points along x 1477 INTEGER(iwp) :: rad_i_l !< possible search radius to the left 1478 INTEGER(iwp) :: rad_i_r !< possible search radius to the right 1479 INTEGER(iwp) :: rad_j !< search radius in grid points along y 1480 INTEGER(iwp) :: rad_j_n !< possible search radius to north 1481 INTEGER(iwp) :: rad_j_s !< possible search radius to south 1482 INTEGER(iwp) :: rad_k !< search radius in grid points along z 1483 INTEGER(iwp) :: rad_k_b !< search radius in grid points along negative z 1484 INTEGER(iwp) :: rad_k_t !< search radius in grid points along positive z 1485 1486 INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE :: vic_yz !< contains a quarter of a single yz-slice of vicinity 1487 1488 INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity !< contains topography information of the vicinity of (i/j/k) 1489 1490 REAL(wp) :: distance_up !< distance of grid box center to its boundary in upper direction 1491 REAL(wp) :: distance_down !< distance of grid box center to its boundary in lower direction 1492 REAL(wp) :: distance_ns !< distance of grid box center to its boundary in y direction 1493 REAL(wp) :: distance_lr !< distance of grid box center to its boundary in x direction 1494 REAL(wp) :: distance_edge_yz_down !< distance of grid box center to its boundary along yz diagonal (down) 1495 REAL(wp) :: distance_edge_yz_up !< distance of grid box center to its boundary along yz diagonal (up) 1496 REAL(wp) :: distance_edge_xz_down !< distance of grid box center to its boundary along xz diagonal 1497 REAL(wp) :: distance_edge_xz_up !< distance of grid box center to its boundary along xz diagonal (up) 1498 REAL(wp) :: distance_edge_xy !< distance of grid box center to its boundary along xy diagonal 1499 REAL(wp) :: distance_corners_down !< distance of grid box center to its upper corners 1500 REAL(wp) :: distance_corners_up !< distance of grid box center to its lower corners 1501 REAL(wp) :: radius !< search radius in meter 1457 1502 1458 1503 REAL(wp), DIMENSION(nzb:nzt+1) :: gridsize_geometric_mean !< geometric mean of grid sizes dx, dy, dz … … 1476 1521 gridsize_geometric_mean(nzt+1) = gridsize_geometric_mean(nzt) 1477 1522 1478 IF ( ANY( gridsize_geometric_mean > 1.5_wp * dx * wall_adjustment_factor ) .OR.&1523 IF ( ANY( gridsize_geometric_mean > 1.5_wp * dx * wall_adjustment_factor ) .OR. & 1479 1524 ANY( gridsize_geometric_mean > 1.5_wp * dy * wall_adjustment_factor ) ) THEN 1480 WRITE( message_string, * ) 'grid anisotropy exceeds threshold', & 1481 ' &starting from height level k = ', k, & 1482 '.' 1525 WRITE( message_string, * ) 'grid anisotropy exceeds threshold', & 1526 ' &starting from height level k = ', k, '.' 1483 1527 CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 ) 1484 1528 ENDIF … … 1491 1535 1492 1536 ! 1493 !-- If Dai et al. (2020) closure is used, the distance to the wall (distance to nearest upward facing surface)1494 !-- must be stored1537 !-- If Dai et al. (2020) closure is used, the distance to the wall (distance to nearest upward 1538 !-- facing surface) must be stored 1495 1539 IF ( les_dai ) THEN 1496 1540 DO i = nxl, nxr … … 1505 1549 IF ( wall_adjustment ) THEN 1506 1550 ! 1507 !-- In case of topography, adjust mixing length if there is any wall at 1508 !-- the surrounding gridboxes:1551 !-- In case of topography, adjust mixing length if there is any wall at the surrounding grid 1552 !-- boxes: 1509 1553 !> @todo check if this is correct also for the ocean case 1510 1554 DO i = nxl, nxr … … 1515 1559 IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) ) THEN 1516 1560 ! 1517 !-- First, check if grid points directly next to current grid point 1518 !-- are surfacegrid points1561 !-- First, check if grid points directly next to current grid point are surface 1562 !-- grid points 1519 1563 !-- Check along... 1520 1564 !-- ...vertical direction, down … … 1533 1577 ! 1534 1578 !-- ...y-direction (north/south) 1535 IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i), 0 ) .OR. &1579 IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i), 0 ) .OR. & 1536 1580 .NOT. BTEST( wall_flags_total_0(k,j+1,i), 0 ) ) THEN 1537 1581 distance_ns = 0.5_wp * dy … … 1541 1585 ! 1542 1586 !-- ...x-direction (left/right) 1543 IF ( .NOT. BTEST( wall_flags_total_0(k,j,i-1), 0 ) .OR. &1587 IF ( .NOT. BTEST( wall_flags_total_0(k,j,i-1), 0 ) .OR. & 1544 1588 .NOT. BTEST( wall_flags_total_0(k,j,i+1), 0 ) ) THEN 1545 1589 distance_lr = 0.5_wp * dx … … 1550 1594 !-- Now, check the edges along... 1551 1595 !-- ...yz-direction (vertical edges, down) 1552 IF ( .NOT. BTEST( wall_flags_total_0(k-1,j-1,i), 0 ) .OR. &1596 IF ( .NOT. BTEST( wall_flags_total_0(k-1,j-1,i), 0 ) .OR. & 1553 1597 .NOT. BTEST( wall_flags_total_0(k-1,j+1,i), 0 ) ) THEN 1554 1598 distance_edge_yz_down = SQRT( 0.25_wp * dy**2 + ( zu(k) - zw(k-1) )**2 ) … … 1558 1602 ! 1559 1603 !-- ...yz-direction (vertical edges, up) 1560 IF ( .NOT. BTEST( wall_flags_total_0(k+1,j-1,i), 0 ) .OR. &1604 IF ( .NOT. BTEST( wall_flags_total_0(k+1,j-1,i), 0 ) .OR. & 1561 1605 .NOT. BTEST( wall_flags_total_0(k+1,j+1,i), 0 ) ) THEN 1562 1606 distance_edge_yz_up = SQRT( 0.25_wp * dy**2 + ( zw(k) - zu(k) )**2 ) … … 1566 1610 ! 1567 1611 !-- ...xz-direction (vertical edges, down) 1568 IF ( .NOT. BTEST( wall_flags_total_0(k-1,j,i-1), 0 ) .OR. &1612 IF ( .NOT. BTEST( wall_flags_total_0(k-1,j,i-1), 0 ) .OR. & 1569 1613 .NOT. BTEST( wall_flags_total_0(k-1,j,i+1), 0 ) ) THEN 1570 1614 distance_edge_xz_down = SQRT( 0.25_wp * dx**2 + ( zu(k) - zw(k-1) )**2 ) … … 1574 1618 ! 1575 1619 !-- ...xz-direction (vertical edges, up) 1576 IF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i-1), 0 ) .OR. &1620 IF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i-1), 0 ) .OR. & 1577 1621 .NOT. BTEST( wall_flags_total_0(k+1,j,i+1), 0 ) ) THEN 1578 1622 distance_edge_xz_up = SQRT( 0.25_wp * dx**2 + ( zw(k) - zu(k) )**2 ) … … 1582 1626 ! 1583 1627 !-- ...xy-direction (horizontal edges) 1584 IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i-1), 0 ) .OR. &1585 .NOT. BTEST( wall_flags_total_0(k,j+1,i-1), 0 ) .OR. &1586 .NOT. BTEST( wall_flags_total_0(k,j-1,i+1), 0 ) .OR. &1628 IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i-1), 0 ) .OR. & 1629 .NOT. BTEST( wall_flags_total_0(k,j+1,i-1), 0 ) .OR. & 1630 .NOT. BTEST( wall_flags_total_0(k,j-1,i+1), 0 ) .OR. & 1587 1631 .NOT. BTEST( wall_flags_total_0(k,j+1,i+1), 0 ) ) THEN 1588 1632 distance_edge_xy = SQRT( 0.25_wp * ( dx**2 + dy**2 ) ) … … 1593 1637 !-- Now, check the corners... 1594 1638 !-- ...lower four corners 1595 IF ( .NOT. BTEST( wall_flags_total_0(k-1,j-1,i-1), 0 ) .OR. &1596 .NOT. BTEST( wall_flags_total_0(k-1,j+1,i-1), 0 ) .OR. &1597 .NOT. BTEST( wall_flags_total_0(k-1,j-1,i+1), 0 ) .OR. &1639 IF ( .NOT. BTEST( wall_flags_total_0(k-1,j-1,i-1), 0 ) .OR. & 1640 .NOT. BTEST( wall_flags_total_0(k-1,j+1,i-1), 0 ) .OR. & 1641 .NOT. BTEST( wall_flags_total_0(k-1,j-1,i+1), 0 ) .OR. & 1598 1642 .NOT. BTEST( wall_flags_total_0(k-1,j+1,i+1), 0 ) ) THEN 1599 distance_corners_down = SQRT( 0.25_wp * ( dx**2 + dy**2 ) &1643 distance_corners_down = SQRT( 0.25_wp * ( dx**2 + dy**2 ) & 1600 1644 + ( zu(k) - zw(k-1) )**2 ) 1601 1645 ELSE … … 1604 1648 ! 1605 1649 !-- ...upper four corners 1606 IF ( .NOT. BTEST( wall_flags_total_0(k+1,j-1,i-1), 0 ) .OR. &1607 .NOT. BTEST( wall_flags_total_0(k+1,j+1,i-1), 0 ) .OR. &1608 .NOT. BTEST( wall_flags_total_0(k+1,j-1,i+1), 0 ) .OR. &1650 IF ( .NOT. BTEST( wall_flags_total_0(k+1,j-1,i-1), 0 ) .OR. & 1651 .NOT. BTEST( wall_flags_total_0(k+1,j+1,i-1), 0 ) .OR. & 1652 .NOT. BTEST( wall_flags_total_0(k+1,j-1,i+1), 0 ) .OR. & 1609 1653 .NOT. BTEST( wall_flags_total_0(k+1,j+1,i+1), 0 ) ) THEN 1610 distance_corners_up = SQRT( 0.25_wp * ( dx**2 + dy**2 ) &1654 distance_corners_up = SQRT( 0.25_wp * ( dx**2 + dy**2 ) & 1611 1655 + ( zw(k) - zu(k) )**2 ) 1612 1656 ELSE … … 1615 1659 1616 1660 ! 1617 !-- Calculate the minimum distance from the wall and store it 1618 !-- temporarily in the array delta 1619 delta(k,j,i) = MIN( & 1620 distance_up, distance_down, distance_ns, distance_lr, & 1621 distance_edge_yz_down, distance_edge_yz_up, & 1622 distance_edge_xz_down, distance_edge_xz_up, & 1623 distance_edge_xy, & 1624 distance_corners_down, distance_corners_up ) 1625 1626 ! 1627 !-- If Dai et al. (2020) closure is used, the distance to the wall 1628 !-- must be permanently stored 1661 !-- Calculate the minimum distance from the wall and store it temporarily in the 1662 !-- array delta 1663 delta(k,j,i) = MIN( distance_up, distance_down, distance_ns, distance_lr, & 1664 distance_edge_yz_down, distance_edge_yz_up, & 1665 distance_edge_xz_down, distance_edge_xz_up, & 1666 distance_edge_xy, distance_corners_down, & 1667 distance_corners_up ) 1668 1669 ! 1670 !-- If Dai et al. (2020) closure is used, the distance to the wall must be 1671 !-- permanently stored 1629 1672 IF ( les_dai .AND. delta(k,j,i) /= 9999999.9_wp ) THEN 1630 1673 distance_to_wall(k,j,i) = delta(k,j,i) … … 1635 1678 delta(k,j,i) = wall_adjustment_factor * delta(k,j,i) 1636 1679 1637 ENDIF ! if grid point belongs to atmosphere1638 1639 1640 1641 ! 1642 !-- The grid size (delta) is defined as the the minimum of the distance to 1643 !-- the nearestwall * 1.8 and the geometric mean grid size.1680 ENDIF ! If grid point belongs to atmosphere 1681 1682 1683 1684 ! 1685 !-- The grid size (delta) is defined as the the minimum of the distance to the nearest 1686 !-- wall * 1.8 and the geometric mean grid size. 1644 1687 delta(k,j,i) = MIN( delta(k,j,i), gridsize_geometric_mean(k) ) 1645 1688 … … 1659 1702 !-- Calculate mixing length according to Blackadar (1962) 1660 1703 IF ( f /= 0.0_wp ) THEN 1661 length_scale_max = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) & 1662 / ABS( f ) + 1.0E-10_wp 1704 length_scale_max = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) / ABS( f ) + 1.0E-10_wp 1663 1705 ELSE 1664 1706 length_scale_max = 30.0_wp … … 1677 1719 DO j = nysg, nyng 1678 1720 DO k = nzb+1, nzt-1 1679 IF ( .NOT. BTEST( wall_flags_total_0(k,j,i), 0 ) .AND. & 1680 k > k_max_topo ) & 1721 IF ( .NOT. BTEST( wall_flags_total_0(k,j,i), 0 ) .AND. k > k_max_topo ) & 1681 1722 k_max_topo = k 1682 1723 ENDDO … … 1687 1728 delta(nzt+1,:,:) = ml_blackadar(nzt+1) 1688 1729 ! 1689 !-- Limit mixing length to either nearest wall or Blackadar mixing length. 1690 !-- For that, analyze each grid point (i/j/k) ("analysed grid point") and 1691 !-- search within its vicinity for the shortest distance to a wall by cal- 1692 !-- culating the distance between the analysed grid point and the "viewed 1693 !-- grid point" if it contains a wall (belongs to topography). 1730 !-- Limit mixing length to either nearest wall or Blackadar mixing length. For that, analyze each 1731 !-- grid point (i/j/k) ("analysed grid point") and search within its vicinity for the shortest 1732 !-- distance to a wall by calculating the distance between the analysed grid point and the 1733 !-- "viewed grid point" if it contains a wall (belongs to topography). 1694 1734 DO k = nzb+1, nzt 1695 1735 … … 1720 1760 rad_k = MAX( rad_k_b, rad_k_t ) 1721 1761 ! 1722 !-- When analysed grid point lies above maximum topography, set search 1723 !-- radius to 0 if the distance between the analysed grid point and max1724 !-- topography height is larger than themaximum search radius1762 !-- When analysed grid point lies above maximum topography, set search radius to 0 if the 1763 !-- distance between the analysed grid point and max topography height is larger than the 1764 !-- maximum search radius 1725 1765 IF ( zu(k-rad_k_b) > zu(k_max_topo) ) rad_k_b = 0 1726 1766 ! … … 1729 1769 1730 1770 !> @note shape of vicinity is larger in z direction 1731 !> Shape of vicinity is two grid points larger than actual search 1732 !> radius in vertical direction. The first and last grid point is 1733 !> always set to 1 to asure correct detection of topography. See 1734 !> function "shortest_distance" for details. 1771 !> Shape of vicinity is two grid points larger than actual search radius in vertical 1772 !> direction. The first and last grid point is always set to 1 to asure correct 1773 !> detection of topography. See function "shortest_distance" for details. 1735 1774 !> 2018-03-16, gronemeier 1736 1775 ALLOCATE( vicinity(-rad_k-1:rad_k+1,-rad_j:rad_j,-rad_i:rad_i) ) … … 1748 1787 vicinity(-rad_k:rad_k,:,:) = 0 1749 1788 ! 1750 !-- Copy area surrounding analysed grid point into vicinity. 1751 !-- First, limit size of data copied to vicinity by the domain 1752 !-- border 1789 !-- Copy area surrounding analysed grid point into vicinity. First, limit size of 1790 !-- data copied to vicinity by the domain border 1753 1791 !> @note limit copied area to 1 grid point in hor. dir. 1754 !> Ignore walls in horizontal direction which are 1755 !> further away than a single grid point. This allows to 1756 !> only search within local subdomain without the need 1757 !> of global topography information. 1758 !> The error made by this assumption are acceptable at 1759 !> the moment. 1792 !> Ignore walls in horizontal direction which are further away than a single 1793 !> grid point. This allows to only search within local subdomain without the 1794 !> need of global topography information. The error made by this assumption 1795 !> are acceptable at the moment. 1760 1796 !> 2018-10-01, gronemeier 1761 1797 rad_i_l = MIN( 1, rad_i, i ) … … 1765 1801 rad_j_n = MIN( 1, rad_j, ny-j ) 1766 1802 1767 CALL copy_into_vicinity( k, j, i, & 1768 -rad_k_b, rad_k_t, & 1769 -rad_j_s, rad_j_n, & 1803 CALL copy_into_vicinity( k, j, i, -rad_k_b, rad_k_t, -rad_j_s, rad_j_n, & 1770 1804 -rad_i_l, rad_i_r ) 1771 !> @note in case of cyclic boundaries, those parts of the 1772 !> topography which lies beyond the domain borders but 1773 !> still within the search radius still needs to be 1774 !> copied into 'vicinity'. As the effective search 1775 !> radius is limited to 1 at the moment, no further 1776 !> copying is needed. Old implementation (prior to 1777 !> 2018-10-01) had this covered but used a global array. 1778 !> 2018-10-01, gronemeier 1805 !> @note in case of cyclic boundaries, those parts of the topography which 1806 !> lies beyond the domain borders but still within the search radius still 1807 !> needs to be copied into 'vicinity'. As the effective search radius is 1808 !> limited to 1 at the moment, no further copying is needed. Old 1809 !> implementation (prior to 2018-10-01) had this covered but used a global 1810 !> array. 1811 !> 2018-10-01, gronemeier 1779 1812 1780 1813 ! … … 1786 1819 DO ii = 0, dist_dx 1787 1820 ! 1788 !-- Search along vertical direction only if below 1789 !-- maximum topography 1790 IF ( rad_k_t > 0 ) THEN 1821 !-- Search along vertical direction only if below maximum topography 1822 IF ( rad_k_t > 0 ) THEN 1791 1823 ! 1792 1824 !-- Search for walls within octant (+++) 1793 1825 vic_yz = vicinity(0:rad_k+1,0:rad_j,ii) 1794 delta(k,j,i) = MIN( delta(k,j,i), &1795 shortest_distance( vic_yz, .TRUE., ii ) )1826 delta(k,j,i) = MIN( delta(k,j,i), & 1827 shortest_distance( vic_yz, .TRUE., ii ) ) 1796 1828 ! 1797 1829 !-- Search for walls within octant (+-+) 1798 !-- Switch order of array so that the analysed grid 1799 !-- point is always located at (0/0) (required by 1800 !-- shortest_distance"). 1830 !-- Switch order of array so that the analysed grid point is always 1831 !-- located at (0/0) (required by shortest_distance"). 1801 1832 vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii) 1802 delta(k,j,i) = MIN( delta(k,j,i), &1803 shortest_distance( vic_yz, .TRUE., ii ) )1833 delta(k,j,i) = MIN( delta(k,j,i), & 1834 shortest_distance( vic_yz, .TRUE., ii ) ) 1804 1835 1805 1836 ENDIF … … 1807 1838 !-- Search for walls within octant (+--) 1808 1839 vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii) 1809 delta(k,j,i) = MIN( delta(k,j,i), &1810 shortest_distance( vic_yz, .FALSE., ii ) )1840 delta(k,j,i) = MIN( delta(k,j,i), & 1841 shortest_distance( vic_yz, .FALSE., ii ) ) 1811 1842 ! 1812 1843 !-- Search for walls within octant (++-) 1813 1844 vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii) 1814 delta(k,j,i) = MIN( delta(k,j,i), &1815 shortest_distance( vic_yz, .FALSE., ii ) )1845 delta(k,j,i) = MIN( delta(k,j,i), & 1846 shortest_distance( vic_yz, .FALSE., ii ) ) 1816 1847 ! 1817 1848 !-- Reduce search along x by already found distance … … 1823 1854 DO ii = 0, -dist_dx, -1 1824 1855 ! 1825 !-- Search along vertical direction only if below 1826 !-- maximum topography 1827 IF ( rad_k_t > 0 ) THEN 1856 !-- Search along vertical direction only if below maximum topography 1857 IF ( rad_k_t > 0 ) THEN 1828 1858 ! 1829 1859 !-- Search for walls within octant (-++) 1830 1860 vic_yz = vicinity(0:rad_k+1,0:rad_j,ii) 1831 delta(k,j,i) = MIN( delta(k,j,i), &1832 shortest_distance( vic_yz, .TRUE., -ii ) )1861 delta(k,j,i) = MIN( delta(k,j,i), & 1862 shortest_distance( vic_yz, .TRUE., -ii ) ) 1833 1863 ! 1834 1864 !-- Search for walls within octant (--+) 1835 !-- Switch order of array so that the analysed grid 1836 !-- point is always located at (0/0) (required by 1837 !-- shortest_distance"). 1865 !-- Switch order of array so that the analysed grid point is always 1866 !-- located at (0/0) (required by shortest_distance"). 1838 1867 vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii) 1839 delta(k,j,i) = MIN( delta(k,j,i), &1840 shortest_distance( vic_yz, .TRUE., -ii ) )1868 delta(k,j,i) = MIN( delta(k,j,i), & 1869 shortest_distance( vic_yz, .TRUE., -ii ) ) 1841 1870 1842 1871 ENDIF … … 1844 1873 !-- Search for walls within octant (---) 1845 1874 vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii) 1846 delta(k,j,i) = MIN( delta(k,j,i), &1847 shortest_distance( vic_yz, .FALSE., -ii ) )1875 delta(k,j,i) = MIN( delta(k,j,i), & 1876 shortest_distance( vic_yz, .FALSE., -ii ) ) 1848 1877 ! 1849 1878 !-- Search for walls within octant (-+-) 1850 1879 vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii) 1851 delta(k,j,i) = MIN( delta(k,j,i), &1852 shortest_distance( vic_yz, .FALSE., -ii ) )1880 delta(k,j,i) = MIN( delta(k,j,i), & 1881 shortest_distance( vic_yz, .FALSE., -ii ) ) 1853 1882 ! 1854 1883 !-- Reduce search along x by already found distance … … 1871 1900 DEALLOCATE( vic_yz ) 1872 1901 1873 ENDIF ! check vertical size of vicinity1902 ENDIF !Check vertical size of vicinity 1874 1903 1875 1904 ENDDO !k loop … … 1892 1921 1893 1922 CONTAINS 1894 !------------------------------------------------------------------------------ !1923 !--------------------------------------------------------------------------------------------------! 1895 1924 ! Description: 1896 1925 ! ------------ 1897 !> Calculate the shortest distance between position (i/j/k)=(0/0/0) and 1898 !> (pos_i/jj/kk), where (jj/kk) is the position of the maximum of 'array' 1899 !> closest to the origin (0/0) of 'array'. 1900 !------------------------------------------------------------------------------! 1926 !> Calculate the shortest distance between position (i/j/k)=(0/0/0) and (pos_i/jj/kk), where (jj/kk) 1927 !> is the position of the maximum of 'array' closest to the origin (0/0) of 'array'. 1928 !--------------------------------------------------------------------------------------------------! 1901 1929 REAL(wp) FUNCTION shortest_distance( array, orientation, pos_i ) 1902 1930 1903 1931 IMPLICIT NONE 1904 1932 1905 LOGICAL, INTENT(IN) :: orientation !< flag if array represents an array oriented upwards (true) or downwards (false)1906 1907 INTEGER(iwp) , INTENT(IN) :: pos_i !< x position of the yz-plane 'array'1908 1909 INTEGER(iwp) :: a!< loop index1910 INTEGER(iwp) :: b !< loop index 1911 INTEGER( iwp) :: jj !< loop index1912 1913 INTEGER( KIND=1) :: maximum !< maximum of array along zdimension1914 1915 INTEGER( iwp), DIMENSION(0:rad_j) :: loc_k !< location of closest wall along vertical dimension1916 1917 INTEGER(KIND=1), DIMENSION(0:rad_k+1,0:rad_j), INTENT(IN) :: array !< array containing a yz-plane at position pos_i1918 1919 ! 1920 !-- Get coordinate of first maximum along vertical dimension 1921 !-- at each y position of array(similar to function maxloc but more stable).1933 INTEGER(iwp), INTENT(IN) :: pos_i !< x position of the yz-plane 'array' 1934 1935 INTEGER(iwp) :: a !< loop index 1936 INTEGER(iwp) :: b !< loop index 1937 INTEGER(iwp) :: jj !< loop index 1938 1939 INTEGER(KIND=1) :: maximum !< maximum of array along z dimension 1940 1941 INTEGER(iwp), DIMENSION(0:rad_j) :: loc_k !< location of closest wall along vertical dimension 1942 1943 INTEGER(KIND=1), DIMENSION(0:rad_k+1,0:rad_j), INTENT(IN) :: array !< array containing a yz-plane at position pos_i 1944 1945 LOGICAL, INTENT(IN) :: orientation !< flag if array represents an array oriented upwards (true) or downwards (false) 1946 1947 ! 1948 !-- Get coordinate of first maximum along vertical dimension at each y position of array 1949 !-- (similar to function maxloc but more stable). 1922 1950 DO a = 0, rad_j 1923 1951 loc_k(a) = rad_k+1 … … 1934 1962 shortest_distance = radius 1935 1963 ! 1936 !-- Calculate distance between position (0/0/0) and 1937 !-- position (pos_i/jj/loc(jj)) and only save theshortest distance.1938 IF ( orientation ) THEN !if array is oriented upwards1964 !-- Calculate distance between position (0/0/0) and position (pos_i/jj/loc(jj)) and only save the 1965 !-- shortest distance. 1966 IF ( orientation ) THEN !if array is oriented upwards 1939 1967 DO jj = 0, rad_j 1940 shortest_distance = &1941 MIN( shortest_distance,&1942 SQRT( MAX(REAL(pos_i, KIND=wp)*dx-0.5_wp*dx, 0.0_wp)**2&1943 + MAX(REAL(jj, KIND=wp)*dy-0.5_wp*dy, 0.0_wp)**2&1944 + MAX(zw(loc_k(jj)+k-1)-zu(k), 0.0_wp)**2&1945 )&1946 1968 shortest_distance = & 1969 MIN( shortest_distance, & 1970 SQRT( MAX( REAL( pos_i, KIND = wp ) * dx - 0.5_wp * dx, 0.0_wp)**2 & 1971 + MAX( REAL( jj, KIND = wp ) * dy - 0.5_wp * dy, 0.0_wp)**2 & 1972 + MAX( zw(loc_k(jj)+k-1) - zu(k), 0.0_wp)**2 & 1973 ) & 1974 ) 1947 1975 ENDDO 1948 1976 ELSE !if array is oriented downwards 1949 1977 !> @note MAX within zw required to circumvent error at domain border 1950 !> At the domain border, if non-cyclic boundary is present, the 1951 !> index for zw could be -1, which will be errorneous (zw(-1) does1952 !> not exist). The MAX function limits theindex to be at least 0.1978 !> At the domain border, if non-cyclic boundary is present, the index for zw could be 1979 !> -1, which will be errorneous (zw(-1) does not exist). The MAX function limits the 1980 !> index to be at least 0. 1953 1981 DO jj = 0, rad_j 1954 shortest_distance = &1955 MIN( shortest_distance, &1956 SQRT( MAX(REAL(pos_i, KIND=wp)*dx-0.5_wp*dx, 0.0_wp)**2 &1957 + MAX(REAL(jj, KIND=wp)*dy-0.5_wp*dy, 0.0_wp)**2 &1958 + MAX(zu(k)-zw(MAX(k-loc_k(jj),0_iwp)), 0.0_wp)**2 &1959 ) &1982 shortest_distance = & 1983 MIN( shortest_distance, & 1984 SQRT( MAX(REAL(pos_i, KIND=wp)*dx-0.5_wp*dx, 0.0_wp)**2 & 1985 + MAX(REAL(jj, KIND=wp)*dy-0.5_wp*dy, 0.0_wp)**2 & 1986 + MAX(zu(k)-zw(MAX(k-loc_k(jj),0_iwp)), 0.0_wp)**2 & 1987 ) & 1960 1988 ) 1961 1989 ENDDO … … 1964 1992 END FUNCTION 1965 1993 1966 !------------------------------------------------------------------------------ !1994 !--------------------------------------------------------------------------------------------------! 1967 1995 ! Description: 1968 1996 ! ------------ 1969 !> Copy a subarray of size (kb:kt,js:jn,il:ir) centered around grid point 1970 !> (kp,jp,ip) containing the first bit of wall_flags_total_0 into the array1971 !> 'vicinity'. Only copy first bit as this indicatesthe presence of topography.1972 !------------------------------------------------------------------------------ !1997 !> Copy a subarray of size (kb:kt,js:jn,il:ir) centered around grid point (kp,jp,ip) containing the 1998 !> first bit of wall_flags_total_0 into the array 'vicinity'. Only copy first bit as this indicates 1999 !> the presence of topography. 2000 !--------------------------------------------------------------------------------------------------! 1973 2001 SUBROUTINE copy_into_vicinity( kp, jp, ip, kb, kt, js, jn, il, ir ) 1974 2002 1975 2003 IMPLICIT NONE 1976 2004 1977 INTEGER(iwp), INTENT(IN) :: il!< left loop boundary1978 INTEGER(iwp), INTENT(IN) :: ip!< center position in x-direction1979 INTEGER(iwp), INTENT(IN) :: ir!< right loop boundary1980 INTEGER(iwp), INTENT(IN) :: jn!< northern loop boundary1981 INTEGER(iwp), INTENT(IN) :: jp!< center position in y-direction1982 INTEGER(iwp), INTENT(IN) :: js!< southern loop boundary1983 INTEGER(iwp), INTENT(IN) :: kb!< bottom loop boundary1984 INTEGER(iwp), INTENT(IN) :: kp!< center position in z-direction1985 INTEGER(iwp), INTENT(IN) :: kt!< top loop boundary1986 1987 INTEGER(iwp) :: i!< loop index1988 INTEGER(iwp) :: j!< loop index1989 INTEGER(iwp) :: k!< loop index2005 INTEGER(iwp), INTENT(IN) :: il !< left loop boundary 2006 INTEGER(iwp), INTENT(IN) :: ip !< center position in x-direction 2007 INTEGER(iwp), INTENT(IN) :: ir !< right loop boundary 2008 INTEGER(iwp), INTENT(IN) :: jn !< northern loop boundary 2009 INTEGER(iwp), INTENT(IN) :: jp !< center position in y-direction 2010 INTEGER(iwp), INTENT(IN) :: js !< southern loop boundary 2011 INTEGER(iwp), INTENT(IN) :: kb !< bottom loop boundary 2012 INTEGER(iwp), INTENT(IN) :: kp !< center position in z-direction 2013 INTEGER(iwp), INTENT(IN) :: kt !< top loop boundary 2014 2015 INTEGER(iwp) :: i !< loop index 2016 INTEGER(iwp) :: j !< loop index 2017 INTEGER(iwp) :: k !< loop index 1990 2018 1991 2019 DO i = il, ir 1992 2020 DO j = js, jn 1993 2021 DO k = kb, kt 1994 vicinity(k,j,i) = MERGE( 0, 1, & 1995 BTEST( wall_flags_total_0(kp+k,jp+j,ip+i), 0 ) ) 2022 vicinity(k,j,i) = MERGE( 0, 1, BTEST( wall_flags_total_0(kp+k,jp+j,ip+i), 0 ) ) 1996 2023 ENDDO 1997 2024 ENDDO … … 2003 2030 2004 2031 2005 !------------------------------------------------------------------------------ !2032 !--------------------------------------------------------------------------------------------------! 2006 2033 ! Description: 2007 2034 ! ------------ 2008 2035 !> Initialize virtual velocities used later in production_e. 2009 !------------------------------------------------------------------------------ !2036 !--------------------------------------------------------------------------------------------------! 2010 2037 SUBROUTINE production_e_init 2011 2038 2012 USE arrays_3d, & 2013 ONLY: drho_air_zw, zu 2014 2015 USE control_parameters, & 2039 USE arrays_3d, & 2040 ONLY: drho_air_zw, & 2041 zu 2042 2043 USE control_parameters, & 2016 2044 ONLY: constant_flux_layer 2017 2045 2018 USE surface_layer_fluxes_mod, &2046 USE surface_layer_fluxes_mod, & 2019 2047 ONLY: phi_m 2020 2048 2021 2049 IMPLICIT NONE 2022 2050 2023 INTEGER(iwp) :: i 2024 INTEGER(iwp) :: j 2025 INTEGER(iwp) :: k 2026 INTEGER(iwp) :: m 2027 2028 REAL(wp) :: km_sfc 2051 INTEGER(iwp) :: i !< grid index x-direction 2052 INTEGER(iwp) :: j !< grid index y-direction 2053 INTEGER(iwp) :: k !< grid index z-direction 2054 INTEGER(iwp) :: m !< running index surface elements 2055 2056 REAL(wp) :: km_sfc !< diffusion coefficient, used to compute virtual velocities 2029 2057 2030 2058 IF ( constant_flux_layer ) THEN 2031 2059 ! 2032 !-- Calculate a virtual velocity at the surface in a way that the 2033 !-- vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the 2034 !-- Prandtl law (-w'u'/km). This gradient is used in the TKE shear 2035 !-- production term at k=1 (see production_e_ij). 2036 !-- The velocity gradient has to be limited in case of too small km 2037 !-- (otherwise the timestep may be significantly reduced by large 2038 !-- surface winds). 2039 !-- not available in case of non-cyclic boundary conditions. 2060 !-- Calculate a virtual velocity at the surface in a way that the vertical velocity gradient at 2061 !-- k = 1 (u(k+1)-u_0) matches the Prandtl law (-w'u'/km). This gradient is used in the TKE shear 2062 !-- production term at k=1 (see production_e_ij). The velocity gradient has to be limited in case 2063 !-- of too small km (otherwise the timestep may be significantly reduced by large surface winds). 2064 !-- Not available in case of non-cyclic boundary conditions. 2040 2065 !-- Default surfaces, upward-facing 2041 2066 !$OMP PARALLEL DO PRIVATE(i,j,k,m) … … 2048 2073 k = surf_def_h(0)%k(m) 2049 2074 ! 2050 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v 2051 !-- and km are not on the same grid. Actually, a further 2052 !-- interpolation of km onto the u/v-grid is necessary. However, the 2075 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v and km are not on the same 2076 !-- grid. Actually, a further interpolation of km onto the u/v-grid is necessary. However, the 2053 2077 !-- effect of this error is negligible. 2054 km_sfc = kappa * surf_def_h(0)%us(m) * surf_def_h(0)%z_mo(m) / &2078 km_sfc = kappa * surf_def_h(0)%us(m) * surf_def_h(0)%z_mo(m) / & 2055 2079 phi_m( surf_def_h(0)%z_mo(m) / surf_def_h(0)%ol(m) ) 2056 2080 2057 surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) * & 2058 drho_air_zw(k-1) * & 2059 ( zu(k+1) - zu(k-1) ) / & 2060 ( km_sfc + 1.0E-20_wp ) 2061 surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) * & 2062 drho_air_zw(k-1) * & 2063 ( zu(k+1) - zu(k-1) ) / & 2064 ( km_sfc + 1.0E-20_wp ) 2065 2066 IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) > & 2067 ABS( u(k+1,j,i) - u(k-1,j,i) ) & 2068 ) surf_def_h(0)%u_0(m) = u(k-1,j,i) 2069 2070 IF ( ABS( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) > & 2071 ABS( v(k+1,j,i) - v(k-1,j,i) ) & 2072 ) surf_def_h(0)%v_0(m) = v(k-1,j,i) 2081 surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) * drho_air_zw(k-1) * & 2082 ( zu(k+1) - zu(k-1) ) / ( km_sfc + 1.0E-20_wp ) 2083 surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) * drho_air_zw(k-1) * & 2084 ( zu(k+1) - zu(k-1) ) / ( km_sfc + 1.0E-20_wp ) 2085 2086 IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) > ABS( u(k+1,j,i) - u(k-1,j,i) ) ) & 2087 surf_def_h(0)%u_0(m) = u(k-1,j,i) 2088 2089 IF ( ABS( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) > ABS( v(k+1,j,i) - v(k-1,j,i) ) ) & 2090 surf_def_h(0)%v_0(m) = v(k-1,j,i) 2073 2091 2074 2092 ENDDO … … 2084 2102 k = surf_def_h(1)%k(m) 2085 2103 ! 2086 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v 2087 !-- and km are not on the same grid. Actually, a further 2088 !-- interpolation of km onto the u/v-grid is necessary. However, the 2104 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v and km are not on the same 2105 !-- grid. Actually, a further interpolation of km onto the u/v-grid is necessary. However, the 2089 2106 !-- effect of this error is negligible. 2090 surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) * & 2091 drho_air_zw(k-1) * & 2092 ( zu(k+1) - zu(k-1) ) / & 2093 ( km(k,j,i) + 1.0E-20_wp ) 2094 surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) * & 2095 drho_air_zw(k-1) * & 2096 ( zu(k+1) - zu(k-1) ) / & 2097 ( km(k,j,i) + 1.0E-20_wp ) 2098 2099 IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) > & 2100 ABS( u(k+1,j,i) - u(k-1,j,i) ) & 2101 ) surf_def_h(1)%u_0(m) = u(k+1,j,i) 2102 2103 IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) > & 2104 ABS( v(k+1,j,i) - v(k-1,j,i) ) & 2105 ) surf_def_h(1)%v_0(m) = v(k+1,j,i) 2107 surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) * drho_air_zw(k-1) * & 2108 ( zu(k+1) - zu(k-1) ) / ( km(k,j,i) + 1.0E-20_wp ) 2109 surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) * drho_air_zw(k-1) * & 2110 ( zu(k+1) - zu(k-1) ) / (km(k,j,i) + 1.0E-20_wp ) 2111 2112 IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) > ABS( u(k+1,j,i) - u(k-1,j,i) ) ) & 2113 surf_def_h(1)%u_0(m) = u(k+1,j,i) 2114 2115 IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) > ABS( v(k+1,j,i) - v(k-1,j,i) ) ) & 2116 surf_def_h(1)%v_0(m) = v(k+1,j,i) 2106 2117 2107 2118 ENDDO … … 2117 2128 k = surf_lsm_h%k(m) 2118 2129 ! 2119 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v 2120 !-- and km are not on the same grid. Actually, a further 2121 !-- interpolation of km onto the u/v-grid is necessary. However, the 2130 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v and km are not on the same 2131 !-- grid. Actually, a further interpolation of km onto the u/v-grid is necessary. However, the 2122 2132 !-- effect of this error is negligible. 2123 km_sfc = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m) / &2133 km_sfc = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m) / & 2124 2134 phi_m( surf_lsm_h%z_mo(m) / surf_lsm_h%ol(m) ) 2125 2135 2126 surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m) * & 2127 drho_air_zw(k-1) * & 2128 ( zu(k+1) - zu(k-1) ) / & 2129 ( km_sfc + 1.0E-20_wp ) 2130 surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m) * & 2131 drho_air_zw(k-1) * & 2132 ( zu(k+1) - zu(k-1) ) / & 2133 ( km_sfc + 1.0E-20_wp ) 2134 2135 IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) ) > & 2136 ABS( u(k+1,j,i) - u(k-1,j,i) ) & 2137 ) surf_lsm_h%u_0(m) = u(k-1,j,i) 2138 2139 IF ( ABS( v(k+1,j,i) - surf_lsm_h%v_0(m) ) > & 2140 ABS( v(k+1,j,i) - v(k-1,j,i) ) & 2141 ) surf_lsm_h%v_0(m) = v(k-1,j,i) 2136 surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m) * drho_air_zw(k-1) * & 2137 ( zu(k+1) - zu(k-1) ) / ( km_sfc + 1.0E-20_wp ) 2138 surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m) * drho_air_zw(k-1) * & 2139 ( zu(k+1) - zu(k-1)) / ( km_sfc + 1.0E-20_wp ) 2140 2141 IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) ) > ABS( u(k+1,j,i) - u(k-1,j,i) ) ) & 2142 surf_lsm_h%u_0(m) = u(k-1,j,i) 2143 2144 IF ( ABS( v(k+1,j,i) - surf_lsm_h%v_0(m) ) > ABS( v(k+1,j,i) - v(k-1,j,i) ) ) & 2145 surf_lsm_h%v_0(m) = v(k-1,j,i) 2142 2146 2143 2147 ENDDO … … 2153 2157 k = surf_usm_h%k(m) 2154 2158 ! 2155 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v 2156 !-- and km are not on the same grid. Actually, a further 2157 !-- interpolation of km onto the u/v-grid is necessary. However, the 2159 !-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v and km are not on the same 2160 !-- grid. Actually, a further interpolation of km onto the u/v-grid is necessary. However, the 2158 2161 !-- effect of this error is negligible. 2159 km_sfc = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m) / &2162 km_sfc = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m) / & 2160 2163 phi_m( surf_usm_h%z_mo(m) / surf_usm_h%ol(m) ) 2161 2164 2162 surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m) * & 2163 drho_air_zw(k-1) * & 2164 ( zu(k+1) - zu(k-1) ) / & 2165 ( km_sfc + 1.0E-20_wp ) 2166 surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m) * & 2167 drho_air_zw(k-1) * & 2168 ( zu(k+1) - zu(k-1) ) / & 2169 ( km_sfc + 1.0E-20_wp ) 2170 2171 IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) ) > & 2172 ABS( u(k+1,j,i) - u(k-1,j,i) ) & 2173 ) surf_usm_h%u_0(m) = u(k-1,j,i) 2174 2175 IF ( ABS( v(k+1,j,i) - surf_usm_h%v_0(m) ) > & 2176 ABS( v(k+1,j,i) - v(k-1,j,i) ) & 2177 ) surf_usm_h%v_0(m) = v(k-1,j,i) 2165 surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m) * drho_air_zw(k-1) * & 2166 ( zu(k+1) - zu(k-1) ) / ( km_sfc + 1.0E-20_wp ) 2167 surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m) * drho_air_zw(k-1) * & 2168 ( zu(k+1) - zu(k-1) ) / ( km_sfc + 1.0E-20_wp ) 2169 2170 IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) ) > ABS( u(k+1,j,i) - u(k-1,j,i) ) ) & 2171 surf_usm_h%u_0(m) = u(k-1,j,i) 2172 2173 IF ( ABS( v(k+1,j,i) - surf_usm_h%v_0(m) ) > ABS( v(k+1,j,i) - v(k-1,j,i) ) ) & 2174 surf_usm_h%v_0(m) = v(k-1,j,i) 2178 2175 2179 2176 ENDDO … … 2192 2189 2193 2190 2194 CHARACTER 2191 CHARACTER(LEN=*) :: location !< 2195 2192 2196 2193 ! INTEGER(iwp) :: i !< … … 2200 2197 ! 2201 2198 !-- Here the module-specific actions follow 2202 !-- No calls for single grid points are allowed at locations before and 2203 !-- after the timestep, since these calls are not within an i,j-loop2199 !-- No calls for single grid points are allowed at locations before and after the timestep, since 2200 !-- these calls are not within an i, j-loop 2204 2201 SELECT CASE ( location ) 2205 2202 … … 2258 2255 2259 2256 2260 CHARACTER (LEN=*) :: location2261 2262 INTEGER(iwp) :: i 2263 INTEGER(iwp) :: j 2257 CHARACTER(LEN=*) :: location !< 2258 2259 INTEGER(iwp) :: i !< 2260 INTEGER(iwp) :: j !< 2264 2261 2265 2262 ! … … 2301 2298 2302 2299 2303 !------------------------------------------------------------------------------ !2300 !--------------------------------------------------------------------------------------------------! 2304 2301 ! Description: 2305 2302 ! ------------ 2306 !> Prognostic equation for subgrid-scale TKE and TKE dissipation rate. 2307 !> Vector-optimized version 2308 !------------------------------------------------------------------------------! 2303 !> Prognostic equation for subgrid-scale TKE and TKE dissipation rate. Vector-optimized version 2304 !--------------------------------------------------------------------------------------------------! 2309 2305 SUBROUTINE tcm_prognostic_equations 2310 2306 2311 USE control_parameters, & 2312 ONLY: scalar_advec, tsc 2307 USE control_parameters, & 2308 ONLY: scalar_advec, & 2309 tsc 2313 2310 2314 2311 IMPLICIT NONE 2315 2312 2316 INTEGER(iwp) :: i !< loop index 2317 INTEGER(iwp) :: j !< loop index 2318 INTEGER(iwp) :: k !< loop index 2319 2320 REAL(wp) :: sbt !< wheighting factor for sub-time step 2321 2322 ! 2323 !-- If required, compute prognostic equation for turbulent kinetic 2324 !-- energy (TKE) 2313 INTEGER(iwp) :: i !< loop index 2314 INTEGER(iwp) :: j !< loop index 2315 INTEGER(iwp) :: k !< loop index 2316 2317 REAL(wp) :: sbt !< wheighting factor for sub-time step 2318 2319 ! 2320 !-- If required, compute prognostic equation for turbulent kinetic energy (TKE) 2325 2321 IF ( .NOT. constant_diffusion ) THEN 2326 2322 … … 2354 2350 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2355 2351 IF ( ws_scheme_sca ) THEN 2356 CALL advec_s_ws( advc_flags_s, e, 'e', & 2357 bc_dirichlet_l .OR. bc_radiation_l, & 2358 bc_dirichlet_n .OR. bc_radiation_n, & 2359 bc_dirichlet_r .OR. bc_radiation_r, & 2360 bc_dirichlet_s .OR. bc_radiation_s ) 2352 CALL advec_s_ws( advc_flags_s, e, 'e', bc_dirichlet_l .OR. bc_radiation_l, & 2353 bc_dirichlet_n .OR. bc_radiation_n, & 2354 bc_dirichlet_r .OR. bc_radiation_r, & 2355 bc_dirichlet_s .OR. bc_radiation_s ) 2361 2356 ELSE 2362 2357 CALL advec_s_pw( e ) … … 2389 2384 ! 2390 2385 !-- Prognostic equation for TKE. 2391 !-- Eliminate negative TKE values, which can occur due to numerical 2392 !-- reasons in the course of the integration. In such cases the old TKE 2393 !-- value is reduced by 90%. 2386 !-- Eliminate negative TKE values, which can occur due to numerical reasons in the course of the 2387 !-- integration. In such cases the old TKE value is reduced by 90%. 2394 2388 !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) & 2395 2389 !$ACC PRESENT(e, tend, te_m, wall_flags_total_0) & … … 2398 2392 DO i = nxl, nxr 2399 2393 DO j = nys, nyn 2400 ! following directive is required to vectorize on Intel192394 !Following directive is required to vectorize on Intel19 2401 2395 !DIR$ IVDEP 2402 2396 DO k = nzb+1, nzt 2403 e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + & 2404 tsc(3) * te_m(k,j,i) ) & 2405 ) & 2406 * MERGE( 1.0_wp, 0.0_wp, & 2407 BTEST( wall_flags_total_0(k,j,i), 0 ) & 2408 ) 2397 e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) ) ) & 2398 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 2409 2399 IF ( e_p(k,j,i) < 0.0_wp ) e_p(k,j,i) = 0.1_wp * e(k,j,i) 2410 2400 ENDDO … … 2425 2415 ENDDO 2426 2416 ENDDO 2427 ELSEIF ( intermediate_timestep_count < & 2428 intermediate_timestep_count_max ) THEN 2417 ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) THEN 2429 2418 !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) & 2430 2419 !$ACC PRESENT(tend, te_m) … … 2432 2421 DO j = nys, nyn 2433 2422 DO k = nzb+1, nzt 2434 te_m(k,j,i) = -9.5625_wp * tend(k,j,i) & 2435 + 5.3125_wp * te_m(k,j,i) 2423 te_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * te_m(k,j,i) 2436 2424 ENDDO 2437 2425 ENDDO … … 2466 2454 2467 2455 ! 2468 !-- dissipation-tendency terms with no communication2456 !-- Dissipation-tendency terms with no communication 2469 2457 IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN 2470 2458 IF ( use_upstream_for_tke ) THEN … … 2475 2463 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2476 2464 IF ( ws_scheme_sca ) THEN 2477 CALL advec_s_ws( advc_flags_s, diss, 'diss', &2478 bc_dirichlet_l .OR. bc_radiation_l, &2479 bc_dirichlet_n .OR. bc_radiation_n, &2480 bc_dirichlet_r .OR. bc_radiation_r, &2465 CALL advec_s_ws( advc_flags_s, diss, 'diss', & 2466 bc_dirichlet_l .OR. bc_radiation_l, & 2467 bc_dirichlet_n .OR. bc_radiation_n, & 2468 bc_dirichlet_r .OR. bc_radiation_r, & 2481 2469 bc_dirichlet_s .OR. bc_radiation_s ) 2482 2470 ELSE … … 2502 2490 ! 2503 2491 !-- Prognostic equation for TKE dissipation. 2504 !-- Eliminate negative dissipation values, which can occur due to numerical 2505 !-- reasons in the course of the integration. In such cases the old 2506 !-- dissipation value is reduced by 90%. 2492 !-- Eliminate negative dissipation values, which can occur due to numerical reasons in the course 2493 !-- of the integration. In such cases the old dissipation value is reduced by 90%. 2507 2494 DO i = nxl, nxr 2508 2495 DO j = nys, nyn 2509 2496 DO k = nzb+1, nzt 2510 diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + & 2511 tsc(3) * tdiss_m(k,j,i) ) & 2512 ) & 2513 * MERGE( 1.0_wp, 0.0_wp, & 2514 BTEST( wall_flags_total_0(k,j,i), 0 ) & 2515 ) 2516 IF ( diss_p(k,j,i) < 0.0_wp ) & 2517 diss_p(k,j,i) = 0.1_wp * diss(k,j,i) 2497 diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + tsc(3) & 2498 * tdiss_m(k,j,i) ) ) & 2499 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 2500 IF ( diss_p(k,j,i) < 0.0_wp ) diss_p(k,j,i) = 0.1_wp * diss(k,j,i) 2518 2501 ENDDO 2519 2502 ENDDO … … 2531 2514 ENDDO 2532 2515 ENDDO 2533 ELSEIF ( intermediate_timestep_count < & 2534 intermediate_timestep_count_max ) THEN 2516 ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) THEN 2535 2517 DO i = nxl, nxr 2536 2518 DO j = nys, nyn 2537 2519 DO k = nzb+1, nzt 2538 tdiss_m(k,j,i) = -9.5625_wp * tend(k,j,i) & 2539 + 5.3125_wp * tdiss_m(k,j,i) 2520 tdiss_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tdiss_m(k,j,i) 2540 2521 ENDDO 2541 2522 ENDDO … … 2551 2532 2552 2533 2553 !------------------------------------------------------------------------------ !2534 !--------------------------------------------------------------------------------------------------! 2554 2535 ! Description: 2555 2536 ! ------------ 2556 !> Prognostic equation for subgrid-scale TKE and TKE dissipation rate. 2557 !> Cache-optimized version 2558 !------------------------------------------------------------------------------! 2537 !> Prognostic equation for subgrid-scale TKE and TKE dissipation rate. Cache-optimized version 2538 !--------------------------------------------------------------------------------------------------! 2559 2539 SUBROUTINE tcm_prognostic_equations_ij( i, j, i_omp, tn ) 2560 2540 2561 USE arrays_3d, & 2562 ONLY: diss_l_diss, diss_l_e, diss_s_diss, diss_s_e, flux_l_diss, & 2563 flux_l_e, flux_s_diss, flux_s_e 2564 2565 USE control_parameters, & 2541 USE arrays_3d, & 2542 ONLY: diss_l_diss, & 2543 diss_l_e, & 2544 diss_s_diss, & 2545 diss_s_e, & 2546 flux_l_diss, & 2547 flux_l_e, & 2548 flux_s_diss, & 2549 flux_s_e 2550 2551 USE control_parameters, & 2566 2552 ONLY: tsc 2567 2553 2568 2554 IMPLICIT NONE 2569 2555 2570 INTEGER(iwp) :: i !< loop index x direction 2571 INTEGER(iwp) :: i_omp !< first loop index of i-loop in prognostic_equations 2572 INTEGER(iwp) :: j !< loop index y direction 2573 INTEGER(iwp) :: k !< loop index z direction 2574 INTEGER(iwp) :: tn !< task number of openmp task 2575 2576 ! 2577 !-- If required, compute prognostic equation for turbulent kinetic 2578 !-- energy (TKE) 2556 INTEGER(iwp) :: i !< loop index x direction 2557 INTEGER(iwp) :: i_omp !< first loop index of i-loop in prognostic_equations 2558 INTEGER(iwp) :: j !< loop index y direction 2559 INTEGER(iwp) :: k !< loop index z direction 2560 INTEGER(iwp) :: tn !< task number of openmp task 2561 2562 ! 2563 !-- If required, compute prognostic equation for turbulent kinetic energy (TKE) 2579 2564 IF ( .NOT. constant_diffusion ) THEN 2580 2565 … … 2582 2567 !-- Tendency-terms for TKE 2583 2568 tend(:,j,i) = 0.0_wp 2584 IF ( timestep_scheme(1:5) == 'runge' & 2585 .AND. .NOT. use_upstream_for_tke ) THEN 2569 IF ( timestep_scheme(1:5) == 'runge' .AND. .NOT. use_upstream_for_tke ) THEN 2586 2570 IF ( ws_scheme_sca ) THEN 2587 CALL advec_s_ws( advc_flags_s, & 2588 i, j, e, 'e', flux_s_e, diss_s_e, & 2589 flux_l_e, diss_l_e , i_omp, tn, & 2590 bc_dirichlet_l .OR. bc_radiation_l, & 2591 bc_dirichlet_n .OR. bc_radiation_n, & 2592 bc_dirichlet_r .OR. bc_radiation_r, & 2571 CALL advec_s_ws( advc_flags_s, i, j, e, 'e', flux_s_e, diss_s_e, flux_l_e, & 2572 diss_l_e , i_omp, tn, & 2573 bc_dirichlet_l .OR. bc_radiation_l, & 2574 bc_dirichlet_n .OR. bc_radiation_n, & 2575 bc_dirichlet_r .OR. bc_radiation_r, & 2593 2576 bc_dirichlet_s .OR. bc_radiation_s ) 2594 2577 ELSE … … 2619 2602 ! 2620 2603 !-- Prognostic equation for TKE. 2621 !-- Eliminate negative TKE values, which can occur due to numerical 2622 !-- reasons in the course of the integration. In such cases the old 2623 !-- TKE value is reduced by 90%. 2604 !-- Eliminate negative TKE values, which can occur due to numerical reasons in the course of the 2605 !-- integration. In such cases the old TKE value is reduced by 90%. 2624 2606 DO k = nzb+1, nzt 2625 e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + & 2626 tsc(3) * te_m(k,j,i) ) & 2627 ) & 2628 * MERGE( 1.0_wp, 0.0_wp, & 2629 BTEST( wall_flags_total_0(k,j,i), 0 ) & 2630 ) 2607 e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) ) ) & 2608 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 2631 2609 IF ( e_p(k,j,i) <= 0.0_wp ) e_p(k,j,i) = 0.1_wp * e(k,j,i) 2632 2610 ENDDO … … 2639 2617 te_m(k,j,i) = tend(k,j,i) 2640 2618 ENDDO 2641 ELSEIF ( intermediate_timestep_count < & 2642 intermediate_timestep_count_max ) THEN 2619 ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) THEN 2643 2620 DO k = nzb+1, nzt 2644 te_m(k,j,i) = -9.5625_wp * tend(k,j,i) + & 2645 5.3125_wp * te_m(k,j,i) 2621 te_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * te_m(k,j,i) 2646 2622 ENDDO 2647 2623 ENDIF … … 2656 2632 !-- Tendency-terms for dissipation 2657 2633 tend(:,j,i) = 0.0_wp 2658 IF ( timestep_scheme(1:5) == 'runge' & 2659 .AND. .NOT. use_upstream_for_tke ) THEN 2634 IF ( timestep_scheme(1:5) == 'runge' .AND. .NOT. use_upstream_for_tke ) THEN 2660 2635 IF ( ws_scheme_sca ) THEN 2661 CALL advec_s_ws( advc_flags_s, & 2662 i, j, diss, 'diss', flux_s_diss, diss_s_diss, & 2663 flux_l_diss, diss_l_diss, i_omp, tn, & 2664 bc_dirichlet_l .OR. bc_radiation_l, & 2665 bc_dirichlet_n .OR. bc_radiation_n, & 2666 bc_dirichlet_r .OR. bc_radiation_r, & 2636 CALL advec_s_ws( advc_flags_s, i, j, diss, 'diss', flux_s_diss, diss_s_diss, & 2637 flux_l_diss, diss_l_diss, i_omp, tn, & 2638 bc_dirichlet_l .OR. bc_radiation_l, & 2639 bc_dirichlet_n .OR. bc_radiation_n, & 2640 bc_dirichlet_r .OR. bc_radiation_r, & 2667 2641 bc_dirichlet_s .OR. bc_radiation_s ) 2668 2642 ELSE … … 2686 2660 ! 2687 2661 !-- Prognostic equation for TKE dissipation 2688 !-- Eliminate negative dissipation values, which can occur due to 2689 !-- numerical reasons in the course of the integration. In such cases 2690 !-- the old dissipation value is reduced by 90%. 2662 !-- Eliminate negative dissipation values, which can occur due to numerical reasons in the course 2663 !-- of the integration. In such cases the old dissipation value is reduced by 90%. 2691 2664 DO k = nzb+1, nzt 2692 diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + & 2693 tsc(3) * tdiss_m(k,j,i) ) & 2694 ) & 2695 * MERGE( 1.0_wp, 0.0_wp, & 2696 BTEST( wall_flags_total_0(k,j,i), 0 )& 2697 ) 2665 diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + tsc(3) & 2666 * tdiss_m(k,j,i) ) ) & 2667 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 2698 2668 ENDDO 2699 2669 … … 2705 2675 tdiss_m(k,j,i) = tend(k,j,i) 2706 2676 ENDDO 2707 ELSEIF ( intermediate_timestep_count < & 2708 intermediate_timestep_count_max ) THEN 2677 ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) THEN 2709 2678 DO k = nzb+1, nzt 2710 tdiss_m(k,j,i) = -9.5625_wp * tend(k,j,i) + & 2711 5.3125_wp * tdiss_m(k,j,i) 2679 tdiss_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tdiss_m(k,j,i) 2712 2680 ENDDO 2713 2681 ENDIF 2714 2682 ENDIF 2715 2683 2716 ENDIF ! dissipation equation2684 ENDIF ! Dissipation equation 2717 2685 2718 2686 END SUBROUTINE tcm_prognostic_equations_ij 2719 2687 2720 2688 2721 !------------------------------------------------------------------------------ !2689 !--------------------------------------------------------------------------------------------------! 2722 2690 ! Description: 2723 2691 ! ------------ 2724 2692 !> Production terms (shear + buoyancy) of the TKE. 2725 2693 !> Vector-optimized version 2726 !> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is 2727 !> not considered well! 2728 !------------------------------------------------------------------------------! 2694 !> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is not considered well! 2695 !--------------------------------------------------------------------------------------------------! 2729 2696 SUBROUTINE production_e( diss_production ) 2730 2697 2731 USE arrays_3d, & 2732 ONLY: ddzw, dd2zu, drho_air_zw, q, ql, d_exner, exner 2733 2734 USE control_parameters, & 2735 ONLY: cloud_droplets, constant_flux_layer, neutral, & 2736 rho_reference, use_single_reference_value, use_surface_fluxes, & 2698 USE arrays_3d, & 2699 ONLY: ddzw, & 2700 dd2zu, & 2701 d_exner, & 2702 drho_air_zw, & 2703 exner, & 2704 q, & 2705 ql 2706 2707 2708 2709 USE control_parameters, & 2710 ONLY: cloud_droplets, & 2711 constant_flux_layer, & 2712 neutral, & 2713 rho_reference, & 2714 use_single_reference_value, & 2715 use_surface_fluxes, & 2737 2716 use_top_fluxes 2738 2717 2739 USE grid_variables, & 2740 ONLY: ddx, dx, ddy, dy 2741 2742 USE bulk_cloud_model_mod, & 2718 USE grid_variables, & 2719 ONLY: ddx, & 2720 dx, & 2721 ddy, & 2722 dy 2723 2724 USE bulk_cloud_model_mod, & 2743 2725 ONLY: bulk_cloud_model 2744 2726 2745 2727 IMPLICIT NONE 2746 2728 2747 LOGICAL :: diss_production2748 2749 INTEGER(iwp) :: i !< running index x-direction2750 INTEGER(iwp) :: j !< running index y-direction2751 INTEGER(iwp) :: k !< running index z-direction2752 INTEGER(iwp) :: l !< running index for different surface type orientation2753 INTEGER(iwp) :: m !< running index surface elements2754 INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position2755 INTEGER(iwp) :: surf_ s !< startindex of surface elements at given i-j position2756 INTEGER(iwp) :: flag_nr !< number of masking flag2729 LOGICAL :: diss_production !< 2730 2731 INTEGER(iwp) :: flag_nr !< number of masking flag 2732 INTEGER(iwp) :: i !< running index x-direction 2733 INTEGER(iwp) :: j !< running index y-direction 2734 INTEGER(iwp) :: k !< running index z-direction 2735 INTEGER(iwp) :: l !< running index for different surface type orientation 2736 INTEGER(iwp) :: m !< running index surface elements 2737 INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position 2738 INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position 2757 2739 2758 2740 REAL(wp) :: def !< ( du_i/dx_j + du_j/dx_i ) * du_i/dx_j … … 2761 2743 REAL(wp) :: k2 !< temporary factor 2762 2744 REAL(wp) :: km_neutral !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces 2745 REAL(wp) :: sign_dir !< sign of wall-tke flux, depending on wall orientation 2763 2746 REAL(wp) :: theta !< virtual potential temperature 2764 2747 REAL(wp) :: temp !< theta * Exner-function 2765 REAL(wp) :: sign_dir !< sign of wall-tke flux, depending on wall orientation2766 2748 REAL(wp) :: usvs !< momentum flux u"v" 2767 2749 REAL(wp) :: vsus !< momentum flux v"u" … … 2769 2751 REAL(wp) :: wsvs !< momentum flux w"v" 2770 2752 2771 REAL(wp), DIMENSION(nzb+1:nzt) :: dudx !< Gradient of u-component in x-direction2772 REAL(wp), DIMENSION(nzb+1:nzt) :: dudy !< Gradient of u-component in y-direction2773 REAL(wp), DIMENSION(nzb+1:nzt) :: dudz !< Gradient of u-component in z-direction2774 REAL(wp), DIMENSION(nzb+1:nzt) :: dvdx !< Gradient of v-component in x-direction2775 REAL(wp), DIMENSION(nzb+1:nzt) :: dvdy !< Gradient of v-component in y-direction2776 REAL(wp), DIMENSION(nzb+1:nzt) :: dvdz !< Gradient of v-component in z-direction2777 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdx !< Gradient of w-component in x-direction2778 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdy !< Gradient of w-component in y-direction2779 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdz !< Gradient of w-component in z-direction2753 REAL(wp), DIMENSION(nzb+1:nzt) :: dudx !< Gradient of u-component in x-direction 2754 REAL(wp), DIMENSION(nzb+1:nzt) :: dudy !< Gradient of u-component in y-direction 2755 REAL(wp), DIMENSION(nzb+1:nzt) :: dudz !< Gradient of u-component in z-direction 2756 REAL(wp), DIMENSION(nzb+1:nzt) :: dvdx !< Gradient of v-component in x-direction 2757 REAL(wp), DIMENSION(nzb+1:nzt) :: dvdy !< Gradient of v-component in y-direction 2758 REAL(wp), DIMENSION(nzb+1:nzt) :: dvdz !< Gradient of v-component in z-direction 2759 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdx !< Gradient of w-component in x-direction 2760 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdy !< Gradient of w-component in y-direction 2761 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdz !< Gradient of w-component in z-direction 2780 2762 REAL(wp), DIMENSION(nzb+1:nzt) :: tmp_flux !< temporary flux-array in z-direction 2781 2763 … … 2783 2765 2784 2766 ! 2785 !-- Calculate TKE production by shear. Calculate gradients at all grid 2786 !-- points first, gradients at surface-bounded grid points will be 2787 !-- overwritten further below. 2767 !-- Calculate TKE production by shear. Calculate gradients at all grid points first, gradients at 2768 !-- surface-bounded grid points will be overwritten further below. 2788 2769 !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, l) & 2789 2770 !$ACC PRIVATE(surf_s, surf_e) & … … 2800 2781 2801 2782 dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx 2802 dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 2803 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 2804 dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 2805 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 2806 2807 dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 2808 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 2783 dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 2784 dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 2785 2786 dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 2809 2787 dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy 2810 dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 2811 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 2812 2813 dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 2814 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 2815 dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 2816 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 2788 dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 2789 2790 dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 2791 dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 2817 2792 dwdz(k) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 2818 2793 … … 2832 2807 !-- 'bottom and wall: use u_0,v_0 and wall functions' 2833 2808 ! 2834 !-- Compute gradients at north- and south-facing surfaces. 2835 !-- First, for default surfaces, then for urban surfaces. 2836 !-- Note, so far no natural vertical surfaces implemented 2809 !-- Compute gradients at north- and south-facing surfaces. First, for default surfaces, 2810 !-- then for urban surfaces. Note, so far no natural vertical surfaces implemented 2837 2811 DO l = 0, 1 2838 2812 surf_s = surf_def_v(l)%start_index(j,i) … … 2844 2818 wsvs = surf_def_v(l)%mom_flux_tke(1,m) 2845 2819 2846 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 2847 * 0.5_wp * dy 2820 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp * 0.5_wp * dy 2848 2821 ! 2849 2822 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 2850 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 2851 BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) ) 2823 sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) ) 2852 2824 dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) 2853 2825 dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) … … 2863 2835 wsvs = surf_lsm_v(l)%mom_flux_tke(1,m) 2864 2836 2865 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 2866 * 0.5_wp * dy 2837 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp * 0.5_wp * dy 2867 2838 ! 2868 2839 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 2869 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 2870 BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) ) 2840 sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) ) 2871 2841 dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) 2872 2842 dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) … … 2882 2852 wsvs = surf_usm_v(l)%mom_flux_tke(1,m) 2883 2853 2884 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 2885 * 0.5_wp * dy 2854 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp * 0.5_wp * dy 2886 2855 ! 2887 2856 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 2888 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 2889 BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) ) 2857 sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) ) 2890 2858 dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) 2891 2859 dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) … … 2903 2871 wsus = surf_def_v(l)%mom_flux_tke(1,m) 2904 2872 2905 km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & 2906 * 0.5_wp * dx 2873 km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp * 0.5_wp * dx 2907 2874 ! 2908 2875 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 2909 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 2910 BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) ) 2876 sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) ) 2911 2877 dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) 2912 2878 dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) … … 2922 2888 wsus = surf_lsm_v(l)%mom_flux_tke(1,m) 2923 2889 2924 km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & 2925 * 0.5_wp * dx 2890 km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp * 0.5_wp * dx 2926 2891 ! 2927 2892 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 2928 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 2929 BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) ) 2893 sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) ) 2930 2894 dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) 2931 2895 dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) … … 2941 2905 wsus = surf_usm_v(l)%mom_flux_tke(1,m) 2942 2906 2943 km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & 2944 * 0.5_wp * dx 2907 km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp * 0.5_wp * dx 2945 2908 ! 2946 2909 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 2947 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 2948 BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) ) 2910 sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) ) 2949 2911 dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) 2950 2912 dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) … … 2959 2921 k = surf_def_h(0)%k(m) 2960 2922 ! 2961 !-- Please note, actually, an interpolation of u_0 and v_0 2962 !-- onto the grid center would be required. However, this 2963 !-- would require several data transfers between 2D-grid and 2964 !-- wall type. The effect of this missing interpolation is 2965 !-- negligible. (See also production_e_init). 2923 !-- Please note, actually, an interpolation of u_0 and v_0 onto the grid center would be 2924 !-- required. However, this would require several data transfers between 2D-grid and 2925 !-- wall type. The effect of this missing interpolation is negligible. 2926 !-- (See also production_e_init). 2966 2927 dudz(k) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k) 2967 2928 dvdz(k) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k) … … 2993 2954 ENDDO 2994 2955 ! 2995 !-- Compute gradients at downward-facing walls, only for 2996 !-- non-natural default surfaces 2956 !-- Compute gradients at downward-facing walls, only for non-natural default surfaces 2997 2957 surf_s = surf_def_h(1)%start_index(j,i) 2998 2958 surf_e = surf_def_h(1)%end_index(j,i) … … 3013 2973 DO k = nzb+1, nzt 3014 2974 3015 def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) + & 3016 dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 + & 3017 dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 + & 3018 2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + & 3019 dwdy(k)*dvdz(k) ) 2975 def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) + dudy(k)**2 + dvdx(k)**2 + & 2976 dwdx(k)**2 + dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 + & 2977 2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) ) 3020 2978 3021 2979 IF ( def < 0.0_wp ) def = 0.0_wp … … 3031 2989 3032 2990 !-- RANS mode: Compute tendency for dissipation-rate-production from shear 3033 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag * &3034 diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_12991 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag * diss(k,j,i) / & 2992 ( e(k,j,i) + 1.0E-20_wp ) * c_1 3035 2993 3036 2994 ENDIF … … 3050 3008 IF ( ocean_mode ) THEN 3051 3009 ! 3052 !-- So far in the ocean no special treatment of density flux 3053 !-- in the bottom and top surfacelayer3010 !-- So far in the ocean no special treatment of density flux in the bottom and top surface 3011 !-- layer 3054 3012 DO i = nxl, nxr 3055 3013 DO j = nys, nyn … … 3059 3017 ENDDO 3060 3018 ! 3061 !-- Treatment of near-surface grid points, at up- and down- 3062 !-- ward facing surfaces 3019 !-- Treatment of near-surface grid points, at up- and down-ward facing surfaces 3063 3020 IF ( use_surface_fluxes ) THEN 3064 3021 DO l = 0, 1 … … 3085 3042 !-- Compute tendency for TKE-production from shear 3086 3043 DO k = nzb+1, nzt 3087 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),0) )3088 tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &3089 MERGE( rho_reference, prho(k,j,i), &3044 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 3045 tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / & 3046 MERGE( rho_reference, prho(k,j,i), & 3090 3047 use_single_reference_value ) ) 3091 3048 ENDDO … … 3095 3052 !-- RANS mode: Compute tendency for dissipation-rate-production from shear 3096 3053 DO k = nzb+1, nzt 3097 flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) ) 3098 tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / & 3099 MERGE( rho_reference, prho(k,j,i), & 3100 use_single_reference_value ) ) * & 3101 diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * & 3102 c_3 3054 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 3055 tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / & 3056 MERGE( rho_reference, prho(k,j,i), & 3057 use_single_reference_value ) ) * & 3058 diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_3 3103 3059 ENDDO 3104 3060 … … 3173 3129 !$ACC LOOP PRIVATE(k, flag) 3174 3130 DO k = nzb+1, nzt 3175 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),0) )3176 tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &3177 MERGE( pt_reference, pt(k,j,i), &3131 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 3132 tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / & 3133 MERGE( pt_reference, pt(k,j,i), & 3178 3134 use_single_reference_value ) ) 3179 3135 ENDDO … … 3183 3139 !-- RANS mode: Compute tendency for dissipation-rate-production from shear 3184 3140 DO k = nzb+1, nzt 3185 flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) ) 3186 tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / & 3187 MERGE( pt_reference, pt(k,j,i), & 3188 use_single_reference_value ) ) * & 3189 diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * & 3190 c_3 3141 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 3142 tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / & 3143 MERGE( pt_reference, pt(k,j,i), & 3144 use_single_reference_value ) ) * & 3145 diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_3 3191 3146 ENDDO 3192 3147 … … 3196 3151 ENDDO 3197 3152 3198 ENDIF ! from IF ( .NOT. ocean_mode )3199 3200 ELSE ! or IF ( humidity ) THEN3153 ENDIF ! From IF ( .NOT. ocean_mode ) 3154 3155 ELSE ! Or IF ( humidity ) THEN 3201 3156 3202 3157 DO i = nxl, nxr … … 3205 3160 DO k = nzb+1, nzt 3206 3161 3207 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN3162 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN 3208 3163 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 3209 3164 k2 = 0.61_wp * pt(k,j,i) 3210 tmp_flux(k) = -1.0_wp * kh(k,j,i) * & 3211 ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) + & 3212 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & 3213 ) * dd2zu(k) 3165 tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) + & 3166 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) ) & 3167 * dd2zu(k) 3214 3168 ELSE IF ( bulk_cloud_model ) THEN 3215 3169 IF ( ql(k,j,i) == 0.0_wp ) THEN … … 3219 3173 theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i) 3220 3174 temp = theta * exner(k) 3221 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 3222 ( q(k,j,i) - ql(k,j,i) ) * & 3223 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / & 3224 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * & 3175 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * ( q(k,j,i) - ql(k,j,i) ) * & 3176 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / & 3177 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * & 3225 3178 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 3226 3179 k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp ) 3227 3180 ENDIF 3228 tmp_flux(k) = -1.0_wp * kh(k,j,i) * & 3229 ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) + & 3230 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & 3231 ) * dd2zu(k) 3181 tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) + & 3182 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) ) & 3183 * dd2zu(k) 3232 3184 ELSE IF ( cloud_droplets ) THEN 3233 3185 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 3234 3186 k2 = 0.61_wp * pt(k,j,i) 3235 tmp_flux(k) = -1.0_wp * kh(k,j,i) * & 3236 ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) + & 3237 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - & 3238 pt(k,j,i) * ( ql(k+1,j,i) - & 3239 ql(k-1,j,i) ) ) * dd2zu(k) 3187 tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) + & 3188 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - & 3189 pt(k,j,i) * ( ql(k+1,j,i) - & 3190 ql(k-1,j,i) ) ) * dd2zu(k) 3240 3191 ENDIF 3241 3192 … … 3251 3202 k = surf_def_h(l)%k(m) 3252 3203 3253 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets )THEN3204 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN 3254 3205 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 3255 3206 k2 = 0.61_wp * pt(k,j,i) … … 3261 3212 theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i) 3262 3213 temp = theta * exner(k) 3263 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 3264 ( q(k,j,i) - ql(k,j,i) ) * & 3265 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / & 3266 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * & 3267 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 3214 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * ( q(k,j,i) - ql(k,j,i) ) * & 3215 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / & 3216 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * & 3217 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 3268 3218 k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp ) 3269 3219 ENDIF … … 3273 3223 ENDIF 3274 3224 3275 tmp_flux(k) = ( k1 * surf_def_h(l)%shf(m) + & 3276 k2 * surf_def_h(l)%qsws(m) & 3277 ) * drho_air_zw(k-1) 3225 tmp_flux(k) = ( k1 * surf_def_h(l)%shf(m) + k2 * surf_def_h(l)%qsws(m) ) & 3226 * drho_air_zw(k-1) 3278 3227 ENDDO 3279 3228 ENDDO … … 3285 3234 k = surf_lsm_h%k(m) 3286 3235 3287 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets )THEN3236 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN 3288 3237 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 3289 3238 k2 = 0.61_wp * pt(k,j,i) … … 3295 3244 theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i) 3296 3245 temp = theta * exner(k) 3297 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 3298 ( q(k,j,i) - ql(k,j,i) ) * & 3299 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / & 3300 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * & 3246 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * ( q(k,j,i) - ql(k,j,i) ) * & 3247 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / & 3248 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * & 3301 3249 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 3302 3250 k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp ) … … 3307 3255 ENDIF 3308 3256 3309 tmp_flux(k) = ( k1 * surf_lsm_h%shf(m) + & 3310 k2 * surf_lsm_h%qsws(m) & 3311 ) * drho_air_zw(k-1) 3257 tmp_flux(k) = ( k1 * surf_lsm_h%shf(m) + k2 * surf_lsm_h%qsws(m) ) & 3258 * drho_air_zw(k-1) 3312 3259 ENDDO 3313 3260 ! … … 3318 3265 k = surf_usm_h%k(m) 3319 3266 3320 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets )THEN3267 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN 3321 3268 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 3322 3269 k2 = 0.61_wp * pt(k,j,i) … … 3328 3275 theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i) 3329 3276 temp = theta * exner(k) 3330 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 3331 ( q(k,j,i) - ql(k,j,i) ) * & 3332 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / & 3333 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * & 3277 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * ( q(k,j,i) - ql(k,j,i) ) * & 3278 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / & 3279 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * & 3334 3280 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 3335 3281 k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp ) … … 3340 3286 ENDIF 3341 3287 3342 tmp_flux(k) = ( k1 * surf_usm_h%shf(m) + & 3343 k2 * surf_usm_h%qsws(m) & 3344 ) * drho_air_zw(k-1) 3288 tmp_flux(k) = ( k1 * surf_usm_h%shf(m) + k2 * surf_usm_h%qsws(m) ) & 3289 * drho_air_zw(k-1) 3345 3290 ENDDO 3346 3291 3347 ENDIF ! from IF ( use_surface_fluxes ) THEN3292 ENDIF ! From IF ( use_surface_fluxes ) THEN 3348 3293 3349 3294 IF ( use_top_fluxes ) THEN … … 3354 3299 k = surf_def_h(2)%k(m) 3355 3300 3356 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets )THEN3301 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN 3357 3302 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 3358 3303 k2 = 0.61_wp * pt(k,j,i) … … 3364 3309 theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i) 3365 3310 temp = theta * exner(k) 3366 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 3367 ( q(k,j,i) - ql(k,j,i) ) * & 3368 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / & 3369 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * & 3311 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * ( q(k,j,i) - ql(k,j,i) ) * & 3312 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / & 3313 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * & 3370 3314 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 3371 3315 k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp ) … … 3376 3320 ENDIF 3377 3321 3378 tmp_flux(k) = ( k1 * surf_def_h(2)%shf(m) + & 3379 k2 * surf_def_h(2)%qsws(m) & 3380 ) * drho_air_zw(k) 3322 tmp_flux(k) = ( k1 * surf_def_h(2)%shf(m) + k2 * surf_def_h(2)%qsws(m) ) & 3323 * drho_air_zw(k) 3381 3324 3382 3325 ENDDO 3383 3326 3384 ENDIF ! from IF ( use_top_fluxes ) THEN3327 ENDIF ! From IF ( use_top_fluxes ) THEN 3385 3328 3386 3329 IF ( .NOT. diss_production ) THEN … … 3388 3331 !-- Compute tendency for TKE-production from shear 3389 3332 DO k = nzb+1, nzt 3390 flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) ) 3391 tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / & 3392 MERGE( vpt_reference, vpt(k,j,i), & 3393 use_single_reference_value ) ) 3333 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 3334 tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / & 3335 MERGE( vpt_reference, vpt(k,j,i), use_single_reference_value ) ) 3394 3336 ENDDO 3395 3337 … … 3398 3340 !-- RANS mode: Compute tendency for dissipation-rate-production from shear 3399 3341 DO k = nzb+1, nzt 3400 flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) ) 3401 tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / & 3402 MERGE( vpt_reference, vpt(k,j,i), & 3403 use_single_reference_value ) ) * & 3404 diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * & 3405 c_3 3342 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 3343 tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * & 3344 ( g / MERGE( vpt_reference, vpt(k,j,i), & 3345 use_single_reference_value ) ) & 3346 * diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_3 3406 3347 ENDDO 3407 3348 … … 3418 3359 3419 3360 3420 !------------------------------------------------------------------------------ !3361 !--------------------------------------------------------------------------------------------------! 3421 3362 ! Description: 3422 3363 ! ------------ 3423 3364 !> Production terms (shear + buoyancy) of the TKE. 3424 3365 !> Cache-optimized version 3425 !> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is 3426 !> not considered well! 3427 !------------------------------------------------------------------------------! 3366 !> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is not considered well! 3367 !--------------------------------------------------------------------------------------------------! 3428 3368 SUBROUTINE production_e_ij( i, j, diss_production ) 3429 3369 3430 USE arrays_3d, & 3431 ONLY: ddzw, dd2zu, drho_air_zw, q, ql, d_exner, exner 3432 3433 USE control_parameters, & 3434 ONLY: cloud_droplets, constant_flux_layer, neutral, & 3435 rho_reference, use_single_reference_value, use_surface_fluxes, & 3370 USE arrays_3d, & 3371 ONLY: ddzw, & 3372 dd2zu, & 3373 drho_air_zw, & 3374 d_exner, & 3375 exner, & 3376 q, & 3377 ql 3378 3379 3380 3381 USE control_parameters, & 3382 ONLY: cloud_droplets, & 3383 constant_flux_layer, & 3384 neutral, & 3385 rho_reference, & 3386 use_single_reference_value, & 3387 use_surface_fluxes, & 3436 3388 use_top_fluxes 3437 3389 3438 USE grid_variables, & 3439 ONLY: ddx, dx, ddy, dy 3440 3441 USE bulk_cloud_model_mod, & 3390 USE grid_variables, & 3391 ONLY: ddx, & 3392 dx, & 3393 ddy, & 3394 dy 3395 3396 USE bulk_cloud_model_mod, & 3442 3397 ONLY: bulk_cloud_model 3443 3398 3444 3399 IMPLICIT NONE 3445 3400 3446 LOGICAL :: diss_production3447 3448 INTEGER(iwp) :: i !< running index x-direction3449 INTEGER(iwp) :: j !< running index y-direction3450 INTEGER(iwp) :: k !< running index z-direction3451 INTEGER(iwp) :: l !< running index for different surface type orientation3452 INTEGER(iwp) :: m !< running index surface elements3453 INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position3454 INTEGER(iwp) :: surf_ s !< startindex of surface elements at given i-j position3455 INTEGER(iwp) :: flag_nr !< number of masking flag3401 LOGICAL :: diss_production !< 3402 3403 INTEGER(iwp) :: flag_nr !< number of masking flag 3404 INTEGER(iwp) :: i !< running index x-direction 3405 INTEGER(iwp) :: j !< running index y-direction 3406 INTEGER(iwp) :: k !< running index z-direction 3407 INTEGER(iwp) :: l !< running index for different surface type orientation 3408 INTEGER(iwp) :: m !< running index surface elements 3409 INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position 3410 INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position 3456 3411 3457 3412 REAL(wp) :: def !< ( du_i/dx_j + du_j/dx_i ) * du_i/dx_j … … 3460 3415 REAL(wp) :: k2 !< temporary factor 3461 3416 REAL(wp) :: km_neutral !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces 3417 REAL(wp) :: sign_dir !< sign of wall-tke flux, depending on wall orientation 3462 3418 REAL(wp) :: theta !< virtual potential temperature 3463 3419 REAL(wp) :: temp !< theta * Exner-function 3464 REAL(wp) :: sign_dir !< sign of wall-tke flux, depending on wall orientation3465 3420 REAL(wp) :: usvs !< momentum flux u"v" 3466 3421 REAL(wp) :: vsus !< momentum flux v"u" … … 3468 3423 REAL(wp) :: wsvs !< momentum flux w"v" 3469 3424 3470 REAL(wp), DIMENSION(nzb+1:nzt) :: dudx !< Gradient of u-component in x-direction3471 REAL(wp), DIMENSION(nzb+1:nzt) :: dudy !< Gradient of u-component in y-direction3472 REAL(wp), DIMENSION(nzb+1:nzt) :: dudz !< Gradient of u-component in z-direction3473 REAL(wp), DIMENSION(nzb+1:nzt) :: dvdx !< Gradient of v-component in x-direction3474 REAL(wp), DIMENSION(nzb+1:nzt) :: dvdy !< Gradient of v-component in y-direction3475 REAL(wp), DIMENSION(nzb+1:nzt) :: dvdz !< Gradient of v-component in z-direction3476 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdx !< Gradient of w-component in x-direction3477 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdy !< Gradient of w-component in y-direction3478 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdz !< Gradient of w-component in z-direction3425 REAL(wp), DIMENSION(nzb+1:nzt) :: dudx !< Gradient of u-component in x-direction 3426 REAL(wp), DIMENSION(nzb+1:nzt) :: dudy !< Gradient of u-component in y-direction 3427 REAL(wp), DIMENSION(nzb+1:nzt) :: dudz !< Gradient of u-component in z-direction 3428 REAL(wp), DIMENSION(nzb+1:nzt) :: dvdx !< Gradient of v-component in x-direction 3429 REAL(wp), DIMENSION(nzb+1:nzt) :: dvdy !< Gradient of v-component in y-direction 3430 REAL(wp), DIMENSION(nzb+1:nzt) :: dvdz !< Gradient of v-component in z-direction 3431 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdx !< Gradient of w-component in x-direction 3432 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdy !< Gradient of w-component in y-direction 3433 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdz !< Gradient of w-component in z-direction 3479 3434 REAL(wp), DIMENSION(nzb+1:nzt) :: tmp_flux !< temporary flux-array in z-direction 3480 3435 … … 3482 3437 3483 3438 ! 3484 !-- Calculate TKE production by shear. Calculate gradients at all grid 3485 !-- points first, gradients at surface-bounded grid points will be 3486 !-- overwritten further below. 3439 !-- Calculate TKE production by shear. Calculate gradients at all grid points first, gradients at 3440 !-- surface-bounded grid points will be overwritten further below. 3487 3441 DO k = nzb+1, nzt 3488 3442 3489 3443 dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx 3490 dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 3491 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 3492 dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 3493 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 3494 3495 dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 3496 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 3444 dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 3445 dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 3446 3447 dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 3497 3448 dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy 3498 dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 3499 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 3500 3501 dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 3502 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 3503 dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 3504 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 3449 dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 3450 3451 dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 3452 dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 3505 3453 dwdz(k) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 3506 3454 … … 3517 3465 !-- 'bottom and wall: use u_0,v_0 and wall functions' 3518 3466 ! 3519 !-- Compute gradients at north- and south-facing surfaces. 3520 !-- First, for default surfaces, then for urban surfaces. 3521 !-- Note, so far no natural vertical surfaces implemented 3467 !-- Compute gradients at north- and south-facing surfaces. First, for default surfaces, then for 3468 !-- urban surfaces. Note, so far no natural vertical surfaces implemented 3522 3469 DO l = 0, 1 3523 3470 surf_s = surf_def_v(l)%start_index(j,i) … … 3528 3475 wsvs = surf_def_v(l)%mom_flux_tke(1,m) 3529 3476 3530 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 3531 * 0.5_wp * dy 3477 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp * 0.5_wp * dy 3532 3478 ! 3533 3479 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 3534 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 3535 BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) ) 3480 sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) ) 3536 3481 dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) 3537 3482 dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) … … 3546 3491 wsvs = surf_lsm_v(l)%mom_flux_tke(1,m) 3547 3492 3548 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 3549 * 0.5_wp * dy 3493 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp * 0.5_wp * dy 3550 3494 ! 3551 3495 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 3552 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 3553 BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) ) 3496 sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) ) 3554 3497 dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) 3555 3498 dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) … … 3564 3507 wsvs = surf_usm_v(l)%mom_flux_tke(1,m) 3565 3508 3566 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 3567 * 0.5_wp * dy 3509 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp * 0.5_wp * dy 3568 3510 ! 3569 3511 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 3570 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 3571 BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) ) 3512 sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) ) 3572 3513 dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) 3573 3514 dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) … … 3584 3525 wsus = surf_def_v(l)%mom_flux_tke(1,m) 3585 3526 3586 km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & 3587 * 0.5_wp * dx 3527 km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp * 0.5_wp * dx 3588 3528 ! 3589 3529 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 3590 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 3591 BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) ) 3530 sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) ) 3592 3531 dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) 3593 3532 dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) … … 3602 3541 wsus = surf_lsm_v(l)%mom_flux_tke(1,m) 3603 3542 3604 km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & 3605 * 0.5_wp * dx 3543 km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp * 0.5_wp * dx 3606 3544 ! 3607 3545 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 3608 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 3609 BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) ) 3546 sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) ) 3610 3547 dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) 3611 3548 dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) … … 3620 3557 wsus = surf_usm_v(l)%mom_flux_tke(1,m) 3621 3558 3622 km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & 3623 * 0.5_wp * dx 3559 km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp * 0.5_wp * dx 3624 3560 ! 3625 3561 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 3626 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 3627 BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) ) 3562 sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) ) 3628 3563 dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) 3629 3564 dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) … … 3637 3572 k = surf_def_h(0)%k(m) 3638 3573 ! 3639 !-- Please note, actually, an interpolation of u_0 and v_0 3640 !-- onto the grid center would be required. However, this 3641 !-- would require several data transfers between 2D-grid and 3642 !-- wall type. The effect of this missing interpolation is 3643 !-- negligible. (See also production_e_init). 3574 !-- Please note, actually, an interpolation of u_0 and v_0 onto the grid center would be 3575 !-- required. However, this would require several data transfers between 2D-grid and wall 3576 !-- type. The effect of this missing interpolation is negligible. (See also production_e_init). 3644 3577 dudz(k) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k) 3645 3578 dvdz(k) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k) … … 3669 3602 ENDDO 3670 3603 ! 3671 !-- Compute gradients at downward-facing walls, only for 3672 !-- non-natural default surfaces 3604 !-- Compute gradients at downward-facing walls, only for non-natural default surfaces 3673 3605 surf_s = surf_def_h(1)%start_index(j,i) 3674 3606 surf_e = surf_def_h(1)%end_index(j,i) … … 3685 3617 DO k = nzb+1, nzt 3686 3618 3687 def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) + & 3688 dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 + & 3689 dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 + & 3690 2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + & 3691 dwdy(k)*dvdz(k) ) 3619 def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) + & 3620 dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 + & 3621 dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 + & 3622 2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) ) 3692 3623 3693 3624 IF ( def < 0.0_wp ) def = 0.0_wp … … 3703 3634 3704 3635 !-- RANS mode: Compute tendency for dissipation-rate-production from shear