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

all OpenACC directives and related parts removed from the code

File:
1 edited

Legend:

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

    r2101 r2118  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! OpenACC version of subroutine removed
    2323!
    2424! Former revisions:
     
    216216!>       are zero at the walls and inside buildings.
    217217!------------------------------------------------------------------------------!
    218 #if ! defined( __openacc )
    219218 SUBROUTINE flow_statistics
    220219 
     
    309308    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
    310309
    311     !$acc update host( km, kh, e, ol, pt, qs, qsws, shf, ts, u, usws, v, vsws, w )
    312310
    313311!
     
    17561754
    17571755 END SUBROUTINE flow_statistics
    1758 
    1759 
    1760 #else
    1761 
    1762 
    1763 !------------------------------------------------------------------------------!
    1764 ! Description:
    1765 ! ------------
    1766 !> flow statistics - accelerator version
    1767 !------------------------------------------------------------------------------!
    1768  SUBROUTINE flow_statistics
    1769 
    1770     USE arrays_3d,                                                             &
    1771         ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
    1772                momentumflux_output_conversion, nr, p, prho, pt, q, qc, ql, qr, &
    1773                qs, qsws, qswst, rho_air, rho_air_zw, rho_ocean, s, sa, saswsb, &
    1774                saswst, shf, ss, ssws, sswst, td_lsa_lpt, td_lsa_q, td_sub_lpt, &
    1775                td_sub_q, time_vert, ts, tswst, u, ug, us, usws, uswst, vsws,   &
    1776                v, vg, vpt, vswst, w, w_subs, waterflux_output_conversion, zw
    1777                  
    1778        
    1779     USE cloud_parameters,                                                      &
    1780         ONLY:  l_d_cp, prr, pt_d_t
    1781        
    1782     USE control_parameters,                                                    &
    1783         ONLY :  average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
    1784                 dt_3d, g, humidity, kappa, large_scale_forcing,                &
    1785                 large_scale_subsidence, max_pr_user, message_string,           &
    1786                 microphysics_seifert, neutral, ocean, passive_scalar,          &
    1787                 simulated_time, use_subsidence_tendencies, use_surface_fluxes, &
    1788                 use_top_fluxes, ws_scheme_mom, ws_scheme_sca
    1789        
    1790     USE cpulog,                                                                &
    1791         ONLY:  cpu_log, log_point
    1792        
    1793     USE grid_variables,                                                        &
    1794         ONLY:  ddx, ddy
    1795        
    1796     USE indices,                                                               &
    1797         ONLY:  ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,       &
    1798                ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner,         &
    1799                nzb_s_inner, nzt, nzt_diff, rflags_invers
    1800        
    1801     USE kinds
    1802    
    1803     USE land_surface_model_mod,                                                &
    1804         ONLY:   ghf_eb, land_surface, m_soil, nzb_soil, nzt_soil,              &
    1805                 qsws_eb, qsws_liq_eb, qsws_soil_eb, qsws_veg_eb, r_a, r_s,     &
    1806                 shf_eb, t_soil
    1807 
    1808     USE netcdf_interface,                                                      &
    1809         ONLY:  dots_rad, dots_soil
    1810 
    1811     USE pegrid
    1812    
    1813     USE radiation_model_mod,                                                   &
    1814         ONLY:  radiation, radiation_scheme, rad_net,                 &
    1815                rad_lw_in, rad_lw_out, rad_sw_in, rad_sw_out
    1816 
    1817 #if defined ( __rrtmg )
    1818     USE radiation_model_mod,                                                   &
    1819         ONLY:  rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rad_lw_cs_hr,   &
    1820                rad_lw_hr,  rad_sw_cs_hr, rad_sw_hr
    1821 #endif
    1822 
    1823     USE statistics
    1824 
    1825     IMPLICIT NONE
    1826 
    1827     INTEGER(iwp) ::  i                   !<
    1828     INTEGER(iwp) ::  j                   !<
    1829     INTEGER(iwp) ::  k                   !<
    1830     INTEGER(iwp) ::  k_surface_level     !<
    1831     INTEGER(iwp) ::  nt                  !<
    1832     INTEGER(iwp) ::  omp_get_thread_num  !<
    1833     INTEGER(iwp) ::  sr                  !<
    1834     INTEGER(iwp) ::  tn                  !<
    1835    
    1836     LOGICAL ::  first  !<
    1837    
    1838     REAL(wp) ::  dptdz_threshold  !<
    1839     REAL(wp) ::  fac              !<
    1840     REAL(wp) ::  height           !<
    1841     REAL(wp) ::  pts              !<
    1842     REAL(wp) ::  sums_l_eper      !<
    1843     REAL(wp) ::  sums_l_etot      !<
    1844     REAL(wp) ::  s1               !<
    1845     REAL(wp) ::  s2               !<
    1846     REAL(wp) ::  s3               !<
    1847     REAL(wp) ::  s4               !<
    1848     REAL(wp) ::  s5               !<
    1849     REAL(wp) ::  s6               !<
    1850     REAL(wp) ::  s7               !<
    1851     REAL(wp) ::  ust              !<
    1852     REAL(wp) ::  ust2             !<
    1853     REAL(wp) ::  u2               !<
    1854     REAL(wp) ::  vst              !<
    1855     REAL(wp) ::  vst2             !<
    1856     REAL(wp) ::  v2               !<
    1857     REAL(wp) ::  w2               !<
    1858     REAL(wp) ::  z_i(2)           !<
    1859 
    1860     REAL(wp) ::  dptdz(nzb+1:nzt+1)    !<
    1861     REAL(wp) ::  sums_ll(nzb:nzt+1,2)  !<
    1862 
    1863     CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
    1864 
    1865 !
    1866 !-- To be on the safe side, check whether flow_statistics has already been
    1867 !-- called once after the current time step
    1868     IF ( flow_statistics_called )  THEN
    1869 
    1870        message_string = 'flow_statistics is called two times within one ' // &
    1871                         'timestep'
    1872        CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
    1873 
    1874     ENDIF
    1875 
    1876     !$acc data create( sums, sums_l )
    1877     !$acc update device( hom )
    1878 
    1879 !
    1880 !-- Compute statistics for each (sub-)region
    1881     DO  sr = 0, statistic_regions
    1882 
    1883 !
    1884 !--    Initialize (local) summation array
    1885        sums_l = 0.0_wp
    1886 
    1887 !
    1888 !--    Store sums that have been computed in other subroutines in summation
    1889 !--    array
    1890        sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
    1891 !--    WARNING: next line still has to be adjusted for OpenMP
    1892        sums_l(:,21,0) = sums_wsts_bc_l(:,sr) *                                 &
    1893                         heatflux_output_conversion  ! heat flux from advec_s_bc
    1894        sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
    1895        sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
    1896 
    1897 !
    1898 !--    When calcuating horizontally-averaged total (resolved- plus subgrid-
    1899 !--    scale) vertical fluxes and velocity variances by using commonly-
    1900 !--    applied Reynolds-based methods ( e.g. <w'pt'> = (w-<w>)*(pt-<pt>) )
    1901 !--    in combination with the 5th order advection scheme, pronounced
    1902 !--    artificial kinks could be observed in the vertical profiles near the
    1903 !--    surface. Please note: these kinks were not related to the model truth,
    1904 !--    i.e. these kinks are just related to an evaluation problem.   
    1905 !--    In order avoid these kinks, vertical fluxes and horizontal as well
    1906 !--    vertical velocity variances are calculated directly within the advection
    1907 !--    routines, according to the numerical discretization, to evaluate the
    1908 !--    statistical quantities as they will appear within the prognostic
    1909 !--    equations.
    1910 !--    Copy the turbulent quantities, evaluated in the advection routines to
    1911 !--    the local array sums_l() for further computations.
    1912        IF ( ws_scheme_mom .AND. sr == 0 )  THEN
    1913 
    1914 !
    1915 !--       According to the Neumann bc for the horizontal velocity components,
    1916 !--       the corresponding fluxes has to satisfiy the same bc.
    1917           IF ( ocean )  THEN
    1918              sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
    1919              sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
    1920           ENDIF
    1921 
    1922           DO  i = 0, threads_per_task-1
    1923 !
    1924 !--          Swap the turbulent quantities evaluated in advec_ws.
    1925              sums_l(:,13,i) = sums_wsus_ws_l(:,i)                              &
    1926                               * momentumflux_output_conversion ! w*u*
    1927              sums_l(:,15,i) = sums_wsvs_ws_l(:,i)                              &
    1928                               * momentumflux_output_conversion ! w*v*
    1929              sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
    1930              sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
    1931              sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
    1932              sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        &
    1933                               ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
    1934                                 sums_ws2_ws_l(:,i) )    ! e*
    1935              DO  k = nzb, nzt
    1936                 sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5_wp * ( &
    1937                                                       sums_us2_ws_l(k,i) +     &
    1938                                                       sums_vs2_ws_l(k,i) +     &
    1939                                                       sums_ws2_ws_l(k,i)     )
    1940              ENDDO
    1941           ENDDO
    1942 
    1943        ENDIF
    1944 
    1945        IF ( ws_scheme_sca .AND. sr == 0 )  THEN
    1946 
    1947           DO  i = 0, threads_per_task-1
    1948              sums_l(:,17,i) = sums_wspts_ws_l(:,i)                             &
    1949                               * heatflux_output_conversion        ! w*pt* from advec_s_ws
    1950              IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
    1951              IF ( humidity       )  sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)      &
    1952                                             * waterflux_output_conversion !w*q*
    1953              IF ( passive_scalar )  sums_l(:,116,i) = sums_wsss_ws_l(:,i) !w*s*
    1954           ENDDO
    1955 
    1956        ENDIF
    1957 !
    1958 !--    Horizontally averaged profiles of horizontal velocities and temperature.
    1959 !--    They must have been computed before, because they are already required
    1960 !--    for other horizontal averages.
    1961        tn = 0
    1962 
    1963        !$OMP PARALLEL PRIVATE( i, j, k, tn )
    1964 !$     tn = omp_get_thread_num()
    1965 
    1966        !$acc update device( sums_l )
    1967 
    1968        !$OMP DO
    1969        !$acc parallel loop gang present( pt, rflags_invers, rmask, sums_l, u, v ) create( s1, s2, s3 )
    1970        DO  k = nzb, nzt+1
    1971           s1 = 0
    1972           s2 = 0
    1973           s3 = 0
    1974           !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
    1975           DO  i = nxl, nxr
    1976              DO  j =  nys, nyn
    1977 !
    1978 !--             k+1 is used in rflags since rflags is set 0 at surface points
    1979                 s1 = s1 + u(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    1980                 s2 = s2 + v(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    1981                 s3 = s3 + pt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    1982              ENDDO
    1983           ENDDO
    1984           sums_l(k,1,tn) = s1
    1985           sums_l(k,2,tn) = s2
    1986           sums_l(k,4,tn) = s3
    1987        ENDDO
    1988        !$acc end parallel loop
    1989 
    1990 !
    1991 !--    Horizontally averaged profile of salinity
    1992        IF ( ocean )  THEN
    1993           !$OMP DO
    1994           !$acc parallel loop gang present( rflags_invers, rmask, sums_l, sa ) create( s1 )
    1995           DO  k = nzb, nzt+1
    1996              s1 = 0
    1997              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    1998              DO  i = nxl, nxr
    1999                 DO  j =  nys, nyn
    2000                    s1 = s1 + sa(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2001                 ENDDO
    2002              ENDDO
    2003              sums_l(k,23,tn) = s1
    2004           ENDDO
    2005           !$acc end parallel loop
    2006        ENDIF
    2007 
    2008 !
    2009 !--    Horizontally averaged profiles of virtual potential temperature,
    2010 !--    total water content, specific humidity and liquid water potential
    2011 !--    temperature
    2012        IF ( humidity )  THEN
    2013 
    2014           !$OMP DO
    2015           !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l, vpt ) create( s1, s2 )
    2016           DO  k = nzb, nzt+1
    2017              s1 = 0
    2018              s2 = 0
    2019              !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
    2020              DO  i = nxl, nxr
    2021                 DO  j =  nys, nyn
    2022                    s1 = s1 + q(k,j,i)   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2023                    s2 = s2 + vpt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2024                 ENDDO
    2025              ENDDO
    2026              sums_l(k,41,tn) = s1
    2027              sums_l(k,44,tn) = s2
    2028           ENDDO
    2029           !$acc end parallel loop
    2030 
    2031           IF ( cloud_physics )  THEN
    2032              !$OMP DO
    2033              !$acc parallel loop gang present( pt, q, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
    2034              DO  k = nzb, nzt+1
    2035                 s1 = 0
    2036                 s2 = 0
    2037                 !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
    2038                 DO  i = nxl, nxr
    2039                    DO  j =  nys, nyn
    2040                       s1 = s1 + ( q(k,j,i) - ql(k,j,i) ) * &
    2041                                 rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2042                       s2 = s2 + ( pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) ) * &
    2043                                 rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2044                    ENDDO
    2045                 ENDDO
    2046                 sums_l(k,42,tn) = s1
    2047                 sums_l(k,43,tn) = s2
    2048              ENDDO
    2049              !$acc end parallel loop
    2050           ENDIF
    2051        ENDIF
    2052 
    2053 !
    2054 !--    Horizontally averaged profiles of passive scalar
    2055        IF ( passive_scalar )  THEN
    2056           !$OMP DO
    2057           !$acc parallel loop gang present( s, rflags_invers, rmask, sums_l ) create( s1 )
    2058           DO  k = nzb, nzt+1
    2059              s1 = 0
    2060              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2061              DO  i = nxl, nxr
    2062                 DO  j =  nys, nyn
    2063                    s1 = s1 + s(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2064                 ENDDO
    2065              ENDDO
    2066              sums_l(k,117,tn) = s1
    2067           ENDDO
    2068           !$acc end parallel loop
    2069        ENDIF
    2070        !$OMP END PARALLEL
    2071 
    2072 !
    2073 !--    Summation of thread sums
    2074        IF ( threads_per_task > 1 )  THEN
    2075           DO  i = 1, threads_per_task-1
    2076              !$acc parallel present( sums_l )
    2077              sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
    2078              sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
    2079              sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
    2080              !$acc end parallel
    2081              IF ( ocean )  THEN
    2082                 !$acc parallel present( sums_l )
    2083                 sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
    2084                 !$acc end parallel
    2085              ENDIF
    2086              IF ( humidity )  THEN
    2087                 !$acc parallel present( sums_l )
    2088                 sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
    2089                 sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
    2090                 !$acc end parallel
    2091                 IF ( cloud_physics )  THEN
    2092                    !$acc parallel present( sums_l )
    2093                    sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
    2094                    sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
    2095                    !$acc end parallel
    2096                 ENDIF
    2097              ENDIF
    2098              IF ( passive_scalar )  THEN
    2099                 !$acc parallel present( sums_l )
    2100                 sums_l(:,117,0) = sums_l(:,117,0) + sums_l(:,117,i)
    2101                 !$acc end parallel
    2102              ENDIF
    2103           ENDDO
    2104        ENDIF
    2105 
    2106 #if defined( __parallel )
    2107 !
    2108 !--    Compute total sum from local sums
    2109        !$acc update host( sums_l )
    2110        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2111        CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL,  &
    2112                            MPI_SUM, comm2d, ierr )
    2113        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2114        CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL,  &
    2115                            MPI_SUM, comm2d, ierr )
    2116        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2117        CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL,  &
    2118                            MPI_SUM, comm2d, ierr )
    2119        IF ( ocean )  THEN
    2120           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2121           CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb,       &
    2122                               MPI_REAL, MPI_SUM, comm2d, ierr )
    2123        ENDIF
    2124        IF ( humidity ) THEN
    2125           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2126           CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb,       &
    2127                               MPI_REAL, MPI_SUM, comm2d, ierr )
    2128           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2129           CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
    2130                               MPI_REAL, MPI_SUM, comm2d, ierr )
    2131           IF ( cloud_physics ) THEN
    2132              IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2133              CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb,    &
    2134                                  MPI_REAL, MPI_SUM, comm2d, ierr )
    2135              IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2136              CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb,    &
    2137                                  MPI_REAL, MPI_SUM, comm2d, ierr )
    2138           ENDIF
    2139        ENDIF
    2140 
    2141        IF ( passive_scalar )  THEN
    2142           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2143           CALL MPI_ALLREDUCE( sums_l(nzb,117,0), sums(nzb,117), nzt+2-nzb,     &
    2144                               MPI_REAL, MPI_SUM, comm2d, ierr )
    2145        ENDIF
    2146        !$acc update device( sums )
    2147 #else
    2148        !$acc parallel present( sums, sums_l )
    2149        sums(:,1) = sums_l(:,1,0)
    2150        sums(:,2) = sums_l(:,2,0)
    2151        sums(:,4) = sums_l(:,4,0)
    2152        !$acc end parallel
    2153        IF ( ocean )  THEN
    2154           !$acc parallel present( sums, sums_l )
    2155           sums(:,23) = sums_l(:,23,0)
    2156           !$acc end parallel
    2157        ENDIF
    2158        IF ( humidity )  THEN
    2159           !$acc parallel present( sums, sums_l )
    2160           sums(:,44) = sums_l(:,44,0)
    2161           sums(:,41) = sums_l(:,41,0)
    2162           !$acc end parallel
    2163           IF ( cloud_physics )  THEN
    2164              !$acc parallel present( sums, sums_l )
    2165              sums(:,42) = sums_l(:,42,0)
    2166              sums(:,43) = sums_l(:,43,0)
    2167              !$acc end parallel
    2168           ENDIF
    2169        ENDIF
    2170        IF ( passive_scalar )  THEN
    2171           !$acc parallel present( sums, sums_l )
    2172           sums(:,117) = sums_l(:,117,0)
    2173           !$acc end parallel
    2174        ENDIF
    2175 #endif
    2176 
    2177 !
    2178 !--    Final values are obtained by division by the total number of grid points
    2179 !--    used for summation. After that store profiles.
    2180        !$acc parallel present( hom, ngp_2dh, ngp_2dh_s_inner, sums )
    2181        sums(:,1) = sums(:,1) / ngp_2dh(sr)
    2182        sums(:,2) = sums(:,2) / ngp_2dh(sr)
    2183        sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
    2184        hom(:,1,1,sr) = sums(:,1)             ! u
    2185        hom(:,1,2,sr) = sums(:,2)             ! v
    2186        hom(:,1,4,sr) = sums(:,4)             ! pt
    2187        !$acc end parallel
    2188 
    2189 !
    2190 !--    Salinity
    2191        IF ( ocean )  THEN
    2192           !$acc parallel present( hom, ngp_2dh_s_inner, sums )
    2193           sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
    2194           hom(:,1,23,sr) = sums(:,23)             ! sa
    2195           !$acc end parallel
    2196        ENDIF
    2197 
    2198 !
    2199 !--    Humidity and cloud parameters
    2200        IF ( humidity ) THEN
    2201           !$acc parallel present( hom, ngp_2dh_s_inner, sums )
    2202           sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
    2203           sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
    2204           hom(:,1,44,sr) = sums(:,44)                ! vpt
    2205           hom(:,1,41,sr) = sums(:,41)                ! qv (q)
    2206           !$acc end parallel
    2207           IF ( cloud_physics ) THEN
    2208              !$acc parallel present( hom, ngp_2dh_s_inner, sums )
    2209              sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
    2210              sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
    2211              hom(:,1,42,sr) = sums(:,42)             ! qv
    2212              hom(:,1,43,sr) = sums(:,43)             ! pt
    2213              !$acc end parallel
    2214           ENDIF
    2215        ENDIF
    2216 
    2217 !
    2218 !--    Passive scalar
    2219        IF ( passive_scalar )  THEN
    2220           !$acc parallel present( hom, ngp_2dh_s_inner, sums )
    2221           sums(:,117)     = sums(:,117) / ngp_2dh_s_inner(:,sr)
    2222           hom(:,1,117,sr) = sums(:,117)                ! s
    2223           !$acc end parallel
    2224        ENDIF
    2225 
    2226 !
    2227 !--    Horizontally averaged profiles of the remaining prognostic variables,
    2228 !--    variances, the total and the perturbation energy (single values in last
    2229 !--    column of sums_l) and some diagnostic quantities.
    2230 !--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
    2231 !--    ----  speaking the following k-loop would have to be split up and
    2232 !--          rearranged according to the staggered grid.
    2233 !--          However, this implies no error since staggered velocity components
    2234 !--          are zero at the walls and inside buildings.
    2235        tn = 0
    2236        !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper,             &
    2237        !$OMP                   sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2,  &
    2238        !$OMP                   w2 )
    2239 !$     tn = omp_get_thread_num()
    2240 
    2241        !$OMP DO
    2242        !$acc parallel loop gang present( e, hom, kh, km, p, pt, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, s5, s6, s7 )
    2243        DO  k = nzb, nzt+1
    2244           s1 = 0
    2245           s2 = 0
    2246           s3 = 0
    2247           s4 = 0
    2248           s5 = 0
    2249           s6 = 0
    2250           s7 = 0
    2251           !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5, s6, s7 )
    2252           DO  i = nxl, nxr
    2253              DO  j =  nys, nyn
    2254 !
    2255 !--             Prognostic and diagnostic variables
    2256                 s1 = s1 + w(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2257                 s2 = s2 + e(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2258                 s3 = s3 + km(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2259                 s4 = s4 + kh(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2260                 s5 = s5 + p(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2261                 s6 = s6 + ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr) * &
    2262                           rflags_invers(j,i,k+1)
    2263 !
    2264 !--             Higher moments
    2265 !--             (Computation of the skewness of w further below)
    2266                 s7 = s7 + w(k,j,i)**3 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2267              ENDDO
    2268           ENDDO
    2269           sums_l(k,3,tn)  = s1
    2270           sums_l(k,8,tn)  = s2
    2271           sums_l(k,9,tn)  = s3
    2272           sums_l(k,10,tn) = s4
    2273           sums_l(k,40,tn) = s5
    2274           sums_l(k,33,tn) = s6
    2275           sums_l(k,38,tn) = s7
    2276        ENDDO
    2277        !$acc end parallel loop
    2278 
    2279        IF ( humidity )  THEN
    2280           !$OMP DO
    2281           !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l ) create( s1 )
    2282           DO  k = nzb, nzt+1
    2283              s1 = 0
    2284              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2285              DO  i = nxl, nxr
    2286                 DO  j =  nys, nyn
    2287                    s1 = s1 + ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr) * &
    2288                              rflags_invers(j,i,k+1)
    2289                 ENDDO
    2290              ENDDO
    2291              sums_l(k,70,tn) = s1
    2292           ENDDO
    2293           !$acc end parallel loop
    2294        ENDIF
    2295 
    2296 !
    2297 !--    Total and perturbation energy for the total domain (being
    2298 !--    collected in the last column of sums_l).
    2299        s1 = 0
    2300        !$OMP DO
    2301        !$acc parallel loop collapse(3) present( rflags_invers, rmask, u, v, w ) reduction(+:s1)
    2302        DO  i = nxl, nxr
    2303           DO  j =  nys, nyn
    2304              DO  k = nzb, nzt+1
    2305                 s1 = s1 + 0.5_wp *                                             &
    2306                           ( u(k,j,i)**2 + v(k,j,i)**2 + w(k,j,i)**2 ) *        &
    2307                           rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2308              ENDDO
    2309           ENDDO
    2310        ENDDO
    2311        !$acc end parallel loop
    2312        !$acc parallel present( sums_l )
    2313        sums_l(nzb+4,pr_palm,tn) = s1
    2314        !$acc end parallel
    2315 
    2316        !$OMP DO
    2317        !$acc parallel present( rmask, sums_l, us, usws, vsws, ts ) create( s1, s2, s3, s4 )
    2318        s1 = 0
    2319        s2 = 0
    2320        s3 = 0
    2321        s4 = 0
    2322        !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
    2323        DO  i = nxl, nxr
    2324           DO  j =  nys, nyn
    2325 !
    2326 !--          2D-arrays (being collected in the last column of sums_l)
    2327              s1 = s1 + us(j,i)   * rmask(j,i,sr)
    2328              s2 = s2 + usws(j,i) * rmask(j,i,sr)
    2329              s3 = s3 + vsws(j,i) * rmask(j,i,sr)
    2330              s4 = s4 + ts(j,i)   * rmask(j,i,sr)
    2331           ENDDO
    2332        ENDDO
    2333        sums_l(nzb,pr_palm,tn)   = s1
    2334        sums_l(nzb+1,pr_palm,tn) = s2
    2335        sums_l(nzb+2,pr_palm,tn) = s3
    2336        sums_l(nzb+3,pr_palm,tn) = s4
    2337        !$acc end parallel
    2338 
    2339        IF ( humidity )  THEN
    2340           !$acc parallel present( qs, rmask, sums_l ) create( s1 )
    2341           s1 = 0
    2342           !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2343           DO  i = nxl, nxr
    2344              DO  j =  nys, nyn
    2345                 s1 = s1 + qs(j,i) * rmask(j,i,sr)
    2346              ENDDO
    2347           ENDDO
    2348           sums_l(nzb+12,pr_palm,tn) = s1
    2349           !$acc end parallel
    2350        ENDIF
    2351 
    2352        IF ( passive_scalar )  THEN
    2353           !$acc parallel present( ss, rmask, sums_l ) create( s1 )
    2354           s1 = 0
    2355           !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2356           DO  i = nxl, nxr
    2357              DO  j =  nys, nyn
    2358                 s1 = s1 + ss(j,i) * rmask(j,i,sr)
    2359              ENDDO
    2360           ENDDO
    2361           sums_l(nzb+13,pr_palm,tn) = s1
    2362           !$acc end parallel
    2363        ENDIF
    2364 
    2365 !
    2366 !--    Computation of statistics when ws-scheme is not used. Else these
    2367 !--    quantities are evaluated in the advection routines.
    2368        IF ( .NOT. ws_scheme_mom .OR. sr /= 0 .OR. simulated_time == 0.0_wp )   &
    2369        THEN
    2370 
    2371           !$OMP DO
    2372           !$acc parallel loop gang present( u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, ust2, vst2, w2 )
    2373           DO  k = nzb, nzt+1
    2374              s1 = 0
    2375              s2 = 0
    2376              s3 = 0
    2377              s4 = 0
    2378              !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
    2379              DO  i = nxl, nxr
    2380                 DO  j =  nys, nyn
    2381                    ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
    2382                    vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
    2383                    w2   = w(k,j,i)**2
    2384 
    2385                    s1 = s1 + ust2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2386                    s2 = s2 + vst2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2387                    s3 = s3 + w2   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2388 !
    2389 !--                Perturbation energy
    2390                    s4 = s4 + 0.5_wp * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) *   &
    2391                              rflags_invers(j,i,k+1)
    2392                 ENDDO
    2393              ENDDO
    2394              sums_l(k,30,tn) = s1
    2395              sums_l(k,31,tn) = s2
    2396              sums_l(k,32,tn) = s3
    2397              sums_l(k,34,tn) = s4
    2398           ENDDO
    2399           !$acc end parallel loop
    2400 !
    2401 !--       Total perturbation TKE
    2402           !$OMP DO
    2403           !$acc parallel present( sums_l ) create( s1 )
    2404           s1 = 0
    2405           !$acc loop reduction( +: s1 )
    2406           DO  k = nzb, nzt+1
    2407              s1 = s1 + sums_l(k,34,tn)
    2408           ENDDO
    2409           sums_l(nzb+5,pr_palm,tn) = s1
    2410           !$acc end parallel
    2411 
    2412        ENDIF
    2413 
    2414 !
    2415 !--    Horizontally averaged profiles of the vertical fluxes
    2416 
    2417 !
    2418 !--    Subgridscale fluxes.
    2419 !--    WARNING: If a Prandtl-layer is used (k=nzb for flat terrain), the fluxes
    2420 !--    -------  should be calculated there in a different way. This is done
    2421 !--             in the next loop further below, where results from this loop are
    2422 !--             overwritten. However, THIS WORKS IN CASE OF FLAT TERRAIN ONLY!
    2423 !--             The non-flat case still has to be handled.
    2424 !--    NOTE: for simplicity, nzb_s_inner is used below, although
    2425 !--    ----  strictly speaking the following k-loop would have to be
    2426 !--          split up according to the staggered grid.
    2427 !--          However, this implies no error since staggered velocity
    2428 !--          components are zero at the walls and inside buildings.
    2429        !$OMP DO
    2430        !$acc parallel loop gang present( ddzu, kh, km, pt, u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
    2431        DO  k = nzb, nzt_diff
    2432           s1 = 0
    2433           s2 = 0
    2434           s3 = 0
    2435           !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
    2436           DO  i = nxl, nxr
    2437              DO  j = nys, nyn
    2438 
    2439 !
    2440 !--             Momentum flux w"u"
    2441                 s1 = s1 - 0.25_wp * (                                          &
    2442                                km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
    2443                                                            ) * (               &
    2444                                    ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
    2445                                  + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
    2446                                                                )               &
    2447                                * rmask(j,i,sr) * rflags_invers(j,i,k+1)        &
    2448                                * rho_air_zw(k)                                 &
    2449                                * momentumflux_output_conversion(k)
    2450 !
    2451 !--             Momentum flux w"v"
    2452                 s2 = s2 - 0.25_wp * (                                          &
    2453                                km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
    2454                                                            ) * (               &
    2455                                    ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
    2456                                  + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
    2457                                                                )               &
    2458                                * rmask(j,i,sr) * rflags_invers(j,i,k+1)        &
    2459                                * rho_air_zw(k)                                 &
    2460                                * momentumflux_output_conversion(k)
    2461 !
    2462 !--             Heat flux w"pt"
    2463                 s3 = s3 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                 &
    2464                                  * ( pt(k+1,j,i) - pt(k,j,i) )                 &
    2465                                  * rho_air_zw(k)                               &
    2466                                  * heatflux_output_conversion(k)               &
    2467                                  * ddzu(k+1) * rmask(j,i,sr)                   &
    2468                                  * rflags_invers(j,i,k+1)
    2469              ENDDO
    2470           ENDDO
    2471           sums_l(k,12,tn) = s1
    2472           sums_l(k,14,tn) = s2
    2473           sums_l(k,16,tn) = s3
    2474        ENDDO
    2475        !$acc end parallel loop
    2476 
    2477 !
    2478 !--    Salinity flux w"sa"
    2479        IF ( ocean )  THEN
    2480           !$acc parallel loop gang present( ddzu, kh, sa, rflags_invers, rmask, sums_l ) create( s1 )
    2481           DO  k = nzb, nzt_diff
    2482              s1 = 0
    2483              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2484              DO  i = nxl, nxr
    2485                 DO  j = nys, nyn
    2486                    s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
    2487                                     * ( sa(k+1,j,i) - sa(k,j,i) )              &
    2488                                     * ddzu(k+1) * rmask(j,i,sr)                &
    2489                                     * rflags_invers(j,i,k+1)
    2490                 ENDDO
    2491              ENDDO
    2492              sums_l(k,65,tn) = s1
    2493           ENDDO
    2494           !$acc end parallel loop
    2495        ENDIF
    2496 
    2497 !
    2498 !--    Buoyancy flux, water flux (humidity flux) w"q"
    2499        IF ( humidity ) THEN
    2500 
    2501           !$acc parallel loop gang present( ddzu, kh, q, vpt, rflags_invers, rmask, sums_l ) create( s1, s2 )
    2502           DO  k = nzb, nzt_diff
    2503              s1 = 0
    2504              s2 = 0
    2505              !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
    2506              DO  i = nxl, nxr
    2507                 DO  j = nys, nyn
    2508                    s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
    2509                                     * ( vpt(k+1,j,i) - vpt(k,j,i) )            &
    2510                                     * rho_air_zw(k)                            &
    2511                                     * heatflux_output_conversion(k)            &
    2512                                     * ddzu(k+1) * rmask(j,i,sr)                &
    2513                                     * rflags_invers(j,i,k+1)
    2514                    s2 = s2 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
    2515                                     * ( q(k+1,j,i) - q(k,j,i) )                &
    2516                                     * rho_air_zw(k)                            &
    2517                                     * waterflux_output_conversion(k)           &
    2518                                     * ddzu(k+1) * rmask(j,i,sr)                &
    2519                                     * rflags_invers(j,i,k+1)
    2520                 ENDDO
    2521              ENDDO
    2522              sums_l(k,45,tn) = s1
    2523              sums_l(k,48,tn) = s2
    2524           ENDDO
    2525           !$acc end parallel loop
    2526 
    2527           IF ( cloud_physics ) THEN
    2528 
    2529              !$acc parallel loop gang present( ddzu, kh, q, ql, rflags_invers, rmask, sums_l ) create( s1 )
    2530              DO  k = nzb, nzt_diff
    2531                 s1 = 0
    2532                 !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2533                 DO  i = nxl, nxr
    2534                    DO  j = nys, nyn
    2535                       s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )           &
    2536                                        *  ( ( q(k+1,j,i) - ql(k+1,j,i) )       &
    2537                                           - ( q(k,j,i) - ql(k,j,i) ) )         &
    2538                                        * rho_air_zw(k)                         &
    2539                                        * waterflux_output_conversion(k)        &
    2540                                        * ddzu(k+1) * rmask(j,i,sr)             &
    2541                                        * rflags_invers(j,i,k+1)
    2542                    ENDDO
    2543                 ENDDO
    2544                 sums_l(k,51,tn) = s1
    2545              ENDDO
    2546              !$acc end parallel loop
    2547 
    2548           ENDIF
    2549 
    2550        ENDIF
    2551 !
    2552 !--    Passive scalar flux
    2553        IF ( passive_scalar )  THEN
    2554 
    2555           !$acc parallel loop gang present( ddzu, kh, s, rflags_invers, rmask, sums_l ) create( s1 )
    2556           DO  k = nzb, nzt_diff
    2557              s1 = 0
    2558              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2559              DO  i = nxl, nxr
    2560                 DO  j = nys, nyn
    2561                    s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
    2562                                     * ( s(k+1,j,i) - s(k,j,i) )                &
    2563                                     * ddzu(k+1) * rmask(j,i,sr)                &
    2564                                     * rflags_invers(j,i,k+1)
    2565                 ENDDO
    2566              ENDDO
    2567              sums_l(k,119,tn) = s1
    2568           ENDDO
    2569           !$acc end parallel loop
    2570 
    2571        ENDIF
    2572 
    2573        IF ( use_surface_fluxes )  THEN
    2574 
    2575           !$OMP DO
    2576           !$acc parallel present( rmask, shf, sums_l, usws, vsws ) create( s1, s2, s3, s4, s5 )
    2577           s1 = 0
    2578           s2 = 0
    2579           s3 = 0
    2580           s4 = 0
    2581           s5 = 0
    2582           !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
    2583           DO  i = nxl, nxr
    2584              DO  j =  nys, nyn
    2585 !
    2586 !--             Subgridscale fluxes in the Prandtl layer
    2587                 s1 = s1 + usws(j,i) * momentumflux_output_conversion(nzb)      &
    2588                                     * rmask(j,i,sr) ! w"u"
    2589                 s2 = s2 + vsws(j,i) * momentumflux_output_conversion(nzb)      &
    2590                                     * rmask(j,i,sr) ! w"v"
    2591                 s3 = s3 + shf(j,i)  * heatflux_output_conversion(nzb)          &
    2592                                     * rmask(j,i,sr) ! w"pt"
    2593                 s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
    2594                 s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
    2595              ENDDO
    2596           ENDDO
    2597           sums_l(nzb,12,tn) = s1
    2598           sums_l(nzb,14,tn) = s2
    2599           sums_l(nzb,16,tn) = s3
    2600           sums_l(nzb,58,tn) = s4
    2601           sums_l(nzb,61,tn) = s5
    2602           !$acc end parallel
    2603 
    2604           IF ( ocean )  THEN
    2605 
    2606              !$OMP DO
    2607              !$acc parallel present( rmask, saswsb, sums_l ) create( s1 )
    2608              s1 = 0
    2609              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2610              DO  i = nxl, nxr
    2611                 DO  j =  nys, nyn
    2612                    s1 = s1 + saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
    2613                 ENDDO
    2614              ENDDO
    2615              sums_l(nzb,65,tn) = s1
    2616              !$acc end parallel
    2617 
    2618           ENDIF
    2619 
    2620           IF ( humidity )  THEN
    2621 
    2622              !$OMP DO
    2623              !$acc parallel present( pt, q, qsws, rmask, shf, sums_l ) create( s1, s2 )
    2624              s1 = 0
    2625              s2 = 0
    2626              !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
    2627              DO  i = nxl, nxr
    2628                 DO  j =  nys, nyn
    2629                    s1 = s1 + qsws(j,i) * waterflux_output_conversion(nzb)      &
    2630                                        * rmask(j,i,sr) ! w"q" (w"qv")
    2631                    s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * shf(j,i)    &
    2632                                + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )           &
    2633                              * heatflux_output_conversion(nzb)
    2634                 ENDDO
    2635              ENDDO
    2636              sums_l(nzb,48,tn) = s1
    2637              sums_l(nzb,45,tn) = s2
    2638              !$acc end parallel
    2639 
    2640              IF ( cloud_droplets )  THEN
    2641 
    2642                 !$OMP DO
    2643                 !$acc parallel present( pt, q, ql, qsws, rmask, shf, sums_l ) create( s1 )
    2644                 s1 = 0
    2645                 !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2646                 DO  i = nxl, nxr
    2647                    DO  j =  nys, nyn
    2648                       s1 = s1 + ( ( 1.0_wp +                                   &
    2649                                     0.61_wp * q(nzb,j,i) - ql(nzb,j,i) ) *     &
    2650                                  shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )&
    2651                                 * heatflux_output_conversion(nzb)
    2652                    ENDDO
    2653                 ENDDO
    2654                 sums_l(nzb,45,tn) = s1
    2655                 !$acc end parallel
    2656 
    2657              ENDIF
    2658 
    2659              IF ( cloud_physics )  THEN
    2660 
    2661                 !$OMP DO
    2662                 !$acc parallel present( qsws, rmask, sums_l ) create( s1 )
    2663                 s1 = 0
    2664                 !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2665                 DO  i = nxl, nxr
    2666                    DO  j =  nys, nyn
    2667 !
    2668 !--                   Formula does not work if ql(nzb) /= 0.0
    2669                       s1 = s1 + qsws(j,i) * waterflux_output_conversion(nzb)   &
    2670                                           * rmask(j,i,sr)   ! w"q" (w"qv")
    2671                    ENDDO
    2672                 ENDDO
    2673                 sums_l(nzb,51,tn) = s1
    2674                 !$acc end parallel
    2675 
    2676              ENDIF
    2677 
    2678           ENDIF
    2679 
    2680           IF ( passive_scalar )  THEN
    2681 
    2682              !$OMP DO
    2683              !$acc parallel present( ssws, rmask, sums_l ) create( s1 )
    2684              s1 = 0
    2685              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2686              DO  i = nxl, nxr
    2687                 DO  j =  nys, nyn
    2688                    s1 = s1 + ssws(j,i) * rmask(j,i,sr)  ! w"s"
    2689                 ENDDO
    2690              ENDDO
    2691              sums_l(nzb,119,tn) = s1
    2692              !$acc end parallel
    2693 
    2694           ENDIF
    2695 
    2696        ENDIF
    2697 
    2698 !
    2699 !--    Subgridscale fluxes at the top surface
    2700        IF ( use_top_fluxes )  THEN
    2701 
    2702           !$OMP DO
    2703           !$acc parallel present( rmask, sums_l, tswst, uswst, vswst ) create( s1, s2, s3, s4, s5 )
    2704           s1 = 0
    2705           s2 = 0
    2706           s3 = 0
    2707           s4 = 0
    2708           s5 = 0
    2709           !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
    2710           DO  i = nxl, nxr
    2711              DO  j =  nys, nyn
    2712                 s1 = s1 + uswst(j,i) * momentumflux_output_conversion(nzt:nzt+1) &
    2713                                      * rmask(j,i,sr)    ! w"u"
    2714                 s2 = s2 + vswst(j,i) * momentumflux_output_conversion(nzt:nzt+1) &
    2715                                      * rmask(j,i,sr)    ! w"v"
    2716                 s3 = s3 + tswst(j,i) * heatflux_output_conversion(nzt:nzt+1)   &
    2717                                      * rmask(j,i,sr)    ! w"pt"
    2718                 s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
    2719                 s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
    2720              ENDDO
    2721           ENDDO
    2722           sums_l(nzt:nzt+1,12,tn) = s1
    2723           sums_l(nzt:nzt+1,14,tn) = s2
    2724           sums_l(nzt:nzt+1,16,tn) = s3
    2725           sums_l(nzt:nzt+1,58,tn) = s4
    2726           sums_l(nzt:nzt+1,61,tn) = s5
    2727           !$acc end parallel
    2728 
    2729           IF ( ocean )  THEN
    2730 
    2731              !$OMP DO
    2732              !$acc parallel present( rmask, saswst, sums_l ) create( s1 )
    2733              s1 = 0
    2734              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2735              DO  i = nxl, nxr
    2736                 DO  j =  nys, nyn
    2737                    s1 = s1 + saswst(j,i) * rmask(j,i,sr)  ! w"sa"
    2738                 ENDDO
    2739              ENDDO
    2740              sums_l(nzt,65,tn) = s1
    2741              !$acc end parallel
    2742 
    2743           ENDIF
    2744 
    2745           IF ( humidity )  THEN
    2746 
    2747              !$OMP DO
    2748              !$acc parallel present( pt, q, qswst, rmask, tswst, sums_l ) create( s1, s2 )
    2749              s1 = 0
    2750              s2 = 0
    2751              !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
    2752              DO  i = nxl, nxr
    2753                 DO  j =  nys, nyn
    2754                    s1 = s1 + qswst(j,i) * waterflux_output_conversion(nzt)     &
    2755                                         * rmask(j,i,sr) ! w"q" (w"qv")
    2756                    s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * tswst(j,i) +&
    2757                                  0.61_wp * pt(nzt,j,i) * qswst(j,i) )          &
    2758                              * heatflux_output_conversion(nzt)
    2759                 ENDDO
    2760              ENDDO
    2761              sums_l(nzt,48,tn) = s1
    2762              sums_l(nzt,45,tn) = s2
    2763              !$acc end parallel
    2764 
    2765              IF ( cloud_droplets )  THEN
    2766 
    2767                 !$OMP DO
    2768                 !$acc parallel present( pt, q, ql, qswst, rmask, tswst, sums_l ) create( s1 )
    2769                 s1 = 0
    2770                 !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2771                 DO  i = nxl, nxr
    2772                    DO  j =  nys, nyn
    2773                       s1 = s1 + ( ( 1.0_wp +                                   &
    2774                                     0.61_wp * q(nzt,j,i) - ql(nzt,j,i) ) *     &
    2775                                   tswst(j,i) +                                 &
    2776                                   0.61_wp * pt(nzt,j,i) * qswst(j,i) )         &
    2777                                 * heatflux_output_conversion(nzt)
    2778                    ENDDO
    2779                 ENDDO
    2780                 sums_l(nzt,45,tn) = s1
    2781                 !$acc end parallel
    2782 
    2783              ENDIF
    2784 
    2785              IF ( cloud_physics )  THEN
    2786 
    2787                 !$OMP DO
    2788                 !$acc parallel present( qswst, rmask, sums_l ) create( s1 )
    2789                 s1 = 0
    2790                 !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2791                 DO  i = nxl, nxr
    2792                    DO  j =  nys, nyn
    2793 !
    2794 !--                   Formula does not work if ql(nzb) /= 0.0
    2795                       s1 = s1 + qswst(j,i) * waterflux_output_conversion(nzt)  &
    2796                                            * rmask(j,i,sr)  ! w"q" (w"qv")
    2797                    ENDDO
    2798                 ENDDO
    2799                 sums_l(nzt,51,tn) = s1
    2800                 !$acc end parallel
    2801 
    2802              ENDIF
    2803 
    2804           ENDIF
    2805 
    2806           IF ( passive_scalar )  THEN
    2807 
    2808              !$OMP DO
    2809              !$acc parallel present( sswst, rmask, sums_l ) create( s1 )
    2810              s1 = 0
    2811              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2812              DO  i = nxl, nxr
    2813                 DO  j =  nys, nyn
    2814                    s1 = s1 + sswst(j,i) * rmask(j,i,sr) ! w"s"
    2815                 ENDDO
    2816              ENDDO
    2817              sums_l(nzt,119,tn) = s1
    2818              !$acc end parallel
    2819 
    2820           ENDIF
    2821 
    2822        ENDIF
    2823 
    2824 !
    2825 !--    Resolved fluxes (can be computed for all horizontal points)
    2826 !--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
    2827 !--    ----  speaking the following k-loop would have to be split up and
    2828 !--          rearranged according to the staggered grid.
    2829        !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2, s3 )
    2830        DO  k = nzb, nzt_diff
    2831           s1 = 0
    2832           s2 = 0
    2833           s3 = 0
    2834           !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
    2835           DO  i = nxl, nxr
    2836              DO  j = nys, nyn
    2837                 ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) + &
    2838                                  u(k+1,j,i) - hom(k+1,1,1,sr) )
    2839                 vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) + &
    2840                                  v(k+1,j,i) - hom(k+1,1,2,sr) )
    2841                 pts = 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) + &
    2842                                  pt(k+1,j,i) - hom(k+1,1,4,sr) )
    2843 !
    2844 !--             Higher moments
    2845                 s1 = s1 + pts * w(k,j,i)**2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2846                 s2 = s2 + pts**2 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2847 !
    2848 !--             Energy flux w*e* (has to be adjusted?)
    2849                 s3 = s3 + w(k,j,i) * 0.5_wp * ( ust**2 + vst**2 + w(k,j,i)**2 )&
    2850                                    * rmask(j,i,sr) * rflags_invers(j,i,k+1)    &
    2851                                    * momentumflux_output_conversion(k)
    2852              ENDDO
    2853           ENDDO
    2854           sums_l(k,35,tn) = s1
    2855           sums_l(k,36,tn) = s2
    2856           sums_l(k,37,tn) = s3
    2857        ENDDO
    2858        !$acc end parallel loop
    2859 
    2860 !
    2861 !--    Salinity flux and density (density does not belong to here,
    2862 !--    but so far there is no other suitable place to calculate)
    2863        IF ( ocean )  THEN
    2864 
    2865           IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
    2866 
    2867              !$acc parallel loop gang present( hom, rflags_invers, rmask, sa, sums_l, w ) create( s1 )
    2868              DO  k = nzb, nzt_diff
    2869                 s1 = 0
    2870                 !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2871                 DO  i = nxl, nxr
    2872                    DO  j = nys, nyn
    2873                       s1 = s1 + 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +      &
    2874                                            sa(k+1,j,i) - hom(k+1,1,23,sr) )    &
    2875                                        * w(k,j,i) * rmask(j,i,sr)              &
    2876                                        * rflags_invers(j,i,k+1)
    2877                    ENDDO
    2878                 ENDDO
    2879                 sums_l(k,66,tn) = s1
    2880              ENDDO
    2881              !$acc end parallel loop
    2882 
    2883           ENDIF
    2884 
    2885           !$acc parallel loop gang present( rflags_invers, rho_ocean, prho, rmask, sums_l ) create( s1, s2 )
    2886           DO  k = nzb, nzt_diff
    2887              s1 = 0
    2888              s2 = 0
    2889              !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
    2890              DO  i = nxl, nxr
    2891                 DO  j = nys, nyn
    2892                    s1 = s1 + rho_ocean(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2893                    s2 = s2 + prho(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2894                 ENDDO
    2895              ENDDO
    2896              sums_l(k,64,tn) = s1
    2897              sums_l(k,71,tn) = s2
    2898           ENDDO
    2899           !$acc end parallel loop
    2900 
    2901        ENDIF
    2902 
    2903 !
    2904 !--    Buoyancy flux, water flux, humidity flux, liquid water
    2905 !--    content, rain drop concentration and rain water content
    2906        IF ( humidity )  THEN
    2907 
    2908           IF ( cloud_physics  .OR.  cloud_droplets )  THEN
    2909 
    2910              !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
    2911              DO  k = nzb, nzt_diff
    2912                 s1 = 0
    2913                 !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2914                 DO  i = nxl, nxr
    2915                    DO  j = nys, nyn
    2916                       s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +     &
    2917                                            vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * &
    2918                                          heatflux_output_conversion(k) *       &
    2919                                          w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2920                    ENDDO
    2921                 ENDDO
    2922                 sums_l(k,46,tn) = s1
    2923              ENDDO
    2924              !$acc end parallel loop
    2925 
    2926              IF ( .NOT. cloud_droplets )  THEN
    2927 
    2928                 !$acc parallel loop gang present( hom, q, ql, rflags_invers, rmask, sums_l, w ) create( s1 )
    2929                 DO  k = nzb, nzt_diff
    2930                    s1 = 0
    2931                    !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2932                    DO  i = nxl, nxr
    2933                       DO  j = nys, nyn
    2934                          s1 = s1 + 0.5_wp * ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) +   &
    2935                                               ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) &
    2936                                           * waterflux_output_conversion(k)                      &
    2937                                           * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2938                       ENDDO
    2939                    ENDDO
    2940                    sums_l(k,52,tn) = s1
    2941                 ENDDO
    2942                 !$acc end parallel loop
    2943 
    2944                 IF ( microphysics_seifert )  THEN
    2945 
    2946                    !$acc parallel loop gang present( qc, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
    2947                    DO  k = nzb, nzt_diff
    2948                       s1 = 0
    2949                       !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
    2950                       DO  i = nxl, nxr
    2951                          DO  j = nys, nyn
    2952                             s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2953                             s2 = s2 + qc(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2954                          ENDDO
    2955                       ENDDO
    2956                       sums_l(k,54,tn) = s1
    2957                       sums_l(k,75,tn) = s2
    2958                    ENDDO
    2959                    !$acc end parallel loop
    2960 
    2961                    !$acc parallel loop gang present( nr, qr, prr, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
    2962                    DO  k = nzb, nzt_diff
    2963                       s1 = 0
    2964                       s2 = 0
    2965                       s3 = 0
    2966                       !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
    2967                       DO  i = nxl, nxr
    2968                          DO  j = nys, nyn
    2969                             s1 = s1 + nr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2970                             s2 = s2 + qr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2971                             s3 = s3 + prr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2972                          ENDDO
    2973                       ENDDO
    2974                       sums_l(k,73,tn) = s1
    2975                       sums_l(k,74,tn) = s2
    2976                       sums_l(k,76,tn) = s3
    2977                    ENDDO
    2978                    !$acc end parallel loop
    2979 
    2980                 ELSE
    2981 
    2982                    !$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
    2983                    DO  k = nzb, nzt_diff
    2984                       s1 = 0
    2985                       !$acc loop vector collapse( 2 ) reduction( +: s1 )
    2986                       DO  i = nxl, nxr
    2987                          DO  j = nys, nyn
    2988                             s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    2989                          ENDDO
    2990                       ENDDO
    2991                       sums_l(k,54,tn) = s1
    2992                    ENDDO
    2993                    !$acc end parallel loop
    2994 
    2995                 ENDIF
    2996 
    2997              ELSE
    2998 
    2999                 !$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
    3000                 DO  k = nzb, nzt_diff
    3001                    s1 = 0
    3002                    !$acc loop vector collapse( 2 ) reduction( +: s1 )
    3003                    DO  i = nxl, nxr
    3004                       DO  j = nys, nyn
    3005                          s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    3006                       ENDDO
    3007                    ENDDO
    3008                    sums_l(k,54,tn) = s1
    3009                 ENDDO
    3010                 !$acc end parallel loop
    3011 
    3012              ENDIF
    3013 
    3014           ELSE
    3015 
    3016              IF( .NOT. ws_scheme_sca  .OR.  sr /= 0 )  THEN
    3017 
    3018                 !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
    3019                 DO  k = nzb, nzt_diff
    3020                    s1 = 0
    3021                    !$acc loop vector collapse( 2 ) reduction( +: s1 )
    3022                    DO  i = nxl, nxr
    3023                       DO  j = nys, nyn
    3024                          s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +   &
    3025                                               vpt(k+1,j,i) - hom(k+1,1,44,sr) ) &
    3026                                           * heatflux_output_conversion(k)       &
    3027                                           * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    3028                       ENDDO
    3029                    ENDDO
    3030                    sums_l(k,46,tn) = s1
    3031                 ENDDO
    3032                 !$acc end parallel loop
    3033 
    3034              ELSEIF ( ws_scheme_sca  .AND.  sr == 0 )  THEN
    3035 
    3036                 !$acc parallel loop present( hom, sums_l )
    3037                 DO  k = nzb, nzt_diff
    3038                    sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp * hom(k,1,41,sr) ) * &
    3039                                        sums_l(k,17,tn) + 0.61_wp *             &
    3040                                        hom(k,1,4,sr) * sums_l(k,49,tn)         &
    3041                                      ) * heatflux_output_conversion(k)
    3042                 ENDDO
    3043                 !$acc end parallel loop
    3044 
    3045              ENDIF
    3046 
    3047           ENDIF
    3048 
    3049        ENDIF
    3050 !
    3051 !--    Passive scalar flux
    3052        IF ( passive_scalar  .AND.  ( .NOT. ws_scheme_sca  .OR.  sr /= 0 ) )  THEN
    3053 
    3054           !$acc parallel loop gang present( hom, s, rflags_invers, rmask, sums_l, w ) create( s1 )
    3055           DO  k = nzb, nzt_diff
    3056              s1 = 0
    3057              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    3058              DO  i = nxl, nxr
    3059                 DO  j = nys, nyn
    3060                    s1 = s1 + 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +          &
    3061                                         s(k+1,j,i) - hom(k+1,1,117,sr) )        &
    3062                                     * w(k,j,i) * rmask(j,i,sr)                 &
    3063                                     * rflags_invers(j,i,k+1)
    3064                 ENDDO
    3065              ENDDO
    3066              sums_l(k,49,tn) = s1
    3067           ENDDO
    3068           !$acc end parallel loop
    3069 
    3070        ENDIF
    3071 
    3072 !
    3073 !--    For speed optimization fluxes which have been computed in part directly
    3074 !--    inside the WS advection routines are treated seperatly
    3075 !--    Momentum fluxes first:
    3076        IF ( .NOT. ws_scheme_mom  .OR.  sr /= 0  )  THEN
    3077 
    3078           !$OMP DO
    3079           !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2 )
    3080           DO  k = nzb, nzt_diff
    3081              s1 = 0
    3082              s2 = 0
    3083              !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
    3084              DO  i = nxl, nxr
    3085                 DO  j = nys, nyn
    3086                    ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +               &
    3087                                     u(k+1,j,i) - hom(k+1,1,1,sr) )
    3088                    vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +               &
    3089                                     v(k+1,j,i) - hom(k+1,1,2,sr) )
    3090 !
    3091 !--                Momentum flux w*u*
    3092                    s1 = s1 + 0.5_wp * ( w(k,j,i-1) + w(k,j,i) )                &
    3093                                     * ust * rmask(j,i,sr)                      &
    3094                                     * momentumflux_output_conversion(k)        &
    3095                                     * rflags_invers(j,i,k+1)
    3096 !
    3097 !--                Momentum flux w*v*
    3098                    s2 = s2 + 0.5_wp * ( w(k,j-1,i) + w(k,j,i) )                &
    3099                                     * vst * rmask(j,i,sr)                      &
    3100                                     * momentumflux_output_conversion(k)        &
    3101                                     * rflags_invers(j,i,k+1)
    3102                 ENDDO
    3103              ENDDO
    3104              sums_l(k,13,tn) = s1
    3105              sums_l(k,15,tn) = s2
    3106           ENDDO
    3107           !$acc end parallel loop
    3108 
    3109        ENDIF
    3110 
    3111        IF ( .NOT. ws_scheme_sca  .OR.  sr /= 0 )  THEN
    3112 
    3113           !$OMP DO
    3114           !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, w ) create( s1 )
    3115           DO  k = nzb, nzt_diff
    3116              s1 = 0
    3117              !$acc loop vector collapse( 2 ) reduction( +: s1 )
    3118              DO  i = nxl, nxr
    3119                 DO  j = nys, nyn
    3120 !
    3121 !--                Vertical heat flux
    3122                    s1 = s1 + 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +          &
    3123                                         pt(k+1,j,i) - hom(k+1,1,4,sr) )        &
    3124                                     * heatflux_output_conversion(k)            &
    3125                                     * w(k,j,i) * rmask(j,i,sr)                 &
    3126                                     * rflags_invers(j,i,k+1)
    3127                 ENDDO
    3128              ENDDO
    3129              sums_l(k,17,tn) = s1
    3130           ENDDO
    3131           !$acc end parallel loop
    3132 
    3133           IF ( humidity )  THEN
    3134 
    3135              !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
    3136              DO  k = nzb, nzt_diff
    3137                 s1 = 0
    3138                 !$acc loop vector collapse( 2 ) reduction( +: s1 )
    3139                 DO  i = nxl, nxr
    3140                    DO  j = nys, nyn
    3141                       s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +       &
    3142                                            q(k+1,j,i) - hom(k+1,1,41,sr) )     &
    3143                                        * waterflux_output_conversion(k)        &
    3144                                        * w(k,j,i) * rmask(j,i,sr)              &
    3145                                        * rflags_invers(j,i,k+1)
    3146                    ENDDO
    3147                 ENDDO
    3148                 sums_l(k,49,tn) = s1
    3149              ENDDO
    3150              !$acc end parallel loop
    3151 
    3152           ENDIF
    3153 
    3154           IF ( passive_scalar )  THEN
    3155 
    3156              !$acc parallel loop gang present( hom, s, rflags_invers, rmask, sums_l, w ) create( s1 )
    3157              DO  k = nzb, nzt_diff
    3158                 s1 = 0
    3159                 !$acc loop vector collapse( 2 ) reduction( +: s1 )
    3160                 DO  i = nxl, nxr
    3161                    DO  j = nys, nyn
    3162                       s1 = s1 + 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +      &
    3163                                            s(k+1,j,i) - hom(k+1,1,117,sr) )    &
    3164                                        * w(k,j,i) * rmask(j,i,sr)              &
    3165                                        * rflags_invers(j,i,k+1)
    3166                    ENDDO
    3167                 ENDDO
    3168                 sums_l(k,116,tn) = s1
    3169              ENDDO
    3170              !$acc end parallel loop
    3171 
    3172           ENDIF
    3173 
    3174        ENDIF
    3175 
    3176 
    3177 !
    3178 !--    Density at top follows Neumann condition
    3179        IF ( ocean )  THEN
    3180           !$acc parallel present( sums_l )
    3181           sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
    3182           sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
    3183           !$acc end parallel
    3184        ENDIF
    3185 
    3186 !
    3187 !--    Divergence of vertical flux of resolved scale energy and pressure
    3188 !--    fluctuations as well as flux of pressure fluctuation itself (68).
    3189 !--    First calculate the products, then the divergence.
    3190 !--    Calculation is time consuming. Do it only, if profiles shall be plotted.
    3191        IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )  THEN
    3192 
    3193           STOP '+++ openACC porting for vertical flux div of resolved scale TKE in flow_statistics is still missing'
    3194           sums_ll = 0.0_wp  ! local array
    3195 
    3196           !$OMP DO
    3197           DO  i = nxl, nxr
    3198              DO  j = nys, nyn
    3199                 DO  k = nzb_s_inner(j,i)+1, nzt
    3200 
    3201                    sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
    3202                   ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1) )  &
    3203                             - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) ) )**2&
    3204                 + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) )  &
    3205                             - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) )**2&
    3206                 + w(k,j,i)**2                                        )
    3207 
    3208                    sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
    3209                                                * ( p(k,j,i) + p(k+1,j,i) )
    3210 
    3211                 ENDDO
    3212              ENDDO
    3213           ENDDO
    3214           sums_ll(0,1)     = 0.0_wp    ! because w is zero at the bottom
    3215           sums_ll(nzt+1,1) = 0.0_wp
    3216           sums_ll(0,2)     = 0.0_wp
    3217           sums_ll(nzt+1,2) = 0.0_wp
    3218 
    3219           DO  k = nzb+1, nzt
    3220              sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
    3221              sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
    3222              sums_l(k,68,tn) = sums_ll(k,2)
    3223           ENDDO
    3224           sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
    3225           sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
    3226           sums_l(nzb,68,tn) = 0.0_wp    ! because w* = 0 at nzb
    3227 
    3228        ENDIF
    3229 
    3230 !
    3231 !--    Divergence of vertical flux of SGS TKE and the flux itself (69)
    3232        IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )  THEN
    3233 
    3234           STOP '+++ openACC porting for vertical flux div of SGS TKE in flow_statistics is still missing'
    3235           !$OMP DO
    3236           DO  i = nxl, nxr
    3237              DO  j = nys, nyn
    3238                 DO  k = nzb_s_inner(j,i)+1, nzt
    3239 
    3240                    sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
    3241                    (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
    3242                  - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
    3243                                                                 ) * ddzw(k)
    3244 
    3245                    sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
    3246                    (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
    3247                                                                 )
    3248 
    3249                 ENDDO
    3250              ENDDO
    3251           ENDDO
    3252           sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
    3253           sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
    3254 
    3255        ENDIF
    3256 
    3257 !
    3258 !--    Horizontal heat fluxes (subgrid, resolved, total).
    3259 !--    Do it only, if profiles shall be plotted.
    3260        IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
    3261 
    3262           STOP '+++ openACC porting for horizontal flux calculation in flow_statistics is still missing'
    3263           !$OMP DO
    3264           DO  i = nxl, nxr
    3265              DO  j = nys, nyn
    3266                 DO  k = nzb_s_inner(j,i)+1, nzt
    3267 !
    3268 !--                Subgrid horizontal heat fluxes u"pt", v"pt"
    3269                    sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
    3270                                                    ( kh(k,j,i) + kh(k,j,i-1) ) &
    3271                                                  * ( pt(k,j,i-1) - pt(k,j,i) ) &
    3272                                                * rho_air_zw(k)                 &
    3273                                                * heatflux_output_conversion(k) &
    3274                                                  * ddx * rmask(j,i,sr)
    3275                    sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
    3276                                                    ( kh(k,j,i) + kh(k,j-1,i) ) &
    3277                                                  * ( pt(k,j-1,i) - pt(k,j,i) ) &
    3278                                                * rho_air_zw(k)                 &
    3279                                                * heatflux_output_conversion(k) &
    3280                                                  * ddy * rmask(j,i,sr)
    3281 !
    3282 !--                Resolved horizontal heat fluxes u*pt*, v*pt*
    3283                    sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
    3284                                      ( u(k,j,i) - hom(k,1,1,sr) ) * 0.5_wp *   &
    3285                                      ( pt(k,j,i-1) - hom(k,1,4,sr) +           &
    3286                                        pt(k,j,i)   - hom(k,1,4,sr) )           &
    3287                                      * heatflux_output_conversion(k)
    3288                    pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
    3289                                     pt(k,j,i)   - hom(k,1,4,sr) )
    3290                    sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
    3291                                      ( v(k,j,i) - hom(k,1,2,sr) ) * 0.5_wp *   &
    3292                                      ( pt(k,j-1,i) - hom(k,1,4,sr) +           &
    3293                                        pt(k,j,i)   - hom(k,1,4,sr) )           &
    3294                                      * heatflux_output_conversion(k)
    3295                 ENDDO
    3296              ENDDO
    3297           ENDDO
    3298 !
    3299 !--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
    3300           sums_l(nzb,58,tn) = 0.0_wp
    3301           sums_l(nzb,59,tn) = 0.0_wp
    3302           sums_l(nzb,60,tn) = 0.0_wp
    3303           sums_l(nzb,61,tn) = 0.0_wp
    3304           sums_l(nzb,62,tn) = 0.0_wp
    3305           sums_l(nzb,63,tn) = 0.0_wp
    3306 
    3307        ENDIF
    3308 
    3309 !
    3310 !--    Collect current large scale advection and subsidence tendencies for
    3311 !--    data output
    3312        IF ( large_scale_forcing  .AND.  ( simulated_time > 0.0_wp ) )  THEN
    3313 !
    3314 !--       Interpolation in time of LSF_DATA
    3315           nt = 1
    3316           DO WHILE ( simulated_time - dt_3d > time_vert(nt) )
    3317              nt = nt + 1
    3318           ENDDO
    3319           IF ( simulated_time - dt_3d /= time_vert(nt) )  THEN
    3320             nt = nt - 1
    3321           ENDIF
    3322 
    3323           fac = ( simulated_time - dt_3d - time_vert(nt) )                     &
    3324                 / ( time_vert(nt+1)-time_vert(nt) )
    3325 
    3326 
    3327           DO  k = nzb, nzt
    3328              sums_ls_l(k,0) = td_lsa_lpt(k,nt)                                 &
    3329                               + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
    3330              sums_ls_l(k,1) = td_lsa_q(k,nt)                                   &
    3331                               + fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
    3332           ENDDO
    3333 
    3334           sums_ls_l(nzt+1,0) = sums_ls_l(nzt,0)
    3335           sums_ls_l(nzt+1,1) = sums_ls_l(nzt,1)
    3336 
    3337           IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
    3338 
    3339              DO  k = nzb, nzt
    3340                 sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac *                      &
    3341                                  ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
    3342                 sums_ls_l(k,3) = td_sub_q(k,nt) + fac *                        &
    3343                                  ( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
    3344              ENDDO
    3345 
    3346              sums_ls_l(nzt+1,2) = sums_ls_l(nzt,2)
    3347              sums_ls_l(nzt+1,3) = sums_ls_l(nzt,3)
    3348 
    3349           ENDIF
    3350 
    3351        ENDIF
    3352 
    3353 
    3354        IF ( land_surface )  THEN
    3355           !$OMP DO
    3356           DO  i = nxl, nxr
    3357              DO  j =  nys, nyn
    3358                 DO  k = nzb_soil, nzt_soil
    3359                    sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i)         &
    3360                                       * rmask(j,i,sr)
    3361                    sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i)         &
    3362                                       * rmask(j,i,sr)
    3363                 ENDDO
    3364              ENDDO
    3365           ENDDO
    3366        ENDIF
    3367 
    3368 
    3369        IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
    3370           !$OMP DO
    3371           DO  i = nxl, nxr
    3372              DO  j =  nys, nyn
    3373                 DO  k = nzb_s_inner(j,i)+1, nzt+1
    3374                    sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i)    &
    3375                                        * rmask(j,i,sr)
    3376                    sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i)   &
    3377                                        * rmask(j,i,sr)
    3378                    sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i)    &
    3379                                        * rmask(j,i,sr)
    3380                    sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i)   &
    3381                                        * rmask(j,i,sr)
    3382 #if defined ( __rrtmg )
    3383                    sums_l(k,106,tn)  = sums_l(k,106,tn)  + rad_lw_cs_hr(k,j,i) &
    3384                                        * rmask(j,i,sr)
    3385                    sums_l(k,107,tn)  = sums_l(k,107,tn)  + rad_lw_hr(k,j,i)    &
    3386                                        * rmask(j,i,sr)
    3387                    sums_l(k,108,tn)  = sums_l(k,108,tn)  + rad_sw_cs_hr(k,j,i) &
    3388                                        * rmask(j,i,sr)
    3389                    sums_l(k,109,tn)  = sums_l(k,109,tn)  + rad_sw_hr(k,j,i)    &
    3390                                        * rmask(j,i,sr)
    3391 #endif
    3392                 ENDDO
    3393              ENDDO
    3394           ENDDO
    3395        ENDIF
    3396 
    3397 !
    3398 !--    Calculate the user-defined profiles
    3399        CALL user_statistics( 'profiles', sr, tn )
    3400        !$OMP END PARALLEL
    3401 
    3402 !
    3403 !--    Summation of thread sums
    3404        IF ( threads_per_task > 1 )  THEN
    3405           STOP '+++ openACC porting for threads_per_task > 1 in flow_statistics is still missing'
    3406           DO  i = 1, threads_per_task-1
    3407              sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
    3408              sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
    3409              sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
    3410                                       sums_l(:,45:pr_palm,i)
    3411              IF ( max_pr_user > 0 )  THEN
    3412                 sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
    3413                                    sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
    3414                                    sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
    3415              ENDIF
    3416           ENDDO
    3417        ENDIF
    3418 
    3419        !$acc update host( hom, sums, sums_l )
    3420 
    3421 #if defined( __parallel )
    3422 
    3423 !
    3424 !--    Compute total sum from local sums
    3425        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    3426        CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
    3427                            MPI_SUM, comm2d, ierr )
    3428        IF ( large_scale_forcing )  THEN
    3429           CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls,     &
    3430                               MPI_REAL, MPI_SUM, comm2d, ierr )
    3431        ENDIF
    3432 #else
    3433        sums = sums_l(:,:,0)
    3434        IF ( large_scale_forcing )  THEN
    3435           sums(:,81:88) = sums_ls_l
    3436        ENDIF
    3437 #endif
    3438 
    3439 !
    3440 !--    Final values are obtained by division by the total number of grid points
    3441 !--    used for summation. After that store profiles.
    3442 !--    Check, if statistical regions do contain at least one grid point at the
    3443 !--    respective k-level, otherwise division by zero will lead to undefined
    3444 !--    values, which may cause e.g. problems with NetCDF output
    3445 !--    Profiles:
    3446        DO  k = nzb, nzt+1
    3447           sums(k,3)             = sums(k,3)             / ngp_2dh(sr)
    3448           sums(k,12:22)         = sums(k,12:22)         / ngp_2dh(sr)
    3449           sums(k,30:32)         = sums(k,30:32)         / ngp_2dh(sr)
    3450           sums(k,35:39)         = sums(k,35:39)         / ngp_2dh(sr)
    3451           sums(k,45:53)         = sums(k,45:53)         / ngp_2dh(sr)
    3452           sums(k,55:63)         = sums(k,55:63)         / ngp_2dh(sr)
    3453           sums(k,81:88)         = sums(k,81:88)         / ngp_2dh(sr)
    3454           sums(k,89:114)        = sums(k,89:114)        / ngp_2dh(sr)
    3455           IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
    3456              sums(k,8:11)          = sums(k,8:11)          / ngp_2dh_s_inner(k,sr)
    3457              sums(k,23:29)         = sums(k,23:29)         / ngp_2dh_s_inner(k,sr)
    3458              sums(k,33:34)         = sums(k,33:34)         / ngp_2dh_s_inner(k,sr)
    3459              sums(k,40)            = sums(k,40)            / ngp_2dh_s_inner(k,sr)
    3460              sums(k,54)            = sums(k,54)            / ngp_2dh_s_inner(k,sr)
    3461              sums(k,64)            = sums(k,64)            / ngp_2dh_s_inner(k,sr)
    3462              sums(k,70:80)         = sums(k,70:80)         / ngp_2dh_s_inner(k,sr)
    3463              sums(k,115:pr_palm-2) = sums(k,115:pr_palm-2) / ngp_2dh_s_inner(k,sr)
    3464           ENDIF
    3465        ENDDO
    3466 
    3467 !--    u* and so on
    3468 !--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
    3469 !--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
    3470 !--    above the topography, they are being divided by ngp_2dh(sr)
    3471        sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
    3472                                     ngp_2dh(sr)
    3473        sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
    3474                                     ngp_2dh(sr)
    3475 !--    eges, e*
    3476        sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
    3477                                     ngp_3d(sr)
    3478 !--    Old and new divergence
    3479        sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
    3480                                     ngp_3d_inner(sr)
    3481 
    3482 !--    User-defined profiles
    3483        IF ( max_pr_user > 0 )  THEN
    3484           DO  k = nzb, nzt+1
    3485              IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
    3486                 sums(k,pr_palm+1:pr_palm+max_pr_user) = &
    3487                                        sums(k,pr_palm+1:pr_palm+max_pr_user) / &
    3488                                        ngp_2dh_s_inner(k,sr)
    3489              ENDIF
    3490           ENDDO
    3491        ENDIF
    3492 
    3493 !
    3494 !--    Collect horizontal average in hom.
    3495 !--    Compute deduced averages (e.g. total heat flux)
    3496        hom(:,1,3,sr)  = sums(:,3)      ! w
    3497        hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
    3498        hom(:,1,9,sr)  = sums(:,9)      ! km
    3499        hom(:,1,10,sr) = sums(:,10)     ! kh
    3500        hom(:,1,11,sr) = sums(:,11)     ! l
    3501        hom(:,1,12,sr) = sums(:,12)     ! w"u"
    3502        hom(:,1,13,sr) = sums(:,13)     ! w*u*
    3503        hom(:,1,14,sr) = sums(:,14)     ! w"v"
    3504        hom(:,1,15,sr) = sums(:,15)     ! w*v*
    3505        hom(:,1,16,sr) = sums(:,16)     ! w"pt"
    3506        hom(:,1,17,sr) = sums(:,17)     ! w*pt*
    3507        hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
    3508        hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
    3509        hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
    3510        hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
    3511        hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
    3512                                        ! profile 24 is initial profile (sa)
    3513                                        ! profiles 25-29 left empty for initial
    3514                                        ! profiles
    3515        hom(:,1,30,sr) = sums(:,30)     ! u*2
    3516        hom(:,1,31,sr) = sums(:,31)     ! v*2
    3517        hom(:,1,32,sr) = sums(:,32)     ! w*2
    3518        hom(:,1,33,sr) = sums(:,33)     ! pt*2
    3519        hom(:,1,34,sr) = sums(:,34)     ! e*
    3520        hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
    3521        hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
    3522        hom(:,1,37,sr) = sums(:,37)     ! w*e*
    3523        hom(:,1,38,sr) = sums(:,38)     ! w*3
    3524        hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
    3525        hom(:,1,40,sr) = sums(:,40)     ! p
    3526        hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
    3527        hom(:,1,46,sr) = sums(:,46)     ! w*vpt*
    3528        hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
    3529        hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
    3530        hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
    3531        hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
    3532        hom(:,1,51,sr) = sums(:,51)     ! w"qv"
    3533        hom(:,1,52,sr) = sums(:,52)     ! w*qv*
    3534        hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
    3535        hom(:,1,54,sr) = sums(:,54)     ! ql
    3536        hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
    3537        hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
    3538        hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho_ocean )/dz
    3539        hom(:,1,58,sr) = sums(:,58)     ! u"pt"
    3540        hom(:,1,59,sr) = sums(:,59)     ! u*pt*
    3541        hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
    3542        hom(:,1,61,sr) = sums(:,61)     ! v"pt"
    3543        hom(:,1,62,sr) = sums(:,62)     ! v*pt*
    3544        hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
    3545        hom(:,1,64,sr) = sums(:,64)     ! rho_ocean
    3546        hom(:,1,65,sr) = sums(:,65)     ! w"sa"
    3547        hom(:,1,66,sr) = sums(:,66)     ! w*sa*
    3548        hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
    3549        hom(:,1,68,sr) = sums(:,68)     ! w*p*
    3550        hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho_ocean
    3551        hom(:,1,70,sr) = sums(:,70)     ! q*2
    3552        hom(:,1,71,sr) = sums(:,71)     ! prho
    3553        hom(:,1,72,sr) = hyp * 1E-4_wp     ! hyp in dbar
    3554        hom(:,1,73,sr) = sums(:,73)     ! nr
    3555        hom(:,1,74,sr) = sums(:,74)     ! qr
    3556        hom(:,1,75,sr) = sums(:,75)     ! qc
    3557        hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
    3558                                        ! 77 is initial density profile
    3559        hom(:,1,78,sr) = ug             ! ug
    3560        hom(:,1,79,sr) = vg             ! vg
    3561        hom(:,1,80,sr) = w_subs         ! w_subs
    3562 
    3563        IF ( large_scale_forcing )  THEN
    3564           hom(:,1,81,sr) = sums_ls_l(:,0)          ! td_lsa_lpt
    3565           hom(:,1,82,sr) = sums_ls_l(:,1)          ! td_lsa_q
    3566           IF ( use_subsidence_tendencies )  THEN
    3567              hom(:,1,83,sr) = sums_ls_l(:,2)       ! td_sub_lpt
    3568              hom(:,1,84,sr) = sums_ls_l(:,3)       ! td_sub_q
    3569           ELSE
    3570              hom(:,1,83,sr) = sums(:,83)           ! td_sub_lpt
    3571              hom(:,1,84,sr) = sums(:,84)           ! td_sub_q
    3572           ENDIF
    3573           hom(:,1,85,sr) = sums(:,85)              ! td_nud_lpt
    3574           hom(:,1,86,sr) = sums(:,86)              ! td_nud_q
    3575           hom(:,1,87,sr) = sums(:,87)              ! td_nud_u
    3576           hom(:,1,88,sr) = sums(:,88)              ! td_nud_v
    3577        END IF
    3578 
    3579        hom(:,1,121,sr) = rho_air       ! rho_air in Kg/m^3
    3580        hom(:,1,122,sr) = rho_air_zw    ! rho_air_zw in Kg/m^3
    3581 
    3582        hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
    3583                                        ! u*, w'u', w'v', t* (in last profile)
    3584 
    3585        IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
    3586           hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
    3587                                sums(:,pr_palm+1:pr_palm+max_pr_user)
    3588        ENDIF
    3589 
    3590 !
    3591 !--    Determine the boundary layer height using two different schemes.
    3592 !--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
    3593 !--    first relative minimum (maximum) of the total heat flux.
    3594 !--    The corresponding height is assumed as the boundary layer height, if it
    3595 !--    is less than 1.5 times the height where the heat flux becomes negative
    3596 !--    (positive) for the first time.
    3597        z_i(1) = 0.0_wp
    3598        first = .TRUE.
    3599 
    3600        IF ( ocean )  THEN
    3601           DO  k = nzt, nzb+1, -1
    3602              IF (  first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
    3603                 first = .FALSE.
    3604                 height = zw(k)
    3605              ENDIF
    3606              IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                           &
    3607                   hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
    3608                 IF ( zw(k) < 1.5_wp * height )  THEN
    3609                    z_i(1) = zw(k)
    3610                 ELSE
    3611                    z_i(1) = height
    3612                 ENDIF
    3613                 EXIT
    3614              ENDIF
    3615           ENDDO
    3616        ELSE
    3617           DO  k = nzb, nzt-1
    3618              IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
    3619                 first = .FALSE.
    3620                 height = zw(k)
    3621              ENDIF
    3622              IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                           &
    3623                   hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
    3624                 IF ( zw(k) < 1.5_wp * height )  THEN
    3625                    z_i(1) = zw(k)
    3626                 ELSE
    3627                    z_i(1) = height
    3628                 ENDIF
    3629                 EXIT
    3630              ENDIF
    3631           ENDDO
    3632        ENDIF
    3633 
    3634 !
    3635 !--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
    3636 !--    by Uhlenbrock(2006). The boundary layer height is the height with the
    3637 !--    maximal local temperature gradient: starting from the second (the last
    3638 !--    but one) vertical gridpoint, the local gradient must be at least
    3639 !--    0.2K/100m and greater than the next four gradients.
    3640 !--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
    3641 !--             ocean case!
    3642        z_i(2) = 0.0_wp
    3643        DO  k = nzb+1, nzt+1
    3644           dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
    3645        ENDDO
    3646        dptdz_threshold = 0.2_wp / 100.0_wp
    3647 
    3648        IF ( ocean )  THEN
    3649           DO  k = nzt+1, nzb+5, -1
    3650              IF ( dptdz(k) > dptdz_threshold  .AND.                           &
    3651                   dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
    3652                   dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
    3653                 z_i(2) = zw(k-1)
    3654                 EXIT
    3655              ENDIF
    3656           ENDDO
    3657        ELSE
    3658           DO  k = nzb+1, nzt-3
    3659              IF ( dptdz(k) > dptdz_threshold  .AND.                           &
    3660                   dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
    3661                   dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
    3662                 z_i(2) = zw(k-1)
    3663                 EXIT
    3664              ENDIF
    3665           ENDDO
    3666        ENDIF
    3667 
    3668        hom(nzb+6,1,pr_palm,sr) = z_i(1)
    3669        hom(nzb+7,1,pr_palm,sr) = z_i(2)
    3670 
    3671 !
    3672 !--    Determine vertical index which is nearest to the mean surface level
    3673 !--    height of the respective statistic region
    3674        DO  k = nzb, nzt
    3675           IF ( zw(k) >= mean_surface_level_height(sr) )  THEN
    3676              k_surface_level = k
    3677              EXIT
    3678           ENDIF
    3679        ENDDO
    3680 
    3681 !
    3682 !--    Computation of both the characteristic vertical velocity and
    3683 !--    the characteristic convective boundary layer temperature.
    3684 !--    The inversion height entering into the equation is defined with respect
    3685 !--    to the mean surface level height of the respective statistic region.
    3686 !--    The horizontal average at surface level index + 1 is input for the
    3687 !--    average temperature.
    3688        IF ( hom(nzb,1,18,sr) > 1.0E-8_wp  .AND.  z_i(1) /= 0.0_wp )  THEN
    3689           hom(nzb+8,1,pr_palm,sr) = &
    3690              ( g / hom(k_surface_level+1,1,4,sr) *                             &
    3691              ( hom(k_surface_level,1,18,sr) / heatflux_output_conversion(nzb) )&
    3692              * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
    3693        ELSE
    3694           hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
    3695        ENDIF
    3696 
    3697 !
    3698 !--    Collect the time series quantities
    3699        ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
    3700        ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
    3701        ts_value(3,sr) = dt_3d
    3702        ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
    3703        ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
    3704        ts_value(6,sr) = u_max
    3705        ts_value(7,sr) = v_max
    3706        ts_value(8,sr) = w_max
    3707        ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
    3708        ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
    3709        ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
    3710        ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
    3711        ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
    3712        ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
    3713        ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
    3714        ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
    3715        ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
    3716        ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
    3717        ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
    3718        ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
    3719        ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
    3720 
    3721        IF ( .NOT. neutral )  THEN
    3722           ts_value(22,sr) = hom(nzb,1,114,sr)          ! L
    3723        ELSE
    3724           ts_value(22,sr) = 1.0E10_wp
    3725        ENDIF
    3726 
    3727        ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
    3728 
    3729 !
    3730 !--    Collect land surface model timeseries
    3731        IF ( land_surface )  THEN
    3732           ts_value(dots_soil  ,sr) = hom(nzb,1,93,sr)           ! ghf_eb
    3733           ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr)           ! shf_eb
    3734           ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr)           ! qsws_eb
    3735           ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr)           ! qsws_liq_eb
    3736           ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr)           ! qsws_soil_eb
    3737           ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr)           ! qsws_veg_eb
    3738           ts_value(dots_soil+6,sr) = hom(nzb,1,99,sr)           ! r_a
    3739           ts_value(dots_soil+7,sr) = hom(nzb,1,100,sr)          ! r_s
    3740        ENDIF
    3741 !
    3742 !--    Collect radiation model timeseries
    3743        IF ( radiation )  THEN
    3744           ts_value(dots_rad,sr)   = hom(nzb,1,101,sr)          ! rad_net
    3745           ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr)          ! rad_lw_in
    3746           ts_value(dots_rad+2,sr) = hom(nzb,1,103,sr)          ! rad_lw_out
    3747           ts_value(dots_rad+3,sr) = hom(nzb,1,104,sr)          ! rad_sw_in
    3748           ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr)          ! rad_sw_out
    3749 
    3750           IF ( radiation_scheme == 'rrtmg' )  THEN
    3751              ts_value(dots_rad+5,sr) = hom(nzb,1,106,sr)          ! rrtm_aldif
    3752              ts_value(dots_rad+6,sr) = hom(nzb,1,107,sr)          ! rrtm_aldir
    3753              ts_value(dots_rad+7,sr) = hom(nzb,1,108,sr)          ! rrtm_asdif
    3754              ts_value(dots_rad+8,sr) = hom(nzb,1,109,sr)          ! rrtm_asdir
    3755           ENDIF
    3756 
    3757        ENDIF
    3758 
    3759 !
    3760 !--    Calculate additional statistics provided by the user interface
    3761        CALL user_statistics( 'time_series', sr, 0 )
    3762 
    3763     ENDDO    ! loop of the subregions
    3764 
    3765     !$acc end data
    3766 
    3767 !
    3768 !-- If required, sum up horizontal averages for subsequent time averaging
    3769 !-- Do not sum, if flow statistics is called before the first initial time step.
    3770     IF ( do_sum  .AND.  simulated_time /= 0.0_wp )  THEN
    3771        IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
    3772        hom_sum = hom_sum + hom(:,1,:,:)
    3773        average_count_pr = average_count_pr + 1
    3774        do_sum = .FALSE.
    3775     ENDIF
    3776 
    3777 !
    3778 !-- Set flag for other UPs (e.g. output routines, but also buoyancy).
    3779 !-- This flag is reset after each time step in time_integration.
    3780     flow_statistics_called = .TRUE.
    3781 
    3782     CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
    3783 
    3784 
    3785  END SUBROUTINE flow_statistics
    3786 #endif
Note: See TracChangeset for help on using the changeset viewer.