Ignore:
Timestamp:
Mar 26, 2014 5:04:47 PM (8 years ago)
Author:
kanani
Message:

REAL constants defined as wp-kind

File:
1 edited

Legend:

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

    r1321 r1342  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants defined as wp-kind
    2323!
    2424! Former revisions:
     
    173173             DO  k = nzb_diff_s_outer(j,i), nzt
    174174
    175                 dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    176                 dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    177                                  u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    178                 dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    179                                  u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    180 
    181                 dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    182                                  v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    183                 dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    184                 dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    185                                  v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    186 
    187                 dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    188                                  w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    189                 dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    190                                  w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    191                 dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
    192 
    193                 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
     175                dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
     176                dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     177                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     178                dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     179                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     180
     181                dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     182                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     183                dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     184                dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     185                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     186
     187                dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     188                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     189                dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     190                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     191                dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
     192
     193                def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
    194194                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
    195                       dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
    196 
    197                 IF ( def < 0.0 )  def = 0.0
     195                      dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
     196
     197                IF ( def < 0.0_wp )  def = 0.0_wp
    198198
    199199                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
     
    210210             DO  j = nys, nyn
    211211
    212                 IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
     212                IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) &
    213213                THEN
    214214
    215215                   k = nzb_diff_s_inner(j,i) - 1
    216216                   dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
    217                    dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    218                                   u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
     217                   dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     218                                     u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
    219219                   dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
    220                    dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    221                                   v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
     220                   dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     221                                     v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
    222222                   dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
    223223
    224                    IF ( wall_e_y(j,i) /= 0.0 )  THEN
     224                   IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
    225225!
    226226!--                   Inconsistency removed: as the thermal stratification is
     
    236236                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
    237237                                          wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
    238                       km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25 * &
    239                                    0.5 * dy
    240                       IF ( km_neutral > 0.0 )  THEN
     238                      km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25_wp * &
     239                                   0.5_wp * dy
     240                      IF ( km_neutral > 0.0_wp )  THEN
    241241                         dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
    242242                         dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
    243243                      ELSE
    244                          dudy = 0.0
    245                          dwdy = 0.0
     244                         dudy = 0.0_wp
     245                         dwdy = 0.0_wp
    246246                      ENDIF
    247247                   ELSE
    248                       dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    249                                       u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    250                       dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    251                                       w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     248                      dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     249                                         u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     250                      dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     251                                         w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    252252                   ENDIF
    253253
    254                    IF ( wall_e_x(j,i) /= 0.0 )  THEN
     254                   IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
    255255!
    256256!--                   Inconsistency removed: as the thermal stratification is
     
    266266                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
    267267                                          wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
    268                       km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25 * &
    269                                    0.5 * dx
    270                       IF ( km_neutral > 0.0 )  THEN
     268                      km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25_wp * &
     269                                   0.5_wp * dx
     270                      IF ( km_neutral > 0.0_wp )  THEN
    271271                         dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
    272272                         dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
    273273                      ELSE
    274                          dvdx = 0.0
    275                          dwdx = 0.0
     274                         dvdx = 0.0_wp
     275                         dwdx = 0.0_wp
    276276                      ENDIF
    277277                   ELSE
    278                       dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    279                                       v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    280                       dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    281                                       w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     278                      dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     279                                         v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     280                      dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     281                                         w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    282282                   ENDIF
    283283
    284                    def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
     284                   def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
    285285                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
    286                          dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
    287 
    288                    IF ( def < 0.0 )  def = 0.0
     286                         dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
     287
     288                   IF ( def < 0.0_wp )  def = 0.0_wp
    289289
    290290                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
     
    301301
    302302                      dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
    303                       dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    304                                      u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    305                       dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    306                       dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    307                                      v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     303                      dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     304                                        u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     305                      dvdy =          ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     306                      dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     307                                        v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    308308                      dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
    309309
    310                       IF ( wall_e_y(j,i) /= 0.0 )  THEN
     310                      IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
    311311!
    312312!--                      Inconsistency removed: as the thermal stratification
     
    319319!--                            validation has been available
    320320                         km_neutral = kappa * ( usvs(k)**2 + &
    321                                                 wsvs(k)**2 )**0.25 * 0.5 * dy
    322                          IF ( km_neutral > 0.0 )  THEN
     321                                                wsvs(k)**2 )**0.25_wp * 0.5_wp * dy
     322                         IF ( km_neutral > 0.0_wp )  THEN
    323323                            dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
    324324                            dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
    325325                         ELSE
    326                             dudy = 0.0
    327                             dwdy = 0.0
     326                            dudy = 0.0_wp
     327                            dwdy = 0.0_wp
    328328                         ENDIF
    329329                      ELSE
    330                          dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    331                                          u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    332                          dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    333                                          w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     330                         dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     331                                            u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     332                         dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     333                                            w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    334334                      ENDIF
    335335
    336                       IF ( wall_e_x(j,i) /= 0.0 )  THEN
     336                      IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
    337337!
    338338!--                      Inconsistency removed: as the thermal stratification
     
    345345!--                            validation has been available
    346346                         km_neutral = kappa * ( vsus(k)**2 + &
    347                                                 wsus(k)**2 )**0.25 * 0.5 * dx
    348                          IF ( km_neutral > 0.0 )  THEN
     347                                                wsus(k)**2 )**0.25_wp * 0.5_wp * dx
     348                         IF ( km_neutral > 0.0_wp )  THEN
    349349                            dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
    350350                            dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
    351351                         ELSE
    352                             dvdx = 0.0
    353                             dwdx = 0.0
     352                            dvdx = 0.0_wp
     353                            dwdx = 0.0_wp
    354354                         ENDIF
    355355                      ELSE
    356                          dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    357                                          v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    358                          dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    359                                          w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     356                         dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     357                                            v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     358                         dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     359                                            w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    360360                      ENDIF
    361361
    362                       def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
     362                      def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
    363363                           dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 +  &
    364                            dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
    365 
    366                       IF ( def < 0.0 )  def = 0.0
     364                           dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
     365
     366                      IF ( def < 0.0_wp )  def = 0.0_wp
    367367
    368368                      tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
     
    379379             DO  j = nys, nyn
    380380
    381                 IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
     381                IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) &
    382382                THEN
    383383
    384384                   k = nzb_diff_s_outer(j,i)-1
    385385
    386                    dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    387                    dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    388                                     u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    389                    dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    390                                     u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    391 
    392                    dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    393                                     v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    394                    dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    395                    dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    396                                     v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    397 
    398                    dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    399                                     w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    400                    dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    401                                     w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    402                    dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
    403 
    404                    def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
     386                   dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
     387                   dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     388                                       u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     389                   dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     390                                       u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     391
     392                   dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     393                                       v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     394                   dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     395                   dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     396                                       v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     397
     398                   dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     399                                       w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     400                   dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     401                                       w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     402                   dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
     403
     404                   def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
    405405                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
    406                          dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
    407 
    408                    IF ( def < 0.0 )  def = 0.0
     406                         dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
     407
     408                   IF ( def < 0.0_wp )  def = 0.0_wp
    409409
    410410                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
     
    420420             DO  j = nys, nyn
    421421
    422                 IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) &
     422                IF ( ( wall_e_x(j,i) == 0.0_wp ) .AND. ( wall_e_y(j,i) == 0.0_wp ) ) &
    423423                THEN
    424424
    425425                   k = nzb_diff_s_inner(j,i)-1
    426426
    427                    dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    428                    dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     427                   dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
     428                   dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     429                                       u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     430                   dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     431                                       u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
     432
     433                   dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     434                                       v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     435                   dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     436                   dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     437                                       v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
     438
     439                   dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     440                                       w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     441                   dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     442                                       w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     443                   dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
     444
     445                   def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
     446                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
     447                         dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
     448
     449                   IF ( def < 0.0_wp )  def = 0.0_wp
     450
     451                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
     452
     453                ENDIF
     454
     455             ENDDO
     456
     457          ELSEIF ( use_surface_fluxes )  THEN
     458
     459             DO  j = nys, nyn
     460
     461                k = nzb_diff_s_outer(j,i)-1
     462
     463                dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
     464                dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    429465                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    430                    dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    431                                     u_0(j,i)   - u_0(j,i+1)  ) * dd2zu(k)
    432 
    433                    dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     466                dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     467                                   u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     468
     469                dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    434470                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    435                    dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    436                    dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    437                                     v_0(j,i)   - v_0(j+1,i)  ) * dd2zu(k)
    438 
    439                    dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     471                dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     472                dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     473                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     474
     475                dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    440476                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    441                    dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     477                dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    442478                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    443                    dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
    444 
    445                    def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
    446                          dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
    447                          dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
    448 
    449                    IF ( def < 0.0 )  def = 0.0
    450 
    451                    tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
    452 
    453                 ENDIF
    454 
    455              ENDDO
    456 
    457           ELSEIF ( use_surface_fluxes )  THEN
    458 
    459              DO  j = nys, nyn
    460 
    461                 k = nzb_diff_s_outer(j,i)-1
    462 
    463                 dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    464                 dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    465                                  u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    466                 dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    467                                  u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    468 
    469                 dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    470                                  v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    471                 dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    472                 dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    473                                  v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    474 
    475                 dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    476                                  w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    477                 dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    478                                  w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    479                 dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
    480 
    481                 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
     479                dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
     480
     481                def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
    482482                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
    483                       dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
    484 
    485                 IF ( def < 0.0 )  def = 0.0
     483                      dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
     484
     485                IF ( def < 0.0_wp )  def = 0.0_wp
    486486
    487487                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
     
    586586
    587587                      IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
    588                          k1 = 1.0 + 0.61 * q(k,j,i)
    589                          k2 = 0.61 * pt(k,j,i)
     588                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     589                         k2 = 0.61_wp * pt(k,j,i)
    590590                         tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
    591591                                         g / vpt(k,j,i) *                      &
     
    594594                                         ) * dd2zu(k)
    595595                      ELSE IF ( cloud_physics )  THEN
    596                          IF ( ql(k,j,i) == 0.0 )  THEN
    597                             k1 = 1.0 + 0.61 * q(k,j,i)
    598                             k2 = 0.61 * pt(k,j,i)
     596                         IF ( ql(k,j,i) == 0.0_wp )  THEN
     597                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     598                            k2 = 0.61_wp * pt(k,j,i)
    599599                         ELSE
    600600                            theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
    601601                            temp  = theta * t_d_pt(k)
    602                             k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
    603                                        ( q(k,j,i) - ql(k,j,i) ) *          &
    604                                  ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
    605                                  ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
     602                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
     603                                          ( q(k,j,i) - ql(k,j,i) ) *          &
     604                                 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
     605                                 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
    606606                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    607                             k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
     607                            k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
    608608                         ENDIF
    609609                         tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
     
    613613                                         ) * dd2zu(k)
    614614                      ELSE IF ( cloud_droplets )  THEN
    615                          k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
    616                          k2 = 0.61 * pt(k,j,i)
     615                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
     616                         k2 = 0.61_wp * pt(k,j,i)
    617617                         tend(k,j,i) = tend(k,j,i) -                          &
    618618                                       kh(k,j,i) * g / vpt(k,j,i) *           &
     
    634634
    635635                      IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
    636                          k1 = 1.0 + 0.61 * q(k,j,i)
    637                          k2 = 0.61 * pt(k,j,i)
     636                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     637                         k2 = 0.61_wp * pt(k,j,i)
    638638                      ELSE IF ( cloud_physics )  THEN
    639                          IF ( ql(k,j,i) == 0.0 )  THEN
    640                             k1 = 1.0 + 0.61 * q(k,j,i)
    641                             k2 = 0.61 * pt(k,j,i)
     639                         IF ( ql(k,j,i) == 0.0_wp )  THEN
     640                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     641                            k2 = 0.61_wp * pt(k,j,i)
    642642                         ELSE
    643643                            theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
    644644                            temp  = theta * t_d_pt(k)
    645                             k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
    646                                        ( q(k,j,i) - ql(k,j,i) ) *          &
    647                                  ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
    648                                  ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
     645                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
     646                                          ( q(k,j,i) - ql(k,j,i) ) *          &
     647                                 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
     648                                 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
    649649                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    650                             k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
     650                            k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
    651651                         ENDIF
    652652                      ELSE IF ( cloud_droplets )  THEN
    653                          k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
    654                          k2 = 0.61 * pt(k,j,i)
     653                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
     654                         k2 = 0.61_wp * pt(k,j,i)
    655655                      ENDIF
    656656
     
    668668
    669669                      IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
    670                          k1 = 1.0 + 0.61 * q(k,j,i)
    671                          k2 = 0.61 * pt(k,j,i)
     670                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     671                         k2 = 0.61_wp * pt(k,j,i)
    672672                      ELSE IF ( cloud_physics )  THEN
    673                          IF ( ql(k,j,i) == 0.0 )  THEN
    674                             k1 = 1.0 + 0.61 * q(k,j,i)
    675                             k2 = 0.61 * pt(k,j,i)
     673                         IF ( ql(k,j,i) == 0.0_wp )  THEN
     674                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     675                            k2 = 0.61_wp * pt(k,j,i)
    676676                         ELSE
    677677                            theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
    678678                            temp  = theta * t_d_pt(k)
    679                             k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
    680                                        ( q(k,j,i) - ql(k,j,i) ) *          &
    681                                  ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
    682                                  ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
     679                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
     680                                          ( q(k,j,i) - ql(k,j,i) ) *          &
     681                                 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
     682                                 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
    683683                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    684                             k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
     684                            k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
    685685                         ENDIF
    686686                      ELSE IF ( cloud_droplets )  THEN
    687                          k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
    688                          k2 = 0.61 * pt(k,j,i)
     687                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
     688                         k2 = 0.61_wp * pt(k,j,i)
    689689                      ENDIF
    690690
     
    783783                IF ( k >= nzb_diff_s_outer(j,i) )  THEN
    784784
    785                    dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    786                    dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    787                                     u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    788                    dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    789                                     u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    790 
    791                    dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    792                                     v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    793                    dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    794                    dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    795                                     v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    796 
    797                    dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    798                                     w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    799                    dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    800                                     w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    801                    dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
    802 
    803                    def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
     785                   dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
     786                   dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     787                                       u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     788                   dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     789                                       u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     790
     791                   dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     792                                       v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     793                   dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     794                   dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     795                                       v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     796
     797                   dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     798                                       w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     799                   dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     800                                       w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     801                   dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
     802
     803                   def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
    804804                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
    805                          dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
    806 
    807                    IF ( def < 0.0 )  def = 0.0
     805                         dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
     806
     807                   IF ( def < 0.0_wp )  def = 0.0_wp
    808808
    809809                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
     
    825825                DO  k = 1, nzt
    826826
    827                    IF ( ( wall_e_x(j,i) /= 0.0 ).OR.( wall_e_y(j,i) /= 0.0 ) ) &
     827                   IF ( ( wall_e_x(j,i) /= 0.0_wp ).OR.( wall_e_y(j,i) /= 0.0_wp ) ) &
    828828                   THEN
    829829
    830830                      IF ( k == nzb_diff_s_inner(j,i) - 1 )  THEN
    831831                         dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
    832                          dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    833                                         u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
     832                         dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     833                                           u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
    834834                         dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
    835                          dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    836                                         v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
     835                         dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     836                                           v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
    837837                         dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
    838838
    839                          IF ( wall_e_y(j,i) /= 0.0 )  THEN
     839                         IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
    840840!
    841841!--                         Inconsistency removed: as the thermal stratification is
     
    852852!                                                wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
    853853                            km_neutral = kappa *                                    &
    854                                         ( usvs(k,j,i)**2 + wsvs(k,j,i)**2 )**0.25 * &
    855                                          0.5 * dy
    856                             IF ( km_neutral > 0.0 )  THEN
     854                                        ( usvs(k,j,i)**2 + wsvs(k,j,i)**2 )**0.25_wp * &
     855                                         0.5_wp * dy
     856                            IF ( km_neutral > 0.0_wp )  THEN
    857857                               dudy = - wall_e_y(j,i) * usvs(k,j,i) / km_neutral
    858858                               dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km_neutral
    859859                            ELSE
    860                                dudy = 0.0
    861                                dwdy = 0.0
     860                               dudy = 0.0_wp
     861                               dwdy = 0.0_wp
    862862                            ENDIF
    863863                         ELSE
    864                             dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    865                                             u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    866                             dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    867                                             w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     864                            dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     865                                               u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     866                            dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     867                                               w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    868868                         ENDIF
    869869
    870                          IF ( wall_e_x(j,i) /= 0.0 )  THEN
     870                         IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
    871871!
    872872!--                         Inconsistency removed: as the thermal stratification is
     
    883883!                                                wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
    884884                            km_neutral = kappa *                                     &
    885                                          ( vsus(k,j,i)**2 + wsus(k,j,i)**2 )**0.25 * &
    886                                          0.5 * dx
    887                             IF ( km_neutral > 0.0 )  THEN
     885                                         ( vsus(k,j,i)**2 + wsus(k,j,i)**2 )**0.25_wp * &
     886                                         0.5_wp * dx
     887                            IF ( km_neutral > 0.0_wp )  THEN
    888888                               dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km_neutral
    889889                               dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km_neutral
    890890                            ELSE
    891                                dvdx = 0.0
    892                                dwdx = 0.0
     891                               dvdx = 0.0_wp
     892                               dwdx = 0.0_wp
    893893                            ENDIF
    894894                         ELSE
    895                             dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    896                                             v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    897                             dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    898                                             w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     895                            dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     896                                               v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     897                            dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     898                                               w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    899899                         ENDIF
    900900
    901                          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
     901                         def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
    902902                               dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
    903                                dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
    904 
    905                          IF ( def < 0.0 )  def = 0.0
     903                               dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
     904
     905                         IF ( def < 0.0_wp )  def = 0.0_wp
    906906
    907907                         tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
     
    919919
    920920                         dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
    921                          dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    922                                         u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    923                          dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    924                          dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    925                                         v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     921                         dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     922                                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     923                         dvdy =          ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     924                         dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     925                                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    926926                         dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
    927927
    928                          IF ( wall_e_y(j,i) /= 0.0 )  THEN
     928                         IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
    929929!
    930930!--                         Inconsistency removed: as the thermal stratification
     
    937937!--                               validation has been available
    938938                            km_neutral = kappa * ( usvs(k,j,i)**2 + &
    939                                                    wsvs(k,j,i)**2 )**0.25 * 0.5 * dy
    940                             IF ( km_neutral > 0.0 )  THEN
     939                                                   wsvs(k,j,i)**2 )**0.25_wp * 0.5_wp * dy
     940                            IF ( km_neutral > 0.0_wp )  THEN
    941941                               dudy = - wall_e_y(j,i) * usvs(k,j,i) / km_neutral
    942942                               dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km_neutral
    943943                            ELSE
    944                                dudy = 0.0
    945                                dwdy = 0.0
     944                               dudy = 0.0_wp
     945                               dwdy = 0.0_wp
    946946                            ENDIF
    947947                         ELSE
    948                             dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    949                                             u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    950                             dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    951                                             w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     948                            dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     949                                               u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     950                            dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     951                                               w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    952952                         ENDIF
    953953
    954                          IF ( wall_e_x(j,i) /= 0.0 )  THEN
     954                         IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
    955955!
    956956!--                         Inconsistency removed: as the thermal stratification
     
    963963!--                               validation has been available
    964964                            km_neutral = kappa * ( vsus(k,j,i)**2 + &
    965                                                    wsus(k,j,i)**2 )**0.25 * 0.5 * dx
    966                             IF ( km_neutral > 0.0 )  THEN
     965                                                   wsus(k,j,i)**2 )**0.25_wp * 0.5_wp * dx
     966                            IF ( km_neutral > 0.0_wp )  THEN
    967967                               dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km_neutral
    968968                               dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km_neutral
    969969                            ELSE
    970                                dvdx = 0.0
    971                                dwdx = 0.0
     970                               dvdx = 0.0_wp
     971                               dwdx = 0.0_wp
    972972                            ENDIF
    973973                         ELSE
    974                             dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    975                                             v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    976                             dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    977                                             w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     974                            dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     975                                               v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     976                            dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     977                                               w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    978978                         ENDIF
    979979
    980                          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
     980                         def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
    981981                              dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 +  &
    982                               dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
    983 
    984                          IF ( def < 0.0 )  def = 0.0
     982                              dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
     983
     984                         IF ( def < 0.0_wp )  def = 0.0_wp
    985985
    986986                         tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
     
    993993                      IF ( k == nzb_diff_s_outer(j,i)-1 )  THEN
    994994
    995                          dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    996                          dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    997                                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    998                          dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    999                                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    1000 
    1001                          dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    1002                                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    1003                          dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    1004                          dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    1005                                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    1006 
    1007                          dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    1008                                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    1009                          dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    1010                                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    1011                          dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
    1012 
    1013                          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
     995                         dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
     996                         dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     997                                             u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     998                         dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     999                                             u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     1000
     1001                         dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     1002                                             v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     1003                         dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     1004                         dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     1005                                             v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     1006
     1007                         dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     1008                                             w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     1009                         dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     1010                                             w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     1011                         dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
     1012
     1013                         def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
    10141014                               dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
    1015                                dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
    1016 
    1017                          IF ( def < 0.0 )  def = 0.0
     1015                               dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
     1016
     1017                         IF ( def < 0.0_wp )  def = 0.0_wp
    10181018
    10191019                         tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
     
    10351035                DO  k = 1, nzt
    10361036
    1037                    IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) &
     1037                   IF ( ( wall_e_x(j,i) == 0.0_wp ) .AND. ( wall_e_y(j,i) == 0.0_wp ) ) &
    10381038                   THEN
    10391039
    10401040                      IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
    10411041
    1042                          dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    1043                          dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    1044                                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    1045                          dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    1046                                           u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
    1047 
    1048                          dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    1049                                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    1050                          dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    1051                          dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    1052                                           v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
    1053 
    1054                          dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    1055                                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    1056                          dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    1057                                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    1058                          dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
    1059 
    1060                          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
     1042                         dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
     1043                         dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     1044                                             u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     1045                         dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     1046                                             u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
     1047
     1048                         dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     1049                                             v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     1050                         dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     1051                         dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     1052                                             v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
     1053
     1054                         dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     1055                                             w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     1056                         dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     1057                                             w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     1058                         dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
     1059
     1060                         def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
    10611061                               dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
    1062                                dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
    1063 
    1064                          IF ( def < 0.0 )  def = 0.0
     1062                               dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
     1063
     1064                         IF ( def < 0.0_wp )  def = 0.0_wp
    10651065
    10661066                         tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
     
    10821082                   IF ( k == nzb_diff_s_outer(j,i)-1 )  THEN
    10831083
    1084                       dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    1085                       dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    1086                                        u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    1087                       dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    1088                                        u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    1089 
    1090                       dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    1091                                        v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    1092                       dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    1093                       dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    1094                                        v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    1095 
    1096                       dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    1097                                        w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    1098                       dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    1099                                        w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    1100                       dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
    1101 
    1102                       def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
     1084                      dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
     1085                      dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     1086                                          u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     1087                      dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     1088                                          u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     1089
     1090                      dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     1091                                          v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     1092                      dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     1093                      dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     1094                                          v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     1095
     1096                      dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     1097                                          w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     1098                      dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     1099                                          w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     1100                      dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
     1101
     1102                      def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
    11031103                            dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
    1104                             dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
    1105 
    1106                       IF ( def < 0.0 )  def = 0.0
     1104                            dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
     1105
     1106                      IF ( def < 0.0_wp )  def = 0.0_wp
    11071107
    11081108                      tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
     
    12321232!
    12331233!                         IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
    1234 !                            k1 = 1.0 + 0.61 * q(k,j,i)
    1235 !                            k2 = 0.61 * pt(k,j,i)
     1234!                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     1235!                            k2 = 0.61_wp * pt(k,j,i)
    12361236!                            tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
    12371237!                                            g / vpt(k,j,i) *                      &
     
    12401240!                                            ) * dd2zu(k)
    12411241!                         ELSE IF ( cloud_physics )  THEN
    1242 !                            IF ( ql(k,j,i) == 0.0 )  THEN
    1243 !                               k1 = 1.0 + 0.61 * q(k,j,i)
    1244 !                               k2 = 0.61 * pt(k,j,i)
     1242!                            IF ( ql(k,j,i) == 0.0_wp )  THEN
     1243!                               k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     1244!                               k2 = 0.61_wp * pt(k,j,i)
    12451245!                            ELSE
    12461246!                               theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
    12471247!                               temp  = theta * t_d_pt(k)
    1248 !                               k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
    1249 !                                          ( q(k,j,i) - ql(k,j,i) ) *          &
    1250 !                                    ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
    1251 !                                    ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
     1248!                               k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
     1249!                                             ( q(k,j,i) - ql(k,j,i) ) *          &
     1250!                                    ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
     1251!                                    ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
    12521252!                                    ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    1253 !                               k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
     1253!                               k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
    12541254!                            ENDIF
    12551255!                            tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
     
    12591259!                                            ) * dd2zu(k)
    12601260!                         ELSE IF ( cloud_droplets )  THEN
    1261 !                            k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
    1262 !                            k2 = 0.61 * pt(k,j,i)
     1261!                            k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
     1262!                            k2 = 0.61_wp * pt(k,j,i)
    12631263!                            tend(k,j,i) = tend(k,j,i) -                          &
    12641264!                                          kh(k,j,i) * g / vpt(k,j,i) *           &
     
    12891289!
    12901290!                            IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
    1291 !                               k1 = 1.0 + 0.61 * q(k,j,i)
    1292 !                               k2 = 0.61 * pt(k,j,i)
     1291!                               k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     1292!                               k2 = 0.61_wp * pt(k,j,i)
    12931293!                            ELSE IF ( cloud_physics )  THEN
    1294 !                               IF ( ql(k,j,i) == 0.0 )  THEN
    1295 !                                  k1 = 1.0 + 0.61 * q(k,j,i)
    1296 !                                  k2 = 0.61 * pt(k,j,i)
     1294!                               IF ( ql(k,j,i) == 0.0_wp )  THEN
     1295!                                  k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     1296!                                  k2 = 0.61_wp * pt(k,j,i)
    12971297!                               ELSE
    12981298!                                  theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
    12991299!                                  temp  = theta * t_d_pt(k)
    1300 !                                  k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
    1301 !                                             ( q(k,j,i) - ql(k,j,i) ) *          &
    1302 !                                       ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
    1303 !                                       ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
     1300!                                  k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *        &
     1301!                                                ( q(k,j,i) - ql(k,j,i) ) *    &
     1302!                                       ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /&
     1303!                                       ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
    13041304!                                       ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    1305 !                                  k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
     1305!                                  k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
    13061306!                               ENDIF
    13071307!                            ELSE IF ( cloud_droplets )  THEN
    1308 !                               k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
    1309 !                               k2 = 0.61 * pt(k,j,i)
     1308!                               k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
     1309!                               k2 = 0.61_wp * pt(k,j,i)
    13101310!                            ENDIF
    13111311!
     
    13301330!
    13311331!                            IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
    1332 !                               k1 = 1.0 + 0.61 * q(k,j,i)
    1333 !                               k2 = 0.61 * pt(k,j,i)
     1332!                               k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     1333!                               k2 = 0.61_wp * pt(k,j,i)
    13341334!                            ELSE IF ( cloud_physics )  THEN
    1335 !                               IF ( ql(k,j,i) == 0.0 )  THEN
    1336 !                                  k1 = 1.0 + 0.61 * q(k,j,i)
    1337 !                                  k2 = 0.61 * pt(k,j,i)
     1335!                               IF ( ql(k,j,i) == 0.0_wp )  THEN
     1336!                                  k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     1337!                                  k2 = 0.61_wp * pt(k,j,i)
    13381338!                               ELSE
    13391339!                                  theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
    13401340!                                  temp  = theta * t_d_pt(k)
    1341 !                                  k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
    1342 !                                             ( q(k,j,i) - ql(k,j,i) ) *          &
    1343 !                                       ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
    1344 !                                       ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
     1341!                                  k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *        &
     1342!                                                ( q(k,j,i) - ql(k,j,i) ) *    &
     1343!                                       ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /&
     1344!                                       ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
    13451345!                                       ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    1346 !                                  k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
     1346!                                  k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
    13471347!                               ENDIF
    13481348!                            ELSE IF ( cloud_droplets )  THEN
    1349 !                               k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
    1350 !                               k2 = 0.61 * pt(k,j,i)
     1349!                               k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
     1350!                               k2 = 0.61_wp * pt(k,j,i)
    13511351!                            ENDIF
    13521352!
     
    14261426       DO  k = nzb_diff_s_outer(j,i), nzt
    14271427
    1428           dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    1429           dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    1430                            u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    1431           dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    1432                            u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    1433 
    1434           dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    1435                            v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    1436           dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    1437           dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    1438                            v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    1439 
    1440           dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    1441                            w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    1442           dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    1443                            w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    1444           dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
    1445 
    1446           def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
    1447                 + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
    1448                 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
    1449 
    1450           IF ( def < 0.0 )  def = 0.0
     1428          dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
     1429          dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     1430                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     1431          dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     1432                              u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     1433
     1434          dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     1435                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     1436          dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     1437          dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     1438                              v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     1439
     1440          dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     1441                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     1442          dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     1443                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     1444          dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
     1445
     1446          def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
     1447                + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2    &
     1448                + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
     1449
     1450          IF ( def < 0.0_wp )  def = 0.0_wp
    14511451
    14521452          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
     
    14561456       IF ( prandtl_layer )  THEN
    14571457
    1458           IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) )  THEN
     1458          IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) )  THEN
    14591459
    14601460!
     
    14651465
    14661466             dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
    1467              dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    1468                             u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
    1469              dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    1470              dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    1471                             v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
     1467             dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     1468                               u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
     1469             dvdy =          ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     1470             dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     1471                               v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
    14721472             dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
    14731473
    1474              IF ( wall_e_y(j,i) /= 0.0 )  THEN
     1474             IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
    14751475!
    14761476!--             Inconsistency removed: as the thermal stratification
     
    14861486                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
    14871487                                    wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
    1488                 km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25 * &
    1489                              0.5 * dy
    1490                 IF ( km_neutral > 0.0 )  THEN
     1488                km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25_wp * &
     1489                             0.5_wp * dy
     1490                IF ( km_neutral > 0.0_wp )  THEN
    14911491                   dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
    14921492                   dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
    14931493                ELSE
    1494                    dudy = 0.0
    1495                    dwdy = 0.0
     1494                   dudy = 0.0_wp
     1495                   dwdy = 0.0_wp
    14961496                ENDIF
    14971497             ELSE
    1498                 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    1499                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    1500                 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    1501                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     1498                dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     1499                                   u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     1500                dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     1501                                   w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    15021502             ENDIF
    15031503
    1504              IF ( wall_e_x(j,i) /= 0.0 )  THEN
     1504             IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
    15051505!
    15061506!--             Inconsistency removed: as the thermal stratification
     
    15161516                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
    15171517                                    wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
    1518                 km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25 * &
    1519                              0.5 * dx
    1520                 IF ( km_neutral > 0.0 )  THEN
     1518                km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25_wp * &
     1519                             0.5_wp * dx
     1520                IF ( km_neutral > 0.0_wp )  THEN
    15211521                   dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
    15221522                   dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
    15231523                ELSE
    1524                    dvdx = 0.0
    1525                    dwdx = 0.0
     1524                   dvdx = 0.0_wp
     1525                   dwdx = 0.0_wp
    15261526                ENDIF
    15271527             ELSE
    1528                 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    1529                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    1530                 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    1531                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     1528                dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     1529                                   v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     1530                dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     1531                                   w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    15321532             ENDIF
    15331533
    1534              def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
     1534             def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
    15351535                   dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
    1536                    dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
    1537 
    1538              IF ( def < 0.0 )  def = 0.0
     1536                   dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
     1537
     1538             IF ( def < 0.0_wp )  def = 0.0_wp
    15391539
    15401540             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
     
    15491549
    15501550                dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
    1551                 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    1552                                u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    1553                 dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    1554                 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    1555                                v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     1551                dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     1552                                  u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     1553                dvdy =          ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     1554                dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     1555                                  v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    15561556                dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
    15571557
    1558                 IF ( wall_e_y(j,i) /= 0.0 )  THEN
     1558                IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
    15591559!
    15601560!--                Inconsistency removed: as the thermal stratification
     
    15671567!--                      validation has been available
    15681568                   km_neutral = kappa * ( usvs(k)**2 + &
    1569                                           wsvs(k)**2 )**0.25 * 0.5 * dy
    1570                    IF ( km_neutral > 0.0 )  THEN
     1569                                          wsvs(k)**2 )**0.25_wp * 0.5_wp * dy
     1570                   IF ( km_neutral > 0.0_wp )  THEN
    15711571                      dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
    15721572                      dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
    15731573                   ELSE
    1574                       dudy = 0.0
    1575                       dwdy = 0.0
     1574                      dudy = 0.0_wp
     1575                      dwdy = 0.0_wp
    15761576                   ENDIF
    15771577                ELSE
    1578                    dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    1579                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    1580                    dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    1581                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     1578                   dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     1579                                      u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     1580                   dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     1581                                      w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    15821582                ENDIF
    15831583
    1584                 IF ( wall_e_x(j,i) /= 0.0 )  THEN
     1584                IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
    15851585!
    15861586!--                Inconsistency removed: as the thermal stratification
     
    15931593!--                      validation has been available
    15941594                   km_neutral = kappa * ( vsus(k)**2 + &
    1595                                           wsus(k)**2 )**0.25 * 0.5 * dx
    1596                    IF ( km_neutral > 0.0 )  THEN
     1595                                          wsus(k)**2 )**0.25_wp * 0.5_wp * dx
     1596                   IF ( km_neutral > 0.0_wp )  THEN
    15971597                      dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
    15981598                      dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
    15991599                   ELSE
    1600                       dvdx = 0.0
    1601                       dwdx = 0.0
     1600                      dvdx = 0.0_wp
     1601                      dwdx = 0.0_wp
    16021602                   ENDIF
    16031603                ELSE
    1604                    dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    1605                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    1606                    dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    1607                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     1604                   dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     1605                                      v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     1606                   dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     1607                                      w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    16081608                ENDIF
    16091609
    1610                 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
     1610                def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
    16111611                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
    1612                       dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
    1613 
    1614                 IF ( def < 0.0 )  def = 0.0
     1612                      dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
     1613
     1614                IF ( def < 0.0_wp )  def = 0.0_wp
    16151615
    16161616                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
     
    16231623             k = nzb_diff_s_outer(j,i)-1
    16241624
    1625              dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    1626              dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    1627                               u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    1628              dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    1629                               u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    1630 
    1631              dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    1632                               v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    1633              dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    1634              dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    1635                               v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    1636 
    1637              dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    1638                               w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    1639              dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    1640                               w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    1641              dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
     1625             dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
     1626             dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     1627                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     1628             dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     1629                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     1630
     1631             dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     1632                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     1633             dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     1634             dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     1635                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     1636
     1637             dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     1638                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     1639             dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     1640                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     1641             dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
    16421642
    16431643             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
     
    16451645                   dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
    16461646
    1647              IF ( def < 0.0 )  def = 0.0
     1647             IF ( def < 0.0_wp )  def = 0.0_wp
    16481648
    16491649             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
     
    16571657             k = nzb_diff_s_inner(j,i)-1
    16581658
    1659              dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    1660              dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     1659             dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
     1660             dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     1661                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     1662             dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     1663                                 u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
     1664
     1665             dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     1666                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     1667             dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     1668             dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     1669                                 v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
     1670
     1671             dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     1672                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     1673             dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     1674                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     1675             dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
     1676
     1677             def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
     1678                   + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
     1679                   + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
     1680
     1681             IF ( def < 0.0_wp )  def = 0.0_wp
     1682
     1683             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
     1684
     1685          ENDIF
     1686
     1687       ELSEIF ( use_surface_fluxes )  THEN
     1688
     1689          k = nzb_diff_s_outer(j,i)-1
     1690
     1691          dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
     1692          dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    16611693                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    1662              dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    1663                               u_0(j,i)   - u_0(j,i+1)  ) * dd2zu(k)
    1664 
    1665              dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     1694          dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
     1695                              u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     1696
     1697          dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    16661698                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    1667              dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    1668              dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    1669                               v_0(j,i)   - v_0(j+1,i)  ) * dd2zu(k)
    1670 
    1671              dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
     1699          dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     1700          dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
     1701                              v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     1702
     1703          dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    16721704                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    1673              dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
     1705          dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    16741706                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    1675              dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
    1676 
    1677              def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
    1678                    + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
    1679                    + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
    1680 
    1681              IF ( def < 0.0 )  def = 0.0
    1682 
    1683              tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
    1684 
    1685           ENDIF
    1686 
    1687        ELSEIF ( use_surface_fluxes )  THEN
    1688 
    1689           k = nzb_diff_s_outer(j,i)-1
    1690 
    1691           dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    1692           dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
    1693                            u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    1694           dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
    1695                            u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    1696 
    1697           dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
    1698                            v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    1699           dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    1700           dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
    1701                            v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    1702 
    1703           dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
    1704                            w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    1705           dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
    1706                            w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    1707           dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
    1708 
    1709           def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
     1707          dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
     1708
     1709          def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
    17101710                dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
    1711                 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
    1712 
    1713           IF ( def < 0.0 )  def = 0.0
     1711                dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
     1712
     1713          IF ( def < 0.0_wp )  def = 0.0_wp
    17141714
    17151715          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
     
    17941794
    17951795                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
    1796                    k1 = 1.0 + 0.61 * q(k,j,i)
    1797                    k2 = 0.61 * pt(k,j,i)
     1796                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     1797                   k2 = 0.61_wp * pt(k,j,i)
    17981798                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *   &
    17991799                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
     
    18011801                                         ) * dd2zu(k)
    18021802                ELSE IF ( cloud_physics )  THEN
    1803                    IF ( ql(k,j,i) == 0.0 )  THEN
    1804                       k1 = 1.0 + 0.61 * q(k,j,i)
    1805                       k2 = 0.61 * pt(k,j,i)
     1803                   IF ( ql(k,j,i) == 0.0_wp )  THEN
     1804                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     1805                      k2 = 0.61_wp * pt(k,j,i)
    18061806                   ELSE
    18071807                      theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
    18081808                      temp  = theta * t_d_pt(k)
    1809                       k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
    1810                                  ( q(k,j,i) - ql(k,j,i) ) *          &
    1811                            ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
    1812                            ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
     1809                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
     1810                                    ( q(k,j,i) - ql(k,j,i) ) *          &
     1811                           ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
     1812                           ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
    18131813                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    1814                       k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
     1814                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
    18151815                   ENDIF
    18161816                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *   &
     
    18191819                                         ) * dd2zu(k)
    18201820                ELSE IF ( cloud_droplets )  THEN
    1821                    k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
    1822                    k2 = 0.61 * pt(k,j,i)
     1821                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
     1822                   k2 = 0.61_wp * pt(k,j,i)
    18231823                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *  &
    18241824                                     ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +    &
     
    18331833
    18341834                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
    1835                    k1 = 1.0 + 0.61 * q(k,j,i)
    1836                    k2 = 0.61 * pt(k,j,i)
     1835                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     1836                   k2 = 0.61_wp * pt(k,j,i)
    18371837                ELSE IF ( cloud_physics )  THEN
    1838                    IF ( ql(k,j,i) == 0.0 )  THEN
    1839                       k1 = 1.0 + 0.61 * q(k,j,i)
    1840                       k2 = 0.61 * pt(k,j,i)
     1838                   IF ( ql(k,j,i) == 0.0_wp )  THEN
     1839                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     1840                      k2 = 0.61_wp * pt(k,j,i)
    18411841                   ELSE
    18421842                      theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
    18431843                      temp  = theta * t_d_pt(k)
    1844                       k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
    1845                                  ( q(k,j,i) - ql(k,j,i) ) *          &
    1846                            ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
    1847                            ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
     1844                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
     1845                                    ( q(k,j,i) - ql(k,j,i) ) *          &
     1846                           ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
     1847                           ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
    18481848                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    1849                       k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
     1849                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
    18501850                   ENDIF
    18511851                ELSE IF ( cloud_droplets )  THEN
    1852                    k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
    1853                    k2 = 0.61 * pt(k,j,i)
     1852                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
     1853                   k2 = 0.61_wp * pt(k,j,i)
    18541854                ENDIF
    18551855
     
    18621862
    18631863                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
    1864                    k1 = 1.0 + 0.61 * q(k,j,i)
    1865                    k2 = 0.61 * pt(k,j,i)
     1864                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     1865                   k2 = 0.61_wp * pt(k,j,i)
    18661866                ELSE IF ( cloud_physics )  THEN
    1867                    IF ( ql(k,j,i) == 0.0 )  THEN
    1868                       k1 = 1.0 + 0.61 * q(k,j,i)
    1869                       k2 = 0.61 * pt(k,j,i)
     1867                   IF ( ql(k,j,i) == 0.0_wp )  THEN
     1868                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     1869                      k2 = 0.61_wp * pt(k,j,i)
    18701870                   ELSE
    18711871                      theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
    18721872                      temp  = theta * t_d_pt(k)
    1873                       k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
    1874                                  ( q(k,j,i) - ql(k,j,i) ) *          &
    1875                            ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
    1876                            ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
     1873                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
     1874                                    ( q(k,j,i) - ql(k,j,i) ) *          &
     1875                           ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
     1876                           ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
    18771877                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    1878                       k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
     1878                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
    18791879                   ENDIF
    18801880                ELSE IF ( cloud_droplets )  THEN
    1881                    k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
    1882                    k2 = 0.61 * pt(k,j,i)
     1881                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
     1882                   k2 = 0.61_wp * pt(k,j,i)
    18831883                ENDIF
    18841884
     
    19171917          IF ( first_call )  THEN
    19181918             ALLOCATE( u_0(nysg:nyng,nxlg:nxrg), v_0(nysg:nyng,nxlg:nxrg) )
    1919              u_0 = 0.0   ! just to avoid access of uninitialized memory
    1920              v_0 = 0.0   ! within exchange_horiz_2d
     1919             u_0 = 0.0_wp   ! just to avoid access of uninitialized memory
     1920             v_0 = 0.0_wp   ! within exchange_horiz_2d
    19211921             first_call = .FALSE.
    19221922          ENDIF
     
    19421942
    19431943                u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / &
    1944                                  ( 0.5 * ( km(ku,j,i) + km(ku,j,i-1) ) +       &
    1945                                    1.0E-20 )
     1944                                 ( 0.5_wp * ( km(ku,j,i) + km(ku,j,i-1) ) +    &
     1945                                   1.0E-20_wp )
    19461946!                                  ( us(j,i) * kappa * zu(1) )
    19471947                v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / &
    1948                                  ( 0.5 * ( km(kv,j,i) + km(kv,j-1,i) ) +       &
    1949                                    1.0E-20 )
     1948                                 ( 0.5_wp * ( km(kv,j,i) + km(kv,j-1,i) ) +    &
     1949                                   1.0E-20_wp )
    19501950!                                  ( us(j,i) * kappa * zu(1) )
    19511951
Note: See TracChangeset for help on using the changeset viewer.