Changeset 1353 for palm/trunk/SOURCE/advec_s_bc.f90
- Timestamp:
- Apr 8, 2014 3:21:23 PM (11 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/advec_s_bc.f90 (modified) (49 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_s_bc.f90
r1347 r1353 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! REAL constants provided with KIND-attribute 23 23 ! 24 24 ! Former revisions: … … 170 170 !-- Array sk_p requires 2 extra elements for each dimension 171 171 ALLOCATE( sk_p(nzb-2:nzt+3,nys-3:nyn+3,nxl-3:nxr+3) ) 172 sk_p = 0.0 172 sk_p = 0.0_wp 173 173 174 174 ! 175 175 !-- Assign reciprocal values in order to avoid divisions later 176 f2 = 0.5 177 f4 = 0.25 178 f8 = 0.125 179 f12 = 0.8333333333333333E-01 180 f24 = 0.4166666666666666E-01 181 f48 = 0.2083333333333333E-01 182 f1920 = 0.5208333333333333E-03 176 f2 = 0.5_wp 177 f4 = 0.25_wp 178 f8 = 0.125_wp 179 f12 = 0.8333333333333333E-01_wp 180 f24 = 0.4166666666666666E-01_wp 181 f48 = 0.2083333333333333E-01_wp 182 f1920 = 0.5208333333333333E-03_wp 183 183 184 184 ! … … 236 236 ! 237 237 !-- Initialise control density 238 d = 0.0 238 d = 0.0_wp 239 239 240 240 ! 241 241 !-- Determine maxima of the first and second derivative in x-direction 242 fmax_l = 0.0 242 fmax_l = 0.0_wp 243 243 DO i = nxl, nxr 244 244 DO j = nys, nyn 245 245 DO k = nzb+1, nzt 246 numera = ABS( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) + sk_p(k,j,i-1) )246 numera = ABS( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) ) 247 247 denomi = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) ) 248 248 fmax_l(1) = MAX( fmax_l(1) , numera ) … … 258 258 #endif 259 259 260 fmax = 0.04 * fmax260 fmax = 0.04_wp * fmax 261 261 262 262 ! … … 271 271 sw(nzb+1:nzt,nxl-1:nxr+1) & 272 272 ) 273 imme = 0.0 ; impe = 0.0; ipme = 0.0; ippe = 0.0273 imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp 274 274 275 275 ! … … 287 287 DO i = nxl-1, nxr+1 288 288 DO k = nzb+1, nzt 289 a12(k,i) = 0.5 * ( sk_p(k,j,i+1) - sk_p(k,j,i-1) )290 a22(k,i) = 0.5 * ( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i)&291 + sk_p(k,j,i-1) )292 a0(k,i) = ( 9.0 * sk_p(k,j,i+2) - 116.0 * sk_p(k,j,i+1)&293 + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k,j,i-1)&294 + 9.0 * sk_p(k,j,i-2) ) * f1920295 a1(k,i) = ( -5.0 * sk_p(k,j,i+2) + 34.0 * sk_p(k,j,i+1)&296 - 34.0 * sk_p(k,j,i-1) + 5.0 * sk_p(k,j,i-2)&289 a12(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - sk_p(k,j,i-1) ) 290 a22(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) & 291 + sk_p(k,j,i-1) ) 292 a0(k,i) = ( 9.0_wp * sk_p(k,j,i+2) - 116.0_wp * sk_p(k,j,i+1) & 293 + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j,i-1) & 294 + 9.0_wp * sk_p(k,j,i-2) ) * f1920 295 a1(k,i) = ( -5.0_wp * sk_p(k,j,i+2) + 34.0_wp * sk_p(k,j,i+1) & 296 - 34.0_wp * sk_p(k,j,i-1) + 5.0_wp * sk_p(k,j,i-2) & 297 297 ) * f48 298 a2(k,i) = ( -3.0 * sk_p(k,j,i+2) + 36.0 * sk_p(k,j,i+1)&299 - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k,j,i-1)&300 - 3.0 * sk_p(k,j,i-2) ) * f48298 a2(k,i) = ( -3.0_wp * sk_p(k,j,i+2) + 36.0_wp * sk_p(k,j,i+1) & 299 - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j,i-1) & 300 - 3.0_wp * sk_p(k,j,i-2) ) * f48 301 301 ENDDO 302 302 ENDDO … … 309 309 cip = MAX( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx ) 310 310 cim = -MIN( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx ) 311 cipf = 1.0 - 2.0* cip312 cimf = 1.0 - 2.0* cim313 ip = a0(k,i) * f2 * ( 1.0 - cipf )&314 + a1(k,i) * f8 * ( 1.0 - cipf*cipf )&315 + a2(k,i) * f24 * ( 1.0 - cipf*cipf*cipf )316 im = a0(k,i+1) * f2 * ( 1.0 - cimf )&317 - a1(k,i+1) * f8 * ( 1.0 - cimf*cimf )&318 + a2(k,i+1) * f24 * ( 1.0 - cimf*cimf*cimf )311 cipf = 1.0_wp - 2.0_wp * cip 312 cimf = 1.0_wp - 2.0_wp * cim 313 ip = a0(k,i) * f2 * ( 1.0_wp - cipf ) & 314 + a1(k,i) * f8 * ( 1.0_wp - cipf*cipf ) & 315 + a2(k,i) * f24 * ( 1.0_wp - cipf*cipf*cipf ) 316 im = a0(k,i+1) * f2 * ( 1.0_wp - cimf ) & 317 - a1(k,i+1) * f8 * ( 1.0_wp - cimf*cimf ) & 318 + a2(k,i+1) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 319 319 ip = MAX( ip, 0.0_wp ) 320 320 im = MAX( im, 0.0_wp ) 321 ippb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15 ) )322 impb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i+1) / (ip+im+1E-15 ) )321 ippb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15_wp) ) 322 impb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i+1) / (ip+im+1E-15_wp) ) 323 323 324 324 cip = MAX( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx ) 325 325 cim = -MIN( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx ) 326 cipf = 1.0 - 2.0* cip327 cimf = 1.0 - 2.0* cim328 ip = a0(k,i-1) * f2 * ( 1.0 - cipf )&329 + a1(k,i-1) * f8 * ( 1.0 - cipf*cipf )&330 + a2(k,i-1) * f24 * ( 1.0 - cipf*cipf*cipf )331 im = a0(k,i) * f2 * ( 1.0 - cimf )&332 - a1(k,i) * f8 * ( 1.0 - cimf*cimf )&333 + a2(k,i) * f24 * ( 1.0 - cimf*cimf*cimf )326 cipf = 1.0_wp - 2.0_wp * cip 327 cimf = 1.0_wp - 2.0_wp * cim 328 ip = a0(k,i-1) * f2 * ( 1.0_wp - cipf ) & 329 + a1(k,i-1) * f8 * ( 1.0_wp - cipf*cipf ) & 330 + a2(k,i-1) * f24 * ( 1.0_wp - cipf*cipf*cipf ) 331 im = a0(k,i) * f2 * ( 1.0_wp - cimf ) & 332 - a1(k,i) * f8 * ( 1.0_wp - cimf*cimf ) & 333 + a2(k,i) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 334 334 ip = MAX( ip, 0.0_wp ) 335 335 im = MAX( im, 0.0_wp ) 336 ipmb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i-1) / (ip+im+1E-15 ) )337 immb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15 ) )336 ipmb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i-1) / (ip+im+1E-15_wp) ) 337 immb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15_wp) ) 338 338 ENDDO 339 339 ENDDO … … 343 343 DO i = nxl-2, nxr+2 344 344 DO k = nzb+1, nzt 345 m1z = ABS( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) + sk_p(k,j,i-1) )345 m1z = ABS( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) ) 346 346 m1n = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) ) 347 IF ( m1n /= 0.0 .AND. m1n >= m1z ) THEN347 IF ( m1n /= 0.0_wp .AND. m1n >= m1z ) THEN 348 348 m1(k,i) = m1z / m1n 349 IF ( m1(k,i) /= 2.0 .AND. m1n < fmax(2) ) m1(k,i) = 0.0349 IF ( m1(k,i) /= 2.0_wp .AND. m1n < fmax(2) ) m1(k,i) = 0.0_wp 350 350 ELSEIF ( m1n < m1z ) THEN 351 m1(k,i) = -1.0 351 m1(k,i) = -1.0_wp 352 352 ELSE 353 m1(k,i) = 0.0 353 m1(k,i) = 0.0_wp 354 354 ENDIF 355 355 ENDDO … … 358 358 ! 359 359 !-- Compute switch sw 360 sw = 0.0 360 sw = 0.0_wp 361 361 DO i = nxl-1, nxr+1 362 362 DO k = nzb+1, nzt 363 m2 = 2.0 * ABS( a1(k,i) - a12(k,i) ) /&363 m2 = 2.0_wp * ABS( a1(k,i) - a12(k,i) ) / & 364 364 MAX( ABS( a1(k,i) + a12(k,i) ), 1E-35_wp ) 365 IF ( ABS( a1(k,i) + a12(k,i) ) < fmax(2) ) m2 = 0.0 366 367 m3 = 2.0 * ABS( a2(k,i) - a22(k,i) ) /&365 IF ( ABS( a1(k,i) + a12(k,i) ) < fmax(2) ) m2 = 0.0_wp 366 367 m3 = 2.0_wp * ABS( a2(k,i) - a22(k,i) ) / & 368 368 MAX( ABS( a2(k,i) + a22(k,i) ), 1E-35_wp ) 369 IF ( ABS( a2(k,i) + a22(k,i) ) < fmax(1) ) m3 = 0.0 370 371 t1 = 0.35 372 t2 = 0.35 373 IF ( m1(k,i) == -1.0 ) t2 = 0.12369 IF ( ABS( a2(k,i) + a22(k,i) ) < fmax(1) ) m3 = 0.0_wp 370 371 t1 = 0.35_wp 372 t2 = 0.35_wp 373 IF ( m1(k,i) == -1.0_wp ) t2 = 0.12_wp 374 374 375 375 !-- *VOCL STMT,IF(10) 376 IF ( m1(k,i-1) == 1.0 .OR. m1(k,i) == 1.0 .OR. m1(k,i+1) == 1.0&377 .OR. m2 > t2 .OR. m3 > t2 .OR.&378 ( m1(k,i) > t1 .AND. m1(k,i-1) /= -1.0 .AND.&379 m1(k,i) /= -1.0 .AND. m1(k,i+1) /= -1.0 )&380 ) sw(k,i) = 1.0 376 IF ( m1(k,i-1) == 1.0_wp .OR. m1(k,i) == 1.0_wp & 377 .OR. m1(k,i+1) == 1.0_wp .OR. m2 > t2 .OR. m3 > t2 .OR. & 378 ( m1(k,i) > t1 .AND. m1(k,i-1) /= -1.0_wp .AND. & 379 m1(k,i) /= -1.0_wp .AND. m1(k,i+1) /= -1.0_wp ) & 380 ) sw(k,i) = 1.0_wp 381 381 ENDDO 382 382 ENDDO … … 389 389 390 390 !-- *VOCL STMT,IF(10) 391 IF ( sw(k,i) == 1.0 ) THEN391 IF ( sw(k,i) == 1.0_wp ) THEN 392 392 snenn = sk_p(k,j,i+1) - sk_p(k,j,i-1) 393 IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9393 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 394 394 sterm = ( sk_p(k,j,i) - sk_p(k,j,i-1) ) / snenn 395 395 sterm = MIN( sterm, 0.9999_wp ) … … 402 402 ippe(k,i) = sk_p(k,j,i-1) * cip + snenn * ( & 403 403 aex(ix) * cip + bex(ix) / dex(ix) * ( & 404 eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) & 404 eex(ix) - & 405 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) ) & 405 406 ) & 406 407 ) 407 IF ( sterm == 0.0001 ) ippe(k,i) = sk_p(k,j,i) * cip408 IF ( sterm == 0.9999 ) ippe(k,i) = sk_p(k,j,i) * cip408 IF ( sterm == 0.0001_wp ) ippe(k,i) = sk_p(k,j,i) * cip 409 IF ( sterm == 0.9999_wp ) ippe(k,i) = sk_p(k,j,i) * cip 409 410 410 411 snenn = sk_p(k,j,i-1) - sk_p(k,j,i+1) 411 IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9412 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 412 413 sterm = ( sk_p(k,j,i) - sk_p(k,j,i+1) ) / snenn 413 414 sterm = MIN( sterm, 0.9999_wp ) … … 420 421 imme(k,i) = sk_p(k,j,i+1) * cim + snenn * ( & 421 422 aex(ix) * cim + bex(ix) / dex(ix) * ( & 422 eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) & 423 eex(ix) - & 424 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) & 423 425 ) & 424 426 ) 425 IF ( sterm == 0.0001 ) imme(k,i) = sk_p(k,j,i) * cim426 IF ( sterm == 0.9999 ) imme(k,i) = sk_p(k,j,i) * cim427 IF ( sterm == 0.0001_wp ) imme(k,i) = sk_p(k,j,i) * cim 428 IF ( sterm == 0.9999_wp ) imme(k,i) = sk_p(k,j,i) * cim 427 429 ENDIF 428 430 429 431 !-- *VOCL STMT,IF(10) 430 IF ( sw(k,i+1) == 1.0 ) THEN432 IF ( sw(k,i+1) == 1.0_wp ) THEN 431 433 snenn = sk_p(k,j,i) - sk_p(k,j,i+2) 432 IF ( ABS( snenn ) .LT. 1E-9 ) snenn = 1E-9434 IF ( ABS( snenn ) .LT. 1E-9_wp ) snenn = 1E-9_wp 433 435 sterm = ( sk_p(k,j,i+1) - sk_p(k,j,i+2) ) / snenn 434 436 sterm = MIN( sterm, 0.9999_wp ) … … 441 443 impe(k,i) = sk_p(k,j,i+2) * cim + snenn * ( & 442 444 aex(ix) * cim + bex(ix) / dex(ix) * ( & 443 eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) & 445 eex(ix) - & 446 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) & 444 447 ) & 445 448 ) … … 449 452 450 453 !-- *VOCL STMT,IF(10) 451 IF ( sw(k,i-1) == 1.0 ) THEN454 IF ( sw(k,i-1) == 1.0_wp ) THEN 452 455 snenn = sk_p(k,j,i) - sk_p(k,j,i-2) 453 IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9456 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 454 457 sterm = ( sk_p(k,j,i-1) - sk_p(k,j,i-2) ) / snenn 455 458 sterm = MIN( sterm, 0.9999_wp ) … … 462 465 ipme(k,i) = sk_p(k,j,i-2) * cip + snenn * ( & 463 466 aex(ix) * cip + bex(ix) / dex(ix) * ( & 464 eex(ix) - EXP( dex(ix)*0.5 * ( 1.0_wp - 2.0 * cip ) ) & 467 eex(ix) - & 468 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) ) & 465 469 ) & 466 470 ) … … 538 542 ! 539 543 !-- Determine the maxima of the first and second derivative in y-direction 540 fmax_l = 0.0 544 fmax_l = 0.0_wp 541 545 DO i = nxl, nxr 542 546 DO j = nys, nyn 543 547 DO k = nzb+1, nzt 544 numera = ABS( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) + sk_p(k,j-1,i) )548 numera = ABS( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j-1,i) ) 545 549 denomi = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) ) 546 550 fmax_l(1) = MAX( fmax_l(1) , numera ) … … 556 560 #endif 557 561 558 fmax = 0.04 * fmax562 fmax = 0.04_wp * fmax 559 563 560 564 ! … … 569 573 sw(nzb+1:nzt,nys-1:nyn+1) & 570 574 ) 571 imme = 0.0 ; impe = 0.0; ipme = 0.0; ippe = 0.0575 imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp 572 576 573 577 ! … … 579 583 DO j = nys-1, nyn+1 580 584 DO k = nzb+1, nzt 581 a12(k,j) = 0.5 * ( sk_p(k,j+1,i) - sk_p(k,j-1,i) )582 a22(k,j) = 0.5 * ( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i)&583 + sk_p(k,j-1,i) )584 a0(k,j) = ( 9.0 * sk_p(k,j+2,i) - 116.0 * sk_p(k,j+1,i)&585 + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k,j-1,i)&586 + 9.0 * sk_p(k,j-2,i) ) * f1920587 a1(k,j) = ( -5.0 * sk_p(k,j+2,i) + 34.0 * sk_p(k,j+1,i)&588 - 34.0 * sk_p(k,j-1,i) + 5.0 * sk_p(k,j-2,i)&585 a12(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - sk_p(k,j-1,i) ) 586 a22(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) & 587 + sk_p(k,j-1,i) ) 588 a0(k,j) = ( 9.0_wp * sk_p(k,j+2,i) - 116.0_wp * sk_p(k,j+1,i) & 589 + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j-1,i) & 590 + 9.0_wp * sk_p(k,j-2,i) ) * f1920 591 a1(k,j) = ( -5.0_wp * sk_p(k,j+2,i) + 34.0_wp * sk_p(k,j+1,i) & 592 - 34.0_wp * sk_p(k,j-1,i) + 5.0_wp * sk_p(k,j-2,i) & 589 593 ) * f48 590 a2(k,j) = ( -3.0 * sk_p(k,j+2,i) + 36.0 * sk_p(k,j+1,i)&591 - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k,j-1,i)&592 - 3.0 * sk_p(k,j-2,i) ) * f48594 a2(k,j) = ( -3.0_wp * sk_p(k,j+2,i) + 36.0_wp * sk_p(k,j+1,i) & 595 - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j-1,i) & 596 - 3.0_wp * sk_p(k,j-2,i) ) * f48 593 597 ENDDO 594 598 ENDDO … … 601 605 cip = MAX( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy ) 602 606 cim = -MIN( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy ) 603 cipf = 1.0 - 2.0* cip604 cimf = 1.0 - 2.0* cim605 ip = a0(k,j) * f2 * ( 1.0 - cipf )&606 + a1(k,j) * f8 * ( 1.0 - cipf*cipf )&607 + a2(k,j) * f24 * ( 1.0 - cipf*cipf*cipf )608 im = a0(k,j+1) * f2 * ( 1.0 - cimf )&609 - a1(k,j+1) * f8 * ( 1.0 - cimf*cimf )&610 + a2(k,j+1) * f24 * ( 1.0 - cimf*cimf*cimf )607 cipf = 1.0_wp - 2.0_wp * cip 608 cimf = 1.0_wp - 2.0_wp * cim 609 ip = a0(k,j) * f2 * ( 1.0_wp - cipf ) & 610 + a1(k,j) * f8 * ( 1.0_wp - cipf*cipf ) & 611 + a2(k,j) * f24 * ( 1.0_wp - cipf*cipf*cipf ) 612 im = a0(k,j+1) * f2 * ( 1.0_wp - cimf ) & 613 - a1(k,j+1) * f8 * ( 1.0_wp - cimf*cimf ) & 614 + a2(k,j+1) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 611 615 ip = MAX( ip, 0.0_wp ) 612 616 im = MAX( im, 0.0_wp ) 613 ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15 ) )614 impb(k,j) = im * MIN( 1.0_wp, sk_p(k,j+1,i) / (ip+im+1E-15 ) )617 ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15_wp) ) 618 impb(k,j) = im * MIN( 1.0_wp, sk_p(k,j+1,i) / (ip+im+1E-15_wp) ) 615 619 616 620 cip = MAX( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy ) 617 621 cim = -MIN( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy ) 618 cipf = 1.0 - 2.0* cip619 cimf = 1.0 - 2.0* cim620 ip = a0(k,j-1) * f2 * ( 1.0 - cipf )&621 + a1(k,j-1) * f8 * ( 1.0 - cipf*cipf )&622 + a2(k,j-1) * f24 * ( 1.0 - cipf*cipf*cipf )623 im = a0(k,j) * f2 * ( 1.0 - cimf )&624 - a1(k,j) * f8 * ( 1.0 - cimf*cimf )&625 + a2(k,j) * f24 * ( 1.0 - cimf*cimf*cimf )622 cipf = 1.0_wp - 2.0_wp * cip 623 cimf = 1.0_wp - 2.0_wp * cim 624 ip = a0(k,j-1) * f2 * ( 1.0_wp - cipf ) & 625 + a1(k,j-1) * f8 * ( 1.0_wp - cipf*cipf ) & 626 + a2(k,j-1) * f24 * ( 1.0_wp - cipf*cipf*cipf ) 627 im = a0(k,j) * f2 * ( 1.0_wp - cimf ) & 628 - a1(k,j) * f8 * ( 1.0_wp - cimf*cimf ) & 629 + a2(k,j) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 626 630 ip = MAX( ip, 0.0_wp ) 627 631 im = MAX( im, 0.0_wp ) 628 ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j-1,i) / (ip+im+1E-15 ) )629 immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15 ) )632 ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j-1,i) / (ip+im+1E-15_wp) ) 633 immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15_wp) ) 630 634 ENDDO 631 635 ENDDO … … 635 639 DO j = nys-2, nyn+2 636 640 DO k = nzb+1, nzt 637 m1z = ABS( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) + sk_p(k,j-1,i) )641 m1z = ABS( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j-1,i) ) 638 642 m1n = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) ) 639 IF ( m1n /= 0.0 .AND. m1n >= m1z ) THEN643 IF ( m1n /= 0.0_wp .AND. m1n >= m1z ) THEN 640 644 m1(k,j) = m1z / m1n 641 IF ( m1(k,j) /= 2.0 .AND. m1n < fmax(2) ) m1(k,j) = 0.0645 IF ( m1(k,j) /= 2.0_wp .AND. m1n < fmax(2) ) m1(k,j) = 0.0_wp 642 646 ELSEIF ( m1n < m1z ) THEN 643 m1(k,j) = -1.0 647 m1(k,j) = -1.0_wp 644 648 ELSE 645 m1(k,j) = 0.0 649 m1(k,j) = 0.0_wp 646 650 ENDIF 647 651 ENDDO … … 650 654 ! 651 655 !-- Compute switch sw 652 sw = 0.0 656 sw = 0.0_wp 653 657 DO j = nys-1, nyn+1 654 658 DO k = nzb+1, nzt 655 m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) / &659 m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) / & 656 660 MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp ) 657 IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0 658 659 m3 = 2.0 * ABS( a2(k,j) - a22(k,j) ) / &661 IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0_wp 662 663 m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) / & 660 664 MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp ) 661 IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0 662 663 t1 = 0.35 664 t2 = 0.35 665 IF ( m1(k,j) == -1.0 ) t2 = 0.12665 IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0_wp 666 667 t1 = 0.35_wp 668 t2 = 0.35_wp 669 IF ( m1(k,j) == -1.0_wp ) t2 = 0.12_wp 666 670 667 671 !-- *VOCL STMT,IF(10) 668 IF ( m1(k,j-1) == 1.0 .OR. m1(k,j) == 1.0 .OR. m1(k,j+1) == 1.0&669 .OR. m2 > t2 .OR. m3 > t2 .OR.&670 ( m1(k,j) > t1 .AND. m1(k,j-1) /= -1.0 .AND.&671 m1(k,j) /= -1.0 .AND. m1(k,j+1) /= -1.0 )&672 ) sw(k,j) = 1.0 672 IF ( m1(k,j-1) == 1.0_wp .OR. m1(k,j) == 1.0_wp & 673 .OR. m1(k,j+1) == 1.0_wp .OR. m2 > t2 .OR. m3 > t2 .OR. & 674 ( m1(k,j) > t1 .AND. m1(k,j-1) /= -1.0_wp .AND. & 675 m1(k,j) /= -1.0_wp .AND. m1(k,j+1) /= -1.0_wp ) & 676 ) sw(k,j) = 1.0_wp 673 677 ENDDO 674 678 ENDDO … … 681 685 682 686 !-- *VOCL STMT,IF(10) 683 IF ( sw(k,j) == 1.0 ) THEN687 IF ( sw(k,j) == 1.0_wp ) THEN 684 688 snenn = sk_p(k,j+1,i) - sk_p(k,j-1,i) 685 IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9689 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 686 690 sterm = ( sk_p(k,j,i) - sk_p(k,j-1,i) ) / snenn 687 691 sterm = MIN( sterm, 0.9999_wp ) … … 694 698 ippe(k,j) = sk_p(k,j-1,i) * cip + snenn * ( & 695 699 aex(ix) * cip + bex(ix) / dex(ix) * ( & 696 eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) & 700 eex(ix) - & 701 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) ) & 697 702 ) & 698 703 ) … … 701 706 702 707 snenn = sk_p(k,j-1,i) - sk_p(k,j+1,i) 703 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9 708 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 704 709 sterm = ( sk_p(k,j,i) - sk_p(k,j+1,i) ) / snenn 705 710 sterm = MIN( sterm, 0.9999_wp ) … … 712 717 imme(k,j) = sk_p(k,j+1,i) * cim + snenn * ( & 713 718 aex(ix) * cim + bex(ix) / dex(ix) * ( & 714 eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) & 719 eex(ix) - & 720 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) & 715 721 ) & 716 722 ) … … 720 726 721 727 !-- *VOCL STMT,IF(10) 722 IF ( sw(k,j+1) == 1.0 ) THEN728 IF ( sw(k,j+1) == 1.0_wp ) THEN 723 729 snenn = sk_p(k,j,i) - sk_p(k,j+2,i) 724 IF ( ABS( snenn ) .LT. 1E-9 ) snenn = 1E-9730 IF ( ABS( snenn ) .LT. 1E-9_wp ) snenn = 1E-9_wp 725 731 sterm = ( sk_p(k,j+1,i) - sk_p(k,j+2,i) ) / snenn 726 732 sterm = MIN( sterm, 0.9999_wp ) … … 733 739 impe(k,j) = sk_p(k,j+2,i) * cim + snenn * ( & 734 740 aex(ix) * cim + bex(ix) / dex(ix) * ( & 735 eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) & 741 eex(ix) - & 742 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) & 736 743 ) & 737 744 ) … … 741 748 742 749 !-- *VOCL STMT,IF(10) 743 IF ( sw(k,j-1) == 1.0 ) THEN750 IF ( sw(k,j-1) == 1.0_wp ) THEN 744 751 snenn = sk_p(k,j,i) - sk_p(k,j-2,i) 745 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9 752 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 746 753 sterm = ( sk_p(k,j-1,i) - sk_p(k,j-2,i) ) / snenn 747 754 sterm = MIN( sterm, 0.9999_wp ) … … 754 761 ipme(k,j) = sk_p(k,j-2,i) * cip + snenn * ( & 755 762 aex(ix) * cip + bex(ix) / dex(ix) * ( & 756 eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) & 763 eex(ix) - & 764 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) ) & 757 765 ) & 758 766 ) … … 769 777 DO j = nys, nyn 770 778 DO k = nzb+1, nzt 771 fplus = ( 1.0 - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j)&772 - ( 1.0 - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j)773 fminus = ( 1.0 - sw(k,j-1) ) * ipmb(k,j) + sw(k,j-1) * ipme(k,j)&774 - ( 1.0 - sw(k,j) ) * immb(k,j) - sw(k,j) * imme(k,j)779 fplus = ( 1.0_wp - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j) & 780 - ( 1.0_wp - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j) 781 fminus = ( 1.0_wp - sw(k,j-1) ) * ipmb(k,j) + sw(k,j-1) * ipme(k,j) & 782 - ( 1.0_wp - sw(k,j) ) * immb(k,j) - sw(k,j) * imme(k,j) 775 783 tendcy = fplus - fminus 776 784 ! … … 781 789 !-- Density correction because of possible remaining divergences 782 790 d_new = d(k,j,i) - ( v(k,j+1,i) - v(k,j,i) ) * dt_3d * ddy 783 sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /&784 ( 1.0 + d_new )791 sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) / & 792 ( 1.0_wp + d_new ) 785 793 d(k,j,i) = d_new 786 794 ENDDO … … 800 808 !-- Initialise for the computation of heat fluxes (see below; required in 801 809 !-- UP flow_statistics) 802 IF ( sk_char == 'pt' ) sums_wsts_bc_l = 0.0 810 IF ( sk_char == 'pt' ) sums_wsts_bc_l = 0.0_wp 803 811 804 812 ! … … 939 947 ! 940 948 !-- Determine the maxima of the first and second derivative in z-direction 941 fmax_l = 0.0 949 fmax_l = 0.0_wp 942 950 DO i = nxl, nxr 943 951 DO j = nys, nyn 944 952 DO k = nzb, nzt+1 945 numera = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )953 numera = ABS( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k-1,j,i) ) 946 954 denomi = ABS( sk_p(k+1,j,i+1) - sk_p(k-1,j,i) ) 947 955 fmax_l(1) = MAX( fmax_l(1) , numera ) … … 957 965 #endif 958 966 959 fmax = 0.04 * fmax967 fmax = 0.04_wp * fmax 960 968 961 969 ! … … 970 978 sw(nzb:nzt+1,nys:nyn) & 971 979 ) 972 imme = 0.0 ; impe = 0.0; ipme = 0.0; ippe = 0.0980 imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp 973 981 974 982 ! … … 980 988 DO j = nys, nyn 981 989 DO k = nzb, nzt+1 982 a12(k,j) = 0.5 * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) )983 a22(k,j) = 0.5 * ( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i)&984 + sk_p(k-1,j,i) )985 a0(k,j) = ( 9.0 * sk_p(k+2,j,i) - 116.0 * sk_p(k+1,j,i)&986 + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k-1,j,i)&987 + 9.0 * sk_p(k-2,j,i) ) * f1920988 a1(k,j) = ( -5.0 * sk_p(k+2,j,i) + 34.0 * sk_p(k+1,j,i)&989 - 34.0 * sk_p(k-1,j,i) + 5.0 * sk_p(k-2,j,i)&990 a12(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) ) 991 a22(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) & 992 + sk_p(k-1,j,i) ) 993 a0(k,j) = ( 9.0_wp * sk_p(k+2,j,i) - 116.0_wp * sk_p(k+1,j,i) & 994 + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k-1,j,i) & 995 + 9.0_wp * sk_p(k-2,j,i) ) * f1920 996 a1(k,j) = ( -5.0_wp * sk_p(k+2,j,i) + 34.0_wp * sk_p(k+1,j,i) & 997 - 34.0_wp * sk_p(k-1,j,i) + 5.0_wp * sk_p(k-2,j,i) & 990 998 ) * f48 991 a2(k,j) = ( -3.0 * sk_p(k+2,j,i) + 36.0 * sk_p(k+1,j,i)&992 - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k-1,j,i)&993 - 3.0 * sk_p(k-2,j,i) ) * f48999 a2(k,j) = ( -3.0_wp * sk_p(k+2,j,i) + 36.0_wp * sk_p(k+1,j,i) & 1000 - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k-1,j,i) & 1001 - 3.0_wp * sk_p(k-2,j,i) ) * f48 994 1002 ENDDO 995 1003 ENDDO … … 1002 1010 cip = MAX( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) ) 1003 1011 cim = -MIN( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) ) 1004 cipf = 1.0 - 2.0* cip1005 cimf = 1.0 - 2.0* cim1006 ip = a0(k,j) * f2 * ( 1.0 - cipf )&1007 + a1(k,j) * f8 * ( 1.0 - cipf*cipf )&1008 + a2(k,j) * f24 * ( 1.0 - cipf*cipf*cipf )1009 im = a0(k+1,j) * f2 * ( 1.0 - cimf )&1010 - a1(k+1,j) * f8 * ( 1.0 - cimf*cimf )&1011 + a2(k+1,j) * f24 * ( 1.0 - cimf*cimf*cimf )1012 cipf = 1.0_wp - 2.0_wp * cip 1013 cimf = 1.0_wp - 2.0_wp * cim 1014 ip = a0(k,j) * f2 * ( 1.0_wp - cipf ) & 1015 + a1(k,j) * f8 * ( 1.0_wp - cipf*cipf ) & 1016 + a2(k,j) * f24 * ( 1.0_wp - cipf*cipf*cipf ) 1017 im = a0(k+1,j) * f2 * ( 1.0_wp - cimf ) & 1018 - a1(k+1,j) * f8 * ( 1.0_wp - cimf*cimf ) & 1019 + a2(k+1,j) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 1012 1020 ip = MAX( ip, 0.0_wp ) 1013 1021 im = MAX( im, 0.0_wp ) 1014 ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15 ) )1015 impb(k,j) = im * MIN( 1.0_wp, sk_p(k+1,j,i) / (ip+im+1E-15 ) )1022 ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15_wp) ) 1023 impb(k,j) = im * MIN( 1.0_wp, sk_p(k+1,j,i) / (ip+im+1E-15_wp) ) 1016 1024 1017 1025 cip = MAX( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) ) 1018 1026 cim = -MIN( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) ) 1019 cipf = 1.0 - 2.0* cip1020 cimf = 1.0 - 2.0* cim1021 ip = a0(k-1,j) * f2 * ( 1.0 - cipf )&1022 + a1(k-1,j) * f8 * ( 1.0 - cipf*cipf )&1023 + a2(k-1,j) * f24 * ( 1.0 - cipf*cipf*cipf )1024 im = a0(k,j) * f2 * ( 1.0 - cimf )&1025 - a1(k,j) * f8 * ( 1.0 - cimf*cimf )&1026 + a2(k,j) * f24 * ( 1.0 - cimf*cimf*cimf )1027 cipf = 1.0_wp - 2.0_wp * cip 1028 cimf = 1.0_wp - 2.0_wp * cim 1029 ip = a0(k-1,j) * f2 * ( 1.0_wp - cipf ) & 1030 + a1(k-1,j) * f8 * ( 1.0_wp - cipf*cipf ) & 1031 + a2(k-1,j) * f24 * ( 1.0_wp - cipf*cipf*cipf ) 1032 im = a0(k,j) * f2 * ( 1.0_wp - cimf ) & 1033 - a1(k,j) * f8 * ( 1.0_wp - cimf*cimf ) & 1034 + a2(k,j) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 1027 1035 ip = MAX( ip, 0.0_wp ) 1028 1036 im = MAX( im, 0.0_wp ) 1029 ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k-1,j,i) / (ip+im+1E-15 ) )1030 immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15 ) )1037 ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k-1,j,i) / (ip+im+1E-15_wp) ) 1038 immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i) / (ip+im+1E-15_wp) ) 1031 1039 ENDDO 1032 1040 ENDDO … … 1036 1044 DO j = nys, nyn 1037 1045 DO k = nzb-1, nzt+2 1038 m1z = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )1046 m1z = ABS( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k-1,j,i) ) 1039 1047 m1n = ABS( sk_p(k+1,j,i) - sk_p(k-1,j,i) ) 1040 IF ( m1n /= 0.0 .AND. m1n >= m1z ) THEN1048 IF ( m1n /= 0.0_wp .AND. m1n >= m1z ) THEN 1041 1049 m1(k,j) = m1z / m1n 1042 IF ( m1(k,j) /= 2.0 .AND. m1n < fmax(2) ) m1(k,j) = 0.01050 IF ( m1(k,j) /= 2.0_wp .AND. m1n < fmax(2) ) m1(k,j) = 0.0_wp 1043 1051 ELSEIF ( m1n < m1z ) THEN 1044 m1(k,j) = -1.0 1052 m1(k,j) = -1.0_wp 1045 1053 ELSE 1046 m1(k,j) = 0.0 1054 m1(k,j) = 0.0_wp 1047 1055 ENDIF 1048 1056 ENDDO … … 1051 1059 ! 1052 1060 !-- Compute switch sw 1053 sw = 0.0 1061 sw = 0.0_wp 1054 1062 DO j = nys, nyn 1055 1063 DO k = nzb, nzt+1 1056 m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) /&1064 m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) / & 1057 1065 MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp ) 1058 IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0 1059 1060 m3 = 2.0 * ABS( a2(k,j) - a22(k,j) ) /&1066 IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0_wp 1067 1068 m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) / & 1061 1069 MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp ) 1062 IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0 1063 1064 t1 = 0.35 1065 t2 = 0.35 1066 IF ( m1(k,j) == -1.0 ) t2 = 0.121070 IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0_wp 1071 1072 t1 = 0.35_wp 1073 t2 = 0.35_wp 1074 IF ( m1(k,j) == -1.0_wp ) t2 = 0.12_wp 1067 1075 1068 1076 !-- *VOCL STMT,IF(10) 1069 IF ( m1(k-1,j) == 1.0 .OR. m1(k,j) == 1.0 .OR. m1(k+1,j) == 1.0&1070 .OR. m2 > t2 .OR. m3 > t2 .OR.&1071 ( m1(k,j) > t1 .AND. m1(k-1,j) /= -1.0 .AND.&1072 m1(k,j) /= -1.0 .AND. m1(k+1,j) /= -1.0 )&1073 ) sw(k,j) = 1.0 1077 IF ( m1(k-1,j) == 1.0_wp .OR. m1(k,j) == 1.0_wp & 1078 .OR. m1(k+1,j) == 1.0_wp .OR. m2 > t2 .OR. m3 > t2 .OR. & 1079 ( m1(k,j) > t1 .AND. m1(k-1,j) /= -1.0_wp .AND. & 1080 m1(k,j) /= -1.0_wp .AND. m1(k+1,j) /= -1.0_wp ) & 1081 ) sw(k,j) = 1.0_wp 1074 1082 ENDDO 1075 1083 ENDDO … … 1082 1090 1083 1091 !-- *VOCL STMT,IF(10) 1084 IF ( sw(k,j) == 1.0 ) THEN1092 IF ( sw(k,j) == 1.0_wp ) THEN 1085 1093 snenn = sk_p(k+1,j,i) - sk_p(k-1,j,i) 1086 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9 1094 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 1087 1095 sterm = ( sk_p(k,j,i) - sk_p(k-1,j,i) ) / snenn 1088 1096 sterm = MIN( sterm, 0.9999_wp ) … … 1095 1103 ippe(k,j) = sk_p(k-1,j,i) * cip + snenn * ( & 1096 1104 aex(ix) * cip + bex(ix) / dex(ix) * ( & 1097 eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) & 1105 eex(ix) - & 1106 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) ) & 1098 1107 ) & 1099 1108 ) 1100 IF ( sterm == 0.0001 ) ippe(k,j) = sk_p(k,j,i) * cip1101 IF ( sterm == 0.9999 ) ippe(k,j) = sk_p(k,j,i) * cip1109 IF ( sterm == 0.0001_wp ) ippe(k,j) = sk_p(k,j,i) * cip 1110 IF ( sterm == 0.9999_wp ) ippe(k,j) = sk_p(k,j,i) * cip 1102 1111 1103 1112 snenn = sk_p(k-1,j,i) - sk_p(k+1,j,i) 1104 IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-91113 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 1105 1114 sterm = ( sk_p(k,j,i) - sk_p(k+1,j,i) ) / snenn 1106 1115 sterm = MIN( sterm, 0.9999_wp ) … … 1113 1122 imme(k,j) = sk_p(k+1,j,i) * cim + snenn * ( & 1114 1123 aex(ix) * cim + bex(ix) / dex(ix) * ( & 1115 eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) & 1124 eex(ix) - & 1125 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) & 1116 1126 ) & 1117 1127 ) … … 1121 1131 1122 1132 !-- *VOCL STMT,IF(10) 1123 IF ( sw(k+1,j) == 1.0 ) THEN1133 IF ( sw(k+1,j) == 1.0_wp ) THEN 1124 1134 snenn = sk_p(k,j,i) - sk_p(k+2,j,i) 1125 IF ( ABS( snenn ) .LT. 1E-9 ) snenn = 1E-91135 IF ( ABS( snenn ) .LT. 1E-9_wp ) snenn = 1E-9_wp 1126 1136 sterm = ( sk_p(k+1,j,i) - sk_p(k+2,j,i) ) / snenn 1127 1137 sterm = MIN( sterm, 0.9999_wp ) … … 1134 1144 impe(k,j) = sk_p(k+2,j,i) * cim + snenn * ( & 1135 1145 aex(ix) * cim + bex(ix) / dex(ix) * ( & 1136 eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) & 1146 eex(ix) - & 1147 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) & 1137 1148 ) & 1138 1149 ) … … 1142 1153 1143 1154 !-- *VOCL STMT,IF(10) 1144 IF ( sw(k-1,j) == 1.0 ) THEN1155 IF ( sw(k-1,j) == 1.0_wp ) THEN 1145 1156 snenn = sk_p(k,j,i) - sk_p(k-2,j,i) 1146 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9 1157 IF ( ABS( snenn ) < 1E-9_wp ) snenn = 1E-9_wp 1147 1158 sterm = ( sk_p(k-1,j,i) - sk_p(k-2,j,i) ) / snenn 1148 1159 sterm = MIN( sterm, 0.9999_wp ) … … 1155 1166 ipme(k,j) = sk_p(k-2,j,i) * cip + snenn * ( & 1156 1167 aex(ix) * cip + bex(ix) / dex(ix) * ( & 1157 eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) & 1168 eex(ix) - & 1169 EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) ) & 1158 1170 ) & 1159 1171 ) … … 1170 1182 DO j = nys, nyn 1171 1183 DO k = nzb+1, nzt 1172 fplus = ( 1.0 - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j)&1173 - ( 1.0 - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j)1174 fminus = ( 1.0 - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j)&1175 - ( 1.0 - sw(k,j) ) * immb(k,j) - sw(k,j) * imme(k,j)1184 fplus = ( 1.0_wp - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j) & 1185 - ( 1.0_wp - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j) 1186 fminus = ( 1.0_wp - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j) & 1187 - ( 1.0_wp - sw(k,j) ) * immb(k,j) - sw(k,j) * imme(k,j) 1176 1188 tendcy = fplus - fminus 1177 1189 ! … … 1182 1194 !-- Density correction because of possible remaining divergences 1183 1195 d_new = d(k,j,i) - ( w(k,j,i) - w(k-1,j,i) ) * dt_3d * ddzw(k) 1184 sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /&1185 ( 1.0 + d_new )1196 sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) / & 1197 ( 1.0_wp + d_new ) 1186 1198 ! 1187 1199 !-- Store heat flux for subsequent statistics output.
Note: See TracChangeset
for help on using the changeset viewer.
