Changeset 1340


Ignore:
Timestamp:
Mar 25, 2014 7:45:13 PM (10 years ago)
Author:
kanani
Message:

REAL constants defined as wp-kind

Location:
palm/trunk/SOURCE
Files:
10 edited

Legend:

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

    r1321 r1340  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants defined as wp-kind
    2323!
    2424! Former revisions:
     
    148148                   dvar_dz = atmos_ocean_sign * &
    149149                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
    150                    IF ( dvar_dz > 0.0 ) THEN
    151                       l_stable = 0.76 * SQRT( e(k,j,i) ) / &
    152                                  SQRT( g / var_reference * dvar_dz ) + 1E-5
     150                   IF ( dvar_dz > 0.0_wp ) THEN
     151                      l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
     152                                    SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
    153153                   ELSE
    154154                      l_stable = l_grid(k)
     
    176176                DO  k = nzb_s_inner(j,i)+1, nzt
    177177
    178                     dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
     178                    dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
    179179                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
    180180
     
    220220                   dvar_dz = atmos_ocean_sign * &
    221221                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
    222                    IF ( dvar_dz > 0.0 ) THEN
    223                       l_stable = 0.76 * SQRT( e(k,j,i) ) / &
    224                                         SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
     222                   IF ( dvar_dz > 0.0_wp ) THEN
     223                      l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
     224                                           SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
    225225                   ELSE
    226226                      l_stable = l_grid(k)
     
    248248                DO  k = nzb_s_inner(j,i)+1, nzt
    249249
    250                     dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
    251                                        e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
     250                    dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
     251                                             e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
    252252
    253253                    tend(k,j,i) = tend(k,j,i)                                  &
     
    356356                      dvar_dz = atmos_ocean_sign * &
    357357                                ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
    358                       IF ( dvar_dz > 0.0 ) THEN
    359                          l_stable = 0.76 * SQRT( e(k,j,i) ) / &
    360                                     SQRT( g / var_reference * dvar_dz ) + 1E-5
     358                      IF ( dvar_dz > 0.0_wp ) THEN
     359                         l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
     360                                       SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
    361361                      ELSE
    362362                         l_stable = l_grid(k)
     
    377377!
    378378!--                   Calculate the tendency terms
    379                       dissipation = ( 0.19 + 0.74 * l / ll ) * &
    380                                     e(k,j,i) * SQRT( e(k,j,i) ) / l
     379                      dissipation = ( 0.19_wp + 0.74_wp * l / ll ) * &
     380                                          e(k,j,i) * SQRT( e(k,j,i) ) / l
    381381
    382382                      tend(k,j,i) = tend(k,j,i)                                  &
     
    423423                      dvar_dz = atmos_ocean_sign * &
    424424                                ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
    425                       IF ( dvar_dz > 0.0 ) THEN
    426                          l_stable = 0.76 * SQRT( e(k,j,i) ) / &
    427                                            SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
     425                      IF ( dvar_dz > 0.0_wp ) THEN
     426                         l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
     427                                              SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
    428428                      ELSE
    429429                         l_stable = l_grid(k)
     
    444444!
    445445!--                   Calculate the tendency terms
    446                       dissipation = ( 0.19 + 0.74 * l / ll ) * &
    447                                     e(k,j,i) * SQRT( e(k,j,i) ) / l
     446                      dissipation = ( 0.19_wp + 0.74_wp * l / ll ) * &
     447                                          e(k,j,i) * SQRT( e(k,j,i) ) / l
    448448
    449449                      tend(k,j,i) = tend(k,j,i)                                &
     
    541541          dvar_dz = atmos_ocean_sign * &
    542542                    ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
    543           IF ( dvar_dz > 0.0 ) THEN
     543          IF ( dvar_dz > 0.0_wp ) THEN
    544544             IF ( use_single_reference_value )  THEN
    545                 l_stable = 0.76 * SQRT( e(k,j,i) ) / &
    546                                   SQRT( g / var_reference * dvar_dz ) + 1E-5
     545                l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
     546                                     SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
    547547             ELSE
    548                 l_stable = 0.76 * SQRT( e(k,j,i) ) / &
    549                                   SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
     548                l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
     549                                     SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
    550550             ENDIF
    551551          ELSE
     
    566566!
    567567!--       Calculate the tendency term
    568           dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
    569                            SQRT( e(k,j,i) ) / l(k)
     568          dissipation(k) = ( 0.19_wp + 0.74_wp * l(k) / ll(k) ) * e(k,j,i) * &
     569                                 SQRT( e(k,j,i) ) / l(k)
    570570
    571571          tend(k,j,i) = tend(k,j,i)                                           &
  • palm/trunk/SOURCE/diffusion_s.f90

    r1321 r1340  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! REAL constants defined as wp-kind
    2323!
    2424! Former revisions:
     
    120120
    121121                tend(k,j,i) = tend(k,j,i)                                      &
    122                                           + 0.5 * (                            &
     122                                          + 0.5_wp * (                         &
    123123                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
    124124                      - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
    125                                                   ) * ddx2                     &
    126                                           + 0.5 * (                            &
     125                                                     ) * ddx2                  &
     126                                          + 0.5_wp * (                         &
    127127                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
    128128                      - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
    129                                                   ) * ddy2
     129                                                     ) * ddy2
    130130             ENDDO
    131131
    132132!
    133133!--          Apply prescribed horizontal wall heatflux where necessary
    134              IF ( ( wall_w_x(j,i) .NE. 0.0 ) .OR. ( wall_w_y(j,i) .NE. 0.0 ) ) &
     134             IF ( ( wall_w_x(j,i) .NE. 0.0_wp ) .OR. ( wall_w_y(j,i) .NE. 0.0_wp ) ) &
    135135             THEN
    136136                DO  k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
    137137
    138138                   tend(k,j,i) = tend(k,j,i)                                   &
    139                                                 + ( fwxp(j,i) * 0.5 *          &
    140                         ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
    141                         + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                 &
    142                                                    -fwxm(j,i) * 0.5 *          &
     139                                                + ( fwxp(j,i) * 0.5_wp *       &
     140                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
     141                        + ( 1.0_wp - fwxp(j,i) ) * wall_s_flux(1)              &
     142                                                   -fwxm(j,i) * 0.5_wp *       &
    143143                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
    144                         + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                 &
     144                        + ( 1.0_wp - fwxm(j,i) ) * wall_s_flux(2)              &
    145145                                                  ) * ddx2                     &
    146                                                 + ( fwyp(j,i) * 0.5 *          &
    147                         ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
    148                         + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                 &
    149                                                    -fwym(j,i) * 0.5 *          &
     146                                                + ( fwyp(j,i) * 0.5_wp *       &
     147                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
     148                        + ( 1.0_wp - fwyp(j,i) ) * wall_s_flux(3)              &
     149                                                   -fwym(j,i) * 0.5_wp *       &
    150150                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
    151                         + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                 &
     151                        + ( 1.0_wp - fwym(j,i) ) * wall_s_flux(4)              &
    152152                                                  ) * ddy2
    153153                ENDDO
     
    161161
    162162                tend(k,j,i) = tend(k,j,i)                                      &
    163                                        + 0.5 * (                               &
     163                                       + 0.5_wp * (                            &
    164164            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
    165165          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
    166                                                ) * ddzw(k)
     166                                                  ) * ddzw(k)
    167167             ENDDO
    168168
     
    175175
    176176                tend(k,j,i) = tend(k,j,i)                                      &
    177                                        + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )     &
    178                                                * ( s(k+1,j,i)-s(k,j,i) )       &
    179                                                * ddzu(k+1)                     &
     177                                       + ( 0.5_wp * ( kh(k,j,i)+kh(k+1,j,i) )  &
     178                                                  * ( s(k+1,j,i)-s(k,j,i) )    &
     179                                                  * ddzu(k+1)                  &
    180180                                           + s_flux_b(j,i)                     &
    181181                                         ) * ddzw(k)
     
    192192                tend(k,j,i) = tend(k,j,i)                                      &
    193193                                       + ( - s_flux_t(j,i)                     &
    194                                            - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )   &
    195                                                  * ( s(k,j,i)-s(k-1,j,i) )     &
    196                                                  * ddzu(k)                     &
     194                                           - 0.5_wp * ( kh(k-1,j,i)+kh(k,j,i) )&
     195                                                    * ( s(k,j,i)-s(k-1,j,i) )  &
     196                                                    * ddzu(k)                  &
    197197                                         ) * ddzw(k)
    198198
     
    251251
    252252                   tend(k,j,i) = tend(k,j,i)                                   &
    253                                           + 0.5 * (                            &
     253                                          + 0.5_wp * (                         &
    254254                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
    255255                      - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
    256                                                   ) * ddx2                     &
    257                                           + 0.5 * (                            &
     256                                                     ) * ddx2                  &
     257                                          + 0.5_wp * (                         &
    258258                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
    259259                      - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
    260                                                   ) * ddy2
     260                                                     ) * ddy2
    261261                ENDIF
    262262             ENDDO
     
    266266             DO  k = 1, nzt
    267267                IF ( k > nzb_s_inner(j,i)  .AND.  k <= nzb_s_outer(j,i)  .AND. &
    268                      ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 ) )    &
     268                     ( wall_w_x(j,i) /= 0.0_wp  .OR.  wall_w_y(j,i) /= 0.0_wp ) )    &
    269269                THEN
    270270                   tend(k,j,i) = tend(k,j,i)                                   &
    271                                                 + ( fwxp(j,i) * 0.5 *          &
    272                         ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
    273                         + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                 &
    274                                                    -fwxm(j,i) * 0.5 *          &
     271                                                + ( fwxp(j,i) * 0.5_wp *       &
     272                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
     273                        + ( 1.0_wp - fwxp(j,i) ) * wall_s_flux(1)              &
     274                                                   -fwxm(j,i) * 0.5_wp *       &
    275275                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
    276                         + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                 &
     276                        + ( 1.0_wp - fwxm(j,i) ) * wall_s_flux(2)              &
    277277                                                  ) * ddx2                     &
    278                                                 + ( fwyp(j,i) * 0.5 *          &
    279                         ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
    280                         + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                 &
    281                                                    -fwym(j,i) * 0.5 *          &
     278                                                + ( fwyp(j,i) * 0.5_wp *       &
     279                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
     280                        + ( 1.0_wp - fwyp(j,i) ) * wall_s_flux(3)              &
     281                                                   -fwym(j,i) * 0.5_wp *       &
    282282                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
    283                         + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                 &
     283                        + ( 1.0_wp - fwym(j,i) ) * wall_s_flux(4)              &
    284284                                                  ) * ddy2
    285285                ENDIF
     
    293293                IF ( k >= nzb_diff_s_inner(j,i) )  THEN
    294294                   tend(k,j,i) = tend(k,j,i)                                   &
    295                                        + 0.5 * (                               &
     295                                       + 0.5_wp * (                            &
    296296            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
    297297          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
    298                                                ) * ddzw(k)
     298                                                  ) * ddzw(k)
    299299                ENDIF
    300300             ENDDO
     
    306306                IF ( use_surface_fluxes  .AND.  k == nzb_s_inner(j,i)+1 )  THEN
    307307                   tend(k,j,i) = tend(k,j,i)                                   &
    308                                           + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )  &
    309                                                   * ( s(k+1,j,i)-s(k,j,i) )    &
    310                                                   * ddzu(k+1)                  &
     308                                          + ( 0.5_wp * ( kh(k,j,i)+kh(k+1,j,i) )&
     309                                                     * ( s(k+1,j,i)-s(k,j,i) ) &
     310                                                     * ddzu(k+1)               &
    311311                                              + s_flux_b(j,i)                  &
    312312                                            ) * ddzw(k)
     
    319319                   tend(k,j,i) = tend(k,j,i)                                   &
    320320                                          + ( - s_flux_t(j,i)                  &
    321                                               - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )&
    322                                                     * ( s(k,j,i)-s(k-1,j,i) )  &
    323                                                     * ddzu(k)                  &
     321                                              - 0.5_wp * ( kh(k-1,j,i)+kh(k,j,i) )&
     322                                                       * ( s(k,j,i)-s(k-1,j,i) )  &
     323                                                       * ddzu(k)                  &
    324324                                            ) * ddzw(k)
    325325                ENDIF
     
    372372
    373373          tend(k,j,i) = tend(k,j,i)                                            &
    374                                           + 0.5 * (                            &
     374                                          + 0.5_wp * (                         &
    375375                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
    376376                      - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
    377                                                   ) * ddx2                     &
    378                                           + 0.5 * (                            &
     377                                                     ) * ddx2                  &
     378                                          + 0.5_wp * (                         &
    379379                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
    380380                      - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
    381                                                   ) * ddy2
     381                                                     ) * ddy2
    382382       ENDDO
    383383
    384384!
    385385!--    Apply prescribed horizontal wall heatflux where necessary
    386        IF ( ( wall_w_x(j,i) .NE. 0.0 ) .OR. ( wall_w_y(j,i) .NE. 0.0 ) )       &
     386       IF ( ( wall_w_x(j,i) .NE. 0.0_wp ) .OR. ( wall_w_y(j,i) .NE. 0.0_wp ) )       &
    387387       THEN
    388388          DO  k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
    389389
    390390             tend(k,j,i) = tend(k,j,i)                                         &
    391                                                 + ( fwxp(j,i) * 0.5 *          &
    392                         ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
    393                         + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                 &
    394                                                    -fwxm(j,i) * 0.5 *          &
     391                                                + ( fwxp(j,i) * 0.5_wp *       &
     392                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
     393                        + ( 1.0_wp - fwxp(j,i) ) * wall_s_flux(1)              &
     394                                                   -fwxm(j,i) * 0.5_wp *       &
    395395                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
    396                         + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                 &
     396                        + ( 1.0_wp - fwxm(j,i) ) * wall_s_flux(2)              &
    397397                                                  ) * ddx2                     &
    398                                                 + ( fwyp(j,i) * 0.5 *          &
    399                         ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
    400                         + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                 &
    401                                                    -fwym(j,i) * 0.5 *          &
     398                                                + ( fwyp(j,i) * 0.5_wp *       &
     399                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
     400                        + ( 1.0_wp - fwyp(j,i) ) * wall_s_flux(3)              &
     401                                                   -fwym(j,i) * 0.5_wp *       &
    402402                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
    403                         + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                 &
     403                        + ( 1.0_wp - fwym(j,i) ) * wall_s_flux(4)              &
    404404                                                  ) * ddy2
    405405          ENDDO
     
    413413
    414414          tend(k,j,i) = tend(k,j,i)                                            &
    415                                        + 0.5 * (                               &
     415                                       + 0.5_wp * (                            &
    416416            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
    417417          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
    418                                                ) * ddzw(k)
     418                                                  ) * ddzw(k)
    419419       ENDDO
    420420
     
    425425          k = nzb_s_inner(j,i)+1
    426426
    427           tend(k,j,i) = tend(k,j,i) + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )        &
    428                                             * ( s(k+1,j,i)-s(k,j,i) )          &
    429                                             * ddzu(k+1)                        &
     427          tend(k,j,i) = tend(k,j,i) + ( 0.5_wp * ( kh(k,j,i)+kh(k+1,j,i) )     &
     428                                               * ( s(k+1,j,i)-s(k,j,i) )       &
     429                                               * ddzu(k+1)                     &
    430430                                        + s_flux_b(j,i)                        &
    431431                                      ) * ddzw(k)
     
    440440
    441441          tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(j,i)                        &
    442                                       - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )        &
    443                                             * ( s(k,j,i)-s(k-1,j,i) )          &
    444                                             * ddzu(k)                          &
     442                                      - 0.5_wp * ( kh(k-1,j,i)+kh(k,j,i) )     &
     443                                               * ( s(k,j,i)-s(k-1,j,i) )       &
     444                                               * ddzu(k)                       &
    445445                                      ) * ddzw(k)
    446446
  • palm/trunk/SOURCE/diffusion_u.f90

    r1321 r1340  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants defined as wp-kind
    2323!
    2424! Former revisions:
     
    131131!
    132132!--             Interpolate eddy diffusivities on staggered gridpoints
    133                 kmyp = 0.25 *                                                  &
     133                kmyp = 0.25_wp *                                               &
    134134                       ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
    135                 kmym = 0.25 *                                                  &
     135                kmym = 0.25_wp *                                               &
    136136                       ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
    137137
    138138                tend(k,j,i) = tend(k,j,i)                                      &
    139                       & + 2.0 * (                                              &
     139                      & + 2.0_wp * (                                           &
    140140                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )    &
    141141                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )    &
    142                       &         ) * ddx2                                       &
     142                      &            ) * ddx2                                    &
    143143                      & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy         &
    144144                      &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx         &
     
    150150!
    151151!--          Wall functions at the north and south walls, respectively
    152              IF ( wall_u(j,i) /= 0.0 )  THEN
     152             IF ( wall_u(j,i) /= 0.0_wp )  THEN
    153153
    154154                DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
    155                    kmyp = 0.25 *                                               &
     155                   kmyp = 0.25_wp *                                            &
    156156                          ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
    157                    kmym = 0.25 *                                               &
     157                   kmym = 0.25_wp *                                            &
    158158                          ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
    159159
    160160                   tend(k,j,i) = tend(k,j,i)                                   &
    161                                  + 2.0 * (                                     &
     161                                 + 2.0_wp * (                                  &
    162162                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
    163163                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
    164                                          ) * ddx2                              &
     164                                            ) * ddx2                           &
    165165                                 + (   fyp(j,i) * (                            &
    166166                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
     
    182182!
    183183!--             Interpolate eddy diffusivities on staggered gridpoints
    184                 kmzp = 0.25 *                                                  &
     184                kmzp = 0.25_wp *                                               &
    185185                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
    186                 kmzm = 0.25 *                                                  &
     186                kmzm = 0.25_wp *                                               &
    187187                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
    188188
     
    211211!
    212212!--             Interpolate eddy diffusivities on staggered gridpoints
    213                 kmzp = 0.25 *                                                  &
     213                kmzp = 0.25_wp *                                               &
    214214                      ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
    215                 kmzm = 0.25 *                                                  &
     215                kmzm = 0.25_wp *                                               &
    216216                      ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
    217217
     
    231231!
    232232!--             Interpolate eddy diffusivities on staggered gridpoints
    233                 kmzp = 0.25 *                                                  &
     233                kmzp = 0.25_wp *                                               &
    234234                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
    235                 kmzm = 0.25 *                                                  &
     235                kmzm = 0.25_wp *                                               &
    236236                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
    237237
     
    303303!
    304304!--                Interpolate eddy diffusivities on staggered gridpoints
    305                    kmyp = 0.25 *                                               &
     305                   kmyp = 0.25_wp *                                            &
    306306                          ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
    307                    kmym = 0.25 *                                               &
     307                   kmym = 0.25_wp *                                            &
    308308                          ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
    309309
    310310                   tend(k,j,i) = tend(k,j,i)                                   &
    311                          & + 2.0 * (                                           &
     311                         & + 2.0_wp * (                                        &
    312312                         &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   ) &
    313313                         &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) ) &
    314                          &         ) * ddx2                                    &
     314                         &            ) * ddx2                                 &
    315315                         & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy      &
    316316                         &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx      &
     
    325325             DO  k = 1, nzt
    326326                IF( k > nzb_u_inner(j,i)  .AND.  k <= nzb_u_outer(j,i)  .AND.  &
    327                     wall_u(j,i) /= 0.0 )  THEN
    328 
    329                    kmyp = 0.25 *                                               &
     327                    wall_u(j,i) /= 0.0_wp )  THEN
     328
     329                   kmyp = 0.25_wp *                                            &
    330330                          ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
    331                    kmym = 0.25 *                                               &
     331                   kmym = 0.25_wp *                                            &
    332332                          ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
    333333
    334334                   tend(k,j,i) = tend(k,j,i)                                   &
    335                                  + 2.0 * (                                     &
     335                                 + 2.0_wp * (                                  &
    336336                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
    337337                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
    338                                          ) * ddx2                              &
     338                                            ) * ddx2                           &
    339339                                 + (   fyp(j,i) * (                            &
    340340                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
     
    357357!
    358358!--                Interpolate eddy diffusivities on staggered gridpoints
    359                    kmzp = 0.25 *                                               &
     359                   kmzp = 0.25_wp *                                            &
    360360                          ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
    361                    kmzm = 0.25 *                                               &
     361                   kmzm = 0.25_wp *                                            &
    362362                          ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
    363363
     
    394394!
    395395!--             Interpolate eddy diffusivities on staggered gridpoints
    396                 kmzp = 0.25 *                                                  &
     396                kmzp = 0.25_wp *                                               &
    397397                      ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
    398                 kmzm = 0.25 *                                                  &
     398                kmzm = 0.25_wp *                                               &
    399399                      ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
    400400
     
    422422!
    423423!--             Interpolate eddy diffusivities on staggered gridpoints
    424                 kmzp = 0.25 *                                                  &
     424                kmzp = 0.25_wp *                                               &
    425425                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
    426                 kmzm = 0.25 *                                                  &
     426                kmzm = 0.25_wp *                                               &
    427427                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
    428428
     
    478478!
    479479!--       Interpolate eddy diffusivities on staggered gridpoints
    480           kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
    481           kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
     480          kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
     481          kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
    482482
    483483          tend(k,j,i) = tend(k,j,i)                                            &
    484                       & + 2.0 * (                                              &
     484                      & + 2.0_wp * (                                           &
    485485                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )    &
    486486                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )    &
    487                       &         ) * ddx2                                       &
     487                      &            ) * ddx2                                    &
    488488                      & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy         &
    489489                      &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx         &
     
    495495!
    496496!--    Wall functions at the north and south walls, respectively
    497        IF ( wall_u(j,i) .NE. 0.0 )  THEN
     497       IF ( wall_u(j,i) .NE. 0.0_wp )  THEN
    498498
    499499!
     
    503503
    504504          DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
    505              kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
    506              kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
     505             kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
     506             kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
    507507
    508508             tend(k,j,i) = tend(k,j,i)                                         &
    509                                  + 2.0 * (                                     &
     509                                 + 2.0_wp * (                                  &
    510510                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
    511511                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
    512                                          ) * ddx2                              &
     512                                            ) * ddx2                           &
    513513                                 + (   fyp(j,i) * (                            &
    514514                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
     
    530530!
    531531!--       Interpolate eddy diffusivities on staggered gridpoints
    532           kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
    533           kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
     532          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
     533          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
    534534
    535535          tend(k,j,i) = tend(k,j,i)                                            &
     
    557557!
    558558!--       Interpolate eddy diffusivities on staggered gridpoints
    559           kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
    560           kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
     559          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
     560          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
    561561
    562562          tend(k,j,i) = tend(k,j,i)                                            &
     
    575575!
    576576!--       Interpolate eddy diffusivities on staggered gridpoints
    577           kmzp = 0.25 *                                                        &
     577          kmzp = 0.25_wp *                                                     &
    578578                 ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
    579           kmzm = 0.25 *                                                        &
     579          kmzm = 0.25_wp *                                                     &
    580580                 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
    581581
  • palm/trunk/SOURCE/diffusion_v.f90

    r1321 r1340  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants defined as wp-kind
    2323!
    2424! Former revisions:
     
    129129!
    130130!--             Interpolate eddy diffusivities on staggered gridpoints
    131                 kmxp = 0.25 * &
     131                kmxp = 0.25_wp * &
    132132                       ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
    133                 kmxm = 0.25 * &
     133                kmxm = 0.25_wp * &
    134134                       ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    135135
     
    140140                      &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
    141141                      &   ) * ddx                                              &
    142                       & + 2.0 * (                                              &
     142                      & + 2.0_wp * (                                           &
    143143                      &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )      &
    144144                      &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )      &
    145                       &         ) * ddy2
     145                      &            ) * ddy2
    146146             ENDDO
    147147
    148148!
    149149!--          Wall functions at the left and right walls, respectively
    150              IF ( wall_v(j,i) /= 0.0 )  THEN
     150             IF ( wall_v(j,i) /= 0.0_wp )  THEN
    151151
    152152                DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
    153                    kmxp = 0.25 *                                               &
     153                   kmxp = 0.25_wp *                                            &
    154154                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
    155                    kmxm = 0.25 *                                               &
     155                   kmxm = 0.25_wp *                                            &
    156156                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    157157                   
    158158                   tend(k,j,i) = tend(k,j,i)                                   &
    159                                  + 2.0 * (                                     &
     159                                 + 2.0_wp * (                                  &
    160160                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
    161161                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
    162                                          ) * ddy2                              &
     162                                            ) * ddy2                           &
    163163                                 + (   fxp(j,i) * (                            &
    164164                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
     
    180180!
    181181!--             Interpolate eddy diffusivities on staggered gridpoints
    182                 kmzp = 0.25 * &
     182                kmzp = 0.25_wp * &
    183183                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    184                 kmzm = 0.25 * &
     184                kmzm = 0.25_wp * &
    185185                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
    186186
     
    209209!
    210210!--             Interpolate eddy diffusivities on staggered gridpoints
    211                 kmzp = 0.25 *                                                  &
     211                kmzp = 0.25_wp *                                               &
    212212                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    213                 kmzm = 0.25 *                                                  &
     213                kmzm = 0.25_wp *                                               &
    214214                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
    215215
     
    229229!
    230230!--             Interpolate eddy diffusivities on staggered gridpoints
    231                 kmzp = 0.25 *                                                  &
     231                kmzp = 0.25_wp *                                               &
    232232                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    233                 kmzm = 0.25 *                                                  &
     233                kmzm = 0.25_wp *                                               &
    234234                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
    235235
     
    301301!
    302302!--                Interpolate eddy diffusivities on staggered gridpoints
    303                    kmxp = 0.25 *                                               &
     303                   kmxp = 0.25_wp *                                            &
    304304                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
    305                    kmxm = 0.25 *                                               &
     305                   kmxm = 0.25_wp *                                            &
    306306                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    307307
     
    312312                         &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy          &
    313313                         &   ) * ddx                                           &
    314                          & + 2.0 * (                                           &
     314                         & + 2.0_wp * (                                        &
    315315                         &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )   &
    316316                         &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )   &
    317                          &         ) * ddy2
     317                         &            ) * ddy2
    318318                ENDIF
    319319             ENDDO
     
    323323             DO  k = 1, nzt
    324324                IF( k > nzb_v_inner(j,i)  .AND.  k <= nzb_v_outer(j,i)  .AND.  &
    325                     wall_v(j,i) /= 0.0 )  THEN
    326 
    327                    kmxp = 0.25 *                                               &
     325                    wall_v(j,i) /= 0.0_wp )  THEN
     326
     327                   kmxp = 0.25_wp *                                            &
    328328                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
    329                    kmxm = 0.25 *                                               &
     329                   kmxm = 0.25_wp *                                            &
    330330                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    331331                   
    332332                   tend(k,j,i) = tend(k,j,i)                                   &
    333                                  + 2.0 * (                                     &
     333                                 + 2.0_wp * (                                  &
    334334                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
    335335                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
    336                                          ) * ddy2                              &
     336                                            ) * ddy2                           &
    337337                                 + (   fxp(j,i) * (                            &
    338338                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
     
    355355!
    356356!--                Interpolate eddy diffusivities on staggered gridpoints
    357                    kmzp = 0.25 *                                               &
     357                   kmzp = 0.25_wp *                                            &
    358358                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    359                    kmzm = 0.25 *                                               &
     359                   kmzm = 0.25_wp *                                            &
    360360                          ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
    361361
     
    392392!
    393393!--             Interpolate eddy diffusivities on staggered gridpoints
    394                 kmzp = 0.25 *                                                  &
     394                kmzp = 0.25_wp *                                               &
    395395                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    396                 kmzm = 0.25 *                                                  &
     396                kmzm = 0.25_wp *                                               &
    397397                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
    398398
     
    420420!
    421421!--             Interpolate eddy diffusivities on staggered gridpoints
    422                 kmzp = 0.25 *                                                  &
     422                kmzp = 0.25_wp *                                               &
    423423                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    424                 kmzm = 0.25 *                                                  &
     424                kmzm = 0.25_wp *                                               &
    425425                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
    426426
     
    476476!
    477477!--       Interpolate eddy diffusivities on staggered gridpoints
    478           kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
    479           kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
     478          kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
     479          kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    480480
    481481          tend(k,j,i) = tend(k,j,i)                                            &
     
    485485                      &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
    486486                      &   ) * ddx                                              &
    487                       & + 2.0 * (                                              &
     487                      & + 2.0_wp * (                                           &
    488488                      &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )      &
    489489                      &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )      &
    490                       &         ) * ddy2
     490                      &            ) * ddy2
    491491       ENDDO
    492492
    493493!
    494494!--    Wall functions at the left and right walls, respectively
    495        IF ( wall_v(j,i) /= 0.0 )  THEN
     495       IF ( wall_v(j,i) /= 0.0_wp )  THEN
    496496
    497497!
     
    501501
    502502          DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
    503              kmxp = 0.25 *                                                     &
     503             kmxp = 0.25_wp *                                                  &
    504504                    ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
    505              kmxm = 0.25 *                                                     &
     505             kmxm = 0.25_wp *                                                  &
    506506                    ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    507507
    508508             tend(k,j,i) = tend(k,j,i)                                         &
    509                                  + 2.0 * (                                     &
     509                                 + 2.0_wp * (                                  &
    510510                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
    511511                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
    512                                          ) * ddy2                              &
     512                                            ) * ddy2                           &
    513513                                 + (   fxp(j,i) * (                            &
    514514                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
     
    530530!
    531531!--       Interpolate eddy diffusivities on staggered gridpoints
    532           kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    533           kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
     532          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
     533          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
    534534
    535535          tend(k,j,i) = tend(k,j,i)                                            &
     
    557557!
    558558!--       Interpolate eddy diffusivities on staggered gridpoints
    559           kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    560           kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
     559          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
     560          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
    561561
    562562          tend(k,j,i) = tend(k,j,i)                                            &
     
    575575!
    576576!--       Interpolate eddy diffusivities on staggered gridpoints
    577           kmzp = 0.25 * &
     577          kmzp = 0.25_wp * &
    578578                 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    579           kmzm = 0.25 * &
     579          kmzm = 0.25_wp * &
    580580                 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
    581581
  • palm/trunk/SOURCE/diffusion_w.f90

    r1323 r1340  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants defined as wp-kind
    2323!
    2424! Former revisions:
     
    158158                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
    159159                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
    160                       &         ) * ddzu(k+1)
     160                      &            ) * ddzu(k+1)
    161161             ENDDO
    162162
    163163!
    164164!--          Wall functions at all vertical walls, where necessary
    165              IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
     165             IF ( wall_w_x(j,i) /= 0.0_wp  .OR.  wall_w_y(j,i) /= 0.0_wp )  THEN
    166166
    167167                DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
     
    198198                                     + wall_w_y(j,i) * wsvs(k,j,i)             &
    199199                                   ) * ddy                                     &
    200                                  + 2.0 * (                                     &
     200                                 + 2.0_wp * (                                  &
    201201                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
    202202                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
    203                                          ) * ddzu(k+1)
     203                                            ) * ddzu(k+1)
    204204                ENDDO
    205205             ENDIF
     
    288288                         &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
    289289                         & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
    290                          &         ) * ddzu(k+1)
     290                         &            ) * ddzu(k+1)
    291291                ENDIF
    292292             ENDDO
     
    297297
    298298                IF ( k > nzb_w_inner(j,i)  .AND.  k <= nzb_w_outer(j,i)  .AND. &
    299                      wall_w_x(j,i) /= 0.0  .AND.  wall_w_y(j,i) /= 0.0 )  THEN
     299                     wall_w_x(j,i) /= 0.0_wp  .AND.  wall_w_y(j,i) /= 0.0_wp )  THEN
    300300!
    301301!--                Interpolate eddy diffusivities on staggered gridpoints
     
    330330                                     + wall_w_y(j,i) * wsvs(k,j,i)             &
    331331                                   ) * ddy                                     &
    332                                  + 2.0 * (                                     &
     332                                 + 2.0_wp * (                                  &
    333333                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
    334334                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
    335                                          ) * ddzu(k+1)
     335                                            ) * ddzu(k+1)
    336336                ENDIF
    337337             ENDDO
     
    400400                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
    401401                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
    402                       &         ) * ddzu(k+1)
     402                      &            ) * ddzu(k+1)
    403403       ENDDO
    404404
    405405!
    406406!--    Wall functions at all vertical walls, where necessary
    407        IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
     407       IF ( wall_w_x(j,i) /= 0.0_wp  .OR.  wall_w_y(j,i) /= 0.0_wp )  THEN
    408408
    409409!
    410410!--       Calculate the horizontal momentum fluxes w'u' and/or w'v'
    411           IF ( wall_w_x(j,i) /= 0.0 )  THEN
     411          IF ( wall_w_x(j,i) /= 0.0_wp )  THEN
    412412             CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i),     &
    413413                               wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
    414414          ELSE
    415              wsus = 0.0
     415             wsus = 0.0_wp
    416416          ENDIF
    417417
    418           IF ( wall_w_y(j,i) /= 0.0 )  THEN
     418          IF ( wall_w_y(j,i) /= 0.0_wp )  THEN
    419419             CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i),     &
    420420                               wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
    421421          ELSE
    422              wsvs = 0.0
     422             wsvs = 0.0_wp
    423423          ENDIF
    424424
     
    452452                                     + wall_w_y(j,i) * wsvs(k)                 &
    453453                                   ) * ddy                                     &
    454                                  + 2.0 * (                                     &
     454                                 + 2.0_wp * (                                  &
    455455                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
    456456                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
    457                                          ) * ddzu(k+1)
     457                                            ) * ddzu(k+1)
    458458          ENDDO
    459459       ENDIF
  • palm/trunk/SOURCE/diffusivities.f90

    r1323 r1340  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants defined as wp-kind
    2323!
    2424! Former revisions:
     
    9999!
    100100!-- Initialization for calculation of the mixing length profile
    101     sums_l_l = 0.0
     101    sums_l_l = 0.0_wp
    102102
    103103!
     
    113113!
    114114!-- Introduce an optional minimum tke
    115     IF ( e_min > 0.0 )  THEN
     115    IF ( e_min > 0.0_wp )  THEN
    116116       !$OMP DO
    117117       !$acc loop
     
    142142                dvar_dz = atmos_ocean_sign * &  ! inverse effect of pt/rho gradient
    143143                          ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
    144                 IF ( dvar_dz > 0.0 ) THEN
     144                IF ( dvar_dz > 0.0_wp ) THEN
    145145                   IF ( use_single_reference_value )  THEN
    146146                      l_stable = 0.76_wp * sqrt_e /                            &
    147                                  SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
     147                                    SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
    148148                   ELSE
    149149                      l_stable = 0.76_wp * sqrt_e /                            &
    150                                  SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
     150                                    SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
    151151                   ENDIF
    152152                ELSE
  • palm/trunk/SOURCE/init_3d_model.f90

    r1323 r1340  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! REAL constants defined as wp-kind
    2323!
    2424! Former revisions:
     
    441441       IF ( humidity_remote )  THEN
    442442          ALLOCATE( qswst_remote(nysg:nyng,nxlg:nxrg))
    443           qswst_remote = 0.0
     443          qswst_remote = 0.0_wp
    444444       ENDIF
    445445    ENDIF
     
    453453    ENDIF
    454454
    455     IF ( dt_dosp /= 9999999.9 )  THEN
     455    IF ( dt_dosp /= 9999999.9_wp )  THEN
    456456       ALLOCATE( spectrum_x( 1:nx/2, 1:10, 1:10 ), &
    457457                 spectrum_y( 1:ny/2, 1:10, 1:10 ) )
    458        spectrum_x = 0.0
    459        spectrum_y = 0.0
     458       spectrum_x = 0.0_wp
     459       spectrum_y = 0.0_wp
    460460    ENDIF
    461461
     
    463463!-- 1D-array for large scale subsidence velocity
    464464    ALLOCATE ( w_subs(nzb:nzt+1) )
    465     w_subs = 0.0
     465    w_subs = 0.0_wp
    466466
    467467
     
    480480       ENDIF
    481481
    482        IF ( cthf /= 0.0 ) THEN
     482       IF ( cthf /= 0.0_wp ) THEN
    483483          ALLOCATE ( lai(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    484484                     canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     
    491491    IF ( topography /= 'flat' )  THEN
    492492       ALLOCATE( rif_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg,1:4) )
    493        rif_wall = 0.0
     493       rif_wall = 0.0_wp
    494494    ENDIF
    495495
     
    578578    ALLOCATE( weight_substep(1:intermediate_timestep_count_max), &
    579579              weight_pres(1:intermediate_timestep_count_max) )
    580     weight_substep = 1.0
    581     weight_pres    = 1.0
     580    weight_substep = 1.0_wp
     581    weight_pres    = 1.0_wp
    582582    intermediate_timestep_count = 1  ! needed when simulated_time = 0.0
    583583       
     
    617617                DO  i = nxlg, nxrg
    618618                   DO  j = nysg, nyng
    619                       qr(:,j,i) = 0.0
    620                       nr(:,j,i) = 0.0
     619                      qr(:,j,i) = 0.0_wp
     620                      nr(:,j,i) = 0.0_wp
    621621                   ENDDO
    622622                ENDDO
     
    640640             IF ( prandtl_layer )  THEN
    641641                rif  = rif1d(nzb+1)
    642                 ts   = 0.0  ! could actually be computed more accurately in the
    643                             ! 1D model. Update when opportunity arises.
     642                ts   = 0.0_wp  ! could actually be computed more accurately in the
     643                               ! 1D model. Update when opportunity arises.
    644644                us   = us1d
    645645                usws = usws1d
    646646                vsws = vsws1d
    647647             ELSE
    648                 ts   = 0.0  ! must be set, because used in
    649                 rif  = 0.0  ! flowste
    650                 us   = 0.0
    651                 usws = 0.0
    652                 vsws = 0.0
     648                ts   = 0.0_wp  ! must be set, because used in
     649                rif  = 0.0_wp  ! flowste
     650                us   = 0.0_wp
     651                usws = 0.0_wp
     652                vsws = 0.0_wp
    653653             ENDIF
    654654
    655655          ELSE
    656              e    = 0.0  ! must be set, because used in
    657              rif  = 0.0  ! flowste
    658              ts   = 0.0
    659              us   = 0.0
    660              usws = 0.0
    661              vsws = 0.0
     656             e    = 0.0_wp  ! must be set, because used in
     657             rif  = 0.0_wp  ! flowste
     658             ts   = 0.0_wp
     659             us   = 0.0_wp
     660             usws = 0.0_wp
     661             vsws = 0.0_wp
    662662          ENDIF
    663663          uswst = top_momentumflux_u
     
    669669!--       Update when opportunity arises!
    670670          IF ( humidity  .OR.  passive_scalar )  THEN
    671              qs = 0.0
     671             qs = 0.0_wp
    672672             IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
    673673                  precipitation )  THEN
    674                 qrs = 0.0
    675                 nrs = 0.0
     674                qrs = 0.0_wp
     675                nrs = 0.0_wp
    676676             ENDIF
    677677          ENDIF
     
    682682             DO  i = nxl-1, nxr+1
    683683                DO  j = nys-1, nyn+1
    684                    u(nzb:nzb_u_inner(j,i),j,i) = 0.0
    685                    v(nzb:nzb_v_inner(j,i),j,i) = 0.0
     684                   u(nzb:nzb_u_inner(j,i),j,i) = 0.0_wp
     685                   v(nzb:nzb_v_inner(j,i),j,i) = 0.0_wp
    686686                ENDDO
    687687             ENDDO
     
    745745          DO  i = nxlg, nxrg
    746746             DO  j = nysg, nyng
    747                 u(nzb:nzb_u_inner(j,i)+1,j,i) = 0.0
    748                 v(nzb:nzb_v_inner(j,i)+1,j,i) = 0.0
     747                u(nzb:nzb_u_inner(j,i)+1,j,i) = 0.0_wp
     748                v(nzb:nzb_v_inner(j,i)+1,j,i) = 0.0_wp
    749749             ENDDO
    750750          ENDDO
     
    764764                   DO  i = nxlg, nxrg
    765765                      DO  j = nysg, nyng
    766                          qr(:,j,i) = 0.0
    767                          nr(:,j,i) = 0.0
     766                         qr(:,j,i) = 0.0_wp
     767                         nr(:,j,i) = 0.0_wp
    768768                      ENDDO
    769769                   ENDDO
     
    784784             km   = km_constant
    785785             kh   = km / prandtl_number
    786              e    = 0.0
    787           ELSEIF ( e_init > 0.0 )  THEN
     786             e    = 0.0_wp
     787          ELSEIF ( e_init > 0.0_wp )  THEN
    788788             DO  k = nzb+1, nzt
    789                 km(k,:,:) = 0.1 * l_grid(k) * SQRT( e_init )
     789                km(k,:,:) = 0.1_wp * l_grid(k) * SQRT( e_init )
    790790             ENDDO
    791791             km(nzb,:,:)   = km(nzb+1,:,:)
     
    795795          ELSE
    796796             IF ( .NOT. ocean )  THEN
    797                 kh   = 0.01   ! there must exist an initial diffusion, because
    798                 km   = 0.01   ! otherwise no TKE would be produced by the
     797                kh   = 0.01_wp   ! there must exist an initial diffusion, because
     798                km   = 0.01_wp   ! otherwise no TKE would be produced by the
    799799                              ! production terms, as long as not yet
    800800                              ! e = (u*/cm)**2 at k=nzb+1
    801801             ELSE
    802                 kh   = 0.00001
    803                 km   = 0.00001
    804              ENDIF
    805              e    = 0.0
    806           ENDIF
    807           rif   = 0.0
    808           ts    = 0.0
    809           us    = 0.0
    810           usws  = 0.0
     802                kh   = 0.00001_wp
     803                km   = 0.00001_wp
     804             ENDIF
     805             e    = 0.0_wp
     806          ENDIF
     807          rif   = 0.0_wp
     808          ts    = 0.0_wp
     809          us    = 0.0_wp
     810          usws  = 0.0_wp
    811811          uswst = top_momentumflux_u
    812           vsws  = 0.0
     812          vsws  = 0.0_wp
    813813          vswst = top_momentumflux_v
    814           IF ( humidity  .OR.  passive_scalar ) qs = 0.0
     814          IF ( humidity  .OR.  passive_scalar ) qs = 0.0_wp
    815815
    816816!
     
    829829!--    Bottom boundary
    830830       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2  )  THEN
    831           u(nzb,:,:) = 0.0
    832           v(nzb,:,:) = 0.0
     831          u(nzb,:,:) = 0.0_wp
     832          v(nzb,:,:) = 0.0_wp
    833833       ENDIF
    834834
     
    836836!--    Apply channel flow boundary condition
    837837       IF ( TRIM( bc_uv_t ) == 'dirichlet_0' )  THEN
    838           u(nzt+1,:,:) = 0.0
    839           v(nzt+1,:,:) = 0.0
     838          u(nzt+1,:,:) = 0.0_wp
     839          v(nzt+1,:,:) = 0.0_wp
    840840       ENDIF
    841841
    842842!
    843843!--    Calculate virtual potential temperature
    844        IF ( humidity ) vpt = pt * ( 1.0 + 0.61 * q )
     844       IF ( humidity ) vpt = pt * ( 1.0_wp + 0.61_wp * q )
    845845
    846846!
     
    849849       hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
    850850       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2)  THEN
    851           hom(nzb,1,5,:) = 0.0   
    852           hom(nzb,1,6,:) = 0.0 
     851          hom(nzb,1,5,:) = 0.0_wp
     852          hom(nzb,1,6,:) = 0.0_wp
    853853       ENDIF
    854854       hom(:,1,7,:)  = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
     
    919919             IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
    920920                  precipitation )  THEN
    921                 qrsws = 0.0
    922                 nrsws = 0.0
     921                qrsws = 0.0_wp
     922                nrsws = 0.0_wp
    923923             ENDIF
    924924             IF ( constant_waterflux )  THEN
     
    954954
    955955             IF ( humidity  .OR.  passive_scalar )  THEN
    956                 qswst = 0.0
     956                qswst = 0.0_wp
    957957                IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
    958958                     precipitation ) THEN
    959                    nrswst = 0.0
    960                    qrswst = 0.0
     959                   nrswst = 0.0_wp
     960                   qrswst = 0.0_wp
    961961                ENDIF
    962962             ENDIF
     
    971971!--       Initialization in case of a coupled model run
    972972          IF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
    973              tswst = 0.0
     973             tswst = 0.0_wp
    974974          ENDIF
    975975
     
    990990!--          to be zero before the first time step. It approaches its correct
    991991!--          value in the course of the first few time steps.
    992              shf   = 0.0
     992             shf   = 0.0_wp
    993993          ENDIF
    994994
    995995          IF ( humidity  .OR.  passive_scalar )  THEN
    996              IF ( .NOT. constant_waterflux )  qsws   = 0.0
     996             IF ( .NOT. constant_waterflux )  qsws   = 0.0_wp
    997997             IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
    998998                  precipitation )  THEN
    999                 qrsws = 0.0
    1000                 nrsws = 0.0
     999                qrsws = 0.0_wp
     1000                nrsws = 0.0_wp
    10011001             ENDIF
    10021002          ENDIF
     
    10231023!
    10241024!--    For the moment, vertical velocity is zero
    1025        w = 0.0
     1025       w = 0.0_wp
    10261026
    10271027!
    10281028!--    Initialize array sums (must be defined in first call of pres)
    1029        sums = 0.0
     1029       sums = 0.0_wp
    10301030
    10311031!
    10321032!--    In case of iterative solvers, p must get an initial value
    1033        IF ( psolver == 'multigrid'  .OR.  psolver == 'sor' )  p = 0.0
     1033       IF ( psolver == 'multigrid'  .OR.  psolver == 'sor' )  p = 0.0_wp
    10341034
    10351035!
     
    10371037!--    are zero at beginning of the simulation
    10381038       IF ( cloud_physics )  THEN
    1039           ql = 0.0
    1040           IF ( precipitation )  precipitation_amount = 0.0
     1039          ql = 0.0_wp
     1040          IF ( precipitation )  precipitation_amount = 0.0_wp
    10411041          IF ( icloud_scheme == 0 )  THEN
    1042              qc = 0.0
     1042             qc = 0.0_wp
    10431043             nc_1d = nc_const
    10441044          ENDIF
     
    10581058!
    10591059!--    If required, change the surface temperature at the start of the 3D run
    1060        IF ( pt_surface_initial_change /= 0.0 )  THEN
     1060       IF ( pt_surface_initial_change /= 0.0_wp )  THEN
    10611061          pt(nzb,:,:) = pt(nzb,:,:) + pt_surface_initial_change
    10621062       ENDIF
     
    10661066!--    run
    10671067       IF ( ( humidity .OR. passive_scalar ) .AND. &
    1068             q_surface_initial_change /= 0.0 )  THEN
     1068            q_surface_initial_change /= 0.0_wp )  THEN
    10691069          q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
    10701070       ENDIF
     
    10751075!
    10761076!--    Initialize old and new time levels.
    1077        te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0
     1077       te_m = 0.0_wp; tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
    10781078       e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
    10791079
    10801080       IF ( humidity  .OR.  passive_scalar )  THEN
    1081           tq_m = 0.0
     1081          tq_m = 0.0_wp
    10821082          q_p = q
    10831083          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
    10841084               precipitation )  THEN
    1085              tqr_m = 0.0
     1085             tqr_m = 0.0_wp
    10861086             qr_p = qr
    1087              tnr_m = 0.0
     1087             tnr_m = 0.0_wp
    10881088             nr_p = nr
    10891089          ENDIF
     
    10911091
    10921092       IF ( ocean )  THEN
    1093           tsa_m = 0.0
     1093          tsa_m = 0.0_wp
    10941094          sa_p  = sa
    10951095       ENDIF
     
    11711171                   u(k,j,nxlg:-1)  = mean_inflow_profiles(k,1)
    11721172                   v(k,j,nxlg:-1)  = mean_inflow_profiles(k,2)
    1173                    w(k,j,nxlg:-1)  = 0.0
     1173                   w(k,j,nxlg:-1)  = 0.0_wp
    11741174                   pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)
    11751175                   e(k,j,nxlg:-1)  = mean_inflow_profiles(k,5)
     
    11831183!--       vertically because otherwise the turbulent inflow layer will grow
    11841184!--       in time.
    1185           IF ( inflow_damping_height == 9999999.9 )  THEN
     1185          IF ( inflow_damping_height == 9999999.9_wp )  THEN
    11861186!
    11871187!--          Default: use the inversion height calculated by the prerun; if
    11881188!--          this is zero, inflow_damping_height must be explicitly
    11891189!--          specified.
    1190              IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0 )  THEN
     1190             IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0_wp )  THEN
    11911191                inflow_damping_height = hom_sum(nzb+6,pr_palm,0)
    11921192             ELSE
     
    11991199          ENDIF
    12001200
    1201           IF ( inflow_damping_width == 9999999.9 )  THEN
     1201          IF ( inflow_damping_width == 9999999.9_wp )  THEN
    12021202!
    12031203!--          Default for the transition range: one tenth of the undamped
    12041204!--          layer
    1205              inflow_damping_width = 0.1 * inflow_damping_height
     1205             inflow_damping_width = 0.1_wp * inflow_damping_height
    12061206
    12071207          ENDIF
     
    12121212
    12131213             IF ( zu(k) <= inflow_damping_height )  THEN
    1214                 inflow_damping_factor(k) = 1.0
     1214                inflow_damping_factor(k) = 1.0_wp
    12151215             ELSEIF ( zu(k) <= ( inflow_damping_height + inflow_damping_width ) )  THEN
    1216                 inflow_damping_factor(k) = 1.0 -                               &
     1216                inflow_damping_factor(k) = 1.0_wp -                            &
    12171217                                           ( zu(k) - inflow_damping_height ) / &
    12181218                                           inflow_damping_width
    12191219             ELSE
    1220                 inflow_damping_factor(k) = 0.0
     1220                inflow_damping_factor(k) = 0.0_wp
    12211221             ENDIF
    12221222
     
    12351235          DO  i = nxlg, nxrg
    12361236             DO  j = nysg, nyng
    1237                 u  (nzb:nzb_u_inner(j,i),j,i)   = 0.0
    1238                 v  (nzb:nzb_v_inner(j,i),j,i)   = 0.0
    1239                 w  (nzb:nzb_w_inner(j,i),j,i)   = 0.0
    1240                 e  (nzb:nzb_w_inner(j,i),j,i)   = 0.0
    1241                 tu_m(nzb:nzb_u_inner(j,i),j,i)  = 0.0
    1242                 tv_m(nzb:nzb_v_inner(j,i),j,i)  = 0.0
    1243                 tw_m(nzb:nzb_w_inner(j,i),j,i)  = 0.0
    1244                 te_m(nzb:nzb_w_inner(j,i),j,i)  = 0.0
    1245                 tpt_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
     1237                u  (nzb:nzb_u_inner(j,i),j,i)   = 0.0_wp
     1238                v  (nzb:nzb_v_inner(j,i),j,i)   = 0.0_wp
     1239                w  (nzb:nzb_w_inner(j,i),j,i)   = 0.0_wp
     1240                e  (nzb:nzb_w_inner(j,i),j,i)   = 0.0_wp
     1241                tu_m(nzb:nzb_u_inner(j,i),j,i)  = 0.0_wp
     1242                tv_m(nzb:nzb_v_inner(j,i),j,i)  = 0.0_wp
     1243                tw_m(nzb:nzb_w_inner(j,i),j,i)  = 0.0_wp
     1244                te_m(nzb:nzb_w_inner(j,i),j,i)  = 0.0_wp
     1245                tpt_m(nzb:nzb_w_inner(j,i),j,i) = 0.0_wp
    12461246             ENDDO
    12471247          ENDDO
     
    12721272!--    have to be predefined here because they are used (but multiplied with 0)
    12731273!--    there before they are set.
    1274        te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0
     1274       te_m = 0.0_wp; tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
    12751275       IF ( humidity  .OR.  passive_scalar )  THEN
    1276           tq_m = 0.0
     1276          tq_m = 0.0_wp
    12771277          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
    12781278               precipitation )  THEN
    1279              tqr_m = 0.0
    1280              tnr_m = 0.0
    1281           ENDIF
    1282        ENDIF
    1283        IF ( ocean )  tsa_m = 0.0
     1279             tqr_m = 0.0_wp
     1280             tnr_m = 0.0_wp
     1281          ENDIF
     1282       ENDIF
     1283       IF ( ocean )  tsa_m = 0.0_wp
    12841284
    12851285    ELSE
     
    13231323       IF ( use_prescribed_profile_data )  THEN
    13241324
    1325           volume_flow_initial_l = 0.0
    1326           volume_flow_area_l    = 0.0
     1325          volume_flow_initial_l = 0.0_wp
     1326          volume_flow_area_l    = 0.0_wp
    13271327
    13281328          IF ( nxr == nx )  THEN
     
    13591359       ELSEIF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
    13601360
    1361           volume_flow_initial_l = 0.0
    1362           volume_flow_area_l    = 0.0
     1361          volume_flow_initial_l = 0.0_wp
     1362          volume_flow_area_l    = 0.0_wp
    13631363
    13641364          IF ( nxr == nx )  THEN
     
    13951395       ELSEIF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
    13961396
    1397           volume_flow_initial_l = 0.0
    1398           volume_flow_area_l    = 0.0
     1397          volume_flow_initial_l = 0.0_wp
     1398          volume_flow_area_l    = 0.0_wp
    13991399
    14001400          IF ( nxr == nx )  THEN
     
    15131513          DO  j = nys, nyn
    15141514             DO  k = nzb, nzt+1
    1515                 IF ( lad_s(k,j,i) > 0.0 )  THEN
     1515                IF ( lad_s(k,j,i) > 0.0_wp )  THEN
    15161516                   lad_u(k,j,i)   = lad_s(k,j,i)
    15171517                   lad_u(k,j,i+1) = lad_s(k,j,i)
     
    15211521             ENDDO
    15221522             DO  k = nzb, nzt
    1523                 lad_w(k,j,i) = 0.5 * ( lad_s(k+1,j,i) + lad_s(k,j,i) )
     1523                lad_w(k,j,i) = 0.5_wp * ( lad_s(k+1,j,i) + lad_s(k,j,i) )
    15241524             ENDDO
    15251525          ENDDO
    15261526       ENDDO
    15271527
    1528        lad_w(pch_index,:,:) = 0.0
     1528       lad_w(pch_index,:,:) = 0.0_wp
    15291529       lad_w(nzt+1,:,:)     = lad_w(nzt,:,:)
    15301530
     
    15351535!
    15361536!--    Initialisation of the canopy heat source distribution
    1537        IF ( cthf /= 0.0 )  THEN
     1537       IF ( cthf /= 0.0_wp )  THEN
    15381538!
    15391539!--       Piecewise evaluation of the leaf area index by
    15401540!--       integration of the leaf area density
    1541           lai(:,:,:) = 0.0
     1541          lai(:,:,:) = 0.0_wp
    15421542          DO  i = nxlg, nxrg
    15431543             DO  j = nysg, nyng
    15441544                DO  k = pch_index-1, 0, -1
    15451545                   lai(k,j,i) = lai(k+1,j,i) +                   &
    1546                                 ( 0.5 * ( lad_w(k+1,j,i) +       &
     1546                                ( 0.5_wp * ( lad_w(k+1,j,i) +    &
    15471547                                          lad_s(k+1,j,i) ) *     &
    15481548                                  ( zw(k+1) - zu(k+1) ) )  +     &
    1549                                 ( 0.5 * ( lad_w(k,j,i)   +       &
     1549                                ( 0.5_wp * ( lad_w(k,j,i)   +    &
    15501550                                          lad_s(k+1,j,i) ) *     &
    15511551                                  ( zu(k+1) - zw(k) ) )
     
    15611561                DO  k = 0, pch_index
    15621562                   canopy_heat_flux(k,j,i) = cthf *                    &
    1563                                              exp( -0.6 * lai(k,j,i) )
     1563                                             exp( -0.6_wp * lai(k,j,i) )
    15641564                ENDDO
    15651565             ENDDO
     
    15771577!
    15781578!-- If required, initialize dvrp-software
    1579     IF ( dt_dvrp /= 9999999.9 )  CALL init_dvrp
     1579    IF ( dt_dvrp /= 9999999.9_wp )  CALL init_dvrp
    15801580
    15811581    IF ( ocean )  THEN
     
    16241624    ELSE                                     ! for Euler-method
    16251625
    1626        weight_substep(1) = 1.0     
    1627        weight_pres(1)    = 1.0                   
     1626       weight_substep(1) = 1.0_wp     
     1627       weight_pres(1)    = 1.0_wp                   
    16281628
    16291629    ENDIF
     
    16311631!
    16321632!-- Initialize Rayleigh damping factors
    1633     rdf    = 0.0
    1634     rdf_sc = 0.0
    1635     IF ( rayleigh_damping_factor /= 0.0 )  THEN
     1633    rdf    = 0.0_wp
     1634    rdf_sc = 0.0_wp
     1635    IF ( rayleigh_damping_factor /= 0.0_wp )  THEN
    16361636       IF ( .NOT. ocean )  THEN
    16371637          DO  k = nzb+1, nzt
    16381638             IF ( zu(k) >= rayleigh_damping_height )  THEN
    16391639                rdf(k) = rayleigh_damping_factor * &
    1640                       ( SIN( pi * 0.5 * ( zu(k) - rayleigh_damping_height )    &
    1641                                       / ( zu(nzt) - rayleigh_damping_height ) )&
     1640                      ( SIN( pi * 0.5_wp * ( zu(k) - rayleigh_damping_height ) &
     1641                                         / ( zu(nzt) - rayleigh_damping_height ) )&
    16421642                      )**2
    16431643             ENDIF
     
    16471647             IF ( zu(k) <= rayleigh_damping_height )  THEN
    16481648                rdf(k) = rayleigh_damping_factor * &
    1649                       ( SIN( pi * 0.5 * ( rayleigh_damping_height - zu(k) )    &
    1650                                       / ( rayleigh_damping_height - zu(nzb+1)))&
     1649                      ( SIN( pi * 0.5_wp * ( rayleigh_damping_height - zu(k) ) &
     1650                                         / ( rayleigh_damping_height - zu(nzb+1)))&
    16511651                      )**2
    16521652             ENDIF
     
    16591659!-- Initialize the starting level and the vertical smoothing factor used for
    16601660!-- the external pressure gradient
    1661     dp_smooth_factor = 1.0
     1661    dp_smooth_factor = 1.0_wp
    16621662    IF ( dp_external )  THEN
    16631663!
     
    16701670       ENDIF
    16711671       IF ( dp_smooth )  THEN
    1672           dp_smooth_factor(:dp_level_ind_b) = 0.0
     1672          dp_smooth_factor(:dp_level_ind_b) = 0.0_wp
    16731673          DO  k = dp_level_ind_b+1, nzt
    1674              dp_smooth_factor(k) = 0.5 * ( 1.0 + SIN( pi * &
    1675                   ( REAL( k - dp_level_ind_b, KIND=wp ) /  &
    1676                     REAL( nzt - dp_level_ind_b, KIND=wp ) - 0.5 ) ) )
     1674             dp_smooth_factor(k) = 0.5_wp * ( 1.0_wp + SIN( pi *              &
     1675                        ( REAL( k - dp_level_ind_b, KIND=wp ) /                &
     1676                          REAL( nzt - dp_level_ind_b, KIND=wp ) - 0.5_wp ) ) )
    16771677          ENDDO
    16781678       ENDIF
     
    16831683!-- non-cyclic lateral boundaries. The damping zone has the maximum value
    16841684!-- at the inflow boundary and decreases to zero at pt_damping_width.
    1685     ptdf_x = 0.0
    1686     ptdf_y = 0.0
     1685    ptdf_x = 0.0_wp
     1686    ptdf_y = 0.0_wp
    16871687    IF ( bc_lr_dirrad )  THEN
    16881688       DO  i = nxl, nxr
    16891689          IF ( ( i * dx ) < pt_damping_width )  THEN
    1690              ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5 *              &
    1691                          REAL( pt_damping_width - i * dx, KIND=wp ) / (                 &
    1692                          REAL( pt_damping_width, KIND=wp )            ) ) )**2
     1690             ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5_wp *              &
     1691                            REAL( pt_damping_width - i * dx, KIND=wp ) / (     &
     1692                            REAL( pt_damping_width, KIND=wp )            ) ) )**2
    16931693          ENDIF
    16941694       ENDDO
     
    16971697          IF ( ( i * dx ) > ( nx * dx - pt_damping_width ) )  THEN
    16981698             ptdf_x(i) = pt_damping_factor *                                   &
    1699                          SIN( pi * 0.5 *                                       &
    1700                               ( ( i - nx ) * dx + pt_damping_width ) /         &
    1701                               REAL( pt_damping_width, KIND=wp ) )**2
     1699                         SIN( pi * 0.5_wp *                                    &
     1700                                 ( ( i - nx ) * dx + pt_damping_width ) /      &
     1701                                 REAL( pt_damping_width, KIND=wp ) )**2
    17021702          ENDIF
    17031703       ENDDO
     
    17061706          IF ( ( j * dy ) > ( ny * dy - pt_damping_width ) )  THEN
    17071707             ptdf_y(j) = pt_damping_factor *                                   &
    1708                          SIN( pi * 0.5 *                                       &
    1709                               ( ( j - ny ) * dy + pt_damping_width ) /         &
    1710                               REAL( pt_damping_width, KIND=wp ) )**2
     1708                         SIN( pi * 0.5_wp *                                    &
     1709                                 ( ( j - ny ) * dy + pt_damping_width ) /      &
     1710                                 REAL( pt_damping_width, KIND=wp ) )**2
    17111711          ENDIF
    17121712       ENDDO
     
    17151715          IF ( ( j * dy ) < pt_damping_width )  THEN
    17161716             ptdf_y(j) = pt_damping_factor *                                   &
    1717                          SIN( pi * 0.5 *                                       &
    1718                              ( pt_damping_width - j * dy ) /                   &
    1719                              REAL( pt_damping_width, KIND=wp ) )**2
     1717                         SIN( pi * 0.5_wp *                                    &
     1718                                ( pt_damping_width - j * dy ) /                &
     1719                                REAL( pt_damping_width, KIND=wp ) )**2
    17201720          ENDIF
    17211721       ENDDO
     
    17271727!-- are called from flow_statistics (or - depending on the chosen model run -
    17281728!-- are never initialized)
    1729     sums_divnew_l      = 0.0
    1730     sums_divold_l      = 0.0
    1731     sums_l_l           = 0.0
    1732     sums_up_fraction_l = 0.0
    1733     sums_wsts_bc_l     = 0.0
     1729    sums_divnew_l      = 0.0_wp
     1730    sums_divold_l      = 0.0_wp
     1731    sums_l_l           = 0.0_wp
     1732    sums_up_fraction_l = 0.0_wp
     1733    sums_wsts_bc_l     = 0.0_wp
    17341734
    17351735!
     
    17371737!-- Ghost points are excluded because counting values at the ghost boundaries
    17381738!-- would bias the statistics
    1739     rmask = 1.0
    1740     rmask(:,nxlg:nxl-1,:) = 0.0;  rmask(:,nxr+1:nxrg,:) = 0.0
    1741     rmask(nysg:nys-1,:,:) = 0.0;  rmask(nyn+1:nyng,:,:) = 0.0
     1739    rmask = 1.0_wp
     1740    rmask(:,nxlg:nxl-1,:) = 0.0_wp;  rmask(:,nxr+1:nxrg,:) = 0.0_wp
     1741    rmask(nysg:nys-1,:,:) = 0.0_wp;  rmask(nyn+1:nyng,:,:) = 0.0_wp
    17421742
    17431743!
     
    17691769    ngp_2dh_l         = 0
    17701770    ngp_2dh           = 0
    1771     ngp_3d_inner_l    = 0.0
     1771    ngp_3d_inner_l    = 0.0_wp
    17721772    ngp_3d_inner      = 0
    17731773    ngp_3d            = 0
     
    17771777       DO  i = nxl, nxr
    17781778          DO  j = nys, nyn
    1779              IF ( rmask(j,i,sr) == 1.0 )  THEN
     1779             IF ( rmask(j,i,sr) == 1.0_wp )  THEN
    17801780!
    17811781!--             All xy-grid points
  • palm/trunk/SOURCE/package_parin.f90

    r1325 r1340  
    2020! Current revisions:
    2121! -----------------
     22! REAL constants defined as wp-kinds
    2223!
    2324! Former revisions:
     
    196197!-- Default setting of dt_dosp here (instead of check_parameters), because its
    197198!-- current value is needed in init_pegrid
    198     IF ( dt_dosp == 9999999.9 )  dt_dosp = dt_data_output
     199    IF ( dt_dosp == 9999999.9_wp )  dt_dosp = dt_data_output
    199200
    200201 30 CONTINUE
  • palm/trunk/SOURCE/plant_canopy_model.f90

    r1321 r1340  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants defined as wp-kind
    2323!
    2424! Former revisions:
     
    100100                                              v(k,j+1,i)      +  &
    101101                                              v(k,j+1,i-1) )     &
    102                                             / 4.0 )**2        +  &
     102                                            / 4.0_wp )**2     +  &
    103103                                          ( ( w(k-1,j,i-1)    +  &
    104104                                              w(k-1,j,i)      +  &
    105105                                              w(k,j,i-1)      +  &
    106106                                              w(k,j,i) )         &
    107                                             / 4.0 )**2 )      *  &
     107                                            / 4.0_wp )**2 )   *  &
    108108                                    u(k,j,i)
    109109                   ENDDO
     
    123123                                              u(k,j,i)        +  &
    124124                                              u(k,j,i+1) )       &
    125                                             / 4.0 )**2        +  &
     125                                            / 4.0_wp )**2     +  &
    126126                                              v(k,j,i)**2     +  &
    127127                                          ( ( w(k-1,j-1,i)    +  &
     
    129129                                              w(k,j-1,i)      +  &
    130130                                              w(k,j,i) )         &
    131                                             / 4.0 )**2 )      *  &
     131                                            / 4.0_wp )**2 )   *  &
    132132                                    v(k,j,i) 
    133133                   ENDDO
     
    147147                                              u(k+1,j,i)      +  &
    148148                                              u(k+1,j,i+1) )     &
    149                                             / 4.0 )**2        +  &
     149                                            / 4.0_wp )**2     +  &
    150150                                          ( ( v(k,j,i)        +  &
    151151                                              v(k,j+1,i)      +  &
    152152                                              v(k+1,j,i)      +  &
    153153                                              v(k+1,j+1,i) )     &
    154                                             / 4.0 )**2        +  &
     154                                            / 4.0_wp )**2     +  &
    155155                                              w(k,j,i)**2 )   *  &
    156156                                    w(k,j,i)
     
    183183                                    SQRT( ( ( u(k,j,i)        +       &
    184184                                              u(k,j,i+1) )            &
    185                                             / 2.0 )**2        +       &
     185                                            / 2.0_wp )**2     +       &
    186186                                          ( ( v(k,j,i)        +       &
    187187                                              v(k,j+1,i) )            &
    188                                             / 2.0 )**2        +       &
     188                                            / 2.0_wp )**2     +       &
    189189                                          ( ( w(k-1,j,i)      +       &
    190190                                              w(k,j,i) )              &
    191                                             / 2.0 )**2 )      *       &
     191                                            / 2.0_wp )**2 )   *       &
    192192                                    ( q(k,j,i) - sls(k,j,i) )
    193193                   ENDDO
     
    202202                   DO  k = nzb_s_inner(j,i)+1, pch_index
    203203                      tend(k,j,i) = tend(k,j,i) -                     &
    204                                     2.0 * cdc(k,j,i) * lad_s(k,j,i) * &
     204                                    2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * &
    205205                                    SQRT( ( ( u(k,j,i)              + &
    206206                                              u(k,j,i+1) )            &
    207                                             / 2.0 )**2              + &
     207                                            / 2.0_wp )**2           + &
    208208                                          ( ( v(k,j,i)              + &
    209209                                              v(k,j+1,i) )            &
    210                                             / 2.0 )**2              + &
     210                                            / 2.0_wp )**2           + &
    211211                                          ( ( w(k,j,i)              + &
    212212                                              w(k+1,j,i) )            &
    213                                             / 2.0 )**2 )            * &
     213                                            / 2.0_wp )**2 )         * &
    214214                                    e(k,j,i)
    215215                   ENDDO
     
    267267                                        v(k,j+1,i)  +        &
    268268                                        v(k,j+1,i-1) )       &
    269                                       / 4.0 )**2    +        &
     269                                      / 4.0_wp )**2 +        &
    270270                                    ( ( w(k-1,j,i-1) +       &
    271271                                        w(k-1,j,i)   +       &
    272272                                        w(k,j,i-1)   +       &
    273273                                        w(k,j,i) )           &
    274                                       / 4.0 )**2 ) *         &
     274                                      / 4.0_wp )**2 ) *      &
    275275                              u(k,j,i)
    276276          ENDDO
     
    286286                                        u(k,j,i)     +       &
    287287                                        u(k,j,i+1) )         &
    288                                       / 4.0 )**2     +       &
     288                                      / 4.0_wp )**2  +       &
    289289                                        v(k,j,i)**2  +       &
    290290                                    ( ( w(k-1,j-1,i) +       &
     
    292292                                        w(k,j-1,i)   +       &
    293293                                        w(k,j,i) )           &
    294                                       / 4.0 )**2 ) *         &
     294                                      / 4.0_wp )**2 ) *      &
    295295                              v(k,j,i)
    296296          ENDDO
     
    306306                                        u(k+1,j,i)  +        &
    307307                                        u(k+1,j,i+1) )       &
    308                                       / 4.0 )**2    +        &
     308                                      / 4.0_wp )**2 +        &
    309309                                    ( ( v(k,j,i)    +        &
    310310                                        v(k,j+1,i)  +        &
    311311                                        v(k+1,j,i)  +        &
    312312                                        v(k+1,j+1,i) )       &
    313                                       / 4.0 )**2    +        &
     313                                      / 4.0_wp )**2 +        &
    314314                                        w(k,j,i)**2 ) *      &
    315315                              w(k,j,i)
     
    336336                              SQRT( ( ( u(k,j,i)        +       &
    337337                                        u(k,j,i+1) )            &
    338                                       / 2.0 )**2        +       &
     338                                      / 2.0_wp )**2     +       &
    339339                                    ( ( v(k,j,i)        +       &
    340340                                        v(k,j+1,i) )            &
    341                                       / 2.0 )**2        +       &
     341                                      / 2.0_wp )**2     +       &
    342342                                    ( ( w(k-1,j,i)      +       &
    343343                                        w(k,j,i) )              &
    344                                       / 2.0 )**2 )      *       &
     344                                      / 2.0_wp )**2 )   *       &
    345345                              ( q(k,j,i) - sls(k,j,i) )
    346346             ENDDO   
     
    351351          DO  k = nzb_s_inner(j,i)+1, pch_index   
    352352             tend(k,j,i) = tend(k,j,i) -                        &
    353                               2.0 * cdc(k,j,i) * lad_s(k,j,i) * &
     353                              2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * &
    354354                              SQRT( ( ( u(k,j,i)           +    &
    355355                                        u(k,j,i+1) )            &
    356                                       / 2.0 )**2           +    & 
     356                                      / 2.0_wp )**2        +    & 
    357357                                    ( ( v(k,j,i)           +    &
    358358                                        v(k,j+1,i) )            &
    359                                       / 2.0 )**2           +    &
     359                                      / 2.0_wp )**2        +    &
    360360                                    ( ( w(k,j,i)           +    &
    361361                                        w(k+1,j,i) )            &
    362                                       / 2.0 )**2 )         *    &
     362                                      / 2.0_wp )**2 )      *    &
    363363                              e(k,j,i)
    364364          ENDDO
  • palm/trunk/SOURCE/prandtl_fluxes.f90

    r1321 r1340  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants defined as wp-kind
    2323!
    2424! Former revisions:
     
    105105       DO  i = nxlg, nxrg
    106106          DO  j = nysg, nyng
    107              ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30 )
     107             ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30_wp )
    108108!
    109109!--          ts must be limited, because otherwise overflow may occur in case of
    110110!--          us=0 when computing rif further below
    111              IF ( ts(j,i) < -1.05E5 )  ts(j,i) = -1.0E5
    112              IF ( ts(j,i) >   1.0E5 )  ts(j,i) =  1.0E5
     111             IF ( ts(j,i) < -1.05E5_wp )  ts(j,i) = -1.0E5_wp
     112             IF ( ts(j,i) >   1.0E5_wp )  ts(j,i) =  1.0E5_wp
    113113          ENDDO
    114114       ENDDO
     
    131131             z_p = zu(k+1) - zw(k)
    132132
    133              IF ( rif(j,i) >= 0.0 )  THEN
     133             IF ( rif(j,i) >= 0.0_wp )  THEN
    134134!
    135135!--             Stable stratification
    136136                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (           &
    137137                                  LOG( z_p / z0h(j,i) ) +                   &
    138                                   5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
     138                                  5.0_wp * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
    139139                                                                )
    140140             ELSE
    141141!
    142142!--             Unstable stratification
    143                 a = SQRT( 1.0 - 16.0 * rif(j,i) )
    144                 b = SQRT( 1.0 - 16.0 * rif(j,i) * z0h(j,i) / z_p )
     143                a = SQRT( 1.0_wp - 16.0_wp * rif(j,i) )
     144                b = SQRT( 1.0_wp - 16.0_wp * rif(j,i) * z0h(j,i) / z_p )
    145145
    146146                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) /  (          &
    147147                          LOG( z_p / z0h(j,i) ) -                           &
    148                           2.0 * LOG( ( 1.0 + a ) / ( 1.0 + b ) ) )
     148                          2.0_wp * LOG( ( 1.0_wp + a ) / ( 1.0_wp + b ) ) )
    149149             ENDIF
    150150
     
    163163             z_p = zu(k+1) - zw(k)
    164164             rif(j,i) = z_p * kappa * g * ts(j,i) / &
    165                         ( pt(k+1,j,i) * ( us(j,i)**2 + 1E-30 ) )
     165                        ( pt(k+1,j,i) * ( us(j,i)**2 + 1E-30_wp ) )
    166166!
    167167!--          Limit the value range of the Richardson numbers.
     
    182182             z_p = zu(k+1) - zw(k)
    183183             rif(j,i) = z_p * kappa * g *                            &
    184                         ( ts(j,i) + 0.61 * pt(k+1,j,i) * qs(j,i) ) / &
    185                         ( vpt(k+1,j,i) * ( us(j,i)**2 + 1E-30 ) )
     184                        ( ts(j,i) + 0.61_wp * pt(k+1,j,i) * qs(j,i) ) / &
     185                        ( vpt(k+1,j,i) * ( us(j,i)**2 + 1E-30_wp ) )
    186186!
    187187!--          Limit the value range of the Richardson numbers.
     
    209209!--       Compute the absolute value of the horizontal velocity
    210210!--       (relative to the surface)
    211           uv_total = SQRT( ( 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1)        &
    212                                    - u(k,j,i)   - u(k,j,i+1) ) )**2 + &
    213                            ( 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i)        &
    214                                    - v(k,j,i)   - v(k,j+1,i) ) )**2 )   
    215 
    216 
    217           IF ( rif(j,i) >= 0.0 )  THEN
     211          uv_total = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)        &
     212                                      - u(k,j,i)   - u(k,j,i+1) ) )**2 + &
     213                           ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)        &
     214                                      - v(k,j,i)   - v(k,j+1,i) ) )**2 )   
     215
     216
     217          IF ( rif(j,i) >= 0.0_wp )  THEN
    218218!
    219219!--          Stable stratification
    220220             us(j,i) = kappa * uv_total / (                                &
    221221                                  LOG( z_p / z0(j,i) ) +                   &
    222                                   5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
     222                                  5.0_wp * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
    223223                                          )
    224224          ELSE
    225225!
    226226!--          Unstable stratification
    227              a = SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
    228              b = SQRT( SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) ) )
     227             a = SQRT( SQRT( 1.0_wp - 16.0_wp * rif(j,i) ) )
     228             b = SQRT( SQRT( 1.0_wp - 16.0_wp * rif(j,i) / z_p * z0(j,i) ) )
    229229
    230230             us(j,i) = kappa * uv_total / (                                  &
    231231                       LOG( z_p / z0(j,i) ) -                                &
    232                        LOG( ( 1.0 + a )**2 * ( 1.0 + a**2 ) / (              &
    233                             ( 1.0 + b )**2 * ( 1.0 + b**2 )   ) ) +          &
    234                             2.0 * ( ATAN( a ) - ATAN( b ) )                  &
     232                       LOG( ( 1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / (              &
     233                            ( 1.0_wp + b )**2 * ( 1.0_wp + b**2 )   ) ) +          &
     234                               2.0_wp * ( ATAN( a ) - ATAN( b ) )                  &
    235235                                           )
    236236          ENDIF
     
    258258!
    259259!--       Compute Richardson-flux number for this point
    260           rifm = 0.5 * ( rif(j,i-1) + rif(j,i) )
    261           IF ( rifm >= 0.0 )  THEN
     260          rifm = 0.5_wp * ( rif(j,i-1) + rif(j,i) )
     261          IF ( rifm >= 0.0_wp )  THEN
    262262!
    263263!--          Stable stratification
    264264             usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) )/ (              &
    265265                                     LOG( z_p / z0(j,i) ) +               &
    266                                      5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
    267                                               )
     266                                     5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p &
     267                                                            )
    268268          ELSE
    269269!
    270270!--          Unstable stratification
    271              a = SQRT( SQRT( 1.0 - 16.0 * rifm ) )
    272              b = SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
     271             a = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm ) )
     272             b = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm / z_p * z0(j,i) ) )
    273273
    274274             usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) ) / (            &
    275275                         LOG( z_p / z0(j,i) ) -                           &
    276                          LOG( (1.0 + a )**2 * ( 1.0 + a**2 ) / (          &
    277                               (1.0 + b )**2 * ( 1.0 + b**2 )   ) ) +      &
    278                               2.0 * ( ATAN( a ) - ATAN( b ) )             &
     276                         LOG( (1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / (          &
     277                              (1.0_wp + b )**2 * ( 1.0_wp + b**2 )   ) ) +      &
     278                                 2.0_wp * ( ATAN( a ) - ATAN( b ) )             &
    279279                                                 )
    280280          ENDIF
    281           usws(j,i) = -usws(j,i) * 0.5 * ( us(j,i-1) + us(j,i) )
     281          usws(j,i) = -usws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )
    282282       ENDDO
    283283    ENDDO
     
    296296!
    297297!--       Compute Richardson-flux number for this point
    298           rifm = 0.5 * ( rif(j-1,i) + rif(j,i) )
    299           IF ( rifm >= 0.0 )  THEN
     298          rifm = 0.5_wp * ( rif(j-1,i) + rif(j,i) )
     299          IF ( rifm >= 0.0_wp )  THEN
    300300!
    301301!--          Stable stratification
    302302             vsws(j,i) = kappa * ( v(k+1,j,i) -  v(k,j,i) ) / (           &
    303303                                     LOG( z_p / z0(j,i) ) +               &
    304                                      5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
    305                                               )
     304                                     5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p &
     305                                                              )
    306306          ELSE
    307307!
    308308!--          Unstable stratification
    309              a = SQRT( SQRT( 1.0 - 16.0 * rifm ) )
    310              b = SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
     309             a = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm ) )
     310             b = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm / z_p * z0(j,i) ) )
    311311
    312312             vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) ) / (            &
    313313                         LOG( z_p / z0(j,i) ) -                           &
    314                          LOG( (1.0 + a )**2 * ( 1.0 + a**2 ) / (          &
    315                               (1.0 + b )**2 * ( 1.0 + b**2 )   ) ) +      &
    316                               2.0 * ( ATAN( a ) - ATAN( b ) )             &
     314                         LOG( (1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / (          &
     315                              (1.0_wp + b )**2 * ( 1.0_wp + b**2 )   ) ) +      &
     316                                 2.0_wp * ( ATAN( a ) - ATAN( b ) )             &
    317317                                                 )
    318318          ENDIF
    319           vsws(j,i) = -vsws(j,i) * 0.5 * ( us(j-1,i) + us(j,i) )
     319          vsws(j,i) = -vsws(j,i) * 0.5_wp * ( us(j-1,i) + us(j,i) )
    320320       ENDDO
    321321    ENDDO
     
    331331          DO  i = nxlg, nxrg
    332332             DO  j = nysg, nyng
    333                 qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30 )
     333                qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30_wp )
    334334             ENDDO
    335335          ENDDO
     
    355355!--             in case of precursor runs)
    356356                IF ( coupled_run )  THEN
    357                    e_q = 6.1 * &
    358                         EXP( 0.07 * ( MIN(pt(0,j,i),pt(1,j,i)) - 273.15 ) )
    359                    q(k,j,i) = 0.622 * e_q / ( surface_pressure - e_q )
     357                   e_q = 6.1_wp * &
     358                           EXP( 0.07_wp * ( MIN(pt(0,j,i),pt(1,j,i)) - 273.15_wp ) )
     359                   q(k,j,i) = 0.622_wp * e_q / ( surface_pressure - e_q )
    360360                ENDIF
    361                 IF ( rif(j,i) >= 0.0 )  THEN
     361                IF ( rif(j,i) >= 0.0_wp )  THEN
    362362!
    363363!--                Stable stratification
    364364                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (          &
    365365                                  LOG( z_p / z0h(j,i) ) +                   &
    366                                   5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
     366                                  5.0_wp * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
    367367                                                                 )
    368368                ELSE
    369369!
    370370!--                Unstable stratification
    371                    a = SQRT( 1.0 - 16.0 * rif(j,i) )
    372                    b = SQRT( 1.0 - 16.0 * rif(j,i) * z0h(j,i) / z_p )
     371                   a = SQRT( 1.0_wp - 16.0_wp * rif(j,i) )
     372                   b = SQRT( 1.0_wp - 16.0_wp * rif(j,i) * z0h(j,i) / z_p )
    373373 
    374374                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) /   (        &
    375375                             LOG( z_p / z0h(j,i) ) -                        &
    376                               2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
     376                              2.0_wp * LOG( (1.0_wp + a ) / ( 1.0_wp + b ) ) )
    377377                ENDIF
    378378
     
    429429          !$acc loop independent
    430430          DO  j = nysg, nyng
    431              e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1 )**2
     431             e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1_wp )**2
    432432!
    433433!--          As a test: cm = 0.4
    434 !            e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.4 )**2
     434!            e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.4_wp )**2
    435435             e(nzb_s_inner(j,i),j,i)   = e(nzb_s_inner(j,i)+1,j,i)
    436436          ENDDO
Note: See TracChangeset for help on using the changeset viewer.