Ignore:
Timestamp:
Jan 17, 2017 4:38:49 PM (5 years ago)
Author:
raasch
Message:

all OpenACC directives and related parts removed from the code

File:
1 edited

Legend:

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

    r2101 r2118  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! OpenACC version of subroutines removed
    2323!
    2424! Former revisions:
     
    181181
    182182    PRIVATE
    183     PUBLIC   advec_s_ws, advec_s_ws_acc, advec_u_ws, advec_u_ws_acc,          &
    184              advec_v_ws, advec_v_ws_acc, advec_w_ws, advec_w_ws_acc,          &
    185              ws_init, ws_init_flags, ws_statistics
     183    PUBLIC   advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, ws_init,          &
     184             ws_init_flags, ws_statistics
    186185
    187186    INTERFACE ws_init
     
    207206    END INTERFACE advec_u_ws
    208207
    209     INTERFACE advec_u_ws_acc
    210        MODULE PROCEDURE advec_u_ws_acc
    211     END INTERFACE advec_u_ws_acc
    212 
    213208    INTERFACE advec_v_ws
    214209       MODULE PROCEDURE advec_v_ws
     
    216211    END INTERFACE advec_v_ws
    217212
    218     INTERFACE advec_v_ws_acc
    219        MODULE PROCEDURE advec_v_ws_acc
    220     END INTERFACE advec_v_ws_acc
    221 
    222213    INTERFACE advec_w_ws
    223214       MODULE PROCEDURE advec_w_ws
    224215       MODULE PROCEDURE advec_w_ws_ij
    225216    END INTERFACE advec_w_ws
    226 
    227     INTERFACE advec_w_ws_acc
    228        MODULE PROCEDURE advec_w_ws_acc
    229     END INTERFACE advec_w_ws_acc
    230217
    231218 CONTAINS
     
    40294016
    40304017
    4031 !------------------------------------------------------------------------------!
    4032 ! Description:
    4033 ! ------------
    4034 !> Scalar advection - Call for all grid points - accelerator version
    4035 !------------------------------------------------------------------------------!
    4036     SUBROUTINE advec_s_ws_acc ( sk, sk_char )
    4037 
    4038        USE arrays_3d,                                                         &
    4039            ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
    4040 
    4041        USE constants,                                                         &
    4042            ONLY:  adv_sca_1, adv_sca_3, adv_sca_5
    4043 
    4044        USE control_parameters,                                                &
    4045            ONLY:  intermediate_timestep_count, monotonic_adjustment, u_gtrans,&
    4046                   v_gtrans
    4047 
    4048        USE grid_variables,                                                    &
    4049            ONLY:  ddx, ddy
    4050 
    4051        USE indices,                                                           &
    4052            ONLY:  i_left, i_right, j_north, j_south, nxlg, nxrg, nyng, nysg,  &
    4053                   nzb, nzb_max, nzt, wall_flags_0
    4054 
    4055        USE kinds
    4056        
    4057 !        USE statistics,                                                       &
    4058 !            ONLY:  sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,          &
    4059 !                   sums_wsqrs_ws_l, sums_wsnrs_ws_l, weight_substep
    4060 
    4061        IMPLICIT NONE
    4062 
    4063        CHARACTER (LEN = *), INTENT(IN)    :: sk_char !<
    4064 
    4065        INTEGER(iwp) ::  i      !<
    4066        INTEGER(iwp) ::  ibit0  !<
    4067        INTEGER(iwp) ::  ibit1  !<
    4068        INTEGER(iwp) ::  ibit2  !<
    4069        INTEGER(iwp) ::  ibit3  !<
    4070        INTEGER(iwp) ::  ibit4  !<
    4071        INTEGER(iwp) ::  ibit5  !<
    4072        INTEGER(iwp) ::  ibit6  !<
    4073        INTEGER(iwp) ::  ibit7  !<
    4074        INTEGER(iwp) ::  ibit8  !<
    4075        INTEGER(iwp) ::  j      !<
    4076        INTEGER(iwp) ::  k      !<
    4077        INTEGER(iwp) ::  k_mm   !<
    4078        INTEGER(iwp) ::  k_mmm  !<
    4079        INTEGER(iwp) ::  k_pp   !<
    4080        INTEGER(iwp) ::  k_ppp  !<
    4081        INTEGER(iwp) ::  tn = 0 !<
    4082 
    4083        REAL(wp)    ::  diss_d !<
    4084        REAL(wp)    ::  diss_l !<
    4085        REAL(wp)    ::  diss_n !<
    4086        REAL(wp)    ::  diss_r !<
    4087        REAL(wp)    ::  diss_s !<
    4088        REAL(wp)    ::  diss_t !<
    4089        REAL(wp)    ::  div    !<
    4090        REAL(wp)    ::  flux_d !<
    4091        REAL(wp)    ::  flux_l !<
    4092        REAL(wp)    ::  flux_n !<
    4093        REAL(wp)    ::  flux_r !<
    4094        REAL(wp)    ::  flux_s !<
    4095        REAL(wp)    ::  flux_t !<
    4096        REAL(wp)    ::  fd_1   !<
    4097        REAL(wp)    ::  fl_1   !<
    4098        REAL(wp)    ::  fn_1   !<
    4099        REAL(wp)    ::  fr_1   !<
    4100        REAL(wp)    ::  fs_1   !<
    4101        REAL(wp)    ::  ft_1   !<
    4102        REAL(wp)    ::  phi_d  !<
    4103        REAL(wp)    ::  phi_l  !<
    4104        REAL(wp)    ::  phi_n  !<
    4105        REAL(wp)    ::  phi_r  !<
    4106        REAL(wp)    ::  phi_s  !<
    4107        REAL(wp)    ::  phi_t  !<
    4108        REAL(wp)    ::  rd     !<
    4109        REAL(wp)    ::  rl     !<
    4110        REAL(wp)    ::  rn     !<
    4111        REAL(wp)    ::  rr     !<
    4112        REAL(wp)    ::  rs     !<
    4113        REAL(wp)    ::  rt     !<
    4114        REAL(wp)    ::  u_comp !<
    4115        REAL(wp)    ::  v_comp !<
    4116 
    4117        REAL(wp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  ::  sk !<
    4118 
    4119 !
    4120 !--    Computation of fluxes and tendency terms
    4121        !$acc kernels present( ddzw, sk, tend, u, v, w, wall_flags_0 )
    4122        DO  i = i_left, i_right
    4123           DO  j = j_south, j_north
    4124              DO  k = nzb+1, nzt
    4125 
    4126                 ibit2 = IBITS(wall_flags_0(k,j,i-1),2,1)
    4127                 ibit1 = IBITS(wall_flags_0(k,j,i-1),1,1)
    4128                 ibit0 = IBITS(wall_flags_0(k,j,i-1),0,1)
    4129 
    4130                 u_comp              = u(k,j,i) - u_gtrans
    4131                 flux_l              = u_comp * (                              &
    4132                                                ( 37.0_wp * ibit2 * adv_sca_5  &
    4133                                             +     7.0_wp * ibit1 * adv_sca_3  &
    4134                                             +              ibit0 * adv_sca_1  &
    4135                                                ) *                            &
    4136                                          ( sk(k,j,i)   + sk(k,j,i-1)    )     &
    4137                                         -      (  8.0_wp * ibit2 * adv_sca_5  &
    4138                                             +              ibit1 * adv_sca_3  &
    4139                                                ) *                            &
    4140                                          ( sk(k,j,i+1) + sk(k,j,i-2)    )     &
    4141                                         +      (           ibit2 * adv_sca_5  &
    4142                                                ) *                            &
    4143                                          ( sk(k,j,i+2) + sk(k,j,i-3)    )     &
    4144                                                )
    4145 
    4146                 diss_l              = -ABS( u_comp ) * (                      &
    4147                                                ( 10.0_wp * ibit2 * adv_sca_5  &
    4148                                             +     3.0_wp * ibit1 * adv_sca_3  &
    4149                                             +              ibit0 * adv_sca_1  &
    4150                                                ) *                            &
    4151                                          ( sk(k,j,i)   - sk(k,j,i-1)    )     &
    4152                                         -      (  5.0_wp * ibit2 * adv_sca_5  &
    4153                                             +              ibit1 * adv_sca_3  &
    4154                                                ) *                            &
    4155                                          ( sk(k,j,i+1) - sk(k,j,i-2)    )     &
    4156                                         +      (           ibit2 * adv_sca_5  &
    4157                                                ) *                            &
    4158                                          ( sk(k,j,i+2) - sk(k,j,i-3)    )     &
    4159                                                         )
    4160 
    4161                 ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
    4162                 ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
    4163                 ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
    4164 
    4165                 u_comp    = u(k,j,i+1) - u_gtrans
    4166                 flux_r    = u_comp * (                                        &
    4167                           ( 37.0_wp * ibit2 * adv_sca_5                       &
    4168                       +      7.0_wp * ibit1 * adv_sca_3                       &
    4169                       +               ibit0 * adv_sca_1                       &
    4170                           ) *                                                 &
    4171                              ( sk(k,j,i+1) + sk(k,j,i)   )                    &
    4172                    -      (  8.0_wp * ibit2 * adv_sca_5                       &
    4173                        +              ibit1 * adv_sca_3                       &
    4174                           ) *                                                 &
    4175                              ( sk(k,j,i+2) + sk(k,j,i-1) )                    &
    4176                    +      (           ibit2 * adv_sca_5                       &
    4177                           ) *                                                 &
    4178                              ( sk(k,j,i+3) + sk(k,j,i-2) )                    &
    4179                                      )
    4180 
    4181                 diss_r    = -ABS( u_comp ) * (                                &
    4182                           ( 10.0_wp * ibit2 * adv_sca_5                       &
    4183                        +     3.0_wp * ibit1 * adv_sca_3                       &
    4184                        +              ibit0 * adv_sca_1                       &
    4185                           ) *                                                 &
    4186                              ( sk(k,j,i+1) - sk(k,j,i)   )                    &
    4187                    -      (  5.0_wp * ibit2 * adv_sca_5                       &
    4188                        +              ibit1 * adv_sca_3                       &
    4189                           ) *                                                 &
    4190                              ( sk(k,j,i+2) - sk(k,j,i-1) )                    &
    4191                    +      (           ibit2 * adv_sca_5                       &
    4192                           ) *                                                 &
    4193                              ( sk(k,j,i+3) - sk(k,j,i-2) )                    &
    4194                                              )
    4195 
    4196                 ibit5 = IBITS(wall_flags_0(k,j-1,i),5,1)
    4197                 ibit4 = IBITS(wall_flags_0(k,j-1,i),4,1)
    4198                 ibit3 = IBITS(wall_flags_0(k,j-1,i),3,1)
    4199 
    4200                 v_comp    = v(k,j,i) - v_gtrans
    4201                 flux_s    = v_comp * (                                        &
    4202                           ( 37.0_wp * ibit5 * adv_sca_5                       &
    4203                        +     7.0_wp * ibit4 * adv_sca_3                       &
    4204                        +              ibit3 * adv_sca_1                       &
    4205                           ) *                                                 &
    4206                              ( sk(k,j,i)  + sk(k,j-1,i)     )                 &
    4207                     -     (  8.0_wp * ibit5 * adv_sca_5                       &
    4208                        +              ibit4 * adv_sca_3                       &
    4209                           ) *                                                 &
    4210                              ( sk(k,j+1,i) + sk(k,j-2,i)    )                 &
    4211                    +      (           ibit5 * adv_sca_5                       &
    4212                           ) *                                                 &
    4213                              ( sk(k,j+2,i) + sk(k,j-3,i)    )                 &
    4214                                      )
    4215 
    4216                 diss_s    = -ABS( v_comp ) * (                                &
    4217                           ( 10.0_wp * ibit5 * adv_sca_5                       &
    4218                        +     3.0_wp * ibit4 * adv_sca_3                       &
    4219                        +              ibit3 * adv_sca_1                       &
    4220                           ) *                                                 &
    4221                              ( sk(k,j,i)   - sk(k,j-1,i)  )                   &
    4222                    -      (  5.0_wp * ibit5 * adv_sca_5                       &
    4223                        +              ibit4 * adv_sca_3                       &
    4224                           ) *                                                 &
    4225                              ( sk(k,j+1,i) - sk(k,j-2,i)  )                   &
    4226                    +      (           ibit5 * adv_sca_5                       &
    4227                           ) *                                                 &
    4228                              ( sk(k,j+2,i) - sk(k,j-3,i)  )                   &
    4229                                              )
    4230 
    4231                 ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
    4232                 ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
    4233                 ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
    4234 
    4235                 v_comp    = v(k,j+1,i) - v_gtrans
    4236                 flux_n    = v_comp * (                                        &
    4237                           ( 37.0_wp * ibit5 * adv_sca_5                       &
    4238                        +     7.0_wp * ibit4 * adv_sca_3                       &
    4239                        +              ibit3 * adv_sca_1                       &
    4240                           ) *                                                 &
    4241                              ( sk(k,j+1,i) + sk(k,j,i)   )                    &
    4242                    -      (  8.0_wp * ibit5 * adv_sca_5                       &
    4243                        +              ibit4 * adv_sca_3                       &
    4244                           ) *                                                 &
    4245                              ( sk(k,j+2,i) + sk(k,j-1,i) )                    &
    4246                    +      (           ibit5 * adv_sca_5                       &
    4247                           ) *                                                 &
    4248                              ( sk(k,j+3,i) + sk(k,j-2,i) )                    &
    4249                                      )
    4250 
    4251                 diss_n    = -ABS( v_comp ) * (                                &
    4252                           ( 10.0_wp * ibit5 * adv_sca_5                       &
    4253                        +     3.0_wp * ibit4 * adv_sca_3                       &
    4254                        +              ibit3 * adv_sca_1                       &
    4255                           ) *                                                 &
    4256                              ( sk(k,j+1,i) - sk(k,j,i)    )                   &
    4257                    -      (  5.0_wp * ibit5 * adv_sca_5                       &
    4258                        +              ibit4 * adv_sca_3                       &
    4259                           ) *                                                 &
    4260                              ( sk(k,j+2,i) - sk(k,j-1,i)  )                   &
    4261                    +      (           ibit5 * adv_sca_5                       &
    4262                           ) *                                                 &
    4263                              ( sk(k,j+3,i) - sk(k,j-2,i)  )                   &
    4264                                              )
    4265 
    4266 !
    4267 !--             indizes k_m, k_mm, ... should be known at these point
    4268                 ibit8 = IBITS(wall_flags_0(k-1,j,i),8,1)
    4269                 ibit7 = IBITS(wall_flags_0(k-1,j,i),7,1)
    4270                 ibit6 = IBITS(wall_flags_0(k-1,j,i),6,1)
    4271 
    4272                 k_pp  = k + 2 * ibit8
    4273                 k_mm  = k - 2 * ( ibit7 + ibit8 )
    4274                 k_mmm = k - 3 * ibit8
    4275 
    4276                 flux_d    = w(k-1,j,i) * (                                    &
    4277                            ( 37.0_wp * ibit8 * adv_sca_5                      &
    4278                         +     7.0_wp * ibit7 * adv_sca_3                      &
    4279                         +              ibit6 * adv_sca_1                      &
    4280                            ) *                                                &
    4281                                    ( sk(k,j,i)    + sk(k-1,j,i)  )            &
    4282                     -      (  8.0_wp * ibit8 * adv_sca_5                      &
    4283                           +            ibit7 * adv_sca_3                      &
    4284                            ) *                                                &
    4285                                    ( sk(k+1,j,i) + sk(k_mm,j,i)  )            &
    4286                     +      (           ibit8 * adv_sca_5                      &
    4287                            ) *     ( sk(k_pp,j,i)+ sk(k_mmm,j,i) )            &
    4288                                          )
    4289 
    4290                 diss_d    = -ABS( w(k-1,j,i) ) * (                            &
    4291                            ( 10.0_wp * ibit8 * adv_sca_5                      &
    4292                         +     3.0_wp * ibit7 * adv_sca_3                      &
    4293                         +              ibit6 * adv_sca_1                      &
    4294                            ) *                                                &
    4295                                    ( sk(k,j,i)    - sk(k-1,j,i)   )           &
    4296                     -      (  5.0_wp * ibit8 * adv_sca_5                      &
    4297                         +              ibit7 * adv_sca_3                      &
    4298                            ) *                                                &
    4299                                    ( sk(k+1,j,i)  - sk(k_mm,j,i)  )           &
    4300                     +      (           ibit8 * adv_sca_5                      &
    4301                            ) *                                                &
    4302                                    ( sk(k_pp,j,i) - sk(k_mmm,j,i) )           &
    4303                                                  )
    4304 
    4305                 ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
    4306                 ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
    4307                 ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
    4308 
    4309                 k_ppp = k + 3 * ibit8
    4310                 k_pp  = k + 2 * ( 1 - ibit6  )
    4311                 k_mm  = k - 2 * ibit8
    4312 
    4313                 flux_t    = w(k,j,i) * rho_air_zw(k) * (                      &
    4314                            ( 37.0_wp * ibit8 * adv_sca_5                      &
    4315                         +     7.0_wp * ibit7 * adv_sca_3                      &
    4316                         +              ibit6 * adv_sca_1                      &
    4317                            ) *                                                &
    4318                                    ( sk(k+1,j,i)  + sk(k,j,i)    )            &
    4319                     -      (  8.0_wp * ibit8 * adv_sca_5                      &
    4320                         +              ibit7 * adv_sca_3                      &
    4321                            ) *                                                &
    4322                                    ( sk(k_pp,j,i) + sk(k-1,j,i)  )            &
    4323                     +      (           ibit8 * adv_sca_5                      &
    4324                            ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )            &
    4325                                        )
    4326 
    4327                 diss_t    = -ABS( w(k,j,i) ) * rho_air_zw(k) * (              &
    4328                            ( 10.0_wp * ibit8 * adv_sca_5                      &
    4329                         +     3.0_wp * ibit7 * adv_sca_3                      &
    4330                         +              ibit6 * adv_sca_1                      &
    4331                            ) *                                                &
    4332                                    ( sk(k+1,j,i)   - sk(k,j,i)    )           &
    4333                     -      (  5.0_wp * ibit8 * adv_sca_5                      &
    4334                         +              ibit7 * adv_sca_3                      &
    4335                            ) *                                                &
    4336                                    ( sk(k_pp,j,i)  - sk(k-1,j,i)  )           &
    4337                     +      (           ibit8 * adv_sca_5                      &
    4338                            ) *                                                &
    4339                                    ( sk(k_ppp,j,i) - sk(k_mm,j,i) )           &
    4340                                          )
    4341 !
    4342 !--             Apply monotonic adjustment.
    4343                 IF ( monotonic_adjustment )  THEN
    4344 !
    4345 !--                At first, calculate first order fluxes.
    4346                    u_comp = u(k,j,i) - u_gtrans
    4347                    fl_1   =  ( u_comp   * ( sk(k,j,i) + sk(k,j,i-1) )         &
    4348                          -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) )         &
    4349                              ) * adv_sca_1
    4350 
    4351                    u_comp = u(k,j,i+1) - u_gtrans
    4352                    fr_1   =  ( u_comp   * ( sk(k,j,i+1) + sk(k,j,i) )         &
    4353                          -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) )         &
    4354                              ) * adv_sca_1
    4355 
    4356                    v_comp = v(k,j,i) - v_gtrans
    4357                    fs_1   =  ( v_comp   * ( sk(k,j,i) + sk(k,j-1,i) )         &
    4358                          -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) )         &
    4359                              ) * adv_sca_1
    4360 
    4361                    v_comp = v(k,j+1,i) - v_gtrans
    4362                    fn_1   =  ( v_comp   * ( sk(k,j+1,i) + sk(k,j,i) )         &
    4363                          -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) )         &
    4364                              ) * adv_sca_1
    4365 
    4366                    fd_1   = (  w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )     &
    4367                         -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )      &
    4368                             ) * adv_sca_1 * rho_air_zw(k)
    4369 
    4370                    ft_1   = (  w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )       &
    4371                         -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )        &
    4372                             ) * adv_sca_1 * rho_air_zw(k)
    4373 !
    4374 !--                Calculate ratio of upwind gradients. Note, Min/Max is just
    4375 !--                to avoid if statements.
    4376                    rl     = ( MAX( 0.0_wp, u(k,j,i) - u_gtrans ) *            &
    4377                                ABS( ( sk(k,j,i-1) - sk(k,j,i-2)            ) /&
    4378                                     ( sk(k,j,i)   - sk(k,j,i-1) + 1E-20_wp )  &
    4379                                   ) +                                         &
    4380                               MIN( 0.0_wp, u(k,j,i) - u_gtrans ) *            &
    4381                                ABS( ( sk(k,j,i)   - sk(k,j,i+1)            ) /&
    4382                                     ( sk(k,j,i-1) - sk(k,j,i)   + 1E-20_wp )  &
    4383                                   )                                           &
    4384                             ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )
    4385 
    4386                    rr     = ( MAX( 0.0_wp, u(k,j,i+1) - u_gtrans ) *          &
    4387                                ABS( ( sk(k,j,i)   - sk(k,j,i-1)            ) /&
    4388                                     ( sk(k,j,i+1) - sk(k,j,i)   + 1E-20_wp )  &
    4389                                   ) +                                         &
    4390                               MIN( 0.0_wp, u(k,j,i+1) - u_gtrans ) *          &
    4391                                ABS( ( sk(k,j,i+1) - sk(k,j,i+2)            ) /&
    4392                                     ( sk(k,j,i)   - sk(k,j,i+1) + 1E-20_wp )  &
    4393                                   )                                           &
    4394                             ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )
    4395 
    4396                    rs     = ( MAX( 0.0_wp, v(k,j,i) - v_gtrans ) *            &
    4397                                ABS( ( sk(k,j-1,i) - sk(k,j-2,i)            ) /&
    4398                                     ( sk(k,j,i)   - sk(k,j-1,i) + 1E-20_wp )  &
    4399                                   ) +                                         &
    4400                               MIN( 0.0_wp, v(k,j,i) - v_gtrans ) *            &
    4401                                ABS( ( sk(k,j,i)   - sk(k,j+1,i)            ) /&
    4402                                     ( sk(k,j-1,i) - sk(k,j,i)   + 1E-20_wp )  &
    4403                                   )                                           &
    4404                             ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )
    4405 
    4406                    rn     = ( MAX( 0.0_wp, v(k,j+1,i) - v_gtrans ) *          &
    4407                                ABS( ( sk(k,j,i)   - sk(k,j-1,i)            ) /&
    4408                                     ( sk(k,j+1,i) - sk(k,j,i)   + 1E-20_wp )  &
    4409                                   ) +                                         &
    4410                               MIN( 0.0_wp, v(k,j+1,i) - v_gtrans ) *          &
    4411                                ABS( ( sk(k,j+1,i) - sk(k,j+2,i)            ) /&
    4412                                     ( sk(k,j,i)   - sk(k,j+1,i) + 1E-20_wp )  &
    4413                                   )                                           &
    4414                             ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )     
    4415 !   
    4416 !--                Reuse k_mm and compute k_mmm for the vertical gradient ratios.
    4417 !--                Note, for vertical advection below the third grid point above
    4418 !--                surface ( or below the model top) rd and rt are set to 0, i.e.
    4419 !--                use of first order scheme is enforced.
    4420                    k_mmm  = k - 3 * ibit8
    4421 
    4422                    rd     = ( MAX( 0.0_wp, w(k-1,j,i) ) *                     &
    4423                             ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i)           ) / &
    4424                                  ( sk(k-1,j,i)  - sk(k_mm,j,i) + 1E-20_wp )   &
    4425                                ) +                                            &
    4426                               MIN( 0.0_wp, w(k-1,j,i) ) *                     &
    4427                             ABS( ( sk(k-1,j,i) - sk(k,j,i)            ) /     &
    4428                                  ( sk(k_mm,j,i) - sk(k-1,j,i)   + 1E-20_wp )  &
    4429                                )                                              &
    4430                             ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp )
    4431  
    4432                    rt     = ( MAX( 0.0_wp, w(k,j,i) ) *                       &
    4433                             ABS( ( sk(k,j,i)   - sk(k-1,j,i)            ) /   &
    4434                                  ( sk(k+1,j,i) - sk(k,j,i)   + 1E-20_wp )     &
    4435                                ) +                                            &
    4436                               MIN( 0.0_wp, w(k,j,i) ) *                       &
    4437                             ABS( ( sk(k+1,j,i) - sk(k_pp,j,i)           ) /   &
    4438                                  ( sk(k,j,i)   - sk(k+1,j,i) + 1E-20_wp )     &
    4439                                )                                              &
    4440                             ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp )
    4441 !
    4442 !--                Calculate empirical limiter function (van Albada2 limiter).
    4443                    phi_l = MIN( 1.0_wp, ( 2.0_wp * ABS( rl ) ) /              &
    4444                                         ( rl**2 + 1.0_wp ) )
    4445                    phi_r = MIN( 1.0_wp, ( 2.0_wp * ABS( rr ) ) /              &
    4446                                         ( rr**2 + 1.0_wp ) )
    4447                    phi_s = MIN( 1.0_wp, ( 2.0_wp * ABS( rs ) ) /              &
    4448                                         ( rs**2 + 1.0_wp ) )
    4449                    phi_n = MIN( 1.0_wp, ( 2.0_wp * ABS( rn ) ) /              &
    4450                                         ( rn**2 + 1.0_wp ) )
    4451                    phi_d = MIN( 1.0_wp, ( 2.0_wp * ABS( rd ) ) /              &
    4452                                         ( rd**2 + 1.0_wp ) )
    4453                    phi_t = MIN( 1.0_wp, ( 2.0_wp * ABS( rt ) ) /              &
    4454                                         ( rt**2 + 1.0_wp ) )
    4455 !
    4456 !--                Calculate the resulting monotone flux.
    4457                    flux_l = fl_1 - phi_l * ( fl_1 - flux_l )
    4458                    flux_r = fr_1 - phi_r * ( fr_1 - flux_r )
    4459                    flux_s = fs_1 - phi_s * ( fs_1 - flux_s )
    4460                    flux_n = fn_1 - phi_n * ( fn_1 - flux_n )
    4461                    flux_d = fd_1 - phi_d * ( fd_1 - flux_d )
    4462                    flux_t = ft_1 - phi_t * ( ft_1 - flux_t )
    4463 !
    4464 !--                Moreover, modify dissipation flux according to the limiter.
    4465                    diss_l = diss_l * phi_l
    4466                    diss_r = diss_r * phi_r
    4467                    diss_s = diss_s * phi_s
    4468                    diss_n = diss_n * phi_n
    4469                    diss_d = diss_d * phi_d
    4470                    diss_t = diss_t * phi_t
    4471 
    4472                 ENDIF
    4473 !
    4474 !--             Calculate the divergence of the velocity field. A respective
    4475 !--             correction is needed to overcome numerical instabilities caused
    4476 !--             by a not sufficient reduction of divergences near topography.
    4477                 div   =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )             &
    4478                           - u(k,j,i)   * ( IBITS(wall_flags_0(k,j,i-1),0,1)    &
    4479                                          + IBITS(wall_flags_0(k,j,i-1),1,1)    &
    4480                                          + IBITS(wall_flags_0(k,j,i-1),2,1)    &
    4481                                          )                                     &
    4482                           ) * rho_air(k) * ddx                                 &
    4483                         + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
    4484                           - v(k,j,i)   * ( IBITS(wall_flags_0(k,j-1,i),3,1)    &
    4485                                          + IBITS(wall_flags_0(k,j-1,i),4,1)    &
    4486                                          + IBITS(wall_flags_0(k,j-1,i),5,1)    &
    4487                                          )                                     &
    4488                           ) * rho_air(k) * ddy                                 &
    4489                         + ( w(k,j,i) * rho_air_zw(k) *                         &
    4490                                          ( ibit6 + ibit7 + ibit8 )             &
    4491                           - w(k-1,j,i) * rho_air_zw(k-1) *                     &
    4492                                          ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
    4493                                          + IBITS(wall_flags_0(k-1,j,i),7,1)    &
    4494                                          + IBITS(wall_flags_0(k-1,j,i),8,1)    &
    4495                                          )                                     &     
    4496                           ) * ddzw(k)
    4497 
    4498 
    4499                 tend(k,j,i) = - (                                             &
    4500                                ( flux_r + diss_r - flux_l - diss_l ) * ddx    &
    4501                              + ( flux_n + diss_n - flux_s - diss_s ) * ddy    &
    4502                              + ( ( flux_t + diss_t ) -                        &
    4503                                  ( flux_d + diss_d )                          &
    4504                                                     ) * drho_air(k) * ddzw(k) &
    4505                                 ) + div * sk(k,j,i)
    4506 
    4507 !++
    4508 !--             Evaluation of statistics
    4509 !                SELECT CASE ( sk_char )
    4510 !
    4511 !                   CASE ( 'pt' )
    4512 !                      sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn)         &
    4513 !                       + ( flux_t + diss_t )                                &
    4514 !                       *   weight_substep(intermediate_timestep_count)
    4515 !                   CASE ( 'sa' )
    4516 !                      sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn)         &
    4517 !                       + ( flux_t + diss_t )                                &
    4518 !                       *   weight_substep(intermediate_timestep_count)
    4519 !                   CASE ( 'q' )
    4520 !                      sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn)           &
    4521 !                      + ( flux_t + diss_t )                                 &
    4522 !                      *   weight_substep(intermediate_timestep_count)
    4523 !                   CASE ( 'qr' )
    4524 !                      sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn)         &
    4525 !                      + ( flux_t + diss_t )                                 &
    4526 !                      *   weight_substep(intermediate_timestep_count)
    4527 !                   CASE ( 'nr' )
    4528 !                      sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn)         &
    4529 !                      + ( flux_t + diss_t )                                 &
    4530 !                      *   weight_substep(intermediate_timestep_count)
    4531 !
    4532 !                END SELECT
    4533 
    4534              ENDDO
    4535          ENDDO
    4536       ENDDO
    4537       !$acc end kernels
    4538 
    4539     END SUBROUTINE advec_s_ws_acc
    4540 
    45414018
    45424019!------------------------------------------------------------------------------!
     
    50394516! Description:
    50404517! ------------
    5041 !> Advection of u - Call for all grid points - accelerator version
    5042 !------------------------------------------------------------------------------!
    5043     SUBROUTINE advec_u_ws_acc
    5044 
    5045        USE arrays_3d,                                                          &
    5046            ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
    5047 
    5048        USE constants,                                                          &
    5049            ONLY:  adv_mom_1, adv_mom_3, adv_mom_5
    5050 
    5051        USE control_parameters,                                                 &
    5052            ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
    5053 
    5054        USE grid_variables,                                                     &
    5055            ONLY:  ddx, ddy
    5056 
    5057        USE indices,                                                            &
    5058            ONLY:  i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb,  &
    5059                   nzb_max, nzt, wall_flags_0
    5060            
    5061        USE kinds
    5062        
    5063 !        USE statistics,                                                       &
    5064 !            ONLY:  hom, sums_us2_ws_l, sums_wsus_ws_l, weight_substep
    5065 
    5066        IMPLICIT NONE
    5067 
    5068        INTEGER(iwp) ::  i      !<
    5069        INTEGER(iwp) ::  ibit9  !<
    5070        INTEGER(iwp) ::  ibit10 !<
    5071        INTEGER(iwp) ::  ibit11 !<
    5072        INTEGER(iwp) ::  ibit12 !<
    5073        INTEGER(iwp) ::  ibit13 !<
    5074        INTEGER(iwp) ::  ibit14 !<
    5075        INTEGER(iwp) ::  ibit15 !<
    5076        INTEGER(iwp) ::  ibit16 !<
    5077        INTEGER(iwp) ::  ibit17 !<
    5078        INTEGER(iwp) ::  j      !<
    5079        INTEGER(iwp) ::  k      !<
    5080        INTEGER(iwp) ::  k_mmm  !<
    5081        INTEGER(iwp) ::  k_mm   !<
    5082        INTEGER(iwp) ::  k_pp   !<
    5083        INTEGER(iwp) ::  k_ppp  !<
    5084        INTEGER(iwp) ::  tn = 0 !<
    5085 
    5086        REAL(wp)    ::  diss_d   !<
    5087        REAL(wp)    ::  diss_l   !<
    5088        REAL(wp)    ::  diss_n   !<
    5089        REAL(wp)    ::  diss_r   !<
    5090        REAL(wp)    ::  diss_s   !<
    5091        REAL(wp)    ::  diss_t   !<
    5092        REAL(wp)    ::  div      !<
    5093        REAL(wp)    ::  flux_d   !<
    5094        REAL(wp)    ::  flux_l   !<
    5095        REAL(wp)    ::  flux_n   !<
    5096        REAL(wp)    ::  flux_r   !<
    5097        REAL(wp)    ::  flux_s   !<
    5098        REAL(wp)    ::  flux_t   !<
    5099        REAL(wp)    ::  gu       !<
    5100        REAL(wp)    ::  gv       !<
    5101        REAL(wp)    ::  u_comp   !<
    5102        REAL(wp)    ::  u_comp_l !<
    5103        REAL(wp)    ::  v_comp   !<
    5104        REAL(wp)    ::  v_comp_s !<
    5105        REAL(wp)    ::  w_comp   !<
    5106 
    5107 
    5108        gu = 2.0_wp * u_gtrans
    5109        gv = 2.0_wp * v_gtrans
    5110 
    5111 !
    5112 !--    Computation of fluxes and tendency terms
    5113        !$acc  kernels present( ddzw, tend, u, v, w, wall_flags_0 )
    5114        DO i = i_left, i_right
    5115           DO  j = j_south, j_north
    5116              DO  k = nzb+1, nzt
    5117 
    5118                 ibit11 = IBITS(wall_flags_0(k,j,i-1),11,1)
    5119                 ibit10 = IBITS(wall_flags_0(k,j,i-1),10,1)
    5120                 ibit9  = IBITS(wall_flags_0(k,j,i-1),9,1)
    5121 
    5122                 u_comp_l           = u(k,j,i) + u(k,j,i-1) - gu
    5123                 flux_l             = u_comp_l * (                          &
    5124                                     ( 37.0_wp * ibit11 * adv_mom_5             &
    5125                                  +     7.0_wp * ibit10 * adv_mom_3             &
    5126                                  +              ibit9  * adv_mom_1             &
    5127                                     ) *                                     &
    5128                                   ( u(k,j,i)   + u(k,j,i-1) )               &
    5129                              -      (  8.0_wp * ibit11 * adv_mom_5             &
    5130                                  +              ibit10 * adv_mom_3             &
    5131                                     ) *                                     &
    5132                                   ( u(k,j,i+1) + u(k,j,i-2) )               &
    5133                              +      (           ibit11 * adv_mom_5             &
    5134                                     ) *                                     &
    5135                                   ( u(k,j,i+2) + u(k,j,i-3) )               &
    5136                                                 )
    5137 
    5138                 diss_l             = - ABS( u_comp_l ) * (                &
    5139                                    ( 10.0_wp * ibit11 * adv_mom_5             &
    5140                                 +     3.0_wp * ibit10 * adv_mom_3             &
    5141                                 +              ibit9  * adv_mom_1             &
    5142                                    ) *                                     &
    5143                                  ( u(k,j,i)   - u(k,j,i-1) )               &
    5144                             -      (  5.0_wp * ibit11 * adv_mom_5             &
    5145                                 +              ibit10 * adv_mom_3             &
    5146                                    ) *                                     &
    5147                                  ( u(k,j,i+1) - u(k,j,i-2) )               &
    5148                             +      (           ibit11 * adv_mom_5             &
    5149                                    ) *                                     &
    5150                                  ( u(k,j,i+2) - u(k,j,i-3) )               &
    5151                                                          )
    5152 
    5153                 ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
    5154                 ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
    5155                 ibit9  = IBITS(wall_flags_0(k,j,i),9,1)
    5156 
    5157                 u_comp    = u(k,j,i+1) + u(k,j,i)
    5158                 flux_r    = ( u_comp   - gu ) * (                           &
    5159                           ( 37.0_wp * ibit11 * adv_mom_5                        &
    5160                        +     7.0_wp * ibit10 * adv_mom_3                        &
    5161                        +              ibit9  * adv_mom_1                        &
    5162                           ) *                                                &
    5163                                  ( u(k,j,i+1) + u(k,j,i)   )                 &
    5164                    -      (  8.0_wp * ibit11 * adv_mom_5                        &
    5165                        +              ibit10 * adv_mom_3                        &
    5166                           ) *                                                &
    5167                                  ( u(k,j,i+2) + u(k,j,i-1) )                 &
    5168                    +      (           ibit11 * adv_mom_5                        &
    5169                           ) *                                                &
    5170                                  ( u(k,j,i+3) + u(k,j,i-2) )                 &
    5171                                                  )
    5172 
    5173                 diss_r    = - ABS( u_comp    - gu ) * (                      &
    5174                           ( 10.0_wp * ibit11 * adv_mom_5                        &
    5175                        +     3.0_wp * ibit10 * adv_mom_3                        &
    5176                        +              ibit9  * adv_mom_1                        &
    5177                           ) *                                                &
    5178                                  ( u(k,j,i+1) - u(k,j,i)  )                  &
    5179                    -      (  5.0_wp * ibit11 * adv_mom_5                        &
    5180                        +              ibit10 * adv_mom_3                        &
    5181                           ) *                                                &
    5182                                  ( u(k,j,i+2) - u(k,j,i-1) )                 &
    5183                    +      (           ibit11 * adv_mom_5                        &
    5184                           ) *                                                &
    5185                                  ( u(k,j,i+3) - u(k,j,i-2) )                 &
    5186                                                      )
    5187 
    5188                 ibit14 = IBITS(wall_flags_0(k,j-1,i),14,1)
    5189                 ibit13 = IBITS(wall_flags_0(k,j-1,i),13,1)
    5190                 ibit12 = IBITS(wall_flags_0(k,j-1,i),12,1)
    5191 
    5192                 v_comp_s                 = v(k,j,i) + v(k,j,i-1) - gv
    5193                 flux_s                   = v_comp_s * (                       &
    5194                                    ( 37.0_wp * ibit14 * adv_mom_5                &
    5195                                 +     7.0_wp * ibit13 * adv_mom_3                &
    5196                                 +              ibit12 * adv_mom_1                &
    5197                                    ) *                                         &
    5198                                      ( u(k,j,i)   + u(k,j-1,i) )              &
    5199                             -      (  8.0_wp * ibit14 * adv_mom_5                &
    5200                             +                  ibit13 * adv_mom_3                    &
    5201                                    ) *                                        &
    5202                                      ( u(k,j+1,i) + u(k,j-2,i) )              &
    5203                         +      (               ibit14 * adv_mom_5                    &
    5204                                ) *                                            &
    5205                                      ( u(k,j+2,i) + u(k,j-3,i) )              &
    5206                                                )
    5207 
    5208                 diss_s                  = - ABS ( v_comp_s ) * (              &
    5209                                    ( 10.0_wp * ibit14 * adv_mom_5                &
    5210                                 +     3.0_wp * ibit13 * adv_mom_3                &
    5211                                 +              ibit12 * adv_mom_1                &
    5212                                    ) *                                        &
    5213                                      ( u(k,j,i)   - u(k,j-1,i) )              &
    5214                             -      (  5.0_wp * ibit14 * adv_mom_5                &
    5215                                 +              ibit13 * adv_mom_3                &
    5216                                    ) *                                        &
    5217                                      ( u(k,j+1,i) - u(k,j-2,i) )              &
    5218                             +      (           ibit14 * adv_mom_5                &
    5219                                    ) *                                        &
    5220                                      ( u(k,j+2,i) - u(k,j-3,i) )              &
    5221                                                          )
    5222 
    5223 
    5224                 ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
    5225                 ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
    5226                 ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
    5227 
    5228                 v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    5229                 flux_n    = v_comp * (                                       &
    5230                           ( 37.0_wp * ibit14 * adv_mom_5                        &
    5231                        +     7.0_wp * ibit13 * adv_mom_3                        &
    5232                        +              ibit12 * adv_mom_1                        &
    5233                           ) *                                                &
    5234                                  ( u(k,j+1,i) + u(k,j,i)   )                 &
    5235                    -      (  8.0_wp * ibit14 * adv_mom_5                        &
    5236                        +              ibit13 * adv_mom_3                        &
    5237                           ) *                                                &
    5238                                  ( u(k,j+2,i) + u(k,j-1,i) )                 &
    5239                    +      (           ibit14 * adv_mom_5                        &
    5240                           ) *                                                &
    5241                                  ( u(k,j+3,i) + u(k,j-2,i) )                 &
    5242                                                  )
    5243 
    5244                 diss_n    = - ABS ( v_comp ) * (                             &
    5245                           ( 10.0_wp * ibit14 * adv_mom_5                        &
    5246                        +     3.0_wp * ibit13 * adv_mom_3                        &
    5247                        +              ibit12 * adv_mom_1                        &
    5248                           ) *                                                &
    5249                                  ( u(k,j+1,i) - u(k,j,i)  )                  &
    5250                    -      (  5.0_wp * ibit14 * adv_mom_5                        &
    5251                        +              ibit13 * adv_mom_3                        &
    5252                           ) *                                                &
    5253                                  ( u(k,j+2,i) - u(k,j-1,i) )                 &
    5254                    +      (           ibit14 * adv_mom_5                        &
    5255                           ) *                                                &
    5256                                  ( u(k,j+3,i) - u(k,j-2,i) )                 &
    5257                                                       )
    5258 
    5259                 ibit17 = IBITS(wall_flags_0(k-1,j,i),17,1)
    5260                 ibit16 = IBITS(wall_flags_0(k-1,j,i),16,1)
    5261                 ibit15 = IBITS(wall_flags_0(k-1,j,i),15,1)
    5262 
    5263                 k_pp  = k + 2 * ibit17
    5264                 k_mm  = k - 2 * ( ibit16 + ibit17 )
    5265                 k_mmm = k - 3 * ibit17
    5266 
    5267                 w_comp    = w(k-1,j,i) + w(k-1,j,i-1)
    5268                 flux_d    = w_comp  * (                                      &
    5269                           ( 37.0_wp * ibit17 * adv_mom_5                        &
    5270                        +     7.0_wp * ibit16 * adv_mom_3                        &
    5271                        +              ibit15 * adv_mom_1                        &
    5272                           ) *                                                &
    5273                              ( u(k,j,i)    + u(k-1,j,i)   )                  &
    5274                    -      (  8.0_wp * ibit17 * adv_mom_5                        &
    5275                        +              ibit16 * adv_mom_3                        &
    5276                           ) *                                                &
    5277                              ( u(k+1,j,i) + u(k_mm,j,i)   )                  &
    5278                    +      (           ibit17 * adv_mom_5                        &
    5279                           ) *                                                 &
    5280                              ( u(k_pp,j,i) + u(k_mmm,j,i) )                  &
    5281                                       )
    5282 
    5283                 diss_d    = - ABS( w_comp ) * (                              &
    5284                           ( 10.0_wp * ibit17 * adv_mom_5                        &
    5285                        +     3.0_wp * ibit16 * adv_mom_3                        &
    5286                        +              ibit15 * adv_mom_1                        &
    5287                           ) *                                                &
    5288                              ( u(k,j,i)     - u(k-1,j,i)  )                  &
    5289                    -      (  5.0_wp * ibit17 * adv_mom_5                        &
    5290                        +              ibit16 * adv_mom_3                        &
    5291                           ) *                                                &
    5292                              ( u(k+1,j,i)  - u(k_mm,j,i)  )                  &
    5293                    +      (           ibit17 * adv_mom_5                        &
    5294                            ) *                                               &
    5295                              ( u(k_pp,j,i) - u(k_mmm,j,i) )                  &
    5296                                               )
    5297 !
    5298 !--             k index has to be modified near bottom and top, else array
    5299 !--             subscripts will be exceeded.
    5300                 ibit17 = IBITS(wall_flags_0(k,j,i),17,1)
    5301                 ibit16 = IBITS(wall_flags_0(k,j,i),16,1)
    5302                 ibit15 = IBITS(wall_flags_0(k,j,i),15,1)
    5303 
    5304                 k_ppp = k + 3 * ibit17
    5305                 k_pp  = k + 2 * ( 1 - ibit15  )
    5306                 k_mm  = k - 2 * ibit17
    5307 
    5308                 w_comp    = w(k,j,i) + w(k,j,i-1)
    5309                 flux_t    = w_comp * rho_air_zw(k) * (                       &
    5310                           ( 37.0_wp * ibit17 * adv_mom_5                        &
    5311                        +     7.0_wp * ibit16 * adv_mom_3                        &
    5312                        +              ibit15 * adv_mom_1                        &
    5313                           ) *                                                &
    5314                              ( u(k+1,j,i)  + u(k,j,i)     )                  &
    5315                    -      (  8.0_wp * ibit17 * adv_mom_5                        &
    5316                        +              ibit16 * adv_mom_3                        &
    5317                           ) *                                                &
    5318                              ( u(k_pp,j,i) + u(k-1,j,i)   )                  &
    5319                    +      (           ibit17 * adv_mom_5                        &
    5320                           ) *                                                &
    5321                              ( u(k_ppp,j,i) + u(k_mm,j,i) )                  &
    5322                                       )
    5323 
    5324                 diss_t    = - ABS( w_comp ) * rho_air_zw(k) * (              &
    5325                           ( 10.0_wp * ibit17 * adv_mom_5                        &
    5326                        +     3.0_wp * ibit16 * adv_mom_3                        &
    5327                        +              ibit15 * adv_mom_1                        &
    5328                           ) *                                                &
    5329                              ( u(k+1,j,i)   - u(k,j,i)    )                  &
    5330                    -      (  5.0_wp * ibit17 * adv_mom_5                        &
    5331                        +              ibit16 * adv_mom_3                        &
    5332                           ) *                                                &
    5333                              ( u(k_pp,j,i)  - u(k-1,j,i)  )                  &
    5334                    +      (           ibit17 * adv_mom_5                        &
    5335                            ) *                                               &
    5336                              ( u(k_ppp,j,i) - u(k_mm,j,i) )                  &
    5337                                               )
    5338 !
    5339 !--             Calculate the divergence of the velocity field. A respective
    5340 !--             correction is needed to overcome numerical instabilities caused
    5341 !--             by a not sufficient reduction of divergences near topography.
    5342                 div = ( ( u_comp    * ( ibit9 + ibit10 + ibit11 )             &
    5343                 - ( u(k,j,i)   + u(k,j,i-1)   )                               &
    5344                                     * ( IBITS(wall_flags_0(k,j,i-1),9,1)      &
    5345                                       + IBITS(wall_flags_0(k,j,i-1),10,1)     &
    5346                                       + IBITS(wall_flags_0(k,j,i-1),11,1)     &
    5347                                       )                                       &
    5348                   ) * rho_air(k) * ddx                                        &
    5349                +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
    5350                   - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
    5351                                     * ( IBITS(wall_flags_0(k,j-1,i),12,1)     &
    5352                                       + IBITS(wall_flags_0(k,j-1,i),13,1)     &
    5353                                       + IBITS(wall_flags_0(k,j-1,i),14,1)     &
    5354                                       )                                       &
    5355                   ) * rho_air(k) * ddy                                        &
    5356                +  ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )     &
    5357                 - ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)             &
    5358                                     * ( IBITS(wall_flags_0(k-1,j,i),15,1)     &
    5359                                       + IBITS(wall_flags_0(k-1,j,i),16,1)     &
    5360                                       + IBITS(wall_flags_0(k-1,j,i),17,1)     &
    5361                                       )                                       &
    5362                   ) * ddzw(k)   &
    5363                 ) * 0.5_wp
    5364 
    5365 
    5366                 tend(k,j,i) = - (                                              &
    5367                                ( flux_r + diss_r - flux_l - diss_l ) * ddx     &
    5368                              + ( flux_n + diss_n - flux_s - diss_s ) * ddy     &
    5369                              + ( ( flux_t + diss_t ) -                         &
    5370                                  ( flux_d + diss_d )                           &
    5371                                                      ) * drho_air(k) * ddzw(k) &
    5372                                 ) + div * u(k,j,i)
    5373 
    5374 !++
    5375 !--             Statistical Evaluation of u'u'. The factor has to be applied
    5376 !--             for right evaluation when gallilei_trans = .T. .
    5377 !                sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                     &
    5378 !                              + ( flux_r    *                                 &
    5379 !                                ( u_comp    - 2.0_wp * hom(k,1,1,0) )            &
    5380 !                              / ( u_comp    - gu + 1.0E-20_wp   )             &
    5381 !                              +   diss_r    *                                 &
    5382 !                                  ABS( u_comp    - 2.0_wp * hom(k,1,1,0) )       &
    5383 !                              / ( ABS( u_comp    - gu ) + 1.0E-20_wp ) )      &
    5384 !                              *   weight_substep(intermediate_timestep_count)
    5385 !
    5386 !--             Statistical Evaluation of w'u'.
    5387 !                sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                   &
    5388 !                              + ( flux_t    + diss_t    )                     &
    5389 !                              *   weight_substep(intermediate_timestep_count)
    5390              ENDDO
    5391           ENDDO
    5392        ENDDO
    5393        !$acc end kernels
    5394 
    5395 !++
    5396 !       sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)
    5397 
    5398     END SUBROUTINE advec_u_ws_acc
    5399 
    5400 
    5401 !------------------------------------------------------------------------------!
    5402 ! Description:
    5403 ! ------------
    54044518!> Advection of v - Call for all grid points
    54054519!------------------------------------------------------------------------------!
     
    59065020   
    59075021   
    5908 !------------------------------------------------------------------------------!
    5909 ! Description:
    5910 ! ------------
    5911 !> Advection of v - Call for all grid points - accelerator version
    5912 !------------------------------------------------------------------------------!
    5913     SUBROUTINE advec_v_ws_acc
    5914 
    5915        USE arrays_3d,                                                          &
    5916            ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
    5917 
    5918        USE constants,                                                          &
    5919            ONLY:  adv_mom_1, adv_mom_3, adv_mom_5
    5920 
    5921        USE control_parameters,                                                 &
    5922            ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
    5923 
    5924        USE grid_variables,                                                     &
    5925            ONLY:  ddx, ddy
    5926 
    5927        USE indices,                                                            &
    5928            ONLY:  i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb,  &
    5929                   nzb_max, nzt, wall_flags_0
    5930            
    5931        USE kinds
    5932        
    5933 !        USE statistics,                                                       &
    5934 !            ONLY:  hom, sums_vs2_ws_l, sums_wsvs_ws_l, weight_substep
    5935 
    5936        IMPLICIT NONE
    5937 
    5938 
    5939        INTEGER(iwp) ::  i      !<
    5940        INTEGER(iwp) ::  ibit18 !<
    5941        INTEGER(iwp) ::  ibit19 !<
    5942        INTEGER(iwp) ::  ibit20 !<
    5943        INTEGER(iwp) ::  ibit21 !<
    5944        INTEGER(iwp) ::  ibit22 !<
    5945        INTEGER(iwp) ::  ibit23 !<
    5946        INTEGER(iwp) ::  ibit24 !<
    5947        INTEGER(iwp) ::  ibit25 !<
    5948        INTEGER(iwp) ::  ibit26 !<
    5949        INTEGER(iwp) ::  j      !<
    5950        INTEGER(iwp) ::  k      !<
    5951        INTEGER(iwp) ::  k_mm   !<
    5952        INTEGER(iwp) ::  k_mmm  !<
    5953        INTEGER(iwp) ::  k_pp   !<
    5954        INTEGER(iwp) ::  k_ppp  !<
    5955        INTEGER(iwp) ::  tn = 0 !<
    5956 
    5957        REAL(wp)    ::  diss_d   !<
    5958        REAL(wp)    ::  diss_l   !<
    5959        REAL(wp)    ::  diss_n   !<
    5960        REAL(wp)    ::  diss_r   !<
    5961        REAL(wp)    ::  diss_s   !<
    5962        REAL(wp)    ::  diss_t   !<
    5963        REAL(wp)    ::  div      !<
    5964        REAL(wp)    ::  flux_d   !<
    5965        REAL(wp)    ::  flux_l   !<
    5966        REAL(wp)    ::  flux_n   !<
    5967        REAL(wp)    ::  flux_r   !<
    5968        REAL(wp)    ::  flux_s   !<
    5969        REAL(wp)    ::  flux_t   !<
    5970        REAL(wp)    ::  gu       !<
    5971        REAL(wp)    ::  gv       !<
    5972        REAL(wp)    ::  u_comp   !<
    5973        REAL(wp)    ::  u_comp_l !<
    5974        REAL(wp)    ::  v_comp   !<
    5975        REAL(wp)    ::  v_comp_s !<
    5976        REAL(wp)    ::  w_comp   !<
    5977 
    5978        gu = 2.0_wp * u_gtrans
    5979        gv = 2.0_wp * v_gtrans
    5980 
    5981 !
    5982 !--    Computation of fluxes and tendency terms
    5983        !$acc kernels present( ddzw, tend, u, v, w, wall_flags_0 )
    5984        DO  i = i_left, i_right
    5985           DO  j = j_south, j_north
    5986              DO  k = nzb+1, nzt
    5987 
    5988                 ibit20 = IBITS(wall_flags_0(k,j,i-1),20,1)
    5989                 ibit19 = IBITS(wall_flags_0(k,j,i-1),19,1)
    5990                 ibit18 = IBITS(wall_flags_0(k,j,i-1),18,1)
    5991 
    5992                 u_comp_l                 = u(k,j-1,i) + u(k,j,i) - gu
    5993                 flux_l                   = u_comp_l * (                          &
    5994                                       ( 37.0_wp * ibit20 * adv_mom_5              &
    5995                                    +     7.0_wp * ibit19 * adv_mom_3              &
    5996                                    +              ibit18 * adv_mom_1              &
    5997                                       ) *                                      &
    5998                                      ( v(k,j,i)   + v(k,j,i-1) )               &
    5999                                -      (  8.0_wp * ibit20 * adv_mom_5              &
    6000                                    +              ibit19 * adv_mom_3              &
    6001                                       ) *                                      &
    6002                                      ( v(k,j,i+1) + v(k,j,i-2) )               &
    6003                                +      (           ibit20 * adv_mom_5              &
    6004                                       ) *                                      &
    6005                                      ( v(k,j,i+2) + v(k,j,i-3) )               &
    6006                                                  )
    6007 
    6008                 diss_l                   = - ABS( u_comp_l ) * (                 &
    6009                                       ( 10.0_wp * ibit20 * adv_mom_5              &
    6010                                    +     3.0_wp * ibit19 * adv_mom_3              &
    6011                                    +              ibit18 * adv_mom_1              &
    6012                                       ) *                                      &
    6013                                      ( v(k,j,i)   - v(k,j,i-1) )               &
    6014                                -      (  5.0_wp * ibit20 * adv_mom_5              &
    6015                                    +              ibit19 * adv_mom_3              &
    6016                                       ) *                                      &
    6017                                      ( v(k,j,i+1) - v(k,j,i-2) )               &
    6018                                +      (           ibit20 * adv_mom_5              &
    6019                                       ) *                                      &
    6020                                      ( v(k,j,i+2) - v(k,j,i-3) )               &
    6021                                                            )
    6022 
    6023                 ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
    6024                 ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
    6025                 ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
    6026 
    6027                 u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
    6028                 flux_r    = u_comp * (                                       &
    6029                           ( 37.0_wp * ibit20 * adv_mom_5                        &
    6030                        +     7.0_wp * ibit19 * adv_mom_3                        &
    6031                        +              ibit18 * adv_mom_1                        &
    6032                           ) *                                                &
    6033                                  ( v(k,j,i+1) + v(k,j,i)   )                 &
    6034                    -      (  8.0_wp * ibit20 * adv_mom_5                        &
    6035                        +              ibit19 * adv_mom_3                        &
    6036                           ) *                                                &
    6037                                  ( v(k,j,i+2) + v(k,j,i-1) )                 &
    6038                    +      (           ibit20 * adv_mom_5                        &
    6039                           ) *                                                &
    6040                                  ( v(k,j,i+3) + v(k,j,i-2) )                 &
    6041                                      )
    6042 
    6043                 diss_r    = - ABS( u_comp ) * (                              &
    6044                           ( 10.0_wp * ibit20 * adv_mom_5                        &
    6045                        +     3.0_wp * ibit19 * adv_mom_3                        &
    6046                        +              ibit18 * adv_mom_1                        &
    6047                           ) *                                                &
    6048                                  ( v(k,j,i+1) - v(k,j,i)  )                  &
    6049                    -      (  5.0_wp * ibit20 * adv_mom_5                        &
    6050                        +              ibit19 * adv_mom_3                        &
    6051                           ) *                                                &
    6052                                  ( v(k,j,i+2) - v(k,j,i-1) )                 &
    6053                    +      (           ibit20 * adv_mom_5                        &
    6054                           ) *                                                &
    6055                                  ( v(k,j,i+3) - v(k,j,i-2) )                 &
    6056                                               )
    6057 
    6058                 ibit23 = IBITS(wall_flags_0(k,j-1,i),23,1)
    6059                 ibit22 = IBITS(wall_flags_0(k,j-1,i),22,1)
    6060                 ibit21 = IBITS(wall_flags_0(k,j-1,i),21,1)
    6061 
    6062 
    6063                 v_comp_s              = v(k,j,i) + v(k,j-1,i) - gv
    6064                 flux_s                = v_comp_s    * (                       &
    6065                                    ( 37.0_wp * ibit23 * adv_mom_5                &
    6066                                 +     7.0_wp * ibit22 * adv_mom_3                &
    6067                                 +              ibit21 * adv_mom_1                &
    6068                                    ) *                                        &
    6069                                      ( v(k,j,i)   + v(k,j-1,i) )              &
    6070                             -      (  8.0_wp * ibit23 * adv_mom_5                &
    6071                                 +              ibit22 * adv_mom_3                &
    6072                                    ) *                                        &
    6073                                      ( v(k,j+1,i) + v(k,j-2,i) )              &
    6074                             +      (           ibit23 * adv_mom_5                &
    6075                                    ) *                                        &
    6076                                      ( v(k,j+2,i) + v(k,j-3,i) )              &
    6077                                                  )
    6078 
    6079                 diss_s                = - ABS( v_comp_s ) * (                 &
    6080                                    ( 10.0_wp * ibit23 * adv_mom_5                &
    6081                                 +     3.0_wp * ibit22 * adv_mom_3                &
    6082                                 +              ibit21 * adv_mom_1                &
    6083                                    ) *                                        &
    6084                                      ( v(k,j,i)   - v(k,j-1,i) )              &
    6085                             -      (  5.0_wp * ibit23 * adv_mom_5                &
    6086                                 +              ibit22 * adv_mom_3                &
    6087                                    ) *                                        &
    6088                                      ( v(k,j+1,i) - v(k,j-2,i) )              &
    6089                             +      (           ibit23 * adv_mom_5                &
    6090                                    ) *                                        &
    6091                                      ( v(k,j+2,i) - v(k,j-3,i) )              &
    6092                                                           )
    6093 
    6094                 ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
    6095                 ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
    6096                 ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
    6097 
    6098                 v_comp = v(k,j+1,i) + v(k,j,i)
    6099                 flux_n = ( v_comp - gv ) * (                                 &
    6100                           ( 37.0_wp * ibit23 * adv_mom_5                        &
    6101                        +     7.0_wp * ibit22 * adv_mom_3                        &
    6102                        +              ibit21 * adv_mom_1                        &
    6103                           ) *                                                &
    6104                                  ( v(k,j+1,i) + v(k,j,i)   )                 &
    6105                    -      (  8.0_wp * ibit23 * adv_mom_5                        &
    6106                        +              ibit22 * adv_mom_3                        &
    6107                           ) *                                                &
    6108                                  ( v(k,j+2,i) + v(k,j-1,i) )                 &
    6109                    +      (           ibit23 * adv_mom_5                        &
    6110                           ) *                                                &
    6111                                  ( v(k,j+3,i) + v(k,j-2,i) )                 &
    6112                                      )
    6113 
    6114                 diss_n = - ABS( v_comp - gv ) * (                         &
    6115                           ( 10.0_wp * ibit23 * adv_mom_5                        &
    6116                        +     3.0_wp * ibit22 * adv_mom_3                        &
    6117                        +              ibit21 * adv_mom_1                        &
    6118                           ) *                                                &
    6119                                  ( v(k,j+1,i) - v(k,j,i)  )                  &
    6120                    -      (  5.0_wp * ibit23 * adv_mom_5                        &
    6121                        +              ibit22 * adv_mom_3                        &
    6122                           ) *                                                &
    6123                                  ( v(k,j+2,i) - v(k,j-1,i) )                 &
    6124                    +      (           ibit23 * adv_mom_5                        &
    6125                           ) *                                                &
    6126                                  ( v(k,j+3,i) - v(k,j-2,i) )                 &
    6127                                                      )
    6128 
    6129                 ibit26 = IBITS(wall_flags_0(k-1,j,i),26,1)
    6130                 ibit25 = IBITS(wall_flags_0(k-1,j,i),25,1)
    6131                 ibit24 = IBITS(wall_flags_0(k-1,j,i),24,1)
    6132 
    6133                 k_pp  = k + 2 * ibit26
    6134                 k_mm  = k - 2 * ( ibit25 + ibit26 )
    6135                 k_mmm = k - 3 * ibit26
    6136 
    6137                 w_comp    = w(k-1,j-1,i) + w(k-1,j,i)
    6138                 flux_d    = w_comp  * (                                      &
    6139                           ( 37.0_wp * ibit26 * adv_mom_5                        &
    6140                        +     7.0_wp * ibit25 * adv_mom_3                        &
    6141                        +              ibit24 * adv_mom_1                        &
    6142                           ) *                                                &
    6143                              ( v(k,j,i)     + v(k-1,j,i)  )                  &
    6144                    -      (  8.0_wp * ibit26 * adv_mom_5                        &
    6145                        +              ibit25 * adv_mom_3                        &
    6146                           ) *                                                &
    6147                              ( v(k+1,j,i)  + v(k_mm,j,i)  )                  &
    6148                    +      (           ibit26 * adv_mom_5                        &
    6149                           ) *                                                &
    6150                              ( v(k_pp,j,i) + v(k_mmm,j,i) )                  &
    6151                                       )
    6152 
    6153                 diss_d    = - ABS( w_comp ) * (                              &
    6154                           ( 10.0_wp * ibit26 * adv_mom_5                        &
    6155                        +     3.0_wp * ibit25 * adv_mom_3                        &
    6156                        +              ibit24 * adv_mom_1                        &
    6157                           ) *                                                &
    6158                              ( v(k,j,i)     - v(k-1,j,i)  )                  &
    6159                    -      (  5.0_wp * ibit26 * adv_mom_5                        &
    6160                        +              ibit25 * adv_mom_3                        &
    6161                           ) *                                                &
    6162                              ( v(k+1,j,i)  - v(k_mm,j,i)  )                  &
    6163                    +      (           ibit26 * adv_mom_5                        &
    6164                           ) *                                                &
    6165                              ( v(k_pp,j,i) - v(k_mmm,j,i) )                  &
    6166                                                )
    6167 !
    6168 !--             k index has to be modified near bottom and top, else array
    6169 !--             subscripts will be exceeded.
    6170                 ibit26 = IBITS(wall_flags_0(k,j,i),26,1)
    6171                 ibit25 = IBITS(wall_flags_0(k,j,i),25,1)
    6172                 ibit24 = IBITS(wall_flags_0(k,j,i),24,1)
    6173 
    6174                 k_ppp = k + 3 * ibit26
    6175                 k_pp  = k + 2 * ( 1 - ibit24  )
    6176                 k_mm  = k - 2 * ibit26
    6177 
    6178                 w_comp    = w(k,j-1,i) + w(k,j,i)
    6179                 flux_t    = w_comp * rho_air_zw(k) * (                       &
    6180                           ( 37.0_wp * ibit26 * adv_mom_5                        &
    6181                        +     7.0_wp * ibit25 * adv_mom_3                        &
    6182                        +              ibit24 * adv_mom_1                        &
    6183                           ) *                                                &
    6184                              ( v(k+1,j,i)   + v(k,j,i)    )                  &
    6185                    -      (  8.0_wp * ibit26 * adv_mom_5                        &
    6186                        +              ibit25 * adv_mom_3                        &
    6187                           ) *                                                &
    6188                              ( v(k_pp,j,i)  + v(k-1,j,i)  )                  &
    6189                    +      (           ibit26 * adv_mom_5                        &
    6190                           ) *                                                &
    6191                              ( v(k_ppp,j,i) + v(k_mm,j,i) )                  &
    6192                                       )
    6193 
    6194                 diss_t    = - ABS( w_comp ) * rho_air_zw(k) * (              &
    6195                           ( 10.0_wp * ibit26 * adv_mom_5                        &
    6196                        +     3.0_wp * ibit25 * adv_mom_3                        &
    6197                        +              ibit24 * adv_mom_1                        &
    6198                           ) *                                                &
    6199                              ( v(k+1,j,i)   - v(k,j,i)    )                  &
    6200                    -      (  5.0_wp * ibit26 * adv_mom_5                        &
    6201                        +              ibit25 * adv_mom_3                        &
    6202                           ) *                                                &
    6203                              ( v(k_pp,j,i)  - v(k-1,j,i)  )                  &
    6204                    +      (           ibit26 * adv_mom_5                        &
    6205                           ) *                                                &
    6206                              ( v(k_ppp,j,i) - v(k_mm,j,i) )                  &
    6207                                                )
    6208 !
    6209 !--             Calculate the divergence of the velocity field. A respective
    6210 !--             correction is needed to overcome numerical instabilities caused
    6211 !--             by a not sufficient reduction of divergences near topography.
    6212                 div = ( ( ( u_comp     + gu )                                 &
    6213                                        * ( ibit18 + ibit19 + ibit20 )         &
    6214                 - ( u(k,j-1,i)   + u(k,j,i) )                                 &
    6215                                        * ( IBITS(wall_flags_0(k,j,i-1),18,1)  &
    6216                                          + IBITS(wall_flags_0(k,j,i-1),19,1)  &
    6217                                          + IBITS(wall_flags_0(k,j,i-1),20,1)  &
    6218                                          )                                    &
    6219                   ) * rho_air(k) * ddx                                        &
    6220                +  ( v_comp                                                    &
    6221                                        * ( ibit21 + ibit22 + ibit23 )         &
    6222                 - ( v(k,j,i)     + v(k,j-1,i) )                               &
    6223                                        * ( IBITS(wall_flags_0(k,j-1,i),21,1)  &
    6224                                          + IBITS(wall_flags_0(k,j-1,i),22,1)  &
    6225                                          + IBITS(wall_flags_0(k,j-1,i),23,1)  &
    6226                                          )                                    &
    6227                   ) * rho_air(k) * ddy                                        &
    6228                +  ( w_comp * rho_air_zw(k)                                    &
    6229                                        * ( ibit24 + ibit25 + ibit26 )         &
    6230                 - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)             &
    6231                                        * ( IBITS(wall_flags_0(k-1,j,i),24,1)  &
    6232                                          + IBITS(wall_flags_0(k-1,j,i),25,1)  &
    6233                                          + IBITS(wall_flags_0(k-1,j,i),26,1)  &
    6234                                          )                                    &
    6235                    ) * ddzw(k)   &
    6236                 ) * 0.5_wp
    6237 
    6238 
    6239                 tend(k,j,i) = - (                                              &
    6240                                ( flux_r + diss_r - flux_l - diss_l ) * ddx     &
    6241                              + ( flux_n + diss_n - flux_s - diss_s ) * ddy     &
    6242                              + ( ( flux_t + diss_t ) -                         &
    6243                                  ( flux_d + diss_d )                           &
    6244                                ) * drho_air(k) * ddzw(k)                       &
    6245                                 ) + div * v(k,j,i)
    6246 
    6247 
    6248 !++
    6249 !--             Statistical Evaluation of v'v'. The factor has to be applied
    6250 !--             for right evaluation when gallilei_trans = .T. .
    6251 !                sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                  &
    6252 !                      + ( flux_n                                           &
    6253 !                      * ( v_comp - 2.0_wp * hom(k,1,2,0) )                    &
    6254 !                      / ( v_comp - gv + 1.0E-20_wp )                       &
    6255 !                      +   diss_n                                           &
    6256 !                      *   ABS( v_comp - 2.0_wp * hom(k,1,2,0) )               &
    6257 !                      / ( ABS( v_comp - gv ) +1.0E-20_wp ) )               &
    6258 !                      *   weight_substep(intermediate_timestep_count)
    6259 !
    6260 !--              Statistical Evaluation of w'v'.
    6261 !                 sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                &
    6262 !                              + ( flux_t + diss_t )                         &
    6263 !                              *   weight_substep(intermediate_timestep_count)
    6264 
    6265              ENDDO
    6266           ENDDO
    6267        ENDDO
    6268        !$acc end kernels
    6269 
    6270 !++
    6271 !       sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
    6272 
    6273     END SUBROUTINE advec_v_ws_acc
    62745022   
    62755023   
     
    67565504
    67575505
    6758 !------------------------------------------------------------------------------!
    6759 ! Description:
    6760 ! ------------
    6761 !> Advection of w - Call for all grid points - accelerator version
    6762 !------------------------------------------------------------------------------!
    6763     SUBROUTINE advec_w_ws_acc
    6764 
    6765        USE arrays_3d,                                                          &
    6766            ONLY:  ddzu, drho_air_zw, tend, u, v, w, rho_air, rho_air_zw
    6767 
    6768        USE constants,                                                          &
    6769            ONLY:  adv_mom_1, adv_mom_3, adv_mom_5
    6770 
    6771        USE control_parameters,                                                 &
    6772            ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
    6773 
    6774        USE grid_variables,                                                     &
    6775            ONLY:  ddx, ddy
    6776 
    6777        USE indices,                                                            &
    6778            ONLY:  i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb,  &
    6779                   nzb_max, nzt, wall_flags_0, wall_flags_00
    6780            
    6781        USE kinds
    6782        
    6783 !        USE statistics,                                                       &
    6784 !            ONLY:  hom, sums_ws2_ws_l, weight_substep
    6785 
    6786        IMPLICIT NONE
    6787 
    6788        INTEGER(iwp) ::  i      !<
    6789        INTEGER(iwp) ::  ibit27 !<
    6790        INTEGER(iwp) ::  ibit28 !<
    6791        INTEGER(iwp) ::  ibit29 !<
    6792        INTEGER(iwp) ::  ibit30 !<
    6793        INTEGER(iwp) ::  ibit31 !<
    6794        INTEGER(iwp) ::  ibit32 !<
    6795        INTEGER(iwp) ::  ibit33 !<
    6796        INTEGER(iwp) ::  ibit34 !<
    6797        INTEGER(iwp) ::  ibit35 !<
    6798        INTEGER(iwp) ::  j      !<
    6799        INTEGER(iwp) ::  k      !<
    6800        INTEGER(iwp) ::  k_mmm  !<
    6801        INTEGER(iwp) ::  k_mm   !<
    6802        INTEGER(iwp) ::  k_pp   !<
    6803        INTEGER(iwp) ::  k_ppp  !<
    6804        INTEGER(iwp) ::  tn = 0 !<
    6805 
    6806        REAL(wp)    ::  diss_d   !<
    6807        REAL(wp)    ::  diss_l   !<
    6808        REAL(wp)    ::  diss_n   !<
    6809        REAL(wp)    ::  diss_r   !<
    6810        REAL(wp)    ::  diss_s   !<
    6811        REAL(wp)    ::  diss_t   !<
    6812        REAL(wp)    ::  div      !<
    6813        REAL(wp)    ::  flux_d   !<
    6814        REAL(wp)    ::  flux_l   !<
    6815        REAL(wp)    ::  flux_n   !<
    6816        REAL(wp)    ::  flux_r   !<
    6817        REAL(wp)    ::  flux_s   !<
    6818        REAL(wp)    ::  flux_t   !<
    6819        REAL(wp)    ::  gu       !<
    6820        REAL(wp)    ::  gv       !<
    6821        REAL(wp)    ::  u_comp   !<
    6822        REAL(wp)    ::  u_comp_l !<
    6823        REAL(wp)    ::  v_comp   !<
    6824        REAL(wp)    ::  v_comp_s !<
    6825        REAL(wp)    ::  w_comp   !<
    6826 
    6827        gu = 2.0_wp * u_gtrans
    6828        gv = 2.0_wp * v_gtrans
    6829 
    6830 
    6831 !
    6832 !--    Computation of fluxes and tendency terms
    6833        !$acc kernels present( ddzu, tend, u, v, w, wall_flags_0, wall_flags_00 )
    6834        DO i = i_left, i_right
    6835           DO  j = j_south, j_north
    6836              DO  k = nzb+1, nzt
    6837 
    6838                 ibit27 = IBITS(wall_flags_0(k,j,i-1),27,1)
    6839                 ibit28 = IBITS(wall_flags_0(k,j,i-1),28,1)
    6840                 ibit29 = IBITS(wall_flags_0(k,j,i-1),29,1)
    6841 
    6842                 u_comp_l                 = u(k+1,j,i) + u(k,j,i) - gu
    6843                 flux_l                   = u_comp_l * (                        &
    6844                                       ( 37.0_wp * ibit29 * adv_mom_5              &
    6845                                    +     7.0_wp * ibit28 * adv_mom_3              &
    6846                                    +              ibit27 * adv_mom_1              &
    6847                                       ) *                                      &
    6848                                      ( w(k,j,i)   + w(k,j,i-1) )               &
    6849                                -      (  8.0_wp * ibit29 * adv_mom_5              &
    6850                                    +              ibit28 * adv_mom_3              &
    6851                                       ) *                                      &
    6852                                      ( w(k,j,i+1) + w(k,j,i-2) )               &
    6853                                +      (           ibit29 * adv_mom_5              &
    6854                                       ) *                                      &
    6855                                      ( w(k,j,i+2) + w(k,j,i-3) )               &
    6856                                                  )
    6857 
    6858                 diss_l                    = - ABS( u_comp_l ) * (              &
    6859                                         ( 10.0_wp * ibit29 * adv_mom_5            &
    6860                                      +     3.0_wp * ibit28 * adv_mom_3            &
    6861                                      +              ibit27 * adv_mom_1            &
    6862                                         ) *                                    &
    6863                                      ( w(k,j,i)   - w(k,j,i-1) )               &
    6864                                  -      (  5.0_wp * ibit29 * adv_mom_5            &
    6865                                      +              ibit28 * adv_mom_3            &
    6866                                         ) *                                    &
    6867                                      ( w(k,j,i+1) - w(k,j,i-2) )               &
    6868                                  +      (           ibit29 * adv_mom_5            &
    6869                                         ) *                                    &
    6870                                      ( w(k,j,i+2) - w(k,j,i-3) )               &
    6871                                                             )
    6872 
    6873                 ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
    6874                 ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
    6875                 ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
    6876 
    6877                 u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    6878                 flux_r    = u_comp * (                                       &
    6879                           ( 37.0_wp * ibit29 * adv_mom_5                        &
    6880                        +     7.0_wp * ibit28 * adv_mom_3                        &
    6881                        +              ibit27 * adv_mom_1                        &
    6882                           ) *                                                &
    6883                                  ( w(k,j,i+1) + w(k,j,i)   )                 &
    6884                    -      (  8.0_wp * ibit29 * adv_mom_5                        &
    6885                        +              ibit28 * adv_mom_3                        &
    6886                           ) *                                                &
    6887                                  ( w(k,j,i+2) + w(k,j,i-1) )                 &
    6888                    +      (           ibit29 * adv_mom_5                        &
    6889                           ) *                                                &
    6890                                  ( w(k,j,i+3) + w(k,j,i-2) )                 &
    6891                                      )
    6892 
    6893                 diss_r    = - ABS( u_comp ) * (                              &
    6894                           ( 10.0_wp * ibit29 * adv_mom_5                        &
    6895                        +     3.0_wp * ibit28 * adv_mom_3                        &
    6896                        +              ibit27 * adv_mom_1                        &
    6897                           ) *                                                &
    6898                                  ( w(k,j,i+1) - w(k,j,i)  )                  &
    6899                    -      (  5.0_wp * ibit29 * adv_mom_5                        &
    6900                        +              ibit28 * adv_mom_3                        &
    6901                           ) *                                                &
    6902                                  ( w(k,j,i+2) - w(k,j,i-1) )                 &
    6903                    +      (           ibit29 * adv_mom_5                        &
    6904                           ) *                                                &
    6905                                  ( w(k,j,i+3) - w(k,j,i-2) )                 &
    6906                                               )
    6907                 ibit32 = IBITS(wall_flags_00(k,j-1,i),0,1)
    6908                 ibit31 = IBITS(wall_flags_0(k,j-1,i),31,1)
    6909                 ibit30 = IBITS(wall_flags_0(k,j-1,i),30,1)
    6910 
    6911                 v_comp_s               = v(k+1,j,i) + v(k,j,i) - gv
    6912                 flux_s                 = v_comp_s * (                         &
    6913                                     ( 37.0_wp * ibit32 * adv_mom_5               &
    6914                                  +     7.0_wp * ibit31 * adv_mom_3               &
    6915                                  +              ibit30 * adv_mom_1               &
    6916                                     ) *                                       &
    6917                                      ( w(k,j,i)   + w(k,j-1,i) )              &
    6918                              -      (  8.0_wp * ibit32 * adv_mom_5               &
    6919                                  +              ibit31 * adv_mom_3               &
    6920                                     ) *                                       &
    6921                                      ( w(k,j+1,i) + w(k,j-2,i) )              &
    6922                              +      (           ibit32 * adv_mom_5               &
    6923                                     ) *                                       &
    6924                                      ( w(k,j+2,i) + w(k,j-3,i) )              &
    6925                                                )
    6926 
    6927                 diss_s                 = - ABS( v_comp_s ) * (                &
    6928                                     ( 10.0_wp * ibit32 * adv_mom_5               &
    6929                                  +     3.0_wp * ibit31 * adv_mom_3               &
    6930                                  +              ibit30 * adv_mom_1               &
    6931                                     ) *                                       &
    6932                                      ( w(k,j,i)   - w(k,j-1,i) )              &
    6933                              -      (  5.0_wp * ibit32 * adv_mom_5               &
    6934                                  +              ibit31 * adv_mom_3               &
    6935                                     ) *                                       &
    6936                                      ( w(k,j+1,i) - w(k,j-2,i) )              &
    6937                              +      (           ibit32 * adv_mom_5               &
    6938                                     ) *                                       &
    6939                                      ( w(k,j+2,i) - w(k,j-3,i) )              &
    6940                                                         )
    6941 
    6942                 ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
    6943                 ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
    6944                 ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
    6945 
    6946                 v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    6947                 flux_n    = v_comp * (                                       &
    6948                           ( 37.0_wp * ibit32 * adv_mom_5                        &
    6949                        +     7.0_wp * ibit31 * adv_mom_3                        &
    6950                        +              ibit30 * adv_mom_1                        &
    6951                           ) *                                                 &
    6952                                  ( w(k,j+1,i) + w(k,j,i)   )                 &
    6953                    -      (  8.0_wp * ibit32 * adv_mom_5                        &
    6954                        +              ibit31 * adv_mom_3                        &
    6955                           ) *                                                &
    6956                                  ( w(k,j+2,i) + w(k,j-1,i) )                 &
    6957                    +      (           ibit32 * adv_mom_5                        &
    6958                           ) *                                                &
    6959                                  ( w(k,j+3,i) + w(k,j-2,i) )                 &
    6960                                      )
    6961 
    6962                 diss_n    = - ABS( v_comp ) * (                              &
    6963                           ( 10.0_wp * ibit32 * adv_mom_5                        &
    6964                        +     3.0_wp * ibit31 * adv_mom_3                        &
    6965                        +              ibit30 * adv_mom_1                        &
    6966                           ) *                                                &
    6967                                  ( w(k,j+1,i) - w(k,j,i)  )                  &
    6968                    -      (  5.0_wp * ibit32 * adv_mom_5                        &
    6969                        +              ibit31 * adv_mom_3                        &
    6970                           ) *                                                &
    6971                                  ( w(k,j+2,i) - w(k,j-1,i) )                 &
    6972                    +      (           ibit32 * adv_mom_5                        &
    6973                           ) *                                                &
    6974                                  ( w(k,j+3,i) - w(k,j-2,i) )                 &
    6975                                               )
    6976 
    6977                 ibit35 = IBITS(wall_flags_00(k-1,j,i),3,1)
    6978                 ibit34 = IBITS(wall_flags_00(k-1,j,i),2,1)
    6979                 ibit33 = IBITS(wall_flags_00(k-1,j,i),1,1)
    6980 
    6981                 k_pp  = k + 2 * ibit35
    6982                 k_mm  = k - 2 * ( ibit34 + ibit35 )
    6983                 k_mmm = k - 3 * ibit35
    6984 
    6985                 w_comp    = w(k,j,i) + w(k-1,j,i)
    6986                 flux_d    = w_comp  * (                                      &
    6987                           ( 37.0_wp * ibit35 * adv_mom_5                        &
    6988                        +     7.0_wp * ibit34 * adv_mom_3                        &
    6989                        +              ibit33 * adv_mom_1                        &
    6990                           ) *                                                &
    6991                              ( w(k,j,i)    + w(k-1,j,i)   )                  &
    6992                    -      (  8.0_wp * ibit35 * adv_mom_5                        &
    6993                        +              ibit34 * adv_mom_3                        &
    6994                           ) *                                                &
    6995                              ( w(k+1,j,i)  + w(k_mm,j,i)  )                  &
    6996                    +      (           ibit35 * adv_mom_5                        &
    6997                           ) *                                                &
    6998                              ( w(k_pp,j,i) + w(k_mmm,j,i) )                  &
    6999                                        )
    7000 
    7001                 diss_d    = - ABS( w_comp ) * (                              &
    7002                           ( 10.0_wp * ibit35 * adv_mom_5                        &
    7003                        +     3.0_wp * ibit34 * adv_mom_3                        &
    7004                        +              ibit33 * adv_mom_1                        &
    7005                           ) *                                                &
    7006                              ( w(k,j,i)    - w(k-1,j,i)   )                  &
    7007                    -      (  5.0_wp * ibit35 * adv_mom_5                        &
    7008                        +              ibit34 * adv_mom_3                        &
    7009                           ) *                                                &
    7010                              ( w(k+1,j,i)  - w(k_mm,j,i)  )                  &
    7011                    +      (           ibit35 * adv_mom_5                        &
    7012                           ) *                                                &
    7013                              ( w(k_pp,j,i) - w(k_mmm,j,i) )                  &
    7014                                                )
    7015 
    7016 !
    7017 !--             k index has to be modified near bottom and top, else array
    7018 !--             subscripts will be exceeded.
    7019                 ibit35 = IBITS(wall_flags_00(k,j,i),3,1)
    7020                 ibit34 = IBITS(wall_flags_00(k,j,i),2,1)
    7021                 ibit33 = IBITS(wall_flags_00(k,j,i),1,1)
    7022 
    7023                 k_ppp = k + 3 * ibit35
    7024                 k_pp  = k + 2 * ( 1 - ibit33  )
    7025                 k_mm  = k - 2 * ibit35
    7026 
    7027                 w_comp    = w(k+1,j,i) + w(k,j,i)
    7028                 flux_t    = w_comp * rho_air(k+1) * (                        &
    7029                           ( 37.0_wp * ibit35 * adv_mom_5                        &
    7030                        +     7.0_wp * ibit34 * adv_mom_3                        &
    7031                        +              ibit33 * adv_mom_1                        &
    7032                           ) *                                                &
    7033                              ( w(k+1,j,i)  + w(k,j,i)     )                  &
    7034                    -      (  8.0_wp * ibit35 * adv_mom_5                        &
    7035                        +              ibit34 * adv_mom_3                        &
    7036                           ) *                                                &
    7037                              ( w(k_pp,j,i)  + w(k-1,j,i)  )                  &
    7038                    +      (           ibit35 * adv_mom_5                        &
    7039                           ) *                                                &
    7040                              ( w(k_ppp,j,i) + w(k_mm,j,i) )                  &
    7041                                        )
    7042 
    7043                 diss_t    = - ABS( w_comp ) * rho_air(k+1) * (               &
    7044                           ( 10.0_wp * ibit35 * adv_mom_5                        &
    7045                        +     3.0_wp * ibit34 * adv_mom_3                        &
    7046                        +              ibit33 * adv_mom_1                        &
    7047                           ) *                                                &
    7048                              ( w(k+1,j,i)   - w(k,j,i)    )                  &
    7049                    -      (  5.0_wp * ibit35 * adv_mom_5                        &
    7050                        +              ibit34 * adv_mom_3                        &
    7051                           ) *                                                &
    7052                              ( w(k_pp,j,i)  - w(k-1,j,i)  )                  &
    7053                    +      (           ibit35 * adv_mom_5                        &
    7054                           ) *                                                &
    7055                              ( w(k_ppp,j,i) - w(k_mm,j,i) )                  &
    7056                                                )
    7057 !
    7058 !--             Calculate the divergence of the velocity field. A respective
    7059 !--             correction is needed to overcome numerical instabilities caused
    7060 !--             by a not sufficient reduction of divergences near topography.
    7061                 div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 )      &
    7062                   - ( u(k+1,j,i) + u(k,j,i)   )                               &
    7063                                     * ( IBITS(wall_flags_0(k,j,i-1),27,1)     &
    7064                                       + IBITS(wall_flags_0(k,j,i-1),28,1)     &
    7065                                       + IBITS(wall_flags_0(k,j,i-1),29,1)     &
    7066                                       )                                       &
    7067                   ) * rho_air_zw(k) * ddx                                     &
    7068               +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            &
    7069                   - ( v(k+1,j,i) + v(k,j,i)   )                               &
    7070                                     * ( IBITS(wall_flags_0(k,j-1,i),30,1)     &
    7071                                       + IBITS(wall_flags_0(k,j-1,i),31,1)     &
    7072                                       + IBITS(wall_flags_00(k,j-1,i),0,1)     &
    7073                                       )                                       &
    7074                   ) * rho_air_zw(k) * ddy                                     &
    7075               +   ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 )      &
    7076                 - ( w(k,j,i)   + w(k-1,j,i)   ) * rho_air(k)                  &
    7077                                     * ( IBITS(wall_flags_00(k-1,j,i),1,1)     &
    7078                                       + IBITS(wall_flags_00(k-1,j,i),2,1)     &
    7079                                       + IBITS(wall_flags_00(k-1,j,i),3,1)     &
    7080                                       )                                       &
    7081                   ) * ddzu(k+1)   &
    7082                 ) * 0.5_wp
    7083 
    7084 
    7085                 tend(k,j,i) = - (                                                &
    7086                                ( flux_r + diss_r - flux_l - diss_l ) * ddx       &
    7087                              + ( flux_n + diss_n - flux_s - diss_s ) * ddy       &
    7088                              + ( ( flux_t + diss_t ) -                           &
    7089                                  ( flux_d + diss_d ) * rho_air(k)                &
    7090                                ) * drho_air_zw(k) * ddzu(k+1)                    &
    7091                                  ) + div * w(k,j,i)
    7092 
    7093 
    7094 !++
    7095 !--             Statistical Evaluation of w'w'.
    7096 !                sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                    &
    7097 !                               + ( flux_t + diss_t )                    &
    7098 !                               *   weight_substep(intermediate_timestep_count)
    7099 
    7100              ENDDO
    7101           ENDDO
    7102        ENDDO
    7103        !$acc end kernels
    7104 
    7105     END SUBROUTINE advec_w_ws_acc
    7106 
    71075506 END MODULE advec_ws
Note: See TracChangeset for help on using the changeset viewer.