Changeset 1342
- Timestamp:
- Mar 26, 2014 5:04:47 PM (10 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 11 edited
-
calc_spectra.f90 (modified) (11 diffs)
-
fft_xy.f90 (modified) (33 diffs)
-
pres.f90 (modified) (12 diffs)
-
production_e.f90 (modified) (48 diffs)
-
random_function.f90 (modified) (2 diffs)
-
random_gauss.f90 (modified) (2 diffs)
-
temperton_fft.f90 (modified) (42 diffs)
-
time_integration.f90 (modified) (8 diffs)
-
time_to_string.f90 (modified) (2 diffs)
-
timestep.f90 (modified) (16 diffs)
-
tridia_solver.f90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/calc_spectra.f90
r1325 r1342 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! REAL constants defined as wp-kinds 22 23 ! 23 24 ! Former revisions: … … 315 316 ! 316 317 !-- Exponent for geometric average 317 exponent = 1.0 / ( ny + 1.0)318 exponent = 1.0_wp / ( ny + 1.0_wp ) 318 319 319 320 ! … … 346 347 DO i = 0, nx/2 347 348 348 sums_spectra_l(i) = 1.0 349 sums_spectra_l(i) = 1.0_wp 349 350 DO j = nys_x, nyn_x 350 351 sums_spectra_l(i) = sums_spectra_l(i) * ddd(i,j,k)**exponent … … 355 356 ELSE 356 357 357 sums_spectra_l = 1.0 358 sums_spectra_l = 1.0_wp 358 359 359 360 ENDIF … … 361 362 ! 362 363 !-- Global sum of spectra on PE0 (from where they are written on file) 363 sums_spectra(:,n) = 0.0 364 sums_spectra(:,n) = 0.0_wp 364 365 #if defined( __parallel ) 365 366 CALL MPI_BARRIER( comm2d, ierr ) ! Necessary? … … 380 381 !-- Temperton fft results need to be normalized 381 382 IF ( fft_method == 'temperton-algorithm' ) THEN 382 fac = nx + 1.0 383 fac = nx + 1.0_wp 383 384 ELSE 384 fac = 1.0 385 fac = 1.0_wp 385 386 ENDIF 386 387 DO i = 1, nx/2 … … 454 455 ! 455 456 !-- Exponent for geometric average 456 exponent = 1.0 / ( nx + 1.0)457 exponent = 1.0_wp / ( nx + 1.0_wp ) 457 458 458 459 ! … … 485 486 DO j = 0, ny/2 486 487 487 sums_spectra_l(j) = 1.0 488 sums_spectra_l(j) = 1.0_wp 488 489 DO i = nxl_yd, nxr_yd 489 490 sums_spectra_l(j) = sums_spectra_l(j) * ddd(j,i,k)**exponent … … 494 495 ELSE 495 496 496 sums_spectra_l = 1.0 497 sums_spectra_l = 1.0_wp 497 498 498 499 ENDIF … … 500 501 ! 501 502 !-- Global sum of spectra on PE0 (from where they are written on file) 502 sums_spectra(:,n) = 0.0 503 sums_spectra(:,n) = 0.0_wp 503 504 #if defined( __parallel ) 504 505 CALL MPI_BARRIER( comm2d, ierr ) ! Necessary? … … 520 521 !-- Temperton fft results need to be normalized 521 522 IF ( fft_method == 'temperton-algorithm' ) THEN 522 fac = ny + 1.0 523 fac = ny + 1.0_wp 523 524 ELSE 524 fac = 1.0 525 fac = 1.0_wp 525 526 ENDIF 526 527 DO j = 1, ny/2 -
palm/trunk/SOURCE/fft_xy.f90
r1323 r1342 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 242 242 IF ( fft_method == 'system-specific' ) THEN 243 243 244 dnx = 1.0 / ( nx + 1.0)245 dny = 1.0 / ( ny + 1.0)244 dnx = 1.0_wp / ( nx + 1.0_wp ) 245 dny = 1.0_wp / ( ny + 1.0_wp ) 246 246 sqr_dnx = SQRT( dnx ) 247 247 sqr_dny = SQRT( dny ) … … 267 267 trig_yb(2*(ny+1)), trig_yf(2*(ny+1)) ) 268 268 269 work_x = 0.0 270 work_y = 0.0 269 work_x = 0.0_wp 270 work_y = 0.0_wp 271 271 nz1 = nz + MOD( nz+1, 2 ) ! odd nz slows down fft significantly 272 272 ! when using the NEC ffts … … 433 433 DO j = nys_x, nyn_x 434 434 435 cwork(0) = CMPLX( ar(0,j,k), 0.0 )435 cwork(0) = CMPLX( ar(0,j,k), 0.0_wp ) 436 436 DO i = 1, (nx+1)/2 - 1 437 437 cwork(i) = CMPLX( ar(i,j,k), -ar(nx+1-i,j,k) ) 438 438 cwork(nx+1-i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) ) 439 439 ENDDO 440 cwork((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )440 cwork((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp ) 441 441 442 442 ishape = SHAPE( cwork ) … … 494 494 work(2*i+1) = ar(nx+1-i,j,k) 495 495 ENDDO 496 work(1) = 0.0 497 work(nx+2) = 0.0 496 work(1) = 0.0_wp 497 work(nx+2) = 0.0_wp 498 498 499 499 CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 ) … … 551 551 IF ( PRESENT( ar_2d ) ) THEN 552 552 553 x_out(0) = CMPLX( ar_2d(0,j), 0.0 )553 x_out(0) = CMPLX( ar_2d(0,j), 0.0_wp ) 554 554 DO i = 1, (nx+1)/2 - 1 555 555 x_out(i) = CMPLX( ar_2d(i,j), ar_2d(nx+1-i,j) ) 556 556 ENDDO 557 x_out((nx+1)/2) = CMPLX( ar_2d((nx+1)/2,j), 0.0 )557 x_out((nx+1)/2) = CMPLX( ar_2d((nx+1)/2,j), 0.0_wp ) 558 558 559 559 ELSE 560 560 561 x_out(0) = CMPLX( ar(0,j,k), 0.0 )561 x_out(0) = CMPLX( ar(0,j,k), 0.0_wp ) 562 562 DO i = 1, (nx+1)/2 - 1 563 563 x_out(i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) ) 564 564 ENDDO 565 x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )565 x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp ) 566 566 567 567 ENDIF … … 614 614 work(2*i+1) = ar(nx+1-i,j,k) 615 615 ENDDO 616 work(1) = 0.0 617 work(nx+2) = 0.0 616 work(1) = 0.0_wp 617 work(nx+2) = 0.0_wp 618 618 619 619 CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, & … … 667 667 work(2*i+1) = ar(nx+1-i,j,k) 668 668 ENDDO 669 work(1) = 0.0 670 work(nx+2) = 0.0 669 work(1) = 0.0_wp 670 work(nx+2) = 0.0_wp 671 671 672 672 CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 ) … … 711 711 DO j = nys_x, nyn_x 712 712 713 ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0 )713 ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0_wp ) 714 714 715 715 DO i = 1, (nx+1)/2 - 1 716 716 ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) ) 717 717 ENDDO 718 ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0 )718 ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp ) 719 719 720 720 ENDDO … … 805 805 ELSE 806 806 807 cwork(0) = CMPLX( ar(0), 0.0 )807 cwork(0) = CMPLX( ar(0), 0.0_wp ) 808 808 DO i = 1, (nx+1)/2 - 1 809 809 cwork(i) = CMPLX( ar(i), -ar(nx+1-i) ) 810 810 cwork(nx+1-i) = CMPLX( ar(i), ar(nx+1-i) ) 811 811 ENDDO 812 cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )812 cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0_wp ) 813 813 814 814 ishape = SHAPE( cwork ) … … 848 848 work(2*i+1) = ar(nx+1-i) 849 849 ENDDO 850 work(1) = 0.0 851 work(nx+2) = 0.0 850 work(1) = 0.0_wp 851 work(nx+2) = 0.0_wp 852 852 853 853 CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 ) … … 873 873 ELSE 874 874 875 x_out(0) = CMPLX( ar(0), 0.0 )875 x_out(0) = CMPLX( ar(0), 0.0_wp ) 876 876 DO i = 1, (nx+1)/2 - 1 877 877 x_out(i) = CMPLX( ar(i), ar(nx+1-i) ) 878 878 ENDDO 879 x_out((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )879 x_out((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0_wp ) 880 880 881 881 CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in) … … 908 908 work(2*i+1) = ar(nx+1-i) 909 909 ENDDO 910 work(1) = 0.0 911 work(nx+2) = 0.0 910 work(1) = 0.0_wp 911 work(nx+2) = 0.0_wp 912 912 913 913 CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, & … … 941 941 work(2*i+1) = ar(nx+1-i) 942 942 ENDDO 943 work(1) = 0.0 944 work(nx+2) = 0.0 943 work(1) = 0.0_wp 944 work(nx+2) = 0.0_wp 945 945 946 946 CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 ) … … 1076 1076 DO i = nxl_y_l, nxr_y_l 1077 1077 1078 cwork(0) = CMPLX( ar_tr(0,i,k), 0.0 )1078 cwork(0) = CMPLX( ar_tr(0,i,k), 0.0_wp ) 1079 1079 DO j = 1, (ny+1)/2 - 1 1080 1080 cwork(j) = CMPLX( ar_tr(j,i,k), -ar_tr(ny+1-j,i,k) ) 1081 1081 cwork(ny+1-j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k) ) 1082 1082 ENDDO 1083 cwork((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0 )1083 cwork((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0_wp ) 1084 1084 1085 1085 jshape = SHAPE( cwork ) … … 1137 1137 work(2*j+1) = ar_tr(ny+1-j,i,k) 1138 1138 ENDDO 1139 work(1) = 0.0 1140 work(ny+2) = 0.0 1139 work(1) = 0.0_wp 1140 work(ny+2) = 0.0_wp 1141 1141 1142 1142 CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 ) … … 1180 1180 DO i = nxl_y_l, nxr_y_l 1181 1181 1182 y_out(0) = CMPLX( ar_tr(0,i,k), 0.0 )1182 y_out(0) = CMPLX( ar_tr(0,i,k), 0.0_wp ) 1183 1183 DO j = 1, (ny+1)/2 - 1 1184 1184 y_out(j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k) ) 1185 1185 ENDDO 1186 y_out((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0 )1186 y_out((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0_wp ) 1187 1187 1188 1188 CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in ) … … 1233 1233 work(2*j+1) = ar_tr(ny+1-j,i,k) 1234 1234 ENDDO 1235 work(1) = 0.0 1236 work(ny+2) = 0.0 1235 work(1) = 0.0_wp 1236 work(ny+2) = 0.0_wp 1237 1237 1238 1238 CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, & … … 1284 1284 work(2*j+1) = ar_tr(ny+1-j,i,k) 1285 1285 ENDDO 1286 work(1) = 0.0 1287 work(ny+2) = 0.0 1286 work(1) = 0.0_wp 1287 work(ny+2) = 0.0_wp 1288 1288 1289 1289 CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 ) … … 1327 1327 DO i = nxl_y, nxr_y 1328 1328 1329 ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0 )1329 ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0_wp ) 1330 1330 1331 1331 DO j = 1, (ny+1)/2 - 1 1332 1332 ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) ) 1333 1333 ENDDO 1334 ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0 )1334 ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0_wp ) 1335 1335 1336 1336 ENDDO … … 1423 1423 ELSE 1424 1424 1425 cwork(0) = CMPLX( ar(0), 0.0 )1425 cwork(0) = CMPLX( ar(0), 0.0_wp ) 1426 1426 DO j = 1, (ny+1)/2 - 1 1427 1427 cwork(j) = CMPLX( ar(j), -ar(ny+1-j) ) 1428 1428 cwork(ny+1-j) = CMPLX( ar(j), ar(ny+1-j) ) 1429 1429 ENDDO 1430 cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 )1430 cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0_wp ) 1431 1431 1432 1432 jshape = SHAPE( cwork ) … … 1466 1466 work(2*j+1) = ar(ny+1-j) 1467 1467 ENDDO 1468 work(1) = 0.0 1469 work(ny+2) = 0.0 1468 work(1) = 0.0_wp 1469 work(ny+2) = 0.0_wp 1470 1470 1471 1471 CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 ) … … 1491 1491 ELSE 1492 1492 1493 y_out(0) = CMPLX( ar(0), 0.0 )1493 y_out(0) = CMPLX( ar(0), 0.0_wp ) 1494 1494 DO j = 1, (ny+1)/2 - 1 1495 1495 y_out(j) = CMPLX( ar(j), ar(ny+1-j) ) 1496 1496 ENDDO 1497 y_out((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 )1497 y_out((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0_wp ) 1498 1498 1499 1499 CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in ) … … 1526 1526 work(2*j+1) = ar(ny+1-j) 1527 1527 ENDDO 1528 work(1) = 0.0 1529 work(ny+2) = 0.0 1528 work(1) = 0.0_wp 1529 work(ny+2) = 0.0_wp 1530 1530 1531 1531 CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, & … … 1559 1559 work(2*j+1) = ar(ny+1-j) 1560 1560 ENDDO 1561 work(1) = 0.0 1562 work(ny+2) = 0.0 1561 work(1) = 0.0_wp 1562 work(ny+2) = 0.0_wp 1563 1563 1564 1564 CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 ) … … 1621 1621 1622 1622 ai(0:nx,1:nz) = ar(0:nx,1:nz) 1623 ai(nx+1:,:) = 0.0 1623 ai(nx+1:,:) = 0.0_wp 1624 1624 1625 1625 CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 ) … … 1643 1643 ai(2*i+1,k) = ar(nx+1-i,k) 1644 1644 ENDDO 1645 ai(1,k) = 0.0 1646 ai(nx+2,k) = 0.0 1645 ai(1,k) = 0.0_wp 1646 ai(nx+2,k) = 0.0_wp 1647 1647 ENDDO 1648 1648 … … 1669 1669 ai(0:nx,1:nz) = ar(0:nx,1:nz) 1670 1670 IF ( nz1 > nz ) THEN 1671 ai(:,nz1) = 0.0 1671 ai(:,nz1) = 0.0_wp 1672 1672 ENDIF 1673 1673 … … 1693 1693 1694 1694 IF ( nz1 > nz ) THEN 1695 work(:,nz1) = 0.0 1695 work(:,nz1) = 0.0_wp 1696 1696 ENDIF 1697 1697 DO k = 1, nz 1698 work(1,k) = CMPLX( ar(0,k), 0.0 )1698 work(1,k) = CMPLX( ar(0,k), 0.0_wp ) 1699 1699 DO i = 1, (nx+1)/2 - 1 1700 1700 work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k) ) 1701 1701 ENDDO 1702 work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0 )1702 work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0_wp ) 1703 1703 ENDDO 1704 1704 … … 1764 1764 1765 1765 ai(0:ny,1:nz) = ar(0:ny,1:nz) 1766 ai(ny+1:,:) = 0.0 1766 ai(ny+1:,:) = 0.0_wp 1767 1767 1768 1768 CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 ) … … 1786 1786 ai(2*j+1,k) = ar(ny+1-j,k) 1787 1787 ENDDO 1788 ai(1,k) = 0.0 1789 ai(ny+2,k) = 0.0 1788 ai(1,k) = 0.0_wp 1789 ai(ny+2,k) = 0.0_wp 1790 1790 ENDDO 1791 1791 … … 1812 1812 ai(0:ny,1:nz) = ar(0:ny,1:nz) 1813 1813 IF ( nz1 > nz ) THEN 1814 ai(:,nz1) = 0.0 1814 ai(:,nz1) = 0.0_wp 1815 1815 ENDIF 1816 1816 … … 1836 1836 1837 1837 IF ( nz1 > nz ) THEN 1838 work(:,nz1) = 0.0 1838 work(:,nz1) = 0.0_wp 1839 1839 ENDIF 1840 1840 DO k = 1, nz 1841 work(1,k) = CMPLX( ar(0,k), 0.0 )1841 work(1,k) = CMPLX( ar(0,k), 0.0_wp ) 1842 1842 DO j = 1, (ny+1)/2 - 1 1843 1843 work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k) ) 1844 1844 ENDDO 1845 work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0 )1845 work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0_wp ) 1846 1846 ENDDO 1847 1847 -
palm/trunk/SOURCE/pres.f90
r1321 r1342 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 135 135 136 136 137 ddt_3d = 1.0 / dt_3d138 d_weight_pres = 1.0 / weight_pres(intermediate_timestep_count)137 ddt_3d = 1.0_wp / dt_3d 138 d_weight_pres = 1.0_wp / weight_pres(intermediate_timestep_count) 139 139 140 140 ! … … 180 180 IF ( conserve_volume_flow .AND. ( outflow_l .OR. outflow_r ) ) THEN 181 181 182 volume_flow(1) = 0.0 183 volume_flow_l(1) = 0.0 182 volume_flow(1) = 0.0_wp 183 volume_flow_l(1) = 0.0_wp 184 184 185 185 IF ( outflow_l ) THEN … … 219 219 IF ( conserve_volume_flow .AND. ( outflow_n .OR. outflow_s ) ) THEN 220 220 221 volume_flow(2) = 0.0 222 volume_flow_l(2) = 0.0 221 volume_flow(2) = 0.0_wp 222 volume_flow_l(2) = 0.0_wp 223 223 224 224 IF ( outflow_s ) THEN … … 257 257 !-- Remove mean vertical velocity 258 258 IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN 259 IF ( simulated_time > 0.0 ) THEN ! otherwise nzb_w_inner not yet known260 w_l = 0.0 ; w_l_l = 0.0259 IF ( simulated_time > 0.0_wp ) THEN ! otherwise nzb_w_inner not yet known 260 w_l = 0.0_wp; w_l_l = 0.0_wp 261 261 DO i = nxl, nxr 262 262 DO j = nys, nyn … … 295 295 DO j = nys-1, nyn+1 296 296 DO k = nzb, nzt+1 297 d(k,j,i) = 0.0 297 d(k,j,i) = 0.0_wp 298 298 ENDDO 299 299 ENDDO … … 305 305 DO j = nys, nyn 306 306 DO k = nzb+1, nzt 307 d(k,j,i) = 0.0 307 d(k,j,i) = 0.0_wp 308 308 ENDDO 309 309 ENDDO … … 312 312 ENDIF 313 313 314 localsum = 0.0 315 threadsum = 0.0 314 localsum = 0.0_wp 315 threadsum = 0.0_wp 316 316 317 317 #if defined( __ibm ) … … 432 432 DO i = nxlg, nxrg 433 433 DO j = nysg, nyng 434 tend(nzb_s_inner(j,i),j,i) = 0.0 434 tend(nzb_s_inner(j,i),j,i) = 0.0_wp 435 435 ENDDO 436 436 ENDDO … … 460 460 DO i = nxlg, nxrg 461 461 DO j = nysg, nyng 462 tend(nzt+1,j,i) = 0.0 462 tend(nzt+1,j,i) = 0.0_wp 463 463 ENDDO 464 464 ENDDO … … 577 577 !-- pressure just computed 578 578 IF ( conserve_volume_flow .AND. ( bc_lr_cyc .OR. bc_ns_cyc ) ) THEN 579 volume_flow_l(1) = 0.0 580 volume_flow_l(2) = 0.0 579 volume_flow_l(1) = 0.0_wp 580 volume_flow_l(2) = 0.0_wp 581 581 ENDIF 582 582 … … 704 704 !-- a possible PE-sum is computed in flow_statistics 705 705 CALL cpu_log( log_point_s(1), 'divergence', 'start' ) 706 sums_divnew_l = 0.0 706 sums_divnew_l = 0.0_wp 707 707 708 708 ! 709 709 !-- d must be reset to zero because it can contain nonzero values below the 710 710 !-- topography 711 IF ( topography /= 'flat' ) d = 0.0 712 713 localsum = 0.0 714 threadsum = 0.0 711 IF ( topography /= 'flat' ) d = 0.0_wp 712 713 localsum = 0.0_wp 714 threadsum = 0.0_wp 715 715 716 716 !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum) -
palm/trunk/SOURCE/production_e.f90
r1321 r1342 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 173 173 DO k = nzb_diff_s_outer(j,i), nzt 174 174 175 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx176 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &177 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy178 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &179 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)180 181 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &182 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx183 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy184 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &185 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)186 187 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &188 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx189 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &190 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy191 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)192 193 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &175 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 176 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 177 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 178 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 179 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 180 181 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 182 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 183 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 184 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 185 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 186 187 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 188 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 189 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 190 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 191 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 192 193 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 194 194 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 195 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )196 197 IF ( def < 0.0 ) def = 0.0195 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 196 197 IF ( def < 0.0_wp ) def = 0.0_wp 198 198 199 199 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 210 210 DO j = nys, nyn 211 211 212 IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0) ) &212 IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) & 213 213 THEN 214 214 215 215 k = nzb_diff_s_inner(j,i) - 1 216 216 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 217 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &218 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)217 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 218 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) 219 219 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 220 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &221 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)220 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 221 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) 222 222 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 223 223 224 IF ( wall_e_y(j,i) /= 0.0 ) THEN224 IF ( wall_e_y(j,i) /= 0.0_wp ) THEN 225 225 ! 226 226 !-- Inconsistency removed: as the thermal stratification is … … 236 236 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 237 237 wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp ) 238 km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25 * &239 0.5 * dy240 IF ( km_neutral > 0.0 ) THEN238 km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25_wp * & 239 0.5_wp * dy 240 IF ( km_neutral > 0.0_wp ) THEN 241 241 dudy = - wall_e_y(j,i) * usvs(k) / km_neutral 242 242 dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral 243 243 ELSE 244 dudy = 0.0 245 dwdy = 0.0 244 dudy = 0.0_wp 245 dwdy = 0.0_wp 246 246 ENDIF 247 247 ELSE 248 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &249 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy250 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &251 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy248 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 249 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 250 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 251 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 252 252 ENDIF 253 253 254 IF ( wall_e_x(j,i) /= 0.0 ) THEN254 IF ( wall_e_x(j,i) /= 0.0_wp ) THEN 255 255 ! 256 256 !-- Inconsistency removed: as the thermal stratification is … … 266 266 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 267 267 wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp ) 268 km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25 * &269 0.5 * dx270 IF ( km_neutral > 0.0 ) THEN268 km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25_wp * & 269 0.5_wp * dx 270 IF ( km_neutral > 0.0_wp ) THEN 271 271 dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral 272 272 dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral 273 273 ELSE 274 dvdx = 0.0 275 dwdx = 0.0 274 dvdx = 0.0_wp 275 dwdx = 0.0_wp 276 276 ENDIF 277 277 ELSE 278 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &279 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx280 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &281 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx278 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 279 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 280 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 281 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 282 282 ENDIF 283 283 284 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &284 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 285 285 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 286 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )287 288 IF ( def < 0.0 ) def = 0.0286 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 287 288 IF ( def < 0.0_wp ) def = 0.0_wp 289 289 290 290 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 301 301 302 302 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 303 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &304 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)305 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy306 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &307 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)303 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 304 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 305 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 306 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 307 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 308 308 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 309 309 310 IF ( wall_e_y(j,i) /= 0.0 ) THEN310 IF ( wall_e_y(j,i) /= 0.0_wp ) THEN 311 311 ! 312 312 !-- Inconsistency removed: as the thermal stratification … … 319 319 !-- validation has been available 320 320 km_neutral = kappa * ( usvs(k)**2 + & 321 wsvs(k)**2 )**0.25 * 0.5* dy322 IF ( km_neutral > 0.0 ) THEN321 wsvs(k)**2 )**0.25_wp * 0.5_wp * dy 322 IF ( km_neutral > 0.0_wp ) THEN 323 323 dudy = - wall_e_y(j,i) * usvs(k) / km_neutral 324 324 dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral 325 325 ELSE 326 dudy = 0.0 327 dwdy = 0.0 326 dudy = 0.0_wp 327 dwdy = 0.0_wp 328 328 ENDIF 329 329 ELSE 330 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &331 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy332 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &333 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy330 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 331 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 332 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 333 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 334 334 ENDIF 335 335 336 IF ( wall_e_x(j,i) /= 0.0 ) THEN336 IF ( wall_e_x(j,i) /= 0.0_wp ) THEN 337 337 ! 338 338 !-- Inconsistency removed: as the thermal stratification … … 345 345 !-- validation has been available 346 346 km_neutral = kappa * ( vsus(k)**2 + & 347 wsus(k)**2 )**0.25 * 0.5* dx348 IF ( km_neutral > 0.0 ) THEN347 wsus(k)**2 )**0.25_wp * 0.5_wp * dx 348 IF ( km_neutral > 0.0_wp ) THEN 349 349 dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral 350 350 dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral 351 351 ELSE 352 dvdx = 0.0 353 dwdx = 0.0 352 dvdx = 0.0_wp 353 dwdx = 0.0_wp 354 354 ENDIF 355 355 ELSE 356 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &357 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx358 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &359 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx356 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 357 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 358 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 359 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 360 360 ENDIF 361 361 362 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &362 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 363 363 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 364 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )365 366 IF ( def < 0.0 ) def = 0.0364 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 365 366 IF ( def < 0.0_wp ) def = 0.0_wp 367 367 368 368 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 379 379 DO j = nys, nyn 380 380 381 IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0) ) &381 IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) & 382 382 THEN 383 383 384 384 k = nzb_diff_s_outer(j,i)-1 385 385 386 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx387 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &388 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy389 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &390 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)391 392 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &393 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx394 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy395 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &396 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)397 398 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &399 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx400 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &401 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy402 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)403 404 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &386 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 387 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 388 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 389 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 390 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 391 392 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 393 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 394 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 395 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 396 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 397 398 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 399 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 400 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 401 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 402 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 403 404 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 405 405 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 406 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )407 408 IF ( def < 0.0 ) def = 0.0406 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 407 408 IF ( def < 0.0_wp ) def = 0.0_wp 409 409 410 410 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 420 420 DO j = nys, nyn 421 421 422 IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0) ) &422 IF ( ( wall_e_x(j,i) == 0.0_wp ) .AND. ( wall_e_y(j,i) == 0.0_wp ) ) & 423 423 THEN 424 424 425 425 k = nzb_diff_s_inner(j,i)-1 426 426 427 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 428 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & 427 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 428 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 429 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 430 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 431 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) 432 433 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 434 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 435 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 436 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 437 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) 438 439 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 440 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 441 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 442 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 443 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 444 445 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 446 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 447 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 448 449 IF ( def < 0.0_wp ) def = 0.0_wp 450 451 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 452 453 ENDIF 454 455 ENDDO 456 457 ELSEIF ( use_surface_fluxes ) THEN 458 459 DO j = nys, nyn 460 461 k = nzb_diff_s_outer(j,i)-1 462 463 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 464 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 429 465 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 430 dudz = 0.5* ( u(k+1,j,i) + u(k+1,j,i+1) - &431 u_0(j,i) - u_0(j,i+1)) * dd2zu(k)432 433 dvdx = 0.25* ( v(k,j,i+1) + v(k,j+1,i+1) - &466 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 467 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 468 469 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 434 470 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 435 dvdy =( v(k,j+1,i) - v(k,j,i) ) * ddy436 dvdz = 0.5* ( v(k+1,j,i) + v(k+1,j+1,i) - &437 v _0(j,i) - v_0(j+1,i)) * dd2zu(k)438 439 dwdx = 0.25* ( w(k,j,i+1) + w(k-1,j,i+1) - &471 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 472 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 473 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 474 475 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 440 476 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 441 dwdy = 0.25* ( w(k,j+1,i) + w(k-1,j+1,i) - &477 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 442 478 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 443 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 444 445 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 446 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 447 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 448 449 IF ( def < 0.0 ) def = 0.0 450 451 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 452 453 ENDIF 454 455 ENDDO 456 457 ELSEIF ( use_surface_fluxes ) THEN 458 459 DO j = nys, nyn 460 461 k = nzb_diff_s_outer(j,i)-1 462 463 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 464 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & 465 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 466 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - & 467 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 468 469 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & 470 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 471 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 472 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - & 473 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 474 475 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - & 476 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 477 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - & 478 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 479 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 480 481 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 479 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 480 481 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 482 482 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 483 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )484 485 IF ( def < 0.0 ) def = 0.0483 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 484 485 IF ( def < 0.0_wp ) def = 0.0_wp 486 486 487 487 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 586 586 587 587 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 588 k1 = 1.0 + 0.61* q(k,j,i)589 k2 = 0.61 * pt(k,j,i)588 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 589 k2 = 0.61_wp * pt(k,j,i) 590 590 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * & 591 591 g / vpt(k,j,i) * & … … 594 594 ) * dd2zu(k) 595 595 ELSE IF ( cloud_physics ) THEN 596 IF ( ql(k,j,i) == 0.0 ) THEN597 k1 = 1.0 + 0.61* q(k,j,i)598 k2 = 0.61 * pt(k,j,i)596 IF ( ql(k,j,i) == 0.0_wp ) THEN 597 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 598 k2 = 0.61_wp * pt(k,j,i) 599 599 ELSE 600 600 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 601 601 temp = theta * t_d_pt(k) 602 k1 = ( 1.0 - q(k,j,i) + 1.61 *&603 ( q(k,j,i) - ql(k,j,i) ) * &604 ( 1.0 + 0.622* l_d_r / temp ) ) / &605 ( 1.0 + 0.622* l_d_r * l_d_cp * &602 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 603 ( q(k,j,i) - ql(k,j,i) ) * & 604 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 605 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 606 606 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 607 k2 = theta * ( l_d_cp / temp * k1 - 1.0 )607 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 608 608 ENDIF 609 609 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * & … … 613 613 ) * dd2zu(k) 614 614 ELSE IF ( cloud_droplets ) THEN 615 k1 = 1.0 + 0.61* q(k,j,i) - ql(k,j,i)616 k2 = 0.61 * pt(k,j,i)615 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 616 k2 = 0.61_wp * pt(k,j,i) 617 617 tend(k,j,i) = tend(k,j,i) - & 618 618 kh(k,j,i) * g / vpt(k,j,i) * & … … 634 634 635 635 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 636 k1 = 1.0 + 0.61* q(k,j,i)637 k2 = 0.61 * pt(k,j,i)636 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 637 k2 = 0.61_wp * pt(k,j,i) 638 638 ELSE IF ( cloud_physics ) THEN 639 IF ( ql(k,j,i) == 0.0 ) THEN640 k1 = 1.0 + 0.61* q(k,j,i)641 k2 = 0.61 * pt(k,j,i)639 IF ( ql(k,j,i) == 0.0_wp ) THEN 640 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 641 k2 = 0.61_wp * pt(k,j,i) 642 642 ELSE 643 643 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 644 644 temp = theta * t_d_pt(k) 645 k1 = ( 1.0 - q(k,j,i) + 1.61 *&646 ( q(k,j,i) - ql(k,j,i) ) * &647 ( 1.0 + 0.622* l_d_r / temp ) ) / &648 ( 1.0 + 0.622* l_d_r * l_d_cp * &645 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 646 ( q(k,j,i) - ql(k,j,i) ) * & 647 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 648 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 649 649 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 650 k2 = theta * ( l_d_cp / temp * k1 - 1.0 )650 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 651 651 ENDIF 652 652 ELSE IF ( cloud_droplets ) THEN 653 k1 = 1.0 + 0.61* q(k,j,i) - ql(k,j,i)654 k2 = 0.61 * pt(k,j,i)653 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 654 k2 = 0.61_wp * pt(k,j,i) 655 655 ENDIF 656 656 … … 668 668 669 669 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 670 k1 = 1.0 + 0.61* q(k,j,i)671 k2 = 0.61 * pt(k,j,i)670 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 671 k2 = 0.61_wp * pt(k,j,i) 672 672 ELSE IF ( cloud_physics ) THEN 673 IF ( ql(k,j,i) == 0.0 ) THEN674 k1 = 1.0 + 0.61* q(k,j,i)675 k2 = 0.61 * pt(k,j,i)673 IF ( ql(k,j,i) == 0.0_wp ) THEN 674 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 675 k2 = 0.61_wp * pt(k,j,i) 676 676 ELSE 677 677 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 678 678 temp = theta * t_d_pt(k) 679 k1 = ( 1.0 - q(k,j,i) + 1.61* &680 ( q(k,j,i) - ql(k,j,i) ) * &681 ( 1.0 + 0.622* l_d_r / temp ) ) / &682 ( 1.0 + 0.622* l_d_r * l_d_cp * &679 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 680 ( q(k,j,i) - ql(k,j,i) ) * & 681 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 682 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 683 683 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 684 k2 = theta * ( l_d_cp / temp * k1 - 1.0 )684 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 685 685 ENDIF 686 686 ELSE IF ( cloud_droplets ) THEN 687 k1 = 1.0 + 0.61* q(k,j,i) - ql(k,j,i)688 k2 = 0.61 * pt(k,j,i)687 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 688 k2 = 0.61_wp * pt(k,j,i) 689 689 ENDIF 690 690 … … 783 783 IF ( k >= nzb_diff_s_outer(j,i) ) THEN 784 784 785 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx786 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &787 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy788 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &789 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)790 791 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &792 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx793 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy794 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &795 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)796 797 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &798 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx799 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &800 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy801 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)802 803 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &785 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 786 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 787 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 788 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 789 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 790 791 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 792 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 793 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 794 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 795 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 796 797 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 798 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 799 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 800 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 801 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 802 803 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 804 804 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 805 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )806 807 IF ( def < 0.0 ) def = 0.0805 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 806 807 IF ( def < 0.0_wp ) def = 0.0_wp 808 808 809 809 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 825 825 DO k = 1, nzt 826 826 827 IF ( ( wall_e_x(j,i) /= 0.0 ).OR.( wall_e_y(j,i) /= 0.0) ) &827 IF ( ( wall_e_x(j,i) /= 0.0_wp ).OR.( wall_e_y(j,i) /= 0.0_wp ) ) & 828 828 THEN 829 829 830 830 IF ( k == nzb_diff_s_inner(j,i) - 1 ) THEN 831 831 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 832 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &833 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)832 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 833 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) 834 834 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 835 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &836 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)835 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 836 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) 837 837 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 838 838 839 IF ( wall_e_y(j,i) /= 0.0 ) THEN839 IF ( wall_e_y(j,i) /= 0.0_wp ) THEN 840 840 ! 841 841 !-- Inconsistency removed: as the thermal stratification is … … 852 852 ! wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp ) 853 853 km_neutral = kappa * & 854 ( usvs(k,j,i)**2 + wsvs(k,j,i)**2 )**0.25 * &855 0.5 * dy856 IF ( km_neutral > 0.0 ) THEN854 ( usvs(k,j,i)**2 + wsvs(k,j,i)**2 )**0.25_wp * & 855 0.5_wp * dy 856 IF ( km_neutral > 0.0_wp ) THEN 857 857 dudy = - wall_e_y(j,i) * usvs(k,j,i) / km_neutral 858 858 dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km_neutral 859 859 ELSE 860 dudy = 0.0 861 dwdy = 0.0 860 dudy = 0.0_wp 861 dwdy = 0.0_wp 862 862 ENDIF 863 863 ELSE 864 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &865 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy866 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &867 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy864 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 865 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 866 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 867 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 868 868 ENDIF 869 869 870 IF ( wall_e_x(j,i) /= 0.0 ) THEN870 IF ( wall_e_x(j,i) /= 0.0_wp ) THEN 871 871 ! 872 872 !-- Inconsistency removed: as the thermal stratification is … … 883 883 ! wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp ) 884 884 km_neutral = kappa * & 885 ( vsus(k,j,i)**2 + wsus(k,j,i)**2 )**0.25 * &886 0.5 * dx887 IF ( km_neutral > 0.0 ) THEN885 ( vsus(k,j,i)**2 + wsus(k,j,i)**2 )**0.25_wp * & 886 0.5_wp * dx 887 IF ( km_neutral > 0.0_wp ) THEN 888 888 dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km_neutral 889 889 dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km_neutral 890 890 ELSE 891 dvdx = 0.0 892 dwdx = 0.0 891 dvdx = 0.0_wp 892 dwdx = 0.0_wp 893 893 ENDIF 894 894 ELSE 895 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &896 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx897 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &898 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx895 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 896 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 897 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 898 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 899 899 ENDIF 900 900 901 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &901 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 902 902 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 903 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )904 905 IF ( def < 0.0 ) def = 0.0903 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 904 905 IF ( def < 0.0_wp ) def = 0.0_wp 906 906 907 907 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 919 919 920 920 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 921 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &922 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)923 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy924 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &925 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)921 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 922 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 923 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 924 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 925 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 926 926 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 927 927 928 IF ( wall_e_y(j,i) /= 0.0 ) THEN928 IF ( wall_e_y(j,i) /= 0.0_wp ) THEN 929 929 ! 930 930 !-- Inconsistency removed: as the thermal stratification … … 937 937 !-- validation has been available 938 938 km_neutral = kappa * ( usvs(k,j,i)**2 + & 939 wsvs(k,j,i)**2 )**0.25 * 0.5* dy940 IF ( km_neutral > 0.0 ) THEN939 wsvs(k,j,i)**2 )**0.25_wp * 0.5_wp * dy 940 IF ( km_neutral > 0.0_wp ) THEN 941 941 dudy = - wall_e_y(j,i) * usvs(k,j,i) / km_neutral 942 942 dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km_neutral 943 943 ELSE 944 dudy = 0.0 945 dwdy = 0.0 944 dudy = 0.0_wp 945 dwdy = 0.0_wp 946 946 ENDIF 947 947 ELSE 948 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &949 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy950 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &951 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy948 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 949 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 950 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 951 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 952 952 ENDIF 953 953 954 IF ( wall_e_x(j,i) /= 0.0 ) THEN954 IF ( wall_e_x(j,i) /= 0.0_wp ) THEN 955 955 ! 956 956 !-- Inconsistency removed: as the thermal stratification … … 963 963 !-- validation has been available 964 964 km_neutral = kappa * ( vsus(k,j,i)**2 + & 965 wsus(k,j,i)**2 )**0.25 * 0.5* dx966 IF ( km_neutral > 0.0 ) THEN965 wsus(k,j,i)**2 )**0.25_wp * 0.5_wp * dx 966 IF ( km_neutral > 0.0_wp ) THEN 967 967 dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km_neutral 968 968 dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km_neutral 969 969 ELSE 970 dvdx = 0.0 971 dwdx = 0.0 970 dvdx = 0.0_wp 971 dwdx = 0.0_wp 972 972 ENDIF 973 973 ELSE 974 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &975 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx976 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &977 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx974 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 975 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 976 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 977 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 978 978 ENDIF 979 979 980 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &980 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 981 981 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 982 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )983 984 IF ( def < 0.0 ) def = 0.0982 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 983 984 IF ( def < 0.0_wp ) def = 0.0_wp 985 985 986 986 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 993 993 IF ( k == nzb_diff_s_outer(j,i)-1 ) THEN 994 994 995 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx996 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &997 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy998 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &999 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)1000 1001 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &1002 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx1003 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy1004 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &1005 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)1006 1007 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &1008 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx1009 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &1010 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy1011 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)1012 1013 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &995 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 996 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 997 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 998 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 999 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 1000 1001 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1002 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1003 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1004 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1005 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 1006 1007 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1008 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1009 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1010 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1011 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1012 1013 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 1014 1014 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 1015 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )1016 1017 IF ( def < 0.0 ) def = 0.01015 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1016 1017 IF ( def < 0.0_wp ) def = 0.0_wp 1018 1018 1019 1019 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 1035 1035 DO k = 1, nzt 1036 1036 1037 IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0) ) &1037 IF ( ( wall_e_x(j,i) == 0.0_wp ) .AND. ( wall_e_y(j,i) == 0.0_wp ) ) & 1038 1038 THEN 1039 1039 1040 1040 IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN 1041 1041 1042 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx1043 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &1044 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy1045 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &1046 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)1047 1048 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &1049 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx1050 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy1051 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &1052 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)1053 1054 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &1055 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx1056 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &1057 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy1058 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)1059 1060 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &1042 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1043 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1044 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1045 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1046 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) 1047 1048 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1049 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1050 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1051 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1052 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) 1053 1054 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1055 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1056 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1057 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1058 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1059 1060 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 1061 1061 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 1062 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )1063 1064 IF ( def < 0.0 ) def = 0.01062 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1063 1064 IF ( def < 0.0_wp ) def = 0.0_wp 1065 1065 1066 1066 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 1082 1082 IF ( k == nzb_diff_s_outer(j,i)-1 ) THEN 1083 1083 1084 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx1085 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &1086 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy1087 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &1088 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)1089 1090 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &1091 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx1092 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy1093 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &1094 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)1095 1096 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &1097 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx1098 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &1099 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy1100 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)1101 1102 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &1084 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1085 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1086 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1087 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1088 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 1089 1090 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1091 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1092 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1093 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1094 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 1095 1096 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1097 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1098 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1099 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1100 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1101 1102 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 1103 1103 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 1104 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )1105 1106 IF ( def < 0.0 ) def = 0.01104 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1105 1106 IF ( def < 0.0_wp ) def = 0.0_wp 1107 1107 1108 1108 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 1232 1232 ! 1233 1233 ! IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1234 ! k1 = 1.0 + 0.61* q(k,j,i)1235 ! k2 = 0.61 * pt(k,j,i)1234 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1235 ! k2 = 0.61_wp * pt(k,j,i) 1236 1236 ! tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * & 1237 1237 ! g / vpt(k,j,i) * & … … 1240 1240 ! ) * dd2zu(k) 1241 1241 ! ELSE IF ( cloud_physics ) THEN 1242 ! IF ( ql(k,j,i) == 0.0 ) THEN1243 ! k1 = 1.0 + 0.61* q(k,j,i)1244 ! k2 = 0.61 * pt(k,j,i)1242 ! IF ( ql(k,j,i) == 0.0_wp ) THEN 1243 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1244 ! k2 = 0.61_wp * pt(k,j,i) 1245 1245 ! ELSE 1246 1246 ! theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1247 1247 ! temp = theta * t_d_pt(k) 1248 ! k1 = ( 1.0 - q(k,j,i) + 1.61* &1249 ! ( q(k,j,i) - ql(k,j,i) ) * &1250 ! ( 1.0 + 0.622* l_d_r / temp ) ) / &1251 ! ( 1.0 + 0.622* l_d_r * l_d_cp * &1248 ! k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1249 ! ( q(k,j,i) - ql(k,j,i) ) * & 1250 ! ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 1251 ! ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1252 1252 ! ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1253 ! k2 = theta * ( l_d_cp / temp * k1 - 1.0 )1253 ! k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1254 1254 ! ENDIF 1255 1255 ! tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * & … … 1259 1259 ! ) * dd2zu(k) 1260 1260 ! ELSE IF ( cloud_droplets ) THEN 1261 ! k1 = 1.0 + 0.61* q(k,j,i) - ql(k,j,i)1262 ! k2 = 0.61 * pt(k,j,i)1261 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1262 ! k2 = 0.61_wp * pt(k,j,i) 1263 1263 ! tend(k,j,i) = tend(k,j,i) - & 1264 1264 ! kh(k,j,i) * g / vpt(k,j,i) * & … … 1289 1289 ! 1290 1290 ! IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1291 ! k1 = 1.0 + 0.61* q(k,j,i)1292 ! k2 = 0.61 * pt(k,j,i)1291 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1292 ! k2 = 0.61_wp * pt(k,j,i) 1293 1293 ! ELSE IF ( cloud_physics ) THEN 1294 ! IF ( ql(k,j,i) == 0.0 ) THEN1295 ! k1 = 1.0 + 0.61* q(k,j,i)1296 ! k2 = 0.61 * pt(k,j,i)1294 ! IF ( ql(k,j,i) == 0.0_wp ) THEN 1295 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1296 ! k2 = 0.61_wp * pt(k,j,i) 1297 1297 ! ELSE 1298 1298 ! theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1299 1299 ! temp = theta * t_d_pt(k) 1300 ! k1 = ( 1.0 - q(k,j,i) + 1.61 *&1301 ! ( q(k,j,i) - ql(k,j,i) ) *&1302 ! ( 1.0 + 0.622 * l_d_r / temp ) ) /&1303 ! ( 1.0 + 0.622 * l_d_r * l_d_cp *&1300 ! k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1301 ! ( q(k,j,i) - ql(k,j,i) ) * & 1302 ! ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /& 1303 ! ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1304 1304 ! ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1305 ! k2 = theta * ( l_d_cp / temp * k1 - 1.0 )1305 ! k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1306 1306 ! ENDIF 1307 1307 ! ELSE IF ( cloud_droplets ) THEN 1308 ! k1 = 1.0 + 0.61* q(k,j,i) - ql(k,j,i)1309 ! k2 = 0.61 * pt(k,j,i)1308 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1309 ! k2 = 0.61_wp * pt(k,j,i) 1310 1310 ! ENDIF 1311 1311 ! … … 1330 1330 ! 1331 1331 ! IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1332 ! k1 = 1.0 + 0.61* q(k,j,i)1333 ! k2 = 0.61 * pt(k,j,i)1332 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1333 ! k2 = 0.61_wp * pt(k,j,i) 1334 1334 ! ELSE IF ( cloud_physics ) THEN 1335 ! IF ( ql(k,j,i) == 0.0 ) THEN1336 ! k1 = 1.0 + 0.61* q(k,j,i)1337 ! k2 = 0.61 * pt(k,j,i)1335 ! IF ( ql(k,j,i) == 0.0_wp ) THEN 1336 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1337 ! k2 = 0.61_wp * pt(k,j,i) 1338 1338 ! ELSE 1339 1339 ! theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1340 1340 ! temp = theta * t_d_pt(k) 1341 ! k1 = ( 1.0 - q(k,j,i) + 1.61 *&1342 ! ( q(k,j,i) - ql(k,j,i) ) *&1343 ! ( 1.0 + 0.622 * l_d_r / temp ) ) /&1344 ! ( 1.0 + 0.622 * l_d_r * l_d_cp *&1341 ! k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1342 ! ( q(k,j,i) - ql(k,j,i) ) * & 1343 ! ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /& 1344 ! ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1345 1345 ! ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1346 ! k2 = theta * ( l_d_cp / temp * k1 - 1.0 )1346 ! k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1347 1347 ! ENDIF 1348 1348 ! ELSE IF ( cloud_droplets ) THEN 1349 ! k1 = 1.0 + 0.61* q(k,j,i) - ql(k,j,i)1350 ! k2 = 0.61 * pt(k,j,i)1349 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1350 ! k2 = 0.61_wp * pt(k,j,i) 1351 1351 ! ENDIF 1352 1352 ! … … 1426 1426 DO k = nzb_diff_s_outer(j,i), nzt 1427 1427 1428 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx1429 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &1430 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy1431 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &1432 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)1433 1434 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &1435 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx1436 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy1437 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &1438 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)1439 1440 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &1441 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx1442 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &1443 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy1444 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)1445 1446 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) &1447 + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &1448 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )1449 1450 IF ( def < 0.0 ) def = 0.01428 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1429 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1430 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1431 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1432 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 1433 1434 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1435 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1436 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1437 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1438 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 1439 1440 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1441 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1442 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1443 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1444 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1445 1446 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) & 1447 + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 & 1448 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1449 1450 IF ( def < 0.0_wp ) def = 0.0_wp 1451 1451 1452 1452 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 1456 1456 IF ( prandtl_layer ) THEN 1457 1457 1458 IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0) ) THEN1458 IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) THEN 1459 1459 1460 1460 ! … … 1465 1465 1466 1466 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1467 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &1468 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)1469 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy1470 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &1471 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)1467 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1468 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) 1469 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1470 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1471 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) 1472 1472 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1473 1473 1474 IF ( wall_e_y(j,i) /= 0.0 ) THEN1474 IF ( wall_e_y(j,i) /= 0.0_wp ) THEN 1475 1475 ! 1476 1476 !-- Inconsistency removed: as the thermal stratification … … 1486 1486 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 1487 1487 wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp ) 1488 km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25 * &1489 0.5 * dy1490 IF ( km_neutral > 0.0 ) THEN1488 km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25_wp * & 1489 0.5_wp * dy 1490 IF ( km_neutral > 0.0_wp ) THEN 1491 1491 dudy = - wall_e_y(j,i) * usvs(k) / km_neutral 1492 1492 dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral 1493 1493 ELSE 1494 dudy = 0.0 1495 dwdy = 0.0 1494 dudy = 0.0_wp 1495 dwdy = 0.0_wp 1496 1496 ENDIF 1497 1497 ELSE 1498 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &1499 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy1500 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &1501 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy1498 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1499 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1500 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1501 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1502 1502 ENDIF 1503 1503 1504 IF ( wall_e_x(j,i) /= 0.0 ) THEN1504 IF ( wall_e_x(j,i) /= 0.0_wp ) THEN 1505 1505 ! 1506 1506 !-- Inconsistency removed: as the thermal stratification … … 1516 1516 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 1517 1517 wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp ) 1518 km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25 * &1519 0.5 * dx1520 IF ( km_neutral > 0.0 ) THEN1518 km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25_wp * & 1519 0.5_wp * dx 1520 IF ( km_neutral > 0.0_wp ) THEN 1521 1521 dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral 1522 1522 dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral 1523 1523 ELSE 1524 dvdx = 0.0 1525 dwdx = 0.0 1524 dvdx = 0.0_wp 1525 dwdx = 0.0_wp 1526 1526 ENDIF 1527 1527 ELSE 1528 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &1529 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx1530 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &1531 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx1528 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1529 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1530 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1531 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1532 1532 ENDIF 1533 1533 1534 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &1534 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 1535 1535 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 1536 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )1537 1538 IF ( def < 0.0 ) def = 0.01536 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1537 1538 IF ( def < 0.0_wp ) def = 0.0_wp 1539 1539 1540 1540 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 1549 1549 1550 1550 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1551 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &1552 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)1553 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy1554 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &1555 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)1551 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1552 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 1553 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1554 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1555 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 1556 1556 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1557 1557 1558 IF ( wall_e_y(j,i) /= 0.0 ) THEN1558 IF ( wall_e_y(j,i) /= 0.0_wp ) THEN 1559 1559 ! 1560 1560 !-- Inconsistency removed: as the thermal stratification … … 1567 1567 !-- validation has been available 1568 1568 km_neutral = kappa * ( usvs(k)**2 + & 1569 wsvs(k)**2 )**0.25 * 0.5* dy1570 IF ( km_neutral > 0.0 ) THEN1569 wsvs(k)**2 )**0.25_wp * 0.5_wp * dy 1570 IF ( km_neutral > 0.0_wp ) THEN 1571 1571 dudy = - wall_e_y(j,i) * usvs(k) / km_neutral 1572 1572 dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral 1573 1573 ELSE 1574 dudy = 0.0 1575 dwdy = 0.0 1574 dudy = 0.0_wp 1575 dwdy = 0.0_wp 1576 1576 ENDIF 1577 1577 ELSE 1578 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &1579 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy1580 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &1581 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy1578 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1579 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1580 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1581 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1582 1582 ENDIF 1583 1583 1584 IF ( wall_e_x(j,i) /= 0.0 ) THEN1584 IF ( wall_e_x(j,i) /= 0.0_wp ) THEN 1585 1585 ! 1586 1586 !-- Inconsistency removed: as the thermal stratification … … 1593 1593 !-- validation has been available 1594 1594 km_neutral = kappa * ( vsus(k)**2 + & 1595 wsus(k)**2 )**0.25 * 0.5* dx1596 IF ( km_neutral > 0.0 ) THEN1595 wsus(k)**2 )**0.25_wp * 0.5_wp * dx 1596 IF ( km_neutral > 0.0_wp ) THEN 1597 1597 dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral 1598 1598 dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral 1599 1599 ELSE 1600 dvdx = 0.0 1601 dwdx = 0.0 1600 dvdx = 0.0_wp 1601 dwdx = 0.0_wp 1602 1602 ENDIF 1603 1603 ELSE 1604 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &1605 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx1606 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &1607 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx1604 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1605 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1606 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1607 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1608 1608 ENDIF 1609 1609 1610 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &1610 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 1611 1611 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 1612 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )1613 1614 IF ( def < 0.0 ) def = 0.01612 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1613 1614 IF ( def < 0.0_wp ) def = 0.0_wp 1615 1615 1616 1616 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 1623 1623 k = nzb_diff_s_outer(j,i)-1 1624 1624 1625 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx1626 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &1627 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy1628 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &1629 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)1630 1631 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &1632 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx1633 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy1634 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &1635 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)1636 1637 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &1638 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx1639 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &1640 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy1641 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)1625 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1626 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1627 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1628 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1629 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 1630 1631 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1632 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1633 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1634 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1635 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 1636 1637 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1638 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1639 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1640 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1641 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1642 1642 1643 1643 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & … … 1645 1645 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1646 1646 1647 IF ( def < 0.0 ) def = 0.01647 IF ( def < 0.0_wp ) def = 0.0_wp 1648 1648 1649 1649 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 1657 1657 k = nzb_diff_s_inner(j,i)-1 1658 1658 1659 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1660 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1659 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1660 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1661 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1662 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1663 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) 1664 1665 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1666 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1667 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1668 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1669 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) 1670 1671 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1672 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1673 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1674 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1675 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1676 1677 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) & 1678 + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 & 1679 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1680 1681 IF ( def < 0.0_wp ) def = 0.0_wp 1682 1683 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 1684 1685 ENDIF 1686 1687 ELSEIF ( use_surface_fluxes ) THEN 1688 1689 k = nzb_diff_s_outer(j,i)-1 1690 1691 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1692 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1661 1693 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1662 dudz = 0.5* ( u(k+1,j,i) + u(k+1,j,i+1) - &1663 u _0(j,i) - u_0(j,i+1)) * dd2zu(k)1664 1665 dvdx = 0.25* ( v(k,j,i+1) + v(k,j+1,i+1) - &1694 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1695 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 1696 1697 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1666 1698 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1667 dvdy =( v(k,j+1,i) - v(k,j,i) ) * ddy1668 dvdz = 0.5* ( v(k+1,j,i) + v(k+1,j+1,i) - &1669 v _0(j,i) - v_0(j+1,i)) * dd2zu(k)1670 1671 dwdx = 0.25* ( w(k,j,i+1) + w(k-1,j,i+1) - &1699 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1700 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1701 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 1702 1703 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1672 1704 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1673 dwdy = 0.25* ( w(k,j+1,i) + w(k-1,j+1,i) - &1705 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1674 1706 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1675 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1676 1677 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) & 1678 + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 & 1679 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1680 1681 IF ( def < 0.0 ) def = 0.0 1682 1683 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 1684 1685 ENDIF 1686 1687 ELSEIF ( use_surface_fluxes ) THEN 1688 1689 k = nzb_diff_s_outer(j,i)-1 1690 1691 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1692 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1693 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1694 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1695 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 1696 1697 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1698 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1699 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1700 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1701 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 1702 1703 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1704 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1705 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1706 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1707 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1708 1709 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 1707 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1708 1709 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 1710 1710 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 1711 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )1712 1713 IF ( def < 0.0 ) def = 0.01711 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1712 1713 IF ( def < 0.0_wp ) def = 0.0_wp 1714 1714 1715 1715 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 1794 1794 1795 1795 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1796 k1 = 1.0 + 0.61* q(k,j,i)1797 k2 = 0.61 * pt(k,j,i)1796 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1797 k2 = 0.61_wp * pt(k,j,i) 1798 1798 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 1799 1799 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & … … 1801 1801 ) * dd2zu(k) 1802 1802 ELSE IF ( cloud_physics ) THEN 1803 IF ( ql(k,j,i) == 0.0 ) THEN1804 k1 = 1.0 + 0.61* q(k,j,i)1805 k2 = 0.61 * pt(k,j,i)1803 IF ( ql(k,j,i) == 0.0_wp ) THEN 1804 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1805 k2 = 0.61_wp * pt(k,j,i) 1806 1806 ELSE 1807 1807 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1808 1808 temp = theta * t_d_pt(k) 1809 k1 = ( 1.0 - q(k,j,i) + 1.61* &1810 ( q(k,j,i) - ql(k,j,i) ) * &1811 ( 1.0 + 0.622* l_d_r / temp ) ) / &1812 ( 1.0 + 0.622* l_d_r * l_d_cp * &1809 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1810 ( q(k,j,i) - ql(k,j,i) ) * & 1811 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 1812 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1813 1813 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1814 k2 = theta * ( l_d_cp / temp * k1 - 1.0 )1814 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1815 1815 ENDIF 1816 1816 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & … … 1819 1819 ) * dd2zu(k) 1820 1820 ELSE IF ( cloud_droplets ) THEN 1821 k1 = 1.0 + 0.61* q(k,j,i) - ql(k,j,i)1822 k2 = 0.61 * pt(k,j,i)1821 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1822 k2 = 0.61_wp * pt(k,j,i) 1823 1823 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 1824 1824 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & … … 1833 1833 1834 1834 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1835 k1 = 1.0 + 0.61* q(k,j,i)1836 k2 = 0.61 * pt(k,j,i)1835 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1836 k2 = 0.61_wp * pt(k,j,i) 1837 1837 ELSE IF ( cloud_physics ) THEN 1838 IF ( ql(k,j,i) == 0.0 ) THEN1839 k1 = 1.0 + 0.61* q(k,j,i)1840 k2 = 0.61 * pt(k,j,i)1838 IF ( ql(k,j,i) == 0.0_wp ) THEN 1839 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1840 k2 = 0.61_wp * pt(k,j,i) 1841 1841 ELSE 1842 1842 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1843 1843 temp = theta * t_d_pt(k) 1844 k1 = ( 1.0 - q(k,j,i) + 1.61* &1845 ( q(k,j,i) - ql(k,j,i) ) * &1846 ( 1.0 + 0.622* l_d_r / temp ) ) / &1847 ( 1.0 + 0.622* l_d_r * l_d_cp * &1844 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1845 ( q(k,j,i) - ql(k,j,i) ) * & 1846 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 1847 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1848 1848 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1849 k2 = theta * ( l_d_cp / temp * k1 - 1.0 )1849 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1850 1850 ENDIF 1851 1851 ELSE IF ( cloud_droplets ) THEN 1852 k1 = 1.0 + 0.61* q(k,j,i) - ql(k,j,i)1853 k2 = 0.61 * pt(k,j,i)1852 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1853 k2 = 0.61_wp * pt(k,j,i) 1854 1854 ENDIF 1855 1855 … … 1862 1862 1863 1863 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1864 k1 = 1.0 + 0.61* q(k,j,i)1865 k2 = 0.61 * pt(k,j,i)1864 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1865 k2 = 0.61_wp * pt(k,j,i) 1866 1866 ELSE IF ( cloud_physics ) THEN 1867 IF ( ql(k,j,i) == 0.0 ) THEN1868 k1 = 1.0 + 0.61* q(k,j,i)1869 k2 = 0.61 * pt(k,j,i)1867 IF ( ql(k,j,i) == 0.0_wp ) THEN 1868 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1869 k2 = 0.61_wp * pt(k,j,i) 1870 1870 ELSE 1871 1871 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1872 1872 temp = theta * t_d_pt(k) 1873 k1 = ( 1.0 - q(k,j,i) + 1.61* &1874 ( q(k,j,i) - ql(k,j,i) ) * &1875 ( 1.0 + 0.622* l_d_r / temp ) ) / &1876 ( 1.0 + 0.622* l_d_r * l_d_cp * &1873 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1874 ( q(k,j,i) - ql(k,j,i) ) * & 1875 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 1876 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1877 1877 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1878 k2 = theta * ( l_d_cp / temp * k1 - 1.0 )1878 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1879 1879 ENDIF 1880 1880 ELSE IF ( cloud_droplets ) THEN 1881 k1 = 1.0 + 0.61* q(k,j,i) - ql(k,j,i)1882 k2 = 0.61 * pt(k,j,i)1881 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1882 k2 = 0.61_wp * pt(k,j,i) 1883 1883 ENDIF 1884 1884 … … 1917 1917 IF ( first_call ) THEN 1918 1918 ALLOCATE( u_0(nysg:nyng,nxlg:nxrg), v_0(nysg:nyng,nxlg:nxrg) ) 1919 u_0 = 0.0 ! just to avoid access of uninitialized memory1920 v_0 = 0.0 ! within exchange_horiz_2d1919 u_0 = 0.0_wp ! just to avoid access of uninitialized memory 1920 v_0 = 0.0_wp ! within exchange_horiz_2d 1921 1921 first_call = .FALSE. 1922 1922 ENDIF … … 1942 1942 1943 1943 u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / & 1944 ( 0.5 * ( km(ku,j,i) + km(ku,j,i-1) ) +&1945 1.0E-20 )1944 ( 0.5_wp * ( km(ku,j,i) + km(ku,j,i-1) ) + & 1945 1.0E-20_wp ) 1946 1946 ! ( us(j,i) * kappa * zu(1) ) 1947 1947 v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / & 1948 ( 0.5 * ( km(kv,j,i) + km(kv,j-1,i) ) +&1949 1.0E-20 )1948 ( 0.5_wp * ( km(kv,j,i) + km(kv,j-1,i) ) + & 1949 1.0E-20_wp ) 1950 1950 ! ( us(j,i) * kappa * zu(1) ) 1951 1951 -
palm/trunk/SOURCE/random_function.f90
r1321 r1342 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 104 104 REAL(wp) :: rnmx !: 105 105 106 PARAMETER ( ia=16807, im=2147483647, am=1.0 /im, iq=127773, ir=2836, &107 ntab=32, ndiv=1+(im-1)/ntab, eps=1.2e-7 , rnmx=1.0-eps )106 PARAMETER ( ia=16807, im=2147483647, am=1.0_wp/im, iq=127773, ir=2836, & 107 ntab=32, ndiv=1+(im-1)/ntab, eps=1.2e-7_wp, rnmx=1.0_wp-eps ) 108 108 109 109 IF ( idum .le. 0 .or. random_iy .eq. 0 ) THEN -
palm/trunk/SOURCE/random_gauss.f90
r1321 r1342 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 78 78 79 79 IF ( iset == 0 ) THEN 80 rsq = 0.0 81 DO WHILE ( rsq >= 1.0 .OR. rsq == 0.0)82 v1 = 2.0 * random_function( idum ) - 1.083 v2 = 2.0 * random_function( idum ) - 1.080 rsq = 0.0_wp 81 DO WHILE ( rsq >= 1.0_wp .OR. rsq == 0.0_wp ) 82 v1 = 2.0_wp * random_function( idum ) - 1.0_wp 83 v2 = 2.0_wp * random_function( idum ) - 1.0_wp 84 84 rsq = v1**2 + v2**2 85 85 ENDDO 86 fac = SQRT( -2.0 * LOG( rsq ) / rsq )86 fac = SQRT( -2.0_wp * LOG( rsq ) / rsq ) 87 87 gset = v1 * fac 88 random_gauss = v2 * fac + 1.0 88 random_gauss = v2 * fac + 1.0_wp 89 89 iset = 1 90 90 ELSE 91 random_gauss = gset + 1.0 91 random_gauss = gset + 1.0_wp 92 92 iset = 0 93 93 ENDIF 94 94 95 IF ( ABS( random_gauss - 1.0 ) <= upper_limit ) EXIT95 IF ( ABS( random_gauss - 1.0_wp ) <= upper_limit ) EXIT 96 96 97 97 ENDDO -
palm/trunk/SOURCE/temperton_fft.f90
r1323 r1342 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! REAL constants defined as wp-kind 7 7 ! 8 8 ! Former revisions: … … 167 167 !OCL NOVREC 168 168 DO j = 1, nvex 169 a(i+inc) = 0.5 *a(i)169 a(i+inc) = 0.5_wp*a(i) 170 170 i = i + jump 171 171 END DO … … 173 173 i = istart + n*inc 174 174 DO j = 1, nvex 175 a(i) = 0.5 *a(i)175 a(i) = 0.5_wp*a(i) 176 176 i = i + jump 177 177 END DO … … 223 223 !OCL NOVREC 224 224 DO j = 1, nvex 225 a(ix) = 0.0 226 a(ix+inc) = 0.0 225 a(ix) = 0.0_wp 226 a(ix+inc) = 0.0_wp 227 227 ix = ix + jump 228 228 END DO … … 285 285 DO j = 1, nvex 286 286 a(ix) = a(ix+inc) 287 a(ix+inc) = 0.0 287 a(ix+inc) = 0.0_wp 288 288 ix = ix + jump 289 289 END DO … … 291 291 iz = istart + (n+1)*inc 292 292 DO j = 1, nvex 293 a(iz) = 0.0 293 a(iz) = 0.0_wp 294 294 iz = iz + jump 295 295 END DO … … 552 552 GO TO 170 553 553 30 CONTINUE 554 z = 1.0_wp/REAL(n )554 z = 1.0_wp/REAL(n,KIND=wp) 555 555 DO l = 1, la 556 556 i = ibase … … 590 590 DO ijk = 1, lot 591 591 c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i)) 592 c(jb+j) = a(ia+i) - 0.5 *(a(ib+i)+a(ic+i))592 c(jb+j) = a(ia+i) - 0.5_wp*(a(ib+i)+a(ic+i)) 593 593 d(jb+j) = sin60*(a(ic+i)-a(ib+i)) 594 594 i = i + inc3 … … 622 622 a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c2*a(ic+i)+s2*b(ic+i)) 623 623 b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c2*b(ic+i)-s2*a(ic+i)) 624 a2 = a(ia+i) - 0.5 *a1625 b2 = b(ia+i) - 0.5 *b1624 a2 = a(ia+i) - 0.5_wp*a1 625 b2 = b(ia+i) - 0.5_wp*b1 626 626 a3 = sin60*((c1*a(ib+i)+s1*b(ib+i))-(c2*a(ic+i)+s2*b(ic+i))) 627 627 b3 = sin60*((c1*b(ib+i)-s1*a(ib+i))-(c2*b(ic+i)-s2*a(ic+i))) … … 653 653 !OCL NOVREC 654 654 DO ijk = 1, lot 655 c(ja+j) = a(ia+i) + 0.5 *(a(ib+i)-a(ic+i))655 c(ja+j) = a(ia+i) + 0.5_wp*(a(ib+i)-a(ic+i)) 656 656 d(ja+j) = -sin60*(a(ib+i)+a(ic+i)) 657 657 c(jb+j) = a(ia+i) - (a(ib+i)-a(ic+i)) … … 665 665 GO TO 170 666 666 60 CONTINUE 667 z = 1.0_wp/REAL(n )667 z = 1.0_wp/REAL(n,KIND=wp) 668 668 zsin60 = z*sin60 669 669 DO l = 1, la … … 675 675 DO ijk = 1, lot 676 676 c(ja+j) = z*(a(ia+i)+(a(ib+i)+a(ic+i))) 677 c(jb+j) = z*(a(ia+i)-0.5 *(a(ib+i)+a(ic+i)))677 c(jb+j) = z*(a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i))) 678 678 d(jb+j) = zsin60*(a(ic+i)-a(ib+i)) 679 679 i = i + inc3 … … 794 794 GO TO 170 795 795 90 CONTINUE 796 z = 1.0_wp/REAL(n )796 z = 1.0_wp/REAL(n,KIND=wp) 797 797 DO l = 1, la 798 798 i = ibase … … 841 841 a2 = a(ic+i) + a(id+i) 842 842 a4 = a(ic+i) - a(id+i) 843 a5 = a(ia+i) - 0.25 *(a1+a2)843 a5 = a(ia+i) - 0.25_wp*(a1+a2) 844 844 a6 = qrt5*(a1-a2) 845 845 c(ja+j) = a(ia+i) + (a1+a2) … … 892 892 b2 = (c2*b(ic+i)-s2*a(ic+i)) + (c3*b(id+i)-s3*a(id+i)) 893 893 b4 = (c2*b(ic+i)-s2*a(ic+i)) - (c3*b(id+i)-s3*a(id+i)) 894 a5 = a(ia+i) - 0.25 *(a1+a2)894 a5 = a(ia+i) - 0.25_wp*(a1+a2) 895 895 a6 = qrt5*(a1-a2) 896 b5 = b(ia+i) - 0.25 *(b1+b2)896 b5 = b(ia+i) - 0.25_wp*(b1+b2) 897 897 b6 = qrt5*(b1-b2) 898 898 a10 = a5 + a6 … … 941 941 a2 = a(ic+i) + a(id+i) 942 942 a4 = a(ic+i) - a(id+i) 943 a5 = a(ia+i) + 0.25 *(a3-a4)943 a5 = a(ia+i) + 0.25_wp*(a3-a4) 944 944 a6 = qrt5*(a3+a4) 945 945 c(ja+j) = a5 + a6 … … 957 957 GO TO 170 958 958 120 CONTINUE 959 z = 1.0_wp/REAL(n )959 z = 1.0_wp/REAL(n,KIND=wp) 960 960 zqrt5 = z*qrt5 961 961 zsin36 = z*sin36 … … 972 972 a2 = a(ic+i) + a(id+i) 973 973 a4 = a(ic+i) - a(id+i) 974 a5 = z*(a(ia+i)-0.25 *(a1+a2))974 a5 = z*(a(ia+i)-0.25_wp*(a1+a2)) 975 975 a6 = zqrt5*(a1-a2) 976 976 c(ja+j) = z*(a(ia+i)+(a1+a2)) … … 1014 1014 a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i)) 1015 1015 c(ja+j) = (a(ia+i)+a(id+i)) + a11 1016 c(jc+j) = (a(ia+i)+a(id+i)-0.5 *a11)1016 c(jc+j) = (a(ia+i)+a(id+i)-0.5_wp*a11) 1017 1017 d(jc+j) = sin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i))) 1018 1018 a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i)) 1019 c(jb+j) = (a(ia+i)-a(id+i)) - 0.5 *a111019 c(jb+j) = (a(ia+i)-a(id+i)) - 0.5_wp*a11 1020 1020 d(jb+j) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i))) 1021 1021 c(jd+j) = (a(ia+i)-a(id+i)) + a11 … … 1071 1071 b5 = c5*b(if+i) - s5*a(if+i) 1072 1072 a11 = (a2+a5) + (a1+a4) 1073 a20 = (a(ia+i)+a3) - 0.5 *a111073 a20 = (a(ia+i)+a3) - 0.5_wp*a11 1074 1074 a21 = sin60*((a2+a5)-(a1+a4)) 1075 1075 b11 = (b2+b5) + (b1+b4) 1076 b20 = (b(ia+i)+b3) - 0.5 *b111076 b20 = (b(ia+i)+b3) - 0.5_wp*b11 1077 1077 b21 = sin60*((b2+b5)-(b1+b4)) 1078 1078 c(ja+j) = (a(ia+i)+a3) + a11 … … 1083 1083 d(je+j) = a21 - b20 1084 1084 a11 = (a2-a5) + (a4-a1) 1085 a20 = (a(ia+i)-a3) - 0.5 *a111085 a20 = (a(ia+i)-a3) - 0.5_wp*a11 1086 1086 a21 = sin60*((a4-a1)-(a2-a5)) 1087 1087 b11 = (b5-b2) - (b4-b1) 1088 b20 = (b3-b(ia+i)) - 0.5 *b111088 b20 = (b3-b(ia+i)) - 0.5_wp*b11 1089 1089 b21 = sin60*((b5-b2)+(b4-b1)) 1090 1090 c(jb+j) = a20 - b21 … … 1118 1118 !OCL NOVREC 1119 1119 DO ijk = 1, lot 1120 c(ja+j) = (a(ia+i)+0.5 *(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i))1121 d(ja+j) = -(a(id+i)+0.5 *(a(ib+i)+a(if+i))) - sin60*(a(ic+i)+a(ie+i))1120 c(ja+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i)) 1121 d(ja+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) - sin60*(a(ic+i)+a(ie+i)) 1122 1122 c(jb+j) = a(ia+i) - (a(ic+i)-a(ie+i)) 1123 1123 d(jb+j) = a(id+i) - (a(ib+i)+a(if+i)) 1124 c(jc+j) = (a(ia+i)+0.5 *(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i))1125 d(jc+j) = -(a(id+i)+0.5 *(a(ib+i)+a(if+i))) + sin60*(a(ic+i)+a(ie+i))1124 c(jc+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i)) 1125 d(jc+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) + sin60*(a(ic+i)+a(ie+i)) 1126 1126 i = i + inc3 1127 1127 j = j + inc4 … … 1133 1133 GO TO 170 1134 1134 150 CONTINUE 1135 z = 1.0_wp/REAL(n )1135 z = 1.0_wp/REAL(n,KIND=wp) 1136 1136 zsin60 = z*sin60 1137 1137 DO l = 1, la … … 1144 1144 a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i)) 1145 1145 c(ja+j) = z*((a(ia+i)+a(id+i))+a11) 1146 c(jc+j) = z*((a(ia+i)+a(id+i))-0.5 *a11)1146 c(jc+j) = z*((a(ia+i)+a(id+i))-0.5_wp*a11) 1147 1147 d(jc+j) = zsin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i))) 1148 1148 a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i)) 1149 c(jb+j) = z*((a(ia+i)-a(id+i))-0.5 *a11)1149 c(jb+j) = z*((a(ia+i)-a(id+i))-0.5_wp*a11) 1150 1150 d(jb+j) = zsin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i))) 1151 1151 c(jd+j) = z*((a(ia+i)-a(id+i))+a11) … … 1176 1176 jd = jc + 2*m*inc2 1177 1177 je = jd + 2*m*inc2 1178 z = 1.0_wp/REAL(n )1178 z = 1.0_wp/REAL(n,KIND=wp) 1179 1179 zsin45 = z*SQRT(0.5_wp) 1180 1180 … … 1419 1419 !OCL NOVREC 1420 1420 DO ijk = 1, lot 1421 c(ja+j) = 2.0 *(a(ia+i)+a(ib+i))1422 c(jb+j) = 2.0 *(a(ia+i)-a(ib+i))1421 c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i)) 1422 c(jb+j) = 2.0_wp*(a(ia+i)-a(ib+i)) 1423 1423 i = i + inc3 1424 1424 j = j + inc4 … … 1449 1449 DO ijk = 1, lot 1450 1450 c(ja+j) = a(ia+i) + a(ib+i) 1451 c(jb+j) = (a(ia+i)-0.5 *a(ib+i)) - (sin60*(b(ib+i)))1452 c(jc+j) = (a(ia+i)-0.5 *a(ib+i)) + (sin60*(b(ib+i)))1451 c(jb+j) = (a(ia+i)-0.5_wp*a(ib+i)) - (sin60*(b(ib+i))) 1452 c(jc+j) = (a(ia+i)-0.5_wp*a(ib+i)) + (sin60*(b(ib+i))) 1453 1453 i = i + inc3 1454 1454 j = j + inc4 … … 1481 1481 c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i)) 1482 1482 d(ja+j) = b(ia+i) + (b(ib+i)-b(ic+i)) 1483 c(jb+j) = c1*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ & 1484 & b(ic+i)))) - s1*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- & 1483 c(jb+j) = c1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ & 1484 & b(ic+i)))) & 1485 & - s1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- & 1485 1486 & a(ic+i)))) 1486 d(jb+j) = s1*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ & 1487 & b(ic+i)))) + c1*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- & 1487 d(jb+j) = s1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ & 1488 & b(ic+i)))) & 1489 & + c1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- & 1488 1490 & a(ic+i)))) 1489 c(jc+j) = c2*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ & 1490 & b(ic+i)))) - s2*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- & 1491 c(jc+j) = c2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ & 1492 & b(ic+i)))) & 1493 & - s2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- & 1491 1494 & a(ic+i)))) 1492 d(jc+j) = s2*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ & 1493 & b(ic+i)))) + c2*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- & 1495 d(jc+j) = s2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ & 1496 & b(ic+i)))) & 1497 & + c2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- & 1494 1498 & a(ic+i)))) 1495 1499 i = i + inc3 … … 1515 1519 DO ijk = 1, lot 1516 1520 c(ja+j) = a(ia+i) + a(ib+i) 1517 c(jb+j) = (0.5 *a(ia+i)-a(ib+i)) - (sin60*b(ia+i))1518 c(jc+j) = -(0.5 *a(ia+i)-a(ib+i)) - (sin60*b(ia+i))1521 c(jb+j) = (0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i)) 1522 c(jc+j) = -(0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i)) 1519 1523 i = i + inc3 1520 1524 j = j + inc4 … … 1526 1530 GO TO 170 1527 1531 60 CONTINUE 1528 ssin60 = 2.0 *sin601529 DO l = 1, la 1530 i = ibase 1531 j = jbase 1532 !DIR$ IVDEP 1533 !CDIR NODEP 1534 !OCL NOVREC 1535 DO ijk = 1, lot 1536 c(ja+j) = 2.0 *(a(ia+i)+a(ib+i))1537 c(jb+j) = (2.0 *a(ia+i)-a(ib+i)) - (ssin60*b(ib+i))1538 c(jc+j) = (2.0 *a(ia+i)-a(ib+i)) + (ssin60*b(ib+i))1532 ssin60 = 2.0_wp*sin60 1533 DO l = 1, la 1534 i = ibase 1535 j = jbase 1536 !DIR$ IVDEP 1537 !CDIR NODEP 1538 !OCL NOVREC 1539 DO ijk = 1, lot 1540 c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i)) 1541 c(jb+j) = (2.0_wp*a(ia+i)-a(ib+i)) - (ssin60*b(ib+i)) 1542 c(jc+j) = (2.0_wp*a(ia+i)-a(ib+i)) + (ssin60*b(ib+i)) 1539 1543 i = i + inc3 1540 1544 j = j + inc4 … … 1659 1663 !OCL NOVREC 1660 1664 DO ijk = 1, lot 1661 c(ja+j) = 2.0 *((a(ia+i)+a(ic+i))+a(ib+i))1662 c(jb+j) = 2.0 *((a(ia+i)-a(ic+i))-b(ib+i))1663 c(jc+j) = 2.0 *((a(ia+i)+a(ic+i))-a(ib+i))1664 c(jd+j) = 2.0 *((a(ia+i)-a(ic+i))+b(ib+i))1665 c(ja+j) = 2.0_wp*((a(ia+i)+a(ic+i))+a(ib+i)) 1666 c(jb+j) = 2.0_wp*((a(ia+i)-a(ic+i))-b(ib+i)) 1667 c(jc+j) = 2.0_wp*((a(ia+i)+a(ic+i))-a(ib+i)) 1668 c(jd+j) = 2.0_wp*((a(ia+i)-a(ic+i))+b(ib+i)) 1665 1669 i = i + inc3 1666 1670 j = j + inc4 … … 1695 1699 DO ijk = 1, lot 1696 1700 c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i)) 1697 c(jb+j) = ((a(ia+i)-0.25 *(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) - &1701 c(jb+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) - & 1698 1702 & (sin72*b(ib+i)+sin36*b(ic+i)) 1699 c(jc+j) = ((a(ia+i)-0.25 *(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) - &1703 c(jc+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) - & 1700 1704 & (sin36*b(ib+i)-sin72*b(ic+i)) 1701 c(jd+j) = ((a(ia+i)-0.25 *(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) + &1705 c(jd+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) + & 1702 1706 & (sin36*b(ib+i)-sin72*b(ic+i)) 1703 c(je+j) = ((a(ia+i)-0.25 *(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) + &1707 c(je+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) + & 1704 1708 & (sin72*b(ib+i)+sin36*b(ic+i)) 1705 1709 i = i + inc3 … … 1740 1744 DO ijk = 1, lot 1741 1745 1742 a10(ijk) = (a(ia+i)-0.25 *((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) + &1746 a10(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) + & 1743 1747 & qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i))) 1744 a20(ijk) = (a(ia+i)-0.25 *((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) - &1748 a20(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) - & 1745 1749 & qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i))) 1746 b10(ijk) = (b(ia+i)-0.25 *((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) + &1750 b10(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) + & 1747 1751 & qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i))) 1748 b20(ijk) = (b(ia+i)-0.25 *((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) - &1752 b20(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) - & 1749 1753 & qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i))) 1750 1754 a11(ijk) = sin72*(b(ib+i)+b(ie+i)) + sin36*(b(ic+i)+b(id+i)) … … 1788 1792 DO ijk = 1, lot 1789 1793 c(ja+j) = (a(ia+i)+a(ib+i)) + a(ic+i) 1790 c(jb+j) = (qrt5*(a(ia+i)-a(ib+i))+(0.25 *(a(ia+i)+a(ib+i))-a(ic+i))) - &1794 c(jb+j) = (qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - & 1791 1795 & (sin36*b(ia+i)+sin72*b(ib+i)) 1792 c(je+j) = -(qrt5*(a(ia+i)-a(ib+i))+(0.25 *(a(ia+i)+a(ib+i))-a(ic+i))) - &1796 c(je+j) = -(qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - & 1793 1797 & (sin36*b(ia+i)+sin72*b(ib+i)) 1794 c(jc+j) = (qrt5*(a(ia+i)-a(ib+i))-(0.25 *(a(ia+i)+a(ib+i))-a(ic+i))) - &1798 c(jc+j) = (qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - & 1795 1799 & (sin72*b(ia+i)-sin36*b(ib+i)) 1796 c(jd+j) = -(qrt5*(a(ia+i)-a(ib+i))-(0.25 *(a(ia+i)+a(ib+i))-a(ic+i))) - &1800 c(jd+j) = -(qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - & 1797 1801 & (sin72*b(ia+i)-sin36*b(ib+i)) 1798 1802 i = i + inc3 … … 1815 1819 !OCL NOVREC 1816 1820 DO ijk = 1, lot 1817 c(ja+j) = 2.0 *(a(ia+i)+(a(ib+i)+a(ic+i)))1818 c(jb+j) = (2.0 *(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &1821 c(ja+j) = 2.0_wp*(a(ia+i)+(a(ib+i)+a(ic+i))) 1822 c(jb+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ & 1819 1823 & i))) - (ssin72*b(ib+i)+ssin36*b(ic+i)) 1820 c(jc+j) = (2.0 *(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &1824 c(jc+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ & 1821 1825 & i))) - (ssin36*b(ib+i)-ssin72*b(ic+i)) 1822 c(jd+j) = (2.0 *(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &1826 c(jd+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ & 1823 1827 & i))) + (ssin36*b(ib+i)-ssin72*b(ic+i)) 1824 c(je+j) = (2.0 *(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &1828 c(je+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ & 1825 1829 & i))) + (ssin72*b(ib+i)+ssin36*b(ic+i)) 1826 1830 i = i + inc3 … … 1859 1863 c(ja+j) = (a(ia+i)+a(id+i)) + (a(ib+i)+a(ic+i)) 1860 1864 c(jd+j) = (a(ia+i)-a(id+i)) - (a(ib+i)-a(ic+i)) 1861 c(jb+j) = ((a(ia+i)-a(id+i))+0.5 *(a(ib+i)-a(ic+i))) - (sin60*(b(ib+ &1865 c(jb+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) - (sin60*(b(ib+ & 1862 1866 & i)+b(ic+i))) 1863 c(jf+j) = ((a(ia+i)-a(id+i))+0.5 *(a(ib+i)-a(ic+i))) + (sin60*(b(ib+ &1867 c(jf+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) + (sin60*(b(ib+ & 1864 1868 & i)+b(ic+i))) 1865 c(jc+j) = ((a(ia+i)+a(id+i))-0.5 *(a(ib+i)+a(ic+i))) - (sin60*(b(ib+ &1869 c(jc+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) - (sin60*(b(ib+ & 1866 1870 & i)-b(ic+i))) 1867 c(je+j) = ((a(ia+i)+a(id+i))-0.5 *(a(ib+i)+a(ic+i))) + (sin60*(b(ib+ &1871 c(je+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) + (sin60*(b(ib+ & 1868 1872 & i)-b(ic+i))) 1869 1873 i = i + inc3 … … 1909 1913 1910 1914 a11(ijk) = (a(ie+i)+a(ib+i)) + (a(ic+i)+a(if+i)) 1911 a20(ijk) = (a(ia+i)+a(id+i)) - 0.5 *a11(ijk)1915 a20(ijk) = (a(ia+i)+a(id+i)) - 0.5_wp*a11(ijk) 1912 1916 a21(ijk) = sin60*((a(ie+i)+a(ib+i))-(a(ic+i)+a(if+i))) 1913 1917 b11(ijk) = (b(ib+i)-b(ie+i)) + (b(ic+i)-b(if+i)) 1914 b20(ijk) = (b(ia+i)-b(id+i)) - 0.5 *b11(ijk)1918 b20(ijk) = (b(ia+i)-b(id+i)) - 0.5_wp*b11(ijk) 1915 1919 b21(ijk) = sin60*((b(ib+i)-b(ie+i))-(b(ic+i)-b(if+i))) 1916 1920 … … 1924 1928 a11(ijk) = (a(ie+i)-a(ib+i)) + (a(ic+i)-a(if+i)) 1925 1929 b11(ijk) = (b(ie+i)+b(ib+i)) - (b(ic+i)+b(if+i)) 1926 a20(ijk) = (a(ia+i)-a(id+i)) - 0.5 *a11(ijk)1930 a20(ijk) = (a(ia+i)-a(id+i)) - 0.5_wp*a11(ijk) 1927 1931 a21(ijk) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i))) 1928 b20(ijk) = (b(ia+i)+b(id+i)) + 0.5 *b11(ijk)1932 b20(ijk) = (b(ia+i)+b(id+i)) + 0.5_wp*b11(ijk) 1929 1933 b21(ijk) = sin60*((b(ie+i)+b(ib+i))+(b(ic+i)+b(if+i))) 1930 1934 … … 1964 1968 c(ja+j) = a(ib+i) + (a(ia+i)+a(ic+i)) 1965 1969 c(jd+j) = b(ib+i) - (b(ia+i)+b(ic+i)) 1966 c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5 *(b(ia+i)+b(ic+i))+b(ib+i))1967 c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5 *(b(ia+i)+b(ic+i))+b(ib+i))1968 c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5 *(a(ia+i)+a(ic+i))-a(ib+i))1969 c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5 *(a(ia+i)+a(ic+i))-a(ib+i))1970 c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i)) 1971 c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i)) 1972 c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i)) 1973 c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i)) 1970 1974 i = i + inc3 1971 1975 j = j + inc4 … … 1977 1981 GO TO 170 1978 1982 150 CONTINUE 1979 ssin60 = 2.0 *sin601980 DO l = 1, la 1981 i = ibase 1982 j = jbase 1983 !DIR$ IVDEP 1984 !CDIR NODEP 1985 !OCL NOVREC 1986 DO ijk = 1, lot 1987 c(ja+j) = (2.0 *(a(ia+i)+a(id+i))) + (2.0*(a(ib+i)+a(ic+i)))1988 c(jd+j) = (2.0 *(a(ia+i)-a(id+i))) - (2.0*(a(ib+i)-a(ic+i)))1989 c(jb+j) = (2.0 *(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) - (ssin60*(b(ib+ &1983 ssin60 = 2.0_wp*sin60 1984 DO l = 1, la 1985 i = ibase 1986 j = jbase 1987 !DIR$ IVDEP 1988 !CDIR NODEP 1989 !OCL NOVREC 1990 DO ijk = 1, lot 1991 c(ja+j) = (2.0_wp*(a(ia+i)+a(id+i))) + (2.0_wp*(a(ib+i)+a(ic+i))) 1992 c(jd+j) = (2.0_wp*(a(ia+i)-a(id+i))) - (2.0_wp*(a(ib+i)-a(ic+i))) 1993 c(jb+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) - (ssin60*(b(ib+ & 1990 1994 & i)+b(ic+i))) 1991 c(jf+j) = (2.0 *(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) + (ssin60*(b(ib+ &1995 c(jf+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) + (ssin60*(b(ib+ & 1992 1996 & i)+b(ic+i))) 1993 c(jc+j) = (2.0 *(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) - (ssin60*(b(ib+ &1997 c(jc+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) - (ssin60*(b(ib+ & 1994 1998 & i)-b(ic+i))) 1995 c(je+j) = (2.0 *(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) + (ssin60*(b(ib+ &1999 c(je+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) + (ssin60*(b(ib+ & 1996 2000 & i)-b(ic+i))) 1997 2001 i = i + inc3 … … 2030 2034 !OCL NOVREC 2031 2035 DO ijk = 1, lot 2032 c(ja+j) = 2.0 *(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i)))2033 c(je+j) = 2.0 *(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i)))2034 c(jc+j) = 2.0 *(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i)))2035 c(jg+j) = 2.0 *(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i)))2036 c(jb+j) = 2.0 *((a(ia+i)-a(ie+i))-b(ic+i)) + ssin45*((a(ib+i)-a(id+ &2037 & i))-(b(ib+i)+b(id+i)))2038 c(jf+j) = 2.0 *((a(ia+i)-a(ie+i))-b(ic+i)) - ssin45*((a(ib+i)-a(id+ &2039 & i))-(b(ib+i)+b(id+i)))2040 c(jd+j) = 2.0 *((a(ia+i)-a(ie+i))+b(ic+i)) - ssin45*((a(ib+i)-a(id+ &2041 & i))+(b(ib+i)+b(id+i)))2042 c(jh+j) = 2.0 *((a(ia+i)-a(ie+i))+b(ic+i)) + ssin45*((a(ib+i)-a(id+ &2043 & i))+(b(ib+i)+b(id+i)))2036 c(ja+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i))) 2037 c(je+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i))) 2038 c(jc+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i))) 2039 c(jg+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i))) 2040 c(jb+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) + ssin45*((a(ib+i)-a(id+ & 2041 & i))-(b(ib+i)+b(id+i))) 2042 c(jf+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) - ssin45*((a(ib+i)-a(id+ & 2043 & i))-(b(ib+i)+b(id+i))) 2044 c(jd+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) - ssin45*((a(ib+i)-a(id+ & 2045 & i))+(b(ib+i)+b(id+i))) 2046 c(jh+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) + ssin45*((a(ib+i)-a(id+ & 2047 & i))+(b(ib+i)+b(id+i))) 2044 2048 i = i + inc3 2045 2049 j = j + inc4 … … 2114 2118 ixxx = 1 2115 2119 2116 del = 4.0_wp*ASIN(1.0_wp)/REAL(n )2120 del = 4.0_wp*ASIN(1.0_wp)/REAL(n,KIND=wp) 2117 2121 nil = 0 2118 2122 nhl = (n/2) - 1 2119 2123 DO k = nil, nhl 2120 angle = REAL(k )*del2124 angle = REAL(k,KIND=wp)*del 2121 2125 trigs(2*k+1) = COS(angle) 2122 2126 trigs(2*k+2) = SIN(angle) -
palm/trunk/SOURCE/time_integration.f90
r1321 r1342 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 205 205 !-- At the beginning of a simulation determine the time step as well as 206 206 !-- determine and print out the run control parameters 207 IF ( simulated_time == 0.0 ) CALL timestep207 IF ( simulated_time == 0.0_wp ) CALL timestep 208 208 209 209 CALL run_control … … 221 221 time_coupling = time_coupling - dt_coupling 222 222 ENDDO 223 IF (time_coupling == 0.0 .AND. time_since_reference_point < dt_coupling)&223 IF (time_coupling == 0.0_wp .AND. time_since_reference_point < dt_coupling)& 224 224 THEN 225 225 time_coupling = time_since_reference_point … … 243 243 ! 244 244 !-- Determine size of next time step 245 IF ( simulated_time /= 0.0 ) THEN245 IF ( simulated_time /= 0.0_wp ) THEN 246 246 CALL timestep 247 247 ENDIF … … 717 717 ! 718 718 !-- If required, call flow_statistics for averaging in time 719 IF ( averaging_interval_pr /= 0.0 .AND. &719 IF ( averaging_interval_pr /= 0.0_wp .AND. & 720 720 ( dt_dopr - time_dopr ) <= averaging_interval_pr .AND. & 721 721 simulated_time >= skip_time_dopr ) THEN … … 731 731 ! 732 732 !-- Sum-up 3d-arrays for later output of time-averaged 2d/3d/masked data 733 IF ( averaging_interval /= 0.0 .AND. &733 IF ( averaging_interval /= 0.0_wp .AND. & 734 734 ( dt_data_output_av - time_do_av ) <= averaging_interval .AND. & 735 735 simulated_time >= skip_time_data_output_av ) & … … 745 745 ! 746 746 !-- Calculate spectra for time averaging 747 IF ( averaging_interval_sp /= 0.0 .AND. &747 IF ( averaging_interval_sp /= 0.0_wp .AND. & 748 748 ( dt_dosp - time_dosp ) <= averaging_interval_sp .AND. & 749 749 simulated_time >= skip_time_dosp ) THEN … … 782 782 IF ( dopr_n /= 0 ) CALL data_output_profiles 783 783 time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) ) 784 time_dopr_av = 0.0 ! due to averaging (see above)784 time_dopr_av = 0.0_wp ! due to averaging (see above) 785 785 ENDIF 786 786 -
palm/trunk/SOURCE/time_to_string.f90
r1321 r1342 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 61 61 ! 62 62 !-- Calculate the number of hours, minutes, and seconds 63 hours = INT( time / 3600.0 )64 rest_time = time - hours * 3600 65 minutes = INT( rest_time / 60.0 )63 hours = INT( time / 3600.0_wp ) 64 rest_time = time - hours * 3600_wp 65 minutes = INT( rest_time / 60.0_wp ) 66 66 seconds = rest_time - minutes * 60 67 67 -
palm/trunk/SOURCE/timestep.f90
r1323 r1342 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 153 153 !-- Horizontal averages already existent, just need to average them 154 154 !-- vertically. 155 u_gtrans = 0.0 156 v_gtrans = 0.0 155 u_gtrans = 0.0_wp 156 v_gtrans = 0.0_wp 157 157 DO k = nzb+1, nzt 158 158 u_gtrans = u_gtrans + hom(k,1,1,0) … … 164 164 ! 165 165 !-- Averaging over the entire model domain. 166 u_gtrans_l = 0.0 167 v_gtrans_l = 0.0 166 u_gtrans_l = 0.0_wp 167 v_gtrans_l = 0.0_wp 168 168 !$acc parallel present( u, v ) 169 169 DO i = nxl, nxr … … 196 196 #if defined( __openacc ) 197 197 IF ( dt_fixed ) THEN ! otherwise do it further below for better cache usage 198 u_max_l = -999999.9 199 u_min_l = 999999.9 200 v_max_l = -999999.9 201 v_min_l = 999999.9 202 w_max_l = -999999.9 203 w_min_l = 999999.9 198 u_max_l = -999999.9_wp 199 u_min_l = 999999.9_wp 200 v_max_l = -999999.9_wp 201 v_min_l = 999999.9_wp 202 w_max_l = -999999.9_wp 203 w_min_l = 999999.9_wp 204 204 !$acc parallel present( u, v, w ) 205 205 DO i = nxl, nxr … … 266 266 !-- Calculate the maximum time step according to the CFL-criterion, 267 267 !-- individually for each velocity component 268 dt_u_l = 999999.9 269 dt_v_l = 999999.9 270 dt_w_l = 999999.9 271 u_max_l = -999999.9 272 u_min_l = 999999.9 273 v_max_l = -999999.9 274 v_min_l = 999999.9 275 w_max_l = -999999.9 276 w_min_l = 999999.9 268 dt_u_l = 999999.9_wp 269 dt_v_l = 999999.9_wp 270 dt_w_l = 999999.9_wp 271 u_max_l = -999999.9_wp 272 u_min_l = 999999.9_wp 273 v_max_l = -999999.9_wp 274 v_min_l = 999999.9_wp 275 w_max_l = -999999.9_wp 276 w_min_l = 999999.9_wp 277 277 !$acc parallel loop collapse(3) present( u, v, w ) 278 278 DO i = nxl, nxr 279 279 DO j = nys, nyn 280 280 DO k = nzb+1, nzt 281 dt_u_l = MIN( dt_u_l, ( dx / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10 ) ) )282 dt_v_l = MIN( dt_v_l, ( dy / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10 ) ) )283 dt_w_l = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) ) + 1.0E-10 ) ) )281 dt_u_l = MIN( dt_u_l, ( dx / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10_wp ) ) ) 282 dt_v_l = MIN( dt_v_l, ( dy / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10_wp ) ) ) 283 dt_w_l = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) ) + 1.0E-10_wp ) ) ) 284 284 u_max_l = MAX( u_max_l, u(k,j,i) ) 285 285 u_min_l = MIN( u_min_l, u(k,j,i) ) … … 346 346 !-- Calculate the maximum time step according to the CFL-criterion, 347 347 !-- individually for each velocity component 348 dt_u_l = 999999.9 349 dt_v_l = 999999.9 350 dt_w_l = 999999.9 348 dt_u_l = 999999.9_wp 349 dt_v_l = 999999.9_wp 350 dt_w_l = 999999.9_wp 351 351 DO i = nxl, nxr 352 352 DO j = nys, nyn 353 353 DO k = nzb+1, nzt 354 dt_u_l = MIN( dt_u_l, ( dx / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10 ) ) )355 dt_v_l = MIN( dt_v_l, ( dy / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10 ) ) )356 dt_w_l = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) ) + 1.0E-10 ) ) )354 dt_u_l = MIN( dt_u_l, ( dx / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10_wp ) ) ) 355 dt_v_l = MIN( dt_v_l, ( dy / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10_wp ) ) ) 356 dt_w_l = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) ) + 1.0E-10_wp ) ) ) 357 357 ENDDO 358 358 ENDDO … … 382 382 !-- in the Prandtl-layer friction term only dz/2 is used. 383 383 !-- Experience from the old model seems to justify this. 384 dt_diff_l = 999999.0 384 dt_diff_l = 999999.0_wp 385 385 386 386 DO k = nzb+1, nzt 387 dxyz2_min(k) = MIN( dx2, dy2, dzw(k)*dzw(k) ) * 0.125 387 dxyz2_min(k) = MIN( dx2, dy2, dzw(k)*dzw(k) ) * 0.125_wp 388 388 ENDDO 389 389 … … 395 395 DO k = nzb+1, nzt 396 396 dt_diff_l = MIN( dt_diff_l, dxyz2_min(k) / & 397 ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20 ) )397 ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20_wp ) ) 398 398 ENDDO 399 399 ENDDO … … 414 414 IF ( plant_canopy ) THEN 415 415 416 dt_plant_canopy_l = 0.0 416 dt_plant_canopy_l = 0.0_wp 417 417 DO i = nxl, nxr 418 418 DO j = nys, nyn … … 424 424 v(k,j+1,i) + & 425 425 v(k,j+1,i-1) ) & 426 / 4.0 )**2+ &426 / 4.0_wp )**2 + & 427 427 ( ( w(k-1,j,i-1) + & 428 428 w(k-1,j,i) + & 429 429 w(k,j,i-1) + & 430 430 w(k,j,i) ) & 431 / 4.0 )**2 )431 / 4.0_wp )**2 ) 432 432 IF ( dt_plant_canopy_u > dt_plant_canopy_l ) THEN 433 433 dt_plant_canopy_l = dt_plant_canopy_u … … 438 438 u(k,j,i) + & 439 439 u(k,j,i+1) ) & 440 / 4.0 )**2+ &440 / 4.0_wp )**2 + & 441 441 v(k,j,i)**2 + & 442 442 ( ( w(k-1,j-1,i) + & … … 444 444 w(k,j-1,i) + & 445 445 w(k,j,i) ) & 446 / 4.0 )**2 )446 / 4.0_wp )**2 ) 447 447 IF ( dt_plant_canopy_v > dt_plant_canopy_l ) THEN 448 448 dt_plant_canopy_l = dt_plant_canopy_v … … 453 453 u(k+1,j,i) + & 454 454 u(k+1,j,i+1) ) & 455 / 4.0 )**2+ &455 / 4.0_wp )**2 + & 456 456 ( ( v(k,j,i) + & 457 457 v(k,j+1,i) + & 458 458 v(k+1,j,i) + & 459 459 v(k+1,j+1,i) ) & 460 / 4.0 )**2 + &460 / 4.0_wp )**2 + & 461 461 w(k,j,i)**2 ) 462 462 IF ( dt_plant_canopy_w > dt_plant_canopy_l ) THEN … … 467 467 ENDDO 468 468 469 IF ( dt_plant_canopy_l > 0.0 ) THEN469 IF ( dt_plant_canopy_l > 0.0_wp ) THEN 470 470 ! 471 471 !-- Invert dt_plant_canopy_l and apply a security timestep factor 0.1 472 dt_plant_canopy_l = 0.1 / dt_plant_canopy_l472 dt_plant_canopy_l = 0.1_wp / dt_plant_canopy_l 473 473 ELSE 474 474 ! … … 515 515 ! 516 516 !-- Set flag if the time step becomes too small. 517 IF ( dt_3d < ( 0.00001 * dt_max ) ) THEN517 IF ( dt_3d < ( 0.00001_wp * dt_max ) ) THEN 518 518 stop_dt = .TRUE. 519 519 … … 553 553 ! 554 554 !-- Ensure a smooth value (two significant digits) of the timestep. 555 div = 1000.0 555 div = 1000.0_wp 556 556 DO WHILE ( dt_3d < div ) 557 div = div / 10.0 557 div = div / 10.0_wp 558 558 ENDDO 559 dt_3d = NINT( dt_3d * 100.0 / div ) * div / 100.0559 dt_3d = NINT( dt_3d * 100.0_wp / div ) * div / 100.0_wp 560 560 561 561 ! -
palm/trunk/SOURCE/tridia_solver.f90
r1323 r1342 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 110 110 ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1) 111 111 ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1) 112 ddzuw(k,3) = -1.0 * &112 ddzuw(k,3) = -1.0_wp * & 113 113 ( ddzu_pres(k+2) * ddzw(k+1) + ddzu_pres(k+1) * ddzw(k+1) ) 114 114 ENDDO … … 172 172 IF ( j >= 0 .AND. j <= nnyh ) THEN 173 173 IF ( i >= 0 .AND. i <= nnxh ) THEN 174 ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0* pi * i ) / &175 176 2.0 * ( 1.0 - COS( ( 2.0* pi * j ) / &177 174 ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / & 175 REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + & 176 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / & 177 REAL( ny+1, KIND=wp ) ) ) / ( dy * dy ) 178 178 ELSE 179 ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0* pi * ( nx+1-i ) ) / &180 181 2.0 * ( 1.0 - COS( ( 2.0* pi * j ) / &182 179 ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / & 180 REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + & 181 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / & 182 REAL( ny+1, KIND=wp ) ) ) / ( dy * dy ) 183 183 ENDIF 184 184 ELSE 185 185 IF ( i >= 0 .AND. i <= nnxh ) THEN 186 ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0* pi * i ) / &187 188 2.0 * ( 1.0 - COS( ( 2.0* pi * ( ny+1-j ) ) / &189 186 ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / & 187 REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + & 188 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / & 189 REAL( ny+1, KIND=wp ) ) ) / ( dy * dy ) 190 190 ELSE 191 ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0* pi * ( nx+1-i ) ) / &192 193 2.0 * ( 1.0 - COS( ( 2.0* pi * ( ny+1-j ) ) / &194 191 ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / & 192 REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + & 193 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / & 194 REAL( ny+1, KIND=wp ) ) ) / ( dy * dy ) 195 195 ENDIF 196 196 ENDIF … … 283 283 284 284 IF ( k == nz-1 ) THEN 285 ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20 )285 ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20_wp ) 286 286 ELSE 287 287 ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) & … … 301 301 !$acc kernels loop present( ar ) 302 302 DO k = 1, nz 303 ar(nxl_z,nys_z,k) = 0.0 303 ar(nxl_z,nys_z,k) = 0.0_wp 304 304 ENDDO 305 305 !$acc end kernels loop … … 367 367 368 368 IF ( k == nz-1 ) THEN 369 ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20 )369 ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20_wp ) 370 370 ELSE 371 371 ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) & … … 385 385 !$acc kernels loop present( ar ) 386 386 DO k = 1, nz 387 ar(nxl_z,nys_z,k) = 0.0 387 ar(nxl_z,nys_z,k) = 0.0_wp 388 388 ENDDO 389 389 ENDIF … … 560 560 DO i = 0, nx 561 561 IF ( i >= 0 .AND. i <= nnxh ) THEN 562 l(i) = 2.0 * ( 1.0 - COS( ( 2.0* pi * i ) / &563 564 2.0 * ( 1.0 - COS( ( 2.0* pi * j ) / &565 562 l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / & 563 REAL( nx+1, KIND=wp ) ) ) * ddx2 + & 564 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / & 565 REAL( ny+1, KIND=wp ) ) ) * ddy2 566 566 ELSE 567 l(i) = 2.0 * ( 1.0 - COS( ( 2.0* pi * ( nx+1-i ) ) / &568 569 2.0 * ( 1.0 - COS( ( 2.0* pi * j ) / &570 567 l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / & 568 REAL( nx+1, KIND=wp ) ) ) * ddx2 + & 569 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / & 570 REAL( ny+1, KIND=wp ) ) ) * ddy2 571 571 ENDIF 572 572 ENDDO … … 574 574 DO k = 0, nz-1 575 575 DO i = 0, nx 576 a = -1.0 * ddzu_pres(k+2) * ddzw(k+1)577 c = -1.0 * ddzu_pres(k+1) * ddzw(k+1)576 a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1) 577 c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1) 578 578 tri_for_1d(1,i,k) = a + c - l(i) 579 579 ENDDO … … 660 660 !-- the model domain. 661 661 DO i = 0, nx 662 ar(i,nz) = ar1(i,nz-1) / ( tri_for_1d(4,i,nz-1) + 1.0E-20 )662 ar(i,nz) = ar1(i,nz-1) / ( tri_for_1d(4,i,nz-1) + 1.0E-20_wp ) 663 663 ENDDO 664 664 DO k = nz-2, 0, -1 … … 676 676 IF ( j == 0 ) THEN 677 677 DO k = 1, nz 678 ar(0,k) = 0.0 678 ar(0,k) = 0.0_wp 679 679 ENDDO 680 680 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
