Changeset 2118 for palm/trunk/SOURCE/flow_statistics.f90
 Timestamp:
 Jan 17, 2017 4:38:49 PM (6 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/flow_statistics.f90
r2101 r2118 20 20 ! Current revisions: 21 21 !  22 ! 22 ! OpenACC version of subroutine removed 23 23 ! 24 24 ! Former revisions: … … 216 216 !> are zero at the walls and inside buildings. 217 217 !! 218 #if ! defined( __openacc )219 218 SUBROUTINE flow_statistics 220 219 … … 309 308 CALL cpu_log( log_point(10), 'flow_statistics', 'start' ) 310 309 311 !$acc update host( km, kh, e, ol, pt, qs, qsws, shf, ts, u, usws, v, vsws, w )312 310 313 311 ! … … 1756 1754 1757 1755 END SUBROUTINE flow_statistics 1758 1759 1760 #else1761 1762 1763 !!1764 ! Description:1765 ! 1766 !> flow statistics  accelerator version1767 !!1768 SUBROUTINE flow_statistics1769 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, zw1777 1778 1779 USE cloud_parameters, &1780 ONLY: l_d_cp, prr, pt_d_t1781 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_sca1789 1790 USE cpulog, &1791 ONLY: cpu_log, log_point1792 1793 USE grid_variables, &1794 ONLY: ddx, ddy1795 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_invers1800 1801 USE kinds1802 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_soil1807 1808 USE netcdf_interface, &1809 ONLY: dots_rad, dots_soil1810 1811 USE pegrid1812 1813 USE radiation_model_mod, &1814 ONLY: radiation, radiation_scheme, rad_net, &1815 rad_lw_in, rad_lw_out, rad_sw_in, rad_sw_out1816 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_hr1821 #endif1822 1823 USE statistics1824 1825 IMPLICIT NONE1826 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 been1867 ! called once after the current time step1868 IF ( flow_statistics_called ) THEN1869 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 ENDIF1875 1876 !$acc data create( sums, sums_l )1877 !$acc update device( hom )1878 1879 !1880 ! Compute statistics for each (sub)region1881 DO sr = 0, statistic_regions1882 1883 !1884 ! Initialize (local) summation array1885 sums_l = 0.0_wp1886 1887 !1888 ! Store sums that have been computed in other subroutines in summation1889 ! array1890 sums_l(:,11,:) = sums_l_l(:,sr,:) ! mixing length from diffusivities1891 ! WARNING: next line still has to be adjusted for OpenMP1892 sums_l(:,21,0) = sums_wsts_bc_l(:,sr) * &1893 heatflux_output_conversion ! heat flux from advec_s_bc1894 sums_l(nzb+9,pr_palm,0) = sums_divold_l(sr) ! old divergence from pres1895 sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr) ! new divergence from pres1896 1897 !1898 ! When calcuating horizontallyaveraged total (resolved plus subgrid1899 ! scale) vertical fluxes and velocity variances by using commonly1900 ! applied Reynoldsbased methods ( e.g. <w'pt'> = (w<w>)*(pt<pt>) )1901 ! in combination with the 5th order advection scheme, pronounced1902 ! artificial kinks could be observed in the vertical profiles near the1903 ! 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 well1906 ! vertical velocity variances are calculated directly within the advection1907 ! routines, according to the numerical discretization, to evaluate the1908 ! statistical quantities as they will appear within the prognostic1909 ! equations.1910 ! Copy the turbulent quantities, evaluated in the advection routines to1911 ! the local array sums_l() for further computations.1912 IF ( ws_scheme_mom .AND. sr == 0 ) THEN1913 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 ) THEN1918 sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)1919 sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)1920 ENDIF1921 1922 DO i = 0, threads_per_task11923 !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*21930 sums_l(:,31,i) = sums_vs2_ws_l(:,i) ! v*21931 sums_l(:,32,i) = sums_ws2_ws_l(:,i) ! w*21932 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, nzt1936 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 ENDDO1941 ENDDO1942 1943 ENDIF1944 1945 IF ( ws_scheme_sca .AND. sr == 0 ) THEN1946 1947 DO i = 0, threads_per_task11948 sums_l(:,17,i) = sums_wspts_ws_l(:,i) &1949 * heatflux_output_conversion ! w*pt* from advec_s_ws1950 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 ENDDO1955 1956 ENDIF1957 !1958 ! Horizontally averaged profiles of horizontal velocities and temperature.1959 ! They must have been computed before, because they are already required1960 ! for other horizontal averages.1961 tn = 01962 1963 !$OMP PARALLEL PRIVATE( i, j, k, tn )1964 !$ tn = omp_get_thread_num()1965 1966 !$acc update device( sums_l )1967 1968 !$OMP DO1969 !$acc parallel loop gang present( pt, rflags_invers, rmask, sums_l, u, v ) create( s1, s2, s3 )1970 DO k = nzb, nzt+11971 s1 = 01972 s2 = 01973 s3 = 01974 !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )1975 DO i = nxl, nxr1976 DO j = nys, nyn1977 !1978 ! k+1 is used in rflags since rflags is set 0 at surface points1979 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 ENDDO1983 ENDDO1984 sums_l(k,1,tn) = s11985 sums_l(k,2,tn) = s21986 sums_l(k,4,tn) = s31987 ENDDO1988 !$acc end parallel loop1989 1990 !1991 ! Horizontally averaged profile of salinity1992 IF ( ocean ) THEN1993 !$OMP DO1994 !$acc parallel loop gang present( rflags_invers, rmask, sums_l, sa ) create( s1 )1995 DO k = nzb, nzt+11996 s1 = 01997 !$acc loop vector collapse( 2 ) reduction( +: s1 )1998 DO i = nxl, nxr1999 DO j = nys, nyn2000 s1 = s1 + sa(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)2001 ENDDO2002 ENDDO2003 sums_l(k,23,tn) = s12004 ENDDO2005 !$acc end parallel loop2006 ENDIF2007 2008 !2009 ! Horizontally averaged profiles of virtual potential temperature,2010 ! total water content, specific humidity and liquid water potential2011 ! temperature2012 IF ( humidity ) THEN2013 2014 !$OMP DO2015 !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l, vpt ) create( s1, s2 )2016 DO k = nzb, nzt+12017 s1 = 02018 s2 = 02019 !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )2020 DO i = nxl, nxr2021 DO j = nys, nyn2022 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 ENDDO2025 ENDDO2026 sums_l(k,41,tn) = s12027 sums_l(k,44,tn) = s22028 ENDDO2029 !$acc end parallel loop2030 2031 IF ( cloud_physics ) THEN2032 !$OMP DO2033 !$acc parallel loop gang present( pt, q, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )2034 DO k = nzb, nzt+12035 s1 = 02036 s2 = 02037 !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )2038 DO i = nxl, nxr2039 DO j = nys, nyn2040 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 ENDDO2045 ENDDO2046 sums_l(k,42,tn) = s12047 sums_l(k,43,tn) = s22048 ENDDO2049 !$acc end parallel loop2050 ENDIF2051 ENDIF2052 2053 !2054 ! Horizontally averaged profiles of passive scalar2055 IF ( passive_scalar ) THEN2056 !$OMP DO2057 !$acc parallel loop gang present( s, rflags_invers, rmask, sums_l ) create( s1 )2058 DO k = nzb, nzt+12059 s1 = 02060 !$acc loop vector collapse( 2 ) reduction( +: s1 )2061 DO i = nxl, nxr2062 DO j = nys, nyn2063 s1 = s1 + s(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)2064 ENDDO2065 ENDDO2066 sums_l(k,117,tn) = s12067 ENDDO2068 !$acc end parallel loop2069 ENDIF2070 !$OMP END PARALLEL2071 2072 !2073 ! Summation of thread sums2074 IF ( threads_per_task > 1 ) THEN2075 DO i = 1, threads_per_task12076 !$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 parallel2081 IF ( ocean ) THEN2082 !$acc parallel present( sums_l )2083 sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)2084 !$acc end parallel2085 ENDIF2086 IF ( humidity ) THEN2087 !$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 parallel2091 IF ( cloud_physics ) THEN2092 !$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 parallel2096 ENDIF2097 ENDIF2098 IF ( passive_scalar ) THEN2099 !$acc parallel present( sums_l )2100 sums_l(:,117,0) = sums_l(:,117,0) + sums_l(:,117,i)2101 !$acc end parallel2102 ENDIF2103 ENDDO2104 ENDIF2105 2106 #if defined( __parallel )2107 !2108 ! Compute total sum from local sums2109 !$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+2nzb, 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+2nzb, 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+2nzb, MPI_REAL, &2118 MPI_SUM, comm2d, ierr )2119 IF ( ocean ) THEN2120 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )2121 CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2nzb, &2122 MPI_REAL, MPI_SUM, comm2d, ierr )2123 ENDIF2124 IF ( humidity ) THEN2125 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )2126 CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2nzb, &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+2nzb, &2130 MPI_REAL, MPI_SUM, comm2d, ierr )2131 IF ( cloud_physics ) THEN2132 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )2133 CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2nzb, &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+2nzb, &2137 MPI_REAL, MPI_SUM, comm2d, ierr )2138 ENDIF2139 ENDIF2140 2141 IF ( passive_scalar ) THEN2142 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )2143 CALL MPI_ALLREDUCE( sums_l(nzb,117,0), sums(nzb,117), nzt+2nzb, &2144 MPI_REAL, MPI_SUM, comm2d, ierr )2145 ENDIF2146 !$acc update device( sums )2147 #else2148 !$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 parallel2153 IF ( ocean ) THEN2154 !$acc parallel present( sums, sums_l )2155 sums(:,23) = sums_l(:,23,0)2156 !$acc end parallel2157 ENDIF2158 IF ( humidity ) THEN2159 !$acc parallel present( sums, sums_l )2160 sums(:,44) = sums_l(:,44,0)2161 sums(:,41) = sums_l(:,41,0)2162 !$acc end parallel2163 IF ( cloud_physics ) THEN2164 !$acc parallel present( sums, sums_l )2165 sums(:,42) = sums_l(:,42,0)2166 sums(:,43) = sums_l(:,43,0)2167 !$acc end parallel2168 ENDIF2169 ENDIF2170 IF ( passive_scalar ) THEN2171 !$acc parallel present( sums, sums_l )2172 sums(:,117) = sums_l(:,117,0)2173 !$acc end parallel2174 ENDIF2175 #endif2176 2177 !2178 ! Final values are obtained by division by the total number of grid points2179 ! 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) ! u2185 hom(:,1,2,sr) = sums(:,2) ! v2186 hom(:,1,4,sr) = sums(:,4) ! pt2187 !$acc end parallel2188 2189 !2190 ! Salinity2191 IF ( ocean ) THEN2192 !$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) ! sa2195 !$acc end parallel2196 ENDIF2197 2198 !2199 ! Humidity and cloud parameters2200 IF ( humidity ) THEN2201 !$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) ! vpt2205 hom(:,1,41,sr) = sums(:,41) ! qv (q)2206 !$acc end parallel2207 IF ( cloud_physics ) THEN2208 !$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) ! qv2212 hom(:,1,43,sr) = sums(:,43) ! pt2213 !$acc end parallel2214 ENDIF2215 ENDIF2216 2217 !2218 ! Passive scalar2219 IF ( passive_scalar ) THEN2220 !$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) ! s2223 !$acc end parallel2224 ENDIF2225 2226 !2227 ! Horizontally averaged profiles of the remaining prognostic variables,2228 ! variances, the total and the perturbation energy (single values in last2229 ! column of sums_l) and some diagnostic quantities.2230 ! NOTE: for simplicity, nzb_s_inner is used below, although strictly2231 !  speaking the following kloop would have to be split up and2232 ! rearranged according to the staggered grid.2233 ! However, this implies no error since staggered velocity components2234 ! are zero at the walls and inside buildings.2235 tn = 02236 !$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 DO2242 !$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+12244 s1 = 02245 s2 = 02246 s3 = 02247 s4 = 02248 s5 = 02249 s6 = 02250 s7 = 02251 !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5, s6, s7 )2252 DO i = nxl, nxr2253 DO j = nys, nyn2254 !2255 ! Prognostic and diagnostic variables2256 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 moments2265 ! (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 ENDDO2268 ENDDO2269 sums_l(k,3,tn) = s12270 sums_l(k,8,tn) = s22271 sums_l(k,9,tn) = s32272 sums_l(k,10,tn) = s42273 sums_l(k,40,tn) = s52274 sums_l(k,33,tn) = s62275 sums_l(k,38,tn) = s72276 ENDDO2277 !$acc end parallel loop2278 2279 IF ( humidity ) THEN2280 !$OMP DO2281 !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l ) create( s1 )2282 DO k = nzb, nzt+12283 s1 = 02284 !$acc loop vector collapse( 2 ) reduction( +: s1 )2285 DO i = nxl, nxr2286 DO j = nys, nyn2287 s1 = s1 + ( q(k,j,i)hom(k,1,41,sr) )**2 * rmask(j,i,sr) * &2288 rflags_invers(j,i,k+1)2289 ENDDO2290 ENDDO2291 sums_l(k,70,tn) = s12292 ENDDO2293 !$acc end parallel loop2294 ENDIF2295 2296 !2297 ! Total and perturbation energy for the total domain (being2298 ! collected in the last column of sums_l).2299 s1 = 02300 !$OMP DO2301 !$acc parallel loop collapse(3) present( rflags_invers, rmask, u, v, w ) reduction(+:s1)2302 DO i = nxl, nxr2303 DO j = nys, nyn2304 DO k = nzb, nzt+12305 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 ENDDO2309 ENDDO2310 ENDDO2311 !$acc end parallel loop2312 !$acc parallel present( sums_l )2313 sums_l(nzb+4,pr_palm,tn) = s12314 !$acc end parallel2315 2316 !$OMP DO2317 !$acc parallel present( rmask, sums_l, us, usws, vsws, ts ) create( s1, s2, s3, s4 )2318 s1 = 02319 s2 = 02320 s3 = 02321 s4 = 02322 !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )2323 DO i = nxl, nxr2324 DO j = nys, nyn2325 !2326 ! 2Darrays (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 ENDDO2332 ENDDO2333 sums_l(nzb,pr_palm,tn) = s12334 sums_l(nzb+1,pr_palm,tn) = s22335 sums_l(nzb+2,pr_palm,tn) = s32336 sums_l(nzb+3,pr_palm,tn) = s42337 !$acc end parallel2338 2339 IF ( humidity ) THEN2340 !$acc parallel present( qs, rmask, sums_l ) create( s1 )2341 s1 = 02342 !$acc loop vector collapse( 2 ) reduction( +: s1 )2343 DO i = nxl, nxr2344 DO j = nys, nyn2345 s1 = s1 + qs(j,i) * rmask(j,i,sr)2346 ENDDO2347 ENDDO2348 sums_l(nzb+12,pr_palm,tn) = s12349 !$acc end parallel2350 ENDIF2351 2352 IF ( passive_scalar ) THEN2353 !$acc parallel present( ss, rmask, sums_l ) create( s1 )2354 s1 = 02355 !$acc loop vector collapse( 2 ) reduction( +: s1 )2356 DO i = nxl, nxr2357 DO j = nys, nyn2358 s1 = s1 + ss(j,i) * rmask(j,i,sr)2359 ENDDO2360 ENDDO2361 sums_l(nzb+13,pr_palm,tn) = s12362 !$acc end parallel2363 ENDIF2364 2365 !2366 ! Computation of statistics when wsscheme is not used. Else these2367 ! quantities are evaluated in the advection routines.2368 IF ( .NOT. ws_scheme_mom .OR. sr /= 0 .OR. simulated_time == 0.0_wp ) &2369 THEN2370 2371 !$OMP DO2372 !$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+12374 s1 = 02375 s2 = 02376 s3 = 02377 s4 = 02378 !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )2379 DO i = nxl, nxr2380 DO j = nys, nyn2381 ust2 = ( u(k,j,i)  hom(k,1,1,sr) )**22382 vst2 = ( v(k,j,i)  hom(k,1,2,sr) )**22383 w2 = w(k,j,i)**22384 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 energy2390 s4 = s4 + 0.5_wp * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) * &2391 rflags_invers(j,i,k+1)2392 ENDDO2393 ENDDO2394 sums_l(k,30,tn) = s12395 sums_l(k,31,tn) = s22396 sums_l(k,32,tn) = s32397 sums_l(k,34,tn) = s42398 ENDDO2399 !$acc end parallel loop2400 !2401 ! Total perturbation TKE2402 !$OMP DO2403 !$acc parallel present( sums_l ) create( s1 )2404 s1 = 02405 !$acc loop reduction( +: s1 )2406 DO k = nzb, nzt+12407 s1 = s1 + sums_l(k,34,tn)2408 ENDDO2409 sums_l(nzb+5,pr_palm,tn) = s12410 !$acc end parallel2411 2412 ENDIF2413 2414 !2415 ! Horizontally averaged profiles of the vertical fluxes2416 2417 !2418 ! Subgridscale fluxes.2419 ! WARNING: If a Prandtllayer is used (k=nzb for flat terrain), the fluxes2420 !  should be calculated there in a different way. This is done2421 ! in the next loop further below, where results from this loop are2422 ! overwritten. However, THIS WORKS IN CASE OF FLAT TERRAIN ONLY!2423 ! The nonflat case still has to be handled.2424 ! NOTE: for simplicity, nzb_s_inner is used below, although2425 !  strictly speaking the following kloop would have to be2426 ! split up according to the staggered grid.2427 ! However, this implies no error since staggered velocity2428 ! components are zero at the walls and inside buildings.2429 !$OMP DO2430 !$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_diff2432 s1 = 02433 s2 = 02434 s3 = 02435 !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )2436 DO i = nxl, nxr2437 DO j = nys, nyn2438 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,i1)+km(k+1,j,i1) &2443 ) * ( &2444 ( u(k+1,j,i)  u(k,j,i) ) * ddzu(k+1) &2445 + ( w(k,j,i)  w(k,j,i1) ) * 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,j1,i)+km(k+1,j1,i) &2454 ) * ( &2455 ( v(k+1,j,i)  v(k,j,i) ) * ddzu(k+1) &2456 + ( w(k,j,i)  w(k,j1,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 ENDDO2470 ENDDO2471 sums_l(k,12,tn) = s12472 sums_l(k,14,tn) = s22473 sums_l(k,16,tn) = s32474 ENDDO2475 !$acc end parallel loop2476 2477 !2478 ! Salinity flux w"sa"2479 IF ( ocean ) THEN2480 !$acc parallel loop gang present( ddzu, kh, sa, rflags_invers, rmask, sums_l ) create( s1 )2481 DO k = nzb, nzt_diff2482 s1 = 02483 !$acc loop vector collapse( 2 ) reduction( +: s1 )2484 DO i = nxl, nxr2485 DO j = nys, nyn2486 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 ENDDO2491 ENDDO2492 sums_l(k,65,tn) = s12493 ENDDO2494 !$acc end parallel loop2495 ENDIF2496 2497 !2498 ! Buoyancy flux, water flux (humidity flux) w"q"2499 IF ( humidity ) THEN2500 2501 !$acc parallel loop gang present( ddzu, kh, q, vpt, rflags_invers, rmask, sums_l ) create( s1, s2 )2502 DO k = nzb, nzt_diff2503 s1 = 02504 s2 = 02505 !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )2506 DO i = nxl, nxr2507 DO j = nys, nyn2508 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 ENDDO2521 ENDDO2522 sums_l(k,45,tn) = s12523 sums_l(k,48,tn) = s22524 ENDDO2525 !$acc end parallel loop2526 2527 IF ( cloud_physics ) THEN2528 2529 !$acc parallel loop gang present( ddzu, kh, q, ql, rflags_invers, rmask, sums_l ) create( s1 )2530 DO k = nzb, nzt_diff2531 s1 = 02532 !$acc loop vector collapse( 2 ) reduction( +: s1 )2533 DO i = nxl, nxr2534 DO j = nys, nyn2535 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 ENDDO2543 ENDDO2544 sums_l(k,51,tn) = s12545 ENDDO2546 !$acc end parallel loop2547 2548 ENDIF2549 2550 ENDIF2551 !2552 ! Passive scalar flux2553 IF ( passive_scalar ) THEN2554 2555 !$acc parallel loop gang present( ddzu, kh, s, rflags_invers, rmask, sums_l ) create( s1 )2556 DO k = nzb, nzt_diff2557 s1 = 02558 !$acc loop vector collapse( 2 ) reduction( +: s1 )2559 DO i = nxl, nxr2560 DO j = nys, nyn2561 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 ENDDO2566 ENDDO2567 sums_l(k,119,tn) = s12568 ENDDO2569 !$acc end parallel loop2570 2571 ENDIF2572 2573 IF ( use_surface_fluxes ) THEN2574 2575 !$OMP DO2576 !$acc parallel present( rmask, shf, sums_l, usws, vsws ) create( s1, s2, s3, s4, s5 )2577 s1 = 02578 s2 = 02579 s3 = 02580 s4 = 02581 s5 = 02582 !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )2583 DO i = nxl, nxr2584 DO j = nys, nyn2585 !2586 ! Subgridscale fluxes in the Prandtl layer2587 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 ENDDO2596 ENDDO2597 sums_l(nzb,12,tn) = s12598 sums_l(nzb,14,tn) = s22599 sums_l(nzb,16,tn) = s32600 sums_l(nzb,58,tn) = s42601 sums_l(nzb,61,tn) = s52602 !$acc end parallel2603 2604 IF ( ocean ) THEN2605 2606 !$OMP DO2607 !$acc parallel present( rmask, saswsb, sums_l ) create( s1 )2608 s1 = 02609 !$acc loop vector collapse( 2 ) reduction( +: s1 )2610 DO i = nxl, nxr2611 DO j = nys, nyn2612 s1 = s1 + saswsb(j,i) * rmask(j,i,sr) ! w"sa"2613 ENDDO2614 ENDDO2615 sums_l(nzb,65,tn) = s12616 !$acc end parallel2617 2618 ENDIF2619 2620 IF ( humidity ) THEN2621 2622 !$OMP DO2623 !$acc parallel present( pt, q, qsws, rmask, shf, sums_l ) create( s1, s2 )2624 s1 = 02625 s2 = 02626 !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )2627 DO i = nxl, nxr2628 DO j = nys, nyn2629 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 ENDDO2635 ENDDO2636 sums_l(nzb,48,tn) = s12637 sums_l(nzb,45,tn) = s22638 !$acc end parallel2639 2640 IF ( cloud_droplets ) THEN2641 2642 !$OMP DO2643 !$acc parallel present( pt, q, ql, qsws, rmask, shf, sums_l ) create( s1 )2644 s1 = 02645 !$acc loop vector collapse( 2 ) reduction( +: s1 )2646 DO i = nxl, nxr2647 DO j = nys, nyn2648 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 ENDDO2653 ENDDO2654 sums_l(nzb,45,tn) = s12655 !$acc end parallel2656 2657 ENDIF2658 2659 IF ( cloud_physics ) THEN2660 2661 !$OMP DO2662 !$acc parallel present( qsws, rmask, sums_l ) create( s1 )2663 s1 = 02664 !$acc loop vector collapse( 2 ) reduction( +: s1 )2665 DO i = nxl, nxr2666 DO j = nys, nyn2667 !2668 ! Formula does not work if ql(nzb) /= 0.02669 s1 = s1 + qsws(j,i) * waterflux_output_conversion(nzb) &2670 * rmask(j,i,sr) ! w"q" (w"qv")2671 ENDDO2672 ENDDO2673 sums_l(nzb,51,tn) = s12674 !$acc end parallel2675 2676 ENDIF2677 2678 ENDIF2679 2680 IF ( passive_scalar ) THEN2681 2682 !$OMP DO2683 !$acc parallel present( ssws, rmask, sums_l ) create( s1 )2684 s1 = 02685 !$acc loop vector collapse( 2 ) reduction( +: s1 )2686 DO i = nxl, nxr2687 DO j = nys, nyn2688 s1 = s1 + ssws(j,i) * rmask(j,i,sr) ! w"s"2689 ENDDO2690 ENDDO2691 sums_l(nzb,119,tn) = s12692 !$acc end parallel2693 2694 ENDIF2695 2696 ENDIF2697 2698 !2699 ! Subgridscale fluxes at the top surface2700 IF ( use_top_fluxes ) THEN2701 2702 !$OMP DO2703 !$acc parallel present( rmask, sums_l, tswst, uswst, vswst ) create( s1, s2, s3, s4, s5 )2704 s1 = 02705 s2 = 02706 s3 = 02707 s4 = 02708 s5 = 02709 !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )2710 DO i = nxl, nxr2711 DO j = nys, nyn2712 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 ENDDO2721 ENDDO2722 sums_l(nzt:nzt+1,12,tn) = s12723 sums_l(nzt:nzt+1,14,tn) = s22724 sums_l(nzt:nzt+1,16,tn) = s32725 sums_l(nzt:nzt+1,58,tn) = s42726 sums_l(nzt:nzt+1,61,tn) = s52727 !$acc end parallel2728 2729 IF ( ocean ) THEN2730 2731 !$OMP DO2732 !$acc parallel present( rmask, saswst, sums_l ) create( s1 )2733 s1 = 02734 !$acc loop vector collapse( 2 ) reduction( +: s1 )2735 DO i = nxl, nxr2736 DO j = nys, nyn2737 s1 = s1 + saswst(j,i) * rmask(j,i,sr) ! w"sa"2738 ENDDO2739 ENDDO2740 sums_l(nzt,65,tn) = s12741 !$acc end parallel2742 2743 ENDIF2744 2745 IF ( humidity ) THEN2746 2747 !$OMP DO2748 !$acc parallel present( pt, q, qswst, rmask, tswst, sums_l ) create( s1, s2 )2749 s1 = 02750 s2 = 02751 !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )2752 DO i = nxl, nxr2753 DO j = nys, nyn2754 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 ENDDO2760 ENDDO2761 sums_l(nzt,48,tn) = s12762 sums_l(nzt,45,tn) = s22763 !$acc end parallel2764 2765 IF ( cloud_droplets ) THEN2766 2767 !$OMP DO2768 !$acc parallel present( pt, q, ql, qswst, rmask, tswst, sums_l ) create( s1 )2769 s1 = 02770 !$acc loop vector collapse( 2 ) reduction( +: s1 )2771 DO i = nxl, nxr2772 DO j = nys, nyn2773 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 ENDDO2779 ENDDO2780 sums_l(nzt,45,tn) = s12781 !$acc end parallel2782 2783 ENDIF2784 2785 IF ( cloud_physics ) THEN2786 2787 !$OMP DO2788 !$acc parallel present( qswst, rmask, sums_l ) create( s1 )2789 s1 = 02790 !$acc loop vector collapse( 2 ) reduction( +: s1 )2791 DO i = nxl, nxr2792 DO j = nys, nyn2793 !2794 ! Formula does not work if ql(nzb) /= 0.02795 s1 = s1 + qswst(j,i) * waterflux_output_conversion(nzt) &2796 * rmask(j,i,sr) ! w"q" (w"qv")2797 ENDDO2798 ENDDO2799 sums_l(nzt,51,tn) = s12800 !$acc end parallel2801 2802 ENDIF2803 2804 ENDIF2805 2806 IF ( passive_scalar ) THEN2807 2808 !$OMP DO2809 !$acc parallel present( sswst, rmask, sums_l ) create( s1 )2810 s1 = 02811 !$acc loop vector collapse( 2 ) reduction( +: s1 )2812 DO i = nxl, nxr2813 DO j = nys, nyn2814 s1 = s1 + sswst(j,i) * rmask(j,i,sr) ! w"s"2815 ENDDO2816 ENDDO2817 sums_l(nzt,119,tn) = s12818 !$acc end parallel2819 2820 ENDIF2821 2822 ENDIF2823 2824 !2825 ! Resolved fluxes (can be computed for all horizontal points)2826 ! NOTE: for simplicity, nzb_s_inner is used below, although strictly2827 !  speaking the following kloop would have to be split up and2828 ! 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_diff2831 s1 = 02832 s2 = 02833 s3 = 02834 !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )2835 DO i = nxl, nxr2836 DO j = nys, nyn2837 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 moments2845 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 ENDDO2853 ENDDO2854 sums_l(k,35,tn) = s12855 sums_l(k,36,tn) = s22856 sums_l(k,37,tn) = s32857 ENDDO2858 !$acc end parallel loop2859 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 ) THEN2864 2865 IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN2866 2867 !$acc parallel loop gang present( hom, rflags_invers, rmask, sa, sums_l, w ) create( s1 )2868 DO k = nzb, nzt_diff2869 s1 = 02870 !$acc loop vector collapse( 2 ) reduction( +: s1 )2871 DO i = nxl, nxr2872 DO j = nys, nyn2873 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 ENDDO2878 ENDDO2879 sums_l(k,66,tn) = s12880 ENDDO2881 !$acc end parallel loop2882 2883 ENDIF2884 2885 !$acc parallel loop gang present( rflags_invers, rho_ocean, prho, rmask, sums_l ) create( s1, s2 )2886 DO k = nzb, nzt_diff2887 s1 = 02888 s2 = 02889 !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )2890 DO i = nxl, nxr2891 DO j = nys, nyn2892 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 ENDDO2895 ENDDO2896 sums_l(k,64,tn) = s12897 sums_l(k,71,tn) = s22898 ENDDO2899 !$acc end parallel loop2900 2901 ENDIF2902 2903 !2904 ! Buoyancy flux, water flux, humidity flux, liquid water2905 ! content, rain drop concentration and rain water content2906 IF ( humidity ) THEN2907 2908 IF ( cloud_physics .OR. cloud_droplets ) THEN2909 2910 !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )2911 DO k = nzb, nzt_diff2912 s1 = 02913 !$acc loop vector collapse( 2 ) reduction( +: s1 )2914 DO i = nxl, nxr2915 DO j = nys, nyn2916 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 ENDDO2921 ENDDO2922 sums_l(k,46,tn) = s12923 ENDDO2924 !$acc end parallel loop2925 2926 IF ( .NOT. cloud_droplets ) THEN2927 2928 !$acc parallel loop gang present( hom, q, ql, rflags_invers, rmask, sums_l, w ) create( s1 )2929 DO k = nzb, nzt_diff2930 s1 = 02931 !$acc loop vector collapse( 2 ) reduction( +: s1 )2932 DO i = nxl, nxr2933 DO j = nys, nyn2934 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 ENDDO2939 ENDDO2940 sums_l(k,52,tn) = s12941 ENDDO2942 !$acc end parallel loop2943 2944 IF ( microphysics_seifert ) THEN2945 2946 !$acc parallel loop gang present( qc, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )2947 DO k = nzb, nzt_diff2948 s1 = 02949 !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )2950 DO i = nxl, nxr2951 DO j = nys, nyn2952 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 ENDDO2955 ENDDO2956 sums_l(k,54,tn) = s12957 sums_l(k,75,tn) = s22958 ENDDO2959 !$acc end parallel loop2960 2961 !$acc parallel loop gang present( nr, qr, prr, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )2962 DO k = nzb, nzt_diff2963 s1 = 02964 s2 = 02965 s3 = 02966 !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )2967 DO i = nxl, nxr2968 DO j = nys, nyn2969 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 ENDDO2973 ENDDO2974 sums_l(k,73,tn) = s12975 sums_l(k,74,tn) = s22976 sums_l(k,76,tn) = s32977 ENDDO2978 !$acc end parallel loop2979 2980 ELSE2981 2982 !$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )2983 DO k = nzb, nzt_diff2984 s1 = 02985 !$acc loop vector collapse( 2 ) reduction( +: s1 )2986 DO i = nxl, nxr2987 DO j = nys, nyn2988 s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)2989 ENDDO2990 ENDDO2991 sums_l(k,54,tn) = s12992 ENDDO2993 !$acc end parallel loop2994 2995 ENDIF2996 2997 ELSE2998 2999 !$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )3000 DO k = nzb, nzt_diff3001 s1 = 03002 !$acc loop vector collapse( 2 ) reduction( +: s1 )3003 DO i = nxl, nxr3004 DO j = nys, nyn3005 s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)3006 ENDDO3007 ENDDO3008 sums_l(k,54,tn) = s13009 ENDDO3010 !$acc end parallel loop3011 3012 ENDIF3013 3014 ELSE3015 3016 IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN3017 3018 !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )3019 DO k = nzb, nzt_diff3020 s1 = 03021 !$acc loop vector collapse( 2 ) reduction( +: s1 )3022 DO i = nxl, nxr3023 DO j = nys, nyn3024 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 ENDDO3029 ENDDO3030 sums_l(k,46,tn) = s13031 ENDDO3032 !$acc end parallel loop3033 3034 ELSEIF ( ws_scheme_sca .AND. sr == 0 ) THEN3035 3036 !$acc parallel loop present( hom, sums_l )3037 DO k = nzb, nzt_diff3038 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 ENDDO3043 !$acc end parallel loop3044 3045 ENDIF3046 3047 ENDIF3048 3049 ENDIF3050 !3051 ! Passive scalar flux3052 IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca .OR. sr /= 0 ) ) THEN3053 3054 !$acc parallel loop gang present( hom, s, rflags_invers, rmask, sums_l, w ) create( s1 )3055 DO k = nzb, nzt_diff3056 s1 = 03057 !$acc loop vector collapse( 2 ) reduction( +: s1 )3058 DO i = nxl, nxr3059 DO j = nys, nyn3060 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 ENDDO3065 ENDDO3066 sums_l(k,49,tn) = s13067 ENDDO3068 !$acc end parallel loop3069 3070 ENDIF3071 3072 !3073 ! For speed optimization fluxes which have been computed in part directly3074 ! inside the WS advection routines are treated seperatly3075 ! Momentum fluxes first:3076 IF ( .NOT. ws_scheme_mom .OR. sr /= 0 ) THEN3077 3078 !$OMP DO3079 !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2 )3080 DO k = nzb, nzt_diff3081 s1 = 03082 s2 = 03083 !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )3084 DO i = nxl, nxr3085 DO j = nys, nyn3086 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,i1) + 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,j1,i) + w(k,j,i) ) &3099 * vst * rmask(j,i,sr) &3100 * momentumflux_output_conversion(k) &3101 * rflags_invers(j,i,k+1)3102 ENDDO3103 ENDDO3104 sums_l(k,13,tn) = s13105 sums_l(k,15,tn) = s23106 ENDDO3107 !$acc end parallel loop3108 3109 ENDIF3110 3111 IF ( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN3112 3113 !$OMP DO3114 !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, w ) create( s1 )3115 DO k = nzb, nzt_diff3116 s1 = 03117 !$acc loop vector collapse( 2 ) reduction( +: s1 )3118 DO i = nxl, nxr3119 DO j = nys, nyn3120 !3121 ! Vertical heat flux3122 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 ENDDO3128 ENDDO3129 sums_l(k,17,tn) = s13130 ENDDO3131 !$acc end parallel loop3132 3133 IF ( humidity ) THEN3134 3135 !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )3136 DO k = nzb, nzt_diff3137 s1 = 03138 !$acc loop vector collapse( 2 ) reduction( +: s1 )3139 DO i = nxl, nxr3140 DO j = nys, nyn3141 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 ENDDO3147 ENDDO3148 sums_l(k,49,tn) = s13149 ENDDO3150 !$acc end parallel loop3151 3152 ENDIF3153 3154 IF ( passive_scalar ) THEN3155 3156 !$acc parallel loop gang present( hom, s, rflags_invers, rmask, sums_l, w ) create( s1 )3157 DO k = nzb, nzt_diff3158 s1 = 03159 !$acc loop vector collapse( 2 ) reduction( +: s1 )3160 DO i = nxl, nxr3161 DO j = nys, nyn3162 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 ENDDO3167 ENDDO3168 sums_l(k,116,tn) = s13169 ENDDO3170 !$acc end parallel loop3171 3172 ENDIF3173 3174 ENDIF3175 3176 3177 !3178 ! Density at top follows Neumann condition3179 IF ( ocean ) THEN3180 !$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 parallel3184 ENDIF3185 3186 !3187 ! Divergence of vertical flux of resolved scale energy and pressure3188 ! 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 ) THEN3192 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 array3195 3196 !$OMP DO3197 DO i = nxl, nxr3198 DO j = nys, nyn3199 DO k = nzb_s_inner(j,i)+1, nzt3200 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 ENDDO3212 ENDDO3213 ENDDO3214 sums_ll(0,1) = 0.0_wp ! because w is zero at the bottom3215 sums_ll(nzt+1,1) = 0.0_wp3216 sums_ll(0,2) = 0.0_wp3217 sums_ll(nzt+1,2) = 0.0_wp3218 3219 DO k = nzb+1, nzt3220 sums_l(k,55,tn) = ( sums_ll(k,1)  sums_ll(k1,1) ) * ddzw(k)3221 sums_l(k,56,tn) = ( sums_ll(k,2)  sums_ll(k1,2) ) * ddzw(k)3222 sums_l(k,68,tn) = sums_ll(k,2)3223 ENDDO3224 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 nzb3227 3228 ENDIF3229 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 ) THEN3233 3234 STOP '+++ openACC porting for vertical flux div of SGS TKE in flow_statistics is still missing'3235 !$OMP DO3236 DO i = nxl, nxr3237 DO j = nys, nyn3238 DO k = nzb_s_inner(j,i)+1, nzt3239 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(k1,j,i)+km(k,j,i)) * (e(k,j,i)e(k1,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 ENDDO3250 ENDDO3251 ENDDO3252 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 ENDIF3256 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 ) THEN3261 3262 STOP '+++ openACC porting for horizontal flux calculation in flow_statistics is still missing'3263 !$OMP DO3264 DO i = nxl, nxr3265 DO j = nys, nyn3266 DO k = nzb_s_inner(j,i)+1, nzt3267 !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,i1) ) &3271 * ( pt(k,j,i1)  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,j1,i) ) &3277 * ( pt(k,j1,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,i1)  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,j1,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,j1,i)  hom(k,1,4,sr) + &3293 pt(k,j,i)  hom(k,1,4,sr) ) &3294 * heatflux_output_conversion(k)3295 ENDDO3296 ENDDO3297 ENDDO3298 !3299 ! Fluxes at the surface must be zero (e.g. due to the Prandtllayer)3300 sums_l(nzb,58,tn) = 0.0_wp3301 sums_l(nzb,59,tn) = 0.0_wp3302 sums_l(nzb,60,tn) = 0.0_wp3303 sums_l(nzb,61,tn) = 0.0_wp3304 sums_l(nzb,62,tn) = 0.0_wp3305 sums_l(nzb,63,tn) = 0.0_wp3306 3307 ENDIF3308 3309 !3310 ! Collect current large scale advection and subsidence tendencies for3311 ! data output3312 IF ( large_scale_forcing .AND. ( simulated_time > 0.0_wp ) ) THEN3313 !3314 ! Interpolation in time of LSF_DATA3315 nt = 13316 DO WHILE ( simulated_time  dt_3d > time_vert(nt) )3317 nt = nt + 13318 ENDDO3319 IF ( simulated_time  dt_3d /= time_vert(nt) ) THEN3320 nt = nt  13321 ENDIF3322 3323 fac = ( simulated_time  dt_3d  time_vert(nt) ) &3324 / ( time_vert(nt+1)time_vert(nt) )3325 3326 3327 DO k = nzb, nzt3328 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 ENDDO3333 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 ) THEN3338 3339 DO k = nzb, nzt3340 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 ENDDO3345 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 ENDIF3350 3351 ENDIF3352 3353 3354 IF ( land_surface ) THEN3355 !$OMP DO3356 DO i = nxl, nxr3357 DO j = nys, nyn3358 DO k = nzb_soil, nzt_soil3359 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 ENDDO3364 ENDDO3365 ENDDO3366 ENDIF3367 3368 3369 IF ( radiation .AND. radiation_scheme == 'rrtmg' ) THEN3370 !$OMP DO3371 DO i = nxl, nxr3372 DO j = nys, nyn3373 DO k = nzb_s_inner(j,i)+1, nzt+13374 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 #endif3392 ENDDO3393 ENDDO3394 ENDDO3395 ENDIF3396 3397 !3398 ! Calculate the userdefined profiles3399 CALL user_statistics( 'profiles', sr, tn )3400 !$OMP END PARALLEL3401 3402 !3403 ! Summation of thread sums3404 IF ( threads_per_task > 1 ) THEN3405 STOP '+++ openACC porting for threads_per_task > 1 in flow_statistics is still missing'3406 DO i = 1, threads_per_task13407 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 ) THEN3412 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 ENDIF3416 ENDDO3417 ENDIF3418 3419 !$acc update host( hom, sums, sums_l )3420 3421 #if defined( __parallel )3422 3423 !3424 ! Compute total sum from local sums3425 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 ) THEN3429 CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls, &3430 MPI_REAL, MPI_SUM, comm2d, ierr )3431 ENDIF3432 #else3433 sums = sums_l(:,:,0)3434 IF ( large_scale_forcing ) THEN3435 sums(:,81:88) = sums_ls_l3436 ENDIF3437 #endif3438 3439 !3440 ! Final values are obtained by division by the total number of grid points3441 ! used for summation. After that store profiles.3442 ! Check, if statistical regions do contain at least one grid point at the3443 ! respective klevel, otherwise division by zero will lead to undefined3444 ! values, which may cause e.g. problems with NetCDF output3445 ! Profiles:3446 DO k = nzb, nzt+13447 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 ) THEN3456 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_palm2) = sums(k,115:pr_palm2) / ngp_2dh_s_inner(k,sr)3464 ENDIF3465 ENDDO3466 3467 ! u* and so on3468 ! As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose3469 ! size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer3470 ! 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) / & ! qs3474 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 divergence3479 sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &3480 ngp_3d_inner(sr)3481 3482 ! Userdefined profiles3483 IF ( max_pr_user > 0 ) THEN3484 DO k = nzb, nzt+13485 IF ( ngp_2dh_s_inner(k,sr) /= 0 ) THEN3486 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 ENDIF3490 ENDDO3491 ENDIF3492 3493 !3494 ! Collect horizontal average in hom.3495 ! Compute deduced averages (e.g. total heat flux)3496 hom(:,1,3,sr) = sums(:,3) ! w3497 hom(:,1,8,sr) = sums(:,8) ! e profiles 57 are initial profiles3498 hom(:,1,9,sr) = sums(:,9) ! km3499 hom(:,1,10,sr) = sums(:,10) ! kh3500 hom(:,1,11,sr) = sums(:,11) ! l3501 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) ! wpt3508 hom(:,1,19,sr) = sums(:,12) + sums(:,13) ! wu3509 hom(:,1,20,sr) = sums(:,14) + sums(:,15) ! wv3510 hom(:,1,21,sr) = sums(:,21) ! w*pt*BC3511 hom(:,1,22,sr) = sums(:,16) + sums(:,21) ! wptBC3512 ! profile 24 is initial profile (sa)3513 ! profiles 2529 left empty for initial3514 ! profiles3515 hom(:,1,30,sr) = sums(:,30) ! u*23516 hom(:,1,31,sr) = sums(:,31) ! v*23517 hom(:,1,32,sr) = sums(:,32) ! w*23518 hom(:,1,33,sr) = sums(:,33) ! pt*23519 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*23522 hom(:,1,37,sr) = sums(:,37) ! w*e*3523 hom(:,1,38,sr) = sums(:,38) ! w*33524 hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E20_wp )**1.5_wp ! Sw3525 hom(:,1,40,sr) = sums(:,40) ! p3526 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) ! wvpt3529 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) ! ql3536 hom(:,1,55,sr) = sums(:,55) ! w*u*u*/dz3537 hom(:,1,56,sr) = sums(:,56) ! w*p*/dz3538 hom(:,1,57,sr) = sums(:,57) ! ( w"e + w"p"/rho_ocean )/dz3539 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_t3542 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_t3545 hom(:,1,64,sr) = sums(:,64) ! rho_ocean3546 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) ! wsa3549 hom(:,1,68,sr) = sums(:,68) ! w*p*3550 hom(:,1,69,sr) = sums(:,69) ! w"e + w"p"/rho_ocean3551 hom(:,1,70,sr) = sums(:,70) ! q*23552 hom(:,1,71,sr) = sums(:,71) ! prho3553 hom(:,1,72,sr) = hyp * 1E4_wp ! hyp in dbar3554 hom(:,1,73,sr) = sums(:,73) ! nr3555 hom(:,1,74,sr) = sums(:,74) ! qr3556 hom(:,1,75,sr) = sums(:,75) ! qc3557 hom(:,1,76,sr) = sums(:,76) ! prr (precipitation rate)3558 ! 77 is initial density profile3559 hom(:,1,78,sr) = ug ! ug3560 hom(:,1,79,sr) = vg ! vg3561 hom(:,1,80,sr) = w_subs ! w_subs3562 3563 IF ( large_scale_forcing ) THEN3564 hom(:,1,81,sr) = sums_ls_l(:,0) ! td_lsa_lpt3565 hom(:,1,82,sr) = sums_ls_l(:,1) ! td_lsa_q3566 IF ( use_subsidence_tendencies ) THEN3567 hom(:,1,83,sr) = sums_ls_l(:,2) ! td_sub_lpt3568 hom(:,1,84,sr) = sums_ls_l(:,3) ! td_sub_q3569 ELSE3570 hom(:,1,83,sr) = sums(:,83) ! td_sub_lpt3571 hom(:,1,84,sr) = sums(:,84) ! td_sub_q3572 ENDIF3573 hom(:,1,85,sr) = sums(:,85) ! td_nud_lpt3574 hom(:,1,86,sr) = sums(:,86) ! td_nud_q3575 hom(:,1,87,sr) = sums(:,87) ! td_nud_u3576 hom(:,1,88,sr) = sums(:,88) ! td_nud_v3577 END IF3578 3579 hom(:,1,121,sr) = rho_air ! rho_air in Kg/m^33580 hom(:,1,122,sr) = rho_air_zw ! rho_air_zw in Kg/m^33581 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 ! userdefined profiles3586 hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &3587 sums(:,pr_palm+1:pr_palm+max_pr_user)3588 ENDIF3589 3590 !3591 ! Determine the boundary layer height using two different schemes.3592 ! First scheme: Starting from the Earth's (Ocean's) surface, look for the3593 ! first relative minimum (maximum) of the total heat flux.3594 ! The corresponding height is assumed as the boundary layer height, if it3595 ! is less than 1.5 times the height where the heat flux becomes negative3596 ! (positive) for the first time.3597 z_i(1) = 0.0_wp3598 first = .TRUE.3599 3600 IF ( ocean ) THEN3601 DO k = nzt, nzb+1, 13602 IF ( first .AND. hom(k,1,18,sr) < 1.0E8_wp ) THEN3603 first = .FALSE.3604 height = zw(k)3605 ENDIF3606 IF ( hom(k,1,18,sr) < 1.0E8_wp .AND. &3607 hom(k1,1,18,sr) > hom(k,1,18,sr) ) THEN3608 IF ( zw(k) < 1.5_wp * height ) THEN3609 z_i(1) = zw(k)3610 ELSE3611 z_i(1) = height3612 ENDIF3613 EXIT3614 ENDIF3615 ENDDO3616 ELSE3617 DO k = nzb, nzt13618 IF ( first .AND. hom(k,1,18,sr) < 1.0E8_wp ) THEN3619 first = .FALSE.3620 height = zw(k)3621 ENDIF3622 IF ( hom(k,1,18,sr) < 1.0E8_wp .AND. &3623 hom(k+1,1,18,sr) > hom(k,1,18,sr) ) THEN3624 IF ( zw(k) < 1.5_wp * height ) THEN3625 z_i(1) = zw(k)3626 ELSE3627 z_i(1) = height3628 ENDIF3629 EXIT3630 ENDIF3631 ENDDO3632 ENDIF3633 3634 !3635 ! Second scheme: Gradient scheme from Sullivan et al. (1998), modified3636 ! by Uhlenbrock(2006). The boundary layer height is the height with the3637 ! maximal local temperature gradient: starting from the second (the last3638 ! but one) vertical gridpoint, the local gradient must be at least3639 ! 0.2K/100m and greater than the next four gradients.3640 ! WARNING: The threshold value of 0.2K/100m must be adjusted for the3641 ! ocean case!3642 z_i(2) = 0.0_wp3643 DO k = nzb+1, nzt+13644 dptdz(k) = ( hom(k,1,4,sr)  hom(k1,1,4,sr) ) * ddzu(k)3645 ENDDO3646 dptdz_threshold = 0.2_wp / 100.0_wp3647 3648 IF ( ocean ) THEN3649 DO k = nzt+1, nzb+5, 13650 IF ( dptdz(k) > dptdz_threshold .AND. &3651 dptdz(k) > dptdz(k1) .AND. dptdz(k) > dptdz(k2) .AND. &3652 dptdz(k) > dptdz(k3) .AND. dptdz(k) > dptdz(k4) ) THEN3653 z_i(2) = zw(k1)3654 EXIT3655 ENDIF3656 ENDDO3657 ELSE3658 DO k = nzb+1, nzt33659 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) ) THEN3662 z_i(2) = zw(k1)3663 EXIT3664 ENDIF3665 ENDDO3666 ENDIF3667 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 level3673 ! height of the respective statistic region3674 DO k = nzb, nzt3675 IF ( zw(k) >= mean_surface_level_height(sr) ) THEN3676 k_surface_level = k3677 EXIT3678 ENDIF3679 ENDDO3680 3681 !3682 ! Computation of both the characteristic vertical velocity and3683 ! the characteristic convective boundary layer temperature.3684 ! The inversion height entering into the equation is defined with respect3685 ! to the mean surface level height of the respective statistic region.3686 ! The horizontal average at surface level index + 1 is input for the3687 ! average temperature.3688 IF ( hom(nzb,1,18,sr) > 1.0E8_wp .AND. z_i(1) /= 0.0_wp ) THEN3689 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_wp3693 ELSE3694 hom(nzb+8,1,pr_palm,sr) = 0.0_wp3695 ENDIF3696 3697 !3698 ! Collect the time series quantities3699 ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr) ! E3700 ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr) ! E*3701 ts_value(3,sr) = dt_3d3702 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_max3705 ts_value(7,sr) = v_max3706 ts_value(8,sr) = w_max3707 ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr) ! new divergence3708 ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr) ! old Divergence3709 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=03713 ts_value(15,sr) = hom(nzb+1,1,16,sr) ! w'pt' at k=13714 ts_value(16,sr) = hom(nzb+1,1,18,sr) ! wpt at k=13715 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=03718 ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr) ! v'w' at k=03719 ts_value(21,sr) = hom(nzb,1,48,sr) ! w"q" at k=03720 3721 IF ( .NOT. neutral ) THEN3722 ts_value(22,sr) = hom(nzb,1,114,sr) ! L3723 ELSE3724 ts_value(22,sr) = 1.0E10_wp3725 ENDIF3726 3727 ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr) ! q*3728 3729 !3730 ! Collect land surface model timeseries3731 IF ( land_surface ) THEN3732 ts_value(dots_soil ,sr) = hom(nzb,1,93,sr) ! ghf_eb3733 ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr) ! shf_eb3734 ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr) ! qsws_eb3735 ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr) ! qsws_liq_eb3736 ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr) ! qsws_soil_eb3737 ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr) ! qsws_veg_eb3738 ts_value(dots_soil+6,sr) = hom(nzb,1,99,sr) ! r_a3739 ts_value(dots_soil+7,sr) = hom(nzb,1,100,sr) ! r_s3740 ENDIF3741 !3742 ! Collect radiation model timeseries3743 IF ( radiation ) THEN3744 ts_value(dots_rad,sr) = hom(nzb,1,101,sr) ! rad_net3745 ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr) ! rad_lw_in3746 ts_value(dots_rad+2,sr) = hom(nzb,1,103,sr) ! rad_lw_out3747 ts_value(dots_rad+3,sr) = hom(nzb,1,104,sr) ! rad_sw_in3748 ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr) ! rad_sw_out3749 3750 IF ( radiation_scheme == 'rrtmg' ) THEN3751 ts_value(dots_rad+5,sr) = hom(nzb,1,106,sr) ! rrtm_aldif3752 ts_value(dots_rad+6,sr) = hom(nzb,1,107,sr) ! rrtm_aldir3753 ts_value(dots_rad+7,sr) = hom(nzb,1,108,sr) ! rrtm_asdif3754 ts_value(dots_rad+8,sr) = hom(nzb,1,109,sr) ! rrtm_asdir3755 ENDIF3756 3757 ENDIF3758 3759 !3760 ! Calculate additional statistics provided by the user interface3761 CALL user_statistics( 'time_series', sr, 0 )3762 3763 ENDDO ! loop of the subregions3764 3765 !$acc end data3766 3767 !3768 ! If required, sum up horizontal averages for subsequent time averaging3769 ! Do not sum, if flow statistics is called before the first initial time step.3770 IF ( do_sum .AND. simulated_time /= 0.0_wp ) THEN3771 IF ( average_count_pr == 0 ) hom_sum = 0.0_wp3772 hom_sum = hom_sum + hom(:,1,:,:)3773 average_count_pr = average_count_pr + 13774 do_sum = .FALSE.3775 ENDIF3776 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_statistics3786 #endif
Note: See TracChangeset
for help on using the changeset viewer.