Changeset 1517 for palm/trunk/SOURCE/advec_s_bc.f90
- Timestamp:
- Jan 7, 2015 7:12:25 PM (9 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/advec_s_bc.f90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_s_bc.f90
r1375 r1517 1 SUBROUTINE advec_s_bc( sk, sk_char ) 1 MODULE advec_s_bc_mod 2 2 3 3 !--------------------------------------------------------------------------------! … … 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! interface added to advec_s_bc 23 23 ! 24 24 ! Former revisions: … … 76 76 !------------------------------------------------------------------------------! 77 77 78 USE advection, & 79 ONLY: aex, bex, dex, eex 80 81 USE arrays_3d, & 82 ONLY: d, ddzw, dzu, dzw, tend, u, v, w 83 84 USE control_parameters, & 85 ONLY: dt_3d, bc_pt_t_val, bc_q_t_val, ibc_pt_b, ibc_pt_t, ibc_q_t, & 86 message_string, pt_slope_offset, sloping_surface, u_gtrans, & 87 v_gtrans 88 89 USE cpulog, & 90 ONLY: cpu_log, log_point_s 91 92 USE grid_variables, & 93 ONLY: ddx, ddy 94 95 USE indices, & 96 ONLY: nx, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt 97 98 USE kinds 99 100 USE pegrid 101 102 USE statistics, & 103 ONLY: rmask, statistic_regions, sums_wsts_bc_l 104 105 106 IMPLICIT NONE 107 108 CHARACTER (LEN=*) :: sk_char !: 109 110 INTEGER(iwp) :: i !: 111 INTEGER(iwp) :: ix !: 112 INTEGER(iwp) :: j !: 113 INTEGER(iwp) :: k !: 114 INTEGER(iwp) :: ngp !: 115 INTEGER(iwp) :: sr !: 116 INTEGER(iwp) :: type_xz_2 !: 117 118 REAL(wp) :: cim !: 119 REAL(wp) :: cimf !: 120 REAL(wp) :: cip !: 121 REAL(wp) :: cipf !: 122 REAL(wp) :: d_new !: 123 REAL(wp) :: denomi !: denominator 124 REAL(wp) :: fminus !: 125 REAL(wp) :: fplus !: 126 REAL(wp) :: f2 !: 127 REAL(wp) :: f4 !: 128 REAL(wp) :: f8 !: 129 REAL(wp) :: f12 !: 130 REAL(wp) :: f24 !: 131 REAL(wp) :: f48 !: 132 REAL(wp) :: f1920 !: 133 REAL(wp) :: im !: 134 REAL(wp) :: ip !: 135 REAL(wp) :: m2 !: 136 REAL(wp) :: m3 !: 137 REAL(wp) :: numera !: numerator 138 REAL(wp) :: snenn !: 139 REAL(wp) :: sterm !: 140 REAL(wp) :: tendcy !: 141 REAL(wp) :: t1 !: 142 REAL(wp) :: t2 !: 143 144 REAL(wp) :: fmax(2) !: 145 REAL(wp) :: fmax_l(2) !: 146 78 PRIVATE 79 PUBLIC advec_s_bc 80 81 INTERFACE advec_s_bc 82 MODULE PROCEDURE advec_s_bc 83 END INTERFACE advec_s_bc 84 85 CONTAINS 86 87 SUBROUTINE advec_s_bc( sk, sk_char ) 88 89 USE advection, & 90 ONLY: aex, bex, dex, eex 91 92 USE arrays_3d, & 93 ONLY: d, ddzw, dzu, dzw, tend, u, v, w 94 95 USE control_parameters, & 96 ONLY: dt_3d, bc_pt_t_val, bc_q_t_val, ibc_pt_b, ibc_pt_t, ibc_q_t, & 97 message_string, pt_slope_offset, sloping_surface, u_gtrans, & 98 v_gtrans 99 100 USE cpulog, & 101 ONLY: cpu_log, log_point_s 102 103 USE grid_variables, & 104 ONLY: ddx, ddy 105 106 USE indices, & 107 ONLY: nx, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt 108 109 USE kinds 110 111 USE pegrid 112 113 USE statistics, & 114 ONLY: rmask, statistic_regions, sums_wsts_bc_l 115 116 117 IMPLICIT NONE 118 119 CHARACTER (LEN=*) :: sk_char !: 120 121 INTEGER(iwp) :: i !: 122 INTEGER(iwp) :: ix !: 123 INTEGER(iwp) :: j !: 124 INTEGER(iwp) :: k !: 125 INTEGER(iwp) :: ngp !: 126 INTEGER(iwp) :: sr !: 127 INTEGER(iwp) :: type_xz_2 !: 128 129 REAL(wp) :: cim !: 130 REAL(wp) :: cimf !: 131 REAL(wp) :: cip !: 132 REAL(wp) :: cipf !: 133 REAL(wp) :: d_new !: 134 REAL(wp) :: denomi !: denominator 135 REAL(wp) :: fminus !: 136 REAL(wp) :: fplus !: 137 REAL(wp) :: f2 !: 138 REAL(wp) :: f4 !: 139 REAL(wp) :: f8 !: 140 REAL(wp) :: f12 !: 141 REAL(wp) :: f24 !: 142 REAL(wp) :: f48 !: 143 REAL(wp) :: f1920 !: 144 REAL(wp) :: im !: 145 REAL(wp) :: ip !: 146 REAL(wp) :: m2 !: 147 REAL(wp) :: m3 !: 148 REAL(wp) :: numera !: numerator 149 REAL(wp) :: snenn !: 150 REAL(wp) :: sterm !: 151 REAL(wp) :: tendcy !: 152 REAL(wp) :: t1 !: 153 REAL(wp) :: t2 !: 154 155 REAL(wp) :: fmax(2) !: 156 REAL(wp) :: fmax_l(2) !: 157 147 158 #if defined( __nopointer ) 148 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !:159 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !: 149 160 #else 150 REAL(wp), DIMENSION(:,:,:), POINTER :: sk161 REAL(wp), DIMENSION(:,:,:), POINTER :: sk 151 162 #endif 152 163 153 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: a0 !:154 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: a1 !:155 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: a12 !:156 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: a2 !:157 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: a22 !:158 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: immb !:159 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: imme !:160 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: impb !:161 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: impe !:162 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ipmb !:163 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ipme !:164 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ippb !:165 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ippe !:166 167 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: sk_p !:164 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: a0 !: 165 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: a1 !: 166 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: a12 !: 167 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: a2 !: 168 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: a22 !: 169 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: immb !: 170 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: imme !: 171 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: impb !: 172 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: impe !: 173 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ipmb !: 174 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ipme !: 175 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ippb !: 176 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ippe !: 177 178 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: sk_p !: 168 179 169 180 #if defined( __nec ) 170 REAL(sp) :: m1n, m1z !Wichtig: Division !:171 REAL(sp), DIMENSION(:,:), ALLOCATABLE :: m1, sw !:181 REAL(sp) :: m1n, m1z !Wichtig: Division !: 182 REAL(sp), DIMENSION(:,:), ALLOCATABLE :: m1, sw !: 172 183 #else 173 REAL(wp) :: m1n, m1z174 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: m1, sw184 REAL(wp) :: m1n, m1z 185 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: m1, sw 175 186 #endif 176 187 177 188 178 189 ! 179 !-- Array sk_p requires 2 extra elements for each dimension 180 ALLOCATE( sk_p(nzb-2:nzt+3,nys-3:nyn+3,nxl-3:nxr+3) ) 181 sk_p = 0.0_wp 182 183 ! 184 !-- Assign reciprocal values in order to avoid divisions later 185 f2 = 0.5_wp 186 f4 = 0.25_wp 187 f8 = 0.125_wp 188 f12 = 0.8333333333333333E-01_wp 189 f24 = 0.4166666666666666E-01_wp 190 f48 = 0.2083333333333333E-01_wp 191 f1920 = 0.5208333333333333E-03_wp 192 193 ! 194 !-- Advection in x-direction: 195 196 ! 197 !-- Save the quantity to be advected in a local array 198 !-- add an enlarged boundary in x-direction 199 DO i = nxl-1, nxr+1 190 !-- Array sk_p requires 2 extra elements for each dimension 191 ALLOCATE( sk_p(nzb-2:nzt+3,nys-3:nyn+3,nxl-3:nxr+3) ) 192 sk_p = 0.0_wp 193 194 ! 195 !-- Assign reciprocal values in order to avoid divisions later 196 f2 = 0.5_wp 197 f4 = 0.25_wp 198 f8 = 0.125_wp 199 f12 = 0.8333333333333333E-01_wp 200 f24 = 0.4166666666666666E-01_wp 201 f48 = 0.2083333333333333E-01_wp 202 f1920 = 0.5208333333333333E-03_wp 203 204 ! 205 !-- Advection in x-direction: 206 207 ! 208 !-- Save the quantity to be advected in a local array 209 !-- add an enlarged boundary in x-direction 210 DO i = nxl-1, nxr+1 211 DO j = nys, nyn 212 DO k = nzb, nzt+1 213 sk_p(k,j,i) = sk(k,j,i) 214 ENDDO 215 ENDDO 216 ENDDO 217 #if defined( __parallel ) 218 ngp = 2 * ( nzt - nzb + 6 ) * ( nyn - nys + 7 ) 219 CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'start' ) 220 ! 221 !-- Send left boundary, receive right boundary 222 CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxl+1), ngp, MPI_REAL, pleft, 0, & 223 sk_p(nzb-2,nys-3,nxr+2), ngp, MPI_REAL, pright, 0, & 224 comm2d, status, ierr ) 225 ! 226 !-- Send right boundary, receive left boundary 227 CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxr-2), ngp, MPI_REAL, pright, 1, & 228 sk_p(nzb-2,nys-3,nxl-3), ngp, MPI_REAL, pleft, 1, & 229 comm2d, status, ierr ) 230 CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' ) 231 #else 232 233 ! 234 !-- Cyclic boundary conditions 235 sk_p(:,nys:nyn,nxl-3) = sk_p(:,nys:nyn,nxr-2) 236 sk_p(:,nys:nyn,nxl-2) = sk_p(:,nys:nyn,nxr-1) 237 sk_p(:,nys:nyn,nxr+2) = sk_p(:,nys:nyn,nxl+1) 238 sk_p(:,nys:nyn,nxr+3) = sk_p(:,nys:nyn,nxl+2) 239 #endif 240 241 ! 242 !-- In case of a sloping surface, the additional gridpoints in x-direction 243 !-- of the temperature field at the left and right boundary of the total 244 !-- domain must be adjusted by the temperature difference between this distance 245 IF ( sloping_surface .AND. sk_char == 'pt' ) THEN 246 IF ( nxl == 0 ) THEN 247 sk_p(:,nys:nyn,nxl-3) = sk_p(:,nys:nyn,nxl-3) - pt_slope_offset 248 sk_p(:,nys:nyn,nxl-2) = sk_p(:,nys:nyn,nxl-2) - pt_slope_offset 249 ENDIF 250 IF ( nxr == nx ) THEN 251 sk_p(:,nys:nyn,nxr+2) = sk_p(:,nys:nyn,nxr+2) + pt_slope_offset 252 sk_p(:,nys:nyn,nxr+3) = sk_p(:,nys:nyn,nxr+3) + pt_slope_offset 253 ENDIF 254 ENDIF 255 256 ! 257 !-- Initialise control density 258 d = 0.0_wp 259 260 ! 261 !-- Determine maxima of the first and second derivative in x-direction 262 fmax_l = 0.0_wp 263 DO i = nxl, nxr 264 DO j = nys, nyn 265 DO k = nzb+1, nzt 266 numera = ABS( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) ) 267 denomi = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) ) 268 fmax_l(1) = MAX( fmax_l(1) , numera ) 269 fmax_l(2) = MAX( fmax_l(2) , denomi ) 270 ENDDO 271 ENDDO 272 ENDDO 273 #if defined( __parallel ) 274 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 275 CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr ) 276 #else 277 fmax = fmax_l 278 #endif 279 280 fmax = 0.04_wp * fmax 281 282 ! 283 !-- Allocate temporary arrays 284 ALLOCATE( a0(nzb+1:nzt,nxl-1:nxr+1), a1(nzb+1:nzt,nxl-1:nxr+1), & 285 a2(nzb+1:nzt,nxl-1:nxr+1), a12(nzb+1:nzt,nxl-1:nxr+1), & 286 a22(nzb+1:nzt,nxl-1:nxr+1), immb(nzb+1:nzt,nxl-1:nxr+1), & 287 imme(nzb+1:nzt,nxl-1:nxr+1), impb(nzb+1:nzt,nxl-1:nxr+1), & 288 impe(nzb+1:nzt,nxl-1:nxr+1), ipmb(nzb+1:nzt,nxl-1:nxr+1), & 289 ipme(nzb+1:nzt,nxl-1:nxr+1), ippb(nzb+1:nzt,nxl-1:nxr+1), & 290 ippe(nzb+1:nzt,nxl-1:nxr+1), m1(nzb+1:nzt,nxl-2:nxr+2), & 291 sw(nzb+1:nzt,nxl-1:nxr+1) & 292 ) 293 imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp 294 295 ! 296 !-- Initialise point of time measuring of the exponential portion (this would 297 !-- not work if done locally within the loop) 298 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'start' ) 299 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' ) 300 301 ! 302 !-- Outer loop of all j 200 303 DO j = nys, nyn 201 DO k = nzb, nzt+1 202 sk_p(k,j,i) = sk(k,j,i) 203 ENDDO 204 ENDDO 205 ENDDO 304 305 ! 306 !-- Compute polynomial coefficients 307 DO i = nxl-1, nxr+1 308 DO k = nzb+1, nzt 309 a12(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - sk_p(k,j,i-1) ) 310 a22(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) & 311 + sk_p(k,j,i-1) ) 312 a0(k,i) = ( 9.0_wp * sk_p(k,j,i+2) - 116.0_wp * sk_p(k,j,i+1) & 313 + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j,i-1) & 314 + 9.0_wp * sk_p(k,j,i-2) ) * f1920 315 a1(k,i) = ( -5.0_wp * sk_p(k,j,i+2) + 34.0_wp * sk_p(k,j,i+1) & 316 - 34.0_wp * sk_p(k,j,i-1) + 5.0_wp * sk_p(k,j,i-2) & 317 ) * f48 318 a2(k,i) = ( -3.0_wp * sk_p(k,j,i+2) + 36.0_wp * sk_p(k,j,i+1) & 319 - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j,i-1) & 320 - 3.0_wp * sk_p(k,j,i-2) ) * f48 321 ENDDO 322 ENDDO 323 324 ! 325 !-- Fluxes using the Bott scheme 326 !-- *VOCL LOOP,UNROLL(2) 327 DO i = nxl, nxr 328 DO k = nzb+1, nzt 329 cip = MAX( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx ) 330 cim = -MIN( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx ) 331 cipf = 1.0_wp - 2.0_wp * cip 332 cimf = 1.0_wp - 2.0_wp * cim 333 ip = a0(k,i) * f2 * ( 1.0_wp - cipf ) & 334 + a1(k,i) * f8 * ( 1.0_wp - cipf*cipf ) & 335 + a2(k,i) * f24 * ( 1.0_wp - cipf*cipf*cipf ) 336 im = a0(k,i+1) * f2 * ( 1.0_wp - cimf ) & 337 - a1(k,i+1) * f8 * ( 1.0_wp - cimf*cimf ) & 338 + a2(k,i+1) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 339 ip = MAX( ip, 0.0_wp ) 340 im = MAX( im, 0.0_wp ) 341 ippb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15_wp) ) 342 impb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i+1) / (ip+im+1E-15_wp) ) 343 344 cip = MAX( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx ) 345 cim = -MIN( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx ) 346 cipf = 1.0_wp - 2.0_wp * cip 347 cimf = 1.0_wp - 2.0_wp * cim 348 ip = a0(k,i-1) * f2 * ( 1.0_wp - cipf ) & 349 + a1(k,i-1) * f8 * ( 1.0_wp - cipf*cipf ) & 350 + a2(k,i-1) * f24 * ( 1.0_wp - cipf*cipf*cipf ) 351 im = a0(k,i) * f2 * ( 1.0_wp - cimf ) & 352 - a1(k,i) * f8 * ( 1.0_wp - cimf*cimf ) & 353 + a2(k,i) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 354 ip = MAX( ip, 0.0_wp ) 355 im = MAX( im, 0.0_wp ) 356 ipmb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i-1) / (ip+im+1E-15_wp) ) 357 immb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15_wp) ) 358 ENDDO 359 ENDDO 360 361 ! 362 !-- Compute monitor function m1 363 DO i = nxl-2, nxr+2 364 DO k = nzb+1, nzt 365 m1z = ABS( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) ) 366 m1n = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) ) 367 IF ( m1n /= 0.0_wp .AND. m1n >= m1z ) THEN 368 m1(k,i) = m1z / m1n 369 IF ( m1(k,i) /= 2.0_wp .AND. m1n < fmax(2) ) m1(k,i) = 0.0_wp 370 ELSEIF ( m1n < m1z ) THEN 371 m1(k,i) = -1.0_wp 372 ELSE 373 m1(k,i) = 0.0_wp 374 ENDIF 375 ENDDO 376 ENDDO 377 378 ! 379 !-- Compute switch sw 380 sw = 0.0_wp 381 DO i = nxl-1, nxr+1 382 DO k = nzb+1, nzt 383 m2 = 2.0_wp * ABS( a1(k,i) - a12(k,i) ) / & 384 MAX( ABS( a1(k,i) + a12(k,i) ), 1E-35_wp ) 385 IF ( ABS( a1(k,i) + a12(k,i) ) < fmax(2) ) m2 = 0.0_wp 386 387 m3 = 2.0_wp * ABS( a2(k,i) - a22(k,i) ) / & 388 MAX( ABS( a2(k,i) + a22(k,i) ), 1E-35_wp ) 389 IF ( ABS( a2(k,i) + a22(k,i) ) < fmax(1) ) m3 = 0.0_wp 390 391 t1 = 0.35_wp 392 t2 = 0.35_wp 393 IF ( m1(k,i) == -1.0_wp ) t2 = 0.12_wp 394 395 !-- *VOCL STMT,IF(10) 396 IF ( m1(k,i-1) == 1.0_wp .OR. m1(k,i) == 1.0_wp & 397 .OR. m1(k,i+1) == 1.0_wp .OR. m2 > t2 .OR. m3 > t2 .OR. & 398 ( m1(k,i) > t1 .AND. m1(k,i-1) /= -1.0_wp .AND. & 399 m1(k,i) /= -1.0_wp .AND. m1(k,i+1) /= -1.0_wp ) & 400 ) sw(k,i) = 1.0_wp 401 ENDDO 402 ENDDO 403 404 ! 405 !-- Fluxes using the exponential scheme 406 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' ) 407 DO i = nxl, nxr 408 DO k = nzb+1, nzt 409 410 !-- *VOCL STMT,IF(10) 411 IF ( sw(k,i) == 1.0_wp ) THEN 412 snenn = sk_p(k,j,i+1) - sk_p(k,j,i-1) 413 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 414 sterm = ( sk_p(k,j,i) - sk_p(k,j,i-1) ) / snenn 415 sterm = MIN( sterm, 0.9999_wp ) 416 sterm = MAX( sterm, 0.0001_wp ) 417 418 ix = INT( sterm * 1000 ) + 1 419 420 cip = MAX( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx ) 421 422 ippe(k,i) = sk_p(k,j,i-1) * cip + snenn * ( & 423 aex(ix) * cip + bex(ix) / dex(ix) * ( & 424 eex(ix) - & 425 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) ) & 426 ) & 427 ) 428 IF ( sterm == 0.0001_wp ) ippe(k,i) = sk_p(k,j,i) * cip 429 IF ( sterm == 0.9999_wp ) ippe(k,i) = sk_p(k,j,i) * cip 430 431 snenn = sk_p(k,j,i-1) - sk_p(k,j,i+1) 432 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 433 sterm = ( sk_p(k,j,i) - sk_p(k,j,i+1) ) / snenn 434 sterm = MIN( sterm, 0.9999_wp ) 435 sterm = MAX( sterm, 0.0001_wp ) 436 437 ix = INT( sterm * 1000 ) + 1 438 439 cim = -MIN( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx ) 440 441 imme(k,i) = sk_p(k,j,i+1) * cim + snenn * ( & 442 aex(ix) * cim + bex(ix) / dex(ix) * ( & 443 eex(ix) - & 444 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) & 445 ) & 446 ) 447 IF ( sterm == 0.0001_wp ) imme(k,i) = sk_p(k,j,i) * cim 448 IF ( sterm == 0.9999_wp ) imme(k,i) = sk_p(k,j,i) * cim 449 ENDIF 450 451 !-- *VOCL STMT,IF(10) 452 IF ( sw(k,i+1) == 1.0_wp ) THEN 453 snenn = sk_p(k,j,i) - sk_p(k,j,i+2) 454 IF ( ABS( snenn ) .LT. 1E-9_wp ) snenn = 1E-9_wp 455 sterm = ( sk_p(k,j,i+1) - sk_p(k,j,i+2) ) / snenn 456 sterm = MIN( sterm, 0.9999_wp ) 457 sterm = MAX( sterm, 0.0001_wp ) 458 459 ix = INT( sterm * 1000 ) + 1 460 461 cim = -MIN( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx ) 462 463 impe(k,i) = sk_p(k,j,i+2) * cim + snenn * ( & 464 aex(ix) * cim + bex(ix) / dex(ix) * ( & 465 eex(ix) - & 466 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) & 467 ) & 468 ) 469 IF ( sterm == 0.0001_wp ) impe(k,i) = sk_p(k,j,i+1) * cim 470 IF ( sterm == 0.9999_wp ) impe(k,i) = sk_p(k,j,i+1) * cim 471 ENDIF 472 473 !-- *VOCL STMT,IF(10) 474 IF ( sw(k,i-1) == 1.0_wp ) THEN 475 snenn = sk_p(k,j,i) - sk_p(k,j,i-2) 476 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 477 sterm = ( sk_p(k,j,i-1) - sk_p(k,j,i-2) ) / snenn 478 sterm = MIN( sterm, 0.9999_wp ) 479 sterm = MAX( sterm, 0.0001_wp ) 480 481 ix = INT( sterm * 1000 ) + 1 482 483 cip = MAX( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx ) 484 485 ipme(k,i) = sk_p(k,j,i-2) * cip + snenn * ( & 486 aex(ix) * cip + bex(ix) / dex(ix) * ( & 487 eex(ix) - & 488 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) ) & 489 ) & 490 ) 491 IF ( sterm == 0.0001_wp ) ipme(k,i) = sk_p(k,j,i-1) * cip 492 IF ( sterm == 0.9999_wp ) ipme(k,i) = sk_p(k,j,i-1) * cip 493 ENDIF 494 495 ENDDO 496 ENDDO 497 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' ) 498 499 ! 500 !-- Prognostic equation 501 DO i = nxl, nxr 502 DO k = nzb+1, nzt 503 fplus = ( 1.0_wp - sw(k,i) ) * ippb(k,i) + sw(k,i) * ippe(k,i) & 504 - ( 1.0_wp - sw(k,i+1) ) * impb(k,i) - sw(k,i+1) * impe(k,i) 505 fminus = ( 1.0_wp - sw(k,i-1) ) * ipmb(k,i) + sw(k,i-1) * ipme(k,i) & 506 - ( 1.0_wp - sw(k,i) ) * immb(k,i) - sw(k,i) * imme(k,i) 507 tendcy = fplus - fminus 508 ! 509 !-- Removed in order to optimize speed 510 ! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35_wp ) 511 ! IF ( ( ABS( tendcy ) / ffmax ) < 1E-7_wp ) tendcy = 0.0 512 ! 513 !-- Density correction because of possible remaining divergences 514 d_new = d(k,j,i) - ( u(k,j,i+1) - u(k,j,i) ) * dt_3d * ddx 515 sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) / & 516 ( 1.0_wp + d_new ) 517 d(k,j,i) = d_new 518 ENDDO 519 ENDDO 520 521 ENDDO ! End of the advection in x-direction 522 523 ! 524 !-- Deallocate temporary arrays 525 DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, & 526 ippb, ippe, m1, sw ) 527 528 529 ! 530 !-- Enlarge boundary of local array cyclically in y-direction 206 531 #if defined( __parallel ) 207 ngp = 2 * ( nzt - nzb + 6 ) * ( nyn - nys + 7 ) 208 CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'start' ) 209 ! 210 !-- Send left boundary, receive right boundary 211 CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxl+1), ngp, MPI_REAL, pleft, 0, & 212 sk_p(nzb-2,nys-3,nxr+2), ngp, MPI_REAL, pright, 0, & 213 comm2d, status, ierr ) 214 ! 215 !-- Send right boundary, receive left boundary 216 CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxr-2), ngp, MPI_REAL, pright, 1, & 217 sk_p(nzb-2,nys-3,nxl-3), ngp, MPI_REAL, pleft, 1, & 218 comm2d, status, ierr ) 219 CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' ) 532 ngp = ( nzt - nzb + 6 ) * ( nyn - nys + 7 ) 533 CALL MPI_TYPE_VECTOR( nxr-nxl+7, 3*(nzt-nzb+6), ngp, MPI_REAL, & 534 type_xz_2, ierr ) 535 CALL MPI_TYPE_COMMIT( type_xz_2, ierr ) 536 ! 537 !-- Send front boundary, receive rear boundary 538 CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' ) 539 CALL MPI_SENDRECV( sk_p(nzb-2,nys,nxl-3), 1, type_xz_2, psouth, 0, & 540 sk_p(nzb-2,nyn+1,nxl-3), 1, type_xz_2, pnorth, 0, & 541 comm2d, status, ierr ) 542 ! 543 !-- Send rear boundary, receive front boundary 544 CALL MPI_SENDRECV( sk_p(nzb-2,nyn-2,nxl-3), 1, type_xz_2, pnorth, 1, & 545 sk_p(nzb-2,nys-3,nxl-3), 1, type_xz_2, psouth, 1, & 546 comm2d, status, ierr ) 547 CALL MPI_TYPE_FREE( type_xz_2, ierr ) 548 CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' ) 220 549 #else 221 222 !223 !-- Cyclic boundary conditions224 sk_p(:,nys:nyn,nxl-3) = sk_p(:,nys:nyn,nxr-2)225 sk_p(:,nys:nyn,nxl-2) = sk_p(:,nys:nyn,nxr-1)226 sk_p(:,nys:nyn,nxr+2) = sk_p(:,nys:nyn,nxl+1)227 sk_p(:,nys:nyn,nxr+3) = sk_p(:,nys:nyn,nxl+2)228 #endif229 230 !231 !-- In case of a sloping surface, the additional gridpoints in x-direction232 !-- of the temperature field at the left and right boundary of the total233 !-- domain must be adjusted by the temperature difference between this distance234 IF ( sloping_surface .AND. sk_char == 'pt' ) THEN235 IF ( nxl == 0 ) THEN236 sk_p(:,nys:nyn,nxl-3) = sk_p(:,nys:nyn,nxl-3) - pt_slope_offset237 sk_p(:,nys:nyn,nxl-2) = sk_p(:,nys:nyn,nxl-2) - pt_slope_offset238 ENDIF239 IF ( nxr == nx ) THEN240 sk_p(:,nys:nyn,nxr+2) = sk_p(:,nys:nyn,nxr+2) + pt_slope_offset241 sk_p(:,nys:nyn,nxr+3) = sk_p(:,nys:nyn,nxr+3) + pt_slope_offset242 ENDIF243 ENDIF244 245 !246 !-- Initialise control density247 d = 0.0_wp248 249 !250 !-- Determine maxima of the first and second derivative in x-direction251 fmax_l = 0.0_wp252 DO i = nxl, nxr253 DO j = nys, nyn254 DO k = nzb+1, nzt255 numera = ABS( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) )256 denomi = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )257 fmax_l(1) = MAX( fmax_l(1) , numera )258 fmax_l(2) = MAX( fmax_l(2) , denomi )259 ENDDO260 ENDDO261 ENDDO262 #if defined( __parallel )263 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )264 CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr )265 #else266 fmax = fmax_l267 #endif268 269 fmax = 0.04_wp * fmax270 271 !272 !-- Allocate temporary arrays273 ALLOCATE( a0(nzb+1:nzt,nxl-1:nxr+1), a1(nzb+1:nzt,nxl-1:nxr+1), &274 a2(nzb+1:nzt,nxl-1:nxr+1), a12(nzb+1:nzt,nxl-1:nxr+1), &275 a22(nzb+1:nzt,nxl-1:nxr+1), immb(nzb+1:nzt,nxl-1:nxr+1), &276 imme(nzb+1:nzt,nxl-1:nxr+1), impb(nzb+1:nzt,nxl-1:nxr+1), &277 impe(nzb+1:nzt,nxl-1:nxr+1), ipmb(nzb+1:nzt,nxl-1:nxr+1), &278 ipme(nzb+1:nzt,nxl-1:nxr+1), ippb(nzb+1:nzt,nxl-1:nxr+1), &279 ippe(nzb+1:nzt,nxl-1:nxr+1), m1(nzb+1:nzt,nxl-2:nxr+2), &280 sw(nzb+1:nzt,nxl-1:nxr+1) &281 )282 imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp283 284 !285 !-- Initialise point of time measuring of the exponential portion (this would286 !-- not work if done locally within the loop)287 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'start' )288 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )289 290 !291 !-- Outer loop of all j292 DO j = nys, nyn293 294 !295 !-- Compute polynomial coefficients296 DO i = nxl-1, nxr+1297 DO k = nzb+1, nzt298 a12(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - sk_p(k,j,i-1) )299 a22(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) &300 + sk_p(k,j,i-1) )301 a0(k,i) = ( 9.0_wp * sk_p(k,j,i+2) - 116.0_wp * sk_p(k,j,i+1) &302 + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j,i-1) &303 + 9.0_wp * sk_p(k,j,i-2) ) * f1920304 a1(k,i) = ( -5.0_wp * sk_p(k,j,i+2) + 34.0_wp * sk_p(k,j,i+1) &305 - 34.0_wp * sk_p(k,j,i-1) + 5.0_wp * sk_p(k,j,i-2) &306 ) * f48307 a2(k,i) = ( -3.0_wp * sk_p(k,j,i+2) + 36.0_wp * sk_p(k,j,i+1) &308 - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j,i-1) &309 - 3.0_wp * sk_p(k,j,i-2) ) * f48310 ENDDO311 ENDDO312 313 !314 !-- Fluxes using the Bott scheme315 !-- *VOCL LOOP,UNROLL(2)316 550 DO i = nxl, nxr 317 551 DO k = nzb+1, nzt 318 cip = MAX( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx ) 319 cim = -MIN( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx ) 320 cipf = 1.0_wp - 2.0_wp * cip 321 cimf = 1.0_wp - 2.0_wp * cim 322 ip = a0(k,i) * f2 * ( 1.0_wp - cipf ) & 323 + a1(k,i) * f8 * ( 1.0_wp - cipf*cipf ) & 324 + a2(k,i) * f24 * ( 1.0_wp - cipf*cipf*cipf ) 325 im = a0(k,i+1) * f2 * ( 1.0_wp - cimf ) & 326 - a1(k,i+1) * f8 * ( 1.0_wp - cimf*cimf ) & 327 + a2(k,i+1) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 328 ip = MAX( ip, 0.0_wp ) 329 im = MAX( im, 0.0_wp ) 330 ippb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15_wp) ) 331 impb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i+1) / (ip+im+1E-15_wp) ) 332 333 cip = MAX( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx ) 334 cim = -MIN( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx ) 335 cipf = 1.0_wp - 2.0_wp * cip 336 cimf = 1.0_wp - 2.0_wp * cim 337 ip = a0(k,i-1) * f2 * ( 1.0_wp - cipf ) & 338 + a1(k,i-1) * f8 * ( 1.0_wp - cipf*cipf ) & 339 + a2(k,i-1) * f24 * ( 1.0_wp - cipf*cipf*cipf ) 340 im = a0(k,i) * f2 * ( 1.0_wp - cimf ) & 341 - a1(k,i) * f8 * ( 1.0_wp - cimf*cimf ) & 342 + a2(k,i) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 343 ip = MAX( ip, 0.0_wp ) 344 im = MAX( im, 0.0_wp ) 345 ipmb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i-1) / (ip+im+1E-15_wp) ) 346 immb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15_wp) ) 552 sk_p(k,nys-1,i) = sk_p(k,nyn,i) 553 sk_p(k,nys-2,i) = sk_p(k,nyn-1,i) 554 sk_p(k,nys-3,i) = sk_p(k,nyn-2,i) 555 sk_p(k,nyn+1,i) = sk_p(k,nys,i) 556 sk_p(k,nyn+2,i) = sk_p(k,nys+1,i) 557 sk_p(k,nyn+3,i) = sk_p(k,nys+2,i) 347 558 ENDDO 348 559 ENDDO 349 350 ! 351 !-- Compute monitor function m1 352 DO i = nxl-2, nxr+2 353 DO k = nzb+1, nzt 354 m1z = ABS( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) ) 355 m1n = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) ) 356 IF ( m1n /= 0.0_wp .AND. m1n >= m1z ) THEN 357 m1(k,i) = m1z / m1n 358 IF ( m1(k,i) /= 2.0_wp .AND. m1n < fmax(2) ) m1(k,i) = 0.0_wp 359 ELSEIF ( m1n < m1z ) THEN 360 m1(k,i) = -1.0_wp 361 ELSE 362 m1(k,i) = 0.0_wp 363 ENDIF 560 #endif 561 562 ! 563 !-- Determine the maxima of the first and second derivative in y-direction 564 fmax_l = 0.0_wp 565 DO i = nxl, nxr 566 DO j = nys, nyn 567 DO k = nzb+1, nzt 568 numera = ABS( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j-1,i) ) 569 denomi = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) ) 570 fmax_l(1) = MAX( fmax_l(1) , numera ) 571 fmax_l(2) = MAX( fmax_l(2) , denomi ) 572 ENDDO 364 573 ENDDO 365 574 ENDDO 366 367 ! 368 !-- Compute switch sw 369 sw = 0.0_wp 370 DO i = nxl-1, nxr+1 371 DO k = nzb+1, nzt 372 m2 = 2.0_wp * ABS( a1(k,i) - a12(k,i) ) / & 373 MAX( ABS( a1(k,i) + a12(k,i) ), 1E-35_wp ) 374 IF ( ABS( a1(k,i) + a12(k,i) ) < fmax(2) ) m2 = 0.0_wp 375 376 m3 = 2.0_wp * ABS( a2(k,i) - a22(k,i) ) / & 377 MAX( ABS( a2(k,i) + a22(k,i) ), 1E-35_wp ) 378 IF ( ABS( a2(k,i) + a22(k,i) ) < fmax(1) ) m3 = 0.0_wp 379 380 t1 = 0.35_wp 381 t2 = 0.35_wp 382 IF ( m1(k,i) == -1.0_wp ) t2 = 0.12_wp 383 384 !-- *VOCL STMT,IF(10) 385 IF ( m1(k,i-1) == 1.0_wp .OR. m1(k,i) == 1.0_wp & 386 .OR. m1(k,i+1) == 1.0_wp .OR. m2 > t2 .OR. m3 > t2 .OR. & 387 ( m1(k,i) > t1 .AND. m1(k,i-1) /= -1.0_wp .AND. & 388 m1(k,i) /= -1.0_wp .AND. m1(k,i+1) /= -1.0_wp ) & 389 ) sw(k,i) = 1.0_wp 390 ENDDO 391 ENDDO 392 393 ! 394 !-- Fluxes using the exponential scheme 395 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' ) 575 #if defined( __parallel ) 576 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 577 CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr ) 578 #else 579 fmax = fmax_l 580 #endif 581 582 fmax = 0.04_wp * fmax 583 584 ! 585 !-- Allocate temporary arrays 586 ALLOCATE( a0(nzb+1:nzt,nys-1:nyn+1), a1(nzb+1:nzt,nys-1:nyn+1), & 587 a2(nzb+1:nzt,nys-1:nyn+1), a12(nzb+1:nzt,nys-1:nyn+1), & 588 a22(nzb+1:nzt,nys-1:nyn+1), immb(nzb+1:nzt,nys-1:nyn+1), & 589 imme(nzb+1:nzt,nys-1:nyn+1), impb(nzb+1:nzt,nys-1:nyn+1), & 590 impe(nzb+1:nzt,nys-1:nyn+1), ipmb(nzb+1:nzt,nys-1:nyn+1), & 591 ipme(nzb+1:nzt,nys-1:nyn+1), ippb(nzb+1:nzt,nys-1:nyn+1), & 592 ippe(nzb+1:nzt,nys-1:nyn+1), m1(nzb+1:nzt,nys-2:nyn+2), & 593 sw(nzb+1:nzt,nys-1:nyn+1) & 594 ) 595 imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp 596 597 ! 598 !-- Outer loop of all i 396 599 DO i = nxl, nxr 397 DO k = nzb+1, nzt 398 399 !-- *VOCL STMT,IF(10) 400 IF ( sw(k,i) == 1.0_wp ) THEN 401 snenn = sk_p(k,j,i+1) - sk_p(k,j,i-1) 402 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 403 sterm = ( sk_p(k,j,i) - sk_p(k,j,i-1) ) / snenn 404 sterm = MIN( sterm, 0.9999_wp ) 405 sterm = MAX( sterm, 0.0001_wp ) 406 407 ix = INT( sterm * 1000 ) + 1 408 409 cip = MAX( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx ) 410 411 ippe(k,i) = sk_p(k,j,i-1) * cip + snenn * ( & 412 aex(ix) * cip + bex(ix) / dex(ix) * ( & 413 eex(ix) - & 414 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) ) & 415 ) & 416 ) 417 IF ( sterm == 0.0001_wp ) ippe(k,i) = sk_p(k,j,i) * cip 418 IF ( sterm == 0.9999_wp ) ippe(k,i) = sk_p(k,j,i) * cip 419 420 snenn = sk_p(k,j,i-1) - sk_p(k,j,i+1) 421 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 422 sterm = ( sk_p(k,j,i) - sk_p(k,j,i+1) ) / snenn 423 sterm = MIN( sterm, 0.9999_wp ) 424 sterm = MAX( sterm, 0.0001_wp ) 425 426 ix = INT( sterm * 1000 ) + 1 427 428 cim = -MIN( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx ) 429 430 imme(k,i) = sk_p(k,j,i+1) * cim + snenn * ( & 431 aex(ix) * cim + bex(ix) / dex(ix) * ( & 432 eex(ix) - & 433 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) & 434 ) & 435 ) 436 IF ( sterm == 0.0001_wp ) imme(k,i) = sk_p(k,j,i) * cim 437 IF ( sterm == 0.9999_wp ) imme(k,i) = sk_p(k,j,i) * cim 438 ENDIF 439 440 !-- *VOCL STMT,IF(10) 441 IF ( sw(k,i+1) == 1.0_wp ) THEN 442 snenn = sk_p(k,j,i) - sk_p(k,j,i+2) 443 IF ( ABS( snenn ) .LT. 1E-9_wp ) snenn = 1E-9_wp 444 sterm = ( sk_p(k,j,i+1) - sk_p(k,j,i+2) ) / snenn 445 sterm = MIN( sterm, 0.9999_wp ) 446 sterm = MAX( sterm, 0.0001_wp ) 447 448 ix = INT( sterm * 1000 ) + 1 449 450 cim = -MIN( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx ) 451 452 impe(k,i) = sk_p(k,j,i+2) * cim + snenn * ( & 453 aex(ix) * cim + bex(ix) / dex(ix) * ( & 454 eex(ix) - & 455 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) & 456 ) & 457 ) 458 IF ( sterm == 0.0001_wp ) impe(k,i) = sk_p(k,j,i+1) * cim 459 IF ( sterm == 0.9999_wp ) impe(k,i) = sk_p(k,j,i+1) * cim 460 ENDIF 461 462 !-- *VOCL STMT,IF(10) 463 IF ( sw(k,i-1) == 1.0_wp ) THEN 464 snenn = sk_p(k,j,i) - sk_p(k,j,i-2) 465 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 466 sterm = ( sk_p(k,j,i-1) - sk_p(k,j,i-2) ) / snenn 467 sterm = MIN( sterm, 0.9999_wp ) 468 sterm = MAX( sterm, 0.0001_wp ) 469 470 ix = INT( sterm * 1000 ) + 1 471 472 cip = MAX( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx ) 473 474 ipme(k,i) = sk_p(k,j,i-2) * cip + snenn * ( & 475 aex(ix) * cip + bex(ix) / dex(ix) * ( & 476 eex(ix) - & 477 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) ) & 478 ) & 479 ) 480 IF ( sterm == 0.0001_wp ) ipme(k,i) = sk_p(k,j,i-1) * cip 481 IF ( sterm == 0.9999_wp ) ipme(k,i) = sk_p(k,j,i-1) * cip 482 ENDIF 483 484 ENDDO 485 ENDDO 486 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' ) 487 488 ! 489 !-- Prognostic equation 490 DO i = nxl, nxr 491 DO k = nzb+1, nzt 492 fplus = ( 1.0_wp - sw(k,i) ) * ippb(k,i) + sw(k,i) * ippe(k,i) & 493 - ( 1.0_wp - sw(k,i+1) ) * impb(k,i) - sw(k,i+1) * impe(k,i) 494 fminus = ( 1.0_wp - sw(k,i-1) ) * ipmb(k,i) + sw(k,i-1) * ipme(k,i) & 495 - ( 1.0_wp - sw(k,i) ) * immb(k,i) - sw(k,i) * imme(k,i) 496 tendcy = fplus - fminus 497 ! 498 !-- Removed in order to optimize speed 499 ! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35_wp ) 500 ! IF ( ( ABS( tendcy ) / ffmax ) < 1E-7_wp ) tendcy = 0.0 501 ! 502 !-- Density correction because of possible remaining divergences 503 d_new = d(k,j,i) - ( u(k,j,i+1) - u(k,j,i) ) * dt_3d * ddx 504 sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) / & 505 ( 1.0_wp + d_new ) 506 d(k,j,i) = d_new 507 ENDDO 508 ENDDO 509 510 ENDDO ! End of the advection in x-direction 511 512 ! 513 !-- Deallocate temporary arrays 514 DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, & 515 ippb, ippe, m1, sw ) 516 517 518 ! 519 !-- Enlarge boundary of local array cyclically in y-direction 520 #if defined( __parallel ) 521 ngp = ( nzt - nzb + 6 ) * ( nyn - nys + 7 ) 522 CALL MPI_TYPE_VECTOR( nxr-nxl+7, 3*(nzt-nzb+6), ngp, MPI_REAL, & 523 type_xz_2, ierr ) 524 CALL MPI_TYPE_COMMIT( type_xz_2, ierr ) 525 ! 526 !-- Send front boundary, receive rear boundary 527 CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' ) 528 CALL MPI_SENDRECV( sk_p(nzb-2,nys,nxl-3), 1, type_xz_2, psouth, 0, & 529 sk_p(nzb-2,nyn+1,nxl-3), 1, type_xz_2, pnorth, 0, & 530 comm2d, status, ierr ) 531 ! 532 !-- Send rear boundary, receive front boundary 533 CALL MPI_SENDRECV( sk_p(nzb-2,nyn-2,nxl-3), 1, type_xz_2, pnorth, 1, & 534 sk_p(nzb-2,nys-3,nxl-3), 1, type_xz_2, psouth, 1, & 535 comm2d, status, ierr ) 536 CALL MPI_TYPE_FREE( type_xz_2, ierr ) 537 CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' ) 538 #else 539 DO i = nxl, nxr 540 DO k = nzb+1, nzt 541 sk_p(k,nys-1,i) = sk_p(k,nyn,i) 542 sk_p(k,nys-2,i) = sk_p(k,nyn-1,i) 543 sk_p(k,nys-3,i) = sk_p(k,nyn-2,i) 544 sk_p(k,nyn+1,i) = sk_p(k,nys,i) 545 sk_p(k,nyn+2,i) = sk_p(k,nys+1,i) 546 sk_p(k,nyn+3,i) = sk_p(k,nys+2,i) 547 ENDDO 548 ENDDO 549 #endif 550 551 ! 552 !-- Determine the maxima of the first and second derivative in y-direction 553 fmax_l = 0.0_wp 554 DO i = nxl, nxr 555 DO j = nys, nyn 556 DO k = nzb+1, nzt 557 numera = ABS( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j-1,i) ) 558 denomi = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) ) 559 fmax_l(1) = MAX( fmax_l(1) , numera ) 560 fmax_l(2) = MAX( fmax_l(2) , denomi ) 561 ENDDO 562 ENDDO 563 ENDDO 564 #if defined( __parallel ) 565 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 566 CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr ) 567 #else 568 fmax = fmax_l 569 #endif 570 571 fmax = 0.04_wp * fmax 572 573 ! 574 !-- Allocate temporary arrays 575 ALLOCATE( a0(nzb+1:nzt,nys-1:nyn+1), a1(nzb+1:nzt,nys-1:nyn+1), & 576 a2(nzb+1:nzt,nys-1:nyn+1), a12(nzb+1:nzt,nys-1:nyn+1), & 577 a22(nzb+1:nzt,nys-1:nyn+1), immb(nzb+1:nzt,nys-1:nyn+1), & 578 imme(nzb+1:nzt,nys-1:nyn+1), impb(nzb+1:nzt,nys-1:nyn+1), & 579 impe(nzb+1:nzt,nys-1:nyn+1), ipmb(nzb+1:nzt,nys-1:nyn+1), & 580 ipme(nzb+1:nzt,nys-1:nyn+1), ippb(nzb+1:nzt,nys-1:nyn+1), & 581 ippe(nzb+1:nzt,nys-1:nyn+1), m1(nzb+1:nzt,nys-2:nyn+2), & 582 sw(nzb+1:nzt,nys-1:nyn+1) & 583 ) 584 imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp 585 586 ! 587 !-- Outer loop of all i 588 DO i = nxl, nxr 589 590 ! 591 !-- Compute polynomial coefficients 592 DO j = nys-1, nyn+1 593 DO k = nzb+1, nzt 594 a12(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - sk_p(k,j-1,i) ) 595 a22(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) & 596 + sk_p(k,j-1,i) ) 597 a0(k,j) = ( 9.0_wp * sk_p(k,j+2,i) - 116.0_wp * sk_p(k,j+1,i) & 598 + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j-1,i) & 599 + 9.0_wp * sk_p(k,j-2,i) ) * f1920 600 a1(k,j) = ( -5.0_wp * sk_p(k,j+2,i) + 34.0_wp * sk_p(k,j+1,i) & 601 - 34.0_wp * sk_p(k,j-1,i) + 5.0_wp * sk_p(k,j-2,i) & 602 ) * f48 603 a2(k,j) = ( -3.0_wp * sk_p(k,j+2,i) + 36.0_wp * sk_p(k,j+1,i) & 604 - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j-1,i) & 605 - 3.0_wp * sk_p(k,j-2,i) ) * f48 606 ENDDO 607 ENDDO 608 609 ! 610 !-- Fluxes using the Bott scheme 611 !-- *VOCL LOOP,UNROLL(2) 612 DO j = nys, nyn 613 DO k = nzb+1, nzt 614 cip = MAX( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy ) 615 cim = -MIN( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy ) 616 cipf = 1.0_wp - 2.0_wp * cip 617 cimf = 1.0_wp - 2.0_wp * cim 618 ip = a0(k,j) * f2 * ( 1.0_wp - cipf ) & 619 + a1(k,j) * f8 * ( 1.0_wp - cipf*cipf ) & 620 + a2(k,j) * f24 * ( 1.0_wp - cipf*cipf*cipf ) 621 im = a0(k,j+1) * f2 * ( 1.0_wp - cimf ) & 622 - a1(k,j+1) * f8 * ( 1.0_wp - cimf*cimf ) & 623 + a2(k,j+1) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 624 ip = MAX( ip, 0.0_wp ) 625 im = MAX( im, 0.0_wp ) 626 ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15_wp) ) 627 impb(k,j) = im * MIN( 1.0_wp, sk_p(k,j+1,i) / (ip+im+1E-15_wp) ) 628 629 cip = MAX( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy ) 630 cim = -MIN( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy ) 631 cipf = 1.0_wp - 2.0_wp * cip 632 cimf = 1.0_wp - 2.0_wp * cim 633 ip = a0(k,j-1) * f2 * ( 1.0_wp - cipf ) & 634 + a1(k,j-1) * f8 * ( 1.0_wp - cipf*cipf ) & 635 + a2(k,j-1) * f24 * ( 1.0_wp - cipf*cipf*cipf ) 636 im = a0(k,j) * f2 * ( 1.0_wp - cimf ) & 637 - a1(k,j) * f8 * ( 1.0_wp - cimf*cimf ) & 638 + a2(k,j) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 639 ip = MAX( ip, 0.0_wp ) 640 im = MAX( im, 0.0_wp ) 641 ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j-1,i) / (ip+im+1E-15_wp) ) 642 immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15_wp) ) 643 ENDDO 644 ENDDO 645 646 ! 647 !-- Compute monitor function m1 648 DO j = nys-2, nyn+2 649 DO k = nzb+1, nzt 650 m1z = ABS( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j-1,i) ) 651 m1n = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) ) 652 IF ( m1n /= 0.0_wp .AND. m1n >= m1z ) THEN 653 m1(k,j) = m1z / m1n 654 IF ( m1(k,j) /= 2.0_wp .AND. m1n < fmax(2) ) m1(k,j) = 0.0_wp 655 ELSEIF ( m1n < m1z ) THEN 656 m1(k,j) = -1.0_wp 657 ELSE 658 m1(k,j) = 0.0_wp 659 ENDIF 660 ENDDO 661 ENDDO 662 663 ! 664 !-- Compute switch sw 665 sw = 0.0_wp 666 DO j = nys-1, nyn+1 667 DO k = nzb+1, nzt 668 m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) / & 669 MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp ) 670 IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0_wp 671 672 m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) / & 673 MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp ) 674 IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0_wp 675 676 t1 = 0.35_wp 677 t2 = 0.35_wp 678 IF ( m1(k,j) == -1.0_wp ) t2 = 0.12_wp 679 680 !-- *VOCL STMT,IF(10) 681 IF ( m1(k,j-1) == 1.0_wp .OR. m1(k,j) == 1.0_wp & 682 .OR. m1(k,j+1) == 1.0_wp .OR. m2 > t2 .OR. m3 > t2 .OR. & 683 ( m1(k,j) > t1 .AND. m1(k,j-1) /= -1.0_wp .AND. & 684 m1(k,j) /= -1.0_wp .AND. m1(k,j+1) /= -1.0_wp ) & 685 ) sw(k,j) = 1.0_wp 686 ENDDO 687 ENDDO 688 689 ! 690 !-- Fluxes using exponential scheme 691 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' ) 692 DO j = nys, nyn 693 DO k = nzb+1, nzt 694 695 !-- *VOCL STMT,IF(10) 696 IF ( sw(k,j) == 1.0_wp ) THEN 697 snenn = sk_p(k,j+1,i) - sk_p(k,j-1,i) 698 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 699 sterm = ( sk_p(k,j,i) - sk_p(k,j-1,i) ) / snenn 700 sterm = MIN( sterm, 0.9999_wp ) 701 sterm = MAX( sterm, 0.0001_wp ) 702 703 ix = INT( sterm * 1000 ) + 1 704 705 cip = MAX( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy ) 706 707 ippe(k,j) = sk_p(k,j-1,i) * cip + snenn * ( & 708 aex(ix) * cip + bex(ix) / dex(ix) * ( & 709 eex(ix) - & 710 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) ) & 711 ) & 712 ) 713 IF ( sterm == 0.0001_wp ) ippe(k,j) = sk_p(k,j,i) * cip 714 IF ( sterm == 0.9999_wp ) ippe(k,j) = sk_p(k,j,i) * cip 715 716 snenn = sk_p(k,j-1,i) - sk_p(k,j+1,i) 717 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 718 sterm = ( sk_p(k,j,i) - sk_p(k,j+1,i) ) / snenn 719 sterm = MIN( sterm, 0.9999_wp ) 720 sterm = MAX( sterm, 0.0001_wp ) 721 722 ix = INT( sterm * 1000 ) + 1 723 724 cim = -MIN( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy ) 725 726 imme(k,j) = sk_p(k,j+1,i) * cim + snenn * ( & 727 aex(ix) * cim + bex(ix) / dex(ix) * ( & 728 eex(ix) - & 729 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) & 730 ) & 731 ) 732 IF ( sterm == 0.0001_wp ) imme(k,j) = sk_p(k,j,i) * cim 733 IF ( sterm == 0.9999_wp ) imme(k,j) = sk_p(k,j,i) * cim 734 ENDIF 735 736 !-- *VOCL STMT,IF(10) 737 IF ( sw(k,j+1) == 1.0_wp ) THEN 738 snenn = sk_p(k,j,i) - sk_p(k,j+2,i) 739 IF ( ABS( snenn ) .LT. 1E-9_wp ) snenn = 1E-9_wp 740 sterm = ( sk_p(k,j+1,i) - sk_p(k,j+2,i) ) / snenn 741 sterm = MIN( sterm, 0.9999_wp ) 742 sterm = MAX( sterm, 0.0001_wp ) 743 744 ix = INT( sterm * 1000 ) + 1 745 746 cim = -MIN( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy ) 747 748 impe(k,j) = sk_p(k,j+2,i) * cim + snenn * ( & 749 aex(ix) * cim + bex(ix) / dex(ix) * ( & 750 eex(ix) - & 751 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) & 752 ) & 753 ) 754 IF ( sterm == 0.0001_wp ) impe(k,j) = sk_p(k,j+1,i) * cim 755 IF ( sterm == 0.9999_wp ) impe(k,j) = sk_p(k,j+1,i) * cim 756 ENDIF 757 758 !-- *VOCL STMT,IF(10) 759 IF ( sw(k,j-1) == 1.0_wp ) THEN 760 snenn = sk_p(k,j,i) - sk_p(k,j-2,i) 761 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 762 sterm = ( sk_p(k,j-1,i) - sk_p(k,j-2,i) ) / snenn 763 sterm = MIN( sterm, 0.9999_wp ) 764 sterm = MAX( sterm, 0.0001_wp ) 765 766 ix = INT( sterm * 1000 ) + 1 767 768 cip = MAX( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy ) 769 770 ipme(k,j) = sk_p(k,j-2,i) * cip + snenn * ( & 771 aex(ix) * cip + bex(ix) / dex(ix) * ( & 772 eex(ix) - & 773 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) ) & 774 ) & 775 ) 776 IF ( sterm == 0.0001_wp ) ipme(k,j) = sk_p(k,j-1,i) * cip 777 IF ( sterm == 0.9999_wp ) ipme(k,j) = sk_p(k,j-1,i) * cip 778 ENDIF 779 780 ENDDO 781 ENDDO 782 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' ) 783 784 ! 785 !-- Prognostic equation 786 DO j = nys, nyn 787 DO k = nzb+1, nzt 788 fplus = ( 1.0_wp - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j) & 789 - ( 1.0_wp - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j) 790 fminus = ( 1.0_wp - sw(k,j-1) ) * ipmb(k,j) + sw(k,j-1) * ipme(k,j) & 791 - ( 1.0_wp - sw(k,j) ) * immb(k,j) - sw(k,j) * imme(k,j) 792 tendcy = fplus - fminus 793 ! 794 !-- Removed in order to optimise speed 795 ! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35_wp ) 796 ! IF ( ( ABS( tendcy ) / ffmax ) < 1E-7_wp ) tendcy = 0.0 797 ! 798 !-- Density correction because of possible remaining divergences 799 d_new = d(k,j,i) - ( v(k,j+1,i) - v(k,j,i) ) * dt_3d * ddy 800 sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) / & 801 ( 1.0_wp + d_new ) 802 d(k,j,i) = d_new 803 ENDDO 804 ENDDO 805 806 ENDDO ! End of the advection in y-direction 807 CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' ) 808 CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'stop' ) 809 810 ! 811 !-- Deallocate temporary arrays 812 DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, & 813 ippb, ippe, m1, sw ) 814 815 816 ! 817 !-- Initialise for the computation of heat fluxes (see below; required in 818 !-- UP flow_statistics) 819 IF ( sk_char == 'pt' ) sums_wsts_bc_l = 0.0_wp 820 821 ! 822 !-- Add top and bottom boundaries according to the relevant boundary conditions 823 IF ( sk_char == 'pt' ) THEN 824 825 ! 826 !-- Temperature boundary condition at the bottom boundary 827 IF ( ibc_pt_b == 0 ) THEN 600 601 ! 602 !-- Compute polynomial coefficients 603 DO j = nys-1, nyn+1 604 DO k = nzb+1, nzt 605 a12(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - sk_p(k,j-1,i) ) 606 a22(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) & 607 + sk_p(k,j-1,i) ) 608 a0(k,j) = ( 9.0_wp * sk_p(k,j+2,i) - 116.0_wp * sk_p(k,j+1,i) & 609 + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j-1,i) & 610 + 9.0_wp * sk_p(k,j-2,i) ) * f1920 611 a1(k,j) = ( -5.0_wp * sk_p(k,j+2,i) + 34.0_wp * sk_p(k,j+1,i) & 612 - 34.0_wp * sk_p(k,j-1,i) + 5.0_wp * sk_p(k,j-2,i) & 613 ) * f48 614 a2(k,j) = ( -3.0_wp * sk_p(k,j+2,i) + 36.0_wp * sk_p(k,j+1,i) & 615 - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j-1,i) & 616 - 3.0_wp * sk_p(k,j-2,i) ) * f48 617 ENDDO 618 ENDDO 619 620 ! 621 !-- Fluxes using the Bott scheme 622 !-- *VOCL LOOP,UNROLL(2) 623 DO j = nys, nyn 624 DO k = nzb+1, nzt 625 cip = MAX( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy ) 626 cim = -MIN( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy ) 627 cipf = 1.0_wp - 2.0_wp * cip 628 cimf = 1.0_wp - 2.0_wp * cim 629 ip = a0(k,j) * f2 * ( 1.0_wp - cipf ) & 630 + a1(k,j) * f8 * ( 1.0_wp - cipf*cipf ) & 631 + a2(k,j) * f24 * ( 1.0_wp - cipf*cipf*cipf ) 632 im = a0(k,j+1) * f2 * ( 1.0_wp - cimf ) & 633 - a1(k,j+1) * f8 * ( 1.0_wp - cimf*cimf ) & 634 + a2(k,j+1) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 635 ip = MAX( ip, 0.0_wp ) 636 im = MAX( im, 0.0_wp ) 637 ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15_wp) ) 638 impb(k,j) = im * MIN( 1.0_wp, sk_p(k,j+1,i) / (ip+im+1E-15_wp) ) 639 640 cip = MAX( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy ) 641 cim = -MIN( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy ) 642 cipf = 1.0_wp - 2.0_wp * cip 643 cimf = 1.0_wp - 2.0_wp * cim 644 ip = a0(k,j-1) * f2 * ( 1.0_wp - cipf ) & 645 + a1(k,j-1) * f8 * ( 1.0_wp - cipf*cipf ) & 646 + a2(k,j-1) * f24 * ( 1.0_wp - cipf*cipf*cipf ) 647 im = a0(k,j) * f2 * ( 1.0_wp - cimf ) & 648 - a1(k,j) * f8 * ( 1.0_wp - cimf*cimf ) & 649 + a2(k,j) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 650 ip = MAX( ip, 0.0_wp ) 651 im = MAX( im, 0.0_wp ) 652 ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j-1,i) / (ip+im+1E-15_wp) ) 653 immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15_wp) ) 654 ENDDO 655 ENDDO 656 657 ! 658 !-- Compute monitor function m1 659 DO j = nys-2, nyn+2 660 DO k = nzb+1, nzt 661 m1z = ABS( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j-1,i) ) 662 m1n = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) ) 663 IF ( m1n /= 0.0_wp .AND. m1n >= m1z ) THEN 664 m1(k,j) = m1z / m1n 665 IF ( m1(k,j) /= 2.0_wp .AND. m1n < fmax(2) ) m1(k,j) = 0.0_wp 666 ELSEIF ( m1n < m1z ) THEN 667 m1(k,j) = -1.0_wp 668 ELSE 669 m1(k,j) = 0.0_wp 670 ENDIF 671 ENDDO 672 ENDDO 673 674 ! 675 !-- Compute switch sw 676 sw = 0.0_wp 677 DO j = nys-1, nyn+1 678 DO k = nzb+1, nzt 679 m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) / & 680 MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp ) 681 IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0_wp 682 683 m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) / & 684 MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp ) 685 IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0_wp 686 687 t1 = 0.35_wp 688 t2 = 0.35_wp 689 IF ( m1(k,j) == -1.0_wp ) t2 = 0.12_wp 690 691 !-- *VOCL STMT,IF(10) 692 IF ( m1(k,j-1) == 1.0_wp .OR. m1(k,j) == 1.0_wp & 693 .OR. m1(k,j+1) == 1.0_wp .OR. m2 > t2 .OR. m3 > t2 .OR. & 694 ( m1(k,j) > t1 .AND. m1(k,j-1) /= -1.0_wp .AND. & 695 m1(k,j) /= -1.0_wp .AND. m1(k,j+1) /= -1.0_wp ) & 696 ) sw(k,j) = 1.0_wp 697 ENDDO 698 ENDDO 699 700 ! 701 !-- Fluxes using exponential scheme 702 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' ) 703 DO j = nys, nyn 704 DO k = nzb+1, nzt 705 706 !-- *VOCL STMT,IF(10) 707 IF ( sw(k,j) == 1.0_wp ) THEN 708 snenn = sk_p(k,j+1,i) - sk_p(k,j-1,i) 709 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 710 sterm = ( sk_p(k,j,i) - sk_p(k,j-1,i) ) / snenn 711 sterm = MIN( sterm, 0.9999_wp ) 712 sterm = MAX( sterm, 0.0001_wp ) 713 714 ix = INT( sterm * 1000 ) + 1 715 716 cip = MAX( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy ) 717 718 ippe(k,j) = sk_p(k,j-1,i) * cip + snenn * ( & 719 aex(ix) * cip + bex(ix) / dex(ix) * ( & 720 eex(ix) - & 721 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) ) & 722 ) & 723 ) 724 IF ( sterm == 0.0001_wp ) ippe(k,j) = sk_p(k,j,i) * cip 725 IF ( sterm == 0.9999_wp ) ippe(k,j) = sk_p(k,j,i) * cip 726 727 snenn = sk_p(k,j-1,i) - sk_p(k,j+1,i) 728 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 729 sterm = ( sk_p(k,j,i) - sk_p(k,j+1,i) ) / snenn 730 sterm = MIN( sterm, 0.9999_wp ) 731 sterm = MAX( sterm, 0.0001_wp ) 732 733 ix = INT( sterm * 1000 ) + 1 734 735 cim = -MIN( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy ) 736 737 imme(k,j) = sk_p(k,j+1,i) * cim + snenn * ( & 738 aex(ix) * cim + bex(ix) / dex(ix) * ( & 739 eex(ix) - & 740 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) & 741 ) & 742 ) 743 IF ( sterm == 0.0001_wp ) imme(k,j) = sk_p(k,j,i) * cim 744 IF ( sterm == 0.9999_wp ) imme(k,j) = sk_p(k,j,i) * cim 745 ENDIF 746 747 !-- *VOCL STMT,IF(10) 748 IF ( sw(k,j+1) == 1.0_wp ) THEN 749 snenn = sk_p(k,j,i) - sk_p(k,j+2,i) 750 IF ( ABS( snenn ) .LT. 1E-9_wp ) snenn = 1E-9_wp 751 sterm = ( sk_p(k,j+1,i) - sk_p(k,j+2,i) ) / snenn 752 sterm = MIN( sterm, 0.9999_wp ) 753 sterm = MAX( sterm, 0.0001_wp ) 754 755 ix = INT( sterm * 1000 ) + 1 756 757 cim = -MIN( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy ) 758 759 impe(k,j) = sk_p(k,j+2,i) * cim + snenn * ( & 760 aex(ix) * cim + bex(ix) / dex(ix) * ( & 761 eex(ix) - & 762 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) & 763 ) & 764 ) 765 IF ( sterm == 0.0001_wp ) impe(k,j) = sk_p(k,j+1,i) * cim 766 IF ( sterm == 0.9999_wp ) impe(k,j) = sk_p(k,j+1,i) * cim 767 ENDIF 768 769 !-- *VOCL STMT,IF(10) 770 IF ( sw(k,j-1) == 1.0_wp ) THEN 771 snenn = sk_p(k,j,i) - sk_p(k,j-2,i) 772 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 773 sterm = ( sk_p(k,j-1,i) - sk_p(k,j-2,i) ) / snenn 774 sterm = MIN( sterm, 0.9999_wp ) 775 sterm = MAX( sterm, 0.0001_wp ) 776 777 ix = INT( sterm * 1000 ) + 1 778 779 cip = MAX( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy ) 780 781 ipme(k,j) = sk_p(k,j-2,i) * cip + snenn * ( & 782 aex(ix) * cip + bex(ix) / dex(ix) * ( & 783 eex(ix) - & 784 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) ) & 785 ) & 786 ) 787 IF ( sterm == 0.0001_wp ) ipme(k,j) = sk_p(k,j-1,i) * cip 788 IF ( sterm == 0.9999_wp ) ipme(k,j) = sk_p(k,j-1,i) * cip 789 ENDIF 790 791 ENDDO 792 ENDDO 793 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' ) 794 795 ! 796 !-- Prognostic equation 797 DO j = nys, nyn 798 DO k = nzb+1, nzt 799 fplus = ( 1.0_wp - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j) & 800 - ( 1.0_wp - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j) 801 fminus = ( 1.0_wp - sw(k,j-1) ) * ipmb(k,j) + sw(k,j-1) * ipme(k,j) & 802 - ( 1.0_wp - sw(k,j) ) * immb(k,j) - sw(k,j) * imme(k,j) 803 tendcy = fplus - fminus 804 ! 805 !-- Removed in order to optimise speed 806 ! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35_wp ) 807 ! IF ( ( ABS( tendcy ) / ffmax ) < 1E-7_wp ) tendcy = 0.0 808 ! 809 !-- Density correction because of possible remaining divergences 810 d_new = d(k,j,i) - ( v(k,j+1,i) - v(k,j,i) ) * dt_3d * ddy 811 sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) / & 812 ( 1.0_wp + d_new ) 813 d(k,j,i) = d_new 814 ENDDO 815 ENDDO 816 817 ENDDO ! End of the advection in y-direction 818 CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' ) 819 CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'stop' ) 820 821 ! 822 !-- Deallocate temporary arrays 823 DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, & 824 ippb, ippe, m1, sw ) 825 826 827 ! 828 !-- Initialise for the computation of heat fluxes (see below; required in 829 !-- UP flow_statistics) 830 IF ( sk_char == 'pt' ) sums_wsts_bc_l = 0.0_wp 831 832 ! 833 !-- Add top and bottom boundaries according to the relevant boundary conditions 834 IF ( sk_char == 'pt' ) THEN 835 836 ! 837 !-- Temperature boundary condition at the bottom boundary 838 IF ( ibc_pt_b == 0 ) THEN 828 839 ! 829 840 !-- Dirichlet (fixed surface temperature) 830 DO i = nxl, nxr 831 DO j = nys, nyn 832 sk_p(nzb-1,j,i) = sk_p(nzb,j,i) 833 sk_p(nzb-2,j,i) = sk_p(nzb,j,i) 834 ENDDO 835 ENDDO 836 837 ELSE 838 ! 839 !-- Neumann (i.e. here zero gradient) 841 DO i = nxl, nxr 842 DO j = nys, nyn 843 sk_p(nzb-1,j,i) = sk_p(nzb,j,i) 844 sk_p(nzb-2,j,i) = sk_p(nzb,j,i) 845 ENDDO 846 ENDDO 847 848 ELSE 849 ! 850 !-- Neumann (i.e. here zero gradient) 851 DO i = nxl, nxr 852 DO j = nys, nyn 853 sk_p(nzb,j,i) = sk_p(nzb+1,j,i) 854 sk_p(nzb-1,j,i) = sk_p(nzb,j,i) 855 sk_p(nzb-2,j,i) = sk_p(nzb,j,i) 856 ENDDO 857 ENDDO 858 859 ENDIF 860 861 ! 862 !-- Temperature boundary condition at the top boundary 863 IF ( ibc_pt_t == 0 .OR. ibc_pt_t == 1 ) THEN 864 ! 865 !-- Dirichlet or Neumann (zero gradient) 866 DO i = nxl, nxr 867 DO j = nys, nyn 868 sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) 869 sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i) 870 ENDDO 871 ENDDO 872 873 ELSEIF ( ibc_pt_t == 2 ) THEN 874 ! 875 !-- Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead 876 DO i = nxl, nxr 877 DO j = nys, nyn 878 sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) + bc_pt_t_val * dzu(nzt+1) 879 sk_p(nzt+3,j,i) = sk_p(nzt+2,j,i) + bc_pt_t_val * dzu(nzt+1) 880 ENDDO 881 ENDDO 882 883 ENDIF 884 885 ELSEIF ( sk_char == 'sa' ) THEN 886 887 ! 888 !-- Salinity boundary condition at the bottom boundary. 889 !-- So far, always Neumann (i.e. here zero gradient) is used 840 890 DO i = nxl, nxr 841 891 DO j = nys, nyn … … 846 896 ENDDO 847 897 848 ENDIF 849 850 ! 851 !-- Temperature boundary condition at the top boundary 852 IF ( ibc_pt_t == 0 .OR. ibc_pt_t == 1 ) THEN 853 ! 898 ! 899 !-- Salinity boundary condition at the top boundary. 854 900 !-- Dirichlet or Neumann (zero gradient) 855 901 DO i = nxl, nxr … … 860 906 ENDDO 861 907 862 ELSEIF ( ibc_pt_t == 2 ) THEN 863 ! 864 !-- Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead 908 ELSEIF ( sk_char == 'q' ) THEN 909 910 ! 911 !-- Specific humidity boundary condition at the bottom boundary. 912 !-- Dirichlet (fixed surface humidity) or Neumann (i.e. zero gradient) 865 913 DO i = nxl, nxr 866 914 DO j = nys, nyn 867 sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) + bc_pt_t_val * dzu(nzt+1) 868 sk_p(nzt+3,j,i) = sk_p(nzt+2,j,i) + bc_pt_t_val * dzu(nzt+1) 869 ENDDO 870 ENDDO 871 872 ENDIF 873 874 ELSEIF ( sk_char == 'sa' ) THEN 875 876 ! 877 !-- Salinity boundary condition at the bottom boundary. 878 !-- So far, always Neumann (i.e. here zero gradient) is used 879 DO i = nxl, nxr 880 DO j = nys, nyn 881 sk_p(nzb,j,i) = sk_p(nzb+1,j,i) 882 sk_p(nzb-1,j,i) = sk_p(nzb,j,i) 883 sk_p(nzb-2,j,i) = sk_p(nzb,j,i) 884 ENDDO 885 ENDDO 886 887 ! 888 !-- Salinity boundary condition at the top boundary. 889 !-- Dirichlet or Neumann (zero gradient) 890 DO i = nxl, nxr 891 DO j = nys, nyn 892 sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) 893 sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i) 894 ENDDO 895 ENDDO 896 897 ELSEIF ( sk_char == 'q' ) THEN 898 899 ! 900 !-- Specific humidity boundary condition at the bottom boundary. 901 !-- Dirichlet (fixed surface humidity) or Neumann (i.e. zero gradient) 902 DO i = nxl, nxr 903 DO j = nys, nyn 904 sk_p(nzb,j,i) = sk_p(nzb+1,j,i) 905 sk_p(nzb-1,j,i) = sk_p(nzb,j,i) 906 sk_p(nzb-2,j,i) = sk_p(nzb,j,i) 907 ENDDO 908 ENDDO 909 910 ! 911 !-- Specific humidity boundary condition at the top boundary 912 IF ( ibc_q_t == 0 ) THEN 913 ! 914 !-- Dirichlet 915 sk_p(nzb,j,i) = sk_p(nzb+1,j,i) 916 sk_p(nzb-1,j,i) = sk_p(nzb,j,i) 917 sk_p(nzb-2,j,i) = sk_p(nzb,j,i) 918 ENDDO 919 ENDDO 920 921 ! 922 !-- Specific humidity boundary condition at the top boundary 923 IF ( ibc_q_t == 0 ) THEN 924 ! 925 !-- Dirichlet 926 DO i = nxl, nxr 927 DO j = nys, nyn 928 sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) 929 sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i) 930 ENDDO 931 ENDDO 932 933 ELSE 934 ! 935 !-- Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead 936 DO i = nxl, nxr 937 DO j = nys, nyn 938 sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) + bc_q_t_val * dzu(nzt+1) 939 sk_p(nzt+3,j,i) = sk_p(nzt+2,j,i) + bc_q_t_val * dzu(nzt+1) 940 ENDDO 941 ENDDO 942 943 ENDIF 944 945 ELSEIF ( sk_char == 'qr' ) THEN 946 947 ! 948 !-- Rain water content boundary condition at the bottom boundary: 949 !-- Dirichlet (fixed surface rain water content). 950 DO i = nxl, nxr 951 DO j = nys, nyn 952 sk_p(nzb,j,i) = sk_p(nzb+1,j,i) 953 sk_p(nzb-1,j,i) = sk_p(nzb,j,i) 954 sk_p(nzb-2,j,i) = sk_p(nzb,j,i) 955 ENDDO 956 ENDDO 957 958 ! 959 !-- Rain water content boundary condition at the top boundary: Dirichlet 915 960 DO i = nxl, nxr 916 961 DO j = nys, nyn … … 920 965 ENDDO 921 966 922 ELSE 923 ! 924 !-- Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead 967 ELSEIF ( sk_char == 'nr' ) THEN 968 969 ! 970 !-- Rain drop concentration boundary condition at the bottom boundary: 971 !-- Dirichlet (fixed surface rain drop concentration). 925 972 DO i = nxl, nxr 926 973 DO j = nys, nyn 927 sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) + bc_q_t_val * dzu(nzt+1) 928 sk_p(nzt+3,j,i) = sk_p(nzt+2,j,i) + bc_q_t_val * dzu(nzt+1) 929 ENDDO 930 ENDDO 931 974 sk_p(nzb,j,i) = sk_p(nzb+1,j,i) 975 sk_p(nzb-1,j,i) = sk_p(nzb,j,i) 976 sk_p(nzb-2,j,i) = sk_p(nzb,j,i) 977 ENDDO 978 ENDDO 979 980 ! 981 !-- Rain drop concentration boundary condition at the top boundary: Dirichlet 982 DO i = nxl, nxr 983 DO j = nys, nyn 984 sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) 985 sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i) 986 ENDDO 987 ENDDO 988 989 ELSEIF ( sk_char == 'e' ) THEN 990 991 ! 992 !-- TKE boundary condition at bottom and top boundary (generally Neumann) 993 DO i = nxl, nxr 994 DO j = nys, nyn 995 sk_p(nzb,j,i) = sk_p(nzb+1,j,i) 996 sk_p(nzb-1,j,i) = sk_p(nzb,j,i) 997 sk_p(nzb-2,j,i) = sk_p(nzb,j,i) 998 sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) 999 sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i) 1000 ENDDO 1001 ENDDO 1002 1003 ELSE 1004 1005 WRITE( message_string, * ) 'no vertical boundary condi', & 1006 'tion for variable "', sk_char, '"' 1007 CALL message( 'advec_s_bc', 'PA0158', 1, 2, 0, 6, 0 ) 1008 932 1009 ENDIF 933 1010 934 ELSEIF ( sk_char == 'qr' ) THEN 935 936 ! 937 !-- Rain water content boundary condition at the bottom boundary: 938 !-- Dirichlet (fixed surface rain water content). 1011 ! 1012 !-- Determine the maxima of the first and second derivative in z-direction 1013 fmax_l = 0.0_wp 939 1014 DO i = nxl, nxr 940 1015 DO j = nys, nyn 941 sk_p(nzb,j,i) = sk_p(nzb+1,j,i) 942 sk_p(nzb-1,j,i) = sk_p(nzb,j,i) 943 sk_p(nzb-2,j,i) = sk_p(nzb,j,i) 1016 DO k = nzb, nzt+1 1017 numera = ABS( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k-1,j,i) ) 1018 denomi = ABS( sk_p(k+1,j,i+1) - sk_p(k-1,j,i) ) 1019 fmax_l(1) = MAX( fmax_l(1) , numera ) 1020 fmax_l(2) = MAX( fmax_l(2) , denomi ) 1021 ENDDO 944 1022 ENDDO 945 1023 ENDDO 946 947 ! 948 !-- Rain water content boundary condition at the top boundary: Dirichlet 1024 #if defined( __parallel ) 1025 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1026 CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr ) 1027 #else 1028 fmax = fmax_l 1029 #endif 1030 1031 fmax = 0.04_wp * fmax 1032 1033 ! 1034 !-- Allocate temporary arrays 1035 ALLOCATE( a0(nzb:nzt+1,nys:nyn), a1(nzb:nzt+1,nys:nyn), & 1036 a2(nzb:nzt+1,nys:nyn), a12(nzb:nzt+1,nys:nyn), & 1037 a22(nzb:nzt+1,nys:nyn), immb(nzb+1:nzt,nys:nyn), & 1038 imme(nzb+1:nzt,nys:nyn), impb(nzb+1:nzt,nys:nyn), & 1039 impe(nzb+1:nzt,nys:nyn), ipmb(nzb+1:nzt,nys:nyn), & 1040 ipme(nzb+1:nzt,nys:nyn), ippb(nzb+1:nzt,nys:nyn), & 1041 ippe(nzb+1:nzt,nys:nyn), m1(nzb-1:nzt+2,nys:nyn), & 1042 sw(nzb:nzt+1,nys:nyn) & 1043 ) 1044 imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp 1045 1046 ! 1047 !-- Outer loop of all i 1048 DO i = nxl, nxr 1049 1050 ! 1051 !-- Compute polynomial coefficients 1052 DO j = nys, nyn 1053 DO k = nzb, nzt+1 1054 a12(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) ) 1055 a22(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) & 1056 + sk_p(k-1,j,i) ) 1057 a0(k,j) = ( 9.0_wp * sk_p(k+2,j,i) - 116.0_wp * sk_p(k+1,j,i) & 1058 + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k-1,j,i) & 1059 + 9.0_wp * sk_p(k-2,j,i) ) * f1920 1060 a1(k,j) = ( -5.0_wp * sk_p(k+2,j,i) + 34.0_wp * sk_p(k+1,j,i) & 1061 - 34.0_wp * sk_p(k-1,j,i) + 5.0_wp * sk_p(k-2,j,i) & 1062 ) * f48 1063 a2(k,j) = ( -3.0_wp * sk_p(k+2,j,i) + 36.0_wp * sk_p(k+1,j,i) & 1064 - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k-1,j,i) & 1065 - 3.0_wp * sk_p(k-2,j,i) ) * f48 1066 ENDDO 1067 ENDDO 1068 1069 ! 1070 !-- Fluxes using the Bott scheme 1071 !-- *VOCL LOOP,UNROLL(2) 1072 DO j = nys, nyn 1073 DO k = nzb+1, nzt 1074 cip = MAX( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) ) 1075 cim = -MIN( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) ) 1076 cipf = 1.0_wp - 2.0_wp * cip 1077 cimf = 1.0_wp - 2.0_wp * cim 1078 ip = a0(k,j) * f2 * ( 1.0_wp - cipf ) & 1079 + a1(k,j) * f8 * ( 1.0_wp - cipf*cipf ) & 1080 + a2(k,j) * f24 * ( 1.0_wp - cipf*cipf*cipf ) 1081 im = a0(k+1,j) * f2 * ( 1.0_wp - cimf ) & 1082 - a1(k+1,j) * f8 * ( 1.0_wp - cimf*cimf ) & 1083 + a2(k+1,j) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 1084 ip = MAX( ip, 0.0_wp ) 1085 im = MAX( im, 0.0_wp ) 1086 ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15_wp) ) 1087 impb(k,j) = im * MIN( 1.0_wp, sk_p(k+1,j,i) / (ip+im+1E-15_wp) ) 1088 1089 cip = MAX( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) ) 1090 cim = -MIN( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) ) 1091 cipf = 1.0_wp - 2.0_wp * cip 1092 cimf = 1.0_wp - 2.0_wp * cim 1093 ip = a0(k-1,j) * f2 * ( 1.0_wp - cipf ) & 1094 + a1(k-1,j) * f8 * ( 1.0_wp - cipf*cipf ) & 1095 + a2(k-1,j) * f24 * ( 1.0_wp - cipf*cipf*cipf ) 1096 im = a0(k,j) * f2 * ( 1.0_wp - cimf ) & 1097 - a1(k,j) * f8 * ( 1.0_wp - cimf*cimf ) & 1098 + a2(k,j) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 1099 ip = MAX( ip, 0.0_wp ) 1100 im = MAX( im, 0.0_wp ) 1101 ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k-1,j,i) / (ip+im+1E-15_wp) ) 1102 immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15_wp) ) 1103 ENDDO 1104 ENDDO 1105 1106 ! 1107 !-- Compute monitor function m1 1108 DO j = nys, nyn 1109 DO k = nzb-1, nzt+2 1110 m1z = ABS( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k-1,j,i) ) 1111 m1n = ABS( sk_p(k+1,j,i) - sk_p(k-1,j,i) ) 1112 IF ( m1n /= 0.0_wp .AND. m1n >= m1z ) THEN 1113 m1(k,j) = m1z / m1n 1114 IF ( m1(k,j) /= 2.0_wp .AND. m1n < fmax(2) ) m1(k,j) = 0.0_wp 1115 ELSEIF ( m1n < m1z ) THEN 1116 m1(k,j) = -1.0_wp 1117 ELSE 1118 m1(k,j) = 0.0_wp 1119 ENDIF 1120 ENDDO 1121 ENDDO 1122 1123 ! 1124 !-- Compute switch sw 1125 sw = 0.0_wp 1126 DO j = nys, nyn 1127 DO k = nzb, nzt+1 1128 m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) / & 1129 MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp ) 1130 IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0_wp 1131 1132 m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) / & 1133 MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp ) 1134 IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0_wp 1135 1136 t1 = 0.35_wp 1137 t2 = 0.35_wp 1138 IF ( m1(k,j) == -1.0_wp ) t2 = 0.12_wp 1139 1140 !-- *VOCL STMT,IF(10) 1141 IF ( m1(k-1,j) == 1.0_wp .OR. m1(k,j) == 1.0_wp & 1142 .OR. m1(k+1,j) == 1.0_wp .OR. m2 > t2 .OR. m3 > t2 .OR. & 1143 ( m1(k,j) > t1 .AND. m1(k-1,j) /= -1.0_wp .AND. & 1144 m1(k,j) /= -1.0_wp .AND. m1(k+1,j) /= -1.0_wp ) & 1145 ) sw(k,j) = 1.0_wp 1146 ENDDO 1147 ENDDO 1148 1149 ! 1150 !-- Fluxes using exponential scheme 1151 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' ) 1152 DO j = nys, nyn 1153 DO k = nzb+1, nzt 1154 1155 !-- *VOCL STMT,IF(10) 1156 IF ( sw(k,j) == 1.0_wp ) THEN 1157 snenn = sk_p(k+1,j,i) - sk_p(k-1,j,i) 1158 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 1159 sterm = ( sk_p(k,j,i) - sk_p(k-1,j,i) ) / snenn 1160 sterm = MIN( sterm, 0.9999_wp ) 1161 sterm = MAX( sterm, 0.0001_wp ) 1162 1163 ix = INT( sterm * 1000 ) + 1 1164 1165 cip = MAX( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) ) 1166 1167 ippe(k,j) = sk_p(k-1,j,i) * cip + snenn * ( & 1168 aex(ix) * cip + bex(ix) / dex(ix) * ( & 1169 eex(ix) - & 1170 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) ) & 1171 ) & 1172 ) 1173 IF ( sterm == 0.0001_wp ) ippe(k,j) = sk_p(k,j,i) * cip 1174 IF ( sterm == 0.9999_wp ) ippe(k,j) = sk_p(k,j,i) * cip 1175 1176 snenn = sk_p(k-1,j,i) - sk_p(k+1,j,i) 1177 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 1178 sterm = ( sk_p(k,j,i) - sk_p(k+1,j,i) ) / snenn 1179 sterm = MIN( sterm, 0.9999_wp ) 1180 sterm = MAX( sterm, 0.0001_wp ) 1181 1182 ix = INT( sterm * 1000 ) + 1 1183 1184 cim = -MIN( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) ) 1185 1186 imme(k,j) = sk_p(k+1,j,i) * cim + snenn * ( & 1187 aex(ix) * cim + bex(ix) / dex(ix) * ( & 1188 eex(ix) - & 1189 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) & 1190 ) & 1191 ) 1192 IF ( sterm == 0.0001_wp ) imme(k,j) = sk_p(k,j,i) * cim 1193 IF ( sterm == 0.9999_wp ) imme(k,j) = sk_p(k,j,i) * cim 1194 ENDIF 1195 1196 !-- *VOCL STMT,IF(10) 1197 IF ( sw(k+1,j) == 1.0_wp ) THEN 1198 snenn = sk_p(k,j,i) - sk_p(k+2,j,i) 1199 IF ( ABS( snenn ) .LT. 1E-9_wp ) snenn = 1E-9_wp 1200 sterm = ( sk_p(k+1,j,i) - sk_p(k+2,j,i) ) / snenn 1201 sterm = MIN( sterm, 0.9999_wp ) 1202 sterm = MAX( sterm, 0.0001_wp ) 1203 1204 ix = INT( sterm * 1000 ) + 1 1205 1206 cim = -MIN( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) ) 1207 1208 impe(k,j) = sk_p(k+2,j,i) * cim + snenn * ( & 1209 aex(ix) * cim + bex(ix) / dex(ix) * ( & 1210 eex(ix) - & 1211 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) & 1212 ) & 1213 ) 1214 IF ( sterm == 0.0001_wp ) impe(k,j) = sk_p(k+1,j,i) * cim 1215 IF ( sterm == 0.9999_wp ) impe(k,j) = sk_p(k+1,j,i) * cim 1216 ENDIF 1217 1218 !-- *VOCL STMT,IF(10) 1219 IF ( sw(k-1,j) == 1.0_wp ) THEN 1220 snenn = sk_p(k,j,i) - sk_p(k-2,j,i) 1221 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 1222 sterm = ( sk_p(k-1,j,i) - sk_p(k-2,j,i) ) / snenn 1223 sterm = MIN( sterm, 0.9999_wp ) 1224 sterm = MAX( sterm, 0.0001_wp ) 1225 1226 ix = INT( sterm * 1000 ) + 1 1227 1228 cip = MAX( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) ) 1229 1230 ipme(k,j) = sk_p(k-2,j,i) * cip + snenn * ( & 1231 aex(ix) * cip + bex(ix) / dex(ix) * ( & 1232 eex(ix) - & 1233 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) ) & 1234 ) & 1235 ) 1236 IF ( sterm == 0.0001_wp ) ipme(k,j) = sk_p(k-1,j,i) * cip 1237 IF ( sterm == 0.9999_wp ) ipme(k,j) = sk_p(k-1,j,i) * cip 1238 ENDIF 1239 1240 ENDDO 1241 ENDDO 1242 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' ) 1243 1244 ! 1245 !-- Prognostic equation 1246 DO j = nys, nyn 1247 DO k = nzb+1, nzt 1248 fplus = ( 1.0_wp - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j) & 1249 - ( 1.0_wp - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j) 1250 fminus = ( 1.0_wp - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j) & 1251 - ( 1.0_wp - sw(k,j) ) * immb(k,j) - sw(k,j) * imme(k,j) 1252 tendcy = fplus - fminus 1253 ! 1254 !-- Removed in order to optimise speed 1255 ! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35_wp ) 1256 ! IF ( ( ABS( tendcy ) / ffmax ) < 1E-7_wp ) tendcy = 0.0 1257 ! 1258 !-- Density correction because of possible remaining divergences 1259 d_new = d(k,j,i) - ( w(k,j,i) - w(k-1,j,i) ) * dt_3d * ddzw(k) 1260 sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) / & 1261 ( 1.0_wp + d_new ) 1262 ! 1263 !-- Store heat flux for subsequent statistics output. 1264 !-- array m1 is here used as temporary storage 1265 m1(k,j) = fplus / dt_3d * dzw(k) 1266 ENDDO 1267 ENDDO 1268 1269 ! 1270 !-- Sum up heat flux in order to order to obtain horizontal averages 1271 IF ( sk_char == 'pt' ) THEN 1272 DO sr = 0, statistic_regions 1273 DO j = nys, nyn 1274 DO k = nzb+1, nzt 1275 sums_wsts_bc_l(k,sr) = sums_wsts_bc_l(k,sr) + & 1276 m1(k,j) * rmask(j,i,sr) 1277 ENDDO 1278 ENDDO 1279 ENDDO 1280 ENDIF 1281 1282 ENDDO ! End of the advection in z-direction 1283 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' ) 1284 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'stop' ) 1285 1286 ! 1287 !-- Deallocate temporary arrays 1288 DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, & 1289 ippb, ippe, m1, sw ) 1290 1291 ! 1292 !-- Store results as tendency and deallocate local array 949 1293 DO i = nxl, nxr 950 1294 DO j = nys, nyn 951 sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) 952 sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i) 1295 DO k = nzb+1, nzt 1296 tend(k,j,i) = tend(k,j,i) + ( sk_p(k,j,i) - sk(k,j,i) ) / dt_3d 1297 ENDDO 953 1298 ENDDO 954 1299 ENDDO 955 1300 956 ELSEIF ( sk_char == 'nr' ) THEN 957 958 ! 959 !-- Rain drop concentration boundary condition at the bottom boundary: 960 !-- Dirichlet (fixed surface rain drop concentration). 961 DO i = nxl, nxr 962 DO j = nys, nyn 963 sk_p(nzb,j,i) = sk_p(nzb+1,j,i) 964 sk_p(nzb-1,j,i) = sk_p(nzb,j,i) 965 sk_p(nzb-2,j,i) = sk_p(nzb,j,i) 966 ENDDO 967 ENDDO 968 969 ! 970 !-- Rain drop concentration boundary condition at the top boundary: Dirichlet 971 DO i = nxl, nxr 972 DO j = nys, nyn 973 sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) 974 sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i) 975 ENDDO 976 ENDDO 977 978 ELSEIF ( sk_char == 'e' ) THEN 979 980 ! 981 !-- TKE boundary condition at bottom and top boundary (generally Neumann) 982 DO i = nxl, nxr 983 DO j = nys, nyn 984 sk_p(nzb,j,i) = sk_p(nzb+1,j,i) 985 sk_p(nzb-1,j,i) = sk_p(nzb,j,i) 986 sk_p(nzb-2,j,i) = sk_p(nzb,j,i) 987 sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) 988 sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i) 989 ENDDO 990 ENDDO 991 992 ELSE 993 994 WRITE( message_string, * ) 'no vertical boundary condi', & 995 'tion for variable "', sk_char, '"' 996 CALL message( 'advec_s_bc', 'PA0158', 1, 2, 0, 6, 0 ) 997 998 ENDIF 999 1000 ! 1001 !-- Determine the maxima of the first and second derivative in z-direction 1002 fmax_l = 0.0_wp 1003 DO i = nxl, nxr 1004 DO j = nys, nyn 1005 DO k = nzb, nzt+1 1006 numera = ABS( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k-1,j,i) ) 1007 denomi = ABS( sk_p(k+1,j,i+1) - sk_p(k-1,j,i) ) 1008 fmax_l(1) = MAX( fmax_l(1) , numera ) 1009 fmax_l(2) = MAX( fmax_l(2) , denomi ) 1010 ENDDO 1011 ENDDO 1012 ENDDO 1013 #if defined( __parallel ) 1014 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1015 CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr ) 1016 #else 1017 fmax = fmax_l 1018 #endif 1019 1020 fmax = 0.04_wp * fmax 1021 1022 ! 1023 !-- Allocate temporary arrays 1024 ALLOCATE( a0(nzb:nzt+1,nys:nyn), a1(nzb:nzt+1,nys:nyn), & 1025 a2(nzb:nzt+1,nys:nyn), a12(nzb:nzt+1,nys:nyn), & 1026 a22(nzb:nzt+1,nys:nyn), immb(nzb+1:nzt,nys:nyn), & 1027 imme(nzb+1:nzt,nys:nyn), impb(nzb+1:nzt,nys:nyn), & 1028 impe(nzb+1:nzt,nys:nyn), ipmb(nzb+1:nzt,nys:nyn), & 1029 ipme(nzb+1:nzt,nys:nyn), ippb(nzb+1:nzt,nys:nyn), & 1030 ippe(nzb+1:nzt,nys:nyn), m1(nzb-1:nzt+2,nys:nyn), & 1031 sw(nzb:nzt+1,nys:nyn) & 1032 ) 1033 imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp 1034 1035 ! 1036 !-- Outer loop of all i 1037 DO i = nxl, nxr 1038 1039 ! 1040 !-- Compute polynomial coefficients 1041 DO j = nys, nyn 1042 DO k = nzb, nzt+1 1043 a12(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) ) 1044 a22(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) & 1045 + sk_p(k-1,j,i) ) 1046 a0(k,j) = ( 9.0_wp * sk_p(k+2,j,i) - 116.0_wp * sk_p(k+1,j,i) & 1047 + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k-1,j,i) & 1048 + 9.0_wp * sk_p(k-2,j,i) ) * f1920 1049 a1(k,j) = ( -5.0_wp * sk_p(k+2,j,i) + 34.0_wp * sk_p(k+1,j,i) & 1050 - 34.0_wp * sk_p(k-1,j,i) + 5.0_wp * sk_p(k-2,j,i) & 1051 ) * f48 1052 a2(k,j) = ( -3.0_wp * sk_p(k+2,j,i) + 36.0_wp * sk_p(k+1,j,i) & 1053 - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k-1,j,i) & 1054 - 3.0_wp * sk_p(k-2,j,i) ) * f48 1055 ENDDO 1056 ENDDO 1057 1058 ! 1059 !-- Fluxes using the Bott scheme 1060 !-- *VOCL LOOP,UNROLL(2) 1061 DO j = nys, nyn 1062 DO k = nzb+1, nzt 1063 cip = MAX( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) ) 1064 cim = -MIN( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) ) 1065 cipf = 1.0_wp - 2.0_wp * cip 1066 cimf = 1.0_wp - 2.0_wp * cim 1067 ip = a0(k,j) * f2 * ( 1.0_wp - cipf ) & 1068 + a1(k,j) * f8 * ( 1.0_wp - cipf*cipf ) & 1069 + a2(k,j) * f24 * ( 1.0_wp - cipf*cipf*cipf ) 1070 im = a0(k+1,j) * f2 * ( 1.0_wp - cimf ) & 1071 - a1(k+1,j) * f8 * ( 1.0_wp - cimf*cimf ) & 1072 + a2(k+1,j) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 1073 ip = MAX( ip, 0.0_wp ) 1074 im = MAX( im, 0.0_wp ) 1075 ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15_wp) ) 1076 impb(k,j) = im * MIN( 1.0_wp, sk_p(k+1,j,i) / (ip+im+1E-15_wp) ) 1077 1078 cip = MAX( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) ) 1079 cim = -MIN( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) ) 1080 cipf = 1.0_wp - 2.0_wp * cip 1081 cimf = 1.0_wp - 2.0_wp * cim 1082 ip = a0(k-1,j) * f2 * ( 1.0_wp - cipf ) & 1083 + a1(k-1,j) * f8 * ( 1.0_wp - cipf*cipf ) & 1084 + a2(k-1,j) * f24 * ( 1.0_wp - cipf*cipf*cipf ) 1085 im = a0(k,j) * f2 * ( 1.0_wp - cimf ) & 1086 - a1(k,j) * f8 * ( 1.0_wp - cimf*cimf ) & 1087 + a2(k,j) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 1088 ip = MAX( ip, 0.0_wp ) 1089 im = MAX( im, 0.0_wp ) 1090 ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k-1,j,i) / (ip+im+1E-15_wp) ) 1091 immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15_wp) ) 1092 ENDDO 1093 ENDDO 1094 1095 ! 1096 !-- Compute monitor function m1 1097 DO j = nys, nyn 1098 DO k = nzb-1, nzt+2 1099 m1z = ABS( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k-1,j,i) ) 1100 m1n = ABS( sk_p(k+1,j,i) - sk_p(k-1,j,i) ) 1101 IF ( m1n /= 0.0_wp .AND. m1n >= m1z ) THEN 1102 m1(k,j) = m1z / m1n 1103 IF ( m1(k,j) /= 2.0_wp .AND. m1n < fmax(2) ) m1(k,j) = 0.0_wp 1104 ELSEIF ( m1n < m1z ) THEN 1105 m1(k,j) = -1.0_wp 1106 ELSE 1107 m1(k,j) = 0.0_wp 1108 ENDIF 1109 ENDDO 1110 ENDDO 1111 1112 ! 1113 !-- Compute switch sw 1114 sw = 0.0_wp 1115 DO j = nys, nyn 1116 DO k = nzb, nzt+1 1117 m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) / & 1118 MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp ) 1119 IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0_wp 1120 1121 m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) / & 1122 MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp ) 1123 IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0_wp 1124 1125 t1 = 0.35_wp 1126 t2 = 0.35_wp 1127 IF ( m1(k,j) == -1.0_wp ) t2 = 0.12_wp 1128 1129 !-- *VOCL STMT,IF(10) 1130 IF ( m1(k-1,j) == 1.0_wp .OR. m1(k,j) == 1.0_wp & 1131 .OR. m1(k+1,j) == 1.0_wp .OR. m2 > t2 .OR. m3 > t2 .OR. & 1132 ( m1(k,j) > t1 .AND. m1(k-1,j) /= -1.0_wp .AND. & 1133 m1(k,j) /= -1.0_wp .AND. m1(k+1,j) /= -1.0_wp ) & 1134 ) sw(k,j) = 1.0_wp 1135 ENDDO 1136 ENDDO 1137 1138 ! 1139 !-- Fluxes using exponential scheme 1140 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' ) 1141 DO j = nys, nyn 1142 DO k = nzb+1, nzt 1143 1144 !-- *VOCL STMT,IF(10) 1145 IF ( sw(k,j) == 1.0_wp ) THEN 1146 snenn = sk_p(k+1,j,i) - sk_p(k-1,j,i) 1147 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 1148 sterm = ( sk_p(k,j,i) - sk_p(k-1,j,i) ) / snenn 1149 sterm = MIN( sterm, 0.9999_wp ) 1150 sterm = MAX( sterm, 0.0001_wp ) 1151 1152 ix = INT( sterm * 1000 ) + 1 1153 1154 cip = MAX( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) ) 1155 1156 ippe(k,j) = sk_p(k-1,j,i) * cip + snenn * ( & 1157 aex(ix) * cip + bex(ix) / dex(ix) * ( & 1158 eex(ix) - & 1159 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) ) & 1160 ) & 1161 ) 1162 IF ( sterm == 0.0001_wp ) ippe(k,j) = sk_p(k,j,i) * cip 1163 IF ( sterm == 0.9999_wp ) ippe(k,j) = sk_p(k,j,i) * cip 1164 1165 snenn = sk_p(k-1,j,i) - sk_p(k+1,j,i) 1166 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 1167 sterm = ( sk_p(k,j,i) - sk_p(k+1,j,i) ) / snenn 1168 sterm = MIN( sterm, 0.9999_wp ) 1169 sterm = MAX( sterm, 0.0001_wp ) 1170 1171 ix = INT( sterm * 1000 ) + 1 1172 1173 cim = -MIN( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) ) 1174 1175 imme(k,j) = sk_p(k+1,j,i) * cim + snenn * ( & 1176 aex(ix) * cim + bex(ix) / dex(ix) * ( & 1177 eex(ix) - & 1178 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) & 1179 ) & 1180 ) 1181 IF ( sterm == 0.0001_wp ) imme(k,j) = sk_p(k,j,i) * cim 1182 IF ( sterm == 0.9999_wp ) imme(k,j) = sk_p(k,j,i) * cim 1183 ENDIF 1184 1185 !-- *VOCL STMT,IF(10) 1186 IF ( sw(k+1,j) == 1.0_wp ) THEN 1187 snenn = sk_p(k,j,i) - sk_p(k+2,j,i) 1188 IF ( ABS( snenn ) .LT. 1E-9_wp ) snenn = 1E-9_wp 1189 sterm = ( sk_p(k+1,j,i) - sk_p(k+2,j,i) ) / snenn 1190 sterm = MIN( sterm, 0.9999_wp ) 1191 sterm = MAX( sterm, 0.0001_wp ) 1192 1193 ix = INT( sterm * 1000 ) + 1 1194 1195 cim = -MIN( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) ) 1196 1197 impe(k,j) = sk_p(k+2,j,i) * cim + snenn * ( & 1198 aex(ix) * cim + bex(ix) / dex(ix) * ( & 1199 eex(ix) - & 1200 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) & 1201 ) & 1202 ) 1203 IF ( sterm == 0.0001_wp ) impe(k,j) = sk_p(k+1,j,i) * cim 1204 IF ( sterm == 0.9999_wp ) impe(k,j) = sk_p(k+1,j,i) * cim 1205 ENDIF 1206 1207 !-- *VOCL STMT,IF(10) 1208 IF ( sw(k-1,j) == 1.0_wp ) THEN 1209 snenn = sk_p(k,j,i) - sk_p(k-2,j,i) 1210 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 1211 sterm = ( sk_p(k-1,j,i) - sk_p(k-2,j,i) ) / snenn 1212 sterm = MIN( sterm, 0.9999_wp ) 1213 sterm = MAX( sterm, 0.0001_wp ) 1214 1215 ix = INT( sterm * 1000 ) + 1 1216 1217 cip = MAX( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) ) 1218 1219 ipme(k,j) = sk_p(k-2,j,i) * cip + snenn * ( & 1220 aex(ix) * cip + bex(ix) / dex(ix) * ( & 1221 eex(ix) - & 1222 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) ) & 1223 ) & 1224 ) 1225 IF ( sterm == 0.0001_wp ) ipme(k,j) = sk_p(k-1,j,i) * cip 1226 IF ( sterm == 0.9999_wp ) ipme(k,j) = sk_p(k-1,j,i) * cip 1227 ENDIF 1228 1229 ENDDO 1230 ENDDO 1231 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' ) 1232 1233 ! 1234 !-- Prognostic equation 1235 DO j = nys, nyn 1236 DO k = nzb+1, nzt 1237 fplus = ( 1.0_wp - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j) & 1238 - ( 1.0_wp - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j) 1239 fminus = ( 1.0_wp - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j) & 1240 - ( 1.0_wp - sw(k,j) ) * immb(k,j) - sw(k,j) * imme(k,j) 1241 tendcy = fplus - fminus 1242 ! 1243 !-- Removed in order to optimise speed 1244 ! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35_wp ) 1245 ! IF ( ( ABS( tendcy ) / ffmax ) < 1E-7_wp ) tendcy = 0.0 1246 ! 1247 !-- Density correction because of possible remaining divergences 1248 d_new = d(k,j,i) - ( w(k,j,i) - w(k-1,j,i) ) * dt_3d * ddzw(k) 1249 sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) / & 1250 ( 1.0_wp + d_new ) 1251 ! 1252 !-- Store heat flux for subsequent statistics output. 1253 !-- array m1 is here used as temporary storage 1254 m1(k,j) = fplus / dt_3d * dzw(k) 1255 ENDDO 1256 ENDDO 1257 1258 ! 1259 !-- Sum up heat flux in order to order to obtain horizontal averages 1260 IF ( sk_char == 'pt' ) THEN 1261 DO sr = 0, statistic_regions 1262 DO j = nys, nyn 1263 DO k = nzb+1, nzt 1264 sums_wsts_bc_l(k,sr) = sums_wsts_bc_l(k,sr) + & 1265 m1(k,j) * rmask(j,i,sr) 1266 ENDDO 1267 ENDDO 1268 ENDDO 1269 ENDIF 1270 1271 ENDDO ! End of the advection in z-direction 1272 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' ) 1273 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'stop' ) 1274 1275 ! 1276 !-- Deallocate temporary arrays 1277 DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, & 1278 ippb, ippe, m1, sw ) 1279 1280 ! 1281 !-- Store results as tendency and deallocate local array 1282 DO i = nxl, nxr 1283 DO j = nys, nyn 1284 DO k = nzb+1, nzt 1285 tend(k,j,i) = tend(k,j,i) + ( sk_p(k,j,i) - sk(k,j,i) ) / dt_3d 1286 ENDDO 1287 ENDDO 1288 ENDDO 1289 1290 DEALLOCATE( sk_p ) 1291 1292 END SUBROUTINE advec_s_bc 1301 DEALLOCATE( sk_p ) 1302 1303 END SUBROUTINE advec_s_bc 1304 1305 END MODULE advec_s_bc_mod
Note: See TracChangeset
for help on using the changeset viewer.
