diff --git a/science/physics_schemes/source/boundary_layer/kmkhz_9c.F90 b/science/physics_schemes/source/boundary_layer/kmkhz_9c.F90 index 24008fac..de50501c 100644 --- a/science/physics_schemes/source/boundary_layer/kmkhz_9c.F90 +++ b/science/physics_schemes/source/boundary_layer/kmkhz_9c.F90 @@ -963,7 +963,7 @@ subroutine kmkhz_9c ( & ! z_uv at k=ntml+1 !Variables for cache-blocking -integer :: jj ! Block index +integer :: ii ! Block index integer(kind=jpim), parameter :: zhook_in = 0 integer(kind=jpim), parameter :: zhook_out = 1 @@ -1011,17 +1011,17 @@ subroutine kmkhz_9c ( & dz_disc_min = one_half * timestep * 1.0e-4_r_bl +j = 1 !Start OpenMP parallel region !$OMP PARALLEL DEFAULT(SHARED) & -!$OMP private (i, j, jj, k, kp, kl, km, i_wt, w_curv_nm, w_del_nm, w_curv, & +!$OMP private (i, ii, k, kp, kl, km, i_wt, w_curv_nm, w_del_nm, w_curv, & !$OMP r_d_eta, rho_dz, z_surf, sl_plume, qw_plume, q_liq_parc, q_liq_env, & !$OMP t_parc, q_vap_parc, t_dens_parc, t_dens_env, dpar_bydz, denv_bydz, & !$OMP z_rad_lim, k_rad_lim, dflw_inv, dfsw_inv, dfsw_top, svl_plume, & !$OMP svl_diff, monotonic_inv, svl_lapse, svl_lapse_base, & !$OMP quad_a, quad_bm, quad_c, dz_disc, dsvl_top, sl_lapse, qw_lapse, & !$OMP kp2, rht_k, rht_kp1, rht_kp2, interp, z_cbase ) - !----------------------------------------------------------------------- ! Index to subroutine KMKHZ9C @@ -1046,32 +1046,32 @@ subroutine kmkhz_9c ( & ! 1.1 Calculate Z_TOP (top of levels) and NTML from previous timestep !----------------------------------------------------------------------- !$OMP do SCHEDULE(STATIC) -do jj = pdims%j_start, pdims%j_end, bl_segment_size - do j = jj, min((jj+bl_segment_size)-1,pdims%j_end) - do i = pdims%i_start, pdims%i_end - ntml_prev(i,j) = 1 - end do +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) + !do i = pdims%i_start, pdims%i_end + ntml_prev(i,j) = 1 + !end do end do do k = 1, bl_levels-1 - do j = jj, min((jj+bl_segment_size)-1,pdims%j_end) - do i = pdims%i_start, pdims%i_end - z_top(i,j,k) = z_uv(i,j,k+1) - !------------------------------------------------------------ - !find NTML from previous TS (for accurate gradient adjustment - !of profiles - also note that NTML le BL_LEVELS-1) - !------------------------------------------------------------ - if ( zh_prev(i,j) >= z_uv(i,j,k+1) ) ntml_prev(i,j)=k - end do + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) + !do i = pdims%i_start, pdims%i_end + z_top(i,j,k) = z_uv(i,j,k+1) + !------------------------------------------------------------ + !find NTML from previous TS (for accurate gradient adjustment + !of profiles - also note that NTML le BL_LEVELS-1) + !------------------------------------------------------------ + if ( zh_prev(i,j) >= z_uv(i,j,k+1) ) ntml_prev(i,j)=k + !end do end do end do k = bl_levels - do j = jj, min((jj+bl_segment_size)-1,pdims%j_end) - do i = pdims%i_start, pdims%i_end - z_top(i,j,k) = z_uv(i,j,k) + dzl(i,j,k) - end do + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) + !do i = pdims%i_start, pdims%i_end + z_top(i,j,k) = z_uv(i,j,k) + dzl(i,j,k) + !end do end do -end do ! jj +end do ! ii !$OMP end do NOWAIT !----------------------------------------------------------------------- @@ -1080,12 +1080,12 @@ subroutine kmkhz_9c ( & !$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 - sl(i,j,k) = tl(i,j,k) + grcp * z_tq(i,j,k) - svl(i,j,k) = sl(i,j,k) * ( one + c_virtual*qw(i,j,k) ) - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + sl(i,j,k) = tl(i,j,k) + grcp * z_tq(i,j,k) + svl(i,j,k) = sl(i,j,k) * ( one + c_virtual*qw(i,j,k) ) end do + !end do end do !$OMP end do NOWAIT @@ -1133,205 +1133,203 @@ subroutine kmkhz_9c ( & do k = 1, bl_levels km = max( 1, k-1 ) kp = min( bl_levels, k+1 ) -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - + !do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end ! If the Intel Compiler vs 12 is used, a job fails with segmentation ! fault when the following loop is vectorised. Thus the following ! directive stop vectorization of this loop when an Intel Compiler ! is used. Other compilers, for example Cray, should vectorise this - ! loop automatically. -#if defined (IFORT_VERSION) -!DIR$ NOVECTOR -#endif - - do i = pdims%i_start, pdims%i_end - w_grad(i,j,k) = (w(i,j,k)-w(i,j,km))*rdz(i,j,k) - w_curv_nm = w(i,j,kp)-2.0_r_bl*w(i,j,k)+w(i,j,km) - w_del_nm = w(i,j,kp)-w(i,j,km) - w_nonmono(i,j,k) = 0 - if ( abs(w_curv_nm) > abs(w_del_nm) ) w_nonmono(i,j,k) = 1 - sls_inc(i,j,k) = zero - qls_inc(i,j,k) = zero - end do ! i - end do ! j -!$OMP end do NOWAIT + ! loop automatically. +! #if defined (IFORT_VERSION) +! !DIR$ NOVECTOR +! #endif + w_grad(i,j,k) = (w(i,j,k)-w(i,j,km))*rdz(i,j,k) + w_curv_nm = w(i,j,kp)-2.0_r_bl*w(i,j,k)+w(i,j,km) + w_del_nm = w(i,j,kp)-w(i,j,km) + w_nonmono(i,j,k) = 0 + if ( abs(w_curv_nm) > abs(w_del_nm) ) w_nonmono(i,j,k) = 1 + sls_inc(i,j,k) = zero + qls_inc(i,j,k) = zero + end do ! i + !$OMP end do NOWAIT + !end do ! j end do ! k !$OMP BARRIER !$OMP do SCHEDULE(DYNAMIC) -do jj = pdims%j_start, pdims%j_end, bl_segment_size +do ii = pdims%i_start, pdims%i_end, bl_segment_size do k = 2, bl_levels-1 - do j = jj, min((jj+bl_segment_size)-1,pdims%j_end) - do i = pdims%i_start, pdims%i_end + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) + !do i = pdims%i_start, pdims%i_end - if ( etadot(i,j,k) < - tiny(one) .and. & - etadot(i,j,k-1)< - tiny(one) ) then - !----------------------------------------------------------- - ! Only needed in subsidence regions - ! Also don't attempt coupling with dynamics if w has - ! significant vertical structure - !----------------------------------------------------------- - w_curv = (w_grad(i,j,k+1)-w_grad(i,j,k))/dzl(i,j,k) - if (abs(w_curv) > 1.0e-6_r_bl .and. w_nonmono(i,j,k) == 1) then - ! large curvature at a turning point - sls_inc(i,j,k-1) = zero ! need to make sure increments in - qls_inc(i,j,k-1) = zero ! level below are also set to zero - etadot(i,j,k-1) = zero - etadot(i,j,k) = zero - etadot(i,j,k+1) = zero - w(i,j,k-1) = zero - w(i,j,k) = zero - w(i,j,k+1) = zero - else - r_d_eta=1.0/(eta_theta_levels(k+1)-eta_theta_levels(k)) - sls_inc(i,j,k) = - etadot(i,j,k) * r_d_eta & - * ( sl(i,j,k+1) - sl(i,j,k) ) - qls_inc(i,j,k) = - etadot(i,j,k) * r_d_eta & - * ( qw(i,j,k+1) - qw(i,j,k) ) - end if ! safe to calculate increments - end if + if ( etadot(i,j,k) < - tiny(one) .and. & + etadot(i,j,k-1)< - tiny(one) ) then + !----------------------------------------------------------- + ! Only needed in subsidence regions + ! Also don't attempt coupling with dynamics if w has + ! significant vertical structure + !----------------------------------------------------------- + w_curv = (w_grad(i,j,k+1)-w_grad(i,j,k))/dzl(i,j,k) + if (abs(w_curv) > 1.0e-6_r_bl .and. w_nonmono(i,j,k) == 1) then + ! large curvature at a turning point + sls_inc(i,j,k-1) = zero ! need to make sure increments in + qls_inc(i,j,k-1) = zero ! level below are also set to zero + etadot(i,j,k-1) = zero + etadot(i,j,k) = zero + etadot(i,j,k+1) = zero + w(i,j,k-1) = zero + w(i,j,k) = zero + w(i,j,k+1) = zero + else + r_d_eta=1.0/(eta_theta_levels(k+1)-eta_theta_levels(k)) + sls_inc(i,j,k) = - etadot(i,j,k) * r_d_eta & + * ( sl(i,j,k+1) - sl(i,j,k) ) + qls_inc(i,j,k) = - etadot(i,j,k) * r_d_eta & + * ( qw(i,j,k+1) - qw(i,j,k) ) + end if ! safe to calculate increments + end if - end do + !end do end do end do -end do !jj +end do !ii !$OMP end do - +! do j = pdims%j_start, pdims%j_end ! Repeat for necessary parts of last 2 loops for water tracers if (l_wtrac) then !$OMP do SCHEDULE(STATIC) do i_wt = 1, n_wtrac do k = 1, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - qls_inc_wtrac(i,j,k,i_wt) = zero - end do ! i - end do ! j + do i = pdims%i_start, pdims%i_end + qls_inc_wtrac(i,j,k,i_wt) = zero + end do ! i end do ! k end do ! i_wt !$OMP end do + !end do ! j -!$OMP do SCHEDULE(STATIC) - do jj = pdims%j_start, pdims%j_end, bl_segment_size +! Convert to dynamic +!$OMP do SCHEDULE(DYNAMIC) + do ii = pdims%i_start, pdims%i_end, bl_segment_size do k = 2, bl_levels-1 - do j = jj, min((jj+bl_segment_size)-1,pdims%j_end) - do i = pdims%i_start, pdims%i_end - - if ( etadot(i,j,k) < - tiny(one) .and. & - etadot(i,j,k-1)< - tiny(one) ) then - !----------------------------------------------------------- - ! Only needed in subsidence regions - ! Also don't attempt coupling with dynamics if w has - ! significant vertical structure - !----------------------------------------------------------- - w_curv = (w_grad(i,j,k+1)-w_grad(i,j,k))/dzl(i,j,k) - if (abs(w_curv) > 1.0e-6_r_bl .and. w_nonmono(i,j,k) == 1) then - ! large curvature at a turning point - do i_wt = 1, n_wtrac - ! level below are also set to zero - qls_inc_wtrac(i,j,k-1,i_wt) = zero - end do - else - kp = k+1 - km = kp-1 - r_d_eta=one/(eta_theta_levels(kp)-eta_theta_levels(km)) - do i_wt = 1, n_wtrac - qls_inc_wtrac(i,j,k,i_wt) = - etadot(i,j,k) * r_d_eta & - * ( wtrac_bl(i_wt)%qw(i,j,kp) - wtrac_bl(i_wt)%qw(i,j,km) ) - end do - end if ! safe to calculate increments - end if - end do + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) + !do i = pdims%i_start, pdims%i_end + + if ( etadot(i,j,k) < - tiny(one) .and. & + etadot(i,j,k-1)< - tiny(one) ) then + !----------------------------------------------------------- + ! Only needed in subsidence regions + ! Also don't attempt coupling with dynamics if w has + ! significant vertical structure + !----------------------------------------------------------- + w_curv = (w_grad(i,j,k+1)-w_grad(i,j,k))/dzl(i,j,k) + if (abs(w_curv) > 1.0e-6_r_bl .and. w_nonmono(i,j,k) == 1) then + ! large curvature at a turning point + do i_wt = 1, n_wtrac + ! level below are also set to zero + qls_inc_wtrac(i,j,k-1,i_wt) = zero + end do + else + kp = k+1 + km = kp-1 + r_d_eta=one/(eta_theta_levels(kp)-eta_theta_levels(km)) + do i_wt = 1, n_wtrac + qls_inc_wtrac(i,j,k,i_wt) = - etadot(i,j,k) * r_d_eta & + * ( wtrac_bl(i_wt)%qw(i,j,kp) - wtrac_bl(i_wt)%qw(i,j,km) ) + end do + end if ! safe to calculate increments + end if + !end do end do - end do - end do !jj + end do !i + end do !ii !$OMP end do end if ! l_wtrac +!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 - ! Non-turbulent fluxes are defined relative to the surface - ! so set them to zero at the surface - frad(i,j,1) = zero - frad_lw(i,j,1) = zero - frad_sw(i,j,1) = zero - fsubs(i,j,1,1) = zero ! for heat - fsubs(i,j,1,2) = zero ! for humidity - fmic(i,j,1,1) = zero ! for heat - fmic(i,j,1,2) = zero ! for humidity - - ! The two expressions are maintained as comment to highlight that - ! ft_nt and fq_net need to be updated if there are any eventual - ! changes in fmic/frad/fsubs at level 1. A job fails with - ! segmentation fault when the Intel compiler is used if the two - ! arrays are not set to 0.0 directly at level 1. - - !ft_nt(i,j,1) = frad(i,j,1) + fmic(i,j,1,1) + fsubs(i,j,1,1) - !fq_nt(i,j,1) = fmic(i,j,1,2) + fsubs(i,j,1,2) - - ft_nt(i,j,1) = zero - fq_nt(i,j,1) = zero +do i = pdims%i_start, pdims%i_end + ! Non-turbulent fluxes are defined relative to the surface + ! so set them to zero at the surface + frad(i,j,1) = zero + frad_lw(i,j,1) = zero + frad_sw(i,j,1) = zero + fsubs(i,j,1,1) = zero ! for heat + fsubs(i,j,1,2) = zero ! for humidity + fmic(i,j,1,1) = zero ! for heat + fmic(i,j,1,2) = zero ! for humidity + + ! The two expressions are maintained as comment to highlight that + ! ft_nt and fq_net need to be updated if there are any eventual + ! changes in fmic/frad/fsubs at level 1. A job fails with + ! segmentation fault when the Intel compiler is used if the two + ! arrays are not set to 0.0 directly at level 1. + + !ft_nt(i,j,1) = frad(i,j,1) + fmic(i,j,1,1) + fsubs(i,j,1,1) + !fq_nt(i,j,1) = fmic(i,j,1,2) + fsubs(i,j,1,2) + + ft_nt(i,j,1) = zero + fq_nt(i,j,1) = zero - end do end do !$OMP end do NOWAIT +!end do ! Repeat for necessary parts of last loop for water tracers + ! do j = pdims%j_start, pdims%j_end if (l_wtrac) then !$OMP do SCHEDULE(STATIC) do i_wt = 1, n_wtrac - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ! Non-turbulent fluxes are defined relative to the surface - ! so set them to zero at the surface - fsubs_wtrac(i,j,1,i_wt) = zero - fmic_wtrac(i,j,1,i_wt) = zero - wtrac_bl(i_wt)%fq_nt(i,j,1) = zero - end do + do i = pdims%i_start, pdims%i_end + ! Non-turbulent fluxes are defined relative to the surface + ! so set them to zero at the surface + fsubs_wtrac(i,j,1,i_wt) = zero + fmic_wtrac(i,j,1,i_wt) = zero + wtrac_bl(i_wt)%fq_nt(i,j,1) = zero end do end do !$OMP end do end if +! end do ! This is the most computational expensive loop of the subroutine. The ! following parallelisation obtains lower times than the ones obtained -! by using the jj/bl_segment_size technique. While it would be +! by using the ii/bl_segment_size technique. While it would be ! possible to remove the "df_over_cp", "dfmic", and "dfsubs" arrays, ! without them some jobs lose bit-comparability with the Intel ! compiler. do k = 1, bl_levels kp = k+1 + !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 - rho_dz = rho_mix_tq(i,j,k) * dzl(i,j,k) - - dflw_over_cp(i,j,k) = - rad_hr(i,j,1,k) * rho_dz - dfsw_over_cp(i,j,k) = - rad_hr(i,j,2,k) * rho_dz - df_over_cp(i,j,k) = dflw_over_cp(i,j,k) + dfsw_over_cp(i,j,k) - - dfmic(i,j,k,1) = - micro_tends(i,j,1,k) * rho_dz - dfmic(i,j,k,2) = - micro_tends(i,j,2,k) * rho_dz - dfsubs(i,j,k,1) = - sls_inc(i,j,k) * rho_dz - dfsubs(i,j,k,2) = - qls_inc(i,j,k) * rho_dz - - frad(i,j,kp) = frad(i,j,k) + df_over_cp(i,j,k) - frad_lw(i,j,kp)= frad_lw(i,j,k) + dflw_over_cp(i,j,k) - frad_sw(i,j,kp)= frad_sw(i,j,k) + dfsw_over_cp(i,j,k) - fsubs(i,j,kp,1)= fsubs(i,j,k,1) + dfsubs(i,j,k,1) - fsubs(i,j,kp,2)= fsubs(i,j,k,2) + dfsubs(i,j,k,2) - fmic(i,j,kp,1) = fmic(i,j,k,1) + dfmic(i,j,k,1) - fmic(i,j,kp,2) = fmic(i,j,k,2) + dfmic(i,j,k,2) - - ft_nt(i,j,kp) = frad(i,j,kp) + fmic(i,j,kp,1) + fsubs(i,j,kp,1) - fq_nt(i,j,kp) = fmic(i,j,kp,2) + fsubs(i,j,kp,2) - end do ! i - end do !j + do i = pdims%i_start, pdims%i_end + rho_dz = rho_mix_tq(i,j,k) * dzl(i,j,k) + + dflw_over_cp(i,j,k) = - rad_hr(i,j,1,k) * rho_dz + dfsw_over_cp(i,j,k) = - rad_hr(i,j,2,k) * rho_dz + df_over_cp(i,j,k) = dflw_over_cp(i,j,k) + dfsw_over_cp(i,j,k) + + dfmic(i,j,k,1) = - micro_tends(i,j,1,k) * rho_dz + dfmic(i,j,k,2) = - micro_tends(i,j,2,k) * rho_dz + dfsubs(i,j,k,1) = - sls_inc(i,j,k) * rho_dz + dfsubs(i,j,k,2) = - qls_inc(i,j,k) * rho_dz + + frad(i,j,kp) = frad(i,j,k) + df_over_cp(i,j,k) + frad_lw(i,j,kp)= frad_lw(i,j,k) + dflw_over_cp(i,j,k) + frad_sw(i,j,kp)= frad_sw(i,j,k) + dfsw_over_cp(i,j,k) + fsubs(i,j,kp,1)= fsubs(i,j,k,1) + dfsubs(i,j,k,1) + fsubs(i,j,kp,2)= fsubs(i,j,k,2) + dfsubs(i,j,k,2) + fmic(i,j,kp,1) = fmic(i,j,k,1) + dfmic(i,j,k,1) + fmic(i,j,kp,2) = fmic(i,j,k,2) + dfmic(i,j,k,2) + + ft_nt(i,j,kp) = frad(i,j,kp) + fmic(i,j,kp,1) + fsubs(i,j,kp,1) + fq_nt(i,j,kp) = fmic(i,j,kp,2) + fsubs(i,j,kp,2) + end do ! i !$OMP end do NOWAIT +!end do !j end do ! k ! Repeat necessary parts of last loop for water tracer @@ -1339,22 +1337,22 @@ subroutine kmkhz_9c ( & !$OMP do SCHEDULE(STATIC) do i_wt = 1, n_wtrac do k = 1, bl_levels - do j = pdims%j_start, pdims%j_end - kp = k+1 - do i = pdims%i_start, pdims%i_end - rho_dz = rho_mix_tq(i,j,k) * dzl(i,j,k) - dfmic_wtrac(i,j,k,i_wt) = & - - wtrac_as(i_wt)%micro_tends(i,j,k) * rho_dz - dfsubs_wtrac(i,j,k,i_wt) = - qls_inc_wtrac(i,j,k,i_wt) * rho_dz - - fsubs_wtrac(i,j,kp,i_wt)= & - fsubs_wtrac(i,j,k,i_wt) + dfsubs_wtrac(i,j,k,i_wt) - fmic_wtrac(i,j,kp,i_wt) = & - fmic_wtrac(i,j,k,i_wt) + dfmic_wtrac(i,j,k,i_wt) - wtrac_bl(i_wt)%fq_nt(i,j,kp) = & - fmic_wtrac(i,j,kp,i_wt) + fsubs_wtrac(i,j,kp,i_wt) - end do ! i - end do !j + !do j = pdims%j_start, pdims%j_end + kp = k+1 + do i = pdims%i_start, pdims%i_end + rho_dz = rho_mix_tq(i,j,k) * dzl(i,j,k) + dfmic_wtrac(i,j,k,i_wt) = & + - wtrac_as(i_wt)%micro_tends(i,j,k) * rho_dz + dfsubs_wtrac(i,j,k,i_wt) = - qls_inc_wtrac(i,j,k,i_wt) * rho_dz + + fsubs_wtrac(i,j,kp,i_wt)= & + fsubs_wtrac(i,j,k,i_wt) + dfsubs_wtrac(i,j,k,i_wt) + fmic_wtrac(i,j,kp,i_wt) = & + fmic_wtrac(i,j,k,i_wt) + dfmic_wtrac(i,j,k,i_wt) + wtrac_bl(i_wt)%fq_nt(i,j,kp) = & + fmic_wtrac(i,j,kp,i_wt) + fsubs_wtrac(i,j,kp,i_wt) + end do ! i + !end do !j end do ! k end do ! i_wt !$OMP end do NOWAIT @@ -1364,12 +1362,12 @@ subroutine kmkhz_9c ( & ! 1.4 Set UNSTABLE flag and find first level above surface layer !----------------------------------------------------------------------- !$OMP do SCHEDULE(DYNAMIC) -do jj = pdims%j_start, pdims%j_end, bl_segment_size - do j = jj, min(jj+bl_segment_size-1,pdims%j_end) - do i = pdims%i_start, pdims%i_end - unstable(i,j) = (fb_surf(i,j) > zero) - k_plume(i,j) = -1 - end do +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min(ii+bl_segment_size-1,pdims%i_end) + !do i = pdims%i_start, pdims%i_end + unstable(i,j) = (fb_surf(i,j) > zero) + k_plume(i,j) = -1 + !end do end do !------------------------------------------------------------ @@ -1381,33 +1379,33 @@ subroutine kmkhz_9c ( & !------------------------------------------------------------ do k = 1, bl_levels-1 - do j = jj, min(jj+bl_segment_size-1,pdims%j_end) - do i = pdims%i_start, pdims%i_end + do i = ii, min(ii+bl_segment_size-1,pdims%i_end) + !do i = pdims%i_start, pdims%i_end - if ( unstable(i,j) ) then - - z_surf = 0.1_r_bl * zh_prev(i,j) - if ( z_tq(i,j,k) >= z_surf .and. k_plume(i,j) == -1 ) then - !reached z_surf - k_plume(i,j)=k - end if - if ( svl(i,j,k+1) >= svl(i,j,k) & - .and. k_plume(i,j) == -1 ) then - !reached inversion - k_plume(i,j)=k - end if + if ( unstable(i,j) ) then + z_surf = 0.1_r_bl * zh_prev(i,j) + if ( z_tq(i,j,k) >= z_surf .and. k_plume(i,j) == -1 ) then + !reached z_surf + k_plume(i,j)=k end if - end do + if ( svl(i,j,k+1) >= svl(i,j,k) & + .and. k_plume(i,j) == -1 ) then + !reached inversion + k_plume(i,j)=k + end if + + end if + !end do end do end do - do j = jj, min(jj+bl_segment_size-1,pdims%j_end) - do i = pdims%i_start, pdims%i_end - if (k_plume(i,j) == -1) k_plume(i,j)=1 - end do + do i = ii, min(ii+bl_segment_size-1,pdims%i_end) + !do i = pdims%i_start, pdims%i_end + if (k_plume(i,j) == -1) k_plume(i,j)=1 + !end do end do -end do ! jj +end do ! ii !$OMP end do NOWAIT !----------------------------------------------------------------------- @@ -1423,147 +1421,147 @@ subroutine kmkhz_9c ( & !----------------------------------------------------------------------- ! Initialise variables +!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 - cloud_base(i,j) = .false. - dsc(i,j) = .false. - coupled(i,j) = .false. - zhsc(i,j) = zero - ntdsc(i,j) = 0 - end do +do i = pdims%i_start, pdims%i_end + cloud_base(i,j) = .false. + dsc(i,j) = .false. + coupled(i,j) = .false. + zhsc(i,j) = zero + ntdsc(i,j) = 0 end do !$OMP end do +!end do if ( .not. sc_diag_opt == sc_diag_all_rh_max ) then ! Use svl-gradient method to diagnose the Sc-top level at non-cumulus points -!$OMP do SCHEDULE(DYNAMIC) - do jj = pdims%j_start, pdims%j_end, bl_segment_size + !$OMP do SCHEDULE(DYNAMIC) + do ii = pdims%i_start, pdims%i_end, bl_segment_size do k = 3, bl_levels-1 - do j = jj, min(jj+bl_segment_size-1,pdims%j_end) - do i = pdims%i_start, pdims%i_end - - !----------------------------------------------------------------- - !..Find cloud-base (where cloud here means CF > SC_CFTOL) - !----------------------------------------------------------------- - - if ( .not. cumulus(i,j) .and. & - z_tq(i,j,k) < zmaxb_for_dsc .and. & - k > ntml(i,j)+1 .and. cf(i,j,k) > sc_cftol & - .and. .not. cloud_base(i,j) & - ! not yet found cloud-base - .and. .not. dsc(i,j) ) then - ! not yet found a Sc layer - cloud_base(i,j) = .true. - end if - if ( cloud_base(i,j) .and. .not. dsc(i,j) .and. & - ! found cloud-base but not yet reached cloud-top - cf(i,j,k+1) < sc_cftol .and. & - z_tq(i,j,k) < zmaxt_for_dsc & - ! got to cloud-top below ZMAXT_FOR_DSC - ) then - cloud_base(i,j) = .false. ! reset CLOUD_BASE - !----------------------------------------------------------- - ! Look to see if at least top of cloud is well mixed: - ! test SVL-gradient for top 2 pairs of levels, in case - ! cloud top extends into the inversion. - ! Parcel descent in Section 4.0 below will determine depth - ! of mixed layer. - !---------------------------------------------------------- - if ( (svl(i,j,k)-svl(i,j,k-1)) & - /(z_tq(i,j,k)-z_tq(i,j,k-1)) & - < max_svl_grad ) then - dsc(i,j) = .true. - ntdsc(i,j) = k - zhsc(i,j) = z_uv(i,j,ntdsc(i,j)+1) - else if ( (svl(i,j,k-1)-svl(i,j,k-2)) & - /(z_tq(i,j,k-1)-z_tq(i,j,k-2)) & + do i = ii, min(ii+bl_segment_size-1,pdims%i_end) + !do i = pdims%i_start, pdims%i_end + + !----------------------------------------------------------------- + !..Find cloud-base (where cloud here means CF > SC_CFTOL) + !----------------------------------------------------------------- + + if ( .not. cumulus(i,j) .and. & + z_tq(i,j,k) < zmaxb_for_dsc .and. & + k > ntml(i,j)+1 .and. cf(i,j,k) > sc_cftol & + .and. .not. cloud_base(i,j) & + ! not yet found cloud-base + .and. .not. dsc(i,j) ) then + ! not yet found a Sc layer + cloud_base(i,j) = .true. + end if + if ( cloud_base(i,j) .and. .not. dsc(i,j) .and. & + ! found cloud-base but not yet reached cloud-top + cf(i,j,k+1) < sc_cftol .and. & + z_tq(i,j,k) < zmaxt_for_dsc & + ! got to cloud-top below ZMAXT_FOR_DSC + ) then + cloud_base(i,j) = .false. ! reset CLOUD_BASE + !----------------------------------------------------------- + ! Look to see if at least top of cloud is well mixed: + ! test SVL-gradient for top 2 pairs of levels, in case + ! cloud top extends into the inversion. + ! Parcel descent in Section 4.0 below will determine depth + ! of mixed layer. + !---------------------------------------------------------- + if ( (svl(i,j,k)-svl(i,j,k-1)) & + /(z_tq(i,j,k)-z_tq(i,j,k-1)) & < max_svl_grad ) then - !--------------------------------------------------------- - ! Well-mixed layer with top at k-1 or k. Check whether - ! there is a buoyancy inversion between levels k-1 and k - ! in a manner similar to the surface-driven plume: compare - ! the buoyancy gradient between levels K-1 and K for an - ! undiluted parcel and the environment - !--------------------------------------------------------- - sl_plume = tl(i,j,k-1) + grcp * z_tq(i,j,k-1) - qw_plume = qw(i,j,k-1) - ! ------------------------------------------------------------ - ! calculate parcel water by linearising qsat about the - ! environmental temperature. - ! ------------------------------------------------------------ - if (t(i,j,k) > tm) then - q_liq_parc = max( zero, ( qw_plume - qs(i,j,k) - & - dqsdt(i,j,k)* & - ( sl_plume-grcp*z_tq(i,j,k)-t(i,j,k) ) & - ) *a_qs(i,j,k) ) - q_liq_env = max( zero, ( qw(i,j,k) - qs(i,j,k) & - -dqsdt(i,j,k)*( tl(i,j,k) - t(i,j,k) ) & - ) *a_qs(i,j,k) ) - ! add on the difference in the environment's ql as - ! calculated by the partial condensation scheme (using - ! some RH_CRIT value) and what it would be if - ! RH_CRIT=1. This then imitates partial condensation in - ! the parcel. - q_liq_parc = q_liq_parc + qcl(i,j,k) + qcf(i,j,k) & - - q_liq_env - t_parc = sl_plume - grcp * z_tq(i,j,k) + & - lcrcp*q_liq_parc - else - q_liq_parc = max( zero, ( qw_plume - qs(i,j,k) - & - dqsdt(i,j,k)* & - ( sl_plume - grcp*z_tq(i,j,k)-t(i,j,k) ) & - ) *a_qs(i,j,k) ) - q_liq_env = max( zero, ( qw(i,j,k) - qs(i,j,k) & - -dqsdt(i,j,k)*( tl(i,j,k) - t(i,j,k) ) & - ) *a_qs(i,j,k) ) - ! add on difference in environment's ql between RH_CRIT and - ! RH_CRIT=1 - q_liq_parc = q_liq_parc + qcl(i,j,k) + qcf(i,j,k) & - - q_liq_env - t_parc = sl_plume - grcp * z_tq(i,j,k) + & - lsrcp*q_liq_parc - end if - q_vap_parc=qw_plume - q_liq_parc - - t_dens_parc=t_parc*(one+c_virtual*q_vap_parc-q_liq_parc) - t_dens_env=t(i,j,k)* & - (one+c_virtual*q(i,j,k)-qcl(i,j,k)-qcf(i,j,k)) - ! find vertical gradients in parcel and environment SVL - ! (using values from level below (K-1)) - env_svl_km1(i,j) = t(i,j,k-1) * ( one+c_virtual*q(i,j,k-1) & - -qcl(i,j,k-1)-qcf(i,j,k-1) ) + grcp*z_tq(i,j,k-1) - dpar_bydz=(t_dens_parc+grcp*z_tq(i,j,k)- & - env_svl_km1(i,j)) / & - (z_tq(i,j,k)-z_tq(i,j,k-1)) - denv_bydz=(t_dens_env+grcp*z_tq(i,j,k)- & - env_svl_km1(i,j))/ & - (z_tq(i,j,k)-z_tq(i,j,k-1)) - - if ( denv_bydz > 1.25_r_bl*dpar_bydz ) then - ! there is an inversion between levels K-1 and K - if ( k >= ntml(i,j)+3 ) then - ! if NTDSC == NTML+1 then assume we're looking - ! at the same inversion and so don't set DSC - ntdsc(i,j) = k-1 - zhsc(i,j) = z_uv(i,j,ntdsc(i,j)+1) - dsc(i,j) = .true. - end if - else - ! no inversion between levels K-1 and K, assume there - ! is an inversion between K and K+1 because of CF change - ntdsc(i,j) = k + dsc(i,j) = .true. + ntdsc(i,j) = k + zhsc(i,j) = z_uv(i,j,ntdsc(i,j)+1) + else if ( (svl(i,j,k-1)-svl(i,j,k-2)) & + /(z_tq(i,j,k-1)-z_tq(i,j,k-2)) & + < max_svl_grad ) then + !--------------------------------------------------------- + ! Well-mixed layer with top at k-1 or k. Check whether + ! there is a buoyancy inversion between levels k-1 and k + ! in a manner similar to the surface-driven plume: compare + ! the buoyancy gradient between levels K-1 and K for an + ! undiluted parcel and the environment + !--------------------------------------------------------- + sl_plume = tl(i,j,k-1) + grcp * z_tq(i,j,k-1) + qw_plume = qw(i,j,k-1) + ! ------------------------------------------------------------ + ! calculate parcel water by linearising qsat about the + ! environmental temperature. + ! ------------------------------------------------------------ + if (t(i,j,k) > tm) then + q_liq_parc = max( zero, ( qw_plume - qs(i,j,k) - & + dqsdt(i,j,k)* & + ( sl_plume-grcp*z_tq(i,j,k)-t(i,j,k) ) & + ) *a_qs(i,j,k) ) + q_liq_env = max( zero, ( qw(i,j,k) - qs(i,j,k) & + -dqsdt(i,j,k)*( tl(i,j,k) - t(i,j,k) ) & + ) *a_qs(i,j,k) ) + ! add on the difference in the environment's ql as + ! calculated by the partial condensation scheme (using + ! some RH_CRIT value) and what it would be if + ! RH_CRIT=1. This then imitates partial condensation in + ! the parcel. + q_liq_parc = q_liq_parc + qcl(i,j,k) + qcf(i,j,k) & + - q_liq_env + t_parc = sl_plume - grcp * z_tq(i,j,k) + & + lcrcp*q_liq_parc + else + q_liq_parc = max( zero, ( qw_plume - qs(i,j,k) - & + dqsdt(i,j,k)* & + ( sl_plume - grcp*z_tq(i,j,k)-t(i,j,k) ) & + ) *a_qs(i,j,k) ) + q_liq_env = max( zero, ( qw(i,j,k) - qs(i,j,k) & + -dqsdt(i,j,k)*( tl(i,j,k) - t(i,j,k) ) & + ) *a_qs(i,j,k) ) + ! add on difference in environment's ql between RH_CRIT and + ! RH_CRIT=1 + q_liq_parc = q_liq_parc + qcl(i,j,k) + qcf(i,j,k) & + - q_liq_env + t_parc = sl_plume - grcp * z_tq(i,j,k) + & + lsrcp*q_liq_parc + end if + q_vap_parc=qw_plume - q_liq_parc + + t_dens_parc=t_parc*(one+c_virtual*q_vap_parc-q_liq_parc) + t_dens_env=t(i,j,k)* & + (one+c_virtual*q(i,j,k)-qcl(i,j,k)-qcf(i,j,k)) + ! find vertical gradients in parcel and environment SVL + ! (using values from level below (K-1)) + env_svl_km1(i,j) = t(i,j,k-1) * ( one+c_virtual*q(i,j,k-1) & + -qcl(i,j,k-1)-qcf(i,j,k-1) ) + grcp*z_tq(i,j,k-1) + dpar_bydz=(t_dens_parc+grcp*z_tq(i,j,k)- & + env_svl_km1(i,j)) / & + (z_tq(i,j,k)-z_tq(i,j,k-1)) + denv_bydz=(t_dens_env+grcp*z_tq(i,j,k)- & + env_svl_km1(i,j))/ & + (z_tq(i,j,k)-z_tq(i,j,k-1)) + + if ( denv_bydz > 1.25_r_bl*dpar_bydz ) then + ! there is an inversion between levels K-1 and K + if ( k >= ntml(i,j)+3 ) then + ! if NTDSC == NTML+1 then assume we're looking + ! at the same inversion and so don't set DSC + ntdsc(i,j) = k-1 zhsc(i,j) = z_uv(i,j,ntdsc(i,j)+1) dsc(i,j) = .true. end if + else + ! no inversion between levels K-1 and K, assume there + ! is an inversion between K and K+1 because of CF change + ntdsc(i,j) = k + zhsc(i,j) = z_uv(i,j,ntdsc(i,j)+1) + dsc(i,j) = .true. end if end if - end do ! i + end if + !end do ! i end do ! j end do ! k - end do ! jj + end do ! ii !$OMP end do end if ! ( .not. sc_diag_opt == sc_diag_all_rh_max ) @@ -1581,141 +1579,143 @@ subroutine kmkhz_9c ( & ! Options that diagnose the Sc-top using the max total-water RH method... ! j-loop outermost to allow parallelisation (k-loop is sequential) -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - - do i = pdims%i_start, pdims%i_end - ! Initialise max RH in the column to zero - rht_max(i,j) = zero - end do - - do k = 1, bl_levels-1 - do i = pdims%i_start, pdims%i_end - if ( ( cumulus(i,j) .and. z_tq(i,j,k) > z_lcl(i,j) ) .or. & - ( sc_diag_opt == sc_diag_all_rh_max .and. (.not. cumulus(i,j)) & - .and. k > ntml(i,j)+1 ) ) then - ! If sc_diag_opt == sc_diag_cu_rh_max, only check cumulus points - ! at heights above the LCL. - ! If sc_diag_opt == sc_diag_all_rh_max, also check non-cumulus points - ! at heights above ntml+1. - - ! Find max of RHt = qw/qsat(Tl) in the column - ! (with extra check that value exceeds next level, to exclude - ! points at bl_levels-1 where a larger value occurs at bl_levels) - if ( cf(i,j,k) > sc_cftol .and. z_tq(i,j,k) < zmaxt_for_dsc ) then - if ( l_mr_physics ) then - call qsat_wat_mix( rht_k, tl(i,j,k), p(i,j,k) ) - call qsat_wat_mix( rht_kp1, tl(i,j,k+1), p(i,j,k+1) ) - else - call qsat_wat( rht_k, tl(i,j,k), p(i,j,k) ) - call qsat_wat( rht_kp1, tl(i,j,k+1), p(i,j,k+1) ) - end if - rht_k = qw(i,j,k) / rht_k - rht_kp1 = qw(i,j,k+1) / rht_kp1 - if ( rht_k > rht_max(i,j) .and. rht_k > rht_kp1 ) then - ntdsc(i,j) = k - rht_max(i,j) = rht_k - end if - end if ! ( cf(i,j,k) > sc_cftol etc ) - end if ! ( cumulus(i,j) etc ) - end do ! i = pdims%i_start, pdims%i_end - end do ! k = 1, bl_levels-1 - + !do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + ! Initialise max RH in the column to zero + rht_max(i,j) = zero + end do + !$OMP end do + do k = 1, bl_levels-1 + !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end - if ( ntdsc(i,j) > 0 .and. ( .not. dsc(i,j) ) ) then - ! If we just found ntdsc above (exclude points where dsc already - ! set true by finding ntdsc at non-cumulus points earlier under - ! the option sc_diag_opt = sc_diag_cu_rh_max)... - ! Set flag indicating we found a potential Sc layer top - dsc(i,j) = .true. - ! Now we assume that theta-level ntdsc is wholly within the Sc-layer, - ! and that theta-level ntdsc+1 is composed of air with RH of - ! level ntdsc up to height zhsc, and air with RH of level ntdsc+2 - ! above that height. Interpolate to find the height of zhsc - ! (fraction of theta-level ntdsc+1 that is below the Sc-top). - k = ntdsc(i,j) - kp2 = min(k+2, bl_levels) ! Avoid out-of-bounds when k = bl_levels-1 - rht_k = rht_max(i,j) - if ( l_mr_physics ) then - call qsat_wat_mix( rht_kp1, tl(i,j,k+1), p(i,j,k+1) ) - call qsat_wat_mix( rht_kp2, tl(i,j,kp2), p(i,j,kp2) ) - else - call qsat_wat( rht_kp1, tl(i,j,k+1), p(i,j,k+1) ) - call qsat_wat( rht_kp2, tl(i,j,kp2), p(i,j,kp2) ) - end if - rht_kp1 = qw(i,j,k+1) / rht_kp1 - rht_kp2 = qw(i,j,kp2) / rht_kp2 - if ( rht_kp2 < rht_kp1 ) then - ! RHt(k+1) lies between RHt(k) and RHt(k+2); compute fraction - interp = (rht_kp1 - rht_kp2) / (rht_k - rht_kp2) - zhsc(i,j) = (one-interp) * z_uv(i,j,k+1) & - + interp * z_uv(i,j,kp2) - else - ! Rht(k+1) is a local minimum; can't construct k+1 as a fraction - ! of properties from k and k+2. Set zhsc to what it would be in - ! the limit RHt(k+2) = RHt(k+1) - zhsc(i,j) = z_uv(i,j,k+1) - end if - end if ! ( ntdsc(i,j) > 0 ) - end do + if ( ( cumulus(i,j) .and. z_tq(i,j,k) > z_lcl(i,j) ) .or. & + ( sc_diag_opt == sc_diag_all_rh_max .and. (.not. cumulus(i,j)) & + .and. k > ntml(i,j)+1 ) ) then + ! If sc_diag_opt == sc_diag_cu_rh_max, only check cumulus points + ! at heights above the LCL. + ! If sc_diag_opt == sc_diag_all_rh_max, also check non-cumulus points + ! at heights above ntml+1. + + ! Find max of RHt = qw/qsat(Tl) in the column + ! (with extra check that value exceeds next level, to exclude + ! points at bl_levels-1 where a larger value occurs at bl_levels) + if ( cf(i,j,k) > sc_cftol .and. z_tq(i,j,k) < zmaxt_for_dsc ) then + if ( l_mr_physics ) then + call qsat_wat_mix( rht_k, tl(i,j,k), p(i,j,k) ) + call qsat_wat_mix( rht_kp1, tl(i,j,k+1), p(i,j,k+1) ) + else + call qsat_wat( rht_k, tl(i,j,k), p(i,j,k) ) + call qsat_wat( rht_kp1, tl(i,j,k+1), p(i,j,k+1) ) + end if + rht_k = qw(i,j,k) / rht_k + rht_kp1 = qw(i,j,k+1) / rht_kp1 + if ( rht_k > rht_max(i,j) .and. rht_k > rht_kp1 ) then + ntdsc(i,j) = k + rht_max(i,j) = rht_k + end if + end if ! ( cf(i,j,k) > sc_cftol etc ) + end if ! ( cumulus(i,j) etc ) + end do ! i = pdims%i_start, pdims%i_end + !$OMP end do + end do ! k = 1, bl_levels-1 + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + if ( ntdsc(i,j) > 0 .and. ( .not. dsc(i,j) ) ) then + ! If we just found ntdsc above (exclude points where dsc already + ! set true by finding ntdsc at non-cumulus points earlier under + ! the option sc_diag_opt = sc_diag_cu_rh_max)... + ! Set flag indicating we found a potential Sc layer top + dsc(i,j) = .true. + ! Now we assume that theta-level ntdsc is wholly within the Sc-layer, + ! and that theta-level ntdsc+1 is composed of air with RH of + ! level ntdsc up to height zhsc, and air with RH of level ntdsc+2 + ! above that height. Interpolate to find the height of zhsc + ! (fraction of theta-level ntdsc+1 that is below the Sc-top). + k = ntdsc(i,j) + kp2 = min(k+2, bl_levels) ! Avoid out-of-bounds when k = bl_levels-1 + rht_k = rht_max(i,j) + if ( l_mr_physics ) then + call qsat_wat_mix( rht_kp1, tl(i,j,k+1), p(i,j,k+1) ) + call qsat_wat_mix( rht_kp2, tl(i,j,kp2), p(i,j,kp2) ) + else + call qsat_wat( rht_kp1, tl(i,j,k+1), p(i,j,k+1) ) + call qsat_wat( rht_kp2, tl(i,j,kp2), p(i,j,kp2) ) + end if + rht_kp1 = qw(i,j,k+1) / rht_kp1 + rht_kp2 = qw(i,j,kp2) / rht_kp2 + if ( rht_kp2 < rht_kp1 ) then + ! RHt(k+1) lies between RHt(k) and RHt(k+2); compute fraction + interp = (rht_kp1 - rht_kp2) / (rht_k - rht_kp2) + zhsc(i,j) = (one-interp) * z_uv(i,j,k+1) & + + interp * z_uv(i,j,kp2) + else + ! Rht(k+1) is a local minimum; can't construct k+1 as a fraction + ! of properties from k and k+2. Set zhsc to what it would be in + ! the limit RHt(k+2) = RHt(k+1) + zhsc(i,j) = z_uv(i,j,k+1) + end if + end if ! ( ntdsc(i,j) > 0 ) end do -!$OMP end do NOWAIT + !$OMP end do NOWAIT + !end do else if ( sc_diag_opt == sc_diag_cu_relax ) then ! Diagnosed simply if significant cloud fraction at ZHPAR ! below the height threshold zmaxt_for_dsc + !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 - k = ntpar(i,j) - if ( cumulus(i,j) .and. k < bl_levels ) then - ! cumulus layer within BL_LEVELS - if ( z_tq(i,j,k) < zmaxt_for_dsc .and. & - ! cloud top below zmaxt_for_dsc - ( max( cf(i,j,k-1),cf(i,j,k),cf(i,j,k+1) ) > sc_cftol ) & - ! cloud-top sufficiently cloudy - ) then - dsc(i,j) = .true. - zhsc(i,j) = zhpar(i,j) - ntdsc(i,j)= ntpar(i,j) - end if + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + k = ntpar(i,j) + if ( cumulus(i,j) .and. k < bl_levels ) then + ! cumulus layer within BL_LEVELS + if ( z_tq(i,j,k) < zmaxt_for_dsc .and. & + ! cloud top below zmaxt_for_dsc + ( max( cf(i,j,k-1),cf(i,j,k),cf(i,j,k+1) ) > sc_cftol ) & + ! cloud-top sufficiently cloudy + ) then + dsc(i,j) = .true. + zhsc(i,j) = zhpar(i,j) + ntdsc(i,j)= ntpar(i,j) end if - end do + end if end do -!$OMP end do NOWAIT + !$OMP end do NOWAIT +!end do else if ( sc_diag_opt == sc_diag_orig ) then ! Original code, only diagnosed if shallow cu or not l_param_conv +!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 ( (l_param_conv .and. & - l_shallow(i,j) .and. ntpar(i,j) < bl_levels ) & - ! shallow cumulus layer within BL_LEVELS - .or. (.not. l_param_conv .and. & - cumulus(i,j) .and. ntpar(i,j) < bl_levels ) ) then - ! cumulus layer and inversion found - if ( cf(i,j,ntpar(i,j)) > sc_cftol .or. & - cf(i,j,ntpar(i,j)+1) > sc_cftol ) then - ! cloudy - dsc(i,j) = .true. - zhsc(i,j) = zhpar(i,j) - ntdsc(i,j)= ntpar(i,j) - end if + do i = pdims%i_start, pdims%i_end + if ( (l_param_conv .and. & + l_shallow(i,j) .and. ntpar(i,j) < bl_levels ) & + ! shallow cumulus layer within BL_LEVELS + .or. (.not. l_param_conv .and. & + cumulus(i,j) .and. ntpar(i,j) < bl_levels ) ) then + ! cumulus layer and inversion found + if ( cf(i,j,ntpar(i,j)) > sc_cftol .or. & + cf(i,j,ntpar(i,j)+1) > sc_cftol ) then + ! cloudy + dsc(i,j) = .true. + zhsc(i,j) = zhpar(i,j) + ntdsc(i,j)= ntpar(i,j) end if - end do + end if end do -!$OMP end do NOWAIT + ! !$OMP end do NOWAIT + !$OMP end do NOWAIT +!end do end if ! test on sc_diag_opt if ( l_use_sml_dsc_fixes ) then ! Need to override "NOWAIT" on the previous blocks if going in here: -!$OMP BARRIER + !$OMP BARRIER ! If conv_diag has diagnosed a SML rising significantly above cloud-base, ! abort any DSC diagnosis higher-up. ! This is because at present, diagnosing an elevated DSC-layer prompts @@ -1727,26 +1727,26 @@ subroutine kmkhz_9c ( & ! The code currently only permits us to have one DSC-layer in the column, ! so if we have a cloudy boundary-layer, we need to reserve the one possible ! DSC for the case where the cloudy SML-top decouples to make a DSC. -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( unstable(i,j) .and. (.not. cumulus(i,j)) .and. dsc(i,j) ) then - if ( zh(i,j) - z_lcl(i,j) >= 400.0 ) then - ! Using cloud-layer depth >= 400m consistent with cumulus_test - ! in conv_diag, as this indicates that conv_diag was prevented - ! from diagnosing cumulus by its test on whether the qw profile - ! looks well-mixed. But it doesn't test on stability, so it can - ! give a very misleading diagnosis of a well-mixed cloudy layer. - ! In particular, it compares the cloud-layer qw gradient with the - ! sub-cloud-layer qw gradient, giving a well-mixed diagnosis if - ! they are similar even if both are large! - dsc(i,j) = .false. - ntdsc(i,j) = 0 - end if +! do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + if ( unstable(i,j) .and. (.not. cumulus(i,j)) .and. dsc(i,j) ) then + if ( zh(i,j) - z_lcl(i,j) >= 400.0 ) then + ! Using cloud-layer depth >= 400m consistent with cumulus_test + ! in conv_diag, as this indicates that conv_diag was prevented + ! from diagnosing cumulus by its test on whether the qw profile + ! looks well-mixed. But it doesn't test on stability, so it can + ! give a very misleading diagnosis of a well-mixed cloudy layer. + ! In particular, it compares the cloud-layer qw gradient with the + ! sub-cloud-layer qw gradient, giving a well-mixed diagnosis if + ! they are similar even if both are large! + dsc(i,j) = .false. + ntdsc(i,j) = 0 end if - end do + end if end do -!$OMP end do NOWAIT + !$OMP end do NOWAIT +! end do end if ! ( l_use_sml_dsc_fixes ) !----------------------------------------------------------------------- @@ -1756,25 +1756,25 @@ subroutine kmkhz_9c ( & !----------------------------------------------------------------------- ! Initialise variables !------------------------------ - +!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 - k_cloud_dsct(i,j) = 0 - df_dsct_over_cp(i,j) = zero - df_inv_dsc(i,j) = zero - end do +do i = pdims%i_start, pdims%i_end + k_cloud_dsct(i,j) = 0 + df_dsct_over_cp(i,j) = zero + df_inv_dsc(i,j) = zero end do !$OMP end do +!end do if (l_new_kcloudtop) then !--------------------------------------------------------------------- ! improved method of finding the k_cloud_dsct, the top of the mixed ! layer as seen by radiation !--------------------------------------------------------------------- -!$OMP do SCHEDULE(DYNAMIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end +!do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(DYNAMIC) + do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) !------------------------------------------------------------- ! Find k_cloud_dsct as the equivalent to ntdsc (the top level ! of the DSC) seen by radiation, which we take as the level two @@ -1810,13 +1810,17 @@ subroutine kmkhz_9c ( & end if ! DSC test separated out end do ! i - !----------------------------------------------------------------- - ! Find bottom grid-level (K_LEVEL) for cloud-top radiative flux - ! divergence: higher of base of LW radiatively cooled layer, - ! ZH and 0.5*ZHSC, since cooling must be in upper part of layer - ! in order to generate turbulence. - !----------------------------------------------------------------- - do i = pdims%i_start, pdims%i_end + end do ! ii + !$OMP end do + !----------------------------------------------------------------- + ! Find bottom grid-level (K_LEVEL) for cloud-top radiative flux + ! divergence: higher of base of LW radiatively cooled layer, + ! ZH and 0.5*ZHSC, since cooling must be in upper part of layer + ! in order to generate turbulence. + !----------------------------------------------------------------- + !$OMP do SCHEDULE(DYNAMIC) + do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) k_level(i,j) = k_cloud_dsct(i,j) if ( k_cloud_dsct(i,j) > 1 ) then k_rad_lim = ntml(i,j)+1 @@ -1831,17 +1835,20 @@ subroutine kmkhz_9c ( & end do ! k end if end do ! i - end do ! j -!$OMP end do + end do ! ii + !$OMP end do + !end do ! j + !##### else ! original method of finding the k_cloud_dsct, the top of the mixed layer ! as seen by radiation, found to be resolution dependent -!$OMP do SCHEDULE(DYNAMIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end + !do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(DYNAMIC) + do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) !------------------------------------------------------------- ! Find the layer with the greatest LW radiative flux jump in ! the upper half of the boundary layer and assume that this @@ -1860,21 +1867,23 @@ subroutine kmkhz_9c ( & ! than half the maximum. DF in level K_CLOUD_DSCT+1 is then ! included as DF_INV_DSC below. if (dflw_over_cp(i,j,k-1) > one_half*dflw_over_cp(i,j,k)) & - k_cloud_dsct(i,j) = k-1 + k_cloud_dsct(i,j) = k-1 df_dsct_over_cp(i,j) = dflw_over_cp(i,j,k) end if end do ! k end do ! i - - !----------------------------------------------------------------- - ! Find bottom grid-level (K_LEVEL) for cloud-top radiative flux - ! divergence: higher of base of LW radiatively cooled layer, - ! ZH and 0.5*ZHSC, since cooling must be in upper part of layer - ! in order to generate turbulence. - !----------------------------------------------------------------- - - do i = pdims%i_start, pdims%i_end + end do ! ii + !$OMP end do + !----------------------------------------------------------------- + ! Find bottom grid-level (K_LEVEL) for cloud-top radiative flux + ! divergence: higher of base of LW radiatively cooled layer, + ! ZH and 0.5*ZHSC, since cooling must be in upper part of layer + ! in order to generate turbulence. + !----------------------------------------------------------------- + !$OMP do SCHEDULE(DYNAMIC) + do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) k_level(i,j) = k_cloud_dsct(i,j) if ( k_cloud_dsct(i,j) > 1 ) then k_rad_lim = ntml(i,j)+1 @@ -1889,8 +1898,9 @@ subroutine kmkhz_9c ( & end do ! k end if end do ! i - end do ! j -!$OMP end do + end do ! ii + !$OMP end do + !end do ! j end if ! test on l_new_kcloudtop @@ -1902,10 +1912,10 @@ subroutine kmkhz_9c ( & ! representative of clear-air rad divergence and so subtract this ! `clear-air' part from the grid-level divergence. !----------------------------------------------------------------- - +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(DYNAMIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) if ( k_cloud_dsct(i,j) > 0 ) then dflw_inv = zero @@ -1914,10 +1924,10 @@ subroutine kmkhz_9c ( & k = k_cloud_dsct(i,j)+1 if ( k < bl_levels ) then dflw_inv = dflw_over_cp(i,j,k) & - - dflw_over_cp(i,j,k+1) & + - dflw_over_cp(i,j,k+1) & * dzl(i,j,k)/dzl(i,j,k+1) dfsw_inv = dfsw_over_cp(i,j,k) & - - dfsw_over_cp(i,j,k+1) & + - dfsw_over_cp(i,j,k+1) & * dzl(i,j,k)/dzl(i,j,k+1) else dflw_inv = dflw_over_cp(i,j,k) @@ -1929,12 +1939,12 @@ subroutine kmkhz_9c ( & df_inv_dsc(i,j) = dflw_inv + dfsw_inv df_dsct_over_cp(i,j) = frad_lw(i,j,k_cloud_dsct(i,j)+1) & - - frad_lw(i,j,k_level(i,j)) & - + dflw_inv + - frad_lw(i,j,k_level(i,j)) & + + dflw_inv dfsw_top = frad_sw(i,j,k_cloud_dsct(i,j)+1) & - - frad_sw(i,j,k_level(i,j)) & - + dfsw_inv + - frad_sw(i,j,k_level(i,j)) & + + dfsw_inv !----------------------------------------------------------- ! Combine SW and LW cloud-top divergences into a net @@ -1946,9 +1956,10 @@ subroutine kmkhz_9c ( & df_dsct_over_cp(i,j) = max( zero, & df_dsct_over_cp(i,j) + dfsw_frac * dfsw_top ) end if - end do -end do + end do !i +end do !ii !$OMP end do +!end do !----------------------------------------------------------------------- ! 2.4 Set NBDSC, the bottom level of the DSC layer. @@ -1967,86 +1978,85 @@ subroutine kmkhz_9c ( & !----------------------------------------------------------------------- !$OMP do SCHEDULE(DYNAMIC) -do jj = pdims%j_start, pdims%j_end, bl_segment_size - do j = jj, min(jj+bl_segment_size-1,pdims%j_end) - do i = pdims%i_start, pdims%i_end - nbdsc(i,j) = ntdsc(i,j)+1 - if (dsc(i,j)) then - ! The depth of the radiatively-cooled layer tends to be less - ! than O(50m) and so RAD_HR will be an underestimate of the - ! cooling tendency there. Compensate by multiplying by - ! DZL/50. (~4) Recall that DF_OVER_CP(I,j,K) = RAD_HR * - ! RHO_MIX_TQ * DZL Thus use cloud-top radiative forcing as - ! follows: - - k = ntdsc(i,j) - rho_dz = rho_mix_tq(i,j,k) * dzl(i,j,k) - svl_plume=svl(i,j,k-1) & - - ct_resid * dzl(i,j,k)*df_dsct_over_cp(i,j) / ( 50.0_r_bl*rho_dz ) - else - svl_plume=zero +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min(ii+bl_segment_size-1,pdims%i_end) + !do i = pdims%i_start, pdims%i_end + nbdsc(i,j) = ntdsc(i,j)+1 + if (dsc(i,j)) then + ! The depth of the radiatively-cooled layer tends to be less + ! than O(50m) and so RAD_HR will be an underestimate of the + ! cooling tendency there. Compensate by multiplying by + ! DZL/50. (~4) Recall that DF_OVER_CP(I,j,K) = RAD_HR * + ! RHO_MIX_TQ * DZL Thus use cloud-top radiative forcing as + ! follows: + + k = ntdsc(i,j) + rho_dz = rho_mix_tq(i,j,k) * dzl(i,j,k) + svl_plume=svl(i,j,k-1) & + - ct_resid * dzl(i,j,k)*df_dsct_over_cp(i,j) / ( 50.0_r_bl*rho_dz ) + else + svl_plume=zero + end if + do k = min(bl_levels-1, ntdsc(i,j)-1), 1, -1 + if (svl_plume < svl(i,j,k) ) then + nbdsc(i,j) = k+1 ! marks lowest level within ML end if - do k = min(bl_levels-1, ntdsc(i,j)-1), 1, -1 - if (svl_plume < svl(i,j,k) ) then - nbdsc(i,j) = k+1 ! marks lowest level within ML - end if - end do ! k - end do ! i - end do ! j -end do !jj + end do ! k + !end do ! j + end do ! i +end do !ii !$OMP end do !---------------------------------------------------------------------- ! 2.5 Tidy up variables associated with decoupled layer ! NOTE that NTDSC ge 3 if non-zero !---------------------------------------------------------------------- - +!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 - ! Note that ZHSC-Z_UV(NTML+2) may = 0, so this test comes first! - if (cumulus(i,j) .and. dsc(i,j)) & - nbdsc(i,j) = max( nbdsc(i,j), ntml(i,j)+2 ) - if ( ntdsc(i,j) >= 1 ) then - if ( nbdsc(i,j) < ntdsc(i,j)+1 ) then - if (sc_diag_opt==sc_diag_orig .or. sc_diag_opt==sc_diag_cu_relax) then - ! Initial zhsc is normally on rho-level ntdsc+1 - dscdepth(i,j) = z_uv(i,j,ntdsc(i,j)+1) - z_uv(i,j,nbdsc(i,j)) - else - ! Initial zhsc was interpolated between levels; use accurate zhsc - dscdepth(i,j) = zhsc(i,j) - z_uv(i,j,nbdsc(i,j)) - end if +do i = pdims%i_start, pdims%i_end + ! Note that ZHSC-Z_UV(NTML+2) may = 0, so this test comes first! + if (cumulus(i,j) .and. dsc(i,j)) & + nbdsc(i,j) = max( nbdsc(i,j), ntml(i,j)+2 ) + if ( ntdsc(i,j) >= 1 ) then + if ( nbdsc(i,j) < ntdsc(i,j)+1 ) then + if (sc_diag_opt==sc_diag_orig .or. sc_diag_opt==sc_diag_cu_relax) then + ! Initial zhsc is normally on rho-level ntdsc+1 + dscdepth(i,j) = z_uv(i,j,ntdsc(i,j)+1) - z_uv(i,j,nbdsc(i,j)) else + ! Initial zhsc was interpolated between levels; use accurate zhsc + dscdepth(i,j) = zhsc(i,j) - z_uv(i,j,nbdsc(i,j)) + end if + else + !---------------------------------------------------------- + ! Indicates a layer of zero depth + !---------------------------------------------------------- + if ( ( sc_diag_opt==sc_diag_orig .or. sc_diag_opt==sc_diag_cu_relax ) & + .and. ntdsc(i,j) == ntpar(i,j) ) then !---------------------------------------------------------- - ! Indicates a layer of zero depth + ! Indicates a Sc layer at the top of Cu: force mixing + ! over single layer. !---------------------------------------------------------- - if ( ( sc_diag_opt==sc_diag_orig .or. sc_diag_opt==sc_diag_cu_relax ) & - .and. ntdsc(i,j) == ntpar(i,j) ) then - !---------------------------------------------------------- - ! Indicates a Sc layer at the top of Cu: force mixing - ! over single layer. - !---------------------------------------------------------- - dscdepth(i,j) = dzl(i,j,ntdsc(i,j)) - else - dsc(i,j)=.false. - ntdsc(i,j)=0 - zhsc(i,j)=zero - df_dsct_over_cp(i,j) = zero - k_cloud_dsct(i,j) = 0 - df_inv_dsc(i,j) = zero - dscdepth(i,j) = zero - end if + dscdepth(i,j) = dzl(i,j,ntdsc(i,j)) + else + dsc(i,j)=.false. + ntdsc(i,j)=0 + zhsc(i,j)=zero + df_dsct_over_cp(i,j) = zero + k_cloud_dsct(i,j) = 0 + df_inv_dsc(i,j) = zero + dscdepth(i,j) = zero end if - else ! ntdsc == 0, just to make sure - dscdepth(i,j)=zero - dsc(i,j)=.false. - zhsc(i,j)=zero - df_dsct_over_cp(i,j) = zero - k_cloud_dsct(i,j) = 0 - df_inv_dsc(i,j) = zero end if - end do + else ! ntdsc == 0, just to make sure + dscdepth(i,j)=zero + dsc(i,j)=.false. + zhsc(i,j)=zero + df_dsct_over_cp(i,j) = zero + k_cloud_dsct(i,j) = 0 + df_inv_dsc(i,j) = zero + end if end do !$OMP end do NOWAIT +!end do !---------------------------------------------------------------------- !2.6 If decoupled cloud-layer found test to see if it is, in fact, @@ -2067,71 +2077,70 @@ subroutine kmkhz_9c ( & ! ON - taper off surface terms to zero for svl_diff between ! svl_coup and svl_coup_max; also ignore cumulus diags !------------------------------------------------------------- -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - coupled(i,j) = .false. - svl_diff_frac(i,j) = zero ! Fully coupled by default - if ( dsc(i,j) ) then - !------------------------------------------------------------- - ! Calculate cloud to surface mixed layer SVL difference - ! - avoid ntdsc as can be within base of inversion - !------------------------------------------------------------- - svl_diff = zero - if ( ntdsc(i,j) >= 2 ) & - svl_diff = svl(i,j,ntdsc(i,j)-1) - svl(i,j,ntml(i,j)) - if ( svl_diff < svl_coup_max ) then - coupled(i,j) = .true. - svl_diff_frac(i,j) = one - max( zero, & - (svl_diff-svl_coup)/(svl_coup_max-svl_coup) ) - ! to give 1 for svl_diff= 2 ) & + svl_diff = svl(i,j,ntdsc(i,j)-1) - svl(i,j,ntml(i,j)) + if ( svl_diff < svl_coup_max ) then + coupled(i,j) = .true. + svl_diff_frac(i,j) = one - max( zero, & + (svl_diff-svl_coup)/(svl_coup_max-svl_coup) ) + ! to give 1 for svl_diff one .and. ntml(i,j)+2 <= bl_levels ) then + dzh(i,j) > one .and. ntml(i,j)+2 <= bl_levels ) then if (zh(i,j)+dzh(i,j) > z_uv(i,j,ntml(i,j)+2) ) res_inv(i,j) = 1 end if @@ -2182,7 +2192,7 @@ subroutine kmkhz_9c ( & dqw_sml(i,j) = qw(i,j,k+1) - qw(i,j,k) if ( .not. cumulus(i,j) .and. .not. coupled(i,j) .and. & - res_inv(i,j) == 0 .and. k > 1 .and. k <= bl_levels-2 ) then + res_inv(i,j) == 0 .and. k > 1 .and. k <= bl_levels-2 ) then ! Require SVL and SL to be monotonically increasing ! and QW to be simply monotonic @@ -2240,7 +2250,7 @@ subroutine kmkhz_9c ( & if ( abs(quad_a) >= rbl_eps ) then ! ...quadratic if QUAD_A /= 0 dz_disc = ( quad_bm - sqrt( quad_bm*quad_bm & - - 4.0_r_bl*quad_a*quad_c ) & + - 4.0_r_bl*quad_a*quad_c ) & ) / (2.0_r_bl*quad_a) else ! ...linear if QUAD_A == 0 @@ -2287,9 +2297,9 @@ subroutine kmkhz_9c ( & qw_lapse = zero if ( k <= bl_levels-3 ) then sl_lapse = max( zero, & - ( sl(i,j,k+3) - sl(i,j,k+2) )*rdz(i,j,k+3) ) + ( sl(i,j,k+3) - sl(i,j,k+2) )*rdz(i,j,k+3) ) qw_lapse = min( zero, & - ( qw(i,j,k+3) - qw(i,j,k+2) )*rdz(i,j,k+3) ) + ( qw(i,j,k+3) - qw(i,j,k+2) )*rdz(i,j,k+3) ) if (l_wtrac) then ! Store method if (qw_lapse >= 0.0) then @@ -2305,12 +2315,12 @@ subroutine kmkhz_9c ( & ! Only reduce 2 level jump by at most half dsl_sml(i,j) = sl(i,j,k+2) - sl(i,j,k) dsl_sml(i,j) = dsl_sml(i,j) - & - min( one_half*dsl_sml(i,j), sl_lapse*dz_disc ) + min( one_half*dsl_sml(i,j), sl_lapse*dz_disc ) !----------------- ! Next QW jump !----------------- if ( qw(i,j,k+2) < qw(i,j,k+1) .and. & - qw(i,j,k+1) < qw(i,j,k) ) then + qw(i,j,k+1) < qw(i,j,k) ) then ! QW monotonically decreasing across inversion ! Only allow for QW lapse rate if both it and the ! 2 grid-level jump are negative (expected sign) @@ -2371,7 +2381,7 @@ subroutine kmkhz_9c ( & if ( monotonic_inv ) then if ( sc_diag_opt == sc_diag_cu_rh_max .or. & - sc_diag_opt == sc_diag_all_rh_max ) then + sc_diag_opt == sc_diag_all_rh_max ) then ! The initial zhsc can be between model-levels rather than exactly ! on a rho-level. ! Store height of base of DSC layer @@ -2403,11 +2413,11 @@ subroutine kmkhz_9c ( & quad_a = one_half*( svl_lapse - svl_lapse_base ) quad_bm = svl(i,j,k+2) - svl(i,j,k) & - - svl_lapse * ( z_tq(i,j,k+2)-z_uv(i,j,k+2) ) & - - svl_lapse_base * ( z_uv(i,j,k+1)-z_tq(i,j,k) + & + - svl_lapse * ( z_tq(i,j,k+2)-z_uv(i,j,k+2) ) & + - svl_lapse_base * ( z_uv(i,j,k+1)-z_tq(i,j,k) + & dzl(i,j,k+1) ) quad_c = dzl(i,j,k+1)*( svl(i,j,k+1) - svl(i,j,k) - & - svl_lapse_base * ( & + svl_lapse_base * ( & z_uv(i,j,k+1)-z_tq(i,j,k) + one_half*dzl(i,j,k+1) ) ) if ( quad_bm > zero ) then @@ -2420,7 +2430,7 @@ subroutine kmkhz_9c ( & if ( abs(quad_a) >= rbl_eps ) then ! ...quadratic if QUAD_A ne 0 dz_disc = ( quad_bm - sqrt( quad_bm*quad_bm & - - 4.0_r_bl*quad_a*quad_c ) & + - 4.0_r_bl*quad_a*quad_c ) & ) / (2.0_r_bl*quad_a) else ! ...linear if QUAD_A = 0 @@ -2465,7 +2475,7 @@ subroutine kmkhz_9c ( & zhsc(i,j) = z_uv(i,j,k+2) - dz_disc if ( sc_diag_opt == sc_diag_cu_rh_max .or. & - sc_diag_opt == sc_diag_all_rh_max ) then + sc_diag_opt == sc_diag_all_rh_max ) then ! The initial zhsc can be between model-levels rather than exactly ! on a rho-level. Assuming height of base of DSC layer stays the ! same, set new depth based on updated DSC top height: @@ -2487,9 +2497,9 @@ subroutine kmkhz_9c ( & qw_lapse = zero if ( k <= bl_levels-3 ) then sl_lapse = max( zero, & - ( sl(i,j,k+3) - sl(i,j,k+2) )*rdz(i,j,k+3) ) + ( sl(i,j,k+3) - sl(i,j,k+2) )*rdz(i,j,k+3) ) qw_lapse = min( zero, & - ( qw(i,j,k+3) - qw(i,j,k+2) )*rdz(i,j,k+3) ) + ( qw(i,j,k+3) - qw(i,j,k+2) )*rdz(i,j,k+3) ) if (l_wtrac) then ! Store method if (qw_lapse >= 0.0) then @@ -2505,12 +2515,12 @@ subroutine kmkhz_9c ( & ! Only reduce 2 level jump by at most half dsl_dsc(i,j) = sl(i,j,k+2) - sl(i,j,k) dsl_dsc(i,j) = dsl_dsc(i,j) - & - min( one_half*dsl_dsc(i,j), sl_lapse*dz_disc ) + min( one_half*dsl_dsc(i,j), sl_lapse*dz_disc ) !----------------- ! Next QW jump !----------------- if ( qw(i,j,k+2) < qw(i,j,k+1) .and. & - qw(i,j,k+1) < qw(i,j,k) ) then + qw(i,j,k+1) < qw(i,j,k) ) then ! QW monotonically decreasing across inversion ! Only allow for QW lapse rate if both it and the ! 2 grid-level jump are negative (expected sign) @@ -2540,9 +2550,10 @@ subroutine kmkhz_9c ( & end if ! monotonic inversion end if ! test on K lt BL_LEVELS-2 end if ! test on DSC - end do -end do + end do !i +end do !ii !$OMP end do +!end do !$OMP end PARALLEL @@ -2590,70 +2601,71 @@ subroutine kmkhz_9c ( & ! 4.1 Calculate gradient adjustment terms !----------------------------------------------------------------------- +!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL DEFAULT(SHARED) & -!$OMP private (i, j, jj, i_wt, k, kl, km, kp, kp2, kmax, wstar3, c_ws, w_m, & +!$OMP private (i, ii, i_wt, k, kl, km, kp, kp2, kmax, wstar3, c_ws, w_m, & !$OMP pr_neut, w_h, k_cff, virt_factor, z_cbase , zdsc_cbase, dsl_ga, & !$OMP dqw_ga, cfl_ml, cff_ml, dqw, dsl, dqcl, dqcf, db_disc, cu_depth_fac, & !$OMP k_rad_lim, z_rad_lim ,dfsw_inv, dflw_inv, dfsw_top, dsldz, cf_for_wb, & -!$OMP grad_t_adj_inv_rdz, grad_q_adj_inv_rdz) - +!$OMP grad_t_adj_inv_rdz, grad_q_adj_inv_rdz, denom) !$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) = zero - grad_q_adj(i,j) = zero - if ( unstable(i,j) ) then - ! Here this is an estimate of the gradient adjustment applied - ! the previous timestep (assumes T1_SD has not changed much, - ! which in turn assumes the surface fluxes have not) - if (flux_grad == Locketal2000) then - grad_t_adj(i,j) = min( max_t_grad , & - a_grad_adj * t1_sd(i,j) / zh_prev(i,j) ) - grad_q_adj(i,j) = zero - else if (flux_grad == HoltBov1993) then - ! Use constants from Holtslag and Boville (1993) - ! Conv limit GAMMA_TH = 10 *FTL1/(wstar*zh) - ! Neut limit GAMMA_TH = 7.2*wstar*FTL1/(ustar^2*zh) - wstar3 = fb_surf(i,j) * zh_prev(i,j) - c_ws = 0.6_r_bl - w_m =( v_s(i,j)**3 + c_ws*wstar3 )**one_third - - grad_t_adj(i,j) = a_ga_hb93*(wstar3**one_third)*ftl(i,j,1) & - / ( rhostar_gb(i,j)*w_m*w_m*zh_prev(i,j) ) - ! GRAD_Q_ADJ(I,j) = A_GA_HB93*(WSTAR3**one_third)*FQW(I,j,1) - ! / ( RHOSTAR_GB(I,j)*W_M*W_M*ZH_PREV(I,j) ) - grad_q_adj(i,j) = zero - else if (flux_grad == LockWhelan2006) then - ! Use constants from LockWhelan2006 - ! Conv limit GAMMA_TH = 10 *FTL1/(wstar*zh) - ! Neut limit GAMMA_TH = 7.5*FTL1/(ustar*zh) - wstar3 = fb_surf(i,j) * zh_prev(i,j) - c_ws = 0.42_r_bl ! = 0.75^3 - pr_neut = 0.75_r_bl - w_h = ( ( v_s(i,j)**3+c_ws*wstar3 )**one_third )/ pr_neut - - grad_t_adj(i,j) = a_ga_lw06 * ftl(i,j,1) & - / ( rhostar_gb(i,j)*w_h*zh_prev(i,j) ) - grad_q_adj(i,j) = a_ga_lw06 * fqw(i,j,1) & - / ( rhostar_gb(i,j)*w_h*zh_prev(i,j) ) - end if - end if ! test on UNSTABLE +do i = pdims%i_start, pdims%i_end + + grad_t_adj(i,j) = zero + grad_q_adj(i,j) = zero + if ( unstable(i,j) ) then + ! Here this is an estimate of the gradient adjustment applied + ! the previous timestep (assumes T1_SD has not changed much, + ! which in turn assumes the surface fluxes have not) + if (flux_grad == Locketal2000) then + grad_t_adj(i,j) = min( max_t_grad , & + a_grad_adj * t1_sd(i,j) / zh_prev(i,j) ) + grad_q_adj(i,j) = zero + else if (flux_grad == HoltBov1993) then + ! Use constants from Holtslag and Boville (1993) + ! Conv limit GAMMA_TH = 10 *FTL1/(wstar*zh) + ! Neut limit GAMMA_TH = 7.2*wstar*FTL1/(ustar^2*zh) + wstar3 = fb_surf(i,j) * zh_prev(i,j) + c_ws = 0.6_r_bl + w_m =( v_s(i,j)**3 + c_ws*wstar3 )**one_third + + grad_t_adj(i,j) = a_ga_hb93*(wstar3**one_third)*ftl(i,j,1) & + / ( rhostar_gb(i,j)*w_m*w_m*zh_prev(i,j) ) + ! GRAD_Q_ADJ(I,j) = A_GA_HB93*(WSTAR3**one_third)*FQW(I,j,1) + ! / ( RHOSTAR_GB(I,j)*W_M*W_M*ZH_PREV(I,j) ) + grad_q_adj(i,j) = zero + else if (flux_grad == LockWhelan2006) then + ! Use constants from LockWhelan2006 + ! Conv limit GAMMA_TH = 10 *FTL1/(wstar*zh) + ! Neut limit GAMMA_TH = 7.5*FTL1/(ustar*zh) + wstar3 = fb_surf(i,j) * zh_prev(i,j) + c_ws = 0.42_r_bl ! = 0.75^3 + pr_neut = 0.75_r_bl + w_h = ( ( v_s(i,j)**3+c_ws*wstar3 )**one_third )/ pr_neut + + grad_t_adj(i,j) = a_ga_lw06 * ftl(i,j,1) & + / ( rhostar_gb(i,j)*w_h*zh_prev(i,j) ) + grad_q_adj(i,j) = a_ga_lw06 * fqw(i,j,1) & + / ( rhostar_gb(i,j)*w_h*zh_prev(i,j) ) + end if + end if ! test on UNSTABLE - end do end do !$OMP end do +!end do ! Water tracers assume flux_grad == Locketal2000 if (l_wtrac) then do i_wt = 1, n_wtrac -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - wtrac_bl(i_wt)%grad_q_adj(i,j) = zero - end do +! !$OMP do SCHEDULE(STATIC) +! do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + wtrac_bl(i_wt)%grad_q_adj(i,j) = zero end do -!$OMP end do + !$OMP end do +! end do +! !$OMP end do end do end if @@ -2673,46 +2685,45 @@ subroutine kmkhz_9c ( & !$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 - if (k <= ntml_prev(i,j)) then - dsldz(i) = -grcp + grad_t_adj(i,j) - else - dsldz(i) = -grcp - end if - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + if (k <= ntml_prev(i,j)) then + dsldz(i) = -grcp + grad_t_adj(i,j) + else + dsldz(i) = -grcp + end if + end do !DIR$ NOFUSION !DIR$ VECTOR ALWAYS - do i = pdims%i_start, pdims%i_end - virt_factor = one + c_virtual*q(i,j,k) - qcl(i,j,k) - & - qcf(i,j,k) - - dqcldz(i,j,k) = -( dsldz(i)*dqsdt(i,j,k) & - + g*qs(i,j,k)/(r*t(i,j,k)*virt_factor) ) & - / ( one + lcrcp*dqsdt(i,j,k) ) - dqcfdz(i,j,k) = -( dsldz(i)*dqsdt(i,j,k) & - + g*qs(i,j,k)/(r*t(i,j,k)*virt_factor) ) * fgf & - / ( one + lsrcp*dqsdt(i,j,k) ) - end do + do i = pdims%i_start, pdims%i_end + virt_factor = one + c_virtual*q(i,j,k) - qcl(i,j,k) - & + qcf(i,j,k) + + dqcldz(i,j,k) = -( dsldz(i)*dqsdt(i,j,k) & + + g*qs(i,j,k)/(r*t(i,j,k)*virt_factor) ) & + / ( one + lcrcp*dqsdt(i,j,k) ) + dqcfdz(i,j,k) = -( dsldz(i)*dqsdt(i,j,k) & + + g*qs(i,j,k)/(r*t(i,j,k)*virt_factor) ) * fgf & + / ( one + lsrcp*dqsdt(i,j,k) ) + end do ! limit calculation to greater than a small cloud fraction !DIR$ NOFUSION - do i = pdims%i_start, pdims%i_end - if ( qcl(i,j,k) + qcf(i,j,k) > zero & - .and. cf(i,j,k) > 1.0e-3_r_bl ) then - cfl(i,j,k) = cf(i,j,k) * qcl(i,j,k) / & - ( qcl(i,j,k) + qcf(i,j,k) ) - cff(i,j,k) = cf(i,j,k) * qcf(i,j,k) / & - ( qcl(i,j,k) + qcf(i,j,k) ) - else - cfl(i,j,k) = zero - cff(i,j,k) = zero - end if + do i = pdims%i_start, pdims%i_end + if ( qcl(i,j,k) + qcf(i,j,k) > zero & + .and. cf(i,j,k) > 1.0e-3_r_bl ) then + cfl(i,j,k) = cf(i,j,k) * qcl(i,j,k) / & + ( qcl(i,j,k) + qcf(i,j,k) ) + cff(i,j,k) = cf(i,j,k) * qcf(i,j,k) / & + ( qcl(i,j,k) + qcf(i,j,k) ) + else + cfl(i,j,k) = zero + cff(i,j,k) = zero + end if - end do end do + !end do end do !$OMP end do @@ -2724,8 +2735,9 @@ subroutine kmkhz_9c ( & ! First the SML !--------------- -!$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end + +! do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end cloud_base(i,j)= .false. zc(i,j) = zero @@ -2742,8 +2754,9 @@ subroutine kmkhz_9c ( & if ( l_check_ntp1 ) cf_sml(i,j) = max( cf_sml(i,j), cf(i,j,k+1) ) end if end do -end do -!$OMP end do NOWAIT + !$OMP end do NOWAIT +! end do + !----------------------------------------------------------------------- ! First find cloud-base as seen by the cloud scheme @@ -2752,258 +2765,276 @@ subroutine kmkhz_9c ( & ! to use as first guess or lower limit !----------------------------------------------------------------------- + +! 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 - k_level(i,j) = ntml(i,j) - if ( cf_sml(i,j) > sc_cftol ) then - if ( .not. l_check_ntp1 ) k_level(i,j) = ntml(i,j)-1 - do while ( cf(i,j,max(k_level(i,j),1)) > sc_cftol & - .and. k_level(i,j) >= 2 ) - k_level(i,j) = k_level(i,j) - 1 - end do - end if - end do +do i = pdims%i_start, pdims%i_end + k_level(i,j) = ntml(i,j) + if ( cf_sml(i,j) > sc_cftol ) then + if ( .not. l_check_ntp1 ) k_level(i,j) = ntml(i,j)-1 + do while ( cf(i,j,max(k_level(i,j),1)) > sc_cftol & + .and. k_level(i,j) >= 2 ) + k_level(i,j) = k_level(i,j) - 1 + end do + end if end do !$OMP end do NOWAIT +! 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 - if ( cf_sml(i,j) > sc_cftol ) then - if ( k_level(i,j) == 1 .and. & - cf(i,j,max(k_level(i,j),1)) > sc_cftol) then - z_cf_base(i,j) = zero - else - z_cf_base(i,j) = z_uv(i,j,k_level(i,j)+1) - end if - zc(i,j) = z_ctop(i,j) - z_cf_base(i,j) +do i = pdims%i_start, pdims%i_end + if ( cf_sml(i,j) > sc_cftol ) then + if ( k_level(i,j) == 1 .and. & + cf(i,j,max(k_level(i,j),1)) > sc_cftol) then + z_cf_base(i,j) = zero + else + z_cf_base(i,j) = z_uv(i,j,k_level(i,j)+1) end if - end do + zc(i,j) = z_ctop(i,j) - z_cf_base(i,j) + end if end do !$OMP end do NOWAIT +! end do + !-------------------------------------------------- ! Find lowest level within ML with max CF !-------------------------------------------------- + +! 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 k = min(bl_levels, ntml(i,j)+p1), 1, -1 - if ( .not. cloud_base(i,j) .and. & - cf_sml(i,j) > sc_cftol ) then - ! within cloudy boundary layer - if ( k == 1) then - cloud_base(i,j) = .true. - else - if ( cf(i,j,k-1) < cf(i,j,k) ) cloud_base(i,j) = .true. - end if - k_cbase(i,j) = k +do i = pdims%i_start, pdims%i_end + do k = min(bl_levels, ntml(i,j)+p1), 1, -1 + if ( .not. cloud_base(i,j) .and. & + cf_sml(i,j) > sc_cftol ) then + ! within cloudy boundary layer + if ( k == 1) then + cloud_base(i,j) = .true. + else + if ( cf(i,j,k-1) < cf(i,j,k) ) cloud_base(i,j) = .true. end if - end do ! K - end do ! I -end do ! j + k_cbase(i,j) = k + end if + end do ! K +end do ! I !$OMP end do NOWAIT +! end do ! j + + +!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 - !-------------------------------------------------- - ! Use adiabatic qcl gradient to estimate cloud-base - ! from in-cloud qcl in level K_CBASE - ! If k_cbase = 0 then it hasn't been initialised - !-------------------------------------------------- + !-------------------------------------------------- + ! Use adiabatic qcl gradient to estimate cloud-base + ! from in-cloud qcl in level K_CBASE + ! If k_cbase = 0 then it hasn't been initialised + !-------------------------------------------------- - if ( cloud_base(i,j) .and. k_cbase(i,j) /= 0 ) then - z_cbase = z_tq(i,j,k_cbase(i,j)) - & - qcl(i,j,k_cbase(i,j)) / & - ( cf(i,j,k_cbase(i,j))*dqcldz(i,j,k_cbase(i,j)) ) - if ( dqcfdz(i,j,k_cbase(i,j)) > zero ) then - z_cbase = min( z_cbase, z_tq(i,j,k_cbase(i,j)) - & - qcf(i,j,k_cbase(i,j)) / & - ( cf(i,j,k_cbase(i,j))*dqcfdz(i,j,k_cbase(i,j)) ) & - ) - else - !--------------------------------------------------------- - ! No adiabatic QCF gradient so find lowest level, K_CFF, - ! with CFF>SC_CFTOL and assume cloud-base within that leve - !--------------------------------------------------------- - !Initialise K_CFF = lowest level with ice cloud - k_cff = k_cbase(i,j) - if (k_cff > 1) then - do while ( cff(i,j,k_cff) > sc_cftol & - .and. k_cff > 1 ) - k_cff = k_cff - 1 - end do - end if - if ( cff(i,j,k_cff) <= sc_cftol .and. & - k_cff < k_cbase(i,j) ) & - k_cff = k_cff + 1 - ! will want to raise K_CFF back up one level unless - ! level 1 is cloudy or no sig frozen cloud at all - z_cbase = min( z_cbase, z_top(i,j,k_cff) - & - dzl(i,j,k_cff) & - * cff(i,j,k_cff)/cf(i,j,k_cff) ) + if ( cloud_base(i,j) .and. k_cbase(i,j) /= 0 ) then + z_cbase = z_tq(i,j,k_cbase(i,j)) - & + qcl(i,j,k_cbase(i,j)) / & + ( cf(i,j,k_cbase(i,j))*dqcldz(i,j,k_cbase(i,j)) ) + if ( dqcfdz(i,j,k_cbase(i,j)) > zero ) then + z_cbase = min( z_cbase, z_tq(i,j,k_cbase(i,j)) - & + qcf(i,j,k_cbase(i,j)) / & + ( cf(i,j,k_cbase(i,j))*dqcfdz(i,j,k_cbase(i,j)) ) & + ) + else + !--------------------------------------------------------- + ! No adiabatic QCF gradient so find lowest level, K_CFF, + ! with CFF>SC_CFTOL and assume cloud-base within that leve + !--------------------------------------------------------- + !Initialise K_CFF = lowest level with ice cloud + k_cff = k_cbase(i,j) + if (k_cff > 1) then + do while ( cff(i,j,k_cff) > sc_cftol & + .and. k_cff > 1 ) + k_cff = k_cff - 1 + end do end if - !------------------------------------------------------ - ! use cloud-base as seen by cloud scheme as lower limit - ! and base of level NTML+1 as upper limit - !------------------------------------------------------ - z_cbase = min( z_uv(i,j,ntml(i,j)+1), & - max( z_cf_base(i,j), z_cbase) ) - - zc(i,j) = z_ctop(i,j) - z_cbase + if ( cff(i,j,k_cff) <= sc_cftol .and. & + k_cff < k_cbase(i,j) ) & + k_cff = k_cff + 1 + ! will want to raise K_CFF back up one level unless + ! level 1 is cloudy or no sig frozen cloud at all + z_cbase = min( z_cbase, z_top(i,j,k_cff) - & + dzl(i,j,k_cff) & + * cff(i,j,k_cff)/cf(i,j,k_cff) ) end if + !------------------------------------------------------ + ! use cloud-base as seen by cloud scheme as lower limit + ! and base of level NTML+1 as upper limit + !------------------------------------------------------ + z_cbase = min( z_uv(i,j,ntml(i,j)+1), & + max( z_cf_base(i,j), z_cbase) ) + + zc(i,j) = z_ctop(i,j) - z_cbase + end if - end do end do !$OMP end do NOWAIT +!end do + !----------------------------------------------------------------------- ! Second DSC layer !----------------------------------------------------------------------- -!$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - cloud_base(i,j) = .false. - zc_dsc(i,j) = zero - k_cbase(i,j) = 0 - z_cf_base(i,j) = zhsc(i,j) - z_ctop(i,j) = zhsc(i,j) - if ( dsc(i,j) ) then - k = ntdsc(i,j) - cf_dsc(i,j) = max( cf(i,j,k), cf(i,j,k-1) ) - if ( l_check_ntp1 ) cf_dsc(i,j) = max( cf_dsc(i,j), cf(i,j,k+1) ) - else - cf_dsc(i,j) = zero - end if - end do +!do j = pdims%j_start, pdims%j_end +!$OMP do SCHEDULE(STATIC) +do i = pdims%i_start, pdims%i_end + cloud_base(i,j) = .false. + zc_dsc(i,j) = zero + k_cbase(i,j) = 0 + z_cf_base(i,j) = zhsc(i,j) + z_ctop(i,j) = zhsc(i,j) + + if ( dsc(i,j) ) then + k = ntdsc(i,j) + cf_dsc(i,j) = max( cf(i,j,k), cf(i,j,k-1) ) + if ( l_check_ntp1 ) cf_dsc(i,j) = max( cf_dsc(i,j), cf(i,j,k+1) ) + else + cf_dsc(i,j) = zero + end if end do !$OMP end do NOWAIT +!end do + !------------------------------------------------------------- ! Find cloud-base as seen by cloud scheme, Z_CF_BASE, ! to use as first guess or lower limit and find cloud top. !------------------------------------------------------------- + +!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 - k_level(i,j) = ntdsc(i,j) - if ( cf_dsc(i,j) > sc_cftol ) then - ! assume level NTDSC is cloudy so start from NTDSC-1 - if ( .not. l_check_ntp1 ) k_level(i,j) = max( 2, ntdsc(i,j) - 1 ) - do while ( cf(i,j,k_level(i,j)) > sc_cftol & - .and. k_level(i,j) >= 2 ) - k_level(i,j) = k_level(i,j) - 1 - end do - end if - end do +do i = pdims%i_start, pdims%i_end + k_level(i,j) = ntdsc(i,j) + if ( cf_dsc(i,j) > sc_cftol ) then + ! assume level NTDSC is cloudy so start from NTDSC-1 + if ( .not. l_check_ntp1 ) k_level(i,j) = max( 2, ntdsc(i,j) - 1 ) + do while ( cf(i,j,k_level(i,j)) > sc_cftol & + .and. k_level(i,j) >= 2 ) + k_level(i,j) = k_level(i,j) - 1 + end do + end if end do !$OMP end do NOWAIT +!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 - if ( cf_dsc(i,j) > sc_cftol ) then - if ( k_level(i,j) == 1 .and. & - cf(i,j,max(k_level(i,j),1)) > sc_cftol) then - z_cf_base(i,j) = zero - else - z_cf_base(i,j) = z_uv(i,j,k_level(i,j)+1) - end if - zc_dsc(i,j) = z_ctop(i,j) - z_cf_base(i,j) ! first guess +do i = pdims%i_start, pdims%i_end + if ( cf_dsc(i,j) > sc_cftol ) then + if ( k_level(i,j) == 1 .and. & + cf(i,j,max(k_level(i,j),1)) > sc_cftol) then + z_cf_base(i,j) = zero + else + z_cf_base(i,j) = z_uv(i,j,k_level(i,j)+1) end if - end do + zc_dsc(i,j) = z_ctop(i,j) - z_cf_base(i,j) ! first guess + end if end do !$OMP end do NOWAIT +!end do + !-------------------------------------------------- ! Find lowest level within ML with max CF !-------------------------------------------------- + +!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 k = min(bl_levels,ntdsc(i,j)+p1), 1, -1 - if ( .not. cloud_base(i,j) .and. & - cf_dsc(i,j) > sc_cftol ) then - ! within cloudy boundary layer - if ( k == 1) then - cloud_base(i,j) = .true. - else - if ( cf(i,j,k-1) < cf(i,j,k) ) cloud_base(i,j) = .true. - end if - k_cbase(i,j) = k +do i = pdims%i_start, pdims%i_end + do k = min(bl_levels,ntdsc(i,j)+p1), 1, -1 + if ( .not. cloud_base(i,j) .and. & + cf_dsc(i,j) > sc_cftol ) then + ! within cloudy boundary layer + if ( k == 1) then + cloud_base(i,j) = .true. + else + if ( cf(i,j,k-1) < cf(i,j,k) ) cloud_base(i,j) = .true. end if - end do ! K - end do ! I -end do ! J + k_cbase(i,j) = k + end if + end do ! K +end do ! I !$OMP end do NOWAIT +!end do ! J -!$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - !-------------------------------------------------- - ! use adiabatic qcl gradient to estimate cloud-base - ! from in-cloud qcl in level K_CBASE - !-------------------------------------------------- - if ( cloud_base(i,j) .and. k_cbase(i,j) /= 0 ) then - z_cbase = z_tq(i,j,k_cbase(i,j)) - & - qcl(i,j,k_cbase(i,j)) / & - ( cf(i,j,k_cbase(i,j))*dqcldz(i,j,k_cbase(i,j)) ) - if ( dqcfdz(i,j,k_cbase(i,j)) > zero ) then - z_cbase = min( z_cbase, z_tq(i,j,k_cbase(i,j)) - & - qcf(i,j,k_cbase(i,j)) / & - ( cf(i,j,k_cbase(i,j))*dqcfdz(i,j,k_cbase(i,j)) ) & - ) - else - ! Initialise K_CFF - k_cff = k_cbase(i,j) - if (k_cff > 1) then - do while ( cff(i,j,k_cff) > sc_cftol & - .and. k_cff > 1) - k_cff = k_cff - 1 - end do - end if - !---------------------------------------------------------- - ! No adiabatic QCF gradient so find lowest level, K_CFF, - ! with CFF>SC_CFTOL and assume cloud-base within that level - !---------------------------------------------------------- - if ( cff(i,j,k_cff) <= sc_cftol .and. & - k_cff < k_cbase(i,j) ) & - k_cff = k_cff + 1 - ! will want to raise K_CFF back up one level unless - ! level 1 is cloudy or no sig frozen cloud at all - z_cbase = min( z_cbase, z_top(i,j,k_cff) - & - dzl(i,j,k_cff) & - * cff(i,j,k_cff)/cf(i,j,k_cff) ) - end if - !------------------------------------------------------ - ! use cloud-base as seen by cloud scheme as lower limit - ! and base of level NTDSC+1 as upper limit - !------------------------------------------------------ - z_cbase = min( z_uv(i,j,ntdsc(i,j)+1), & - max( z_cf_base(i,j) , z_cbase) ) - zc_dsc(i,j) = z_ctop(i,j) - z_cbase +! do j = pdims%j_start, pdims%j_end +!$OMP do SCHEDULE(STATIC) +do i = pdims%i_start, pdims%i_end + + !-------------------------------------------------- + ! use adiabatic qcl gradient to estimate cloud-base + ! from in-cloud qcl in level K_CBASE + !-------------------------------------------------- + if ( cloud_base(i,j) .and. k_cbase(i,j) /= 0 ) then + z_cbase = z_tq(i,j,k_cbase(i,j)) - & + qcl(i,j,k_cbase(i,j)) / & + ( cf(i,j,k_cbase(i,j))*dqcldz(i,j,k_cbase(i,j)) ) + if ( dqcfdz(i,j,k_cbase(i,j)) > zero ) then + z_cbase = min( z_cbase, z_tq(i,j,k_cbase(i,j)) - & + qcf(i,j,k_cbase(i,j)) / & + ( cf(i,j,k_cbase(i,j))*dqcfdz(i,j,k_cbase(i,j)) ) & + ) + else + ! Initialise K_CFF + k_cff = k_cbase(i,j) + if (k_cff > 1) then + do while ( cff(i,j,k_cff) > sc_cftol & + .and. k_cff > 1) + k_cff = k_cff - 1 + end do + end if + !---------------------------------------------------------- + ! No adiabatic QCF gradient so find lowest level, K_CFF, + ! with CFF>SC_CFTOL and assume cloud-base within that level + !---------------------------------------------------------- + if ( cff(i,j,k_cff) <= sc_cftol .and. & + k_cff < k_cbase(i,j) ) & + k_cff = k_cff + 1 + ! will want to raise K_CFF back up one level unless + ! level 1 is cloudy or no sig frozen cloud at all + z_cbase = min( z_cbase, z_top(i,j,k_cff) - & + dzl(i,j,k_cff) & + * cff(i,j,k_cff)/cf(i,j,k_cff) ) end if + !------------------------------------------------------ + ! use cloud-base as seen by cloud scheme as lower limit + ! and base of level NTDSC+1 as upper limit + !------------------------------------------------------ + z_cbase = min( z_uv(i,j,ntdsc(i,j)+1), & + max( z_cf_base(i,j) , z_cbase) ) + + zc_dsc(i,j) = z_ctop(i,j) - z_cbase + end if - !----------------------------------------------------------------- - ! Layer cloud depth cannot be > the layer depth itself. - !----------------------------------------------------------------- + !----------------------------------------------------------------- + ! Layer cloud depth cannot be > the layer depth itself. + !----------------------------------------------------------------- - zc_dsc(i,j) = min( zc_dsc(i,j), dscdepth(i,j) ) + zc_dsc(i,j) = min( zc_dsc(i,j), dscdepth(i,j) ) - end do !I -end do !J +end do !I !$OMP end do +! end do !J + !----------------------------------------------------------------------- ! 6. Calculate buoyancy flux factor used in the diagnosis of decoupling @@ -3020,62 +3051,62 @@ subroutine kmkhz_9c ( & !$OMP do SCHEDULE(DYNAMIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end + !do j = pdims%j_start, pdims%j_end - ! This is to help vectorization - do i = pdims%i_start, pdims%i_end - if ( k <= ntml_prev(i,j) .or. l_converge_ga ) then - ! If using more accurate treatment of gradient adjustment in the - ! buoyancy-flux integration, it needs to be set on all levels. - grad_t_adj_inv_rdz(i) = grad_t_adj(i,j)/rdz(i,j,k) - grad_q_adj_inv_rdz(i) = grad_q_adj(i,j)/rdz(i,j,k) - else - grad_t_adj_inv_rdz(i) = zero - grad_q_adj_inv_rdz(i) = zero - end if + ! This is to help vectorization + do i = pdims%i_start, pdims%i_end + if ( k <= ntml_prev(i,j) .or. l_converge_ga ) then + ! If using more accurate treatment of gradient adjustment in the + ! buoyancy-flux integration, it needs to be set on all levels. + grad_t_adj_inv_rdz(i) = grad_t_adj(i,j)/rdz(i,j,k) + grad_q_adj_inv_rdz(i) = grad_q_adj(i,j)/rdz(i,j,k) + else + grad_t_adj_inv_rdz(i) = zero + grad_q_adj_inv_rdz(i) = zero + end if - !---------------------------------------------------------- - ! CF_FOR_WB is uniform `bl' CF for use within cloud layers - !---------------------------------------------------------- - cf_for_wb(i) = zero - z_cbase = zh(i,j)-zc(i,j) - zdsc_cbase = zhsc(i,j)-zc_dsc(i,j) - if ( z_tq(i,j,k) <= zh(i,j) .and. & - z_tq(i,j,k) >= z_cbase) cf_for_wb(i) = cf_sml(i,j) - if ( z_tq(i,j,k) <= zhsc(i,j) .and. & - z_tq(i,j,k) >= zdsc_cbase) cf_for_wb(i) = cf_dsc(i,j) - end do + !---------------------------------------------------------- + ! CF_FOR_WB is uniform `bl' CF for use within cloud layers + !---------------------------------------------------------- + cf_for_wb(i) = zero + z_cbase = zh(i,j)-zc(i,j) + zdsc_cbase = zhsc(i,j)-zc_dsc(i,j) + if ( z_tq(i,j,k) <= zh(i,j) .and. & + z_tq(i,j,k) >= z_cbase) cf_for_wb(i) = cf_sml(i,j) + if ( z_tq(i,j,k) <= zhsc(i,j) .and. & + z_tq(i,j,k) >= zdsc_cbase) cf_for_wb(i) = cf_dsc(i,j) + end do !DIR$ NOFUSION !DIR$ VECTOR ALWAYS - do i = pdims%i_start, pdims%i_end - dqw = qw(i,j,k) - qw(i,j,k-1) - dsl = sl(i,j,k) - sl(i,j,k-1) - dsl_ga = dsl - grad_t_adj_inv_rdz(i) - dqw_ga = dqw - grad_q_adj_inv_rdz(i) + do i = pdims%i_start, pdims%i_end + dqw = qw(i,j,k) - qw(i,j,k-1) + dsl = sl(i,j,k) - sl(i,j,k-1) + dsl_ga = dsl - grad_t_adj_inv_rdz(i) + dqw_ga = dqw - grad_q_adj_inv_rdz(i) - !---------------------------------------------------------- - ! WB = -K_SURF*(DB/DZ - gamma_buoy) - K_TOP*DB/DZ - ! This is integrated in EXCF_NL, iterating the K profiles. - ! Here the relevant integrated DB/DZ factors are calculated - !---------------------------------------------------------- - db_ga_dry(i,j,k) = - g * & - ( btm(i,j,k-1)*dsl_ga + bqm(i,j,k-1)*dqw_ga ) - db_noga_dry(i,j,k) = - g * & - ( btm(i,j,k-1)*dsl + bqm(i,j,k-1)*dqw ) - db_ga_cld(i,j,k) = - g * & - ( btm_cld(i,j,k-1)*dsl_ga + bqm_cld(i,j,k-1)*dqw_ga ) - db_noga_cld(i,j,k) = - g * & - ( btm_cld(i,j,k-1)*dsl + bqm_cld(i,j,k-1)*dqw ) - !------------------------------------------------------- - ! Weight cloud layer factors with cloud fraction - !------------------------------------------------------- - db_ga_cld(i,j,k) = db_ga_dry(i,j,k)*(one-cf_for_wb(i)) + & - db_ga_cld(i,j,k)*cf_for_wb(i) - db_noga_cld(i,j,k) = db_noga_dry(i,j,k)*(one-cf_for_wb(i)) + & - db_noga_cld(i,j,k)*cf_for_wb(i) - end do + !---------------------------------------------------------- + ! WB = -K_SURF*(DB/DZ - gamma_buoy) - K_TOP*DB/DZ + ! This is integrated in EXCF_NL, iterating the K profiles. + ! Here the relevant integrated DB/DZ factors are calculated + !---------------------------------------------------------- + db_ga_dry(i,j,k) = - g * & + ( btm(i,j,k-1)*dsl_ga + bqm(i,j,k-1)*dqw_ga ) + db_noga_dry(i,j,k) = - g * & + ( btm(i,j,k-1)*dsl + bqm(i,j,k-1)*dqw ) + db_ga_cld(i,j,k) = - g * & + ( btm_cld(i,j,k-1)*dsl_ga + bqm_cld(i,j,k-1)*dqw_ga ) + db_noga_cld(i,j,k) = - g * & + ( btm_cld(i,j,k-1)*dsl + bqm_cld(i,j,k-1)*dqw ) + !------------------------------------------------------- + ! Weight cloud layer factors with cloud fraction + !------------------------------------------------------- + db_ga_cld(i,j,k) = db_ga_dry(i,j,k)*(one-cf_for_wb(i)) + & + db_ga_cld(i,j,k)*cf_for_wb(i) + db_noga_cld(i,j,k) = db_noga_dry(i,j,k)*(one-cf_for_wb(i)) + & + db_noga_cld(i,j,k)*cf_for_wb(i) end do + !end do end do !$OMP end do NOWAIT @@ -3083,60 +3114,65 @@ subroutine kmkhz_9c ( & ! 7. Calculate inputs for the top of b.l. entrainment parametrization !----------------------------------------------------------------------- + +!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 - zeta_r_dsc(i,j) = zero - chi_s_dsct(i,j) = zero - d_siems_dsc(i,j) = zero - br_fback_dsc(i,j)= zero - cld_factor_dsc(i,j) = zero - bt_dsct(i,j) = zero - btt_dsct(i,j) = zero - db_dsct(i,j) = zero - db_dsct_cld(i,j) = zero - chi_s_top(i,j) = zero - d_siems(i,j) = zero - br_fback(i,j)= zero - cld_factor(i,j) = zero - bt_top(i,j) = zero - btt_top(i,j) = zero - btc_top(i,j) = zero - db_top(i,j) = zero - db_top_cld(i,j) = zero ! default required if COUPLED - z_cld(i,j) = zero - z_cld_dsc(i,j) = zero - end do +do i = pdims%i_start, pdims%i_end + zeta_r_dsc(i,j) = zero + chi_s_dsct(i,j) = zero + d_siems_dsc(i,j) = zero + br_fback_dsc(i,j)= zero + cld_factor_dsc(i,j) = zero + bt_dsct(i,j) = zero + btt_dsct(i,j) = zero + db_dsct(i,j) = zero + db_dsct_cld(i,j) = zero + chi_s_top(i,j) = zero + d_siems(i,j) = zero + br_fback(i,j)= zero + cld_factor(i,j) = zero + bt_top(i,j) = zero + btt_top(i,j) = zero + btc_top(i,j) = zero + db_top(i,j) = zero + db_top_cld(i,j) = zero ! default required if COUPLED + z_cld(i,j) = zero + z_cld_dsc(i,j) = zero end do !$OMP end do +!end do + !----------------------------------------------------------------------- ! 7.1 Calculate surface buoyancy flux !----------------------------------------------------------------------- -!$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ! use mixed-layer average of buoyancy parameters - bflux_surf(i,j) = one_half * g * ( & - (btm(i,j,1)+btm(i,j,ntml(i,j)))*ftl(i,j,1) + & - (bqm(i,j,1)+bqm(i,j,ntml(i,j)))*fqw(i,j,1) ) - - if ( bflux_surf(i,j) > zero ) then - bflux_surf_sat(i,j) = one_half * g * ( & - (btm_cld(i,j,1)+btm_cld(i,j,ntml(i,j)))*ftl(i,j,1) + & - (bqm_cld(i,j,1)+bqm_cld(i,j,ntml(i,j)))*fqw(i,j,1) ) - if ( coupled(i,j) ) bflux_surf_sat(i,j) = one_half * g * ( & - (btm_cld(i,j,1)+btm_cld(i,j,ntdsc(i,j)))*ftl(i,j,1) + & - (bqm_cld(i,j,1)+bqm_cld(i,j,ntdsc(i,j)))*fqw(i,j,1) ) - else - bflux_surf_sat(i,j) = zero - end if +!do j = pdims%j_start, pdims%j_end +!$OMP do SCHEDULE(STATIC) +do i = pdims%i_start, pdims%i_end + + ! use mixed-layer average of buoyancy parameters + bflux_surf(i,j) = one_half * g * ( & + (btm(i,j,1)+btm(i,j,ntml(i,j)))*ftl(i,j,1) + & + (bqm(i,j,1)+bqm(i,j,ntml(i,j)))*fqw(i,j,1) ) + + if ( bflux_surf(i,j) > zero ) then + bflux_surf_sat(i,j) = one_half * g * ( & + (btm_cld(i,j,1)+btm_cld(i,j,ntml(i,j)))*ftl(i,j,1) + & + (bqm_cld(i,j,1)+bqm_cld(i,j,ntml(i,j)))*fqw(i,j,1) ) + if ( coupled(i,j) ) bflux_surf_sat(i,j) = one_half * g * ( & + (btm_cld(i,j,1)+btm_cld(i,j,ntdsc(i,j)))*ftl(i,j,1) + & + (bqm_cld(i,j,1)+bqm_cld(i,j,ntdsc(i,j)))*fqw(i,j,1) ) + else + bflux_surf_sat(i,j) = zero + end if - end do end do !$OMP end do NOWAIT +! No wait removal due to dimension change? +!end do + !----------------------------------------------------------------------- ! 7.2 Calculation of cloud fraction weighted thickness of @@ -3144,44 +3180,44 @@ subroutine kmkhz_9c ( & !----------------------------------------------------------------------- !$OMP do SCHEDULE(DYNAMIC) -do jj = pdims%j_start, pdims%j_end, bl_segment_size +do ii = pdims%i_start, pdims%i_end, bl_segment_size do k = 1, bl_levels - do j = jj, min((jj+bl_segment_size)-1,pdims%j_end) - do i = pdims%i_start, pdims%i_end - if ( k <= ntml(i,j)+1 ) then - z_cld(i,j) = z_cld(i,j) + & - cf(i,j,k) * one_half * dzl(i,j,k) + & - min( cfl(i,j,k) * one_half * dzl(i,j,k) , & - qcl(i,j,k) / dqcldz(i,j,k) ) - if ( dqcfdz(i,j,k) > zero) then - z_cld(i,j) = z_cld(i,j) + & - min( cff(i,j,k) * one_half * dzl(i,j,k) , & - qcf(i,j,k) / dqcfdz(i,j,k) ) - else - z_cld(i,j) = z_cld(i,j) + cff(i,j,k) * one_half * dzl(i,j,k) - end if + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) + !do i = pdims%i_start, pdims%i_end + if ( k <= ntml(i,j)+1 ) then + z_cld(i,j) = z_cld(i,j) + & + cf(i,j,k) * one_half * dzl(i,j,k) + & + min( cfl(i,j,k) * one_half * dzl(i,j,k) , & + qcl(i,j,k) / dqcldz(i,j,k) ) + if ( dqcfdz(i,j,k) > zero) then + z_cld(i,j) = z_cld(i,j) + & + min( cff(i,j,k) * one_half * dzl(i,j,k) , & + qcf(i,j,k) / dqcfdz(i,j,k) ) + else + z_cld(i,j) = z_cld(i,j) + cff(i,j,k) * one_half * dzl(i,j,k) end if + end if - if ( dsc(i,j) .and. k <= ntdsc(i,j)+1 .and. & - ( coupled(i,j) .or. & - z_top(i,j,k) >= zhsc(i,j)-zc_dsc(i,j) ) ) then - z_cld_dsc(i,j) = z_cld_dsc(i,j) + & - cf(i,j,k) * one_half * dzl(i,j,k) + & - min( cfl(i,j,k) * one_half * dzl(i,j,k) , & - qcl(i,j,k) / dqcldz(i,j,k) ) - if ( dqcfdz(i,j,k) > zero) then - z_cld_dsc(i,j) = z_cld_dsc(i,j) + & - min( cff(i,j,k) * one_half * dzl(i,j,k) , & - qcf(i,j,k) / dqcfdz(i,j,k) ) - else - z_cld_dsc(i,j) = z_cld_dsc(i,j) + & - cff(i,j,k) * one_half * dzl(i,j,k) - end if + if ( dsc(i,j) .and. k <= ntdsc(i,j)+1 .and. & + ( coupled(i,j) .or. & + z_top(i,j,k) >= zhsc(i,j)-zc_dsc(i,j) ) ) then + z_cld_dsc(i,j) = z_cld_dsc(i,j) + & + cf(i,j,k) * one_half * dzl(i,j,k) + & + min( cfl(i,j,k) * one_half * dzl(i,j,k) , & + qcl(i,j,k) / dqcldz(i,j,k) ) + if ( dqcfdz(i,j,k) > zero) then + z_cld_dsc(i,j) = z_cld_dsc(i,j) + & + min( cff(i,j,k) * one_half * dzl(i,j,k) , & + qcf(i,j,k) / dqcfdz(i,j,k) ) + else + z_cld_dsc(i,j) = z_cld_dsc(i,j) + & + cff(i,j,k) * one_half * dzl(i,j,k) end if - end do + end if + !end do end do end do -end do ! jj +end do ! ii !$OMP end do NOWAIT !----------------------------------------------------------------------- @@ -3191,9 +3227,11 @@ subroutine kmkhz_9c ( & !..assuming a discontinuous subgrid inversion structure. !---------------------------------------------------------------------- + +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(DYNAMIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) !-------------------------- ! First the SML !-------------------------- @@ -3215,14 +3253,14 @@ subroutine kmkhz_9c ( & if ( l_check_ntp1 .and. cf(i,j,km+1) > cf(i,j,kmax) ) kmax = km+1 cfl_ml = cf_sml(i,j)*cfl(i,j,kmax) & - /(cfl(i,j,kmax)+cff(i,j,kmax)+rbl_eps) + /(cfl(i,j,kmax)+cff(i,j,kmax)+rbl_eps) cff_ml = cf_sml(i,j)*cff(i,j,kmax) & - /(cfl(i,j,kmax)+cff(i,j,kmax)+rbl_eps) + /(cfl(i,j,kmax)+cff(i,j,kmax)+rbl_eps) if (cfl_ml > 0.01_r_bl) qcl_ic_top(i,j) = qcl(i,j,kmax)/cfl_ml & - + ( zh(i,j)-z_tq(i,j,km) )*dqcldz(i,j,km) + + ( zh(i,j)-z_tq(i,j,km) )*dqcldz(i,j,km) if (cff_ml > 0.01_r_bl) qcf_ic_top(i,j) = qcf(i,j,kmax)/cff_ml & - + ( zh(i,j)-z_tq(i,j,km) )*dqcfdz(i,j,km) + + ( zh(i,j)-z_tq(i,j,km) )*dqcfdz(i,j,km) end if dqw = dqw_sml(i,j) @@ -3232,8 +3270,8 @@ subroutine kmkhz_9c ( & dqcf = - cff_ml*qcf_ic_top(i,j) db_disc = g * ( btm(i,j,km)*dsl + bqm(i,j,km)*dqw + & - (lcrcp*btm(i,j,km) - etar*bqm(i,j,km)) * dqcl + & - (lsrcp*btm(i,j,km) - etar*bqm(i,j,km)) * dqcf ) + (lcrcp*btm(i,j,km) - etar*bqm(i,j,km)) * dqcl + & + (lsrcp*btm(i,j,km) - etar*bqm(i,j,km)) * dqcf ) if ( db_disc > 0.03_r_bl ) then ! Diagnosed inversion statically stable and at least ~1K @@ -3247,7 +3285,7 @@ subroutine kmkhz_9c ( & end if ! disc inversion diagnosed if ( sml_disc_inv(i,j) == 0 ) then - ! Calculate using simple grid-level differences + ! Calculate using simple grid-level differences kp = km+1 dqw = qw(i,j,kp) - qw(i,j,km) dsl = sl(i,j,kp) - sl(i,j,km) @@ -3277,17 +3315,22 @@ subroutine kmkhz_9c ( & chi_s_top(i,j) = zero ! term to zero end if - end do -end do + end do !i +end do !ii !$OMP end do +!end do + !-------------------------- ! Then the DSC layer !-------------------------- + +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(DYNAMIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) + ! do i = pdims%i_start, pdims%i_end if ( dsc(i,j) ) then qcl_ic_top(i,j) = zero @@ -3297,11 +3340,11 @@ subroutine kmkhz_9c ( & km = ntdsc(i,j) if ( dsc_disc_inv(i,j) == 1 ) then - !----------------------------------------------------- - ! Extrapolate water contents from level with max CF, - ! out of NTDSC and NTDSC-1 (ie near top of DSC), - ! to the top of the mixed layer - !----------------------------------------------------- + !----------------------------------------------------- + ! Extrapolate water contents from level with max CF, + ! out of NTDSC and NTDSC-1 (ie near top of DSC), + ! to the top of the mixed layer + !----------------------------------------------------- if (cf_dsc(i,j) > zero) then kmax = km if (cf(i,j,km-1) > cf(i,j,km)) kmax = km-1 @@ -3320,7 +3363,7 @@ subroutine kmkhz_9c ( & dqw = dqw_dsc(i,j) dsl = dsl_dsc(i,j) - ! ignore any cloud above the inversion + ! ignore any cloud above the inversion dqcl = - cfl_ml*qcl_ic_top(i,j) dqcf = - cff_ml*qcf_ic_top(i,j) @@ -3340,7 +3383,7 @@ subroutine kmkhz_9c ( & end if ! disc inversion diagnosed if ( dsc_disc_inv(i,j) == 0 ) then - ! Calculate using simple grid-level differences + ! Calculate using simple grid-level differences kp = km+1 dqw = qw(i,j,kp) - qw(i,j,km) dsl = sl(i,j,kp) - sl(i,j,km) @@ -3354,7 +3397,7 @@ subroutine kmkhz_9c ( & end if ! no disc inversion diagnosed db_dsct_cld(i,j) = g * ( btm_cld(i,j,km)*dsl & - + bqm_cld(i,j,km)*dqw ) + + bqm_cld(i,j,km)*dqw ) denom = a_qsm(i,j,km)*dqw - a_dqsdtm(i,j,km)*dsl if (abs(denom) > rbl_eps) then chi_s_dsct(i,j) = -qcl_ic_top(i,j) / denom @@ -3362,9 +3405,9 @@ subroutine kmkhz_9c ( & end if if ( db_dsct(i,j) < 0.003_r_bl ) then - ! Diagnosed inversion statically unstable: - ! ensure DB>0 so that entrainment is non-zero and - ! instability can be removed. + ! Diagnosed inversion statically unstable: + ! ensure DB>0 so that entrainment is non-zero and + ! instability can be removed. db_dsct(i,j) = 0.003_r_bl db_dsct_cld(i,j) = zero ! set buoyancy reversal chi_s_dsct(i,j) = zero ! term to zero @@ -3399,14 +3442,18 @@ subroutine kmkhz_9c ( & end do end do !$OMP end do +!end do + !----------------------------------------------------------------------- ! 7.4 Next those terms which depend on the presence of buoyancy reversal !----------------------------------------------------------------------- + +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(DYNAMIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) z_cld(i,j) = min( z_cld(i,j), zh(i,j) ) z_cld_dsc(i,j) = min( z_cld_dsc(i,j), zhsc(i,j) ) !--------------------------------------------------------------- @@ -3434,7 +3481,7 @@ subroutine kmkhz_9c ( & !---------------------------- db_top_cld(i,j) = -db_top_cld(i,j) * cld_factor(i,j) d_siems(i,j) = max( zero, & - chi_s_top(i,j) * db_top_cld(i,j) / (db_top(i,j)+rbl_eps) ) + chi_s_top(i,j) * db_top_cld(i,j) / (db_top(i,j)+rbl_eps) ) ! Linear feedback dependence for D<0.1 br_fback(i,j)= min( one, 10.0_r_bl*d_siems(i,j) ) zeta_r(i,j) = zeta_r(i,j) + (one-zeta_r(i,j))*br_fback(i,j) @@ -3468,9 +3515,9 @@ subroutine kmkhz_9c ( & br_fback_dsc(i,j) = min( one, 10.0_r_bl*d_siems_dsc(i,j) ) if ( entr_enhance_by_cu == Buoyrev_feedback & - .and. cumulus(i,j) & - .and. d_siems_dsc(i,j) < 0.1_r_bl & - .and. d_siems_dsc(i,j) > rbl_eps ) then + .and. cumulus(i,j) & + .and. d_siems_dsc(i,j) < 0.1_r_bl & + .and. d_siems_dsc(i,j) > rbl_eps ) then ! Assume mixing from cumulus can enhance the ! buoyancy reversal feedback in regime 01 for Cu>1000. br_fback_dsc(i,j) = cu_depth_fac + & - (one-cu_depth_fac)*br_fback_dsc(i,j) + (one-cu_depth_fac)*br_fback_dsc(i,j) end if zeta_r_dsc(i,j) = zeta_r_dsc(i,j) + & - (one-zeta_r_dsc(i,j))*br_fback_dsc(i,j) + (one-zeta_r_dsc(i,j))*br_fback_dsc(i,j) end if end if - end do -end do + end do !i +end do !ii !$OMP end do +! end do + !----------------------------------------------------------------------- ! 8. Calculate the radiative flux change across cloud top for mixed- @@ -3498,15 +3547,18 @@ subroutine kmkhz_9c ( & ! Initialise variables !------------------------------ + +!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 - k_cloud_top(i,j) = 0 - df_top_over_cp(i,j) = zero - df_inv_sml(i,j) = zero - end do +do i = pdims%i_start, pdims%i_end + + k_cloud_top(i,j) = 0 + df_top_over_cp(i,j) = zero + df_inv_sml(i,j) = zero end do !$OMP end do +! end do + if (l_new_kcloudtop) then !--------------------------------------------------------------------- @@ -3520,9 +3572,11 @@ subroutine kmkhz_9c ( & ! First find the level with maximum LW cooling, below z_rad_lim and ! in the upper half of the BL (ie, restrict search to `close' to ZH) !--------------------------------------------------------------------- -!$OMP do SCHEDULE(DYNAMIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end + + ! do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(DYNAMIC) + do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) k_rad_lim = ntml(i,j)+1 z_rad_lim = max( z_tq(i,j,k_rad_lim)+0.1_r_bl, 1.2_r_bl*zh(i,j) ) k = 1 @@ -3548,12 +3602,16 @@ subroutine kmkhz_9c ( & k_cloud_top(i,j) = k-1 end if end do ! i - !----------------------------------------------------------------- - ! Find bottom grid-level (K_LEVEL) for cloud-top radiative fux - ! divergence: higher of base of LW radiatively cooled layer, - ! 0.5*ZH, since cooling must be in upper part of layer - !----------------------------------------------------------------- - do i = pdims%i_start, pdims%i_end + end do ! ii + !$OMP end do + !----------------------------------------------------------------- + ! Find bottom grid-level (K_LEVEL) for cloud-top radiative fux + ! divergence: higher of base of LW radiatively cooled layer, + ! 0.5*ZH, since cooling must be in upper part of layer + !----------------------------------------------------------------- + !$OMP do SCHEDULE(DYNAMIC) + do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) k_level(i,j) = k_cloud_top(i,j) if ( k_cloud_top(i,j) > 1 ) then k_rad_lim = 1 @@ -3569,17 +3627,21 @@ subroutine kmkhz_9c ( & end do ! k end if end do ! i - end do ! j -!$OMP end do + end do ! ii + !$OMP end do + !end do ! j + else ! original method of finding the k_cloud_dsct, the top of the mixed layer ! as seen by radiation, found to be resolution dependent -!$OMP do SCHEDULE(DYNAMIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end + + !do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(DYNAMIC) + do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) ! restrict search to `close' to ZH k_rad_lim = ntml(i,j)+1 do k = 1, min(bl_levels,k_rad_lim) @@ -3603,13 +3665,17 @@ subroutine kmkhz_9c ( & end do ! k end do ! i + end do ! ii + !$OMP end do - !----------------------------------------------------------------- - ! Find bottom grid-level (K_LEVEL) for cloud-top radiative fux - ! divergence: higher of base of LW radiatively cooled layer, - ! 0.5*ZH, since cooling must be in upper part of layer - !----------------------------------------------------------------- - do i = pdims%i_start, pdims%i_end + !----------------------------------------------------------------- + ! Find bottom grid-level (K_LEVEL) for cloud-top radiative fux + ! divergence: higher of base of LW radiatively cooled layer, + ! 0.5*ZH, since cooling must be in upper part of layer + !----------------------------------------------------------------- + !$OMP do SCHEDULE(DYNAMIC) + do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) k_level(i,j) = k_cloud_top(i,j) if ( k_cloud_top(i,j) > 1 ) then k_rad_lim = 1 @@ -3625,8 +3691,10 @@ subroutine kmkhz_9c ( & end do end if end do ! i - end do ! j -!$OMP end do + end do ! ii + !$OMP end do +!end do ! j + end if ! test on l_new_kcloudtop @@ -3639,9 +3707,12 @@ subroutine kmkhz_9c ( & ! `clear-air' part from the grid-level divergence. !----------------------------------------------------------------- + +!do j = pdims%j_start, pdims%j_end !$OMP do SCHEDULE(DYNAMIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min((ii+bl_segment_size)-1,pdims%i_end) +! do i = pdims%i_start, pdims%i_end if ( k_cloud_top(i,j) > 0 ) then dflw_inv = zero @@ -3650,10 +3721,10 @@ subroutine kmkhz_9c ( & k = k_cloud_top(i,j)+1 if ( k < bl_levels ) then dflw_inv = dflw_over_cp(i,j,k) & - - dflw_over_cp(i,j,k+1) & + - dflw_over_cp(i,j,k+1) & * dzl(i,j,k)/dzl(i,j,k+1) dfsw_inv = dfsw_over_cp(i,j,k) & - - dfsw_over_cp(i,j,k+1) & + - dfsw_over_cp(i,j,k+1) & * dzl(i,j,k)/dzl(i,j,k+1) else dflw_inv = dflw_over_cp(i,j,k) @@ -3665,12 +3736,12 @@ subroutine kmkhz_9c ( & df_inv_sml(i,j) = dflw_inv + dfsw_inv df_top_over_cp(i,j) = frad_lw(i,j,k_cloud_top(i,j)+1) & - - frad_lw(i,j,k_level(i,j)) & - + dflw_inv + - frad_lw(i,j,k_level(i,j)) & + + dflw_inv dfsw_top = frad_sw(i,j,k_cloud_top(i,j)+1) & - - frad_sw(i,j,k_level(i,j)) & - + dfsw_inv + - frad_sw(i,j,k_level(i,j)) & + + dfsw_inv !----------------------------------------------------------- ! Combine SW and LW cloud-top divergences into a net @@ -3682,9 +3753,11 @@ subroutine kmkhz_9c ( & df_top_over_cp(i,j) = max( zero, & df_top_over_cp(i,j) + dfsw_frac * dfsw_top ) end if - end do -end do + end do !i +end do !ii !$OMP end do +!end do + ! ------------------------------------------------------------------ ! 9. Calculate the non-turbulent fluxes at the layer boundaries. @@ -3703,51 +3776,55 @@ subroutine kmkhz_9c ( & ! K_CLOUD_TOP and K_CLOUD_DSCT to top level in mixed layer ! ------------------------------------------------------------------ + +! 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 ( k_cloud_top(i,j) == 0 ) k_cloud_top(i,j) = ntml(i,j) - - ft_nt_zh(i,j) = frad(i,j,k_cloud_top(i,j)+1) & - + df_inv_sml(i,j) - ft_nt_zh(i,j) = ft_nt_zh(i,j) + fmic(i,j,ntml(i,j)+2,1) & - + fsubs(i,j,ntml(i,j),1) - fq_nt_zh(i,j) = fmic(i,j,ntml(i,j)+2,2) & - + fsubs(i,j,ntml(i,j),2) - - ft_nt_zhsc(i,j) = zero - ft_nt_zhsc(i,j) = zero - fq_nt_zhsc(i,j) = zero - if ( dsc(i,j) ) then - if ( k_cloud_dsct(i,j) == 0 ) k_cloud_dsct(i,j) = ntdsc(i,j) - ft_nt_zhsc(i,j) = frad(i,j,k_cloud_dsct(i,j)+1) & - + df_inv_dsc(i,j) - ft_nt_zhsc(i,j) = ft_nt_zhsc(i,j) + fmic(i,j,ntdsc(i,j)+2,1) & - + fsubs(i,j,ntdsc(i,j),1) - fq_nt_zhsc(i,j) = fmic(i,j,ntdsc(i,j)+2,2) & - + fsubs(i,j,ntdsc(i,j),2) - end if - end do +do i = pdims%i_start, pdims%i_end + if ( k_cloud_top(i,j) == 0 ) k_cloud_top(i,j) = ntml(i,j) + + ft_nt_zh(i,j) = frad(i,j,k_cloud_top(i,j)+1) & + + df_inv_sml(i,j) + ft_nt_zh(i,j) = ft_nt_zh(i,j) + fmic(i,j,ntml(i,j)+2,1) & + + fsubs(i,j,ntml(i,j),1) + fq_nt_zh(i,j) = fmic(i,j,ntml(i,j)+2,2) & + + fsubs(i,j,ntml(i,j),2) + + ft_nt_zhsc(i,j) = zero + ft_nt_zhsc(i,j) = zero + fq_nt_zhsc(i,j) = zero + if ( dsc(i,j) ) then + if ( k_cloud_dsct(i,j) == 0 ) k_cloud_dsct(i,j) = ntdsc(i,j) + ft_nt_zhsc(i,j) = frad(i,j,k_cloud_dsct(i,j)+1) & + + df_inv_dsc(i,j) + ft_nt_zhsc(i,j) = ft_nt_zhsc(i,j) + fmic(i,j,ntdsc(i,j)+2,1) & + + fsubs(i,j,ntdsc(i,j),1) + fq_nt_zhsc(i,j) = fmic(i,j,ntdsc(i,j)+2,2) & + + fsubs(i,j,ntdsc(i,j),2) + end if end do !$OMP end do NOWAIT +! end do + ! Repeat for water tracers if (l_wtrac) then do i_wt = 1, n_wtrac -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - fq_nt_zh_wtrac(i,j,i_wt) = fmic_wtrac(i,j,ntml(i,j)+2,i_wt) & - + fsubs_wtrac(i,j,ntml(i,j),i_wt) - fq_nt_zhsc_wtrac(i,j,i_wt) = zero - if ( dsc(i,j) ) then - fq_nt_zhsc_wtrac(i,j,i_wt) = fmic_wtrac(i,j,ntdsc(i,j)+2,i_wt) & - + fsubs_wtrac(i,j,ntdsc(i,j),i_wt) - end if - end do + ! do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + fq_nt_zh_wtrac(i,j,i_wt) = fmic_wtrac(i,j,ntml(i,j)+2,i_wt) & + + fsubs_wtrac(i,j,ntml(i,j),i_wt) + + fq_nt_zhsc_wtrac(i,j,i_wt) = zero + if ( dsc(i,j) ) then + fq_nt_zhsc_wtrac(i,j,i_wt) = fmic_wtrac(i,j,ntdsc(i,j)+2,i_wt) & + + fsubs_wtrac(i,j,ntdsc(i,j),i_wt) + end if end do -!$OMP end do NOWAIT + !$OMP end do NOWAIT + ! end do + end do end if @@ -3755,15 +3832,17 @@ subroutine kmkhz_9c ( & ! 10. Subroutine EXCF_NL. !----------------------------------------------------------------------- + +!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_save(i,j) = ntml(i,j) ! needed to identify changes - dsc_save(i,j) = dsc(i,j) ! in excf_nl - dsc_removed(i,j) = 0 - end do +do i = pdims%i_start, pdims%i_end + ntml_save(i,j) = ntml(i,j) ! needed to identify changes + dsc_save(i,j) = dsc(i,j) ! in excf_nl + dsc_removed(i,j) = 0 end do !$OMP end do NOWAIT +!end do + !$OMP end PARALLEL @@ -3790,70 +3869,124 @@ subroutine kmkhz_9c ( & rhof2,rhofsc,f_ngstress,tke_nl,zdsc_base,nbdsc & ) +!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL DEFAULT(SHARED) & -!$OMP private (i, j, i_wt, k, kl, kp, c_ws, c_tke, w_m, tothf_efl, totqf_efl, & +!$OMP private (i, i_wt, k, kl, kp, c_ws, c_tke, w_m, tothf_efl, totqf_efl, & !$OMP ml_tend, fa_tend, inv_tend, Prandtl, svl_lapse_rho, & !$OMP recip_svl_lapse, rhok_inv, svl_lapse, svl_target, svl_flux) - !----------------------------------------------------------------------- !-adjust SML/DSC properties depending on diagnoses in EXCF_NL ! Note that the non-turbulent fluxes at inversions will have been ! swapped in EXCF_NL (ie. FT/Q_NT_ZH/ZHSC) !----------------------------------------------------------------------- !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( dsc(i,j) .and. .not. dsc_save(i,j) ) then - !..decoupling diagnosed in EXCF_NL - change parameters around - dsl_dsc(i,j) = dsl_sml(i,j) - dqw_dsc(i,j) = dqw_sml(i,j) - db_dsct(i,j) = db_top(i,j) ! copy diagnostics across - zc_dsc(i,j) = zc(i,j) ! - dsc_disc_inv(i,j) = sml_disc_inv(i,j) - sml_disc_inv(i,j) = 0 - dsl_sml(i,j) = zero - dqw_sml(i,j) = zero - df_inv_dsc(i,j) = df_inv_sml(i,j) - df_inv_sml(i,j) = zero - res_inv(i,j) = 0 ! not clear what to do at dsc top (dzh diagnosed - ! assuming well-mixed) so do nothing for now! - dzh(i,j) = rmdi ! don't want to associate dzh with new zh - end if - end do +do i = pdims%i_start, pdims%i_end + if ( dsc(i,j) .and. .not. dsc_save(i,j) ) then + !..decoupling diagnosed in EXCF_NL - change parameters around + dsl_dsc(i,j) = dsl_sml(i,j) + dqw_dsc(i,j) = dqw_sml(i,j) + db_dsct(i,j) = db_top(i,j) ! copy diagnostics across + zc_dsc(i,j) = zc(i,j) ! + dsc_disc_inv(i,j) = sml_disc_inv(i,j) + sml_disc_inv(i,j) = 0 + dsl_sml(i,j) = zero + dqw_sml(i,j) = zero + df_inv_dsc(i,j) = df_inv_sml(i,j) + df_inv_sml(i,j) = zero + res_inv(i,j) = 0 ! not clear what to do at dsc top (dzh diagnosed + ! assuming well-mixed) so do nothing for now! + dzh(i,j) = rmdi ! don't want to associate dzh with new zh + end if end do !$OMP end do +!end do + ! Repeat for water tracers if (l_wtrac) then do i_wt = 1, n_wtrac -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( dsc(i,j) .and. .not. dsc_save(i,j) ) then - !..decoupling diagnosed in EXCF_NL - change parameters around - dqw_dsc_wtrac(i,j,i_wt) = dqw_sml_wtrac(i,j,i_wt) - dqw_sml_wtrac(i,j,i_wt) = zero - end if - end do + ! do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + if ( dsc(i,j) .and. .not. dsc_save(i,j) ) then + !..decoupling diagnosed in EXCF_NL - change parameters around + dqw_dsc_wtrac(i,j,i_wt) = dqw_sml_wtrac(i,j,i_wt) + dqw_sml_wtrac(i,j,i_wt) = zero + end if end do -!$OMP end do + !$OMP end do + ! end do end do end if if ( l_use_sml_dsc_fixes ) then -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( dsc_removed(i,j) == 1 ) then - ! decoupled layer removed in EXCF_NL because it had no - ! turbulence forcing + + ! do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + if ( dsc_removed(i,j) == 1 ) then + ! decoupled layer removed in EXCF_NL because it had no + ! turbulence forcing + dsl_dsc(i,j) = zero + dqw_dsc(i,j) = zero + dsc_disc_inv(i,j) = 0 + df_inv_dsc(i,j) = zero + else if ( dsc_removed(i,j) == 2 ) then + ! decoupled layer recoupled with surface layer + dsl_sml(i,j) = dsl_dsc(i,j) + dqw_sml(i,j) = dqw_dsc(i,j) + dsl_dsc(i,j) = zero + dqw_dsc(i,j) = zero + sml_disc_inv(i,j) = dsc_disc_inv(i,j) + dsc_disc_inv(i,j) = 0 + df_inv_sml(i,j) = df_inv_dsc(i,j) + df_inv_dsc(i,j) = zero + res_inv(i,j) = 0 ! dzh was diagnosed for a capping inversion + dzh(i,j) = rmdi ! subsumed into mixed layer so remove its flags + end if + end do + !$OMP end do + ! end do + + + ! Repeat for water tracers + if (l_wtrac) then + do i_wt = 1, n_wtrac + !do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + if ( dsc_removed(i,j) == 1 ) then + ! decoupled layer removed in EXCF_NL because it had no + ! turbulence forcing + dqw_dsc_wtrac(i,j,i_wt) = zero + else if ( dsc_removed(i,j) == 2 ) then + ! decoupled layer recoupled with surface layer + dqw_sml_wtrac(i,j,i_wt) = dqw_dsc_wtrac(i,j,i_wt) + dqw_dsc_wtrac(i,j,i_wt) = zero + end if + end do + !$OMP end do + !end do + end do + end if + +else ! not l_use_sml_dsc_fixes + + + !do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + if ( .not. dsc(i,j) .and. dsc_save(i,j) ) then + !..decoupled layer removed in EXCF_NL; either... + if ( ntml_save(i,j) == ntml(i,j) ) then + !...had no turbulence forcing dsl_dsc(i,j) = zero dqw_dsc(i,j) = zero dsc_disc_inv(i,j) = 0 df_inv_dsc(i,j) = zero - else if ( dsc_removed(i,j) == 2 ) then - ! decoupled layer recoupled with surface layer + else + !...recoupled with surface layer dsl_sml(i,j) = dsl_dsc(i,j) dqw_sml(i,j) = dqw_dsc(i,j) dsl_dsc(i,j) = zero @@ -3865,82 +3998,34 @@ subroutine kmkhz_9c ( & res_inv(i,j) = 0 ! dzh was diagnosed for a capping inversion dzh(i,j) = rmdi ! subsumed into mixed layer so remove its flags end if - end do + end if end do -!$OMP end do + !$OMP end do + !end do + ! Repeat for water tracers if (l_wtrac) then do i_wt = 1, n_wtrac -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( dsc_removed(i,j) == 1 ) then - ! decoupled layer removed in EXCF_NL because it had no - ! turbulence forcing + + !do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + if ( .not. dsc(i,j) .and. dsc_save(i,j) ) then + !..decoupled layer removed in EXCF_NL; either... + if ( ntml_save(i,j) == ntml(i,j) ) then + !...had no turbulence forcing dqw_dsc_wtrac(i,j,i_wt) = zero - else if ( dsc_removed(i,j) == 2 ) then - ! decoupled layer recoupled with surface layer + else + !...recoupled with surface layer dqw_sml_wtrac(i,j,i_wt) = dqw_dsc_wtrac(i,j,i_wt) dqw_dsc_wtrac(i,j,i_wt) = zero end if - end do - end do -!$OMP end do - end do - end if - -else ! not l_use_sml_dsc_fixes - -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( .not. dsc(i,j) .and. dsc_save(i,j) ) then - !..decoupled layer removed in EXCF_NL; either... - if ( ntml_save(i,j) == ntml(i,j) ) then - !...had no turbulence forcing - dsl_dsc(i,j) = zero - dqw_dsc(i,j) = zero - dsc_disc_inv(i,j) = 0 - df_inv_dsc(i,j) = zero - else - !...recoupled with surface layer - dsl_sml(i,j) = dsl_dsc(i,j) - dqw_sml(i,j) = dqw_dsc(i,j) - dsl_dsc(i,j) = zero - dqw_dsc(i,j) = zero - sml_disc_inv(i,j) = dsc_disc_inv(i,j) - dsc_disc_inv(i,j) = 0 - df_inv_sml(i,j) = df_inv_dsc(i,j) - df_inv_dsc(i,j) = zero - res_inv(i,j) = 0 ! dzh was diagnosed for a capping inversion - dzh(i,j) = rmdi ! subsumed into mixed layer so remove its flags end if - end if - end do - end do -!$OMP end do - - ! Repeat for water tracers - if (l_wtrac) then - do i_wt = 1, n_wtrac -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( .not. dsc(i,j) .and. dsc_save(i,j) ) then - !..decoupled layer removed in EXCF_NL; either... - if ( ntml_save(i,j) == ntml(i,j) ) then - !...had no turbulence forcing - dqw_dsc_wtrac(i,j,i_wt) = zero - else - !...recoupled with surface layer - dqw_sml_wtrac(i,j,i_wt) = dqw_dsc_wtrac(i,j,i_wt) - dqw_dsc_wtrac(i,j,i_wt) = zero - end if - end if - end do end do -!$OMP end do + !$OMP end do + !end do + end do end if @@ -3952,22 +4037,24 @@ subroutine kmkhz_9c ( & ! Calculate the non-turbulent fluxes at the DSC base ! ------------------------------------------------------------------ + +! 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 - ft_nt_dscb(i,j) = ft_nt(i,j,1) - if ( nbdsc(i,j) > 1 ) then - k = nbdsc(i,j) ! NBDSC marks the lowest flux-level - ! within the DSC layer - ! Interpolate non-turb flux to base - ! of DSC layer: - ft_nt_dscb(i,j) = ft_nt(i,j,k-1) + & - (ft_nt(i,j,k)-ft_nt(i,j,k-1)) & - *(zdsc_base(i,j)-z_uv(i,j,k-1))/dzl(i,j,k-1) - end if - end do +do i = pdims%i_start, pdims%i_end + ft_nt_dscb(i,j) = ft_nt(i,j,1) + if ( nbdsc(i,j) > 1 ) then + k = nbdsc(i,j) ! NBDSC marks the lowest flux-level + ! within the DSC layer + ! Interpolate non-turb flux to base + ! of DSC layer: + ft_nt_dscb(i,j) = ft_nt(i,j,k-1) + & + (ft_nt(i,j,k)-ft_nt(i,j,k-1)) & + *(zdsc_base(i,j)-z_uv(i,j,k-1))/dzl(i,j,k-1) + end if end do !$OMP end do NOWAIT +! end do + !----------------------------------------------------------------------- !..Specify entrainment fluxes at NTML+1 and NTDSC+1 directly through FTL @@ -3993,103 +4080,105 @@ subroutine kmkhz_9c ( & !..ie. the inversion is well-defined) !----------------------------------------------------------------------- -!$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - - zh_np1(i,j) = zero - t_frac(i,j) = zero - zrzi(i,j) = zero - we_rho(i,j) = zero - tothf_zh(i,j) = zero - totqf_zh(i,j) = zero - k=ntml(i,j)+1 - ! Only RHOKH_ENT is passed out of EXCFNL so recalculate WE: - we_parm(i,j) = rdz(i,j,k)* & - ( rhokh_top_ent(i,j)+rhokh_surf_ent(i,j) ) & - / rho_mix(i,j,k) - - if ( sml_disc_inv(i,j) == 1 .and. .not. coupled(i,j) .and. & - (rhokh_top_ent(i,j)+rhokh_surf_ent(i,j)) > zero ) then +! do j = pdims%j_start, pdims%j_end +!$OMP do SCHEDULE(STATIC) +do i = pdims%i_start, pdims%i_end - !----------------------------------------------------------------- - !..Calculate ZH at end of timestep, ZH_NP1 - !----------------------------------------------------------------- - !..linearly interpolate vertical velocity to ZH - if ( zh(i,j) >= z_tq(i,j,k) ) then - w_ls(i,j) = w(i,j,k) + ( w(i,j,k+1) - w(i,j,k) ) & - * (zh(i,j)-z_tq(i,j,k)) * rdz(i,j,k+1) - else - w_ls(i,j) = w(i,j,k) + ( w(i,j,k) - w(i,j,k-1) ) & - * (zh(i,j)-z_tq(i,j,k)) * rdz(i,j,k) - end if - w_ls(i,j) = min ( w_ls(i,j), zero ) - ! only interested in subsidence - - zh_np1(i,j) = zh(i,j) + & - timestep * ( we_parm(i,j) + w_ls(i,j) ) - zh_np1(i,j) = max( zh_np1(i,j), z_uv(i,j,k-1) ) - if ( zh_np1(i,j) > z_top(i,j,k+1) ) then - ! limit ZH and W_e (and therefore the entraiment fluxes) - ! because the inversion cannot rise more than one level - ! in a timestep. - zh_np1(i,j) = z_top(i,j,k+1) - we_parm(i,j) = & - (z_top(i,j,k+1) - zh(i,j))/timestep - w_ls(i,j) - end if - !----------------------------------------------------------------- - !..Decide on which grid-level to apply entrainment flux - !----------------------------------------------------------------- - if ( zh_np1(i,j) > z_uv(i,j,ntml(i,j)+2) ) then - ! ZH risen above level K+1 so specify appropriate flux - ! at this level and raise NTML by one (this means - ! gradient-adjustment is also applied at half-level - ! old_NTML+1). Note KH profiles should already be - ! calculated at level NTML+1 because ZH is above this level. - ntml(i,j) = ntml(i,j) + 1 - k=ntml(i,j)+1 - sml_disc_inv(i,j)=2 + zh_np1(i,j) = zero + t_frac(i,j) = zero + zrzi(i,j) = zero + we_rho(i,j) = zero + tothf_zh(i,j) = zero + totqf_zh(i,j) = zero + k=ntml(i,j)+1 - ! T_FRAC is fraction of timestep inversion is above - ! the entrainment flux grid-level (at Z_UV(K)) - t_frac(i,j) = (zh_np1(i,j)-z_uv(i,j,k)) / & - (zh_np1(i,j)-zh(i,j)) - ! ZH_FRAC is the timestep-average fraction of mixed layer - ! air in the inversion grid-level, level NTML+1 - zh_frac(i,j) = one_half*t_frac(i,j)*(zh_np1(i,j)-z_uv(i,j,k) ) & - / dzl(i,j,k) + ! Only RHOKH_ENT is passed out of EXCFNL so recalculate WE: + we_parm(i,j) = rdz(i,j,k)* & + ( rhokh_top_ent(i,j)+rhokh_surf_ent(i,j) ) & + / rho_mix(i,j,k) - else if ( zh_np1(i,j) >= z_uv(i,j,ntml(i,j)+1) ) then - ! ZH always between half-levels NTML+1 and NTML+2 + if ( sml_disc_inv(i,j) == 1 .and. .not. coupled(i,j) .and. & + (rhokh_top_ent(i,j)+rhokh_surf_ent(i,j)) > zero ) then - t_frac(i,j) = one - zh_frac(i,j) = ( one_half*(zh(i,j)+zh_np1(i,j)) - z_uv(i,j,k) ) & - / dzl(i,j,k) + !----------------------------------------------------------------- + !..Calculate ZH at end of timestep, ZH_NP1 + !----------------------------------------------------------------- + !..linearly interpolate vertical velocity to ZH + if ( zh(i,j) >= z_tq(i,j,k) ) then + w_ls(i,j) = w(i,j,k) + ( w(i,j,k+1) - w(i,j,k) ) & + * (zh(i,j)-z_tq(i,j,k)) * rdz(i,j,k+1) + else + w_ls(i,j) = w(i,j,k) + ( w(i,j,k) - w(i,j,k-1) ) & + * (zh(i,j)-z_tq(i,j,k)) * rdz(i,j,k) + end if + w_ls(i,j) = min ( w_ls(i,j), zero ) + ! only interested in subsidence + + zh_np1(i,j) = zh(i,j) + & + timestep * ( we_parm(i,j) + w_ls(i,j) ) + zh_np1(i,j) = max( zh_np1(i,j), z_uv(i,j,k-1) ) + if ( zh_np1(i,j) > z_top(i,j,k+1) ) then + ! limit ZH and W_e (and therefore the entraiment fluxes) + ! because the inversion cannot rise more than one level + ! in a timestep. + zh_np1(i,j) = z_top(i,j,k+1) + we_parm(i,j) = & + (z_top(i,j,k+1) - zh(i,j))/timestep - w_ls(i,j) + end if + !----------------------------------------------------------------- + !..Decide on which grid-level to apply entrainment flux + !----------------------------------------------------------------- + if ( zh_np1(i,j) > z_uv(i,j,ntml(i,j)+2) ) then + ! ZH risen above level K+1 so specify appropriate flux + ! at this level and raise NTML by one (this means + ! gradient-adjustment is also applied at half-level + ! old_NTML+1). Note KH profiles should already be + ! calculated at level NTML+1 because ZH is above this level. + ntml(i,j) = ntml(i,j) + 1 + k=ntml(i,j)+1 + sml_disc_inv(i,j)=2 + + ! T_FRAC is fraction of timestep inversion is above + ! the entrainment flux grid-level (at Z_UV(K)) + t_frac(i,j) = (zh_np1(i,j)-z_uv(i,j,k)) / & + (zh_np1(i,j)-zh(i,j)) + ! ZH_FRAC is the timestep-average fraction of mixed layer + ! air in the inversion grid-level, level NTML+1 + zh_frac(i,j) = one_half*t_frac(i,j)*(zh_np1(i,j)-z_uv(i,j,k) ) & + / dzl(i,j,k) + + else if ( zh_np1(i,j) >= z_uv(i,j,ntml(i,j)+1) ) then + ! ZH always between half-levels NTML+1 and NTML+2 + + t_frac(i,j) = one + zh_frac(i,j) = ( one_half*(zh(i,j)+zh_np1(i,j)) - z_uv(i,j,k) ) & + / dzl(i,j,k) - else - ! ZH falls below half-level NTML+1 - ! Keep implicit (diffusive) entrainment but apply - ! at the level below - if (ntml(i,j) >= 2) then ! ftl(k=1) is surface flux - ntml(i,j) = ntml(i,j) - 1 - k=ntml(i,j)+1 - rhokh_top(i,j,k+1) = zero ! also need to remove diffusion - rhokh(i,j,k+1) = zero ! at old entrainment grid-level - end if - t_frac(i,j) = zero - zh_frac(i,j) = zero - sml_disc_inv(i,j) = 0 + else + ! ZH falls below half-level NTML+1 + ! Keep implicit (diffusive) entrainment but apply + ! at the level below + if (ntml(i,j) >= 2) then ! ftl(k=1) is surface flux + ntml(i,j) = ntml(i,j) - 1 + k=ntml(i,j)+1 + rhokh_top(i,j,k+1) = zero ! also need to remove diffusion + rhokh(i,j,k+1) = zero ! at old entrainment grid-level + end if + t_frac(i,j) = zero + zh_frac(i,j) = zero + sml_disc_inv(i,j) = 0 - end if ! test on where to apply entrainment flux + end if ! test on where to apply entrainment flux - we_rho(i,j) = rho_mix(i,j,k) * we_parm(i,j) - zrzi(i,j) = z_uv(i,j,k)*2.0_r_bl/(zh(i,j)+zh_np1(i,j)) + we_rho(i,j) = rho_mix(i,j,k) * we_parm(i,j) + zrzi(i,j) = z_uv(i,j,k)*2.0_r_bl/(zh(i,j)+zh_np1(i,j)) - end if ! test on SML_DISC_INV, etc - end do + end if ! test on SML_DISC_INV, etc end do !$OMP end do NOWAIT +!end do + !----------------------------------------------------------------------- !..Linearly interpolate between the known total (turb+rad+subs+micro) @@ -4104,313 +4193,319 @@ subroutine kmkhz_9c ( & end if c_tke = 1.33_r_bl/(vkman*c_ws**two_thirds) -!$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ! Entrainment flux applied to level NTML+1 which is the - ! flux-level above the top of the SML - k=ntml(i,j)+1 - - if ( t_frac(i,j) > zero ) then - - rhokh_top(i,j,k) = zero ! apply entrainment explicitly - rhokh(i,j,k) = zero ! " - - tothf_zh(i,j) = - we_rho(i,j)*dsl_sml(i,j) + ft_nt_zh(i,j) - ! Linearly interpolate to entrainment flux grid-level - tothf_efl = ft_nt(i,j,1) + ftl(i,j,1) + & - ( tothf_zh(i,j)-ft_nt(i,j,1)-ftl(i,j,1) )*zrzi(i,j) - ! Ensure total heat flux gradient in inversion grid-level is - ! consistent with inversion rising (ie. implies cooling in - ! level K relative to the mixed layer) or falling - ! (implies warming) - - ml_tend = -( tothf_zh(i,j)-ft_nt(i,j,1)-ftl(i,j,1) ) / zh(i,j) - fa_tend = zero - if ( k+1 <= bl_levels ) & - fa_tend = - ( ft_nt(i,j,k+2) - ft_nt(i,j,k+1) ) & - / dzl(i,j,k+1) - inv_tend = zh_frac(i,j) * ml_tend & - + (one-zh_frac(i,j)) * fa_tend - if (we_parm(i,j)+w_ls(i,j) >= zero) then - ! Inversion moving up so inversion level should cool - ! Ensure it does cool relative to ML - tothf_efl = min( tothf_efl, & - ft_nt(i,j,k+1)+inv_tend*dzl(i,j,k) ) - ! Ensure inversion level won't end up colder than - ! NTML by end of timestep. - ! Set INV_TEND to max allowable cooling rate, also - ! allowing for change in ML_TEND arising from this change - ! to TOTHF_EFL: - inv_tend = (sl(i,j,k-1)-sl(i,j,k))/timestep & - + (ft_nt(i,j,1)+ftl(i,j,1))/z_uv(i,j,k) - tothf_efl = max( tothf_efl, & - (ft_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & - /(one+ dzl(i,j,k)/z_uv(i,j,k)) ) - else ! WE_PARM+W_LS < 0 - ! Ensure inversion level does warm relative to ML - tothf_efl = max( tothf_efl, & - ft_nt(i,j,k+1)+inv_tend*dzl(i,j,k) ) +! do j = pdims%j_start, pdims%j_end +!$OMP do SCHEDULE(STATIC) +do i = pdims%i_start, pdims%i_end + + ! Entrainment flux applied to level NTML+1 which is the + ! flux-level above the top of the SML + k=ntml(i,j)+1 + + if ( t_frac(i,j) > zero ) then + + rhokh_top(i,j,k) = zero ! apply entrainment explicitly + rhokh(i,j,k) = zero ! " + + tothf_zh(i,j) = - we_rho(i,j)*dsl_sml(i,j) + ft_nt_zh(i,j) + ! Linearly interpolate to entrainment flux grid-level + tothf_efl = ft_nt(i,j,1) + ftl(i,j,1) + & + ( tothf_zh(i,j)-ft_nt(i,j,1)-ftl(i,j,1) )*zrzi(i,j) + ! Ensure total heat flux gradient in inversion grid-level is + ! consistent with inversion rising (ie. implies cooling in + ! level K relative to the mixed layer) or falling + ! (implies warming) + + ml_tend = -( tothf_zh(i,j)-ft_nt(i,j,1)-ftl(i,j,1) ) / zh(i,j) + fa_tend = zero + if ( k+1 <= bl_levels ) & + fa_tend = - ( ft_nt(i,j,k+2) - ft_nt(i,j,k+1) ) & + / dzl(i,j,k+1) + inv_tend = zh_frac(i,j) * ml_tend & + + (one-zh_frac(i,j)) * fa_tend + if (we_parm(i,j)+w_ls(i,j) >= zero) then + ! Inversion moving up so inversion level should cool + ! Ensure it does cool relative to ML + tothf_efl = min( tothf_efl, & + ft_nt(i,j,k+1)+inv_tend*dzl(i,j,k) ) + ! Ensure inversion level won't end up colder than + ! NTML by end of timestep. + ! Set INV_TEND to max allowable cooling rate, also + ! allowing for change in ML_TEND arising from this change + ! to TOTHF_EFL: + inv_tend = (sl(i,j,k-1)-sl(i,j,k))/timestep & + + (ft_nt(i,j,1)+ftl(i,j,1))/z_uv(i,j,k) + tothf_efl = max( tothf_efl, & + (ft_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & + /(one+ dzl(i,j,k)/z_uv(i,j,k)) ) + else ! WE_PARM+W_LS < 0 + ! Ensure inversion level does warm relative to ML + tothf_efl = max( tothf_efl, & + ft_nt(i,j,k+1)+inv_tend*dzl(i,j,k) ) + end if + ! Turbulent entrainment flux is then the residual of the total + ! flux and the net flux from other processes + ftl(i,j,k) = t_frac(i,j) * ( tothf_efl - ft_nt(i,j,k) ) + else ! not specifying entrainment flux but KH + ! Include entrainment KH in K-profiles, if greater + ! (for COUPLED layers these will be zero) + rhokh_top(i,j,k) = max( rhokh_top(i,j,k), rhokh_top_ent(i,j) ) + rhokh(i,j,k) = max( rhokh(i,j,k), rhokh_surf_ent(i,j) ) + + if (res_inv(i,j) == 1) then + Prandtl = min( rhokm(i,j,k)/(rbl_eps+rhokh_surf_ent(i,j)), & + pr_max ) + if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then + ! need velocity scale for TKE diagnostic + w_m = ( v_s(i,j)*v_s(i,j)*v_s(i,j) + & + c_ws * zh(i,j) * fb_surf(i,j) ) ** one_third end if - ! Turbulent entrainment flux is then the residual of the total - ! flux and the net flux from other processes - ftl(i,j,k) = t_frac(i,j) * ( tothf_efl - ft_nt(i,j,k) ) - else ! not specifying entrainment flux but KH - ! Include entrainment KH in K-profiles, if greater - ! (for COUPLED layers these will be zero) - rhokh_top(i,j,k) = max( rhokh_top(i,j,k), rhokh_top_ent(i,j) ) - rhokh(i,j,k) = max( rhokh(i,j,k), rhokh_surf_ent(i,j) ) - - if (res_inv(i,j) == 1) then - Prandtl = min( rhokm(i,j,k)/(rbl_eps+rhokh_surf_ent(i,j)), & - pr_max ) - if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then - ! need velocity scale for TKE diagnostic - w_m = ( v_s(i,j)*v_s(i,j)*v_s(i,j) + & - c_ws * zh(i,j) * fb_surf(i,j) ) ** one_third - end if - if (bl_res_inv == cosine_inv_flux) then - svl_lapse_rho = (svl(i,j,k)-svl(i,j,k-1)) / & - ( (z_tq(i,j,k)-z_tq(i,j,k-1))*rho_mix(i,j,k) ) - kl=k+1 - do while ( z_uv(i,j,kl) < zh(i,j)+dzh(i,j) .and. & - kl <= bl_levels ) - recip_svl_lapse = (z_tq(i,j,kl)-z_tq(i,j,kl-1))/ & - max( 0.01_r_bl, svl(i,j,kl)-svl(i,j,kl-1) ) - rhok_inv = rhokh_surf_ent(i,j) * svl_lapse_rho * & - rho_mix(i,j,kl) * recip_svl_lapse * & - cos(one_half*pi*(z_uv(i,j,kl)-zh(i,j))/dzh(i,j)) - rhok_inv = min( rhok_inv, 1000.0_r_bl ) - rhokh(i,j,kl) = max( rhokh(i,j,kl), rhok_inv ) - ! rescale for KM on staggered grid - rhok_inv = Prandtl * rhok_inv & - * rdz(i,j,kl) * (z_uv(i,j,kl)-z_uv(i,j,kl-1)) & - * rho_wet_tq(i,j,kl-1) / rho_mix(i,j,kl) - rhokm(i,j,kl) = max( rhokm(i,j,kl), rhok_inv ) - if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then - ! save Km/timescale for TKE diag, completed in bdy_expl2 - tke_nl(i,j,kl) = max( tke_nl(i,j,kl), rhok_inv*c_tke*w_m/zh(i,j)) - end if - kl=kl+1 - end do - else if (bl_res_inv == target_inv_profile) then - svl_lapse = (svl(i,j,k)-svl(i,j,k-1)) / & - ( (z_tq(i,j,k)-z_tq(i,j,k-1)) ) - kp=k+1 ! kp marks the lowest level above the inversion - do while ( z_uv(i,j,kp) < zh(i,j)+dzh(i,j) .and. & - kp <= bl_levels ) - kp=kp+1 - end do - svl_flux(k) = - rhokh_surf_ent(i,j) * svl_lapse - kl=k+1 - do while ( z_uv(i,j,kl) < zh(i,j)+dzh(i,j) .and. & - kl <= bl_levels ) - ! assume a linear target svl profile within inversion - svl_target = svl(i,j,k-1) + (svl(i,j,kp)-svl(i,j,k-1)) * & - (z_uv(i,j,kl)-zh(i,j)) / dzh(i,j) - rho_dz = rho_mix_tq(i,j,kl) * dzl(i,j,kl) - svl_flux(kl) = svl_flux(kl-1) - & - (svl_target-svl(i,j,kl))*rho_dz/timestep - kl=kl+1 - end do - ! linearly extrapolate flux to inversion top - svl_flux(kp)=svl_flux(kp-1) + (svl_flux(kp-1)-svl_flux(kp-2))* & - (zh(i,j)+dzh(i,j)-z_uv(i,j,kp-1))*rdz(i,j,kp-1) - kl=k+1 - do while ( z_uv(i,j,kl) < zh(i,j)+dzh(i,j) .and. & - kl <= bl_levels ) - ! rescale svl_flux so as to have zero flux at the inversion top - ! ie so svl_flux(kp)=0 - svl_flux(kl) = svl_flux(k)*( one - & - (svl_flux(kl)-svl_flux(k))/ & - (svl_flux(kp)-svl_flux(k)) ) - recip_svl_lapse = (z_tq(i,j,kl)-z_tq(i,j,kl-1))/ & - max( 0.01_r_bl, svl(i,j,kl)-svl(i,j,kl-1) ) - rhok_inv = - svl_flux(kl) * recip_svl_lapse - - rhok_inv = min( rhok_inv, 1000.0_r_bl ) - rhokh(i,j,kl) = max( rhokh(i,j,kl), rhok_inv ) - ! rescale for KM on staggered grid - rhok_inv = Prandtl * rhok_inv & - * rdz(i,j,kl) * (z_uv(i,j,kl)-z_uv(i,j,kl-1)) & - * rho_wet_tq(i,j,kl-1) / rho_mix(i,j,kl) - rhokm(i,j,kl) = max( rhokm(i,j,kl), rhok_inv ) - if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then - ! save Km/timescale for TKE diag, completed in bdy_expl2 - tke_nl(i,j,kl) = max( tke_nl(i,j,kl), rhok_inv*c_tke*w_m/zh(i,j)) - end if - kl=kl+1 - end do - end if ! bl_res_inv option - end if ! res_inv - end if ! test on T_FRAC gt 0 + if (bl_res_inv == cosine_inv_flux) then + svl_lapse_rho = (svl(i,j,k)-svl(i,j,k-1)) / & + ( (z_tq(i,j,k)-z_tq(i,j,k-1))*rho_mix(i,j,k) ) + kl=k+1 + do while ( z_uv(i,j,kl) < zh(i,j)+dzh(i,j) .and. & + kl <= bl_levels ) + recip_svl_lapse = (z_tq(i,j,kl)-z_tq(i,j,kl-1))/ & + max( 0.01_r_bl, svl(i,j,kl)-svl(i,j,kl-1) ) + rhok_inv = rhokh_surf_ent(i,j) * svl_lapse_rho * & + rho_mix(i,j,kl) * recip_svl_lapse * & + cos(one_half*pi*(z_uv(i,j,kl)-zh(i,j))/dzh(i,j)) + rhok_inv = min( rhok_inv, 1000.0_r_bl ) + rhokh(i,j,kl) = max( rhokh(i,j,kl), rhok_inv ) + ! rescale for KM on staggered grid + rhok_inv = Prandtl * rhok_inv & + * rdz(i,j,kl) * (z_uv(i,j,kl)-z_uv(i,j,kl-1)) & + * rho_wet_tq(i,j,kl-1) / rho_mix(i,j,kl) + rhokm(i,j,kl) = max( rhokm(i,j,kl), rhok_inv ) + if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then + ! save Km/timescale for TKE diag, completed in bdy_expl2 + tke_nl(i,j,kl) = max( tke_nl(i,j,kl), rhok_inv*c_tke*w_m/zh(i,j)) + end if + kl=kl+1 + end do + else if (bl_res_inv == target_inv_profile) then + svl_lapse = (svl(i,j,k)-svl(i,j,k-1)) / & + ( (z_tq(i,j,k)-z_tq(i,j,k-1)) ) + kp=k+1 ! kp marks the lowest level above the inversion + do while ( z_uv(i,j,kp) < zh(i,j)+dzh(i,j) .and. & + kp <= bl_levels ) + kp=kp+1 + end do + svl_flux(k) = - rhokh_surf_ent(i,j) * svl_lapse + kl=k+1 + do while ( z_uv(i,j,kl) < zh(i,j)+dzh(i,j) .and. & + kl <= bl_levels ) + ! assume a linear target svl profile within inversion + svl_target = svl(i,j,k-1) + (svl(i,j,kp)-svl(i,j,k-1)) * & + (z_uv(i,j,kl)-zh(i,j)) / dzh(i,j) + rho_dz = rho_mix_tq(i,j,kl) * dzl(i,j,kl) + svl_flux(kl) = svl_flux(kl-1) - & + (svl_target-svl(i,j,kl))*rho_dz/timestep + kl=kl+1 + end do + ! linearly extrapolate flux to inversion top + svl_flux(kp)=svl_flux(kp-1) + (svl_flux(kp-1)-svl_flux(kp-2))* & + (zh(i,j)+dzh(i,j)-z_uv(i,j,kp-1))*rdz(i,j,kp-1) + kl=k+1 + do while ( z_uv(i,j,kl) < zh(i,j)+dzh(i,j) .and. & + kl <= bl_levels ) + ! rescale svl_flux so as to have zero flux at the inversion top + ! ie so svl_flux(kp)=0 + svl_flux(kl) = svl_flux(k)*( one - & + (svl_flux(kl)-svl_flux(k))/ & + (svl_flux(kp)-svl_flux(k)) ) + recip_svl_lapse = (z_tq(i,j,kl)-z_tq(i,j,kl-1))/ & + max( 0.01_r_bl, svl(i,j,kl)-svl(i,j,kl-1) ) + rhok_inv = - svl_flux(kl) * recip_svl_lapse + + rhok_inv = min( rhok_inv, 1000.0_r_bl ) + rhokh(i,j,kl) = max( rhokh(i,j,kl), rhok_inv ) + ! rescale for KM on staggered grid + rhok_inv = Prandtl * rhok_inv & + * rdz(i,j,kl) * (z_uv(i,j,kl)-z_uv(i,j,kl-1)) & + * rho_wet_tq(i,j,kl-1) / rho_mix(i,j,kl) + rhokm(i,j,kl) = max( rhokm(i,j,kl), rhok_inv ) + if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then + ! save Km/timescale for TKE diag, completed in bdy_expl2 + tke_nl(i,j,kl) = max( tke_nl(i,j,kl), rhok_inv*c_tke*w_m/zh(i,j)) + end if + kl=kl+1 + end do + end if ! bl_res_inv option + end if ! res_inv + end if ! test on T_FRAC gt 0 - end do end do !$OMP end do NOWAIT +!end do + !------------------------------------------------- !..Second the decoupled mixed layer, if entraining !------------------------------------------------- + +!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 - zhsc_np1(i,j) = zero - t_frac_dsc(i,j) = zero - zrzi_dsc(i,j) = zero - we_rho_dsc(i,j) = zero - tothf_zhsc(i,j) = zero - totqf_zhsc(i,j) = zero + zhsc_np1(i,j) = zero + t_frac_dsc(i,j) = zero + zrzi_dsc(i,j) = zero + we_rho_dsc(i,j) = zero + tothf_zhsc(i,j) = zero + totqf_zhsc(i,j) = zero - k=ntdsc(i,j)+1 - we_dsc_parm(i,j) = rdz(i,j,k)*rhokh_dsct_ent(i,j) & - / rho_mix(i,j,k) + k=ntdsc(i,j)+1 + we_dsc_parm(i,j) = rdz(i,j,k)*rhokh_dsct_ent(i,j) & + / rho_mix(i,j,k) - if ( dsc_disc_inv(i,j) == 1 & - .and. rhokh_dsct_ent(i,j) > zero ) then + if ( dsc_disc_inv(i,j) == 1 & + .and. rhokh_dsct_ent(i,j) > zero ) then - !----------------------------------------------------------------- - !..Calculate ZHSC at end of timestep, ZHSC_NP1 - !----------------------------------------------------------------- - !..interpolate vertical velocity to ZH - if ( zhsc(i,j) >= z_tq(i,j,k) ) then - w_ls_dsc(i,j) = w(i,j,k) + ( w(i,j,k+1) - w(i,j,k) ) * & - (zhsc(i,j)-z_tq(i,j,k)) * rdz(i,j,k+1) - else - w_ls_dsc(i,j) = w(i,j,k) + ( w(i,j,k) - w(i,j,k-1) ) * & - (zhsc(i,j)-z_tq(i,j,k)) * rdz(i,j,k) - end if - w_ls_dsc(i,j) = min ( w_ls_dsc(i,j), zero ) - ! only interested in subsidence - - zhsc_np1(i,j) = zhsc(i,j) + & - timestep * ( we_dsc_parm(i,j) + w_ls_dsc(i,j) ) - zhsc_np1(i,j) = max( zhsc_np1(i,j), z_uv(i,j,k-1) ) - if ( zhsc_np1(i,j) > z_top(i,j,k+1) ) then - ! limit ZHSC and W_e (and therefore the entrainment fluxes) - ! because the inversion cannot rise more than one level - ! in a timestep. - zhsc_np1(i,j) = z_top(i,j,k+1) - we_dsc_parm(i,j) = & - (z_top(i,j,k+1) - zhsc(i,j))/timestep - w_ls_dsc(i,j) - end if - !----------------------------------------------------------------- - !..Decide on which grid-level to apply entrainment flux - !----------------------------------------------------------------- - if ( zhsc_np1(i,j) > z_uv(i,j,ntdsc(i,j)+2) ) then - ! ZHSC risen above level K+1 so specify appropriate - ! flux at this level and raise NTDSC by one + !----------------------------------------------------------------- + !..Calculate ZHSC at end of timestep, ZHSC_NP1 + !----------------------------------------------------------------- + !..interpolate vertical velocity to ZH + if ( zhsc(i,j) >= z_tq(i,j,k) ) then + w_ls_dsc(i,j) = w(i,j,k) + ( w(i,j,k+1) - w(i,j,k) ) * & + (zhsc(i,j)-z_tq(i,j,k)) * rdz(i,j,k+1) + else + w_ls_dsc(i,j) = w(i,j,k) + ( w(i,j,k) - w(i,j,k-1) ) * & + (zhsc(i,j)-z_tq(i,j,k)) * rdz(i,j,k) + end if + w_ls_dsc(i,j) = min ( w_ls_dsc(i,j), zero ) + ! only interested in subsidence + + zhsc_np1(i,j) = zhsc(i,j) + & + timestep * ( we_dsc_parm(i,j) + w_ls_dsc(i,j) ) + zhsc_np1(i,j) = max( zhsc_np1(i,j), z_uv(i,j,k-1) ) + if ( zhsc_np1(i,j) > z_top(i,j,k+1) ) then + ! limit ZHSC and W_e (and therefore the entrainment fluxes) + ! because the inversion cannot rise more than one level + ! in a timestep. + zhsc_np1(i,j) = z_top(i,j,k+1) + we_dsc_parm(i,j) = & + (z_top(i,j,k+1) - zhsc(i,j))/timestep - w_ls_dsc(i,j) + end if + !----------------------------------------------------------------- + !..Decide on which grid-level to apply entrainment flux + !----------------------------------------------------------------- + if ( zhsc_np1(i,j) > z_uv(i,j,ntdsc(i,j)+2) ) then + ! ZHSC risen above level K+1 so specify appropriate + ! flux at this level and raise NTDSC by one - ntdsc(i,j) = ntdsc(i,j) + 1 - k = ntdsc(i,j)+1 - dsc_disc_inv(i,j) = 2 - t_frac_dsc(i,j) = (zhsc_np1(i,j)-z_uv(i,j,k)) / & - (zhsc_np1(i,j)-zhsc(i,j)) + ntdsc(i,j) = ntdsc(i,j) + 1 + k = ntdsc(i,j)+1 + dsc_disc_inv(i,j) = 2 + t_frac_dsc(i,j) = (zhsc_np1(i,j)-z_uv(i,j,k)) / & + (zhsc_np1(i,j)-zhsc(i,j)) - zhsc_frac(i,j) = one_half*t_frac_dsc(i,j)* & - ( zhsc_np1(i,j)-z_uv(i,j,k) )/ dzl(i,j,k) + zhsc_frac(i,j) = one_half*t_frac_dsc(i,j)* & + ( zhsc_np1(i,j)-z_uv(i,j,k) )/ dzl(i,j,k) - else if ( zhsc_np1(i,j) > z_uv(i,j,ntdsc(i,j)+1) ) then - ! ZHSC always between half-levels NTDSC+1 and NTDSC+2 + else if ( zhsc_np1(i,j) > z_uv(i,j,ntdsc(i,j)+1) ) then + ! ZHSC always between half-levels NTDSC+1 and NTDSC+2 - t_frac_dsc(i,j) = one - zhsc_frac(i,j) = ( one_half*(zhsc(i,j)+zhsc_np1(i,j)) & - - z_uv(i,j,k) )/ dzl(i,j,k) + t_frac_dsc(i,j) = one + zhsc_frac(i,j) = ( one_half*(zhsc(i,j)+zhsc_np1(i,j)) & + - z_uv(i,j,k) )/ dzl(i,j,k) - else - ! ZHSC falls below half-level NTDSC+1 - ! Keep implicit (diffusive) entrainment but apply - ! at the level below - ntdsc(i,j) = ntdsc(i,j) - 1 ! could reduce NTDSC to 1 - k = ntdsc(i,j)+1 - rhokh_top(i,j,k+1) = zero - rhokh(i,j,k+1) = zero - - t_frac_dsc(i,j) = zero - zhsc_frac(i,j) = zero - dsc_disc_inv(i,j) = 0 + else + ! ZHSC falls below half-level NTDSC+1 + ! Keep implicit (diffusive) entrainment but apply + ! at the level below + ntdsc(i,j) = ntdsc(i,j) - 1 ! could reduce NTDSC to 1 + k = ntdsc(i,j)+1 + rhokh_top(i,j,k+1) = zero + rhokh(i,j,k+1) = zero - end if ! test on where to apply entrainment flux + t_frac_dsc(i,j) = zero + zhsc_frac(i,j) = zero + dsc_disc_inv(i,j) = 0 - we_rho_dsc(i,j) = rho_mix(i,j,k) * we_dsc_parm(i,j) - ! for z'/z_i' assume height of DSC base is fixed in time - zrzi_dsc(i,j) =( z_uv(i,j,k)-(zhsc(i,j)-dscdepth(i,j)) ) & - /( dscdepth(i,j)+one_half*(zhsc_np1(i,j)-zhsc(i,j)) ) + end if ! test on where to apply entrainment flux - end if ! test on DSC_DISC_INV, etc - end do + we_rho_dsc(i,j) = rho_mix(i,j,k) * we_dsc_parm(i,j) + ! for z'/z_i' assume height of DSC base is fixed in time + zrzi_dsc(i,j) =( z_uv(i,j,k)-(zhsc(i,j)-dscdepth(i,j)) ) & + /( dscdepth(i,j)+one_half*(zhsc_np1(i,j)-zhsc(i,j)) ) + + end if ! test on DSC_DISC_INV, etc end do !$OMP end do NOWAIT +! end do + !----------------------------------------------------------------------- !..Linearly interpolate between the known total (turb+rad+subs+micro) !..flux at the DSC base and the parametrized flux at the inversion !----------------------------------------------------------------------- -!$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - - ! Entrainment flux applied to level NTDSC+1 which is the - ! flux-level above the top of the DSC layer - k=ntdsc(i,j)+1 - - if ( t_frac_dsc(i,j) > zero ) then - - rhokh_top(i,j,k) = zero ! apply entrainment explicitly - rhokh(i,j,k) = zero ! " - - tothf_zhsc(i,j) = - we_rho_dsc(i,j)*dsl_dsc(i,j) & - + ft_nt_zhsc(i,j) - tothf_efl = ft_nt_dscb(i,j) + & - ( tothf_zhsc(i,j)-ft_nt_dscb(i,j) )*zrzi_dsc(i,j) - ! Ensure total heat flux gradient in inversion grid-level is - ! consistent with inversion rising (implies cooling in - ! level K, relative to mixed layer) or falling - ! (implies warming) - ml_tend = - ( tothf_zhsc(i,j)-ft_nt_dscb(i,j) )/ dscdepth(i,j) - fa_tend = zero - if ( k+1 <= bl_levels ) & - fa_tend = - ( ft_nt(i,j,k+2) - ft_nt(i,j,k+1) ) & - / dzl(i,j,k+1) - inv_tend = zhsc_frac(i,j) * ml_tend & - + (one-zhsc_frac(i,j)) * fa_tend - if (we_dsc_parm(i,j)+w_ls_dsc(i,j) >= zero) then - ! Inversion moving up so inversion level should cool - ! Ensure it does cool relative to ML - tothf_efl = min( tothf_efl, & - ft_nt(i,j,k+1)+inv_tend*dzl(i,j,k) ) - ! Ensure inversion level won't end up colder than - ! NTDSC by end of timestep. - inv_tend = (sl(i,j,k-1)-sl(i,j,k))/timestep & - + ft_nt_dscb(i,j)/dscdepth(i,j) - tothf_efl = max( tothf_efl, & - (ft_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & - /(one+ dzl(i,j,k)/dscdepth(i,j)) ) - else ! WE_DSC_PARM+W_LS_DSC < 0 - tothf_efl = max( tothf_efl, & - ft_nt(i,j,k+1)+inv_tend*dzl(i,j,k) ) - end if - ! Turbulent entrainment flux is then the residual of the total - ! flux and the net flux from other processes - ftl(i,j,k) = t_frac_dsc(i,j) * ( tothf_efl - ft_nt(i,j,k) ) +! do j = pdims%j_start, pdims%j_end +!$OMP do SCHEDULE(STATIC) +do i = pdims%i_start, pdims%i_end + + ! Entrainment flux applied to level NTDSC+1 which is the + ! flux-level above the top of the DSC layer + k=ntdsc(i,j)+1 + + if ( t_frac_dsc(i,j) > zero ) then + + rhokh_top(i,j,k) = zero ! apply entrainment explicitly + rhokh(i,j,k) = zero ! " + + tothf_zhsc(i,j) = - we_rho_dsc(i,j)*dsl_dsc(i,j) & + + ft_nt_zhsc(i,j) + tothf_efl = ft_nt_dscb(i,j) + & + ( tothf_zhsc(i,j)-ft_nt_dscb(i,j) )*zrzi_dsc(i,j) + ! Ensure total heat flux gradient in inversion grid-level is + ! consistent with inversion rising (implies cooling in + ! level K, relative to mixed layer) or falling + ! (implies warming) + ml_tend = - ( tothf_zhsc(i,j)-ft_nt_dscb(i,j) )/ dscdepth(i,j) + fa_tend = zero + if ( k+1 <= bl_levels ) & + fa_tend = - ( ft_nt(i,j,k+2) - ft_nt(i,j,k+1) ) & + / dzl(i,j,k+1) + inv_tend = zhsc_frac(i,j) * ml_tend & + + (one-zhsc_frac(i,j)) * fa_tend + + if (we_dsc_parm(i,j)+w_ls_dsc(i,j) >= zero) then + ! Inversion moving up so inversion level should cool + ! Ensure it does cool relative to ML + tothf_efl = min( tothf_efl, & + ft_nt(i,j,k+1)+inv_tend*dzl(i,j,k) ) + ! Ensure inversion level won't end up colder than + ! NTDSC by end of timestep. + inv_tend = (sl(i,j,k-1)-sl(i,j,k))/timestep & + + ft_nt_dscb(i,j)/dscdepth(i,j) + tothf_efl = max( tothf_efl, & + (ft_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & + /(one+ dzl(i,j,k)/dscdepth(i,j)) ) + else ! WE_DSC_PARM+W_LS_DSC < 0 + tothf_efl = max( tothf_efl, & + ft_nt(i,j,k+1)+inv_tend*dzl(i,j,k) ) + end if + ! Turbulent entrainment flux is then the residual of the total + ! flux and the net flux from other processes + ftl(i,j,k) = t_frac_dsc(i,j) * ( tothf_efl - ft_nt(i,j,k) ) - else if ( dsc(i,j) ) then + else if ( dsc(i,j) ) then - ! Not specifying entrainment flux but KH - ! Include entrainment KH in K-profile, if greater - rhokh_top(i,j,k) = max( rhokh_top(i,j,k),rhokh_dsct_ent(i,j) ) + ! Not specifying entrainment flux but KH + ! Include entrainment KH in K-profile, if greater + rhokh_top(i,j,k) = max( rhokh_top(i,j,k),rhokh_dsct_ent(i,j) ) - end if ! if not DSC + end if ! if not DSC - end do end do !$OMP end do NOWAIT +!end do + !----------------------------------------------------------------------- ! Specify QW entrainment fluxes @@ -4424,42 +4519,45 @@ subroutine kmkhz_9c ( & ! the inversion grid-level is physically within the BL) ! ------------------------------------------------------------------ + +! 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 - fq_nt_dscb(i,j) = fq_nt(i,j,1) - if ( nbdsc(i,j) > 1 ) then - k = nbdsc(i,j) ! NBDSC marks the lowest flux-level - ! within the DSC layer - ! Interpolate non-turb flux to base - ! of DSC layer: - fq_nt_dscb(i,j) = fq_nt(i,j,k-1) + & - (fq_nt(i,j,k)-fq_nt(i,j,k-1)) & - *(zdsc_base(i,j)-z_uv(i,j,k-1))/dzl(i,j,k-1) - end if - end do +do i = pdims%i_start, pdims%i_end + fq_nt_dscb(i,j) = fq_nt(i,j,1) + if ( nbdsc(i,j) > 1 ) then + k = nbdsc(i,j) ! NBDSC marks the lowest flux-level + ! within the DSC layer + ! Interpolate non-turb flux to base + ! of DSC layer: + fq_nt_dscb(i,j) = fq_nt(i,j,k-1) + & + (fq_nt(i,j,k)-fq_nt(i,j,k-1)) & + *(zdsc_base(i,j)-z_uv(i,j,k-1))/dzl(i,j,k-1) + end if end do !$OMP end do NOWAIT +! end do + ! Repeat for water tracers if (l_wtrac) then do i_wt = 1, n_wtrac -!$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_dscb(i,j) = wtrac_bl(i_wt)%fq_nt(i,j,1) - if ( nbdsc(i,j) > 1 ) then - k = nbdsc(i,j) ! NBDSC marks the lowest flux-level - ! within the DSC layer - ! Interpolate non-turb flux to base - ! of DSC layer: - wtrac_bl(i_wt)%fq_nt_dscb(i,j) = wtrac_bl(i_wt)%fq_nt(i,j,k-1) + & - ( wtrac_bl(i_wt)%fq_nt(i,j,k)- wtrac_bl(i_wt)%fq_nt(i,j,k-1)) & - *(zdsc_base(i,j)-z_uv(i,j,k-1))/dzl(i,j,k-1) - end if - 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_dscb(i,j) = wtrac_bl(i_wt)%fq_nt(i,j,1) + if ( nbdsc(i,j) > 1 ) then + k = nbdsc(i,j) ! NBDSC marks the lowest flux-level + ! within the DSC layer + ! Interpolate non-turb flux to base + ! of DSC layer: + wtrac_bl(i_wt)%fq_nt_dscb(i,j) = wtrac_bl(i_wt)%fq_nt(i,j,k-1) + & + ( wtrac_bl(i_wt)%fq_nt(i,j,k)- wtrac_bl(i_wt)%fq_nt(i,j,k-1)) & + *(zdsc_base(i,j)-z_uv(i,j,k-1))/dzl(i,j,k-1) + end if end do -!$OMP end do + !$OMP end do + ! end do + end do end if ! l_wtrac @@ -4468,101 +4566,103 @@ subroutine kmkhz_9c ( & ! microphysical and subsidence fluxes are correctly coupled. !----------------------------------------------------------------------- -!$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - moisten(i,j) = .false. - if (l_wtrac) totqf_efl_meth1(i,j) = 0 - if (l_wtrac) totqf_efl_meth2(i,j) = 0 - k = ntml(i,j)+1 - - if ( t_frac(i,j) > zero ) then - ! Calculate total (turb+micro+subs) QW flux at subgrid - ! inversion height - totqf_zh(i,j) = - we_rho(i,j)*dqw_sml(i,j) + fq_nt_zh(i,j) - ! Interpolate to entrainment flux-level below - totqf_efl = fq_nt(i,j,1) + fqw(i,j,1) + zrzi(i,j) * & - ( totqf_zh(i,j) - fq_nt(i,j,1) - fqw(i,j,1) ) - ! Need to ensure the total QW flux gradient in inversion - ! grid-level is consistent with inversion rising or falling. - ! If QW(K) is drier than mixed layer then inversion rising - ! implies moistening in level K relative to mixed layer - ! while falling would imply relative drying of level K. - ! If QW(K) is moister than ML then want opposite tendencies. - ml_tend = - ( totqf_zh(i,j)-fq_nt(i,j,1)-fqw(i,j,1) ) /zh(i,j) - fa_tend = zero - if ( k+1 <= bl_levels ) & - fa_tend = - ( fq_nt(i,j,k+2)-fq_nt(i,j,k+1) ) & - / dzl(i,j,k+1) - inv_tend = zh_frac(i,j) * ml_tend & - + (one-zh_frac(i,j)) * fa_tend - if (we_parm(i,j)+w_ls(i,j) >= zero) then - ! inversion moving up so inversion will moisten/dry - ! depending on relative QW in level below - moisten(i,j) = ( qw(i,j,k) <= qw(i,j,k-1) ) - else - ! inversion moving down so inversion will moisten/dry - ! depending on relative QW in level above - moisten(i,j) = ( qw(i,j,k) <= qw(i,j,k+1) ) - end if +! do j = pdims%j_start, pdims%j_end +!$OMP do SCHEDULE(STATIC) +do i = pdims%i_start, pdims%i_end + moisten(i,j) = .false. + if (l_wtrac) totqf_efl_meth1(i,j) = 0 + if (l_wtrac) totqf_efl_meth2(i,j) = 0 + k = ntml(i,j)+1 + + if ( t_frac(i,j) > zero ) then + ! Calculate total (turb+micro+subs) QW flux at subgrid + ! inversion height + totqf_zh(i,j) = - we_rho(i,j)*dqw_sml(i,j) + fq_nt_zh(i,j) + ! Interpolate to entrainment flux-level below + totqf_efl = fq_nt(i,j,1) + fqw(i,j,1) + zrzi(i,j) * & + ( totqf_zh(i,j) - fq_nt(i,j,1) - fqw(i,j,1) ) + ! Need to ensure the total QW flux gradient in inversion + ! grid-level is consistent with inversion rising or falling. + ! If QW(K) is drier than mixed layer then inversion rising + ! implies moistening in level K relative to mixed layer + ! while falling would imply relative drying of level K. + ! If QW(K) is moister than ML then want opposite tendencies. + ml_tend = - ( totqf_zh(i,j)-fq_nt(i,j,1)-fqw(i,j,1) ) /zh(i,j) + fa_tend = zero + if ( k+1 <= bl_levels ) & + fa_tend = - ( fq_nt(i,j,k+2)-fq_nt(i,j,k+1) ) & + / dzl(i,j,k+1) + inv_tend = zh_frac(i,j) * ml_tend & + + (one-zh_frac(i,j)) * fa_tend + + if (we_parm(i,j)+w_ls(i,j) >= zero) then + ! inversion moving up so inversion will moisten/dry + ! depending on relative QW in level below + moisten(i,j) = ( qw(i,j,k) <= qw(i,j,k-1) ) + else + ! inversion moving down so inversion will moisten/dry + ! depending on relative QW in level above + moisten(i,j) = ( qw(i,j,k) <= qw(i,j,k+1) ) + end if - if ( moisten(i,j) ) then - ! Ensure inversion level does moisten relative to ML + if ( moisten(i,j) ) then + ! Ensure inversion level does moisten relative to ML - if (l_wtrac .and. totqf_efl < (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) ) & - totqf_efl_meth1(i,j) = 1 ! Store method + if (l_wtrac .and. totqf_efl < (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) ) & + totqf_efl_meth1(i,j) = 1 ! Store method - totqf_efl = max( totqf_efl, & - fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k) ) - if (we_parm(i,j)+w_ls(i,j) >= zero) then - ! Ensure inversion level won't end up more moist than - ! NTML by end of timestep. - ! Set INV_TEND to max allowable moistening rate, also - ! allowing for change in ML_TEND arising from this change - ! to TOTQF_EFL: - inv_tend = (qw(i,j,k-1)-qw(i,j,k))/timestep & - + (fq_nt(i,j,1)+fqw(i,j,1))/z_uv(i,j,k) - - if (l_wtrac .and. totqf_efl > & - ( (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & - /(one+ dzl(i,j,k)/z_uv(i,j,k)) ) ) & - totqf_efl_meth2(i,j) = 1 ! Store method - - totqf_efl = min( totqf_efl, & - (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & - /(one+ dzl(i,j,k)/z_uv(i,j,k)) ) - end if - else - if (l_wtrac .and. totqf_efl > (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) ) & - totqf_efl_meth1(i,j) = 1 ! Store method + totqf_efl = max( totqf_efl, & + fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k) ) + if (we_parm(i,j)+w_ls(i,j) >= zero) then + ! Ensure inversion level won't end up more moist than + ! NTML by end of timestep. + ! Set INV_TEND to max allowable moistening rate, also + ! allowing for change in ML_TEND arising from this change + ! to TOTQF_EFL: + inv_tend = (qw(i,j,k-1)-qw(i,j,k))/timestep & + + (fq_nt(i,j,1)+fqw(i,j,1))/z_uv(i,j,k) + + if (l_wtrac .and. totqf_efl > & + ( (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & + /(one+ dzl(i,j,k)/z_uv(i,j,k)) ) ) & + totqf_efl_meth2(i,j) = 1 ! Store method + + totqf_efl = min( totqf_efl, & + (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & + /(one+ dzl(i,j,k)/z_uv(i,j,k)) ) + end if + else + if (l_wtrac .and. totqf_efl > (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) ) & + totqf_efl_meth1(i,j) = 1 ! Store method - totqf_efl = min( totqf_efl, & - fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k) ) - if (we_parm(i,j)+w_ls(i,j) >= zero) then - ! Ensure inversion level won't end up drier than - ! NTML by end of timestep. - ! Set INV_TEND to max allowable drying rate: - inv_tend = (qw(i,j,k-1)-qw(i,j,k))/timestep & - + (fq_nt(i,j,1)+fqw(i,j,1))/z_uv(i,j,k) - - if (l_wtrac .and. totqf_efl < & - ( (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & - /(one+ dzl(i,j,k)/z_uv(i,j,k)) ) ) & - totqf_efl_meth2(i,j) = 1 ! Store method - - totqf_efl = max( totqf_efl, & - (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & - /(one+ dzl(i,j,k)/z_uv(i,j,k)) ) - end if + totqf_efl = min( totqf_efl, & + fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k) ) + if (we_parm(i,j)+w_ls(i,j) >= zero) then + ! Ensure inversion level won't end up drier than + ! NTML by end of timestep. + ! Set INV_TEND to max allowable drying rate: + inv_tend = (qw(i,j,k-1)-qw(i,j,k))/timestep & + + (fq_nt(i,j,1)+fqw(i,j,1))/z_uv(i,j,k) + + if (l_wtrac .and. totqf_efl < & + ( (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & + /(one+ dzl(i,j,k)/z_uv(i,j,k)) ) ) & + totqf_efl_meth2(i,j) = 1 ! Store method + + totqf_efl = max( totqf_efl, & + (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & + /(one+ dzl(i,j,k)/z_uv(i,j,k)) ) end if - fqw(i,j,k) = t_frac(i,j) * & - ( totqf_efl - fq_nt(i,j,k) ) end if + fqw(i,j,k) = t_frac(i,j) * & + ( totqf_efl - fq_nt(i,j,k) ) + end if - end do end do !$OMP end do +!end do + !$OMP end PARALLEL ! Repeat the last block of code for water tracers @@ -4571,14 +4671,16 @@ subroutine kmkhz_9c ( & ! call to calc_fqw_inv_wtrac allocate(z_uv_ntmlp1(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end)) -!$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(SHARED) private(i,j,k) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - k = ntml(i,j) + 1 - z_uv_ntmlp1(i,j) = z_uv(i,j,ntml(i,j)+1) - end do + + !do j = pdims%j_start, pdims%j_end + !$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(SHARED) private(i,k) + do i = pdims%i_start, pdims%i_end + k = ntml(i,j) + 1 + z_uv_ntmlp1(i,j) = z_uv(i,j,ntml(i,j)+1) end do -!$OMP end PARALLEL do + !$OMP end PARALLEL do + !end do + call calc_fqw_inv_wtrac(bl_levels, ntml, totqf_efl_meth1, & totqf_efl_meth2, t_frac, zh, zh_frac, & @@ -4593,95 +4695,93 @@ subroutine kmkhz_9c ( & !----------------------------------------------------------------------- ! Now decoupled layer !----------------------------------------------------------------------- - +!do j = pdims%j_start, pdims%j_end !$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(SHARED) & -!$OMP private ( i, j, k, totqf_efl, ml_tend, fa_tend, inv_tend) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - - moisten(i,j) = .false. - if (l_wtrac) totqf_efl_meth1(i,j) = 0 - if (l_wtrac) totqf_efl_meth2(i,j) = 0 - if ( t_frac_dsc(i,j) > zero ) then - - k = ntdsc(i,j)+1 +!$OMP private ( i, k, totqf_efl, ml_tend, fa_tend, inv_tend) +do i = pdims%i_start, pdims%i_end + moisten(i,j) = .false. + if (l_wtrac) totqf_efl_meth1(i,j) = 0 + if (l_wtrac) totqf_efl_meth2(i,j) = 0 + if ( t_frac_dsc(i,j) > zero ) then + + k = ntdsc(i,j)+1 + + ! Calculate total (turb+micro) QW flux at subgrid inversion + totqf_zhsc(i,j) = - we_rho_dsc(i,j)*dqw_dsc(i,j) & + + fq_nt_zhsc(i,j) + ! Interpolate to entrainment flux-level + totqf_efl = fq_nt_dscb(i,j) + & + ( totqf_zhsc(i,j) - fq_nt_dscb(i,j) )*zrzi_dsc(i,j) + + ml_tend = - ( totqf_zhsc(i,j)-fq_nt_dscb(i,j) )/dscdepth(i,j) + fa_tend = zero + if ( k+1 <= bl_levels ) & + fa_tend = - ( fq_nt(i,j,k+2)-fq_nt(i,j,k+1) ) & + / dzl(i,j,k+1) + inv_tend = zhsc_frac(i,j) * ml_tend & + + (one-zhsc_frac(i,j)) * fa_tend - ! Calculate total (turb+micro) QW flux at subgrid inversion - totqf_zhsc(i,j) = - we_rho_dsc(i,j)*dqw_dsc(i,j) & - + fq_nt_zhsc(i,j) - ! Interpolate to entrainment flux-level - totqf_efl = fq_nt_dscb(i,j) + & - ( totqf_zhsc(i,j) - fq_nt_dscb(i,j) )*zrzi_dsc(i,j) - - ml_tend = - ( totqf_zhsc(i,j)-fq_nt_dscb(i,j) )/dscdepth(i,j) - fa_tend = zero - if ( k+1 <= bl_levels ) & - fa_tend = - ( fq_nt(i,j,k+2)-fq_nt(i,j,k+1) ) & - / dzl(i,j,k+1) - inv_tend = zhsc_frac(i,j) * ml_tend & - + (one-zhsc_frac(i,j)) * fa_tend - - if (we_dsc_parm(i,j)+w_ls_dsc(i,j) >= zero) then - ! inversion moving up so inversion will moisten/dry - ! depending on relative QW in level below - moisten(i,j) = ( qw(i,j,k) <= qw(i,j,k-1) ) - else - ! inversion moving down so inversion will moisten/dry - ! depending on relative QW in level above - moisten(i,j) = ( qw(i,j,k) <= qw(i,j,k+1) ) - end if + if (we_dsc_parm(i,j)+w_ls_dsc(i,j) >= zero) then + ! inversion moving up so inversion will moisten/dry + ! depending on relative QW in level below + moisten(i,j) = ( qw(i,j,k) <= qw(i,j,k-1) ) + else + ! inversion moving down so inversion will moisten/dry + ! depending on relative QW in level above + moisten(i,j) = ( qw(i,j,k) <= qw(i,j,k+1) ) + end if - if ( moisten(i,j) ) then + if ( moisten(i,j) ) then - if (l_wtrac .and. (totqf_efl < fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) ) & - totqf_efl_meth1(i,j) = 1 ! Store method + if (l_wtrac .and. (totqf_efl < fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) ) & + totqf_efl_meth1(i,j) = 1 ! Store method - totqf_efl = max( totqf_efl, & - fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k) ) - if (we_dsc_parm(i,j)+w_ls_dsc(i,j) >= zero) then - ! Ensure inversion level won't end up more moist than - ! NTDSC by end of timestep. - inv_tend = (qw(i,j,k-1)-qw(i,j,k))/timestep & - + fq_nt_dscb(i,j)/dscdepth(i,j) + totqf_efl = max( totqf_efl, & + fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k) ) + if (we_dsc_parm(i,j)+w_ls_dsc(i,j) >= zero) then + ! Ensure inversion level won't end up more moist than + ! NTDSC by end of timestep. + inv_tend = (qw(i,j,k-1)-qw(i,j,k))/timestep & + + fq_nt_dscb(i,j)/dscdepth(i,j) - if (l_wtrac .and. totqf_efl > & - ( (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & - /(one+ dzl(i,j,k)/dscdepth(i,j)) ) ) & - totqf_efl_meth2(i,j) = 1 ! Store method + if (l_wtrac .and. totqf_efl > & + ( (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & + /(one+ dzl(i,j,k)/dscdepth(i,j)) ) ) & + totqf_efl_meth2(i,j) = 1 ! Store method - totqf_efl = min( totqf_efl, & - (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & - /(one+ dzl(i,j,k)/dscdepth(i,j)) ) - end if - else - if (l_wtrac .and. totqf_efl > (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) ) & - totqf_efl_meth1(i,j) = 1 ! Store method - - totqf_efl = min( totqf_efl, & - fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k) ) - if (we_dsc_parm(i,j)+w_ls_dsc(i,j) >= zero) then - ! Ensure inversion level won't end up drier than - ! NTDSC by end of timestep. - ! Set INV_TEND to max allowable drying rate: - inv_tend = (qw(i,j,k-1)-qw(i,j,k))/timestep & - + fq_nt_dscb(i,j)/dscdepth(i,j) - - if (l_wtrac .and. totqf_efl < & - ( (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & - /(one+ dzl(i,j,k)/dscdepth(i,j))) ) & - totqf_efl_meth2(i,j) = 1 ! Store method - - totqf_efl = max( totqf_efl, & - (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & - /(one+ dzl(i,j,k)/dscdepth(i,j)) ) - end if + totqf_efl = min( totqf_efl, & + (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & + /(one+ dzl(i,j,k)/dscdepth(i,j)) ) end if - fqw(i,j,k) = t_frac_dsc(i,j) * ( totqf_efl - fq_nt(i,j,k) ) + else + if (l_wtrac .and. totqf_efl > (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) ) & + totqf_efl_meth1(i,j) = 1 ! Store method + totqf_efl = min( totqf_efl, & + fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k) ) + if (we_dsc_parm(i,j)+w_ls_dsc(i,j) >= zero) then + ! Ensure inversion level won't end up drier than + ! NTDSC by end of timestep. + ! Set INV_TEND to max allowable drying rate: + inv_tend = (qw(i,j,k-1)-qw(i,j,k))/timestep & + + fq_nt_dscb(i,j)/dscdepth(i,j) + + if (l_wtrac .and. totqf_efl < & + ( (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & + /(one+ dzl(i,j,k)/dscdepth(i,j))) ) & + totqf_efl_meth2(i,j) = 1 ! Store method + + totqf_efl = max( totqf_efl, & + (fq_nt(i,j,k+1)+inv_tend*dzl(i,j,k)) & + /(one+ dzl(i,j,k)/dscdepth(i,j)) ) + end if end if - end do + fqw(i,j,k) = t_frac_dsc(i,j) * ( totqf_efl - fq_nt(i,j,k) ) + + end if end do !$OMP end PARALLEL do +!end do ! Repeat last block of code for water tracers if (l_wtrac) then @@ -4695,81 +4795,82 @@ subroutine kmkhz_9c ( & end if ! l_wtrac !$OMP PARALLEL DEFAULT(SHARED) & -!$OMP private (i, j, k, kp, w_var_inv, weight, tke_nl_rh, delta_tke, & +!$OMP private (i, k, kp, w_var_inv, weight, tke_nl_rh, delta_tke, & !$OMP w_s_ent, w_s_cubed, w_m, wstar3, w_h) - !----------------------------------------------------------------------- ! Estimate turbulent w-variance scale at discontinuous inversions !----------------------------------------------------------------------- if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ! We expect the turbulent temperature perturbation scale close - ! to the inversion, Tl' ~ ftl / w_scale, to be at most 1/2 dsl_inv, - ! where ftl is the entrainment heat flux, and w_scale = sqrt(w_var) - ! => w_scale >= 2 ftl / dsl_inv - ! and then square that to get w_var. - - ! Want to set non-local TKE to max of existing value and this - ! w_scale^2, but this is complicated by the fact that - ! tke_nl is on theta-levels, while w_scale is on rho-levels. - ! Compare w_scale^2 with the existing non-local TKE - ! interpolated to the same rho-level; increase TKE on the neighbouring - ! theta-levels if needed to increase interpolated rho-level value - ! to the desired inversion value. - - ! Note: heights z_tq are on actual theta-levels (surface at k=0), - ! but tke_nl is offset (surface at k=1) - - ! Do the above for discontinuous inversion at top of surface mixed-layer - if ( sml_disc_inv(i,j) >= 1 ) then - k = ntml(i,j) + 1 - kp = min( k+1, bl_levels ) - - w_var_inv = 2.0_r_bl * (ftl(i,j,k)/rho_mix(i,j,k)) / dsl_sml(i,j) - w_var_inv = w_var_inv * w_var_inv * rho_mix(i,j,k) - - weight = ( z_uv(i,j,k) - z_tq(i,j,k-1) ) & - / ( z_tq(i,j,k) - z_tq(i,j,k-1) ) - tke_nl_rh = (one-weight) * tke_nl(i,j,k) & - + weight * tke_nl(i,j,kp) - - delta_tke = w_var_inv - tke_nl_rh - if ( delta_tke > zero ) then - tke_nl(i,j,k) = tke_nl(i,j,k) + delta_tke - tke_nl(i,j,kp) = tke_nl(i,j,kp) + delta_tke - end if + ! do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + + ! We expect the turbulent temperature perturbation scale close + ! to the inversion, Tl' ~ ftl / w_scale, to be at most 1/2 dsl_inv, + ! where ftl is the entrainment heat flux, and w_scale = sqrt(w_var) + ! => w_scale >= 2 ftl / dsl_inv + ! and then square that to get w_var. + + ! Want to set non-local TKE to max of existing value and this + ! w_scale^2, but this is complicated by the fact that + ! tke_nl is on theta-levels, while w_scale is on rho-levels. + ! Compare w_scale^2 with the existing non-local TKE + ! interpolated to the same rho-level; increase TKE on the neighbouring + ! theta-levels if needed to increase interpolated rho-level value + ! to the desired inversion value. + + ! Note: heights z_tq are on actual theta-levels (surface at k=0), + ! but tke_nl is offset (surface at k=1) + + ! Do the above for discontinuous inversion at top of surface mixed-layer + if ( sml_disc_inv(i,j) >= 1 ) then + k = ntml(i,j) + 1 + kp = min( k+1, bl_levels ) + + w_var_inv = 2.0_r_bl * (ftl(i,j,k)/rho_mix(i,j,k)) / dsl_sml(i,j) + w_var_inv = w_var_inv * w_var_inv * rho_mix(i,j,k) + weight = ( z_uv(i,j,k) - z_tq(i,j,k-1) ) & + / ( z_tq(i,j,k) - z_tq(i,j,k-1) ) + tke_nl_rh = (one-weight) * tke_nl(i,j,k) & + + weight * tke_nl(i,j,kp) + + delta_tke = w_var_inv - tke_nl_rh + if ( delta_tke > zero ) then + tke_nl(i,j,k) = tke_nl(i,j,k) + delta_tke + tke_nl(i,j,kp) = tke_nl(i,j,kp) + delta_tke end if - ! Repeat for discontinuous inversion at top of decoupled Sc layer - if ( dsc_disc_inv(i,j) >= 1 ) then - k = ntdsc(i,j) + 1 - kp = min( k+1, bl_levels ) + end if - w_var_inv = 2.0_r_bl * (ftl(i,j,k)/rho_mix(i,j,k)) / dsl_dsc(i,j) - w_var_inv = w_var_inv * w_var_inv * rho_mix(i,j,k) + ! Repeat for discontinuous inversion at top of decoupled Sc layer + if ( dsc_disc_inv(i,j) >= 1 ) then + k = ntdsc(i,j) + 1 + kp = min( k+1, bl_levels ) - weight = ( z_uv(i,j,k) - z_tq(i,j,k-1) ) & - / ( z_tq(i,j,k) - z_tq(i,j,k-1) ) - tke_nl_rh = (one-weight) * tke_nl(i,j,k) & - + weight * tke_nl(i,j,kp) + w_var_inv = 2.0_r_bl * (ftl(i,j,k)/rho_mix(i,j,k)) / dsl_dsc(i,j) + w_var_inv = w_var_inv * w_var_inv * rho_mix(i,j,k) - delta_tke = w_var_inv - tke_nl_rh + weight = ( z_uv(i,j,k) - z_tq(i,j,k-1) ) & + / ( z_tq(i,j,k) - z_tq(i,j,k-1) ) + tke_nl_rh = (one-weight) * tke_nl(i,j,k) & + + weight * tke_nl(i,j,kp) - if ( delta_tke > zero ) then - tke_nl(i,j,k) = tke_nl(i,j,k) + delta_tke - tke_nl(i,j,kp) = tke_nl(i,j,kp) + delta_tke - end if + delta_tke = w_var_inv - tke_nl_rh + if ( delta_tke > zero ) then + tke_nl(i,j,k) = tke_nl(i,j,k) + delta_tke + tke_nl(i,j,kp) = tke_nl(i,j,kp) + delta_tke end if - end do + end if + end do -!$OMP end do NOWAIT + !$OMP end do NOWAIT + ! end do + end if ! (BL_diag%l_tke) !----------------------------------------------------------------------- @@ -4781,119 +4882,125 @@ subroutine kmkhz_9c ( & ! grid-level so only one element of these 3D arrays is used. !----------------------------------------------------------------------- + +! 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 - we_lim(i,j,1) = zero - t_frac_tr(i,j,1) = zero - zrzi_tr(i,j,1) = zero - we_lim_dsc(i,j,1) = zero - t_frac_dsc_tr(i,j,1) = zero - zrzi_dsc_tr(i,j,1) = zero - we_lim(i,j,3) = zero - t_frac_tr(i,j,3) = zero - zrzi_tr(i,j,3) = zero - we_lim_dsc(i,j,3) = zero - t_frac_dsc_tr(i,j,3) = zero - zrzi_dsc_tr(i,j,3) = zero - end do ! i -end do ! j +do i = pdims%i_start, pdims%i_end + we_lim(i,j,1) = zero + t_frac_tr(i,j,1) = zero + zrzi_tr(i,j,1) = zero + we_lim_dsc(i,j,1) = zero + t_frac_dsc_tr(i,j,1) = zero + zrzi_dsc_tr(i,j,1) = zero + we_lim(i,j,3) = zero + t_frac_tr(i,j,3) = zero + zrzi_tr(i,j,3) = zero + we_lim_dsc(i,j,3) = zero + t_frac_dsc_tr(i,j,3) = zero + zrzi_dsc_tr(i,j,3) = zero +end do ! i !$OMP end do NOWAIT +! end do ! j + + +! 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 - kent(i,j) = ntml(i,j)+1 - t_frac_tr(i,j,2) = t_frac(i,j) - zrzi_tr(i,j,2) = zrzi(i,j) - if ( t_frac(i,j) > zero ) then - w_s_ent = zero - k = ntml(i,j) - if ( abs( dsl_sml(i,j) ) >= rbl_eps ) w_s_ent = & - min( zero, -sls_inc(i,j,k) * dzl(i,j,k) /dsl_sml(i,j) ) - ! Only allow w_e to be reduced to zero! - we_lim(i,j,2) = rho_mix(i,j,k+1) * & - max( zero, we_parm(i,j) + w_s_ent ) - else - we_lim(i,j,2) = zero - end if - kent_dsc(i,j) = ntdsc(i,j)+1 - t_frac_dsc_tr(i,j,2) = t_frac_dsc(i,j) - zrzi_dsc_tr(i,j,2) = zrzi_dsc(i,j) - if ( t_frac_dsc(i,j) > zero ) then - w_s_ent = zero - k = ntdsc(i,j) - if ( abs( dsl_dsc(i,j) ) >= rbl_eps ) w_s_ent = & - min( zero, -sls_inc(i,j,k) * dzl(i,j,k) /dsl_dsc(i,j) ) - ! Only allow w_e to be reduced to zero! - we_lim_dsc(i,j,2) = rho_mix(i,j,k) * & - max( zero, we_dsc_parm(i,j) + w_s_ent ) - else - we_lim_dsc(i,j,2) = zero - end if - end do +do i = pdims%i_start, pdims%i_end + kent(i,j) = ntml(i,j)+1 + t_frac_tr(i,j,2) = t_frac(i,j) + zrzi_tr(i,j,2) = zrzi(i,j) + if ( t_frac(i,j) > zero ) then + w_s_ent = zero + k = ntml(i,j) + if ( abs( dsl_sml(i,j) ) >= rbl_eps ) w_s_ent = & + min( zero, -sls_inc(i,j,k) * dzl(i,j,k) /dsl_sml(i,j) ) + ! Only allow w_e to be reduced to zero! + we_lim(i,j,2) = rho_mix(i,j,k+1) * & + max( zero, we_parm(i,j) + w_s_ent ) + else + we_lim(i,j,2) = zero + end if + kent_dsc(i,j) = ntdsc(i,j)+1 + t_frac_dsc_tr(i,j,2) = t_frac_dsc(i,j) + zrzi_dsc_tr(i,j,2) = zrzi_dsc(i,j) + if ( t_frac_dsc(i,j) > zero ) then + w_s_ent = zero + k = ntdsc(i,j) + if ( abs( dsl_dsc(i,j) ) >= rbl_eps ) w_s_ent = & + min( zero, -sls_inc(i,j,k) * dzl(i,j,k) /dsl_dsc(i,j) ) + ! Only allow w_e to be reduced to zero! + we_lim_dsc(i,j,2) = rho_mix(i,j,k) * & + max( zero, we_dsc_parm(i,j) + w_s_ent ) + else + we_lim_dsc(i,j,2) = zero + end if end do !$OMP end do NOWAIT +! end do + !----------------------------------------------------------------------- ! 12. Update standard deviations and gradient adjustment to use this ! timestep's ZH (code from SF_EXCH) !----------------------------------------------------------------------- + +! 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 ( unstable(i,j) ) then - if (flux_grad == Locketal2000) then - w_s_cubed = 0.25_r_bl * zh(i,j) * fb_surf(i,j) - if (w_s_cubed > zero) then - w_m = & - ( w_s_cubed + v_s(i,j) * v_s(i,j) * v_s(i,j) ) ** one_third - t1_sd(i,j) = 1.93_r_bl * ftl(i,j,1) / (rhostar_gb(i,j) * w_m) - q1_sd(i,j) = 1.93_r_bl * fqw(i,j,1) / (rhostar_gb(i,j) * w_m) - tv1_sd(i,j) = t(i,j,1) * & - ( one + c_virtual*q(i,j,1) - qcl(i,j,1) - qcf(i,j,1) ) * & - ( bt(i,j,1)*t1_sd(i,j) + bq(i,j,1)*q1_sd(i,j) ) - t1_sd(i,j) = max ( zero , t1_sd(i,j) ) - q1_sd(i,j) = max ( zero , q1_sd(i,j) ) - if (tv1_sd(i,j) <= zero) then - tv1_sd(i,j) = zero - t1_sd(i,j) = zero - q1_sd(i,j) = zero - end if +do i = pdims%i_start, pdims%i_end + if ( unstable(i,j) ) then + if (flux_grad == Locketal2000) then + w_s_cubed = 0.25_r_bl * zh(i,j) * fb_surf(i,j) + if (w_s_cubed > zero) then + w_m = & + ( w_s_cubed + v_s(i,j) * v_s(i,j) * v_s(i,j) ) ** one_third + t1_sd(i,j) = 1.93_r_bl * ftl(i,j,1) / (rhostar_gb(i,j) * w_m) + q1_sd(i,j) = 1.93_r_bl * fqw(i,j,1) / (rhostar_gb(i,j) * w_m) + tv1_sd(i,j) = t(i,j,1) * & + ( one + c_virtual*q(i,j,1) - qcl(i,j,1) - qcf(i,j,1) ) * & + ( bt(i,j,1)*t1_sd(i,j) + bq(i,j,1)*q1_sd(i,j) ) + t1_sd(i,j) = max ( zero , t1_sd(i,j) ) + q1_sd(i,j) = max ( zero , q1_sd(i,j) ) + if (tv1_sd(i,j) <= zero) then + tv1_sd(i,j) = zero + t1_sd(i,j) = zero + q1_sd(i,j) = zero end if - grad_t_adj(i,j) = min( max_t_grad , & - a_grad_adj * t1_sd(i,j) / zh(i,j) ) - grad_q_adj(i,j) = zero - else if (flux_grad == HoltBov1993) then - ! Use constants from Holtslag and Boville (1993) - ! Conv limit GAMMA_TH = 10 *FTL1/(wstar*zh) - ! Neut limit GAMMA_TH = 7.2*wstar*FTL1/(ustar^2*zh) - wstar3 = fb_surf(i,j) * zh(i,j) - w_m =( v_s(i,j)**3 + 0.6_r_bl*wstar3 )**one_third - - grad_t_adj(i,j) = a_ga_hb93*(wstar3**one_third)*ftl(i,j,1) & - / ( rhostar_gb(i,j)*w_m*w_m*zh(i,j) ) - ! GRAD_Q_ADJ(I,j) = A_GA_HB93*(WSTAR3**one_third)*FQW(I,j,1) - ! / ( RHOSTAR_GB(I,j)*W_M*W_M*ZH(I,j) ) - ! Set q term to zero for same empirical reasons as Lock et al - grad_q_adj(i,j) = zero - else if (flux_grad == LockWhelan2006) then - ! Use constants LockWhelan2006 - ! Conv limit GAMMA_TH = 10 *FTL1/(wstar*zh) - ! Neut limit GAMMA_TH = 7.5*FTL1/(ustar*zh) - wstar3 = fb_surf(i,j) * zh(i,j) - w_h =( ((4.0_r_bl/3.0_r_bl)*v_s(i,j))**3 + wstar3 )**one_third - - grad_t_adj(i,j) = a_ga_lw06 * ftl(i,j,1) & - / ( rhostar_gb(i,j)*w_h*zh(i,j) ) - grad_q_adj(i,j) = a_ga_lw06 * fqw(i,j,1) & - / ( rhostar_gb(i,j)*w_h*zh(i,j) ) end if - end if ! test on UNSTABLE - end do + grad_t_adj(i,j) = min( max_t_grad , & + a_grad_adj * t1_sd(i,j) / zh(i,j) ) + grad_q_adj(i,j) = zero + else if (flux_grad == HoltBov1993) then + ! Use constants from Holtslag and Boville (1993) + ! Conv limit GAMMA_TH = 10 *FTL1/(wstar*zh) + ! Neut limit GAMMA_TH = 7.2*wstar*FTL1/(ustar^2*zh) + wstar3 = fb_surf(i,j) * zh(i,j) + w_m =( v_s(i,j)**3 + 0.6_r_bl*wstar3 )**one_third + + grad_t_adj(i,j) = a_ga_hb93*(wstar3**one_third)*ftl(i,j,1) & + / ( rhostar_gb(i,j)*w_m*w_m*zh(i,j) ) + ! GRAD_Q_ADJ(I,j) = A_GA_HB93*(WSTAR3**one_third)*FQW(I,j,1) + ! / ( RHOSTAR_GB(I,j)*W_M*W_M*ZH(I,j) ) + ! Set q term to zero for same empirical reasons as Lock et al + grad_q_adj(i,j) = zero + else if (flux_grad == LockWhelan2006) then + ! Use constants LockWhelan2006 + ! Conv limit GAMMA_TH = 10 *FTL1/(wstar*zh) + ! Neut limit GAMMA_TH = 7.5*FTL1/(ustar*zh) + wstar3 = fb_surf(i,j) * zh(i,j) + w_h =( ((4.0_r_bl/3.0_r_bl)*v_s(i,j))**3 + wstar3 )**one_third + + grad_t_adj(i,j) = a_ga_lw06 * ftl(i,j,1) & + / ( rhostar_gb(i,j)*w_h*zh(i,j) ) + grad_q_adj(i,j) = a_ga_lw06 * fqw(i,j,1) & + / ( rhostar_gb(i,j)*w_h*zh(i,j) ) + end if + end if ! test on UNSTABLE end do !$OMP end do NOWAIT +! end do + ! (Note, water tracers assume flux_grad = Locketal2000 so no need to ! update wtrac_bl%grad_q_adj as it is always zero) @@ -4903,23 +5010,25 @@ subroutine kmkhz_9c ( & !----------------------------------------------------------------------- if (BL_diag%l_dzh) then -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ! fill unset values (rmdi<0) with zero - ! and in cumulus set dzh=zh-z_lcl - if (kprof_cu >= on .and. cumulus(i,j)) then - BL_diag%dzh(i,j)= max( zero, zh(i,j)-z_lcl(i,j) ) - else - BL_diag%dzh(i,j)= max( zero, dzh(i,j) ) - end if - end do + + !do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + ! fill unset values (rmdi<0) with zero + ! and in cumulus set dzh=zh-z_lcl + if (kprof_cu >= on .and. cumulus(i,j)) then + BL_diag%dzh(i,j)= max( zero, zh(i,j)-z_lcl(i,j) ) + else + BL_diag%dzh(i,j)= max( zero, dzh(i,j) ) + end if end do -!$OMP end do NOWAIT + !$OMP end do NOWAIT + !end do end if if (BL_diag%l_dscbase) then -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end + + ! do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end if ( dsc(i,j) ) then BL_diag%dscbase(i,j)= zhsc(i,j)-dscdepth(i,j) @@ -4927,43 +5036,50 @@ subroutine kmkhz_9c ( & BL_diag%dscbase(i,j)= rmdi end if end do - end do -!$OMP end do NOWAIT + !$OMP end do NOWAIT + ! end do + end if if (BL_diag%l_cldbase) then -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( dsc(i,j) ) then - BL_diag%cldbase(i,j)= zhsc(i,j)-zc_dsc(i,j) - else - BL_diag%cldbase(i,j)= zh(i,j)-zc(i,j) - end if - end do + + ! do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + if ( dsc(i,j) ) then + BL_diag%cldbase(i,j)= zhsc(i,j)-zc_dsc(i,j) + else + BL_diag%cldbase(i,j)= zh(i,j)-zc(i,j) + end if end do -!$OMP end do NOWAIT + !$OMP end do NOWAIT + ! end do + end if if (BL_diag%l_weparm_dsc) then -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( dsc(i,j) ) then - BL_diag%weparm_dsc(i,j)= we_dsc_parm(i,j) - else - BL_diag%weparm_dsc(i,j)= we_parm(i,j) - end if - end do + + ! do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + if ( dsc(i,j) ) then + BL_diag%weparm_dsc(i,j)= we_dsc_parm(i,j) + else + BL_diag%weparm_dsc(i,j)= we_parm(i,j) + end if end do -!$OMP end do NOWAIT + !$OMP end do NOWAIT + ! end do + end if if (BL_diag%l_weparm) then -!$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - BL_diag%weparm(i,j)= we_parm(i,j) - end do + + !do j = pdims%j_start, pdims%j_end + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + BL_diag%weparm(i,j)= we_parm(i,j) end do -!$OMP end do NOWAIT + !$OMP end do NOWAIT + !end do + end if !-----------------------------------------------------------------------