Changeset 4472 for palm/trunk/SOURCE/flow_statistics.f90
- Timestamp:
- Mar 24, 2020 12:21:00 PM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/flow_statistics.f90 (modified) (47 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/flow_statistics.f90
r4464 r4472 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Calculations of the Kolmogorov lengt scale eta implemented 28 ! 29 ! 4464 2020-03-17 11:08:46Z Giersch 27 30 ! Reset last change (r4463) 28 ! 31 ! 29 32 ! 4463 2020-03-17 09:27:36Z Giersch 30 33 ! Calculate horizontally averaged profiles of all velocity components at the … … 33 36 ! 4444 2020-03-05 15:59:50Z raasch 34 37 ! bugfix: cpp-directives for serial mode added 35 ! 38 ! 36 39 ! 4442 2020-03-04 19:21:13Z suehring 37 ! Change order of dimension in surface array %frac to allow for better 40 ! Change order of dimension in surface array %frac to allow for better 38 41 ! vectorization. 39 ! 42 ! 40 43 ! 4441 2020-03-04 19:20:35Z suehring 41 44 ! Introduction of wall_flags_total_0, which currently sets bits based on static 42 45 ! topography information used in wall_flags_static_0 43 ! 46 ! 44 47 ! 4329 2019-12-10 15:46:36Z motisi 45 48 ! Renamed wall_flags_0 to wall_flags_static_0 46 ! 49 ! 47 50 ! 4182 2019-08-22 15:20:23Z scharf 48 51 ! Corrected "Former revisions" section 49 ! 52 ! 50 53 ! 4131 2019-08-02 11:06:18Z monakurppa 51 54 ! Allow profile output for salsa variables. 52 ! 55 ! 53 56 ! 4039 2019-06-18 10:32:41Z suehring 54 ! Correct conversion to kinematic scalar fluxes in case of pw-scheme and 57 ! Correct conversion to kinematic scalar fluxes in case of pw-scheme and 55 58 ! statistic regions 56 ! 59 ! 57 60 ! 3828 2019-03-27 19:36:23Z raasch 58 61 ! unused variables removed 59 ! 62 ! 60 63 ! 3676 2019-01-16 15:07:05Z knoop 61 64 ! Bugfix, terminate OMP Parallel block … … 71 74 !> 72 75 !> @note For simplicity, nzb_s_inner and nzb_diff_s_inner are being used as a 73 !> lower vertical index for k-loops for all variables, although strictly 74 !> speaking the k-loops would have to be split up according to the staggered 75 !> grid. However, this implies no error since staggered velocity components 76 !> lower vertical index for k-loops for all variables, although strictly 77 !> speaking the k-loops would have to be split up according to the staggered 78 !> grid. However, this implies no error since staggered velocity components 76 79 !> are zero at the walls and inside buildings. 77 80 !------------------------------------------------------------------------------! … … 97 100 USE control_parameters, & 98 101 ONLY: air_chemistry, average_count_pr, cloud_droplets, do_sum, & 99 dt_3d, humidity, initializing_actions, land_surface,&100 la rge_scale_forcing, large_scale_subsidence, max_pr_user,&101 m essage_string, neutral, ocean_mode, passive_scalar,&102 simulated_time, simulated_time_at_begin,&102 dt_3d, humidity, initializing_actions, kolmogorov_length_scale,& 103 land_surface, large_scale_forcing, large_scale_subsidence, & 104 max_pr_user, message_string, neutral, ocean_mode, & 105 passive_scalar, simulated_time, simulated_time_at_begin, & 103 106 use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, & 104 107 ws_scheme_mom, ws_scheme_sca, salsa, max_pr_salsa … … 109 112 USE grid_variables, & 110 113 ONLY: ddx, ddy 111 114 112 115 USE indices, & 113 116 ONLY: ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, nxl, nxr, nyn, & … … 118 121 ONLY: ngp_sums, ngp_sums_ls 119 122 #endif 120 123 121 124 USE kinds 122 125 123 126 USE land_surface_model_mod, & 124 127 ONLY: m_soil_h, nzb_soil, nzt_soil, t_soil_h … … 151 154 INTEGER(iwp) :: j !< 152 155 INTEGER(iwp) :: k !< 153 INTEGER(iwp) :: ki !< 156 INTEGER(iwp) :: ki !< 154 157 INTEGER(iwp) :: k_surface_level !< 155 INTEGER(iwp) :: m !< loop variable over all horizontal wall elements 158 INTEGER(iwp) :: m !< loop variable over all horizontal wall elements 156 159 INTEGER(iwp) :: l !< loop variable over surface facing -- up- or downward-facing 157 160 INTEGER(iwp) :: nt !< … … 161 164 162 165 LOGICAL :: first !< 163 164 REAL(wp) :: dptdz_threshold !< 166 167 REAL(wp) :: dissipation !< dissipation rate 168 REAL(wp) :: dptdz_threshold !< 169 REAL(wp) :: du_dx !< Derivative of u fluctuations with respect to x 170 REAL(wp) :: du_dy !< Derivative of u fluctuations with respect to y 171 REAL(wp) :: du_dz !< Derivative of u fluctuations with respect to z 172 REAL(wp) :: dv_dx !< Derivative of v fluctuations with respect to x 173 REAL(wp) :: dv_dy !< Derivative of v fluctuations with respect to y 174 REAL(wp) :: dv_dz !< Derivative of v fluctuations with respect to z 175 REAL(wp) :: dw_dx !< Derivative of w fluctuations with respect to x 176 REAL(wp) :: dw_dy !< Derivative of w fluctuations with respect to y 177 REAL(wp) :: dw_dz !< Derivative of w fluctuations with respect to z 178 REAL(wp) :: eta !< Kolmogorov length scale 165 179 REAL(wp) :: fac !< 166 180 REAL(wp) :: flag !< 167 181 REAL(wp) :: height !< 168 182 REAL(wp) :: pts !< 183 REAL(wp) :: s11 !< fluctuating rate-of-strain tensor component 11 184 REAL(wp) :: s21 !< fluctuating rate-of-strain tensor component 21 185 REAL(wp) :: s31 !< fluctuating rate-of-strain tensor component 31 186 REAL(wp) :: s12 !< fluctuating rate-of-strain tensor component 12 187 REAL(wp) :: s22 !< fluctuating rate-of-strain tensor component 22 188 REAL(wp) :: s32 !< fluctuating rate-of-strain tensor component 32 189 REAL(wp) :: s13 !< fluctuating rate-of-strain tensor component 13 190 REAL(wp) :: s23 !< fluctuating rate-of-strain tensor component 23 191 REAL(wp) :: s33 !< fluctuating rate-of-strain tensor component 33 169 192 REAL(wp) :: sums_l_etot !< 170 193 REAL(wp) :: ust !< … … 175 198 REAL(wp) :: v2 !< 176 199 REAL(wp) :: w2 !< 177 200 178 201 REAL(wp) :: dptdz(nzb+1:nzt+1) !< 179 202 REAL(wp) :: sums_ll(nzb:nzt+1,2) !< … … 210 233 !-- array 211 234 sums_l(:,11,:) = sums_l_l(:,sr,:) ! mixing length from diffusivities 212 !-- WARNING: next line still has to be adjusted for OpenMP 235 !-- WARNING: next line still has to be adjusted for OpenMP 213 236 sums_l(:,21,0) = sums_wsts_bc_l(:,sr) * & 214 237 heatflux_output_conversion ! heat flux from advec_s_bc … … 220 243 !-- scale) vertical fluxes and velocity variances by using commonly- 221 244 !-- applied Reynolds-based methods ( e.g. <w'pt'> = (w-<w>)*(pt-<pt>) ) 222 !-- in combination with the 5th order advection scheme, pronounced 223 !-- artificial kinks could be observed in the vertical profiles near the 224 !-- surface. Please note: these kinks were not related to the model truth, 225 !-- i.e. these kinks are just related to an evaluation problem. 226 !-- In order avoid these kinks, vertical fluxes and horizontal as well 245 !-- in combination with the 5th order advection scheme, pronounced 246 !-- artificial kinks could be observed in the vertical profiles near the 247 !-- surface. Please note: these kinks were not related to the model truth, 248 !-- i.e. these kinks are just related to an evaluation problem. 249 !-- In order avoid these kinks, vertical fluxes and horizontal as well 227 250 !-- vertical velocity variances are calculated directly within the advection 228 !-- routines, according to the numerical discretization, to evaluate the 229 !-- statistical quantities as they will appear within the prognostic 251 !-- routines, according to the numerical discretization, to evaluate the 252 !-- statistical quantities as they will appear within the prognostic 230 253 !-- equations. 231 !-- Copy the turbulent quantities, evaluated in the advection routines to 254 !-- Copy the turbulent quantities, evaluated in the advection routines to 232 255 !-- the local array sums_l() for further computations. 233 256 IF ( ws_scheme_mom .AND. sr == 0 ) THEN … … 248 271 sums_l(:,15,i) = sums_wsvs_ws_l(:,i) & 249 272 * momentumflux_output_conversion ! w*v* 250 sums_l(:,30,i) = sums_us2_ws_l(:,i) ! u*2 251 sums_l(:,31,i) = sums_vs2_ws_l(:,i) ! v*2 252 sums_l(:,32,i) = sums_ws2_ws_l(:,i) ! w*2 253 sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp * & 273 sums_l(:,30,i) = sums_us2_ws_l(:,i) ! u*2 274 sums_l(:,31,i) = sums_vs2_ws_l(:,i) ! v*2 275 sums_l(:,32,i) = sums_ws2_ws_l(:,i) ! w*2 276 sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp * & 254 277 ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) + & 255 278 sums_ws2_ws_l(:,i) ) ! e* … … 270 293 271 294 ENDIF 272 ! 295 ! 273 296 !-- Horizontally averaged profiles of horizontal velocities and temperature. 274 297 !-- They must have been computed before, because they are already required … … 446 469 447 470 ! 448 !-- Final values are obtained by division by the total number of grid points 471 !-- Final values are obtained by division by the total number of grid points 449 472 !-- used for summation. After that store profiles. 450 473 sums(:,1) = sums(:,1) / ngp_2dh(sr) … … 482 505 !-- Passive scalar 483 506 IF ( passive_scalar ) hom(:,1,115,sr) = sums(:,115) / & 484 ngp_2dh_s_inner(:,sr) ! s 507 ngp_2dh_s_inner(:,sr) ! s 485 508 486 509 ! … … 488 511 !-- variances, the total and the perturbation energy (single values in last 489 512 !-- column of sums_l) and some diagnostic quantities. 490 !-- NOTE: for simplicity, nzb_s_inner is used below, although strictly 491 !-- ---- speaking the following k-loop would have to be split up and 513 !-- NOTE: for simplicity, nzb_s_inner is used below, although strictly 514 !-- ---- speaking the following k-loop would have to be split up and 492 515 !-- rearranged according to the staggered grid. 493 516 !-- However, this implies no error since staggered velocity components … … 500 523 !$OMP DO 501 524 !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, m) & 502 !$ACC PRIVATE(sums_l_etot, flag) & 525 !$ACC PRIVATE(sums_l_etot, flag, du_dx, du_dy, du_dz) & 526 !$ACC PRIVATE(dv_dx, dv_dy, dv_dz, dw_dx, dw_dy, dw_dz) & 527 !$ACC PRIVATE(s11, s21, s31, s12, s22, s32, s13, s23, s33) & 528 !$ACC PRIVATE(dissipation, eta) & 503 529 !$ACC PRESENT(wall_flags_total_0, rmask, momentumflux_output_conversion) & 504 !$ACC PRESENT(hom(:,1, 4,sr)) &530 !$ACC PRESENT(hom(:,1,1:2,sr), hom(:,1,4,sr)) & 505 531 !$ACC PRESENT(e, u, v, w, km, kh, p, pt) & 532 !$ACC PRESENT(ddx, ddy, ddzu, ddzw) & 506 533 !$ACC PRESENT(surf_def_h(0), surf_lsm_h, surf_usm_h) & 507 534 !$ACC PRESENT(sums_l) … … 557 584 w(k,j,i)**2 ) * rmask(j,i,sr)& 558 585 * flag 559 ENDDO 586 587 ! 588 !-- Computation of the Kolmogorov length scale. Calculation is based 589 !-- on gradients of the deviations from the horizontal mean. 590 !-- Kolmogorov scale at the boundaries (k=0/z=0m and k=nzt+1) is set to zero. 591 IF ( kolmogorov_length_scale .AND. k /= nzb .AND. k /= nzt+1) THEN 592 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 22 ) ) 593 594 ! 595 !-- Calculate components of the fluctuating rate-of-strain tensor 596 !-- (0.5*(del u'_i/del x_j + del u'_j/del x_i)) defined in the 597 !-- center of each grid box. 598 du_dx = ( ( u(k,j,i+1) - hom(k,1,1,sr) ) - & 599 ( u(k,j,i) - hom(k,1,1,sr) ) ) * ddx 600 du_dy = 0.25_wp * ddy * & 601 ( ( u(k,j+1,i) - hom(k,1,1,sr) ) - & 602 ( u(k,j-1,i) - hom(k,1,1,sr) ) + & 603 ( u(k,j+1,i+1) - hom(k,1,1,sr) ) - & 604 ( u(k,j-1,i+1) - hom(k,1,1,sr) ) ) 605 du_dz = 0.25_wp * ( ( ( u(k+1,j,i) - hom(k+1,1,1,sr) ) - & 606 ( u(k,j,i) - hom(k,1,1,sr) ) ) * & 607 ddzu(k+1) + & 608 ( ( u(k,j,i) - hom(k,1,1,sr) ) - & 609 ( u(k-1,j,i) - hom(k-1,1,1,sr) ) )* & 610 ddzu(k) + & 611 ( ( u(k+1,j,i+1) - hom(k+1,1,1,sr) )- & 612 ( u(k,j,i+1) - hom(k,1,1,sr) ) ) * & 613 ddzu(k+1) + & 614 ( ( u(k,j,i+1) - hom(k,1,1,sr) ) - & 615 ( u(k-1,j,i+1) - hom(k-1,1,1,sr) ) ) *& 616 ddzu(k) ) 617 618 dv_dx = 0.25_wp * ddx * & 619 ( ( v(k,j,i+1) - hom(k,1,2,sr) ) - & 620 ( v(k,j,i-1) - hom(k,1,2,sr) ) + & 621 ( v(k,j+1,i+1) - hom(k,1,2,sr) ) - & 622 ( v(k,j+1,i-1) - hom(k,1,2,sr) ) ) 623 dv_dy = ( ( v(k,j+1,i) - hom(k,1,2,sr) ) - & 624 ( v(k,j,i) - hom(k,1,2,sr) ) ) * ddy 625 dv_dz = 0.25_wp * ( ( ( v(k+1,j,i) - hom(k+1,1,2,sr) ) - & 626 ( v(k,j,i) - hom(k,1,2,sr) ) ) * & 627 ddzu(k+1) + & 628 ( ( v(k,j,i) - hom(k,1,2,sr) ) - & 629 ( v(k-1,j,i) - hom(k-1,1,2,sr) ) ) * & 630 ddzu(k) + & 631 ( ( v(k+1,j+1,i) - hom(k+1,1,2,sr) ) - & 632 ( v(k,j+1,i) - hom(k,1,2,sr) ) ) * & 633 ddzu(k+1) + & 634 ( ( v(k,j+1,i) - hom(k,1,2,sr) ) - & 635 ( v(k-1,j+1,i) - hom(k-1,1,2,sr) ) ) *& 636 ddzu(k) ) 637 638 dw_dx = 0.25_wp * ddx * ( w(k,j,i+1) - w(k,j,i-1) + & 639 w(k-1,j,i+1) - w(k-1,j,i-1) ) 640 dw_dy = 0.25_wp * ddy * ( w(k,j+1,i) - w(k,j-1,i) + & 641 w(k-1,j+1,i) - w(k-1,j-1,i) ) 642 dw_dz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 643 644 s11 = 0.5_wp * ( du_dx + du_dx ) 645 s21 = 0.5_wp * ( dv_dx + du_dy ) 646 s31 = 0.5_wp * ( dw_dx + du_dz ) 647 648 s12 = s21 649 s22 = 0.5 * ( dv_dy + dv_dy ) 650 s32 = 0.5 * ( dw_dy + dv_dz ) 651 652 s13 = s31 653 s23 = s32 654 s33 = 0.5_wp * ( dw_dz + dw_dz ) 655 656 !-- Calculate 3D instantaneous energy dissipation rate after 657 !-- Pope (2000): Turbulent flows, p.259. It is defined in the center 658 !-- of each grid volume. 659 dissipation = 2.0_wp * km(k,j,i) * & 660 ( s11*s11 + s21*s21 + s31*s31 + & 661 s12*s12 + s22*s22 + s32*s32 + & 662 s13*s13 + s23*s23 + s33*s33 ) 663 eta = ( km(k,j,i)**3.0_wp / ( dissipation+1.0E-12 ) )**(1.0_wp/4.0_wp) 664 665 !$ACC ATOMIC 666 sums_l(k,121,tn) = sums_l(k,121,tn) + eta * rmask(j,i,sr) & 667 * flag 668 669 670 ENDIF !Kolmogorov length scale 671 672 ENDDO !k-loop 560 673 ! 561 674 !-- Total and perturbation energy for the total domain (being … … 661 774 rmask(j,i,sr) 662 775 ENDIF 663 ENDDO 664 ENDDO 776 ENDDO !j-loop 777 ENDDO !i-loop 665 778 !$ACC UPDATE & 666 779 !$ACC HOST(sums_l(:,3,tn), sums_l(:,8,tn), sums_l(:,9,tn)) & 667 780 !$ACC HOST(sums_l(:,10,tn), sums_l(:,40,tn), sums_l(:,33,tn)) & 668 !$ACC HOST(sums_l(:,38,tn) ) &781 !$ACC HOST(sums_l(:,38,tn), sums_l(:,121,tn)) & 669 782 !$ACC HOST(sums_l(nzb:nzb+4,pr_palm,tn), sums_l(nzb+14:nzb+14,pr_palm,tn)) 670 783 671 784 ! 672 !-- Computation of statistics when ws-scheme is not used. Else these 785 !-- Computation of statistics when ws-scheme is not used. Else these 673 786 !-- quantities are evaluated in the advection routines. 674 787 IF ( .NOT. ws_scheme_mom .OR. sr /= 0 .OR. simulated_time == 0.0_wp ) & … … 703 816 ENDIF 704 817 ! 705 !-- Computaion of domain-averaged perturbation energy. Please note, 706 !-- to prevent that perturbation energy is larger (even if only slightly) 818 !-- Computaion of domain-averaged perturbation energy. Please note, 819 !-- to prevent that perturbation energy is larger (even if only slightly) 707 820 !-- than the total kinetic energy, calculation is based on deviations from 708 821 !-- the horizontal mean, instead of spatial descretization of the advection 709 !-- term. 822 !-- term. 710 823 !$OMP DO 711 824 !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k, flag, w2, ust2, vst2) & … … 727 840 * rmask(j,i,sr) & 728 841 * flag 842 729 843 ENDDO 730 844 ENDDO … … 746 860 DO j = nys, nyn 747 861 ! 748 !-- Subgridscale fluxes (without Prandtl layer from k=nzb, 862 !-- Subgridscale fluxes (without Prandtl layer from k=nzb, 749 863 !-- oterwise from k=nzb+1) 750 864 !-- NOTE: for simplicity, nzb_diff_s_inner is used below, although 751 865 !-- ---- strictly speaking the following k-loop would have to be 752 866 !-- split up according to the staggered grid. 753 !-- However, this implies no error since staggered velocity 867 !-- However, this implies no error since staggered velocity 754 868 !-- components are zero at the walls and inside buildings. 755 869 !-- Flag 23 is used to mask surface fluxes as well as model-top fluxes, 756 !-- which are added further below. 870 !-- which are added further below. 757 871 DO k = nzb, nzt 758 872 flag = MERGE( 1.0_wp, 0.0_wp, & … … 1184 1298 0.0_wp * rmask(j,i,sr) ! u"pt" 1185 1299 sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + & 1186 0.0_wp * rmask(j,i,sr) ! v"pt" 1300 0.0_wp * rmask(j,i,sr) ! v"pt" 1187 1301 #endif 1188 1302 #ifndef _OPENACC … … 1227 1341 ! 1228 1342 !-- Resolved fluxes (can be computed for all horizontal points) 1229 !-- NOTE: for simplicity, nzb_s_inner is used below, although strictly 1230 !-- ---- speaking the following k-loop would have to be split up and 1343 !-- NOTE: for simplicity, nzb_s_inner is used below, although strictly 1344 !-- ---- speaking the following k-loop would have to be split up and 1231 1345 !-- rearranged according to the staggered grid. 1232 1346 DO k = nzb, nzt … … 1303 1417 sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) * & 1304 1418 rmask(j,i,sr) *& 1305 flag 1419 flag 1306 1420 sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) * & 1307 1421 rmask(j,i,sr) *& … … 1320 1434 flag 1321 1435 ELSE IF ( ws_scheme_sca .AND. sr == 0 ) THEN 1322 sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp * & 1436 sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp * & 1323 1437 hom(k,1,41,sr) ) * & 1324 1438 sums_l(k,17,tn) + & … … 1359 1473 !$ACC HOST(sums_l(:,35,tn), sums_l(:,36,tn), sums_l(:,37,tn)) 1360 1474 ! 1361 !-- Treat land-surface quantities according to new wall model structure. 1475 !-- Treat land-surface quantities according to new wall model structure. 1362 1476 IF ( land_surface ) THEN 1363 1477 tn = 0 … … 1368 1482 i = surf_lsm_h%i(m) 1369 1483 j = surf_lsm_h%j(m) 1370 1484 1371 1485 IF ( i >= nxl .AND. i <= nxr .AND. & 1372 j >= nys .AND. j <= nyn ) THEN 1486 j >= nys .AND. j <= nyn ) THEN 1373 1487 sums_l(nzb,93,tn) = sums_l(nzb,93,tn) + surf_lsm_h%ghf(m) 1374 1488 sums_l(nzb,94,tn) = sums_l(nzb,94,tn) + surf_lsm_h%qsws_liq(m) … … 1387 1501 DO m = 1, surf_lsm_h%ns 1388 1502 1389 i = surf_lsm_h%i(m) 1503 i = surf_lsm_h%i(m) 1390 1504 j = surf_lsm_h%j(m) 1391 1505 1392 1506 IF ( i >= nxl .AND. i <= nxr .AND. & 1393 j >= nys .AND. j <= nyn ) THEN 1507 j >= nys .AND. j <= nyn ) THEN 1394 1508 1395 1509 DO k = nzb_soil, nzt_soil … … 1417 1531 DO k = nzb, nzt 1418 1532 ! 1419 !-- Flag 23 is used to mask surface fluxes as well as model-top 1420 !-- fluxes, which are added further below. 1533 !-- Flag 23 is used to mask surface fluxes as well as model-top 1534 !-- fluxes, which are added further below. 1421 1535 flag = MERGE( 1.0_wp, 0.0_wp, & 1422 1536 BTEST( wall_flags_total_0(k,j,i), 23 ) ) * & … … 1569 1683 1570 1684 ! 1571 !-- Horizontal heat fluxes (subgrid, resolved, total). 1572 !-- Do it only, if profiles shall be plotted. 1685 !-- Horizontal heat fluxes (subgrid, resolved, total). 1686 !-- Do it only, if profiles shall be plotted. 1573 1687 IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN 1574 1688 … … 1628 1742 IF ( large_scale_forcing .AND. ( simulated_time > 0.0_wp ) ) THEN 1629 1743 ! 1630 !-- Interpolation in time of LSF_DATA 1744 !-- Interpolation in time of LSF_DATA 1631 1745 nt = 1 1632 1746 DO WHILE ( simulated_time - dt_3d > time_vert(nt) ) … … 1669 1783 tn = 0 1670 1784 !$OMP PARALLEL PRIVATE( i, j, k, tn ) 1671 !$ tn = omp_get_thread_num() 1785 !$ tn = omp_get_thread_num() 1672 1786 IF ( radiation .AND. radiation_scheme == 'rrtmg' ) THEN 1673 1787 !$OMP DO … … 1718 1832 1719 1833 IF ( air_chemistry ) THEN 1720 IF ( max_pr_cs > 0 ) THEN 1834 IF ( max_pr_cs > 0 ) THEN 1721 1835 sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+ max_pr_cs,0) = & 1722 1836 sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+max_pr_cs,0) + & … … 1774 1888 1775 1889 ! 1776 !-- Final values are obtained by division by the total number of grid points 1890 !-- Final values are obtained by division by the total number of grid points 1777 1891 !-- used for summation. After that store profiles. 1778 1892 !-- Check, if statistical regions do contain at least one grid point at the … … 1834 1948 1835 1949 IF ( air_chemistry ) THEN 1836 IF ( max_pr_cs > 0 ) THEN 1950 IF ( max_pr_cs > 0 ) THEN 1837 1951 DO k = nzb, nzt+1 1838 1952 sums(k, pr_palm+1:pr_palm+max_pr_user+max_pr_cs) = & … … 1840 1954 ngp_2dh_s_inner(k,sr) 1841 1955 ENDDO 1842 ENDIF 1956 ENDIF 1843 1957 ENDIF 1844 1958 … … 1850 1964 / ngp_2dh_s_inner(k,sr) 1851 1965 ENDDO 1852 ENDIF 1966 ENDIF 1853 1967 ENDIF 1854 1968 … … 1873 1987 hom(:,1,22,sr) = sums(:,16) + sums(:,21) ! wptBC 1874 1988 ! profile 24 is initial profile (sa) 1875 ! profiles 25-29 left empty for initial 1989 ! profiles 25-29 left empty for initial 1876 1990 ! profiles 1877 1991 hom(:,1,30,sr) = sums(:,30) ! u*2 … … 1887 2001 hom(:,1,40,sr) = sums(:,40) ! p 1888 2002 hom(:,1,45,sr) = sums(:,45) ! w"vpt" 1889 hom(:,1,46,sr) = sums(:,46) ! w*vpt* 2003 hom(:,1,46,sr) = sums(:,46) ! w*vpt* 1890 2004 hom(:,1,47,sr) = sums(:,45) + sums(:,46) ! wvpt 1891 2005 hom(:,1,48,sr) = sums(:,48) ! w"q" (w"qv") … … 1893 2007 hom(:,1,50,sr) = sums(:,48) + sums(:,49) ! wq (wqv) 1894 2008 hom(:,1,51,sr) = sums(:,51) ! w"qv" 1895 hom(:,1,52,sr) = sums(:,52) ! w*qv* 2009 hom(:,1,52,sr) = sums(:,52) ! w*qv* 1896 2010 hom(:,1,53,sr) = sums(:,52) + sums(:,51) ! wq (wqv) 1897 2011 hom(:,1,54,sr) = sums(:,54) ! ql … … 1979 2093 hom(:,1,117,sr) = sums(:,117) ! w"s" 1980 2094 hom(:,1,114,sr) = sums(:,114) ! w*s* 1981 hom(:,1,118,sr) = sums(:,117) + sums(:,114) ! ws 2095 hom(:,1,118,sr) = sums(:,117) + sums(:,114) ! ws 1982 2096 hom(:,1,116,sr) = sums(:,116) ! s*2 1983 2097 ENDIF … … 1985 2099 hom(:,1,119,sr) = rho_air ! rho_air in Kg/m^3 1986 2100 hom(:,1,120,sr) = rho_air_zw ! rho_air_zw in Kg/m^3 2101 2102 IF ( kolmogorov_length_scale ) THEN 2103 hom(:,1,121,sr) = sums(:,121) * 1E3_wp ! eta in mm 2104 ENDIF 2105 1987 2106 1988 2107 hom(:,1,pr_palm,sr) = sums(:,pr_palm) … … 1995 2114 1996 2115 IF ( air_chemistry ) THEN 1997 IF ( max_pr_cs > 0 ) THEN ! chem_spcs profiles 2116 IF ( max_pr_cs > 0 ) THEN ! chem_spcs profiles 1998 2117 hom(:, 1, pr_palm+max_pr_user+1:pr_palm + max_pr_user+max_pr_cs, sr) = & 1999 2118 sums(:, pr_palm+max_pr_user+1:pr_palm+max_pr_user+max_pr_cs) … … 2013 2132 !-- The corresponding height is assumed as the boundary layer height, if it 2014 2133 !-- is less than 1.5 times the height where the heat flux becomes negative 2015 !-- (positive) for the first time. Attention: the resolved vertical sensible 2134 !-- (positive) for the first time. Attention: the resolved vertical sensible 2016 2135 !-- heat flux (hom(:,1,17,sr) = w*pt*) is not known at the beginning because 2017 !-- the calculation happens in advec_s_ws which is called after 2018 !-- flow_statistics. Therefore z_i is directly taken from restart data at 2019 !-- the beginning of restart runs. 2136 !-- the calculation happens in advec_s_ws which is called after 2137 !-- flow_statistics. Therefore z_i is directly taken from restart data at 2138 !-- the beginning of restart runs. 2020 2139 IF ( TRIM( initializing_actions ) /= 'read_restart_data' .OR. & 2021 2140 simulated_time_at_begin /= simulated_time ) THEN … … 2059 2178 2060 2179 ! 2061 !-- Second scheme: Gradient scheme from Sullivan et al. (1998), modified 2062 !-- by Uhlenbrock(2006). The boundary layer height is the height with the 2180 !-- Second scheme: Gradient scheme from Sullivan et al. (1998), modified 2181 !-- by Uhlenbrock(2006). The boundary layer height is the height with the 2063 2182 !-- maximal local temperature gradient: starting from the second (the 2064 !-- last but one) vertical gridpoint, the local gradient must be at least 2183 !-- last but one) vertical gridpoint, the local gradient must be at least 2065 2184 !-- 0.2K/100m and greater than the next four gradients. 2066 2185 !-- WARNING: The threshold value of 0.2K/100m must be adjusted for the 2067 !-- ocean case! 2186 !-- ocean case! 2068 2187 z_i(2) = 0.0_wp 2069 2188 DO k = nzb+1, nzt+1 … … 2126 2245 2127 2246 ! 2128 !-- Collect the time series quantities. Please note, timeseries quantities 2129 !-- which are collected from horizontally averaged profiles, e.g. wpt 2130 !-- or pt(zp), are treated specially. In case of elevated model surfaces, 2247 !-- Collect the time series quantities. Please note, timeseries quantities 2248 !-- which are collected from horizontally averaged profiles, e.g. wpt 2249 !-- or pt(zp), are treated specially. In case of elevated model surfaces, 2131 2250 !-- index nzb+1 might be within topography and data will be zero. Therefore, 2132 !-- take value for the first atmosphere index, which is topo_min_level+1. 2251 !-- take value for the first atmosphere index, which is topo_min_level+1. 2133 2252 ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr) ! E 2134 2253 ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr) ! E*
Note: See TracChangeset
for help on using the changeset viewer.
