Changeset 940 for palm/trunk/SOURCE
- Timestamp:
- Jul 9, 2012 2:31:00 PM (12 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 9 edited
-
check_parameters.f90 (modified) (3 diffs)
-
header.f90 (modified) (3 diffs)
-
modules.f90 (modified) (2 diffs)
-
parin.f90 (modified) (2 diffs)
-
poisfft.f90 (modified) (9 diffs)
-
production_e.f90 (modified) (6 diffs)
-
prognostic_equations.f90 (modified) (16 diffs)
-
read_var_list.f90 (modified) (2 diffs)
-
write_var_list.f90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/check_parameters.f90
r925 r940 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! checks for parameter neutral 7 7 ! 8 8 ! Former revisions: … … 1039 1039 ! 1040 1040 !-- Compute initial temperature profile using the given temperature gradients 1041 i = 1 1042 gradient = 0.0 1043 1044 IF ( .NOT. ocean ) THEN 1045 1046 pt_vertical_gradient_level_ind(1) = 0 1047 DO k = 1, nzt+1 1048 IF ( i < 11 ) THEN 1049 IF ( pt_vertical_gradient_level(i) < zu(k) .AND. & 1050 pt_vertical_gradient_level(i) >= 0.0 ) THEN 1051 gradient = pt_vertical_gradient(i) / 100.0 1052 pt_vertical_gradient_level_ind(i) = k - 1 1053 i = i + 1 1041 IF ( .NOT. neutral ) THEN 1042 1043 i = 1 1044 gradient = 0.0 1045 1046 IF ( .NOT. ocean ) THEN 1047 1048 pt_vertical_gradient_level_ind(1) = 0 1049 DO k = 1, nzt+1 1050 IF ( i < 11 ) THEN 1051 IF ( pt_vertical_gradient_level(i) < zu(k) .AND. & 1052 pt_vertical_gradient_level(i) >= 0.0 ) THEN 1053 gradient = pt_vertical_gradient(i) / 100.0 1054 pt_vertical_gradient_level_ind(i) = k - 1 1055 i = i + 1 1056 ENDIF 1054 1057 ENDIF 1055 ENDIF 1056 IF ( gradient /= 0.0 ) THEN 1057 IF ( k /= 1 ) THEN 1058 pt_init(k) = pt_init(k-1) + dzu(k) * gradient 1058 IF ( gradient /= 0.0 ) THEN 1059 IF ( k /= 1 ) THEN 1060 pt_init(k) = pt_init(k-1) + dzu(k) * gradient 1061 ELSE 1062 pt_init(k) = pt_surface + 0.5 * dzu(k) * gradient 1063 ENDIF 1059 1064 ELSE 1060 pt_init(k) = pt_ surface + 0.5 * dzu(k) * gradient1065 pt_init(k) = pt_init(k-1) 1061 1066 ENDIF 1062 ELSE 1063 pt_init(k) = pt_init(k-1) 1064 ENDIF 1065 ENDDO 1066 1067 ELSE 1068 1069 pt_vertical_gradient_level_ind(1) = nzt+1 1070 DO k = nzt, 0, -1 1071 IF ( i < 11 ) THEN 1072 IF ( pt_vertical_gradient_level(i) > zu(k) .AND. & 1073 pt_vertical_gradient_level(i) <= 0.0 ) THEN 1074 gradient = pt_vertical_gradient(i) / 100.0 1075 pt_vertical_gradient_level_ind(i) = k + 1 1076 i = i + 1 1067 ENDDO 1068 1069 ELSE 1070 1071 pt_vertical_gradient_level_ind(1) = nzt+1 1072 DO k = nzt, 0, -1 1073 IF ( i < 11 ) THEN 1074 IF ( pt_vertical_gradient_level(i) > zu(k) .AND. & 1075 pt_vertical_gradient_level(i) <= 0.0 ) THEN 1076 gradient = pt_vertical_gradient(i) / 100.0 1077 pt_vertical_gradient_level_ind(i) = k + 1 1078 i = i + 1 1079 ENDIF 1077 1080 ENDIF 1078 ENDIF 1079 IF ( gradient /= 0.0 ) THEN 1080 IF ( k /= nzt ) THEN 1081 pt_init(k) = pt_init(k+1) - dzu(k+1) * gradient 1081 IF ( gradient /= 0.0 ) THEN 1082 IF ( k /= nzt ) THEN 1083 pt_init(k) = pt_init(k+1) - dzu(k+1) * gradient 1084 ELSE 1085 pt_init(k) = pt_surface - 0.5 * dzu(k+1) * gradient 1086 pt_init(k+1) = pt_surface + 0.5 * dzu(k+1) * gradient 1087 ENDIF 1082 1088 ELSE 1083 pt_init(k) = pt_surface - 0.5 * dzu(k+1) * gradient 1084 pt_init(k+1) = pt_surface + 0.5 * dzu(k+1) * gradient 1089 pt_init(k) = pt_init(k+1) 1085 1090 ENDIF 1086 ELSE 1087 pt_init(k) = pt_init(k+1) 1088 ENDIF 1089 ENDDO 1091 ENDDO 1092 1093 ENDIF 1090 1094 1091 1095 ENDIF … … 1505 1509 IF ( surface_heatflux == 9999999.9 ) constant_heatflux = .FALSE. 1506 1510 IF ( top_heatflux == 9999999.9 ) constant_top_heatflux = .FALSE. 1511 1512 IF ( neutral ) THEN 1513 1514 IF ( surface_heatflux /= 0.0 .AND. surface_heatflux /= 9999999.9 ) & 1515 THEN 1516 message_string = 'heatflux must not be set for pure neutral flow' 1517 CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 ) 1518 ENDIF 1519 1520 IF ( top_heatflux /= 0.0 .AND. top_heatflux /= 9999999.9 ) & 1521 THEN 1522 message_string = 'heatflux must not be set for pure neutral flow' 1523 CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 ) 1524 ENDIF 1525 1526 ENDIF 1527 1507 1528 IF ( top_momentumflux_u /= 9999999.9 .AND. & 1508 1529 top_momentumflux_v /= 9999999.9 ) THEN -
palm/trunk/SOURCE/header.f90
r928 r940 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! Output in case of simulations for pure neutral stratification (no pt-equation 7 ! solved) 7 8 ! 8 9 ! Former revisions: … … 349 350 ENDIF 350 351 ENDIF 352 IF ( neutral ) WRITE ( io, 131 ) pt_surface 351 353 IF ( humidity ) THEN 352 354 IF ( .NOT. cloud_physics ) THEN … … 1682 1684 129 FORMAT (' --> Additional prognostic equation for the specific humidity') 1683 1685 130 FORMAT (' --> Additional prognostic equation for the total water content') 1686 131 FORMAT (' --> No pt-equation solved. Neutral stratification with pt = ', & 1687 F6.2, ' K assumed') 1684 1688 132 FORMAT (' Parameterization of long-wave radiation processes via'/ & 1685 1689 ' effective emissivity scheme') -
palm/trunk/SOURCE/modules.f90
r937 r940 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! +neutral 7 7 ! 8 8 ! Former revisions: … … 591 591 iso2d_output = .FALSE., large_scale_subsidence = .FALSE., & 592 592 masking_method = .FALSE., mg_switch_to_pe0 = .FALSE., & 593 netcdf_output = .FALSE., ocean = .FALSE., &593 netcdf_output = .FALSE., neutral = .FALSE., ocean = .FALSE., & 594 594 outflow_l = .FALSE., outflow_n = .FALSE., outflow_r = .FALSE., & 595 595 outflow_s = .FALSE., passive_scalar = .FALSE., & -
palm/trunk/SOURCE/parin.f90
r937 r940 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! +neutral in inipar 7 7 ! 8 8 ! Former revisions: … … 165 165 leaf_surface_concentration, long_filter_factor, & 166 166 loop_optimization, masking_method, mg_cycles, & 167 mg_switch_to_pe0_level, &168 mixing_length_1d, momentum_advec, netcdf_precision, ngsrb, nsor, &167 mg_switch_to_pe0_level, mixing_length_1d, momentum_advec, & 168 netcdf_precision, neutral, ngsrb, nsor, & 169 169 nsor_ini, nx, ny, nz, ocean, omega, omega_sor, & 170 170 outflow_damping_width, overshoot_limit_e, overshoot_limit_pt, & -
palm/trunk/SOURCE/poisfft.f90
r878 r940 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! special handling of tri-array as an argument in tridia_1dd routines switched 7 ! off because it caused segmentation faults with intel 12.1 compiler 6 8 ! 7 9 ! Former revisions: … … 1381 1383 ! tridia) 1382 1384 ! 1383 ! Attention: when using the intel compiler, array tri must be passed as an 1384 ! argument to the contained subroutines. Otherwise addres faults 1385 ! will occur. 1385 ! Attention: when using the intel compilers older than 12.0, array tri must 1386 ! be passed as an argument to the contained subroutines. Otherwise 1387 ! addres faults will occur. This feature can be activated with 1388 ! cpp-switch __intel11 1386 1389 ! On NEC, tri should not be passed (except for routine substi_1dd) 1387 1390 ! because this causes very bad performance. … … 1425 1428 1426 1429 IF ( j <= nnyh ) THEN 1427 #if defined( __ lc)1430 #if defined( __intel11 ) 1428 1431 CALL maketri_1dd( j, tri ) 1429 1432 #else … … 1431 1434 #endif 1432 1435 ELSE 1433 #if defined( __ lc)1436 #if defined( __intel11 ) 1434 1437 CALL maketri_1dd( ny+1-j, tri ) 1435 1438 #else … … 1437 1440 #endif 1438 1441 ENDIF 1439 #if defined( __ lc)1442 #if defined( __intel11 ) 1440 1443 CALL split_1dd( tri ) 1441 1444 #else … … 1446 1449 CONTAINS 1447 1450 1448 #if defined( __ lc)1451 #if defined( __intel11 ) 1449 1452 SUBROUTINE maketri_1dd( j, tri ) 1450 1453 #else … … 1465 1468 REAL, DIMENSION(0:nx) :: l 1466 1469 1467 #if defined( __ lc)1470 #if defined( __intel11 ) 1468 1471 REAL, DIMENSION(5,0:nx,0:nz-1) :: tri 1469 1472 #endif … … 1512 1515 1513 1516 1514 #if defined( __ lc)1517 #if defined( __intel11 ) 1515 1518 SUBROUTINE split_1dd( tri ) 1516 1519 #else … … 1526 1529 INTEGER :: i, k 1527 1530 1528 #if defined( __ lc)1531 #if defined( __intel11 ) 1529 1532 REAL, DIMENSION(5,0:nx,0:nz-1) :: tri 1530 1533 #endif -
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 -
palm/trunk/SOURCE/prognostic_equations.f90
r786 r940 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! temperature equation can be switched off 7 7 ! 8 8 ! Former revisions: … … 159 159 !-- Calculate those variables needed in the tendency terms which need 160 160 !-- global communication 161 CALL calc_mean_profile( pt, 4 )162 IF ( ocean ) CALL calc_mean_profile( rho, 64 )163 IF ( humidity ) CALL calc_mean_profile( vpt, 44 )161 IF ( .NOT. neutral ) CALL calc_mean_profile( pt, 4 ) 162 IF ( ocean ) CALL calc_mean_profile( rho, 64 ) 163 IF ( humidity ) CALL calc_mean_profile( vpt, 44 ) 164 164 IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. & 165 165 intermediate_timestep_count == 1 ) CALL ws_statistics … … 205 205 ENDIF 206 206 CALL coriolis( i, j, 1 ) 207 IF ( sloping_surface ) CALL buoyancy( i, j, pt, pt_reference, 1, 4 ) 207 IF ( sloping_surface .AND. .NOT. neutral ) THEN 208 CALL buoyancy( i, j, pt, pt_reference, 1, 4 ) 209 ENDIF 208 210 209 211 ! … … 375 377 ENDIF 376 378 CALL coriolis( i, j, 3 ) 377 IF ( ocean ) THEN 378 CALL buoyancy( i, j, rho, rho_reference, 3, 64 ) 379 ELSE 380 IF ( .NOT. humidity ) THEN 381 CALL buoyancy( i, j, pt, pt_reference, 3, 4 ) 382 ELSE 383 CALL buoyancy( i, j, vpt, pt_reference, 3, 44 ) 379 380 IF ( .NOT. neutral ) THEN 381 IF ( ocean ) THEN 382 CALL buoyancy( i, j, rho, rho_reference, 3, 64 ) 383 ELSE 384 IF ( .NOT. humidity ) THEN 385 CALL buoyancy( i, j, pt, pt_reference, 3, 4 ) 386 ELSE 387 CALL buoyancy( i, j, vpt, pt_reference, 3, 44 ) 388 ENDIF 384 389 ENDIF 385 390 ENDIF … … 422 427 423 428 ! 424 !-- potential temperature 425 CALL cpu_log( log_point(13), 'pt-equation', 'start' ) 426 427 ! 428 !-- pt-tendency terms with communication 429 sat = tsc(1) 430 sbt = tsc(2) 431 IF ( scalar_advec == 'bc-scheme' ) THEN 432 433 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 434 ! 435 !-- Bott-Chlond scheme always uses Euler time step when leapfrog is 436 !-- switched on. Thus: 437 sat = 1.0 438 sbt = 1.0 439 ENDIF 440 tend = 0.0 441 CALL advec_s_bc( pt, 'pt' ) 442 ELSE 443 IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN 429 !-- If required, compute prognostic equation for potential temperature 430 IF ( .NOT. neutral ) THEN 431 432 CALL cpu_log( log_point(13), 'pt-equation', 'start' ) 433 434 ! 435 !-- pt-tendency terms with communication 436 sat = tsc(1) 437 sbt = tsc(2) 438 IF ( scalar_advec == 'bc-scheme' ) THEN 439 440 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 441 ! 442 !-- Bott-Chlond scheme always uses Euler time step when leapfrog is 443 !-- switched on. Thus: 444 sat = 1.0 445 sbt = 1.0 446 ENDIF 444 447 tend = 0.0 445 CALL advec_s_ups( pt, 'pt' ) 446 ENDIF 447 ENDIF 448 449 ! 450 !-- pt-tendency terms with no communication 451 DO i = nxl, nxr 452 DO j = nys, nyn 453 ! 454 !-- Tendency terms 455 IF ( scalar_advec == 'bc-scheme' ) THEN 456 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, & 457 wall_heatflux, tend ) 458 ELSE 459 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 460 tend(:,j,i) = 0.0 461 IF ( ws_scheme_sca ) THEN 462 CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, & 463 diss_s_pt, flux_l_pt, diss_l_pt, i_omp_start, tn ) 464 ELSE 465 CALL advec_s_pw( i, j, pt ) 466 ENDIF 467 ELSE 468 IF ( scalar_advec /= 'ups-scheme' ) THEN 469 tend(:,j,i) = 0.0 470 CALL advec_s_up( i, j, pt ) 471 ENDIF 472 ENDIF 473 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 474 THEN 475 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, & 476 tswst_m, wall_heatflux, tend ) 477 ELSE 448 CALL advec_s_bc( pt, 'pt' ) 449 ELSE 450 IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN 451 tend = 0.0 452 CALL advec_s_ups( pt, 'pt' ) 453 ENDIF 454 ENDIF 455 456 ! 457 !-- pt-tendency terms with no communication 458 DO i = nxl, nxr 459 DO j = nys, nyn 460 ! 461 !-- Tendency terms 462 IF ( scalar_advec == 'bc-scheme' ) THEN 478 463 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, & 479 464 wall_heatflux, tend ) 480 ENDIF 481 ENDIF 482 483 ! 484 !-- If required compute heating/cooling due to long wave radiation 485 !-- processes 486 IF ( radiation ) THEN 487 CALL calc_radiation( i, j ) 488 ENDIF 489 490 ! 491 !-- If required compute impact of latent heat due to precipitation 492 IF ( precipitation ) THEN 493 CALL impact_of_latent_heat( i, j ) 494 ENDIF 495 496 ! 497 !-- Consideration of heat sources within the plant canopy 498 IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN 499 CALL plant_canopy_model( i, j, 4 ) 500 ENDIF 501 502 ! 503 !-- If required compute influence of large-scale subsidence/ascent 504 IF ( large_scale_subsidence ) THEN 505 CALL subsidence ( i, j, tend, pt, pt_init ) 506 ENDIF 507 508 CALL user_actions( i, j, 'pt-tendency' ) 509 510 ! 511 !-- Prognostic equation for potential temperature 512 DO k = nzb_s_inner(j,i)+1, nzt 513 pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + & 514 dt_3d * ( & 515 sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) & 516 ) - & 517 tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) ) 465 ELSE 466 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) & 467 THEN 468 tend(:,j,i) = 0.0 469 IF ( ws_scheme_sca ) THEN 470 CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, & 471 diss_s_pt, flux_l_pt, diss_l_pt, & 472 i_omp_start, tn ) 473 ELSE 474 CALL advec_s_pw( i, j, pt ) 475 ENDIF 476 ELSE 477 IF ( scalar_advec /= 'ups-scheme' ) THEN 478 tend(:,j,i) = 0.0 479 CALL advec_s_up( i, j, pt ) 480 ENDIF 481 ENDIF 482 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 483 THEN 484 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, & 485 tswst_m, wall_heatflux, tend ) 486 ELSE 487 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, & 488 wall_heatflux, tend ) 489 ENDIF 490 ENDIF 491 492 ! 493 !-- If required compute heating/cooling due to long wave radiation 494 !-- processes 495 IF ( radiation ) THEN 496 CALL calc_radiation( i, j ) 497 ENDIF 498 499 ! 500 !-- If required compute impact of latent heat due to precipitation 501 IF ( precipitation ) THEN 502 CALL impact_of_latent_heat( i, j ) 503 ENDIF 504 505 ! 506 !-- Consideration of heat sources within the plant canopy 507 IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN 508 CALL plant_canopy_model( i, j, 4 ) 509 ENDIF 510 511 ! 512 !-- If required compute influence of large-scale subsidence/ascent 513 IF ( large_scale_subsidence ) THEN 514 CALL subsidence( i, j, tend, pt, pt_init ) 515 ENDIF 516 517 CALL user_actions( i, j, 'pt-tendency' ) 518 519 ! 520 !-- Prognostic equation for potential temperature 521 DO k = nzb_s_inner(j,i)+1, nzt 522 pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + & 523 dt_3d * ( & 524 sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) & 525 ) - & 526 tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) ) 527 ENDDO 528 529 ! 530 !-- Calculate tendencies for the next Runge-Kutta step 531 IF ( timestep_scheme(1:5) == 'runge' ) THEN 532 IF ( intermediate_timestep_count == 1 ) THEN 533 DO k = nzb_s_inner(j,i)+1, nzt 534 tpt_m(k,j,i) = tend(k,j,i) 535 ENDDO 536 ELSEIF ( intermediate_timestep_count < & 537 intermediate_timestep_count_max ) THEN 538 DO k = nzb_s_inner(j,i)+1, nzt 539 tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + & 540 5.3125 * tpt_m(k,j,i) 541 ENDDO 542 ENDIF 543 ENDIF 544 518 545 ENDDO 519 520 !521 !-- Calculate tendencies for the next Runge-Kutta step522 IF ( timestep_scheme(1:5) == 'runge' ) THEN523 IF ( intermediate_timestep_count == 1 ) THEN524 DO k = nzb_s_inner(j,i)+1, nzt525 tpt_m(k,j,i) = tend(k,j,i)526 ENDDO527 ELSEIF ( intermediate_timestep_count < &528 intermediate_timestep_count_max ) THEN529 DO k = nzb_s_inner(j,i)+1, nzt530 tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)531 ENDDO532 ENDIF533 ENDIF534 535 546 ENDDO 536 ENDDO 537 538 CALL cpu_log( log_point(13), 'pt-equation', 'stop' ) 547 548 CALL cpu_log( log_point(13), 'pt-equation', 'stop' ) 549 550 ENDIF 539 551 540 552 ! … … 714 726 ! 715 727 !-- If required compute influence of large-scale subsidence/ascent 716 IF ( large_scale_subsidence ) THEN717 CALL subsidence 728 IF ( large_scale_subsidence ) THEN 729 CALL subsidence( i, j, tend, q, q_init ) 718 730 ENDIF 719 731 … … 936 948 !-- Calculate those variables needed in the tendency terms which need 937 949 !-- global communication 938 CALL calc_mean_profile( pt, 4 )939 IF ( ocean ) CALL calc_mean_profile( rho, 64 )940 IF ( humidity ) CALL calc_mean_profile( vpt, 44 )950 IF ( .NOT. neutral ) CALL calc_mean_profile( pt, 4 ) 951 IF ( ocean ) CALL calc_mean_profile( rho, 64 ) 952 IF ( humidity ) CALL calc_mean_profile( vpt, 44 ) 941 953 IF ( .NOT. constant_diffusion ) CALL production_e_init 942 954 IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. & … … 989 1001 ENDIF 990 1002 CALL coriolis( i, j, 1 ) 991 IF ( sloping_surface ) CALL buoyancy( i, j, pt, pt_reference, 1, & 992 4 ) 1003 IF ( sloping_surface .AND. .NOT. neutral ) THEN 1004 CALL buoyancy( i, j, pt, pt_reference, 1, 4 ) 1005 ENDIF 993 1006 994 1007 ! … … 1121 1134 ENDIF 1122 1135 CALL coriolis( i, j, 3 ) 1123 IF ( ocean ) THEN 1124 CALL buoyancy( i, j, rho, rho_reference, 3, 64 ) 1125 ELSE 1126 IF ( .NOT. humidity ) THEN 1127 CALL buoyancy( i, j, pt, pt_reference, 3, 4 ) 1128 ELSE 1129 CALL buoyancy( i, j, vpt, pt_reference, 3, 44 ) 1136 1137 IF ( .NOT. neutral ) THEN 1138 IF ( ocean ) THEN 1139 CALL buoyancy( i, j, rho, rho_reference, 3, 64 ) 1140 ELSE 1141 IF ( .NOT. humidity ) THEN 1142 CALL buoyancy( i, j, pt, pt_reference, 3, 4 ) 1143 ELSE 1144 CALL buoyancy( i, j, vpt, pt_reference, 3, 44 ) 1145 ENDIF 1130 1146 ENDIF 1131 1147 ENDIF … … 1163 1179 1164 1180 ! 1165 !-- Tendency terms for potential temperature 1166 tend(:,j,i) = 0.0 1167 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 1168 IF ( ws_scheme_sca ) THEN 1169 ! CALL local_diss( i, j, pt ) 1170 CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, & 1171 diss_s_pt, flux_l_pt, diss_l_pt, i_omp_start, tn ) 1172 ELSE 1173 CALL advec_s_pw( i, j, pt ) 1174 ENDIF 1175 ELSE 1176 CALL advec_s_up( i, j, pt ) 1177 ENDIF 1178 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 1179 THEN 1180 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, & 1181 tswst_m, wall_heatflux, tend ) 1182 ELSE 1183 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, & 1184 wall_heatflux, tend ) 1185 ENDIF 1186 1187 ! 1188 !-- If required compute heating/cooling due to long wave radiation 1189 !-- processes 1190 IF ( radiation ) THEN 1191 CALL calc_radiation( i, j ) 1192 ENDIF 1193 1194 ! 1195 !-- If required compute impact of latent heat due to precipitation 1196 IF ( precipitation ) THEN 1197 CALL impact_of_latent_heat( i, j ) 1198 ENDIF 1199 1200 ! 1201 !-- Consideration of heat sources within the plant canopy 1202 IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN 1203 CALL plant_canopy_model( i, j, 4 ) 1204 ENDIF 1205 1206 1207 !-- If required compute influence of large-scale subsidence/ascent 1208 IF ( large_scale_subsidence ) THEN 1209 CALL subsidence ( i, j, tend, pt, pt_init ) 1210 ENDIF 1211 1212 1213 CALL user_actions( i, j, 'pt-tendency' ) 1214 1215 ! 1216 !-- Prognostic equation for potential temperature 1217 DO k = nzb_s_inner(j,i)+1, nzt 1218 pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +& 1219 dt_3d * ( & 1220 tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) & 1221 ) - & 1222 tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) ) 1223 ENDDO 1224 1225 ! 1226 !-- Calculate tendencies for the next Runge-Kutta step 1227 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1228 IF ( intermediate_timestep_count == 1 ) THEN 1229 DO k = nzb_s_inner(j,i)+1, nzt 1230 tpt_m(k,j,i) = tend(k,j,i) 1231 ENDDO 1232 ELSEIF ( intermediate_timestep_count < & 1233 intermediate_timestep_count_max ) THEN 1234 DO k = nzb_s_inner(j,i)+1, nzt 1235 tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1236 5.3125 * tpt_m(k,j,i) 1237 ENDDO 1238 ENDIF 1181 !-- If required, compute prognostic equation for potential temperature 1182 IF ( .NOT. neutral ) THEN 1183 ! 1184 !-- Tendency terms for potential temperature 1185 tend(:,j,i) = 0.0 1186 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 1187 IF ( ws_scheme_sca ) THEN 1188 CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, diss_s_pt, & 1189 flux_l_pt, diss_l_pt, i_omp_start, tn ) 1190 ELSE 1191 CALL advec_s_pw( i, j, pt ) 1192 ENDIF 1193 ELSE 1194 CALL advec_s_up( i, j, pt ) 1195 ENDIF 1196 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 1197 THEN 1198 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, & 1199 tswst_m, wall_heatflux, tend ) 1200 ELSE 1201 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, & 1202 wall_heatflux, tend ) 1203 ENDIF 1204 1205 ! 1206 !-- If required compute heating/cooling due to long wave radiation 1207 !-- processes 1208 IF ( radiation ) THEN 1209 CALL calc_radiation( i, j ) 1210 ENDIF 1211 1212 ! 1213 !-- If required compute impact of latent heat due to precipitation 1214 IF ( precipitation ) THEN 1215 CALL impact_of_latent_heat( i, j ) 1216 ENDIF 1217 1218 ! 1219 !-- Consideration of heat sources within the plant canopy 1220 IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN 1221 CALL plant_canopy_model( i, j, 4 ) 1222 ENDIF 1223 1224 ! 1225 !-- If required, compute influence of large-scale subsidence/ascent 1226 IF ( large_scale_subsidence ) THEN 1227 CALL subsidence( i, j, tend, pt, pt_init ) 1228 ENDIF 1229 1230 1231 CALL user_actions( i, j, 'pt-tendency' ) 1232 1233 ! 1234 !-- Prognostic equation for potential temperature 1235 DO k = nzb_s_inner(j,i)+1, nzt 1236 pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +& 1237 dt_3d * ( & 1238 tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) & 1239 ) - & 1240 tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) ) 1241 ENDDO 1242 1243 ! 1244 !-- Calculate tendencies for the next Runge-Kutta step 1245 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1246 IF ( intermediate_timestep_count == 1 ) THEN 1247 DO k = nzb_s_inner(j,i)+1, nzt 1248 tpt_m(k,j,i) = tend(k,j,i) 1249 ENDDO 1250 ELSEIF ( intermediate_timestep_count < & 1251 intermediate_timestep_count_max ) THEN 1252 DO k = nzb_s_inner(j,i)+1, nzt 1253 tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1254 5.3125 * tpt_m(k,j,i) 1255 ENDDO 1256 ENDIF 1257 ENDIF 1258 1239 1259 ENDIF 1240 1260 … … 1337 1357 1338 1358 !-- If required compute influence of large-scale subsidence/ascent 1339 IF ( large_scale_subsidence ) THEN1340 CALL subsidence 1359 IF ( large_scale_subsidence ) THEN 1360 CALL subsidence( i, j, tend, q, q_init ) 1341 1361 ENDIF 1342 1362 … … 1484 1504 !-- Calculate those variables needed in the tendency terms which need 1485 1505 !-- global communication 1486 CALL calc_mean_profile( pt, 4 )1487 IF ( ocean ) CALL calc_mean_profile( rho, 64 )1488 IF ( humidity ) CALL calc_mean_profile( vpt, 44 )1506 IF ( .NOT. neutral ) CALL calc_mean_profile( pt, 4 ) 1507 IF ( ocean ) CALL calc_mean_profile( rho, 64 ) 1508 IF ( humidity ) CALL calc_mean_profile( vpt, 44 ) 1489 1509 IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. & 1490 1510 intermediate_timestep_count == 1 ) CALL ws_statistics … … 1523 1543 ENDIF 1524 1544 CALL coriolis( 1 ) 1525 IF ( sloping_surface ) CALL buoyancy( pt, pt_reference, 1, 4 ) 1545 IF ( sloping_surface .AND. .NOT. neutral ) THEN 1546 CALL buoyancy( pt, pt_reference, 1, 4 ) 1547 ENDIF 1526 1548 1527 1549 ! … … 1706 1728 ENDIF 1707 1729 CALL coriolis( 3 ) 1708 IF ( ocean ) THEN 1709 CALL buoyancy( rho, rho_reference, 3, 64 ) 1710 ELSE 1711 IF ( .NOT. humidity ) THEN 1712 CALL buoyancy( pt, pt_reference, 3, 4 ) 1730 1731 IF ( .NOT. neutral ) THEN 1732 IF ( ocean ) THEN 1733 CALL buoyancy( rho, rho_reference, 3, 64 ) 1713 1734 ELSE 1714 CALL buoyancy( vpt, pt_reference, 3, 44 ) 1735 IF ( .NOT. humidity ) THEN 1736 CALL buoyancy( pt, pt_reference, 3, 4 ) 1737 ELSE 1738 CALL buoyancy( vpt, pt_reference, 3, 44 ) 1739 ENDIF 1715 1740 ENDIF 1716 1741 ENDIF … … 1761 1786 CALL cpu_log( log_point(7), 'w-equation', 'stop' ) 1762 1787 1763 ! 1764 !-- potential temperature 1765 CALL cpu_log( log_point(13), 'pt-equation', 'start' ) 1766 1767 ! 1768 !-- pt-tendency terms with communication 1769 sat = tsc(1) 1770 sbt = tsc(2) 1771 IF ( scalar_advec == 'bc-scheme' ) THEN 1772 1773 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 1774 ! 1775 !-- Bott-Chlond scheme always uses Euler time step when leapfrog is 1776 !-- switched on. Thus: 1777 sat = 1.0 1778 sbt = 1.0 1779 ENDIF 1780 tend = 0.0 1781 CALL advec_s_bc( pt, 'pt' ) 1782 ELSE 1783 IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN 1788 1789 ! 1790 !-- If required, compute prognostic equation for potential temperature 1791 IF ( .NOT. neutral ) THEN 1792 1793 CALL cpu_log( log_point(13), 'pt-equation', 'start' ) 1794 1795 ! 1796 !-- pt-tendency terms with communication 1797 sat = tsc(1) 1798 sbt = tsc(2) 1799 IF ( scalar_advec == 'bc-scheme' ) THEN 1800 1801 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 1802 ! 1803 !-- Bott-Chlond scheme always uses Euler time step when leapfrog is 1804 !-- switched on. Thus: 1805 sat = 1.0 1806 sbt = 1.0 1807 ENDIF 1784 1808 tend = 0.0 1785 CALL advec_s_ups( pt, 'pt' ) 1786 ENDIF 1787 ENDIF 1788 1789 ! 1790 !-- pt-tendency terms with no communication 1791 IF ( scalar_advec == 'bc-scheme' ) THEN 1792 CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, & 1793 tend ) 1794 ELSE 1795 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 1796 tend = 0.0 1797 IF ( ws_scheme_sca ) THEN 1798 CALL advec_s_ws( pt, 'pt' ) 1799 ELSE 1800 CALL advec_s_pw( pt ) 1801 ENDIF 1809 CALL advec_s_bc( pt, 'pt' ) 1802 1810 ELSE 1803 IF ( scalar_advec /= 'ups-scheme' ) THEN1811 IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN 1804 1812 tend = 0.0 1805 CALL advec_s_up ( pt)1806 ENDIF 1807 ENDIF 1808 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 1809 CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tswst_m, & 1810 wall_heatflux, tend ) 1811 ELSE1813 CALL advec_s_ups( pt, 'pt' ) 1814 ENDIF 1815 ENDIF 1816 1817 ! 1818 !-- pt-tendency terms with no communication 1819 IF ( scalar_advec == 'bc-scheme' ) THEN 1812 1820 CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, & 1813 1821 tend ) 1814 ENDIF 1815 ENDIF 1816 1817 ! 1818 !-- If required compute heating/cooling due to long wave radiation 1819 !-- processes 1820 IF ( radiation ) THEN 1821 CALL calc_radiation 1822 ENDIF 1823 1824 ! 1825 !-- If required compute impact of latent heat due to precipitation 1826 IF ( precipitation ) THEN 1827 CALL impact_of_latent_heat 1828 ENDIF 1829 1830 ! 1831 !-- Consideration of heat sources within the plant canopy 1832 IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN 1833 CALL plant_canopy_model( 4 ) 1834 ENDIF 1835 1836 !--If required compute influence of large-scale subsidence/ascent 1837 IF ( large_scale_subsidence ) THEN 1838 CALL subsidence ( tend, pt, pt_init ) 1839 ENDIF 1840 1841 CALL user_actions( 'pt-tendency' ) 1842 1843 ! 1844 !-- Prognostic equation for potential temperature 1845 DO i = nxl, nxr 1846 DO j = nys, nyn 1847 DO k = nzb_s_inner(j,i)+1, nzt 1848 pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + & 1849 dt_3d * ( & 1850 sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) & 1851 ) - & 1852 tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) ) 1822 ELSE 1823 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 1824 tend = 0.0 1825 IF ( ws_scheme_sca ) THEN 1826 CALL advec_s_ws( pt, 'pt' ) 1827 ELSE 1828 CALL advec_s_pw( pt ) 1829 ENDIF 1830 ELSE 1831 IF ( scalar_advec /= 'ups-scheme' ) THEN 1832 tend = 0.0 1833 CALL advec_s_up( pt ) 1834 ENDIF 1835 ENDIF 1836 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 1837 CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tswst_m, & 1838 wall_heatflux, tend ) 1839 ELSE 1840 CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, & 1841 tend ) 1842 ENDIF 1843 ENDIF 1844 1845 ! 1846 !-- If required compute heating/cooling due to long wave radiation processes 1847 IF ( radiation ) THEN 1848 CALL calc_radiation 1849 ENDIF 1850 1851 ! 1852 !-- If required compute impact of latent heat due to precipitation 1853 IF ( precipitation ) THEN 1854 CALL impact_of_latent_heat 1855 ENDIF 1856 1857 ! 1858 !-- Consideration of heat sources within the plant canopy 1859 IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN 1860 CALL plant_canopy_model( 4 ) 1861 ENDIF 1862 1863 ! 1864 !-- If required compute influence of large-scale subsidence/ascent 1865 IF ( large_scale_subsidence ) THEN 1866 CALL subsidence( tend, pt, pt_init ) 1867 ENDIF 1868 1869 CALL user_actions( 'pt-tendency' ) 1870 1871 ! 1872 !-- Prognostic equation for potential temperature 1873 DO i = nxl, nxr 1874 DO j = nys, nyn 1875 DO k = nzb_s_inner(j,i)+1, nzt 1876 pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + & 1877 dt_3d * ( & 1878 sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) & 1879 ) - & 1880 tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) ) 1881 ENDDO 1853 1882 ENDDO 1854 1883 ENDDO 1855 ENDDO 1856 1857 ! 1858 !-- Calculate tendencies for the next Runge-Kutta step 1859 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1860 IF ( intermediate_timestep_count == 1 ) THEN 1861 DO i = nxl, nxr 1862 DO j = nys, nyn 1863 DO k = nzb_s_inner(j,i)+1, nzt 1864 tpt_m(k,j,i) = tend(k,j,i) 1865 ENDDO 1866 ENDDO 1867 ENDDO 1868 ELSEIF ( intermediate_timestep_count < & 1869 intermediate_timestep_count_max ) THEN 1870 DO i = nxl, nxr 1871 DO j = nys, nyn 1872 DO k = nzb_s_inner(j,i)+1, nzt 1873 tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i) 1874 ENDDO 1875 ENDDO 1876 ENDDO 1877 ENDIF 1878 ENDIF 1879 1880 CALL cpu_log( log_point(13), 'pt-equation', 'stop' ) 1884 1885 ! 1886 !-- Calculate tendencies for the next Runge-Kutta step 1887 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1888 IF ( intermediate_timestep_count == 1 ) THEN 1889 DO i = nxl, nxr 1890 DO j = nys, nyn 1891 DO k = nzb_s_inner(j,i)+1, nzt 1892 tpt_m(k,j,i) = tend(k,j,i) 1893 ENDDO 1894 ENDDO 1895 ENDDO 1896 ELSEIF ( intermediate_timestep_count < & 1897 intermediate_timestep_count_max ) THEN 1898 DO i = nxl, nxr 1899 DO j = nys, nyn 1900 DO k = nzb_s_inner(j,i)+1, nzt 1901 tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1902 5.3125 * tpt_m(k,j,i) 1903 ENDDO 1904 ENDDO 1905 ENDDO 1906 ENDIF 1907 ENDIF 1908 1909 CALL cpu_log( log_point(13), 'pt-equation', 'stop' ) 1910 1911 ENDIF 1881 1912 1882 1913 ! … … 2054 2085 ! 2055 2086 !-- If required compute influence of large-scale subsidence/ascent 2056 IF ( large_scale_subsidence ) THEN2057 CALL subsidence 2087 IF ( large_scale_subsidence ) THEN 2088 CALL subsidence( tend, q, q_init ) 2058 2089 ENDIF 2059 2090 -
palm/trunk/SOURCE/read_var_list.f90
r928 r940 4 4 ! Current revisions: 5 5 ! ------------------ 6 ! 6 ! +neutral 7 7 ! 8 8 ! Former revisions: … … 421 421 CASE ( 'netcdf_precision' ) 422 422 READ ( 13 ) netcdf_precision 423 CASE ( 'neutral' ) 424 READ ( 13 ) neutral 423 425 CASE ( 'ngsrb' ) 424 426 READ ( 13 ) ngsrb -
palm/trunk/SOURCE/write_var_list.f90
r928 r940 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! +neutral 7 7 ! 8 8 ! Former revisions: … … 335 335 WRITE ( 14 ) 'netcdf_precision ' 336 336 WRITE ( 14 ) netcdf_precision 337 WRITE ( 14 ) 'neutral ' 338 WRITE ( 14 ) neutral 337 339 WRITE ( 14 ) 'ngsrb ' 338 340 WRITE ( 14 ) ngsrb
Note: See TracChangeset
for help on using the changeset viewer.
