From ae344fb0afdfe55b5232666c157e3c7fb4037cb7 Mon Sep 17 00:00:00 2001 From: MetBenjaminWent Date: Wed, 21 Jan 2026 13:15:08 +0000 Subject: [PATCH 1/5] Hand adjusted, remove j, good build --- .../source/boundary_layer/bdy_expl2.F90 | 2627 +++++++++-------- 1 file changed, 1359 insertions(+), 1268 deletions(-) diff --git a/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 b/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 index af0136c0..cbf26cb0 100644 --- a/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 +++ b/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 @@ -19,13 +19,12 @@ module bdy_expl2_mod use fm_drag_mod, only: fm_drag -use tuning_segments_mod, only: & - l_autotune_segments, & - bl_segment_size - +use tuning_segments_mod, only: bl_segment_size use um_types, only: r_bl +use timer_mod, only: timer + implicit none !Automatic segment size tuning @@ -634,8 +633,9 @@ subroutine bdy_expl2 ( & ! (qw - qsat(Tl))/(1 + Lc/cp dqsat/dT) ! gradient of this is used for estimating ! saturated fraction on rho-levels - qs_tl(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end) + qs_tl(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,1:bl_levels) ! qsat(Tl) on current level + ! Vectorised over levels (k) for open-mp ease real(kind=r_bl) :: & frac_sat, frac_dry, frac_edg, frac_lev, qc_tot, bt_rh, bq_rh @@ -781,7 +781,7 @@ subroutine bdy_expl2 ( & real(kind=r_bl) :: & b2, sh, exner, root6, delta_x, var_fac, sl_var, qw_var, sl_qw, & sgm(tdims%i_start:tdims%i_end), & - qsw_arr(tdims%i_start:tdims%i_end), & + qsw_arr(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,1:bl_levels-1), & max_rhcpt(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end), & min_rhcpt(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end) @@ -800,10 +800,17 @@ subroutine bdy_expl2 ( & integer :: & - omp_block, & - ! for open mp blocking - jj - ! for indexing over open mp block + pdims_omp_block, pdims_seg_block, & + ! for pdims open-mp blocking + tdims_omp_block, tdims_seg_block, & + ! for tdims open mp blocking + ii, & + ! for indexing over open-mp block + max_threads, & + ! Hold the number of threads for open-mp + seg_slice_start, seg_slice_end, & + ! Segmentaion of calls into qstat, start and end + bl_segment_range ! Range of segmentaion, derived from start and end real(kind=r_bl), parameter :: max_abs_obkhov = 1.0e6_r_bl ! Maximum permitted magnitude of the Obukhov @@ -824,40 +831,43 @@ subroutine bdy_expl2 ( & real(kind=jprb) :: zhook_handle if (lhook) call dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) +call timer('bdy_expl2') -!Set up automatic segment tuning +! Set up blocking for OMP sections +max_threads = 1 +!$ max_threads = omp_get_max_threads() +j = 1 !----------------------------------------------------------------------- ! 1) Set various diagnostics and switches !----------------------------------------------------------------------- ! checking of nl_bl_levels has been moved to readsize/scm_shell ! so that it is only executed at initialisation - +!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL DEFAULT(none) & -!$OMP SHARED( pdims, dynamic_bl_diag, ntml_save, ntml, bl_levels, & +!$OMP SHARED( j, pdims, dynamic_bl_diag, ntml_save, ntml, bl_levels, & !$OMP weight_1dbl, weight_1dbl_rho, BL_diag, u_s, fb_surf ) & -!$OMP private( i, j, k ) +!$OMP private( i, k ) !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end do i = pdims%i_start, pdims%i_end dynamic_bl_diag(i,j) = .false. ntml_save(i,j) = ntml(i,j) end do -end do !$OMP end do NOWAIT + !end do !------------------------------------------------------------------ ! Initialize weighting applied to 1d BL scheme ! (used to blend with 3D Smagorinsky scheme) !------------------------------------------------------------------ + !do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - weight_1dbl(i,j,k) = one - weight_1dbl_rho(i,j,k) = one - end do + do i = pdims%i_start, pdims%i_end + weight_1dbl(i,j,k) = one + weight_1dbl_rho(i,j,k) = one end do + ! end do end do !$OMP end do NOWAIT @@ -865,42 +875,42 @@ subroutine bdy_expl2 ( & ! Set surface scaling diagnostics !----------------------------------------------------------------------- ! Obukhov length + ! do j = pdims%j_start, pdims%j_end if (BL_diag%l_oblen) then !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ! Limit the magnitude of the Obukhov length to avoid - ! problems with packing. - BL_diag%oblen(i,j)= u_s(i,j)*u_s(i,j)*u_s(i,j) - if ( BL_diag%oblen(i,j) < & - max_abs_obkhov*abs(vkman*fb_surf(i,j)) ) then - BL_diag%oblen(i,j)=-BL_diag%oblen(i,j)/(vkman*fb_surf(i,j)) - else - BL_diag%oblen(i,j)=-sign(max_abs_obkhov, fb_surf(i,j)) - end if - end do + do i = pdims%i_start, pdims%i_end + ! Limit the magnitude of the Obukhov length to avoid + ! problems with packing. + BL_diag%oblen(i,j)= u_s(i,j)*u_s(i,j)*u_s(i,j) + if ( BL_diag%oblen(i,j) < & + max_abs_obkhov*abs(vkman*fb_surf(i,j)) ) then + BL_diag%oblen(i,j)=-BL_diag%oblen(i,j)/(vkman*fb_surf(i,j)) + else + BL_diag%oblen(i,j)=-sign(max_abs_obkhov, fb_surf(i,j)) + end if + ! end do end do !$OMP end do NOWAIT end if ! Ustar + !do j = pdims%j_start, pdims%j_end if (BL_diag%l_ustar) then !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - BL_diag%ustar(i,j)=u_s(i,j) - end do + do i = pdims%i_start, pdims%i_end + BL_diag%ustar(i,j)=u_s(i,j) + !end do end do !$OMP end do NOWAIT end if ! Surface buoyancy flux + !do j = pdims%j_start, pdims%j_end if (BL_diag%l_wbsurf) then !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - BL_diag%wbsurf(i,j)=fb_surf(i,j) - end do + do i = pdims%i_start, pdims%i_end + BL_diag%wbsurf(i,j)=fb_surf(i,j) + !end do end do !$OMP end do end if @@ -922,32 +932,40 @@ subroutine bdy_expl2 ( & !----------------------------------------------------------------------- ! Calculate lapse rates !----------------------------------------------------------------------- -!$OMP PARALLEL DEFAULT(SHARED) private(i, j, k, weight1, weight2, & +! Sub divide pdims domain len over the number of threads +pdims_omp_block = ceiling(real(pdims%i_len)/max_threads) +! Pick the smallest option. +! Work should be split accross the threads, this provides a cap. +! In the occurence where the pdims domain len is very small, it will be seleted. +! Otherwise bl_segment_size will be selected. +!do j = pdims%j_start, pdims%j_end +pdims_seg_block = min(bl_segment_size, pdims_omp_block, pdims%i_len) +!$OMP PARALLEL DEFAULT(SHARED) private(ii, i, k, weight1, weight2, & !$OMP weight3, zpr, dzv, dzu, l, slope, dsldzm_ga, & -!$OMP qs_tl, frac_sat, frac_dry, frac_edg, frac_lev, qc_tot, bt_rh, bq_rh ) +!$OMP qs_tl, frac_sat, frac_dry, frac_edg, frac_lev, qc_tot, bt_rh, bq_rh, & +!$OMP seg_slice_start, seg_slice_end, bl_segment_range) !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - grad_t_adj(i,j) = min( max_t_grad, & - a_grad_adj * t1_sd(i,j) / zh_prev(i,j) ) - end do +do i = pdims%i_start, pdims%i_end + grad_t_adj(i,j) = min( max_t_grad, & + a_grad_adj * t1_sd(i,j) / zh_prev(i,j) ) end do +!end do !$OMP end do + ! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - dsldz(i,j,k) = ( tl(i,j,k) - tl(i,j,k-1) ) & - *rdz_charney_grid(i,j,k) + grcp - dsldz_ga(i,j,k) = dsldz(i,j,k) - if ( z_tq(i,j,k) <= zh_prev(i,j) ) then - dsldz_ga(i,j,k) = dsldz_ga(i,j,k) - grad_t_adj(i,j) - end if - dqwdz(i,j,k) = ( qw(i,j,k) - qw(i,j,k-1) ) & - * rdz_charney_grid(i,j,k) - end do + do i = pdims%i_start, pdims%i_end + dsldz(i,j,k) = ( tl(i,j,k) - tl(i,j,k-1) ) & + *rdz_charney_grid(i,j,k) + grcp + dsldz_ga(i,j,k) = dsldz(i,j,k) + if ( z_tq(i,j,k) <= zh_prev(i,j) ) then + dsldz_ga(i,j,k) = dsldz_ga(i,j,k) - grad_t_adj(i,j) + end if + dqwdz(i,j,k) = ( qw(i,j,k) - qw(i,j,k-1) ) & + * rdz_charney_grid(i,j,k) end do + !end do end do !$OMP end do @@ -963,51 +981,71 @@ subroutine bdy_expl2 ( & ! extrapolate dbdz itself from level 2, but the sl and qw gradients are ! used in the newer variance calculations and with i_interp_local_cf_dbdz ! so extrapolate them here +! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - dsldz(i,j,1) = dsldz(i,j,2) - dsldz_ga(i,j,1) = dsldz_ga(i,j,2) - dqwdz(i,j,1) = dqwdz(i,j,2) - end do + do i = pdims%i_start, pdims%i_end + dsldz(i,j,1) = dsldz(i,j,2) + dsldz_ga(i,j,1) = dsldz_ga(i,j,2) + dqwdz(i,j,1) = dqwdz(i,j,2) end do + !end do !$OMP end do else ! l_use_surf_in_ri = true -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end +!! $OMP do SCHEDULE(STATIC) + !do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) + do ii = pdims%i_start, pdims%i_end, pdims_seg_block + seg_slice_start = ii + seg_slice_end = MIN(((ii+pdims_seg_block)-1), pdims%i_end) + bl_segment_range = (seg_slice_end - seg_slice_start) + 1 + if (l_noice_in_turb) then ! use qsat_wat if ( l_mr_physics ) then - call qsat_wat_mix(qssurf(:,j),tstar(:,j),pstar(:,j),pdims%i_len) - ! No halos + call qsat_wat_mix(qssurf(seg_slice_start:seg_slice_end,j), & + tstar(seg_slice_start:seg_slice_end,j), & + pstar(seg_slice_start:seg_slice_end,j), & + bl_segment_range) ! No halos else - call qsat_wat(qssurf(:,j),tstar(:,j),pstar(:,j),pdims%i_len) ! No halos + call qsat_wat(qssurf(seg_slice_start:seg_slice_end,j), & + tstar(seg_slice_start:seg_slice_end,j), & + pstar(seg_slice_start:seg_slice_end,j), & + bl_segment_range) ! No halos end if else ! l_noice_in_turb ! use qsat if ( l_mr_physics ) then - call qsat_mix(qssurf(:,j),tstar(:,j),pstar(:,j),pdims%i_len) ! No halos + call qsat_mix(qssurf(seg_slice_start:seg_slice_end,j), & + tstar(seg_slice_start:seg_slice_end,j), & + pstar(seg_slice_start:seg_slice_end,j), & + bl_segment_range) ! No halos else - call qsat(qssurf(:,j),tstar(:,j),pstar(:,j),pdims%i_len) ! No halos + call qsat(qssurf(seg_slice_start:seg_slice_end,j), & + tstar(seg_slice_start:seg_slice_end,j), & + pstar(seg_slice_start:seg_slice_end,j), & + bl_segment_range) ! No halos end if end if ! l_noice_in_turb - end do ! j -!$OMP end do NOWAIT + end do ! ii + !$OMP end do NOWAIT + !end do ! j +!! $OMP end do NOWAIT k=1 -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - dsldz(i,j,k) = ( tl(i,j,k) - tstar(i,j) ) & - *rdz_charney_grid(i,j,k) + grcp - dsldz_ga(i,j,k) = dsldz(i,j,k) ! no GA below level 1 - if ( flandg(i,j) < 0.2_r_bl) then - dqwdz(i,j,k) = ( qw(i,j,k) - qssurf(i,j) ) & - * rdz_charney_grid(i,j,k) - else - dqwdz(i,j,k) = dqwdz(i,j,2) ! extrapolate qw if mainly land - end if - end do + + ! do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + dsldz(i,j,k) = ( tl(i,j,k) - tstar(i,j) ) & + *rdz_charney_grid(i,j,k) + grcp + dsldz_ga(i,j,k) = dsldz(i,j,k) ! no GA below level 1 + if ( flandg(i,j) < 0.2_r_bl) then + dqwdz(i,j,k) = ( qw(i,j,k) - qssurf(i,j) ) & + * rdz_charney_grid(i,j,k) + else + dqwdz(i,j,k) = dqwdz(i,j,2) ! extrapolate qw if mainly land + end if + !end do end do !$OMP end do @@ -1022,33 +1060,32 @@ subroutine bdy_expl2 ( & ! Interpolate gradients to theta-levels, then ! calculate `buoyancy' gradient, DBDZ, on theta-levels ! NOTE: DBDZ(K) is on theta-level K-1 - + ! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - weight1 = r_rho_levels(i,j,k) - & - r_rho_levels(i,j,k-1) - weight2 = r_theta_levels(i,j,k-1)- & - r_rho_levels(i,j,k-1) - weight3 = r_rho_levels(i,j,k) - & - r_theta_levels(i,j,k-1) - dsldzm(i,j,k) = weight2 * dsldz(i,j,k) & - + weight3 * dsldz(i,j,k-1) - dqwdzm(i,j,k) = weight2 * dqwdz(i,j,k) & - + weight3 * dqwdz(i,j,k-1) - dbdz(i,j,k) = g*( bt_gb(i,j,k-1)*dsldzm(i,j,k) + & - bq_gb(i,j,k-1)*dqwdzm(i,j,k) )/weight1 - ! ! Now with gradient adjustment - dsldzm_ga = weight2 * dsldz_ga(i,j,k) & - + weight3 * dsldz_ga(i,j,k-1) - dbdz_ga(i,j,k) = g*( bt_gb(i,j,k-1)*dsldzm_ga + & - bq_gb(i,j,k-1)*dqwdzm(i,j,k) )/weight1 - ! Complete calculation of interpolated gradients - dsldzm(i,j,k) = dsldzm(i,j,k) / weight1 - dqwdzm(i,j,k) = dqwdzm(i,j,k) / weight1 - end do - end do + do i = pdims%i_start, pdims%i_end + weight1 = r_rho_levels(i,j,k) - & + r_rho_levels(i,j,k-1) + weight2 = r_theta_levels(i,j,k-1)- & + r_rho_levels(i,j,k-1) + weight3 = r_rho_levels(i,j,k) - & + r_theta_levels(i,j,k-1) + dsldzm(i,j,k) = weight2 * dsldz(i,j,k) & + + weight3 * dsldz(i,j,k-1) + dqwdzm(i,j,k) = weight2 * dqwdz(i,j,k) & + + weight3 * dqwdz(i,j,k-1) + dbdz(i,j,k) = g*( bt_gb(i,j,k-1)*dsldzm(i,j,k) + & + bq_gb(i,j,k-1)*dqwdzm(i,j,k) )/weight1 + ! ! Now with gradient adjustment + dsldzm_ga = weight2 * dsldz_ga(i,j,k) & + + weight3 * dsldz_ga(i,j,k-1) + dbdz_ga(i,j,k) = g*( bt_gb(i,j,k-1)*dsldzm_ga + & + bq_gb(i,j,k-1)*dqwdzm(i,j,k) )/weight1 + ! Complete calculation of interpolated gradients + dsldzm(i,j,k) = dsldzm(i,j,k) / weight1 + dqwdzm(i,j,k) = dqwdzm(i,j,k) / weight1 + end do + !end do end do !$OMP end do @@ -1056,15 +1093,15 @@ subroutine bdy_expl2 ( & ! (ie instead of using centred value that uses surface parameters) if ( .not. l_use_surf_in_ri ) then k = 2 + ! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - dbdz(i,j,k) = g*( bt_gb(i,j,k-1)*dsldz(i,j,k) + & - bq_gb(i,j,k-1)*dqwdz(i,j,k) ) - dbdz_ga(i,j,k) = g*( bt_gb(i,j,k-1)*dsldz_ga(i,j,k) + & - bq_gb(i,j,k-1)*dqwdz(i,j,k) ) - end do + do i = pdims%i_start, pdims%i_end + dbdz(i,j,k) = g*( bt_gb(i,j,k-1)*dsldz(i,j,k) + & + bq_gb(i,j,k-1)*dqwdz(i,j,k) ) + dbdz_ga(i,j,k) = g*( bt_gb(i,j,k-1)*dsldz_ga(i,j,k) + & + bq_gb(i,j,k-1)*dqwdz(i,j,k) ) end do + ! end do !$OMP end do end if @@ -1084,134 +1121,147 @@ subroutine bdy_expl2 ( & ! Calculate total-water supersaturation qw - qsat(Tl), on theta-levels !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels - ! Calculate qsat(Tl)... - if ( l_mr_physics ) then - call qsat_mix( qs_tl, tl(:,:,k), p_theta_levels(:,:,k), & - tdims%i_len, tdims%j_len ) - else - call qsat( qs_tl, tl(:,:,k), p_theta_levels(:,:,k), & - tdims%i_len, tdims%j_len ) - end if + do ii = tdims%i_start, tdims%i_end, bl_segment_size + seg_slice_start = ii + seg_slice_end = MIN(((ii+bl_segment_size)-1), tdims%i_end) + bl_segment_range = (seg_slice_end - seg_slice_start) + 1 + + ! Calculate qsat(Tl)... + if ( l_mr_physics ) then + call qsat_mix(qs_tl(seg_slice_start:seg_slice_end,:,k), & + tl(seg_slice_start:seg_slice_end,:,k), & + p_theta_levels(seg_slice_start:seg_slice_end,:,k), & + bl_segment_range, & + tdims%j_len ) + else + call qsat(qs_tl(seg_slice_start:seg_slice_end,:,k), & + tl(seg_slice_start:seg_slice_end,:,k), & + p_theta_levels(seg_slice_start:seg_slice_end,:,k), & + bl_segment_range, & + tdims%j_len ) + end if + end do ! ii ! ...then subtract from qw to get supersaturation, and multiply by ! 1/(1 + Lc/cp dqsat/dT) ! in order to compare with values of qcl+qcf. - do j = tdims%j_start, tdims%j_end - do i = tdims%i_start, tdims%i_end - supersat(i,j,k) = a_qs(i,j,k) * ( qw(i,j,k) - qs_tl(i,j) ) - end do - end do - end do + !do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end + supersat(i,j,k) = a_qs(i,j,k) * ( qw(i,j,k) - qs_tl(i,j,k) ) + end do !i + !end do !j + end do !k !$OMP end do ! Calculate dbdz on rho-levels, so between th-levels k-1 and k !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = tdims%j_start, tdims%j_end - do i = tdims%i_start, tdims%i_end + !do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end - ! Interpolate the buoyancy coefficients onto rho-levels... - ! Split the interpolation into weighted contributions from: - ! a) Area that is cloudy at k-1 and at k - ! (just use saturated coefficients). - ! b) Area that is cloudy at one height but not at the other - ! (interpolate qw-qsat to find fraction of level that is saturated). - ! c) Area that is cloud-free at k-1 and k - ! (just use cloud-free coefficients). - - ! Area (a) is min of cloud-fractions from both levels (max overlap) - frac_sat = min( cf_bulk(i,j,k-1), cf_bulk(i,j,k) ) - ! Area (c) is 1 - max of both levels - frac_dry = one - max( cf_bulk(i,j,k-1), cf_bulk(i,j,k) ) - ! Area (b) (cloud-edge between k-1 and k) is the remainder - frac_edg = one - (frac_sat + frac_dry) - - ! Find volume fraction of the model-level that is saturated within - ! the "edge" area - if ( cf_bulk(i,j,k-1) > cf_bulk(i,j,k) ) then - ! There is less cloud area at k than k-1 - qc_tot = (qcl(i,j,k-1)+qcf(i,j,k-1))/cf_bulk(i,j,k-1) - if ( supersat(i,j,k-1) - supersat(i,j,k) > qc_tot ) then - frac_lev = qc_tot / ( supersat(i,j,k-1) - supersat(i,j,k) ) - else - frac_lev = one - end if - else if ( cf_bulk(i,j,k) > cf_bulk(i,j,k-1) ) then - ! There is more cloud area at k than k-1 - qc_tot = (qcl(i,j,k)+qcf(i,j,k))/cf_bulk(i,j,k) - if ( supersat(i,j,k) - supersat(i,j,k-1) > qc_tot ) then - frac_lev = qc_tot / ( supersat(i,j,k) - supersat(i,j,k-1) ) - else - frac_lev = one - end if + ! Interpolate the buoyancy coefficients onto rho-levels... + ! Split the interpolation into weighted contributions from: + ! a) Area that is cloudy at k-1 and at k + ! (just use saturated coefficients). + ! b) Area that is cloudy at one height but not at the other + ! (interpolate qw-qsat to find fraction of level that is saturated). + ! c) Area that is cloud-free at k-1 and k + ! (just use cloud-free coefficients). + + ! Area (a) is min of cloud-fractions from both levels (max overlap) + frac_sat = min( cf_bulk(i,j,k-1), cf_bulk(i,j,k) ) + ! Area (c) is 1 - max of both levels + frac_dry = one - max( cf_bulk(i,j,k-1), cf_bulk(i,j,k) ) + ! Area (b) (cloud-edge between k-1 and k) is the remainder + frac_edg = one - (frac_sat + frac_dry) + + ! Find volume fraction of the model-level that is saturated within + ! the "edge" area + if ( cf_bulk(i,j,k-1) > cf_bulk(i,j,k) ) then + ! There is less cloud area at k than k-1 + qc_tot = (qcl(i,j,k-1)+qcf(i,j,k-1))/cf_bulk(i,j,k-1) + if ( supersat(i,j,k-1) - supersat(i,j,k) > qc_tot ) then + frac_lev = qc_tot / ( supersat(i,j,k-1) - supersat(i,j,k) ) else - ! CF equal on both levels (edge contribution will be zero) - frac_lev = zero + frac_lev = one end if + else if ( cf_bulk(i,j,k) > cf_bulk(i,j,k-1) ) then + ! There is more cloud area at k than k-1 + qc_tot = (qcl(i,j,k)+qcf(i,j,k))/cf_bulk(i,j,k) + if ( supersat(i,j,k) - supersat(i,j,k-1) > qc_tot ) then + frac_lev = qc_tot / ( supersat(i,j,k) - supersat(i,j,k-1) ) + else + frac_lev = one + end if + else + ! CF equal on both levels (edge contribution will be zero) + frac_lev = zero + end if - ! Find saturated volume-fraction of the rho-level - weight2 = frac_sat + frac_lev * frac_edg - - ! Find vertical interpolation weight - weight1 = ( z_uv(i,j,k) - z_tq(i,j,k-1) ) & - / ( z_tq(i,j,k) - z_tq(i,j,k-1) ) - - ! Interpolate to find combined volume-average buoyancy coefficients - bt_rh = weight2 * ( (one-weight1) * bt_cld(i,j,k-1) & - + weight1 * bt_cld(i,j,k) ) & - + (one-weight2) * ( (one-weight1) * bt(i,j,k-1) & - + weight1 * bt(i,j,k) ) - bq_rh = weight2 * ( (one-weight1) * bq_cld(i,j,k-1) & - + weight1 * bq_cld(i,j,k) ) & - + (one-weight2) * ( (one-weight1) * bq(i,j,k-1) & - + weight1 * bq(i,j,k) ) - - ! Compute dbdz - dbdz_rh(i,j,k) = g * ( bt_rh * dsldz(i,j,k) & - + bq_rh * dqwdz(i,j,k) ) - ! Also compute gradient-adjusted version - dbdz_ga_rh(i,j,k) = g * ( bt_rh * dsldz_ga(i,j,k) & - + bq_rh * dqwdz(i,j,k) ) + ! Find saturated volume-fraction of the rho-level + weight2 = frac_sat + frac_lev * frac_edg - end do + ! Find vertical interpolation weight + weight1 = ( z_uv(i,j,k) - z_tq(i,j,k-1) ) & + / ( z_tq(i,j,k) - z_tq(i,j,k-1) ) + + ! Interpolate to find combined volume-average buoyancy coefficients + bt_rh = weight2 * ( (one-weight1) * bt_cld(i,j,k-1) & + + weight1 * bt_cld(i,j,k) ) & + + (one-weight2) * ( (one-weight1) * bt(i,j,k-1) & + + weight1 * bt(i,j,k) ) + bq_rh = weight2 * ( (one-weight1) * bq_cld(i,j,k-1) & + + weight1 * bq_cld(i,j,k) ) & + + (one-weight2) * ( (one-weight1) * bq(i,j,k-1) & + + weight1 * bq(i,j,k) ) + + ! Compute dbdz + dbdz_rh(i,j,k) = g * ( bt_rh * dsldz(i,j,k) & + + bq_rh * dqwdz(i,j,k) ) + ! Also compute gradient-adjusted version + dbdz_ga_rh(i,j,k) = g * ( bt_rh * dsldz_ga(i,j,k) & + + bq_rh * dqwdz(i,j,k) ) + + !end do end do end do !$OMP end do NOWAIT k = 1 -!$OMP do SCHEDULE(STATIC) - do j = tdims%j_start, tdims%j_end - do i = tdims%i_start, tdims%i_end - ! At surface, just use buoyancy coefficients from the first theta-level - dbdz_rh(i,j,k) = g * ( bt_gb(i,j,k) * dsldz(i,j,k) & - + bq_gb(i,j,k) * dqwdz(i,j,k) ) - dbdz_ga_rh(i,j,k) = g * ( bt_gb(i,j,k) * dsldz_ga(i,j,k) & - + bq_gb(i,j,k) * dqwdz(i,j,k) ) - end do + + !do j = tdims%j_start, tdims%j_end + !$OMP do SCHEDULE(STATIC) + do i = tdims%i_start, tdims%i_end + ! At surface, just use buoyancy coefficients from the first theta-level + dbdz_rh(i,j,k) = g * ( bt_gb(i,j,k) * dsldz(i,j,k) & + + bq_gb(i,j,k) * dqwdz(i,j,k) ) + dbdz_ga_rh(i,j,k) = g * ( bt_gb(i,j,k) * dsldz_ga(i,j,k) & + + bq_gb(i,j,k) * dqwdz(i,j,k) ) + !end do end do -!$OMP end do + !$OMP end do + !do j = tdims%j_start, tdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = tdims%j_start, tdims%j_end - do i = tdims%i_start, tdims%i_end + do i = tdims%i_start, tdims%i_end - ! Interpolate dbdz to theta-levels - ! Note: dbdz(k) is defined on theta-level k-1 - weight1 = ( z_tq(i,j,k-1) - z_uv(i,j,k-1) ) & - / ( z_uv(i,j,k) - z_uv(i,j,k-1) ) - dbdz(i,j,k) = (one-weight1) * dbdz_rh(i,j,k-1) & - + weight1 * dbdz_rh(i,j,k) - dbdz_ga(i,j,k) = (one-weight1) * dbdz_ga_rh(i,j,k-1) & - + weight1 * dbdz_ga_rh(i,j,k) - - ! Also need to interpolate dsldz, dqwdz onto theta-levels for use - ! in the RHcrit calculation - dsldzm(i,j,k) = (one-weight1) * dsldz(i,j,k-1) & - + weight1 * dsldz(i,j,k) - dqwdzm(i,j,k) = (one-weight1) * dqwdz(i,j,k-1) & - + weight1 * dqwdz(i,j,k) + ! Interpolate dbdz to theta-levels + ! Note: dbdz(k) is defined on theta-level k-1 + weight1 = ( z_tq(i,j,k-1) - z_uv(i,j,k-1) ) & + / ( z_uv(i,j,k) - z_uv(i,j,k-1) ) + dbdz(i,j,k) = (one-weight1) * dbdz_rh(i,j,k-1) & + + weight1 * dbdz_rh(i,j,k) + dbdz_ga(i,j,k) = (one-weight1) * dbdz_ga_rh(i,j,k-1) & + + weight1 * dbdz_ga_rh(i,j,k) + + ! Also need to interpolate dsldz, dqwdz onto theta-levels for use + ! in the RHcrit calculation + dsldzm(i,j,k) = (one-weight1) * dsldz(i,j,k-1) & + + weight1 * dsldz(i,j,k) + dqwdzm(i,j,k) = (one-weight1) * dqwdz(i,j,k-1) & + + weight1 * dqwdz(i,j,k) - end do end do + !end do end do !$OMP end do @@ -1224,15 +1274,15 @@ subroutine bdy_expl2 ( & !-------------------------------------------------- if (.not. l_subfilter_vert) then +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - dzu = u_p(i,j,k) - u_p(i,j,k-1) - dzv = v_p(i,j,k) - v_p(i,j,k-1) - dvdzm(i,j,k) = max( 1.0e-12_r_bl, & - sqrt(dzu*dzu + dzv*dzv) * rdz(i,j,k) ) - end do + do i = pdims%i_start, pdims%i_end + dzu = u_p(i,j,k) - u_p(i,j,k-1) + dzv = v_p(i,j,k) - v_p(i,j,k-1) + dvdzm(i,j,k) = max( 1.0e-12_r_bl, & + sqrt(dzu*dzu + dzv*dzv) * rdz(i,j,k) ) + !end do end do end do !$OMP end do @@ -1242,31 +1292,30 @@ subroutine bdy_expl2 ( & if ( model_type == mt_single_column ) then ! visc_m, visc_h need to be the 1D shear(k) on theta-level(k) +! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - dzu = u_p(i,j,k) - u_p(i,j,k-1) - dzv = v_p(i,j,k) - v_p(i,j,k-1) - dvdzm(i,j,k) = max( 1.0e-12_r_bl, & - sqrt(dzu*dzu + dzv*dzv) * rdz(i,j,k) ) - visc_m(i,j,k-1) = dvdzm(i,j,k) - visc_h(i,j,k-1) = dvdzm(i,j,k) - end do + do i = pdims%i_start, pdims%i_end + dzu = u_p(i,j,k) - u_p(i,j,k-1) + dzv = v_p(i,j,k) - v_p(i,j,k-1) + dvdzm(i,j,k) = max( 1.0e-12_r_bl, & + sqrt(dzu*dzu + dzv*dzv) * rdz(i,j,k) ) + visc_m(i,j,k-1) = dvdzm(i,j,k) + visc_h(i,j,k-1) = dvdzm(i,j,k) + !end do end do end do !$OMP end do else - ! On entry, visc_m is 3D shear(k) on theta-level(k) +! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - dvdzm(i,j,k) = max( 1.0e-12_r_bl, visc_m(i,j,k-1) ) - end do + do i = pdims%i_start, pdims%i_end + dvdzm(i,j,k) = max( 1.0e-12_r_bl, visc_m(i,j,k-1) ) + ! end do end do end do !$OMP end do @@ -1278,13 +1327,13 @@ subroutine bdy_expl2 ( & if (smag_l_calc == smag_l_calc_use_geo) then ! use 3d grid geometric mean +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - rmlmax2(i,j,k) = ( delta_smag(i,j) * delta_smag(i,j) * & - dzl_charney(i,j,k) )**one_third - end do + do i = pdims%i_start, pdims%i_end + rmlmax2(i,j,k) = ( delta_smag(i,j) * delta_smag(i,j) * & + dzl_charney(i,j,k) )**one_third + !end do end do end do !$OMP end do @@ -1293,10 +1342,10 @@ subroutine bdy_expl2 ( & ! only use horizontal grid !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - rmlmax2(i,j,k) = delta_smag(i,j) - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + rmlmax2(i,j,k) = delta_smag(i,j) + !end do end do end do !$OMP end do @@ -1305,38 +1354,38 @@ subroutine bdy_expl2 ( & if ( l_rp2 .and. i_rp_scheme == i_rp2b ) then +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - rmlmax2(i,j,k) = ( cs_rp(rp_idx) * rmlmax2(i,j,k) )**2 - end do + do i = pdims%i_start, pdims%i_end + rmlmax2(i,j,k) = ( cs_rp(rp_idx) * rmlmax2(i,j,k) )**2 + !end do end do end do !$OMP end do else +! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - rmlmax2(i,j,k) = ( mix_factor * rmlmax2(i,j,k) )**2 - end do + do i = pdims%i_start, pdims%i_end + rmlmax2(i,j,k) = ( mix_factor * rmlmax2(i,j,k) )**2 + !end do end do end do !$OMP end do end if +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - rneutml_sq(i,j,k) = one / ( & - one/( vkman*(z_tq(i,j,k) + z0m_eff_gb(i,j)) )**2 & - + one/rmlmax2(i,j,k) ) - end do + do i = pdims%i_start, pdims%i_end + rneutml_sq(i,j,k) = one / ( & + one/( vkman*(z_tq(i,j,k) + z0m_eff_gb(i,j)) )**2 & + + one/rmlmax2(i,j,k) ) + !end do end do end do !$OMP end do @@ -1347,11 +1396,11 @@ subroutine bdy_expl2 ( & !----------------------------------------------------------------------- ! Set-up 2D array for standard deviation of subgrid orography. !----------------------------------------------------------------------- +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - sigma_h(i,j) = zero - end do +do i = pdims%i_start, pdims%i_end + sigma_h(i,j) = zero +!end do end do !$OMP end do @@ -1368,77 +1417,77 @@ subroutine bdy_expl2 ( & if (sg_orog_mixing == sg_shear .or. & sg_orog_mixing == sg_shear_enh_lambda) then +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end + do i = pdims%i_start, pdims%i_end - if (sigma_h(i,j) > one ) then - zpr = z_tq(i,j,k-1)/sigma_h(i,j) - ! Height dependence, to reduce effect to zero with height - ! gives z_scale~[1,0.95,0.5,0] at zpr=[0,0.6,1,1.7] - weight1 = one_half*( one - tanh(4.0_r_bl*(zpr-one) ) ) + if (sigma_h(i,j) > one ) then + zpr = z_tq(i,j,k-1)/sigma_h(i,j) + ! Height dependence, to reduce effect to zero with height + ! gives z_scale~[1,0.95,0.5,0] at zpr=[0,0.6,1,1.7] + weight1 = one_half*( one - tanh(4.0_r_bl*(zpr-one) ) ) - ! Take slope ~ sd/h_scale for small sd; - ! tends to 0.2 for large sd - slope = one / sqrt( 25.0_r_bl + (h_scale/sigma_h(i,j))**2 ) + ! Take slope ~ sd/h_scale for small sd; + ! tends to 0.2 for large sd + slope = one / sqrt( 25.0_r_bl + (h_scale/sigma_h(i,j))**2 ) - dvdzm(i,j,k) = max ( dvdzm(i,j,k), & - weight1*slope*t_drain*dbdz(i,j,k) ) + dvdzm(i,j,k) = max ( dvdzm(i,j,k), & + weight1*slope*t_drain*dbdz(i,j,k) ) - if (k==2 .and. BL_diag%l_dvdzm) & - BL_diag%dvdzm(i,j,1)=weight1*slope*t_drain*dbdz(i,j,k) + if (k==2 .and. BL_diag%l_dvdzm) & + BL_diag%dvdzm(i,j,1)=weight1*slope*t_drain*dbdz(i,j,k) - end if - end do + end if + !end do end do end do !$OMP end do end if ! sg_orog_mixing +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ri(i,j,k) = dbdz(i,j,k) / ( dvdzm(i,j,k)*dvdzm(i,j,k) ) - ri(i,j,k) = max(min(ri(i,j,k),max_ri),-max_ri) - ri_ga(i,j,k) = dbdz_ga(i,j,k) / ( dvdzm(i,j,k)*dvdzm(i,j,k) ) - end do + do i = pdims%i_start, pdims%i_end + ri(i,j,k) = dbdz(i,j,k) / ( dvdzm(i,j,k)*dvdzm(i,j,k) ) + ri(i,j,k) = max(min(ri(i,j,k),max_ri),-max_ri) + ri_ga(i,j,k) = dbdz_ga(i,j,k) / ( dvdzm(i,j,k)*dvdzm(i,j,k) ) + ! end do end do end do !$OMP end do +!do j = pdims%j_start, pdims%j_end if (BL_diag%l_gradrich) then !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - BL_diag%gradrich(i,j,k)=ri(i,j,k) - end do + do i = pdims%i_start, pdims%i_end + BL_diag%gradrich(i,j,k)=ri(i,j,k) + !end do end do end do !$OMP end do end if +!do j = pdims%j_start, pdims%j_end if (BL_diag%l_dbdz) then !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - BL_diag%dbdz(i,j,k)=dbdz(i,j,k) - end do + do i = pdims%i_start, pdims%i_end + BL_diag%dbdz(i,j,k)=dbdz(i,j,k) + !end do end do end do !$OMP end do end if +!do j = pdims%j_start, pdims%j_end if (BL_diag%l_dvdzm) then !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - BL_diag%dvdzm(i,j,k)=dvdzm(i,j,k) - end do + do i = pdims%i_start, pdims%i_end + BL_diag%dvdzm(i,j,k)=dvdzm(i,j,k) + !end do end do end do !$OMP end do @@ -1467,32 +1516,32 @@ subroutine bdy_expl2 ( & ! Orographic stress diagnostics !------------------------------------------------------------------ if (BL_diag%l_ostressx) then +!do j = tdims%j_start, tdims%j_end !$OMP PARALLEL do & !$OMP SCHEDULE(STATIC) & !$OMP DEFAULT(none) & -!$OMP private(k,j,i) & -!$OMP SHARED(bl_levels,tdims,BL_diag,tau_fd_x) +!$OMP private(k,i) & +!$OMP SHARED(j,bl_levels,tdims,BL_diag,tau_fd_x) do k = 1, bl_levels - do j = tdims%j_start, tdims%j_end - do i = tdims%i_start, tdims%i_end - BL_diag%ostressx(i,j,k)=tau_fd_x(i,j,k) - end do + do i = tdims%i_start, tdims%i_end + BL_diag%ostressx(i,j,k)=tau_fd_x(i,j,k) + !end do end do end do !$OMP end PARALLEL do end if + ! do j = tdims%j_start, tdims%j_end if (BL_diag%l_ostressy) then !$OMP PARALLEL do & !$OMP SCHEDULE(STATIC) & !$OMP DEFAULT(none) & -!$OMP private(k,j,i) & -!$OMP SHARED(bl_levels,tdims,BL_diag,tau_fd_y) +!$OMP private(k,i) & +!$OMP SHARED(j, bl_levels,tdims,BL_diag,tau_fd_y) do k = 1, bl_levels - do j = tdims%j_start, tdims%j_end - do i = tdims%i_start, tdims%i_end - BL_diag%ostressy(i,j,k)=tau_fd_y(i,j,k) - end do + do i = tdims%i_start, tdims%i_end + BL_diag%ostressy(i,j,k)=tau_fd_y(i,j,k) + !end do end do end do !$OMP end PARALLEL do @@ -1511,25 +1560,21 @@ subroutine bdy_expl2 ( & ! of importance mainly at sea-points, to avoid complications ! with coastal tiling, the scheme operates only at points ! where the land fraction is below 0.5. - -omp_block = pdims%j_end -!$ omp_block = ceiling(real(pdims%j_end)/omp_get_max_threads()) - -!$OMP PARALLEL DEFAULT(SHARED) private(i, j, k, jj, ntop, z_scale) - -!$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ! First override the provisional cumulus diagnosis if the - ! actual surface buoyancy flux is stable - if ( fb_surf(i,j) < zero ) then - cumulus(i,j) = .false. - l_shallow(i,j) = .false. - ntml(i,j) = 1 - zh(i,j) = z_uv(i,j,2) - dzh(i,j) = rmdi - end if - end do +!do j = pdims%j_start, pdims%j_end +!$OMP PARALLEL DEFAULT(SHARED) private(i, k, ii, ntop, z_scale) + +!$OMP do SCHEDULE(STATIC) +do i = pdims%i_start, pdims%i_end + ! First override the provisional cumulus diagnosis if the + ! actual surface buoyancy flux is stable + if ( fb_surf(i,j) < zero ) then + cumulus(i,j) = .false. + l_shallow(i,j) = .false. + ntml(i,j) = 1 + zh(i,j) = z_uv(i,j,2) + dzh(i,j) = rmdi + end if +!end do end do !$OMP end do @@ -1539,82 +1584,79 @@ subroutine bdy_expl2 ( & ! Original version - causes spuriously deep boundary layers if ! BL_LEVELS is >> 3km +! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end + do i = pdims%i_start, pdims%i_end - if ( flandg(i,j) < one_half ) then - ntop = min(ntpar(i,j),bl_levels-1) - if ( -z_uv(i,j,ntop+1) * recip_l_mo_sea(i,j) < near_neut_z_on_l ) then - cumulus(i,j) = .false. - l_shallow(i,j) = .false. - ntml(i,j) = ntop - zh(i,j) = z_uv(i,j,ntml(i,j)+1) - end if + if ( flandg(i,j) < one_half ) then + ntop = min(ntpar(i,j),bl_levels-1) + if ( -z_uv(i,j,ntop+1) * recip_l_mo_sea(i,j) < near_neut_z_on_l ) then + cumulus(i,j) = .false. + l_shallow(i,j) = .false. + ntml(i,j) = ntop + zh(i,j) = z_uv(i,j,ntml(i,j)+1) end if + end if - end do + !end do end do !$OMP end do else if (idyndiag == DynDiag_ZL_corrn) then + ! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end + do i = pdims%i_start, pdims%i_end - !------------------------------------------------------------ - ! As original code above, except restrict depth scale to <3km - ! and, if near-neutral, ignore the BL depth diagnosed by the - ! adiabatic parcel (ie ZH, NTML) completely. - !------------------------------------------------------------ - if ( flandg(i,j) < one_half ) then - ntop = min(ntpar(i,j),bl_levels-1) - z_scale = min( 3000.0_r_bl, z_uv(i,j,ntop+1) ) - if ( -z_scale*recip_l_mo_sea(i,j) < near_neut_z_on_l ) then - cumulus(i,j) = .false. - l_shallow(i,j) = .false. - ntml(i,j) = 1 - zh(i,j) = z_uv(i,j,ntml(i,j)+1) - end if + !------------------------------------------------------------ + ! As original code above, except restrict depth scale to <3km + ! and, if near-neutral, ignore the BL depth diagnosed by the + ! adiabatic parcel (ie ZH, NTML) completely. + !------------------------------------------------------------ + if ( flandg(i,j) < one_half ) then + ntop = min(ntpar(i,j),bl_levels-1) + z_scale = min( 3000.0_r_bl, z_uv(i,j,ntop+1) ) + if ( -z_scale*recip_l_mo_sea(i,j) < near_neut_z_on_l ) then + cumulus(i,j) = .false. + l_shallow(i,j) = .false. + ntml(i,j) = 1 + zh(i,j) = z_uv(i,j,ntml(i,j)+1) end if - - end do + end if + !end do end do !$OMP end do else if (idyndiag == DynDiag_ZL_CuOnly) then - + ! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end + do i = pdims%i_start, pdims%i_end - !------------------------------------------------------------ - ! As DynDiag_ZL_corrn but only affects cumulus points and - ! points that are entirely sea. Note that DynDiag_ZL_corrn - ! has been found to switch off non-local mixing in - ! stratocumulus, where surface fluxes are typically small - !------------------------------------------------------------ - if ( cumulus(i,j) .and. flandg(i,j) < 0.01_r_bl ) then - ntop = min(ntpar(i,j),bl_levels-1) - z_scale = min( 3000.0_r_bl, z_uv(i,j,ntop+1) ) - if ( -z_scale*recip_l_mo_sea(i,j) < near_neut_z_on_l ) then - ! - ZH/L indicates BL close to neutral - dynamic_bl_diag(i,j) = .true. - cumulus(i,j) = .false. - l_shallow(i,j) = .false. - if (l_fix_dyndiag) then - ntml(i,j) = 1 - else - ! note that ntop can be very high making a spuriously - ! deep well-mixed PBL - hence depricated - ntml(i,j) = ntop - end if - zh(i,j) = z_uv(i,j,ntml(i,j)+1) + !------------------------------------------------------------ + ! As DynDiag_ZL_corrn but only affects cumulus points and + ! points that are entirely sea. Note that DynDiag_ZL_corrn + ! has been found to switch off non-local mixing in + ! stratocumulus, where surface fluxes are typically small + !------------------------------------------------------------ + if ( cumulus(i,j) .and. flandg(i,j) < 0.01_r_bl ) then + ntop = min(ntpar(i,j),bl_levels-1) + z_scale = min( 3000.0_r_bl, z_uv(i,j,ntop+1) ) + if ( -z_scale*recip_l_mo_sea(i,j) < near_neut_z_on_l ) then + ! - ZH/L indicates BL close to neutral + dynamic_bl_diag(i,j) = .true. + cumulus(i,j) = .false. + l_shallow(i,j) = .false. + if (l_fix_dyndiag) then + ntml(i,j) = 1 + else + ! note that ntop can be very high making a spuriously + ! deep well-mixed PBL - hence depricated + ntml(i,j) = ntop end if + zh(i,j) = z_uv(i,j,ntml(i,j)+1) end if - - end do + end if + !end do end do !$OMP end do @@ -1623,29 +1665,29 @@ subroutine bdy_expl2 ( & ! As DynDiag_ZL_CuOnly but also allow ZH(Ri) to overrule the ! Cumulus diagnosis !------------------------------------------------------------ + !do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - topbl(i,j) = .false. - end do + do i = pdims%i_start, pdims%i_end + topbl(i,j) = .false. + ! end do end do !$OMP end do !--------------------------------------------------------------- ! Loop over levels to find Ri > RiCrit_sharp (=0.25) to find ! level to which Ri really is close to neutral !--------------------------------------------------------------- - +! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) - do jj = pdims%j_start, pdims%j_end, omp_block + do ii = pdims%i_start, pdims%i_end, pdims_omp_block do k = 2, bl_levels - do j = jj, min(jj+omp_block-1,pdims%j_end) - do i = pdims%i_start, pdims%i_end - if ( .not. topbl(i,j) .and. & - (ri_ga(i,j,k) > RiCrit_sharp .or. k > bl_levels-1) ) then - topbl(i,j) = .true. - zh_local(i,j) = z_uv(i,j,k) - end if - end do ! Loop over points + + do i = ii, min(((ii+pdims_omp_block)-1),pdims%i_end) + if ( .not. topbl(i,j) .and. & + (ri_ga(i,j,k) > RiCrit_sharp .or. k > bl_levels-1) ) then + topbl(i,j) = .true. + zh_local(i,j) = z_uv(i,j,k) + end if + ! end do ! Loop over points end do ! Loop over points end do ! Loop over levels end do @@ -1661,63 +1703,63 @@ subroutine bdy_expl2 ( & zhloc_depth_fac = zhloc_depth_fac_rp(rp_idx) end if +! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end + do i = pdims%i_start, pdims%i_end - if ( cumulus(i,j) .and. flandg(i,j) < 0.01_r_bl ) then - ntop = min(ntpar(i,j),bl_levels-1) - z_scale = min( 3000.0_r_bl, z_uv(i,j,ntop+1) ) - if ( zh_local(i,j) & - > zh(i,j)+zhloc_depth_fac*(zhpar(i,j)-zh(i,j)) ) then - ! ZH(Ri>RiCrit) more than zhloc_depth_fac up the - ! cloud layer, indicating significant shear disruption - dynamic_bl_diag(i,j) = .true. - cumulus(i,j) = .false. - l_shallow(i,j) = .false. - ntml(i,j) = ntop - zh(i,j) = z_uv(i,j,ntml(i,j)+1) - else if ( -z_scale*recip_l_mo_sea(i,j) < near_neut_z_on_l ) then - ! - ZH/L indicates BL close to neutral so overrule - ! cumulus diagnosis and leave mixing to the local - ! Ri-based scheme, given no better information on where - ! unstable bl top should be - dynamic_bl_diag(i,j) = .true. - cumulus(i,j) = .false. - l_shallow(i,j) = .false. - ntml(i,j) = 1 - zh(i,j) = z_uv(i,j,ntml(i,j)+1) - end if - end if ! Cu over sea + if ( cumulus(i,j) .and. flandg(i,j) < 0.01_r_bl ) then + ntop = min(ntpar(i,j),bl_levels-1) + z_scale = min( 3000.0_r_bl, z_uv(i,j,ntop+1) ) + if ( zh_local(i,j) & + > zh(i,j)+zhloc_depth_fac*(zhpar(i,j)-zh(i,j)) ) then + ! ZH(Ri>RiCrit) more than zhloc_depth_fac up the + ! cloud layer, indicating significant shear disruption + dynamic_bl_diag(i,j) = .true. + cumulus(i,j) = .false. + l_shallow(i,j) = .false. + ntml(i,j) = ntop + zh(i,j) = z_uv(i,j,ntml(i,j)+1) + else if ( -z_scale*recip_l_mo_sea(i,j) < near_neut_z_on_l ) then + ! - ZH/L indicates BL close to neutral so overrule + ! cumulus diagnosis and leave mixing to the local + ! Ri-based scheme, given no better information on where + ! unstable bl top should be + dynamic_bl_diag(i,j) = .true. + cumulus(i,j) = .false. + l_shallow(i,j) = .false. + ntml(i,j) = 1 + zh(i,j) = z_uv(i,j,ntml(i,j)+1) + end if + end if ! Cu over sea - end do + !end do end do !$OMP end do else ! old version: not good to set ntml=ntop for small zh/L +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end + do i = pdims%i_start, pdims%i_end - if ( cumulus(i,j) .and. flandg(i,j) < 0.01_r_bl ) then - ntop = min(ntpar(i,j),bl_levels-1) - z_scale = min( 3000.0_r_bl, z_uv(i,j,ntop+1) ) - if ( -z_scale*recip_l_mo_sea(i,j) < near_neut_z_on_l .or. & - ! - ZH/L indicates BL close to neutral - zh_local(i,j) > zh(i,j)+zhloc_depth_fac*(z_scale-zh(i,j)) & - ! ZH(Ri>RiCrit) more than zhloc_depth_fac up the - ! cloud layer, indicating significant shear disruption - ) then - dynamic_bl_diag(i,j) = .true. - cumulus(i,j) = .false. - l_shallow(i,j) = .false. - ntml(i,j) = ntop - zh(i,j) = z_uv(i,j,ntml(i,j)+1) - end if - end if ! Cu over sea + if ( cumulus(i,j) .and. flandg(i,j) < 0.01_r_bl ) then + ntop = min(ntpar(i,j),bl_levels-1) + z_scale = min( 3000.0_r_bl, z_uv(i,j,ntop+1) ) + if ( -z_scale*recip_l_mo_sea(i,j) < near_neut_z_on_l .or. & + ! - ZH/L indicates BL close to neutral + zh_local(i,j) > zh(i,j)+zhloc_depth_fac*(z_scale-zh(i,j)) & + ! ZH(Ri>RiCrit) more than zhloc_depth_fac up the + ! cloud layer, indicating significant shear disruption + ) then + dynamic_bl_diag(i,j) = .true. + cumulus(i,j) = .false. + l_shallow(i,j) = .false. + ntml(i,j) = ntop + zh(i,j) = z_uv(i,j,ntml(i,j)+1) + end if + end if ! Cu over sea - end do + !end do end do !$OMP end do @@ -1738,16 +1780,16 @@ subroutine bdy_expl2 ( & ! Set NTML_NL to NTML as passed in from initial diagnosis routine !----------------------------------------------------------------------- +!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL do & !$OMP SCHEDULE(STATIC) & !$OMP DEFAULT(none) & -!$OMP private(j,i) & -!$OMP SHARED(pdims,ntml_nl,ntml) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ntml_nl(i,j) = ntml(i,j) - end do +!$OMP private(i) & +!$OMP SHARED(j,pdims,ntml_nl,ntml) +do i = pdims%i_start, pdims%i_end + ntml_nl(i,j) = ntml(i,j) end do +!end do !$OMP end PARALLEL do if (nl_bl_levels < bl_levels) then @@ -1755,18 +1797,18 @@ subroutine bdy_expl2 ( & zmaxb_for_dsc = 1.0e10_r_bl zmaxt_for_dsc = zmaxb_for_dsc +! do j = pdims%j_start, pdims%j_end !$OMP PARALLEL do & !$OMP SCHEDULE(STATIC) & !$OMP DEFAULT(none) & -!$OMP private(j,i) & -!$OMP SHARED(pdims,ntml_nl,nl_bl_levels,zh,z_uv) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( ntml_nl(i,j) > nl_bl_levels-1 ) then - ntml_nl(i,j) = nl_bl_levels-1 - zh(i,j) = z_uv(i,j,ntml_nl(i,j)+1) - end if - end do +!$OMP private(i) & +!$OMP SHARED(j,pdims,ntml_nl,nl_bl_levels,zh,z_uv) + do i = pdims%i_start, pdims%i_end + if ( ntml_nl(i,j) > nl_bl_levels-1 ) then + ntml_nl(i,j) = nl_bl_levels-1 + zh(i,j) = z_uv(i,j,ntml_nl(i,j)+1) + end if + !end do end do !$OMP end PARALLEL do else @@ -1777,49 +1819,49 @@ subroutine bdy_expl2 ( & ! Initialise non-local K and fluxes to zero; necessary for levels ! above NL_BL_LEVELS !----------------------------------------------------------------------- + ! do j = pdims%j_start, pdims%j_end !$OMP PARALLEL DEFAULT(none) & -!$OMP SHARED( bl_levels, pdims, ftl, fqw, rhokmz, rhokhz, rhokm_top, & +!$OMP SHARED( j, bl_levels, pdims, ftl, fqw, rhokmz, rhokhz, rhokm_top, & !$OMP rhokh_top, tke_nl, f_ngstress, rhof2, rhofsc, ft_nt, fq_nt, & !$OMP tothf_zh, tothf_zhsc, totqf_zh, totqf_zhsc, ft_nt_dscb, & -!$OMP zhnl, zh, fq_nt_dscb, tke_loc ) private( i, j, k ) +!$OMP zhnl, zh, fq_nt_dscb, tke_loc ) private( i, k ) !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ftl(i,j,k) = zero - fqw(i,j,k) = zero - rhokmz(i,j,k) = zero - rhokhz(i,j,k) = zero - rhokm_top(i,j,k) = zero - rhokh_top(i,j,k) = zero - tke_nl(i,j,k) = zero - tke_loc(i,j,k) = zero - f_ngstress(i,j,k) = zero - rhof2(i,j,k) = zero - rhofsc(i,j,k) = zero - ! Initialise Lock-Whelan non-gradient terms to zero - ! Only calculated in KMKHZ9C - ft_nt(i,j,k) = zero - fq_nt(i,j,k) = zero - end do + do i = pdims%i_start, pdims%i_end + ftl(i,j,k) = zero + fqw(i,j,k) = zero + rhokmz(i,j,k) = zero + rhokhz(i,j,k) = zero + rhokm_top(i,j,k) = zero + rhokh_top(i,j,k) = zero + tke_nl(i,j,k) = zero + tke_loc(i,j,k) = zero + f_ngstress(i,j,k) = zero + rhof2(i,j,k) = zero + rhofsc(i,j,k) = zero + ! Initialise Lock-Whelan non-gradient terms to zero + ! Only calculated in KMKHZ9C + ft_nt(i,j,k) = zero + fq_nt(i,j,k) = zero + !end do end do end do !$OMP end do NOWAIT -!$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ft_nt(i,j,1) = zero - fq_nt(i,j,1) = zero - tothf_zh(i,j) = zero - tothf_zhsc(i,j) = zero - totqf_zh(i,j) = zero - totqf_zhsc(i,j) = zero - ft_nt_dscb(i,j) = zero - fq_nt_dscb(i,j) = zero - zhnl(i,j) = zh(i,j) ! initialise non-local PBL depth - ! to that diagnosed in conv_diag - end do +!do j = pdims%j_start, pdims%j_end +!$OMP do SCHEDULE(STATIC) +do i = pdims%i_start, pdims%i_end + ft_nt(i,j,1) = zero + fq_nt(i,j,1) = zero + tothf_zh(i,j) = zero + tothf_zhsc(i,j) = zero + totqf_zh(i,j) = zero + totqf_zhsc(i,j) = zero + ft_nt_dscb(i,j) = zero + fq_nt_dscb(i,j) = zero + zhnl(i,j) = zh(i,j) ! initialise non-local PBL depth + ! to that diagnosed in conv_diag +!end do end do !$OMP end do !$OMP end PARALLEL @@ -1828,27 +1870,27 @@ subroutine bdy_expl2 ( & if (l_wtrac) then do i_wt = 1, n_wtrac +!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL DEFAULT(none) & -!$OMP SHARED( i_wt, bl_levels, pdims, wtrac_bl) private( i, j, k) +!$OMP SHARED( j, i_wt, bl_levels, pdims, wtrac_bl) private( i, k) !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - wtrac_bl(i_wt)%fqw(i,j,k) = zero - wtrac_bl(i_wt)%fq_nt(i,j,k) = zero - end do + do i = pdims%i_start, pdims%i_end + wtrac_bl(i_wt)%fqw(i,j,k) = zero + wtrac_bl(i_wt)%fq_nt(i,j,k) = zero + ! end do end do end do !$OMP end do +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - wtrac_bl(i_wt)%fq_nt(i,j,1) = zero - wtrac_bl(i_wt)%fq_nt_dscb(i,j) = zero - wtrac_bl(i_wt)%totqf_zh(i,j) = zero - wtrac_bl(i_wt)%totqf_zhsc(i,j) = zero - end do + do i = pdims%i_start, pdims%i_end + wtrac_bl(i_wt)%fq_nt(i,j,1) = zero + wtrac_bl(i_wt)%fq_nt_dscb(i,j) = zero + wtrac_bl(i_wt)%totqf_zh(i,j) = zero + wtrac_bl(i_wt)%totqf_zhsc(i,j) = zero + !end do end do !$OMP end do !$OMP end PARALLEL @@ -1893,105 +1935,106 @@ subroutine bdy_expl2 ( & ! Set all variables from the non-local scheme to zero or "off" ! - reset all fluxes and K's arising from the non-local scheme !------------------------------------------------------------- - -!$OMP PARALLEL DEFAULT(none) private(j,i,ient) & -!$OMP SHARED(pdims,unstable,fb_surf,cumulus,l_shallow,sml_disc_inv,ntpar, & +!do j = pdims%j_start, pdims%j_end +!$OMP PARALLEL DEFAULT(none) private(i,ient) & +!$OMP SHARED(j,pdims,unstable,fb_surf,cumulus,l_shallow,sml_disc_inv,ntpar, & !$OMP ntml_nl,zhnl,grad_t_adj,grad_q_adj,dsc,dsc_disc_inv,ntdsc,nbdsc, & !$OMP zhsc,dzh,coupled,kent,kent_dsc,t_frac,zrzi,we_lim,t_frac_dsc, & !$OMP zrzi_dsc,we_lim_dsc,kplume) !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ! surface mixed layer - unstable(i,j) = (fb_surf(i,j) > zero) - cumulus(i,j) = .false. - l_shallow(i,j) = .false. - sml_disc_inv(i,j) = 0 - ntpar(i,j) = 0 - ntml_nl(i,j) = -1 ! to ensure correct diagnostics - zhnl(i,j) = zero - dzh(i,j) = zero - grad_t_adj(i,j) = zero - grad_q_adj(i,j) = zero - ! decoupled mixed layer - dsc(i,j) = .false. - dsc_disc_inv(i,j) = 0 - ntdsc(i,j) = 0 - nbdsc(i,j) = 0 - zhsc(i,j) = zero - coupled(i,j) = .false. - ! entrainment variables for non-local tracer mixing - kent(i,j) = 2 - kent_dsc(i,j) = 2 - ! kplume for bi-modal cloud scheme - kplume(i,j) = 1 - end do + do i = pdims%i_start, pdims%i_end + ! surface mixed layer + unstable(i,j) = (fb_surf(i,j) > zero) + cumulus(i,j) = .false. + l_shallow(i,j) = .false. + sml_disc_inv(i,j) = 0 + ntpar(i,j) = 0 + ntml_nl(i,j) = -1 ! to ensure correct diagnostics + zhnl(i,j) = zero + dzh(i,j) = zero + grad_t_adj(i,j) = zero + grad_q_adj(i,j) = zero + ! decoupled mixed layer + dsc(i,j) = .false. + dsc_disc_inv(i,j) = 0 + ntdsc(i,j) = 0 + nbdsc(i,j) = 0 + zhsc(i,j) = zero + coupled(i,j) = .false. + ! entrainment variables for non-local tracer mixing + kent(i,j) = 2 + kent_dsc(i,j) = 2 + ! kplume for bi-modal cloud scheme + kplume(i,j) = 1 end do + ! end do !$OMP end do NOWAIT do ient = 1, 3 +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - t_frac(i,j,ient) = zero - zrzi(i,j,ient) = zero - we_lim(i,j,ient) = zero - t_frac_dsc(i,j,ient) = zero - zrzi_dsc(i,j,ient) = zero - we_lim_dsc(i,j,ient) = zero - end do + + do i = pdims%i_start, pdims%i_end + t_frac(i,j,ient) = zero + zrzi(i,j,ient) = zero + we_lim(i,j,ient) = zero + t_frac_dsc(i,j,ient) = zero + zrzi_dsc(i,j,ient) = zero + we_lim_dsc(i,j,ient) = zero + ! end do end do !$OMP end do NOWAIT end do !$OMP end PARALLEL end if ! test on NON_LOCAL_BL +!do j = pdims%j_start, pdims%j_end ! default values !$OMP PARALLEL do & !$OMP SCHEDULE(STATIC) & !$OMP DEFAULT(none) & -!$OMP private(j,i) & -!$OMP SHARED(pdims,l_shallow_cth) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - l_shallow_cth(i,j) = .false. - end do +!$OMP private(i) & +!$OMP SHARED(j,pdims,l_shallow_cth) +do i = pdims%i_start, pdims%i_end + l_shallow_cth(i,j) = .false. end do +! end do !$OMP end PARALLEL do if (blending_option == blend_cth_shcu_only) then +! do j = pdims%j_start, pdims%j_end ! only going to use the parcel top as the length scale in blending ! if the convection is shallow, where we define shallow convection here ! as the resolved cloud top (from cf_bluk) being below shallow_cu_maxtop !$OMP PARALLEL & !$OMP DEFAULT(none) & -!$OMP private(i,j,k) & -!$OMP SHARED(pdims,bl_levels,z_uv,l_shallow_cth,cumulus,ntml,cf_bulk, & +!$OMP private(i,k) & +!$OMP SHARED(j,pdims,bl_levels,z_uv,l_shallow_cth,cumulus,ntml,cf_bulk, & !$OMP shallow_cu_maxtop,cloud_base_found,cloud_top_found) !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - cloud_base_found(i,j) = .false. - cloud_top_found(i,j) = .false. - end do + do i = pdims%i_start, pdims%i_end + cloud_base_found(i,j) = .false. + cloud_top_found(i,j) = .false. + !end do end do !$OMP end do + +!do j = pdims%j_start, pdims%j_end do k = 3, bl_levels-1 !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( cumulus(i,j) .and. k > ntml(i,j)+1 .and. & - .not. cloud_base_found(i,j) .and. cf_bulk(i,j,k) > sc_cftol ) then - ! found a cumulus cloud base - cloud_base_found(i,j) = .true. - end if - if ( cloud_base_found(i,j) .and. .not. cloud_top_found(i,j) .and. & - ! found cloud-base but not yet reached cloud-top - cf_bulk(i,j,k+1) < sc_cftol ) then - ! got to cloud-top - cloud_top_found(i,j) = .true. - l_shallow_cth(i,j) = z_uv(i,j,k+1) < shallow_cu_maxtop - end if - end do + do i = pdims%i_start, pdims%i_end + if ( cumulus(i,j) .and. k > ntml(i,j)+1 .and. & + .not. cloud_base_found(i,j) .and. cf_bulk(i,j,k) > sc_cftol ) then + ! found a cumulus cloud base + cloud_base_found(i,j) = .true. + end if + if ( cloud_base_found(i,j) .and. .not. cloud_top_found(i,j) .and. & + ! found cloud-base but not yet reached cloud-top + cf_bulk(i,j,k+1) < sc_cftol ) then + ! got to cloud-top + cloud_top_found(i,j) = .true. + l_shallow_cth(i,j) = z_uv(i,j,k+1) < shallow_cu_maxtop + end if + !end do end do !$OMP end do end do @@ -2019,144 +2062,144 @@ subroutine bdy_expl2 ( & !----------------------------------------------------------------------- ! interpolate rhokh_th to rho levels 2 to bl_levels +! do j = pdims%j_start, pdims%j_end !$OMP PARALLEL DEFAULT(SHARED) & -!$OMP private(i, j, k, weight1, weight2, weight3, lambdah, z_scale, & +!$OMP private(i, k, weight1, weight2, weight3, lambdah, z_scale, & !$OMP vkz, f_log) !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - weight1 = r_theta_levels(i,j,k) - & - r_theta_levels(i,j, k-1) - weight2 = r_theta_levels(i,j,k) - & - r_rho_levels(i,j,k) - weight3 = r_rho_levels(i,j,k) - & - r_theta_levels(i,j,k-1) - if ( k == bl_levels ) then - ! assume rhokh_th(BL_LEVELS+1) is zero - rhokh(i,j,k) = ( weight2/weight1 ) * rhokh_th(i,j,k) - if (blending_option /= off) weight_1dbl_rho(i,j,k) = & - (weight2/weight1) * weight_1dbl(i,j,k) - else - rhokh(i,j,k) = (weight3/weight1) * rhokh_th(i,j,k+1) & - + (weight2/weight1) * rhokh_th(i,j,k) - if (blending_option /= off) weight_1dbl_rho(i,j,k) = & - (weight3/weight1)*weight_1dbl(i,j,k+1) & - + (weight2/weight1)*weight_1dbl(i,j,k) - end if - - if (local_fa == free_trop_layers) then - ! elh already included in rhokh_th so no need to calculate - ! here, but interpolate elh separately for diagnostic - if (BL_diag%l_elh3d) then - if ( k == bl_levels ) then + do i = pdims%i_start, pdims%i_end + weight1 = r_theta_levels(i,j,k) - & + r_theta_levels(i,j, k-1) + weight2 = r_theta_levels(i,j,k) - & + r_rho_levels(i,j,k) + weight3 = r_rho_levels(i,j,k) - & + r_theta_levels(i,j,k-1) + if ( k == bl_levels ) then ! assume rhokh_th(BL_LEVELS+1) is zero - elh_rho(i,j,k) = ( weight2/weight1 ) * elh(i,j,k) + rhokh(i,j,k) = ( weight2/weight1 ) * rhokh_th(i,j,k) + if (blending_option /= off) weight_1dbl_rho(i,j,k) = & + (weight2/weight1) * weight_1dbl(i,j,k) + else + rhokh(i,j,k) = (weight3/weight1) * rhokh_th(i,j,k+1) & + + (weight2/weight1) * rhokh_th(i,j,k) + if (blending_option /= off) weight_1dbl_rho(i,j,k) = & + (weight3/weight1)*weight_1dbl(i,j,k+1) & + + (weight2/weight1)*weight_1dbl(i,j,k) + end if + + if (local_fa == free_trop_layers) then + ! elh already included in rhokh_th so no need to calculate + ! here, but interpolate elh separately for diagnostic + if (BL_diag%l_elh3d) then + if ( k == bl_levels ) then + ! assume rhokh_th(BL_LEVELS+1) is zero + elh_rho(i,j,k) = ( weight2/weight1 ) * elh(i,j,k) + else + elh_rho(i,j,k) = & + weight3/weight1 * & + elh(i,j,k+1) & + +weight2/weight1 * & + elh(i,j,k) + end if + BL_diag%elh3d(i,j,k)=elh_rho(i,j,k) + end if + else + if ((sbl_op/=equilibrium_sbl) .or. (fb_surf(i,j) > zero)) then + !------------------------------------------------------------ + ! Include mixing length, ELH, in RHOKH. + ! Code moved from EX_COEF to avoid interpolation + !------------------------------------------------------------ + if ( k >= ntml_local(i,j)+2 .and. l_full_lambdas .and. & + local_fa == to_sharp_across_1km ) then + ! Assuming only LOCAL_FA = "to_sharp_across_1km" option + ! will have L_FULL_LAMBDAS. + ! If other LOCAL_FA options are coded here then + ! changes must be included in section 2.1 of ex_coef + if (l_rp2) then + lambdah = max ( lambda_min , & + par_mezcla_rp(rp_idx)*zh_local(i,j) ) else - elh_rho(i,j,k) = & - weight3/weight1 * & - elh(i,j,k+1) & - +weight2/weight1 * & - elh(i,j,k) + lambdah = max ( lambda_min , 0.15_r_bl*zh_local(i,j) ) end if - BL_diag%elh3d(i,j,k)=elh_rho(i,j,k) + z_scale = 1000.0_r_bl + weight1 = one_half*( one - & + tanh(3.0_r_bl*((z_uv(i,j,k)/z_scale )-one) ) ) + lambdah = lambdah * weight1 & + + lambda_min*( one - weight1) + ! no need to do log profile correction as klog_layr eq 2 + vkz = vkman * ( z_uv(i,j,k) + z0m_eff_gb(i,j) ) + elh_rho(i,j,k) = vkz / (one + vkz/lambdah ) end if - else - if ((sbl_op/=equilibrium_sbl) .or. (fb_surf(i,j) > zero)) then - !------------------------------------------------------------ - ! Include mixing length, ELH, in RHOKH. - ! Code moved from EX_COEF to avoid interpolation - !------------------------------------------------------------ - if ( k >= ntml_local(i,j)+2 .and. l_full_lambdas .and. & - local_fa == to_sharp_across_1km ) then - ! Assuming only LOCAL_FA = "to_sharp_across_1km" option - ! will have L_FULL_LAMBDAS. - ! If other LOCAL_FA options are coded here then - ! changes must be included in section 2.1 of ex_coef - if (l_rp2) then - lambdah = max ( lambda_min , & - par_mezcla_rp(rp_idx)*zh_local(i,j) ) - else - lambdah = max ( lambda_min , 0.15_r_bl*zh_local(i,j) ) - end if - z_scale = 1000.0_r_bl - weight1 = one_half*( one - & - tanh(3.0_r_bl*((z_uv(i,j,k)/z_scale )-one) ) ) - lambdah = lambdah * weight1 & - + lambda_min*( one - weight1) - ! no need to do log profile correction as klog_layr eq 2 + ! Reinstate UKV drainage flow bug here, where lambdah was not + ! enhanced as intended (and as was done in ex_coef)! + if (sg_orog_mixing == sg_shear_enh_lambda) then + if (l_rp2) then + lambdah = max ( lambda_min , & + par_mezcla_rp(rp_idx)*zh_local(i,j) ) + else + lambdah = max ( lambda_min , 0.15_r_bl*zh_local(i,j) ) + end if + if (k >= ntml_local(i,j)+2 .and. .not. l_full_lambdas) then + lambdah = lambda_min + end if + if (k <= k_log_layr) then + vkz = vkman * ( z_tq(i,j,k) - z_tq(i,j,k-1) ) + f_log = log( ( z_tq(i,j,k) + z0m_eff_gb(i,j) ) / & + ( z_tq(i,j,k-1) + z0m_eff_gb(i,j) ) ) + elh_rho(i,j,k) = vkz / ( f_log + vkz / lambdah ) + else vkz = vkman * ( z_uv(i,j,k) + z0m_eff_gb(i,j) ) elh_rho(i,j,k) = vkz / (one + vkz/lambdah ) end if - ! Reinstate UKV drainage flow bug here, where lambdah was not - ! enhanced as intended (and as was done in ex_coef)! - if (sg_orog_mixing == sg_shear_enh_lambda) then - if (l_rp2) then - lambdah = max ( lambda_min , & - par_mezcla_rp(rp_idx)*zh_local(i,j) ) - else - lambdah = max ( lambda_min , 0.15_r_bl*zh_local(i,j) ) - end if - if (k >= ntml_local(i,j)+2 .and. .not. l_full_lambdas) then - lambdah = lambda_min - end if - if (k <= k_log_layr) then - vkz = vkman * ( z_tq(i,j,k) - z_tq(i,j,k-1) ) - f_log = log( ( z_tq(i,j,k) + z0m_eff_gb(i,j) ) / & - ( z_tq(i,j,k-1) + z0m_eff_gb(i,j) ) ) - elh_rho(i,j,k) = vkz / ( f_log + vkz / lambdah ) - else - vkz = vkman * ( z_uv(i,j,k) + z0m_eff_gb(i,j) ) - elh_rho(i,j,k) = vkz / (one + vkz/lambdah ) - end if - end if - ! End of UKV bug! + end if + ! End of UKV bug! - if (BL_diag%l_elh3d) BL_diag%elh3d(i,j,k)=elh_rho(i,j,k) + if (BL_diag%l_elh3d) BL_diag%elh3d(i,j,k)=elh_rho(i,j,k) - rhokh(i,j,k) = elh_rho(i,j,k) * rhokh(i,j,k) + rhokh(i,j,k) = elh_rho(i,j,k) * rhokh(i,j,k) - end if ! test on sbl_op - end if ! test on local_fa = free_trop_layers + end if ! test on sbl_op + end if ! test on local_fa = free_trop_layers - ! Finally multiply RHOKH by dry density - if (l_mr_physics) & - rhokh(i,j,k) = rho_mix(i,j,k) * rhokh(i,j,k) + ! Finally multiply RHOKH by dry density + if (l_mr_physics) & + rhokh(i,j,k) = rho_mix(i,j,k) * rhokh(i,j,k) - end do + !end do end do end do !$OMP end do NOWAIT -!$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - !---------------------------------------------------------- - ! Use local NTML if significantly higher (to allow for - ! local averaging) than the non-local or if the non-local - ! is on the ground (=1) - !---------------------------------------------------------- - if ( .not. cumulus(i,j) .and. & - ( ntml_local(i,j) > ntml_nl(i,j)+1 & - .or. ntml_nl(i,j) == 1 ) ) then - ntml(i,j) = ntml_local(i,j) - sml_disc_inv(i,j) = 0 ! reset flag for subgrid inversion - else - ntml(i,j) = ntml_nl(i,j) - end if - !---------------------------------------------------------- - ! If local NTML is higher than NTDSC then ignore DSC layer - ! for diagnostics but keep mixing associated with it - !---------------------------------------------------------- - if ( ntml_local(i,j) > ntdsc(i,j)+1 ) then - dsc_disc_inv(i,j) = 0 - ntdsc(i,j) = 0 - nbdsc(i,j) = 0 - zhsc(i,j) = zero - dsc(i,j) = .false. - coupled(i,j) = .false. - end if - end do +!do j = pdims%j_start, pdims%j_end +!$OMP do SCHEDULE(STATIC) +do i = pdims%i_start, pdims%i_end + !---------------------------------------------------------- + ! Use local NTML if significantly higher (to allow for + ! local averaging) than the non-local or if the non-local + ! is on the ground (=1) + !---------------------------------------------------------- + if ( .not. cumulus(i,j) .and. & + ( ntml_local(i,j) > ntml_nl(i,j)+1 & + .or. ntml_nl(i,j) == 1 ) ) then + ntml(i,j) = ntml_local(i,j) + sml_disc_inv(i,j) = 0 ! reset flag for subgrid inversion + else + ntml(i,j) = ntml_nl(i,j) + end if + !---------------------------------------------------------- + ! If local NTML is higher than NTDSC then ignore DSC layer + ! for diagnostics but keep mixing associated with it + !---------------------------------------------------------- + if ( ntml_local(i,j) > ntdsc(i,j)+1 ) then + dsc_disc_inv(i,j) = 0 + ntdsc(i,j) = 0 + nbdsc(i,j) = 0 + zhsc(i,j) = zero + dsc(i,j) = .false. + coupled(i,j) = .false. + end if +!end do end do !$OMP end do !$OMP end PARALLEL @@ -2171,50 +2214,51 @@ subroutine bdy_expl2 ( & rhokm,rhokh,rhokmz,rhokhz,rhokm_top,rhokh_top,tke_loc & ) +!do j = pdims%j_start, pdims%j_end if (BL_diag%l_weight1d) then !$OMP PARALLEL do & !$OMP SCHEDULE(STATIC) & !$OMP DEFAULT(none) & -!$OMP private(k,j,i) & -!$OMP SHARED(bl_levels,pdims,BL_diag,weight_1dbl) +!$OMP private(k,i) & +!$OMP SHARED(j,bl_levels,pdims,BL_diag,weight_1dbl) do k = 1, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - BL_diag%weight1d(i,j,k)=weight_1dbl(i,j,k) - end do + do i = pdims%i_start, pdims%i_end + BL_diag%weight1d(i,j,k)=weight_1dbl(i,j,k) end do + !end do end do !$OMP end PARALLEL do end if + +!do j = pdims%j_start, pdims%j_end if (BL_diag%l_rhokm) then !$OMP PARALLEL do & !$OMP SCHEDULE(STATIC) & !$OMP DEFAULT(none) & -!$OMP private(k,j,i) & -!$OMP SHARED(bl_levels,pdims,BL_diag,rhokm) +!$OMP private(k,i) & +!$OMP SHARED(j,bl_levels,pdims,BL_diag,rhokm) do k = 1, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - BL_diag%rhokm(i,j,k)=rhokm(i,j,k) - end do + do i = pdims%i_start, pdims%i_end + BL_diag%rhokm(i,j,k)=rhokm(i,j,k) end do + !end do end do !$OMP end PARALLEL do end if +!do j = pdims%j_start, pdims%j_end if (BL_diag%l_rhokh) then !$OMP PARALLEL do & !$OMP SCHEDULE(STATIC) & !$OMP DEFAULT(none) & -!$OMP private(k,j,i) & -!$OMP SHARED(bl_levels,pdims,BL_diag,rhokh) +!$OMP private(k,i) & +!$OMP SHARED(j,bl_levels,pdims,BL_diag,rhokh) do k = 1, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - BL_diag%rhokh(i,j,k)=rhokh(i,j,k) - end do + do i = pdims%i_start, pdims%i_end + BL_diag%rhokh(i,j,k)=rhokh(i,j,k) end do +! end do end do !$OMP end PARALLEL do end if @@ -2230,23 +2274,23 @@ subroutine bdy_expl2 ( & ! BL_diag%tke currently contains (the reciprocal of) the non-local ! (SML and DSC) mixed layer timescale, calculated in excf_nl +!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(none) & -!$OMP SHARED(BL_diag, rho_wet_tq, rhokm, dbdz, bl_levels, pdims) & -!$OMP private(i, j, k, recip_time_sbl, recip_time_cbl) +!$OMP SHARED(j, BL_diag, rho_wet_tq, rhokm, dbdz, bl_levels, pdims) & +!$OMP private(i, k, recip_time_sbl, recip_time_cbl) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end + do i = pdims%i_start, pdims%i_end - ! stable timescale - recip_time_sbl = sqrt( max(dbdz(i,j,k),1.0e-5_r_bl) )/0.7_r_bl - recip_time_cbl = BL_diag%tke(i,j,k) + ! stable timescale + recip_time_sbl = sqrt( max(dbdz(i,j,k),1.0e-5_r_bl) )/0.7_r_bl + recip_time_cbl = BL_diag%tke(i,j,k) - ! TKE diagnostic - taking 5 m2/s2 as a suitable maximum - BL_diag%tke(i,j,k)= min( max_tke, & - ( rhokm(i,j,k) / rho_wet_tq(i,j,k-1) )*( & - recip_time_cbl + recip_time_sbl ) ) + ! TKE diagnostic - taking 5 m2/s2 as a suitable maximum + BL_diag%tke(i,j,k)= min( max_tke, & + ( rhokm(i,j,k) / rho_wet_tq(i,j,k-1) )*( & + recip_time_cbl + recip_time_sbl ) ) - end do + !end do end do end do !$OMP end PARALLEL do @@ -2254,57 +2298,59 @@ subroutine bdy_expl2 ( & else if (var_diags_opt == split_tke_and_inv) then ! Combine the, separately calculated, local and non-local TKE diagnostics +!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(none) & -!$OMP SHARED(BL_diag, tke_nl, tke_loc, rho_wet_tq, weight_1dbl, & +!$OMP SHARED(j, BL_diag, tke_nl, tke_loc, rho_wet_tq, weight_1dbl, & !$OMP tke_diag_fac, bl_levels, pdims) & -!$OMP private(i, j, k) +!$OMP private(i, k) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end + do i = pdims%i_start, pdims%i_end - ! TKE diagnostic - ! Currently tke_nl is strictly rho*sigma_w^2 while tke_loc is TKE - ! Assume isotropic turb so TKE_nl = 3/2 sigma_w^2 - tke_nl(i,j,k) = weight_1dbl(i,j,k) * 1.5_r_bl * & - tke_nl(i,j,k)/rho_wet_tq(i,j,k-1) - BL_diag%tke(i,j,k) = max( tke_loc(i,j,k), tke_nl(i,j,k) ) + ! TKE diagnostic + ! Currently tke_nl is strictly rho*sigma_w^2 while tke_loc is TKE + ! Assume isotropic turb so TKE_nl = 3/2 sigma_w^2 + tke_nl(i,j,k) = weight_1dbl(i,j,k) * 1.5_r_bl * & + tke_nl(i,j,k)/rho_wet_tq(i,j,k-1) + BL_diag%tke(i,j,k) = max( tke_loc(i,j,k), tke_nl(i,j,k) ) - ! Multiply by tuning factor - BL_diag%tke(i,j,k) = tke_diag_fac * BL_diag%tke(i,j,k) + ! Multiply by tuning factor + BL_diag%tke(i,j,k) = tke_diag_fac * BL_diag%tke(i,j,k) - ! Keep TKE below a sensible max value of max_tke - BL_diag%tke(i,j,k) = min( max_tke, BL_diag%tke(i,j,k) ) - ! Applying this limit can occasionally cause the length-scale - ! Km / sqrt(w_var) to become unrealistically large, since no - ! equivalent limiting is done on Km. + ! Keep TKE below a sensible max value of max_tke + BL_diag%tke(i,j,k) = min( max_tke, BL_diag%tke(i,j,k) ) + ! Applying this limit can occasionally cause the length-scale + ! Km / sqrt(w_var) to become unrealistically large, since no + ! equivalent limiting is done on Km. - end do end do + !end do end do !$OMP end PARALLEL do end if ! var_diags_opt + +!do j = pdims%j_start, pdims%j_end if ( i_bm_ez_opt == i_bm_ez_entpar ) then ! Calculate mixing-length to pass to bimodal cloud scheme, ! using Km and TKE. -!$OMP PARALLEL DEFAULT(none) private( i, j, k, weight1 ) & -!$OMP SHARED( bl_diag, mix_len_bm, rhokm, rho_wet_tq, pdims, tdims, bl_levels, & -!$OMP sml_disc_inv, ntml_nl, zhnl, dsc_disc_inv, ntdsc, zhsc, z_uv ) +!$OMP PARALLEL DEFAULT(none) private( i, k, weight1 ) & +!$OMP SHARED( j, bl_diag, mix_len_bm, rhokm, rho_wet_tq, pdims, tdims, & +!$OMP bl_levels, sml_disc_inv, ntml_nl, zhnl, dsc_disc_inv, ntdsc, & +!$OMP zhsc, z_uv ) !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels-1 - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ! TKE and rhokm are on boundary-layer theta-levels (surface at k=1), - ! mix_len_bm and rho_wet_tq are on UM theta-levels (surface at k=0). - if ( BL_diag%tke(i,j,k+1) > 0.0 ) then - mix_len_bm(i,j,k) = max( ( rhokm(i,j,k+1) / rho_wet_tq(i,j,k) ) & - / sqrt( two_thirds * BL_diag%tke(i,j,k+1) ),& - bm_tiny ) - else - mix_len_bm(i,j,k) = bm_tiny - end if - end do + do i = pdims%i_start, pdims%i_end + ! TKE and rhokm are on boundary-layer theta-levels (surface at k=1), + ! mix_len_bm and rho_wet_tq are on UM theta-levels (surface at k=0). + if ( BL_diag%tke(i,j,k+1) > 0.0 ) then + mix_len_bm(i,j,k) = max( ( rhokm(i,j,k+1) / rho_wet_tq(i,j,k) ) & + / sqrt( two_thirds * BL_diag%tke(i,j,k+1) ),& + bm_tiny ) + else + mix_len_bm(i,j,k) = bm_tiny + end if + !end do end do end do !$OMP end do @@ -2314,38 +2360,38 @@ subroutine bdy_expl2 ( & ! increment ntml/ntdsc due to the inversion rise or fall over the ! current timestep, but doesn't increment zhnl/zhsc to keep it consistent. ! We need to check for this to get sensible interpolation weights... +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( sml_disc_inv(i,j) > 0 .and. ntml_nl(i,j) > 1 ) then - k = ntml_nl(i,j) + 1 - if ( zhnl(i,j) < z_uv(i,j,k) ) k = k - 1 - if ( zhnl(i,j) > z_uv(i,j,k+1) ) k = k + 1 - weight1 = ( zhnl(i,j) - z_uv(i,j,k) )/( z_uv(i,j,k+1) - z_uv(i,j,k) ) - mix_len_bm(i,j,k) = mix_len_bm(i,j,k) & - + weight1 * max( 0.5 * mix_len_bm(i,j,k-1) & - - mix_len_bm(i,j,k), 0.0 ) - end if - if ( dsc_disc_inv(i,j) > 0 .and. ntdsc(i,j) > 1 ) then - k = ntdsc(i,j) + 1 - if ( zhsc(i,j) < z_uv(i,j,k) ) k = k - 1 - if ( zhsc(i,j) > z_uv(i,j,k+1) ) k = k + 1 - weight1 = ( zhsc(i,j) - z_uv(i,j,k) )/( z_uv(i,j,k+1) - z_uv(i,j,k) ) - mix_len_bm(i,j,k) = mix_len_bm(i,j,k) & - + weight1 * max( 0.5 * mix_len_bm(i,j,k-1) & - - mix_len_bm(i,j,k), 0.0 ) - end if - end do + do i = pdims%i_start, pdims%i_end + if ( sml_disc_inv(i,j) > 0 .and. ntml_nl(i,j) > 1 ) then + k = ntml_nl(i,j) + 1 + if ( zhnl(i,j) < z_uv(i,j,k) ) k = k - 1 + if ( zhnl(i,j) > z_uv(i,j,k+1) ) k = k + 1 + weight1 = ( zhnl(i,j) - z_uv(i,j,k) )/( z_uv(i,j,k+1) - z_uv(i,j,k) ) + mix_len_bm(i,j,k) = mix_len_bm(i,j,k) & + + weight1 * max( 0.5 * mix_len_bm(i,j,k-1) & + - mix_len_bm(i,j,k), 0.0 ) + end if + if ( dsc_disc_inv(i,j) > 0 .and. ntdsc(i,j) > 1 ) then + k = ntdsc(i,j) + 1 + if ( zhsc(i,j) < z_uv(i,j,k) ) k = k - 1 + if ( zhsc(i,j) > z_uv(i,j,k+1) ) k = k + 1 + weight1 = ( zhsc(i,j) - z_uv(i,j,k) )/( z_uv(i,j,k+1) - z_uv(i,j,k) ) + mix_len_bm(i,j,k) = mix_len_bm(i,j,k) & + + weight1 * max( 0.5 * mix_len_bm(i,j,k-1) & + - mix_len_bm(i,j,k), 0.0 ) + end if + !end do end do !$OMP end do NOWAIT ! Set values above bl_levels +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = bl_levels, tdims%k_end - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - mix_len_bm(i,j,k) = bm_tiny - end do + do i = pdims%i_start, pdims%i_end + mix_len_bm(i,j,k) = bm_tiny end do + !end do end do !$OMP end do NOWAIT !$OMP end PARALLEL @@ -2361,42 +2407,42 @@ subroutine bdy_expl2 ( & ! or any very low values being passed through ! to the microphysics turbulence call -!$OMP PARALLEL DEFAULT(none) private(k,j,i) & -!$OMP SHARED(tdims, bl_levels,pdims,BL_diag,bl_w_var,var_diags_opt) +!$OMP PARALLEL DEFAULT(none) private(k,i) & +!$OMP SHARED(j,tdims, bl_levels,pdims,BL_diag,bl_w_var,var_diags_opt) +!do j = tdims%j_start, tdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, tdims%k_end+1 - do j = tdims%j_start, tdims%j_end - do i = tdims%i_start, tdims%i_end - bl_w_var(i,j,k) = 1.0e-12_r_bl - end do + do i = tdims%i_start, tdims%i_end + bl_w_var(i,j,k) = 1.0e-12_r_bl + !end do end do end do !$OMP end do if (var_diags_opt == original_vars) then ! Original version: w_var = TKE +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if (BL_diag%tke(i,j,k) > 1.0e-12_r_bl) then - bl_w_var(i,j,k) = BL_diag%tke(i,j,k) - end if - end do + do i = pdims%i_start, pdims%i_end + if (BL_diag%tke(i,j,k) > 1.0e-12_r_bl) then + bl_w_var(i,j,k) = BL_diag%tke(i,j,k) + end if + !end do end do end do !$OMP end do else if (var_diags_opt == split_tke_and_inv) then ! Improved version: w_var = 2/3 TKE +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if (BL_diag%tke(i,j,k) > 1.0e-12_r_bl) then - bl_w_var(i,j,k) = two_thirds * BL_diag%tke(i,j,k) - end if - end do + do i = pdims%i_start, pdims%i_end + if (BL_diag%tke(i,j,k) > 1.0e-12_r_bl) then + bl_w_var(i,j,k) = two_thirds * BL_diag%tke(i,j,k) + end if + !end do end do end do !$OMP end do @@ -2409,14 +2455,14 @@ subroutine bdy_expl2 ( & ! at this point, BL_diag%tke really contains 1.5*sigma_w^2. To make it look ! a bit more like TKE near the surface, we will keep it constant below ! the max of rhokm_surf + !do j = pdims%j_start, pdims%j_end !$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(none) & -!$OMP SHARED( pdims, rhokmz, BL_diag ) private(i, j, k, kmax ) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - kmax = maxloc(rhokmz(i,j,:), dim=1) - do k = 2, kmax - BL_diag%tke(i,j,k) = BL_diag%tke(i,j,kmax) - end do +!$OMP SHARED( j, pdims, rhokmz, BL_diag ) private(i, k, kmax ) + do i = pdims%i_start, pdims%i_end + kmax = maxloc(rhokmz(i,j,:), dim=1) + do k = 2, kmax + BL_diag%tke(i,j,k) = BL_diag%tke(i,j,kmax) + !end do end do end do !$OMP end PARALLEL do @@ -2425,19 +2471,19 @@ subroutine bdy_expl2 ( & ! a bit more like TKE near the surface, we will keep it constant below ! the max of rhokm_surf (here we find the first local maximum in case there ! is a larger rhokmz in a resolved inversion +!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(none) & -!$OMP SHARED( pdims, bl_levels, rhokmz, BL_diag, tke_loc, tke_nl ) & -!$OMP private(i, j, k, kp) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - kp=2 - do while ( rhokmz(i,j,kp+1) > rhokmz(i,j,kp) .and. kp < bl_levels ) - kp = kp +1 - end do - do k = 2, kp-1 - BL_diag%tke(i,j,k) = max( tke_loc(i,j,k), tke_nl(i,j,kp) ) - end do +!$OMP SHARED( j, pdims, bl_levels, rhokmz, BL_diag, tke_loc, tke_nl ) & +!$OMP private(i, k, kp) + do i = pdims%i_start, pdims%i_end + kp=2 + do while ( rhokmz(i,j,kp+1) > rhokmz(i,j,kp) .and. kp < bl_levels ) + kp = kp +1 end do + do k = 2, kp-1 + BL_diag%tke(i,j,k) = max( tke_loc(i,j,k), tke_nl(i,j,kp) ) + end do + !end do end do !$OMP end PARALLEL do end if ! test on var_diags_opt @@ -2449,62 +2495,61 @@ subroutine bdy_expl2 ( & ! ! so copy to visc_m,h for horizontal diffusion too. ! ! Need to interpolate rhokh back to theta levels for visc_h -!$OMP PARALLEL DEFAULT(SHARED) private(i, j, k, weight1, weight2, & +!$OMP PARALLEL DEFAULT(SHARED) private(i, k, weight1, weight2, & !$OMP weight3) - +! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels-1 - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( blending_option == blend_except_cu .and. & - cumulus(i,j) .and. ntdsc(i,j) == 0) then - ! pure cumulus layer so revert to Smag scheme - visc_m(i,j,k) = visc_m(i,j,k) & - *rneutml_sq(i,j,k)*fm_3d(i,j,k+1) - visc_h(i,j,k) = visc_h(i,j,k) & - *rneutml_sq(i,j,k)*fh_3d(i,j,k+1) - if (k == bl_levels-1) then - ! need to do bl_levels too apparently - ! (see non-blend code below) - visc_m(i,j,k+1) = visc_m(i,j,k+1)*rneutml_sq(i,j,k+1) - visc_h(i,j,k+1) = visc_h(i,j,k+1)*rneutml_sq(i,j,k+1) - end if - else - visc_m(i,j,k) = rhokm(i,j,k+1)/rho_wet_tq(i,j,k) + do i = pdims%i_start, pdims%i_end + if ( blending_option == blend_except_cu .and. & + cumulus(i,j) .and. ntdsc(i,j) == 0) then + ! pure cumulus layer so revert to Smag scheme + visc_m(i,j,k) = visc_m(i,j,k) & + *rneutml_sq(i,j,k)*fm_3d(i,j,k+1) + visc_h(i,j,k) = visc_h(i,j,k) & + *rneutml_sq(i,j,k)*fh_3d(i,j,k+1) + if (k == bl_levels-1) then + ! need to do bl_levels too apparently + ! (see non-blend code below) + visc_m(i,j,k+1) = visc_m(i,j,k+1)*rneutml_sq(i,j,k+1) + visc_h(i,j,k+1) = visc_h(i,j,k+1)*rneutml_sq(i,j,k+1) end if - end do + else + visc_m(i,j,k) = rhokm(i,j,k+1)/rho_wet_tq(i,j,k) + end if + !end do end do end do !$OMP end do NOWAIT !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels-1 if ( k > 1 ) then - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if (.not. ( blending_option == blend_except_cu .and. & - cumulus(i,j) .and. ntdsc(i,j) == 0)) then - ! not a pure cumulus layer or blending_option ne - ! blend_except_cu, so use standard blending - weight1 = r_rho_levels(i,j,k+1) - r_rho_levels(i,j,k) - weight2 = r_theta_levels(i,j,k) - r_rho_levels(i,j,k) - weight3 = r_rho_levels(i,j,k+1) - r_theta_levels(i,j,k) - visc_h(i,j,k) = ( weight2*(rhokh(i,j,k+1)/rho_mix(i,j,k+1)) & - + weight3*(rhokh(i,j,k) /rho_mix(i,j,k) )) & - / weight1 - end if - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + if (.not. ( blending_option == blend_except_cu .and. & + cumulus(i,j) .and. ntdsc(i,j) == 0)) then + ! not a pure cumulus layer or blending_option ne + ! blend_except_cu, so use standard blending + weight1 = r_rho_levels(i,j,k+1) - r_rho_levels(i,j,k) + weight2 = r_theta_levels(i,j,k) - r_rho_levels(i,j,k) + weight3 = r_rho_levels(i,j,k+1) - r_theta_levels(i,j,k) + visc_h(i,j,k) = ( weight2*(rhokh(i,j,k+1)/rho_mix(i,j,k+1)) & + + weight3*(rhokh(i,j,k) /rho_mix(i,j,k) )) & + / weight1 + end if end do + !end do else - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if (.not. ( blending_option == blend_except_cu .and. & - cumulus(i,j) .and. ntdsc(i,j) == 0)) then - ! not a pure cumulus layer or blending_option ne - ! blend_except_cu, so use standard blending - visc_h(i,j,k) = rhokh_th(i,j,k+1)/rho_wet_tq(i,j,k) - end if - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + if (.not. ( blending_option == blend_except_cu .and. & + cumulus(i,j) .and. ntdsc(i,j) == 0)) then + ! not a pure cumulus layer or blending_option ne + ! blend_except_cu, so use standard blending + visc_h(i,j,k) = rhokh_th(i,j,k+1)/rho_wet_tq(i,j,k) + end if end do + !end do end if end do !$OMP end do @@ -2514,31 +2559,31 @@ subroutine bdy_expl2 ( & ! visc_m,h on in are just S and visc_m,h(k) are co-located with w(k) +!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL do SCHEDULE(STATIC) & !$OMP DEFAULT(none) & -!$OMP private(k,j,i) & -!$OMP SHARED(bl_levels,pdims,visc_m,rneutml_sq,visc_h,fh_3d,fm_3d,max_diff) +!$OMP private(k,i) & +!$OMP SHARED(j,bl_levels,pdims,visc_m,rneutml_sq,visc_h,fh_3d,fm_3d,max_diff) do k = 1, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - visc_m(i,j,k) = visc_m(i,j,k)*rneutml_sq(i,j,k) - visc_h(i,j,k) = visc_h(i,j,k)*rneutml_sq(i,j,k) - end do + do i = pdims%i_start, pdims%i_end + visc_m(i,j,k) = visc_m(i,j,k)*rneutml_sq(i,j,k) + visc_h(i,j,k) = visc_h(i,j,k)*rneutml_sq(i,j,k) + !end do end do - +!do j = pdims%j_start, pdims%j_end if ( k < bl_levels ) then - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ! stability functions are indexed with Ri, fm(k) on w(k-1) - visc_h(i,j,k) = visc_h(i,j,k)*fh_3d(i,j,k+1) - visc_m(i,j,k) = visc_m(i,j,k)*fm_3d(i,j,k+1) - ! APL why apply this cap here for implicit vertical diffusion? - ! (also applied in atm_step_phys_init for horiz diffn, - ! that actually needs it)? - visc_h(i,j,k) = min(visc_h(i,j,k),max_diff(i,j)) - visc_m(i,j,k) = min(visc_m(i,j,k),max_diff(i,j)) - end do + + do i = pdims%i_start, pdims%i_end + ! stability functions are indexed with Ri, fm(k) on w(k-1) + visc_h(i,j,k) = visc_h(i,j,k)*fh_3d(i,j,k+1) + visc_m(i,j,k) = visc_m(i,j,k)*fm_3d(i,j,k+1) + ! APL why apply this cap here for implicit vertical diffusion? + ! (also applied in atm_step_phys_init for horiz diffn, + ! that actually needs it)? + visc_h(i,j,k) = min(visc_h(i,j,k),max_diff(i,j)) + visc_m(i,j,k) = min(visc_m(i,j,k),max_diff(i,j)) + !end do end do end if end do @@ -2552,25 +2597,26 @@ subroutine bdy_expl2 ( & allocate (visc_h_rho(pdims%i_start:pdims%i_end, & pdims%j_start:pdims%j_end, bl_levels)) + +!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(none) & -!$OMP SHARED(bl_levels, pdims, r_theta_levels, r_rho_levels, visc_h_rho, & +!$OMP SHARED(j,bl_levels, pdims, r_theta_levels, r_rho_levels, visc_h_rho, & !$OMP visc_h, turb_startlev_vert, turb_endlev_vert, rhokm, visc_m, & !$OMP rho_wet_tq, rhokh, rho_mix) & -!$OMP private(i, j, k, weight1, weight2, weight3) +!$OMP private(i, k, weight1, weight2, weight3) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - weight1 = r_theta_levels(i,j,k) - r_theta_levels(i,j,k-1) - weight2 = r_theta_levels(i,j,k) - r_rho_levels(i,j,k) - weight3 = r_rho_levels(i,j,k) - r_theta_levels(i,j,k-1) - if ( k == bl_levels ) then - ! assume visc_h(bl_levels) is zero (Ri and thence f_h not defined) - visc_h_rho(i,j,k) = ( weight2/weight1 ) * visc_h(i,j,k-1) - else - visc_h_rho(i,j,k) = ( weight3/weight1 ) * visc_h(i,j,k) & - + ( weight2/weight1 ) * visc_h(i,j,k-1) - end if - end do + do i = pdims%i_start, pdims%i_end + weight1 = r_theta_levels(i,j,k) - r_theta_levels(i,j,k-1) + weight2 = r_theta_levels(i,j,k) - r_rho_levels(i,j,k) + weight3 = r_rho_levels(i,j,k) - r_theta_levels(i,j,k-1) + if ( k == bl_levels ) then + ! assume visc_h(bl_levels) is zero (Ri and thence f_h not defined) + visc_h_rho(i,j,k) = ( weight2/weight1 ) * visc_h(i,j,k-1) + else + visc_h_rho(i,j,k) = ( weight3/weight1 ) * visc_h(i,j,k) & + + ( weight2/weight1 ) * visc_h(i,j,k-1) + end if + !end do end do ! Overwrite the diffusion coefficients from the BL scheme @@ -2578,19 +2624,19 @@ subroutine bdy_expl2 ( & if (k >= turb_startlev_vert .and. & k <= turb_endlev_vert) then - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - rhokm(i,j,k) = visc_m(i,j,k-1)*rho_wet_tq(i,j,k-1) - rhokh(i,j,k) = visc_h_rho(i,j,k)*rho_mix(i,j,k) - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + rhokm(i,j,k) = visc_m(i,j,k-1)*rho_wet_tq(i,j,k-1) + rhokh(i,j,k) = visc_h_rho(i,j,k)*rho_mix(i,j,k) end do + !end do else - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - rhokm(i,j,k) = zero - rhokh(i,j,k) = zero - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + rhokm(i,j,k) = zero + rhokh(i,j,k) = zero end do + !end do end if end do !$OMP end PARALLEL do @@ -2625,95 +2671,95 @@ subroutine bdy_expl2 ( & BL_diag%zhpar=zhpar end if -!$OMP PARALLEL DEFAULT(SHARED) private(i, j, k) +!$OMP PARALLEL DEFAULT(SHARED) private(i, k) +!do j = pdims%j_start, pdims%j_end if (BL_diag%l_zht) then !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - bl_diag%zht(i,j) = max( zhnl(i,j) , zhsc(i,j) ) - if ( ntml(i,j) > ntml_nl(i,j) ) then - bl_diag%zht(i,j) = max( bl_diag%zht(i,j), zh_local(i,j) ) - end if - end do + do i = pdims%i_start, pdims%i_end + bl_diag%zht(i,j) = max( zhnl(i,j) , zhsc(i,j) ) + if ( ntml(i,j) > ntml_nl(i,j) ) then + bl_diag%zht(i,j) = max( bl_diag%zht(i,j), zh_local(i,j) ) + end if + !end do end do !$OMP end do NOWAIT end if -!$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ! Max height of BL turbulent mixing - ! PBL depth diagnostic: start with top of non-local BL - zh(i,j) = zhnl(i,j) - if ( ntml(i,j) > ntml_nl(i,j) ) then - ! Higher local K allowed so reset ZH, ZHT diagnostics - zh(i,j) = max( zh(i,j) , zh_local(i,j) ) - if (.not. l_fix_zh) then - if (forced_cu == off) then - ! going to use zhnl in tr_mix and ex_flux_uv which - ! was formerly spuriously set to zh - zhnl(i,j)=zh(i,j) - else - ! spuriously overwriting zh with zhnl - zh(i,j)=zhnl(i,j) - end if +!do j = pdims%j_start, pdims%j_end +!$OMP do SCHEDULE(STATIC) +do i = pdims%i_start, pdims%i_end + ! Max height of BL turbulent mixing + ! PBL depth diagnostic: start with top of non-local BL + zh(i,j) = zhnl(i,j) + if ( ntml(i,j) > ntml_nl(i,j) ) then + ! Higher local K allowed so reset ZH, ZHT diagnostics + zh(i,j) = max( zh(i,j) , zh_local(i,j) ) + if (.not. l_fix_zh) then + if (forced_cu == off) then + ! going to use zhnl in tr_mix and ex_flux_uv which + ! was formerly spuriously set to zh + zhnl(i,j)=zh(i,j) + else + ! spuriously overwriting zh with zhnl + zh(i,j)=zhnl(i,j) end if end if - bl_type_1(i,j) = zero - bl_type_2(i,j) = zero - bl_type_3(i,j) = zero - bl_type_4(i,j) = zero - bl_type_5(i,j) = zero - bl_type_6(i,j) = zero - bl_type_7(i,j) = zero - end do + end if + bl_type_1(i,j) = zero + bl_type_2(i,j) = zero + bl_type_3(i,j) = zero + bl_type_4(i,j) = zero + bl_type_5(i,j) = zero + bl_type_6(i,j) = zero + bl_type_7(i,j) = zero end do +!end do !$OMP end do NOWAIT +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if (dynamic_bl_diag(i,j)) then - ! ! shear-dominated, via iDynDiag option - bl_type_7(i,j) = one - else - if (.not. unstable(i,j) .and. .not. dsc(i,j) .and. & - .not. cumulus(i,j)) then - ! ! Stable b.l. - bl_type_1(i,j) = one - else if (.not. unstable(i,j) .and. dsc(i,j) .and. & - .not. cumulus(i,j)) then - ! ! Stratocumulus over a stable surface layer - bl_type_2(i,j) = one - else if (unstable(i,j) .and. .not. cumulus(i,j) .and. & - .not. dsc(i,j) ) then - ! ! Well mixed b.l. (possibly with stratocumulus) - if ( ntml(i,j) > ntml_nl(i,j) ) then - ! shear-dominated - currently identified - ! by local NTML overriding non-local - bl_type_7(i,j) = one - else - ! buoyancy-dominated - bl_type_3(i,j) = one - end if - else if (unstable(i,j) .and. dsc(i,j) .and. & - .not. cumulus(i,j)) then - ! ! Decoupled stratocumulus (not over cumulus) - bl_type_4(i,j) = one - else if (dsc(i,j) .and. cumulus(i,j)) then - ! ! Decoupled stratocumulus over cumulus - bl_type_5(i,j) = one - else if (.not. dsc(i,j) .and. cumulus(i,j)) then - ! ! Cumulus capped b.l. - bl_type_6(i,j) = one - ! Label this shallow regime as 2.0_r_bl here, to be able to identify it - ! in diagnostics_bl, but the "cumulus" stash output will still be 1.0 - if (l_shallow_cth(i,j)) bl_type_6(i,j) = 2.0_r_bl +do i = pdims%i_start, pdims%i_end + if (dynamic_bl_diag(i,j)) then + ! ! shear-dominated, via iDynDiag option + bl_type_7(i,j) = one + else + if (.not. unstable(i,j) .and. .not. dsc(i,j) .and. & + .not. cumulus(i,j)) then + ! ! Stable b.l. + bl_type_1(i,j) = one + else if (.not. unstable(i,j) .and. dsc(i,j) .and. & + .not. cumulus(i,j)) then + ! ! Stratocumulus over a stable surface layer + bl_type_2(i,j) = one + else if (unstable(i,j) .and. .not. cumulus(i,j) .and. & + .not. dsc(i,j) ) then + ! ! Well mixed b.l. (possibly with stratocumulus) + if ( ntml(i,j) > ntml_nl(i,j) ) then + ! shear-dominated - currently identified + ! by local NTML overriding non-local + bl_type_7(i,j) = one + else + ! buoyancy-dominated + bl_type_3(i,j) = one end if + else if (unstable(i,j) .and. dsc(i,j) .and. & + .not. cumulus(i,j)) then + ! ! Decoupled stratocumulus (not over cumulus) + bl_type_4(i,j) = one + else if (dsc(i,j) .and. cumulus(i,j)) then + ! ! Decoupled stratocumulus over cumulus + bl_type_5(i,j) = one + else if (.not. dsc(i,j) .and. cumulus(i,j)) then + ! ! Cumulus capped b.l. + bl_type_6(i,j) = one + ! Label this shallow regime as 2.0_r_bl here, to be able to identify it + ! in diagnostics_bl, but the "cumulus" stash output will still be 1.0 + if (l_shallow_cth(i,j)) bl_type_6(i,j) = 2.0_r_bl end if + end if - end do +! end do end do !$OMP end do NOWAIT !$OMP end PARALLEL @@ -2733,14 +2779,14 @@ subroutine bdy_expl2 ( & ftl,fqw,wtrac_bl & ) -!$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(none) private(k, j, i) & -!$OMP SHARED(bl_levels, pdims, f_ngstress, weight_1dbl) +!do j = pdims%j_start, pdims%j_end +!$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(none) private(k, i) & +!$OMP SHARED(j, bl_levels, pdims, f_ngstress, weight_1dbl) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - f_ngstress(i,j,k) = weight_1dbl(i,j,k) * f_ngstress(i,j,k) - end do + do i = pdims%i_start, pdims%i_end + f_ngstress(i,j,k) = weight_1dbl(i,j,k) * f_ngstress(i,j,k) end do + !end do end do !$OMP end PARALLEL do @@ -2751,31 +2797,31 @@ subroutine bdy_expl2 ( & ! 5.6.1 Calculate explicit surface fluxes of U and V on ! P-grid for convection scheme !----------------------------------------------------------------------- -!$OMP PARALLEL DEFAULT(none) private(j,i,l) & -!$OMP SHARED(pdims,uw0,rhokm,u_p,u_0_px,vw0,v_p,v_0_px,wstar,wthvs, & +!$OMP PARALLEL DEFAULT(none) private(i,l) & +!$OMP SHARED(j,pdims,uw0,rhokm,u_p,u_0_px,vw0,v_p,v_0_px,wstar,wthvs, & !$OMP non_local_bl,cu_over_orog,cumulus,fb_surf,zh,g,bt,l_param_conv, & !$OMP ntml,ntml_nl,ntml_local,formdrag,tau_fd_x,tau_fd_y,ntdsc,BL_diag, & !$OMP land_pts,land_index,ho2r2_orog,l_shallow,bl_type_5,bl_type_6, & !$OMP ntml_save,shallowc) !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end +!do j = pdims%j_start, pdims%j_end do i = pdims%i_start, pdims%i_end uw0(i,j) = -rhokm(i,j,1) * & ( u_p(i,j,1) - u_0_px(i,j) ) vw0(i,j) = -rhokm(i,j,1) * & ( v_p(i,j,1) - v_0_px(i,j) ) end do -end do +!end do !$OMP end do NOWAIT if (formdrag == explicit_stress) then +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - uw0(i,j) = uw0(i,j) - tau_fd_x(i,j,1) - vw0(i,j) = vw0(i,j) - tau_fd_y(i,j,1) - end do + do i = pdims%i_start, pdims%i_end + uw0(i,j) = uw0(i,j) - tau_fd_x(i,j,1) + vw0(i,j) = vw0(i,j) - tau_fd_y(i,j,1) end do + ! end do !$OMP end do NOWAIT end if @@ -2787,62 +2833,62 @@ subroutine bdy_expl2 ( & if (non_local_bl == on) then +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( cumulus(i,j) ) then - if (.not. l_param_conv) then - ntml(i,j) = max( 2, ntml_nl(i,j) - 1 ) - end if - else - ntml(i,j) = max( ntml_nl(i,j) , ntdsc(i,j) ) + do i = pdims%i_start, pdims%i_end + if ( cumulus(i,j) ) then + if (.not. l_param_conv) then + ntml(i,j) = max( 2, ntml_nl(i,j) - 1 ) end if - end do + else + ntml(i,j) = max( ntml_nl(i,j) , ntdsc(i,j) ) + end if + !end do end do !$OMP end do NOWAIT else ! non-local scheme not being used so always use local scheme values +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ntml(i,j) = ntml_local(i,j) - end do + do i = pdims%i_start, pdims%i_end + ntml(i,j) = ntml_local(i,j) end do + !end do !$OMP end do NOWAIT end if +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - wstar(i,j) = zero - wthvs(i,j) = zero - cu_over_orog(i,j) = zero - if ( cumulus(i,j) ) then - if ( fb_surf(i,j) > zero ) then - wstar(i,j) = ( zh(i,j)*fb_surf(i,j) )**one_third - wthvs(i,j) = fb_surf(i,j) / ( g * bt(i,j,1) ) - end if - wstar(i,j) = max( 0.1_r_bl, wstar(i,j) ) +do i = pdims%i_start, pdims%i_end + wstar(i,j) = zero + wthvs(i,j) = zero + cu_over_orog(i,j) = zero + if ( cumulus(i,j) ) then + if ( fb_surf(i,j) > zero ) then + wstar(i,j) = ( zh(i,j)*fb_surf(i,j) )**one_third + wthvs(i,j) = fb_surf(i,j) / ( g * bt(i,j,1) ) end if - ! Limit explicitly calculated surface stresses - ! to a physically plausible level. - if ( uw0(i,j) >= 5.0_r_bl ) then - uw0(i,j) = 5.0_r_bl - else if ( uw0(i,j) <= -5.0_r_bl ) then - uw0(i,j) = -5.0_r_bl - end if - if ( vw0(i,j) >= 5.0_r_bl ) then - vw0(i,j) = 5.0_r_bl - else if ( vw0(i,j) <= -5.0_r_bl ) then - vw0(i,j) = -5.0_r_bl - end if - if (BL_diag%l_wstar .and. (fb_surf(i,j) >zero)) then - BL_diag%wstar(i,j)= (zh(i,j)*fb_surf(i,j))**one_third - end if - end do + wstar(i,j) = max( 0.1_r_bl, wstar(i,j) ) + end if + ! Limit explicitly calculated surface stresses + ! to a physically plausible level. + if ( uw0(i,j) >= 5.0_r_bl ) then + uw0(i,j) = 5.0_r_bl + else if ( uw0(i,j) <= -5.0_r_bl ) then + uw0(i,j) = -5.0_r_bl + end if + if ( vw0(i,j) >= 5.0_r_bl ) then + vw0(i,j) = 5.0_r_bl + else if ( vw0(i,j) <= -5.0_r_bl ) then + vw0(i,j) = -5.0_r_bl + end if + if (BL_diag%l_wstar .and. (fb_surf(i,j) >zero)) then + BL_diag%wstar(i,j)= (zh(i,j)*fb_surf(i,j))**one_third + end if end do +!end do !$OMP end do if (l_param_conv) then @@ -2872,13 +2918,13 @@ subroutine bdy_expl2 ( & ! if cumulus is still true in order to trigger convection ! at the correct level +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( .not. cumulus(i,j) ) l_shallow(i,j) = .false. - if ( cumulus(i,j) ) ntml(i,j) = ntml_save(i,j) - end do + do i = pdims%i_start, pdims%i_end + if ( .not. cumulus(i,j) ) l_shallow(i,j) = .false. + if ( cumulus(i,j) ) ntml(i,j) = ntml_save(i,j) end do + !end do !$OMP end do NOWAIT end if ! (l_param_conv) @@ -2887,16 +2933,16 @@ subroutine bdy_expl2 ( & ! 0.0 if .not. CUMULUS !----------------------------------------------------------------------- +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( cumulus(i,j) .and. l_shallow(i,j) ) then - shallowc(i,j) = one - else - shallowc(i,j) = zero - end if - end do +do i = pdims%i_start, pdims%i_end + if ( cumulus(i,j) .and. l_shallow(i,j) ) then + shallowc(i,j) = one + else + shallowc(i,j) = zero + end if end do +!end do !$OMP end do !$OMP end PARALLEL @@ -2928,42 +2974,198 @@ subroutine bdy_expl2 ( & ! ftl is on rho-level k ! rhokm is on theta-level K-1 + ! Sub divide tdims domain len over the number of threads. + tdims_omp_block = ceiling(real(tdims%i_len)/max_threads) + ! Pick the smallest option. + ! Work should be split accross the threads, this provides a cap. + ! In the occurence where the tdims domain len is very small, it will be seleted. + ! Otherwise bl_segment_size will be selected. + tdims_seg_block = min(bl_segment_size, tdims_omp_block, tdims%i_len) ! level 2 needs special treatment because of the surface - -!$OMP PARALLEL DEFAULT(SHARED) & -!$OMP private(k,j,i,km,kp,sh,exner,sgm,qsw_arr,weight1,weight2,weight3,var_fac,& -!$OMP sl_var,qw_var,sl_qw,delta_x) + !do j = tdims%j_start, tdims%j_end + !$OMP PARALLEL DEFAULT(SHARED) & + !$OMP private(k,i,ii,km,kp,sh,exner,sgm,qsw_arr,weight1,weight2,weight3, & + !$OMP var_fac,sl_var,qw_var,sl_qw,delta_x, & + !$OMP seg_slice_start, seg_slice_end, bl_segment_range) k = 2 -!$OMP do SCHEDULE(STATIC) - do j = tdims%j_start, tdims%j_end +! !$OMP do SCHEDULE(STATIC) + !$OMP do SCHEDULE(STATIC) + do ii = tdims%i_start, tdims%i_end, tdims_seg_block + seg_slice_start = ii + seg_slice_end = MIN(((ii+tdims_seg_block)-1), tdims%i_end) + bl_segment_range = (seg_slice_end - seg_slice_start) + 1 if ( l_mr_physics ) then - call qsat_wat_mix(qsw_arr,tl(:,j,k-1),p_theta_levels(:,j,k-1),tdims%i_len) + call qsat_wat_mix(qsw_arr(seg_slice_start:seg_slice_end, j, k-1), & + tl(seg_slice_start:seg_slice_end,j,k-1), & + p_theta_levels(seg_slice_start:seg_slice_end,j,k-1), & + bl_segment_range ) else - call qsat_wat(qsw_arr,tl(:,j,k-1),p_theta_levels(:,j,k-1),tdims%i_len) + call qsat_wat(qsw_arr(seg_slice_start:seg_slice_end, j, k-1), & + tl(seg_slice_start:seg_slice_end,j,k-1), & + p_theta_levels(seg_slice_start:seg_slice_end,j,k-1), & + bl_segment_range ) end if - if (var_diags_opt == original_vars) then + end do ! ii + !$OMP end do + + if (var_diags_opt == original_vars) then + !$OMP do SCHEDULE(STATIC) + do i = tdims%i_start, tdims%i_end + ! calculate sh, don't interpolate with the surface flux or divide by 0 + if ( BL_diag%tke(i,j,k) > 1.0e-10_r_bl ) then + sh = rhokh(i,j,k) & + / ( rho_wet_tq(i,j,k-1) * elm(i,j,k) & + * sqrt( 2.0_r_bl * BL_diag%tke(i,j,k) ) ) + else + sh = small_sh + end if + if ( sh < small_sh ) sh = small_sh + ! calculate exner + exner = ( p_theta_levels(i,j,k-1) / pref )**kappa + ! calculate the variance, use gradient interpolated between levs 1 and 2 + sgm(i) = a_dqsdt(i,j,k-1)**2 * exner**2 * b2 * sh & + * elm(i,j,k)**2 * dsldz(i,j,k)**2 & + + a_qs(i,j,k-1)**2 * b2 * sh & + * elm(i,j,k)**2 * dqwdz(i,j,k)**2 & + - 2.0_r_bl * a_qs(i,j,k-1) * a_dqsdt(i,j,k-1) * exner * b2 & + * sh * elm(i,j,k)**2 * dsldz(i,j,k) * dqwdz(i,j,k) + if (BL_diag%l_slvar) then + BL_diag%slvar(i,j,k-1) = b2 * sh * elm(i,j,k)**2 * dsldz(i,j,k)**2 + end if + if (BL_diag%l_qwvar) then + BL_diag%qwvar(i,j,k-1) = b2 * sh * elm(i,j,k)**2 * dqwdz(i,j,k)**2 + end if + if (BL_diag%l_slqw) then + BL_diag%slqw(i,j,k-1) = b2 * sh * elm(i,j,k)**2 * dsldz(i,j,k) & + * dqwdz(i,j,k) + end if + ! do this for safety, not sure if it's really needed + sgm(i) = sqrt ( max( sgm(i), zero ) ) + end do + !$OMP end do + + else ! var_diags_opt + !$OMP do SCHEDULE(STATIC) + do i = tdims%i_start, tdims%i_end + ! calculate the variance + sl_var = zero + qw_var = zero + sl_qw = zero + if ( BL_diag%tke(i,j,k) > small_tke ) then + ! vertical interpolation weights + weight1 = r_rho_levels(i,j,k) - r_rho_levels(i,j,k-1) + weight2 = r_theta_levels(i,j,k-1) - r_rho_levels(i,j,k-1) + weight3 = r_rho_levels(i,j,k) - r_theta_levels(i,j,k-1) + ! var_fac=b2*L/sqrt(TKE) + var_fac = b2 * rhokm(i,j,k) / ( weight1 * BL_diag%tke(i,j,k) * & + rho_wet_tq(i,j,k-1) * rho_mix_tq(i,j,k-1) ) + ! Note that flux*gradient can be negative so the absolute values + ! are used + qw_var= abs( -var_fac*( weight2 * fqw(i,j,k) * dqwdz(i,j,k) + & + weight3 * fqw(i,j,k-1)* dqwdz(i,j,k-1) ) ) + sl_var= abs( -var_fac*( weight2 * ftl(i,j,k) * dsldz(i,j,k) + & + weight3 * ftl(i,j,k-1)* dsldz(i,j,k-1) ) ) + sl_qw = - one_half * var_fac*( & + weight2*( ftl(i,j,k)*dqwdz(i,j,k) + fqw(i,j,k)*dsldz(i,j,k) ) + & + weight3*( ftl(i,j,k-1)*dqwdz(i,j,k-1) + fqw(i,j,k-1)*dsldz(i,j,k-1))& + ) + if (BL_diag%l_slvar) then + BL_diag%slvar(i,j,k-1) = sl_var + end if + if (BL_diag%l_qwvar) then + BL_diag%qwvar(i,j,k-1) = qw_var + end if + if (BL_diag%l_slqw) then + BL_diag%slqw(i,j,k-1) = sl_qw + end if + end if + exner = ( p_theta_levels(i,j,k-1) / pref )**kappa + sgm(i) = a_dqsdt(i,j,k-1)**2 * exner**2 * sl_var & + + a_qs(i,j,k-1)**2 * qw_var & + - 2.0_r_bl * a_qs(i,j,k-1) * a_dqsdt(i,j,k-1) * exner * sl_qw + ! do this for safety, not sure if it's really needed + sgm(i) = sqrt ( max( sgm(i), zero ) ) + end do !i + !$OMP end do + end if ! var_diags_opt + if (i_rhcpt == rhcpt_tke_based) then + !$OMP do SCHEDULE(STATIC) + do i = tdims%i_start, tdims%i_end + ! calculate rhcrit, with appropriate limits + ! calculate grid-box size, just take surface for simplicity + delta_x = sqrt( r_theta_levels(i,j,k-1) * delta_lambda * & + r_theta_levels(i,j,k-1) * delta_phi * & + cos_theta_latitude(i,j) ) + ! max limit, based on curve fitted to aircraft observations + max_rhcpt(i,j) = min( 0.99_r_bl, 0.997_r_bl - 0.0078_r_bl * & + log( 0.001_r_bl * delta_x ) ) + ! min limit, based on curve fitted to aircraft observations + min_rhcpt(i,j) = max( 0.6_r_bl, 0.846_r_bl - 0.065_r_bl * & + log( 0.001_r_bl * delta_x ) ) + ! full expression + rhcpt(i,j,k-1) = min( max_rhcpt(i,j), max( min_rhcpt(i,j), & + one - root6 * sgm(i) / (a_qs(i,j,k-1) * qsw_arr(i,j,k-1)))) + end do !i + !$OMP end do + end if + !end do !j +!!$OMP end do + + ! remaining bl levels use proper interpolation + + do k = 3, bl_levels-1 +!! $OMP do SCHEDULE(STATIC) + !do j = tdims%j_start, tdims%j_end + !$OMP do SCHEDULE(STATIC) + do ii = tdims%i_start, tdims%i_end, tdims_seg_block + seg_slice_start = ii + seg_slice_end = MIN(((ii+tdims_seg_block)-1), tdims%i_end) + bl_segment_range = (seg_slice_end - seg_slice_start) + 1 + + if ( l_mr_physics ) then + !dir$ safe_address + call qsat_wat_mix(qsw_arr(seg_slice_start:seg_slice_end, j, k-1), & + tl(seg_slice_start:seg_slice_end,j,k-1), & + p_theta_levels(seg_slice_start:seg_slice_end,j,k-1), & + bl_segment_range) + else + !dir$ safe_address + call qsat_wat(qsw_arr(seg_slice_start:seg_slice_end, j, k-1), & + tl(seg_slice_start:seg_slice_end,j,k-1), & + p_theta_levels(seg_slice_start:seg_slice_end,j,k-1), & + bl_segment_range) + end if + + end do ! ii + !$OMP end do NOWAIT + + if (var_diags_opt == original_vars) then + !$OMP do SCHEDULE(STATIC) do i = tdims%i_start, tdims%i_end - ! calculate sh, don't interpolate with the surface flux or divide by 0 + ! calculate sh, don't divide by 0 if ( BL_diag%tke(i,j,k) > 1.0e-10_r_bl ) then - sh = rhokh(i,j,k) & - / ( rho_wet_tq(i,j,k-1) * elm(i,j,k) & - * sqrt( 2.0_r_bl * BL_diag%tke(i,j,k) ) ) + weight1 = r_rho_levels(i,j,k) - r_rho_levels(i,j,k-1) + weight2 = r_theta_levels(i,j,k-1)-r_rho_levels(i,j,k-1) + weight3 = r_rho_levels(i,j,k) - r_theta_levels(i,j,k-1) + sh = ( weight2 * rhokh(i,j,k) + weight3 * rhokh(i,j,k-1) ) & + / ( weight1 * rho_wet_tq(i,j,k-1) * elm(i,j,k) & + * sqrt( 2.0_r_bl * BL_diag%tke(i,j,k) ) ) else sh = small_sh end if if ( sh < small_sh ) sh = small_sh ! calculate exner exner = ( p_theta_levels(i,j,k-1) / pref )**kappa - ! calculate the variance, use gradient interpolated between levs 1 and 2 - sgm(i) = a_dqsdt(i,j,k-1)**2 * exner**2 * b2 * sh & - * elm(i,j,k)**2 * dsldz(i,j,k)**2 & - + a_qs(i,j,k-1)**2 * b2 * sh & - * elm(i,j,k)**2 * dqwdz(i,j,k)**2 & - - 2.0_r_bl * a_qs(i,j,k-1) * a_dqsdt(i,j,k-1) * exner * b2 & - * sh * elm(i,j,k)**2 * dsldz(i,j,k) * dqwdz(i,j,k) + ! calculate the variance + sgm(i) = a_dqsdt(i,j,k-1)**2 * exner**2 * b2 * sh & + * elm(i,j,k)**2 * dsldzm(i,j,k)**2 & + + a_qs(i,j,k-1)**2 * b2 * sh & + * elm(i,j,k)**2 * dqwdzm(i,j,k)**2 & + - 2.0_r_bl * a_qs(i,j,k-1) * a_dqsdt(i,j,k-1) * exner * b2 & + * sh * elm(i,j,k)**2 * dsldzm(i,j,k) * dqwdzm(i,j,k) if (BL_diag%l_slvar) then BL_diag%slvar(i,j,k-1) = b2 * sh * elm(i,j,k)**2 * dsldz(i,j,k)**2 end if @@ -2971,37 +3173,55 @@ subroutine bdy_expl2 ( & BL_diag%qwvar(i,j,k-1) = b2 * sh * elm(i,j,k)**2 * dqwdz(i,j,k)**2 end if if (BL_diag%l_slqw) then - BL_diag%slqw(i,j,k-1) = b2 * sh * elm(i,j,k)**2 * dsldz(i,j,k) & - * dqwdz(i,j,k) + BL_diag%slqw(i,j,k-1) = b2 * sh * elm(i,j,k)**2 * dsldzm(i,j,k) & + * dqwdz(i,j,k) end if - ! do this for safety, not sure if it's really needed - sgm(i) = sqrt ( max( sgm(i), zero ) ) - end do + end do !i + !$OMP end do NOWAIT + if (i_rhcpt == rhcpt_tke_based) then + !$OMP do SCHEDULE(STATIC) + do i = tdims%i_start, tdims%i_end + ! do this for safety, not sure if it's really needed + sgm(i) = sqrt ( max( sgm(i), zero ) ) + ! calculate rhcrit, with appropriate limits + rhcpt(i,j,k-1) = min( max_rhcpt(i,j), max( min_rhcpt(i,j), & + one - root6 * sgm(i) / (a_qs(i,j,k-1) * qsw_arr(i,j,k-1)))) + end do !i + !$OMP end do NOWAIT + end if else ! var_diags_opt - + !$OMP do SCHEDULE(STATIC) do i = tdims%i_start, tdims%i_end ! calculate the variance sl_var = zero qw_var = zero - sl_qw = zero + sl_qw = zero if ( BL_diag%tke(i,j,k) > small_tke ) then ! vertical interpolation weights weight1 = r_rho_levels(i,j,k) - r_rho_levels(i,j,k-1) weight2 = r_theta_levels(i,j,k-1) - r_rho_levels(i,j,k-1) weight3 = r_rho_levels(i,j,k) - r_theta_levels(i,j,k-1) - ! var_fac=b2*L/sqrt(TKE) - var_fac = b2 * rhokm(i,j,k) / ( weight1 * BL_diag%tke(i,j,k) * & - rho_wet_tq(i,j,k-1) * rho_mix_tq(i,j,k-1) ) + var_fac = b2 * rhokm(i,j,k) / ( weight1*BL_diag%tke(i,j,k) * & + rho_wet_tq(i,j,k-1) * rho_mix_tq(i,j,k-1)) + kp=k + km=k-1 + ! Don't use level ntml (or ntdsc) if disc_inv=2 as this indicates + ! the inversion has just risen and the gradients between ntml and + ! ntml-1 are likely to give excessive variances + if ( (kp == ntml(i,j) .and. sml_disc_inv(i,j) == 2) .or. & + (kp == ntdsc(i,j) .and. dsc_disc_inv(i,j) == 2) ) kp = km + if ( (km == ntml(i,j) .and. sml_disc_inv(i,j) == 2) .or. & + (km == ntdsc(i,j) .and. dsc_disc_inv(i,j) == 2) ) km = kp ! Note that flux*gradient can be negative so the absolute values ! are used - qw_var= abs( -var_fac*( weight2 * fqw(i,j,k) * dqwdz(i,j,k) + & - weight3 * fqw(i,j,k-1)* dqwdz(i,j,k-1) ) ) - sl_var= abs( -var_fac*( weight2 * ftl(i,j,k) * dsldz(i,j,k) + & - weight3 * ftl(i,j,k-1)* dsldz(i,j,k-1) ) ) - sl_qw = - one_half * var_fac*( & - weight2*( ftl(i,j,k)*dqwdz(i,j,k) + fqw(i,j,k)*dsldz(i,j,k) ) + & - weight3*( ftl(i,j,k-1)*dqwdz(i,j,k-1) + fqw(i,j,k-1)*dsldz(i,j,k-1))& + qw_var= abs( -var_fac*( weight2 * fqw(i,j,kp) * dqwdz(i,j,kp) + & + weight3 * fqw(i,j,km) * dqwdz(i,j,km) ) ) + sl_var= abs( -var_fac*( weight2 * ftl(i,j,kp) * dsldz(i,j,kp) + & + weight3 * ftl(i,j,km) * dsldz(i,j,km) ) ) + sl_qw = - one_half * var_fac*( & + weight2*( ftl(i,j,k)*dqwdz(i,j,k) + fqw(i,j,k)*dsldz(i,j,k) ) + & + weight3*( ftl(i,j,k-1)*dqwdz(i,j,k-1) + fqw(i,j,k-1)*dsldz(i,j,k-1))& ) if (BL_diag%l_slvar) then BL_diag%slvar(i,j,k-1) = sl_var @@ -3010,176 +3230,47 @@ subroutine bdy_expl2 ( & BL_diag%qwvar(i,j,k-1) = qw_var end if if (BL_diag%l_slqw) then - BL_diag%slqw(i,j,k-1) = sl_qw + BL_diag%slqw(i,j,k-1) = sl_qw end if end if exner = ( p_theta_levels(i,j,k-1) / pref )**kappa - sgm(i) = a_dqsdt(i,j,k-1)**2 * exner**2 * sl_var & - + a_qs(i,j,k-1)**2 * qw_var & - - 2.0_r_bl * a_qs(i,j,k-1) * a_dqsdt(i,j,k-1) * exner * sl_qw - ! do this for safety, not sure if it's really needed - sgm(i) = sqrt ( max( sgm(i), zero ) ) - + sgm(i) = a_dqsdt(i,j,k-1)**2 * exner**2 * sl_var & + + a_qs(i,j,k-1)**2 * qw_var & + - 2.0_r_bl * a_qs(i,j,k-1) * a_dqsdt(i,j,k-1) * exner * sl_qw end do !i - end if ! var_diags_opt - - if (i_rhcpt == rhcpt_tke_based) then - do i = tdims%i_start, tdims%i_end - ! calculate rhcrit, with appropriate limits - ! calculate grid-box size, just take surface for simplicity - delta_x = sqrt( r_theta_levels(i,j,k-1) * delta_lambda * & - r_theta_levels(i,j,k-1) * delta_phi * & - cos_theta_latitude(i,j) ) - ! max limit, based on curve fitted to aircraft observations - max_rhcpt(i,j) = min( 0.99_r_bl, 0.997_r_bl - 0.0078_r_bl * & - log( 0.001_r_bl * delta_x ) ) - ! min limit, based on curve fitted to aircraft observations - min_rhcpt(i,j) = max( 0.6_r_bl, 0.846_r_bl - 0.065_r_bl * & - log( 0.001_r_bl * delta_x ) ) - ! full expression - rhcpt(i,j,k-1) = min( max_rhcpt(i,j), max( min_rhcpt(i,j), & - one - root6 * sgm(i) / (a_qs(i,j,k-1) * qsw_arr(i)))) - end do !i - end if - end do !j -!$OMP end do - - ! remaining bl levels use proper interpolation - - do k = 3, bl_levels-1 -!$OMP do SCHEDULE(STATIC) - do j = tdims%j_start, tdims%j_end - - if ( l_mr_physics ) then - call qsat_wat_mix(qsw_arr,tl(:,j,k-1),p_theta_levels(:,j,k-1), & - tdims%i_len) - else - call qsat_wat(qsw_arr,tl(:,j,k-1),p_theta_levels(:,j,k-1),tdims%i_len) - end if - - if (var_diags_opt == original_vars) then - - do i = tdims%i_start, tdims%i_end - ! calculate sh, don't divide by 0 - if ( BL_diag%tke(i,j,k) > 1.0e-10_r_bl ) then - weight1 = r_rho_levels(i,j,k) - r_rho_levels(i,j,k-1) - weight2 = r_theta_levels(i,j,k-1)-r_rho_levels(i,j,k-1) - weight3 = r_rho_levels(i,j,k) - r_theta_levels(i,j,k-1) - sh = ( weight2 * rhokh(i,j,k) + weight3 * rhokh(i,j,k-1) ) & - / ( weight1 * rho_wet_tq(i,j,k-1) * elm(i,j,k) & - * sqrt( 2.0_r_bl * BL_diag%tke(i,j,k) ) ) - else - sh = small_sh - end if - if ( sh < small_sh ) sh = small_sh - ! calculate exner - exner = ( p_theta_levels(i,j,k-1) / pref )**kappa - ! calculate the variance - sgm(i) = a_dqsdt(i,j,k-1)**2 * exner**2 * b2 * sh & - * elm(i,j,k)**2 * dsldzm(i,j,k)**2 & - + a_qs(i,j,k-1)**2 * b2 * sh & - * elm(i,j,k)**2 * dqwdzm(i,j,k)**2 & - - 2.0_r_bl * a_qs(i,j,k-1) * a_dqsdt(i,j,k-1) * exner * b2 & - * sh * elm(i,j,k)**2 * dsldzm(i,j,k) * dqwdzm(i,j,k) - if (BL_diag%l_slvar) then - BL_diag%slvar(i,j,k-1) = b2 * sh * elm(i,j,k)**2 * dsldz(i,j,k)**2 - end if - if (BL_diag%l_qwvar) then - BL_diag%qwvar(i,j,k-1) = b2 * sh * elm(i,j,k)**2 * dqwdz(i,j,k)**2 - end if - if (BL_diag%l_slqw) then - BL_diag%slqw(i,j,k-1) = b2 * sh * elm(i,j,k)**2 * dsldzm(i,j,k) & - * dqwdz(i,j,k) - end if - end do !i - if (i_rhcpt == rhcpt_tke_based) then - do i = tdims%i_start, tdims%i_end - ! do this for safety, not sure if it's really needed - sgm(i) = sqrt ( max( sgm(i), zero ) ) - ! calculate rhcrit, with appropriate limits - rhcpt(i,j,k-1) = min( max_rhcpt(i,j), max( min_rhcpt(i,j), & - one - root6 * sgm(i) / (a_qs(i,j,k-1) * qsw_arr(i)))) - end do !i - end if - - else ! var_diags_opt - + !$OMP end do NOWAIT + if (i_rhcpt == rhcpt_tke_based) then + !$OMP do SCHEDULE(STATIC) do i = tdims%i_start, tdims%i_end - ! calculate the variance - sl_var = zero - qw_var = zero - sl_qw = zero - if ( BL_diag%tke(i,j,k) > small_tke ) then - ! vertical interpolation weights - weight1 = r_rho_levels(i,j,k) - r_rho_levels(i,j,k-1) - weight2 = r_theta_levels(i,j,k-1) - r_rho_levels(i,j,k-1) - weight3 = r_rho_levels(i,j,k) - r_theta_levels(i,j,k-1) - var_fac = b2 * rhokm(i,j,k) / ( weight1*BL_diag%tke(i,j,k) * & - rho_wet_tq(i,j,k-1) * rho_mix_tq(i,j,k-1)) - kp=k - km=k-1 - ! Don't use level ntml (or ntdsc) if disc_inv=2 as this indicates - ! the inversion has just risen and the gradients between ntml and - ! ntml-1 are likely to give excessive variances - if ( (kp == ntml(i,j) .and. sml_disc_inv(i,j) == 2) .or. & - (kp == ntdsc(i,j) .and. dsc_disc_inv(i,j) == 2) ) kp = km - if ( (km == ntml(i,j) .and. sml_disc_inv(i,j) == 2) .or. & - (km == ntdsc(i,j) .and. dsc_disc_inv(i,j) == 2) ) km = kp - ! Note that flux*gradient can be negative so the absolute values - ! are used - qw_var= abs( -var_fac*( weight2 * fqw(i,j,kp) * dqwdz(i,j,kp) + & - weight3 * fqw(i,j,km) * dqwdz(i,j,km) ) ) - sl_var= abs( -var_fac*( weight2 * ftl(i,j,kp) * dsldz(i,j,kp) + & - weight3 * ftl(i,j,km) * dsldz(i,j,km) ) ) - sl_qw = - one_half * var_fac*( & - weight2*( ftl(i,j,k)*dqwdz(i,j,k) + fqw(i,j,k)*dsldz(i,j,k) ) + & - weight3*( ftl(i,j,k-1)*dqwdz(i,j,k-1) + fqw(i,j,k-1)*dsldz(i,j,k-1))& - ) - if (BL_diag%l_slvar) then - BL_diag%slvar(i,j,k-1) = sl_var - end if - if (BL_diag%l_qwvar) then - BL_diag%qwvar(i,j,k-1) = qw_var - end if - if (BL_diag%l_slqw) then - BL_diag%slqw(i,j,k-1) = sl_qw - end if - end if - exner = ( p_theta_levels(i,j,k-1) / pref )**kappa - sgm(i) = a_dqsdt(i,j,k-1)**2 * exner**2 * sl_var & - + a_qs(i,j,k-1)**2 * qw_var & - - 2.0_r_bl * a_qs(i,j,k-1) * a_dqsdt(i,j,k-1) * exner * sl_qw + ! do this for safety, not sure if it's really needed + sgm(i) = sqrt ( max( sgm(i), zero ) ) + ! calculate rhcrit, with appropriate limits + rhcpt(i,j,k-1) = min( max_rhcpt(i,j), max( min_rhcpt(i,j), & + one - root6 * sgm(i) / (a_qs(i,j,k-1) * qsw_arr(i,j,k-1)))) end do !i - if (i_rhcpt == rhcpt_tke_based) then - do i = tdims%i_start, tdims%i_end - ! do this for safety, not sure if it's really needed - sgm(i) = sqrt ( max( sgm(i), zero ) ) - ! calculate rhcrit, with appropriate limits - rhcpt(i,j,k-1) = min( max_rhcpt(i,j), max( min_rhcpt(i,j), & - one - root6 * sgm(i) / (a_qs(i,j,k-1) * qsw_arr(i)))) - end do !i - end if + !$OMP end do NOWAIT + end if - end if ! var_diags_opt + end if ! var_diags_opt - end do !j -!$OMP end do NOWAIT + !end do !j +!!$OMP end do NOWAIT end do !k -!$OMP BARRIER + !$OMP BARRIER if (i_rhcpt == rhcpt_tke_based) then ! just use rhcpt=0.8 above bl_levels -!$OMP do SCHEDULE(STATIC) + !$OMP do SCHEDULE(STATIC) do k = bl_levels-1, tdims%k_end - do j = tdims%j_start, tdims%j_end - do i = tdims%i_start, tdims%i_end - rhcpt(i,j,k) = 0.8_r_bl - end do + !do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end + rhcpt(i,j,k) = 0.8_r_bl end do + !end do end do -!$OMP end do + !$OMP end do end if -!$OMP end PARALLEL + !$OMP end PARALLEL end if !i_rhcpt or variance diagnostics ! ------ liquid potential temperature gradient for bi-modal cloud scheme ------- @@ -3188,42 +3279,42 @@ subroutine bdy_expl2 ( & if ( i_cld_vn == i_cld_bimodal .or. & (i_cld_vn == i_cld_pc2 .and. i_pc2_init_method == pc2init_bimodal) ) then -!$OMP PARALLEL DEFAULT(none) private(k,j,i) & -!$OMP SHARED(tdims,tgrad_bm,bl_levels,kplume,dsldzm,z_tq,zhnl,zh,zhsc, & +!$OMP PARALLEL DEFAULT(none) private(k,i) & +!$OMP SHARED(j,tdims,tgrad_bm,bl_levels,kplume,dsldzm,z_tq,zhnl,zh,zhsc, & !$OMP i_bm_ez_opt,ez_max_bm) !$OMP do SCHEDULE(STATIC) do k = 1, tdims%k_end - do j = tdims%j_start, tdims%j_end - do i = tdims%i_start, tdims%i_end - ! Initialise liquid potential temperature gradient - tgrad_bm(i,j,k) = bm_negative_init - end do + ! do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end + ! Initialise liquid potential temperature gradient + tgrad_bm(i,j,k) = bm_negative_init end do + !end do end do !$OMP end do !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels-1 - do j = tdims%j_start, tdims%j_end - do i = tdims%i_start, tdims%i_end - ! Calculate liquid potential temperature gradient for bimodal cloud - ! scheme entrainment zone identification. Only calculate this gradient - ! if the level is within the top half of the boundary layer and resides - ! below 3000 m (low cloud), or if the level is within ez_max_bm above - ! the top of the mixed layer (defined by zh or zhsc). Else keep the - ! gradient at its initialised zero-value. This change avoids unrealistc - ! free-tropospheric deep entrainment zones. Note that dsldzm is held at - ! k+1 - if ( k-1 > kplume(i,j) .and. k > 2 .and. & - z_tq(i,j,k-1) > one_half*zhnl(i,j) .and. & - (i_bm_ez_opt==i_bm_ez_subcrit .or. & - (z_tq(i,j,k-1) < 3.0e3_r_bl) .or. & - ((z_tq(i,j,k-1)-zh(i,j)) < ez_max_bm) .or. & - ((z_tq(i,j,k-1)-zhsc(i,j)) < ez_max_bm))) then - tgrad_bm(i,j,k-1) = dsldzm(i,j,k) - end if - end do + !do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end + ! Calculate liquid potential temperature gradient for bimodal cloud + ! scheme entrainment zone identification. Only calculate this gradient + ! if the level is within the top half of the boundary layer and resides + ! below 3000 m (low cloud), or if the level is within ez_max_bm above + ! the top of the mixed layer (defined by zh or zhsc). Else keep the + ! gradient at its initialised zero-value. This change avoids unrealistc + ! free-tropospheric deep entrainment zones. Note that dsldzm is held at + ! k+1 + if ( k-1 > kplume(i,j) .and. k > 2 .and. & + z_tq(i,j,k-1) > one_half*zhnl(i,j) .and. & + (i_bm_ez_opt==i_bm_ez_subcrit .or. & + (z_tq(i,j,k-1) < 3.0e3_r_bl) .or. & + ((z_tq(i,j,k-1)-zh(i,j)) < ez_max_bm) .or. & + ((z_tq(i,j,k-1)-zhsc(i,j)) < ez_max_bm))) then + tgrad_bm(i,j,k-1) = dsldzm(i,j,k) + end if end do + !end do end do !$OMP end do !$OMP end PARALLEL @@ -3232,7 +3323,7 @@ subroutine bdy_expl2 ( & !If autotuning is active, decide what to do with the !trial segment size and report the current status. - +call timer('bdy_expl2') if (lhook) call dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) return end subroutine bdy_expl2 From 643f3d9b74990ce74bf05d6d0e36751d9dd90e0e Mon Sep 17 00:00:00 2001 From: MetBenjaminWent Date: Mon, 2 Feb 2026 13:09:31 +0000 Subject: [PATCH 2/5] Add kgo findings for full, revert psyclone workaround for better perfomance --- .../site/meto/groups/groups_lfric_atm.cylc | 5 ++++ ...ca_1T-C48_MG_ex1a_cce_full-debug-32bit.txt | 9 +++++++ ...ca_4T-C48_MG_ex1a_cce_full-debug-32bit.txt | 9 +++++++ .../source/boundary_layer/bdy_expl2.F90 | 26 +++++++++---------- 4 files changed, 36 insertions(+), 13 deletions(-) create mode 100644 rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_noukca_1T-C48_MG_ex1a_cce_full-debug-32bit.txt create mode 100644 rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_noukca_4T-C48_MG_ex1a_cce_full-debug-32bit.txt diff --git a/rose-stem/site/meto/groups/groups_lfric_atm.cylc b/rose-stem/site/meto/groups/groups_lfric_atm.cylc index 44467096..5566752f 100644 --- a/rose-stem/site/meto/groups/groups_lfric_atm.cylc +++ b/rose-stem/site/meto/groups/groups_lfric_atm.cylc @@ -222,6 +222,11 @@ "lfric_atm_clim_gal9_chem_1T-C12_ex1a_cce_fast-debug-32bit", "lfric_atm_clim_gal9_chem_2T-C12_ex1a_cce_fast-debug-32bit", ], + "ex1a_omp_C48_cce_full": [ + "lfric_atm_nwp_gal9_noukca_1T-C48_MG_ex1a_cce_full-debug-32bit", + "lfric_atm_nwp_gal9_noukca_2T-C48_MG_ex1a_cce_full-debug-32bit", + "lfric_atm_nwp_gal9_noukca_4T-C48_MG_ex1a_cce_full-debug-32bit", + ], "ex1a_omp_C48_cce": [ "lfric_atm_nwp_gal9_noukca_1T-C48_MG_ex1a_cce_fast-debug-32bit", "lfric_atm_nwp_gal9_noukca_2T-C48_MG_ex1a_cce_fast-debug-32bit", diff --git a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_noukca_1T-C48_MG_ex1a_cce_full-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_noukca_1T-C48_MG_ex1a_cce_full-debug-32bit.txt new file mode 100644 index 00000000..ecdac22f --- /dev/null +++ b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_noukca_1T-C48_MG_ex1a_cce_full-debug-32bit.txt @@ -0,0 +1,9 @@ +Inner product checksum rho = 48D7FB47 +Inner product checksum theta = 5392A6D6 +Inner product checksum u = 6A97B5F7 +Inner product checksum mr1 = 41CCED78 +Inner product checksum mr2 = 39CD58A2 +Inner product checksum mr3 = 37A6C0C1 +Inner product checksum mr4 = 3970B330 +Inner product checksum mr5 = 0 +Inner product checksum mr6 = 0 diff --git a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_noukca_4T-C48_MG_ex1a_cce_full-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_noukca_4T-C48_MG_ex1a_cce_full-debug-32bit.txt new file mode 100644 index 00000000..4dce9788 --- /dev/null +++ b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_noukca_4T-C48_MG_ex1a_cce_full-debug-32bit.txt @@ -0,0 +1,9 @@ +Inner product checksum rho = 48D7FB59 +Inner product checksum theta = 5392A6CA +Inner product checksum u = 6A97B3E6 +Inner product checksum mr1 = 41CCEE63 +Inner product checksum mr2 = 39D05F5C +Inner product checksum mr3 = 37A7FEC3 +Inner product checksum mr4 = 3970D685 +Inner product checksum mr5 = 0 +Inner product checksum mr6 = 0 diff --git a/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 b/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 index cbf26cb0..0c971a55 100644 --- a/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 +++ b/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 @@ -633,7 +633,7 @@ subroutine bdy_expl2 ( & ! (qw - qsat(Tl))/(1 + Lc/cp dqsat/dT) ! gradient of this is used for estimating ! saturated fraction on rho-levels - qs_tl(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,1:bl_levels) + qs_tl(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end) ! qsat(Tl) on current level ! Vectorised over levels (k) for open-mp ease @@ -781,7 +781,7 @@ subroutine bdy_expl2 ( & real(kind=r_bl) :: & b2, sh, exner, root6, delta_x, var_fac, sl_var, qw_var, sl_qw, & sgm(tdims%i_start:tdims%i_end), & - qsw_arr(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,1:bl_levels-1), & + qsw_arr(tdims%i_start:tdims%i_end), & max_rhcpt(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end), & min_rhcpt(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end) @@ -1128,13 +1128,13 @@ subroutine bdy_expl2 ( & ! Calculate qsat(Tl)... if ( l_mr_physics ) then - call qsat_mix(qs_tl(seg_slice_start:seg_slice_end,:,k), & + call qsat_mix(qs_tl(seg_slice_start:seg_slice_end,:), & tl(seg_slice_start:seg_slice_end,:,k), & p_theta_levels(seg_slice_start:seg_slice_end,:,k), & bl_segment_range, & tdims%j_len ) else - call qsat(qs_tl(seg_slice_start:seg_slice_end,:,k), & + call qsat(qs_tl(seg_slice_start:seg_slice_end,:), & tl(seg_slice_start:seg_slice_end,:,k), & p_theta_levels(seg_slice_start:seg_slice_end,:,k), & bl_segment_range, & @@ -1146,7 +1146,7 @@ subroutine bdy_expl2 ( & ! in order to compare with values of qcl+qcf. !do j = tdims%j_start, tdims%j_end do i = tdims%i_start, tdims%i_end - supersat(i,j,k) = a_qs(i,j,k) * ( qw(i,j,k) - qs_tl(i,j,k) ) + supersat(i,j,k) = a_qs(i,j,k) * ( qw(i,j,k) - qs_tl(i,j) ) end do !i !end do !j end do !k @@ -2984,7 +2984,7 @@ subroutine bdy_expl2 ( & ! level 2 needs special treatment because of the surface !do j = tdims%j_start, tdims%j_end !$OMP PARALLEL DEFAULT(SHARED) & - !$OMP private(k,i,ii,km,kp,sh,exner,sgm,qsw_arr,weight1,weight2,weight3, & + !$OMP private(k,i,ii,km,kp,sh,exner,sgm,weight1,weight2,weight3, & !$OMP var_fac,sl_var,qw_var,sl_qw,delta_x, & !$OMP seg_slice_start, seg_slice_end, bl_segment_range) k = 2 @@ -2996,12 +2996,12 @@ subroutine bdy_expl2 ( & bl_segment_range = (seg_slice_end - seg_slice_start) + 1 if ( l_mr_physics ) then - call qsat_wat_mix(qsw_arr(seg_slice_start:seg_slice_end, j, k-1), & + call qsat_wat_mix(qsw_arr(seg_slice_start:seg_slice_end), & tl(seg_slice_start:seg_slice_end,j,k-1), & p_theta_levels(seg_slice_start:seg_slice_end,j,k-1), & bl_segment_range ) else - call qsat_wat(qsw_arr(seg_slice_start:seg_slice_end, j, k-1), & + call qsat_wat(qsw_arr(seg_slice_start:seg_slice_end), & tl(seg_slice_start:seg_slice_end,j,k-1), & p_theta_levels(seg_slice_start:seg_slice_end,j,k-1), & bl_segment_range ) @@ -3107,7 +3107,7 @@ subroutine bdy_expl2 ( & log( 0.001_r_bl * delta_x ) ) ! full expression rhcpt(i,j,k-1) = min( max_rhcpt(i,j), max( min_rhcpt(i,j), & - one - root6 * sgm(i) / (a_qs(i,j,k-1) * qsw_arr(i,j,k-1)))) + one - root6 * sgm(i) / (a_qs(i,j,k-1) * qsw_arr(i)))) end do !i !$OMP end do end if @@ -3127,13 +3127,13 @@ subroutine bdy_expl2 ( & if ( l_mr_physics ) then !dir$ safe_address - call qsat_wat_mix(qsw_arr(seg_slice_start:seg_slice_end, j, k-1), & + call qsat_wat_mix(qsw_arr(seg_slice_start:seg_slice_end), & tl(seg_slice_start:seg_slice_end,j,k-1), & p_theta_levels(seg_slice_start:seg_slice_end,j,k-1), & bl_segment_range) else !dir$ safe_address - call qsat_wat(qsw_arr(seg_slice_start:seg_slice_end, j, k-1), & + call qsat_wat(qsw_arr(seg_slice_start:seg_slice_end), & tl(seg_slice_start:seg_slice_end,j,k-1), & p_theta_levels(seg_slice_start:seg_slice_end,j,k-1), & bl_segment_range) @@ -3185,7 +3185,7 @@ subroutine bdy_expl2 ( & sgm(i) = sqrt ( max( sgm(i), zero ) ) ! calculate rhcrit, with appropriate limits rhcpt(i,j,k-1) = min( max_rhcpt(i,j), max( min_rhcpt(i,j), & - one - root6 * sgm(i) / (a_qs(i,j,k-1) * qsw_arr(i,j,k-1)))) + one - root6 * sgm(i) / (a_qs(i,j,k-1) * qsw_arr(i)))) end do !i !$OMP end do NOWAIT end if @@ -3246,7 +3246,7 @@ subroutine bdy_expl2 ( & sgm(i) = sqrt ( max( sgm(i), zero ) ) ! calculate rhcrit, with appropriate limits rhcpt(i,j,k-1) = min( max_rhcpt(i,j), max( min_rhcpt(i,j), & - one - root6 * sgm(i) / (a_qs(i,j,k-1) * qsw_arr(i,j,k-1)))) + one - root6 * sgm(i) / (a_qs(i,j,k-1) * qsw_arr(i)))) end do !i !$OMP end do NOWAIT end if From 3b12566b0aab6f430a8d1d7016026ce74e5ceef5 Mon Sep 17 00:00:00 2001 From: MetBenjaminWent Date: Tue, 3 Feb 2026 16:00:33 +0000 Subject: [PATCH 3/5] Remove j loop formatting - omp_blocksize to do --- .../source/boundary_layer/bdy_expl2.F90 | 140 +----------------- 1 file changed, 2 insertions(+), 138 deletions(-) diff --git a/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 b/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 index 0c971a55..0add1d00 100644 --- a/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 +++ b/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 @@ -23,8 +23,6 @@ module bdy_expl2_mod use um_types, only: r_bl -use timer_mod, only: timer - implicit none !Automatic segment size tuning @@ -831,7 +829,6 @@ subroutine bdy_expl2 ( & real(kind=jprb) :: zhook_handle if (lhook) call dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) -call timer('bdy_expl2') ! Set up blocking for OMP sections max_threads = 1 @@ -843,7 +840,6 @@ subroutine bdy_expl2 ( & !----------------------------------------------------------------------- ! checking of nl_bl_levels has been moved to readsize/scm_shell ! so that it is only executed at initialisation -!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL DEFAULT(none) & !$OMP SHARED( j, pdims, dynamic_bl_diag, ntml_save, ntml, bl_levels, & !$OMP weight_1dbl, weight_1dbl_rho, BL_diag, u_s, fb_surf ) & @@ -854,20 +850,17 @@ subroutine bdy_expl2 ( & ntml_save(i,j) = ntml(i,j) end do !$OMP end do NOWAIT - !end do !------------------------------------------------------------------ ! Initialize weighting applied to 1d BL scheme ! (used to blend with 3D Smagorinsky scheme) !------------------------------------------------------------------ - !do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels do i = pdims%i_start, pdims%i_end weight_1dbl(i,j,k) = one weight_1dbl_rho(i,j,k) = one end do - ! end do end do !$OMP end do NOWAIT @@ -875,7 +868,6 @@ subroutine bdy_expl2 ( & ! Set surface scaling diagnostics !----------------------------------------------------------------------- ! Obukhov length - ! do j = pdims%j_start, pdims%j_end if (BL_diag%l_oblen) then !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end @@ -899,18 +891,15 @@ subroutine bdy_expl2 ( & !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end BL_diag%ustar(i,j)=u_s(i,j) - !end do end do !$OMP end do NOWAIT end if ! Surface buoyancy flux - !do j = pdims%j_start, pdims%j_end if (BL_diag%l_wbsurf) then !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end BL_diag%wbsurf(i,j)=fb_surf(i,j) - !end do end do !$OMP end do end if @@ -938,7 +927,6 @@ subroutine bdy_expl2 ( & ! Work should be split accross the threads, this provides a cap. ! In the occurence where the pdims domain len is very small, it will be seleted. ! Otherwise bl_segment_size will be selected. -!do j = pdims%j_start, pdims%j_end pdims_seg_block = min(bl_segment_size, pdims_omp_block, pdims%i_len) !$OMP PARALLEL DEFAULT(SHARED) private(ii, i, k, weight1, weight2, & !$OMP weight3, zpr, dzv, dzu, l, slope, dsldzm_ga, & @@ -949,10 +937,8 @@ subroutine bdy_expl2 ( & grad_t_adj(i,j) = min( max_t_grad, & a_grad_adj * t1_sd(i,j) / zh_prev(i,j) ) end do -!end do !$OMP end do - ! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels do i = pdims%i_start, pdims%i_end @@ -965,7 +951,6 @@ subroutine bdy_expl2 ( & dqwdz(i,j,k) = ( qw(i,j,k) - qw(i,j,k-1) ) & * rdz_charney_grid(i,j,k) end do - !end do end do !$OMP end do @@ -981,19 +966,15 @@ subroutine bdy_expl2 ( & ! extrapolate dbdz itself from level 2, but the sl and qw gradients are ! used in the newer variance calculations and with i_interp_local_cf_dbdz ! so extrapolate them here -! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end dsldz(i,j,1) = dsldz(i,j,2) dsldz_ga(i,j,1) = dsldz_ga(i,j,2) dqwdz(i,j,1) = dqwdz(i,j,2) end do - !end do !$OMP end do else ! l_use_surf_in_ri = true -!! $OMP do SCHEDULE(STATIC) - !do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do ii = pdims%i_start, pdims%i_end, pdims_seg_block seg_slice_start = ii @@ -1029,11 +1010,8 @@ subroutine bdy_expl2 ( & end if ! l_noice_in_turb end do ! ii !$OMP end do NOWAIT - !end do ! j -!! $OMP end do NOWAIT k=1 - ! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end dsldz(i,j,k) = ( tl(i,j,k) - tstar(i,j) ) & @@ -1045,7 +1023,6 @@ subroutine bdy_expl2 ( & else dqwdz(i,j,k) = dqwdz(i,j,2) ! extrapolate qw if mainly land end if - !end do end do !$OMP end do @@ -1060,7 +1037,6 @@ subroutine bdy_expl2 ( & ! Interpolate gradients to theta-levels, then ! calculate `buoyancy' gradient, DBDZ, on theta-levels ! NOTE: DBDZ(K) is on theta-level K-1 - ! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels do i = pdims%i_start, pdims%i_end @@ -1085,7 +1061,6 @@ subroutine bdy_expl2 ( & dsldzm(i,j,k) = dsldzm(i,j,k) / weight1 dqwdzm(i,j,k) = dqwdzm(i,j,k) / weight1 end do - !end do end do !$OMP end do @@ -1093,7 +1068,6 @@ subroutine bdy_expl2 ( & ! (ie instead of using centred value that uses surface parameters) if ( .not. l_use_surf_in_ri ) then k = 2 - ! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end dbdz(i,j,k) = g*( bt_gb(i,j,k-1)*dsldz(i,j,k) + & @@ -1101,7 +1075,6 @@ subroutine bdy_expl2 ( & dbdz_ga(i,j,k) = g*( bt_gb(i,j,k-1)*dsldz_ga(i,j,k) + & bq_gb(i,j,k-1)*dqwdz(i,j,k) ) end do - ! end do !$OMP end do end if @@ -1144,20 +1117,16 @@ subroutine bdy_expl2 ( & ! ...then subtract from qw to get supersaturation, and multiply by ! 1/(1 + Lc/cp dqsat/dT) ! in order to compare with values of qcl+qcf. - !do j = tdims%j_start, tdims%j_end do i = tdims%i_start, tdims%i_end supersat(i,j,k) = a_qs(i,j,k) * ( qw(i,j,k) - qs_tl(i,j) ) end do !i - !end do !j end do !k !$OMP end do ! Calculate dbdz on rho-levels, so between th-levels k-1 and k !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - !do j = tdims%j_start, tdims%j_end do i = tdims%i_start, tdims%i_end - ! Interpolate the buoyancy coefficients onto rho-levels... ! Split the interpolation into weighted contributions from: ! a) Area that is cloudy at k-1 and at k @@ -1221,13 +1190,11 @@ subroutine bdy_expl2 ( & dbdz_ga_rh(i,j,k) = g * ( bt_rh * dsldz_ga(i,j,k) & + bq_rh * dqwdz(i,j,k) ) - !end do end do end do !$OMP end do NOWAIT k = 1 - !do j = tdims%j_start, tdims%j_end !$OMP do SCHEDULE(STATIC) do i = tdims%i_start, tdims%i_end ! At surface, just use buoyancy coefficients from the first theta-level @@ -1235,11 +1202,9 @@ subroutine bdy_expl2 ( & + bq_gb(i,j,k) * dqwdz(i,j,k) ) dbdz_ga_rh(i,j,k) = g * ( bt_gb(i,j,k) * dsldz_ga(i,j,k) & + bq_gb(i,j,k) * dqwdz(i,j,k) ) - !end do end do !$OMP end do - !do j = tdims%j_start, tdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels do i = tdims%i_start, tdims%i_end @@ -1261,7 +1226,6 @@ subroutine bdy_expl2 ( & + weight1 * dqwdz(i,j,k) end do - !end do end do !$OMP end do @@ -1274,7 +1238,6 @@ subroutine bdy_expl2 ( & !-------------------------------------------------- if (.not. l_subfilter_vert) then -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels do i = pdims%i_start, pdims%i_end @@ -1282,7 +1245,6 @@ subroutine bdy_expl2 ( & dzv = v_p(i,j,k) - v_p(i,j,k-1) dvdzm(i,j,k) = max( 1.0e-12_r_bl, & sqrt(dzu*dzu + dzv*dzv) * rdz(i,j,k) ) - !end do end do end do !$OMP end do @@ -1292,7 +1254,6 @@ subroutine bdy_expl2 ( & if ( model_type == mt_single_column ) then ! visc_m, visc_h need to be the 1D shear(k) on theta-level(k) -! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels do i = pdims%i_start, pdims%i_end @@ -1302,7 +1263,6 @@ subroutine bdy_expl2 ( & sqrt(dzu*dzu + dzv*dzv) * rdz(i,j,k) ) visc_m(i,j,k-1) = dvdzm(i,j,k) visc_h(i,j,k-1) = dvdzm(i,j,k) - !end do end do end do !$OMP end do @@ -1310,12 +1270,10 @@ subroutine bdy_expl2 ( & else ! On entry, visc_m is 3D shear(k) on theta-level(k) -! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels do i = pdims%i_start, pdims%i_end dvdzm(i,j,k) = max( 1.0e-12_r_bl, visc_m(i,j,k-1) ) - ! end do end do end do !$OMP end do @@ -1327,13 +1285,11 @@ subroutine bdy_expl2 ( & if (smag_l_calc == smag_l_calc_use_geo) then ! use 3d grid geometric mean -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels do i = pdims%i_start, pdims%i_end rmlmax2(i,j,k) = ( delta_smag(i,j) * delta_smag(i,j) * & dzl_charney(i,j,k) )**one_third - !end do end do end do !$OMP end do @@ -1342,10 +1298,8 @@ subroutine bdy_expl2 ( & ! only use horizontal grid !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels - !do j = pdims%j_start, pdims%j_end do i = pdims%i_start, pdims%i_end rmlmax2(i,j,k) = delta_smag(i,j) - !end do end do end do !$OMP end do @@ -1354,38 +1308,32 @@ subroutine bdy_expl2 ( & if ( l_rp2 .and. i_rp_scheme == i_rp2b ) then -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels do i = pdims%i_start, pdims%i_end rmlmax2(i,j,k) = ( cs_rp(rp_idx) * rmlmax2(i,j,k) )**2 - !end do end do end do !$OMP end do else -! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels do i = pdims%i_start, pdims%i_end rmlmax2(i,j,k) = ( mix_factor * rmlmax2(i,j,k) )**2 - !end do end do end do !$OMP end do end if -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels do i = pdims%i_start, pdims%i_end rneutml_sq(i,j,k) = one / ( & one/( vkman*(z_tq(i,j,k) + z0m_eff_gb(i,j)) )**2 & + one/rmlmax2(i,j,k) ) - !end do end do end do !$OMP end do @@ -1396,11 +1344,9 @@ subroutine bdy_expl2 ( & !----------------------------------------------------------------------- ! Set-up 2D array for standard deviation of subgrid orography. !----------------------------------------------------------------------- -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end sigma_h(i,j) = zero -!end do end do !$OMP end do @@ -1417,7 +1363,6 @@ subroutine bdy_expl2 ( & if (sg_orog_mixing == sg_shear .or. & sg_orog_mixing == sg_shear_enh_lambda) then -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels do i = pdims%i_start, pdims%i_end @@ -1439,55 +1384,46 @@ subroutine bdy_expl2 ( & BL_diag%dvdzm(i,j,1)=weight1*slope*t_drain*dbdz(i,j,k) end if - !end do end do end do !$OMP end do end if ! sg_orog_mixing -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels do i = pdims%i_start, pdims%i_end ri(i,j,k) = dbdz(i,j,k) / ( dvdzm(i,j,k)*dvdzm(i,j,k) ) ri(i,j,k) = max(min(ri(i,j,k),max_ri),-max_ri) ri_ga(i,j,k) = dbdz_ga(i,j,k) / ( dvdzm(i,j,k)*dvdzm(i,j,k) ) - ! end do end do end do !$OMP end do -!do j = pdims%j_start, pdims%j_end if (BL_diag%l_gradrich) then !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels do i = pdims%i_start, pdims%i_end BL_diag%gradrich(i,j,k)=ri(i,j,k) - !end do end do end do !$OMP end do end if -!do j = pdims%j_start, pdims%j_end if (BL_diag%l_dbdz) then !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels do i = pdims%i_start, pdims%i_end BL_diag%dbdz(i,j,k)=dbdz(i,j,k) - !end do end do end do !$OMP end do end if -!do j = pdims%j_start, pdims%j_end if (BL_diag%l_dvdzm) then !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels do i = pdims%i_start, pdims%i_end BL_diag%dvdzm(i,j,k)=dvdzm(i,j,k) - !end do end do end do !$OMP end do @@ -1516,7 +1452,6 @@ subroutine bdy_expl2 ( & ! Orographic stress diagnostics !------------------------------------------------------------------ if (BL_diag%l_ostressx) then -!do j = tdims%j_start, tdims%j_end !$OMP PARALLEL do & !$OMP SCHEDULE(STATIC) & !$OMP DEFAULT(none) & @@ -1525,7 +1460,6 @@ subroutine bdy_expl2 ( & do k = 1, bl_levels do i = tdims%i_start, tdims%i_end BL_diag%ostressx(i,j,k)=tau_fd_x(i,j,k) - !end do end do end do !$OMP end PARALLEL do @@ -1541,7 +1475,6 @@ subroutine bdy_expl2 ( & do k = 1, bl_levels do i = tdims%i_start, tdims%i_end BL_diag%ostressy(i,j,k)=tau_fd_y(i,j,k) - !end do end do end do !$OMP end PARALLEL do @@ -1574,7 +1507,6 @@ subroutine bdy_expl2 ( & zh(i,j) = z_uv(i,j,2) dzh(i,j) = rmdi end if -!end do end do !$OMP end do @@ -1584,7 +1516,6 @@ subroutine bdy_expl2 ( & ! Original version - causes spuriously deep boundary layers if ! BL_LEVELS is >> 3km -! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end @@ -1598,13 +1529,11 @@ subroutine bdy_expl2 ( & end if end if - !end do end do !$OMP end do else if (idyndiag == DynDiag_ZL_corrn) then - ! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end @@ -1623,12 +1552,10 @@ subroutine bdy_expl2 ( & zh(i,j) = z_uv(i,j,ntml(i,j)+1) end if end if - !end do end do !$OMP end do else if (idyndiag == DynDiag_ZL_CuOnly) then - ! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end @@ -1656,7 +1583,6 @@ subroutine bdy_expl2 ( & zh(i,j) = z_uv(i,j,ntml(i,j)+1) end if end if - !end do end do !$OMP end do @@ -1676,7 +1602,6 @@ subroutine bdy_expl2 ( & ! Loop over levels to find Ri > RiCrit_sharp (=0.25) to find ! level to which Ri really is close to neutral !--------------------------------------------------------------- -! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do ii = pdims%i_start, pdims%i_end, pdims_omp_block do k = 2, bl_levels @@ -1703,7 +1628,6 @@ subroutine bdy_expl2 ( & zhloc_depth_fac = zhloc_depth_fac_rp(rp_idx) end if -! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end @@ -1732,7 +1656,6 @@ subroutine bdy_expl2 ( & end if end if ! Cu over sea - !end do end do !$OMP end do @@ -1759,7 +1682,6 @@ subroutine bdy_expl2 ( & end if end if ! Cu over sea - !end do end do !$OMP end do @@ -1789,7 +1711,6 @@ subroutine bdy_expl2 ( & do i = pdims%i_start, pdims%i_end ntml_nl(i,j) = ntml(i,j) end do -!end do !$OMP end PARALLEL do if (nl_bl_levels < bl_levels) then @@ -1797,7 +1718,6 @@ subroutine bdy_expl2 ( & zmaxb_for_dsc = 1.0e10_r_bl zmaxt_for_dsc = zmaxb_for_dsc -! do j = pdims%j_start, pdims%j_end !$OMP PARALLEL do & !$OMP SCHEDULE(STATIC) & !$OMP DEFAULT(none) & @@ -1808,7 +1728,6 @@ subroutine bdy_expl2 ( & ntml_nl(i,j) = nl_bl_levels-1 zh(i,j) = z_uv(i,j,ntml_nl(i,j)+1) end if - !end do end do !$OMP end PARALLEL do else @@ -1819,7 +1738,6 @@ subroutine bdy_expl2 ( & ! Initialise non-local K and fluxes to zero; necessary for levels ! above NL_BL_LEVELS !----------------------------------------------------------------------- - ! do j = pdims%j_start, pdims%j_end !$OMP PARALLEL DEFAULT(none) & !$OMP SHARED( j, bl_levels, pdims, ftl, fqw, rhokmz, rhokhz, rhokm_top, & !$OMP rhokh_top, tke_nl, f_ngstress, rhof2, rhofsc, ft_nt, fq_nt, & @@ -1843,7 +1761,6 @@ subroutine bdy_expl2 ( & ! Only calculated in KMKHZ9C ft_nt(i,j,k) = zero fq_nt(i,j,k) = zero - !end do end do end do !$OMP end do NOWAIT @@ -1861,7 +1778,6 @@ subroutine bdy_expl2 ( & fq_nt_dscb(i,j) = zero zhnl(i,j) = zh(i,j) ! initialise non-local PBL depth ! to that diagnosed in conv_diag -!end do end do !$OMP end do !$OMP end PARALLEL @@ -1890,7 +1806,6 @@ subroutine bdy_expl2 ( & wtrac_bl(i_wt)%fq_nt_dscb(i,j) = zero wtrac_bl(i_wt)%totqf_zh(i,j) = zero wtrac_bl(i_wt)%totqf_zhsc(i,j) = zero - !end do end do !$OMP end do !$OMP end PARALLEL @@ -2001,7 +1916,6 @@ subroutine bdy_expl2 ( & !$OMP end PARALLEL do if (blending_option == blend_cth_shcu_only) then -! do j = pdims%j_start, pdims%j_end ! only going to use the parcel top as the length scale in blending ! if the convection is shallow, where we define shallow convection here ! as the resolved cloud top (from cf_bluk) being below shallow_cu_maxtop @@ -2014,7 +1928,6 @@ subroutine bdy_expl2 ( & do i = pdims%i_start, pdims%i_end cloud_base_found(i,j) = .false. cloud_top_found(i,j) = .false. - !end do end do !$OMP end do @@ -2034,7 +1947,6 @@ subroutine bdy_expl2 ( & cloud_top_found(i,j) = .true. l_shallow_cth(i,j) = z_uv(i,j,k+1) < shallow_cu_maxtop end if - !end do end do !$OMP end do end do @@ -2062,7 +1974,6 @@ subroutine bdy_expl2 ( & !----------------------------------------------------------------------- ! interpolate rhokh_th to rho levels 2 to bl_levels -! do j = pdims%j_start, pdims%j_end !$OMP PARALLEL DEFAULT(SHARED) & !$OMP private(i, k, weight1, weight2, weight3, lambdah, z_scale, & !$OMP vkz, f_log) @@ -2166,7 +2077,6 @@ subroutine bdy_expl2 ( & if (l_mr_physics) & rhokh(i,j,k) = rho_mix(i,j,k) * rhokh(i,j,k) - !end do end do end do !$OMP end do NOWAIT @@ -2199,7 +2109,6 @@ subroutine bdy_expl2 ( & dsc(i,j) = .false. coupled(i,j) = .false. end if -!end do end do !$OMP end do !$OMP end PARALLEL @@ -2225,7 +2134,6 @@ subroutine bdy_expl2 ( & do i = pdims%i_start, pdims%i_end BL_diag%weight1d(i,j,k)=weight_1dbl(i,j,k) end do - !end do end do !$OMP end PARALLEL do end if @@ -2242,7 +2150,6 @@ subroutine bdy_expl2 ( & do i = pdims%i_start, pdims%i_end BL_diag%rhokm(i,j,k)=rhokm(i,j,k) end do - !end do end do !$OMP end PARALLEL do end if @@ -2289,8 +2196,6 @@ subroutine bdy_expl2 ( & BL_diag%tke(i,j,k)= min( max_tke, & ( rhokm(i,j,k) / rho_wet_tq(i,j,k-1) )*( & recip_time_cbl + recip_time_sbl ) ) - - !end do end do end do !$OMP end PARALLEL do @@ -2323,7 +2228,6 @@ subroutine bdy_expl2 ( & ! equivalent limiting is done on Km. end do - !end do end do !$OMP end PARALLEL do @@ -2350,7 +2254,6 @@ subroutine bdy_expl2 ( & else mix_len_bm(i,j,k) = bm_tiny end if - !end do end do end do !$OMP end do @@ -2381,7 +2284,6 @@ subroutine bdy_expl2 ( & + weight1 * max( 0.5 * mix_len_bm(i,j,k-1) & - mix_len_bm(i,j,k), 0.0 ) end if - !end do end do !$OMP end do NOWAIT ! Set values above bl_levels @@ -2391,7 +2293,6 @@ subroutine bdy_expl2 ( & do i = pdims%i_start, pdims%i_end mix_len_bm(i,j,k) = bm_tiny end do - !end do end do !$OMP end do NOWAIT !$OMP end PARALLEL @@ -2410,12 +2311,10 @@ subroutine bdy_expl2 ( & !$OMP PARALLEL DEFAULT(none) private(k,i) & !$OMP SHARED(j,tdims, bl_levels,pdims,BL_diag,bl_w_var,var_diags_opt) -!do j = tdims%j_start, tdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, tdims%k_end+1 do i = tdims%i_start, tdims%i_end bl_w_var(i,j,k) = 1.0e-12_r_bl - !end do end do end do !$OMP end do @@ -2429,7 +2328,6 @@ subroutine bdy_expl2 ( & if (BL_diag%tke(i,j,k) > 1.0e-12_r_bl) then bl_w_var(i,j,k) = BL_diag%tke(i,j,k) end if - !end do end do end do !$OMP end do @@ -2442,7 +2340,6 @@ subroutine bdy_expl2 ( & if (BL_diag%tke(i,j,k) > 1.0e-12_r_bl) then bl_w_var(i,j,k) = two_thirds * BL_diag%tke(i,j,k) end if - !end do end do end do !$OMP end do @@ -2462,7 +2359,6 @@ subroutine bdy_expl2 ( & kmax = maxloc(rhokmz(i,j,:), dim=1) do k = 2, kmax BL_diag%tke(i,j,k) = BL_diag%tke(i,j,kmax) - !end do end do end do !$OMP end PARALLEL do @@ -2483,7 +2379,6 @@ subroutine bdy_expl2 ( & do k = 2, kp-1 BL_diag%tke(i,j,k) = max( tke_loc(i,j,k), tke_nl(i,j,kp) ) end do - !end do end do !$OMP end PARALLEL do end if ! test on var_diags_opt @@ -2497,7 +2392,6 @@ subroutine bdy_expl2 ( & !$OMP PARALLEL DEFAULT(SHARED) private(i, k, weight1, weight2, & !$OMP weight3) -! do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels-1 do i = pdims%i_start, pdims%i_end @@ -2517,7 +2411,6 @@ subroutine bdy_expl2 ( & else visc_m(i,j,k) = rhokm(i,j,k+1)/rho_wet_tq(i,j,k) end if - !end do end do end do !$OMP end do NOWAIT @@ -2538,7 +2431,6 @@ subroutine bdy_expl2 ( & / weight1 end if end do - !end do else !do j = pdims%j_start, pdims%j_end do i = pdims%i_start, pdims%i_end @@ -2549,7 +2441,6 @@ subroutine bdy_expl2 ( & visc_h(i,j,k) = rhokh_th(i,j,k+1)/rho_wet_tq(i,j,k) end if end do - !end do end if end do !$OMP end do @@ -2569,7 +2460,6 @@ subroutine bdy_expl2 ( & do i = pdims%i_start, pdims%i_end visc_m(i,j,k) = visc_m(i,j,k)*rneutml_sq(i,j,k) visc_h(i,j,k) = visc_h(i,j,k)*rneutml_sq(i,j,k) - !end do end do !do j = pdims%j_start, pdims%j_end if ( k < bl_levels ) then @@ -2583,7 +2473,6 @@ subroutine bdy_expl2 ( & ! that actually needs it)? visc_h(i,j,k) = min(visc_h(i,j,k),max_diff(i,j)) visc_m(i,j,k) = min(visc_m(i,j,k),max_diff(i,j)) - !end do end do end if end do @@ -2616,7 +2505,6 @@ subroutine bdy_expl2 ( & visc_h_rho(i,j,k) = ( weight3/weight1 ) * visc_h(i,j,k) & + ( weight2/weight1 ) * visc_h(i,j,k-1) end if - !end do end do ! Overwrite the diffusion coefficients from the BL scheme @@ -2629,14 +2517,12 @@ subroutine bdy_expl2 ( & rhokm(i,j,k) = visc_m(i,j,k-1)*rho_wet_tq(i,j,k-1) rhokh(i,j,k) = visc_h_rho(i,j,k)*rho_mix(i,j,k) end do - !end do else !do j = pdims%j_start, pdims%j_end do i = pdims%i_start, pdims%i_end rhokm(i,j,k) = zero rhokh(i,j,k) = zero end do - !end do end if end do !$OMP end PARALLEL do @@ -2681,7 +2567,6 @@ subroutine bdy_expl2 ( & if ( ntml(i,j) > ntml_nl(i,j) ) then bl_diag%zht(i,j) = max( bl_diag%zht(i,j), zh_local(i,j) ) end if - !end do end do !$OMP end do NOWAIT end if @@ -2714,7 +2599,6 @@ subroutine bdy_expl2 ( & bl_type_6(i,j) = zero bl_type_7(i,j) = zero end do -!end do !$OMP end do NOWAIT !do j = pdims%j_start, pdims%j_end @@ -2786,7 +2670,6 @@ subroutine bdy_expl2 ( & do i = pdims%i_start, pdims%i_end f_ngstress(i,j,k) = weight_1dbl(i,j,k) * f_ngstress(i,j,k) end do - !end do end do !$OMP end PARALLEL do @@ -2812,7 +2695,6 @@ subroutine bdy_expl2 ( & vw0(i,j) = -rhokm(i,j,1) * & ( v_p(i,j,1) - v_0_px(i,j) ) end do -!end do !$OMP end do NOWAIT if (formdrag == explicit_stress) then !do j = pdims%j_start, pdims%j_end @@ -2843,7 +2725,6 @@ subroutine bdy_expl2 ( & else ntml(i,j) = max( ntml_nl(i,j) , ntdsc(i,j) ) end if - !end do end do !$OMP end do NOWAIT @@ -2854,7 +2735,6 @@ subroutine bdy_expl2 ( & do i = pdims%i_start, pdims%i_end ntml(i,j) = ntml_local(i,j) end do - !end do !$OMP end do NOWAIT end if @@ -2888,7 +2768,6 @@ subroutine bdy_expl2 ( & BL_diag%wstar(i,j)= (zh(i,j)*fb_surf(i,j))**one_third end if end do -!end do !$OMP end do if (l_param_conv) then @@ -2924,7 +2803,6 @@ subroutine bdy_expl2 ( & if ( .not. cumulus(i,j) ) l_shallow(i,j) = .false. if ( cumulus(i,j) ) ntml(i,j) = ntml_save(i,j) end do - !end do !$OMP end do NOWAIT end if ! (l_param_conv) @@ -2942,7 +2820,6 @@ subroutine bdy_expl2 ( & shallowc(i,j) = zero end if end do -!end do !$OMP end do !$OMP end PARALLEL @@ -2982,7 +2859,6 @@ subroutine bdy_expl2 ( & ! Otherwise bl_segment_size will be selected. tdims_seg_block = min(bl_segment_size, tdims_omp_block, tdims%i_len) ! level 2 needs special treatment because of the surface - !do j = tdims%j_start, tdims%j_end !$OMP PARALLEL DEFAULT(SHARED) & !$OMP private(k,i,ii,km,kp,sh,exner,sgm,weight1,weight2,weight3, & !$OMP var_fac,sl_var,qw_var,sl_qw,delta_x, & @@ -3111,14 +2987,10 @@ subroutine bdy_expl2 ( & end do !i !$OMP end do end if - !end do !j -!!$OMP end do ! remaining bl levels use proper interpolation do k = 3, bl_levels-1 -!! $OMP do SCHEDULE(STATIC) - !do j = tdims%j_start, tdims%j_end !$OMP do SCHEDULE(STATIC) do ii = tdims%i_start, tdims%i_end, tdims_seg_block seg_slice_start = ii @@ -3253,20 +3125,16 @@ subroutine bdy_expl2 ( & end if ! var_diags_opt - !end do !j -!!$OMP end do NOWAIT end do !k !$OMP BARRIER if (i_rhcpt == rhcpt_tke_based) then ! just use rhcpt=0.8 above bl_levels !$OMP do SCHEDULE(STATIC) - do k = bl_levels-1, tdims%k_end - !do j = tdims%j_start, tdims%j_end + do k = bl_levels-1, tdims%k_end do i = tdims%i_start, tdims%i_end rhcpt(i,j,k) = 0.8_r_bl - end do - !end do + end do end do !$OMP end do end if @@ -3289,13 +3157,11 @@ subroutine bdy_expl2 ( & ! Initialise liquid potential temperature gradient tgrad_bm(i,j,k) = bm_negative_init end do - !end do end do !$OMP end do !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels-1 - !do j = tdims%j_start, tdims%j_end do i = tdims%i_start, tdims%i_end ! Calculate liquid potential temperature gradient for bimodal cloud ! scheme entrainment zone identification. Only calculate this gradient @@ -3314,7 +3180,6 @@ subroutine bdy_expl2 ( & tgrad_bm(i,j,k-1) = dsldzm(i,j,k) end if end do - !end do end do !$OMP end do !$OMP end PARALLEL @@ -3323,7 +3188,6 @@ subroutine bdy_expl2 ( & !If autotuning is active, decide what to do with the !trial segment size and report the current status. -call timer('bdy_expl2') if (lhook) call dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) return end subroutine bdy_expl2 From 5f76d3e6c2b2ffa3d9ab1b82ee2b8d5c12ef2a6f Mon Sep 17 00:00:00 2001 From: MetBenjaminWent Date: Tue, 3 Feb 2026 16:05:20 +0000 Subject: [PATCH 4/5] remove addtional j loops --- .../source/boundary_layer/bdy_expl2.F90 | 52 ++----------------- 1 file changed, 4 insertions(+), 48 deletions(-) diff --git a/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 b/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 index 0add1d00..2cf7e12c 100644 --- a/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 +++ b/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 @@ -798,9 +798,9 @@ subroutine bdy_expl2 ( & integer :: & - pdims_omp_block, pdims_seg_block, & + pdims_omp_block, pdims_seg_block, & ! for pdims open-mp blocking - tdims_omp_block, tdims_seg_block, & + tdims_omp_block, tdims_seg_block, & ! for tdims open mp blocking ii, & ! for indexing over open-mp block @@ -886,7 +886,6 @@ subroutine bdy_expl2 ( & end if ! Ustar - !do j = pdims%j_start, pdims%j_end if (BL_diag%l_ustar) then !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end @@ -1493,7 +1492,6 @@ subroutine bdy_expl2 ( & ! of importance mainly at sea-points, to avoid complications ! with coastal tiling, the scheme operates only at points ! where the land fraction is below 0.5. -!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL DEFAULT(SHARED) private(i, k, ii, ntop, z_scale) !$OMP do SCHEDULE(STATIC) @@ -1591,7 +1589,6 @@ subroutine bdy_expl2 ( & ! As DynDiag_ZL_CuOnly but also allow ZH(Ri) to overrule the ! Cumulus diagnosis !------------------------------------------------------------ - !do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end topbl(i,j) = .false. @@ -1660,8 +1657,6 @@ subroutine bdy_expl2 ( & !$OMP end do else ! old version: not good to set ntml=ntop for small zh/L - -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end @@ -1702,7 +1697,6 @@ subroutine bdy_expl2 ( & ! Set NTML_NL to NTML as passed in from initial diagnosis routine !----------------------------------------------------------------------- -!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL do & !$OMP SCHEDULE(STATIC) & !$OMP DEFAULT(none) & @@ -1765,7 +1759,6 @@ subroutine bdy_expl2 ( & end do !$OMP end do NOWAIT -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end ft_nt(i,j,1) = zero @@ -1786,7 +1779,6 @@ subroutine bdy_expl2 ( & if (l_wtrac) then do i_wt = 1, n_wtrac -!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL DEFAULT(none) & !$OMP SHARED( j, i_wt, bl_levels, pdims, wtrac_bl) private( i, k) !$OMP do SCHEDULE(STATIC) @@ -1799,7 +1791,6 @@ subroutine bdy_expl2 ( & end do !$OMP end do -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end wtrac_bl(i_wt)%fq_nt(i,j,1) = zero @@ -1850,7 +1841,6 @@ subroutine bdy_expl2 ( & ! Set all variables from the non-local scheme to zero or "off" ! - reset all fluxes and K's arising from the non-local scheme !------------------------------------------------------------- -!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL DEFAULT(none) private(i,ient) & !$OMP SHARED(j,pdims,unstable,fb_surf,cumulus,l_shallow,sml_disc_inv,ntpar, & !$OMP ntml_nl,zhnl,grad_t_adj,grad_q_adj,dsc,dsc_disc_inv,ntdsc,nbdsc, & @@ -1885,7 +1875,6 @@ subroutine bdy_expl2 ( & ! end do !$OMP end do NOWAIT do ient = 1, 3 -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end @@ -1902,7 +1891,6 @@ subroutine bdy_expl2 ( & !$OMP end PARALLEL end if ! test on NON_LOCAL_BL -!do j = pdims%j_start, pdims%j_end ! default values !$OMP PARALLEL do & !$OMP SCHEDULE(STATIC) & @@ -1931,7 +1919,6 @@ subroutine bdy_expl2 ( & end do !$OMP end do -!do j = pdims%j_start, pdims%j_end do k = 3, bl_levels-1 !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end @@ -2081,7 +2068,6 @@ subroutine bdy_expl2 ( & end do !$OMP end do NOWAIT -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end !---------------------------------------------------------- @@ -2123,7 +2109,6 @@ subroutine bdy_expl2 ( & rhokm,rhokh,rhokmz,rhokhz,rhokm_top,rhokh_top,tke_loc & ) -!do j = pdims%j_start, pdims%j_end if (BL_diag%l_weight1d) then !$OMP PARALLEL do & !$OMP SCHEDULE(STATIC) & @@ -2139,7 +2124,6 @@ subroutine bdy_expl2 ( & end if -!do j = pdims%j_start, pdims%j_end if (BL_diag%l_rhokm) then !$OMP PARALLEL do & !$OMP SCHEDULE(STATIC) & @@ -2154,7 +2138,6 @@ subroutine bdy_expl2 ( & !$OMP end PARALLEL do end if -!do j = pdims%j_start, pdims%j_end if (BL_diag%l_rhokh) then !$OMP PARALLEL do & !$OMP SCHEDULE(STATIC) & @@ -2181,7 +2164,6 @@ subroutine bdy_expl2 ( & ! BL_diag%tke currently contains (the reciprocal of) the non-local ! (SML and DSC) mixed layer timescale, calculated in excf_nl -!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(none) & !$OMP SHARED(j, BL_diag, rho_wet_tq, rhokm, dbdz, bl_levels, pdims) & !$OMP private(i, k, recip_time_sbl, recip_time_cbl) @@ -2203,7 +2185,6 @@ subroutine bdy_expl2 ( & else if (var_diags_opt == split_tke_and_inv) then ! Combine the, separately calculated, local and non-local TKE diagnostics -!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(none) & !$OMP SHARED(j, BL_diag, tke_nl, tke_loc, rho_wet_tq, weight_1dbl, & !$OMP tke_diag_fac, bl_levels, pdims) & @@ -2234,7 +2215,6 @@ subroutine bdy_expl2 ( & end if ! var_diags_opt -!do j = pdims%j_start, pdims%j_end if ( i_bm_ez_opt == i_bm_ez_entpar ) then ! Calculate mixing-length to pass to bimodal cloud scheme, ! using Km and TKE. @@ -2263,7 +2243,6 @@ subroutine bdy_expl2 ( & ! increment ntml/ntdsc due to the inversion rise or fall over the ! current timestep, but doesn't increment zhnl/zhsc to keep it consistent. ! We need to check for this to get sensible interpolation weights... -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end if ( sml_disc_inv(i,j) > 0 .and. ntml_nl(i,j) > 1 ) then @@ -2287,7 +2266,6 @@ subroutine bdy_expl2 ( & end do !$OMP end do NOWAIT ! Set values above bl_levels -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = bl_levels, tdims%k_end do i = pdims%i_start, pdims%i_end @@ -2321,7 +2299,6 @@ subroutine bdy_expl2 ( & if (var_diags_opt == original_vars) then ! Original version: w_var = TKE -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels do i = pdims%i_start, pdims%i_end @@ -2333,7 +2310,6 @@ subroutine bdy_expl2 ( & !$OMP end do else if (var_diags_opt == split_tke_and_inv) then ! Improved version: w_var = 2/3 TKE -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels do i = pdims%i_start, pdims%i_end @@ -2351,8 +2327,7 @@ subroutine bdy_expl2 ( & if (var_diags_opt == original_vars) then ! at this point, BL_diag%tke really contains 1.5*sigma_w^2. To make it look ! a bit more like TKE near the surface, we will keep it constant below - ! the max of rhokm_surf - !do j = pdims%j_start, pdims%j_end + ! the max of rhokm_surf !$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(none) & !$OMP SHARED( j, pdims, rhokmz, BL_diag ) private(i, k, kmax ) do i = pdims%i_start, pdims%i_end @@ -2367,7 +2342,6 @@ subroutine bdy_expl2 ( & ! a bit more like TKE near the surface, we will keep it constant below ! the max of rhokm_surf (here we find the first local maximum in case there ! is a larger rhokmz in a resolved inversion -!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(none) & !$OMP SHARED( j, pdims, bl_levels, rhokmz, BL_diag, tke_loc, tke_nl ) & !$OMP private(i, k, kp) @@ -2417,7 +2391,6 @@ subroutine bdy_expl2 ( & !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels-1 if ( k > 1 ) then - !do j = pdims%j_start, pdims%j_end do i = pdims%i_start, pdims%i_end if (.not. ( blending_option == blend_except_cu .and. & cumulus(i,j) .and. ntdsc(i,j) == 0)) then @@ -2432,7 +2405,6 @@ subroutine bdy_expl2 ( & end if end do else - !do j = pdims%j_start, pdims%j_end do i = pdims%i_start, pdims%i_end if (.not. ( blending_option == blend_except_cu .and. & cumulus(i,j) .and. ntdsc(i,j) == 0)) then @@ -2450,7 +2422,6 @@ subroutine bdy_expl2 ( & ! visc_m,h on in are just S and visc_m,h(k) are co-located with w(k) -!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL do SCHEDULE(STATIC) & !$OMP DEFAULT(none) & !$OMP private(k,i) & @@ -2461,7 +2432,6 @@ subroutine bdy_expl2 ( & visc_m(i,j,k) = visc_m(i,j,k)*rneutml_sq(i,j,k) visc_h(i,j,k) = visc_h(i,j,k)*rneutml_sq(i,j,k) end do -!do j = pdims%j_start, pdims%j_end if ( k < bl_levels ) then do i = pdims%i_start, pdims%i_end @@ -2486,8 +2456,7 @@ subroutine bdy_expl2 ( & allocate (visc_h_rho(pdims%i_start:pdims%i_end, & pdims%j_start:pdims%j_end, bl_levels)) - -!do j = pdims%j_start, pdims%j_end + !$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(none) & !$OMP SHARED(j,bl_levels, pdims, r_theta_levels, r_rho_levels, visc_h_rho, & !$OMP visc_h, turb_startlev_vert, turb_endlev_vert, rhokm, visc_m, & @@ -2512,13 +2481,11 @@ subroutine bdy_expl2 ( & if (k >= turb_startlev_vert .and. & k <= turb_endlev_vert) then - !do j = pdims%j_start, pdims%j_end do i = pdims%i_start, pdims%i_end rhokm(i,j,k) = visc_m(i,j,k-1)*rho_wet_tq(i,j,k-1) rhokh(i,j,k) = visc_h_rho(i,j,k)*rho_mix(i,j,k) end do else - !do j = pdims%j_start, pdims%j_end do i = pdims%i_start, pdims%i_end rhokm(i,j,k) = zero rhokh(i,j,k) = zero @@ -2559,7 +2526,6 @@ subroutine bdy_expl2 ( & !$OMP PARALLEL DEFAULT(SHARED) private(i, k) -!do j = pdims%j_start, pdims%j_end if (BL_diag%l_zht) then !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end @@ -2571,7 +2537,6 @@ subroutine bdy_expl2 ( & !$OMP end do NOWAIT end if -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end ! Max height of BL turbulent mixing @@ -2601,7 +2566,6 @@ subroutine bdy_expl2 ( & end do !$OMP end do NOWAIT -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end if (dynamic_bl_diag(i,j)) then @@ -2663,7 +2627,6 @@ subroutine bdy_expl2 ( & ftl,fqw,wtrac_bl & ) -!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(none) private(k, i) & !$OMP SHARED(j, bl_levels, pdims, f_ngstress, weight_1dbl) do k = 2, bl_levels @@ -2688,7 +2651,6 @@ subroutine bdy_expl2 ( & !$OMP ntml_save,shallowc) !$OMP do SCHEDULE(STATIC) -!do j = pdims%j_start, pdims%j_end do i = pdims%i_start, pdims%i_end uw0(i,j) = -rhokm(i,j,1) * & ( u_p(i,j,1) - u_0_px(i,j) ) @@ -2697,7 +2659,6 @@ subroutine bdy_expl2 ( & end do !$OMP end do NOWAIT if (formdrag == explicit_stress) then -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end uw0(i,j) = uw0(i,j) - tau_fd_x(i,j,1) @@ -2715,7 +2676,6 @@ subroutine bdy_expl2 ( & if (non_local_bl == on) then -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end if ( cumulus(i,j) ) then @@ -2730,7 +2690,6 @@ subroutine bdy_expl2 ( & else ! non-local scheme not being used so always use local scheme values -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end ntml(i,j) = ntml_local(i,j) @@ -2739,7 +2698,6 @@ subroutine bdy_expl2 ( & end if -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end wstar(i,j) = zero @@ -2797,7 +2755,6 @@ subroutine bdy_expl2 ( & ! if cumulus is still true in order to trigger convection ! at the correct level -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end if ( .not. cumulus(i,j) ) l_shallow(i,j) = .false. @@ -2811,7 +2768,6 @@ subroutine bdy_expl2 ( & ! 0.0 if .not. CUMULUS !----------------------------------------------------------------------- -!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end if ( cumulus(i,j) .and. l_shallow(i,j) ) then From b30fa49c92f82e97a12dc46af274a714e3f117d6 Mon Sep 17 00:00:00 2001 From: MetBenjaminWent Date: Tue, 3 Feb 2026 16:20:48 +0000 Subject: [PATCH 5/5] Tidy groups - keep for the moment, before removing in PR --- rose-stem/site/meto/groups/groups_lfric_atm.cylc | 6 +++--- ...p_gal9_noukca_1T-C48_MG_ex1a_cce_full-debug-32bit.txt | 9 --------- ...p_gal9_noukca_4T-C48_MG_ex1a_cce_full-debug-32bit.txt | 9 --------- 3 files changed, 3 insertions(+), 21 deletions(-) delete mode 100644 rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_noukca_1T-C48_MG_ex1a_cce_full-debug-32bit.txt delete mode 100644 rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_noukca_4T-C48_MG_ex1a_cce_full-debug-32bit.txt diff --git a/rose-stem/site/meto/groups/groups_lfric_atm.cylc b/rose-stem/site/meto/groups/groups_lfric_atm.cylc index cc82b7ba..be3d7267 100644 --- a/rose-stem/site/meto/groups/groups_lfric_atm.cylc +++ b/rose-stem/site/meto/groups/groups_lfric_atm.cylc @@ -223,9 +223,9 @@ "lfric_atm_clim_gal9_chem_2T-C12_ex1a_cce_fast-debug-32bit", ], "ex1a_omp_C48_cce_full": [ - "lfric_atm_nwp_gal9_noukca_1T-C48_MG_ex1a_cce_full-debug-32bit", - "lfric_atm_nwp_gal9_noukca_2T-C48_MG_ex1a_cce_full-debug-32bit", - "lfric_atm_nwp_gal9_noukca_4T-C48_MG_ex1a_cce_full-debug-32bit", + "lfric_atm_nwp_gal9_1T-C48_MG_ex1a_cce_full-debug-32bit", + "lfric_atm_nwp_gal9_2T-C48_MG_ex1a_cce_full-debug-32bit", + "lfric_atm_nwp_gal9_4T-C48_MG_ex1a_cce_full-debug-32bit", ], "ex1a_omp_C48_cce": [ "lfric_atm_nwp_gal9_1T-C48_MG_ex1a_cce_fast-debug-32bit", diff --git a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_noukca_1T-C48_MG_ex1a_cce_full-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_noukca_1T-C48_MG_ex1a_cce_full-debug-32bit.txt deleted file mode 100644 index ecdac22f..00000000 --- a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_noukca_1T-C48_MG_ex1a_cce_full-debug-32bit.txt +++ /dev/null @@ -1,9 +0,0 @@ -Inner product checksum rho = 48D7FB47 -Inner product checksum theta = 5392A6D6 -Inner product checksum u = 6A97B5F7 -Inner product checksum mr1 = 41CCED78 -Inner product checksum mr2 = 39CD58A2 -Inner product checksum mr3 = 37A6C0C1 -Inner product checksum mr4 = 3970B330 -Inner product checksum mr5 = 0 -Inner product checksum mr6 = 0 diff --git a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_noukca_4T-C48_MG_ex1a_cce_full-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_noukca_4T-C48_MG_ex1a_cce_full-debug-32bit.txt deleted file mode 100644 index 4dce9788..00000000 --- a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_noukca_4T-C48_MG_ex1a_cce_full-debug-32bit.txt +++ /dev/null @@ -1,9 +0,0 @@ -Inner product checksum rho = 48D7FB59 -Inner product checksum theta = 5392A6CA -Inner product checksum u = 6A97B3E6 -Inner product checksum mr1 = 41CCEE63 -Inner product checksum mr2 = 39D05F5C -Inner product checksum mr3 = 37A7FEC3 -Inner product checksum mr4 = 3970D685 -Inner product checksum mr5 = 0 -Inner product checksum mr6 = 0