Ignore:
Timestamp:
Sep 27, 2012 9:23:24 AM (9 years ago)
Author:
raasch
Message:

Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added.
Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed.
mask array is set zero for ghost boundaries

File:
1 edited

Legend:

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

    r1011 r1015  
    44! Current revisions:
    55! ------------------
    6 !
     6! accelerator versions (*_acc) added
    77!
    88! Former revisions:
     
    8888
    8989    PRIVATE
    90     PUBLIC   advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, &
     90    PUBLIC   advec_s_ws, advec_s_ws_acc, advec_u_ws, advec_u_ws_acc, &
     91             advec_v_ws, advec_v_ws_acc, advec_w_ws, advec_w_ws_acc, &
    9192             ws_init, ws_statistics
    9293
     
    109110    END INTERFACE advec_u_ws
    110111
     112    INTERFACE advec_u_ws_acc
     113       MODULE PROCEDURE advec_u_ws_acc
     114    END INTERFACE advec_u_ws_acc
     115
    111116    INTERFACE advec_v_ws
    112117       MODULE PROCEDURE advec_v_ws
     
    114119    END INTERFACE advec_v_ws
    115120
     121    INTERFACE advec_v_ws_acc
     122       MODULE PROCEDURE advec_v_ws_acc
     123    END INTERFACE advec_v_ws_acc
     124
    116125    INTERFACE advec_w_ws
    117126       MODULE PROCEDURE advec_w_ws
    118127       MODULE PROCEDURE advec_w_ws_ij
    119128    END INTERFACE advec_w_ws
     129
     130    INTERFACE advec_w_ws_acc
     131       MODULE PROCEDURE advec_w_ws_acc
     132    END INTERFACE advec_w_ws_acc
    120133
    121134 CONTAINS
     
    23352348
    23362349!------------------------------------------------------------------------------!
     2350! Scalar advection - Call for all grid points - accelerator version
     2351!------------------------------------------------------------------------------!
     2352    SUBROUTINE advec_s_ws_acc ( sk, sk_char )
     2353
     2354       USE arrays_3d
     2355       USE constants
     2356       USE control_parameters
     2357       USE grid_variables
     2358       USE indices
     2359       USE statistics
     2360
     2361       IMPLICIT NONE
     2362
     2363       CHARACTER (LEN = *), INTENT(IN)    :: sk_char
     2364
     2365       INTEGER ::  i, ibit0, ibit1, ibit2, ibit3, ibit4, ibit5, ibit6,        &
     2366                   ibit7, ibit8, j, k, k_mm, k_mmm, k_pp, k_ppp, tn = 0
     2367
     2368       REAL    :: diss_d, diss_l, diss_n, diss_r, diss_s, diss_t, div, flux_d, &
     2369                  flux_l, flux_n, flux_r, flux_s, flux_t, u_comp, v_comp
     2370
     2371       REAL, INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  ::  sk
     2372
     2373!
     2374!--    Computation of fluxes and tendency terms
     2375       !$acc kernels present( ddzw, sk, tend, u, v, w, wall_flags_0 )
     2376       !$acc loop
     2377       DO  i = nxl, nxr
     2378          DO  j = nys, nyn
     2379             !$acc loop vector( 32 )
     2380             DO  k = nzb+1, nzt
     2381
     2382                ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
     2383                ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
     2384                ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
     2385
     2386                u_comp              = u(k,j,i) - u_gtrans
     2387                flux_l              = u_comp * (                           &
     2388                                               ( 37.0 * ibit2 * adv_sca_5  &
     2389                                            +     7.0 * ibit1 * adv_sca_3  &
     2390                                            +           ibit0 * adv_sca_1  &
     2391                                               ) *                         &
     2392                                         ( sk(k,j,i)   + sk(k,j,i-1)    )  &
     2393                                        -      (  8.0 * ibit2 * adv_sca_5  &
     2394                                            +           ibit1 * adv_sca_3  &
     2395                                               ) *                         &
     2396                                         ( sk(k,j,i+1) + sk(k,j,i-2)    )  &
     2397                                        +      (        ibit2 * adv_sca_5  &
     2398                                               ) *                         &
     2399                                         ( sk(k,j,i+2) + sk(k,j,i-3)    )  &
     2400                                            )
     2401
     2402                diss_l              = -ABS( u_comp ) * (                   &
     2403                                               ( 10.0 * ibit2 * adv_sca_5  &
     2404                                            +     3.0 * ibit1 * adv_sca_3  &
     2405                                            +           ibit0 * adv_sca_1  &
     2406                                               ) *                         &
     2407                                         ( sk(k,j,i)   - sk(k,j,i-1)    )  &
     2408                                        -      (  5.0 * ibit2 * adv_sca_5  &
     2409                                            +           ibit1 * adv_sca_3  &
     2410                                               ) *                         &
     2411                                         ( sk(k,j,i+1) - sk(k,j,i-2)  )    &
     2412                                        +      (        ibit2 * adv_sca_5  &
     2413                                               ) *                         &
     2414                                         ( sk(k,j,i+2) - sk(k,j,i-3) )     &
     2415                                                    )
     2416
     2417                u_comp    = u(k,j,i+1) - u_gtrans
     2418                flux_r    = u_comp * (                                      &
     2419                          ( 37.0 * ibit2 * adv_sca_5                        &
     2420                      +      7.0 * ibit1 * adv_sca_3                        &
     2421                      +            ibit0 * adv_sca_1                        &
     2422                          ) *                                               &
     2423                             ( sk(k,j,i+1) + sk(k,j,i)   )                  &
     2424                   -      (  8.0 * ibit2 * adv_sca_5                        &
     2425                       +           ibit1 * adv_sca_3                        &
     2426                          ) *                                               &
     2427                             ( sk(k,j,i+2) + sk(k,j,i-1) )                  &
     2428                   +      (        ibit2 * adv_sca_5                        &
     2429                          ) *                                               &
     2430                             ( sk(k,j,i+3) + sk(k,j,i-2) )                  &
     2431                                     )
     2432
     2433                diss_r    = -ABS( u_comp ) * (                              &
     2434                          ( 10.0 * ibit2 * adv_sca_5                        &
     2435                       +     3.0 * ibit1 * adv_sca_3                        &
     2436                       +           ibit0 * adv_sca_1                        &
     2437                          ) *                                               &
     2438                             ( sk(k,j,i+1) - sk(k,j,i)  )                   &
     2439                   -      (  5.0 * ibit2 * adv_sca_5                        &
     2440                       +           ibit1 * adv_sca_3                        &
     2441                          ) *                                               &
     2442                             ( sk(k,j,i+2) - sk(k,j,i-1) )                  &
     2443                   +      (        ibit2 * adv_sca_5                        &
     2444                          ) *                                               &
     2445                             ( sk(k,j,i+3) - sk(k,j,i-2) )                  &
     2446                                             )
     2447
     2448                ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
     2449                ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
     2450                ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
     2451
     2452                v_comp               = v(k,j,i) - v_gtrans
     2453                flux_s               = v_comp * (                             &
     2454                                                  ( 37.0 * ibit5 * adv_sca_5  &
     2455                                               +     7.0 * ibit4 * adv_sca_3  &
     2456                                               +           ibit3 * adv_sca_1  &
     2457                                                  ) *                         &
     2458                                           ( sk(k,j,i)  + sk(k,j-1,i)     )   &
     2459                                            -     (  8.0 * ibit5 * adv_sca_5  &
     2460                                               +           ibit4 * adv_sca_3  &
     2461                                                   ) *                        &
     2462                                           ( sk(k,j+1,i) + sk(k,j-2,i)    )   &
     2463                                           +      (        ibit5 * adv_sca_5  &
     2464                                                  ) *                         &
     2465                                           ( sk(k,j+2,i) + sk(k,j-3,i)    )   &
     2466                                             )
     2467
     2468                diss_s               = -ABS( v_comp ) * (                     &
     2469                                                  ( 10.0 * ibit5 * adv_sca_5  &
     2470                                               +     3.0 * ibit4 * adv_sca_3  &
     2471                                               +           ibit3 * adv_sca_1  &
     2472                                                  ) *                         &
     2473                                            ( sk(k,j,i)   - sk(k,j-1,i)    )  &
     2474                                           -      (  5.0 * ibit5 * adv_sca_5  &
     2475                                               +           ibit4 * adv_sca_3  &
     2476                                            ) *                               &
     2477                                            ( sk(k,j+1,i) - sk(k,j-2,i)  )    &
     2478                                           +      (        ibit5 * adv_sca_5  &
     2479                                                  ) *                         &
     2480                                            ( sk(k,j+2,i) - sk(k,j-3,i) )     &
     2481                                                     )
     2482
     2483
     2484                v_comp    = v(k,j+1,i) - v_gtrans
     2485                flux_n    = v_comp * (                                      &
     2486                          ( 37.0 * ibit5 * adv_sca_5                        &
     2487                       +     7.0 * ibit4 * adv_sca_3                        &
     2488                       +           ibit3 * adv_sca_1                        &
     2489                          ) *                                               &
     2490                             ( sk(k,j+1,i) + sk(k,j,i)   )                  &
     2491                   -      (  8.0 * ibit5 * adv_sca_5                        &
     2492                       +           ibit4 * adv_sca_3                        &
     2493                          ) *                                               &
     2494                             ( sk(k,j+2,i) + sk(k,j-1,i) )                  &
     2495                   +      (        ibit5 * adv_sca_5                        &
     2496                          ) *                                               &
     2497                             ( sk(k,j+3,i) + sk(k,j-2,i) )                  &
     2498                                     )
     2499
     2500                diss_n    = -ABS( v_comp ) * (                              &
     2501                          ( 10.0 * ibit5 * adv_sca_5                        &
     2502                       +     3.0 * ibit4 * adv_sca_3                        &
     2503                       +           ibit3 * adv_sca_1                        &
     2504                          ) *                                               &
     2505                             ( sk(k,j+1,i) - sk(k,j,i)    )                 &
     2506                   -      (  5.0 * ibit5 * adv_sca_5                        &
     2507                       +           ibit4 * adv_sca_3                        &
     2508                          ) *                                               &
     2509                             ( sk(k,j+2,i) - sk(k,j-1,i)  )                 &
     2510                   +      (        ibit5 * adv_sca_5                        &
     2511                          ) *                                               &
     2512                             ( sk(k,j+3,i) - sk(k,j-2,i) )                  &
     2513                                             )
     2514
     2515!
     2516!--             indizes k_m, k_mm, ... should be known at these point
     2517                ibit8 = IBITS(wall_flags_0(k-1,j,i),8,1)
     2518                ibit7 = IBITS(wall_flags_0(k-1,j,i),7,1)
     2519                ibit6 = IBITS(wall_flags_0(k-1,j,i),6,1)
     2520
     2521                k_pp  = k + 2 * ( 1 - ibit6 )
     2522                k_mm  = k - 2 * ( 1 - ibit6 )
     2523                k_mmm = k - 3 * ibit8
     2524
     2525                flux_d    = w(k-1,j,i) * (                                    &
     2526                           ( 37.0 * ibit8 * adv_sca_5                         &
     2527                        +     7.0 * ibit7 * adv_sca_3                         &
     2528                        +           ibit6 * adv_sca_1                         &
     2529                           ) *                                                &
     2530                                   ( sk(k,j,i)    + sk(k-1,j,i) )             &
     2531                          -      (  8.0 * ibit8 * adv_sca_5                   &
     2532                          +               ibit7 * adv_sca_3                   &
     2533                           ) *                                                &
     2534                                   ( sk(k+1,j,i) + sk(k_mm,j,i) )             &
     2535                    +      (        ibit8 * adv_sca_5                         &
     2536                           ) *     ( sk(k_pp,j,i)+ sk(k_mmm,j,i) )            &
     2537                                       )
     2538
     2539                diss_d    = -ABS( w(k-1,j,i) ) * (                            &
     2540                           ( 10.0 * ibit8 * adv_sca_5                         &
     2541                        +     3.0 * ibit7 * adv_sca_3                         &
     2542                        +           ibit6 * adv_sca_1                         &
     2543                           ) *                                                &
     2544                                   ( sk(k,j,i)    - sk(k-1,j,i)   )           &
     2545                    -      (  5.0 * ibit8 * adv_sca_5                         &
     2546                        +           ibit7 * adv_sca_3                         &
     2547                           ) *                                                &
     2548                                   ( sk(k+1,j,i)  - sk(k_mm,j,i)  )           &
     2549                    +      (        ibit8 * adv_sca_5                         &
     2550                           ) *                                                &
     2551                                   ( sk(k_pp,j,i) - sk(k_mmm,j,i) )           &
     2552                                         )
     2553
     2554                ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
     2555                ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
     2556                ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
     2557
     2558                k_ppp = k + 3 * ibit8
     2559                k_pp  = k + 2 * ( 1 - ibit6  )
     2560                k_mm  = k - 2 * ibit8
     2561
     2562                flux_t    = w(k,j,i) * (                                      &
     2563                           ( 37.0 * ibit8 * adv_sca_5                         &
     2564                        +     7.0 * ibit7 * adv_sca_3                         &
     2565                        +           ibit6 * adv_sca_1                         &
     2566                           ) *                                                &
     2567                                   ( sk(k+1,j,i)  + sk(k,j,i)   )             &
     2568                          -      (  8.0 * ibit8 * adv_sca_5                   &
     2569                        +                 ibit7 * adv_sca_3                   &
     2570                           ) *                                                &
     2571                                   ( sk(k_pp,j,i) + sk(k-1,j,i) )             &
     2572                    +      (        ibit8 * adv_sca_5                         &
     2573                           ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )            &
     2574                                       )
     2575
     2576                diss_t    = -ABS( w(k,j,i) ) * (                              &
     2577                           ( 10.0 * ibit8 * adv_sca_5                         &
     2578                        +     3.0 * ibit7 * adv_sca_3                         &
     2579                        +           ibit6 * adv_sca_1                         &
     2580                           ) *                                                &
     2581                                   ( sk(k+1,j,i)   - sk(k,j,i)    )           &
     2582                    -      (  5.0 * ibit8 * adv_sca_5                         &
     2583                        +           ibit7 * adv_sca_3                         &
     2584                           ) *                                                &
     2585                                   ( sk(k_pp,j,i)  - sk(k-1,j,i)  )           &
     2586                    +      (        ibit8 * adv_sca_5                         &
     2587                           ) *                                                &
     2588                                   ( sk(k_ppp,j,i) - sk(k_mm,j,i) )           &
     2589                                         )
     2590!
     2591!--             Calculate the divergence of the velocity field. A respective
     2592!--             correction is needed to overcome numerical instabilities caused
     2593!--             by a not sufficient reduction of divergences near topography.
     2594                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
     2595                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
     2596                              + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
     2597
     2598                tend(k,j,i) = - (                                             &
     2599                               ( flux_r + diss_r - flux_l - diss_l ) * ddx    &
     2600                             + ( flux_n + diss_n - flux_s - diss_s ) * ddy    &
     2601                             + ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k)&
     2602                                ) + div * sk(k,j,i)
     2603
     2604!++
     2605!--             Evaluation of statistics
     2606!                SELECT CASE ( sk_char )
     2607!
     2608!                   CASE ( 'pt' )
     2609!                      sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn)         &
     2610!                       + ( flux_t + diss_t )                                &
     2611!                       *   weight_substep(intermediate_timestep_count)
     2612!                   CASE ( 'sa' )
     2613!                      sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn)         &
     2614!                       + ( flux_t + diss_t )                                &
     2615!                       *   weight_substep(intermediate_timestep_count)
     2616!                   CASE ( 'q' )
     2617!                      sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn)           &
     2618!                      + ( flux_t + diss_t )                                &
     2619!                      *   weight_substep(intermediate_timestep_count)
     2620!
     2621!                END SELECT
     2622
     2623             ENDDO
     2624         ENDDO
     2625      ENDDO
     2626      !$acc end kernels
     2627
     2628    END SUBROUTINE advec_s_ws_acc
     2629
     2630
     2631!------------------------------------------------------------------------------!
    23372632! Advection of u - Call for all grid points
    23382633!------------------------------------------------------------------------------!
     
    27523047   
    27533048   
     3049!------------------------------------------------------------------------------!
     3050! Advection of u - Call for all grid points - accelerator version
     3051!------------------------------------------------------------------------------!
     3052    SUBROUTINE advec_u_ws_acc
     3053
     3054       USE arrays_3d
     3055       USE constants
     3056       USE control_parameters
     3057       USE grid_variables
     3058       USE indices
     3059       USE statistics
     3060
     3061       IMPLICIT NONE
     3062
     3063       INTEGER ::  i, ibit9, ibit10, ibit11, ibit12, ibit13, ibit14, ibit15,   &
     3064                   ibit16, ibit17, j, k, k_mmm, k_mm, k_pp, k_ppp, tn = 0
     3065
     3066       REAL    ::  diss_d, diss_l, diss_n, diss_r, diss_s, diss_t, div,    &
     3067                   flux_d, flux_l, flux_n, flux_r, flux_s, flux_t, gu, gv, &
     3068                   u_comp, u_comp_l, v_comp, v_comp_s, w_comp
     3069
     3070
     3071       gu = 2.0 * u_gtrans
     3072       gv = 2.0 * v_gtrans
     3073
     3074!
     3075!--    Computation of fluxes and tendency terms
     3076       !$acc  kernels present( ddzw, tend, u, v, w, wall_flags_0 )
     3077       !$acc  loop
     3078       DO i = nxlu, nxr
     3079          DO  j = nys, nyn
     3080             !$acc  loop vector( 32 )
     3081             DO  k = nzb+1, nzt
     3082
     3083                ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
     3084                ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
     3085                ibit9  = IBITS(wall_flags_0(k,j,i),9,1)
     3086
     3087                u_comp_l           = u(k,j,i) + u(k,j,i-1) - gu
     3088                flux_l             = u_comp_l * (                          &
     3089                                    ( 37.0 * ibit11 * adv_mom_5             &
     3090                                 +     7.0 * ibit10 * adv_mom_3             &
     3091                                 +           ibit9  * adv_mom_1             &
     3092                                    ) *                                     &
     3093                                  ( u(k,j,i)   + u(k,j,i-1) )               &
     3094                             -      (  8.0 * ibit11 * adv_mom_5             &
     3095                                 +           ibit10 * adv_mom_3             &
     3096                                    ) *                                     &
     3097                                  ( u(k,j,i+1) + u(k,j,i-2) )               &
     3098                             +      (        ibit11 * adv_mom_5             &
     3099                                    ) *                                     &
     3100                                  ( u(k,j,i+2) + u(k,j,i-3) )               &
     3101                                                )
     3102
     3103                diss_l             = - ABS( u_comp_l ) * (                &
     3104                                   ( 10.0 * ibit11 * adv_mom_5             &
     3105                                +     3.0 * ibit10 * adv_mom_3             &
     3106                                +           ibit9  * adv_mom_1             &
     3107                                   ) *                                     &
     3108                                 ( u(k,j,i)   - u(k,j,i-1) )               &
     3109                            -      (  5.0 * ibit11 * adv_mom_5             &
     3110                                +           ibit10 * adv_mom_3             &
     3111                                   ) *                                     &
     3112                                 ( u(k,j,i+1) - u(k,j,i-2) )               &
     3113                            +      (        ibit11 * adv_mom_5             &
     3114                                   ) *                                     &
     3115                                 ( u(k,j,i+2) - u(k,j,i-3) )               &
     3116                                                         )
     3117
     3118                u_comp    = u(k,j,i+1) + u(k,j,i)
     3119                flux_r    = ( u_comp   - gu ) * (                           &
     3120                          ( 37.0 * ibit11 * adv_mom_5                        &
     3121                       +     7.0 * ibit10 * adv_mom_3                        &
     3122                       +           ibit9  * adv_mom_1                        &
     3123                          ) *                                                &
     3124                                 ( u(k,j,i+1) + u(k,j,i)   )                 &
     3125                   -      (  8.0 * ibit11 * adv_mom_5                        &
     3126                       +           ibit10 * adv_mom_3                        &
     3127                          ) *                                                &
     3128                                 ( u(k,j,i+2) + u(k,j,i-1) )                 &
     3129                   +      (        ibit11 * adv_mom_5                        &
     3130                          ) *                                                &
     3131                                 ( u(k,j,i+3) + u(k,j,i-2) )                 &
     3132                                                 )
     3133
     3134                diss_r    = - ABS( u_comp    - gu ) * (                      &
     3135                          ( 10.0 * ibit11 * adv_mom_5                        &
     3136                       +     3.0 * ibit10 * adv_mom_3                        &
     3137                       +           ibit9  * adv_mom_1                        &
     3138                          ) *                                                &
     3139                                 ( u(k,j,i+1) - u(k,j,i)  )                  &
     3140                   -      (  5.0 * ibit11 * adv_mom_5                        &
     3141                       +           ibit10 * adv_mom_3                        &
     3142                          ) *                                                &
     3143                                 ( u(k,j,i+2) - u(k,j,i-1) )                 &
     3144                   +      (        ibit11 * adv_mom_5                        &
     3145                          ) *                                                &
     3146                                 ( u(k,j,i+3) - u(k,j,i-2) )                 &
     3147                                                     )
     3148
     3149                ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
     3150                ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
     3151                ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
     3152
     3153                v_comp_s                 = v(k,j,i) + v(k,j,i-1) - gv
     3154                flux_s                   = v_comp_s * (                       &
     3155                                   ( 37.0 * ibit14 * adv_mom_5                &
     3156                                +     7.0 * ibit13 * adv_mom_3                &
     3157                                +           ibit12 * adv_mom_1                &
     3158                                   ) *                                        &
     3159                                     ( u(k,j,i)   + u(k,j-1,i) )              &
     3160                            -      (  8.0 * ibit14 * adv_mom_5                &
     3161                            +           ibit13 * adv_mom_3                    &
     3162                                   ) *                                        &
     3163                                     ( u(k,j+1,i) + u(k,j-2,i) )              &
     3164                        +      (        ibit14 * adv_mom_5                    &
     3165                               ) *                                            &
     3166                                     ( u(k,j+2,i) + u(k,j-3,i) )              &
     3167                                               )
     3168
     3169                diss_s                  = - ABS ( v_comp_s ) * (              &
     3170                                   ( 10.0 * ibit14 * adv_mom_5                &
     3171                                +     3.0 * ibit13 * adv_mom_3                &
     3172                                +           ibit12 * adv_mom_1                &
     3173                                   ) *                                        &
     3174                                     ( u(k,j,i)   - u(k,j-1,i) )              &
     3175                            -      (  5.0 * ibit14 * adv_mom_5                &
     3176                                +           ibit13 * adv_mom_3                &
     3177                                   ) *                                        &
     3178                                     ( u(k,j+1,i) - u(k,j-2,i) )              &
     3179                            +      (        ibit14 * adv_mom_5                &
     3180                                   ) *                                        &
     3181                                     ( u(k,j+2,i) - u(k,j-3,i) )              &
     3182                                                         )
     3183
     3184
     3185                v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
     3186                flux_n    = v_comp * (                                       &
     3187                          ( 37.0 * ibit14 * adv_mom_5                        &
     3188                       +     7.0 * ibit13 * adv_mom_3                        &
     3189                       +           ibit12 * adv_mom_1                        &
     3190                          ) *                                                &
     3191                                 ( u(k,j+1,i) + u(k,j,i)   )                 &
     3192                   -      (  8.0 * ibit14 * adv_mom_5                        &
     3193                       +           ibit13 * adv_mom_3                        &
     3194                          ) *                                                &
     3195                                 ( u(k,j+2,i) + u(k,j-1,i) )                 &
     3196                   +      (        ibit14 * adv_mom_5                        &
     3197                          ) *                                                &
     3198                                 ( u(k,j+3,i) + u(k,j-2,i) )                 &
     3199                                                 )
     3200
     3201                diss_n    = - ABS ( v_comp ) * (                             &
     3202                          ( 10.0 * ibit14 * adv_mom_5                        &
     3203                       +     3.0 * ibit13 * adv_mom_3                        &
     3204                       +           ibit12 * adv_mom_1                        &
     3205                          ) *                                                &
     3206                                 ( u(k,j+1,i) - u(k,j,i)  )                  &
     3207                   -      (  5.0 * ibit14 * adv_mom_5                        &
     3208                       +           ibit13 * adv_mom_3                        &
     3209                          ) *                                                &
     3210                                 ( u(k,j+2,i) - u(k,j-1,i) )                 &
     3211                   +      (        ibit14 * adv_mom_5                        &
     3212                          ) *                                                &
     3213                                 ( u(k,j+3,i) - u(k,j-2,i) )                 &
     3214                                                      )
     3215
     3216                ibit17 = IBITS(wall_flags_0(k-1,j,i),17,1)
     3217                ibit16 = IBITS(wall_flags_0(k-1,j,i),16,1)
     3218                ibit15 = IBITS(wall_flags_0(k-1,j,i),15,1)
     3219
     3220                k_pp  = k + 2 * ( 1 - ibit15 )
     3221                k_mm  = k - 2 * ( 1 - ibit15 )
     3222                k_mmm = k - 3 * ibit17
     3223
     3224                w_comp    = w(k-1,j,i) + w(k-1,j,i-1)
     3225                flux_d    = w_comp  * (                                      &
     3226                          ( 37.0 * ibit17 * adv_mom_5                        &
     3227                       +     7.0 * ibit16 * adv_mom_3                        &
     3228                       +           ibit15 * adv_mom_1                        &
     3229                          ) *                                                &
     3230                             ( u(k,j,i)    + u(k-1,j,i)   )                  &
     3231                   -      (  8.0 * ibit17 * adv_mom_5                        &
     3232                       +           ibit16 * adv_mom_3                        &
     3233                          ) *                                                &
     3234                             ( u(k+1,j,i) + u(k_mm,j,i)   )                  &
     3235                   +      (        ibit17 * adv_mom_5                        &
     3236                          ) *                                                &
     3237                             ( u(k_pp,j,i) + u(k_mmm,j,i) )                  &
     3238                                      )
     3239
     3240                diss_d    = - ABS( w_comp ) * (                              &
     3241                          ( 10.0 * ibit17 * adv_mom_5                        &
     3242                       +     3.0 * ibit16 * adv_mom_3                        &
     3243                       +           ibit15 * adv_mom_1                        &
     3244                          ) *                                                &
     3245                             ( u(k,j,i)     - u(k-1,j,i)  )                  &
     3246                   -      (  5.0 * ibit17 * adv_mom_5                        &
     3247                       +           ibit16 * adv_mom_3                        &
     3248                          ) *                                                &
     3249                             ( u(k+1,j,i)  - u(k_mm,j,i)  )                  &
     3250                   +      (        ibit17 * adv_mom_5                        &
     3251                           ) *                                               &
     3252                             ( u(k_pp,j,i) - u(k_mmm,j,i) )                  &
     3253                                              )
     3254!
     3255!--             k index has to be modified near bottom and top, else array
     3256!--             subscripts will be exceeded.
     3257                ibit17 = IBITS(wall_flags_0(k,j,i),17,1)
     3258                ibit16 = IBITS(wall_flags_0(k,j,i),16,1)
     3259                ibit15 = IBITS(wall_flags_0(k,j,i),15,1)
     3260
     3261                k_ppp = k + 3 * ibit17
     3262                k_pp  = k + 2 * ( 1 - ibit15  )
     3263                k_mm  = k - 2 * ibit17
     3264
     3265                w_comp    = w(k,j,i) + w(k,j,i-1)
     3266                flux_t    = w_comp  * (                                      &
     3267                          ( 37.0 * ibit17 * adv_mom_5                        &
     3268                       +     7.0 * ibit16 * adv_mom_3                        &
     3269                       +           ibit15 * adv_mom_1                        &
     3270                          ) *                                                &
     3271                             ( u(k+1,j,i)  + u(k,j,i)     )                  &
     3272                   -      (  8.0 * ibit17 * adv_mom_5                        &
     3273                       +           ibit16 * adv_mom_3                        &
     3274                          ) *                                                &
     3275                             ( u(k_pp,j,i) + u(k-1,j,i)   )                  &
     3276                   +      (        ibit17 * adv_mom_5                        &
     3277                          ) *                                                &
     3278                             ( u(k_ppp,j,i) + u(k_mm,j,i) )                  &
     3279                                      )
     3280
     3281                diss_t    = - ABS( w_comp ) * (                              &
     3282                          ( 10.0 * ibit17 * adv_mom_5                        &
     3283                       +     3.0 * ibit16 * adv_mom_3                        &
     3284                       +           ibit15 * adv_mom_1                        &
     3285                          ) *                                                &
     3286                             ( u(k+1,j,i)   - u(k,j,i)    )                  &
     3287                   -      (  5.0 * ibit17 * adv_mom_5                        &
     3288                       +           ibit16 * adv_mom_3                        &
     3289                          ) *                                                &
     3290                             ( u(k_pp,j,i)  - u(k-1,j,i)  )                  &
     3291                   +      (        ibit17 * adv_mom_5                        &
     3292                           ) *                                               &
     3293                             ( u(k_ppp,j,i) - u(k_mm,j,i) )                  &
     3294                                              )
     3295!
     3296!--             Calculate the divergence of the velocity field. A respective
     3297!--             correction is needed to overcome numerical instabilities caused
     3298!--             by a not sufficient reduction of divergences near topography.
     3299                div = ( ( u_comp      - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx  &
     3300                     +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy  &
     3301                     +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) )        &
     3302                                                                    * ddzw(k)  &
     3303                      ) * 0.5
     3304
     3305                tend(k,j,i) = - (                                              &
     3306                               ( flux_r + diss_r - flux_l - diss_l ) * ddx     &
     3307                             + ( flux_n + diss_n - flux_s - diss_s ) * ddy     &
     3308                             + ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k) &
     3309                                ) + div * u(k,j,i)
     3310
     3311!++
     3312!--             Statistical Evaluation of u'u'. The factor has to be applied
     3313!--             for right evaluation when gallilei_trans = .T. .
     3314!                sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                     &
     3315!                              + ( flux_r    *                                 &
     3316!                                ( u_comp    - 2.0 * hom(k,1,1,0) )            &
     3317!                              / ( u_comp    - gu + 1.0E-20      )             &
     3318!                              +   diss_r    *                                 &
     3319!                                  ABS( u_comp    - 2.0 * hom(k,1,1,0) )       &
     3320!                              / ( ABS( u_comp    - gu ) + 1.0E-20 ) )         &
     3321!                              *   weight_substep(intermediate_timestep_count)
     3322!
     3323!--             Statistical Evaluation of w'u'.
     3324!                sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                   &
     3325!                              + ( flux_t    + diss_t    )                     &
     3326!                              *   weight_substep(intermediate_timestep_count)
     3327             ENDDO
     3328          ENDDO
     3329       ENDDO
     3330       !$acc end kernels
     3331
     3332!++
     3333!       sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)
     3334
     3335    END SUBROUTINE advec_u_ws_acc
     3336
     3337
    27543338!------------------------------------------------------------------------------!
    27553339! Advection of v - Call for all grid points
     
    31813765   
    31823766!------------------------------------------------------------------------------!
     3767! Advection of v - Call for all grid points - accelerator version
     3768!------------------------------------------------------------------------------!
     3769    SUBROUTINE advec_v_ws_acc
     3770
     3771       USE arrays_3d
     3772       USE constants
     3773       USE control_parameters
     3774       USE grid_variables
     3775       USE indices
     3776       USE statistics
     3777
     3778       IMPLICIT NONE
     3779
     3780
     3781       INTEGER ::  i, ibit18, ibit19, ibit20, ibit21, ibit22, ibit23, ibit24, &
     3782                    ibit25, ibit26, j, k, k_mm, k_mmm, k_pp, k_ppp, tn = 0
     3783
     3784       REAL    ::  diss_d, diss_l, diss_n, diss_r, diss_s, diss_t, div,    &
     3785                   flux_d, flux_l, flux_n, flux_r, flux_s, flux_t, gu, gv, &
     3786                   u_comp, u_comp_l, v_comp, v_comp_s, w_comp
     3787
     3788       gu = 2.0 * u_gtrans
     3789       gv = 2.0 * v_gtrans
     3790
     3791!
     3792!--    Computation of fluxes and tendency terms
     3793       !$acc kernels present( ddzw, tend, u, v, w, wall_flags_0 )
     3794       !$acc loop
     3795       DO  i = nxl, nxr
     3796          DO  j = nysv, nyn
     3797             !$acc loop vector( 32 )
     3798             DO  k = nzb+1, nzt
     3799
     3800                ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
     3801                ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
     3802                ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
     3803
     3804                u_comp_l                 = u(k,j-1,i) + u(k,j,i) - gu
     3805                flux_l                   = u_comp_l * (                          &
     3806                                      ( 37.0 * ibit20 * adv_mom_5              &
     3807                                   +     7.0 * ibit19 * adv_mom_3              &
     3808                                   +           ibit18 * adv_mom_1              &
     3809                                      ) *                                      &
     3810                                     ( v(k,j,i)   + v(k,j,i-1) )               &
     3811                               -      (  8.0 * ibit20 * adv_mom_5              &
     3812                                   +           ibit19 * adv_mom_3              &
     3813                                      ) *                                      &
     3814                                     ( v(k,j,i+1) + v(k,j,i-2) )               &
     3815                               +      (        ibit20 * adv_mom_5              &
     3816                                      ) *                                      &
     3817                                     ( v(k,j,i+2) + v(k,j,i-3) )               &
     3818                                                 )
     3819
     3820                diss_l                   = - ABS( u_comp_l ) * (                 &
     3821                                      ( 10.0 * ibit20 * adv_mom_5              &
     3822                                   +     3.0 * ibit19 * adv_mom_3              &
     3823                                   +           ibit18 * adv_mom_1              &
     3824                                      ) *                                      &
     3825                                     ( v(k,j,i)   - v(k,j,i-1) )               &
     3826                               -      (  5.0 * ibit20 * adv_mom_5              &
     3827                                   +           ibit19 * adv_mom_3              &
     3828                                      ) *                                      &
     3829                                     ( v(k,j,i+1) - v(k,j,i-2) )               &
     3830                               +      (        ibit20 * adv_mom_5              &
     3831                                      ) *                                      &
     3832                                     ( v(k,j,i+2) - v(k,j,i-3) )               &
     3833                                                           )
     3834
     3835                u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
     3836                flux_r    = u_comp * (                                       &
     3837                          ( 37.0 * ibit20 * adv_mom_5                        &
     3838                       +     7.0 * ibit19 * adv_mom_3                        &
     3839                       +           ibit18 * adv_mom_1                        &
     3840                          ) *                                                &
     3841                                 ( v(k,j,i+1) + v(k,j,i)   )                 &
     3842                   -      (  8.0 * ibit20 * adv_mom_5                        &
     3843                       +           ibit19 * adv_mom_3                        &
     3844                          ) *                                                &
     3845                                 ( v(k,j,i+2) + v(k,j,i-1) )                 &
     3846                   +      (        ibit20 * adv_mom_5                        &
     3847                          ) *                                                &
     3848                                 ( v(k,j,i+3) + v(k,j,i-2) )                 &
     3849                                     )
     3850
     3851                diss_r    = - ABS( u_comp ) * (                              &
     3852                          ( 10.0 * ibit20 * adv_mom_5                        &
     3853                       +     3.0 * ibit19 * adv_mom_3                        &
     3854                       +           ibit18 * adv_mom_1                        &
     3855                          ) *                                                &
     3856                                 ( v(k,j,i+1) - v(k,j,i)  )                  &
     3857                   -      (  5.0 * ibit20 * adv_mom_5                        &
     3858                       +           ibit19 * adv_mom_3                        &
     3859                          ) *                                                &
     3860                                 ( v(k,j,i+2) - v(k,j,i-1) )                 &
     3861                   +      (        ibit20 * adv_mom_5                        &
     3862                          ) *                                                &
     3863                                 ( v(k,j,i+3) - v(k,j,i-2) )                 &
     3864                                              )
     3865
     3866                ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
     3867                ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
     3868                ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
     3869
     3870
     3871                v_comp_s              = v(k,j,i) + v(k,j-1,i) - gv
     3872                flux_s                = v_comp_s    * (                       &
     3873                                   ( 37.0 * ibit23 * adv_mom_5                &
     3874                                +     7.0 * ibit22 * adv_mom_3                &
     3875                                +           ibit21 * adv_mom_1                &
     3876                                   ) *                                        &
     3877                                     ( v(k,j,i)   + v(k,j-1,i) )              &
     3878                            -      (  8.0 * ibit23 * adv_mom_5                &
     3879                                +           ibit22 * adv_mom_3                &
     3880                                   ) *                                        &
     3881                                     ( v(k,j+1,i) + v(k,j-2,i) )              &
     3882                            +      (        ibit23 * adv_mom_5                &
     3883                                   ) *                                        &
     3884                                     ( v(k,j+2,i) + v(k,j-3,i) )              &
     3885                                                 )
     3886
     3887                diss_s                = - ABS( v_comp_s ) * (                 &
     3888                                   ( 10.0 * ibit23 * adv_mom_5                &
     3889                                +     3.0 * ibit22 * adv_mom_3                &
     3890                                +           ibit21 * adv_mom_1                &
     3891                                   ) *                                        &
     3892                                     ( v(k,j,i)   - v(k,j-1,i) )              &
     3893                            -      (  5.0 * ibit23 * adv_mom_5                &
     3894                                +           ibit22 * adv_mom_3                &
     3895                                   ) *                                        &
     3896                                     ( v(k,j+1,i) - v(k,j-2,i) )              &
     3897                            +      (        ibit23 * adv_mom_5                &
     3898                                   ) *                                        &
     3899                                     ( v(k,j+2,i) - v(k,j-3,i) )              &
     3900                                                          )
     3901
     3902                v_comp = v(k,j+1,i) + v(k,j,i)
     3903                flux_n = ( v_comp - gv ) * (                                 &
     3904                          ( 37.0 * ibit23 * adv_mom_5                        &
     3905                       +     7.0 * ibit22 * adv_mom_3                        &
     3906                       +           ibit21 * adv_mom_1                        &
     3907                          ) *                                                &
     3908                                 ( v(k,j+1,i) + v(k,j,i)   )                 &
     3909                   -      (  8.0 * ibit23 * adv_mom_5                        &
     3910                       +           ibit22 * adv_mom_3                        &
     3911                          ) *                                                &
     3912                                 ( v(k,j+2,i) + v(k,j-1,i) )                 &
     3913                   +      (        ibit23 * adv_mom_5                        &
     3914                          ) *                                                &
     3915                                 ( v(k,j+3,i) + v(k,j-2,i) )                 &
     3916                                     )
     3917
     3918                diss_n = - ABS( v_comp - gv ) * (                         &
     3919                          ( 10.0 * ibit23 * adv_mom_5                        &
     3920                       +     3.0 * ibit22 * adv_mom_3                        &
     3921                       +           ibit21 * adv_mom_1                        &
     3922                          ) *                                                &
     3923                                 ( v(k,j+1,i) - v(k,j,i)  )                  &
     3924                   -      (  5.0 * ibit23 * adv_mom_5                        &
     3925                       +           ibit22 * adv_mom_3                        &
     3926                          ) *                                                &
     3927                                 ( v(k,j+2,i) - v(k,j-1,i) )                 &
     3928                   +      (        ibit23 * adv_mom_5                        &
     3929                          ) *                                                &
     3930                                 ( v(k,j+3,i) - v(k,j-2,i) )                 &
     3931                                                     )
     3932
     3933                ibit26 = IBITS(wall_flags_0(k-1,j,i),26,1)
     3934                ibit25 = IBITS(wall_flags_0(k-1,j,i),25,1)
     3935                ibit24 = IBITS(wall_flags_0(k-1,j,i),24,1)
     3936
     3937                k_pp  = k + 2 * ( 1 - ibit24 )
     3938                k_mm  = k - 2 * ( 1 - ibit24 )
     3939                k_mmm = k - 3 * ibit26
     3940
     3941                w_comp    = w(k-1,j-1,i) + w(k-1,j,i)
     3942                flux_d    = w_comp  * (                                      &
     3943                          ( 37.0 * ibit26 * adv_mom_5                        &
     3944                       +     7.0 * ibit25 * adv_mom_3                        &
     3945                       +           ibit24 * adv_mom_1                        &
     3946                          ) *                                                &
     3947                             ( v(k,j,i)     + v(k-1,j,i)  )                  &
     3948                   -      (  8.0 * ibit26 * adv_mom_5                        &
     3949                       +           ibit25 * adv_mom_3                        &
     3950                          ) *                                                &
     3951                             ( v(k+1,j,i)  + v(k_mm,j,i)  )                  &
     3952                   +      (        ibit26 * adv_mom_5                        &
     3953                          ) *                                                &
     3954                             ( v(k_pp,j,i) + v(k_mmm,j,i) )                  &
     3955                                      )
     3956
     3957                diss_d    = - ABS( w_comp ) * (                              &
     3958                          ( 10.0 * ibit26 * adv_mom_5                        &
     3959                       +     3.0 * ibit25 * adv_mom_3                        &
     3960                       +           ibit24 * adv_mom_1                        &
     3961                          ) *                                                &
     3962                             ( v(k,j,i)     - v(k-1,j,i)  )                  &
     3963                   -      (  5.0 * ibit26 * adv_mom_5                        &
     3964                       +           ibit25 * adv_mom_3                        &
     3965                          ) *                                                &
     3966                             ( v(k+1,j,i)  - v(k_mm,j,i)  )                  &
     3967                   +      (        ibit26 * adv_mom_5                        &
     3968                          ) *                                                &
     3969                             ( v(k_pp,j,i) - v(k_mmm,j,i) )                  &
     3970                                               )
     3971!
     3972!--             k index has to be modified near bottom and top, else array
     3973!--             subscripts will be exceeded.
     3974                ibit26 = IBITS(wall_flags_0(k,j,i),26,1)
     3975                ibit25 = IBITS(wall_flags_0(k,j,i),25,1)
     3976                ibit24 = IBITS(wall_flags_0(k,j,i),24,1)
     3977
     3978                k_ppp = k + 3 * ibit26
     3979                k_pp  = k + 2 * ( 1 - ibit24  )
     3980                k_mm  = k - 2 * ibit26
     3981
     3982                w_comp    = w(k,j-1,i) + w(k,j,i)
     3983                flux_t    = w_comp  * (                                      &
     3984                          ( 37.0 * ibit26 * adv_mom_5                        &
     3985                       +     7.0 * ibit25 * adv_mom_3                        &
     3986                       +           ibit24 * adv_mom_1                        &
     3987                          ) *                                                &
     3988                             ( v(k+1,j,i)   + v(k,j,i)    )                  &
     3989                   -      (  8.0 * ibit26 * adv_mom_5                        &
     3990                       +           ibit25 * adv_mom_3                        &
     3991                          ) *                                                &
     3992                             ( v(k_pp,j,i)  + v(k-1,j,i)  )                  &
     3993                   +      (        ibit26 * adv_mom_5                        &
     3994                          ) *                                                &
     3995                             ( v(k_ppp,j,i) + v(k_mm,j,i) )                  &
     3996                                      )
     3997
     3998                diss_t    = - ABS( w_comp ) * (                              &
     3999                          ( 10.0 * ibit26 * adv_mom_5                        &
     4000                       +     3.0 * ibit25 * adv_mom_3                        &
     4001                       +           ibit24 * adv_mom_1                        &
     4002                          ) *                                                &
     4003                             ( v(k+1,j,i)   - v(k,j,i)    )                  &
     4004                   -      (  5.0 * ibit26 * adv_mom_5                        &
     4005                       +           ibit25 * adv_mom_3                        &
     4006                          ) *                                                &
     4007                             ( v(k_pp,j,i)  - v(k-1,j,i)  )                  &
     4008                   +      (        ibit26 * adv_mom_5                        &
     4009                          ) *                                                &
     4010                             ( v(k_ppp,j,i) - v(k_mm,j,i) )                  &
     4011                                               )
     4012!
     4013!--             Calculate the divergence of the velocity field. A respective
     4014!--             correction is needed to overcome numerical instabilities caused
     4015!--             by a not sufficient reduction of divergences near topography.
     4016                div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx &
     4017                     +  ( v_comp      - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy &
     4018                     +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) )         &
     4019                                                                  ) * ddzw(k) &
     4020                      ) * 0.5
     4021
     4022                tend(k,j,i) = - (                                              &
     4023                               ( flux_r + diss_r - flux_l - diss_l ) * ddx     &
     4024                             + ( flux_n + diss_n - flux_s - diss_s ) * ddy     &
     4025                             + ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k) &
     4026                                ) + div * v(k,j,i)
     4027
     4028
     4029!++
     4030!--             Statistical Evaluation of v'v'. The factor has to be applied
     4031!--             for right evaluation when gallilei_trans = .T. .
     4032!                sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                  &
     4033!                      + ( flux_n                                           &
     4034!                      * ( v_comp - 2.0 * hom(k,1,2,0) )                    &
     4035!                      / ( v_comp - gv + 1.0E-20 )                          &
     4036!                      +   diss_n                                           &
     4037!                      *   ABS( v_comp - 2.0 * hom(k,1,2,0) )               &
     4038!                      / ( ABS( v_comp - gv ) +1.0E-20 ) )                  &
     4039!                      *   weight_substep(intermediate_timestep_count)
     4040!
     4041!--              Statistical Evaluation of w'v'.
     4042!                 sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                &
     4043!                              + ( flux_t + diss_t )                         &
     4044!                              *   weight_substep(intermediate_timestep_count)
     4045
     4046             ENDDO
     4047          ENDDO
     4048       ENDDO
     4049       !$acc end kernels
     4050
     4051!++
     4052!       sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
     4053
     4054    END SUBROUTINE advec_v_ws_acc
     4055   
     4056   
     4057!------------------------------------------------------------------------------!
    31834058! Advection of w - Call for all grid points
    31844059!------------------------------------------------------------------------------!
     
    35844459    END SUBROUTINE advec_w_ws
    35854460
     4461
     4462!------------------------------------------------------------------------------!
     4463! Advection of w - Call for all grid points - accelerator version
     4464!------------------------------------------------------------------------------!
     4465    SUBROUTINE advec_w_ws_acc
     4466
     4467       USE arrays_3d
     4468       USE constants
     4469       USE control_parameters
     4470       USE grid_variables
     4471       USE indices
     4472       USE statistics
     4473
     4474       IMPLICIT NONE
     4475
     4476       INTEGER ::  i, ibit27, ibit28, ibit29, ibit30, ibit31, ibit32, ibit33, &
     4477                   ibit34, ibit35, j, k, k_mmm, k_mm, k_pp, k_ppp, tn = 0
     4478
     4479       REAL    ::  diss_d, diss_l, diss_n, diss_r, diss_s, diss_t, div,    &
     4480                   flux_d, flux_l, flux_n, flux_r, flux_s, flux_t, gu, gv, &
     4481                   u_comp, u_comp_l, v_comp, v_comp_s, w_comp
     4482
     4483       gu = 2.0 * u_gtrans
     4484       gv = 2.0 * v_gtrans
     4485
     4486!
     4487!--    Computation of fluxes and tendency terms
     4488       !$acc kernels present( ddzu, tend, u, v, w, wall_flags_0 )
     4489       !$acc loop
     4490       DO i = nxl, nxr
     4491          DO  j = nys, nyn
     4492             !$acc loop vector( 32 )
     4493             DO  k = nzb+1, nzt
     4494
     4495                ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
     4496                ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
     4497                ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
     4498
     4499
     4500                u_comp_l                 = u(k+1,j,i) + u(k,j,i) - gu
     4501                flux_l                   = u_comp_l * (                        &
     4502                                      ( 37.0 * ibit29 * adv_mom_5              &
     4503                                   +     7.0 * ibit28 * adv_mom_3              &
     4504                                   +           ibit27 * adv_mom_1              &
     4505                                      ) *                                      &
     4506                                     ( w(k,j,i)   + w(k,j,i-1) )               &
     4507                               -      (  8.0 * ibit29 * adv_mom_5              &
     4508                                   +           ibit28 * adv_mom_3              &
     4509                                      ) *                                      &
     4510                                     ( w(k,j,i+1) + w(k,j,i-2) )               &
     4511                               +      (        ibit29 * adv_mom_5              &
     4512                                      ) *                                      &
     4513                                     ( w(k,j,i+2) + w(k,j,i-3) )               &
     4514                                                 )
     4515
     4516                diss_l                    = - ABS( u_comp_l ) * (              &
     4517                                        ( 10.0 * ibit29 * adv_mom_5            &
     4518                                     +     3.0 * ibit28 * adv_mom_3            &
     4519                                     +           ibit27 * adv_mom_1            &
     4520                                        ) *                                    &
     4521                                     ( w(k,j,i)   - w(k,j,i-1) )               &
     4522                                 -      (  5.0 * ibit29 * adv_mom_5            &
     4523                                     +           ibit28 * adv_mom_3            &
     4524                                        ) *                                    &
     4525                                     ( w(k,j,i+1) - w(k,j,i-2) )               &
     4526                                 +      (        ibit29 * adv_mom_5            &
     4527                                        ) *                                    &
     4528                                     ( w(k,j,i+2) - w(k,j,i-3) )               &
     4529                                                            )
     4530
     4531                u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
     4532                flux_r    = u_comp * (                                       &
     4533                          ( 37.0 * ibit29 * adv_mom_5                        &
     4534                       +     7.0 * ibit28 * adv_mom_3                        &
     4535                       +           ibit27 * adv_mom_1                        &
     4536                          ) *                                                &
     4537                                 ( w(k,j,i+1) + w(k,j,i)   )                 &
     4538                   -      (  8.0 * ibit29 * adv_mom_5                        &
     4539                       +           ibit28 * adv_mom_3                        &
     4540                          ) *                                                &
     4541                                 ( w(k,j,i+2) + w(k,j,i-1) )                 &
     4542                   +      (        ibit29 * adv_mom_5                        &
     4543                          ) *                                                &
     4544                                 ( w(k,j,i+3) + w(k,j,i-2) )                 &
     4545                                     )
     4546
     4547                diss_r    = - ABS( u_comp ) * (                              &
     4548                          ( 10.0 * ibit29 * adv_mom_5                        &
     4549                       +     3.0 * ibit28 * adv_mom_3                        &
     4550                       +           ibit27 * adv_mom_1                        &
     4551                          ) *                                                &
     4552                                 ( w(k,j,i+1) - w(k,j,i)  )                  &
     4553                   -      (  5.0 * ibit29 * adv_mom_5                        &
     4554                       +           ibit28 * adv_mom_3                        &
     4555                          ) *                                                &
     4556                                 ( w(k,j,i+2) - w(k,j,i-1) )                 &
     4557                   +      (        ibit29 * adv_mom_5                        &
     4558                          ) *                                                &
     4559                                 ( w(k,j,i+3) - w(k,j,i-2) )                 &
     4560                                              )
     4561
     4562                ibit32 = IBITS(wall_flags_0(k,j,i),32,1)
     4563                ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
     4564                ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
     4565
     4566                v_comp_s               = v(k+1,j,i) + v(k,j,i) - gv
     4567                flux_s                 = v_comp_s * (                         &
     4568                                    ( 37.0 * ibit32 * adv_mom_5               &
     4569                                 +     7.0 * ibit31 * adv_mom_3               &
     4570                                 +           ibit30 * adv_mom_1               &
     4571                                    ) *                                       &
     4572                                     ( w(k,j,i)   + w(k,j-1,i) )              &
     4573                             -      (  8.0 * ibit32 * adv_mom_5               &
     4574                                 +           ibit31 * adv_mom_3               &
     4575                                    ) *                                       &
     4576                                     ( w(k,j+1,i) + w(k,j-2,i) )              &
     4577                             +      (        ibit32 * adv_mom_5               &
     4578                                    ) *                                       &
     4579                                     ( w(k,j+2,i) + w(k,j-3,i) )              &
     4580                                               )
     4581
     4582                diss_s                 = - ABS( v_comp_s ) * (                &
     4583                                    ( 10.0 * ibit32 * adv_mom_5               &
     4584                                 +     3.0 * ibit31 * adv_mom_3               &
     4585                                 +           ibit30 * adv_mom_1               &
     4586                                    ) *                                       &
     4587                                     ( w(k,j,i)   - w(k,j-1,i) )              &
     4588                             -      (  5.0 * ibit32 * adv_mom_5               &
     4589                                 +           ibit31 * adv_mom_3               &
     4590                                    ) *                                       &
     4591                                     ( w(k,j+1,i) - w(k,j-2,i) )              &
     4592                             +      (        ibit32 * adv_mom_5               &
     4593                                    ) *                                       &
     4594                                     ( w(k,j+2,i) - w(k,j-3,i) )              &
     4595                                                        )
     4596
     4597                v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
     4598                flux_n    = v_comp * (                                       &
     4599                          ( 37.0 * ibit32 * adv_mom_5                        &
     4600                       +     7.0 * ibit31 * adv_mom_3                        &
     4601                       +           ibit30 * adv_mom_1                        &
     4602                          ) *                                                &
     4603                                 ( w(k,j+1,i) + w(k,j,i)   )                 &
     4604                   -      (  8.0 * ibit32 * adv_mom_5                        &
     4605                       +           ibit31 * adv_mom_3                        &
     4606                          ) *                                                &
     4607                                 ( w(k,j+2,i) + w(k,j-1,i) )                 &
     4608                   +      (        ibit32 * adv_mom_5                        &
     4609                          ) *                                                &
     4610                                 ( w(k,j+3,i) + w(k,j-2,i) )                 &
     4611                                     )
     4612
     4613                diss_n    = - ABS( v_comp ) * (                              &
     4614                          ( 10.0 * ibit32 * adv_mom_5                        &
     4615                       +     3.0 * ibit31 * adv_mom_3                        &
     4616                       +           ibit30 * adv_mom_1                        &
     4617                          ) *                                                &
     4618                                 ( w(k,j+1,i) - w(k,j,i)  )                  &
     4619                   -      (  5.0 * ibit32 * adv_mom_5                        &
     4620                       +           ibit31 * adv_mom_3                        &
     4621                          ) *                                                &
     4622                                 ( w(k,j+2,i) - w(k,j-1,i) )                 &
     4623                   +      (        ibit32 * adv_mom_5                        &
     4624                          ) *                                                &
     4625                                 ( w(k,j+3,i) - w(k,j-2,i) )                 &
     4626                                              )
     4627
     4628
     4629                ibit35 = IBITS(wall_flags_0(k-1,j,i),35,1)
     4630                ibit34 = IBITS(wall_flags_0(k-1,j,i),34,1)
     4631                ibit33 = IBITS(wall_flags_0(k-1,j,i),33,1)
     4632
     4633                k_pp  = k + 2 * ( 1 - ibit33 )
     4634                k_mm  = k - 2 * ( 1 - ibit33 )
     4635                k_mmm = k - 3 * ibit35
     4636
     4637                w_comp    = w(k,j,i) + w(k-1,j,i)
     4638                flux_d    = w_comp  * (                                      &
     4639                          ( 37.0 * ibit35 * adv_mom_5                        &
     4640                       +     7.0 * ibit34 * adv_mom_3                        &
     4641                       +           ibit33 * adv_mom_1                        &
     4642                          ) *                                                &
     4643                             ( w(k,j,i)    + w(k-1,j,i)   )                  &
     4644                   -      (  8.0 * ibit35 * adv_mom_5                        &
     4645                       +           ibit34 * adv_mom_3                        &
     4646                          ) *                                                &
     4647                             ( w(k+1,j,i)  + w(k_mm,j,i)  )                  &
     4648                   +      (        ibit35 * adv_mom_5                        &
     4649                          ) *                                                &
     4650                             ( w(k_pp,j,i) + w(k_mmm,j,i) )                  &
     4651                                       )
     4652
     4653                diss_d    = - ABS( w_comp ) * (                              &
     4654                          ( 10.0 * ibit35 * adv_mom_5                        &
     4655                       +     3.0 * ibit34 * adv_mom_3                        &
     4656                       +           ibit33 * adv_mom_1                        &
     4657                          ) *                                                &
     4658                             ( w(k,j,i)    - w(k-1,j,i)   )                  &
     4659                   -      (  5.0 * ibit35 * adv_mom_5                        &
     4660                       +           ibit34 * adv_mom_3                        &
     4661                          ) *                                                &
     4662                             ( w(k+1,j,i)  - w(k_mm,j,i)  )                  &
     4663                   +      (        ibit35 * adv_mom_5                        &
     4664                          ) *                                                &
     4665                             ( w(k_pp,j,i) - w(k_mmm,j,i) )                  &
     4666                                               )
     4667
     4668!
     4669!--             k index has to be modified near bottom and top, else array
     4670!--             subscripts will be exceeded.
     4671                ibit35 = IBITS(wall_flags_0(k,j,i),35,1)
     4672                ibit34 = IBITS(wall_flags_0(k,j,i),34,1)
     4673                ibit33 = IBITS(wall_flags_0(k,j,i),33,1)
     4674
     4675                k_ppp = k + 3 * ibit35
     4676                k_pp  = k + 2 * ( 1 - ibit33  )
     4677                k_mm  = k - 2 * ibit35
     4678
     4679                w_comp    = w(k+1,j,i) + w(k,j,i)
     4680                flux_t    = w_comp  * (                                      &
     4681                          ( 37.0 * ibit35 * adv_mom_5                        &
     4682                       +     7.0 * ibit34 * adv_mom_3                        &
     4683                       +           ibit33 * adv_mom_1                        &
     4684                          ) *                                                &
     4685                             ( w(k+1,j,i)  + w(k,j,i)     )                  &
     4686                   -      (  8.0 * ibit35 * adv_mom_5                        &
     4687                       +           ibit34 * adv_mom_3                        &
     4688                          ) *                                                &
     4689                             ( w(k_pp,j,i)  + w(k-1,j,i)  )                  &
     4690                   +      (        ibit35 * adv_mom_5                        &
     4691                          ) *                                                &
     4692                             ( w(k_ppp,j,i) + w(k_mm,j,i) )                  &
     4693                                       )
     4694
     4695                diss_t    = - ABS( w_comp ) * (                              &
     4696                          ( 10.0 * ibit35 * adv_mom_5                        &
     4697                       +     3.0 * ibit34 * adv_mom_3                        &
     4698                       +           ibit33 * adv_mom_1                        &
     4699                          ) *                                                &
     4700                             ( w(k+1,j,i)   - w(k,j,i)    )                  &
     4701                   -      (  5.0 * ibit35 * adv_mom_5                        &
     4702                       +           ibit34 * adv_mom_3                        &
     4703                          ) *                                                &
     4704                             ( w(k_pp,j,i)  - w(k-1,j,i)  )                  &
     4705                   +      (        ibit35 * adv_mom_5                        &
     4706                          ) *                                                &
     4707                             ( w(k_ppp,j,i) - w(k_mm,j,i) )                  &
     4708                                               )
     4709!
     4710!--             Calculate the divergence of the velocity field. A respective
     4711!--             correction is needed to overcome numerical instabilities caused
     4712!--             by a not sufficient reduction of divergences near topography.
     4713                div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx  &
     4714                    +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy  &
     4715                    +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) )        &
     4716                                                                 * ddzu(k+1) &
     4717                      ) * 0.5
     4718
     4719                tend(k,j,i) = - (                                                &
     4720                               ( flux_r + diss_r - flux_l - diss_l ) * ddx       &
     4721                             + ( flux_n + diss_n - flux_s - diss_s ) * ddy       &
     4722                             + ( flux_t + diss_t - flux_d - diss_d ) * ddzu(k+1) &
     4723                                 ) + div * w(k,j,i)
     4724
     4725
     4726!++
     4727!--             Statistical Evaluation of w'w'.
     4728!                sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                    &
     4729!                               + ( flux_t + diss_t )                    &
     4730!                               *   weight_substep(intermediate_timestep_count)
     4731
     4732             ENDDO
     4733          ENDDO
     4734       ENDDO
     4735       !$acc end kernels
     4736
     4737    END SUBROUTINE advec_w_ws_acc
     4738
    35864739 END MODULE advec_ws
Note: See TracChangeset for help on using the changeset viewer.