Changeset 1342


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

REAL constants defined as wp-kind

Location:
palm/trunk/SOURCE
Files:
11 edited

Legend:

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

    r1325 r1342  
    2020! Current revisions:
    2121! -----------------
     22! REAL constants defined as wp-kinds
    2223!
    2324! Former revisions:
     
    315316!
    316317!-- Exponent for geometric average
    317     exponent = 1.0 / ( ny + 1.0 )
     318    exponent = 1.0_wp / ( ny + 1.0_wp )
    318319
    319320!
     
    346347          DO  i = 0, nx/2
    347348
    348              sums_spectra_l(i) = 1.0
     349             sums_spectra_l(i) = 1.0_wp
    349350             DO  j = nys_x, nyn_x
    350351                sums_spectra_l(i) = sums_spectra_l(i) * ddd(i,j,k)**exponent
     
    355356       ELSE
    356357
    357           sums_spectra_l = 1.0
     358          sums_spectra_l = 1.0_wp
    358359
    359360       ENDIF
     
    361362!
    362363!--    Global sum of spectra on PE0 (from where they are written on file)
    363        sums_spectra(:,n) = 0.0
     364       sums_spectra(:,n) = 0.0_wp
    364365#if defined( __parallel )   
    365366       CALL MPI_BARRIER( comm2d, ierr )  ! Necessary?
     
    380381!--    Temperton fft results need to be normalized
    381382       IF ( fft_method == 'temperton-algorithm' )  THEN
    382           fac = nx + 1.0
     383          fac = nx + 1.0_wp
    383384       ELSE
    384           fac = 1.0
     385          fac = 1.0_wp
    385386       ENDIF
    386387       DO  i = 1, nx/2
     
    454455!
    455456!-- Exponent for geometric average
    456     exponent = 1.0 / ( nx + 1.0 )
     457    exponent = 1.0_wp / ( nx + 1.0_wp )
    457458
    458459!
     
    485486          DO  j = 0, ny/2
    486487
    487              sums_spectra_l(j) = 1.0
     488             sums_spectra_l(j) = 1.0_wp
    488489             DO  i = nxl_yd, nxr_yd
    489490                sums_spectra_l(j) = sums_spectra_l(j) * ddd(j,i,k)**exponent
     
    494495       ELSE
    495496
    496           sums_spectra_l = 1.0
     497          sums_spectra_l = 1.0_wp
    497498
    498499       ENDIF
     
    500501!
    501502!--    Global sum of spectra on PE0 (from where they are written on file)
    502        sums_spectra(:,n) = 0.0
     503       sums_spectra(:,n) = 0.0_wp
    503504#if defined( __parallel )   
    504505       CALL MPI_BARRIER( comm2d, ierr )  ! Necessary?
     
    520521!--    Temperton fft results need to be normalized
    521522       IF ( fft_method == 'temperton-algorithm' )  THEN
    522           fac = ny + 1.0
     523          fac = ny + 1.0_wp
    523524       ELSE
    524           fac = 1.0
     525          fac = 1.0_wp
    525526       ENDIF
    526527       DO  j = 1, ny/2
  • palm/trunk/SOURCE/fft_xy.f90

    r1323 r1342  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants defined as wp-kind
    2323!
    2424! Former revisions:
     
    242242       IF ( fft_method == 'system-specific' )  THEN
    243243
    244           dnx = 1.0 / ( nx + 1.0 )
    245           dny = 1.0 / ( ny + 1.0 )
     244          dnx = 1.0_wp / ( nx + 1.0_wp )
     245          dny = 1.0_wp / ( ny + 1.0_wp )
    246246          sqr_dnx = SQRT( dnx )
    247247          sqr_dny = SQRT( dny )
     
    267267                    trig_yb(2*(ny+1)), trig_yf(2*(ny+1)) )
    268268
    269           work_x = 0.0
    270           work_y = 0.0
     269          work_x = 0.0_wp
     270          work_y = 0.0_wp
    271271          nz1  = nz + MOD( nz+1, 2 )  ! odd nz slows down fft significantly
    272272                                      ! when using the NEC ffts
     
    433433                DO  j = nys_x, nyn_x
    434434
    435                    cwork(0) = CMPLX( ar(0,j,k), 0.0 )
     435                   cwork(0) = CMPLX( ar(0,j,k), 0.0_wp )
    436436                   DO  i = 1, (nx+1)/2 - 1
    437437                      cwork(i)      = CMPLX( ar(i,j,k), -ar(nx+1-i,j,k) )
    438438                      cwork(nx+1-i) = CMPLX( ar(i,j,k),  ar(nx+1-i,j,k) )
    439439                   ENDDO
    440                    cwork((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
     440                   cwork((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp )
    441441
    442442                   ishape = SHAPE( cwork )
     
    494494                      work(2*i+1) = ar(nx+1-i,j,k)
    495495                   ENDDO
    496                    work(1)    = 0.0
    497                    work(nx+2) = 0.0
     496                   work(1)    = 0.0_wp
     497                   work(nx+2) = 0.0_wp
    498498
    499499                   CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
     
    551551                   IF ( PRESENT( ar_2d ) )  THEN
    552552
    553                       x_out(0) = CMPLX( ar_2d(0,j), 0.0 )
     553                      x_out(0) = CMPLX( ar_2d(0,j), 0.0_wp )
    554554                      DO  i = 1, (nx+1)/2 - 1
    555555                         x_out(i) = CMPLX( ar_2d(i,j), ar_2d(nx+1-i,j) )
    556556                      ENDDO
    557                       x_out((nx+1)/2) = CMPLX( ar_2d((nx+1)/2,j), 0.0 )
     557                      x_out((nx+1)/2) = CMPLX( ar_2d((nx+1)/2,j), 0.0_wp )
    558558
    559559                   ELSE
    560560
    561                       x_out(0) = CMPLX( ar(0,j,k), 0.0 )
     561                      x_out(0) = CMPLX( ar(0,j,k), 0.0_wp )
    562562                      DO  i = 1, (nx+1)/2 - 1
    563563                         x_out(i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
    564564                      ENDDO
    565                       x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
     565                      x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp )
    566566
    567567                   ENDIF
     
    614614                      work(2*i+1) = ar(nx+1-i,j,k)
    615615                   ENDDO
    616                    work(1) = 0.0
    617                    work(nx+2) = 0.0
     616                   work(1) = 0.0_wp
     617                   work(nx+2) = 0.0_wp
    618618
    619619                   CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx,      &
     
    667667                      work(2*i+1) = ar(nx+1-i,j,k)
    668668                   ENDDO
    669                    work(1) = 0.0
    670                    work(nx+2) = 0.0
     669                   work(1) = 0.0_wp
     670                   work(nx+2) = 0.0_wp
    671671
    672672                   CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
     
    711711                DO  j = nys_x, nyn_x
    712712
    713                    ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0 )
     713                   ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0_wp )
    714714
    715715                   DO  i = 1, (nx+1)/2 - 1
    716716                      ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
    717717                   ENDDO
    718                    ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
     718                   ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp )
    719719
    720720                ENDDO
     
    805805          ELSE
    806806
    807              cwork(0) = CMPLX( ar(0), 0.0 )
     807             cwork(0) = CMPLX( ar(0), 0.0_wp )
    808808             DO  i = 1, (nx+1)/2 - 1
    809809                cwork(i)      = CMPLX( ar(i), -ar(nx+1-i) )
    810810                cwork(nx+1-i) = CMPLX( ar(i),  ar(nx+1-i) )
    811811             ENDDO
    812              cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )
     812             cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0_wp )
    813813
    814814             ishape = SHAPE( cwork )
     
    848848                work(2*i+1) = ar(nx+1-i)
    849849             ENDDO
    850              work(1)    = 0.0
    851              work(nx+2) = 0.0
     850             work(1)    = 0.0_wp
     851             work(nx+2) = 0.0_wp
    852852
    853853             CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
     
    873873         ELSE
    874874
    875              x_out(0) = CMPLX( ar(0), 0.0 )
     875             x_out(0) = CMPLX( ar(0), 0.0_wp )
    876876             DO  i = 1, (nx+1)/2 - 1
    877877                x_out(i) = CMPLX( ar(i), ar(nx+1-i) )
    878878             ENDDO
    879              x_out((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )
     879             x_out((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0_wp )
    880880
    881881             CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in)
     
    908908                work(2*i+1) = ar(nx+1-i)
    909909             ENDDO
    910              work(1) = 0.0
    911              work(nx+2) = 0.0
     910             work(1) = 0.0_wp
     911             work(nx+2) = 0.0_wp
    912912
    913913             CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
     
    941941                work(2*i+1) = ar(nx+1-i)
    942942             ENDDO
    943              work(1) = 0.0
    944              work(nx+2) = 0.0
     943             work(1) = 0.0_wp
     944             work(nx+2) = 0.0_wp
    945945
    946946             CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
     
    10761076                DO  i = nxl_y_l, nxr_y_l
    10771077
    1078                    cwork(0) = CMPLX( ar_tr(0,i,k), 0.0 )
     1078                   cwork(0) = CMPLX( ar_tr(0,i,k), 0.0_wp )
    10791079                   DO  j = 1, (ny+1)/2 - 1
    10801080                      cwork(j)      = CMPLX( ar_tr(j,i,k), -ar_tr(ny+1-j,i,k) )
    10811081                      cwork(ny+1-j) = CMPLX( ar_tr(j,i,k),  ar_tr(ny+1-j,i,k) )
    10821082                   ENDDO
    1083                    cwork((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0 )
     1083                   cwork((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0_wp )
    10841084
    10851085                   jshape = SHAPE( cwork )
     
    11371137                      work(2*j+1) = ar_tr(ny+1-j,i,k)
    11381138                   ENDDO
    1139                    work(1)    = 0.0
    1140                    work(ny+2) = 0.0
     1139                   work(1)    = 0.0_wp
     1140                   work(ny+2) = 0.0_wp
    11411141
    11421142                   CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
     
    11801180                DO  i = nxl_y_l, nxr_y_l
    11811181
    1182                    y_out(0) = CMPLX( ar_tr(0,i,k), 0.0 )
     1182                   y_out(0) = CMPLX( ar_tr(0,i,k), 0.0_wp )
    11831183                   DO  j = 1, (ny+1)/2 - 1
    11841184                      y_out(j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k) )
    11851185                   ENDDO
    1186                    y_out((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0 )
     1186                   y_out((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0_wp )
    11871187
    11881188                   CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in )
     
    12331233                      work(2*j+1) = ar_tr(ny+1-j,i,k)
    12341234                   ENDDO
    1235                    work(1)    = 0.0
    1236                    work(ny+2) = 0.0
     1235                   work(1)    = 0.0_wp
     1236                   work(ny+2) = 0.0_wp
    12371237
    12381238                   CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny,      &
     
    12841284                      work(2*j+1) = ar_tr(ny+1-j,i,k)
    12851285                   ENDDO
    1286                    work(1) = 0.0
    1287                    work(ny+2) = 0.0
     1286                   work(1) = 0.0_wp
     1287                   work(ny+2) = 0.0_wp
    12881288
    12891289                   CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
     
    13271327                DO  i = nxl_y, nxr_y
    13281328
    1329                    ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0 )
     1329                   ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0_wp )
    13301330
    13311331                   DO  j = 1, (ny+1)/2 - 1
    13321332                      ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
    13331333                   ENDDO
    1334                    ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
     1334                   ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0_wp )
    13351335
    13361336                ENDDO
     
    14231423          ELSE
    14241424
    1425              cwork(0) = CMPLX( ar(0), 0.0 )
     1425             cwork(0) = CMPLX( ar(0), 0.0_wp )
    14261426             DO  j = 1, (ny+1)/2 - 1
    14271427                cwork(j)      = CMPLX( ar(j), -ar(ny+1-j) )
    14281428                cwork(ny+1-j) = CMPLX( ar(j),  ar(ny+1-j) )
    14291429             ENDDO
    1430              cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 )
     1430             cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0_wp )
    14311431
    14321432             jshape = SHAPE( cwork )
     
    14661466                work(2*j+1) = ar(ny+1-j)
    14671467             ENDDO
    1468              work(1)    = 0.0
    1469              work(ny+2) = 0.0
     1468             work(1)    = 0.0_wp
     1469             work(ny+2) = 0.0_wp
    14701470
    14711471             CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
     
    14911491          ELSE
    14921492
    1493              y_out(0) = CMPLX( ar(0), 0.0 )
     1493             y_out(0) = CMPLX( ar(0), 0.0_wp )
    14941494             DO  j = 1, (ny+1)/2 - 1
    14951495                y_out(j) = CMPLX( ar(j), ar(ny+1-j) )
    14961496             ENDDO
    1497              y_out((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 )
     1497             y_out((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0_wp )
    14981498
    14991499             CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in )
     
    15261526                work(2*j+1) = ar(ny+1-j)
    15271527             ENDDO
    1528              work(1)    = 0.0
    1529              work(ny+2) = 0.0
     1528             work(1)    = 0.0_wp
     1529             work(ny+2) = 0.0_wp
    15301530
    15311531             CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3,      &
     
    15591559                work(2*j+1) = ar(ny+1-j)
    15601560             ENDDO
    1561              work(1) = 0.0
    1562              work(ny+2) = 0.0
     1561             work(1) = 0.0_wp
     1562             work(ny+2) = 0.0_wp
    15631563
    15641564             CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
     
    16211621
    16221622             ai(0:nx,1:nz) = ar(0:nx,1:nz)
    1623              ai(nx+1:,:)   = 0.0
     1623             ai(nx+1:,:)   = 0.0_wp
    16241624
    16251625             CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
     
    16431643                   ai(2*i+1,k) = ar(nx+1-i,k)
    16441644                ENDDO
    1645                 ai(1,k) = 0.0
    1646                 ai(nx+2,k) = 0.0
     1645                ai(1,k) = 0.0_wp
     1646                ai(nx+2,k) = 0.0_wp
    16471647             ENDDO
    16481648
     
    16691669             ai(0:nx,1:nz) = ar(0:nx,1:nz)
    16701670             IF ( nz1 > nz )  THEN
    1671                 ai(:,nz1) = 0.0
     1671                ai(:,nz1) = 0.0_wp
    16721672             ENDIF
    16731673
     
    16931693
    16941694             IF ( nz1 > nz )  THEN
    1695                 work(:,nz1) = 0.0
     1695                work(:,nz1) = 0.0_wp
    16961696             ENDIF
    16971697             DO  k = 1, nz
    1698                 work(1,k) = CMPLX( ar(0,k), 0.0 )
     1698                work(1,k) = CMPLX( ar(0,k), 0.0_wp )
    16991699                DO  i = 1, (nx+1)/2 - 1
    17001700                   work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k) )
    17011701                ENDDO
    1702                 work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0 )
     1702                work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0_wp )
    17031703             ENDDO
    17041704
     
    17641764
    17651765             ai(0:ny,1:nz) = ar(0:ny,1:nz)
    1766              ai(ny+1:,:)   = 0.0
     1766             ai(ny+1:,:)   = 0.0_wp
    17671767
    17681768             CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
     
    17861786                   ai(2*j+1,k) = ar(ny+1-j,k)
    17871787                ENDDO
    1788                 ai(1,k) = 0.0
    1789                 ai(ny+2,k) = 0.0
     1788                ai(1,k) = 0.0_wp
     1789                ai(ny+2,k) = 0.0_wp
    17901790             ENDDO
    17911791
     
    18121812             ai(0:ny,1:nz) = ar(0:ny,1:nz)
    18131813             IF ( nz1 > nz )  THEN
    1814                 ai(:,nz1) = 0.0
     1814                ai(:,nz1) = 0.0_wp
    18151815             ENDIF
    18161816
     
    18361836
    18371837             IF ( nz1 > nz )  THEN
    1838                 work(:,nz1) = 0.0
     1838                work(:,nz1) = 0.0_wp
    18391839             ENDIF
    18401840             DO  k = 1, nz
    1841                 work(1,k) = CMPLX( ar(0,k), 0.0 )
     1841                work(1,k) = CMPLX( ar(0,k), 0.0_wp )
    18421842                DO  j = 1, (ny+1)/2 - 1
    18431843                   work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k) )
    18441844                ENDDO
    1845                 work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0 )
     1845                work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0_wp )
    18461846             ENDDO
    18471847
  • palm/trunk/SOURCE/pres.f90

    r1321 r1342  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! REAL constants defined as wp-kind
    2323!
    2424! Former revisions:
     
    135135
    136136
    137     ddt_3d = 1.0 / dt_3d
    138     d_weight_pres = 1.0 / weight_pres(intermediate_timestep_count)
     137    ddt_3d = 1.0_wp / dt_3d
     138    d_weight_pres = 1.0_wp / weight_pres(intermediate_timestep_count)
    139139
    140140!
     
    180180    IF ( conserve_volume_flow  .AND.  ( outflow_l .OR. outflow_r ) )  THEN
    181181
    182        volume_flow(1)   = 0.0
    183        volume_flow_l(1) = 0.0
     182       volume_flow(1)   = 0.0_wp
     183       volume_flow_l(1) = 0.0_wp
    184184
    185185       IF ( outflow_l )  THEN
     
    219219    IF ( conserve_volume_flow  .AND.  ( outflow_n .OR. outflow_s ) )  THEN
    220220
    221        volume_flow(2)   = 0.0
    222        volume_flow_l(2) = 0.0
     221       volume_flow(2)   = 0.0_wp
     222       volume_flow_l(2) = 0.0_wp
    223223
    224224       IF ( outflow_s )  THEN
     
    257257!-- Remove mean vertical velocity
    258258    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
    259        IF ( simulated_time > 0.0 )  THEN ! otherwise nzb_w_inner not yet known
    260           w_l = 0.0;  w_l_l = 0.0
     259       IF ( simulated_time > 0.0_wp )  THEN ! otherwise nzb_w_inner not yet known
     260          w_l = 0.0_wp;  w_l_l = 0.0_wp
    261261          DO  i = nxl, nxr
    262262             DO  j = nys, nyn
     
    295295          DO  j = nys-1, nyn+1
    296296             DO  k = nzb, nzt+1
    297                 d(k,j,i) = 0.0
     297                d(k,j,i) = 0.0_wp
    298298             ENDDO
    299299          ENDDO
     
    305305          DO  j = nys, nyn
    306306             DO  k = nzb+1, nzt
    307                 d(k,j,i) = 0.0
     307                d(k,j,i) = 0.0_wp
    308308             ENDDO
    309309          ENDDO
     
    312312    ENDIF
    313313
    314     localsum  = 0.0
    315     threadsum = 0.0
     314    localsum  = 0.0_wp
     315    threadsum = 0.0_wp
    316316
    317317#if defined( __ibm )
     
    432432          DO  i = nxlg, nxrg
    433433             DO  j = nysg, nyng
    434                 tend(nzb_s_inner(j,i),j,i) = 0.0
     434                tend(nzb_s_inner(j,i),j,i) = 0.0_wp
    435435             ENDDO
    436436          ENDDO
     
    460460          DO  i = nxlg, nxrg
    461461             DO  j = nysg, nyng
    462                 tend(nzt+1,j,i) = 0.0
     462                tend(nzt+1,j,i) = 0.0_wp
    463463             ENDDO
    464464          ENDDO
     
    577577!-- pressure just computed
    578578    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc .OR. bc_ns_cyc ) )  THEN
    579        volume_flow_l(1) = 0.0
    580        volume_flow_l(2) = 0.0
     579       volume_flow_l(1) = 0.0_wp
     580       volume_flow_l(2) = 0.0_wp
    581581    ENDIF
    582582
     
    704704!-- a possible PE-sum is computed in flow_statistics
    705705    CALL cpu_log( log_point_s(1), 'divergence', 'start' )
    706     sums_divnew_l = 0.0
     706    sums_divnew_l = 0.0_wp
    707707
    708708!
    709709!-- d must be reset to zero because it can contain nonzero values below the
    710710!-- topography
    711     IF ( topography /= 'flat' )  d = 0.0
    712 
    713     localsum  = 0.0
    714     threadsum = 0.0
     711    IF ( topography /= 'flat' )  d = 0.0_wp
     712
     713    localsum  = 0.0_wp
     714    threadsum = 0.0_wp
    715715
    716716    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
  • 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
  • palm/trunk/SOURCE/random_function.f90

    r1321 r1342  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants defined as wp-kind
    2323!
    2424! Former revisions:
     
    104104       REAL(wp)     ::  rnmx             !:
    105105
    106        PARAMETER ( ia=16807, im=2147483647, am=1.0/im, iq=127773, ir=2836, &
    107                    ntab=32, ndiv=1+(im-1)/ntab, eps=1.2e-7, rnmx=1.0-eps )
     106       PARAMETER ( ia=16807, im=2147483647, am=1.0_wp/im, iq=127773, ir=2836, &
     107                   ntab=32, ndiv=1+(im-1)/ntab, eps=1.2e-7_wp, rnmx=1.0_wp-eps )
    108108
    109109       IF ( idum .le. 0  .or.  random_iy .eq. 0 )  THEN
  • palm/trunk/SOURCE/random_gauss.f90

    r1321 r1342  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants defined as wp-kind
    2323!
    2424! Former revisions:
     
    7878
    7979       IF ( iset == 0 )  THEN
    80           rsq = 0.0
    81           DO  WHILE ( rsq >= 1.0  .OR.  rsq == 0.0 )
    82              v1  = 2.0 * random_function( idum ) - 1.0
    83              v2  = 2.0 * random_function( idum ) - 1.0
     80          rsq = 0.0_wp
     81          DO  WHILE ( rsq >= 1.0_wp  .OR.  rsq == 0.0_wp )
     82             v1  = 2.0_wp * random_function( idum ) - 1.0_wp
     83             v2  = 2.0_wp * random_function( idum ) - 1.0_wp
    8484             rsq = v1**2 + v2**2
    8585          ENDDO
    86           fac          = SQRT( -2.0 * LOG( rsq ) / rsq )
     86          fac          = SQRT( -2.0_wp * LOG( rsq ) / rsq )
    8787          gset         = v1 * fac
    88           random_gauss = v2 * fac + 1.0
     88          random_gauss = v2 * fac + 1.0_wp
    8989          iset         = 1
    9090       ELSE
    91           random_gauss = gset + 1.0
     91          random_gauss = gset + 1.0_wp
    9292          iset         = 0
    9393       ENDIF
    9494
    95        IF ( ABS( random_gauss - 1.0 ) <= upper_limit )  EXIT
     95       IF ( ABS( random_gauss - 1.0_wp ) <= upper_limit )  EXIT
    9696
    9797    ENDDO
  • palm/trunk/SOURCE/temperton_fft.f90

    r1323 r1342  
    44! Current revisions:
    55! -----------------
    6 !
     6! REAL constants defined as wp-kind
    77!
    88! Former revisions:
     
    167167!OCL NOVREC
    168168       DO j = 1, nvex
    169           a(i+inc) = 0.5*a(i)
     169          a(i+inc) = 0.5_wp*a(i)
    170170          i = i + jump
    171171       END DO
     
    173173       i = istart + n*inc
    174174       DO j = 1, nvex
    175           a(i) = 0.5*a(i)
     175          a(i) = 0.5_wp*a(i)
    176176          i = i + jump
    177177       END DO
     
    223223!OCL NOVREC
    224224       DO j = 1, nvex
    225           a(ix) = 0.0
    226           a(ix+inc) = 0.0
     225          a(ix) = 0.0_wp
     226          a(ix+inc) = 0.0_wp
    227227          ix = ix + jump
    228228       END DO
     
    285285       DO j = 1, nvex
    286286          a(ix) = a(ix+inc)
    287           a(ix+inc) = 0.0
     287          a(ix+inc) = 0.0_wp
    288288          ix = ix + jump
    289289       END DO
     
    291291       iz = istart + (n+1)*inc
    292292       DO j = 1, nvex
    293           a(iz) = 0.0
     293          a(iz) = 0.0_wp
    294294          iz = iz + jump
    295295       END DO
     
    552552    GO TO 170
    55355330  CONTINUE
    554     z = 1.0_wp/REAL(n)
     554    z = 1.0_wp/REAL(n,KIND=wp)
    555555    DO l = 1, la
    556556       i = ibase
     
    590590       DO ijk = 1, lot
    591591          c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
    592           c(jb+j) = a(ia+i) - 0.5*(a(ib+i)+a(ic+i))
     592          c(jb+j) = a(ia+i) - 0.5_wp*(a(ib+i)+a(ic+i))
    593593          d(jb+j) = sin60*(a(ic+i)-a(ib+i))
    594594          i = i + inc3
     
    622622             a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c2*a(ic+i)+s2*b(ic+i))
    623623             b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c2*b(ic+i)-s2*a(ic+i))
    624              a2 = a(ia+i) - 0.5*a1
    625              b2 = b(ia+i) - 0.5*b1
     624             a2 = a(ia+i) - 0.5_wp*a1
     625             b2 = b(ia+i) - 0.5_wp*b1
    626626             a3 = sin60*((c1*a(ib+i)+s1*b(ib+i))-(c2*a(ic+i)+s2*b(ic+i)))
    627627             b3 = sin60*((c1*b(ib+i)-s1*a(ib+i))-(c2*b(ic+i)-s2*a(ic+i)))
     
    653653!OCL NOVREC
    654654       DO ijk = 1, lot
    655           c(ja+j) = a(ia+i) + 0.5*(a(ib+i)-a(ic+i))
     655          c(ja+j) = a(ia+i) + 0.5_wp*(a(ib+i)-a(ic+i))
    656656          d(ja+j) = -sin60*(a(ib+i)+a(ic+i))
    657657          c(jb+j) = a(ia+i) - (a(ib+i)-a(ic+i))
     
    665665    GO TO 170
    66666660  CONTINUE
    667     z = 1.0_wp/REAL(n)
     667    z = 1.0_wp/REAL(n,KIND=wp)
    668668    zsin60 = z*sin60
    669669    DO l = 1, la
     
    675675       DO ijk = 1, lot
    676676          c(ja+j) = z*(a(ia+i)+(a(ib+i)+a(ic+i)))
    677           c(jb+j) = z*(a(ia+i)-0.5*(a(ib+i)+a(ic+i)))
     677          c(jb+j) = z*(a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))
    678678          d(jb+j) = zsin60*(a(ic+i)-a(ib+i))
    679679          i = i + inc3
     
    794794    GO TO 170
    79579590  CONTINUE
    796     z = 1.0_wp/REAL(n)
     796    z = 1.0_wp/REAL(n,KIND=wp)
    797797    DO l = 1, la
    798798       i = ibase
     
    841841          a2 = a(ic+i) + a(id+i)
    842842          a4 = a(ic+i) - a(id+i)
    843           a5 = a(ia+i) - 0.25*(a1+a2)
     843          a5 = a(ia+i) - 0.25_wp*(a1+a2)
    844844          a6 = qrt5*(a1-a2)
    845845          c(ja+j) = a(ia+i) + (a1+a2)
     
    892892             b2 = (c2*b(ic+i)-s2*a(ic+i)) + (c3*b(id+i)-s3*a(id+i))
    893893             b4 = (c2*b(ic+i)-s2*a(ic+i)) - (c3*b(id+i)-s3*a(id+i))
    894              a5 = a(ia+i) - 0.25*(a1+a2)
     894             a5 = a(ia+i) - 0.25_wp*(a1+a2)
    895895             a6 = qrt5*(a1-a2)
    896              b5 = b(ia+i) - 0.25*(b1+b2)
     896             b5 = b(ia+i) - 0.25_wp*(b1+b2)
    897897             b6 = qrt5*(b1-b2)
    898898             a10 = a5 + a6
     
    941941          a2 = a(ic+i) + a(id+i)
    942942          a4 = a(ic+i) - a(id+i)
    943           a5 = a(ia+i) + 0.25*(a3-a4)
     943          a5 = a(ia+i) + 0.25_wp*(a3-a4)
    944944          a6 = qrt5*(a3+a4)
    945945          c(ja+j) = a5 + a6
     
    957957    GO TO 170
    958958120 CONTINUE
    959     z = 1.0_wp/REAL(n)
     959    z = 1.0_wp/REAL(n,KIND=wp)
    960960    zqrt5 = z*qrt5
    961961    zsin36 = z*sin36
     
    972972          a2 = a(ic+i) + a(id+i)
    973973          a4 = a(ic+i) - a(id+i)
    974           a5 = z*(a(ia+i)-0.25*(a1+a2))
     974          a5 = z*(a(ia+i)-0.25_wp*(a1+a2))
    975975          a6 = zqrt5*(a1-a2)
    976976          c(ja+j) = z*(a(ia+i)+(a1+a2))
     
    10141014          a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
    10151015          c(ja+j) = (a(ia+i)+a(id+i)) + a11
    1016           c(jc+j) = (a(ia+i)+a(id+i)-0.5*a11)
     1016          c(jc+j) = (a(ia+i)+a(id+i)-0.5_wp*a11)
    10171017          d(jc+j) = sin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
    10181018          a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
    1019           c(jb+j) = (a(ia+i)-a(id+i)) - 0.5*a11
     1019          c(jb+j) = (a(ia+i)-a(id+i)) - 0.5_wp*a11
    10201020          d(jb+j) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
    10211021          c(jd+j) = (a(ia+i)-a(id+i)) + a11
     
    10711071             b5 = c5*b(if+i) - s5*a(if+i)
    10721072             a11 = (a2+a5) + (a1+a4)
    1073              a20 = (a(ia+i)+a3) - 0.5*a11
     1073             a20 = (a(ia+i)+a3) - 0.5_wp*a11
    10741074             a21 = sin60*((a2+a5)-(a1+a4))
    10751075             b11 = (b2+b5) + (b1+b4)
    1076              b20 = (b(ia+i)+b3) - 0.5*b11
     1076             b20 = (b(ia+i)+b3) - 0.5_wp*b11
    10771077             b21 = sin60*((b2+b5)-(b1+b4))
    10781078             c(ja+j) = (a(ia+i)+a3) + a11
     
    10831083             d(je+j) = a21 - b20
    10841084             a11 = (a2-a5) + (a4-a1)
    1085              a20 = (a(ia+i)-a3) - 0.5*a11
     1085             a20 = (a(ia+i)-a3) - 0.5_wp*a11
    10861086             a21 = sin60*((a4-a1)-(a2-a5))
    10871087             b11 = (b5-b2) - (b4-b1)
    1088              b20 = (b3-b(ia+i)) - 0.5*b11
     1088             b20 = (b3-b(ia+i)) - 0.5_wp*b11
    10891089             b21 = sin60*((b5-b2)+(b4-b1))
    10901090             c(jb+j) = a20 - b21
     
    11181118!OCL NOVREC
    11191119       DO ijk = 1, lot
    1120           c(ja+j) = (a(ia+i)+0.5*(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i))
    1121           d(ja+j) = -(a(id+i)+0.5*(a(ib+i)+a(if+i))) - sin60*(a(ic+i)+a(ie+i))
     1120          c(ja+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i))
     1121          d(ja+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) - sin60*(a(ic+i)+a(ie+i))
    11221122          c(jb+j) = a(ia+i) - (a(ic+i)-a(ie+i))
    11231123          d(jb+j) = a(id+i) - (a(ib+i)+a(if+i))
    1124           c(jc+j) = (a(ia+i)+0.5*(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i))
    1125           d(jc+j) = -(a(id+i)+0.5*(a(ib+i)+a(if+i))) + sin60*(a(ic+i)+a(ie+i))
     1124          c(jc+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i))
     1125          d(jc+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) + sin60*(a(ic+i)+a(ie+i))
    11261126          i = i + inc3
    11271127          j = j + inc4
     
    11331133    GO TO 170
    11341134150 CONTINUE
    1135     z = 1.0_wp/REAL(n)
     1135    z = 1.0_wp/REAL(n,KIND=wp)
    11361136    zsin60 = z*sin60
    11371137    DO l = 1, la
     
    11441144          a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
    11451145          c(ja+j) = z*((a(ia+i)+a(id+i))+a11)
    1146           c(jc+j) = z*((a(ia+i)+a(id+i))-0.5*a11)
     1146          c(jc+j) = z*((a(ia+i)+a(id+i))-0.5_wp*a11)
    11471147          d(jc+j) = zsin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
    11481148          a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
    1149           c(jb+j) = z*((a(ia+i)-a(id+i))-0.5*a11)
     1149          c(jb+j) = z*((a(ia+i)-a(id+i))-0.5_wp*a11)
    11501150          d(jb+j) = zsin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
    11511151          c(jd+j) = z*((a(ia+i)-a(id+i))+a11)
     
    11761176    jd = jc + 2*m*inc2
    11771177    je = jd + 2*m*inc2
    1178     z = 1.0_wp/REAL(n)
     1178    z = 1.0_wp/REAL(n,KIND=wp)
    11791179    zsin45 = z*SQRT(0.5_wp)
    11801180
     
    14191419!OCL NOVREC
    14201420       DO ijk = 1, lot
    1421           c(ja+j) = 2.0*(a(ia+i)+a(ib+i))
    1422           c(jb+j) = 2.0*(a(ia+i)-a(ib+i))
     1421          c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
     1422          c(jb+j) = 2.0_wp*(a(ia+i)-a(ib+i))
    14231423          i = i + inc3
    14241424          j = j + inc4
     
    14491449       DO ijk = 1, lot
    14501450          c(ja+j) = a(ia+i) + a(ib+i)
    1451           c(jb+j) = (a(ia+i)-0.5*a(ib+i)) - (sin60*(b(ib+i)))
    1452           c(jc+j) = (a(ia+i)-0.5*a(ib+i)) + (sin60*(b(ib+i)))
     1451          c(jb+j) = (a(ia+i)-0.5_wp*a(ib+i)) - (sin60*(b(ib+i)))
     1452          c(jc+j) = (a(ia+i)-0.5_wp*a(ib+i)) + (sin60*(b(ib+i)))
    14531453          i = i + inc3
    14541454          j = j + inc4
     
    14811481             c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
    14821482             d(ja+j) = b(ia+i) + (b(ib+i)-b(ic+i))
    1483              c(jb+j) = c1*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
    1484                   &            b(ic+i)))) - s1*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
     1483             c(jb+j) = c1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
     1484                  &            b(ic+i))))                                      &
     1485                  &    - s1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
    14851486                  &            a(ic+i))))
    1486              d(jb+j) = s1*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
    1487                   &            b(ic+i)))) + c1*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
     1487             d(jb+j) = s1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
     1488                  &            b(ic+i))))                                      &
     1489                  &    + c1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
    14881490                  &            a(ic+i))))
    1489              c(jc+j) = c2*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
    1490                   &            b(ic+i)))) - s2*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
     1491             c(jc+j) = c2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
     1492                  &            b(ic+i))))                                      &
     1493                  &    - s2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
    14911494                  &            a(ic+i))))
    1492              d(jc+j) = s2*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
    1493                   &            b(ic+i)))) + c2*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
     1495             d(jc+j) = s2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
     1496                  &            b(ic+i))))                                      &
     1497                  &    + c2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
    14941498                  &            a(ic+i))))
    14951499             i = i + inc3
     
    15151519       DO ijk = 1, lot
    15161520          c(ja+j) = a(ia+i) + a(ib+i)
    1517           c(jb+j) = (0.5*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
    1518           c(jc+j) = -(0.5*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
     1521          c(jb+j) = (0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
     1522          c(jc+j) = -(0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
    15191523          i = i + inc3
    15201524          j = j + inc4
     
    15261530    GO TO 170
    1527153160  CONTINUE
    1528     ssin60 = 2.0*sin60
    1529     DO l = 1, la
    1530        i = ibase
    1531        j = jbase
    1532 !DIR$ IVDEP
    1533 !CDIR NODEP
    1534 !OCL NOVREC
    1535        DO ijk = 1, lot
    1536           c(ja+j) = 2.0*(a(ia+i)+a(ib+i))
    1537           c(jb+j) = (2.0*a(ia+i)-a(ib+i)) - (ssin60*b(ib+i))
    1538           c(jc+j) = (2.0*a(ia+i)-a(ib+i)) + (ssin60*b(ib+i))
     1532    ssin60 = 2.0_wp*sin60
     1533    DO l = 1, la
     1534       i = ibase
     1535       j = jbase
     1536!DIR$ IVDEP
     1537!CDIR NODEP
     1538!OCL NOVREC
     1539       DO ijk = 1, lot
     1540          c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
     1541          c(jb+j) = (2.0_wp*a(ia+i)-a(ib+i)) - (ssin60*b(ib+i))
     1542          c(jc+j) = (2.0_wp*a(ia+i)-a(ib+i)) + (ssin60*b(ib+i))
    15391543          i = i + inc3
    15401544          j = j + inc4
     
    16591663!OCL NOVREC
    16601664       DO ijk = 1, lot
    1661           c(ja+j) = 2.0*((a(ia+i)+a(ic+i))+a(ib+i))
    1662           c(jb+j) = 2.0*((a(ia+i)-a(ic+i))-b(ib+i))
    1663           c(jc+j) = 2.0*((a(ia+i)+a(ic+i))-a(ib+i))
    1664           c(jd+j) = 2.0*((a(ia+i)-a(ic+i))+b(ib+i))
     1665          c(ja+j) = 2.0_wp*((a(ia+i)+a(ic+i))+a(ib+i))
     1666          c(jb+j) = 2.0_wp*((a(ia+i)-a(ic+i))-b(ib+i))
     1667          c(jc+j) = 2.0_wp*((a(ia+i)+a(ic+i))-a(ib+i))
     1668          c(jd+j) = 2.0_wp*((a(ia+i)-a(ic+i))+b(ib+i))
    16651669          i = i + inc3
    16661670          j = j + inc4
     
    16951699       DO ijk = 1, lot
    16961700          c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
    1697           c(jb+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) - &
     1701          c(jb+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) - &
    16981702               &          (sin72*b(ib+i)+sin36*b(ic+i))
    1699           c(jc+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) - &
     1703          c(jc+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) - &
    17001704               &          (sin36*b(ib+i)-sin72*b(ic+i))
    1701           c(jd+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) + &
     1705          c(jd+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) + &
    17021706               &          (sin36*b(ib+i)-sin72*b(ic+i))
    1703           c(je+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) + &
     1707          c(je+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) + &
    17041708               &          (sin72*b(ib+i)+sin36*b(ic+i))
    17051709          i = i + inc3
     
    17401744          DO ijk = 1, lot
    17411745
    1742              a10(ijk) = (a(ia+i)-0.25*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) + &
     1746             a10(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) + &
    17431747                  &            qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
    1744              a20(ijk) = (a(ia+i)-0.25*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) - &
     1748             a20(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) - &
    17451749                  &            qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
    1746              b10(ijk) = (b(ia+i)-0.25*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) + &
     1750             b10(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) + &
    17471751                  &            qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
    1748              b20(ijk) = (b(ia+i)-0.25*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) - &
     1752             b20(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) - &
    17491753                  &            qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
    17501754             a11(ijk) = sin72*(b(ib+i)+b(ie+i)) + sin36*(b(ic+i)+b(id+i))
     
    17881792       DO ijk = 1, lot
    17891793          c(ja+j) = (a(ia+i)+a(ib+i)) + a(ic+i)
    1790           c(jb+j) = (qrt5*(a(ia+i)-a(ib+i))+(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
     1794          c(jb+j) = (qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
    17911795               &          (sin36*b(ia+i)+sin72*b(ib+i))
    1792           c(je+j) = -(qrt5*(a(ia+i)-a(ib+i))+(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
     1796          c(je+j) = -(qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
    17931797               &          (sin36*b(ia+i)+sin72*b(ib+i))
    1794           c(jc+j) = (qrt5*(a(ia+i)-a(ib+i))-(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
     1798          c(jc+j) = (qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
    17951799               &          (sin72*b(ia+i)-sin36*b(ib+i))
    1796           c(jd+j) = -(qrt5*(a(ia+i)-a(ib+i))-(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
     1800          c(jd+j) = -(qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
    17971801               &          (sin72*b(ia+i)-sin36*b(ib+i))
    17981802          i = i + inc3
     
    18151819!OCL NOVREC
    18161820       DO ijk = 1, lot
    1817           c(ja+j) = 2.0*(a(ia+i)+(a(ib+i)+a(ic+i)))
    1818           c(jb+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
     1821          c(ja+j) = 2.0_wp*(a(ia+i)+(a(ib+i)+a(ic+i)))
     1822          c(jb+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
    18191823               &          i))) - (ssin72*b(ib+i)+ssin36*b(ic+i))
    1820           c(jc+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
     1824          c(jc+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
    18211825               &          i))) - (ssin36*b(ib+i)-ssin72*b(ic+i))
    1822           c(jd+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
     1826          c(jd+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
    18231827               &          i))) + (ssin36*b(ib+i)-ssin72*b(ic+i))
    1824           c(je+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
     1828          c(je+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
    18251829               &          i))) + (ssin72*b(ib+i)+ssin36*b(ic+i))
    18261830          i = i + inc3
     
    18591863          c(ja+j) = (a(ia+i)+a(id+i)) + (a(ib+i)+a(ic+i))
    18601864          c(jd+j) = (a(ia+i)-a(id+i)) - (a(ib+i)-a(ic+i))
    1861           c(jb+j) = ((a(ia+i)-a(id+i))+0.5*(a(ib+i)-a(ic+i))) - (sin60*(b(ib+ &
     1865          c(jb+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) - (sin60*(b(ib+ &
    18621866               &          i)+b(ic+i)))
    1863           c(jf+j) = ((a(ia+i)-a(id+i))+0.5*(a(ib+i)-a(ic+i))) + (sin60*(b(ib+ &
     1867          c(jf+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) + (sin60*(b(ib+ &
    18641868               &          i)+b(ic+i)))
    1865           c(jc+j) = ((a(ia+i)+a(id+i))-0.5*(a(ib+i)+a(ic+i))) - (sin60*(b(ib+ &
     1869          c(jc+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) - (sin60*(b(ib+ &
    18661870               &          i)-b(ic+i)))
    1867           c(je+j) = ((a(ia+i)+a(id+i))-0.5*(a(ib+i)+a(ic+i))) + (sin60*(b(ib+ &
     1871          c(je+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) + (sin60*(b(ib+ &
    18681872               &          i)-b(ic+i)))
    18691873          i = i + inc3
     
    19091913
    19101914             a11(ijk) = (a(ie+i)+a(ib+i)) + (a(ic+i)+a(if+i))
    1911              a20(ijk) = (a(ia+i)+a(id+i)) - 0.5*a11(ijk)
     1915             a20(ijk) = (a(ia+i)+a(id+i)) - 0.5_wp*a11(ijk)
    19121916             a21(ijk) = sin60*((a(ie+i)+a(ib+i))-(a(ic+i)+a(if+i)))
    19131917             b11(ijk) = (b(ib+i)-b(ie+i)) + (b(ic+i)-b(if+i))
    1914              b20(ijk) = (b(ia+i)-b(id+i)) - 0.5*b11(ijk)
     1918             b20(ijk) = (b(ia+i)-b(id+i)) - 0.5_wp*b11(ijk)
    19151919             b21(ijk) = sin60*((b(ib+i)-b(ie+i))-(b(ic+i)-b(if+i)))
    19161920
     
    19241928             a11(ijk) = (a(ie+i)-a(ib+i)) + (a(ic+i)-a(if+i))
    19251929             b11(ijk) = (b(ie+i)+b(ib+i)) - (b(ic+i)+b(if+i))
    1926              a20(ijk) = (a(ia+i)-a(id+i)) - 0.5*a11(ijk)
     1930             a20(ijk) = (a(ia+i)-a(id+i)) - 0.5_wp*a11(ijk)
    19271931             a21(ijk) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
    1928              b20(ijk) = (b(ia+i)+b(id+i)) + 0.5*b11(ijk)
     1932             b20(ijk) = (b(ia+i)+b(id+i)) + 0.5_wp*b11(ijk)
    19291933             b21(ijk) = sin60*((b(ie+i)+b(ib+i))+(b(ic+i)+b(if+i)))
    19301934
     
    19641968          c(ja+j) = a(ib+i) + (a(ia+i)+a(ic+i))
    19651969          c(jd+j) = b(ib+i) - (b(ia+i)+b(ic+i))
    1966           c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5*(b(ia+i)+b(ic+i))+b(ib+i))
    1967           c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5*(b(ia+i)+b(ic+i))+b(ib+i))
    1968           c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5*(a(ia+i)+a(ic+i))-a(ib+i))
    1969           c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5*(a(ia+i)+a(ic+i))-a(ib+i))
     1970          c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i))
     1971          c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i))
     1972          c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i))
     1973          c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i))
    19701974          i = i + inc3
    19711975          j = j + inc4
     
    19771981    GO TO 170
    19781982150 CONTINUE
    1979     ssin60 = 2.0*sin60
    1980     DO l = 1, la
    1981        i = ibase
    1982        j = jbase
    1983 !DIR$ IVDEP
    1984 !CDIR NODEP
    1985 !OCL NOVREC
    1986        DO ijk = 1, lot
    1987           c(ja+j) = (2.0*(a(ia+i)+a(id+i))) + (2.0*(a(ib+i)+a(ic+i)))
    1988           c(jd+j) = (2.0*(a(ia+i)-a(id+i))) - (2.0*(a(ib+i)-a(ic+i)))
    1989           c(jb+j) = (2.0*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) - (ssin60*(b(ib+ &
     1983    ssin60 = 2.0_wp*sin60
     1984    DO l = 1, la
     1985       i = ibase
     1986       j = jbase
     1987!DIR$ IVDEP
     1988!CDIR NODEP
     1989!OCL NOVREC
     1990       DO ijk = 1, lot
     1991          c(ja+j) = (2.0_wp*(a(ia+i)+a(id+i))) + (2.0_wp*(a(ib+i)+a(ic+i)))
     1992          c(jd+j) = (2.0_wp*(a(ia+i)-a(id+i))) - (2.0_wp*(a(ib+i)-a(ic+i)))
     1993          c(jb+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) - (ssin60*(b(ib+ &
    19901994               &          i)+b(ic+i)))
    1991           c(jf+j) = (2.0*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) + (ssin60*(b(ib+ &
     1995          c(jf+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) + (ssin60*(b(ib+ &
    19921996               &          i)+b(ic+i)))
    1993           c(jc+j) = (2.0*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) - (ssin60*(b(ib+ &
     1997          c(jc+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) - (ssin60*(b(ib+ &
    19941998               &          i)-b(ic+i)))
    1995           c(je+j) = (2.0*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) + (ssin60*(b(ib+ &
     1999          c(je+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) + (ssin60*(b(ib+ &
    19962000               &          i)-b(ic+i)))
    19972001          i = i + inc3
     
    20302034!OCL NOVREC
    20312035       DO ijk = 1, lot
    2032           c(ja+j) = 2.0*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i)))
    2033           c(je+j) = 2.0*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i)))
    2034           c(jc+j) = 2.0*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i)))
    2035           c(jg+j) = 2.0*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i)))
    2036           c(jb+j) = 2.0*((a(ia+i)-a(ie+i))-b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
    2037                &          i))-(b(ib+i)+b(id+i)))
    2038           c(jf+j) = 2.0*((a(ia+i)-a(ie+i))-b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
    2039                &          i))-(b(ib+i)+b(id+i)))
    2040           c(jd+j) = 2.0*((a(ia+i)-a(ie+i))+b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
    2041                &          i))+(b(ib+i)+b(id+i)))
    2042           c(jh+j) = 2.0*((a(ia+i)-a(ie+i))+b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
    2043                &          i))+(b(ib+i)+b(id+i)))
     2036          c(ja+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i)))
     2037          c(je+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i)))
     2038          c(jc+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i)))
     2039          c(jg+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i)))
     2040          c(jb+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
     2041               &              i))-(b(ib+i)+b(id+i)))
     2042          c(jf+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
     2043               &              i))-(b(ib+i)+b(id+i)))
     2044          c(jd+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
     2045               &              i))+(b(ib+i)+b(id+i)))
     2046          c(jh+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
     2047               &              i))+(b(ib+i)+b(id+i)))
    20442048          i = i + inc3
    20452049          j = j + inc4
     
    21142118    ixxx = 1
    21152119
    2116     del = 4.0_wp*ASIN(1.0_wp)/REAL(n)
     2120    del = 4.0_wp*ASIN(1.0_wp)/REAL(n,KIND=wp)
    21172121    nil = 0
    21182122    nhl = (n/2) - 1
    21192123    DO k = nil, nhl
    2120        angle = REAL(k)*del
     2124       angle = REAL(k,KIND=wp)*del
    21212125       trigs(2*k+1) = COS(angle)
    21222126       trigs(2*k+2) = SIN(angle)
  • palm/trunk/SOURCE/time_integration.f90

    r1321 r1342  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! REAL constants defined as wp-kind
    2323!
    2424! Former revisions:
     
    205205!-- At the beginning of a simulation determine the time step as well as
    206206!-- determine and print out the run control parameters
    207     IF ( simulated_time == 0.0 )  CALL timestep
     207    IF ( simulated_time == 0.0_wp )  CALL timestep
    208208
    209209    CALL run_control
     
    221221          time_coupling = time_coupling - dt_coupling
    222222       ENDDO
    223        IF (time_coupling == 0.0 .AND. time_since_reference_point < dt_coupling)&
     223       IF (time_coupling == 0.0_wp .AND. time_since_reference_point < dt_coupling)&
    224224       THEN
    225225          time_coupling = time_since_reference_point
     
    243243!
    244244!--    Determine size of next time step
    245        IF ( simulated_time /= 0.0 )  THEN
     245       IF ( simulated_time /= 0.0_wp )  THEN
    246246          CALL timestep
    247247       ENDIF
     
    717717!
    718718!--    If required, call flow_statistics for averaging in time
    719        IF ( averaging_interval_pr /= 0.0  .AND.  &
     719       IF ( averaging_interval_pr /= 0.0_wp  .AND.  &
    720720            ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.  &
    721721            simulated_time >= skip_time_dopr )  THEN
     
    731731!
    732732!--    Sum-up 3d-arrays for later output of time-averaged 2d/3d/masked data
    733        IF ( averaging_interval /= 0.0  .AND.                                &
     733       IF ( averaging_interval /= 0.0_wp  .AND.                                &
    734734            ( dt_data_output_av - time_do_av ) <= averaging_interval  .AND. &
    735735            simulated_time >= skip_time_data_output_av )                    &
     
    745745!
    746746!--    Calculate spectra for time averaging
    747        IF ( averaging_interval_sp /= 0.0  .AND.  &
     747       IF ( averaging_interval_sp /= 0.0_wp  .AND.  &
    748748            ( dt_dosp - time_dosp ) <= averaging_interval_sp  .AND.  &
    749749            simulated_time >= skip_time_dosp )  THEN
     
    782782          IF ( dopr_n /= 0 )  CALL data_output_profiles
    783783          time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
    784           time_dopr_av = 0.0    ! due to averaging (see above)
     784          time_dopr_av = 0.0_wp    ! due to averaging (see above)
    785785       ENDIF
    786786
  • palm/trunk/SOURCE/time_to_string.f90

    r1321 r1342  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants defined as wp-kind
    2323!
    2424! Former revisions:
     
    6161!
    6262!-- Calculate the number of hours, minutes, and seconds
    63     hours     = INT( time / 3600.0 )
    64     rest_time = time - hours * 3600
    65     minutes   = INT( rest_time / 60.0 )
     63    hours     = INT( time / 3600.0_wp )
     64    rest_time = time - hours * 3600_wp
     65    minutes   = INT( rest_time / 60.0_wp )
    6666    seconds   = rest_time - minutes * 60
    6767
  • palm/trunk/SOURCE/timestep.f90

    r1323 r1342  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! REAL constants defined as wp-kind
    2323!
    2424! Former revisions:
     
    153153!--       Horizontal averages already existent, just need to average them
    154154!--       vertically.
    155           u_gtrans = 0.0
    156           v_gtrans = 0.0
     155          u_gtrans = 0.0_wp
     156          v_gtrans = 0.0_wp
    157157          DO  k = nzb+1, nzt
    158158             u_gtrans = u_gtrans + hom(k,1,1,0)
     
    164164!
    165165!--       Averaging over the entire model domain.
    166           u_gtrans_l = 0.0
    167           v_gtrans_l = 0.0
     166          u_gtrans_l = 0.0_wp
     167          v_gtrans_l = 0.0_wp
    168168          !$acc parallel present( u, v )
    169169          DO  i = nxl, nxr
     
    196196#if defined( __openacc )
    197197    IF ( dt_fixed )  THEN  ! otherwise do it further below for better cache usage
    198        u_max_l = -999999.9
    199        u_min_l =  999999.9
    200        v_max_l = -999999.9
    201        v_min_l =  999999.9
    202        w_max_l = -999999.9
    203        w_min_l =  999999.9
     198       u_max_l = -999999.9_wp
     199       u_min_l =  999999.9_wp
     200       v_max_l = -999999.9_wp
     201       v_min_l =  999999.9_wp
     202       w_max_l = -999999.9_wp
     203       w_min_l =  999999.9_wp
    204204       !$acc parallel present( u, v, w )
    205205       DO  i = nxl, nxr
     
    266266!--    Calculate the maximum time step according to the CFL-criterion,
    267267!--    individually for each velocity component
    268        dt_u_l  =  999999.9
    269        dt_v_l  =  999999.9
    270        dt_w_l  =  999999.9
    271        u_max_l = -999999.9
    272        u_min_l =  999999.9
    273        v_max_l = -999999.9
    274        v_min_l =  999999.9
    275        w_max_l = -999999.9
    276        w_min_l =  999999.9
     268       dt_u_l  =  999999.9_wp
     269       dt_v_l  =  999999.9_wp
     270       dt_w_l  =  999999.9_wp
     271       u_max_l = -999999.9_wp
     272       u_min_l =  999999.9_wp
     273       v_max_l = -999999.9_wp
     274       v_min_l =  999999.9_wp
     275       w_max_l = -999999.9_wp
     276       w_min_l =  999999.9_wp
    277277       !$acc parallel loop collapse(3) present( u, v, w )
    278278       DO  i = nxl, nxr
    279279          DO  j = nys, nyn
    280280             DO  k = nzb+1, nzt
    281                 dt_u_l  = MIN( dt_u_l, ( dx     / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10 ) ) )
    282                 dt_v_l  = MIN( dt_v_l, ( dy     / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10 ) ) )
    283                 dt_w_l  = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) )            + 1.0E-10 ) ) )
     281                dt_u_l  = MIN( dt_u_l, ( dx     / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10_wp ) ) )
     282                dt_v_l  = MIN( dt_v_l, ( dy     / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10_wp ) ) )
     283                dt_w_l  = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) )            + 1.0E-10_wp ) ) )
    284284                u_max_l = MAX( u_max_l, u(k,j,i) )
    285285                u_min_l = MIN( u_min_l, u(k,j,i) )
     
    346346!--    Calculate the maximum time step according to the CFL-criterion,
    347347!--    individually for each velocity component
    348        dt_u_l = 999999.9
    349        dt_v_l = 999999.9
    350        dt_w_l = 999999.9
     348       dt_u_l = 999999.9_wp
     349       dt_v_l = 999999.9_wp
     350       dt_w_l = 999999.9_wp
    351351       DO  i = nxl, nxr
    352352          DO  j = nys, nyn
    353353             DO  k = nzb+1, nzt
    354                 dt_u_l = MIN( dt_u_l, ( dx     / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10 ) ) )
    355                 dt_v_l = MIN( dt_v_l, ( dy     / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10 ) ) )
    356                 dt_w_l = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) )            + 1.0E-10 ) ) )
     354                dt_u_l = MIN( dt_u_l, ( dx     / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10_wp ) ) )
     355                dt_v_l = MIN( dt_v_l, ( dy     / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10_wp ) ) )
     356                dt_w_l = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) )            + 1.0E-10_wp ) ) )
    357357             ENDDO
    358358          ENDDO
     
    382382!--          in the Prandtl-layer friction term only dz/2 is used.
    383383!--          Experience from the old model seems to justify this.
    384        dt_diff_l = 999999.0
     384       dt_diff_l = 999999.0_wp
    385385
    386386       DO  k = nzb+1, nzt
    387            dxyz2_min(k) = MIN( dx2, dy2, dzw(k)*dzw(k) ) * 0.125
     387           dxyz2_min(k) = MIN( dx2, dy2, dzw(k)*dzw(k) ) * 0.125_wp
    388388       ENDDO
    389389
     
    395395             DO  k = nzb+1, nzt
    396396                dt_diff_l = MIN( dt_diff_l, dxyz2_min(k) / &
    397                                        ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20 ) )
     397                                       ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20_wp ) )
    398398             ENDDO
    399399          ENDDO
     
    414414       IF ( plant_canopy ) THEN
    415415
    416           dt_plant_canopy_l = 0.0
     416          dt_plant_canopy_l = 0.0_wp
    417417          DO  i = nxl, nxr
    418418             DO  j = nys, nyn
     
    424424                                                 v(k,j+1,i)      +  &
    425425                                                 v(k,j+1,i-1) )     &
    426                                                / 4.0 )**2        +  &
     426                                               / 4.0_wp )**2     +  &
    427427                                             ( ( w(k-1,j,i-1)    +  &
    428428                                                 w(k-1,j,i)      +  &
    429429                                                 w(k,j,i-1)      +  &
    430430                                                 w(k,j,i) )         &
    431                                                  / 4.0 )**2 )
     431                                                 / 4.0_wp )**2 )
    432432                   IF ( dt_plant_canopy_u > dt_plant_canopy_l ) THEN
    433433                      dt_plant_canopy_l = dt_plant_canopy_u 
     
    438438                                                 u(k,j,i)        +  &
    439439                                                 u(k,j,i+1) )       &
    440                                                / 4.0 )**2        +  &
     440                                               / 4.0_wp )**2     +  &
    441441                                                 v(k,j,i)**2     +  &
    442442                                             ( ( w(k-1,j-1,i)    +  &
     
    444444                                                 w(k,j-1,i)      +  &
    445445                                                 w(k,j,i) )         &
    446                                                  / 4.0 )**2 )
     446                                                 / 4.0_wp )**2 )
    447447                   IF ( dt_plant_canopy_v > dt_plant_canopy_l ) THEN
    448448                      dt_plant_canopy_l = dt_plant_canopy_v
     
    453453                                                 u(k+1,j,i)      +  &
    454454                                                 u(k+1,j,i+1) )     &
    455                                                / 4.0 )**2        +  &
     455                                               / 4.0_wp )**2     +  &
    456456                                             ( ( v(k,j,i)        +  &
    457457                                                 v(k,j+1,i)      +  &
    458458                                                 v(k+1,j,i)      +  &
    459459                                                 v(k+1,j+1,i) )     &
    460                                                / 4.0 )**2        +  &
     460                                               / 4.0_wp )**2        +  &
    461461                                                 w(k,j,i)**2 )     
    462462                   IF ( dt_plant_canopy_w > dt_plant_canopy_l ) THEN
     
    467467          ENDDO
    468468
    469           IF ( dt_plant_canopy_l > 0.0 ) THEN
     469          IF ( dt_plant_canopy_l > 0.0_wp ) THEN
    470470!
    471471!--          Invert dt_plant_canopy_l and apply a security timestep factor 0.1
    472              dt_plant_canopy_l = 0.1 / dt_plant_canopy_l
     472             dt_plant_canopy_l = 0.1_wp / dt_plant_canopy_l
    473473          ELSE
    474474!
     
    515515!
    516516!--    Set flag if the time step becomes too small.
    517        IF ( dt_3d < ( 0.00001 * dt_max ) )  THEN
     517       IF ( dt_3d < ( 0.00001_wp * dt_max ) )  THEN
    518518          stop_dt = .TRUE.
    519519
     
    553553!
    554554!--    Ensure a smooth value (two significant digits) of the timestep.
    555        div = 1000.0
     555       div = 1000.0_wp
    556556       DO  WHILE ( dt_3d < div )
    557           div = div / 10.0
     557          div = div / 10.0_wp
    558558       ENDDO
    559        dt_3d = NINT( dt_3d * 100.0 / div ) * div / 100.0
     559       dt_3d = NINT( dt_3d * 100.0_wp / div ) * div / 100.0_wp
    560560
    561561!
  • palm/trunk/SOURCE/tridia_solver.f90

    r1323 r1342  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! REAL constants defined as wp-kind
    2323!
    2424! Former revisions:
     
    110110          ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1)
    111111          ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1)
    112           ddzuw(k,3) = -1.0 * &
     112          ddzuw(k,3) = -1.0_wp * &
    113113                       ( ddzu_pres(k+2) * ddzw(k+1) + ddzu_pres(k+1) * ddzw(k+1) )
    114114       ENDDO
     
    172172                IF ( j >= 0  .AND.  j <= nnyh )  THEN
    173173                   IF ( i >= 0  .AND.  i <= nnxh )  THEN
    174                       ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
    175                                                   REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
    176                                 2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
    177                                                   REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
     174                      ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
     175                                            REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
     176                                2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
     177                                            REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
    178178                   ELSE
    179                       ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
    180                                                   REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
    181                                 2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
    182                                                   REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
     179                      ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
     180                                            REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
     181                                2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
     182                                            REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
    183183                   ENDIF
    184184                ELSE
    185185                   IF ( i >= 0  .AND.  i <= nnxh )  THEN
    186                       ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
    187                                                   REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
    188                                 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( ny+1-j ) ) / &
    189                                                   REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
     186                      ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
     187                                            REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
     188                                2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / &
     189                                            REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
    190190                   ELSE
    191                       ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
    192                                                   REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
    193                                 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( ny+1-j ) ) / &
    194                                                   REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
     191                      ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
     192                                            REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
     193                                2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / &
     194                                            REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
    195195                   ENDIF
    196196                ENDIF
     
    283283
    284284                   IF ( k == nz-1 )  THEN
    285                       ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20 )
     285                      ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20_wp )
    286286                   ELSE
    287287                      ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
     
    301301                !$acc kernels loop present( ar )
    302302                DO  k = 1, nz
    303                    ar(nxl_z,nys_z,k) = 0.0
     303                   ar(nxl_z,nys_z,k) = 0.0_wp
    304304                ENDDO
    305305                !$acc end kernels loop
     
    367367
    368368                   IF ( k == nz-1 )  THEN
    369                       ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20 )
     369                      ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20_wp )
    370370                   ELSE
    371371                      ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
     
    385385                !$acc kernels loop present( ar )
    386386                DO  k = 1, nz
    387                    ar(nxl_z,nys_z,k) = 0.0
     387                   ar(nxl_z,nys_z,k) = 0.0_wp
    388388                ENDDO
    389389             ENDIF
     
    560560          DO  i = 0, nx
    561561             IF ( i >= 0 .AND. i <= nnxh ) THEN
    562                 l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
    563                                          REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
    564                        2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
    565                                          REAL( ny+1, KIND=wp ) ) ) * ddy2
     562                l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
     563                                   REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
     564                       2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
     565                                   REAL( ny+1, KIND=wp ) ) ) * ddy2
    566566             ELSE
    567                 l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
    568                                          REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
    569                        2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
    570                                          REAL( ny+1, KIND=wp ) ) ) * ddy2
     567                l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
     568                                   REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
     569                       2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
     570                                   REAL( ny+1, KIND=wp ) ) ) * ddy2
    571571             ENDIF
    572572          ENDDO
     
    574574          DO  k = 0, nz-1
    575575             DO  i = 0, nx
    576                 a = -1.0 * ddzu_pres(k+2) * ddzw(k+1)
    577                 c = -1.0 * ddzu_pres(k+1) * ddzw(k+1)
     576                a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1)
     577                c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1)
    578578                tri_for_1d(1,i,k) = a + c - l(i)
    579579             ENDDO
     
    660660!--       the model domain.
    661661          DO  i = 0, nx
    662              ar(i,nz) = ar1(i,nz-1) / ( tri_for_1d(4,i,nz-1) + 1.0E-20 )
     662             ar(i,nz) = ar1(i,nz-1) / ( tri_for_1d(4,i,nz-1) + 1.0E-20_wp )
    663663          ENDDO
    664664          DO  k = nz-2, 0, -1
     
    676676             IF ( j == 0 )  THEN
    677677                DO  k = 1, nz
    678                    ar(0,k) = 0.0
     678                   ar(0,k) = 0.0_wp
    679679                ENDDO
    680680             ENDIF
Note: See TracChangeset for help on using the changeset viewer.