Changeset 2118 for palm/trunk/SOURCE/advec_ws.f90
- Timestamp:
- Jan 17, 2017 4:38:49 PM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/advec_ws.f90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_ws.f90
r2101 r2118 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! OpenACC version of subroutines removed 23 23 ! 24 24 ! Former revisions: … … 181 181 182 182 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 186 185 187 186 INTERFACE ws_init … … 207 206 END INTERFACE advec_u_ws 208 207 209 INTERFACE advec_u_ws_acc210 MODULE PROCEDURE advec_u_ws_acc211 END INTERFACE advec_u_ws_acc212 213 208 INTERFACE advec_v_ws 214 209 MODULE PROCEDURE advec_v_ws … … 216 211 END INTERFACE advec_v_ws 217 212 218 INTERFACE advec_v_ws_acc219 MODULE PROCEDURE advec_v_ws_acc220 END INTERFACE advec_v_ws_acc221 222 213 INTERFACE advec_w_ws 223 214 MODULE PROCEDURE advec_w_ws 224 215 MODULE PROCEDURE advec_w_ws_ij 225 216 END INTERFACE advec_w_ws 226 227 INTERFACE advec_w_ws_acc228 MODULE PROCEDURE advec_w_ws_acc229 END INTERFACE advec_w_ws_acc230 217 231 218 CONTAINS … … 4029 4016 4030 4017 4031 !------------------------------------------------------------------------------!4032 ! Description:4033 ! ------------4034 !> Scalar advection - Call for all grid points - accelerator version4035 !------------------------------------------------------------------------------!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_zw4040 4041 USE constants, &4042 ONLY: adv_sca_1, adv_sca_3, adv_sca_54043 4044 USE control_parameters, &4045 ONLY: intermediate_timestep_count, monotonic_adjustment, u_gtrans,&4046 v_gtrans4047 4048 USE grid_variables, &4049 ONLY: ddx, ddy4050 4051 USE indices, &4052 ONLY: i_left, i_right, j_north, j_south, nxlg, nxrg, nyng, nysg, &4053 nzb, nzb_max, nzt, wall_flags_04054 4055 USE kinds4056 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_substep4060 4061 IMPLICIT NONE4062 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 terms4121 !$acc kernels present( ddzw, sk, tend, u, v, w, wall_flags_0 )4122 DO i = i_left, i_right4123 DO j = j_south, j_north4124 DO k = nzb+1, nzt4125 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_gtrans4131 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_gtrans4166 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_gtrans4201 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_gtrans4236 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 point4268 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 * ibit84273 k_mm = k - 2 * ( ibit7 + ibit8 )4274 k_mmm = k - 3 * ibit84275 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 * ibit84310 k_pp = k + 2 * ( 1 - ibit6 )4311 k_mm = k - 2 * ibit84312 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 ) THEN4344 !4345 !-- At first, calculate first order fluxes.4346 u_comp = u(k,j,i) - u_gtrans4347 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_14350 4351 u_comp = u(k,j,i+1) - u_gtrans4352 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_14355 4356 v_comp = v(k,j,i) - v_gtrans4357 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_14360 4361 v_comp = v(k,j+1,i) - v_gtrans4362 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_14365 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 just4375 !-- 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 above4418 !-- 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 * ibit84421 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_l4466 diss_r = diss_r * phi_r4467 diss_s = diss_s * phi_s4468 diss_n = diss_n * phi_n4469 diss_d = diss_d * phi_d4470 diss_t = diss_t * phi_t4471 4472 ENDIF4473 !4474 !-- Calculate the divergence of the velocity field. A respective4475 !-- correction is needed to overcome numerical instabilities caused4476 !-- 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 statistics4509 ! 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 SELECT4533 4534 ENDDO4535 ENDDO4536 ENDDO4537 !$acc end kernels4538 4539 END SUBROUTINE advec_s_ws_acc4540 4541 4018 4542 4019 !------------------------------------------------------------------------------! … … 5039 4516 ! Description: 5040 4517 ! ------------ 5041 !> Advection of u - Call for all grid points - accelerator version5042 !------------------------------------------------------------------------------!5043 SUBROUTINE advec_u_ws_acc5044 5045 USE arrays_3d, &5046 ONLY: ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw5047 5048 USE constants, &5049 ONLY: adv_mom_1, adv_mom_3, adv_mom_55050 5051 USE control_parameters, &5052 ONLY: intermediate_timestep_count, u_gtrans, v_gtrans5053 5054 USE grid_variables, &5055 ONLY: ddx, ddy5056 5057 USE indices, &5058 ONLY: i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb, &5059 nzb_max, nzt, wall_flags_05060 5061 USE kinds5062 5063 ! USE statistics, &5064 ! ONLY: hom, sums_us2_ws_l, sums_wsus_ws_l, weight_substep5065 5066 IMPLICIT NONE5067 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_gtrans5109 gv = 2.0_wp * v_gtrans5110 5111 !5112 !-- Computation of fluxes and tendency terms5113 !$acc kernels present( ddzw, tend, u, v, w, wall_flags_0 )5114 DO i = i_left, i_right5115 DO j = j_south, j_north5116 DO k = nzb+1, nzt5117 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) - gu5123 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) - gv5193 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) - gv5229 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 * ibit175264 k_mm = k - 2 * ( ibit16 + ibit17 )5265 k_mmm = k - 3 * ibit175266 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 array5299 !-- 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 * ibit175305 k_pp = k + 2 * ( 1 - ibit15 )5306 k_mm = k - 2 * ibit175307 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 respective5340 !-- correction is needed to overcome numerical instabilities caused5341 !-- 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_wp5364 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 applied5376 !-- 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 ENDDO5391 ENDDO5392 ENDDO5393 !$acc end kernels5394 5395 !++5396 ! sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)5397 5398 END SUBROUTINE advec_u_ws_acc5399 5400 5401 !------------------------------------------------------------------------------!5402 ! Description:5403 ! ------------5404 4518 !> Advection of v - Call for all grid points 5405 4519 !------------------------------------------------------------------------------! … … 5906 5020 5907 5021 5908 !------------------------------------------------------------------------------!5909 ! Description:5910 ! ------------5911 !> Advection of v - Call for all grid points - accelerator version5912 !------------------------------------------------------------------------------!5913 SUBROUTINE advec_v_ws_acc5914 5915 USE arrays_3d, &5916 ONLY: ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw5917 5918 USE constants, &5919 ONLY: adv_mom_1, adv_mom_3, adv_mom_55920 5921 USE control_parameters, &5922 ONLY: intermediate_timestep_count, u_gtrans, v_gtrans5923 5924 USE grid_variables, &5925 ONLY: ddx, ddy5926 5927 USE indices, &5928 ONLY: i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb, &5929 nzb_max, nzt, wall_flags_05930 5931 USE kinds5932 5933 ! USE statistics, &5934 ! ONLY: hom, sums_vs2_ws_l, sums_wsvs_ws_l, weight_substep5935 5936 IMPLICIT NONE5937 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_gtrans5979 gv = 2.0_wp * v_gtrans5980 5981 !5982 !-- Computation of fluxes and tendency terms5983 !$acc kernels present( ddzw, tend, u, v, w, wall_flags_0 )5984 DO i = i_left, i_right5985 DO j = j_south, j_north5986 DO k = nzb+1, nzt5987 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) - gu5993 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) - gu6028 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) - gv6064 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 * ibit266134 k_mm = k - 2 * ( ibit25 + ibit26 )6135 k_mmm = k - 3 * ibit266136 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 array6169 !-- 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 * ibit266175 k_pp = k + 2 * ( 1 - ibit24 )6176 k_mm = k - 2 * ibit266177 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 respective6210 !-- correction is needed to overcome numerical instabilities caused6211 !-- 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_wp6237 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 applied6250 !-- 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 ENDDO6266 ENDDO6267 ENDDO6268 !$acc end kernels6269 6270 !++6271 ! sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)6272 6273 END SUBROUTINE advec_v_ws_acc6274 5022 6275 5023 … … 6756 5504 6757 5505 6758 !------------------------------------------------------------------------------!6759 ! Description:6760 ! ------------6761 !> Advection of w - Call for all grid points - accelerator version6762 !------------------------------------------------------------------------------!6763 SUBROUTINE advec_w_ws_acc6764 6765 USE arrays_3d, &6766 ONLY: ddzu, drho_air_zw, tend, u, v, w, rho_air, rho_air_zw6767 6768 USE constants, &6769 ONLY: adv_mom_1, adv_mom_3, adv_mom_56770 6771 USE control_parameters, &6772 ONLY: intermediate_timestep_count, u_gtrans, v_gtrans6773 6774 USE grid_variables, &6775 ONLY: ddx, ddy6776 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_006780 6781 USE kinds6782 6783 ! USE statistics, &6784 ! ONLY: hom, sums_ws2_ws_l, weight_substep6785 6786 IMPLICIT NONE6787 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_gtrans6828 gv = 2.0_wp * v_gtrans6829 6830 6831 !6832 !-- Computation of fluxes and tendency terms6833 !$acc kernels present( ddzu, tend, u, v, w, wall_flags_0, wall_flags_00 )6834 DO i = i_left, i_right6835 DO j = j_south, j_north6836 DO k = nzb+1, nzt6837 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) - gu6843 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) - gu6878 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) - gv6912 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) - gv6947 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 * ibit356982 k_mm = k - 2 * ( ibit34 + ibit35 )6983 k_mmm = k - 3 * ibit356984 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 array7018 !-- 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 * ibit357024 k_pp = k + 2 * ( 1 - ibit33 )7025 k_mm = k - 2 * ibit357026 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 respective7059 !-- correction is needed to overcome numerical instabilities caused7060 !-- 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_wp7083 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 ENDDO7101 ENDDO7102 ENDDO7103 !$acc end kernels7104 7105 END SUBROUTINE advec_w_ws_acc7106 7107 5506 END MODULE advec_ws
Note: See TracChangeset
for help on using the changeset viewer.
