diff --git a/science/physics_schemes/source/boundary_layer/excf_nl_9c.F90 b/science/physics_schemes/source/boundary_layer/excf_nl_9c.F90 index 4525be05d..1246882cf 100644 --- a/science/physics_schemes/source/boundary_layer/excf_nl_9c.F90 +++ b/science/physics_schemes/source/boundary_layer/excf_nl_9c.F90 @@ -700,7 +700,7 @@ subroutine excf_nl_9c ( & integer :: i1, j1, l -integer :: jj !Cache blocking - loop index +integer :: jj, ii !Cache blocking - loop index(s) integer(kind=jpim), parameter :: zhook_in = 0 integer(kind=jpim), parameter :: zhook_out = 1 @@ -765,9 +765,9 @@ subroutine excf_nl_9c ( & !----------------------------------------------------------------------- ! 0. Calculate top-of-b.l. velocity scales and Prandtl number. !----------------------------------------------------------------------- - +j = 1 !$OMP PARALLEL DEFAULT(SHARED) & -!$OMP private(i, j, k, jj, i_wt, c_ws, wstar3, pr_neut, pr_conv, w_m_neut, & +!$OMP private(i, k, ii, jj, i_wt, c_ws, wstar3, pr_neut, pr_conv, w_m_neut, & !$OMP zeta_s_fac, sf_term, sf_shear_term, zeta_r_sq, ir_term, zr, & !$OMP evap_term, dz_inv, l_rad, alpha_t, dr_term, zil_corr, & !$OMP rhokh_ent, frac_top, zh_pr, factor, rhokm_dsct, & @@ -780,8 +780,10 @@ subroutine excf_nl_9c ( & !cdir collapse !$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 +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 rhokh_surf_ent(i,j) = zero rhokh_top_ent(i,j) = zero @@ -811,7 +813,7 @@ subroutine excf_nl_9c ( & if (coupled(i,j)) then if ( entr_smooth_dec == entr_taper_zh ) then wstar3 = ( svl_diff_frac(i,j) * zhsc(i,j) & - + (one-svl_diff_frac(i,j)) * zh(i,j) )* fb_surf(i,j) + + (one-svl_diff_frac(i,j)) * zh(i,j) )* fb_surf(i,j) else wstar3 = zhsc(i,j) * fb_surf(i,j) end if @@ -830,10 +832,10 @@ subroutine excf_nl_9c ( & ! Turbulent Prandtl number and velocity scale for scalars ! gives 0.375 0 @@ -963,20 +968,20 @@ subroutine excf_nl_9c ( & dz_inv = min( v_sum(i,j)*v_sum(i,j) / db_top(i,j) ,100.0_r_bl ) l_rad = 15.0_r_bl * max( one, 200.0_r_bl/(zc(i,j)+rbl_eps) ) alpha_t = one - exp(-one_half*dz_inv/l_rad) - ! Make enhancement due to buoyancy reversal feedback + ! Make enhancement due to buoyancy reversal feedback alpha_t = alpha_t + br_fback(i,j)*(one-alpha_t) dr_term = btc_top(i,j) * alpha_t * df_top_over_cp(i,j) - !---------------------------------------------------------- - ! Combine terms to calculate the entrainment - ! mixing coefficients - !---------------------------------------------------------- + !---------------------------------------------------------- + ! Combine terms to calculate the entrainment + ! mixing coefficients + !---------------------------------------------------------- zil_corr = c_t * ( (sf_term + sf_shear_term + & ir_term + evap_term) / & (rho_mix(i,j,k) * sqrt(zh(i,j))) )**two_thirds rho_we_sml(i,j) = (sf_term + sf_shear_term & + ir_term + evap_term + dr_term) & - / ( db_top(i,j) + zil_corr ) + / ( db_top(i,j) + zil_corr ) rhokh_ent = rho_we_sml(i,j)/ rdz(i,j,k) @@ -989,8 +994,8 @@ subroutine excf_nl_9c ( & ! flux at the top and then added the entrainment flux onto ! this, which is double-counting). l_apply_surf_ent = ( .not. ( ( kprof_cu==buoy_integ .or. & - kprof_cu==buoy_integ_low ) .and. & - cumulus(i,j) ) ) + kprof_cu==buoy_integ_low ) .and. & + cumulus(i,j) ) ) else ! If bug-fix is off, always add on surface-driven entrainment flux l_apply_surf_ent = .true. @@ -1004,8 +1009,8 @@ subroutine excf_nl_9c ( & end if rhokh_top_ent(i,j) = rhokh_ent * frac_top rhokm_top_ent(i,j) = 0.75_r_bl * rhokh_top_ent(i,j) & - * rdz(i,j,k) * (z_uv(i,j,k)-z_uv(i,j,k-1)) & - * rho_wet_tq(i,j,k-1) / rho_mix(i,j,k) + * rdz(i,j,k) * (z_uv(i,j,k)-z_uv(i,j,k-1)) & + * rho_wet_tq(i,j,k-1) / rho_mix(i,j,k) if (.not. l_use_var_fixes) then rhokm(i,j,k) = rhokm_surf_ent(i,j) rhokm_top(i,j,k) = rhokm_top_ent(i,j) @@ -1071,28 +1076,28 @@ subroutine excf_nl_9c ( & ir_term + evap_term) & * dscdepth(i,j) / (a_ent_1*rho_mix(i,j,k)) )**one_third v_top_dsc(i,j) =( (ir_term + evap_term) * dscdepth(i,j) / & - (a_ent_1*rho_mix(i,j,k)) )**one_third + (a_ent_1*rho_mix(i,j,k)) )**one_third !----------------------------------------------------------- ! Calculate the direct radiative term !----------------------------------------------------------- if (db_dsct(i,j) > zero) then dz_inv = min( v_sum_dsc(i,j)*v_sum_dsc(i,j)/db_dsct(i,j), & - 100.0_r_bl ) + 100.0_r_bl ) l_rad = 15.0_r_bl * max( one, 200.0_r_bl/(zc_dsc(i,j)+one) ) alpha_t = one - exp(-one_half*dz_inv/l_rad) - ! Make enhancement due to buoyancy reversal feedback + ! Make enhancement due to buoyancy reversal feedback alpha_t = alpha_t + br_fback_dsc(i,j)*(one-alpha_t) dr_term = btc_top(i,j) * alpha_t * df_dsct_over_cp(i,j) - !---------------------------------------------------------- - ! Finally combine terms to calculate the entrainment - ! rate and mixing coefficients - !---------------------------------------------------------- + !---------------------------------------------------------- + ! Finally combine terms to calculate the entrainment + ! rate and mixing coefficients + !---------------------------------------------------------- zil_corr = c_t * ( (sf_term + sf_shear_term + & ir_term + evap_term) / & - (rho_mix(i,j,k) * sqrt(dscdepth(i,j))) )**two_thirds + (rho_mix(i,j,k) * sqrt(dscdepth(i,j))) )**two_thirds rho_we_dsc(i,j) = ( sf_term + sf_shear_term & + ir_term + evap_term + dr_term ) & - / ( db_dsct(i,j) + zil_corr ) + / ( db_dsct(i,j) + zil_corr ) rhokh_dsct_ent(i,j) = rho_we_dsc(i,j)/ rdz(i,j,k) rhokm_dsct_ent(i,j) = 0.75_r_bl * rhokh_dsct_ent(i,j) & * rdz(i,j,k) * (z_uv(i,j,k)-z_uv(i,j,k-1)) & @@ -1100,27 +1105,28 @@ subroutine excf_nl_9c ( & if (.not. l_use_var_fixes) rhokm_top(i,j,k) = rhokm_dsct_ent(i,j) end if ! test on DB_DSCT gt 0 end if - end do -end do + end do ! I +end do ! II +!end do !$OMP end do ! If there is no turbulence generation in DSC layer, ignore it. !cdir collapse !$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. v_top_dsc(i,j) <= zero ) then - dsc_removed(i,j) = 1 - dsc(i,j) = .false. - ntdsc(i,j) = 0 - zhsc(i,j) = zero - zc_dsc(i,j) = zero - dscdepth(i,j) = zero - coupled(i,j) = .false. - end if - end do +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + if ( dsc(i,j) .and. v_top_dsc(i,j) <= zero ) then + dsc_removed(i,j) = 1 + dsc(i,j) = .false. + ntdsc(i,j) = 0 + zhsc(i,j) = zero + zc_dsc(i,j) = zero + dscdepth(i,j) = zero + coupled(i,j) = .false. + end if end do +!end do !$OMP end do ! ---------------------------------------------------------------------- @@ -1144,13 +1150,13 @@ subroutine excf_nl_9c ( & ! Need to use pdims%k_end(=bl_levels-1) here because bl_levels ! in this routine is actually nl_bl_levels (< bl_levels) do k = pdims%k_start, pdims%k_end - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - wbmix(i,j,k) = zero ! WB if were diag as well-mixed - wbend(i,j,k) = zero ! WB after dec diag - wbend_sml(i,j,k) = zero ! WB after ksurf_iterate - end do + ! do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + wbmix(i,j,k) = zero ! WB if were diag as well-mixed + wbend(i,j,k) = zero ! WB after dec diag + wbend_sml(i,j,k) = zero ! WB after ksurf_iterate end do + ! end do end do !$OMP end do end if ! model_type @@ -1160,40 +1166,40 @@ subroutine excf_nl_9c ( & !cdir collapse !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - zsml_top(i,j) = zh(i,j) - zsml_base(i,j) = 0.1_r_bl * zh(i,j) - zdsc_base(i,j) = 0.1_r_bl * zhsc(i,j) - z_top_lim(i,j) = zero ! initialising - z_inv(i,j) = zero ! inversion height (top of K profiles) - zwb0(i,j) = zero ! height at which WB goes to zero - wbp_int(i,j) = zero - wbn_int(i,j) = zero - wb_surf_int(i,j) = zero - wb_dzrad_int(i,j) = zero - dzrad(i,j) = 100.0_r_bl - tothf_zi(i,j) = zero - totqf_zi(i,j) = zero - kstatus(i,j)= .true. - kwb0(i,j) = 2 - ntop(i,j) = -1 - ksurf(i,j) = 1 - dec_thres(i,j) = dec_thres_cloud ! use cloudy by default - end do +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + zsml_top(i,j) = zh(i,j) + zsml_base(i,j) = 0.1_r_bl * zh(i,j) + zdsc_base(i,j) = 0.1_r_bl * zhsc(i,j) + z_top_lim(i,j) = zero ! initialising + z_inv(i,j) = zero ! inversion height (top of K profiles) + zwb0(i,j) = zero ! height at which WB goes to zero + wbp_int(i,j) = zero + wbn_int(i,j) = zero + wb_surf_int(i,j) = zero + wb_dzrad_int(i,j) = zero + dzrad(i,j) = 100.0_r_bl + tothf_zi(i,j) = zero + totqf_zi(i,j) = zero + kstatus(i,j)= .true. + kwb0(i,j) = 2 + ntop(i,j) = -1 + ksurf(i,j) = 1 + dec_thres(i,j) = dec_thres_cloud ! use cloudy by default end do +!end do !$OMP end do ! Find KSURF, the first theta-level above the surface layer !$OMP do SCHEDULE(STATIC) -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 !cdir collapse - do j = jj, min(jj+bl_segment_size-1, pdims%j_end) - do i = pdims%i_start, pdims%i_end - if ( z_tq(i,j,k-1) < 0.1_r_bl*zh(i,j) ) ksurf(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 + if ( z_tq(i,j,k-1) < 0.1_r_bl*zh(i,j) ) ksurf(i,j) = k + !end do end do end do end do @@ -1205,22 +1211,24 @@ subroutine excf_nl_9c ( & !cdir collapse !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - test_well_mixed(i,j) = .false. - ksurf_iterate(i,j)= .false. - ktop_iterate(i,j) = .false. - if ( ntdsc(i,j) > 2 ) then - ktop_iterate(i,j) = .true. - end if - end do ! I -end do ! J +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + test_well_mixed(i,j) = .false. + ksurf_iterate(i,j)= .false. + ktop_iterate(i,j) = .false. + if ( ntdsc(i,j) > 2 ) then + ktop_iterate(i,j) = .true. + end if +end do ! I +!end do ! J !$OMP end do !cdir collapse !$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 +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 ( bflux_surf(i,j) > zero) then ! can only be coupled to an unstable SML if ( ntdsc(i,j) > 2 ) then @@ -1271,51 +1279,52 @@ subroutine excf_nl_9c ( & end if end if end do ! I -end do ! J +end do ! II +!end do ! J !$OMP end do if ( l_converge_ga ) then ! Set gradient adjustment terms consistent with zh = z_inv !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( test_well_mixed(i,j) ) then - ! The gradient adjustment scales with 1/zh. The terms were calculated - ! in kmkhz_9c using zh_prev (mixed-layer height from previous timestep), - ! but here we are setting the mixed-layer height to z_inv, - ! so need to scale the gradient adjustment by zh_prev/z_inv - ga_fac(i,j) = zh_prev(i,j) / z_inv(i,j) - else - ! Set to 1 where not used. - ga_fac(i,j) = one - end if - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + if ( test_well_mixed(i,j) ) then + ! The gradient adjustment scales with 1/zh. The terms were calculated + ! in kmkhz_9c using zh_prev (mixed-layer height from previous timestep), + ! but here we are setting the mixed-layer height to z_inv, + ! so need to scale the gradient adjustment by zh_prev/z_inv + ga_fac(i,j) = zh_prev(i,j) / z_inv(i,j) + else + ! Set to 1 where not used. + ga_fac(i,j) = one + end if end do + !end do !$OMP end do !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ! Updated gradient-adjusted db is non-adjusted value + ga_fac times - ! the original adjustment. - db_ga_dry_n(i,j,k) = db_noga_dry(i,j,k) & - + ( db_ga_dry(i,j,k) - db_noga_dry(i,j,k) ) * ga_fac(i,j) - db_ga_cld_n(i,j,k) = db_noga_cld(i,j,k) & - + ( db_ga_cld(i,j,k) - db_noga_cld(i,j,k) ) * ga_fac(i,j) - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + ! Updated gradient-adjusted db is non-adjusted value + ga_fac times + ! the original adjustment. + db_ga_dry_n(i,j,k) = db_noga_dry(i,j,k) & + + ( db_ga_dry(i,j,k) - db_noga_dry(i,j,k) ) * ga_fac(i,j) + db_ga_cld_n(i,j,k) = db_noga_cld(i,j,k) & + + ( db_ga_cld(i,j,k) - db_noga_cld(i,j,k) ) * ga_fac(i,j) end do + ! end do end do !$OMP end do else ! ( .not. l_converge_ga ) ! Use original gradient adjustment set consistent with zh_prev !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - db_ga_dry_n(i,j,k) = db_ga_dry(i,j,k) - db_ga_cld_n(i,j,k) = db_ga_cld(i,j,k) - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + db_ga_dry_n(i,j,k) = db_ga_dry(i,j,k) + db_ga_cld_n(i,j,k) = db_ga_cld(i,j,k) end do + !end do end do !$OMP end do end if ! ( l_converge_ga ) @@ -1323,18 +1332,18 @@ subroutine excf_nl_9c ( & ! Find kwb0, level with lowest positive cloud-free buoyancy gradient !$OMP do SCHEDULE(STATIC) -do jj = pdims%j_start, pdims%j_end, bl_segment_size +do ii = pdims%j_start, pdims%i_end, bl_segment_size do k = 2, bl_levels - do j = jj, min(jj+bl_segment_size-1, pdims%j_end) - do i = pdims%i_start, pdims%i_end - if (kstatus(i,j)) then - if ( (db_ga_dry_n(i,j,k) <= zero) .or. & - (k >= ntml(i,j)) ) then - kstatus(i,j)=.false. - kwb0(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 (kstatus(i,j)) then + if ( (db_ga_dry_n(i,j,k) <= zero) .or. & + (k >= ntml(i,j)) ) then + kstatus(i,j)=.false. + kwb0(i,j)=k end if - end do + end if + !end do end do end do end do @@ -1353,8 +1362,10 @@ subroutine excf_nl_9c ( & ! ---------------------------------------------------------------------- !cdir collapse !$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 +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 ( test_well_mixed(i,j) ) then dzrad(i,j) = 100.0_r_bl @@ -1369,7 +1380,7 @@ subroutine excf_nl_9c ( & ntop(i,j) = ntop(i,j) - 1 if ( ntop(i,j) > 2 ) then inner_loop1: do while ( z_tq(i,j,ntop(i,j)) > z_inv(i,j)-dzrad(i,j)& - .and. z_tq(i,j,ntop(i,j)-1) > z_inv(i,j)-z_cbase(i,j)) + .and. z_tq(i,j,ntop(i,j)-1) > z_inv(i,j)-z_cbase(i,j)) ntop(i,j) = ntop(i,j) - 1 if (ntop(i,j) == 2 ) then exit inner_loop1 @@ -1391,14 +1402,14 @@ subroutine excf_nl_9c ( & k_inv = ntop(i,j) if ( z_inv(i,j) < z_tq(i,j,k_inv+1) ) then interp = ( z_inv(i,j) - z_uv(i,j,k_inv+1) ) & - / ( z_tq(i,j,k_inv+1) - z_uv(i,j,k_inv+1) ) + / ( z_tq(i,j,k_inv+1) - z_uv(i,j,k_inv+1) ) z_pr = (one-interp) * z_tq(i,j,k_inv-1) & - + interp * z_uv(i,j,k_inv) + + interp * z_uv(i,j,k_inv) else interp = ( z_inv(i,j) - z_tq(i,j,k_inv+1) ) & - / ( z_uv(i,j,k_inv+2) - z_tq(i,j,k_inv+1) ) + / ( z_uv(i,j,k_inv+2) - z_tq(i,j,k_inv+1) ) z_pr = (one-interp) * z_uv(i,j,k_inv) & - + interp * z_tq(i,j,k_inv) + + interp * z_tq(i,j,k_inv) end if ! Go down to 100m below z_inv (or cloud-base) if that is lower, ! but don't allow lower than theta-level 1 @@ -1413,7 +1424,7 @@ subroutine excf_nl_9c ( & ! Fraction of rho-level ntop below the dzrad layer interp = ( z_pr - z_tq(i,j,ntop(i,j)-1) ) & - / ( z_tq(i,j,ntop(i,j)) - z_tq(i,j,ntop(i,j)-1) ) + / ( z_tq(i,j,ntop(i,j)) - z_tq(i,j,ntop(i,j)-1) ) ! Add on SL and qw differences between the discontinuous inversion base ! and the base of the dzrad layer dsl(i,j) = dsl(i,j) + sl(i,j,k_inv) - (one-interp)*sl(i,j,ntop(i,j)-1) & @@ -1445,7 +1456,8 @@ subroutine excf_nl_9c ( & end if end do ! I -end do ! J +end do ! II +!end do ! J !$OMP end do ! For WB diagnostics, convert integrated WB to uniform profile @@ -1453,19 +1465,19 @@ subroutine excf_nl_9c ( & !$OMP do SCHEDULE(STATIC) do k=2, bl_levels - do j=pdims%j_start, pdims%j_end - do i=pdims%i_start, pdims%i_end - if ( test_well_mixed(i,j) .and. k >= ntop(i,j)+1 ) then - if ( k <= ntdsc(i,j)+1 ) then - wbend(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) - wbmix(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) - else if ( k <= ntml(i,j)+1 ) then - wbend(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) - wbmix(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) - end if + !do j=pdims%j_start, pdims%j_end + do i=pdims%i_start, pdims%i_end + if ( test_well_mixed(i,j) .and. k >= ntop(i,j)+1 ) then + if ( k <= ntdsc(i,j)+1 ) then + wbend(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) + wbmix(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) + else if ( k <= ntml(i,j)+1 ) then + wbend(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) + wbmix(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) end if - end do ! i - end do ! j + end if + end do ! i + !end do ! j end do ! k !$OMP end do @@ -1480,8 +1492,10 @@ subroutine excf_nl_9c ( & !cdir collapse !$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 +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 ( test_well_mixed(i,j) ) then if ( kwb0(i,j) == ntml(i,j) ) then @@ -1509,19 +1523,20 @@ subroutine excf_nl_9c ( & wbp_int(i,j) = wbp_int(i,j) + wb_surf_int(i,j) ! must be >0 end do ! I -end do ! J +end do ! I +!end do ! J !$OMP end do if (model_type == mt_single_column) then !$OMP do SCHEDULE(STATIC) ! Save surface and bl-top layer integral for diagnostics - do j=pdims%j_start, pdims%j_end - do i=pdims%i_start, pdims%i_end - wbmix(i,j,ksurf(i,j)) = wb_surf_int(i,j) - wbend(i,j,ksurf(i,j)) = wb_surf_int(i,j) - end do ! i - end do ! j + !do j=pdims%j_start, pdims%j_end + do i=pdims%i_start, pdims%i_end + wbmix(i,j,ksurf(i,j)) = wb_surf_int(i,j) + wbend(i,j,ksurf(i,j)) = wb_surf_int(i,j) + end do ! i + !end do ! j !$OMP end do end if ! model_type @@ -1531,135 +1546,135 @@ subroutine excf_nl_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 = 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 ( test_well_mixed(i,j) ) then - ! ---------------------------------------------- - ! worth testing layer as well-mixed to cloud-top - ! ---------------------------------------------- - zb_ktop = 0.1_r_bl*z_inv(i,j) - zinv_pr(i,j) = z_inv(i,j) - zb_ktop - ! DB(K)is the K to K-1 difference and already - ! integrated up to KSURF, so start this loop at KSURF+1 - if ( (k >= ksurf(i,j)+1) .and. (k <= ntop(i,j)) ) then - - khtop(i,j) = zero - khsurf(i,j)= zero - f2 = zero - fsc = zero - - z_pr = z_uv(i,j,k) - zb_ktop - if (z_pr > zero .and. z_pr < zinv_pr(i,j)) then - z_ratio = z_pr/zinv_pr(i,j) + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) + !do i = pdims%i_start, pdims%i_end - if (flux_grad == LockWhelan2006) then - khtop(i,j) = 3.6_r_bl * vkman * rho_mix(i,j,k) *v_ktop(i,j) & - * zinv_pr(i,j) * (z_ratio**3) & - * ( (one-z_ratio)*(one-z_ratio) ) - f2 = rho_mix(i,j,k) * one_half * z_ratio & - * 2.0_r_bl**( z_ratio**4 ) - if ( v_ksum(i,j) > zero ) then - fsc = rho_mix(i,j,k) * 3.5_r_bl*(v_ktop(i,j)/v_ksum(i,j)) & - * (z_ratio**3) * (one-z_ratio) - end if - else - ! max to avoid rounding errors giving small negative numbers - khtop(i,j) = g1 * vkman * rho_mix(i,j,k) * v_ktop(i,j) & - * ((max(one - z_ratio,zero))**0.8_r_bl) & - * z_pr * z_ratio - end if - end if - - z_pr = z_uv(i,j,k) - if ( z_pr < z_inv(i,j)) then - !-------------------------------- - ! include surface-driven profile - !-------------------------------- - khsurf(i,j) = vkman * rho_mix(i,j,k) * & - w_h_top(i,j)*z_pr*( one - z_pr/z_inv(i,j) ) & - *( one - z_pr/z_inv(i,j) ) - end if + if ( test_well_mixed(i,j) ) then + ! ---------------------------------------------- + ! worth testing layer as well-mixed to cloud-top + ! ---------------------------------------------- + zb_ktop = 0.1_r_bl*z_inv(i,j) + zinv_pr(i,j) = z_inv(i,j) - zb_ktop + ! DB(K)is the K to K-1 difference and already + ! integrated up to KSURF, so start this loop at KSURF+1 + if ( (k >= ksurf(i,j)+1) .and. (k <= ntop(i,j)) ) then + + khtop(i,j) = zero + khsurf(i,j)= zero + f2 = zero + fsc = zero + + z_pr = z_uv(i,j,k) - zb_ktop + if (z_pr > zero .and. z_pr < zinv_pr(i,j)) then + z_ratio = z_pr/zinv_pr(i,j) if (flux_grad == LockWhelan2006) then - wslng = (f2+fsc)*tothf_zi(i,j) - ft_nt(i,j,k) - wqwng = (f2+fsc)*totqf_zi(i,j) - fq_nt(i,j,k) - end if - - if ( z_tq(i,j,k) <= z_cbase(i,j) ) then - ! Completely below cloud-base so use cloud-free formula - wb_scld = khsurf(i,j) * db_ga_dry_n(i,j,k) + & - khtop(i,j) * db_noga_dry(i,j,k) - if (flux_grad == LockWhelan2006) then - wb_scld = wb_scld + ( g/rdz(i,j,k) ) * & - ( btm(i,j,k-1)*wslng + bqm(i,j,k-1)*wqwng ) - end if - wb_cld = zero - else if (z_tq(i,j,k-1) >= z_cbase(i,j)) then - ! Completely above cloud-base so use cloudy formula - wb_cld = ( khsurf(i,j) * db_ga_cld_n(i,j,k) + & - khtop(i,j) * db_noga_cld(i,j,k) ) - if (flux_grad == LockWhelan2006) then - wb_cld = wb_cld + ( g/rdz(i,j,k) ) * ( & - ( btm(i,j,k-1)*(one-cf_ml(i,j)) + & - btm_cld(i,j,k-1)*cf_ml(i,j) )*wslng + & - ( bqm(i,j,k-1)*(one-cf_ml(i,j)) + & - bqm_cld(i,j,k-1)*cf_ml(i,j) )*wqwng ) + khtop(i,j) = 3.6_r_bl * vkman * rho_mix(i,j,k) *v_ktop(i,j) & + * zinv_pr(i,j) * (z_ratio**3) & + * ( (one-z_ratio)*(one-z_ratio) ) + f2 = rho_mix(i,j,k) * one_half * z_ratio & + * 2.0_r_bl**( z_ratio**4 ) + if ( v_ksum(i,j) > zero ) then + fsc = rho_mix(i,j,k) * 3.5_r_bl*(v_ktop(i,j)/v_ksum(i,j)) & + * (z_ratio**3) * (one-z_ratio) end if - wb_scld = zero else - ! cloud-base within this integration range - ! so treat cloud and sub-cloud layer wb separately - wb_scld = khsurf(i,j) * db_ga_dry_n(i,j,k) + & - khtop(i,j) * db_noga_dry(i,j,k) - wb_cld = ( khsurf(i,j) * db_ga_cld_n(i,j,k) + & - khtop(i,j) * db_noga_cld(i,j,k) ) - if (flux_grad == LockWhelan2006) then - wb_scld = wb_scld + ( g/rdz(i,j,k) ) * & - ( btm(i,j,k-1)*wslng + bqm(i,j,k-1)*wqwng ) - wb_cld = wb_cld + ( g/rdz(i,j,k) ) * ( & - ( btm(i,j,k-1)*(one-cf_ml(i,j)) + & - btm_cld(i,j,k-1)*cf_ml(i,j) )*wslng + & - ( bqm(i,j,k-1)*(one-cf_ml(i,j)) + & - bqm_cld(i,j,k-1)*cf_ml(i,j) )*wqwng ) - end if - cld_frac = (z_tq(i,j,k)-z_cbase(i,j)) & - /(z_tq(i,j,k)-z_tq(i,j,k-1)) - wb_cld = cld_frac * wb_cld - wb_scld = (one-cld_frac) * wb_scld + ! max to avoid rounding errors giving small negative numbers + khtop(i,j) = g1 * vkman * rho_mix(i,j,k) * v_ktop(i,j) & + * ((max(one - z_ratio,zero))**0.8_r_bl) & + * z_pr * z_ratio end if + end if - if ( dzrad_disc_opt == dzrad_1p5dz .and. k == ntop(i,j) ) then - ! At top of Sc layer, only include the part of the integral - ! that is below the radiatively-cooled cloud-top layer. - interp = ( z_inv(i,j) - dzrad(i,j) - z_tq(i,j,k-1) ) & - / ( z_tq(i,j,k) - z_tq(i,j,k-1) ) - wb_cld = wb_cld * interp - wb_scld = wb_scld * interp - end if + z_pr = z_uv(i,j,k) + if ( z_pr < z_inv(i,j)) then + !-------------------------------- + ! include surface-driven profile + !-------------------------------- + khsurf(i,j) = vkman * rho_mix(i,j,k) * & + w_h_top(i,j)*z_pr*( one - z_pr/z_inv(i,j) ) & + *( one - z_pr/z_inv(i,j) ) + end if - if (wb_cld >= zero) then - wbp_int(i,j) = wbp_int(i,j) + wb_cld - else - wbn_int(i,j) = wbn_int(i,j) - wb_cld + if (flux_grad == LockWhelan2006) then + wslng = (f2+fsc)*tothf_zi(i,j) - ft_nt(i,j,k) + wqwng = (f2+fsc)*totqf_zi(i,j) - fq_nt(i,j,k) + end if + + if ( z_tq(i,j,k) <= z_cbase(i,j) ) then + ! Completely below cloud-base so use cloud-free formula + wb_scld = khsurf(i,j) * db_ga_dry_n(i,j,k) + & + khtop(i,j) * db_noga_dry(i,j,k) + if (flux_grad == LockWhelan2006) then + wb_scld = wb_scld + ( g/rdz(i,j,k) ) * & + ( btm(i,j,k-1)*wslng + bqm(i,j,k-1)*wqwng ) end if - if (wb_scld >= zero) then - wbp_int(i,j) = wbp_int(i,j) + wb_scld - else - wbn_int(i,j) = wbn_int(i,j) - wb_scld + wb_cld = zero + else if (z_tq(i,j,k-1) >= z_cbase(i,j)) then + ! Completely above cloud-base so use cloudy formula + wb_cld = ( khsurf(i,j) * db_ga_cld_n(i,j,k) + & + khtop(i,j) * db_noga_cld(i,j,k) ) + if (flux_grad == LockWhelan2006) then + wb_cld = wb_cld + ( g/rdz(i,j,k) ) * ( & + ( btm(i,j,k-1)*(one-cf_ml(i,j)) + & + btm_cld(i,j,k-1)*cf_ml(i,j) )*wslng + & + ( bqm(i,j,k-1)*(one-cf_ml(i,j)) + & + bqm_cld(i,j,k-1)*cf_ml(i,j) )*wqwng ) end if - - if (model_type == mt_single_column) then - wbmix(i,j,k) = wb_cld+wb_scld - wbend(i,j,k) = wb_cld+wb_scld + wb_scld = zero + else + ! cloud-base within this integration range + ! so treat cloud and sub-cloud layer wb separately + wb_scld = khsurf(i,j) * db_ga_dry_n(i,j,k) + & + khtop(i,j) * db_noga_dry(i,j,k) + wb_cld = ( khsurf(i,j) * db_ga_cld_n(i,j,k) + & + khtop(i,j) * db_noga_cld(i,j,k) ) + if (flux_grad == LockWhelan2006) then + wb_scld = wb_scld + ( g/rdz(i,j,k) ) * & + ( btm(i,j,k-1)*wslng + bqm(i,j,k-1)*wqwng ) + wb_cld = wb_cld + ( g/rdz(i,j,k) ) * ( & + ( btm(i,j,k-1)*(one-cf_ml(i,j)) + & + btm_cld(i,j,k-1)*cf_ml(i,j) )*wslng + & + ( bqm(i,j,k-1)*(one-cf_ml(i,j)) + & + bqm_cld(i,j,k-1)*cf_ml(i,j) )*wqwng ) end if + cld_frac = (z_tq(i,j,k)-z_cbase(i,j)) & + /(z_tq(i,j,k)-z_tq(i,j,k-1)) + wb_cld = cld_frac * wb_cld + wb_scld = (one-cld_frac) * wb_scld + end if - end if ! K + if ( dzrad_disc_opt == dzrad_1p5dz .and. k == ntop(i,j) ) then + ! At top of Sc layer, only include the part of the integral + ! that is below the radiatively-cooled cloud-top layer. + interp = ( z_inv(i,j) - dzrad(i,j) - z_tq(i,j,k-1) ) & + / ( z_tq(i,j,k) - z_tq(i,j,k-1) ) + wb_cld = wb_cld * interp + wb_scld = wb_scld * interp + end if - end if ! TEST_WELL_MIXED - end do ! I + if (wb_cld >= zero) then + wbp_int(i,j) = wbp_int(i,j) + wb_cld + else + wbn_int(i,j) = wbn_int(i,j) - wb_cld + end if + if (wb_scld >= zero) then + wbp_int(i,j) = wbp_int(i,j) + wb_scld + else + wbn_int(i,j) = wbn_int(i,j) - wb_scld + end if + + if (model_type == mt_single_column) then + wbmix(i,j,k) = wb_cld+wb_scld + wbend(i,j,k) = wb_cld+wb_scld + end if + + end if ! K + + end if ! TEST_WELL_MIXED + !end do ! I end do ! J end do ! K end do @@ -1671,8 +1686,10 @@ subroutine excf_nl_9c ( & !cdir collapse !$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 +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 (l_wtrac) iset_wtrac(i,j) = 0 if ( test_well_mixed(i,j) ) then @@ -1692,14 +1709,14 @@ subroutine excf_nl_9c ( & if ( db_top(i,j) > zero .and. db_dsct(i,j) > 0.01_r_bl ) then ! can't calc Zil. term rhokh_surf_ent(i,j) = rhokh_surf_ent(i,j) * & - ( rho_mix(i,j,ntdsc(i,j)+1) * db_top(i,j) * & - rdz(i,j,ntml(i,j)+1) ) / & - ( rho_mix(i,j,ntml(i,j)+1) * db_dsct(i,j) * & + ( rho_mix(i,j,ntdsc(i,j)+1) * db_top(i,j) * & + rdz(i,j,ntml(i,j)+1) ) / & + ( rho_mix(i,j,ntml(i,j)+1) * db_dsct(i,j) * & rdz(i,j,ntdsc(i,j)+1) ) rhokm_surf_ent(i,j) = rhokm_surf_ent(i,j) * & - ( rho_wet_tq(i,j,ntdsc(i,j)) * db_top(i,j) * & - rdz(i,j,ntml(i,j)+1) ) / & - ( rho_wet_tq(i,j,ntml(i,j)) * db_dsct(i,j) * & + ( rho_wet_tq(i,j,ntdsc(i,j)) * db_top(i,j) * & + rdz(i,j,ntml(i,j)+1) ) / & + ( rho_wet_tq(i,j,ntml(i,j)) * db_dsct(i,j) * & rdz(i,j,ntdsc(i,j)+1) ) if (.not. l_use_var_fixes) & rhokm(i,j,ntdsc(i,j)+1) = rhokm_surf_ent(i,j) @@ -1755,7 +1772,7 @@ subroutine excf_nl_9c ( & ! turn on ktop_iterate to find the new DSC's base ktop_iterate(i,j) = .true. if ( (.not. cumulus(i,j)) .or. kprof_cu == buoy_integ .or. & - kprof_cu == buoy_integ_low ) then + kprof_cu == buoy_integ_low ) then ! Switch to using b-flux integration to find new lower SML-top, ! but not in Cu layers unless using appropriate kprof_cu option. ksurf_iterate(i,j) = .true. @@ -1774,11 +1791,11 @@ subroutine excf_nl_9c ( & ! Extent of mixing must be reduced !--------------------------------- if ( .not. cumulus(i,j) .or. kprof_cu == buoy_integ .or. & - kprof_cu == buoy_integ_low) & - ksurf_iterate(i,j) = .true. + kprof_cu == buoy_integ_low) & + ksurf_iterate(i,j) = .true. if ( cumulus(i,j) .and. ( kprof_cu == buoy_integ .or. & kprof_cu == buoy_integ_low) ) & - dec_thres(i,j) = dec_thres_cu + dec_thres(i,j) = dec_thres_cu ktop_iterate(i,j) = .true. end if ! ( l_use_sml_dsc_fixes ) @@ -1821,7 +1838,7 @@ subroutine excf_nl_9c ( & rhokm_surf_ent(i,j) = zero if (.not. l_use_var_fixes) then rhokm_top(i,j,ntml(i,j)+1) = rhokm_top(i,j,ntml(i,j)+1) & - + rhokm(i,j,ntml(i,j)+1) + + rhokm(i,j,ntml(i,j)+1) rhokm(i,j,ntml(i,j)+1) = zero end if end if ! test on l_use_var_fixes and v_top <= zero @@ -1830,23 +1847,24 @@ subroutine excf_nl_9c ( & end if ! testing for well-mixed layer (TEST_WELL_MIXED) end do ! I -end do ! J +end do ! I +!end do ! J !$OMP end do ! Set water tracer fields in the same manner as water in last loop 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 (iset_wtrac(i,j) == 1) then - fq_nt_zh_wtrac(i,j,i_wt) = fq_nt_zhsc_wtrac(i,j,i_wt) - fq_nt_zhsc_wtrac(i,j,i_wt) = zero - else if (iset_wtrac(i,j) == 2) then - fq_nt_zhsc_wtrac(i,j,i_wt) = fq_nt_zh_wtrac(i,j,i_wt) - end if - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + if (iset_wtrac(i,j) == 1) then + fq_nt_zh_wtrac(i,j,i_wt) = fq_nt_zhsc_wtrac(i,j,i_wt) + fq_nt_zhsc_wtrac(i,j,i_wt) = zero + else if (iset_wtrac(i,j) == 2) then + fq_nt_zhsc_wtrac(i,j,i_wt) = fq_nt_zh_wtrac(i,j,i_wt) + end if end do + !end do !$OMP end do end do @@ -1865,57 +1883,59 @@ subroutine excf_nl_9c ( & !cdir collapse !$OMP do SCHEDULE(DYNAMIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - - if ( ksurf_iterate(i,j) ) then - !----------------------------------------------------------- - ! Mixing must extend just above surface layer - ! (not clear precisely how to define this here: for now use - ! KSURF calculated from ZH) - !----------------------------------------------------------- - z_bot_lim(i,j)=z_uv(i,j,ksurf(i,j)+1) & - + 0.1_r_bl * (z_uv(i,j,ksurf(i,j)+2)-z_uv(i,j,ksurf(i,j)+1)) - if ( l_converge_ga ) then - ! Allow K-surf to go up to the Sc-top (increasing zsml_top within - ! the cloud-top layer can still make the buoyancy fluxes below more - ! negative because this reduces the gradient adjustment) - z_top_lim(i,j)=max( z_top_lim(i,j), & - max( z_bot_lim(i,j), zhsc(i,j) ) ) - ! ksurf_iterate currently only true if test_well_mixed was true and - ! decoupling was diagnosed - if (cumulus(i,j)) then - ! ignore cloudy component of wb for cumulus sml iteration as the - ! cloudy flux is carried by the convection scheme - z_cbase(i,j) = z_top_lim(i,j) - else - ! sml has decoupled from the cloud layer above. For consistency with - ! the test_well_mixed calculation, continue to use saturated wb above - ! the now decoupled layer's cloud base. - z_cbase(i,j) = zhsc(i,j) - zc_dsc(i,j) - end if - else ! ( .not. l_converge_ga ) - ! limit K-surf to below cloud-top radiatively cooled layer - ! or original zh (held in z_top_lim) if no top-driven turbulence - z_top_lim(i,j)=max( z_top_lim(i,j), & - max( z_bot_lim(i,j), zhsc(i,j) - dzrad(i,j) ) ) - ! Always use saturated wb above decoupled-layer cloud-base +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) +!do j = pdims%j_start, pdims%j_end +!do i = pdims%i_start, pdims%i_end + + if ( ksurf_iterate(i,j) ) then + !----------------------------------------------------------- + ! Mixing must extend just above surface layer + ! (not clear precisely how to define this here: for now use + ! KSURF calculated from ZH) + !----------------------------------------------------------- + z_bot_lim(i,j)=z_uv(i,j,ksurf(i,j)+1) & + + 0.1_r_bl * (z_uv(i,j,ksurf(i,j)+2)-z_uv(i,j,ksurf(i,j)+1)) + if ( l_converge_ga ) then + ! Allow K-surf to go up to the Sc-top (increasing zsml_top within + ! the cloud-top layer can still make the buoyancy fluxes below more + ! negative because this reduces the gradient adjustment) + z_top_lim(i,j)=max( z_top_lim(i,j), & + max( z_bot_lim(i,j), zhsc(i,j) ) ) + ! ksurf_iterate currently only true if test_well_mixed was true and + ! decoupling was diagnosed + if (cumulus(i,j)) then + ! ignore cloudy component of wb for cumulus sml iteration as the + ! cloudy flux is carried by the convection scheme + z_cbase(i,j) = z_top_lim(i,j) + else + ! sml has decoupled from the cloud layer above. For consistency with + ! the test_well_mixed calculation, continue to use saturated wb above + ! the now decoupled layer's cloud base. z_cbase(i,j) = zhsc(i,j) - zc_dsc(i,j) - end if ! ( l_converge_ga ) - !----------------------------------------------------- - ! Initial increment to ZSML_TOP found by dividing - ! up depth of layer within which it is allowed: - ! Start with ZSML_TOP at lower limit and work upwards - !----------------------------------------------------- - z_inc(i,j)=(z_top_lim(i,j)-z_bot_lim(i,j)) / real(n_steps, r_bl) - zsml_top(i,j) = z_bot_lim(i,j) - - wb_ratio(i,j) = dec_thres(i,j) - one ! to be < DEC_THRES - - end if ! KSURF_ITERATE - + end if + else ! ( .not. l_converge_ga ) + ! limit K-surf to below cloud-top radiatively cooled layer + ! or original zh (held in z_top_lim) if no top-driven turbulence + z_top_lim(i,j)=max( z_top_lim(i,j), & + max( z_bot_lim(i,j), zhsc(i,j) - dzrad(i,j) ) ) + ! Always use saturated wb above decoupled-layer cloud-base + z_cbase(i,j) = zhsc(i,j) - zc_dsc(i,j) + end if ! ( l_converge_ga ) + !----------------------------------------------------- + ! Initial increment to ZSML_TOP found by dividing + ! up depth of layer within which it is allowed: + ! Start with ZSML_TOP at lower limit and work upwards + !----------------------------------------------------- + z_inc(i,j)=(z_top_lim(i,j)-z_bot_lim(i,j)) / real(n_steps, r_bl) + zsml_top(i,j) = z_bot_lim(i,j) + + wb_ratio(i,j) = dec_thres(i,j) - one ! to be < DEC_THRES + + end if ! KSURF_ITERATE end do end do +!end do !$OMP end do ! ---------------------------------------------------------------------- @@ -1927,53 +1947,56 @@ subroutine excf_nl_9c ( & if (kprof_cu == buoy_integ_low) lcl_fac = one_half !$OMP do SCHEDULE(DYNAMIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( cumulus(i,j) .and. .not. dsc(i,j) .and. & - .not. ksurf_iterate(i,j)) then - ! pure cumulus layer and necessary parameters not already set up - ksurf_iterate(i,j) = .true. - dec_thres(i,j) = dec_thres_cu ! use cu threshold - ! limit top of K_surf to be between... - z_bot_lim(i,j) = lcl_fac*z_lcl(i,j) ! ...a fraction of the LCL and... - z_top_lim(i,j) = z_lcl(i,j) + 1000.0_r_bl ! ...1km above the LCL - ntop(i,j)=ntml(i,j) - do while ( z_uv(i,j,ntop(i,j)+1) <= z_top_lim(i,j) .and. & - ntop(i,j)+1 < bl_levels-2 ) - ntop(i,j) = ntop(i,j) + 1 ! z_uv(ntop+1) > z_top_lim - end do - z_top_lim(i,j)=z_uv(i,j,ntop(i,j)) ! highest z_uv below z_top_lim - z_cbase(i,j) = z_top_lim(i,j) - !----------------------------------------------------- - ! Initial increment to ZSML_TOP found by dividing - ! up depth of layer within which it is allowed: - ! Start with ZSML_TOP at lower limit and work upwards - !----------------------------------------------------- - z_inc(i,j)=(z_top_lim(i,j)-z_bot_lim(i,j)) & - / real(n_steps, r_bl) - zsml_top(i,j) = z_bot_lim(i,j) - wb_ratio(i,j) = dec_thres(i,j) - one ! to be < DEC_THRES - end if - end do ! I - end do ! J + !do j = pdims%j_start, pdims%j_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 ( cumulus(i,j) .and. .not. dsc(i,j) .and. & + .not. ksurf_iterate(i,j)) then + ! pure cumulus layer and necessary parameters not already set up + ksurf_iterate(i,j) = .true. + dec_thres(i,j) = dec_thres_cu ! use cu threshold + ! limit top of K_surf to be between... + z_bot_lim(i,j) = lcl_fac*z_lcl(i,j) ! ...a fraction of the LCL and... + z_top_lim(i,j) = z_lcl(i,j) + 1000.0_r_bl ! ...1km above the LCL + ntop(i,j)=ntml(i,j) + do while ( z_uv(i,j,ntop(i,j)+1) <= z_top_lim(i,j) .and. & + ntop(i,j)+1 < bl_levels-2 ) + ntop(i,j) = ntop(i,j) + 1 ! z_uv(ntop+1) > z_top_lim + end do + z_top_lim(i,j)=z_uv(i,j,ntop(i,j)) ! highest z_uv below z_top_lim + z_cbase(i,j) = z_top_lim(i,j) + !----------------------------------------------------- + ! Initial increment to ZSML_TOP found by dividing + ! up depth of layer within which it is allowed: + ! Start with ZSML_TOP at lower limit and work upwards + !----------------------------------------------------- + z_inc(i,j)=(z_top_lim(i,j)-z_bot_lim(i,j)) & + / real(n_steps, r_bl) + zsml_top(i,j) = z_bot_lim(i,j) + wb_ratio(i,j) = dec_thres(i,j) - one ! to be < DEC_THRES + end if + end do ! I +end do ! II + !end do ! J !$OMP end do end if ! test in kprof_cu ! ---------------------------------------------------------------------- !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - l = i - pdims%i_start + 1 + pdims%i_len * (j - pdims%j_start) - ind_todo(l) = l - up(l) = 1 - if (ksurf_iterate(i,j)) then - to_do(l) = .true. - else - to_do(l) = .false. - end if - end do +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + l = i - pdims%i_start + 1 + pdims%i_len * (j - pdims%j_start) + ind_todo(l) = l + up(l) = 1 + if (ksurf_iterate(i,j)) then + to_do(l) = .true. + else + to_do(l) = .false. + end if end do +!end do !$OMP end do !$OMP MASTER @@ -2180,56 +2203,56 @@ subroutine excf_nl_9c ( & end do ! n_sweep !$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 - ntml_new(i,j) = 2 - status_ntml(i,j)=.true. - 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_new(i,j) = 2 + status_ntml(i,j)=.true. + !end do end do do k = 2, bl_levels-2 !cdir collapse - do j = jj, min(jj+bl_segment_size-1, pdims%j_end) - do i = pdims%i_start, pdims%i_end - if ( ksurf_iterate(i,j) .and. status_ntml(i,j) ) then - ! ------------- - ! find new NTML - ! ------------- - if (z_uv(i,j,k+1) < zsml_top(i,j)) then - ntml_new(i,j) = k+1 + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) + !do i = pdims%i_start, pdims%i_end + if ( ksurf_iterate(i,j) .and. status_ntml(i,j) ) then + ! ------------- + ! find new NTML + ! ------------- + if (z_uv(i,j,k+1) < zsml_top(i,j)) then + ntml_new(i,j) = k+1 + else + status_ntml(i,j)=.false. + end if + ! -------------------------------------------------------- + ! Rounding error previously found to give + ! ZSML_TOP > Z_TOP_LIM = ZHSC + ! Test on ZSML_TOP hitting thresholds consequently changed + ! but also include the following failsafe tests here. + ! -------------------------------------------------------- + if (dsc(i,j)) then + ntml(i,j) = min( ntdsc(i,j), ntml_new(i,j)-1 ) + zh(i,j) = min( zhsc(i,j), zsml_top(i,j) ) + else if (.not. status_ntml(i,j)) then + if (kprof_cu == buoy_integ_low) then + ! just use zh from buoyancy flux integration + ntml(i,j) = ntml_new(i,j)-1 + zh(i,j) = zsml_top(i,j) else - status_ntml(i,j)=.false. - end if - ! -------------------------------------------------------- - ! Rounding error previously found to give - ! ZSML_TOP > Z_TOP_LIM = ZHSC - ! Test on ZSML_TOP hitting thresholds consequently changed - ! but also include the following failsafe tests here. - ! -------------------------------------------------------- - if (dsc(i,j)) then - ntml(i,j) = min( ntdsc(i,j), ntml_new(i,j)-1 ) - zh(i,j) = min( zhsc(i,j), zsml_top(i,j) ) - else if (.not. status_ntml(i,j)) then - if (kprof_cu == buoy_integ_low) then - ! just use zh from buoyancy flux integration - ntml(i,j) = ntml_new(i,j)-1 - zh(i,j) = zsml_top(i,j) + ! for kprof_cu eq buoy_integ use max of ksurf top calculated + ! from buoyancy flux integration and original LCL based values + ntml(i,j) = max( ntml(i,j), ntml_new(i,j)-1 ) + if (l_use_var_fixes) then + zh(i,j) = max( zh(i,j), zsml_top(i,j) ) else - ! for kprof_cu eq buoy_integ use max of ksurf top calculated - ! from buoyancy flux integration and original LCL based values - ntml(i,j) = max( ntml(i,j), ntml_new(i,j)-1 ) - if (l_use_var_fixes) then - zh(i,j) = max( zh(i,j), zsml_top(i,j) ) - else - zh(i,j) = zsml_top(i,j) - end if + zh(i,j) = zsml_top(i,j) end if end if + end if - end if ! KSURF_ITERATE true + end if ! KSURF_ITERATE true - end do + !end do end do end do end do @@ -2241,8 +2264,10 @@ subroutine excf_nl_9c ( & ! ---------------------------------------------------------------------- !cdir collapse !$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 +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 ( ktop_iterate(i,j) ) then @@ -2263,7 +2288,7 @@ subroutine excf_nl_9c ( & !..Divide up depth of layer within which ZDSC_BASE is allowed z_inc(i,j) = (z_top_lim(i,j)-z_bot_lim(i,j)) / real(n_steps, r_bl) zdsc_base(i,j) = z_bot_lim(i,j) - ! will start at Z_BOT_LIM+Z_INC + ! will start at Z_BOT_LIM+Z_INC if (l_reset_dec_thres) then dec_thres(i,j) = dec_thres_cloud ! reset after potential use @@ -2275,8 +2300,13 @@ subroutine excf_nl_9c ( & end if ! KTOP_ITERATE end do +end do +!$OMP end do - do i = pdims%i_start, pdims%i_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) +! do i = pdims%i_start, pdims%i_end if ( ktop_iterate(i,j) .and. wb_dzrad_int(i,j) < zero ) then !------------------------------------------------------------- ! Estimation of wb integral over radiatively cooled cloud-top @@ -2305,7 +2335,7 @@ subroutine excf_nl_9c ( & ntop(i,j) = ntop(i,j) - 1 if ( ntop(i,j) > 1 ) then inner_loop2: do while ( z_tq(i,j,ntop(i,j)) > z_inv(i,j)-dzrad(i,j)& - .and. z_tq(i,j,ntop(i,j)-1) > z_inv(i,j)-z_cbase(i,j)) + .and. z_tq(i,j,ntop(i,j)-1) > z_inv(i,j)-z_cbase(i,j)) ntop(i,j) = ntop(i,j) - 1 if (ntop(i,j) == 1 ) then exit inner_loop2 @@ -2327,14 +2357,14 @@ subroutine excf_nl_9c ( & k_inv = ntop(i,j) if ( z_inv(i,j) < z_tq(i,j,k_inv+1) ) then interp = ( z_inv(i,j) - z_uv(i,j,k_inv+1) ) & - / ( z_tq(i,j,k_inv+1) - z_uv(i,j,k_inv+1) ) + / ( z_tq(i,j,k_inv+1) - z_uv(i,j,k_inv+1) ) z_pr = (one-interp) * z_tq(i,j,k_inv-1) & - + interp * z_uv(i,j,k_inv) + + interp * z_uv(i,j,k_inv) else interp = ( z_inv(i,j) - z_tq(i,j,k_inv+1) ) & - / ( z_uv(i,j,k_inv+2) - z_tq(i,j,k_inv+1) ) + / ( z_uv(i,j,k_inv+2) - z_tq(i,j,k_inv+1) ) z_pr = (one-interp) * z_uv(i,j,k_inv) & - + interp * z_tq(i,j,k_inv) + + interp * z_tq(i,j,k_inv) end if ! Go down to 100m below z_inv (or cloud-base) if that is lower, ! but don't allow lower than theta-level 1 @@ -2349,7 +2379,7 @@ subroutine excf_nl_9c ( & ! Fraction of rho-level ntop below the dzrad layer interp = ( z_pr - z_tq(i,j,ntop(i,j)-1) ) & - / ( z_tq(i,j,ntop(i,j)) - z_tq(i,j,ntop(i,j)-1) ) + / ( z_tq(i,j,ntop(i,j)) - z_tq(i,j,ntop(i,j)-1) ) ! Add on SL and qw differences between the discontinuous inversion base ! and the base of the dzrad layer dsl(i,j) = dsl(i,j) + sl(i,j,k_inv) - (one-interp)*sl(i,j,ntop(i,j)-1) & @@ -2381,24 +2411,24 @@ subroutine excf_nl_9c ( & end if wb_dzrad_int(i,j) = max( rbl_eps, wb_dzrad_int(i,j) ) - end do ! I -end do ! J +end do ! II +!end do ! J !$OMP end do !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - l = i - pdims%i_start + 1 + pdims%i_len * (j - pdims%j_start) - ind_todo(l) = l - up(l) = 1 - if (ktop_iterate(i,j)) then - to_do(l) = .true. - else - to_do(l) = .false. - end if - end do +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + l = i - pdims%i_start + 1 + pdims%i_len * (j - pdims%j_start) + ind_todo(l) = l + up(l) = 1 + if (ktop_iterate(i,j)) then + to_do(l) = .true. + else + to_do(l) = .false. + end if end do +!end do !$OMP end do !$OMP MASTER @@ -2679,60 +2709,62 @@ subroutine excf_nl_9c ( & !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels-1 - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ! Include sml - if ( .not. zdsc_base(i,j) < zsml_top(i,j) ) then - ! Note: wbend already includes the SML buoyancy flux - ! in the case where the SML and DSC mixing profiles overlap - ! (since the DSC b-flux integration sums both the SML and DSC - ! mixing contributions and stores this in wbend). - ! So only add the SML contribution back on if they do not overlap. - wbend(i,j,k) = wbend(i,j,k) + wbend_sml(i,j,k) - end if - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + ! Include sml + if ( .not. zdsc_base(i,j) < zsml_top(i,j) ) then + ! Note: wbend already includes the SML buoyancy flux + ! in the case where the SML and DSC mixing profiles overlap + ! (since the DSC b-flux integration sums both the SML and DSC + ! mixing contributions and stores this in wbend). + ! So only add the SML contribution back on if they do not overlap. + wbend(i,j,k) = wbend(i,j,k) + wbend_sml(i,j,k) + end if end do + !end do end do !$OMP end do -!$OMP do SCHEDULE(STATIC) + ! Note parallelised over j as k isn't independent - do j = pdims%j_start, pdims%j_end - do k = 1, bl_levels-1 - do i = pdims%i_start, pdims%i_end - ! convert to m2/s-3 - if ( k <= ntop(i,j) ) then - if ( k <= ksurf(i,j) ) then - gamma_wbs = ( (wbmix(i,j,ksurf(i,j))/z_tq(i,j,ksurf(i,j))) & - - bflux_surf(i,j) )*2.0_r_bl/z_tq(i,j,ksurf(i,j)) - wbmix(i,j,k) = bflux_surf(i,j) + gamma_wbs*z_uv(i,j,k) - - gamma_wbs = ( (wbend(i,j,ksurf(i,j))/z_tq(i,j,ksurf(i,j))) & - - bflux_surf(i,j) )*2.0_r_bl/z_tq(i,j,ksurf(i,j)) - wbend(i,j,k) = bflux_surf(i,j) + gamma_wbs*z_uv(i,j,k) - else - wbmix(i,j,k)=wbmix(i,j,k)*rdz(i,j,k) - wbend(i,j,k)=wbend(i,j,k)*rdz(i,j,k) - end if +!do j = pdims%j_start, pdims%j_end + do k = 1, bl_levels-1 + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + ! convert to m2/s-3 + if ( k <= ntop(i,j) ) then + if ( k <= ksurf(i,j) ) then + gamma_wbs = ( (wbmix(i,j,ksurf(i,j))/z_tq(i,j,ksurf(i,j))) & + - bflux_surf(i,j) )*2.0_r_bl/z_tq(i,j,ksurf(i,j)) + wbmix(i,j,k) = bflux_surf(i,j) + gamma_wbs*z_uv(i,j,k) + + gamma_wbs = ( (wbend(i,j,ksurf(i,j))/z_tq(i,j,ksurf(i,j))) & + - bflux_surf(i,j) )*2.0_r_bl/z_tq(i,j,ksurf(i,j)) + wbend(i,j,k) = bflux_surf(i,j) + gamma_wbs*z_uv(i,j,k) + else + wbmix(i,j,k)=wbmix(i,j,k)*rdz(i,j,k) + wbend(i,j,k)=wbend(i,j,k)*rdz(i,j,k) end if - end do + end if end do + !$OMP end do end do -!$OMP end do + !end do + !$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 - ! convert to m2/s-3 - if ( k >= ntop(i,j)+1 .and. k <= ntdsc(i,j)+1 ) then - wbend(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) - end if - if ( k >= ntop(i,j)+1 .and. k <= ntml(i,j)+1 ) then - wbmix(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) - end if - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + ! convert to m2/s-3 + if ( k >= ntop(i,j)+1 .and. k <= ntdsc(i,j)+1 ) then + wbend(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) + end if + if ( k >= ntop(i,j)+1 .and. k <= ntml(i,j)+1 ) then + wbmix(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) + end if end do + !end do end do !$OMP end do end if ! model_type @@ -2748,74 +2780,74 @@ subroutine excf_nl_9c ( & !cdir collapse !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( cumulus(i,j) .or. & - ( dsc(i,j) .and. zdsc_base(i,j) < zh(i,j) ) ) then - ! ignore SML `cloud-top' driven mixing - zsml_base(i,j) = zh(i,j) - v_top(i,j) = zero - else - zsml_base(i,j) = 0.1_r_bl*zh(i,j) - end if - end do ! loop over j - end do ! loop over I + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + if ( cumulus(i,j) .or. & + ( dsc(i,j) .and. zdsc_base(i,j) < zh(i,j) ) ) then + ! ignore SML `cloud-top' driven mixing + zsml_base(i,j) = zh(i,j) + v_top(i,j) = zero + else + zsml_base(i,j) = 0.1_r_bl*zh(i,j) + end if + end do ! loop over j + !end do ! loop over I !$OMP end do else ! entr_smooth_dec off !cdir collapse !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( cumulus(i,j) .or. coupled(i,j) ) then - ! ignore SML `cloud-top' driven mixing - zsml_base(i,j) = zh(i,j) - v_top(i,j) = zero - else - zsml_base(i,j) = 0.1_r_bl*zh(i,j) - end if - end do ! loop over j - end do ! loop over I + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + if ( cumulus(i,j) .or. coupled(i,j) ) then + ! ignore SML `cloud-top' driven mixing + zsml_base(i,j) = zh(i,j) + v_top(i,j) = zero + else + zsml_base(i,j) = 0.1_r_bl*zh(i,j) + end if + end do ! loop over j + !end do ! loop over I !$OMP end do end if ! test on entr_smooth_dec !cdir collapse !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - zsml_top(i,j) = zh(i,j) - if (bl_res_inv /= off .and. dzh(i,j) > zero) & - zsml_top(i,j) = zh(i,j)+dzh(i,j) - end do ! loop over j -end do ! loop over I +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + zsml_top(i,j) = zh(i,j) + if (bl_res_inv /= off .and. dzh(i,j) > zero) & + zsml_top(i,j) = zh(i,j)+dzh(i,j) +end do ! loop over j +!end do ! loop over I !$OMP end do if ( kprof_cu == klcl_entr ) then !cdir collapse !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( cumulus(i,j) ) then - zsml_top(i,j) = zh(i,j)+max_cu_depth - ! Keep top of Kprof for cu below DSC base: - if (dsc(i,j)) then - zsml_top(i,j) = min( zsml_top(i,j), zdsc_base(i,j) ) - end if - ! ...but at least 1.1*z_lcl: - zsml_top(i,j) = max( 1.1_r_bl*z_lcl(i,j), zsml_top(i,j) ) - ! ...but no higher than parcel top: - zsml_top(i,j) = min( zhpar(i,j), zsml_top(i,j) ) + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + if ( cumulus(i,j) ) then + zsml_top(i,j) = zh(i,j)+max_cu_depth + ! Keep top of Kprof for cu below DSC base: + if (dsc(i,j)) then + zsml_top(i,j) = min( zsml_top(i,j), zdsc_base(i,j) ) + end if + ! ...but at least 1.1*z_lcl: + zsml_top(i,j) = max( 1.1_r_bl*z_lcl(i,j), zsml_top(i,j) ) + ! ...but no higher than parcel top: + zsml_top(i,j) = min( zhpar(i,j), zsml_top(i,j) ) - cu_depth_scale(i,j) = (zsml_top(i,j)-zh(i,j))/3.0_r_bl + cu_depth_scale(i,j) = (zsml_top(i,j)-zh(i,j))/3.0_r_bl - rhokh_lcl(i,j) = zero - ! Use the BL entrainment parametrization as calculated above - rhokh_lcl(i,j) = min( rhokh_surf_ent(i,j), 5.0_r_bl) - end if - end do ! loop over j - end do ! loop over I + rhokh_lcl(i,j) = zero + ! Use the BL entrainment parametrization as calculated above + rhokh_lcl(i,j) = min( rhokh_surf_ent(i,j), 5.0_r_bl) + end if + end do ! loop over j + !end do ! loop over I !$OMP end do end if ! test on kprof_cu !----------------------------------------------------------------------- @@ -2826,18 +2858,20 @@ subroutine excf_nl_9c ( & !cdir collapse !$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 +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 k=ntml(i,j)+1 kh_top_factor(i,j) = max( 0.7_r_bl , one - sqrt( & - rhokh_surf_ent(i,j) / & - ( rho_mix(i,j,k)*w_h_top(i,j)*vkman*zh(i,j) ) ) ) + rhokh_surf_ent(i,j) / & + ( rho_mix(i,j,k)*w_h_top(i,j)*vkman*zh(i,j) ) ) ) if (l_use_var_fixes) then km_top_factor(i,j) = max( 0.7_r_bl, one-sqrt( rhokm_surf_ent(i,j) / & - ( rho_wet_tq(i,j,k-1)*w_m_top(i,j)*vkman*zh(i,j) ) ) ) + ( rho_wet_tq(i,j,k-1)*w_m_top(i,j)*vkman*zh(i,j) ) ) ) else km_top_factor(i,j) = max( 0.7_r_bl , one - sqrt( rhokm(i,j,k) / & - ( rho_wet_tq(i,j,k-1)*w_m_top(i,j)*vkman*zh(i,j) ) ) ) + ( rho_wet_tq(i,j,k-1)*w_m_top(i,j)*vkman*zh(i,j) ) ) ) end if scdepth(i,j) = zh(i,j) - zsml_base(i,j) factor = g1 * rho_mix(i,j,k) * v_top(i,j) *vkman *scdepth(i,j) @@ -2852,11 +2886,11 @@ subroutine excf_nl_9c ( & if ( factor > zero) then if (l_use_var_fixes) then km_sct_factor(i,j) = one - & - ( rhokm_top_ent(i,j) / factor )**1.25_r_bl + ( rhokm_top_ent(i,j) / factor )**1.25_r_bl ! 1.25=1/0.8 else km_sct_factor(i,j) = one - & - ( rhokm_top(i,j,k) / factor )**1.25_r_bl + ( rhokm_top(i,j,k) / factor )**1.25_r_bl end if else km_sct_factor(i,j) = one @@ -2882,22 +2916,23 @@ subroutine excf_nl_9c ( & end if factor = 0.75_r_bl * g1 * rho_wet_tq(i,j,k-1) * v_top_dsc(i,j) * & - vkman * dscdepth(i,j) + vkman * dscdepth(i,j) if ( factor > zero) then if (l_use_var_fixes) then km_dsct_factor(i,j) = one - & - ( rhokm_dsct_ent(i,j) / factor )**1.25_r_bl + ( rhokm_dsct_ent(i,j) / factor )**1.25_r_bl ! 1.25=1/0.8 else km_dsct_factor(i,j) = one - & - ( rhokm_top(i,j,k) / factor )**1.25_r_bl + ( rhokm_top(i,j,k) / factor )**1.25_r_bl end if else km_dsct_factor(i,j) = one end if end if - end do -end do + end do !I +end do !II +! end do !$OMP end do !----------------------------------------------------------------------- @@ -2910,12 +2945,12 @@ subroutine excf_nl_9c ( & !$OMP MASTER !cdir collapse -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - scbase(i,j) = .false. - nbdsc(i,j) = 0 - end do +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + scbase(i,j) = .false. + nbdsc(i,j) = 0 end do +!end do !$OMP end MASTER !$OMP BARRIER @@ -2932,150 +2967,150 @@ subroutine excf_nl_9c ( & c_tke = 1.33_r_bl/(vkman*g1) !$OMP do SCHEDULE(STATIC) -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 !cdir collapse - do j = jj, min(jj+bl_segment_size-1, pdims%j_end) - do i = pdims%i_start, pdims%i_end - - ! Calculate the height of u,v-level above the surface - ! *APL: z0m removed from z in K(z) - zk_uv = z_uv(i,j,k) - - ! Calculate the height of T,q-level above the surface - - zk_tq = z_tq(i,j,k-1) - - if ( zk_uv < zh(i,j) .and. & - zk_uv > zsml_base(i,j) ) then - z_pr = zk_uv - zsml_base(i,j) - zh_pr = zh(i,j) - zsml_base(i,j) - z_ratio = z_pr/zh_pr - if (flux_grad == LockWhelan2006) then - - rhokh_top(i,j,k) = 3.6_r_bl*vkman * rho_mix(i,j,k) * v_top(i,j) & - * zh_pr * (z_ratio**3) & - * (( one - z_ratio )**2) - - if ( .not. coupled(i,j) ) then - rhof2(i,j,k) = rho_mix(i,j,k) * one_half * z_ratio & - * 2.0_r_bl**( z_ratio**4 ) - if ( v_sum(i,j) > zero ) then - rhofsc(i,j,k) = 3.5_r_bl * rho_mix(i,j,k) & - * (v_top(i,j)/v_sum(i,j)) & - * (z_ratio**3) * (one-z_ratio) - end if + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) + !do i = pdims%i_start, pdims%i_end + + ! Calculate the height of u,v-level above the surface + ! *APL: z0m removed from z in K(z) + zk_uv = z_uv(i,j,k) + + ! Calculate the height of T,q-level above the surface + + zk_tq = z_tq(i,j,k-1) + + if ( zk_uv < zh(i,j) .and. & + zk_uv > zsml_base(i,j) ) then + z_pr = zk_uv - zsml_base(i,j) + zh_pr = zh(i,j) - zsml_base(i,j) + z_ratio = z_pr/zh_pr + if (flux_grad == LockWhelan2006) then + + rhokh_top(i,j,k) = 3.6_r_bl*vkman * rho_mix(i,j,k) * v_top(i,j) & + * zh_pr * (z_ratio**3) & + * (( one - z_ratio )**2) + + if ( .not. coupled(i,j) ) then + rhof2(i,j,k) = rho_mix(i,j,k) * one_half * z_ratio & + * 2.0_r_bl**( z_ratio**4 ) + if ( v_sum(i,j) > zero ) then + rhofsc(i,j,k) = 3.5_r_bl * rho_mix(i,j,k) & + * (v_top(i,j)/v_sum(i,j)) & + * (z_ratio**3) * (one-z_ratio) end if - - else ! Not LockWhelan2006 - ! max to avoid rounding errors giving small negative numbers - rhokh_top(i,j,k) = rho_mix(i,j,k) * v_top(i,j) * g1 * & - vkman * ((max(one - kh_sct_factor(i,j)*z_ratio,zero))**0.8_r_bl )& - * z_pr * z_ratio end if + else ! Not LockWhelan2006 + ! max to avoid rounding errors giving small negative numbers + rhokh_top(i,j,k) = rho_mix(i,j,k) * v_top(i,j) * g1 * & + vkman * ((max(one - kh_sct_factor(i,j)*z_ratio,zero))**0.8_r_bl )& + * z_pr * z_ratio end if - !------------------------------------------------------- - ! For LockWhelan2006, KM_TOP could be changed to match - ! the shape of KH_TOP. This has not been done on the - ! grounds that the change in shape arises with the - ! inclusion of the other non-gradient terms. - !------------------------------------------------------- - if ( zk_tq < zh(i,j) .and. & - zk_tq > zsml_base(i,j) ) then - z_pr = zk_tq - zsml_base(i,j) - zh_pr = zh(i,j) - zsml_base(i,j) - rhokm_top(i,j,k) = 0.75_r_bl * rho_wet_tq(i,j,k-1) * v_top(i,j) * & - g1 * vkman * & - ( (max(one - km_sct_factor(i,j)*z_pr/zh_pr, zero))**0.8_r_bl ) & - * z_pr * z_pr / zh_pr - ! PRANDTL=0.75 - 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,k)=rhokm_top(i,j,k)*c_tke*v_top(i,j)/zh(i,j) - end if + + end if + !------------------------------------------------------- + ! For LockWhelan2006, KM_TOP could be changed to match + ! the shape of KH_TOP. This has not been done on the + ! grounds that the change in shape arises with the + ! inclusion of the other non-gradient terms. + !------------------------------------------------------- + if ( zk_tq < zh(i,j) .and. & + zk_tq > zsml_base(i,j) ) then + z_pr = zk_tq - zsml_base(i,j) + zh_pr = zh(i,j) - zsml_base(i,j) + rhokm_top(i,j,k) = 0.75_r_bl * rho_wet_tq(i,j,k-1) * v_top(i,j) * & + g1 * vkman * & + ( (max(one - km_sct_factor(i,j)*z_pr/zh_pr, zero))**0.8_r_bl ) & + * z_pr * z_pr / zh_pr + ! PRANDTL=0.75 + 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,k)=rhokm_top(i,j,k)*c_tke*v_top(i,j)/zh(i,j) end if - !------------------------------------------------------------- - ! Add contribution to top-down mixing coefficient - ! profiles for decoupled stratocumulus layers when - ! one exists - !------------------------------------------------------------- - if ( zk_uv < zhsc(i,j) .and. & - zk_uv > zdsc_base(i,j) ) then - if (.not. scbase(i,j) ) then - scbase(i,j) = .true. - ! identifies lowest layer below which there is mixing - nbdsc(i,j) = k + end if + !------------------------------------------------------------- + ! Add contribution to top-down mixing coefficient + ! profiles for decoupled stratocumulus layers when + ! one exists + !------------------------------------------------------------- + if ( zk_uv < zhsc(i,j) .and. & + zk_uv > zdsc_base(i,j) ) then + if (.not. scbase(i,j) ) then + scbase(i,j) = .true. + ! identifies lowest layer below which there is mixing + nbdsc(i,j) = k + end if + !----------------------------------------------------------- + ! Calculate RHOK(H/M)_TOP, top-down turbulent mixing + ! profiles and add to any generated in the surface mixing + ! layer. + ! This is a variation on an up-side-down version of the + ! cubic surface-forced profiles above. Implement between + ! at least the top of the `surface layer' (at Z=0.1*ZH) and + ! ZHSC. + !----------------------------------------------------------- + z_pr = zk_uv - zdsc_base(i,j) + zh_pr = zhsc(i,j) - zdsc_base(i,j) + z_ratio = z_pr/zh_pr + + if (flux_grad == LockWhelan2006) then + + rhokh_top(i,j,k) = 3.6_r_bl*vkman * rho_mix(i,j,k) & + * v_top_dsc(i,j) * zh_pr * (z_ratio**3) & + * (( one - z_ratio )**2) + + rhof2(i,j,k) = rho_mix(i,j,k) * one_half * z_ratio & + * 2.0_r_bl**( z_ratio**4 ) + if ( v_sum_dsc(i,j) > zero ) then + rhofsc(i,j,k) = 3.5_r_bl * rho_mix(i,j,k) & + * (v_top_dsc(i,j)/v_sum_dsc(i,j)) & + * (z_ratio**3) * (one-z_ratio) end if - !----------------------------------------------------------- - ! Calculate RHOK(H/M)_TOP, top-down turbulent mixing - ! profiles and add to any generated in the surface mixing - ! layer. - ! This is a variation on an up-side-down version of the - ! cubic surface-forced profiles above. Implement between - ! at least the top of the `surface layer' (at Z=0.1*ZH) and - ! ZHSC. - !----------------------------------------------------------- - z_pr = zk_uv - zdsc_base(i,j) - zh_pr = zhsc(i,j) - zdsc_base(i,j) - z_ratio = z_pr/zh_pr - - if (flux_grad == LockWhelan2006) then - - rhokh_top(i,j,k) = 3.6_r_bl*vkman * rho_mix(i,j,k) & - * v_top_dsc(i,j) * zh_pr * (z_ratio**3) & - * (( one - z_ratio )**2) - rhof2(i,j,k) = rho_mix(i,j,k) * one_half * z_ratio & - * 2.0_r_bl**( z_ratio**4 ) - if ( v_sum_dsc(i,j) > zero ) then - rhofsc(i,j,k) = 3.5_r_bl * rho_mix(i,j,k) & - * (v_top_dsc(i,j)/v_sum_dsc(i,j)) & - * (z_ratio**3) * (one-z_ratio) - end if - - else ! Not LockWhelan2006 - ! max to avoid rounding errors giving small negative numbers - rhokh_top(i,j,k) = rhokh_top(i,j,k) + & - rho_mix(i,j,k)*v_top_dsc(i,j)*g1*vkman* & - ( (max(one - kh_dsct_factor(i,j)*z_ratio,zero))**0.8_r_bl ) & - * z_pr * z_ratio - end if - end if - !------------------------------------------------------------- - ! Now momentum - !------------------------------------------------------------- - if ( zk_tq < zhsc(i,j) .and. & - zk_tq > zdsc_base(i,j) ) then - z_pr = zk_tq - zdsc_base(i,j) - zh_pr = zhsc(i,j) - zdsc_base(i,j) + else ! Not LockWhelan2006 ! max to avoid rounding errors giving small negative numbers - rhokm_dsct = & - 0.75_r_bl*rho_wet_tq(i,j,k-1)*v_top_dsc(i,j)*g1*vkman* & - ( (max(one - km_dsct_factor(i,j)*z_pr/zh_pr,zero))**0.8_r_bl ) & - * z_pr * z_pr / zh_pr - 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,k) = tke_nl(i,j,k) + & - rhokm_dsct*c_tke*v_top_dsc(i,j)/dscdepth(i,j) - end if - rhokm_top(i,j,k) = rhokm_top(i,j,k) + rhokm_dsct + rhokh_top(i,j,k) = rhokh_top(i,j,k) + & + rho_mix(i,j,k)*v_top_dsc(i,j)*g1*vkman* & + ( (max(one - kh_dsct_factor(i,j)*z_ratio,zero))**0.8_r_bl ) & + * z_pr * z_ratio + end if + end if + !------------------------------------------------------------- + ! Now momentum + !------------------------------------------------------------- + if ( zk_tq < zhsc(i,j) .and. & + zk_tq > zdsc_base(i,j) ) then + z_pr = zk_tq - zdsc_base(i,j) + zh_pr = zhsc(i,j) - zdsc_base(i,j) + ! max to avoid rounding errors giving small negative numbers + rhokm_dsct = & + 0.75_r_bl*rho_wet_tq(i,j,k-1)*v_top_dsc(i,j)*g1*vkman* & + ( (max(one - km_dsct_factor(i,j)*z_pr/zh_pr,zero))**0.8_r_bl ) & + * z_pr * z_pr / zh_pr + 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,k) = tke_nl(i,j,k) + & + rhokm_dsct*c_tke*v_top_dsc(i,j)/dscdepth(i,j) end if - if (BL_diag%l_tke .and. var_diags_opt == original_vars) then + rhokm_top(i,j,k) = rhokm_top(i,j,k) + rhokm_dsct + end if + if (BL_diag%l_tke .and. var_diags_opt == original_vars) then + ! save 1/timescale for TKE diag, completed in bdy_expl2 + if ( zk_tq < zsml_top(i,j) .and. & + zk_tq > zsml_base(i,j) ) then + BL_diag%tke(i,j,k) = c_tke*v_top(i,j)/zh(i,j) + end if + if ( zk_tq < zhsc(i,j) .and. & + zk_tq > zdsc_base(i,j) ) then ! save 1/timescale for TKE diag, completed in bdy_expl2 - if ( zk_tq < zsml_top(i,j) .and. & - zk_tq > zsml_base(i,j) ) then - BL_diag%tke(i,j,k) = c_tke*v_top(i,j)/zh(i,j) - end if - if ( zk_tq < zhsc(i,j) .and. & - zk_tq > zdsc_base(i,j) ) then - ! save 1/timescale for TKE diag, completed in bdy_expl2 - BL_diag%tke(i,j,k) = max( BL_diag%tke(i,j,k), & - c_tke*v_top_dsc(i,j)/dscdepth(i,j) ) - end if + BL_diag%tke(i,j,k) = max( BL_diag%tke(i,j,k), & + c_tke*v_top_dsc(i,j)/dscdepth(i,j) ) end if + end if - end do + !end do end do end do end do @@ -3095,109 +3130,109 @@ subroutine excf_nl_9c ( & c_tke = 1.33_r_bl/(vkman*c_ws**two_thirds) !$OMP do SCHEDULE(STATIC) - 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 !cdir collapse - 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 - zk_uv = z_uv(i,j,k) - zk_tq = z_tq(i,j,k-1) + zk_uv = z_uv(i,j,k) + zk_tq = z_tq(i,j,k-1) - if (fb_surf(i,j) >= zero) then + if (fb_surf(i,j) >= zero) then - ! Calculate the free-convective scaling velocity at z(k) + ! Calculate the free-convective scaling velocity at z(k) - if (coupled(i,j)) then ! coupled - if ( entr_smooth_dec == entr_taper_zh ) then - wstar3 = ( svl_diff_frac(i,j) * zhsc(i,j) & - + (one-svl_diff_frac(i,j)) * zh(i,j) ) *fb_surf(i,j) - else - wstar3 = zhsc(i,j) * fb_surf(i,j) - end if + if (coupled(i,j)) then ! coupled + if ( entr_smooth_dec == entr_taper_zh ) then + wstar3 = ( svl_diff_frac(i,j) * zhsc(i,j) & + + (one-svl_diff_frac(i,j)) * zh(i,j) ) *fb_surf(i,j) else - wstar3 = zh(i,j) * fb_surf(i,j) + wstar3 = zhsc(i,j) * fb_surf(i,j) end if + else + wstar3 = zh(i,j) * fb_surf(i,j) + end if - if (zk_uv <= 0.1_r_bl*zh(i,j)) then - ! Surface layer calculation - w_s_cubed_uv = 10.0_r_bl*c_ws * zk_uv * fb_surf(i,j) - else - ! Outer layer calculation - w_s_cubed_uv = c_ws * wstar3 - end if + if (zk_uv <= 0.1_r_bl*zh(i,j)) then + ! Surface layer calculation + w_s_cubed_uv = 10.0_r_bl*c_ws * zk_uv * fb_surf(i,j) + else + ! Outer layer calculation + w_s_cubed_uv = c_ws * wstar3 + end if - if (zk_tq <= 0.1_r_bl*zh(i,j)) then - ! Surface layer calculation - w_s_cubed_tq = 10.0_r_bl*c_ws * zk_tq * fb_surf(i,j) - else - ! Outer layer calculation - w_s_cubed_tq = c_ws * wstar3 - end if + if (zk_tq <= 0.1_r_bl*zh(i,j)) then + ! Surface layer calculation + w_s_cubed_tq = 10.0_r_bl*c_ws * zk_tq * fb_surf(i,j) + else + ! Outer layer calculation + w_s_cubed_tq = c_ws * wstar3 + end if - ! Turbulent velocity scale for scalars + ! Turbulent velocity scale for scalars - w_m_neut = ( v_s(i,j)*v_s(i,j)*v_s(i,j) + w_s_cubed_uv ) & - **one_third - w_h_uv = w_m_neut/pr_neut + w_m_neut = ( v_s(i,j)*v_s(i,j)*v_s(i,j) + w_s_cubed_uv ) & + **one_third + w_h_uv = w_m_neut/pr_neut - ! Also calc on TQ levels for W_M_TQ - w_m_neut = ( v_s(i,j)*v_s(i,j)*v_s(i,j) + w_s_cubed_tq ) & - **one_third - w_h_tq = w_m_neut/pr_neut + ! Also calc on TQ levels for W_M_TQ + w_m_neut = ( v_s(i,j)*v_s(i,j)*v_s(i,j) + w_s_cubed_tq ) & + **one_third + w_h_tq = w_m_neut/pr_neut - ! The calculations below involve v_s**4. In very rare circumstances - ! v_s can be order(10^-11), therefore v_s**4 is order(10^-44) - ! which is outside the range of single precision calculations - ! Hence we create a double precision version to enforce the correct - ! calculation - v_s_dbl(i,j) = v_s(i,j) + ! The calculations below involve v_s**4. In very rare circumstances + ! v_s can be order(10^-11), therefore v_s**4 is order(10^-44) + ! which is outside the range of single precision calculations + ! Hence we create a double precision version to enforce the correct + ! calculation + v_s_dbl(i,j) = v_s(i,j) - ! Turbulent Prandtl number and velocity scale for scalars + ! Turbulent Prandtl number and velocity scale for scalars - Prandtl = pr_neut* & - ( v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j) + & - (one/(c_ws*25.0_r_bl))*w_s_cubed_tq*w_m_neut ) / & - ( v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j) + & - (one/(c_ws*25.0_r_bl))*(pr_neut/pr_conv)*w_s_cubed_tq* & - w_m_neut ) + Prandtl = pr_neut* & + ( v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j) + & + (one/(c_ws*25.0_r_bl))*w_s_cubed_tq*w_m_neut ) / & + ( v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j) + & + (one/(c_ws*25.0_r_bl))*(pr_neut/pr_conv)*w_s_cubed_tq* & + w_m_neut ) - w_m_tq = Prandtl * w_h_tq + w_m_tq = Prandtl * w_h_tq - if ( zk_uv < zh(i,j) ) then - !--------------------------------------------------------- - ! Calculate RHOKH(w_h,z/z_h) - !--------------------------------------------------------- + if ( zk_uv < zh(i,j) ) then + !--------------------------------------------------------- + ! Calculate RHOKH(w_h,z/z_h) + !--------------------------------------------------------- - rhokh(i,j,k) = rho_mix(i,j,k) * w_h_uv * vkman * zk_uv * & - ( one - ( zk_uv / zh(i,j) ) ) * & - ( one - ( zk_uv / zh(i,j) ) ) + rhokh(i,j,k) = rho_mix(i,j,k) * w_h_uv * vkman * zk_uv * & + ( one - ( zk_uv / zh(i,j) ) ) * & + ( one - ( zk_uv / zh(i,j) ) ) - end if - if ( zk_tq < zh(i,j) ) then - !--------------------------------------------------------- - ! Calculate RHOKM(w_m,z/z_h) - !--------------------------------------------------------- + end if + if ( zk_tq < zh(i,j) ) then + !--------------------------------------------------------- + ! Calculate RHOKM(w_m,z/z_h) + !--------------------------------------------------------- - rhokm(i,j,k) = rho_wet_tq(i,j,k-1)*w_m_tq*vkman*zk_tq * & - ( one - ( zk_tq / zh(i,j) ) ) * & - ( one - ( zk_tq / zh(i,j) ) ) + rhokm(i,j,k) = rho_wet_tq(i,j,k-1)*w_m_tq*vkman*zk_tq * & + ( one - ( zk_tq / zh(i,j) ) ) * & + ( one - ( zk_tq / zh(i,j) ) ) - 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,k) = tke_nl(i,j,k) + & - rhokm(i,j,k)*c_tke*w_m_tq/zh(i,j) - end if + 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,k) = tke_nl(i,j,k) + & + rhokm(i,j,k)*c_tke*w_m_tq/zh(i,j) end if - if (BL_diag%l_tke .and. var_diags_opt == original_vars) then - ! save 1/timescale for TKE diag, completed in bdy_expl2 - if ( zk_tq < zsml_top(i,j) ) then - BL_diag%tke(i,j,k) = max( BL_diag%tke(i,j,k), & - c_tke*w_m_tq/zh(i,j) ) - end if + end if + if (BL_diag%l_tke .and. var_diags_opt == original_vars) then + ! save 1/timescale for TKE diag, completed in bdy_expl2 + if ( zk_tq < zsml_top(i,j) ) then + BL_diag%tke(i,j,k) = max( BL_diag%tke(i,j,k), & + c_tke*w_m_tq/zh(i,j) ) end if end if - end do + end if + !end do end do end do end do @@ -3220,148 +3255,148 @@ subroutine excf_nl_9c ( & c_tke = 1.33_r_bl/(vkman*c_ws**two_thirds) !$OMP do SCHEDULE(STATIC) - 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 !cdir collapse - 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 - ! Calculate the height of u,v-level above the surface - ! *APL: z0m removed from z in K(z) - zk_uv = z_uv(i,j,k) + ! Calculate the height of u,v-level above the surface + ! *APL: z0m removed from z in K(z) + zk_uv = z_uv(i,j,k) - ! Calculate the height of T,q-level above the surface + ! Calculate the height of T,q-level above the surface - zk_tq = z_tq(i,j,k-1) + zk_tq = z_tq(i,j,k-1) - if (fb_surf(i,j) >= zero) then + if (fb_surf(i,j) >= zero) then - ! Calculate the free-convective scaling velocity at z(k) + ! Calculate the free-convective scaling velocity at z(k) - if (zk_uv <= 0.1_r_bl*zh(i,j)) then + if (zk_uv <= 0.1_r_bl*zh(i,j)) then - ! Surface layer calculation + ! Surface layer calculation - w_s_cubed_uv = 10.0_r_bl*c_ws * zk_uv * fb_surf(i,j) - else + w_s_cubed_uv = 10.0_r_bl*c_ws * zk_uv * fb_surf(i,j) + else - ! Outer layer calculation + ! Outer layer calculation - if (coupled(i,j)) then ! coupled and cloudy - if ( entr_smooth_dec == entr_taper_zh ) then - w_s_cubed_uv = c_ws & - * ( svl_diff_frac(i,j) * zhsc(i,j) & - + (one-svl_diff_frac(i,j)) * zh(i,j) ) * fb_surf(i,j) - else - w_s_cubed_uv = c_ws * zhsc(i,j) * fb_surf(i,j) - end if + if (coupled(i,j)) then ! coupled and cloudy + if ( entr_smooth_dec == entr_taper_zh ) then + w_s_cubed_uv = c_ws & + * ( svl_diff_frac(i,j) * zhsc(i,j) & + + (one-svl_diff_frac(i,j)) * zh(i,j) ) * fb_surf(i,j) else - w_s_cubed_uv = c_ws * zh(i,j) * fb_surf(i,j) + w_s_cubed_uv = c_ws * zhsc(i,j) * fb_surf(i,j) end if + else + w_s_cubed_uv = c_ws * zh(i,j) * fb_surf(i,j) end if + end if - if (zk_tq <= 0.1_r_bl*zh(i,j)) then + if (zk_tq <= 0.1_r_bl*zh(i,j)) then - ! Surface layer calculation + ! Surface layer calculation - w_s_cubed_tq = 10.0_r_bl*c_ws * zk_tq * fb_surf(i,j) - else + w_s_cubed_tq = 10.0_r_bl*c_ws * zk_tq * fb_surf(i,j) + else - ! Outer layer calculation + ! Outer layer calculation - if (coupled(i,j)) then ! coupled and cloudy - if ( entr_smooth_dec == entr_taper_zh ) then - w_s_cubed_tq = c_ws & - * ( svl_diff_frac(i,j) * zhsc(i,j) & - + (one-svl_diff_frac(i,j)) * zh(i,j) ) * fb_surf(i,j) - else - w_s_cubed_tq = c_ws * zhsc(i,j) * fb_surf(i,j) - end if + if (coupled(i,j)) then ! coupled and cloudy + if ( entr_smooth_dec == entr_taper_zh ) then + w_s_cubed_tq = c_ws & + * ( svl_diff_frac(i,j) * zhsc(i,j) & + + (one-svl_diff_frac(i,j)) * zh(i,j) ) * fb_surf(i,j) else - w_s_cubed_tq = c_ws * zh(i,j) * fb_surf(i,j) + w_s_cubed_tq = c_ws * zhsc(i,j) * fb_surf(i,j) end if + else + w_s_cubed_tq = c_ws * zh(i,j) * fb_surf(i,j) end if + end if - ! Turbulent velocity scale for momentum - - w_m_uv = (v_s(i,j)*v_s(i,j)*v_s(i,j) + w_s_cubed_uv) & - **one_third - - w_m_tq = (v_s(i,j)*v_s(i,j)*v_s(i,j) + w_s_cubed_tq) & - **one_third - - ! The calculations below involve v_s**4. In very rare circumstances - ! v_s can be order(10^-11), therefore v_s**4 is order(10^-44) - ! which is outside the range of single precision calculations - ! Hence we create a double precision version to enforce the correct - ! calculation - v_s_dbl(i,j) = v_s(i,j) - - ! Turbulent Prandtl number and velocity scale for scalars - - Prandtl = pr_neut* & - ( v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j) + & - (one/(c_ws*25.0_r_bl))*w_s_cubed_uv*w_m_uv ) / & - ( v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j) + & - (one/(c_ws*25.0_r_bl))*(pr_neut/pr_conv)*w_s_cubed_uv* & - w_m_uv ) - w_h_uv = w_m_uv / Prandtl - - if ( zk_uv < zh(i,j) ) then - !--------------------------------------------------------- - ! Calculate RHOKH(w_h,z/z_h) - !--------------------------------------------------------- - - rhokh(i,j,k) = rho_mix(i,j,k) * w_h_uv * vkman * zk_uv * & - ( one - kh_top_factor(i,j) * ( zk_uv / zh(i,j) ) ) * & - ( one - kh_top_factor(i,j) * ( zk_uv / zh(i,j) ) ) - else if ( kprof_cu == klcl_entr .and. cumulus(i,j) ) then - if ( zk_uv < zsml_top(i,j) ) then - ! Exponential decay from ZH but tends to zero - ! at zsml_top - rhokh(i,j,k) = rhokh_lcl(i,j) * & - exp(-(zk_uv-zh(i,j))/cu_depth_scale(i,j)) * & - (one-(zk_uv-zh(i,j))/(zsml_top(i,j)-zh(i,j))) - end if + ! Turbulent velocity scale for momentum + + w_m_uv = (v_s(i,j)*v_s(i,j)*v_s(i,j) + w_s_cubed_uv) & + **one_third + + w_m_tq = (v_s(i,j)*v_s(i,j)*v_s(i,j) + w_s_cubed_tq) & + **one_third + + ! The calculations below involve v_s**4. In very rare circumstances + ! v_s can be order(10^-11), therefore v_s**4 is order(10^-44) + ! which is outside the range of single precision calculations + ! Hence we create a double precision version to enforce the correct + ! calculation + v_s_dbl(i,j) = v_s(i,j) + + ! Turbulent Prandtl number and velocity scale for scalars + + Prandtl = pr_neut* & + ( v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j) + & + (one/(c_ws*25.0_r_bl))*w_s_cubed_uv*w_m_uv ) / & + ( v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j) + & + (one/(c_ws*25.0_r_bl))*(pr_neut/pr_conv)*w_s_cubed_uv* & + w_m_uv ) + w_h_uv = w_m_uv / Prandtl + + if ( zk_uv < zh(i,j) ) then + !--------------------------------------------------------- + ! Calculate RHOKH(w_h,z/z_h) + !--------------------------------------------------------- + + rhokh(i,j,k) = rho_mix(i,j,k) * w_h_uv * vkman * zk_uv * & + ( one - kh_top_factor(i,j) * ( zk_uv / zh(i,j) ) ) * & + ( one - kh_top_factor(i,j) * ( zk_uv / zh(i,j) ) ) + else if ( kprof_cu == klcl_entr .and. cumulus(i,j) ) then + if ( zk_uv < zsml_top(i,j) ) then + ! Exponential decay from ZH but tends to zero + ! at zsml_top + rhokh(i,j,k) = rhokh_lcl(i,j) * & + exp(-(zk_uv-zh(i,j))/cu_depth_scale(i,j)) * & + (one-(zk_uv-zh(i,j))/(zsml_top(i,j)-zh(i,j))) end if - if ( zk_tq < zh(i,j) ) then - !--------------------------------------------------------- - ! Calculate RHOKM(w_m,z/z_h) - !--------------------------------------------------------- + end if + if ( zk_tq < zh(i,j) ) then + !--------------------------------------------------------- + ! Calculate RHOKM(w_m,z/z_h) + !--------------------------------------------------------- - rhokm(i,j,k) = rho_wet_tq(i,j,k-1)*w_m_tq*vkman* zk_tq * & - ( one - km_top_factor(i,j) * ( zk_tq / zh(i,j) ) ) * & - ( one - km_top_factor(i,j) * ( zk_tq / zh(i,j) ) ) + rhokm(i,j,k) = rho_wet_tq(i,j,k-1)*w_m_tq*vkman* zk_tq * & + ( one - km_top_factor(i,j) * ( zk_tq / zh(i,j) ) ) * & + ( one - km_top_factor(i,j) * ( zk_tq / zh(i,j) ) ) - if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then + 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,k) = tke_nl(i,j,k) + & + rhokm(i,j,k)*c_tke*w_m_tq/zh(i,j) + end if + else if ( kprof_cu == klcl_entr .and. cumulus(i,j) ) then + if ( zk_tq < zsml_top(i,j) ) then + ! Exponential decay from ZH but tends to zero + ! at zsml_top + rhokm(i,j,k) = prandtl_top(i,j) * rhokh_lcl(i,j) * & + exp(-(zk_tq-zh(i,j))/cu_depth_scale(i,j)) * & + (one-(zk_tq-zh(i,j))/(zsml_top(i,j)-zh(i,j))) + 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,k) = tke_nl(i,j,k) + & - rhokm(i,j,k)*c_tke*w_m_tq/zh(i,j) - end if - else if ( kprof_cu == klcl_entr .and. cumulus(i,j) ) then - if ( zk_tq < zsml_top(i,j) ) then - ! Exponential decay from ZH but tends to zero - ! at zsml_top - rhokm(i,j,k) = prandtl_top(i,j) * rhokh_lcl(i,j) * & - exp(-(zk_tq-zh(i,j))/cu_depth_scale(i,j)) * & - (one-(zk_tq-zh(i,j))/(zsml_top(i,j)-zh(i,j))) - 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,k) = tke_nl(i,j,k) + & - rhokm(i,j,k)*c_tke*w_m_tq/zh(i,j) - end if + tke_nl(i,j,k) = tke_nl(i,j,k) + & + rhokm(i,j,k)*c_tke*w_m_tq/zh(i,j) end if end if - if (BL_diag%l_tke .and. var_diags_opt == original_vars) then - ! save 1/timescale for TKE diag, completed in bdy_expl2 - if ( zk_tq < zsml_top(i,j) ) then - BL_diag%tke(i,j,k) = max( BL_diag%tke(i,j,k), & - c_tke*w_m_tq/zh(i,j) ) - end if + end if + if (BL_diag%l_tke .and. var_diags_opt == original_vars) then + ! save 1/timescale for TKE diag, completed in bdy_expl2 + if ( zk_tq < zsml_top(i,j) ) then + BL_diag%tke(i,j,k) = max( BL_diag%tke(i,j,k), & + c_tke*w_m_tq/zh(i,j) ) end if - end if - end do + + end if + !end do end do end do end do @@ -3374,59 +3409,59 @@ subroutine excf_nl_9c ( & ng_stress == BrownGrant97_original ) then !$OMP do SCHEDULE(STATIC) - 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 !cdir collapse - do j = jj, min(jj+bl_segment_size-1, pdims%j_end) - do i = pdims%i_start, pdims%i_end - zk_tq = z_tq(i,j,k-1) ! stresses are calc on theta-levs - if ( fb_surf(i,j) > zero .and. zk_tq < zh(i,j) ) then - !--------------------------------------------------------- - ! Calculate non-gradient stress function - ! (Brown and Grant 1997) - ! Shape function chosen such that non-gradient stress - ! either goes to zero at 0.1*ZH and ZH (ng_stress==BrownGrant97 - ! and ng_stress==BrownGrant97_limited), or goes to zero at - ! at the ground surface and ZH (ng_stress==BrownGrant97_original) - !--------------------------------------------------------- - ng_stress_calculate = .false. - - if (ng_stress == BrownGrant97 .or. & - ng_stress == BrownGrant97_limited) then - if ( zk_tq > 0.1_r_bl*zh(i,j) ) then - z_pr = zk_tq - 0.1_r_bl*zh(i,j) - zh_pr = 0.9_r_bl*zh(i,j) - ng_stress_calculate = .true. - end if - else if (ng_stress==BrownGrant97_original) then - z_pr = zk_tq - zh_pr = zh(i,j) + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) + !do i = pdims%i_start, pdims%i_end + zk_tq = z_tq(i,j,k-1) ! stresses are calc on theta-levs + if ( fb_surf(i,j) > zero .and. zk_tq < zh(i,j) ) then + !--------------------------------------------------------- + ! Calculate non-gradient stress function + ! (Brown and Grant 1997) + ! Shape function chosen such that non-gradient stress + ! either goes to zero at 0.1*ZH and ZH (ng_stress==BrownGrant97 + ! and ng_stress==BrownGrant97_limited), or goes to zero at + ! at the ground surface and ZH (ng_stress==BrownGrant97_original) + !--------------------------------------------------------- + ng_stress_calculate = .false. + + if (ng_stress == BrownGrant97 .or. & + ng_stress == BrownGrant97_limited) then + if ( zk_tq > 0.1_r_bl*zh(i,j) ) then + z_pr = zk_tq - 0.1_r_bl*zh(i,j) + zh_pr = 0.9_r_bl*zh(i,j) ng_stress_calculate = .true. end if + else if (ng_stress==BrownGrant97_original) then + z_pr = zk_tq + zh_pr = zh(i,j) + ng_stress_calculate = .true. + end if - if ( ng_stress_calculate ) then - ! Outer layer calculation - if (coupled(i,j)) then ! coupled and cloudy - if ( entr_smooth_dec == entr_taper_zh ) then - wstar3 = ( svl_diff_frac(i,j) * zhsc(i,j) & - + (one-svl_diff_frac(i,j)) * zh(i,j) ) *fb_surf(i,j) - else - wstar3 = zhsc(i,j) * fb_surf(i,j) - end if + if ( ng_stress_calculate ) then + ! Outer layer calculation + if (coupled(i,j)) then ! coupled and cloudy + if ( entr_smooth_dec == entr_taper_zh ) then + wstar3 = ( svl_diff_frac(i,j) * zhsc(i,j) & + + (one-svl_diff_frac(i,j)) * zh(i,j) ) *fb_surf(i,j) else - wstar3 = zh(i,j) * fb_surf(i,j) + wstar3 = zhsc(i,j) * fb_surf(i,j) end if + else + wstar3 = zh(i,j) * fb_surf(i,j) + end if - ! Use the Holtslag and Boville velocity scale for - ! non-gradient stress stability dependence, as in BG97 - w_m_hb_3 = v_s(i,j)*v_s(i,j)*v_s(i,j) + 0.6_r_bl*wstar3 - f_ngstress(i,j,k) =(rho_wet_tq(i,j,k-1)/rhostar_gb(i,j)) & - * s_m * ( a_ngs * wstar3 / w_m_hb_3 ) & - * ( z_pr / zh_pr ) * ( one - ( z_pr / zh_pr ) ) * & - ( one - ( z_pr / zh_pr ) ) - end if ! ng_stress_calculate - end if ! fb_surf>0 and z0 and z