Changeset 2508 for palm/trunk/SOURCE/production_e.f90
- Timestamp:
- Oct 2, 2017 8:57:09 AM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/production_e.f90 (modified) (29 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/production_e.f90
r2478 r2508 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - Bugfix in buoyancy production term, wrong base state was set. 28 ! - Consider use_single_reference_value case if humidity is used. 29 ! - In case of use_top_fluxes, use correct inverse density at model top 30 ! - Consider use_surface_fluxes and use_top_fluxes in ocean case 31 ! 32 ! 2478 2017-09-18 13:37:24Z suehring 27 33 ! Bugfix, consider case where no constant-flux layer and no surfaces fluxes 28 34 ! are used … … 152 158 humidity, kappa, neutral, ocean, pt_reference, & 153 159 rho_reference, use_single_reference_value, & 154 use_surface_fluxes, use_top_fluxes 160 use_surface_fluxes, use_top_fluxes, vpt_reference 155 161 156 162 USE grid_variables, & … … 372 378 DO m = surf_s, surf_e 373 379 k = surf_lsm_h%k(m) 374 ! 375 !-- Please note, actually, an interpolation of u_0 and v_0 376 !-- onto the grid center would be required. However, this 377 !-- would require several data transfers between 2D-grid and 378 !-- wall type. The effect of this missing interpolation is 379 !-- negligible. (See also production_e_init). 380 380 381 dudz(k,j) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k) 381 382 dvdz(k,j) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k) … … 388 389 DO m = surf_s, surf_e 389 390 k = surf_usm_h%k(m) 390 ! 391 !-- Please note, actually, an interpolation of u_0 and v_0 392 !-- onto the grid center would be required. However, this 393 !-- would require several data transfers between 2D-grid and 394 !-- wall type. The effect of this missing interpolation is 395 !-- negligible. (See also production_e_init). 391 396 392 dudz(k,j) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k) 397 393 dvdz(k,j) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k) … … 405 401 DO m = surf_s, surf_e 406 402 k = surf_def_h(1)%k(m) 407 ! 408 !-- Please note, actually, an interpolation of u_0 and v_0 409 !-- onto the grid center would be required. However, this 410 !-- would require several data transfers between 2D-grid and 411 !-- wall type. The effect of this missing interpolation is 412 !-- negligible. (See also production_e_init). 403 413 404 dudz(k,j) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k) 414 405 dvdz(k,j) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k) 415 406 416 407 ENDDO 417 418 408 ENDDO 419 409 … … 493 483 IF ( .NOT. humidity ) THEN 494 484 495 IF ( use_single_reference_value ) THEN 496 497 IF ( ocean ) THEN 498 ! 499 !-- So far in the ocean no special treatment of density flux 500 !-- in the bottom and top surface layer 501 DO j = nys, nyn 502 DO k = nzb+1, nzt 503 tend(k,j,i) = tend(k,j,i) + & 504 kh(k,j,i) * g / rho_reference * & 505 ( prho(k+1,j,i) - prho(k-1,j,i) ) * & 506 dd2zu(k) * & 507 MERGE( 1.0_wp, 0.0_wp, & 508 BTEST( wall_flags_0(k,j,i), 0 ) & 509 ) 485 IF ( ocean ) THEN 486 ! 487 !-- So far in the ocean no special treatment of density flux 488 !-- in the bottom and top surface layer 489 DO j = nys, nyn 490 DO k = nzb+1, nzt 491 tend(k,j,i) = tend(k,j,i) + & 492 kh(k,j,i) * g / & 493 MERGE( rho_reference, prho(k,j,i), & 494 use_single_reference_value ) * & 495 ( prho(k+1,j,i) - prho(k-1,j,i) ) * & 496 dd2zu(k) * & 497 MERGE( 1.0_wp, 0.0_wp, & 498 BTEST( wall_flags_0(k,j,i), 30 ) & 499 ) * & 500 MERGE( 1.0_wp, 0.0_wp, & 501 BTEST( wall_flags_0(k,j,i), 9 ) & 502 ) 503 ENDDO 504 ! 505 !-- Treatment of near-surface grid points, at up- and down- 506 !-- ward facing surfaces 507 IF ( use_surface_fluxes ) THEN 508 DO l = 0, 1 509 surf_s = surf_def_h(l)%start_index(j,i) 510 surf_e = surf_def_h(l)%end_index(j,i) 511 DO m = surf_s, surf_e 512 k = surf_def_h(l)%k(m) 513 tend(k,j,i) = tend(k,j,i) + g / & 514 MERGE( rho_reference, prho(k,j,i), & 515 use_single_reference_value ) * & 516 drho_air_zw(k-1) * & 517 surf_def_h(l)%shf(m) 518 ENDDO 510 519 ENDDO 520 521 ENDIF 522 523 IF ( use_top_fluxes ) THEN 524 surf_s = surf_def_h(2)%start_index(j,i) 525 surf_e = surf_def_h(2)%end_index(j,i) 526 DO m = surf_s, surf_e 527 k = surf_def_h(2)%k(m) 528 tend(k,j,i) = tend(k,j,i) + g / & 529 MERGE( rho_reference, prho(k,j,i), & 530 use_single_reference_value ) * & 531 drho_air_zw(k) * & 532 surf_def_h(2)%shf(m) 533 ENDDO 534 ENDIF 535 536 ENDDO 537 538 ELSE 539 540 DO j = nys, nyn 541 DO k = nzb+1, nzt 542 ! 543 !-- Flag 9 is used to mask top fluxes, flag 30 to mask 544 !-- surface fluxes 545 tend(k,j,i) = tend(k,j,i) - & 546 kh(k,j,i) * g / & 547 MERGE( pt_reference, pt(k,j,i), & 548 use_single_reference_value ) * & 549 ( pt(k+1,j,i) - pt(k-1,j,i) ) * & 550 dd2zu(k) * & 551 MERGE( 1.0_wp, 0.0_wp, & 552 BTEST( wall_flags_0(k,j,i), 30 ) & 553 ) * & 554 MERGE( 1.0_wp, 0.0_wp, & 555 BTEST( wall_flags_0(k,j,i), 9 ) & 556 ) 511 557 ENDDO 512 558 513 ELSE 514 515 DO j = nys, nyn 516 DO k = nzb+1, nzt 517 ! 518 !-- Flag 9 is used to mask top fluxes, flag 30 to mask 519 !-- surface fluxes 520 tend(k,j,i) = tend(k,j,i) - & 521 kh(k,j,i) * g / pt_reference * & 522 ( pt(k+1,j,i) - pt(k-1,j,i) ) * & 523 dd2zu(k) * & 524 MERGE( 1.0_wp, 0.0_wp, & 525 BTEST( wall_flags_0(k,j,i), 30 ) & 526 ) * & 527 MERGE( 1.0_wp, 0.0_wp, & 528 BTEST( wall_flags_0(k,j,i), 9 ) & 529 ) 559 IF ( use_surface_fluxes ) THEN 560 ! 561 !-- Default surfaces, up- and downward-facing 562 DO l = 0, 1 563 surf_s = surf_def_h(l)%start_index(j,i) 564 surf_e = surf_def_h(l)%end_index(j,i) 565 DO m = surf_s, surf_e 566 k = surf_def_h(l)%k(m) 567 tend(k,j,i) = tend(k,j,i) + g / & 568 MERGE( pt_reference, pt(k,j,i), & 569 use_single_reference_value ) & 570 * drho_air_zw(k-1) & 571 * surf_def_h(l)%shf(m) 572 ENDDO 530 573 ENDDO 531 532 IF ( use_surface_fluxes ) THEN 533 ! 534 !-- Default surfaces, up- and downward-facing 535 DO l = 0, 1 536 surf_s = surf_def_h(l)%start_index(j,i) 537 surf_e = surf_def_h(l)%end_index(j,i) 538 DO m = surf_s, surf_e 539 k = surf_def_h(l)%k(m) 540 tend(k,j,i) = tend(k,j,i) + g / pt_reference & 541 * drho_air_zw(k-1) & 542 * surf_def_h(l)%shf(m) 543 ENDDO 544 ENDDO 545 ! 546 !-- Natural surfaces 547 surf_s = surf_lsm_h%start_index(j,i) 548 surf_e = surf_lsm_h%end_index(j,i) 549 DO m = surf_s, surf_e 550 k = surf_lsm_h%k(m) 551 tend(k,j,i) = tend(k,j,i) + g / pt_reference & 552 * drho_air_zw(k-1) & 553 * surf_lsm_h%shf(m) 554 ENDDO 555 ! 556 !-- Urban surfaces 557 surf_s = surf_usm_h%start_index(j,i) 558 surf_e = surf_usm_h%end_index(j,i) 559 DO m = surf_s, surf_e 560 k = surf_usm_h%k(m) 561 tend(k,j,i) = tend(k,j,i) + g / pt_reference & 562 * drho_air_zw(k-1) & 563 * surf_usm_h%shf(m) 564 ENDDO 565 ENDIF 566 567 IF ( use_top_fluxes ) THEN 568 surf_s = surf_def_h(2)%start_index(j,i) 569 surf_e = surf_def_h(2)%end_index(j,i) 570 DO m = surf_s, surf_e 571 k = surf_def_h(2)%k(m) 572 tend(k,j,i) = tend(k,j,i) + g / pt_reference & 573 * drho_air_zw(k-1) & 574 * surf_def_h(2)%shf(m) 575 ENDDO 576 ENDIF 577 ENDDO 578 579 ENDIF 580 581 ELSE 582 583 IF ( ocean ) THEN 584 ! 585 !-- So far in the ocean no special treatment of density flux 586 !-- in the bottom and top surface layer 587 DO j = nys, nyn 588 DO k = nzb+1, nzt 589 tend(k,j,i) = tend(k,j,i) + & 590 kh(k,j,i) * g / prho(k,j,i) * & 591 ( prho(k+1,j,i) - prho(k-1,j,i) ) * & 592 dd2zu(k) * & 593 MERGE( 1.0_wp, 0.0_wp, & 594 BTEST( wall_flags_0(k,j,i), 0 ) & 595 ) 574 ! 575 !-- Natural surfaces 576 surf_s = surf_lsm_h%start_index(j,i) 577 surf_e = surf_lsm_h%end_index(j,i) 578 DO m = surf_s, surf_e 579 k = surf_lsm_h%k(m) 580 tend(k,j,i) = tend(k,j,i) + g / & 581 MERGE( pt_reference, pt(k,j,i), & 582 use_single_reference_value ) & 583 * drho_air_zw(k-1) & 584 * surf_lsm_h%shf(m) 596 585 ENDDO 597 ENDDO 598 599 ELSE 600 601 DO j = nys, nyn 602 DO k = nzb+1, nzt 603 ! 604 !-- Flag 9 is used to mask top fluxes, flag 30 to mask 605 !-- surface fluxes 606 tend(k,j,i) = tend(k,j,i) - & 607 kh(k,j,i) * g / pt(k,j,i) * & 608 ( pt(k+1,j,i) - pt(k-1,j,i) ) * & 609 dd2zu(k) * & 610 MERGE( 1.0_wp, 0.0_wp, & 611 BTEST( wall_flags_0(k,j,i), 30 ) & 612 ) * & 613 MERGE( 1.0_wp, 0.0_wp, & 614 BTEST( wall_flags_0(k,j,i), 9 ) & 615 ) 586 ! 587 !-- Urban surfaces 588 surf_s = surf_usm_h%start_index(j,i) 589 surf_e = surf_usm_h%end_index(j,i) 590 DO m = surf_s, surf_e 591 k = surf_usm_h%k(m) 592 tend(k,j,i) = tend(k,j,i) + g / & 593 MERGE( pt_reference, pt(k,j,i), & 594 use_single_reference_value ) & 595 * drho_air_zw(k-1) & 596 * surf_usm_h%shf(m) 597 ENDDO 598 ENDIF 599 600 IF ( use_top_fluxes ) THEN 601 surf_s = surf_def_h(2)%start_index(j,i) 602 surf_e = surf_def_h(2)%end_index(j,i) 603 DO m = surf_s, surf_e 604 k = surf_def_h(2)%k(m) 605 tend(k,j,i) = tend(k,j,i) + g / & 606 MERGE( pt_reference, pt(k,j,i), & 607 use_single_reference_value ) & 608 * drho_air_zw(k) & 609 * surf_def_h(2)%shf(m) 616 610 ENDDO 617 618 IF ( use_surface_fluxes ) THEN 619 ! 620 !-- Default surfaces, up- and downwrd-facing 621 DO l = 0, 1 622 surf_s = surf_def_h(l)%start_index(j,i) 623 surf_e = surf_def_h(l)%end_index(j,i) 624 DO m = surf_s, surf_e 625 k = surf_def_h(l)%k(m) 626 tend(k,j,i) = tend(k,j,i) + g / pt_reference & 627 * drho_air_zw(k-1) & 628 * surf_def_h(l)%shf(m) 629 ENDDO 630 ENDDO 631 ! 632 !-- Natural surfaces 633 surf_s = surf_lsm_h%start_index(j,i) 634 surf_e = surf_lsm_h%end_index(j,i) 635 DO m = surf_s, surf_e 636 k = surf_lsm_h%k(m) 637 tend(k,j,i) = tend(k,j,i) + g / pt_reference & 638 * drho_air_zw(k-1) & 639 * surf_lsm_h%shf(m) 640 ENDDO 641 ! 642 !-- Urban surfaces 643 surf_s = surf_usm_h%start_index(j,i) 644 surf_e = surf_usm_h%end_index(j,i) 645 DO m = surf_s, surf_e 646 k = surf_usm_h%k(m) 647 tend(k,j,i) = tend(k,j,i) + g / pt_reference & 648 * drho_air_zw(k-1) & 649 * surf_usm_h%shf(m) 650 ENDDO 651 ENDIF 652 653 IF ( use_top_fluxes ) THEN 654 surf_s = surf_def_h(2)%start_index(j,i) 655 surf_e = surf_def_h(2)%end_index(j,i) 656 DO m = surf_s, surf_e 657 k = surf_def_h(2)%k(m) 658 tend(k,j,i) = tend(k,j,i) + g / pt_reference & 659 * drho_air_zw(k-1) & 660 * surf_def_h(2)%shf(m) 661 ENDDO 662 ENDIF 663 ENDDO 664 665 ENDIF 611 ENDIF 612 ENDDO 666 613 667 614 ENDIF … … 679 626 k2 = 0.61_wp * pt(k,j,i) 680 627 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * & 681 g / vpt(k,j,i) * & 628 g / & 629 MERGE( vpt_reference, vpt(k,j,i), & 630 use_single_reference_value ) * & 682 631 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 683 632 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & … … 704 653 ENDIF 705 654 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * & 706 g / vpt(k,j,i) * & 655 g / & 656 MERGE( vpt_reference, vpt(k,j,i), & 657 use_single_reference_value ) * & 707 658 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 708 659 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & … … 718 669 k2 = 0.61_wp * pt(k,j,i) 719 670 tend(k,j,i) = tend(k,j,i) - & 720 kh(k,j,i) * g / vpt(k,j,i) * & 671 kh(k,j,i) * g / & 672 MERGE( vpt_reference, vpt(k,j,i), & 673 use_single_reference_value ) * & 721 674 ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) + & 722 675 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - & … … 739 692 DO j = nys, nyn 740 693 ! 741 !-- Treat horizontal default surfaces , up- and downward-facing694 !-- Treat horizontal default surfaces 742 695 DO l = 0, 1 743 696 surf_s = surf_def_h(l)%start_index(j,i) … … 768 721 ENDIF 769 722 770 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 723 tend(k,j,i) = tend(k,j,i) + g / & 724 MERGE( vpt_reference, vpt(k,j,i), & 725 use_single_reference_value ) * & 771 726 ( k1 * surf_def_h(l)%shf(m) + & 772 k2 * surf_def_h(l)%qsws(m) 727 k2 * surf_def_h(l)%qsws(m) & 773 728 ) * drho_air_zw(k-1) 774 729 ENDDO … … 803 758 ENDIF 804 759 805 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 760 tend(k,j,i) = tend(k,j,i) + g / & 761 MERGE( vpt_reference, vpt(k,j,i), & 762 use_single_reference_value ) * & 806 763 ( k1 * surf_lsm_h%shf(m) + & 807 764 k2 * surf_lsm_h%qsws(m) & … … 837 794 ENDIF 838 795 839 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 796 tend(k,j,i) = tend(k,j,i) + g / & 797 MERGE( vpt_reference, vpt(k,j,i), & 798 use_single_reference_value ) * & 840 799 ( k1 * surf_usm_h%shf(m) + & 841 800 k2 * surf_usm_h%qsws(m) & … … 878 837 ENDIF 879 838 880 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 839 tend(k,j,i) = tend(k,j,i) + g / & 840 MERGE( vpt_reference, vpt(k,j,i), & 841 use_single_reference_value ) * & 881 842 ( k1 * surf_def_h(2)%shf(m) + & 882 k2 * surf_def_h(2)%qsws(m) 883 ) * drho_air_zw(k -1)843 k2 * surf_def_h(2)%qsws(m) & 844 ) * drho_air_zw(k) 884 845 885 846 ENDDO … … 916 877 humidity, kappa, neutral, ocean, pt_reference, & 917 878 rho_reference, use_single_reference_value, & 918 use_surface_fluxes, use_top_fluxes 879 use_surface_fluxes, use_top_fluxes, vpt_reference 919 880 920 881 USE grid_variables, & … … 1065 1026 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 1066 1027 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 1067 BTEST( wall_flags_0(k,j,i-1), 0 ) ) 1028 BTEST( wall_flags_0(k,j,i-1), 0 ) ) 1068 1029 dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) 1069 1030 dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) … … 1129 1090 DO m = surf_s, surf_e 1130 1091 k = surf_lsm_h%k(m) 1131 ! 1132 !-- Please note, actually, an interpolation of u_0 and v_0 1133 !-- onto the grid center would be required. However, this 1134 !-- would require several data transfers between 2D-grid and 1135 !-- wall type. The effect of this missing interpolation is 1136 !-- negligible. (See also production_e_init). 1092 1137 1093 dudz(k) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k) 1138 1094 dvdz(k) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k) … … 1144 1100 DO m = surf_s, surf_e 1145 1101 k = surf_usm_h%k(m) 1146 ! 1147 !-- Please note, actually, an interpolation of u_0 and v_0 1148 !-- onto the grid center would be required. However, this 1149 !-- would require several data transfers between 2D-grid and 1150 !-- wall type. The effect of this missing interpolation is 1151 !-- negligible. (See also production_e_init). 1102 1152 1103 dudz(k) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k) 1153 1104 dvdz(k) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k) … … 1160 1111 DO m = surf_s, surf_e 1161 1112 k = surf_def_h(1)%k(m) 1162 ! 1163 !-- Please note, actually, an interpolation of u_0 and v_0 1164 !-- onto the grid center would be required. However, this 1165 !-- would require several data transfers between 2D-grid and 1166 !-- wall type. The effect of this missing interpolation is 1167 !-- negligible. (See also production_e_init). 1113 1168 1114 dudz(k) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k) 1169 1115 dvdz(k) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k) … … 1232 1178 IF ( .NOT. humidity ) THEN 1233 1179 1234 IF ( use_single_reference_value ) THEN 1235 1236 IF ( ocean ) THEN 1237 ! 1238 !-- So far in the ocean no special treatment of density flux in 1239 !-- the bottom and top surface layer 1240 DO k = nzb+1, nzt 1180 IF ( ocean ) THEN 1181 ! 1182 !-- So far in the ocean no special treatment of density flux in 1183 !-- the bottom and top surface layer 1184 DO k = nzb+1, nzt 1241 1185 1242 tend(k,j,i) = tend(k,j,i) + & 1243 kh(k,j,i) * g / rho_reference * & 1244 ( prho(k+1,j,i) - prho(k-1,j,i) ) * & 1245 dd2zu(k) * & 1246 MERGE( 1.0_wp, 0.0_wp, & 1247 BTEST( wall_flags_0(k,j,i), 0 ) & 1248 ) 1186 tend(k,j,i) = tend(k,j,i) + & 1187 kh(k,j,i) * g / & 1188 MERGE( rho_reference, prho(k,j,i), & 1189 use_single_reference_value ) * & 1190 ( prho(k+1,j,i) - prho(k-1,j,i) ) * & 1191 dd2zu(k) * & 1192 MERGE( 1.0_wp, 0.0_wp, & 1193 BTEST( wall_flags_0(k,j,i), 30 ) & 1194 ) * & 1195 MERGE( 1.0_wp, 0.0_wp, & 1196 BTEST( wall_flags_0(k,j,i), 9 ) & 1197 ) 1198 ENDDO 1199 1200 IF ( use_surface_fluxes ) THEN 1201 ! 1202 !-- Default surfaces, up- and downward-facing 1203 DO l = 0, 1 1204 surf_s = surf_def_h(l)%start_index(j,i) 1205 surf_e = surf_def_h(l)%end_index(j,i) 1206 DO m = surf_s, surf_e 1207 k = surf_def_h(l)%k(m) 1208 tend(k,j,i) = tend(k,j,i) + g / & 1209 MERGE( rho_reference, prho(k,j,i), & 1210 use_single_reference_value ) * & 1211 drho_air_zw(k-1) * & 1212 surf_def_h(l)%shf(m) 1213 ENDDO 1249 1214 ENDDO 1250 1215 1251 ELSE 1252 1253 DO k = nzb+1, nzt 1254 ! 1255 !-- Flag 9 is used to mask top fluxes, flag 30 to mask 1256 !-- surface fluxes 1257 tend(k,j,i) = tend(k,j,i) - & 1258 kh(k,j,i) * g / pt_reference * & 1259 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) * & 1260 MERGE( 1.0_wp, 0.0_wp, & 1261 BTEST( wall_flags_0(k,j,i), 30 ) & 1262 ) * & 1263 MERGE( 1.0_wp, 0.0_wp, & 1264 BTEST( wall_flags_0(k,j,i), 9 ) & 1265 ) 1266 1216 ENDIF 1217 1218 IF ( use_top_fluxes ) THEN 1219 surf_s = surf_def_h(2)%start_index(j,i) 1220 surf_e = surf_def_h(2)%end_index(j,i) 1221 DO m = surf_s, surf_e 1222 k = surf_def_h(2)%k(m) 1223 tend(k,j,i) = tend(k,j,i) + g / & 1224 MERGE( rho_reference, prho(k,j,i), & 1225 use_single_reference_value ) * & 1226 drho_air_zw(k) * & 1227 surf_def_h(2)%shf(m) 1267 1228 ENDDO 1268 1269 IF ( use_surface_fluxes ) THEN 1270 ! 1271 !-- Default surfaces, up- and downward-facing 1272 DO l = 0, 1 1273 surf_s = surf_def_h(l)%start_index(j,i) 1274 surf_e = surf_def_h(l)%end_index(j,i) 1275 DO m = surf_s, surf_e 1276 k = surf_def_h(l)%k(m) 1277 tend(k,j,i) = tend(k,j,i) + g / pt_reference * & 1278 drho_air_zw(k-1) * & 1279 surf_def_h(l)%shf(m) 1280 ENDDO 1229 ENDIF 1230 1231 ELSE 1232 1233 DO k = nzb+1, nzt 1234 ! 1235 !-- Flag 9 is used to mask top fluxes, flag 30 to mask 1236 !-- surface fluxes 1237 tend(k,j,i) = tend(k,j,i) - & 1238 kh(k,j,i) * g / & 1239 MERGE( pt_reference, pt(k,j,i), & 1240 use_single_reference_value ) * & 1241 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) * & 1242 MERGE( 1.0_wp, 0.0_wp, & 1243 BTEST( wall_flags_0(k,j,i), 30 ) & 1244 ) * & 1245 MERGE( 1.0_wp, 0.0_wp, & 1246 BTEST( wall_flags_0(k,j,i), 9 ) & 1247 ) 1248 1249 ENDDO 1250 1251 IF ( use_surface_fluxes ) THEN 1252 ! 1253 !-- Default surfaces, up- and downward-facing 1254 DO l = 0, 1 1255 surf_s = surf_def_h(l)%start_index(j,i) 1256 surf_e = surf_def_h(l)%end_index(j,i) 1257 DO m = surf_s, surf_e 1258 k = surf_def_h(l)%k(m) 1259 tend(k,j,i) = tend(k,j,i) + g / & 1260 MERGE( pt_reference, pt(k,j,i), & 1261 use_single_reference_value ) * & 1262 drho_air_zw(k-1) * & 1263 surf_def_h(l)%shf(m) 1281 1264 ENDDO 1282 ! 1283 !-- Natural surfaces 1284 surf_s = surf_lsm_h%start_index(j,i) 1285 surf_e = surf_lsm_h%end_index(j,i) 1286 DO m = surf_s, surf_e 1287 k = surf_lsm_h%k(m) 1288 tend(k,j,i) = tend(k,j,i) + g / pt_reference * & 1289 drho_air_zw(k-1) * & 1290 surf_lsm_h%shf(m) 1291 ENDDO 1292 ! 1293 !-- Urban surfaces 1294 surf_s = surf_usm_h%start_index(j,i) 1295 surf_e = surf_usm_h%end_index(j,i) 1296 DO m = surf_s, surf_e 1297 k = surf_usm_h%k(m) 1298 tend(k,j,i) = tend(k,j,i) + g / pt_reference * & 1299 drho_air_zw(k-1) * & 1300 surf_usm_h%shf(m) 1301 ENDDO 1302 ENDIF 1303 1304 IF ( use_top_fluxes ) THEN 1305 surf_s = surf_def_h(2)%start_index(j,i) 1306 surf_e = surf_def_h(2)%end_index(j,i) 1307 DO m = surf_s, surf_e 1308 k = surf_def_h(2)%k(m) 1309 tend(k,j,i) = tend(k,j,i) + g / pt_reference * & 1310 drho_air_zw(k-1) * & 1311 surf_def_h(2)%shf(m) 1312 ENDDO 1313 ENDIF 1314 1265 ENDDO 1266 ! 1267 !-- Natural surfaces 1268 surf_s = surf_lsm_h%start_index(j,i) 1269 surf_e = surf_lsm_h%end_index(j,i) 1270 DO m = surf_s, surf_e 1271 k = surf_lsm_h%k(m) 1272 tend(k,j,i) = tend(k,j,i) + g / & 1273 MERGE( pt_reference, pt(k,j,i), & 1274 use_single_reference_value ) * & 1275 drho_air_zw(k-1) * & 1276 surf_lsm_h%shf(m) 1277 ENDDO 1278 ! 1279 !-- Urban surfaces 1280 surf_s = surf_usm_h%start_index(j,i) 1281 surf_e = surf_usm_h%end_index(j,i) 1282 DO m = surf_s, surf_e 1283 k = surf_usm_h%k(m) 1284 tend(k,j,i) = tend(k,j,i) + g / & 1285 MERGE( pt_reference, pt(k,j,i), & 1286 use_single_reference_value ) * & 1287 drho_air_zw(k-1) * & 1288 surf_usm_h%shf(m) 1289 ENDDO 1315 1290 ENDIF 1316 1291 1317 ELSE 1318 1319 IF ( ocean ) THEN 1320 ! 1321 !-- So far in the ocean no special treatment of density flux in 1322 !-- the bottom and top surface layer 1323 DO k = nzb+1, nzt 1324 tend(k,j,i) = tend(k,j,i) + & 1325 kh(k,j,i) * g / prho(k,j,i) * & 1326 ( prho(k+1,j,i) - prho(k-1,j,i) ) * & 1327 dd2zu(k) * & 1328 MERGE( 1.0_wp, 0.0_wp, & 1329 BTEST( wall_flags_0(k,j,i), 0 ) & 1330 ) 1292 IF ( use_top_fluxes ) THEN 1293 surf_s = surf_def_h(2)%start_index(j,i) 1294 surf_e = surf_def_h(2)%end_index(j,i) 1295 DO m = surf_s, surf_e 1296 k = surf_def_h(2)%k(m) 1297 tend(k,j,i) = tend(k,j,i) + g / & 1298 MERGE( pt_reference, pt(k,j,i), & 1299 use_single_reference_value ) * & 1300 drho_air_zw(k) * & 1301 surf_def_h(2)%shf(m) 1331 1302 ENDDO 1332 1333 ELSE1334 1335 DO k = nzb+1, nzt1336 !1337 !-- Flag 9 is used to mask top fluxes, flag 30 to mask1338 !-- surface fluxes1339 tend(k,j,i) = tend(k,j,i) - &1340 kh(k,j,i) * g / pt(k,j,i) * &1341 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) * &1342 MERGE( 1.0_wp, 0.0_wp, &1343 BTEST( wall_flags_0(k,j,i), 30 ) &1344 ) * &1345 MERGE( 1.0_wp, 0.0_wp, &1346 BTEST( wall_flags_0(k,j,i), 9 ) &1347 )1348 ENDDO1349 1350 IF ( use_surface_fluxes ) THEN1351 !1352 !-- Default surfaces, up- and downward-facing1353 DO l = 0, 11354 surf_s = surf_def_h(l)%start_index(j,i)1355 surf_e = surf_def_h(l)%end_index(j,i)1356 DO m = surf_s, surf_e1357 k = surf_def_h(l)%k(m)1358 tend(k,j,i) = tend(k,j,i) + g / pt_reference &1359 * drho_air_zw(k-1) &1360 * surf_def_h(l)%shf(m)1361 ENDDO1362 ENDDO1363 !1364 !-- Natural surfaces1365 surf_s = surf_lsm_h%start_index(j,i)1366 surf_e = surf_lsm_h%end_index(j,i)1367 DO m = surf_s, surf_e1368 k = surf_lsm_h%k(m)1369 tend(k,j,i) = tend(k,j,i) + g / pt_reference &1370 * drho_air_zw(k-1) &1371 * surf_lsm_h%shf(m)1372 ENDDO1373 !1374 !-- Urban surfaces1375 surf_s = surf_usm_h%start_index(j,i)1376 surf_e = surf_usm_h%end_index(j,i)1377 DO m = surf_s, surf_e1378 k = surf_usm_h%k(m)1379 tend(k,j,i) = tend(k,j,i) + g / pt_reference &1380 * drho_air_zw(k-1) &1381 * surf_usm_h%shf(m)1382 ENDDO1383 ENDIF1384 1385 IF ( use_top_fluxes ) THEN1386 surf_s = surf_def_h(2)%start_index(j,i)1387 surf_e = surf_def_h(2)%end_index(j,i)1388 DO m = surf_s, surf_e1389 k = surf_def_h(2)%k(m)1390 tend(k,j,i) = tend(k,j,i) + g / pt_reference * &1391 drho_air_zw(k-1) * &1392 surf_def_h(2)%shf(m)1393 ENDDO1394 ENDIF1395 1396 1303 ENDIF 1397 1304 … … 1406 1313 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1407 1314 k2 = 0.61_wp * pt(k,j,i) 1408 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 1315 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / & 1316 MERGE( vpt_reference, vpt(k,j,i), & 1317 use_single_reference_value ) * & 1409 1318 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 1410 1319 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & … … 1431 1340 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1432 1341 ENDIF 1433 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 1342 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / & 1343 MERGE( vpt_reference, vpt(k,j,i), & 1344 use_single_reference_value ) * & 1434 1345 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 1435 1346 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & … … 1444 1355 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1445 1356 k2 = 0.61_wp * pt(k,j,i) 1446 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 1357 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / & 1358 MERGE( vpt_reference, vpt(k,j,i), & 1359 use_single_reference_value ) * & 1447 1360 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 1448 1361 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - & … … 1489 1402 ENDIF 1490 1403 1491 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 1404 tend(k,j,i) = tend(k,j,i) + g / & 1405 MERGE( vpt_reference, vpt(k,j,i), & 1406 use_single_reference_value ) * & 1492 1407 ( k1 * surf_def_h(l)%shf(m) + & 1493 1408 k2 * surf_def_h(l)%qsws(m) & … … 1524 1439 ENDIF 1525 1440 1526 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 1441 tend(k,j,i) = tend(k,j,i) + g / & 1442 MERGE( vpt_reference, vpt(k,j,i), & 1443 use_single_reference_value ) * & 1527 1444 ( k1 * surf_lsm_h%shf(m) + & 1528 1445 k2 * surf_lsm_h%qsws(m) & … … 1558 1475 ENDIF 1559 1476 1560 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 1477 tend(k,j,i) = tend(k,j,i) + g / & 1478 MERGE( vpt_reference, vpt(k,j,i), & 1479 use_single_reference_value ) * & 1561 1480 ( k1 * surf_usm_h%shf(m) + & 1562 1481 k2 * surf_usm_h%qsws(m) & … … 1596 1515 ENDIF 1597 1516 1598 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 1517 tend(k,j,i) = tend(k,j,i) + g / & 1518 MERGE( vpt_reference, vpt(k,j,i), & 1519 use_single_reference_value ) * & 1599 1520 ( k1* surf_def_h(2)%shf(m) + & 1600 1521 k2 * surf_def_h(2)%qsws(m) & 1601 ) * drho_air_zw(k -1)1522 ) * drho_air_zw(k) 1602 1523 ENDDO 1603 1524 … … 1682 1603 ENDDO 1683 1604 ! 1684 !-- Default surfaces, downward-facing 1605 !-- Default surfaces, downward-facing surfaces 1685 1606 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 1686 1607 DO m = 1, surf_def_h(1)%ns … … 1689 1610 j = surf_def_h(1)%j(m) 1690 1611 k = surf_def_h(1)%k(m) 1691 ! 1692 !-- Note, calculatione of u_0 and v_0 is not fully accurate, as u/v 1693 !-- and km are not on the same grid. Actually, a further 1694 !-- interpolation of km onto the u/v-grid is necessary. However, the 1695 !-- effect of this error is negligible. 1696 !-- In case of downward-facing surfaces, gradient is calculated 1697 !-- between u_0 and u(k-1). 1698 surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) * & 1699 drho_air_zw(k-1) * & 1700 ( zu(k+1) - zu(k-1) ) / & 1612 1613 surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) * & 1614 drho_air_zw(k-1) * & 1615 ( zu(k+1) - zu(k-1) ) / & 1701 1616 ( km(k,j,i) + 1.0E-20_wp ) 1702 surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) * 1703 drho_air_zw(k-1) * 1704 ( zu(k+1) - zu(k-1) ) / 1617 surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) * & 1618 drho_air_zw(k-1) * & 1619 ( zu(k+1) - zu(k-1) ) / & 1705 1620 ( km(k,j,i) + 1.0E-20_wp ) 1706 1621 1707 IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) > 1708 ABS( u(k+1,j,i) - u(k-1,j,i) ) 1622 IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) > & 1623 ABS( u(k+1,j,i) - u(k-1,j,i) ) & 1709 1624 ) surf_def_h(1)%u_0(m) = u(k+1,j,i) 1710 1625 1711 IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) > 1712 ABS( v(k+1,j,i) - v(k-1,j,i) ) 1626 IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) > & 1627 ABS( v(k+1,j,i) - v(k-1,j,i) ) & 1713 1628 ) surf_def_h(1)%v_0(m) = v(k+1,j,i) 1714 1629
Note: See TracChangeset
for help on using the changeset viewer.
