Changeset 2118


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

Location:
palm/trunk/SOURCE
Files:
1 deleted
32 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile

    r2051 r2118  
    2020# Current revisions:
    2121# ------------------
    22 #
     22# -cuda_fft_interfaces_mod
    2323#
    2424# Former revisions:
     
    312312        check_for_restart.f90 check_open.f90 check_parameters.f90 \
    313313        close_file.f90 compute_vpt.f90 coriolis.f90 cpulog_mod.f90 \
    314         cuda_fft_interfaces_mod.f90 data_log.f90 data_output_dvrp.f90 \
     314        data_log.f90 data_output_dvrp.f90 \
    315315        data_output_mask.f90 data_output_profiles.f90 \
    316316        data_output_ptseries.f90 data_output_spectra.f90 data_output_flight.f90\
     
    426426cpulog_mod.o: modules.o mod_kinds.o
    427427cpu_statistics.o: modules.o mod_kinds.o
    428 cuda_fft_interfaces_mod.o: cuda_fft_interfaces_mod.f90 modules.o mod_kinds.o
    429428data_log.o: modules.o mod_kinds.o
    430429data_output_dvrp.o: modules.o cpulog_mod.o mod_kinds.o
     
    456455exchange_horiz.o: modules.o cpulog_mod.o mod_kinds.o
    457456exchange_horiz_2d.o: modules.o cpulog_mod.o mod_kinds.o pmc_interface_mod.o
    458 fft_xy_mod.o: cuda_fft_interfaces_mod.o modules.o mod_kinds.o singleton_mod.o temperton_fft_mod.o
     457fft_xy_mod.o: modules.o mod_kinds.o singleton_mod.o temperton_fft_mod.o
    459458flow_statistics.o: modules.o cpulog_mod.o mod_kinds.o land_surface_model_mod.o \
    460459   netcdf_interface_mod.o radiation_model_mod.o
  • 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
  • palm/trunk/SOURCE/boundary_conds.f90

    r2101 r2118  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! OpenACC directives removed
    2323!
    2424! Former revisions:
     
    188188!-- Bottom boundary
    189189    IF ( ibc_uv_b == 1 )  THEN
    190        !$acc kernels present( u_p, v_p )
    191190       u_p(nzb,:,:) = u_p(nzb+1,:,:)
    192191       v_p(nzb,:,:) = v_p(nzb+1,:,:)
    193        !$acc end kernels
    194     ENDIF
    195 
    196     !$acc kernels present( nzb_w_inner, w_p )
     192    ENDIF
     193
    197194    DO  i = nxlg, nxrg
    198195       DO  j = nysg, nyng
     
    200197       ENDDO
    201198    ENDDO
    202     !$acc end kernels
    203199
    204200!
    205201!-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings.
    206202    IF ( ibc_uv_t == 0 )  THEN
    207        !$acc kernels present( u_init, u_p, v_init, v_p )
    208203        u_p(nzt+1,:,:) = u_init(nzt+1)
    209204        v_p(nzt+1,:,:) = v_init(nzt+1)
    210        !$acc end kernels
    211205    ELSEIF ( ibc_uv_t == 1 )  THEN
    212        !$acc kernels present( u_p, v_p )
    213206        u_p(nzt+1,:,:) = u_p(nzt,:,:)
    214207        v_p(nzt+1,:,:) = v_p(nzt,:,:)
    215        !$acc end kernels
    216208    ENDIF
    217209
    218210    IF ( .NOT. nest_domain )  THEN
    219        !$acc kernels present( w_p )
    220211       w_p(nzt:nzt+1,:,:) = 0.0_wp  ! nzt is not a prognostic level (but cf. pres)
    221        !$acc end kernels
    222212    ENDIF
    223213
     
    227217!-- the sea surface temperature of the coupled ocean model.
    228218    IF ( ibc_pt_b == 0 )  THEN
    229        !$acc kernels present( nzb_s_inner, pt, pt_p )
    230        !$acc loop independent
    231219       DO  i = nxlg, nxrg
    232           !$acc loop independent
    233220          DO  j = nysg, nyng
    234221             pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
    235222          ENDDO
    236223       ENDDO
    237        !$acc end kernels
    238224    ELSEIF ( ibc_pt_b == 1 )  THEN
    239        !$acc kernels present( nzb_s_inner, pt_p )
    240        !$acc loop independent
    241225       DO  i = nxlg, nxrg
    242           !$acc loop independent
    243226          DO  j = nysg, nyng
    244227             pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
    245228          ENDDO
    246229       ENDDO
    247       !$acc end kernels
    248230    ENDIF
    249231
     
    251233!-- Temperature at top boundary
    252234    IF ( ibc_pt_t == 0 )  THEN
    253        !$acc kernels present( pt, pt_p )
    254235        pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
    255236!
     
    259240           pt_p(nzt+1,:,:) = pt_init(nzt+1)
    260241        ENDIF
    261        !$acc end kernels
    262242    ELSEIF ( ibc_pt_t == 1 )  THEN
    263        !$acc kernels present( pt_p )
    264243        pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
    265        !$acc end kernels
    266244    ELSEIF ( ibc_pt_t == 2 )  THEN
    267        !$acc kernels present( dzu, pt_p )
    268245        pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
    269        !$acc end kernels
    270246    ENDIF
    271247
     
    274250!-- Generally Neumann conditions with de/dz=0 are assumed
    275251    IF ( .NOT. constant_diffusion )  THEN
    276        !$acc kernels present( e_p, nzb_s_inner )
    277        !$acc loop independent
    278252       DO  i = nxlg, nxrg
    279           !$acc loop independent
    280253          DO  j = nysg, nyng
    281254             e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
     
    285258          e_p(nzt+1,:,:) = e_p(nzt,:,:)
    286259       ENDIF
    287        !$acc end kernels
    288260    ENDIF
    289261
  • palm/trunk/SOURCE/buoyancy.f90

    r2101 r2118  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! OpenACC version of subroutine removed
    2323!
    2424! Former revisions:
     
    102102
    103103    PRIVATE
    104     PUBLIC buoyancy, buoyancy_acc
     104    PUBLIC buoyancy
    105105
    106106    INTERFACE buoyancy
     
    108108       MODULE PROCEDURE buoyancy_ij
    109109    END INTERFACE buoyancy
    110 
    111     INTERFACE buoyancy_acc
    112        MODULE PROCEDURE buoyancy_acc
    113     END INTERFACE buoyancy_acc
    114110
    115111 CONTAINS
     
    212208
    213209    END SUBROUTINE buoyancy
    214 
    215 
    216 !------------------------------------------------------------------------------!
    217 ! Description:
    218 ! ------------
    219 !> Call for all grid points - accelerator version
    220 !------------------------------------------------------------------------------!
    221     SUBROUTINE buoyancy_acc( var, wind_component )
    222 
    223        USE arrays_3d,                                                          &
    224            ONLY:  pt, pt_slope_ref, ref_state, tend
    225 
    226        USE control_parameters,                                                 &
    227            ONLY:  atmos_ocean_sign, cos_alpha_surface, g, message_string,      &
    228                   pt_surface, sin_alpha_surface, sloping_surface
    229 
    230        USE indices,                                                            &
    231            ONLY:  i_left, i_right, j_north, j_south, nxl, nxlg, nxlu, nxr,     &
    232                   nxrg, nyn, nyng, nys, nysg, nzb, nzb_s_inner, nzt
    233 
    234        USE kinds
    235 
    236        USE pegrid
    237 
    238 
    239        IMPLICIT NONE
    240 
    241        INTEGER(iwp) ::  i              !<
    242        INTEGER(iwp) ::  j              !<
    243        INTEGER(iwp) ::  k              !<
    244        INTEGER(iwp) ::  wind_component !<
    245        
    246 #if defined( __nopointer )
    247        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var !<
    248 #else
    249        REAL(wp), DIMENSION(:,:,:), POINTER ::  var
    250 #endif
    251 
    252 
    253        IF ( .NOT. sloping_surface )  THEN
    254 !
    255 !--       Normal case: horizontal surface
    256           !$acc kernels present( nzb_s_inner, ref_state, tend, var )
    257           !$acc loop
    258           DO  i = i_left, i_right
    259              DO  j = j_south, j_north
    260                 !$acc loop independent vector
    261                 DO  k = nzb_s_inner(j,i)+1, nzt-1
    262                    tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5_wp * &
    263                           (                                                    &
    264                              ( var(k,j,i)   - ref_state(k) )   / ref_state(k) + &
    265                              ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1) &
    266                           )
    267                 ENDDO
    268              ENDDO
    269           ENDDO
    270           !$acc end kernels
    271 
    272        ELSE
    273 !
    274 !--       Buoyancy term for a surface with a slope in x-direction. The equations
    275 !--       for both the u and w velocity-component contain proportionate terms.
    276 !--       Temperature field at time t=0 serves as environmental temperature.
    277 !--       Reference temperature (pt_surface) is the one at the lower left corner
    278 !--       of the total domain.
    279           IF ( wind_component == 1 )  THEN
    280 
    281              DO  i = nxlu, nxr
    282                 DO  j = nys, nyn
    283                    DO  k = nzb_s_inner(j,i)+1, nzt-1
    284                       tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *      &
    285                            0.5_wp * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
    286                                     - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
    287                                     ) / pt_surface
    288                    ENDDO
    289                 ENDDO
    290              ENDDO
    291 
    292           ELSEIF ( wind_component == 3 )  THEN
    293 
    294              DO  i = nxl, nxr
    295                 DO  j = nys, nyn
    296                    DO  k = nzb_s_inner(j,i)+1, nzt-1
    297                       tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *      &
    298                            0.5_wp * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
    299                                     - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
    300                                     ) / pt_surface
    301                    ENDDO
    302                 ENDDO
    303             ENDDO
    304 
    305           ELSE
    306 
    307              WRITE( message_string, * ) 'no term for component "',             &
    308                                        wind_component,'"'
    309              CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
    310 
    311           ENDIF
    312 
    313        ENDIF
    314 
    315     END SUBROUTINE buoyancy_acc
    316210
    317211
  • palm/trunk/SOURCE/check_parameters.f90

    r2101 r2118  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! OpenACC related parts of code removed
    2323!
    2424! Former revisions:
     
    518518    IF ( transpose_compute_overlap )  THEN
    519519       IF ( numprocs == 1 )  STOP '+++ transpose-compute-overlap not implemented for single PE runs'
    520 #if defined( __openacc )
    521        STOP '+++ transpose-compute-overlap not implemented for GPU usage'
    522 #endif
    523520    ENDIF
    524521
     
    774771    SELECT CASE ( TRIM( loop_optimization ) )
    775772
    776        CASE ( 'acc', 'cache', 'vector' )
     773       CASE ( 'cache', 'vector' )
    777774          CONTINUE
    778775
  • palm/trunk/SOURCE/coriolis.f90

    r2101 r2118  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! OpenACC version of subroutine removed
    2323!
    2424! Former revisions:
     
    3535! 1850 2016-04-08 13:29:27Z maronga
    3636! Module renamed
    37 !
    3837!
    3938! 1682 2015-10-07 23:56:08Z knoop
     
    7675
    7776    PRIVATE
    78     PUBLIC coriolis, coriolis_acc
     77    PUBLIC coriolis
    7978
    8079    INTERFACE coriolis
     
    8281       MODULE PROCEDURE coriolis_ij
    8382    END INTERFACE coriolis
    84 
    85     INTERFACE coriolis_acc
    86        MODULE PROCEDURE coriolis_acc
    87     END INTERFACE coriolis_acc
    8883
    8984 CONTAINS
     
    177172! Description:
    178173! ------------
    179 !> Call for all grid points - accelerator version
    180 !------------------------------------------------------------------------------!
    181     SUBROUTINE coriolis_acc( component )
    182 
    183        USE arrays_3d,                                                          &
    184            ONLY:  tend, u, ug, v, vg, w
    185            
    186        USE control_parameters,                                                 &
    187            ONLY:  f, fs, message_string
    188            
    189        USE indices,                                                            &
    190            ONLY:  i_left, i_right, j_north, j_south, nzb_u_inner,              &
    191                   nzb_v_inner, nzb_w_inner, nzt
    192                    
    193        USE kinds
    194 
    195        IMPLICIT NONE
    196 
    197        INTEGER(iwp) ::  component  !<
    198        INTEGER(iwp) ::  i          !<
    199        INTEGER(iwp) ::  j          !<
    200        INTEGER(iwp) ::  k          !<
    201 
    202 
    203 !
    204 !--    Compute Coriolis terms for the three velocity components
    205        SELECT CASE ( component )
    206 
    207 !
    208 !--       u-component
    209           CASE ( 1 )
    210              !$acc  kernels present( nzb_u_inner, tend, v, vg, w )
    211              DO  i = i_left, i_right
    212                 DO  j = j_south, j_north
    213                    DO  k = 1, nzt
    214                       IF  ( k > nzb_u_inner(j,i) )  THEN
    215                          tend(k,j,i) = tend(k,j,i) + f  *    ( 0.25_wp *       &
    216                                       ( v(k,j,i-1) + v(k,j,i) + v(k,j+1,i-1) + &
    217                                         v(k,j+1,i) ) - vg(k) )                 &
    218                                                    - fs *    ( 0.25_wp *       &
    219                                       ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) &
    220                                         + w(k,j,i)   )                         &
    221                                                              )
    222                       ENDIF
    223                    ENDDO
    224                 ENDDO
    225              ENDDO
    226              !$acc end kernels
    227 
    228 !
    229 !--       v-component
    230           CASE ( 2 )
    231              !$acc  kernels present( nzb_v_inner, tend, u, ug )
    232              DO  i = i_left, i_right
    233                 DO  j = j_south, j_north
    234                    DO  k = 1, nzt
    235                       IF  ( k > nzb_v_inner(j,i) )  THEN
    236                          tend(k,j,i) = tend(k,j,i) - f *     ( 0.25_wp *       &
    237                                       ( u(k,j-1,i) + u(k,j,i) + u(k,j-1,i+1) + &
    238                                         u(k,j,i+1) ) - ug(k) )
    239                       ENDIF
    240                    ENDDO
    241                 ENDDO
    242              ENDDO
    243              !$acc end kernels
    244 
    245 !
    246 !--       w-component
    247           CASE ( 3 )
    248              !$acc  kernels present( nzb_w_inner, tend, u )
    249              DO  i = i_left, i_right
    250                 DO  j = j_south, j_north
    251                    DO  k = 1, nzt
    252                       IF  ( k > nzb_w_inner(j,i) )  THEN
    253                          tend(k,j,i) = tend(k,j,i) + fs * 0.25_wp *            &
    254                                       ( u(k,j,i) + u(k+1,j,i) + u(k,j,i+1) +   &
    255                                         u(k+1,j,i+1) )
    256                       ENDIF
    257                    ENDDO
    258                 ENDDO
    259              ENDDO
    260              !$acc end kernels
    261 
    262           CASE DEFAULT
    263 
    264              WRITE( message_string, * ) ' wrong component: ', component
    265              CALL message( 'coriolis', 'PA0173', 1, 2, 0, 6, 0 )
    266 
    267        END SELECT
    268 
    269     END SUBROUTINE coriolis_acc
    270 
    271 
    272 !------------------------------------------------------------------------------!
    273 ! Description:
    274 ! ------------
    275174!> Call for grid point i,j
    276175!------------------------------------------------------------------------------!
  • palm/trunk/SOURCE/cpulog_mod.f90

    r2101 r2118  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! OpenACC relevant code removed
    2323!
    2424! Former revisions:
     
    401401                             average_cputime
    402402
    403           IF ( num_acc_per_node /= 0 )  WRITE ( 18, 108 )  num_acc_per_node
    404403          WRITE ( 18, 110 )
    405404#else
     
    409408                             average_cputime
    410409
    411           IF ( num_acc_per_node /= 0 )  WRITE ( 18, 109 )  num_acc_per_node
    412410          WRITE ( 18, 110 )
    413411#endif
     
    565563   106 FORMAT (/'Exchange of ghostpoints via MPI_ISEND/MPI_IRECV')
    566564   107 FORMAT (//)
    567    108 FORMAT ('Accelerator boards per node: ',14X,I2)
    568    109 FORMAT ('Accelerator boards: ',23X,I2)
    569565   110 FORMAT ('-------------------------------------------------------------',     &
    570566               &'---------'//&
  • palm/trunk/SOURCE/diffusion_e.f90

    r2101 r2118  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! OpenACC version of subroutine removed
    2323!
    2424! Former revisions:
     
    3434! 1873 2016-04-18 14:50:06Z maronga
    3535! Module renamed (removed _mod)
    36 !
    3736!
    3837! 1850 2016-04-08 13:29:27Z maronga
     
    108107
    109108    PRIVATE
    110     PUBLIC diffusion_e, diffusion_e_acc
     109    PUBLIC diffusion_e
    111110   
    112111
     
    116115    END INTERFACE diffusion_e
    117116 
    118     INTERFACE diffusion_e_acc
    119        MODULE PROCEDURE diffusion_e_acc
    120     END INTERFACE diffusion_e_acc
    121 
    122117 CONTAINS
    123118
     
    337332
    338333    END SUBROUTINE diffusion_e
    339 
    340 
    341 !------------------------------------------------------------------------------!
    342 ! Description:
    343 ! ------------
    344 !> Call for all grid points - accelerator version
    345 !------------------------------------------------------------------------------!
    346     SUBROUTINE diffusion_e_acc( var, var_reference )
    347 
    348        USE arrays_3d,                                                          &
    349            ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw,        &
    350                   drho_air, rho_air_zw
    351          
    352        USE control_parameters,                                                 &
    353            ONLY:  atmos_ocean_sign, g, use_single_reference_value,             &
    354                   wall_adjustment, wall_adjustment_factor
    355 
    356        USE grid_variables,                                                     &
    357            ONLY:  ddx2, ddy2
    358            
    359        USE indices,                                                            &
    360            ONLY:  i_left, i_right, j_north, j_south, nxlg, nxrg, nyng, nysg,   &
    361                   nzb, nzb_s_inner, nzt
    362            
    363        USE kinds
    364 
    365        USE microphysics_mod,                                                   &
    366            ONLY:  collision_turbulence
    367 
    368        USE particle_attributes,                                                &
    369            ONLY:  use_sgs_for_particles, wang_kernel
    370 
    371        IMPLICIT NONE
    372 
    373        INTEGER(iwp) ::  i              !<
    374        INTEGER(iwp) ::  j              !<
    375        INTEGER(iwp) ::  k              !<
    376        REAL(wp)     ::  dissipation    !<
    377        REAL(wp)     ::  dvar_dz        !<
    378        REAL(wp)     ::  l              !<
    379        REAL(wp)     ::  ll             !<
    380        REAL(wp)     ::  l_stable       !<
    381        REAL(wp)     ::  var_reference  !<
    382 
    383 #if defined( __nopointer )
    384        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !<
    385 #else
    386        REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !<
    387 #endif
    388 
    389 
    390 !
    391 !--    This if clause must be outside the k-loop because otherwise
    392 !--    runtime errors occur with -C hopt on NEC
    393        IF ( use_single_reference_value )  THEN
    394 
    395           !$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
    396           !$acc         present( nzb_s_inner, tend, var, zu, zw )
    397           DO  i = i_left, i_right
    398              DO  j = j_south, j_north
    399                 DO  k = 1, nzt
    400 
    401                    IF ( k > nzb_s_inner(j,i) )  THEN
    402 !
    403 !--                   Calculate the mixing length (for dissipation)
    404                       dvar_dz = atmos_ocean_sign * &
    405                                 ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
    406                       IF ( dvar_dz > 0.0_wp ) THEN
    407                          l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
    408                                        SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
    409                       ELSE
    410                          l_stable = l_grid(k)
    411                       ENDIF
    412 !
    413 !--                   Adjustment of the mixing length
    414                       IF ( wall_adjustment )  THEN
    415                          l  = MIN( wall_adjustment_factor *          &
    416                                    ( zu(k) - zw(nzb_s_inner(j,i)) ), &
    417                                    l_grid(k), l_stable )
    418                          ll = MIN( wall_adjustment_factor *          &
    419                                    ( zu(k) - zw(nzb_s_inner(j,i)) ), &
    420                                    l_grid(k) )
    421                       ELSE
    422                          l  = MIN( l_grid(k), l_stable )
    423                          ll = l_grid(k)
    424                       ENDIF
    425 !
    426 !--                   Calculate the tendency terms
    427                       dissipation = ( 0.19_wp + 0.74_wp * l / ll ) * &
    428                                           e(k,j,i) * SQRT( e(k,j,i) ) / l
    429 
    430                       tend(k,j,i) = tend(k,j,i)                                  &
    431                                         + (                                    &
    432                           ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
    433                         - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
    434                                           ) * ddx2                             &
    435                                         + (                                    &
    436                           ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
    437                         - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
    438                                           ) * ddy2                             &
    439                                         + (                                    &
    440                ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
    441                                                              * rho_air_zw(k)   &
    442              - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
    443                                                              * rho_air_zw(k-1) &
    444                                           ) * ddzw(k) * drho_air(k)            &
    445                                   - dissipation
    446 
    447 !
    448 !--                   Store dissipation if needed for calculating the sgs particle
    449 !--                   velocities
    450                       IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.      &
    451                            collision_turbulence )  THEN
    452                          diss(k,j,i) = dissipation
    453                       ENDIF
    454 
    455                    ENDIF
    456 
    457                 ENDDO
    458              ENDDO
    459           ENDDO
    460           !$acc end kernels
    461 
    462        ELSE
    463 
    464           !$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
    465           !$acc         present( nzb_s_inner, tend, var, zu, zw )
    466           DO  i = i_left, i_right
    467              DO  j = j_south, j_north
    468                 DO  k = 1, nzt
    469 
    470                    IF ( k > nzb_s_inner(j,i) )  THEN
    471 !
    472 !--                   Calculate the mixing length (for dissipation)
    473                       dvar_dz = atmos_ocean_sign * &
    474                                 ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
    475                       IF ( dvar_dz > 0.0_wp ) THEN
    476                          l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
    477                                               SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
    478                       ELSE
    479                          l_stable = l_grid(k)
    480                       ENDIF
    481 !
    482 !--                   Adjustment of the mixing length
    483                       IF ( wall_adjustment )  THEN
    484                          l  = MIN( wall_adjustment_factor *          &
    485                                    ( zu(k) - zw(nzb_s_inner(j,i)) ), &
    486                                      l_grid(k), l_stable )
    487                          ll = MIN( wall_adjustment_factor *          &
    488                                    ( zu(k) - zw(nzb_s_inner(j,i)) ), &
    489                                    l_grid(k) )
    490                       ELSE
    491                          l  = MIN( l_grid(k), l_stable )
    492                          ll = l_grid(k)
    493                       ENDIF
    494 !
    495 !--                   Calculate the tendency terms
    496                       dissipation = ( 0.19_wp + 0.74_wp * l / ll ) * &
    497                                           e(k,j,i) * SQRT( e(k,j,i) ) / l
    498 
    499                       tend(k,j,i) = tend(k,j,i)                                &
    500                                         + (                                    &
    501                           ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
    502                         - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
    503                                           ) * ddx2                             &
    504                                         + (                                    &
    505                           ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
    506                         - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
    507                                           ) * ddy2                             &
    508                                         + (                                    &
    509                ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
    510                                                              * rho_air_zw(k)   &
    511              - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
    512                                                              * rho_air_zw(k-1) &
    513                                           ) * ddzw(k) * drho_air(k)            &
    514                                   - dissipation
    515 
    516 !
    517 !--                   Store dissipation if needed for calculating the sgs
    518 !--                   particle  velocities
    519                       IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.      &
    520                            collision_turbulence )  THEN
    521                          diss(k,j,i) = dissipation
    522                       ENDIF
    523 
    524                    ENDIF
    525 
    526                 ENDDO
    527              ENDDO
    528           ENDDO
    529           !$acc end kernels
    530 
    531        ENDIF
    532 
    533 !
    534 !--    Boundary condition for dissipation
    535        IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.                     &
    536             collision_turbulence )  THEN
    537           !$acc kernels present( diss, nzb_s_inner )
    538           DO  i = i_left, i_right
    539              DO  j = j_south, j_north
    540                 diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
    541              ENDDO
    542           ENDDO
    543           !$acc end kernels
    544        ENDIF
    545 
    546     END SUBROUTINE diffusion_e_acc
    547334
    548335
  • palm/trunk/SOURCE/diffusion_s.f90

    r2101 r2118  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! OpenACC version of subroutine removed
    2323!
    2424! Former revisions:
     
    3434! 1873 2016-04-18 14:50:06Z maronga
    3535! Module renamed (removed _mod)
    36 !
    37 !
     36!
    3837! 1850 2016-04-08 13:29:27Z maronga
    3938! Module renamed
    40 !
    4139!
    4240! 1691 2015-10-26 16:17:44Z maronga
     
    9492
    9593    PRIVATE
    96     PUBLIC diffusion_s, diffusion_s_acc
     94    PUBLIC diffusion_s
    9795
    9896    INTERFACE diffusion_s
     
    10098       MODULE PROCEDURE diffusion_s_ij
    10199    END INTERFACE diffusion_s
    102 
    103     INTERFACE diffusion_s_acc
    104        MODULE PROCEDURE diffusion_s_acc
    105     END INTERFACE diffusion_s_acc
    106100
    107101 CONTAINS
     
    242236! Description:
    243237! ------------
    244 !> Call for all grid points - accelerator version
    245 !------------------------------------------------------------------------------!
    246     SUBROUTINE diffusion_s_acc( s, s_flux_b, s_flux_t, wall_s_flux )
    247 
    248        USE arrays_3d,                                                          &
    249            ONLY:  ddzu, ddzw, kh, tend, drho_air, rho_air_zw
    250            
    251        USE control_parameters,                                                 &
    252            ONLY: use_surface_fluxes, use_top_fluxes
    253        
    254        USE grid_variables,                                                     &
    255            ONLY:  ddx2, ddy2, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y
    256        
    257        USE indices, &
    258            ONLY: i_left, i_right, j_north, j_south, nxlg, nxrg, nyng, nysg,    &
    259                  nzb, nzb_diff_s_inner, nzb_s_inner, nzb_s_outer, nzt, nzt_diff
    260            
    261        USE kinds
    262 
    263        IMPLICIT NONE
    264 
    265        INTEGER(iwp) ::  i                 !<
    266        INTEGER(iwp) ::  j                 !<
    267        INTEGER(iwp) ::  k                 !<
    268        REAL(wp)     ::  wall_s_flux(0:4)  !<
    269        REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_b !<
    270        REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_t !<
    271 #if defined( __nopointer )
    272        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s  !<
    273 #else
    274        REAL(wp), DIMENSION(:,:,:), POINTER ::  s  !<
    275 #endif
    276 
    277        !$acc kernels present( ddzu, ddzw, fwxm, fwxp, fwym, fwyp, kh )        &
    278        !$acc         present( nzb_diff_s_inner, nzb_s_inner, nzb_s_outer, s ) &
    279        !$acc         present( s_flux_b, s_flux_t, tend, wall_s_flux )         &
    280        !$acc         present( wall_w_x, wall_w_y )
    281        DO  i = i_left, i_right
    282           DO  j = j_south, j_north
    283 !
    284 !--          Compute horizontal diffusion
    285              DO  k = 1, nzt
    286                 IF ( k > nzb_s_outer(j,i) )  THEN
    287 
    288                    tend(k,j,i) = tend(k,j,i)                                   &
    289                                           + 0.5_wp * (                         &
    290                         ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
    291                       - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
    292                                                      ) * ddx2                  &
    293                                           + 0.5_wp * (                         &
    294                         ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
    295                       - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
    296                                                      ) * ddy2
    297                 ENDIF
    298              ENDDO
    299 
    300 !
    301 !--          Apply prescribed horizontal wall heatflux where necessary
    302              DO  k = 1, nzt
    303                 IF ( k > nzb_s_inner(j,i)  .AND.  k <= nzb_s_outer(j,i)  .AND. &
    304                      ( wall_w_x(j,i) /= 0.0_wp  .OR.  wall_w_y(j,i) /= 0.0_wp ) )    &
    305                 THEN
    306                    tend(k,j,i) = tend(k,j,i)                                   &
    307                                                 + ( fwxp(j,i) * 0.5_wp *       &
    308                         ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
    309                         + ( 1.0_wp - fwxp(j,i) ) * wall_s_flux(1)              &
    310                                                    -fwxm(j,i) * 0.5_wp *       &
    311                         ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
    312                         + ( 1.0_wp - fwxm(j,i) ) * wall_s_flux(2)              &
    313                                                   ) * ddx2                     &
    314                                                 + ( fwyp(j,i) * 0.5_wp *       &
    315                         ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
    316                         + ( 1.0_wp - fwyp(j,i) ) * wall_s_flux(3)              &
    317                                                    -fwym(j,i) * 0.5_wp *       &
    318                         ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
    319                         + ( 1.0_wp - fwym(j,i) ) * wall_s_flux(4)              &
    320                                                   ) * ddy2
    321                 ENDIF
    322              ENDDO
    323 
    324 !
    325 !--          Compute vertical diffusion. In case that surface fluxes have been
    326 !--          prescribed or computed at bottom and/or top, index k starts/ends at
    327 !--          nzb+2 or nzt-1, respectively.
    328              DO  k = 1, nzt_diff
    329                 IF ( k >= nzb_diff_s_inner(j,i) )  THEN
    330                    tend(k,j,i) = tend(k,j,i)                                   &
    331                                        + 0.5_wp * (                            &
    332             ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
    333                                                             * rho_air_zw(k)    &
    334           - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
    335                                                             * rho_air_zw(k-1)  &
    336                                                   ) * ddzw(k) * drho_air(k)
    337                 ENDIF
    338              ENDDO
    339 
    340 !
    341 !--          Vertical diffusion at the first computational gridpoint along
    342 !--          z-direction
    343              DO  k = 1, nzt
    344                 IF ( use_surface_fluxes  .AND.  k == nzb_s_inner(j,i)+1 )  THEN
    345                    tend(k,j,i) = tend(k,j,i)                                   &
    346                                           + ( 0.5_wp * ( kh(k,j,i)+kh(k+1,j,i) )&
    347                                                      * ( s(k+1,j,i)-s(k,j,i) ) &
    348                                                      * ddzu(k+1)               &
    349                                                      * rho_air_zw(k)           &
    350                                               + s_flux_b(j,i)                  &
    351                                             ) * ddzw(k) * drho_air(k)
    352                 ENDIF
    353 
    354 !
    355 !--             Vertical diffusion at the last computational gridpoint along
    356 !--             z-direction
    357                 IF ( use_top_fluxes  .AND.  k == nzt )  THEN
    358                    tend(k,j,i) = tend(k,j,i)                                   &
    359                                           + ( - s_flux_t(j,i)                  &
    360                                               - 0.5_wp * ( kh(k-1,j,i)+kh(k,j,i) )&
    361                                                        * ( s(k,j,i)-s(k-1,j,i) )  &
    362                                                        * ddzu(k)                  &
    363                                                        * rho_air_zw(k-1)          &
    364                                             ) * ddzw(k) * drho_air(k)
    365                 ENDIF
    366              ENDDO
    367 
    368           ENDDO
    369        ENDDO
    370        !$acc end kernels
    371 
    372     END SUBROUTINE diffusion_s_acc
    373 
    374 
    375 !------------------------------------------------------------------------------!
    376 ! Description:
    377 ! ------------
    378238!> Call for grid point i,j
    379239!------------------------------------------------------------------------------!
  • palm/trunk/SOURCE/diffusion_u.f90

    r2101 r2118  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! OpenACC version of subroutine removed
    2323!
    2424! Former revisions:
     
    3535! Module renamed (removed _mod)
    3636!
    37 !
    3837! 1850 2016-04-08 13:29:27Z maronga
    3938! Module renamed
    40 !
    4139!
    4240! 1740 2016-01-13 08:19:40Z raasch
     
    9795
    9896    PRIVATE
    99     PUBLIC diffusion_u, diffusion_u_acc
     97    PUBLIC diffusion_u
    10098
    10199    INTERFACE diffusion_u
     
    103101       MODULE PROCEDURE diffusion_u_ij
    104102    END INTERFACE diffusion_u
    105 
    106     INTERFACE diffusion_u_acc
    107        MODULE PROCEDURE diffusion_u_acc
    108     END INTERFACE diffusion_u_acc
    109103
    110104 CONTAINS
     
    280274! Description:
    281275! ------------
    282 !> Call for all grid points - accelerator version
    283 !------------------------------------------------------------------------------!
    284     SUBROUTINE diffusion_u_acc
    285 
    286        USE arrays_3d,                                                          &
    287            ONLY:  ddzu, ddzw, km, tend, u, usws, uswst, v, w,                  &
    288                   drho_air, rho_air_zw
    289        
    290        USE control_parameters,                                                 &
    291            ONLY:  constant_top_momentumflux, topography, use_surface_fluxes,   &
    292                   use_top_fluxes
    293        
    294        USE grid_variables,                                                     &
    295            ONLY:  ddx, ddx2, ddy, fym, fyp, wall_u
    296        
    297        USE indices,                                                            &
    298            ONLY:  i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb,  &
    299                   nzb_diff_u, nzb_u_inner, nzb_u_outer, nzt, nzt_diff
    300        
    301        USE kinds
    302 
    303        IMPLICIT NONE
    304 
    305        INTEGER(iwp) ::  i     !<
    306        INTEGER(iwp) ::  j     !<
    307        INTEGER(iwp) ::  k     !<
    308        REAL(wp)     ::  kmym  !<
    309        REAL(wp)     ::  kmyp  !<
    310        REAL(wp)     ::  kmzm  !<
    311        REAL(wp)     ::  kmzp  !<
    312 
    313        REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs  !<
    314        !$acc declare create ( usvs )
    315 
    316 !
    317 !--    First calculate horizontal momentum flux u'v' at vertical walls,
    318 !--    if neccessary
    319        IF ( topography /= 'flat' )  THEN
    320           CALL wall_fluxes_acc( usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,          &
    321                                 nzb_u_inner, nzb_u_outer, wall_u )
    322        ENDIF
    323 
    324        !$acc kernels present ( u, v, w, km, tend, usws, uswst )                &
    325        !$acc         present ( ddzu, ddzw, fym, fyp, wall_u )                  &
    326        !$acc         present ( nzb_u_inner, nzb_u_outer, nzb_diff_u )
    327        DO  i = i_left, i_right
    328           DO  j = j_south, j_north
    329 !
    330 !--          Compute horizontal diffusion
    331              DO  k = 1, nzt
    332                 IF ( k > nzb_u_outer(j,i) )  THEN
    333 !
    334 !--                Interpolate eddy diffusivities on staggered gridpoints
    335                    kmyp = 0.25_wp *                                            &
    336                           ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
    337                    kmym = 0.25_wp *                                            &
    338                           ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
    339 
    340                    tend(k,j,i) = tend(k,j,i)                                   &
    341                          & + 2.0_wp * (                                        &
    342                          &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   ) &
    343                          &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) ) &
    344                          &            ) * ddx2                                 &
    345                          & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy      &
    346                          &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx      &
    347                          &   - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy          &
    348                          &   - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx          &
    349                          &   ) * ddy
    350                 ENDIF
    351              ENDDO
    352 
    353 !
    354 !--          Wall functions at the north and south walls, respectively
    355              DO  k = 1, nzt
    356                 IF( k > nzb_u_inner(j,i)  .AND.  k <= nzb_u_outer(j,i)  .AND.  &
    357                     wall_u(j,i) /= 0.0_wp )  THEN
    358 
    359                    kmyp = 0.25_wp *                                            &
    360                           ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
    361                    kmym = 0.25_wp *                                            &
    362                           ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
    363 
    364                    tend(k,j,i) = tend(k,j,i)                                   &
    365                                  + 2.0_wp * (                                  &
    366                                        km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
    367                                      - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
    368                                             ) * ddx2                           &
    369                                  + (   fyp(j,i) * (                            &
    370                                   kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
    371                                 + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx   &
    372                                                   )                            &
    373                                      - fym(j,i) * (                            &
    374                                   kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
    375                                 + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
    376                                                   )                            &
    377                                      + wall_u(j,i) * usvs(k,j,i)               &
    378                                    ) * ddy
    379                 ENDIF
    380              ENDDO
    381 
    382 !
    383 !--          Compute vertical diffusion. In case of simulating a Prandtl layer,
    384 !--          index k starts at nzb_u_inner+2.
    385              DO  k = 1, nzt_diff
    386                 IF ( k >= nzb_diff_u(j,i) )  THEN
    387 !
    388 !--                Interpolate eddy diffusivities on staggered gridpoints
    389                    kmzp = 0.25_wp *                                            &
    390                           ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
    391                    kmzm = 0.25_wp *                                            &
    392                           ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
    393 
    394                    tend(k,j,i) = tend(k,j,i)                                   &
    395                          & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)&
    396                          &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx      &
    397                          &            ) * rho_air_zw(k)                        &
    398                          &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)&
    399                          &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx    &
    400                          &            ) * rho_air_zw(k-1)                      &
    401                          &   ) * ddzw(k) * drho_air(k)
    402                 ENDIF
    403              ENDDO
    404 
    405           ENDDO
    406        ENDDO
    407 
    408 !
    409 !--    Vertical diffusion at the first grid point above the surface,
    410 !--    if the momentum flux at the bottom is given by the Prandtl law or
    411 !--    if it is prescribed by the user.
    412 !--    Difference quotient of the momentum flux is not formed over half
    413 !--    of the grid spacing (2.0*ddzw(k)) any more, since the comparison
    414 !--    with other (LES) models showed that the values of the momentum
    415 !--    flux becomes too large in this case.
    416 !--    The term containing w(k-1,..) (see above equation) is removed here
    417 !--    because the vertical velocity is assumed to be zero at the surface.
    418        IF ( use_surface_fluxes )  THEN
    419 
    420           DO  i = i_left, i_right
    421              DO  j = j_south, j_north
    422          
    423                 k = nzb_u_inner(j,i)+1
    424 !
    425 !--             Interpolate eddy diffusivities on staggered gridpoints
    426                 kmzp = 0.25_wp *                                               &
    427                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
    428 
    429                 tend(k,j,i) = tend(k,j,i)                                      &
    430                       & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
    431                       &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
    432                       &            ) * rho_air_zw(k)                           &
    433                       &   - ( -usws(j,i) )                                     &
    434                       &   ) * ddzw(k) * drho_air(k)
    435              ENDDO
    436           ENDDO
    437 
    438        ENDIF
    439 
    440 !
    441 !--    Vertical diffusion at the first gridpoint below the top boundary,
    442 !--    if the momentum flux at the top is prescribed by the user
    443        IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
    444 
    445           k = nzt
    446 
    447           DO  i = i_left, i_right
    448              DO  j = j_south, j_north
    449 
    450 !
    451 !--             Interpolate eddy diffusivities on staggered gridpoints
    452                 kmzm = 0.25_wp *                                               &
    453                        ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
    454 
    455                 tend(k,j,i) = tend(k,j,i)                                      &
    456                       & + ( ( -uswst(j,i) )                                    &
    457                       &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
    458                       &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
    459                       &            ) * rho_air_zw(k-1)                         &
    460                       &   ) * ddzw(k) * drho_air(k)
    461              ENDDO
    462           ENDDO
    463 
    464        ENDIF
    465        !$acc end kernels
    466 
    467     END SUBROUTINE diffusion_u_acc
    468 
    469 
    470 !------------------------------------------------------------------------------!
    471 ! Description:
    472 ! ------------
    473276!> Call for grid point i,j
    474277!------------------------------------------------------------------------------!
  • palm/trunk/SOURCE/diffusion_v.f90

    r2101 r2118  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! OpenACC version of subroutine removed
    2323!
    2424! Former revisions:
     
    3535! Module renamed (removed _mod)
    3636!
    37 !
    3837! 1850 2016-04-08 13:29:27Z maronga
    3938! Module renamed
    40 !
    4139!
    4240! 1740 2016-01-13 08:19:40Z raasch
     
    9290
    9391    PRIVATE
    94     PUBLIC diffusion_v, diffusion_v_acc
     92    PUBLIC diffusion_v
    9593
    9694    INTERFACE diffusion_v
     
    9896       MODULE PROCEDURE diffusion_v_ij
    9997    END INTERFACE diffusion_v
    100 
    101     INTERFACE diffusion_v_acc
    102        MODULE PROCEDURE diffusion_v_acc
    103     END INTERFACE diffusion_v_acc
    10498
    10599 CONTAINS
     
    275269! Description:
    276270! ------------
    277 !> Call for all grid points - accelerator version
    278 !------------------------------------------------------------------------------!
    279     SUBROUTINE diffusion_v_acc
    280 
    281        USE arrays_3d,                                                          &
    282            ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w,                  &
    283                   drho_air, rho_air_zw
    284        
    285        USE control_parameters,                                                 &
    286            ONLY:  constant_top_momentumflux, topography, use_surface_fluxes,   &
    287                   use_top_fluxes
    288        
    289        USE grid_variables,                                                     &
    290            ONLY:  ddx, ddy, ddy2, fxm, fxp, wall_v
    291        
    292        USE indices,                                                            &
    293            ONLY:  i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb,  &
    294                   nzb_diff_v, nzb_v_inner, nzb_v_outer, nzt, nzt_diff
    295        
    296        USE kinds
    297 
    298        IMPLICIT NONE
    299 
    300        INTEGER(iwp) ::  i     !<
    301        INTEGER(iwp) ::  j     !<
    302        INTEGER(iwp) ::  k     !<
    303        REAL(wp)     ::  kmxm  !<
    304        REAL(wp)     ::  kmxp  !<
    305        REAL(wp)     ::  kmzm  !<
    306        REAL(wp)     ::  kmzp  !<
    307 
    308        REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus  !<
    309        !$acc declare create ( vsus )
    310 
    311 !
    312 !--    First calculate horizontal momentum flux v'u' at vertical walls,
    313 !--    if neccessary
    314        IF ( topography /= 'flat' )  THEN
    315           CALL wall_fluxes_acc( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp,          &
    316                                 nzb_v_inner, nzb_v_outer, wall_v )
    317        ENDIF
    318 
    319        !$acc kernels present ( u, v, w, km, tend, vsws, vswst )                &
    320        !$acc         present ( ddzu, ddzw, fxm, fxp, wall_v )                  &
    321        !$acc         present ( nzb_v_inner, nzb_v_outer, nzb_diff_v )
    322        DO  i = i_left, i_right
    323           DO  j = j_south, j_north
    324 !
    325 !--          Compute horizontal diffusion
    326              DO  k = 1, nzt
    327                 IF ( k > nzb_v_outer(j,i) )  THEN
    328 !
    329 !--                Interpolate eddy diffusivities on staggered gridpoints
    330                    kmxp = 0.25_wp *                                            &
    331                           ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
    332                    kmxm = 0.25_wp *                                            &
    333                           ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    334 
    335                    tend(k,j,i) = tend(k,j,i)                                   &
    336                          & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx      &
    337                          &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy      &
    338                          &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx          &
    339                          &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy          &
    340                          &   ) * ddx                                           &
    341                          & + 2.0_wp * (                                        &
    342                          &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )   &
    343                          &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )   &
    344                          &            ) * ddy2
    345                 ENDIF
    346              ENDDO
    347 
    348 !
    349 !--          Wall functions at the left and right walls, respectively
    350              DO  k = 1, nzt
    351                 IF( k > nzb_v_inner(j,i)  .AND.  k <= nzb_v_outer(j,i)  .AND.  &
    352                     wall_v(j,i) /= 0.0_wp )  THEN
    353 
    354                    kmxp = 0.25_wp *                                            &
    355                           ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
    356                    kmxm = 0.25_wp *                                            &
    357                           ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    358                    
    359                    tend(k,j,i) = tend(k,j,i)                                   &
    360                                  + 2.0_wp * (                                  &
    361                                        km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
    362                                      - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
    363                                             ) * ddy2                           &
    364                                  + (   fxp(j,i) * (                            &
    365                                   kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
    366                                 + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
    367                                                   )                            &
    368                                      - fxm(j,i) * (                            &
    369                                   kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
    370                                 + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
    371                                                   )                            &
    372                                      + wall_v(j,i) * vsus(k,j,i)               &
    373                                    ) * ddx
    374                 ENDIF
    375              ENDDO
    376 
    377 !
    378 !--          Compute vertical diffusion. In case of simulating a Prandtl
    379 !--          layer, index k starts at nzb_v_inner+2.
    380              DO  k = 1, nzt_diff
    381                 IF ( k >= nzb_diff_v(j,i) )  THEN
    382 !
    383 !--                Interpolate eddy diffusivities on staggered gridpoints
    384                    kmzp = 0.25_wp *                                            &
    385                           ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    386                    kmzm = 0.25_wp *                                            &
    387                           ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
    388 
    389                    tend(k,j,i) = tend(k,j,i)                                   &
    390                          & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)&
    391                          &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy      &
    392                          &            ) * rho_air_zw(k)                        &
    393                          &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)&
    394                          &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy    &
    395                          &            ) * rho_air_zw(k-1)                      &
    396                          &   ) * ddzw(k) * drho_air(k)
    397                 ENDIF
    398              ENDDO
    399 
    400           ENDDO
    401        ENDDO
    402 
    403 !
    404 !--    Vertical diffusion at the first grid point above the surface,
    405 !--    if the momentum flux at the bottom is given by the Prandtl law
    406 !--    or if it is prescribed by the user.
    407 !--    Difference quotient of the momentum flux is not formed over
    408 !--    half of the grid spacing (2.0*ddzw(k)) any more, since the
    409 !--    comparison with other (LES) models showed that the values of
    410 !--    the momentum flux becomes too large in this case.
    411 !--    The term containing w(k-1,..) (see above equation) is removed here
    412 !--    because the vertical velocity is assumed to be zero at the surface.
    413        IF ( use_surface_fluxes )  THEN
    414 
    415           DO  i = i_left, i_right
    416              DO  j = j_south, j_north
    417          
    418                 k = nzb_v_inner(j,i)+1
    419 !
    420 !--             Interpolate eddy diffusivities on staggered gridpoints
    421                 kmzp = 0.25_wp *                                               &
    422                        ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    423 
    424                 tend(k,j,i) = tend(k,j,i)                                      &
    425                       & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
    426                       &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
    427                       &            ) * rho_air_zw(k)                           &
    428                       &   - ( -vsws(j,i) )                                     &
    429                       &   ) * ddzw(k) * drho_air(k)
    430              ENDDO
    431           ENDDO
    432 
    433        ENDIF
    434 
    435 !
    436 !--    Vertical diffusion at the first gridpoint below the top boundary,
    437 !--    if the momentum flux at the top is prescribed by the user
    438        IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
    439 
    440           k = nzt
    441 
    442           DO  i = i_left, i_right
    443              DO  j = j_south, j_north
    444 
    445 !
    446 !--             Interpolate eddy diffusivities on staggered gridpoints
    447                 kmzm = 0.25_wp *                                               &
    448                        ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
    449 
    450                 tend(k,j,i) = tend(k,j,i)                                      &
    451                       & + ( ( -vswst(j,i) )                                    &
    452                       &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
    453                       &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
    454                       &            ) * rho_air_zw(k-1)                         &
    455                       &   ) * ddzw(k) * drho_air(k)
    456              ENDDO
    457           ENDDO
    458 
    459        ENDIF
    460        !$acc end kernels
    461 
    462     END SUBROUTINE diffusion_v_acc
    463 
    464 
    465 !------------------------------------------------------------------------------!
    466 ! Description:
    467 ! ------------
    468271!> Call for grid point i,j
    469272!------------------------------------------------------------------------------!
  • palm/trunk/SOURCE/diffusion_w.f90

    r2101 r2118  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! OpenACC version of subroutine removed
    2323!
    2424! Former revisions:
     
    3535! Module renamed (removed _mod)
    3636!
    37 !
    3837! 1850 2016-04-08 13:29:27Z maronga
    3938! Module renamed
    40 !
    4139!
    4240! 1682 2015-10-07 23:56:08Z knoop
     
    9795
    9896    USE wall_fluxes_mod,                                                       &
    99         ONLY :  wall_fluxes, wall_fluxes_acc
     97        ONLY :  wall_fluxes
    10098
    10199    PRIVATE
    102     PUBLIC diffusion_w, diffusion_w_acc
     100    PUBLIC diffusion_w
    103101
    104102    INTERFACE diffusion_w
     
    106104       MODULE PROCEDURE diffusion_w_ij
    107105    END INTERFACE diffusion_w
    108 
    109     INTERFACE diffusion_w_acc
    110        MODULE PROCEDURE diffusion_w_acc
    111     END INTERFACE diffusion_w_acc
    112106
    113107 CONTAINS
     
    248242! Description:
    249243! ------------
    250 !> Call for all grid points - accelerator version
    251 !------------------------------------------------------------------------------!
    252     SUBROUTINE diffusion_w_acc
    253 
    254        USE arrays_3d,                                                          &
    255            ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
    256            
    257        USE control_parameters,                                                 &
    258            ONLY :  topography
    259            
    260        USE grid_variables,                                                     &
    261            ONLY : ddx, ddy, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y
    262            
    263        USE indices,                                                            &
    264            ONLY :  i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb, &
    265                    nzb_w_inner, nzb_w_outer, nzt
    266                    
    267        USE kinds
    268 
    269        IMPLICIT NONE
    270 
    271        INTEGER(iwp) ::  i     !<
    272        INTEGER(iwp) ::  j     !<
    273        INTEGER(iwp) ::  k     !<
    274        
    275        REAL(wp) ::  kmxm  !<
    276        REAL(wp) ::  kmxp  !<
    277        REAL(wp) ::  kmym  !<
    278        REAL(wp) ::  kmyp  !<
    279 
    280        REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus  !<
    281        REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsvs  !<
    282        !$acc declare create ( wsus, wsvs )
    283 
    284 !
    285 !--    First calculate horizontal momentum flux w'u' and/or w'v' at vertical
    286 !--    walls, if neccessary
    287        IF ( topography /= 'flat' )  THEN
    288           CALL wall_fluxes_acc( wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp,          &
    289                                 nzb_w_inner, nzb_w_outer, wall_w_x )
    290           CALL wall_fluxes_acc( wsvs, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp,          &
    291                                 nzb_w_inner, nzb_w_outer, wall_w_y )
    292        ENDIF
    293 
    294        !$acc kernels present ( u, v, w, km, tend )                             &
    295        !$acc         present ( ddzu, ddzw, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y )           &
    296        !$acc         present ( nzb_w_inner, nzb_w_outer )
    297        DO  i = i_left, i_right
    298           DO  j = j_south, j_north
    299              DO  k = 1, nzt
    300                 IF ( k > nzb_w_outer(j,i) )  THEN
    301 !
    302 !--                Interpolate eddy diffusivities on staggered gridpoints
    303                    kmxp = 0.25_wp *                                            &
    304                           ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
    305                    kmxm = 0.25_wp *                                            &
    306                           ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
    307                    kmyp = 0.25_wp *                                            &
    308                           ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
    309                    kmym = 0.25_wp *                                            &
    310                           ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    311 
    312                    tend(k,j,i) = tend(k,j,i)                                     &
    313                          & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx        &
    314                          &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)  &
    315                          &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx          &
    316                          &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)    &
    317                          &   ) * ddx                                             &
    318                          & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy        &
    319                          &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)  &
    320                          &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy          &
    321                          &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)    &
    322                          &   ) * ddy                                             &
    323                          & + 2.0_wp * (                                          &
    324                          &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
    325                          &               * rho_air(k+1)                          &
    326                          & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
    327                          &               * rho_air(k)                            &
    328                          &            ) * ddzu(k+1) * drho_air_zw(k)
    329                 ENDIF
    330              ENDDO
    331 
    332 !
    333 !--          Wall functions at all vertical walls, where necessary
    334              DO  k = 1,nzt
    335 
    336                 IF ( k > nzb_w_inner(j,i)  .AND.  k <= nzb_w_outer(j,i)  .AND. &
    337                      wall_w_x(j,i) /= 0.0_wp  .AND.  wall_w_y(j,i) /= 0.0_wp )  THEN
    338 !
    339 !--                Interpolate eddy diffusivities on staggered gridpoints
    340                    kmxp = 0.25_wp *                                            &
    341                           ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
    342                    kmxm = 0.25_wp *                                            &
    343                           ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
    344                    kmyp = 0.25_wp *                                            &
    345                           ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
    346                    kmym = 0.25_wp *                                            &
    347                           ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    348 
    349                    tend(k,j,i) = tend(k,j,i)                                   &
    350                                  + (   fwxp(j,i) * (                           &
    351                             kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
    352                           + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
    353                                                    )                           &
    354                                      - fwxm(j,i) * (                           &
    355                             kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
    356                           + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
    357                                                    )                           &
    358                                      + wall_w_x(j,i) * wsus(k,j,i)             &
    359                                    ) * ddx                                     &
    360                                  + (   fwyp(j,i) * (                           &
    361                             kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
    362                           + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
    363                                                    )                           &
    364                                      - fwym(j,i) * (                           &
    365                             kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
    366                           + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
    367                                                    )                           &
    368                                      + wall_w_y(j,i) * wsvs(k,j,i)             &
    369                                    ) * ddy                                     &
    370                                  + 2.0_wp * (                                  &
    371                            km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
    372                                        * rho_air(k+1)                          &
    373                          - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
    374                                        * rho_air(k)                            &
    375                                             ) * ddzu(k+1) * drho_air_zw(k)
    376                 ENDIF
    377              ENDDO
    378 
    379           ENDDO
    380        ENDDO
    381        !$acc end kernels
    382 
    383     END SUBROUTINE diffusion_w_acc
    384 
    385 
    386 !------------------------------------------------------------------------------!
    387 ! Description:
    388 ! ------------
    389244!> Call for grid point i,j
    390245!------------------------------------------------------------------------------!
  • palm/trunk/SOURCE/diffusivities.f90

    r2101 r2118  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! OpenACC directives removed
    2323!
    2424! Former revisions:
     
    130130
    131131!
    132 !-- Data declerations for accelerators
    133     !$acc data present( dd2zu, e, km, kh, l_grid, l_wall, nzb_s_inner, var )
    134     !$acc kernels
    135 
    136 !
    137132!-- Introduce an optional minimum tke
    138133    IF ( e_min > 0.0_wp )  THEN
    139134       !$OMP DO
    140        !$acc loop
    141135       DO  i = nxlg, nxrg
    142136          DO  j = nysg, nyng
    143              !$acc loop vector( 32 )
    144137             DO  k = 1, nzt
    145138                IF ( k > nzb_s_inner(j,i) )  THEN
     
    152145
    153146    !$OMP DO
    154     !$acc loop
    155147    DO  i = nxlg, nxrg
    156148       DO  j = nysg, nyng
    157           !$acc loop vector( 32 )
    158149          DO  k = 1, nzt
    159150
     
    191182                kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / ll ) * km(k,j,i)
    192183
    193 #if ! defined( __openacc )
    194 !
    195 !++             Statistics still have to be realized for accelerators
     184!
    196185!--             Summation for averaged profile (cf. flow_statistics)
    197186                DO  sr = 0, statistic_regions
    198187                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr)
    199188                ENDDO
    200 #endif
     189
    201190             ENDIF
    202191
     
    205194    ENDDO
    206195
    207 #if ! defined( __openacc )
    208 !
    209 !++ Statistics still have to be realized for accelerators
    210196    sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
    211                                                   ! data output
    212 #endif
    213     !$OMP END PARALLEL
     197                                                ! data output
     198!$OMP END PARALLEL
    214199
    215200!
     
    219204!-- values of the diffusivities are not needed
    220205    !$OMP PARALLEL DO
    221     !$acc loop
    222206    DO  i = nxlg, nxrg
    223207       DO  j = nysg, nyng
     
    249233    ENDIF
    250234
    251     !$acc end kernels
    252     !$acc end data
    253 
    254235 END SUBROUTINE diffusivities
  • palm/trunk/SOURCE/exchange_horiz.f90

    r2101 r2118  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! OpenACC directives and related code removed
    2323!
    2424! Former revisions:
     
    8686    USE control_parameters,                                                    &
    8787        ONLY:  bc_lr, bc_lr_cyc, bc_ns, bc_ns_cyc, grid_level,                 &
    88                mg_switch_to_pe0, on_device, synchronous_exchange
     88               mg_switch_to_pe0, synchronous_exchange
    8989               
    9090    USE cpulog,                                                                &
     
    254254!-- with array syntax, explicit loops are used.
    255255    IF ( bc_lr == 'cyclic' )  THEN
    256        IF ( on_device )  THEN
    257           !$acc kernels present( ar )
    258           !$acc loop independent
    259           DO  i = 0, nbgp_local-1
    260              DO  j = nys-nbgp_local, nyn+nbgp_local
    261                 DO  k = nzb, nzt+1
    262                    ar(k,j,nxl-nbgp_local+i) = ar(k,j,nxr-nbgp_local+1+i)
    263                    ar(k,j,nxr+1+i)          = ar(k,j,nxl+i)
    264                 ENDDO
    265              ENDDO
    266           ENDDO
    267           !$acc end kernels
    268        ELSE
    269           ar(:,:,nxl-nbgp_local:nxl-1) = ar(:,:,nxr-nbgp_local+1:nxr)
    270           ar(:,:,nxr+1:nxr+nbgp_local) = ar(:,:,nxl:nxl+nbgp_local-1)
    271        ENDIF
     256       ar(:,:,nxl-nbgp_local:nxl-1) = ar(:,:,nxr-nbgp_local+1:nxr)
     257       ar(:,:,nxr+1:nxr+nbgp_local) = ar(:,:,nxl:nxl+nbgp_local-1)
    272258    ENDIF
    273259
    274260    IF ( bc_ns == 'cyclic' )  THEN
    275        IF ( on_device )  THEN
    276           !$acc kernels present( ar )
    277           DO  i = nxl-nbgp_local, nxr+nbgp_local
    278              !$acc loop independent
    279              DO  j = 0, nbgp_local-1
    280                 !$acc loop independent
    281                 DO  k = nzb, nzt+1
    282                    ar(k,nys-nbgp_local+j,i) = ar(k,nyn-nbgp_local+1+j,i)
    283                    ar(k,nyn+1+j,i)          = ar(k,nys+j,i)
    284                 ENDDO
    285              ENDDO
    286           ENDDO
    287           !$acc end kernels
    288        ELSE
    289           ar(:,nys-nbgp_local:nys-1,:) = ar(:,nyn-nbgp_local+1:nyn,:)
    290           ar(:,nyn+1:nyn+nbgp_local,:) = ar(:,nys:nys+nbgp_local-1,:)
    291        ENDIF
     261       ar(:,nys-nbgp_local:nys-1,:) = ar(:,nyn-nbgp_local+1:nyn,:)
     262       ar(:,nyn+1:nyn+nbgp_local,:) = ar(:,nys:nys+nbgp_local-1,:)
    292263    ENDIF
    293264
  • palm/trunk/SOURCE/fft_xy_mod.f90

    r2101 r2118  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! OpenACC directives and CUDA-fft related code removed
    2323!
    2424! Former revisions:
     
    3131! 1850 2016-04-08 13:29:27Z maronga
    3232! Module renamed
    33 !
    3433!
    3534! 1815 2016-04-06 13:49:59Z raasch
     
    139138        ONLY:  nx, ny, nz
    140139       
    141 #if defined( __cuda_fft )
    142     USE ISO_C_BINDING
    143 #elif defined( __fftw )
     140#if defined( __fftw )
    144141    USE, INTRINSIC ::  ISO_C_BINDING
    145142#endif
     
    192189    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trig_yf  !<
    193190   
    194 #elif defined( __cuda_fft )
    195     INTEGER(C_INT), SAVE ::  plan_xf  !<
    196     INTEGER(C_INT), SAVE ::  plan_xi  !<
    197     INTEGER(C_INT), SAVE ::  plan_yf  !<
    198     INTEGER(C_INT), SAVE ::  plan_yi  !<
    199    
    200     INTEGER(iwp), SAVE   ::  total_points_x_transpo  !<
    201     INTEGER(iwp), SAVE   ::  total_points_y_transpo  !<
    202191#endif
    203192
     
    261250    SUBROUTINE fft_init
    262251
    263        USE cuda_fft_interfaces
    264 
    265252       IMPLICIT NONE
    266253
     
    338325          CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4,      &
    339326                       trig_yb, worky, 0 )
    340 #elif defined( __cuda_fft )
    341           total_points_x_transpo = (nx+1) * (nyn_x-nys_x+1) * (nzt_x-nzb_x+1)
    342           total_points_y_transpo = (ny+1) * (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1)
    343           CALL CUFFTPLAN1D( plan_xf, nx+1, CUFFT_D2Z, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
    344           CALL CUFFTPLAN1D( plan_xi, nx+1, CUFFT_Z2D, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
    345           CALL CUFFTPLAN1D( plan_yf, ny+1, CUFFT_D2Z, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
    346           CALL CUFFTPLAN1D( plan_yi, ny+1, CUFFT_Z2D, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
    347327#else
    348328          message_string = 'no system-specific fft-call available'
     
    403383
    404384
    405        USE cuda_fft_interfaces
    406 #if defined( __cuda_fft )
    407        USE ISO_C_BINDING
    408 #endif
    409 
    410385       IMPLICIT NONE
    411386
     
    429404#elif defined( __nec )
    430405       REAL(wp), DIMENSION(6*(nx+1)) ::  work2  !<
    431 #elif defined( __cuda_fft )
    432        COMPLEX(dp), DIMENSION(0:(nx+1)/2,nys_x:nyn_x,nzb_x:nzt_x) ::           &
    433           ar_tmp  !<
    434        ! following does not work for PGI 14.1 -> to be removed later
    435        ! !$acc declare create( ar_tmp )
    436406#endif
    437407
     
    737707
    738708          ENDIF
    739 
    740 #elif defined( __cuda_fft )
    741 
    742           !$acc data create( ar_tmp )
    743           IF ( forward_fft )  THEN
    744 
    745              !$acc data present( ar )
    746              CALL CUFFTEXECD2Z( plan_xf, ar, ar_tmp )
    747 
    748              !$acc kernels
    749              DO  k = nzb_x, nzt_x
    750                 DO  j = nys_x, nyn_x
    751 
    752                    DO  i = 0, (nx+1)/2
    753                       ar(i,j,k)      = REAL( ar_tmp(i,j,k), KIND=wp )  * dnx
    754                    ENDDO
    755 
    756                    DO  i = 1, (nx+1)/2 - 1
    757                       ar(nx+1-i,j,k) = AIMAG( ar_tmp(i,j,k) ) * dnx
    758                    ENDDO
    759 
    760                 ENDDO
    761              ENDDO
    762              !$acc end kernels
    763              !$acc end data
    764 
    765           ELSE
    766 
    767              !$acc data present( ar )
    768              !$acc kernels
    769              DO  k = nzb_x, nzt_x
    770                 DO  j = nys_x, nyn_x
    771 
    772                    ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0_wp, KIND=wp )
    773 
    774                    DO  i = 1, (nx+1)/2 - 1
    775                       ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k),        &
    776                                              KIND=wp )
    777                    ENDDO
    778                    ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp,     &
    779                                                  KIND=wp )
    780 
    781                 ENDDO
    782              ENDDO
    783              !$acc end kernels
    784 
    785              CALL CUFFTEXECZ2D( plan_xi, ar_tmp, ar )
    786              !$acc end data
    787 
    788           ENDIF
    789           !$acc end data
    790709
    791710#else
     
    1052971
    1053972
    1054        USE cuda_fft_interfaces
    1055 #if defined( __cuda_fft )
    1056        USE ISO_C_BINDING
    1057 #endif
    1058 
    1059973       IMPLICIT NONE
    1060974
     
    1082996#elif defined( __nec )
    1083997       REAL(wp), DIMENSION(6*(ny+1)) ::  work2  !<
    1084 #elif defined( __cuda_fft )
    1085        COMPLEX(dp), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) ::           &
    1086           ar_tmp  !<
    1087        ! following does not work for PGI 14.1 -> to be removed later
    1088        ! !$acc declare create( ar_tmp )
    1089998#endif
    1090999
     
    13641273
    13651274          ENDIF
    1366 #elif defined( __cuda_fft )
    1367 
    1368           !$acc data create( ar_tmp )
    1369           IF ( forward_fft )  THEN
    1370 
    1371              !$acc data present( ar )
    1372              CALL CUFFTEXECD2Z( plan_yf, ar, ar_tmp )
    1373 
    1374              !$acc kernels
    1375              DO  k = nzb_y, nzt_y
    1376                 DO  i = nxl_y, nxr_y
    1377 
    1378                    DO  j = 0, (ny+1)/2
    1379                       ar(j,i,k)      = REAL( ar_tmp(j,i,k), KIND=wp )  * dny
    1380                    ENDDO
    1381 
    1382                    DO  j = 1, (ny+1)/2 - 1
    1383                       ar(ny+1-j,i,k) = AIMAG( ar_tmp(j,i,k) ) * dny
    1384                    ENDDO
    1385 
    1386                 ENDDO
    1387              ENDDO
    1388              !$acc end kernels
    1389              !$acc end data
    1390 
    1391           ELSE
    1392 
    1393              !$acc data present( ar )
    1394              !$acc kernels
    1395              DO  k = nzb_y, nzt_y
    1396                 DO  i = nxl_y, nxr_y
    1397 
    1398                    ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0_wp, KIND=wp )
    1399 
    1400                    DO  j = 1, (ny+1)/2 - 1
    1401                       ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k),        &
    1402                                              KIND=wp )
    1403                    ENDDO
    1404                    ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0_wp,     &
    1405                                                  KIND=wp )
    1406 
    1407                 ENDDO
    1408              ENDDO
    1409              !$acc end kernels
    1410 
    1411              CALL CUFFTEXECZ2D( plan_yi, ar_tmp, ar )
    1412              !$acc end data
    1413 
    1414           ENDIF
    1415           !$acc end data
    1416 
    14171275#else
    14181276          message_string = 'no system-specific fft-call available'
  • palm/trunk/SOURCE/flow_statistics.f90

    r2101 r2118  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! OpenACC version of subroutine removed
    2323!
    2424! Former revisions:
     
    216216!>       are zero at the walls and inside buildings.
    217217!------------------------------------------------------------------------------!
    218 #if ! defined( __openacc )
    219218 SUBROUTINE flow_statistics
    220219 
     
    309308    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
    310309
    311     !$acc update host( km, kh, e, ol, pt, qs, qsws, shf, ts, u, usws, v, vsws, w )
    312310
    313311!
     
    17561754
    17571755 END SUBROUTINE flow_statistics
    1758 
    1759 
    1760 #else
    1761 
    1762 
    1763 !------------------------------------------------------------------------------!
    1764 ! Description:
    1765 ! ------------
    1766 !> flow statistics - accelerator version
    1767 !------------------------------------------------------------------------------!
    1768  SUBROUTINE flow_statistics
    1769 
    1770     USE arrays_3d,                                                             &
    1771         ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
    1772                momentumflux_output_conversion, nr, p, prho, pt, q, qc, ql, qr, &
    1773                qs, qsws, qswst, rho_air, rho_air_zw, rho_ocean, s, sa, saswsb, &
    1774                saswst, shf, ss, ssws, sswst, td_lsa_lpt, td_lsa_q, td_sub_lpt, &
    1775                td_sub_q, time_vert, ts, tswst, u, ug, us, usws, uswst, vsws,   &
    1776                v, vg, vpt, vswst, w, w_subs, waterflux_output_conversion, zw
    1777                  
    1778        
    1779     USE cloud_parameters,                                                      &
    1780         ONLY:  l_d_cp, prr, pt_d_t
    1781        
    1782     USE control_parameters,                                                    &
    1783         ONLY :  average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
    1784                 dt_3d, g, humidity, kappa, large_scale_forcing,                &
    1785                 large_scale_subsidence, max_pr_user, message_string,           &
    1786                 microphysics_seifert, neutral, ocean, passive_scalar,          &
    1787                 simulated_time, use_subsidence_tendencies, use_surface_fluxes, &
    1788                 use_top_fluxes, ws_scheme_mom, ws_scheme_sca
    1789        
    1790     USE cpulog,                                                                &
    1791         ONLY:  cpu_log, log_point
    1792        
    1793     USE grid_variables,                                                        &
    1794         ONLY:  ddx, ddy
    1795        
    1796     USE indices,                                                               &
    1797         ONLY:  ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,       &
    1798                ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner,         &
    1799                nzb_s_inner, nzt, nzt_diff, rflags_invers
    1800        
    1801     USE kinds
    1802    
    1803     USE land_surface_model_mod,                                                &
    1804         ONLY:   ghf_eb, land_surface, m_soil, nzb_soil, nzt_soil,              &
    1805                 qsws_eb, qsws_liq_eb, qsws_soil_eb, qsws_veg_eb, r_a, r_s,     &
    1806                 shf_eb, t_soil
    1807 
    1808     USE netcdf_interface,                                                      &
    1809         ONLY:  dots_rad, dots_soil
    1810 
    1811     USE pegrid
    1812    
    1813     USE radiation_model_mod,                                                   &
    1814         ONLY:  radiation, radiation_scheme, rad_net,                 &
    1815                rad_lw_in, rad_lw_out, rad_sw_in, rad_sw_out
    1816 
    1817 #if defined ( __rrtmg )
    1818     USE radiation_model_mod,                                                   &
    1819         ONLY:  rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rad_lw_cs_hr,   &
    1820                rad_lw_hr,  rad_sw_cs_hr, rad_sw_hr
    1821 #endif
    1822 
    1823     USE statistics
    1824 
    1825     IMPLICIT NONE
    1826 
    1827     INTEGER(iwp) ::  i                   !<
    1828     INTEGER(iwp) ::  j                   !<
    1829     INTEGER(iwp) ::  k                   !<
    1830     INTEGER(iwp) ::  k_surface_level     !<
    1831     INTEGER(iwp) ::  nt                  !<
    1832     INTEGER(iwp) ::  omp_get_thread_num  !<
    1833     INTEGER(iwp) ::  sr                  !<
    1834     INTEGER(iwp) ::  tn                  !<
    1835    
    1836     LOGICAL ::  first  !<
    1837    
    1838     REAL(wp) ::  dptdz_threshold  !<
    1839     REAL(wp) ::  fac              !<
    1840     REAL(wp) ::  height           !<
    1841     REAL(wp) ::  pts              !<
    1842     REAL(wp) ::  sums_l_eper      !<
    1843     REAL(wp) ::  sums_l_etot      !<
    1844     REAL(wp) ::  s1               !<
    1845     REAL(wp) ::  s2               !<
    1846     REAL(wp) ::  s3               !<
    1847     REAL(wp) ::  s4               !<
    1848     REAL(wp) ::  s5               !<
    1849     REAL(wp) ::  s6               !<
    1850     REAL(wp) ::  s7               !<
    1851     REAL(wp) ::  ust              !<
    1852     REAL(wp) ::  ust2             !<
    1853     REAL(wp) ::  u2               !<
    1854     REAL(wp) ::  vst              !<
    1855     REAL(wp) ::  vst2             !<
    1856     REAL(wp) ::  v2               !<
    1857     REAL(wp) ::  w2               !<
    1858     REAL(wp) ::  z_i(2)           !<
    1859 
    1860     REAL(wp) ::  dptdz(nzb+1:nzt+1)    !<
    1861     REAL(wp) ::  sums_ll(nzb:nzt+1,2)  !<
    1862 
    1863     CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
    1864 
    1865 !
    1866 !-- To be on the safe side, check whether flow_statistics has already been
    1867 !-- called once after the current time step
    1868     IF ( flow_statistics_called )  THEN
    1869 
    1870        message_string = 'flow_statistics is called two times within one ' // &
    1871                         'timestep'
    1872        CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
    1873 
    1874     ENDIF
    1875 
    1876     !$acc data create( sums, sums_l )
    1877     !$acc update device( hom )
    1878 
    1879 !
    1880 !-- Compute statistics for each (sub-)region
    1881     DO  sr = 0, statistic_regions
    1882 
    1883 !
    1884 !--    Initialize (local) summation array
    1885        sums_l = 0.0_wp
    1886 
    1887 !
    1888 !--    Store sums that have been computed in other subroutines in summation
    1889 !--    array
    1890        sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
    1891 !--    WARNING: next line still has to be adjusted for OpenMP
    1892        sums_l(:,21,0) = sums_wsts_bc_l(:,sr) *                                 &
    1893                         heatflux_output_conversion  ! heat flux from advec_s_bc
    1894        sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
    1895        sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
    1896 
    1897 !
    1898 !--    When calcuating horizontally-averaged total (resolved- plus subgrid-
    1899 !--    scale) vertical fluxes and velocity variances by using commonly-
    1900 !--    applied Reynolds-based methods ( e.g. <w'pt'> = (w-<w>)*(pt-<pt>) )
    1901 !--    in combination with the 5th order advection scheme, pronounced
    1902 !--    artificial kinks could be observed in the vertical profiles near the
    1903 !--    surface. Please note: these kinks were not related to the model truth,
    1904 !--    i.e. these kinks are just related to an evaluation problem.   
    1905 !--    In order avoid these kinks, vertical fluxes and horizontal as well
    1906 !--    vertical velocity variances are calculated directly within the advection
    1907 !--    routines, according to the numerical discretization, to evaluate the
    1908 !--    statistical quantities as they will appear within the prognostic
    1909 !--    equations.
    1910 !--    Copy the turbulent quantities, evaluated in the advection routines to
    1911 !--    the local array sums_l() for further computations.
    1912        IF ( ws_scheme_mom .AND. sr == 0 )  THEN
    1913 
    1914 !
    1915 !--       According to the Neumann bc for the horizontal velocity components,
    1916 !--       the corresponding fluxes has to satisfiy the same bc.
    1917           IF ( ocean )  THEN
    1918              sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
    1919              sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
    1920           ENDIF
    1921 
    1922           DO  i = 0, threads_per_task-1
    1923 !
    1924 !--          Swap the turbulent quantities evaluated in advec_ws.
    1925              sums_l(:,13,i) = sums_wsus_ws_l(:,i)                              &
    1926                               * momentumflux_output_conversion ! w*u*
    1927              sums_l(:,15,i) = sums_wsvs_ws_l(:,i)                              &
    1928                               * momentumflux_output_conversion ! w*v*
    1929              sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
    1930              sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
    1931              sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
    1932              sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        &
    1933                               ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
    1934                                 sums_ws2_ws_l(:,i) )    ! e*
    1935              DO  k = nzb, nzt
    1936                 sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5_wp * ( &
    1937                                                       sums_us2_ws_l(k,i) +     &
    1938                                                       sums_vs2_ws_l(k,i) +     &
    1939                                                       sums_ws2_ws_l(k,i)     )
    1940              ENDDO
    1941           ENDDO
    1942 
    1943        ENDIF
    1944 
    1945        IF ( ws_scheme_sca .AND. sr == 0 )  THEN
    1946 
    1947           DO  i = 0, threads_per_task-1
    1948              sums_l(:,17,i) = sums_wspts_ws_l(:,i)                             &
    1949                               * heatflux_output_conversion        ! w*pt* from advec_s_ws
    1950              IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
    1951              IF ( humidity       )  sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)      &
    1952                                             * waterflux_output_conversion !w*q*
    1953              IF ( passive_scalar )  sums_l(:,116,i) = sums_wsss_ws_l(:,i) !w*s*
    1954           ENDDO
    1955 
    1956        ENDIF
    1957 !
    1958 !--    Horizontally averaged profiles of horizontal velocities and temperature.
    1959 !--    They must have been computed before, because they are already required
    1960 !--    for other horizontal averages.
    1961        tn = 0
    1962 
    1963        !$OMP PARALLEL PRIVATE( i, j, k, tn )
    1964 !$     tn = omp_get_thread_num()
    1965 
    1966        !$acc update device( sums_l )
    1967 
    1968        !$OMP DO
    1969        !$acc parallel loop gang present( pt, rflags_invers, rmask, sums_l, u, v ) create( s1, s2, s3 )
    1970        DO  k = nzb, nzt+1
    1971           s1 = 0
    1972           s2 = 0
    1973           s3 = 0
    1974           !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
    1975           DO  i = nxl, nxr
    1976              DO  j =  nys, nyn
    1977 !
    1978 !--             k+1 is used in rflags since rflags is set 0 at surface points
    1979                 s1 = s1 + u(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    1980                 s2 = s2 + v(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    1981                 s3 = s3 + pt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    1982              ENDDO
    1983           ENDDO
    1984           sums_l(k,1,tn) = s1
    1985           sums_l(k,2,tn) = s2
    1986           sums_l(k,4,tn) = s3
    1987        ENDDO
    1988        !$acc end parallel loop
    1989 
    1990 !
    1991 !--    Horizontally averaged profile of salinity
    1992        IF ( ocean )  THEN
    1993           !$OMP DO
    1994           !$acc parallel loop gang present( rflags_invers, rmask, sums_l, sa ) create( s1 )
    1995           DO  k = nzb, nzt+1
    1996              s1 = 0
    1997              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    1998              DO  i = nxl, nxr
    1999                 DO  j =  nys, nyn
    2000                    s1 = s1 + sa(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2001                 ENDDO
    2002              ENDDO
    2003              sums_l(k,23,tn) = s1
    2004           ENDDO
    2005           !$acc end parallel loop
    2006        ENDIF
    2007 
    2008 !
    2009 !--    Horizontally averaged profiles of virtual potential temperature,
    2010 !--    total water content, specific humidity and liquid water potential
    2011 !--    temperature
    2012        IF ( humidity )  THEN
    2013 
    2014           !$OMP DO
    2015           !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l, vpt ) create( s1, s2 )
    2016           DO  k = nzb, nzt+1
    2017              s1 = 0
    2018              s2 = 0
    2019              !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
    2020              DO  i = nxl, nxr
    2021                 DO  j =  nys, nyn
    2022                    s1 = s1 + q(k,j,i)   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2023                    s2 = s2 + vpt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2024                 ENDDO
    2025              ENDDO
    2026              sums_l(k,41,tn) = s1
    2027              sums_l(k,44,tn) = s2
    2028           ENDDO
    2029           !$acc end parallel loop
    2030 
    2031           IF ( cloud_physics )  THEN
    2032              !$OMP DO
    2033              !$acc parallel loop gang present( pt, q, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
    2034              DO  k = nzb, nzt+1
    2035                 s1 = 0
    2036                 s2 = 0
    2037                 !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
    2038                 DO  i = nxl, nxr
    2039                    DO  j =  nys, nyn
    2040                       s1 = s1 + ( q(k,j,i) - ql(k,j,i) ) * &
    2041                                 rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2042                       s2 = s2 + ( pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) ) * &
    2043                                 rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2044                    ENDDO
    2045                 ENDDO
    2046                 sums_l(k,42,tn) = s1
    2047                 sums_l(k,43,tn) = s2
    2048              ENDDO
    2049              !$acc end parallel loop
    2050           ENDIF
    2051        ENDIF
    2052 
    2053 !
    2054 !--    Horizontally averaged profiles of passive scalar
    2055        IF ( passive_scalar )  THEN
    2056           !$OMP DO
    2057           !$acc parallel loop gang present( s, rflags_invers, rmask, sums_l ) create( s1 )
    2058           DO  k = nzb, nzt+1
    2059              s1 = 0
    2060              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2061              DO  i = nxl, nxr
    2062                 DO  j =  nys, nyn
    2063                    s1 = s1 + s(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2064                 ENDDO
    2065              ENDDO
    2066              sums_l(k,117,tn) = s1
    2067           ENDDO
    2068           !$acc end parallel loop
    2069        ENDIF
    2070        !$OMP END PARALLEL
    2071 
    2072 !
    2073 !--    Summation of thread sums
    2074        IF ( threads_per_task > 1 )  THEN
    2075           DO  i = 1, threads_per_task-1
    2076              !$acc parallel present( sums_l )
    2077              sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
    2078              sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
    2079              sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
    2080              !$acc end parallel
    2081              IF ( ocean )  THEN
    2082                 !$acc parallel present( sums_l )
    2083                 sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
    2084                 !$acc end parallel
    2085              ENDIF
    2086              IF ( humidity )  THEN
    2087                 !$acc parallel present( sums_l )
    2088                 sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
    2089                 sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
    2090                 !$acc end parallel
    2091                 IF ( cloud_physics )  THEN
    2092                    !$acc parallel present( sums_l )
    2093                    sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
    2094                    sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
    2095                    !$acc end parallel
    2096                 ENDIF
    2097              ENDIF
    2098              IF ( passive_scalar )  THEN
    2099                 !$acc parallel present( sums_l )
    2100                 sums_l(:,117,0) = sums_l(:,117,0) + sums_l(:,117,i)
    2101                 !$acc end parallel
    2102              ENDIF
    2103           ENDDO
    2104        ENDIF
    2105 
    2106 #if defined( __parallel )
    2107 !
    2108 !--    Compute total sum from local sums
    2109        !$acc update host( sums_l )
    2110        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2111        CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL,  &
    2112                            MPI_SUM, comm2d, ierr )
    2113        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2114        CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL,  &
    2115                            MPI_SUM, comm2d, ierr )
    2116        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2117        CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL,  &
    2118                            MPI_SUM, comm2d, ierr )
    2119        IF ( ocean )  THEN
    2120           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2121           CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb,       &
    2122                               MPI_REAL, MPI_SUM, comm2d, ierr )
    2123        ENDIF
    2124        IF ( humidity ) THEN
    2125           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2126           CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb,       &
    2127                               MPI_REAL, MPI_SUM, comm2d, ierr )
    2128           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2129           CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
    2130                               MPI_REAL, MPI_SUM, comm2d, ierr )
    2131           IF ( cloud_physics ) THEN
    2132              IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2133              CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb,    &
    2134                                  MPI_REAL, MPI_SUM, comm2d, ierr )
    2135              IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2136              CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb,    &
    2137                                  MPI_REAL, MPI_SUM, comm2d, ierr )
    2138           ENDIF
    2139        ENDIF
    2140 
    2141        IF ( passive_scalar )  THEN
    2142           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2143           CALL MPI_ALLREDUCE( sums_l(nzb,117,0), sums(nzb,117), nzt+2-nzb,     &
    2144                               MPI_REAL, MPI_SUM, comm2d, ierr )
    2145        ENDIF
    2146        !$acc update device( sums )
    2147 #else
    2148        !$acc parallel present( sums, sums_l )
    2149        sums(:,1) = sums_l(:,1,0)
    2150        sums(:,2) = sums_l(:,2,0)
    2151        sums(:,4) = sums_l(:,4,0)
    2152        !$acc end parallel
    2153        IF ( ocean )  THEN
    2154           !$acc parallel present( sums, sums_l )
    2155           sums(:,23) = sums_l(:,23,0)
    2156           !$acc end parallel
    2157        ENDIF
    2158        IF ( humidity )  THEN
    2159           !$acc parallel present( sums, sums_l )
    2160           sums(:,44) = sums_l(:,44,0)
    2161           sums(:,41) = sums_l(:,41,0)
    2162           !$acc end parallel
    2163           IF ( cloud_physics )  THEN
    2164              !$acc parallel present( sums, sums_l )
    2165              sums(:,42) = sums_l(:,42,0)
    2166              sums(:,43) = sums_l(:,43,0)
    2167              !$acc end parallel
    2168           ENDIF
    2169        ENDIF
    2170        IF ( passive_scalar )  THEN
    2171           !$acc parallel present( sums, sums_l )
    2172           sums(:,117) = sums_l(:,117,0)
    2173           !$acc end parallel
    2174        ENDIF
    2175 #endif
    2176 
    2177 !
    2178 !--    Final values are obtained by division by the total number of grid points
    2179 !--    used for summation. After that store profiles.
    2180        !$acc parallel present( hom, ngp_2dh, ngp_2dh_s_inner, sums )
    2181        sums(:,1) = sums(:,1) / ngp_2dh(sr)
    2182        sums(:,2) = sums(:,2) / ngp_2dh(sr)
    2183        sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
    2184        hom(:,1,1,sr) = sums(:,1)             ! u
    2185        hom(:,1,2,sr) = sums(:,2)             ! v
    2186        hom(:,1,4,sr) = sums(:,4)             ! pt
    2187        !$acc end parallel
    2188 
    2189 !
    2190 !--    Salinity
    2191        IF ( ocean )  THEN
    2192           !$acc parallel present( hom, ngp_2dh_s_inner, sums )
    2193           sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
    2194           hom(:,1,23,sr) = sums(:,23)             ! sa
    2195           !$acc end parallel
    2196        ENDIF
    2197 
    2198 !
    2199 !--    Humidity and cloud parameters
    2200        IF ( humidity ) THEN
    2201           !$acc parallel present( hom, ngp_2dh_s_inner, sums )
    2202           sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
    2203           sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
    2204           hom(:,1,44,sr) = sums(:,44)                ! vpt
    2205           hom(:,1,41,sr) = sums(:,41)                ! qv (q)
    2206           !$acc end parallel
    2207           IF ( cloud_physics ) THEN
    2208              !$acc parallel present( hom, ngp_2dh_s_inner, sums )
    2209              sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
    2210              sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
    2211              hom(:,1,42,sr) = sums(:,42)             ! qv
    2212              hom(:,1,43,sr) = sums(:,43)             ! pt
    2213              !$acc end parallel
    2214           ENDIF
    2215        ENDIF
    2216 
    2217 !
    2218 !--    Passive scalar
    2219        IF ( passive_scalar )  THEN
    2220           !$acc parallel present( hom, ngp_2dh_s_inner, sums )
    2221           sums(:,117)     = sums(:,117) / ngp_2dh_s_inner(:,sr)
    2222           hom(:,1,117,sr) = sums(:,117)                ! s
    2223           !$acc end parallel
    2224        ENDIF
    2225 
    2226 !
    2227 !--    Horizontally averaged profiles of the remaining prognostic variables,
    2228 !--    variances, the total and the perturbation energy (single values in last
    2229 !--    column of sums_l) and some diagnostic quantities.
    2230 !--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
    2231 !--    ----  speaking the following k-loop would have to be split up and
    2232 !--          rearranged according to the staggered grid.
    2233 !--          However, this implies no error since staggered velocity components
    2234 !--          are zero at the walls and inside buildings.
    2235        tn = 0
    2236        !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper,             &
    2237        !$OMP                   sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2,  &
    2238        !$OMP                   w2 )
    2239 !$     tn = omp_get_thread_num()
    2240 
    2241        !$OMP DO
    2242        !$acc parallel loop gang present( e, hom, kh, km, p, pt, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, s5, s6, s7 )
    2243        DO  k = nzb, nzt+1
    2244           s1 = 0
    2245           s2 = 0
    2246           s3 = 0
    2247           s4 = 0
    2248           s5 = 0
    2249           s6 = 0
    2250           s7 = 0
    2251           !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5, s6, s7 )
    2252           DO  i = nxl, nxr
    2253              DO  j =  nys, nyn
    2254 !
    2255 !--             Prognostic and diagnostic variables
    2256                 s1 = s1 + w(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2257                 s2 = s2 + e(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2258                 s3 = s3 + km(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2259                 s4 = s4 + kh(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2260                 s5 = s5 + p(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2261                 s6 = s6 + ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr) * &
    2262                           rflags_invers(j,i,k+1)
    2263 !
    2264 !--             Higher moments
    2265 !--             (Computation of the skewness of w further below)
    2266                 s7 = s7 + w(k,j,i)**3 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2267              ENDDO
    2268           ENDDO
    2269           sums_l(k,3,tn)  = s1
    2270           sums_l(k,8,tn)  = s2
    2271           sums_l(k,9,tn)  = s3
    2272           sums_l(k,10,tn) = s4
    2273           sums_l(k,40,tn) = s5
    2274           sums_l(k,33,tn) = s6
    2275           sums_l(k,38,tn) = s7
    2276        ENDDO
    2277        !$acc end parallel loop
    2278 
    2279        IF ( humidity )  THEN
    2280           !$OMP DO
    2281           !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l ) create( s1 )
    2282           DO  k = nzb, nzt+1
    2283              s1 = 0
    2284              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2285              DO  i = nxl, nxr
    2286                 DO  j =  nys, nyn
    2287                    s1 = s1 + ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr) * &
    2288                              rflags_invers(j,i,k+1)
    2289                 ENDDO
    2290              ENDDO
    2291              sums_l(k,70,tn) = s1
    2292           ENDDO
    2293           !$acc end parallel loop
    2294        ENDIF
    2295 
    2296 !
    2297 !--    Total and perturbation energy for the total domain (being
    2298 !--    collected in the last column of sums_l).
    2299        s1 = 0
    2300        !$OMP DO
    2301        !$acc parallel loop collapse(3) present( rflags_invers, rmask, u, v, w ) reduction(+:s1)
    2302        DO  i = nxl, nxr
    2303           DO  j =  nys, nyn
    2304              DO  k = nzb, nzt+1
    2305                 s1 = s1 + 0.5_wp *                                             &
    2306                           ( u(k,j,i)**2 + v(k,j,i)**2 + w(k,j,i)**2 ) *        &
    2307                           rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2308              ENDDO
    2309           ENDDO
    2310        ENDDO
    2311        !$acc end parallel loop
    2312        !$acc parallel present( sums_l )
    2313        sums_l(nzb+4,pr_palm,tn) = s1
    2314        !$acc end parallel
    2315 
    2316        !$OMP DO
    2317        !$acc parallel present( rmask, sums_l, us, usws, vsws, ts ) create( s1, s2, s3, s4 )
    2318        s1 = 0
    2319        s2 = 0
    2320        s3 = 0
    2321        s4 = 0
    2322        !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
    2323        DO  i = nxl, nxr
    2324           DO  j =  nys, nyn
    2325 !
    2326 !--          2D-arrays (being collected in the last column of sums_l)
    2327              s1 = s1 + us(j,i)   * rmask(j,i,sr)
    2328              s2 = s2 + usws(j,i) * rmask(j,i,sr)
    2329              s3 = s3 + vsws(j,i) * rmask(j,i,sr)
    2330              s4 = s4 + ts(j,i)   * rmask(j,i,sr)
    2331           ENDDO
    2332        ENDDO
    2333        sums_l(nzb,pr_palm,tn)   = s1
    2334        sums_l(nzb+1,pr_palm,tn) = s2
    2335        sums_l(nzb+2,pr_palm,tn) = s3
    2336        sums_l(nzb+3,pr_palm,tn) = s4
    2337        !$acc end parallel
    2338 
    2339        IF ( humidity )  THEN
    2340           !$acc parallel present( qs, rmask, sums_l ) create( s1 )
    2341           s1 = 0
    2342           !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2343           DO  i = nxl, nxr
    2344              DO  j =  nys, nyn
    2345                 s1 = s1 + qs(j,i) * rmask(j,i,sr)
    2346              ENDDO
    2347           ENDDO
    2348           sums_l(nzb+12,pr_palm,tn) = s1
    2349           !$acc end parallel
    2350        ENDIF
    2351 
    2352        IF ( passive_scalar )  THEN
    2353           !$acc parallel present( ss, rmask, sums_l ) create( s1 )
    2354           s1 = 0
    2355           !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2356           DO  i = nxl, nxr
    2357              DO  j =  nys, nyn
    2358                 s1 = s1 + ss(j,i) * rmask(j,i,sr)
    2359              ENDDO
    2360           ENDDO
    2361           sums_l(nzb+13,pr_palm,tn) = s1
    2362           !$acc end parallel
    2363        ENDIF
    2364 
    2365 !
    2366 !--    Computation of statistics when ws-scheme is not used. Else these
    2367 !--    quantities are evaluated in the advection routines.
    2368        IF ( .NOT. ws_scheme_mom .OR. sr /= 0 .OR. simulated_time == 0.0_wp )   &
    2369        THEN
    2370 
    2371           !$OMP DO
    2372           !$acc parallel loop gang present( u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, ust2, vst2, w2 )
    2373           DO  k = nzb, nzt+1
    2374              s1 = 0
    2375              s2 = 0
    2376              s3 = 0
    2377              s4 = 0
    2378              !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
    2379              DO  i = nxl, nxr
    2380                 DO  j =  nys, nyn
    2381                    ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
    2382                    vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
    2383                    w2   = w(k,j,i)**2
    2384 
    2385                    s1 = s1 + ust2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2386                    s2 = s2 + vst2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2387                    s3 = s3 + w2   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2388 !
    2389 !--                Perturbation energy
    2390                    s4 = s4 + 0.5_wp * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) *   &
    2391                              rflags_invers(j,i,k+1)
    2392                 ENDDO
    2393              ENDDO
    2394              sums_l(k,30,tn) = s1
    2395              sums_l(k,31,tn) = s2
    2396              sums_l(k,32,tn) = s3
    2397              sums_l(k,34,tn) = s4
    2398           ENDDO
    2399           !$acc end parallel loop
    2400 !
    2401 !--       Total perturbation TKE
    2402           !$OMP DO
    2403           !$acc parallel present( sums_l ) create( s1 )
    2404           s1 = 0
    2405           !$acc loop reduction( +: s1 )
    2406           DO  k = nzb, nzt+1
    2407              s1 = s1 + sums_l(k,34,tn)
    2408           ENDDO
    2409           sums_l(nzb+5,pr_palm,tn) = s1
    2410           !$acc end parallel
    2411 
    2412        ENDIF
    2413 
    2414 !
    2415 !--    Horizontally averaged profiles of the vertical fluxes
    2416 
    2417 !
    2418 !--    Subgridscale fluxes.
    2419 !--    WARNING: If a Prandtl-layer is used (k=nzb for flat terrain), the fluxes
    2420 !--    -------  should be calculated there in a different way. This is done
    2421 !--             in the next loop further below, where results from this loop are
    2422 !--             overwritten. However, THIS WORKS IN CASE OF FLAT TERRAIN ONLY!
    2423 !--             The non-flat case still has to be handled.
    2424 !--    NOTE: for simplicity, nzb_s_inner is used below, although
    2425 !--    ----  strictly speaking the following k-loop would have to be
    2426 !--          split up according to the staggered grid.
    2427 !--          However, this implies no error since staggered velocity
    2428 !--          components are zero at the walls and inside buildings.
    2429        !$OMP DO
    2430        !$acc parallel loop gang present( ddzu, kh, km, pt, u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
    2431        DO  k = nzb, nzt_diff
    2432           s1 = 0
    2433           s2 = 0
    2434           s3 = 0
    2435           !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
    2436           DO  i = nxl, nxr
    2437              DO  j = nys, nyn
    2438 
    2439 !
    2440 !--             Momentum flux w"u"
    2441                 s1 = s1 - 0.25_wp * (                                          &
    2442                                km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
    2443                                                            ) * (               &
    2444                                    ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
    2445                                  + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
    2446                                                                )               &
    2447                                * rmask(j,i,sr) * rflags_invers(j,i,k+1)        &
    2448                                * rho_air_zw(k)                                 &
    2449                                * momentumflux_output_conversion(k)
    2450 !
    2451 !--             Momentum flux w"v"
    2452                 s2 = s2 - 0.25_wp * (                                          &
    2453                                km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
    2454                                                            ) * (               &
    2455                                    ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
    2456                                  + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
    2457                                                                )               &
    2458                                * rmask(j,i,sr) * rflags_invers(j,i,k+1)        &
    2459                                * rho_air_zw(k)                                 &
    2460                                * momentumflux_output_conversion(k)
    2461 !
    2462 !--             Heat flux w"pt"
    2463                 s3 = s3 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                 &
    2464                                  * ( pt(k+1,j,i) - pt(k,j,i) )                 &
    2465                                  * rho_air_zw(k)                               &
    2466                                  * heatflux_output_conversion(k)               &
    2467                                  * ddzu(k+1) * rmask(j,i,sr)                   &
    2468                                  * rflags_invers(j,i,k+1)
    2469              ENDDO
    2470           ENDDO
    2471           sums_l(k,12,tn) = s1
    2472           sums_l(k,14,tn) = s2
    2473           sums_l(k,16,tn) = s3
    2474        ENDDO
    2475        !$acc end parallel loop
    2476 
    2477 !
    2478 !--    Salinity flux w"sa"
    2479        IF ( ocean )  THEN
    2480           !$acc parallel loop gang present( ddzu, kh, sa, rflags_invers, rmask, sums_l ) create( s1 )
    2481           DO  k = nzb, nzt_diff
    2482              s1 = 0
    2483              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2484              DO  i = nxl, nxr
    2485                 DO  j = nys, nyn
    2486                    s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
    2487                                     * ( sa(k+1,j,i) - sa(k,j,i) )              &
    2488                                     * ddzu(k+1) * rmask(j,i,sr)                &
    2489                                     * rflags_invers(j,i,k+1)
    2490                 ENDDO
    2491              ENDDO
    2492              sums_l(k,65,tn) = s1
    2493           ENDDO
    2494           !$acc end parallel loop
    2495        ENDIF
    2496 
    2497 !
    2498 !--    Buoyancy flux, water flux (humidity flux) w"q"
    2499        IF ( humidity ) THEN
    2500 
    2501           !$acc parallel loop gang present( ddzu, kh, q, vpt, rflags_invers, rmask, sums_l ) create( s1, s2 )
    2502           DO  k = nzb, nzt_diff
    2503              s1 = 0
    2504              s2 = 0
    2505              !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
    2506              DO  i = nxl, nxr
    2507                 DO  j = nys, nyn
    2508                    s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
    2509                                     * ( vpt(k+1,j,i) - vpt(k,j,i) )            &
    2510                                     * rho_air_zw(k)                            &
    2511                                     * heatflux_output_conversion(k)            &
    2512                                     * ddzu(k+1) * rmask(j,i,sr)                &
    2513                                     * rflags_invers(j,i,k+1)
    2514                    s2 = s2 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
    2515                                     * ( q(k+1,j,i) - q(k,j,i) )                &
    2516                                     * rho_air_zw(k)                            &
    2517                                     * waterflux_output_conversion(k)           &
    2518                                     * ddzu(k+1) * rmask(j,i,sr)                &
    2519                                     * rflags_invers(j,i,k+1)
    2520                 ENDDO
    2521              ENDDO
    2522              sums_l(k,45,tn) = s1
    2523              sums_l(k,48,tn) = s2
    2524           ENDDO
    2525           !$acc end parallel loop
    2526 
    2527           IF ( cloud_physics ) THEN
    2528 
    2529              !$acc parallel loop gang present( ddzu, kh, q, ql, rflags_invers, rmask, sums_l ) create( s1 )
    2530              DO  k = nzb, nzt_diff
    2531                 s1 = 0
    2532                 !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2533                 DO  i = nxl, nxr
    2534                    DO  j = nys, nyn
    2535                       s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )           &
    2536                                        *  ( ( q(k+1,j,i) - ql(k+1,j,i) )       &
    2537                                           - ( q(k,j,i) - ql(k,j,i) ) )         &
    2538                                        * rho_air_zw(k)                         &
    2539                                        * waterflux_output_conversion(k)        &
    2540                                        * ddzu(k+1) * rmask(j,i,sr)             &
    2541                                        * rflags_invers(j,i,k+1)
    2542                    ENDDO
    2543                 ENDDO
    2544                 sums_l(k,51,tn) = s1
    2545              ENDDO
    2546              !$acc end parallel loop
    2547 
    2548           ENDIF
    2549 
    2550        ENDIF
    2551 !
    2552 !--    Passive scalar flux
    2553        IF ( passive_scalar )  THEN
    2554 
    2555           !$acc parallel loop gang present( ddzu, kh, s, rflags_invers, rmask, sums_l ) create( s1 )
    2556           DO  k = nzb, nzt_diff
    2557              s1 = 0
    2558              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2559              DO  i = nxl, nxr
    2560                 DO  j = nys, nyn
    2561                    s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
    2562                                     * ( s(k+1,j,i) - s(k,j,i) )                &
    2563                                     * ddzu(k+1) * rmask(j,i,sr)                &
    2564                                     * rflags_invers(j,i,k+1)
    2565                 ENDDO
    2566              ENDDO
    2567              sums_l(k,119,tn) = s1
    2568           ENDDO
    2569           !$acc end parallel loop
    2570 
    2571        ENDIF
    2572 
    2573        IF ( use_surface_fluxes )  THEN
    2574 
    2575           !$OMP DO
    2576           !$acc parallel present( rmask, shf, sums_l, usws, vsws ) create( s1, s2, s3, s4, s5 )
    2577           s1 = 0
    2578           s2 = 0
    2579           s3 = 0
    2580           s4 = 0
    2581           s5 = 0
    2582           !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
    2583           DO  i = nxl, nxr
    2584              DO  j =  nys, nyn
    2585 !
    2586 !--             Subgridscale fluxes in the Prandtl layer
    2587                 s1 = s1 + usws(j,i) * momentumflux_output_conversion(nzb)      &
    2588                                     * rmask(j,i,sr) ! w"u"
    2589                 s2 = s2 + vsws(j,i) * momentumflux_output_conversion(nzb)      &
    2590                                     * rmask(j,i,sr) ! w"v"
    2591                 s3 = s3 + shf(j,i)  * heatflux_output_conversion(nzb)          &
    2592                                     * rmask(j,i,sr) ! w"pt"
    2593                 s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
    2594                 s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
    2595              ENDDO
    2596           ENDDO
    2597           sums_l(nzb,12,tn) = s1
    2598           sums_l(nzb,14,tn) = s2
    2599           sums_l(nzb,16,tn) = s3
    2600           sums_l(nzb,58,tn) = s4
    2601           sums_l(nzb,61,tn) = s5
    2602           !$acc end parallel
    2603 
    2604           IF ( ocean )  THEN
    2605 
    2606              !$OMP DO
    2607              !$acc parallel present( rmask, saswsb, sums_l ) create( s1 )
    2608              s1 = 0
    2609              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2610              DO  i = nxl, nxr
    2611                 DO  j =  nys, nyn
    2612                    s1 = s1 + saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
    2613                 ENDDO
    2614              ENDDO
    2615              sums_l(nzb,65,tn) = s1
    2616              !$acc end parallel
    2617 
    2618           ENDIF
    2619 
    2620           IF ( humidity )  THEN
    2621 
    2622              !$OMP DO
    2623              !$acc parallel present( pt, q, qsws, rmask, shf, sums_l ) create( s1, s2 )
    2624              s1 = 0
    2625              s2 = 0
    2626              !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
    2627              DO  i = nxl, nxr
    2628                 DO  j =  nys, nyn
    2629                    s1 = s1 + qsws(j,i) * waterflux_output_conversion(nzb)      &
    2630                                        * rmask(j,i,sr) ! w"q" (w"qv")
    2631                    s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * shf(j,i)    &
    2632                                + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )           &
    2633                              * heatflux_output_conversion(nzb)
    2634                 ENDDO
    2635              ENDDO
    2636              sums_l(nzb,48,tn) = s1
    2637              sums_l(nzb,45,tn) = s2
    2638              !$acc end parallel
    2639 
    2640              IF ( cloud_droplets )  THEN
    2641 
    2642                 !$OMP DO
    2643                 !$acc parallel present( pt, q, ql, qsws, rmask, shf, sums_l ) create( s1 )
    2644                 s1 = 0
    2645                 !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2646                 DO  i = nxl, nxr
    2647                    DO  j =  nys, nyn
    2648                       s1 = s1 + ( ( 1.0_wp +                                   &
    2649                                     0.61_wp * q(nzb,j,i) - ql(nzb,j,i) ) *     &
    2650                                  shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )&
    2651                                 * heatflux_output_conversion(nzb)
    2652                    ENDDO
    2653                 ENDDO
    2654                 sums_l(nzb,45,tn) = s1
    2655                 !$acc end parallel
    2656 
    2657              ENDIF
    2658 
    2659              IF ( cloud_physics )  THEN
    2660 
    2661                 !$OMP DO
    2662                 !$acc parallel present( qsws, rmask, sums_l ) create( s1 )
    2663                 s1 = 0
    2664                 !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2665                 DO  i = nxl, nxr
    2666                    DO  j =  nys, nyn
    2667 !
    2668 !--                   Formula does not work if ql(nzb) /= 0.0
    2669                       s1 = s1 + qsws(j,i) * waterflux_output_conversion(nzb)   &
    2670                                           * rmask(j,i,sr)   ! w"q" (w"qv")
    2671                    ENDDO
    2672                 ENDDO
    2673                 sums_l(nzb,51,tn) = s1
    2674                 !$acc end parallel
    2675 
    2676              ENDIF
    2677 
    2678           ENDIF
    2679 
    2680           IF ( passive_scalar )  THEN
    2681 
    2682              !$OMP DO
    2683              !$acc parallel present( ssws, rmask, sums_l ) create( s1 )
    2684              s1 = 0
    2685              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2686              DO  i = nxl, nxr
    2687                 DO  j =  nys, nyn
    2688                    s1 = s1 + ssws(j,i) * rmask(j,i,sr)  ! w"s"
    2689                 ENDDO
    2690              ENDDO
    2691              sums_l(nzb,119,tn) = s1
    2692              !$acc end parallel
    2693 
    2694           ENDIF
    2695 
    2696        ENDIF
    2697 
    2698 !
    2699 !--    Subgridscale fluxes at the top surface
    2700        IF ( use_top_fluxes )  THEN
    2701 
    2702           !$OMP DO
    2703           !$acc parallel present( rmask, sums_l, tswst, uswst, vswst ) create( s1, s2, s3, s4, s5 )
    2704           s1 = 0
    2705           s2 = 0
    2706           s3 = 0
    2707           s4 = 0
    2708           s5 = 0
    2709           !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
    2710           DO  i = nxl, nxr
    2711              DO  j =  nys, nyn
    2712                 s1 = s1 + uswst(j,i) * momentumflux_output_conversion(nzt:nzt+1) &
    2713                                      * rmask(j,i,sr)    ! w"u"
    2714                 s2 = s2 + vswst(j,i) * momentumflux_output_conversion(nzt:nzt+1) &
    2715                                      * rmask(j,i,sr)    ! w"v"
    2716                 s3 = s3 + tswst(j,i) * heatflux_output_conversion(nzt:nzt+1)   &
    2717                                      * rmask(j,i,sr)    ! w"pt"
    2718                 s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
    2719                 s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
    2720              ENDDO
    2721           ENDDO
    2722           sums_l(nzt:nzt+1,12,tn) = s1
    2723           sums_l(nzt:nzt+1,14,tn) = s2
    2724           sums_l(nzt:nzt+1,16,tn) = s3
    2725           sums_l(nzt:nzt+1,58,tn) = s4
    2726           sums_l(nzt:nzt+1,61,tn) = s5
    2727           !$acc end parallel
    2728 
    2729           IF ( ocean )  THEN
    2730 
    2731              !$OMP DO
    2732              !$acc parallel present( rmask, saswst, sums_l ) create( s1 )
    2733              s1 = 0
    2734              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2735              DO  i = nxl, nxr
    2736                 DO  j =  nys, nyn
    2737                    s1 = s1 + saswst(j,i) * rmask(j,i,sr)  ! w"sa"
    2738                 ENDDO
    2739              ENDDO
    2740              sums_l(nzt,65,tn) = s1
    2741              !$acc end parallel
    2742 
    2743           ENDIF
    2744 
    2745           IF ( humidity )  THEN
    2746 
    2747              !$OMP DO
    2748              !$acc parallel present( pt, q, qswst, rmask, tswst, sums_l ) create( s1, s2 )
    2749              s1 = 0
    2750              s2 = 0
    2751              !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
    2752              DO  i = nxl, nxr
    2753                 DO  j =  nys, nyn
    2754                    s1 = s1 + qswst(j,i) * waterflux_output_conversion(nzt)     &
    2755                                         * rmask(j,i,sr) ! w"q" (w"qv")
    2756                    s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * tswst(j,i) +&
    2757                                  0.61_wp * pt(nzt,j,i) * qswst(j,i) )          &
    2758                              * heatflux_output_conversion(nzt)
    2759                 ENDDO
    2760              ENDDO
    2761              sums_l(nzt,48,tn) = s1
    2762              sums_l(nzt,45,tn) = s2
    2763              !$acc end parallel
    2764 
    2765              IF ( cloud_droplets )  THEN
    2766 
    2767                 !$OMP DO
    2768                 !$acc parallel present( pt, q, ql, qswst, rmask, tswst, sums_l ) create( s1 )
    2769                 s1 = 0
    2770                 !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2771                 DO  i = nxl, nxr
    2772                    DO  j =  nys, nyn
    2773                       s1 = s1 + ( ( 1.0_wp +                                   &
    2774                                     0.61_wp * q(nzt,j,i) - ql(nzt,j,i) ) *     &
    2775                                   tswst(j,i) +                                 &
    2776                                   0.61_wp * pt(nzt,j,i) * qswst(j,i) )         &
    2777                                 * heatflux_output_conversion(nzt)
    2778                    ENDDO
    2779                 ENDDO
    2780                 sums_l(nzt,45,tn) = s1
    2781                 !$acc end parallel
    2782 
    2783              ENDIF
    2784 
    2785              IF ( cloud_physics )  THEN
    2786 
    2787                 !$OMP DO
    2788                 !$acc parallel present( qswst, rmask, sums_l ) create( s1 )
    2789                 s1 = 0
    2790                 !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2791                 DO  i = nxl, nxr
    2792                    DO  j =  nys, nyn
    2793 !
    2794 !--                   Formula does not work if ql(nzb) /= 0.0
    2795                       s1 = s1 + qswst(j,i) * waterflux_output_conversion(nzt)  &
    2796                                            * rmask(j,i,sr)  ! w"q" (w"qv")
    2797                    ENDDO
    2798                 ENDDO
    2799                 sums_l(nzt,51,tn) = s1
    2800                 !$acc end parallel
    2801 
    2802              ENDIF
    2803 
    2804           ENDIF
    2805 
    2806           IF ( passive_scalar )  THEN
    2807 
    2808              !$OMP DO
    2809              !$acc parallel present( sswst, rmask, sums_l ) create( s1 )
    2810              s1 = 0
    2811              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2812              DO  i = nxl, nxr
    2813                 DO  j =  nys, nyn
    2814                    s1 = s1 + sswst(j,i) * rmask(j,i,sr) ! w"s"
    2815                 ENDDO
    2816              ENDDO
    2817              sums_l(nzt,119,tn) = s1
    2818              !$acc end parallel
    2819 
    2820           ENDIF
    2821 
    2822        ENDIF
    2823 
    2824 !
    2825 !--    Resolved fluxes (can be computed for all horizontal points)
    2826 !--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
    2827 !--    ----  speaking the following k-loop would have to be split up and
    2828 !--          rearranged according to the staggered grid.
    2829        !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2, s3 )
    2830        DO  k = nzb, nzt_diff
    2831           s1 = 0
    2832           s2 = 0
    2833           s3 = 0
    2834           !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
    2835           DO  i = nxl, nxr
    2836              DO  j = nys, nyn
    2837                 ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) + &
    2838                                  u(k+1,j,i) - hom(k+1,1,1,sr) )
    2839                 vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) + &
    2840                                  v(k+1,j,i) - hom(k+1,1,2,sr) )
    2841                 pts = 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) + &
    2842                                  pt(k+1,j,i) - hom(k+1,1,4,sr) )
    2843 !
    2844 !--             Higher moments
    2845                 s1 = s1 + pts * w(k,j,i)**2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2846                 s2 = s2 + pts**2 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2847 !
    2848 !--             Energy flux w*e* (has to be adjusted?)
    2849                 s3 = s3 + w(k,j,i) * 0.5_wp * ( ust**2 + vst**2 + w(k,j,i)**2 )&
    2850                                    * rmask(j,i,sr) * rflags_invers(j,i,k+1)    &
    2851