Changeset 940


Ignore:
Timestamp:
Jul 9, 2012 2:31:00 PM (12 years ago)
Author:
raasch
Message:

temperature equation can be switched off; bugfix of tridia_1dd for current Intel (12.1) compilers

Location:
palm/trunk/SOURCE
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/check_parameters.f90

    r925 r940  
    44! Current revisions:
    55! -----------------
    6 !
     6! checks for parameter neutral
    77!
    88! Former revisions:
     
    10391039!
    10401040!--    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
    10541057                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
    10591064                ELSE
    1060                    pt_init(k) = pt_surface   + 0.5 * dzu(k) * gradient
     1065                   pt_init(k) = pt_init(k-1)
    10611066                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
    10771080                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
    10821088                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)
    10851090                ENDIF
    1086              ELSE
    1087                 pt_init(k) = pt_init(k+1)
    1088              ENDIF
    1089           ENDDO
     1091             ENDDO
     1092
     1093          ENDIF
    10901094
    10911095       ENDIF
     
    15051509    IF ( surface_heatflux == 9999999.9 )  constant_heatflux     = .FALSE.
    15061510    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
    15071528    IF ( top_momentumflux_u /= 9999999.9  .AND.  &
    15081529         top_momentumflux_v /= 9999999.9 )  THEN
  • palm/trunk/SOURCE/header.f90

    r928 r940  
    44! Current revisions:
    55! -----------------
    6 !
     6! Output in case of simulations for pure neutral stratification (no pt-equation
     7! solved)
    78!
    89! Former revisions:
     
    349350       ENDIF
    350351    ENDIF
     352    IF ( neutral )  WRITE ( io, 131 )  pt_surface
    351353    IF ( humidity )  THEN
    352354       IF ( .NOT. cloud_physics )  THEN
     
    16821684129 FORMAT (' --> Additional prognostic equation for the specific humidity')
    16831685130 FORMAT (' --> Additional prognostic equation for the total water content')
     1686131 FORMAT (' --> No pt-equation solved. Neutral stratification with pt = ', &
     1687                  F6.2, ' K assumed')
    16841688132 FORMAT ('     Parameterization of long-wave radiation processes via'/ &
    16851689            '     effective emissivity scheme')
  • palm/trunk/SOURCE/modules.f90

    r937 r940  
    44! Current revisions:
    55! -----------------
    6 !
     6! +neutral
    77!
    88! Former revisions:
     
    591591                iso2d_output = .FALSE., large_scale_subsidence = .FALSE., &
    592592                masking_method = .FALSE., mg_switch_to_pe0 = .FALSE., &
    593                 netcdf_output = .FALSE., ocean = .FALSE., &
     593                netcdf_output = .FALSE., neutral = .FALSE., ocean = .FALSE., &
    594594                outflow_l = .FALSE., outflow_n = .FALSE., outflow_r = .FALSE., &
    595595                outflow_s = .FALSE., passive_scalar = .FALSE., &
  • palm/trunk/SOURCE/parin.f90

    r937 r940  
    44! Current revisions:
    55! -----------------
    6 !
     6! +neutral in inipar
    77!
    88! Former revisions:
     
    165165             leaf_surface_concentration, long_filter_factor, &
    166166             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, &
    169169             nsor_ini, nx, ny, nz, ocean, omega, omega_sor, &
    170170             outflow_damping_width, overshoot_limit_e, overshoot_limit_pt, &
  • palm/trunk/SOURCE/poisfft.f90

    r878 r940  
    44! Current revisions:
    55! -----------------
     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
    68!
    79! Former revisions:
     
    13811383! tridia)
    13821384!
    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
    13861389!             On NEC, tri should not be passed (except for routine substi_1dd)
    13871390!             because this causes very bad performance.
     
    14251428
    14261429       IF ( j <= nnyh )  THEN
    1427 #if defined( __lc )
     1430#if defined( __intel11 )
    14281431          CALL maketri_1dd( j, tri )
    14291432#else
     
    14311434#endif
    14321435       ELSE
    1433 #if defined( __lc )
     1436#if defined( __intel11 )
    14341437          CALL maketri_1dd( ny+1-j, tri )
    14351438#else
     
    14371440#endif
    14381441       ENDIF
    1439 #if defined( __lc )
     1442#if defined( __intel11 )
    14401443       CALL split_1dd( tri )
    14411444#else
     
    14461449    CONTAINS
    14471450
    1448 #if defined( __lc )
     1451#if defined( __intel11 )
    14491452       SUBROUTINE maketri_1dd( j, tri )
    14501453#else
     
    14651468          REAL, DIMENSION(0:nx) ::  l
    14661469
    1467 #if defined( __lc )
     1470#if defined( __intel11 )
    14681471          REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
    14691472#endif
     
    15121515
    15131516
    1514 #if defined( __lc )
     1517#if defined( __intel11 )
    15151518       SUBROUTINE split_1dd( tri )
    15161519#else
     
    15261529          INTEGER ::  i, k
    15271530
    1528 #if defined( __lc )
     1531#if defined( __intel11 )
    15291532          REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
    15301533#endif
  • palm/trunk/SOURCE/production_e.f90

    r760 r940  
    44! Current revisions:
    55! -----------------
    6 !
     6! TKE production by buoyancy can be switched off in case of runs with pure
     7! neutral stratification
    78!
    89! Former revisions:
     
    127128!       ENDIF
    128129
    129 !
    130 !--    Calculate TKE production by shear
     130
    131131       DO  i = nxl, nxr
    132132
     133!
     134!--       Calculate TKE production by shear
    133135          DO  j = nys, nyn
    134136             DO  k = nzb_diff_s_outer(j,i), nzt
     
    453455
    454456!
    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
    469475                      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
    471501
    472502                ELSE
    473503
    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
    479515                      ENDDO
    480516
    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
    492541
    493542                ENDIF
     
    495544             ELSE
    496545
    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
    544551                         k1 = 1.0 + 0.61 * q(k,j,i)
    545552                         k2 = 0.61 * pt(k,j,i)
    546553                      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
    555567                      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
    562575                ENDDO
    563576
    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
    577584                         k1 = 1.0 + 0.61 * q(k,j,i)
    578585                         k2 = 0.61 * pt(k,j,i)
    579586                      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
    588600                      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
    608615                         k1 = 1.0 + 0.61 * q(k,j,i)
    609616                         k2 = 0.61 * pt(k,j,i)
    610617                      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
    619631                      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
    625638
    626639             ENDIF
     
    950963
    951964!
    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
    9651001
    9661002             ELSE
    9671003
    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
    9811032                ENDIF
    9821033
     
    9851036          ELSE
    9861037
    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
    10261041                   k1 = 1.0 + 0.61 * q(k,j,i)
    10271042                   k2 = 0.61 * pt(k,j,i)
    10281043                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                 ENDIF
    1038              ENDIF
    1039 
    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           ENDDO
    1045 
    1046           IF ( use_surface_fluxes )  THEN
    1047              k = nzb_diff_s_inner(j,i)-1
    1048 
    1049              IF ( .NOT. cloud_physics ) THEN
    1050                 k1 = 1.0 + 0.61 * q(k,j,i)
    1051                 k2 = 0.61 * pt(k,j,i)
    1052              ELSE
    1053                 IF ( ql(k,j,i) == 0.0 ) THEN
     1044                   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
    10541069                   k1 = 1.0 + 0.61 * q(k,j,i)
    10551070                   k2 = 0.61 * pt(k,j,i)
    10561071                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) )
    10661089             ENDIF
    10671090
    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
    10801095                   k1 = 1.0 + 0.61 * q(k,j,i)
    10811096                   k2 = 0.61 * pt(k,j,i)
    10821097                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) )
    10921115             ENDIF
    10931116
    1094              tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
    1095                                          ( k1* tswst(j,i) + k2 * qswst(j,i) )
    10961117          ENDIF
    10971118
  • palm/trunk/SOURCE/prognostic_equations.f90

    r786 r940  
    44! Current revisions:
    55! -----------------
    6 !
     6! temperature equation can be switched off
    77!
    88! Former revisions:
     
    159159!-- Calculate those variables needed in the tendency terms which need
    160160!-- 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 )
    164164    IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
    165165         intermediate_timestep_count == 1 )  CALL ws_statistics
     
    205205          ENDIF
    206206          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
    208210
    209211!
     
    375377          ENDIF
    376378          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
    384389             ENDIF
    385390          ENDIF
     
    422427
    423428!
    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
    444447          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
    478463                CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
    479464                                  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
    518545          ENDDO
    519 
    520 !
    521 !--       Calculate tendencies for the next Runge-Kutta step
    522           IF ( timestep_scheme(1:5) == 'runge' )  THEN
    523              IF ( intermediate_timestep_count == 1 )  THEN
    524                 DO  k = nzb_s_inner(j,i)+1, nzt
    525                    tpt_m(k,j,i) = tend(k,j,i)
    526                 ENDDO
    527              ELSEIF ( intermediate_timestep_count < &
    528                       intermediate_timestep_count_max )  THEN
    529                 DO  k = nzb_s_inner(j,i)+1, nzt
    530                    tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
    531                 ENDDO
    532              ENDIF
    533           ENDIF
    534 
    535546       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
    539551
    540552!
     
    714726!
    715727!--          If required compute influence of large-scale subsidence/ascent
    716              IF ( large_scale_subsidence ) THEN
    717                 CALL subsidence ( i, j, tend, q, q_init )
     728             IF ( large_scale_subsidence )  THEN
     729                CALL subsidence( i, j, tend, q, q_init )
    718730             ENDIF
    719731
     
    936948!-- Calculate those variables needed in the tendency terms which need
    937949!-- 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 )
    941953    IF ( .NOT. constant_diffusion )  CALL production_e_init
    942954    IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
     
    9891001             ENDIF
    9901002             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
    9931006
    9941007!
     
    11211134          ENDIF
    11221135          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
    11301146             ENDIF
    11311147          ENDIF
     
    11631179
    11641180!
    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
    12391259          ENDIF
    12401260
     
    13371357
    13381358!--          If required compute influence of large-scale subsidence/ascent
    1339              IF ( large_scale_subsidence ) THEN
    1340                 CALL subsidence ( i, j, tend, q, q_init )
     1359             IF ( large_scale_subsidence )  THEN
     1360                CALL subsidence( i, j, tend, q, q_init )
    13411361             ENDIF
    13421362
     
    14841504!-- Calculate those variables needed in the tendency terms which need
    14851505!-- 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 )
    14891509    IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
    14901510         intermediate_timestep_count == 1 )  CALL ws_statistics
     
    15231543    ENDIF
    15241544    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
    15261548
    15271549!
     
    17061728    ENDIF
    17071729    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 )
    17131734       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
    17151740       ENDIF
    17161741    ENDIF
     
    17611786    CALL cpu_log( log_point(7), 'w-equation', 'stop' )
    17621787
    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
    17841808          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' )
    18021810       ELSE
    1803           IF ( scalar_advec /= 'ups-scheme' )  THEN
     1811          IF ( tsc(2) /= 2.0  .AND.  scalar_advec == 'ups-scheme' )  THEN
    18041812             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        ELSE
     1813             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
    18121820          CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, &
    18131821                            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
    18531882          ENDDO
    18541883       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
    18811912
    18821913!
     
    20542085!
    20552086!--    If required compute influence of large-scale subsidence/ascent
    2056        IF ( large_scale_subsidence ) THEN
    2057          CALL subsidence ( tend, q, q_init )
     2087       IF ( large_scale_subsidence )  THEN
     2088         CALL subsidence( tend, q, q_init )
    20582089       ENDIF
    20592090
  • palm/trunk/SOURCE/read_var_list.f90

    r928 r940  
    44! Current revisions:
    55! ------------------
    6 !
     6! +neutral
    77!
    88! Former revisions:
     
    421421          CASE ( 'netcdf_precision' )
    422422             READ ( 13 )  netcdf_precision
     423          CASE ( 'neutral' )
     424             READ ( 13 )  neutral
    423425          CASE ( 'ngsrb' )
    424426             READ ( 13 )  ngsrb
  • palm/trunk/SOURCE/write_var_list.f90

    r928 r940  
    44! Current revisions:
    55! -----------------
    6 !
     6! +neutral
    77!
    88! Former revisions:
     
    335335    WRITE ( 14 )  'netcdf_precision              '
    336336    WRITE ( 14 )  netcdf_precision
     337    WRITE ( 14 )  'neutral                       '
     338    WRITE ( 14 )  neutral
    337339    WRITE ( 14 )  'ngsrb                         '
    338340    WRITE ( 14 )  ngsrb
Note: See TracChangeset for help on using the changeset viewer.