Changeset 940 for palm/trunk/SOURCE/production_e.f90
- Timestamp:
- Jul 9, 2012 2:31:00 PM (12 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/production_e.f90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/production_e.f90
r760 r940 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! TKE production by buoyancy can be switched off in case of runs with pure 7 ! neutral stratification 7 8 ! 8 9 ! Former revisions: … … 127 128 ! ENDIF 128 129 129 ! 130 !-- Calculate TKE production by shear 130 131 131 DO i = nxl, nxr 132 132 133 ! 134 !-- Calculate TKE production by shear 133 135 DO j = nys, nyn 134 136 DO k = nzb_diff_s_outer(j,i), nzt … … 453 455 454 456 ! 455 !-- Calculate TKE production by buoyancy 456 IF ( .NOT. humidity ) THEN 457 458 IF ( use_reference ) THEN 459 460 IF ( ocean ) THEN 461 ! 462 !-- So far in the ocean no special treatment of density flux in 463 !-- the bottom and top surface layer 464 DO j = nys, nyn 465 DO k = nzb_s_inner(j,i)+1, nzt 466 tend(k,j,i) = tend(k,j,i) + & 467 kh(k,j,i) * g / rho_reference * & 468 ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k) 457 !-- If required, calculate TKE production by buoyancy 458 IF ( .NOT. neutral ) THEN 459 460 IF ( .NOT. humidity ) THEN 461 462 IF ( use_reference ) THEN 463 464 IF ( ocean ) THEN 465 ! 466 !-- So far in the ocean no special treatment of density flux 467 !-- in the bottom and top surface layer 468 DO j = nys, nyn 469 DO k = nzb_s_inner(j,i)+1, nzt 470 tend(k,j,i) = tend(k,j,i) + & 471 kh(k,j,i) * g / rho_reference * & 472 ( rho(k+1,j,i) - rho(k-1,j,i) ) * & 473 dd2zu(k) 474 ENDDO 469 475 ENDDO 470 ENDDO 476 477 ELSE 478 479 DO j = nys, nyn 480 DO k = nzb_diff_s_inner(j,i), nzt_diff 481 tend(k,j,i) = tend(k,j,i) - & 482 kh(k,j,i) * g / pt_reference * & 483 ( pt(k+1,j,i) - pt(k-1,j,i) ) * & 484 dd2zu(k) 485 ENDDO 486 487 IF ( use_surface_fluxes ) THEN 488 k = nzb_diff_s_inner(j,i)-1 489 tend(k,j,i) = tend(k,j,i) + g / pt_reference * & 490 shf(j,i) 491 ENDIF 492 493 IF ( use_top_fluxes ) THEN 494 k = nzt 495 tend(k,j,i) = tend(k,j,i) + g / pt_reference * & 496 tswst(j,i) 497 ENDIF 498 ENDDO 499 500 ENDIF 471 501 472 502 ELSE 473 503 474 DO j = nys, nyn 475 DO k = nzb_diff_s_inner(j,i), nzt_diff 476 tend(k,j,i) = tend(k,j,i) - & 477 kh(k,j,i) * g / pt_reference * & 478 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) 504 IF ( ocean ) THEN 505 ! 506 !-- So far in the ocean no special treatment of density flux 507 !-- in the bottom and top surface layer 508 DO j = nys, nyn 509 DO k = nzb_s_inner(j,i)+1, nzt 510 tend(k,j,i) = tend(k,j,i) + & 511 kh(k,j,i) * g / rho(k,j,i) * & 512 ( rho(k+1,j,i) - rho(k-1,j,i) ) * & 513 dd2zu(k) 514 ENDDO 479 515 ENDDO 480 516 481 IF ( use_surface_fluxes ) THEN 482 k = nzb_diff_s_inner(j,i)-1 483 tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i) 484 ENDIF 485 486 IF ( use_top_fluxes ) THEN 487 k = nzt 488 tend(k,j,i) = tend(k,j,i) + g / pt_reference * & 489 tswst(j,i) 490 ENDIF 491 ENDDO 517 ELSE 518 519 DO j = nys, nyn 520 DO k = nzb_diff_s_inner(j,i), nzt_diff 521 tend(k,j,i) = tend(k,j,i) - & 522 kh(k,j,i) * g / pt(k,j,i) * & 523 ( pt(k+1,j,i) - pt(k-1,j,i) ) * & 524 dd2zu(k) 525 ENDDO 526 527 IF ( use_surface_fluxes ) THEN 528 k = nzb_diff_s_inner(j,i)-1 529 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * & 530 shf(j,i) 531 ENDIF 532 533 IF ( use_top_fluxes ) THEN 534 k = nzt 535 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * & 536 tswst(j,i) 537 ENDIF 538 ENDDO 539 540 ENDIF 492 541 493 542 ENDIF … … 495 544 ELSE 496 545 497 IF ( ocean ) THEN 498 ! 499 !-- So far in the ocean no special treatment of density flux in 500 !-- the bottom and top surface layer 501 DO j = nys, nyn 502 DO k = nzb_s_inner(j,i)+1, nzt 503 tend(k,j,i) = tend(k,j,i) + & 504 kh(k,j,i) * g / rho(k,j,i) * & 505 ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k) 506 ENDDO 507 ENDDO 508 509 ELSE 510 511 DO j = nys, nyn 512 DO k = nzb_diff_s_inner(j,i), nzt_diff 513 tend(k,j,i) = tend(k,j,i) - & 514 kh(k,j,i) * g / pt(k,j,i) * & 515 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) 516 ENDDO 517 518 IF ( use_surface_fluxes ) THEN 519 k = nzb_diff_s_inner(j,i)-1 520 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i) 521 ENDIF 522 523 IF ( use_top_fluxes ) THEN 524 k = nzt 525 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i) 526 ENDIF 527 ENDDO 528 529 ENDIF 530 531 ENDIF 532 533 ELSE 534 535 DO j = nys, nyn 536 537 DO k = nzb_diff_s_inner(j,i), nzt_diff 538 539 IF ( .NOT. cloud_physics ) THEN 540 k1 = 1.0 + 0.61 * q(k,j,i) 541 k2 = 0.61 * pt(k,j,i) 542 ELSE 543 IF ( ql(k,j,i) == 0.0 ) THEN 546 DO j = nys, nyn 547 548 DO k = nzb_diff_s_inner(j,i), nzt_diff 549 550 IF ( .NOT. cloud_physics ) THEN 544 551 k1 = 1.0 + 0.61 * q(k,j,i) 545 552 k2 = 0.61 * pt(k,j,i) 546 553 ELSE 547 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 548 temp = theta * t_d_pt(k) 549 k1 = ( 1.0 - q(k,j,i) + 1.61 * & 550 ( q(k,j,i) - ql(k,j,i) ) * & 551 ( 1.0 + 0.622 * l_d_r / temp ) ) / & 552 ( 1.0 + 0.622 * l_d_r * l_d_cp * & 553 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 554 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 554 IF ( ql(k,j,i) == 0.0 ) THEN 555 k1 = 1.0 + 0.61 * q(k,j,i) 556 k2 = 0.61 * pt(k,j,i) 557 ELSE 558 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 559 temp = theta * t_d_pt(k) 560 k1 = ( 1.0 - q(k,j,i) + 1.61 * & 561 ( q(k,j,i) - ql(k,j,i) ) * & 562 ( 1.0 + 0.622 * l_d_r / temp ) ) / & 563 ( 1.0 + 0.622 * l_d_r * l_d_cp * & 564 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 565 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 566 ENDIF 555 567 ENDIF 556 ENDIF 557 558 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 559 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 560 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & 561 ) * dd2zu(k) 568 569 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 570 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 571 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & 572 ) * dd2zu(k) 573 ENDDO 574 562 575 ENDDO 563 576 564 ENDDO 565 566 IF ( use_surface_fluxes ) THEN 567 568 DO j = nys, nyn 569 570 k = nzb_diff_s_inner(j,i)-1 571 572 IF ( .NOT. cloud_physics ) THEN 573 k1 = 1.0 + 0.61 * q(k,j,i) 574 k2 = 0.61 * pt(k,j,i) 575 ELSE 576 IF ( ql(k,j,i) == 0.0 ) THEN 577 IF ( use_surface_fluxes ) THEN 578 579 DO j = nys, nyn 580 581 k = nzb_diff_s_inner(j,i)-1 582 583 IF ( .NOT. cloud_physics ) THEN 577 584 k1 = 1.0 + 0.61 * q(k,j,i) 578 585 k2 = 0.61 * pt(k,j,i) 579 586 ELSE 580 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 581 temp = theta * t_d_pt(k) 582 k1 = ( 1.0 - q(k,j,i) + 1.61 * & 583 ( q(k,j,i) - ql(k,j,i) ) * & 584 ( 1.0 + 0.622 * l_d_r / temp ) ) / & 585 ( 1.0 + 0.622 * l_d_r * l_d_cp * & 586 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 587 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 587 IF ( ql(k,j,i) == 0.0 ) THEN 588 k1 = 1.0 + 0.61 * q(k,j,i) 589 k2 = 0.61 * pt(k,j,i) 590 ELSE 591 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 592 temp = theta * t_d_pt(k) 593 k1 = ( 1.0 - q(k,j,i) + 1.61 * & 594 ( q(k,j,i) - ql(k,j,i) ) * & 595 ( 1.0 + 0.622 * l_d_r / temp ) ) / & 596 ( 1.0 + 0.622 * l_d_r * l_d_cp * & 597 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 598 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 599 ENDIF 588 600 ENDIF 589 ENDIF 590 591 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 592 ( k1* shf(j,i) + k2 * qsws(j,i) ) 593 ENDDO 594 595 ENDIF 596 597 IF ( use_top_fluxes ) THEN 598 599 DO j = nys, nyn 600 601 k = nzt 602 603 IF ( .NOT. cloud_physics ) THEN 604 k1 = 1.0 + 0.61 * q(k,j,i) 605 k2 = 0.61 * pt(k,j,i) 606 ELSE 607 IF ( ql(k,j,i) == 0.0 ) THEN 601 602 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 603 ( k1* shf(j,i) + k2 * qsws(j,i) ) 604 ENDDO 605 606 ENDIF 607 608 IF ( use_top_fluxes ) THEN 609 610 DO j = nys, nyn 611 612 k = nzt 613 614 IF ( .NOT. cloud_physics ) THEN 608 615 k1 = 1.0 + 0.61 * q(k,j,i) 609 616 k2 = 0.61 * pt(k,j,i) 610 617 ELSE 611 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 612 temp = theta * t_d_pt(k) 613 k1 = ( 1.0 - q(k,j,i) + 1.61 * & 614 ( q(k,j,i) - ql(k,j,i) ) * & 615 ( 1.0 + 0.622 * l_d_r / temp ) ) / & 616 ( 1.0 + 0.622 * l_d_r * l_d_cp * & 617 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 618 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 618 IF ( ql(k,j,i) == 0.0 ) THEN 619 k1 = 1.0 + 0.61 * q(k,j,i) 620 k2 = 0.61 * pt(k,j,i) 621 ELSE 622 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 623 temp = theta * t_d_pt(k) 624 k1 = ( 1.0 - q(k,j,i) + 1.61 * & 625 ( q(k,j,i) - ql(k,j,i) ) * & 626 ( 1.0 + 0.622 * l_d_r / temp ) ) / & 627 ( 1.0 + 0.622 * l_d_r * l_d_cp * & 628 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 629 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 630 ENDIF 619 631 ENDIF 620 ENDIF 621 622 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 623 ( k1* tswst(j,i) + k2 * qswst(j,i) ) 624 ENDDO 632 633 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 634 ( k1* tswst(j,i) + k2 * qswst(j,i) ) 635 ENDDO 636 637 ENDIF 625 638 626 639 ENDIF … … 950 963 951 964 ! 952 !-- Calculate TKE production by buoyancy 953 IF ( .NOT. humidity ) THEN 954 955 IF ( use_reference ) THEN 956 957 IF ( ocean ) THEN 958 ! 959 !-- So far in the ocean no special treatment of density flux in the 960 !-- bottom and top surface layer 961 DO k = nzb_s_inner(j,i)+1, nzt 962 tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / rho_reference * & 963 ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k) 964 ENDDO 965 !-- If required, calculate TKE production by buoyancy 966 IF ( .NOT. neutral ) THEN 967 968 IF ( .NOT. humidity ) THEN 969 970 IF ( use_reference ) THEN 971 972 IF ( ocean ) THEN 973 ! 974 !-- So far in the ocean no special treatment of density flux in 975 !-- the bottom and top surface layer 976 DO k = nzb_s_inner(j,i)+1, nzt 977 tend(k,j,i) = tend(k,j,i) + & 978 kh(k,j,i) * g / rho_reference * & 979 ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k) 980 ENDDO 981 982 ELSE 983 984 DO k = nzb_diff_s_inner(j,i), nzt_diff 985 tend(k,j,i) = tend(k,j,i) - & 986 kh(k,j,i) * g / pt_reference * & 987 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) 988 ENDDO 989 990 IF ( use_surface_fluxes ) THEN 991 k = nzb_diff_s_inner(j,i)-1 992 tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i) 993 ENDIF 994 995 IF ( use_top_fluxes ) THEN 996 k = nzt 997 tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i) 998 ENDIF 999 1000 ENDIF 965 1001 966 1002 ELSE 967 1003 968 DO k = nzb_diff_s_inner(j,i), nzt_diff 969 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt_reference * & 970 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) 971 ENDDO 972 973 IF ( use_surface_fluxes ) THEN 974 k = nzb_diff_s_inner(j,i)-1 975 tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i) 976 ENDIF 977 978 IF ( use_top_fluxes ) THEN 979 k = nzt 980 tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i) 1004 IF ( ocean ) THEN 1005 ! 1006 !-- So far in the ocean no special treatment of density flux in 1007 !-- the bottom and top surface layer 1008 DO k = nzb_s_inner(j,i)+1, nzt 1009 tend(k,j,i) = tend(k,j,i) + & 1010 kh(k,j,i) * g / rho(k,j,i) * & 1011 ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k) 1012 ENDDO 1013 1014 ELSE 1015 1016 DO k = nzb_diff_s_inner(j,i), nzt_diff 1017 tend(k,j,i) = tend(k,j,i) - & 1018 kh(k,j,i) * g / pt(k,j,i) * & 1019 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) 1020 ENDDO 1021 1022 IF ( use_surface_fluxes ) THEN 1023 k = nzb_diff_s_inner(j,i)-1 1024 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i) 1025 ENDIF 1026 1027 IF ( use_top_fluxes ) THEN 1028 k = nzt 1029 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i) 1030 ENDIF 1031 981 1032 ENDIF 982 1033 … … 985 1036 ELSE 986 1037 987 IF ( ocean ) THEN 988 ! 989 !-- So far in the ocean no special treatment of density flux in the 990 !-- bottom and top surface layer 991 DO k = nzb_s_inner(j,i)+1, nzt 992 tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / rho(k,j,i) * & 993 ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k) 994 ENDDO 995 996 ELSE 997 998 DO k = nzb_diff_s_inner(j,i), nzt_diff 999 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * & 1000 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) 1001 ENDDO 1002 1003 IF ( use_surface_fluxes ) THEN 1004 k = nzb_diff_s_inner(j,i)-1 1005 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i) 1006 ENDIF 1007 1008 IF ( use_top_fluxes ) THEN 1009 k = nzt 1010 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i) 1011 ENDIF 1012 1013 ENDIF 1014 1015 ENDIF 1016 1017 ELSE 1018 1019 DO k = nzb_diff_s_inner(j,i), nzt_diff 1020 1021 IF ( .NOT. cloud_physics ) THEN 1022 k1 = 1.0 + 0.61 * q(k,j,i) 1023 k2 = 0.61 * pt(k,j,i) 1024 ELSE 1025 IF ( ql(k,j,i) == 0.0 ) THEN 1038 DO k = nzb_diff_s_inner(j,i), nzt_diff 1039 1040 IF ( .NOT. cloud_physics ) THEN 1026 1041 k1 = 1.0 + 0.61 * q(k,j,i) 1027 1042 k2 = 0.61 * pt(k,j,i) 1028 1043 ELSE 1029 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)1030 temp = theta * t_d_pt(k)1031 k1 = ( 1.0 - q(k,j,i) + 1.61 * &1032 ( q(k,j,i) - ql(k,j,i) ) * &1033 ( 1.0 + 0.622 * l_d_r / temp ) ) / &1034 ( 1.0 + 0.622 * l_d_r * l_d_cp * &1035 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )1036 k2 = theta * ( l_d_cp / temp * k1 - 1.0 )1037 ENDIF1038 ENDIF1039 1040 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &1041 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &1042 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &1043 ) * dd2zu(k) 1044 ENDDO1045 1046 IF ( use_surface_fluxes ) THEN1047 k = nzb_diff_s_inner(j,i)-11048 1049 IF ( .NOT. cloud_physics ) THEN 1050 k1 = 1.0 + 0.61 * q(k,j,i)1051 k 2 = 0.61 * pt(k,j,i)1052 ELSE 1053 IF ( ql(k,j,i) == 0.0 )THEN1044 IF ( ql(k,j,i) == 0.0 ) THEN 1045 k1 = 1.0 + 0.61 * q(k,j,i) 1046 k2 = 0.61 * pt(k,j,i) 1047 ELSE 1048 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1049 temp = theta * t_d_pt(k) 1050 k1 = ( 1.0 - q(k,j,i) + 1.61 * & 1051 ( q(k,j,i) - ql(k,j,i) ) * & 1052 ( 1.0 + 0.622 * l_d_r / temp ) ) / & 1053 ( 1.0 + 0.622 * l_d_r * l_d_cp * & 1054 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1055 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 1056 ENDIF 1057 ENDIF 1058 1059 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 1060 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 1061 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & 1062 ) * dd2zu(k) 1063 ENDDO 1064 1065 IF ( use_surface_fluxes ) THEN 1066 k = nzb_diff_s_inner(j,i)-1 1067 1068 IF ( .NOT. cloud_physics ) THEN 1054 1069 k1 = 1.0 + 0.61 * q(k,j,i) 1055 1070 k2 = 0.61 * pt(k,j,i) 1056 1071 ELSE 1057 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1058 temp = theta * t_d_pt(k) 1059 k1 = ( 1.0 - q(k,j,i) + 1.61 * & 1060 ( q(k,j,i) - ql(k,j,i) ) * & 1061 ( 1.0 + 0.622 * l_d_r / temp ) ) / & 1062 ( 1.0 + 0.622 * l_d_r * l_d_cp * & 1063 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1064 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 1065 ENDIF 1072 IF ( ql(k,j,i) == 0.0 ) THEN 1073 k1 = 1.0 + 0.61 * q(k,j,i) 1074 k2 = 0.61 * pt(k,j,i) 1075 ELSE 1076 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1077 temp = theta * t_d_pt(k) 1078 k1 = ( 1.0 - q(k,j,i) + 1.61 * & 1079 ( q(k,j,i) - ql(k,j,i) ) * & 1080 ( 1.0 + 0.622 * l_d_r / temp ) ) / & 1081 ( 1.0 + 0.622 * l_d_r * l_d_cp * & 1082 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1083 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 1084 ENDIF 1085 ENDIF 1086 1087 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 1088 ( k1* shf(j,i) + k2 * qsws(j,i) ) 1066 1089 ENDIF 1067 1090 1068 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 1069 ( k1* shf(j,i) + k2 * qsws(j,i) ) 1070 ENDIF 1071 1072 IF ( use_top_fluxes ) THEN 1073 k = nzt 1074 1075 IF ( .NOT. cloud_physics ) THEN 1076 k1 = 1.0 + 0.61 * q(k,j,i) 1077 k2 = 0.61 * pt(k,j,i) 1078 ELSE 1079 IF ( ql(k,j,i) == 0.0 ) THEN 1091 IF ( use_top_fluxes ) THEN 1092 k = nzt 1093 1094 IF ( .NOT. cloud_physics ) THEN 1080 1095 k1 = 1.0 + 0.61 * q(k,j,i) 1081 1096 k2 = 0.61 * pt(k,j,i) 1082 1097 ELSE 1083 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1084 temp = theta * t_d_pt(k) 1085 k1 = ( 1.0 - q(k,j,i) + 1.61 * & 1086 ( q(k,j,i) - ql(k,j,i) ) * & 1087 ( 1.0 + 0.622 * l_d_r / temp ) ) / & 1088 ( 1.0 + 0.622 * l_d_r * l_d_cp * & 1089 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1090 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 1091 ENDIF 1098 IF ( ql(k,j,i) == 0.0 ) THEN 1099 k1 = 1.0 + 0.61 * q(k,j,i) 1100 k2 = 0.61 * pt(k,j,i) 1101 ELSE 1102 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1103 temp = theta * t_d_pt(k) 1104 k1 = ( 1.0 - q(k,j,i) + 1.61 * & 1105 ( q(k,j,i) - ql(k,j,i) ) * & 1106 ( 1.0 + 0.622 * l_d_r / temp ) ) / & 1107 ( 1.0 + 0.622 * l_d_r * l_d_cp * & 1108 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1109 k2 = theta * ( l_d_cp / temp * k1 - 1.0 ) 1110 ENDIF 1111 ENDIF 1112 1113 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 1114 ( k1* tswst(j,i) + k2 * qswst(j,i) ) 1092 1115 ENDIF 1093 1116 1094 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &1095 ( k1* tswst(j,i) + k2 * qswst(j,i) )1096 1117 ENDIF 1097 1118
Note: See TracChangeset
for help on using the changeset viewer.
