Ignore:
Timestamp:
Jan 7, 2015 7:12:25 PM (7 years ago)
Author:
hoffmann
Message:

bug-fix: Bott-Chlond scheme

File:
1 edited

Legend:

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

    r1375 r1517  
    1  SUBROUTINE advec_s_bc( sk, sk_char )
     1MODULE advec_s_bc_mod
    22
    33!--------------------------------------------------------------------------------!
     
    2020! Current revisions:
    2121! -----------------
    22 !
     22! interface added to advec_s_bc
    2323!
    2424! Former revisions:
     
    7676!------------------------------------------------------------------------------!
    7777
    78     USE advection,                                                             &
    79         ONLY:  aex, bex, dex, eex
    80 
    81     USE arrays_3d,                                                             &
    82         ONLY:  d, ddzw, dzu, dzw, tend, u, v, w
    83 
    84     USE control_parameters,                                                    &
    85         ONLY:  dt_3d, bc_pt_t_val, bc_q_t_val, ibc_pt_b, ibc_pt_t, ibc_q_t,    &
    86                message_string, pt_slope_offset, sloping_surface, u_gtrans,     &
    87                v_gtrans
    88 
    89     USE cpulog,                                                                &
    90         ONLY:  cpu_log, log_point_s
    91 
    92     USE grid_variables,                                                        &
    93         ONLY:  ddx, ddy
    94 
    95     USE indices,                                                               &
    96         ONLY:  nx, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
    97 
    98     USE kinds
    99 
    100     USE pegrid
    101 
    102     USE statistics,                                                            &
    103         ONLY:  rmask, statistic_regions, sums_wsts_bc_l
    104 
    105 
    106     IMPLICIT NONE
    107 
    108     CHARACTER (LEN=*) ::  sk_char !:
    109 
    110     INTEGER(iwp) ::  i         !:
    111     INTEGER(iwp) ::  ix        !:
    112     INTEGER(iwp) ::  j         !:
    113     INTEGER(iwp) ::  k         !:
    114     INTEGER(iwp) ::  ngp       !:
    115     INTEGER(iwp) ::  sr        !:
    116     INTEGER(iwp) ::  type_xz_2 !:
    117 
    118     REAL(wp) ::  cim    !:
    119     REAL(wp) ::  cimf   !:
    120     REAL(wp) ::  cip    !:
    121     REAL(wp) ::  cipf   !:
    122     REAL(wp) ::  d_new  !:
    123     REAL(wp) ::  denomi !: denominator
    124     REAL(wp) ::  fminus !:
    125     REAL(wp) ::  fplus  !:
    126     REAL(wp) ::  f2     !:
    127     REAL(wp) ::  f4     !:
    128     REAL(wp) ::  f8     !:
    129     REAL(wp) ::  f12    !:
    130     REAL(wp) ::  f24    !:
    131     REAL(wp) ::  f48    !:
    132     REAL(wp) ::  f1920  !:
    133     REAL(wp) ::  im     !:
    134     REAL(wp) ::  ip     !:
    135     REAL(wp) ::  m2     !:
    136     REAL(wp) ::  m3     !:
    137     REAL(wp) ::  numera !: numerator
    138     REAL(wp) ::  snenn  !:
    139     REAL(wp) ::  sterm  !:
    140     REAL(wp) ::  tendcy !:
    141     REAL(wp) ::  t1     !:
    142     REAL(wp) ::  t2     !:
    143 
    144     REAL(wp) ::  fmax(2)   !:
    145     REAL(wp) ::  fmax_l(2) !:
    146    
     78    PRIVATE
     79    PUBLIC advec_s_bc
     80
     81    INTERFACE advec_s_bc
     82       MODULE PROCEDURE advec_s_bc
     83    END INTERFACE advec_s_bc
     84
     85 CONTAINS
     86
     87    SUBROUTINE advec_s_bc( sk, sk_char )
     88
     89       USE advection,                                                             &
     90           ONLY:  aex, bex, dex, eex
     91
     92       USE arrays_3d,                                                             &
     93           ONLY:  d, ddzw, dzu, dzw, tend, u, v, w
     94
     95       USE control_parameters,                                                    &
     96           ONLY:  dt_3d, bc_pt_t_val, bc_q_t_val, ibc_pt_b, ibc_pt_t, ibc_q_t,    &
     97                  message_string, pt_slope_offset, sloping_surface, u_gtrans,     &
     98                  v_gtrans
     99
     100       USE cpulog,                                                                &
     101           ONLY:  cpu_log, log_point_s
     102
     103       USE grid_variables,                                                        &
     104           ONLY:  ddx, ddy
     105
     106       USE indices,                                                               &
     107           ONLY:  nx, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
     108
     109       USE kinds
     110
     111       USE pegrid
     112
     113       USE statistics,                                                            &
     114           ONLY:  rmask, statistic_regions, sums_wsts_bc_l
     115
     116
     117       IMPLICIT NONE
     118
     119       CHARACTER (LEN=*) ::  sk_char !:
     120
     121       INTEGER(iwp) ::  i         !:
     122       INTEGER(iwp) ::  ix        !:
     123       INTEGER(iwp) ::  j         !:
     124       INTEGER(iwp) ::  k         !:
     125       INTEGER(iwp) ::  ngp       !:
     126       INTEGER(iwp) ::  sr        !:
     127       INTEGER(iwp) ::  type_xz_2 !:
     128
     129       REAL(wp) ::  cim    !:
     130       REAL(wp) ::  cimf   !:
     131       REAL(wp) ::  cip    !:
     132       REAL(wp) ::  cipf   !:
     133       REAL(wp) ::  d_new  !:
     134       REAL(wp) ::  denomi !: denominator
     135       REAL(wp) ::  fminus !:
     136       REAL(wp) ::  fplus  !:
     137       REAL(wp) ::  f2     !:
     138       REAL(wp) ::  f4     !:
     139       REAL(wp) ::  f8     !:
     140       REAL(wp) ::  f12    !:
     141       REAL(wp) ::  f24    !:
     142       REAL(wp) ::  f48    !:
     143       REAL(wp) ::  f1920  !:
     144       REAL(wp) ::  im     !:
     145       REAL(wp) ::  ip     !:
     146       REAL(wp) ::  m2     !:
     147       REAL(wp) ::  m3     !:
     148       REAL(wp) ::  numera !: numerator
     149       REAL(wp) ::  snenn  !:
     150       REAL(wp) ::  sterm  !:
     151       REAL(wp) ::  tendcy !:
     152       REAL(wp) ::  t1     !:
     153       REAL(wp) ::  t2     !:
     154
     155       REAL(wp) ::  fmax(2)   !:
     156       REAL(wp) ::  fmax_l(2) !:
     157       
    147158#if defined( __nopointer )
    148     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !:
     159       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !:
    149160#else
    150     REAL(wp), DIMENSION(:,:,:), POINTER ::  sk
     161       REAL(wp), DIMENSION(:,:,:), POINTER ::  sk
    151162#endif
    152163
    153     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  a0   !:
    154     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  a1   !:
    155     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  a12  !:
    156     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  a2   !:
    157     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  a22  !:
    158     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  immb !:
    159     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  imme !:
    160     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  impb !:
    161     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  impe !:
    162     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ipmb !:
    163     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ipme !:
    164     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ippb !:
    165     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ippe !:
    166    
    167     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  sk_p !:
     164       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  a0   !:
     165       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  a1   !:
     166       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  a12  !:
     167       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  a2   !:
     168       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  a22  !:
     169       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  immb !:
     170       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  imme !:
     171       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  impb !:
     172       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  impe !:
     173       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ipmb !:
     174       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ipme !:
     175       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ippb !:
     176       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ippe !:
     177       
     178       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  sk_p !:
    168179
    169180#if defined( __nec )
    170     REAL(sp) ::  m1n, m1z  !Wichtig: Division !:
    171     REAL(sp), DIMENSION(:,:), ALLOCATABLE :: m1, sw !:
     181       REAL(sp) ::  m1n, m1z  !Wichtig: Division !:
     182       REAL(sp), DIMENSION(:,:), ALLOCATABLE :: m1, sw !:
    172183#else
    173     REAL(wp) ::  m1n, m1z
    174     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: m1, sw
     184       REAL(wp) ::  m1n, m1z
     185       REAL(wp), DIMENSION(:,:), ALLOCATABLE :: m1, sw
    175186#endif
    176187
    177188
    178189!
    179 !-- Array sk_p requires 2 extra elements for each dimension
    180     ALLOCATE( sk_p(nzb-2:nzt+3,nys-3:nyn+3,nxl-3:nxr+3) )
    181     sk_p = 0.0_wp
    182 
    183 !
    184 !-- Assign reciprocal values in order to avoid divisions later
    185     f2    = 0.5_wp
    186     f4    = 0.25_wp
    187     f8    = 0.125_wp
    188     f12   = 0.8333333333333333E-01_wp
    189     f24   = 0.4166666666666666E-01_wp
    190     f48   = 0.2083333333333333E-01_wp
    191     f1920 = 0.5208333333333333E-03_wp
    192 
    193 !
    194 !-- Advection in x-direction:
    195 
    196 !
    197 !-- Save the quantity to be advected in a local array
    198 !-- add an enlarged boundary in x-direction
    199     DO  i = nxl-1, nxr+1
     190!--    Array sk_p requires 2 extra elements for each dimension
     191       ALLOCATE( sk_p(nzb-2:nzt+3,nys-3:nyn+3,nxl-3:nxr+3) )
     192       sk_p = 0.0_wp
     193
     194!
     195!--    Assign reciprocal values in order to avoid divisions later
     196       f2    = 0.5_wp
     197       f4    = 0.25_wp
     198       f8    = 0.125_wp
     199       f12   = 0.8333333333333333E-01_wp
     200       f24   = 0.4166666666666666E-01_wp
     201       f48   = 0.2083333333333333E-01_wp
     202       f1920 = 0.5208333333333333E-03_wp
     203
     204!
     205!--    Advection in x-direction:
     206
     207!
     208!--    Save the quantity to be advected in a local array
     209!--    add an enlarged boundary in x-direction
     210       DO  i = nxl-1, nxr+1
     211          DO  j = nys, nyn
     212             DO  k = nzb, nzt+1
     213                sk_p(k,j,i) = sk(k,j,i)
     214             ENDDO
     215          ENDDO
     216       ENDDO
     217#if defined( __parallel )
     218       ngp = 2 * ( nzt - nzb + 6 ) * ( nyn - nys + 7 )
     219       CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'start' )
     220!
     221!--    Send left boundary, receive right boundary
     222       CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxl+1), ngp, MPI_REAL, pleft,  0,      &
     223                          sk_p(nzb-2,nys-3,nxr+2), ngp, MPI_REAL, pright, 0,      &
     224                          comm2d, status, ierr )
     225!
     226!--    Send right boundary, receive left boundary
     227       CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxr-2), ngp, MPI_REAL, pright, 1,      &
     228                          sk_p(nzb-2,nys-3,nxl-3), ngp, MPI_REAL, pleft,  1,      &
     229                          comm2d, status, ierr )
     230       CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' )
     231#else
     232
     233!
     234!--    Cyclic boundary conditions
     235       sk_p(:,nys:nyn,nxl-3) = sk_p(:,nys:nyn,nxr-2)
     236       sk_p(:,nys:nyn,nxl-2) = sk_p(:,nys:nyn,nxr-1)
     237       sk_p(:,nys:nyn,nxr+2) = sk_p(:,nys:nyn,nxl+1)
     238       sk_p(:,nys:nyn,nxr+3) = sk_p(:,nys:nyn,nxl+2)
     239#endif
     240
     241!
     242!--    In case of a sloping surface, the additional gridpoints in x-direction
     243!--    of the temperature field at the left and right boundary of the total
     244!--    domain must be adjusted by the temperature difference between this distance
     245       IF ( sloping_surface  .AND.  sk_char == 'pt' )  THEN
     246          IF ( nxl ==  0 )  THEN
     247             sk_p(:,nys:nyn,nxl-3) = sk_p(:,nys:nyn,nxl-3) - pt_slope_offset
     248             sk_p(:,nys:nyn,nxl-2) = sk_p(:,nys:nyn,nxl-2) - pt_slope_offset
     249          ENDIF
     250          IF ( nxr == nx )  THEN
     251             sk_p(:,nys:nyn,nxr+2) = sk_p(:,nys:nyn,nxr+2) + pt_slope_offset
     252             sk_p(:,nys:nyn,nxr+3) = sk_p(:,nys:nyn,nxr+3) + pt_slope_offset
     253          ENDIF
     254       ENDIF
     255
     256!
     257!--    Initialise control density
     258       d = 0.0_wp
     259
     260!
     261!--    Determine maxima of the first and second derivative in x-direction
     262       fmax_l = 0.0_wp
     263       DO  i = nxl, nxr
     264          DO  j = nys, nyn
     265             DO  k = nzb+1, nzt
     266                numera = ABS( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) )
     267                denomi  = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
     268                fmax_l(1) = MAX( fmax_l(1) , numera )
     269                fmax_l(2) = MAX( fmax_l(2) , denomi  )
     270             ENDDO
     271          ENDDO
     272       ENDDO
     273#if defined( __parallel )
     274       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     275       CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr )
     276#else
     277       fmax = fmax_l
     278#endif
     279
     280       fmax = 0.04_wp * fmax
     281
     282!
     283!--    Allocate temporary arrays
     284       ALLOCATE( a0(nzb+1:nzt,nxl-1:nxr+1),   a1(nzb+1:nzt,nxl-1:nxr+1),          &
     285                 a2(nzb+1:nzt,nxl-1:nxr+1),   a12(nzb+1:nzt,nxl-1:nxr+1),         &
     286                 a22(nzb+1:nzt,nxl-1:nxr+1),  immb(nzb+1:nzt,nxl-1:nxr+1),        &
     287                 imme(nzb+1:nzt,nxl-1:nxr+1), impb(nzb+1:nzt,nxl-1:nxr+1),        &
     288                 impe(nzb+1:nzt,nxl-1:nxr+1), ipmb(nzb+1:nzt,nxl-1:nxr+1),        &
     289                 ipme(nzb+1:nzt,nxl-1:nxr+1), ippb(nzb+1:nzt,nxl-1:nxr+1),        &
     290                 ippe(nzb+1:nzt,nxl-1:nxr+1), m1(nzb+1:nzt,nxl-2:nxr+2),          &
     291                 sw(nzb+1:nzt,nxl-1:nxr+1)                                        &
     292               )
     293       imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp
     294
     295!
     296!--    Initialise point of time measuring of the exponential portion (this would
     297!--    not work if done locally within the loop)
     298       CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'start' )
     299       CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
     300
     301!
     302!--    Outer loop of all j
    200303       DO  j = nys, nyn
    201           DO  k = nzb, nzt+1
    202              sk_p(k,j,i) = sk(k,j,i)
    203           ENDDO
    204        ENDDO
    205     ENDDO
     304
     305!
     306!--       Compute polynomial coefficients
     307          DO  i = nxl-1, nxr+1
     308             DO  k = nzb+1, nzt
     309                a12(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
     310                a22(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i)        &
     311                                                    + sk_p(k,j,i-1) )
     312                a0(k,i) = ( 9.0_wp * sk_p(k,j,i+2)    - 116.0_wp * sk_p(k,j,i+1)  &
     313                            + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j,i-1)  &
     314                            + 9.0_wp * sk_p(k,j,i-2) ) * f1920
     315                a1(k,i) = ( -5.0_wp * sk_p(k,j,i+2)   + 34.0_wp * sk_p(k,j,i+1)   &
     316                            - 34.0_wp * sk_p(k,j,i-1) + 5.0_wp * sk_p(k,j,i-2)    &
     317                          ) * f48
     318                a2(k,i) = ( -3.0_wp * sk_p(k,j,i+2) + 36.0_wp * sk_p(k,j,i+1)     &
     319                            - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j,i-1)     &
     320                            - 3.0_wp * sk_p(k,j,i-2) ) * f48
     321             ENDDO
     322          ENDDO
     323
     324!
     325!--       Fluxes using the Bott scheme
     326!--       *VOCL LOOP,UNROLL(2)
     327          DO  i = nxl, nxr
     328             DO  k = nzb+1, nzt
     329                cip  =  MAX( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
     330                cim  = -MIN( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
     331                cipf = 1.0_wp - 2.0_wp * cip
     332                cimf = 1.0_wp - 2.0_wp * cim
     333                ip   =   a0(k,i)   * f2  * ( 1.0_wp - cipf )                      &
     334                       + a1(k,i)   * f8  * ( 1.0_wp - cipf*cipf )                 &
     335                       + a2(k,i)   * f24 * ( 1.0_wp - cipf*cipf*cipf )
     336                im   =   a0(k,i+1) * f2  * ( 1.0_wp - cimf )                      &
     337                       - a1(k,i+1) * f8  * ( 1.0_wp - cimf*cimf )                 &
     338                       + a2(k,i+1) * f24 * ( 1.0_wp - cimf*cimf*cimf )
     339                ip   = MAX( ip, 0.0_wp )
     340                im   = MAX( im, 0.0_wp )
     341                ippb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
     342                impb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i+1) / (ip+im+1E-15_wp) )
     343
     344                cip  =  MAX( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
     345                cim  = -MIN( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
     346                cipf = 1.0_wp - 2.0_wp * cip
     347                cimf = 1.0_wp - 2.0_wp * cim
     348                ip   =   a0(k,i-1) * f2  * ( 1.0_wp - cipf )                      &
     349                       + a1(k,i-1) * f8  * ( 1.0_wp - cipf*cipf )                 &
     350                       + a2(k,i-1) * f24 * ( 1.0_wp - cipf*cipf*cipf )
     351                im   =   a0(k,i)   * f2  * ( 1.0_wp - cimf )                      &
     352                       - a1(k,i)   * f8  * ( 1.0_wp - cimf*cimf )                 &
     353                       + a2(k,i)   * f24 * ( 1.0_wp - cimf*cimf*cimf )
     354                ip   = MAX( ip, 0.0_wp )
     355                im   = MAX( im, 0.0_wp )
     356                ipmb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i-1) / (ip+im+1E-15_wp) )
     357                immb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
     358             ENDDO
     359          ENDDO
     360
     361!
     362!--       Compute monitor function m1
     363          DO  i = nxl-2, nxr+2
     364             DO  k = nzb+1, nzt
     365                m1z = ABS( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) )
     366                m1n = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
     367                IF ( m1n /= 0.0_wp  .AND.  m1n >= m1z )  THEN
     368                   m1(k,i) = m1z / m1n
     369                   IF ( m1(k,i) /= 2.0_wp  .AND.  m1n < fmax(2) )  m1(k,i) = 0.0_wp
     370                ELSEIF ( m1n < m1z )  THEN
     371                   m1(k,i) = -1.0_wp
     372                ELSE
     373                   m1(k,i) = 0.0_wp
     374                ENDIF
     375             ENDDO
     376          ENDDO
     377
     378!
     379!--       Compute switch sw
     380          sw = 0.0_wp
     381          DO  i = nxl-1, nxr+1
     382             DO  k = nzb+1, nzt
     383                m2 = 2.0_wp * ABS( a1(k,i) - a12(k,i) ) /                         &
     384                     MAX( ABS( a1(k,i) + a12(k,i) ), 1E-35_wp )
     385                IF ( ABS( a1(k,i) + a12(k,i) ) < fmax(2) )  m2 = 0.0_wp
     386
     387                m3 = 2.0_wp * ABS( a2(k,i) - a22(k,i) ) /                         &
     388                     MAX( ABS( a2(k,i) + a22(k,i) ), 1E-35_wp )
     389                IF ( ABS( a2(k,i) + a22(k,i) ) < fmax(1) )  m3 = 0.0_wp
     390
     391                t1 = 0.35_wp
     392                t2 = 0.35_wp
     393                IF ( m1(k,i) == -1.0_wp )  t2 = 0.12_wp
     394
     395!--             *VOCL STMT,IF(10)
     396                IF ( m1(k,i-1) == 1.0_wp .OR. m1(k,i) == 1.0_wp                   &
     397                     .OR. m1(k,i+1) == 1.0_wp .OR.  m2 > t2  .OR.  m3 > t2  .OR.  &
     398                     ( m1(k,i) > t1  .AND.  m1(k,i-1) /= -1.0_wp  .AND.           &
     399                       m1(k,i) /= -1.0_wp  .AND.  m1(k,i+1) /= -1.0_wp )          &
     400                   )  sw(k,i) = 1.0_wp
     401             ENDDO
     402          ENDDO
     403
     404!
     405!--       Fluxes using the exponential scheme
     406          CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
     407          DO  i = nxl, nxr
     408             DO  k = nzb+1, nzt
     409
     410!--             *VOCL STMT,IF(10)
     411                IF ( sw(k,i) == 1.0_wp )  THEN
     412                   snenn = sk_p(k,j,i+1) - sk_p(k,j,i-1)
     413                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
     414                   sterm = ( sk_p(k,j,i) - sk_p(k,j,i-1) ) / snenn
     415                   sterm = MIN( sterm, 0.9999_wp )
     416                   sterm = MAX( sterm, 0.0001_wp )
     417
     418                   ix = INT( sterm * 1000 ) + 1
     419
     420                   cip =  MAX( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
     421
     422                   ippe(k,i) = sk_p(k,j,i-1) * cip + snenn * (                    &
     423                               aex(ix) * cip + bex(ix) / dex(ix) * (              &
     424                               eex(ix) -                                          &
     425                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
     426                                                                   )              &
     427                                                             )
     428                   IF ( sterm == 0.0001_wp )  ippe(k,i) = sk_p(k,j,i) * cip
     429                   IF ( sterm == 0.9999_wp )  ippe(k,i) = sk_p(k,j,i) * cip
     430
     431                   snenn = sk_p(k,j,i-1) - sk_p(k,j,i+1)
     432                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
     433                   sterm = ( sk_p(k,j,i) - sk_p(k,j,i+1) ) / snenn
     434                   sterm = MIN( sterm, 0.9999_wp )
     435                   sterm = MAX( sterm, 0.0001_wp )
     436
     437                   ix = INT( sterm * 1000 ) + 1
     438
     439                   cim = -MIN( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
     440
     441                   imme(k,i) = sk_p(k,j,i+1) * cim + snenn * (                    &
     442                               aex(ix) * cim + bex(ix) / dex(ix) * (              &
     443                               eex(ix) -                                          &
     444                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
     445                                                                   )              &
     446                                                             )
     447                   IF ( sterm == 0.0001_wp )  imme(k,i) = sk_p(k,j,i) * cim
     448                   IF ( sterm == 0.9999_wp )  imme(k,i) = sk_p(k,j,i) * cim
     449                ENDIF
     450
     451!--             *VOCL STMT,IF(10)
     452                IF ( sw(k,i+1) == 1.0_wp )  THEN
     453                   snenn = sk_p(k,j,i) - sk_p(k,j,i+2)
     454                   IF ( ABS( snenn ) .LT. 1E-9_wp )  snenn = 1E-9_wp
     455                   sterm = ( sk_p(k,j,i+1) - sk_p(k,j,i+2) ) / snenn
     456                   sterm = MIN( sterm, 0.9999_wp )
     457                   sterm = MAX( sterm, 0.0001_wp )
     458
     459                   ix = INT( sterm * 1000 ) + 1
     460
     461                   cim = -MIN( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
     462
     463                   impe(k,i) = sk_p(k,j,i+2) * cim + snenn * (                    &
     464                               aex(ix) * cim + bex(ix) / dex(ix) * (              &
     465                               eex(ix) -                                          &
     466                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
     467                                                                   )              &
     468                                                             )
     469                   IF ( sterm == 0.0001_wp )  impe(k,i) = sk_p(k,j,i+1) * cim
     470                   IF ( sterm == 0.9999_wp )  impe(k,i) = sk_p(k,j,i+1) * cim
     471                ENDIF
     472
     473!--             *VOCL STMT,IF(10)
     474                IF ( sw(k,i-1) == 1.0_wp )  THEN
     475                   snenn = sk_p(k,j,i) - sk_p(k,j,i-2)
     476                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
     477                   sterm = ( sk_p(k,j,i-1) - sk_p(k,j,i-2) ) / snenn
     478                   sterm = MIN( sterm, 0.9999_wp )
     479                   sterm = MAX( sterm, 0.0001_wp )
     480
     481                   ix = INT( sterm * 1000 ) + 1
     482
     483                   cip = MAX( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
     484
     485                   ipme(k,i) = sk_p(k,j,i-2) * cip + snenn * (                    &
     486                               aex(ix) * cip + bex(ix) / dex(ix) * (              &
     487                               eex(ix) -                                          &
     488                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
     489                                                                   )              &
     490                                                             )
     491                   IF ( sterm == 0.0001_wp )  ipme(k,i) = sk_p(k,j,i-1) * cip
     492                   IF ( sterm == 0.9999_wp )  ipme(k,i) = sk_p(k,j,i-1) * cip
     493                ENDIF
     494
     495             ENDDO
     496          ENDDO
     497          CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
     498
     499!
     500!--       Prognostic equation
     501          DO  i = nxl, nxr
     502             DO  k = nzb+1, nzt
     503                fplus  = ( 1.0_wp - sw(k,i)   ) * ippb(k,i) + sw(k,i)   * ippe(k,i)  &
     504                       - ( 1.0_wp - sw(k,i+1) ) * impb(k,i) - sw(k,i+1) * impe(k,i)
     505                fminus = ( 1.0_wp - sw(k,i-1) ) * ipmb(k,i) + sw(k,i-1) * ipme(k,i)  &
     506                       - ( 1.0_wp - sw(k,i)   ) * immb(k,i) - sw(k,i)   * imme(k,i)
     507                tendcy = fplus - fminus
     508!
     509!--             Removed in order to optimize speed
     510!                ffmax   = MAX( ABS( fplus ), ABS( fminus ), 1E-35_wp )
     511!                IF ( ( ABS( tendcy ) / ffmax ) < 1E-7_wp )  tendcy = 0.0
     512!
     513!--             Density correction because of possible remaining divergences
     514                d_new = d(k,j,i) - ( u(k,j,i+1) - u(k,j,i) ) * dt_3d * ddx
     515                sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /    &
     516                              ( 1.0_wp + d_new )
     517                d(k,j,i)  = d_new
     518             ENDDO
     519          ENDDO
     520
     521       ENDDO   ! End of the advection in x-direction
     522
     523!
     524!--    Deallocate temporary arrays
     525       DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme,      &
     526                   ippb, ippe, m1, sw )
     527
     528
     529!
     530!--    Enlarge boundary of local array cyclically in y-direction
    206531#if defined( __parallel )
    207     ngp = 2 * ( nzt - nzb + 6 ) * ( nyn - nys + 7 )
    208     CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'start' )
    209 !
    210 !-- Send left boundary, receive right boundary
    211     CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxl+1), ngp, MPI_REAL, pleft,  0,      &
    212                        sk_p(nzb-2,nys-3,nxr+2), ngp, MPI_REAL, pright, 0,      &
    213                        comm2d, status, ierr )
    214 !
    215 !-- Send right boundary, receive left boundary
    216     CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxr-2), ngp, MPI_REAL, pright, 1,      &
    217                        sk_p(nzb-2,nys-3,nxl-3), ngp, MPI_REAL, pleft,  1,      &
    218                        comm2d, status, ierr )
    219     CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' )
     532       ngp = ( nzt - nzb + 6 ) * ( nyn - nys + 7 )
     533       CALL MPI_TYPE_VECTOR( nxr-nxl+7, 3*(nzt-nzb+6), ngp, MPI_REAL,             &
     534                             type_xz_2, ierr )
     535       CALL MPI_TYPE_COMMIT( type_xz_2, ierr )
     536!
     537!--    Send front boundary, receive rear boundary
     538       CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' )
     539       CALL MPI_SENDRECV( sk_p(nzb-2,nys,nxl-3),   1, type_xz_2, psouth, 0,       &
     540                          sk_p(nzb-2,nyn+1,nxl-3), 1, type_xz_2, pnorth, 0,       &
     541                          comm2d, status, ierr )
     542!
     543!--    Send rear boundary, receive front boundary
     544       CALL MPI_SENDRECV( sk_p(nzb-2,nyn-2,nxl-3), 1, type_xz_2, pnorth, 1,       &
     545                          sk_p(nzb-2,nys-3,nxl-3), 1, type_xz_2, psouth, 1,       &
     546                          comm2d, status, ierr )
     547       CALL MPI_TYPE_FREE( type_xz_2, ierr )
     548       CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' )
    220549#else
    221 
    222 !
    223 !-- Cyclic boundary conditions
    224     sk_p(:,nys:nyn,nxl-3) = sk_p(:,nys:nyn,nxr-2)
    225     sk_p(:,nys:nyn,nxl-2) = sk_p(:,nys:nyn,nxr-1)
    226     sk_p(:,nys:nyn,nxr+2) = sk_p(:,nys:nyn,nxl+1)
    227     sk_p(:,nys:nyn,nxr+3) = sk_p(:,nys:nyn,nxl+2)
    228 #endif
    229 
    230 !
    231 !-- In case of a sloping surface, the additional gridpoints in x-direction
    232 !-- of the temperature field at the left and right boundary of the total
    233 !-- domain must be adjusted by the temperature difference between this distance
    234     IF ( sloping_surface  .AND.  sk_char == 'pt' )  THEN
    235        IF ( nxl ==  0 )  THEN
    236           sk_p(:,nys:nyn,nxl-3) = sk_p(:,nys:nyn,nxl-3) - pt_slope_offset
    237           sk_p(:,nys:nyn,nxl-2) = sk_p(:,nys:nyn,nxl-2) - pt_slope_offset
    238        ENDIF
    239        IF ( nxr == nx )  THEN
    240           sk_p(:,nys:nyn,nxr+2) = sk_p(:,nys:nyn,nxr+2) + pt_slope_offset
    241           sk_p(:,nys:nyn,nxr+3) = sk_p(:,nys:nyn,nxr+3) + pt_slope_offset
    242        ENDIF
    243     ENDIF
    244 
    245 !
    246 !-- Initialise control density
    247     d = 0.0_wp
    248 
    249 !
    250 !-- Determine maxima of the first and second derivative in x-direction
    251     fmax_l = 0.0_wp
    252     DO  i = nxl, nxr
    253        DO  j = nys, nyn
    254           DO  k = nzb+1, nzt
    255              numera = ABS( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) )
    256              denomi  = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
    257              fmax_l(1) = MAX( fmax_l(1) , numera )
    258              fmax_l(2) = MAX( fmax_l(2) , denomi  )
    259           ENDDO
    260        ENDDO
    261     ENDDO
    262 #if defined( __parallel )
    263     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    264     CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr )
    265 #else
    266     fmax = fmax_l
    267 #endif
    268 
    269     fmax = 0.04_wp * fmax
    270 
    271 !
    272 !-- Allocate temporary arrays
    273     ALLOCATE( a0(nzb+1:nzt,nxl-1:nxr+1),   a1(nzb+1:nzt,nxl-1:nxr+1),          &
    274               a2(nzb+1:nzt,nxl-1:nxr+1),   a12(nzb+1:nzt,nxl-1:nxr+1),         &
    275               a22(nzb+1:nzt,nxl-1:nxr+1),  immb(nzb+1:nzt,nxl-1:nxr+1),        &
    276               imme(nzb+1:nzt,nxl-1:nxr+1), impb(nzb+1:nzt,nxl-1:nxr+1),        &
    277               impe(nzb+1:nzt,nxl-1:nxr+1), ipmb(nzb+1:nzt,nxl-1:nxr+1),        &
    278               ipme(nzb+1:nzt,nxl-1:nxr+1), ippb(nzb+1:nzt,nxl-1:nxr+1),        &
    279               ippe(nzb+1:nzt,nxl-1:nxr+1), m1(nzb+1:nzt,nxl-2:nxr+2),          &
    280               sw(nzb+1:nzt,nxl-1:nxr+1)                                        &
    281             )
    282     imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp
    283 
    284 !
    285 !-- Initialise point of time measuring of the exponential portion (this would
    286 !-- not work if done locally within the loop)
    287     CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'start' )
    288     CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
    289 
    290 !
    291 !-- Outer loop of all j
    292     DO  j = nys, nyn
    293 
    294 !
    295 !--    Compute polynomial coefficients
    296        DO  i = nxl-1, nxr+1
    297           DO  k = nzb+1, nzt
    298              a12(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
    299              a22(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i)        &
    300                                                  + sk_p(k,j,i-1) )
    301              a0(k,i) = ( 9.0_wp * sk_p(k,j,i+2)    - 116.0_wp * sk_p(k,j,i+1)  &
    302                          + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j,i-1)  &
    303                          + 9.0_wp * sk_p(k,j,i-2) ) * f1920
    304              a1(k,i) = ( -5.0_wp * sk_p(k,j,i+2)   + 34.0_wp * sk_p(k,j,i+1)   &
    305                          - 34.0_wp * sk_p(k,j,i-1) + 5.0_wp * sk_p(k,j,i-2)    &
    306                        ) * f48
    307              a2(k,i) = ( -3.0_wp * sk_p(k,j,i+2) + 36.0_wp * sk_p(k,j,i+1)     &
    308                          - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j,i-1)     &
    309                          - 3.0_wp * sk_p(k,j,i-2) ) * f48
    310           ENDDO
    311        ENDDO
    312 
    313 !
    314 !--    Fluxes using the Bott scheme
    315 !--    *VOCL LOOP,UNROLL(2)
    316550       DO  i = nxl, nxr
    317551          DO  k = nzb+1, nzt
    318              cip  =  MAX( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
    319              cim  = -MIN( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
    320              cipf = 1.0_wp - 2.0_wp * cip
    321              cimf = 1.0_wp - 2.0_wp * cim
    322              ip   =   a0(k,i)   * f2  * ( 1.0_wp - cipf )                      &
    323                     + a1(k,i)   * f8  * ( 1.0_wp - cipf*cipf )                 &
    324                     + a2(k,i)   * f24 * ( 1.0_wp - cipf*cipf*cipf )
    325              im   =   a0(k,i+1) * f2  * ( 1.0_wp - cimf )                      &
    326                     - a1(k,i+1) * f8  * ( 1.0_wp - cimf*cimf )                 &
    327                     + a2(k,i+1) * f24 * ( 1.0_wp - cimf*cimf*cimf )
    328              ip   = MAX( ip, 0.0_wp )
    329              im   = MAX( im, 0.0_wp )
    330              ippb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
    331              impb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i+1) / (ip+im+1E-15_wp) )
    332 
    333              cip  =  MAX( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
    334              cim  = -MIN( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
    335              cipf = 1.0_wp - 2.0_wp * cip
    336              cimf = 1.0_wp - 2.0_wp * cim
    337              ip   =   a0(k,i-1) * f2  * ( 1.0_wp - cipf )                      &
    338                     + a1(k,i-1) * f8  * ( 1.0_wp - cipf*cipf )                 &
    339                     + a2(k,i-1) * f24 * ( 1.0_wp - cipf*cipf*cipf )
    340              im   =   a0(k,i)   * f2  * ( 1.0_wp - cimf )                      &
    341                     - a1(k,i)   * f8  * ( 1.0_wp - cimf*cimf )                 &
    342                     + a2(k,i)   * f24 * ( 1.0_wp - cimf*cimf*cimf )
    343              ip   = MAX( ip, 0.0_wp )
    344              im   = MAX( im, 0.0_wp )
    345              ipmb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i-1) / (ip+im+1E-15_wp) )
    346              immb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
     552             sk_p(k,nys-1,i) = sk_p(k,nyn,i)
     553             sk_p(k,nys-2,i) = sk_p(k,nyn-1,i)
     554             sk_p(k,nys-3,i) = sk_p(k,nyn-2,i)
     555             sk_p(k,nyn+1,i) = sk_p(k,nys,i)
     556             sk_p(k,nyn+2,i) = sk_p(k,nys+1,i)
     557             sk_p(k,nyn+3,i) = sk_p(k,nys+2,i)
    347558          ENDDO
    348559       ENDDO
    349 
    350 !
    351 !--    Compute monitor function m1
    352        DO  i = nxl-2, nxr+2
    353           DO  k = nzb+1, nzt
    354              m1z = ABS( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) )
    355              m1n = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
    356              IF ( m1n /= 0.0_wp  .AND.  m1n >= m1z )  THEN
    357                 m1(k,i) = m1z / m1n
    358                 IF ( m1(k,i) /= 2.0_wp  .AND.  m1n < fmax(2) )  m1(k,i) = 0.0_wp
    359              ELSEIF ( m1n < m1z )  THEN
    360                 m1(k,i) = -1.0_wp
    361              ELSE
    362                 m1(k,i) = 0.0_wp
    363              ENDIF
     560#endif
     561
     562!
     563!--    Determine the maxima of the first and second derivative in y-direction
     564       fmax_l = 0.0_wp
     565       DO  i = nxl, nxr
     566          DO  j = nys, nyn
     567             DO  k = nzb+1, nzt
     568                numera = ABS( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j-1,i) )
     569                denomi  = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
     570                fmax_l(1) = MAX( fmax_l(1) , numera )
     571                fmax_l(2) = MAX( fmax_l(2) , denomi  )
     572             ENDDO
    364573          ENDDO
    365574       ENDDO
    366 
    367 !
    368 !--    Compute switch sw
    369        sw = 0.0_wp
    370        DO  i = nxl-1, nxr+1
    371           DO  k = nzb+1, nzt
    372              m2 = 2.0_wp * ABS( a1(k,i) - a12(k,i) ) /                         &
    373                   MAX( ABS( a1(k,i) + a12(k,i) ), 1E-35_wp )
    374              IF ( ABS( a1(k,i) + a12(k,i) ) < fmax(2) )  m2 = 0.0_wp
    375 
    376              m3 = 2.0_wp * ABS( a2(k,i) - a22(k,i) ) /                         &
    377                   MAX( ABS( a2(k,i) + a22(k,i) ), 1E-35_wp )
    378              IF ( ABS( a2(k,i) + a22(k,i) ) < fmax(1) )  m3 = 0.0_wp
    379 
    380              t1 = 0.35_wp
    381              t2 = 0.35_wp
    382              IF ( m1(k,i) == -1.0_wp )  t2 = 0.12_wp
    383 
    384 !--          *VOCL STMT,IF(10)
    385              IF ( m1(k,i-1) == 1.0_wp .OR. m1(k,i) == 1.0_wp                   &
    386                   .OR. m1(k,i+1) == 1.0_wp .OR.  m2 > t2  .OR.  m3 > t2  .OR.  &
    387                   ( m1(k,i) > t1  .AND.  m1(k,i-1) /= -1.0_wp  .AND.           &
    388                     m1(k,i) /= -1.0_wp  .AND.  m1(k,i+1) /= -1.0_wp )          &
    389                 )  sw(k,i) = 1.0_wp
    390           ENDDO
    391        ENDDO
    392 
    393 !
    394 !--    Fluxes using the exponential scheme
    395        CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
     575#if defined( __parallel )
     576       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     577       CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr )
     578#else
     579       fmax = fmax_l
     580#endif
     581
     582       fmax = 0.04_wp * fmax
     583
     584!
     585!--    Allocate temporary arrays
     586       ALLOCATE( a0(nzb+1:nzt,nys-1:nyn+1),   a1(nzb+1:nzt,nys-1:nyn+1),          &
     587                 a2(nzb+1:nzt,nys-1:nyn+1),   a12(nzb+1:nzt,nys-1:nyn+1),         &
     588                 a22(nzb+1:nzt,nys-1:nyn+1),  immb(nzb+1:nzt,nys-1:nyn+1),        &
     589                 imme(nzb+1:nzt,nys-1:nyn+1), impb(nzb+1:nzt,nys-1:nyn+1),        &
     590                 impe(nzb+1:nzt,nys-1:nyn+1), ipmb(nzb+1:nzt,nys-1:nyn+1),        &
     591                 ipme(nzb+1:nzt,nys-1:nyn+1), ippb(nzb+1:nzt,nys-1:nyn+1),        &
     592                 ippe(nzb+1:nzt,nys-1:nyn+1), m1(nzb+1:nzt,nys-2:nyn+2),          &
     593                 sw(nzb+1:nzt,nys-1:nyn+1)                                        &
     594               )
     595       imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp
     596
     597!
     598!--    Outer loop of all i
    396599       DO  i = nxl, nxr
    397           DO  k = nzb+1, nzt
    398 
    399 !--          *VOCL STMT,IF(10)
    400              IF ( sw(k,i) == 1.0_wp )  THEN
    401                 snenn = sk_p(k,j,i+1) - sk_p(k,j,i-1)
    402                 IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
    403                 sterm = ( sk_p(k,j,i) - sk_p(k,j,i-1) ) / snenn
    404                 sterm = MIN( sterm, 0.9999_wp )
    405                 sterm = MAX( sterm, 0.0001_wp )
    406 
    407                 ix = INT( sterm * 1000 ) + 1
    408 
    409                 cip =  MAX( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
    410 
    411                 ippe(k,i) = sk_p(k,j,i-1) * cip + snenn * (                    &
    412                             aex(ix) * cip + bex(ix) / dex(ix) * (              &
    413                             eex(ix) -                                          &
    414                             EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
    415                                                                 )              &
    416                                                           )
    417                 IF ( sterm == 0.0001_wp )  ippe(k,i) = sk_p(k,j,i) * cip
    418                 IF ( sterm == 0.9999_wp )  ippe(k,i) = sk_p(k,j,i) * cip
    419 
    420                 snenn = sk_p(k,j,i-1) - sk_p(k,j,i+1)
    421                 IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
    422                 sterm = ( sk_p(k,j,i) - sk_p(k,j,i+1) ) / snenn
    423                 sterm = MIN( sterm, 0.9999_wp )
    424                 sterm = MAX( sterm, 0.0001_wp )
    425 
    426                 ix = INT( sterm * 1000 ) + 1
    427 
    428                 cim = -MIN( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
    429 
    430                 imme(k,i) = sk_p(k,j,i+1) * cim + snenn * (                    &
    431                             aex(ix) * cim + bex(ix) / dex(ix) * (              &
    432                             eex(ix) -                                          &
    433                             EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
    434                                                                 )              &
    435                                                           )
    436                 IF ( sterm == 0.0001_wp )  imme(k,i) = sk_p(k,j,i) * cim
    437                 IF ( sterm == 0.9999_wp )  imme(k,i) = sk_p(k,j,i) * cim
    438              ENDIF
    439 
    440 !--          *VOCL STMT,IF(10)
    441              IF ( sw(k,i+1) == 1.0_wp )  THEN
    442                 snenn = sk_p(k,j,i) - sk_p(k,j,i+2)
    443                 IF ( ABS( snenn ) .LT. 1E-9_wp )  snenn = 1E-9_wp
    444                 sterm = ( sk_p(k,j,i+1) - sk_p(k,j,i+2) ) / snenn
    445                 sterm = MIN( sterm, 0.9999_wp )
    446                 sterm = MAX( sterm, 0.0001_wp )
    447 
    448                 ix = INT( sterm * 1000 ) + 1
    449 
    450                 cim = -MIN( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
    451 
    452                 impe(k,i) = sk_p(k,j,i+2) * cim + snenn * (                    &
    453                             aex(ix) * cim + bex(ix) / dex(ix) * (              &
    454                             eex(ix) -                                          &
    455                             EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
    456                                                                 )              &
    457                                                           )
    458                 IF ( sterm == 0.0001_wp )  impe(k,i) = sk_p(k,j,i+1) * cim
    459                 IF ( sterm == 0.9999_wp )  impe(k,i) = sk_p(k,j,i+1) * cim
    460              ENDIF
    461 
    462 !--          *VOCL STMT,IF(10)
    463              IF ( sw(k,i-1) == 1.0_wp )  THEN
    464                 snenn = sk_p(k,j,i) - sk_p(k,j,i-2)
    465                 IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
    466                 sterm = ( sk_p(k,j,i-1) - sk_p(k,j,i-2) ) / snenn
    467                 sterm = MIN( sterm, 0.9999_wp )
    468                 sterm = MAX( sterm, 0.0001_wp )
    469 
    470                 ix = INT( sterm * 1000 ) + 1
    471 
    472                 cip = MAX( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
    473 
    474                 ipme(k,i) = sk_p(k,j,i-2) * cip + snenn * (                    &
    475                             aex(ix) * cip + bex(ix) / dex(ix) * (              &
    476                             eex(ix) -                                          &
    477                             EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
    478                                                                 )              &
    479                                                           )
    480                 IF ( sterm == 0.0001_wp )  ipme(k,i) = sk_p(k,j,i-1) * cip
    481                 IF ( sterm == 0.9999_wp )  ipme(k,i) = sk_p(k,j,i-1) * cip
    482              ENDIF
    483 
    484           ENDDO
    485        ENDDO
    486        CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
    487 
    488 !
    489 !--    Prognostic equation
    490        DO  i = nxl, nxr
    491           DO  k = nzb+1, nzt
    492              fplus  = ( 1.0_wp - sw(k,i)   ) * ippb(k,i) + sw(k,i)   * ippe(k,i)  &
    493                     - ( 1.0_wp - sw(k,i+1) ) * impb(k,i) - sw(k,i+1) * impe(k,i)
    494              fminus = ( 1.0_wp - sw(k,i-1) ) * ipmb(k,i) + sw(k,i-1) * ipme(k,i)  &
    495                     - ( 1.0_wp - sw(k,i)   ) * immb(k,i) - sw(k,i)   * imme(k,i)
    496              tendcy = fplus - fminus
    497 !
    498 !--           Removed in order to optimize speed
    499 !             ffmax   = MAX( ABS( fplus ), ABS( fminus ), 1E-35_wp )
    500 !             IF ( ( ABS( tendcy ) / ffmax ) < 1E-7_wp )  tendcy = 0.0
    501 !
    502 !--          Density correction because of possible remaining divergences
    503              d_new = d(k,j,i) - ( u(k,j,i+1) - u(k,j,i) ) * dt_3d * ddx
    504              sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /    &
    505                            ( 1.0_wp + d_new )
    506              d(k,j,i)  = d_new
    507           ENDDO
    508        ENDDO
    509 
    510     ENDDO   ! End of the advection in x-direction
    511 
    512 !
    513 !-- Deallocate temporary arrays
    514     DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme,      &
    515                 ippb, ippe, m1, sw )
    516 
    517 
    518 !
    519 !-- Enlarge boundary of local array cyclically in y-direction
    520 #if defined( __parallel )
    521     ngp = ( nzt - nzb + 6 ) * ( nyn - nys + 7 )
    522     CALL MPI_TYPE_VECTOR( nxr-nxl+7, 3*(nzt-nzb+6), ngp, MPI_REAL,             &
    523                           type_xz_2, ierr )
    524     CALL MPI_TYPE_COMMIT( type_xz_2, ierr )
    525 !
    526 !-- Send front boundary, receive rear boundary
    527     CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' )
    528     CALL MPI_SENDRECV( sk_p(nzb-2,nys,nxl-3),   1, type_xz_2, psouth, 0,       &
    529                        sk_p(nzb-2,nyn+1,nxl-3), 1, type_xz_2, pnorth, 0,       &
    530                        comm2d, status, ierr )
    531 !
    532 !-- Send rear boundary, receive front boundary
    533     CALL MPI_SENDRECV( sk_p(nzb-2,nyn-2,nxl-3), 1, type_xz_2, pnorth, 1,       &
    534                        sk_p(nzb-2,nys-3,nxl-3), 1, type_xz_2, psouth, 1,       &
    535                        comm2d, status, ierr )
    536     CALL MPI_TYPE_FREE( type_xz_2, ierr )
    537     CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' )
    538 #else
    539     DO  i = nxl, nxr
    540        DO  k = nzb+1, nzt
    541           sk_p(k,nys-1,i) = sk_p(k,nyn,i)
    542           sk_p(k,nys-2,i) = sk_p(k,nyn-1,i)
    543           sk_p(k,nys-3,i) = sk_p(k,nyn-2,i)
    544           sk_p(k,nyn+1,i) = sk_p(k,nys,i)
    545           sk_p(k,nyn+2,i) = sk_p(k,nys+1,i)
    546           sk_p(k,nyn+3,i) = sk_p(k,nys+2,i)
    547        ENDDO
    548     ENDDO
    549 #endif
    550 
    551 !
    552 !-- Determine the maxima of the first and second derivative in y-direction
    553     fmax_l = 0.0_wp
    554     DO  i = nxl, nxr
    555        DO  j = nys, nyn
    556           DO  k = nzb+1, nzt
    557              numera = ABS( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j-1,i) )
    558              denomi  = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
    559              fmax_l(1) = MAX( fmax_l(1) , numera )
    560              fmax_l(2) = MAX( fmax_l(2) , denomi  )
    561           ENDDO
    562        ENDDO
    563     ENDDO
    564 #if defined( __parallel )
    565     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    566     CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr )
    567 #else
    568     fmax = fmax_l
    569 #endif
    570 
    571     fmax = 0.04_wp * fmax
    572 
    573 !
    574 !-- Allocate temporary arrays
    575     ALLOCATE( a0(nzb+1:nzt,nys-1:nyn+1),   a1(nzb+1:nzt,nys-1:nyn+1),          &
    576               a2(nzb+1:nzt,nys-1:nyn+1),   a12(nzb+1:nzt,nys-1:nyn+1),         &
    577               a22(nzb+1:nzt,nys-1:nyn+1),  immb(nzb+1:nzt,nys-1:nyn+1),        &
    578               imme(nzb+1:nzt,nys-1:nyn+1), impb(nzb+1:nzt,nys-1:nyn+1),        &
    579               impe(nzb+1:nzt,nys-1:nyn+1), ipmb(nzb+1:nzt,nys-1:nyn+1),        &
    580               ipme(nzb+1:nzt,nys-1:nyn+1), ippb(nzb+1:nzt,nys-1:nyn+1),        &
    581               ippe(nzb+1:nzt,nys-1:nyn+1), m1(nzb+1:nzt,nys-2:nyn+2),          &
    582               sw(nzb+1:nzt,nys-1:nyn+1)                                        &
    583             )
    584     imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp
    585 
    586 !
    587 !-- Outer loop of all i
    588     DO  i = nxl, nxr
    589 
    590 !
    591 !--    Compute polynomial coefficients
    592        DO  j = nys-1, nyn+1
    593           DO  k = nzb+1, nzt
    594              a12(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
    595              a22(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i)        &
    596                                                  + sk_p(k,j-1,i) )
    597              a0(k,j) = ( 9.0_wp * sk_p(k,j+2,i)    - 116.0_wp * sk_p(k,j+1,i)  &
    598                          + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j-1,i)  &
    599                          + 9.0_wp * sk_p(k,j-2,i) ) * f1920
    600              a1(k,j) = ( -5.0_wp   * sk_p(k,j+2,i) + 34.0_wp * sk_p(k,j+1,i)   &
    601                          - 34.0_wp * sk_p(k,j-1,i) + 5.0_wp  * sk_p(k,j-2,i)   &
    602                        ) * f48
    603              a2(k,j) = ( -3.0_wp * sk_p(k,j+2,i) + 36.0_wp * sk_p(k,j+1,i)     &
    604                          - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j-1,i)     &
    605                          - 3.0_wp * sk_p(k,j-2,i) ) * f48
    606           ENDDO
    607        ENDDO
    608 
    609 !
    610 !--    Fluxes using the Bott scheme
    611 !--    *VOCL LOOP,UNROLL(2)
    612        DO  j = nys, nyn
    613           DO  k = nzb+1, nzt
    614              cip  =  MAX( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
    615              cim  = -MIN( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
    616              cipf = 1.0_wp - 2.0_wp * cip
    617              cimf = 1.0_wp - 2.0_wp * cim
    618              ip   =   a0(k,j)   * f2  * ( 1.0_wp - cipf )                      &
    619                     + a1(k,j)   * f8  * ( 1.0_wp - cipf*cipf )                 &
    620                     + a2(k,j)   * f24 * ( 1.0_wp - cipf*cipf*cipf )
    621              im   =   a0(k,j+1) * f2  * ( 1.0_wp - cimf )                      &
    622                     - a1(k,j+1) * f8  * ( 1.0_wp - cimf*cimf )                 &
    623                     + a2(k,j+1) * f24 * ( 1.0_wp - cimf*cimf*cimf )
    624              ip   = MAX( ip, 0.0_wp )
    625              im   = MAX( im, 0.0_wp )
    626              ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
    627              impb(k,j) = im * MIN( 1.0_wp, sk_p(k,j+1,i) / (ip+im+1E-15_wp) )
    628 
    629              cip  =  MAX( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
    630              cim  = -MIN( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
    631              cipf = 1.0_wp - 2.0_wp * cip
    632              cimf = 1.0_wp - 2.0_wp * cim
    633              ip   =   a0(k,j-1) * f2  * ( 1.0_wp - cipf )                      &
    634                     + a1(k,j-1) * f8  * ( 1.0_wp - cipf*cipf )                 &
    635                     + a2(k,j-1) * f24 * ( 1.0_wp - cipf*cipf*cipf )
    636              im   =   a0(k,j)   * f2  * ( 1.0_wp - cimf )                      &
    637                     - a1(k,j)   * f8  * ( 1.0_wp - cimf*cimf )                 &
    638                     + a2(k,j)   * f24 * ( 1.0_wp - cimf*cimf*cimf )
    639              ip   = MAX( ip, 0.0_wp )
    640              im   = MAX( im, 0.0_wp )
    641              ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j-1,i) / (ip+im+1E-15_wp) )
    642              immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
    643           ENDDO
    644        ENDDO
    645 
    646 !
    647 !--    Compute monitor function m1
    648        DO  j = nys-2, nyn+2
    649           DO  k = nzb+1, nzt
    650              m1z = ABS( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j-1,i) )
    651              m1n = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
    652              IF ( m1n /= 0.0_wp  .AND.  m1n >= m1z )  THEN
    653                 m1(k,j) = m1z / m1n
    654                 IF ( m1(k,j) /= 2.0_wp  .AND.  m1n < fmax(2) )  m1(k,j) = 0.0_wp
    655              ELSEIF ( m1n < m1z )  THEN
    656                 m1(k,j) = -1.0_wp
    657              ELSE
    658                 m1(k,j) = 0.0_wp
    659              ENDIF
    660           ENDDO
    661        ENDDO
    662 
    663 !
    664 !--    Compute switch sw
    665        sw = 0.0_wp
    666        DO  j = nys-1, nyn+1
    667           DO  k = nzb+1, nzt
    668              m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) /                            &
    669                   MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp )
    670              IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) )  m2 = 0.0_wp
    671 
    672              m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) /                            &
    673                   MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp )
    674              IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) )  m3 = 0.0_wp
    675 
    676              t1 = 0.35_wp
    677              t2 = 0.35_wp
    678              IF ( m1(k,j) == -1.0_wp )  t2 = 0.12_wp
    679 
    680 !--          *VOCL STMT,IF(10)
    681              IF ( m1(k,j-1) == 1.0_wp .OR. m1(k,j) == 1.0_wp                   &
    682                   .OR. m1(k,j+1) == 1.0_wp .OR.  m2 > t2  .OR.  m3 > t2  .OR.  &
    683                   ( m1(k,j) > t1  .AND.  m1(k,j-1) /= -1.0_wp  .AND.           &
    684                     m1(k,j) /= -1.0_wp  .AND.  m1(k,j+1) /= -1.0_wp )          &
    685                 )  sw(k,j) = 1.0_wp
    686           ENDDO
    687        ENDDO
    688 
    689 !
    690 !--    Fluxes using exponential scheme
    691        CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
    692        DO  j = nys, nyn
    693           DO  k = nzb+1, nzt
    694 
    695 !--          *VOCL STMT,IF(10)
    696              IF ( sw(k,j) == 1.0_wp )  THEN
    697                 snenn = sk_p(k,j+1,i) - sk_p(k,j-1,i)
    698                 IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
    699                 sterm = ( sk_p(k,j,i) - sk_p(k,j-1,i) ) / snenn
    700                 sterm = MIN( sterm, 0.9999_wp )
    701                 sterm = MAX( sterm, 0.0001_wp )
    702 
    703                 ix = INT( sterm * 1000 ) + 1
    704 
    705                 cip =  MAX( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
    706 
    707                 ippe(k,j) = sk_p(k,j-1,i) * cip + snenn * (                    &
    708                             aex(ix) * cip + bex(ix) / dex(ix) * (              &
    709                             eex(ix) -                                          &
    710                             EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
    711                                                                 )              &
    712                                                           )
    713                 IF ( sterm == 0.0001_wp )  ippe(k,j) = sk_p(k,j,i) * cip
    714                 IF ( sterm == 0.9999_wp )  ippe(k,j) = sk_p(k,j,i) * cip
    715 
    716                 snenn = sk_p(k,j-1,i) - sk_p(k,j+1,i)
    717                 IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
    718                 sterm = ( sk_p(k,j,i) - sk_p(k,j+1,i) ) / snenn
    719                 sterm = MIN( sterm, 0.9999_wp )
    720                 sterm = MAX( sterm, 0.0001_wp )
    721 
    722                 ix = INT( sterm * 1000 ) + 1
    723 
    724                 cim = -MIN( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
    725 
    726                 imme(k,j) = sk_p(k,j+1,i) * cim + snenn * (                    &
    727                             aex(ix) * cim + bex(ix) / dex(ix) * (              &
    728                             eex(ix) -                                          &
    729                             EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
    730                                                                 )              &
    731                                                           )
    732                 IF ( sterm == 0.0001_wp )  imme(k,j) = sk_p(k,j,i) * cim
    733                 IF ( sterm == 0.9999_wp )  imme(k,j) = sk_p(k,j,i) * cim
    734              ENDIF
    735 
    736 !--          *VOCL STMT,IF(10)
    737              IF ( sw(k,j+1) == 1.0_wp )  THEN
    738                 snenn = sk_p(k,j,i) - sk_p(k,j+2,i)
    739                 IF ( ABS( snenn ) .LT. 1E-9_wp )  snenn = 1E-9_wp
    740                 sterm = ( sk_p(k,j+1,i) - sk_p(k,j+2,i) ) / snenn
    741                 sterm = MIN( sterm, 0.9999_wp )
    742                 sterm = MAX( sterm, 0.0001_wp )
    743 
    744                 ix = INT( sterm * 1000 ) + 1
    745 
    746                 cim = -MIN( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
    747 
    748                 impe(k,j) = sk_p(k,j+2,i) * cim + snenn * (                    &
    749                             aex(ix) * cim + bex(ix) / dex(ix) * (              &
    750                             eex(ix) -                                          &
    751                             EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
    752                                                                 )              &
    753                                                           )
    754                 IF ( sterm == 0.0001_wp )  impe(k,j) = sk_p(k,j+1,i) * cim
    755                 IF ( sterm == 0.9999_wp )  impe(k,j) = sk_p(k,j+1,i) * cim
    756              ENDIF
    757 
    758 !--          *VOCL STMT,IF(10)
    759              IF ( sw(k,j-1) == 1.0_wp )  THEN
    760                 snenn = sk_p(k,j,i) - sk_p(k,j-2,i)
    761                 IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
    762                 sterm = ( sk_p(k,j-1,i) - sk_p(k,j-2,i) ) / snenn
    763                 sterm = MIN( sterm, 0.9999_wp )
    764                 sterm = MAX( sterm, 0.0001_wp )
    765 
    766                 ix = INT( sterm * 1000 ) + 1
    767 
    768                 cip = MAX( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
    769 
    770                 ipme(k,j) = sk_p(k,j-2,i) * cip + snenn * (                    &
    771                             aex(ix) * cip + bex(ix) / dex(ix) * (              &
    772                             eex(ix) -                                          &
    773                             EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
    774                                                                 )              &
    775                                                           )
    776                 IF ( sterm == 0.0001_wp )  ipme(k,j) = sk_p(k,j-1,i) * cip
    777                 IF ( sterm == 0.9999_wp )  ipme(k,j) = sk_p(k,j-1,i) * cip
    778              ENDIF
    779 
    780           ENDDO
    781        ENDDO
    782        CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
    783 
    784 !
    785 !--    Prognostic equation
    786        DO  j = nys, nyn
    787           DO  k = nzb+1, nzt
    788              fplus  = ( 1.0_wp - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
    789                     - ( 1.0_wp - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j)
    790              fminus = ( 1.0_wp - sw(k,j-1) ) * ipmb(k,j) + sw(k,j-1) * ipme(k,j) &
    791                     - ( 1.0_wp - sw(k,j)   ) * immb(k,j) - sw(k,j)   * imme(k,j)
    792              tendcy = fplus - fminus
    793 !
    794 !--           Removed in order to optimise speed
    795 !             ffmax   = MAX( ABS( fplus ), ABS( fminus ), 1E-35_wp )
    796 !             IF ( ( ABS( tendcy ) / ffmax ) < 1E-7_wp )  tendcy = 0.0
    797 !
    798 !--          Density correction because of possible remaining divergences
    799              d_new = d(k,j,i) - ( v(k,j+1,i) - v(k,j,i) ) * dt_3d * ddy
    800              sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /  &
    801                            ( 1.0_wp + d_new )
    802              d(k,j,i)  = d_new
    803           ENDDO
    804        ENDDO
    805 
    806     ENDDO   ! End of the advection in y-direction
    807     CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' )
    808     CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'stop' )
    809 
    810 !
    811 !-- Deallocate temporary arrays
    812     DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme,      &
    813                 ippb, ippe, m1, sw )
    814 
    815 
    816 !
    817 !-- Initialise for the computation of heat fluxes (see below; required in
    818 !-- UP flow_statistics)
    819     IF ( sk_char == 'pt' )  sums_wsts_bc_l = 0.0_wp
    820 
    821 !
    822 !-- Add top and bottom boundaries according to the relevant boundary conditions
    823     IF ( sk_char == 'pt' )  THEN
    824 
    825 !
    826 !--    Temperature boundary condition at the bottom boundary
    827        IF ( ibc_pt_b == 0 )  THEN
     600
     601!
     602!--       Compute polynomial coefficients
     603          DO  j = nys-1, nyn+1
     604             DO  k = nzb+1, nzt
     605                a12(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
     606                a22(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i)        &
     607                                                    + sk_p(k,j-1,i) )
     608                a0(k,j) = ( 9.0_wp * sk_p(k,j+2,i)    - 116.0_wp * sk_p(k,j+1,i)  &
     609                            + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j-1,i)  &
     610                            + 9.0_wp * sk_p(k,j-2,i) ) * f1920
     611                a1(k,j) = ( -5.0_wp   * sk_p(k,j+2,i) + 34.0_wp * sk_p(k,j+1,i)   &
     612                            - 34.0_wp * sk_p(k,j-1,i) + 5.0_wp  * sk_p(k,j-2,i)   &
     613                          ) * f48
     614                a2(k,j) = ( -3.0_wp * sk_p(k,j+2,i) + 36.0_wp * sk_p(k,j+1,i)     &
     615                            - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j-1,i)     &
     616                            - 3.0_wp * sk_p(k,j-2,i) ) * f48
     617             ENDDO
     618          ENDDO
     619
     620!
     621!--       Fluxes using the Bott scheme
     622!--       *VOCL LOOP,UNROLL(2)
     623          DO  j = nys, nyn
     624             DO  k = nzb+1, nzt
     625                cip  =  MAX( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
     626                cim  = -MIN( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
     627                cipf = 1.0_wp - 2.0_wp * cip
     628                cimf = 1.0_wp - 2.0_wp * cim
     629                ip   =   a0(k,j)   * f2  * ( 1.0_wp - cipf )                      &
     630                       + a1(k,j)   * f8  * ( 1.0_wp - cipf*cipf )                 &
     631                       + a2(k,j)   * f24 * ( 1.0_wp - cipf*cipf*cipf )
     632                im   =   a0(k,j+1) * f2  * ( 1.0_wp - cimf )                      &
     633                       - a1(k,j+1) * f8  * ( 1.0_wp - cimf*cimf )                 &
     634                       + a2(k,j+1) * f24 * ( 1.0_wp - cimf*cimf*cimf )
     635                ip   = MAX( ip, 0.0_wp )
     636                im   = MAX( im, 0.0_wp )
     637                ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
     638                impb(k,j) = im * MIN( 1.0_wp, sk_p(k,j+1,i) / (ip+im+1E-15_wp) )
     639
     640                cip  =  MAX( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
     641                cim  = -MIN( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
     642                cipf = 1.0_wp - 2.0_wp * cip
     643                cimf = 1.0_wp - 2.0_wp * cim
     644                ip   =   a0(k,j-1) * f2  * ( 1.0_wp - cipf )                      &
     645                       + a1(k,j-1) * f8  * ( 1.0_wp - cipf*cipf )                 &
     646                       + a2(k,j-1) * f24 * ( 1.0_wp - cipf*cipf*cipf )
     647                im   =   a0(k,j)   * f2  * ( 1.0_wp - cimf )                      &
     648                       - a1(k,j)   * f8  * ( 1.0_wp - cimf*cimf )                 &
     649                       + a2(k,j)   * f24 * ( 1.0_wp - cimf*cimf*cimf )
     650                ip   = MAX( ip, 0.0_wp )
     651                im   = MAX( im, 0.0_wp )
     652                ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j-1,i) / (ip+im+1E-15_wp) )
     653                immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
     654             ENDDO
     655          ENDDO
     656
     657!
     658!--       Compute monitor function m1
     659          DO  j = nys-2, nyn+2
     660             DO  k = nzb+1, nzt
     661                m1z = ABS( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j-1,i) )
     662                m1n = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
     663                IF ( m1n /= 0.0_wp  .AND.  m1n >= m1z )  THEN
     664                   m1(k,j) = m1z / m1n
     665                   IF ( m1(k,j) /= 2.0_wp  .AND.  m1n < fmax(2) )  m1(k,j) = 0.0_wp
     666                ELSEIF ( m1n < m1z )  THEN
     667                   m1(k,j) = -1.0_wp
     668                ELSE
     669                   m1(k,j) = 0.0_wp
     670                ENDIF
     671             ENDDO
     672          ENDDO
     673
     674!
     675!--       Compute switch sw
     676          sw = 0.0_wp
     677          DO  j = nys-1, nyn+1
     678             DO  k = nzb+1, nzt
     679                m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) /                            &
     680                     MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp )
     681                IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) )  m2 = 0.0_wp
     682
     683                m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) /                            &
     684                     MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp )
     685                IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) )  m3 = 0.0_wp
     686
     687                t1 = 0.35_wp
     688                t2 = 0.35_wp
     689                IF ( m1(k,j) == -1.0_wp )  t2 = 0.12_wp
     690
     691!--             *VOCL STMT,IF(10)
     692                IF ( m1(k,j-1) == 1.0_wp .OR. m1(k,j) == 1.0_wp                   &
     693                     .OR. m1(k,j+1) == 1.0_wp .OR.  m2 > t2  .OR.  m3 > t2  .OR.  &
     694                     ( m1(k,j) > t1  .AND.  m1(k,j-1) /= -1.0_wp  .AND.           &
     695                       m1(k,j) /= -1.0_wp  .AND.  m1(k,j+1) /= -1.0_wp )          &
     696                   )  sw(k,j) = 1.0_wp
     697             ENDDO
     698          ENDDO
     699
     700!
     701!--       Fluxes using exponential scheme
     702          CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
     703          DO  j = nys, nyn
     704             DO  k = nzb+1, nzt
     705
     706!--             *VOCL STMT,IF(10)
     707                IF ( sw(k,j) == 1.0_wp )  THEN
     708                   snenn = sk_p(k,j+1,i) - sk_p(k,j-1,i)
     709                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
     710                   sterm = ( sk_p(k,j,i) - sk_p(k,j-1,i) ) / snenn
     711                   sterm = MIN( sterm, 0.9999_wp )
     712                   sterm = MAX( sterm, 0.0001_wp )
     713
     714                   ix = INT( sterm * 1000 ) + 1
     715
     716                   cip =  MAX( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
     717
     718                   ippe(k,j) = sk_p(k,j-1,i) * cip + snenn * (                    &
     719                               aex(ix) * cip + bex(ix) / dex(ix) * (              &
     720                               eex(ix) -                                          &
     721                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
     722                                                                   )              &
     723                                                             )
     724                   IF ( sterm == 0.0001_wp )  ippe(k,j) = sk_p(k,j,i) * cip
     725                   IF ( sterm == 0.9999_wp )  ippe(k,j) = sk_p(k,j,i) * cip
     726
     727                   snenn = sk_p(k,j-1,i) - sk_p(k,j+1,i)
     728                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
     729                   sterm = ( sk_p(k,j,i) - sk_p(k,j+1,i) ) / snenn
     730                   sterm = MIN( sterm, 0.9999_wp )
     731                   sterm = MAX( sterm, 0.0001_wp )
     732
     733                   ix = INT( sterm * 1000 ) + 1
     734
     735                   cim = -MIN( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
     736
     737                   imme(k,j) = sk_p(k,j+1,i) * cim + snenn * (                    &
     738                               aex(ix) * cim + bex(ix) / dex(ix) * (              &
     739                               eex(ix) -                                          &
     740                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
     741                                                                   )              &
     742                                                             )
     743                   IF ( sterm == 0.0001_wp )  imme(k,j) = sk_p(k,j,i) * cim
     744                   IF ( sterm == 0.9999_wp )  imme(k,j) = sk_p(k,j,i) * cim
     745                ENDIF
     746
     747!--             *VOCL STMT,IF(10)
     748                IF ( sw(k,j+1) == 1.0_wp )  THEN
     749                   snenn = sk_p(k,j,i) - sk_p(k,j+2,i)
     750                   IF ( ABS( snenn ) .LT. 1E-9_wp )  snenn = 1E-9_wp
     751                   sterm = ( sk_p(k,j+1,i) - sk_p(k,j+2,i) ) / snenn
     752                   sterm = MIN( sterm, 0.9999_wp )
     753                   sterm = MAX( sterm, 0.0001_wp )
     754
     755                   ix = INT( sterm * 1000 ) + 1
     756
     757                   cim = -MIN( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
     758
     759                   impe(k,j) = sk_p(k,j+2,i) * cim + snenn * (                    &
     760                               aex(ix) * cim + bex(ix) / dex(ix) * (              &
     761                               eex(ix) -                                          &
     762                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
     763                                                                   )              &
     764                                                             )
     765                   IF ( sterm == 0.0001_wp )  impe(k,j) = sk_p(k,j+1,i) * cim
     766                   IF ( sterm == 0.9999_wp )  impe(k,j) = sk_p(k,j+1,i) * cim
     767                ENDIF
     768
     769!--             *VOCL STMT,IF(10)
     770                IF ( sw(k,j-1) == 1.0_wp )  THEN
     771                   snenn = sk_p(k,j,i) - sk_p(k,j-2,i)
     772                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
     773                   sterm = ( sk_p(k,j-1,i) - sk_p(k,j-2,i) ) / snenn
     774                   sterm = MIN( sterm, 0.9999_wp )
     775                   sterm = MAX( sterm, 0.0001_wp )
     776
     777                   ix = INT( sterm * 1000 ) + 1
     778
     779                   cip = MAX( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
     780
     781                   ipme(k,j) = sk_p(k,j-2,i) * cip + snenn * (                    &
     782                               aex(ix) * cip + bex(ix) / dex(ix) * (              &
     783                               eex(ix) -                                          &
     784                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
     785                                                                   )              &
     786                                                             )
     787                   IF ( sterm == 0.0001_wp )  ipme(k,j) = sk_p(k,j-1,i) * cip
     788                   IF ( sterm == 0.9999_wp )  ipme(k,j) = sk_p(k,j-1,i) * cip
     789                ENDIF
     790
     791             ENDDO
     792          ENDDO
     793          CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
     794
     795!
     796!--       Prognostic equation
     797          DO  j = nys, nyn
     798             DO  k = nzb+1, nzt
     799                fplus  = ( 1.0_wp - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
     800                       - ( 1.0_wp - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j)
     801                fminus = ( 1.0_wp - sw(k,j-1) ) * ipmb(k,j) + sw(k,j-1) * ipme(k,j) &
     802                       - ( 1.0_wp - sw(k,j)   ) * immb(k,j) - sw(k,j)   * imme(k,j)
     803                tendcy = fplus - fminus
     804!
     805!--             Removed in order to optimise speed
     806!                ffmax   = MAX( ABS( fplus ), ABS( fminus ), 1E-35_wp )
     807!                 IF ( ( ABS( tendcy ) / ffmax ) < 1E-7_wp )  tendcy = 0.0
     808!
     809!--             Density correction because of possible remaining divergences
     810                d_new = d(k,j,i) - ( v(k,j+1,i) - v(k,j,i) ) * dt_3d * ddy
     811                sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /  &
     812                              ( 1.0_wp + d_new )
     813                d(k,j,i)  = d_new
     814             ENDDO
     815          ENDDO
     816
     817       ENDDO   ! End of the advection in y-direction
     818       CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' )
     819       CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'stop' )
     820
     821!
     822!--    Deallocate temporary arrays
     823       DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme,      &
     824                   ippb, ippe, m1, sw )
     825
     826
     827!
     828!--    Initialise for the computation of heat fluxes (see below; required in
     829!--    UP flow_statistics)
     830       IF ( sk_char == 'pt' )  sums_wsts_bc_l = 0.0_wp
     831
     832!
     833!--    Add top and bottom boundaries according to the relevant boundary conditions
     834       IF ( sk_char == 'pt' )  THEN
     835
     836!
     837!--       Temperature boundary condition at the bottom boundary
     838          IF ( ibc_pt_b == 0 )  THEN
    828839!
    829840!--       Dirichlet (fixed surface temperature)
    830           DO  i = nxl, nxr
    831              DO  j = nys, nyn
    832                 sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
    833                 sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
    834              ENDDO
    835           ENDDO
    836 
    837        ELSE
    838 !
    839 !--       Neumann (i.e. here zero gradient)
     841             DO  i = nxl, nxr
     842                DO  j = nys, nyn
     843                   sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     844                   sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
     845                ENDDO
     846             ENDDO
     847
     848          ELSE
     849!
     850!--          Neumann (i.e. here zero gradient)
     851             DO  i = nxl, nxr
     852                DO  j = nys, nyn
     853                   sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
     854                   sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     855                   sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
     856                ENDDO
     857             ENDDO
     858
     859          ENDIF
     860
     861!
     862!--       Temperature boundary condition at the top boundary
     863          IF ( ibc_pt_t == 0  .OR.  ibc_pt_t == 1 )  THEN
     864!
     865!--          Dirichlet or Neumann (zero gradient)
     866             DO  i = nxl, nxr
     867                DO  j = nys, nyn
     868                   sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
     869                   sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
     870                ENDDO
     871             ENDDO
     872
     873          ELSEIF ( ibc_pt_t == 2 )  THEN
     874!
     875!--          Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead
     876             DO  i = nxl, nxr
     877                DO  j = nys, nyn
     878                   sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i) + bc_pt_t_val * dzu(nzt+1)
     879                   sk_p(nzt+3,j,i)   = sk_p(nzt+2,j,i) + bc_pt_t_val * dzu(nzt+1)
     880                ENDDO
     881             ENDDO
     882
     883          ENDIF
     884
     885       ELSEIF ( sk_char == 'sa' )  THEN
     886
     887!
     888!--       Salinity boundary condition at the bottom boundary.
     889!--       So far, always Neumann (i.e. here zero gradient) is used
    840890          DO  i = nxl, nxr
    841891             DO  j = nys, nyn
     
    846896          ENDDO
    847897
    848        ENDIF
    849 
    850 !
    851 !--    Temperature boundary condition at the top boundary
    852        IF ( ibc_pt_t == 0  .OR.  ibc_pt_t == 1 )  THEN
    853 !
     898!
     899!--       Salinity boundary condition at the top boundary.
    854900!--       Dirichlet or Neumann (zero gradient)
    855901          DO  i = nxl, nxr
     
    860906          ENDDO
    861907
    862        ELSEIF ( ibc_pt_t == 2 )  THEN
    863 !
    864 !--       Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead
     908       ELSEIF ( sk_char == 'q' )  THEN
     909
     910!
     911!--       Specific humidity boundary condition at the bottom boundary.
     912!--       Dirichlet (fixed surface humidity) or Neumann (i.e. zero gradient)
    865913          DO  i = nxl, nxr
    866914             DO  j = nys, nyn
    867                 sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i) + bc_pt_t_val * dzu(nzt+1)
    868                 sk_p(nzt+3,j,i)   = sk_p(nzt+2,j,i) + bc_pt_t_val * dzu(nzt+1)
    869              ENDDO
    870           ENDDO
    871 
    872        ENDIF
    873 
    874     ELSEIF ( sk_char == 'sa' )  THEN
    875 
    876 !
    877 !--    Salinity boundary condition at the bottom boundary.
    878 !--    So far, always Neumann (i.e. here zero gradient) is used
    879        DO  i = nxl, nxr
    880           DO  j = nys, nyn
    881              sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
    882              sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
    883              sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
    884           ENDDO
    885        ENDDO
    886 
    887 !
    888 !--    Salinity boundary condition at the top boundary.
    889 !--    Dirichlet or Neumann (zero gradient)
    890        DO  i = nxl, nxr
    891           DO  j = nys, nyn
    892              sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
    893              sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
    894           ENDDO
    895        ENDDO
    896 
    897     ELSEIF ( sk_char == 'q' )  THEN
    898 
    899 !
    900 !--    Specific humidity boundary condition at the bottom boundary.
    901 !--    Dirichlet (fixed surface humidity) or Neumann (i.e. zero gradient)
    902        DO  i = nxl, nxr
    903           DO  j = nys, nyn
    904              sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
    905              sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
    906              sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
    907           ENDDO
    908        ENDDO
    909 
    910 !
    911 !--    Specific humidity boundary condition at the top boundary
    912        IF ( ibc_q_t == 0 )  THEN
    913 !
    914 !--       Dirichlet
     915                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
     916                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     917                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
     918             ENDDO
     919          ENDDO
     920
     921!
     922!--       Specific humidity boundary condition at the top boundary
     923          IF ( ibc_q_t == 0 )  THEN
     924!
     925!--          Dirichlet
     926             DO  i = nxl, nxr
     927                DO  j = nys, nyn
     928                   sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
     929                   sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
     930                ENDDO
     931             ENDDO
     932
     933          ELSE
     934!
     935!--          Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead
     936             DO  i = nxl, nxr
     937                DO  j = nys, nyn
     938                   sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i) + bc_q_t_val * dzu(nzt+1)
     939                   sk_p(nzt+3,j,i)   = sk_p(nzt+2,j,i) + bc_q_t_val * dzu(nzt+1)
     940                ENDDO
     941             ENDDO
     942
     943          ENDIF
     944
     945       ELSEIF ( sk_char == 'qr' )  THEN
     946
     947!
     948!--       Rain water content boundary condition at the bottom boundary:
     949!--       Dirichlet (fixed surface rain water content).
     950          DO  i = nxl, nxr
     951             DO  j = nys, nyn
     952                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
     953                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     954                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
     955             ENDDO
     956          ENDDO
     957
     958!
     959!--       Rain water content boundary condition at the top boundary: Dirichlet
    915960          DO  i = nxl, nxr
    916961             DO  j = nys, nyn
     
    920965          ENDDO
    921966
    922        ELSE
    923 !
    924 !--       Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead
     967       ELSEIF ( sk_char == 'nr' )  THEN
     968
     969!
     970!--       Rain drop concentration boundary condition at the bottom boundary:
     971!--       Dirichlet (fixed surface rain drop concentration).
    925972          DO  i = nxl, nxr
    926973             DO  j = nys, nyn
    927                 sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i) + bc_q_t_val * dzu(nzt+1)
    928                 sk_p(nzt+3,j,i)   = sk_p(nzt+2,j,i) + bc_q_t_val * dzu(nzt+1)
    929              ENDDO
    930           ENDDO
    931 
     974                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
     975                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     976                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
     977             ENDDO
     978          ENDDO
     979
     980!
     981!--       Rain drop concentration boundary condition at the top boundary: Dirichlet
     982          DO  i = nxl, nxr
     983             DO  j = nys, nyn
     984                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
     985                sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
     986             ENDDO
     987          ENDDO
     988
     989       ELSEIF ( sk_char == 'e' )  THEN
     990
     991!
     992!--       TKE boundary condition at bottom and top boundary (generally Neumann)
     993          DO  i = nxl, nxr
     994             DO  j = nys, nyn
     995                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
     996                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     997                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
     998                sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i)
     999                sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i)
     1000             ENDDO
     1001          ENDDO
     1002
     1003       ELSE
     1004
     1005          WRITE( message_string, * ) 'no vertical boundary condi',                &
     1006                                   'tion for variable "', sk_char, '"'
     1007          CALL message( 'advec_s_bc', 'PA0158', 1, 2, 0, 6, 0 )
     1008         
    9321009       ENDIF
    9331010
    934     ELSEIF ( sk_char == 'qr' )  THEN
    935 
    936 !
    937 !--    Rain water content boundary condition at the bottom boundary:
    938 !--    Dirichlet (fixed surface rain water content).
     1011!
     1012!--    Determine the maxima of the first and second derivative in z-direction
     1013       fmax_l = 0.0_wp
    9391014       DO  i = nxl, nxr
    9401015          DO  j = nys, nyn
    941              sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
    942              sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
    943              sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
     1016             DO  k = nzb, nzt+1
     1017                numera = ABS( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k-1,j,i) )
     1018                denomi  = ABS( sk_p(k+1,j,i+1) - sk_p(k-1,j,i) )
     1019                fmax_l(1) = MAX( fmax_l(1) , numera )
     1020                fmax_l(2) = MAX( fmax_l(2) , denomi  )
     1021             ENDDO
    9441022          ENDDO
    9451023       ENDDO
    946 
    947 !
    948 !--    Rain water content boundary condition at the top boundary: Dirichlet
     1024#if defined( __parallel )
     1025       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     1026       CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr )
     1027#else
     1028       fmax = fmax_l
     1029#endif
     1030
     1031       fmax = 0.04_wp * fmax
     1032
     1033!
     1034!--    Allocate temporary arrays
     1035       ALLOCATE( a0(nzb:nzt+1,nys:nyn),   a1(nzb:nzt+1,nys:nyn),                  &
     1036                 a2(nzb:nzt+1,nys:nyn),   a12(nzb:nzt+1,nys:nyn),                 &
     1037                 a22(nzb:nzt+1,nys:nyn),  immb(nzb+1:nzt,nys:nyn),                &
     1038                 imme(nzb+1:nzt,nys:nyn), impb(nzb+1:nzt,nys:nyn),                &
     1039                 impe(nzb+1:nzt,nys:nyn), ipmb(nzb+1:nzt,nys:nyn),                &
     1040                 ipme(nzb+1:nzt,nys:nyn), ippb(nzb+1:nzt,nys:nyn),                &
     1041                 ippe(nzb+1:nzt,nys:nyn), m1(nzb-1:nzt+2,nys:nyn),                &
     1042                 sw(nzb:nzt+1,nys:nyn)                                            &
     1043               )
     1044       imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp
     1045
     1046!
     1047!--    Outer loop of all i
     1048       DO  i = nxl, nxr
     1049
     1050!
     1051!--       Compute polynomial coefficients
     1052          DO  j = nys, nyn
     1053             DO  k = nzb, nzt+1
     1054                a12(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
     1055                a22(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i)        &
     1056                                                    + sk_p(k-1,j,i) )
     1057                a0(k,j) = ( 9.0_wp * sk_p(k+2,j,i)    - 116.0_wp * sk_p(k+1,j,i)  &
     1058                            + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k-1,j,i)  &
     1059                            + 9.0_wp * sk_p(k-2,j,i) ) * f1920
     1060                a1(k,j) = ( -5.0_wp   * sk_p(k+2,j,i) + 34.0_wp * sk_p(k+1,j,i)   &
     1061                            - 34.0_wp * sk_p(k-1,j,i) + 5.0_wp  * sk_p(k-2,j,i)   &
     1062                          ) * f48
     1063                a2(k,j) = ( -3.0_wp * sk_p(k+2,j,i) + 36.0_wp * sk_p(k+1,j,i)     &
     1064                            - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k-1,j,i)     &
     1065                            - 3.0_wp * sk_p(k-2,j,i) ) * f48
     1066             ENDDO
     1067          ENDDO
     1068
     1069!
     1070!--       Fluxes using the Bott scheme
     1071!--       *VOCL LOOP,UNROLL(2)
     1072          DO  j = nys, nyn
     1073             DO  k = nzb+1, nzt
     1074                cip  =  MAX( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) )
     1075                cim  = -MIN( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) )
     1076                cipf = 1.0_wp - 2.0_wp * cip
     1077                cimf = 1.0_wp - 2.0_wp * cim
     1078                ip   =   a0(k,j)   * f2  * ( 1.0_wp - cipf )                      &
     1079                       + a1(k,j)   * f8  * ( 1.0_wp - cipf*cipf )                 &
     1080                       + a2(k,j)   * f24 * ( 1.0_wp - cipf*cipf*cipf )
     1081                im   =   a0(k+1,j) * f2  * ( 1.0_wp - cimf )                      &
     1082                       - a1(k+1,j) * f8  * ( 1.0_wp - cimf*cimf )                 &
     1083                       + a2(k+1,j) * f24 * ( 1.0_wp - cimf*cimf*cimf )
     1084                ip   = MAX( ip, 0.0_wp )
     1085                im   = MAX( im, 0.0_wp )
     1086                ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
     1087                impb(k,j) = im * MIN( 1.0_wp, sk_p(k+1,j,i) / (ip+im+1E-15_wp) )
     1088
     1089                cip  =  MAX( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) )
     1090                cim  = -MIN( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) )
     1091                cipf = 1.0_wp - 2.0_wp * cip
     1092                cimf = 1.0_wp - 2.0_wp * cim
     1093                ip   =   a0(k-1,j) * f2  * ( 1.0_wp - cipf )                      &
     1094                       + a1(k-1,j) * f8  * ( 1.0_wp - cipf*cipf )                 &
     1095                       + a2(k-1,j) * f24 * ( 1.0_wp - cipf*cipf*cipf )
     1096                im   =   a0(k,j)   * f2  * ( 1.0_wp - cimf )                      &
     1097                       - a1(k,j)   * f8  * ( 1.0_wp - cimf*cimf )                 &
     1098                       + a2(k,j)   * f24 * ( 1.0_wp - cimf*cimf*cimf )
     1099                ip   = MAX( ip, 0.0_wp )
     1100                im   = MAX( im, 0.0_wp )
     1101                ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k-1,j,i) / (ip+im+1E-15_wp) )
     1102                immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
     1103             ENDDO
     1104          ENDDO
     1105
     1106!
     1107!--       Compute monitor function m1
     1108          DO  j = nys, nyn
     1109             DO  k = nzb-1, nzt+2
     1110                m1z = ABS( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k-1,j,i) )
     1111                m1n = ABS( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
     1112                IF ( m1n /= 0.0_wp  .AND.  m1n >= m1z )  THEN
     1113                   m1(k,j) = m1z / m1n
     1114                   IF ( m1(k,j) /= 2.0_wp  .AND.  m1n < fmax(2) )  m1(k,j) = 0.0_wp
     1115                ELSEIF ( m1n < m1z )  THEN
     1116                   m1(k,j) = -1.0_wp
     1117                ELSE
     1118                   m1(k,j) = 0.0_wp
     1119                ENDIF
     1120             ENDDO
     1121          ENDDO
     1122
     1123!
     1124!--       Compute switch sw
     1125          sw = 0.0_wp
     1126          DO  j = nys, nyn
     1127             DO  k = nzb, nzt+1
     1128                m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) /                         &
     1129                     MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp )
     1130                IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) )  m2 = 0.0_wp
     1131
     1132                m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) /                         &
     1133                     MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp )
     1134                IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) )  m3 = 0.0_wp
     1135
     1136                t1 = 0.35_wp
     1137                t2 = 0.35_wp
     1138                IF ( m1(k,j) == -1.0_wp )  t2 = 0.12_wp
     1139
     1140!--             *VOCL STMT,IF(10)
     1141                IF ( m1(k-1,j) == 1.0_wp .OR. m1(k,j) == 1.0_wp                   &
     1142                     .OR. m1(k+1,j) == 1.0_wp .OR.  m2 > t2  .OR.  m3 > t2  .OR.  &
     1143                     ( m1(k,j) > t1  .AND.  m1(k-1,j) /= -1.0_wp  .AND.           &
     1144                       m1(k,j) /= -1.0_wp  .AND.  m1(k+1,j) /= -1.0_wp )          &
     1145                   )  sw(k,j) = 1.0_wp
     1146             ENDDO
     1147          ENDDO
     1148
     1149!
     1150!--       Fluxes using exponential scheme
     1151          CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
     1152          DO  j = nys, nyn
     1153             DO  k = nzb+1, nzt
     1154
     1155!--             *VOCL STMT,IF(10)
     1156                IF ( sw(k,j) == 1.0_wp )  THEN
     1157                   snenn = sk_p(k+1,j,i) - sk_p(k-1,j,i)
     1158                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
     1159                   sterm = ( sk_p(k,j,i) - sk_p(k-1,j,i) ) / snenn
     1160                   sterm = MIN( sterm, 0.9999_wp )
     1161                   sterm = MAX( sterm, 0.0001_wp )
     1162
     1163                   ix = INT( sterm * 1000 ) + 1
     1164
     1165                   cip =  MAX( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) )
     1166
     1167                   ippe(k,j) = sk_p(k-1,j,i) * cip + snenn * (                    &
     1168                               aex(ix) * cip + bex(ix) / dex(ix) * (              &
     1169                               eex(ix) -                                          &
     1170                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
     1171                                                                   )              &
     1172                                                             )
     1173                   IF ( sterm == 0.0001_wp )  ippe(k,j) = sk_p(k,j,i) * cip
     1174                   IF ( sterm == 0.9999_wp )  ippe(k,j) = sk_p(k,j,i) * cip
     1175
     1176                   snenn = sk_p(k-1,j,i) - sk_p(k+1,j,i)
     1177                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
     1178                   sterm = ( sk_p(k,j,i) - sk_p(k+1,j,i) ) / snenn
     1179                   sterm = MIN( sterm, 0.9999_wp )
     1180                   sterm = MAX( sterm, 0.0001_wp )
     1181
     1182                   ix = INT( sterm * 1000 ) + 1
     1183
     1184                   cim = -MIN( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) )
     1185
     1186                   imme(k,j) = sk_p(k+1,j,i) * cim + snenn * (                    &
     1187                               aex(ix) * cim + bex(ix) / dex(ix) * (              &
     1188                               eex(ix) -                                          &
     1189                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) &
     1190                                                                   )              &
     1191                                                             )
     1192                   IF ( sterm == 0.0001_wp )  imme(k,j) = sk_p(k,j,i) * cim
     1193                   IF ( sterm == 0.9999_wp )  imme(k,j) = sk_p(k,j,i) * cim
     1194                ENDIF
     1195
     1196!--             *VOCL STMT,IF(10)
     1197                IF ( sw(k+1,j) == 1.0_wp )  THEN
     1198                   snenn = sk_p(k,j,i) - sk_p(k+2,j,i)
     1199                   IF ( ABS( snenn ) .LT. 1E-9_wp )  snenn = 1E-9_wp
     1200                   sterm = ( sk_p(k+1,j,i) - sk_p(k+2,j,i) ) / snenn
     1201                   sterm = MIN( sterm, 0.9999_wp )
     1202                   sterm = MAX( sterm, 0.0001_wp )
     1203
     1204                   ix = INT( sterm * 1000 ) + 1
     1205
     1206                   cim = -MIN( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) )
     1207
     1208                   impe(k,j) = sk_p(k+2,j,i) * cim + snenn * (                    &
     1209                               aex(ix) * cim + bex(ix) / dex(ix) * (              &
     1210                               eex(ix) -                                          &
     1211                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
     1212                                                                   )              &
     1213                                                             )
     1214                   IF ( sterm == 0.0001_wp )  impe(k,j) = sk_p(k+1,j,i) * cim
     1215                   IF ( sterm == 0.9999_wp )  impe(k,j) = sk_p(k+1,j,i) * cim
     1216                ENDIF
     1217
     1218!--             *VOCL STMT,IF(10)
     1219                IF ( sw(k-1,j) == 1.0_wp )  THEN
     1220                   snenn = sk_p(k,j,i) - sk_p(k-2,j,i)
     1221                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
     1222                   sterm = ( sk_p(k-1,j,i) - sk_p(k-2,j,i) ) / snenn
     1223                   sterm = MIN( sterm, 0.9999_wp )
     1224                   sterm = MAX( sterm, 0.0001_wp )
     1225
     1226                   ix = INT( sterm * 1000 ) + 1
     1227
     1228                   cip = MAX( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) )
     1229
     1230                   ipme(k,j) = sk_p(k-2,j,i) * cip + snenn * (                    &
     1231                               aex(ix) * cip + bex(ix) / dex(ix) * (              &
     1232                               eex(ix) -                                          &
     1233                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
     1234                                                                   )              &
     1235                                                             )
     1236                   IF ( sterm == 0.0001_wp )  ipme(k,j) = sk_p(k-1,j,i) * cip
     1237                   IF ( sterm == 0.9999_wp )  ipme(k,j) = sk_p(k-1,j,i) * cip
     1238                ENDIF
     1239
     1240             ENDDO
     1241          ENDDO
     1242          CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
     1243
     1244!
     1245!--       Prognostic equation
     1246          DO  j = nys, nyn
     1247             DO  k = nzb+1, nzt
     1248                fplus  = ( 1.0_wp - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
     1249                       - ( 1.0_wp - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j)
     1250                fminus = ( 1.0_wp - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j) &
     1251                       - ( 1.0_wp - sw(k,j)   ) * immb(k,j) - sw(k,j)   * imme(k,j)
     1252                tendcy = fplus - fminus
     1253!
     1254!--              Removed in order to optimise speed
     1255!                ffmax   = MAX( ABS( fplus ), ABS( fminus ), 1E-35_wp )
     1256!                IF ( ( ABS( tendcy ) / ffmax ) < 1E-7_wp )  tendcy = 0.0
     1257!
     1258!--             Density correction because of possible remaining divergences
     1259                d_new = d(k,j,i) - ( w(k,j,i) - w(k-1,j,i) ) * dt_3d * ddzw(k)
     1260                sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /  &
     1261                              ( 1.0_wp + d_new )
     1262!
     1263!--             Store heat flux for subsequent statistics output.
     1264!--             array m1 is here used as temporary storage
     1265                m1(k,j) = fplus / dt_3d * dzw(k)
     1266             ENDDO
     1267          ENDDO
     1268
     1269!
     1270!--       Sum up heat flux in order to order to obtain horizontal averages
     1271          IF ( sk_char == 'pt' )  THEN
     1272             DO  sr = 0, statistic_regions
     1273                DO  j = nys, nyn
     1274                   DO  k = nzb+1, nzt
     1275                      sums_wsts_bc_l(k,sr) = sums_wsts_bc_l(k,sr) +               &
     1276                                             m1(k,j) * rmask(j,i,sr)
     1277                   ENDDO
     1278                ENDDO
     1279             ENDDO
     1280          ENDIF
     1281
     1282       ENDDO   ! End of the advection in z-direction
     1283       CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
     1284       CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'stop' )
     1285
     1286!
     1287!--    Deallocate temporary arrays
     1288       DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme,      &
     1289                   ippb, ippe, m1, sw )
     1290
     1291!
     1292!--    Store results as tendency and deallocate local array
    9491293       DO  i = nxl, nxr
    9501294          DO  j = nys, nyn
    951              sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
    952              sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
     1295             DO  k = nzb+1, nzt
     1296                tend(k,j,i) = tend(k,j,i) + ( sk_p(k,j,i) - sk(k,j,i) ) / dt_3d
     1297             ENDDO
    9531298          ENDDO
    9541299       ENDDO
    9551300
    956     ELSEIF ( sk_char == 'nr' )  THEN
    957 
    958 !
    959 !--    Rain drop concentration boundary condition at the bottom boundary:
    960 !--    Dirichlet (fixed surface rain drop concentration).
    961        DO  i = nxl, nxr
    962           DO  j = nys, nyn
    963              sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
    964              sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
    965              sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
    966           ENDDO
    967        ENDDO
    968 
    969 !
    970 !--    Rain drop concentration boundary condition at the top boundary: Dirichlet
    971        DO  i = nxl, nxr
    972           DO  j = nys, nyn
    973              sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
    974              sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
    975           ENDDO
    976        ENDDO
    977 
    978     ELSEIF ( sk_char == 'e' )  THEN
    979 
    980 !
    981 !--    TKE boundary condition at bottom and top boundary (generally Neumann)
    982        DO  i = nxl, nxr
    983           DO  j = nys, nyn
    984              sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
    985              sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
    986              sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
    987              sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i)
    988              sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i)
    989           ENDDO
    990        ENDDO
    991 
    992     ELSE
    993 
    994        WRITE( message_string, * ) 'no vertical boundary condi',                &
    995                                 'tion for variable "', sk_char, '"'
    996        CALL message( 'advec_s_bc', 'PA0158', 1, 2, 0, 6, 0 )
    997      
    998     ENDIF
    999 
    1000 !
    1001 !-- Determine the maxima of the first and second derivative in z-direction
    1002     fmax_l = 0.0_wp
    1003     DO  i = nxl, nxr
    1004        DO  j = nys, nyn
    1005           DO  k = nzb, nzt+1
    1006              numera = ABS( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k-1,j,i) )
    1007              denomi  = ABS( sk_p(k+1,j,i+1) - sk_p(k-1,j,i) )
    1008              fmax_l(1) = MAX( fmax_l(1) , numera )
    1009              fmax_l(2) = MAX( fmax_l(2) , denomi  )
    1010           ENDDO
    1011        ENDDO
    1012     ENDDO
    1013 #if defined( __parallel )
    1014     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1015     CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr )
    1016 #else
    1017     fmax = fmax_l
    1018 #endif
    1019 
    1020     fmax = 0.04_wp * fmax
    1021 
    1022 !
    1023 !-- Allocate temporary arrays
    1024     ALLOCATE( a0(nzb:nzt+1,nys:nyn),   a1(nzb:nzt+1,nys:nyn),                  &
    1025               a2(nzb:nzt+1,nys:nyn),   a12(nzb:nzt+1,nys:nyn),                 &
    1026               a22(nzb:nzt+1,nys:nyn),  immb(nzb+1:nzt,nys:nyn),                &
    1027               imme(nzb+1:nzt,nys:nyn), impb(nzb+1:nzt,nys:nyn),                &
    1028               impe(nzb+1:nzt,nys:nyn), ipmb(nzb+1:nzt,nys:nyn),                &
    1029               ipme(nzb+1:nzt,nys:nyn), ippb(nzb+1:nzt,nys:nyn),                &
    1030               ippe(nzb+1:nzt,nys:nyn), m1(nzb-1:nzt+2,nys:nyn),                &
    1031               sw(nzb:nzt+1,nys:nyn)                                            &
    1032             )
    1033     imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp
    1034 
    1035 !
    1036 !-- Outer loop of all i
    1037     DO  i = nxl, nxr
    1038 
    1039 !
    1040 !--    Compute polynomial coefficients
    1041        DO  j = nys, nyn
    1042           DO  k = nzb, nzt+1
    1043              a12(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
    1044              a22(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i)        &
    1045                                                  + sk_p(k-1,j,i) )
    1046              a0(k,j) = ( 9.0_wp * sk_p(k+2,j,i)    - 116.0_wp * sk_p(k+1,j,i)  &
    1047                          + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k-1,j,i)  &
    1048                          + 9.0_wp * sk_p(k-2,j,i) ) * f1920
    1049              a1(k,j) = ( -5.0_wp   * sk_p(k+2,j,i) + 34.0_wp * sk_p(k+1,j,i)   &
    1050                          - 34.0_wp * sk_p(k-1,j,i) + 5.0_wp  * sk_p(k-2,j,i)   &
    1051                        ) * f48
    1052              a2(k,j) = ( -3.0_wp * sk_p(k+2,j,i) + 36.0_wp * sk_p(k+1,j,i)     &
    1053                          - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k-1,j,i)     &
    1054                          - 3.0_wp * sk_p(k-2,j,i) ) * f48
    1055           ENDDO
    1056        ENDDO
    1057 
    1058 !
    1059 !--    Fluxes using the Bott scheme
    1060 !--    *VOCL LOOP,UNROLL(2)
    1061        DO  j = nys, nyn
    1062           DO  k = nzb+1, nzt
    1063              cip  =  MAX( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) )
    1064              cim  = -MIN( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) )
    1065              cipf = 1.0_wp - 2.0_wp * cip
    1066              cimf = 1.0_wp - 2.0_wp * cim
    1067              ip   =   a0(k,j)   * f2  * ( 1.0_wp - cipf )                      &
    1068                     + a1(k,j)   * f8  * ( 1.0_wp - cipf*cipf )                 &
    1069                     + a2(k,j)   * f24 * ( 1.0_wp - cipf*cipf*cipf )
    1070              im   =   a0(k+1,j) * f2  * ( 1.0_wp - cimf )                      &
    1071                     - a1(k+1,j) * f8  * ( 1.0_wp - cimf*cimf )                 &
    1072                     + a2(k+1,j) * f24 * ( 1.0_wp - cimf*cimf*cimf )
    1073              ip   = MAX( ip, 0.0_wp )
    1074              im   = MAX( im, 0.0_wp )
    1075              ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
    1076              impb(k,j) = im * MIN( 1.0_wp, sk_p(k+1,j,i) / (ip+im+1E-15_wp) )
    1077 
    1078              cip  =  MAX( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) )
    1079              cim  = -MIN( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) )
    1080              cipf = 1.0_wp - 2.0_wp * cip
    1081              cimf = 1.0_wp - 2.0_wp * cim
    1082              ip   =   a0(k-1,j) * f2  * ( 1.0_wp - cipf )                      &
    1083                     + a1(k-1,j) * f8  * ( 1.0_wp - cipf*cipf )                 &
    1084                     + a2(k-1,j) * f24 * ( 1.0_wp - cipf*cipf*cipf )
    1085              im   =   a0(k,j)   * f2  * ( 1.0_wp - cimf )                      &
    1086                     - a1(k,j)   * f8  * ( 1.0_wp - cimf*cimf )                 &
    1087                     + a2(k,j)   * f24 * ( 1.0_wp - cimf*cimf*cimf )
    1088              ip   = MAX( ip, 0.0_wp )
    1089              im   = MAX( im, 0.0_wp )
    1090              ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k-1,j,i) / (ip+im+1E-15_wp) )
    1091              immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
    1092           ENDDO
    1093        ENDDO
    1094 
    1095 !
    1096 !--    Compute monitor function m1
    1097        DO  j = nys, nyn
    1098           DO  k = nzb-1, nzt+2
    1099              m1z = ABS( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k-1,j,i) )
    1100              m1n = ABS( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
    1101              IF ( m1n /= 0.0_wp  .AND.  m1n >= m1z )  THEN
    1102                 m1(k,j) = m1z / m1n
    1103                 IF ( m1(k,j) /= 2.0_wp  .AND.  m1n < fmax(2) )  m1(k,j) = 0.0_wp
    1104              ELSEIF ( m1n < m1z )  THEN
    1105                 m1(k,j) = -1.0_wp
    1106              ELSE
    1107                 m1(k,j) = 0.0_wp
    1108              ENDIF
    1109           ENDDO
    1110        ENDDO
    1111 
    1112 !
    1113 !--    Compute switch sw
    1114        sw = 0.0_wp
    1115        DO  j = nys, nyn
    1116           DO  k = nzb, nzt+1
    1117              m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) /                         &
    1118                   MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp )
    1119              IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) )  m2 = 0.0_wp
    1120 
    1121              m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) /                         &
    1122                   MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp )
    1123              IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) )  m3 = 0.0_wp
    1124 
    1125              t1 = 0.35_wp
    1126              t2 = 0.35_wp
    1127              IF ( m1(k,j) == -1.0_wp )  t2 = 0.12_wp
    1128 
    1129 !--          *VOCL STMT,IF(10)
    1130              IF ( m1(k-1,j) == 1.0_wp .OR. m1(k,j) == 1.0_wp                   &
    1131                   .OR. m1(k+1,j) == 1.0_wp .OR.  m2 > t2  .OR.  m3 > t2  .OR.  &
    1132                   ( m1(k,j) > t1  .AND.  m1(k-1,j) /= -1.0_wp  .AND.           &
    1133                     m1(k,j) /= -1.0_wp  .AND.  m1(k+1,j) /= -1.0_wp )          &
    1134                 )  sw(k,j) = 1.0_wp
    1135           ENDDO
    1136        ENDDO
    1137 
    1138 !
    1139 !--    Fluxes using exponential scheme
    1140        CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
    1141        DO  j = nys, nyn
    1142           DO  k = nzb+1, nzt
    1143 
    1144 !--          *VOCL STMT,IF(10)
    1145              IF ( sw(k,j) == 1.0_wp )  THEN
    1146                 snenn = sk_p(k+1,j,i) - sk_p(k-1,j,i)
    1147                 IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
    1148                 sterm = ( sk_p(k,j,i) - sk_p(k-1,j,i) ) / snenn
    1149                 sterm = MIN( sterm, 0.9999_wp )
    1150                 sterm = MAX( sterm, 0.0001_wp )
    1151 
    1152                 ix = INT( sterm * 1000 ) + 1
    1153 
    1154                 cip =  MAX( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) )
    1155 
    1156                 ippe(k,j) = sk_p(k-1,j,i) * cip + snenn * (                    &
    1157                             aex(ix) * cip + bex(ix) / dex(ix) * (              &
    1158                             eex(ix) -                                          &
    1159                             EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
    1160                                                                 )              &
    1161                                                           )
    1162                 IF ( sterm == 0.0001_wp )  ippe(k,j) = sk_p(k,j,i) * cip
    1163                 IF ( sterm == 0.9999_wp )  ippe(k,j) = sk_p(k,j,i) * cip
    1164 
    1165                 snenn = sk_p(k-1,j,i) - sk_p(k+1,j,i)
    1166                 IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
    1167                 sterm = ( sk_p(k,j,i) - sk_p(k+1,j,i) ) / snenn
    1168                 sterm = MIN( sterm, 0.9999_wp )
    1169                 sterm = MAX( sterm, 0.0001_wp )
    1170 
    1171                 ix = INT( sterm * 1000 ) + 1
    1172 
    1173                 cim = -MIN( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) )
    1174 
    1175                 imme(k,j) = sk_p(k+1,j,i) * cim + snenn * (                    &
    1176                             aex(ix) * cim + bex(ix) / dex(ix) * (              &
    1177                             eex(ix) -                                          &
    1178                             EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) &
    1179                                                                 )              &
    1180                                                           )
    1181                 IF ( sterm == 0.0001_wp )  imme(k,j) = sk_p(k,j,i) * cim
    1182                 IF ( sterm == 0.9999_wp )  imme(k,j) = sk_p(k,j,i) * cim
    1183              ENDIF
    1184 
    1185 !--          *VOCL STMT,IF(10)
    1186              IF ( sw(k+1,j) == 1.0_wp )  THEN
    1187                 snenn = sk_p(k,j,i) - sk_p(k+2,j,i)
    1188                 IF ( ABS( snenn ) .LT. 1E-9_wp )  snenn = 1E-9_wp
    1189                 sterm = ( sk_p(k+1,j,i) - sk_p(k+2,j,i) ) / snenn
    1190                 sterm = MIN( sterm, 0.9999_wp )
    1191                 sterm = MAX( sterm, 0.0001_wp )
    1192 
    1193                 ix = INT( sterm * 1000 ) + 1
    1194 
    1195                 cim = -MIN( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) )
    1196 
    1197                 impe(k,j) = sk_p(k+2,j,i) * cim + snenn * (                    &
    1198                             aex(ix) * cim + bex(ix) / dex(ix) * (              &
    1199                             eex(ix) -                                          &
    1200                             EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
    1201                                                                 )              &
    1202                                                           )
    1203                 IF ( sterm == 0.0001_wp )  impe(k,j) = sk_p(k+1,j,i) * cim
    1204                 IF ( sterm == 0.9999_wp )  impe(k,j) = sk_p(k+1,j,i) * cim
    1205              ENDIF
    1206 
    1207 !--          *VOCL STMT,IF(10)
    1208              IF ( sw(k-1,j) == 1.0_wp )  THEN
    1209                 snenn = sk_p(k,j,i) - sk_p(k-2,j,i)
    1210                 IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
    1211                 sterm = ( sk_p(k-1,j,i) - sk_p(k-2,j,i) ) / snenn
    1212                 sterm = MIN( sterm, 0.9999_wp )
    1213                 sterm = MAX( sterm, 0.0001_wp )
    1214 
    1215                 ix = INT( sterm * 1000 ) + 1
    1216 
    1217                 cip = MAX( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) )
    1218 
    1219                 ipme(k,j) = sk_p(k-2,j,i) * cip + snenn * (                    &
    1220                             aex(ix) * cip + bex(ix) / dex(ix) * (              &
    1221                             eex(ix) -                                          &
    1222                             EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
    1223                                                                 )              &
    1224                                                           )
    1225                 IF ( sterm == 0.0001_wp )  ipme(k,j) = sk_p(k-1,j,i) * cip
    1226                 IF ( sterm == 0.9999_wp )  ipme(k,j) = sk_p(k-1,j,i) * cip
    1227              ENDIF
    1228 
    1229           ENDDO
    1230        ENDDO
    1231        CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
    1232 
    1233 !
    1234 !--    Prognostic equation
    1235        DO  j = nys, nyn
    1236           DO  k = nzb+1, nzt
    1237              fplus  = ( 1.0_wp - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
    1238                     - ( 1.0_wp - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j)
    1239              fminus = ( 1.0_wp - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j) &
    1240                     - ( 1.0_wp - sw(k,j)   ) * immb(k,j) - sw(k,j)   * imme(k,j)
    1241              tendcy = fplus - fminus
    1242 !
    1243 !--           Removed in order to optimise speed
    1244 !             ffmax   = MAX( ABS( fplus ), ABS( fminus ), 1E-35_wp )
    1245 !             IF ( ( ABS( tendcy ) / ffmax ) < 1E-7_wp )  tendcy = 0.0
    1246 !
    1247 !--          Density correction because of possible remaining divergences
    1248              d_new = d(k,j,i) - ( w(k,j,i) - w(k-1,j,i) ) * dt_3d * ddzw(k)
    1249              sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /  &
    1250                            ( 1.0_wp + d_new )
    1251 !
    1252 !--          Store heat flux for subsequent statistics output.
    1253 !--          array m1 is here used as temporary storage
    1254              m1(k,j) = fplus / dt_3d * dzw(k)
    1255           ENDDO
    1256        ENDDO
    1257 
    1258 !
    1259 !--    Sum up heat flux in order to order to obtain horizontal averages
    1260        IF ( sk_char == 'pt' )  THEN
    1261           DO  sr = 0, statistic_regions
    1262              DO  j = nys, nyn
    1263                 DO  k = nzb+1, nzt
    1264                    sums_wsts_bc_l(k,sr) = sums_wsts_bc_l(k,sr) +               &
    1265                                           m1(k,j) * rmask(j,i,sr)
    1266                 ENDDO
    1267              ENDDO
    1268           ENDDO
    1269        ENDIF
    1270 
    1271     ENDDO   ! End of the advection in z-direction
    1272     CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
    1273     CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'stop' )
    1274 
    1275 !
    1276 !-- Deallocate temporary arrays
    1277     DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme,      &
    1278                 ippb, ippe, m1, sw )
    1279 
    1280 !
    1281 !-- Store results as tendency and deallocate local array
    1282     DO  i = nxl, nxr
    1283        DO  j = nys, nyn
    1284           DO  k = nzb+1, nzt
    1285              tend(k,j,i) = tend(k,j,i) + ( sk_p(k,j,i) - sk(k,j,i) ) / dt_3d
    1286           ENDDO
    1287        ENDDO
    1288     ENDDO
    1289 
    1290     DEALLOCATE( sk_p )
    1291 
    1292  END SUBROUTINE advec_s_bc
     1301       DEALLOCATE( sk_p )
     1302
     1303    END SUBROUTINE advec_s_bc
     1304
     1305 END MODULE advec_s_bc_mod
Note: See TracChangeset for help on using the changeset viewer.