From cf2453bd5a76458b50587f37b039971ce5216676 Mon Sep 17 00:00:00 2001 From: William Putman Date: Mon, 31 Aug 2026 08:43:29 -0400 Subject: [PATCH 1/2] 015e just updated the logfile output for Rayleigh damping timescales --- model/fv_dynamics.F90 | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/model/fv_dynamics.F90 b/model/fv_dynamics.F90 index 667049fec..91c74945f 100644 --- a/model/fv_dynamics.F90 +++ b/model/fv_dynamics.F90 @@ -1053,11 +1053,15 @@ subroutine Rayleigh_fast(dt, npx, npy, npz, pfull, tau, u, v, w, & do k=1, npz if ( pfull(k) < rf_cutoff ) then rff(k) = dt/tau0*sin(0.5*pi*log(rf_cutoff/pfull(k))/log(rf_cutoff/ptop))**2 -! Re-FACTOR rf - if( is_master() ) write(6,*) k, 0.01*pfull(k), dt/(rff(k)*sday) kmax = k rff(k) = 1.d0 / (1.0d0+rff(k)) rf(k) = rff(k) + ! Added safety check to prevent a log(1.0) divide-by-zero crash + if( is_master() ) then + if (rff(k) < 1.0d0) then + write(6,*) k, 0.01*pfull(k), (-dt / log(rff(k))) / sday + endif + endif else exit endif From 8f6a2d65f1ed839400c65237d473531511b463a6 Mon Sep 17 00:00:00 2001 From: William Putman Date: Wed, 30 Sep 2026 08:18:16 -0400 Subject: [PATCH 2/2] important hires tracer optimizations --- model/dyn_core.F90 | 60 +++- model/fv_control.F90 | 3 + model/fv_dynamics.F90 | 114 +++---- model/fv_mapz.F90 | 722 +++++++++++++++++++++++++++++++++++++++++- model/fv_tracer2d.F90 | 600 ++++++++++++++++++++++------------- model/nh_utils.F90 | 46 ++- tools/fv_mp_mod.F90 | 52 ++- 7 files changed, 1268 insertions(+), 329 deletions(-) diff --git a/model/dyn_core.F90 b/model/dyn_core.F90 index 452d75c48..24f52a785 100644 --- a/model/dyn_core.F90 +++ b/model/dyn_core.F90 @@ -110,7 +110,7 @@ module dyn_core_mod use mpp_parameter_mod, only: CORNER use fv_mp_mod, only: is_master use fv_mp_mod, only: start_group_halo_update, complete_group_halo_update - use fv_mp_mod, only: group_halo_update_type + use fv_mp_mod, only: group_halo_update_type, mp_ireduce_max use sw_core_mod, only: c_sw, d_sw, d2a2c_vect use a2b_edge_mod, only: a2b_ord2, a2b_ord4 use nh_core_mod, only: Riem_Solver3, Riem_Solver_C, update_dz_c, update_dz_d, nest_halo_nh @@ -162,9 +162,9 @@ module dyn_core_mod subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, k_split, n_split, zvir, cp, akap, cappa, grav, hydrostatic, & u, v, w, delz, pt, q, delp, pe, pk, phis, varflt, ws, omga, ptop, pfull, ua, va, & - dudt_rf, dvdt_rf, dwdt_rf, uc, vc, mfx, mfy, cx, cy, pkz, peln, q_con, ak, bk, dpx, & + dudt_rf, dvdt_rf, dwdt_rf, uc, vc, mfx, mfy, cx, cy, cmax_z, imax_req, q_split, pkz, peln, q_con, ak, bk, dpx, & ks, gridstruct, flagstruct, neststruct, idiag, bd, domain, & - init_step, i_pack, end_step, diss_est,time_total) + init_step, i_pack, end_step, diss_est, time_total) integer, intent(IN) :: npx integer, intent(IN) :: npy @@ -225,6 +225,10 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, k_split, n_split, zvir, c real, intent(inout), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: ua, va real, intent(inout):: q_con(bd%isd:, bd%jsd:, 1:) + real, intent(out):: cmax_z(npz) + integer, intent(out):: imax_req + integer, intent(in):: q_split + ! The Flux capacitors: accumulated Mass flux arrays real(kind=8), intent(inout):: mfx(bd%is:bd%ie+1, bd%js:bd%je, npz) real(kind=8), intent(inout):: mfy(bd%is:bd%ie , bd%js:bd%je+1, npz) @@ -242,6 +246,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, k_split, n_split, zvir, c real, allocatable, dimension(:,:,:):: pem, heat_source ! Auto 1D & 2D arrays: real, dimension(bd%isd:bd%ied,bd%jsd:bd%jed):: ws3, z_rat + real:: cmax real:: dp_ref(npz) real:: zs(bd%isd:bd%ied,bd%jsd:bd%jed) !< surface height (m) real:: p1d(bd%is:bd%ie) @@ -393,8 +398,6 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, k_split, n_split, zvir, c endif endif - - !----------------------------------------------------- do it=1,n_split !----------------------------------------------------- @@ -826,12 +829,48 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, k_split, n_split, zvir, c enddo enddo endif + enddo ! end openMP k-loop if ( (.not. hydrostatic) .and. flagstruct%fv_debug ) & call prt_mxm('W_dsw ', w, is, ie , js, je , ng, npz, 1., gridstruct%area_64, domain) call timing_off('d_sw') + + ! Initialize the request handle so the later WAIT call doesn't crash + ! if the reduction is never started. + if ( (it==n_split) .and. nq > 0 .and. .not. flagstruct%inline_q ) then + + ! --- 1. NON-BLOCKING CMAX REDUCTION START --- + ! A. Compute and store local cmax_z for this processor + !$OMP parallel do default(none) shared(npz,is,ie,js,je,cx,cy,cmax_z,gridstruct) & + !$OMP private(i,j,k,cmax) + do k = 1, npz + ! Compute local cmax_z (cx and cy are already hot in cache) + cmax_z(k) = 0. + do j=js,je + do i=is,ie + cmax = max(abs(cx(i,j,k)),abs(cy(i,j,k)))+1.-gridstruct%sin_sg(i,j,5) + cmax_z(k) = max( cmax_z(k), cmax ) + enddo + enddo + enddo + if ( q_split == 0 ) then + ! B. Initiate the non-blocking global reduction + call timing_on('COMM_TOTAL') + call timing_on('COMM_TRACER_MAX') + call mp_ireduce_max(cmax_z, npz, imax_req) + call timing_off('COMM_TRACER_MAX') + call timing_off('COMM_TOTAL') + endif + + call timing_on('COMM_TOTAL') + call timing_on('COMM_TRACER') + call start_group_halo_update(i_pack(10), q, domain) + call timing_off('COMM_TRACER') + call timing_off('COMM_TOTAL') + endif + if( flagstruct%fill_dp ) call mix_dp(hydrostatic, w, delp, pt, npz, ak, bk, .false., flagstruct%fv_debug, bd) if ( flagstruct%fv_debug ) & @@ -1042,6 +1081,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, k_split, n_split, zvir, c ! *** Inline Rayleigh friction here if( flagstruct%RF_fast .and. flagstruct%tau > 0. ) then + call timing_on('RAY_FAST') dudt_rf = u(is:ie,js:je,:) dvdt_rf = v(is:ie,js:je,:) if (.not. hydrostatic) dwdt_rf = w(is:ie,js:je,:) @@ -1050,6 +1090,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, k_split, n_split, zvir, c dudt_rf = ( u(is:ie,js:je,:) - dudt_rf)/dt dvdt_rf = ( v(is:ie,js:je,:) - dvdt_rf)/dt if (.not. hydrostatic) dwdt_rf = ( w(is:ie,js:je,:) - dwdt_rf)/dt + call timing_off('RAY_FAST') endif ! *** Inline Beljaars turbulent-orographic-form-drag here @@ -1195,20 +1236,11 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, k_split, n_split, zvir, c !----------------------------------------------------- enddo ! time split loop !----------------------------------------------------- - if ( nq > 0 .and. .not. flagstruct%inline_q ) then - call timing_on('COMM_TOTAL') - call timing_on('COMM_TRACER') - call start_group_halo_update(i_pack(10), q, domain) - call timing_off('COMM_TRACER') - call timing_off('COMM_TOTAL') - endif - if ( flagstruct%fv_debug ) then if(is_master()) write(*,*) 'End of n_split loop' endif - if ( n_con/=0 .and. flagstruct%d_con > 1.e-5 ) then nf_ke = min(3, flagstruct%nord+1) call del2_cubed(heat_source, cnst_0p20*gridstruct%da_min, gridstruct, domain, npx, npy, npz, nf_ke, bd) diff --git a/model/fv_control.F90 b/model/fv_control.F90 index 17c474b9f..f909aaa92 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -834,6 +834,9 @@ subroutine run_setup(Atm, dt_atmos, grids_on_this_pe, p_split) write(*,199) 'Using p_split = ', p_split endif + ! Safety check on q_split and z_tracer + if ( q_split /= 0) z_tracer = .false. + if (Atm(n)%neststruct%nested) then do i=1,n-1 if (Atm(i)%grid_number == parent_grid_num) then diff --git a/model/fv_dynamics.F90 b/model/fv_dynamics.F90 index 91c74945f..ebad58464 100644 --- a/model/fv_dynamics.F90 +++ b/model/fv_dynamics.F90 @@ -245,20 +245,15 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, & type(fv_diag_type), intent(IN) :: idiag ! Local Arrays -#ifdef SINGLE_FV ! R8 Mass flux arrays: the "Flux Capacitor" real(kind=8) :: mfxR8(bd%is:bd%ie+1, bd%js:bd%je, npz) real(kind=8) :: mfyR8(bd%is:bd%ie , bd%js:bd%je+1, npz) ! R8 Courant number arrays real(kind=8) :: cxR8(bd%is:bd%ie+1, bd%jsd:bd%jed, npz) real(kind=8) :: cyR8(bd%isd:bd%ied ,bd%js:bd%je+1, npz) -#endif -! Local Mass flux arrays: the "Flux Capacitor" - real :: mfxL(bd%is:bd%ie+1, bd%js:bd%je, npz) - real :: mfyL(bd%is:bd%ie , bd%js:bd%je+1, npz) -! Local Courant number arrays - real :: cxL(bd%is:bd%ie+1, bd%jsd:bd%jed, npz) - real :: cyL(bd%isd:bd%ied ,bd%js:bd%je+1, npz) + + integer :: imax_req + real :: cmax_z(npz) ! More Local arrays real:: ws(bd%is:bd%ie,bd%js:bd%je) real(kind=8):: te_2d(bd%is:bd%ie,bd%js:bd%je) @@ -270,8 +265,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, & real, dimension(bd%is:bd%ie):: cvm real, allocatable :: dp1(:,:,:), dtdt_m(:,:,:), cappa(:,:,:) real(kind=8), allocatable :: psx(:,:) - !real(kind=8), allocatable :: dpx(:,:) - real(kind=8), allocatable :: dpx(:,:,:) !needed for OpenMP + real(kind=8), allocatable :: dpx(:,:,:) real:: akap, rdg, ph1, ph2, mdt, gam, amdt, u0 integer:: kord_tracer(ncnst) integer :: i,j,k, n, iq, n_map, nq, nwat @@ -604,14 +598,18 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, & endif !DryMassRoundoffControl - allocate(psx(isd:ied,jsd:jed),dpx(is:ie,js:je,npz)) - psx(:,:) = 0.0 - do j=js,je - do i=is,ie - psx(i,j) = pe(i,npz+1,j) - dpx(i,j,:) = 0.0 - enddo - enddo + allocate(dpx(is:ie,js:je,npz)) + dpx(:,:,:) = 0.0 + if (hydrostatic) then + allocate(psx(isd:ied,jsd:jed)) + psx(:,:) = 0.0 + do j=js,je + do i=is,ie + psx(i,j) = pe(i,npz+1,j) + enddo + enddo + endif + call timing_on('FV_DYN_LOOP') do n_map=1, k_split ! first level of time-split call timing_on('COMM_TOTAL') @@ -645,30 +643,12 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, & call dyn_core(npx, npy, npz, ng, sphum, nq, mdt, k_split, n_split, zvir, cp_air, akap, cappa, grav, hydrostatic, & u, v, w, delz, pt, q, delp, pe, pk, phis, varflt, ws, omga, ptop, pfull, ua, va, & dudt_rf, dvdt_rf, dwdt_rf, uc, vc, & -#ifdef SINGLE_FV - mfxR8, mfyR8, cxR8, cyR8, & -#else - mfxL, mfyL, cxL, cyL, & -#endif + mfxR8, mfyR8, cxR8, cyR8, cmax_z, imax_req, q_split, & pkz, peln, q_con, ak, bk, dpx, ks, & gridstruct, flagstruct, neststruct, idiag, bd, & domain, n_map==1, i_pack, last_step, diss_est,time_total) call timing_off('DYN_CORE') -!MassFluxRoundoffControl -#ifdef SINGLE_FV - mfxL=mfxR8 - mfyL=mfyR8 - cxL= cxR8 - cyL= cyR8 -#endif -! if ( flagstruct%range_warn ) then -! call range_check('CX_dyn', cxL(is:ie,js:je,:)/real(n_split), is, ie, js, je, 0, npz, gridstruct%agrid, & -! -0.5, 0.5, bad_range) -! call range_check('CY_dyn', cyL(is:ie,js:je,:)/real(n_split), is, ie, js, je, 0, npz, gridstruct%agrid, & -! -0.5, 0.5, bad_range) -! endif - if ( flagstruct%fv_debug ) then if(is_master()) write(6,*) '' if(is_master()) write(6,*) '------------------------------------' @@ -688,24 +668,25 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, & !DryMassRoundoffControl if(last_step) then if (hydrostatic) then - do k = 2, npz - dpx(:,:,1) = dpx(:,:,1) + dpx(:,:,k) - enddo - do j=js,je - do i=is,ie - psx(i,j) = psx(i,j) + dpx(i,j,1) - enddo - enddo - call timing_on('COMM_TOTAL') - call mpp_update_domains(psx, domain) - call timing_off('COMM_TOTAL') - do j=js-1,je+1 - do i=is-1,ie+1 - pe(i,npz+1,j) = psx(i,j) - enddo - enddo + do k = 2, npz + dpx(:,:,1) = dpx(:,:,1) + dpx(:,:,k) + enddo + do j=js,je + do i=is,ie + psx(i,j) = psx(i,j) + dpx(i,j,1) + enddo + enddo + call timing_on('COMM_TOTAL') + call mpp_update_domains(psx, domain) + call timing_off('COMM_TOTAL') + do j=js-1,je+1 + do i=is-1,ie+1 + pe(i,npz+1,j) = psx(i,j) + enddo + enddo + deallocate(psx) endif - deallocate(psx,dpx) + deallocate(dpx) end if #ifdef SW_DYNAMICS @@ -723,18 +704,18 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, & call timing_on('tracer_2d') !!! CLEANUP: merge these two calls? if (gridstruct%nested) then - call tracer_2d_nested(q, dp1, mfxL, mfyL, cxL, cyL, gridstruct, bd, domain, npx, npy, npz, nq, & + call tracer_2d_nested(q, dp1, mfxR8, mfyR8, cxR8, cyR8, cmax_z, imax_req, gridstruct, bd, domain, npx, npy, npz, nq, & flagstruct%hord_tr, q_split, mdt, idiag%id_divg, i_pack(10), & flagstruct%nord_tr, flagstruct%trdm2, & k_split, neststruct, parent_grid, flagstruct%lim_fac) else if ( flagstruct%z_tracer .and. q_split == 0 ) then - call tracer_2d_1L(q, dp1, mfxL, mfyL, cxL, cyL, gridstruct, bd, domain, npx, npy, npz, nq, & + call tracer_2d_1L(q, dp1, mfxR8, mfyR8, cxR8, cyR8, cmax_z, imax_req, gridstruct, bd, domain, npx, npy, npz, nq, & flagstruct%hord_tr, mdt, idiag%id_divg, i_pack(10), & flagstruct%nord_tr, flagstruct%trdm2, flagstruct%lim_fac) else - call tracer_2d(q, dp1, mfxL, mfyL, cxL, cyL, gridstruct, bd, domain, npx, npy, npz, nq, & - flagstruct%hord_tr, q_split, mdt, idiag%id_divg, i_pack(10), & + call tracer_2d(q, dp1, mfxR8, mfyR8, cxR8, cyR8, cmax_z, imax_req, gridstruct, bd, domain, npx, npy, npz, nq, & + flagstruct%hord_tr, flagstruct%n_sponge, q_split, mdt, idiag%id_divg, i_pack(10), & flagstruct%nord_tr, flagstruct%trdm2, flagstruct%lim_fac) endif endif @@ -793,7 +774,6 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, & flagstruct%adiabatic, do_adiabatic_init, & flagstruct%remap_option, flagstruct%gmao_remap, & flagstruct%gmao_top_bc, flagstruct%gmao_bot_bc) -!!! mfx=mfxL, mfy=mfyL, cx=cxL, cy=cyL) #ifdef AVEC_TIMERS call avec_timer_stop(6) @@ -808,10 +788,10 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, & #endif ! Accumulate the total Mass flux and Courant numbers for export - mfx = mfx + mfxL - mfy = mfy + mfyL - cx = cx + cxL - cy = cy + cyL + mfx = mfx + mfxR8 + mfy = mfy + mfyR8 + cx = cx + cxR8 + cy = cy + cyR8 if ( flagstruct%fv_debug ) then if(is_master()) write(6,*) '' @@ -1053,15 +1033,11 @@ subroutine Rayleigh_fast(dt, npx, npy, npz, pfull, tau, u, v, w, & do k=1, npz if ( pfull(k) < rf_cutoff ) then rff(k) = dt/tau0*sin(0.5*pi*log(rf_cutoff/pfull(k))/log(rf_cutoff/ptop))**2 +! Re-FACTOR rf + if( is_master() ) write(6,*) k, 0.01*pfull(k), dt/(rff(k)*sday) kmax = k rff(k) = 1.d0 / (1.0d0+rff(k)) rf(k) = rff(k) - ! Added safety check to prevent a log(1.0) divide-by-zero crash - if( is_master() ) then - if (rff(k) < 1.0d0) then - write(6,*) k, 0.01*pfull(k), (-dt / log(rff(k))) / sday - endif - endif else exit endif diff --git a/model/fv_mapz.F90 b/model/fv_mapz.F90 index 6ca2c4059..9e75262de 100644 --- a/model/fv_mapz.F90 +++ b/model/fv_mapz.F90 @@ -1543,6 +1543,177 @@ subroutine map_scalar( km, pe1, q1, & logical, optional, intent(in):: optional_bot !< optional GMAO bottom BC logical, optional, intent(in):: optional_top !< optional GMAO top BC +! DESCRIPTION: +! IV = 0: constituents +! pe1: pressure at layer edges (from model top to bottom surface) +! in the original vertical coordinate +! pe2: pressure at layer edges (from model top to bottom surface) +! in the new vertical coordinate +! LOCAL VARIABLES: + real, allocatable :: q4(:,:,:) + real qsum, pl, pr, pfac0, pfac1, pfac2, dp, esl + integer i, k, l, m, k0 + integer LM1,LP0,LP1 + logical gmao_bot, gmao_top + integer LP0_start ! OPTIMIZATION: State tracking variable + + ! Assuming these constants from context (1/3 and 2/3) + real, parameter :: r3 = 1.0/3.0 + real, parameter :: r23 = 2.0/3.0 + + gmao_bot=.false. + if (present(optional_bot)) gmao_bot=optional_bot + + gmao_top=.false. + if (present(optional_top)) gmao_top=optional_top + + allocate ( q4(4,i1:i2,km) ) + + do k=1,km + do i=i1,i2 + q4(1,i,k) = q1(i,j,k) + enddo + enddo + if ( kord > 7 ) then + if (present(q_min)) then + call scalar_profile( q4(1,i1,1), dpe1, km, i1, i2, iv, kord, q_min, qs=q_bot) + else + call scalar_profile( q4(1,i1,1), dpe1, km, i1, i2, iv, kord, no_min, qs=q_bot) + endif + else + call ppm_profile( q4(1,i1,1), dpe1, km, i1, i2, iv, kord) + endif + +! Interpolate field onto target Pressures +! --------------------------------------- + do i=i1,i2 + k0 = 1 + LP0_start = 1 ! OPTIMIZATION: Track the starting pressure index + + k_loop: do k=1,kn + LM1 = 1 + LP0 = LP0_start ! Start from the previous k's index instead of 1 + + do while( LP0.le.km ) + if (pe1(i,LP0).lt.pe2(i,k)) then + LP0 = LP0+1 + else + exit + endif + enddo + + LP0_start = LP0 ! Save the state for the next vertical level (k+1) + + LM1 = max(LP0-1,1) + LP0 = min(LP0, km) + +! Entire grid below old ps +! ---------------------------------------------------- + if( gmao_bot .and. (pe2(i,k) .ge. pe1(i,km+1)) ) then + q2(i,k) = q1(i,j,km) +! Extrapolate Linearly below last model level +! ---------------------------------------------------- + else if( gmao_bot .and. (LM1.eq.km .and. LP0.eq.km) ) then + pfac0=(pe2(i,k )-pe1(i,km ))/(pe1(i,km)-pe1(i,km-1)) + q2(i,k) = q1(i,j,km) + ( q1(i,j,km)-q1(i,j,km-1) )*pfac0 +! Interpolate Linearly between levels 1 => 2 and km-1 => km +! ----------------------------------------------------------------- + else if( gmao_bot .and. (LP0.eq.km) ) then + pfac0=(pe2(i,k )-pe1(i,LP0))/(pe1(i,LM1)-pe1(i,LP0)) + q2(i,k) = q1(i,j,LP0) + ( q1(i,j,LM1)-q1(i,j,LP0) )*pfac0 +! Extrapolate Linearly above first model level +! ---------------------------------------------------- + else if( gmao_top .and. (LM1.eq.1 .and. LP0.eq.1) ) then + pfac0=(pe2(i,k)-pe1(i,1))/(pe1(i,2)-pe1(i,1)) + q2(i,k) = q1(i,j,1) + ( q1(i,j,2)-q1(i,j,1) )*pfac0 +! Interpolate Linearly between levels 1 => 2 and km-1 => km +! ----------------------------------------------------------------- + else if( gmao_top .and. (LM1.eq.1) ) then + pfac0=(pe2(i,k )-pe1(i,LP0))/(pe1(i,LM1)-pe1(i,LP0)) + q2(i,k) = q1(i,j,LP0) + ( q1(i,j,LM1)-q1(i,j,LP0) )*pfac0 +! Use standard FV3 remapping in between +! ---------------------------------------------------- + else + l_loop: do l=k0,km + ! locate the top edge: pe2(i,k) + if( pe2(i,k) >= pe1(i,l) .and. pe2(i,k) <= pe1(i,l+1) ) then + pl = (pe2(i,k)-pe1(i,l)) / dpe1(i,l) + if( pe2(i,k+1) <= pe1(i,l+1) ) then + ! entire new grid is within the original grid + pr = (pe2(i,k+1)-pe1(i,l)) / dpe1(i,l) + pfac1 = (pr+pl) + pfac2 = r3*(pr*pfac1+pl**2) + q2(i,k) = q4(2,i,l) + 0.5*(q4(4,i,l)+q4(3,i,l)-q4(2,i,l)) & + *pfac1-q4(4,i,l)*pfac2 + k0 = l + cycle k_loop ! Replaces goto 555 + else + ! Fractional area... + pfac0 = pe1(i,l+1)-pe2(i,k) + pfac1 = (1.+pl) + pfac2 = r3*(1.+pl*pfac1) + qsum = pfac0*(q4(2,i,l)+0.5*(q4(4,i,l)+ & + q4(3,i,l)-q4(2,i,l))*pfac1-q4(4,i,l)*pfac2) + + m_loop: do m=l+1,km + ! locate the bottom edge: pe2(i,k+1) + if( pe2(i,k+1) > pe1(i,m+1) ) then + ! Whole layer + qsum = qsum + dpe1(i,m)*q4(1,i,m) + else + dp = pe2(i,k+1)-pe1(i,m) + esl = (dp / dpe1(i,m)) + pfac0 = (1.-r23*esl) + qsum = qsum + dp*(q4(2,i,m)+0.5*esl* & + (q4(3,i,m)-q4(2,i,m)+q4(4,i,m)*pfac0)) + k0 = m + exit l_loop ! Replaces goto 123 + endif + enddo m_loop + + exit l_loop ! Replaces goto 123 + endif + endif + enddo l_loop + + ! Equivalent to 123 continue + q2(i,k) = qsum / dpe2(i,k) + + endif + enddo k_loop + enddo + + deallocate( q4 ) + + end subroutine map_scalar + + subroutine map_scalar_old( km, pe1, q1, & + kn, pe2, q2, & + dpe1, dpe2, i1, i2, & + j, ibeg, iend, jbeg, jend, iv, kord, q_min, & + q_bot, optional_bot, optional_top) +! iv=1 + integer, intent(in) :: i1 !< Starting longitude + integer, intent(in) :: i2 !< Finishing longitude + integer, intent(in) :: iv !< Mode: 0 == constituents 1 == temp 2 == remap temp with cs scheme + integer, intent(in) :: kord !< Method order + integer, intent(in) :: j !< Current latitude + integer, intent(in) :: ibeg, iend, jbeg, jend + integer, intent(in) :: km !< Original vertical dimension + integer, intent(in) :: kn !< Target vertical dimension + real, intent(in) :: pe1(i1:i2,km+1) !< pressure at layer edges from model top to bottom surface in the original vertical coordinate + real, intent(in) :: pe2(i1:i2,kn+1) !< pressure at layer edges from model top to bottom surface in the new vertical coordinate + real, intent(in) :: dpe1(i1:i2,km) !< pressure thickness in the original vertical coordinate + real, intent(in) :: dpe2(i1:i2,kn) !< pressure thickness in the new vertical coordinate + real, intent(in) :: q1(ibeg:iend,jbeg:jend,1:km) !< Field input +! INPUT/OUTPUT PARAMETERS: + real, intent(inout):: q2(i1:i2,1:kn) !< Field output +! Optional aruguments: + real, optional, intent(in):: q_min !< minimum for scheme + real, optional, intent(in):: q_bot(i1:i2) !< bottom BC + logical, optional, intent(in):: optional_bot !< optional GMAO bottom BC + logical, optional, intent(in):: optional_top !< optional GMAO top BC + ! DESCRIPTION: ! IV = 0: constituents ! pe1: pressure at layer edges (from model top to bottom surface) @@ -1668,7 +1839,7 @@ subroutine map_scalar( km, pe1, q1, & deallocate( q4 ) - end subroutine map_scalar + end subroutine map_scalar_old subroutine mapn_tracer(nq, km, pe1, pe2, q1, dp2, kord, j, & i1, i2, isd, ied, jsd, jed, q_min, fill, qs) @@ -1709,6 +1880,124 @@ subroutine mapn_tracer(nq, km, pe1, pe2, q1, dp2, kord, j, & call scalar_profile( q4(1,i1,1,iq), dp1, km, i1, i2, 0, kord(iq), q_min, qs=qs) enddo +! Mapping +! OPTIMIZATION: Replaced GOTOs with modern named loop exits/cycles for vectorization + do i=i1,i2 + k0 = 1 + + k_loop: do k=1,km + l_loop: do l=k0,km + ! locate the top edge: pe2(i,k) + if(pe2(i,k) >= pe1(i,l) .and. pe2(i,k) <= pe1(i,l+1)) then + pl = (pe2(i,k)-pe1(i,l)) / dp1(i,l) + + if(pe2(i,k+1) <= pe1(i,l+1)) then + ! entire new grid is within the original grid + pr = (pe2(i,k+1)-pe1(i,l)) / dp1(i,l) + fac1 = pr + pl + fac2 = r3*(pr*fac1 + pl*pl) + fac1 = 0.5*fac1 + do iq=1,nq + q2(i,k,iq) = q4(2,i,l,iq) + (q4(4,i,l,iq)+q4(3,i,l,iq)-q4(2,i,l,iq))*fac1 & + - q4(4,i,l,iq)*fac2 + enddo + k0 = l + cycle k_loop ! Replaces goto 555 + else + ! Fractional area... + dp = pe1(i,l+1) - pe2(i,k) + fac1 = 1. + pl + fac2 = r3*(1.+pl*fac1) + fac1 = 0.5*fac1 + do iq=1,nq + qsum(iq) = dp*(q4(2,i,l,iq) + (q4(4,i,l,iq)+ & + q4(3,i,l,iq) - q4(2,i,l,iq))*fac1 - q4(4,i,l,iq)*fac2) + enddo + + m_loop: do m=l+1,km + ! locate the bottom edge: pe2(i,k+1) + if(pe2(i,k+1) > pe1(i,m+1) ) then + ! Whole layer.. + do iq=1,nq + qsum(iq) = qsum(iq) + dp1(i,m)*q4(1,i,m,iq) + enddo + else + dp = pe2(i,k+1)-pe1(i,m) + esl = dp / dp1(i,m) + fac1 = 0.5*esl + fac2 = 1.-r23*esl + do iq=1,nq + qsum(iq) = qsum(iq) + dp*( q4(2,i,m,iq) + fac1*( & + q4(3,i,m,iq)-q4(2,i,m,iq)+q4(4,i,m,iq)*fac2 ) ) + enddo + k0 = m + exit l_loop ! Replaces goto 123 + endif + enddo m_loop + + exit l_loop ! Replaces goto 123 (if m_loop finished without exit) + endif + endif + enddo l_loop + + ! 123 continue equivalent + do iq=1,nq + q2(i,k,iq) = qsum(iq) / dp2(i,k) + enddo + enddo k_loop + enddo + + if (fill) call fillz(i2-i1+1, km, nq, q2, dp2) + + do iq=1,nq + do k=1,km + do i=i1,i2 + q1(i,j,k,iq) = q2(i,k,iq) + enddo + enddo + enddo + + end subroutine mapn_tracer + + subroutine mapn_tracer_old(nq, km, pe1, pe2, q1, dp2, kord, j, & + i1, i2, isd, ied, jsd, jed, q_min, fill, qs) +! INPUT PARAMETERS: + integer, intent(in):: km !< vertical dimension + integer, intent(in):: j, nq, i1, i2 + integer, intent(in):: isd, ied, jsd, jed + integer, intent(in):: kord(nq) + real, intent(in):: pe1(i1:i2,km+1) !< pressure at layer edges from model top to bottom surface in the original vertical coordinate + real, intent(in):: pe2(i1:i2,km+1) !< pressure at layer edges from model top to bottom surface in the new vertical coordinate + real, intent(in):: dp2(i1:i2,km) + real, intent(in):: q_min + logical, intent(in):: fill + real, intent(inout):: q1(isd:ied,jsd:jed,km,nq) ! Field input + + real, optional, intent(in) :: qs(i1:i2) + +! LOCAL VARIABLES: + real:: q4(4,i1:i2,km,nq) + real:: q2(i1:i2,km,nq) !< Field output + real:: qsum(nq) + real:: dp1(i1:i2,km) + real:: pl, pr, dp, esl, fac1, fac2 + integer:: i, k, l, m, k0, iq + + do k=1,km + do i=i1,i2 + dp1(i,k) = pe1(i,k+1) - pe1(i,k) + enddo + enddo + + do iq=1,nq + do k=1,km + do i=i1,i2 + q4(1,i,k,iq) = q1(i,j,k,iq) + enddo + enddo + call scalar_profile( q4(1,i1,1,iq), dp1, km, i1, i2, 0, kord(iq), q_min, qs=qs) + enddo + ! Mapping do 4000 i=i1,i2 k0 = 1 @@ -1781,7 +2070,7 @@ subroutine mapn_tracer(nq, km, pe1, pe2, q1, dp2, kord, j, & enddo enddo - end subroutine mapn_tracer + end subroutine mapn_tracer_old subroutine remap_2d(km, pe1, q1, & kn, pe2, q2, & @@ -1876,9 +2165,434 @@ subroutine remap_2d(km, pe1, q1, & end subroutine remap_2d + subroutine scalar_profile(a4, delp, km, i1, i2, iv, kord, qmin, qs) +! Optimized vertical profile reconstruction: +! Latest: Apr 2008 S.-J. Lin, NOAA/GFDL + integer, intent(in):: i1, i2 + integer, intent(in):: km !< vertical dimension + integer, intent(in):: iv !< iv =-1: winds iv = 0: positive definite scalars iv = 1: others + integer, intent(in):: kord + real, intent(in) :: delp(i1:i2,km) !< Layer pressure thickness + real, intent(inout):: a4(4,i1:i2,km) !< Interpolated values + real, intent(in):: qmin + real, optional, intent(in) :: qs(i1:i2) +!----------------------------------------------------------------------- + logical, dimension(i1:i2,km):: extm, ext5, ext6 + real gam(i1:i2,km) + real q(i1:i2,km+1) + real d4(i1:i2) + real bet, a_bot, grat + real pmp_1, lac_1, pmp_2, lac_2, x0, x1 + integer i, k, im + + if ( iv .eq. -2 ) then + if (.not. present(qs)) call mpp_error (FATAL, 'fv_mapz::scalar_profile - qs is not present') + do i=i1,i2 + gam(i,2) = 0.5 + q(i,1) = 1.5*a4(1,i,1) + enddo + do k=2,km-1 + do i=i1, i2 + grat = delp(i,k-1) / delp(i,k) + bet = 2. + grat + grat - gam(i,k) + q(i,k) = (3.*(a4(1,i,k-1)+a4(1,i,k)) - q(i,k-1))/bet + gam(i,k+1) = grat / bet + enddo + enddo + do i=i1,i2 + grat = delp(i,km-1) / delp(i,km) + q(i,km) = (3.*(a4(1,i,km-1)+a4(1,i,km)) - grat*qs(i) - q(i,km-1)) / & + (2. + grat + grat - gam(i,km)) + q(i,km+1) = qs(i) + enddo + do k=km-1,1,-1 + do i=i1,i2 + q(i,k) = q(i,k) - gam(i,k+1)*q(i,k+1) + enddo + enddo + else + do i=i1,i2 + grat = delp(i,2) / delp(i,1) ! grid ratio + bet = grat*(grat+0.5) + q(i,1) = ( (grat+grat)*(grat+1.)*a4(1,i,1) + a4(1,i,2) ) / bet + gam(i,1) = ( 1. + grat*(grat+1.5) ) / bet + enddo + + do k=2,km + do i=i1,i2 + d4(i) = delp(i,k-1) / delp(i,k) + bet = 2. + d4(i) + d4(i) - gam(i,k-1) + q(i,k) = ( 3.*(a4(1,i,k-1)+d4(i)*a4(1,i,k)) - q(i,k-1) )/bet + gam(i,k) = d4(i) / bet + enddo + enddo + + do i=i1,i2 + a_bot = 1. + d4(i)*(d4(i)+1.5) + q(i,km+1) = (2.*d4(i)*(d4(i)+1.)*a4(1,i,km)+a4(1,i,km-1)-a_bot*q(i,km)) & + / ( d4(i)*(d4(i)+0.5) - a_bot*gam(i,km) ) + enddo + + do k=km,1,-1 + do i=i1,i2 + q(i,k) = q(i,k) - gam(i,k)*q(i,k+1) + enddo + enddo + endif + +!----- Perfectly linear scheme -------------------------------- + if ( abs(kord) > 16 ) then + do k=1,km + do i=i1,i2 + a4(2,i,k) = q(i,k ) + a4(3,i,k) = q(i,k+1) + a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))) + enddo + enddo + return + endif +!----- Perfectly linear scheme -------------------------------- + +!------------------ +! Apply constraints +!------------------ + im = i2 - i1 + 1 + +! Apply *large-scale* constraints + do i=i1,i2 + q(i,2) = min( q(i,2), max(a4(1,i,1), a4(1,i,2)) ) + q(i,2) = max( q(i,2), min(a4(1,i,1), a4(1,i,2)) ) + enddo + + do k=2,km + do i=i1,i2 + gam(i,k) = a4(1,i,k) - a4(1,i,k-1) + enddo + enddo + +! Interior: + do k=3,km-1 + do i=i1,i2 + if ( gam(i,k-1)*gam(i,k+1)>0. ) then +! Apply large-scale constraint to ALL fields if not local max/min + q(i,k) = min( q(i,k), max(a4(1,i,k-1),a4(1,i,k)) ) + q(i,k) = max( q(i,k), min(a4(1,i,k-1),a4(1,i,k)) ) + else + if ( gam(i,k-1) > 0. ) then +! There exists a local max + q(i,k) = max(q(i,k), min(a4(1,i,k-1),a4(1,i,k))) + else +! There exists a local min + q(i,k) = min(q(i,k), max(a4(1,i,k-1),a4(1,i,k))) + if ( iv==0 ) q(i,k) = max(0., q(i,k)) + endif + endif + enddo + enddo + +! Bottom: + do i=i1,i2 + q(i,km) = min( q(i,km), max(a4(1,i,km-1), a4(1,i,km)) ) + q(i,km) = max( q(i,km), min(a4(1,i,km-1), a4(1,i,km)) ) + enddo + + do k=1,km + do i=i1,i2 + a4(2,i,k) = q(i,k ) + a4(3,i,k) = q(i,k+1) + enddo + enddo + + do k=1,km + if ( k==1 .or. k==km ) then + do i=i1,i2 + extm(i,k) = (a4(2,i,k)-a4(1,i,k)) * (a4(3,i,k)-a4(1,i,k)) > 0. + enddo + else + do i=i1,i2 + extm(i,k) = gam(i,k)*gam(i,k+1) < 0. + enddo + endif + if ( abs(kord) > 9 ) then + do i=i1,i2 + x0 = 2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)) + x1 = abs(a4(2,i,k)-a4(3,i,k)) + a4(4,i,k) = 3.*x0 + ext5(i,k) = abs(x0) > x1 + ext6(i,k) = abs(a4(4,i,k)) > x1 + enddo + endif + enddo + +!--------------------------- +! Apply subgrid constraints: +!--------------------------- +! f(s) = AL + s*[(AR-AL) + A6*(1-s)] ( 0 <= s <= 1 ) +! Top 2 and bottom 2 layers always use monotonic mapping + + if ( iv==0 ) then + do i=i1,i2 + a4(2,i,1) = max(0., a4(2,i,1)) + enddo + elseif ( iv==-1 ) then + do i=i1,i2 + if ( a4(2,i,1)*a4(1,i,1) <= 0. ) a4(2,i,1) = 0. + enddo + elseif ( iv==2 ) then + do i=i1,i2 + a4(2,i,1) = a4(1,i,1) + a4(3,i,1) = a4(1,i,1) + a4(4,i,1) = 0. + enddo + endif + + if ( iv/=2 ) then + do i=i1,i2 + a4(4,i,1) = 3.*(2.*a4(1,i,1) - (a4(2,i,1)+a4(3,i,1))) + enddo + call cs_limiters(im, extm(i1,1), a4(1,i1,1), 1) + endif + +! k=2 + do i=i1,i2 + a4(4,i,2) = 3.*(2.*a4(1,i,2) - (a4(2,i,2)+a4(3,i,2))) + enddo + call cs_limiters(im, extm(i1,2), a4(1,i1,2), 2) + + +!------------------------------------- +! Huynh's 2nd constraint for interior: +!------------------------------------- +! OPTIMIZATION: Hoisted kord select block outside the vertical loop. Loop Fusion applied where possible. + select case(abs(kord)) + + case (0:8) + do k=3,km-2 + do i=i1,i2 + ! Left edges + pmp_1 = a4(1,i,k) - 2.*gam(i,k+1) + lac_1 = pmp_1 + 1.5*gam(i,k+2) + a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), & + max(a4(1,i,k), pmp_1, lac_1) ) + ! Right edges + pmp_2 = a4(1,i,k) + 2.*gam(i,k) + lac_2 = pmp_2 - 1.5*gam(i,k-1) + a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), & + max(a4(1,i,k), pmp_2, lac_2) ) + ! Fused + a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))) + enddo + enddo + + case (9) + do k=3,km-2 + do i=i1,i2 + if ( extm(i,k) .and. extm(i,k-1) ) then + a4(2,i,k) = a4(1,i,k) + a4(3,i,k) = a4(1,i,k) + a4(4,i,k) = 0. + else if ( extm(i,k) .and. extm(i,k+1) ) then + a4(2,i,k) = a4(1,i,k) + a4(3,i,k) = a4(1,i,k) + a4(4,i,k) = 0. + else if ( extm(i,k) .and. a4(1,i,k) abs(a4(2,i,k)-a4(3,i,k)) ) then + pmp_1 = a4(1,i,k) - 2.*gam(i,k+1) + lac_1 = pmp_1 + 1.5*gam(i,k+2) + a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), & + max(a4(1,i,k), pmp_1, lac_1) ) + pmp_2 = a4(1,i,k) + 2.*gam(i,k) + lac_2 = pmp_2 - 1.5*gam(i,k-1) + a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), & + max(a4(1,i,k), pmp_2, lac_2) ) + a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))) + endif + endif + enddo + enddo + + case (10) + do k=3,km-2 + do i=i1,i2 + if( ext5(i,k) ) then + if( ext5(i,k-1) .or. ext5(i,k+1) ) then + a4(2,i,k) = a4(1,i,k) + a4(3,i,k) = a4(1,i,k) + elseif ( ext6(i,k-1) .or. ext6(i,k+1) ) then + pmp_1 = a4(1,i,k) - 2.*gam(i,k+1) + lac_1 = pmp_1 + 1.5*gam(i,k+2) + a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), & + max(a4(1,i,k), pmp_1, lac_1) ) + pmp_2 = a4(1,i,k) + 2.*gam(i,k) + lac_2 = pmp_2 - 1.5*gam(i,k-1) + a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), & + max(a4(1,i,k), pmp_2, lac_2) ) + endif + elseif( ext6(i,k) ) then + if( ext5(i,k-1) .or. ext5(i,k+1) ) then + pmp_1 = a4(1,i,k) - 2.*gam(i,k+1) + lac_1 = pmp_1 + 1.5*gam(i,k+2) + a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), & + max(a4(1,i,k), pmp_1, lac_1) ) + pmp_2 = a4(1,i,k) + 2.*gam(i,k) + lac_2 = pmp_2 - 1.5*gam(i,k-1) + a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), & + max(a4(1,i,k), pmp_2, lac_2) ) + endif + endif + ! Fused + a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))) + enddo + enddo + + case (12) + do k=3,km-2 + do i=i1,i2 + if( extm(i,k) ) then + a4(2,i,k) = a4(1,i,k) + a4(3,i,k) = a4(1,i,k) + a4(4,i,k) = 0. + else + a4(4,i,k) = 6.*a4(1,i,k) - 3.*(a4(2,i,k)+a4(3,i,k)) + if( abs(a4(4,i,k)) > abs(a4(2,i,k)-a4(3,i,k)) ) then + pmp_1 = a4(1,i,k) - 2.*gam(i,k+1) + lac_1 = pmp_1 + 1.5*gam(i,k+2) + a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), & + max(a4(1,i,k), pmp_1, lac_1) ) + pmp_2 = a4(1,i,k) + 2.*gam(i,k) + lac_2 = pmp_2 - 1.5*gam(i,k-1) + a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), & + max(a4(1,i,k), pmp_2, lac_2) ) + a4(4,i,k) = 6.*a4(1,i,k) - 3.*(a4(2,i,k)+a4(3,i,k)) + endif + endif + enddo + enddo + + case (13) + do k=3,km-2 + do i=i1,i2 + if( ext6(i,k) ) then + if ( ext6(i,k-1) .and. ext6(i,k+1) ) then + a4(2,i,k) = a4(1,i,k) + a4(3,i,k) = a4(1,i,k) + endif + endif + ! Fused + a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))) + enddo + enddo + + case (14) + do k=3,km-2 + do i=i1,i2 + a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))) + enddo + enddo + + case (15) + do k=3,km-2 + do i=i1,i2 + if ( ext5(i,k) .and. ext5(i,k-1) ) then + a4(2,i,k) = a4(1,i,k) + a4(3,i,k) = a4(1,i,k) + else if ( ext5(i,k) .and. ext5(i,k+1) ) then + a4(2,i,k) = a4(1,i,k) + a4(3,i,k) = a4(1,i,k) + else if ( ext5(i,k) .and. a4(1,i,k)@brief Optimized vertical profile reconstruction: !> Latest: Apr 2008 S.-J. Lin, NOAA/GFDL - subroutine scalar_profile(a4, delp, km, i1, i2, iv, kord, qmin, qs) + subroutine scalar_profile_old(a4, delp, km, i1, i2, iv, kord, qmin, qs) ! Optimized vertical profile reconstruction: ! Latest: Apr 2008 S.-J. Lin, NOAA/GFDL integer, intent(in):: i1, i2 @@ -2285,7 +2999,7 @@ subroutine scalar_profile(a4, delp, km, i1, i2, iv, kord, qmin, qs) if(k== km ) call cs_limiters(im, extm(i1,k), a4(1,i1,k), 1) enddo - end subroutine scalar_profile + end subroutine scalar_profile_old !>@brief The subroutine 'cs_profile' performs the optimized vertical profile reconstruction: !>@date April 2008 diff --git a/model/fv_tracer2d.F90 b/model/fv_tracer2d.F90 index c807a98be..78436cc49 100644 --- a/model/fv_tracer2d.F90 +++ b/model/fv_tracer2d.F90 @@ -63,7 +63,8 @@ module fv_tracer2d_mod use tp_core_mod, only: fv_tp_2d, copy_corners - use fv_mp_mod, only: mp_reduce_max + use fv_mp_mod, only: mp_barrier, mp_reduce_max, mp_reduce_min + use fv_mp_mod, only: mp_ireduce_max, mp_wait_reduce_max use fv_mp_mod, only: ng, mp_gather, is_master use fv_mp_mod, only: group_halo_update_type use fv_mp_mod, only: start_group_halo_update, complete_group_halo_update @@ -76,6 +77,9 @@ module fv_tracer2d_mod use fv_grid_utils_mod, only: g_sum_r8 implicit none + +#include "mpif.h" + private public :: tracer_2d, tracer_2d_nested, tracer_2d_1L, offline_tracer_advection @@ -89,7 +93,7 @@ module fv_tracer2d_mod !! It modifies 'tracer_2d' so that each layer uses a different diagnosed number !! of split tracer timesteps. This potentially accelerates tracer advection when there !! is a large difference in layer-maximum wind speeds (cf. polar night jet). -subroutine tracer_2d_1L(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, npz, & +subroutine tracer_2d_1L(q, dp1, mfx, mfy, cx, cy, cmax, imax_req, gridstruct, bd, domain, npx, npy, npz, & nq, hord, dt, id_divg, q_pack, nord_tr, trdm, lim_fac, dpA) type(fv_grid_bounds_type), intent(IN) :: bd @@ -104,10 +108,15 @@ subroutine tracer_2d_1L(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, n type(group_halo_update_type), intent(inout) :: q_pack real , intent(INOUT) :: q(bd%isd:bd%ied,bd%jsd:bd%jed,npz,nq) !< Tracers real , intent(INOUT) :: dp1(bd%isd:bd%ied,bd%jsd:bd%jed,npz) !< DELP before dyn_core - real , intent(IN ) :: mfx(bd%is:bd%ie+1,bd%js:bd%je, npz) !< Mass Flux X-Dir - real , intent(IN ) :: mfy(bd%is:bd%ie ,bd%js:bd%je+1,npz) !< Mass Flux Y-Dir - real , intent(IN ) :: cx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz) !< Courant Number X-Dir - real , intent(IN ) :: cy(bd%isd:bd%ied,bd%js :bd%je +1,npz) !< Courant Number Y-Dir + + real(kind=8) , intent(IN ) :: mfx(bd%is:bd%ie+1,bd%js:bd%je, npz) !< Mass Flux X-Dir + real(kind=8) , intent(IN ) :: mfy(bd%is:bd%ie ,bd%js:bd%je+1,npz) !< Mass Flux Y-Dir + real(kind=8) , intent(IN ) :: cx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz) !< Courant Number X-Dir + real(kind=8) , intent(IN ) :: cy(bd%isd:bd%ied,bd%js :bd%je +1,npz) !< Courant Number Y-Dir + + real , intent(INOUT) :: cmax(npz) + integer, intent(INOUT) :: imax_req + real , optional, intent(OUT) :: dpA(bd%is:bd%ie,bd%js:bd%je,npz) ! DELP after advection type(fv_grid_type), intent(IN), target :: gridstruct type(domain2d), intent(INOUT) :: domain @@ -126,8 +135,6 @@ subroutine tracer_2d_1L(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, n real :: cx2(bd%is:bd%ie+1,bd%jsd:bd%jed, npz) real :: cy2(bd%isd:bd%ied,bd%js :bd%je +1, npz) real :: c2d(bd%is:bd%ie,bd%js:bd%je) - real :: cmax(npz), cmin(npz) - real :: qmax(npz*(nq+1)), qmin(npz*(nq+1)) integer :: icount(npz,nq) real :: frac integer :: nsplt @@ -137,6 +144,10 @@ subroutine tracer_2d_1L(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, n real, pointer, dimension(:,:,:) :: sin_sg real, pointer, dimension(:,:) :: dxa, dya, dx, dy + real, parameter :: TRACER_EPS = tiny(1.0) ! Adjust threshold as needed + + real(kind=8) :: t_start, t_wait, t_min_wait, t_max_wait + integer :: is, ie, js, je integer :: isd, ied, jsd, jed @@ -161,6 +172,7 @@ subroutine tracer_2d_1L(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, n !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx,dxa,dy, & !$OMP sin_sg,cy,yfx,dya,dx) private(i,j,k) do k=1,npz + ! 1. Compute xfx do j=jsd,jed do i=is,ie+1 if (cx(i,j,k) > 0.) then @@ -170,6 +182,8 @@ subroutine tracer_2d_1L(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, n endif enddo enddo + + ! 2. Compute yfx do j=js,je+1 do i=isd,ied if (cy(i,j,k) > 0.) then @@ -181,58 +195,22 @@ subroutine tracer_2d_1L(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, n enddo enddo ! k-loop + call timing_on('COMM_TOTAL') + call timing_on('COMM_TRACER_WAIT') + call mp_wait_reduce_max(cmax, npz, imax_req) + call timing_off('COMM_TRACER_WAIT') + call timing_off('COMM_TOTAL') + call timing_on('COMM_TOTAL') call timing_on('COMM_TRACER') call complete_group_halo_update(q_pack, domain) call timing_off('COMM_TRACER') call timing_off('COMM_TOTAL') -! Check for levels where Q does not need to be advected -!$OMP parallel do default(none) shared(nq,npz,is,ie,js,je,q,qmax) private(iq,i,j,k,n) - do iq=1,nq - do k=1,npz - n=(iq-1)*npz + k - qmax(n) = 0.0 - do j=js,je - do i=is,ie - qmax(n) = max(qmax(n),q(i,j,k,iq)) - enddo - enddo - enddo - enddo -!$OMP parallel do default(none) shared(is,ie,js,je,npz,nq,& -!$OMP sin_sg,cx,cy,cmax,qmax) private(i,j,k,n) - do k=1,npz - cmax(k) = 0. - if ( k < npz/6 ) then - do j=js,je - do i=is,ie - cmax(k) = max( cmax(k), abs(cx(i,j,k)), abs(cy(i,j,k)) ) - enddo - enddo - else - do j=js,je - do i=is,ie - cmax(k) = max( cmax(k), max(abs(cx(i,j,k)),abs(cy(i,j,k)))+1.-sin_sg(i,j,5) ) - enddo - enddo - endif - !!! if ( is_master() ) write(*,*) 'tracer_2d_1L: k, nsplt, cmax =', k, int(1. + cmax(k)), cmax(k) - ! add to qmax for allreduce - n=nq*npz + k - qmax(n) = cmax(k) - enddo ! k-loop - - call timing_on('COMM_TOTAL') - call timing_on('COMM_TRACER_MAX') - call mp_reduce_max(qmax,(nq+1)*npz) - call timing_off('COMM_TRACER_MAX') - call timing_off('COMM_TOTAL') - ! get cmax from allreduce array + ! check cmax from allreduce array do k=1,npz - n=nq*npz + k - cmax(k) = qmax(n) + if ( is_master() .and. (cmax(k) > 1.0) ) write(*,*) 'tracer_2d_1L: k, qsplt =', k, int(1. + cmax(k)) enddo ! k-loop !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx, & @@ -362,8 +340,8 @@ subroutine tracer_2d_1L(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, n end subroutine tracer_2d_1L !>@brief The subroutine 'tracer_2d' is the standard routine for sub-cycled tracer advection. -subroutine tracer_2d(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, npz, & - nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, lim_fac, dpA) +subroutine tracer_2d(q, dp1, mfx, mfy, cx, cy, cmax, imax_req, gridstruct, bd, domain, npx, npy, npz, & + nq, hord, n_sponge, q_split, dt, id_divg, q_pack, nord_tr, trdm, lim_fac, dpA) type(fv_grid_bounds_type), intent(IN) :: bd integer, intent(IN) :: npx @@ -371,17 +349,22 @@ subroutine tracer_2d(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, integer, intent(IN) :: npz integer, intent(IN) :: nq !< number of tracers to be advected integer, intent(IN) :: hord, nord_tr - integer, intent(IN) :: q_split + integer, intent(IN) :: n_sponge, q_split integer, intent(IN) :: id_divg real , intent(IN) :: dt, trdm real , intent(IN) :: lim_fac type(group_halo_update_type), intent(inout) :: q_pack real , intent(INOUT) :: q(bd%isd:bd%ied,bd%jsd:bd%jed,npz,nq) !< Tracers real , intent(INOUT) :: dp1(bd%isd:bd%ied,bd%jsd:bd%jed,npz) !< DELP before dyn_core - real , intent(IN ) :: mfx(bd%is:bd%ie+1,bd%js:bd%je, npz) !< Mass Flux X-Dir - real , intent(IN ) :: mfy(bd%is:bd%ie ,bd%js:bd%je+1,npz) !< Mass Flux Y-Dir - real , intent(IN ) :: cx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz) !< Courant Number X-Dir - real , intent(IN ) :: cy(bd%isd:bd%ied,bd%js :bd%je +1,npz) !< Courant Number Y-Dir + + real(kind=8) , intent(IN ) :: mfx(bd%is:bd%ie+1,bd%js:bd%je, npz) !< Mass Flux X-Dir + real(kind=8) , intent(IN ) :: mfy(bd%is:bd%ie ,bd%js:bd%je+1,npz) !< Mass Flux Y-Dir + real(kind=8) , intent(IN ) :: cx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz) !< Courant Number X-Dir + real(kind=8) , intent(IN ) :: cy(bd%isd:bd%ied,bd%js :bd%je +1,npz) !< Courant Number Y-Dir + + real , intent(INOUT) :: cmax(npz) + integer, intent(INOUT) :: imax_req + real , optional, intent(OUT) :: dpA(bd%is:bd%ie,bd%js:bd%je,npz)! DELP after advection type(fv_grid_type), intent(IN), target :: gridstruct type(domain2d), intent(INOUT) :: domain @@ -398,7 +381,11 @@ subroutine tracer_2d(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, real :: mfy2(bd%is:bd%ie,bd%js:bd%je+1,npz) real :: cx2(bd%is:bd%ie+1,bd%jsd:bd%jed, npz) real :: cy2(bd%isd:bd%ied,bd%js :bd%je +1, npz) - real :: cmax(npz) + + real, parameter :: cfl_tol = 1.25 ! 25% + logical :: severe_violation + character(len=128) :: error_msg + real :: frac, rdt integer :: maxsplt, ksplt(npz) integer :: i,j,k,n,it,iq @@ -429,7 +416,7 @@ subroutine tracer_2d(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, dy => gridstruct%dy !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx,dxa,dy, & -!$OMP sin_sg,cy,yfx,dya,dx,cmax,q_split) +!$OMP sin_sg,cy,yfx,dya,dx) do k=1,npz do j=jsd,jed do i=is,ie+1 @@ -449,41 +436,58 @@ subroutine tracer_2d(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, endif enddo enddo - if ( q_split == 0 ) then - cmax(k) = 0. - if ( k < npz/6 ) then - do j=js,je - do i=is,ie - cmax(k) = max( cmax(k), abs(cx(i,j,k)), abs(cy(i,j,k)) ) - enddo - enddo - else - do j=js,je - do i=is,ie - cmax(k) = max( cmax(k), max(abs(cx(i,j,k)),abs(cy(i,j,k)))+1.-sin_sg(i,j,5) ) - enddo - enddo - endif - endif enddo !-------------------------------------------------------------------------------- -! Determine global cmax on levels: if ( q_split == 0 ) then - call mp_reduce_max(cmax,npz) + !------------------------------------------------------- + ! ORIGINAL FV3 PATH (UNCHANGED) + !------------------------------------------------------- + call timing_on('COMM_TOTAL') + call timing_on('COMM_TRACER_MAX') + call mp_wait_reduce_max(cmax, npz, imax_req) + call timing_off('COMM_TRACER_MAX') + call timing_off('COMM_TOTAL') + maxsplt = 0 + do k=1,npz + ksplt(k) = int(1. + cmax(k)) + maxsplt = max(maxsplt, ksplt(k)) + enddo else - cmax = q_split-1. + !------------------------------------------------------- + ! USER-SUPPLIED q_split PATH (NO GLOBAL REDUCTION) + !------------------------------------------------------- + maxsplt = 0 + severe_violation = .false. + do k=1,npz + ksplt(k) = q_split + !! if (k <= n_sponge) ksplt(k) = ksplt(k) + 1 + maxsplt = max(maxsplt, ksplt(k)) + if (cmax(k) > real(ksplt(k))) then + !------------------------------- + ! mild violation → allow + !------------------------------- + if (cmax(k) <= cfl_tol * real(ksplt(k))) then + write(*,'(A,I4,A,F10.5,A,F10.5)') & + 'Tracer CFL warning (mild): k=', k, & + ' cmax=', cmax(k), & + ' ksplt(k)=', real(ksplt(k)) + !------------------------------- + ! severe violation → fail + !------------------------------- + else + severe_violation = .true. + write(error_msg,'(A,I4,A,F10.5,A,F10.5)') & + 'FATAL tracer_2d CFL violation at k=', k, & + ' cmax=', cmax(k), & + ' ksplt(k)=', real(ksplt(k)) + call mpp_error(FATAL, trim(error_msg)) + endif + endif + enddo endif -! Determine maxsplt for outer iteration loop. -! Advection will be done base on levels using ksplt - maxsplt = 0 - do k=1,npz - ksplt(k) = int(1. + cmax(k)) - maxsplt = max(maxsplt,ksplt(k)) - enddo - !-------------------------------------------------------------------------------- !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,cx2,xfx,mfx,mfx2,cy,cy2,yfx,mfy,mfy2,cmax,ksplt) & @@ -493,31 +497,40 @@ subroutine tracer_2d(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, mfy2(:,:,k)=mfy(:,:,k) cx2(:,:,k)=cx(:,:,k) cy2(:,:,k)=cy(:,:,k) - if( ksplt(k) /= 1 ) then - frac = 1. / real(ksplt(k)) - do j=jsd,jed - do i=is,ie+1 - cx2(i,j,k) = cx(i,j,k) * frac - xfx(i,j,k) = xfx(i,j,k) * frac - enddo - enddo - do j=js,je - do i=is,ie+1 - mfx2(i,j,k) = mfx(i,j,k) * frac - enddo - enddo - do j=js,je+1 - do i=isd,ied - cy2(i,j,k) = cy(i,j,k) * frac - yfx(i,j,k) = yfx(i,j,k) * frac - enddo - enddo - do j=js,je+1 - do i=is,ie - mfy2(i,j,k) = mfy(i,j,k) * frac - enddo - enddo - endif + frac = 1.0 / real(ksplt(k)) + do j=jsd,jed + do i=is,ie+1 + cx2(i,j,k) = cx(i,j,k) * frac + if (abs(cx2(i,j,k)) > 1.0) then + cx2(i,j,k) = sign(1.0, cx2(i,j,k)) + xfx(i,j,k) = xfx(i,j,k) * frac / abs(cx(i,j,k)*frac) + else + xfx(i,j,k) = xfx(i,j,k) * frac + endif + enddo + enddo + do j=js,je + do i=is,ie+1 + mfx2(i,j,k) = mfx(i,j,k) * frac + enddo + enddo + + do j=js,je+1 + do i=isd,ied + cy2(i,j,k) = cy(i,j,k) * frac + if (abs(cy2(i,j,k)) > 1.0) then + cy2(i,j,k) = sign(1.0, cy2(i,j,k)) + yfx(i,j,k) = yfx(i,j,k) * frac / abs(cy(i,j,k)*frac) + else + yfx(i,j,k) = yfx(i,j,k) * frac + endif + enddo + enddo + do j=js,je+1 + do i=is,ie + mfy2(i,j,k) = mfy(i,j,k) * frac + enddo + enddo enddo do it=1,maxsplt @@ -543,13 +556,11 @@ subroutine tracer_2d(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, do j=jsd,jed do i=is,ie ra_x(i,j) = area(i,j) + (xfx(i,j,k) - xfx(i+1,j,k)) - if (cx2(i,j,k) > 1.0) write(*,*) 'cx2(i,j,k) > 1.0 : ', cx2(i,j,k), i, j, k enddo enddo do j=js,je do i=isd,ied ra_y(i,j) = area(i,j) + (yfx(i,j,k) - yfx(i,j+1,k)) - if (cy2(i,j,k) > 1.0) write(*,*) 'cy2(i,j,k) > 1.0 : ', cy2(i,j,k), i, j, k enddo enddo @@ -602,7 +613,7 @@ subroutine tracer_2d(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, end subroutine tracer_2d -subroutine tracer_2d_nested(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, npz, & +subroutine tracer_2d_nested(q, dp1, mfx, mfy, cx, cy, cmax, imax_req, gridstruct, bd, domain, npx, npy, npz, & nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, & k_split, neststruct, parent_grid, lim_fac) @@ -619,10 +630,15 @@ subroutine tracer_2d_nested(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, np type(group_halo_update_type), intent(inout) :: q_pack real , intent(INOUT) :: q(bd%isd:bd%ied,bd%jsd:bd%jed,npz,nq) !< Tracers real , intent(INOUT) :: dp1(bd%isd:bd%ied,bd%jsd:bd%jed,npz) !< DELP before dyn_core - real , intent(INOUT) :: mfx(bd%is:bd%ie+1,bd%js:bd%je, npz) !< Mass Flux X-Dir - real , intent(INOUT) :: mfy(bd%is:bd%ie ,bd%js:bd%je+1,npz) !< Mass Flux Y-Dir - real , intent(INOUT) :: cx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz) !< Courant Number X-Dir - real , intent(INOUT) :: cy(bd%isd:bd%ied,bd%js :bd%je +1,npz) !< Courant Number Y-Dir + + real(kind=8) , intent(IN ) :: mfx(bd%is:bd%ie+1,bd%js:bd%je, npz) !< Mass Flux X-Dir + real(kind=8) , intent(IN ) :: mfy(bd%is:bd%ie ,bd%js:bd%je+1,npz) !< Mass Flux Y-Dir + real(kind=8) , intent(IN ) :: cx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz) !< Courant Number X-Dir + real(kind=8) , intent(IN ) :: cy(bd%isd:bd%ied,bd%js :bd%je +1,npz) !< Courant Number Y-Dir + + real , intent(INOUT) :: cmax(npz) + integer, intent(INOUT) :: imax_req + type(fv_grid_type), intent(IN), target :: gridstruct type(fv_nest_type), intent(INOUT) :: neststruct type(fv_atmos_type), intent(INOUT) :: parent_grid @@ -636,8 +652,12 @@ subroutine tracer_2d_nested(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, np real :: ra_y(bd%isd:bd%ied,bd%js:bd%je) real :: xfx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz) real :: yfx(bd%isd:bd%ied,bd%js: bd%je+1, npz) - real :: cmax(npz) - real :: cmax_t + + real :: mfxL(bd%is:bd%ie+1,bd%js:bd%je, npz) !< Mass Flux X-Dir + real :: mfyL(bd%is:bd%ie ,bd%js:bd%je+1,npz) !< Mass Flux Y-Dir + real :: cxL(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz) !< Courant Number X-Dir + real :: cyL(bd%isd:bd%ied,bd%js :bd%je +1,npz) !< Courant Number Y-Dir + real :: c_global real :: frac, rdt integer :: nsplt, nsplt_parent, msg_split_steps = 1 @@ -695,28 +715,11 @@ subroutine tracer_2d_nested(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, np if ( q_split == 0 ) then ! Determine nsplt -!$OMP parallel do default(none) shared(is,ie,js,je,npz,cmax,cx,cy,sin_sg) & -!$OMP private(cmax_t ) - do k=1,npz - cmax(k) = 0. - if ( k < 4 ) then -! Top layers: C < max( abs(c_x), abs(c_y) ) - do j=js,je - do i=is,ie - cmax_t = max( abs(cx(i,j,k)), abs(cy(i,j,k)) ) - cmax(k) = max( cmax_t, cmax(k) ) - enddo - enddo - else - do j=js,je - do i=is,ie - cmax_t = max(abs(cx(i,j,k)), abs(cy(i,j,k))) + 1.-sin_sg(i,j,5) - cmax(k) = max( cmax_t, cmax(k) ) - enddo - enddo - endif - enddo - call mp_reduce_max(cmax,npz) + call timing_on('COMM_TOTAL') + call timing_on('COMM_TRACER_WAIT') + call mp_wait_reduce_max(cmax, npz, imax_req) + call timing_off('COMM_TRACER_WAIT') + call timing_off('COMM_TOTAL') ! find global max courant number and define nsplt to scale cx,cy,mfx,mfy c_global = cmax(1) @@ -737,33 +740,38 @@ subroutine tracer_2d_nested(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, np frac = 1. / real(nsplt) if( nsplt /= 1 ) then -!$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,frac,xfx,mfx,cy,yfx,mfy) +!$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,cxL,frac,xfx,mfx,mfxL,cy,cyL,yfx,mfy,mfyL) do k=1,npz do j=jsd,jed do i=is,ie+1 - cx(i,j,k) = cx(i,j,k) * frac + cxL(i,j,k) = cx(i,j,k) * frac xfx(i,j,k) = xfx(i,j,k) * frac enddo enddo do j=js,je do i=is,ie+1 - mfx(i,j,k) = mfx(i,j,k) * frac + mfxL(i,j,k) = mfx(i,j,k) * frac enddo enddo do j=js,je+1 do i=isd,ied - cy(i,j,k) = cy(i,j,k) * frac + cyL(i,j,k) = cy(i,j,k) * frac yfx(i,j,k) = yfx(i,j,k) * frac enddo enddo do j=js,je+1 do i=is,ie - mfy(i,j,k) = mfy(i,j,k) * frac + mfyL(i,j,k) = mfy(i,j,k) * frac enddo enddo enddo + else + cxL= cx + cyL= cy + mfxL=mfx + mfyL=mfy endif @@ -787,14 +795,14 @@ subroutine tracer_2d_nested(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, np endif -!$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,mfx,mfy,rarea,nq, & -!$OMP area,xfx,yfx,q,cx,cy,npx,npy,hord,gridstruct,bd,it,nsplt,nord_tr,trdm,lim_fac) & +!$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,mfxL,mfyL,rarea,nq, & +!$OMP area,xfx,yfx,q,cxL,cyL,npx,npy,hord,gridstruct,bd,it,nsplt,nord_tr,trdm,lim_fac) & !$OMP private(dp2, ra_x, ra_y, fx, fy) do k=1,npz do j=js,je do i=is,ie - dp2(i,j) = dp1(i,j,k) + ((mfx(i,j,k)-mfx(i+1,j,k))+(mfy(i,j,k)-mfy(i,j+1,k)))*rarea(i,j) + dp2(i,j) = dp1(i,j,k) + ((mfxL(i,j,k)-mfxL(i+1,j,k))+(mfyL(i,j,k)-mfyL(i,j+1,k)))*rarea(i,j) enddo enddo @@ -811,14 +819,14 @@ subroutine tracer_2d_nested(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, np do iq=1,nq if ( it==1 .and. trdm>1.e-4 ) then - call fv_tp_2d(q(isd,jsd,k,iq), cx(is,jsd,k), cy(isd,js,k), & + call fv_tp_2d(q(isd,jsd,k,iq), cxL(is,jsd,k), cyL(isd,js,k), & npx, npy, hord, fx, fy, xfx(is,jsd,k), yfx(isd,js,k), & - gridstruct, bd, ra_x, ra_y, lim_fac, mfx=mfx(is,js,k), mfy=mfy(is,js,k), & + gridstruct, bd, ra_x, ra_y, lim_fac, mfx=mfxL(is,js,k), mfy=mfyL(is,js,k), & mass=dp1(isd,jsd,k), nord=nord_tr, damp_c=trdm) else - call fv_tp_2d(q(isd,jsd,k,iq), cx(is,jsd,k), cy(isd,js,k), & + call fv_tp_2d(q(isd,jsd,k,iq), cxL(is,jsd,k), cyL(isd,js,k), & npx, npy, hord, fx, fy, xfx(is,jsd,k), yfx(isd,js,k), & - gridstruct, bd, ra_x, ra_y, lim_fac, mfx=mfx(is,js,k), mfy=mfy(is,js,k)) + gridstruct, bd, ra_x, ra_y, lim_fac, mfx=mfxL(is,js,k), mfy=mfyL(is,js,k)) endif do j=js,je do i=is,ie @@ -896,13 +904,20 @@ subroutine offline_tracer_advection(q, pleB, pleA, mfx, mfy, cx, cy, & ! Local Arrays real :: xL(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz) ! X-Dir for MPP Updates real :: yL(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz) ! Y-Dir for MPP Updates - real :: cxL(bd%is :bd%ie +1,bd%jsd:bd%jed ,npz) ! Courant Number X-Dir - real :: cyL(bd%isd:bd%ied ,bd%js :bd%je +1,npz) ! Courant Number Y-Dir - real :: mfxL(bd%is :bd%ie +1,bd%js :bd%je ,npz) ! Mass Flux X-Dir - real :: mfyL(bd%is :bd%ie ,bd%js :bd%je +1,npz) ! Mass Flux Y-Dir + + real(kind=8) :: cxL(bd%is :bd%ie +1,bd%jsd:bd%jed ,npz) ! Courant Number X-Dir + real(kind=8) :: cyL(bd%isd:bd%ied ,bd%js :bd%je +1,npz) ! Courant Number Y-Dir + real(kind=8) :: mfxL(bd%is :bd%ie +1,bd%js :bd%je ,npz) ! Mass Flux X-Dir + real(kind=8) :: mfyL(bd%is :bd%ie ,bd%js :bd%je +1,npz) ! Mass Flux Y-Dir + + real :: cmax + real :: cmax_z(npz) + integer :: imax_req + real :: dpL(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz) ! Pressure Thickness real :: dpA(bd%is :bd%ie ,bd%js :bd%je ,npz) ! Pressure Thickness ! Local Tracer Arrays + real :: m1(bd%isd:bd%ied,npz,nq) real :: q3(bd%isd:bd%ied,bd%jsd:bd%jed, npz,nq)! Ghosted 3D Tracers real :: q2(bd%is :bd%ie , npz )! 2D Tmp integer :: kord_tracers(nq) @@ -914,14 +929,13 @@ subroutine offline_tracer_advection(q, pleB, pleA, mfx, mfy, cx, cy, & ! Local Remap Arrays real pe1(bd%is:bd%ie,npz+1) real pe2(bd%is:bd%ie,npz+1) - real dp1(bd%is:bd%ie,npz) real dp2(bd%is:bd%ie,npz) +! Local Array for scalings + real :: scalingFactors(nq) ! Local indices integer :: i,j,k,n,iq - real :: scalingFactor - type(group_halo_update_type), save :: i_pack integer :: is, ie, js, je @@ -950,6 +964,29 @@ subroutine offline_tracer_advection(q, pleB, pleA, mfx, mfy, cx, cy, & cxL(is:ie+1,jsd:jed,:) = xL(is:ie+1,jsd:jed,:) cyL(isd:ied,js:je+1,:) = yL(isd:ied,js:je+1,:) + ! --- 1. NON-BLOCKING CMAX REDUCTION START --- + if ( q_split == 0 ) then + ! A. Compute local cmax_z for this processor +!$OMP parallel do default(none) shared(npz,is,ie,js,je,cx,cy,cmax_z,gridstruct) & +!$OMP private(i,j,k,cmax) + do k = 1, npz + ! Compute local cmax_z (cx and cy are already hot in cache) + cmax_z(k) = 0. + do j=js,je + do i=is,ie + cmax = max(abs(cx(i,j,k)),abs(cy(i,j,k)))+1.-gridstruct%sin_sg(i,j,5) + cmax_z(k) = max( cmax_z(k), cmax ) + enddo + enddo + enddo + ! B. Initiate the non-blocking global reduction + call timing_on('COMM_TOTAL') + call timing_on('COMM_TRACER_MAX') + call mp_ireduce_max(cmax_z, npz, imax_req) + call timing_off('COMM_TRACER_MAX') + call timing_off('COMM_TOTAL') + endif + ! Fill MFX/MFY C-Grid boundaries xL(is:ie,js:je,:) = mfx(:,:,:) yL(is:ie,js:je,:) = mfy(:,:,:) @@ -962,19 +999,23 @@ subroutine offline_tracer_advection(q, pleB, pleA, mfx, mfy, cx, cy, & mfxL(is:ie+1,js:je,:) = xL(is:ie+1,js:je,:) mfyL(is:ie,js:je+1,:) = yL(is:ie,js:je+1,:) -! Fill local tracers and pressure thickness - dpL(bd%is:bd%ie,bd%js:bd%je,:) = pleB(:,:,2:npz+1) - pleB(:,:,1:npz) - q3(is:ie,js:je,:,:) = q(is:ie,js:je,:,:) +! Fill local tracer mass and pressure thickness + dpL(is:ie,js:je,:) = max(pleB(:,:,2:npz+1) - pleB(:,:,1:npz), 1.0) + do iq=1,nq + do k=1,npz + q3(is:ie,js:je,k,iq) = q(is:ie,js:je,k,iq) + enddo + enddo call start_group_halo_update(i_pack, q3, domain) if ( flagstruct%z_tracer .and. q_split==0 ) then - call tracer_2d_1L(q3, dpL, mfxL, mfyL, cxL, cyL, & + call tracer_2d_1L(q3, dpL, mfxL, mfyL, cxL, cyL, cmax_z, imax_req, & gridstruct, bd, domain, npx, npy, npz, nq, & flagstruct%hord_tr, dt, 0, i_pack, & flagstruct%nord_tr, flagstruct%trdm2, flagstruct%lim_fac, dpA=dpA) else - call tracer_2d(q3, dpL, mfxL, mfyL, cxL, cyL, gridstruct, bd, domain, npx, npy, npz, nq, & - flagstruct%hord_tr, q_split, dt, 0, i_pack, & + call tracer_2d(q3, dpL, mfxL, mfyL, cxL, cyL, cmax_z, imax_req, gridstruct, bd, domain, npx, npy, npz, nq, & + flagstruct%hord_tr, flagstruct%n_sponge, q_split, dt, 0, i_pack, & flagstruct%nord_tr, flagstruct%trdm2, flagstruct%lim_fac, dpA=dpA) endif @@ -984,88 +1025,201 @@ subroutine offline_tracer_advection(q, pleB, pleA, mfx, mfy, cx, cy, & kord_tracers = flagstruct%kord_tr !$OMP parallel do default(none) shared(is,ie,isd,ied,js,je,jsd,jed,npz,nq, & !$OMP pleB,dpA,pleA,q3,flagstruct,kord_tracers) & -!$OMP private(i,j,k,pe1,dp1,pe2,dp2) +!$OMP private(i,j,k,iq,pe1,pe2,dp2,m1) do j=js,je ! pressures mapping from (dpA is new delp after tracer_2d) pe1(:,1) = pleB(:,j,1) do k=2,npz+1 pe1(:,k) = pe1(:,k-1) + dpA(:,j,k-1) enddo - do k=1,npz - dp1(:,k) = pe1(:,k+1) - pe1(:,k) - enddo ! pressures mapping to pe2 = pleA(:,j,:) do k=1,npz - dp2(:,k) = pe2(:,k+1) - pe2(:,k) + dp2(:,k) = max(pe2(:,k+1) - pe2(:,k), 1.0) enddo call mapn_tracer(nq, npz, pe1, pe2, q3, dp2, kord_tracers, j, & is, ie, isd, ied, jsd, jed, 0., flagstruct%fill) enddo - ! Rescale tracers based on pleA at destination timestep - !------------------------------------------------------ + !------------------------------------------------------------------ + ! Rescale tracers globally + !------------------------------------------------------------------ + ! Calculate all scaling factors in one MPI pass + scalingFactors = calcScalingFactors(q(is:ie,js:je,1:npz,1:nq), q3(is:ie,js:je,1:npz,1:nq), & + pleB, pleA, npz, nq, gridstruct, bd) + !! Calculate bit-reproducible scaling factors (with verbose logging ON) + !scalingFactors = calcScalingFactors_BitReproducible(q(is:ie,js:je,1:npz,1:nq), q3(is:ie,js:je,1:npz,1:nq), & + ! pleB, pleA, npz, nq, gridstruct, bd, domain) + ! Apply scaling to each tracer do iq=1,nq - ! Scale tracers - !--------------- - scalingFactor = calcScalingFactor(q(is:ie,js:je,1:npz,iq), q3(is:ie,js:je,1:npz,iq), pleB, pleA, & - npz, domain, gridstruct, flagstruct, bd) - q(is:ie,js:je,1:npz,iq) = q3(is:ie,js:je,1:npz,iq) * scalingFactor + do k=1,npz + do j=js,je + do i=is,ie + q(i,j,k,iq) = q3(i,j,k,iq) * scalingFactors(iq) + enddo + enddo + enddo enddo end subroutine offline_tracer_advection !------------------------------------------------------------------------------------ - function calcScalingFactor(q1, q2, ple1, ple2, npz, domain, gridstruct, flagstruct, bd) result(scaling) - integer, intent(in) :: npz - type(fv_grid_bounds_type), intent(IN ) :: bd - real, intent(in) :: q1(bd%is:bd%ie,bd%js:bd%je,npz) - real, intent(in) :: q2(bd%is:bd%ie,bd%js:bd%je,npz) + function calcScalingFactors(q1, q2, ple1, ple2, npz, nq, gridstruct, bd, verbose_opt) result(scalings) + use mpp_mod, only: mpp_sum, mpp_pe, mpp_root_pe ! Added mpp_pe to check for master node + + integer, intent(in) :: npz, nq + type(fv_grid_bounds_type), intent(IN) :: bd + real, intent(in) :: q1(bd%is:bd%ie,bd%js:bd%je,npz,nq) + real, intent(in) :: q2(bd%is:bd%ie,bd%js:bd%je,npz,nq) real, intent(in) :: ple1(bd%is:bd%ie,bd%js:bd%je,npz+1) real, intent(in) :: ple2(bd%is:bd%ie,bd%js:bd%je,npz+1) + type(fv_grid_type), intent(IN) :: gridstruct + logical, intent(in), optional :: verbose_opt + + real :: scalings(nq) + + integer :: i, j, k, iq + real(REAL8) :: local_mass(2*nq) + real(REAL8) :: dp1, dp2, area_wt + real(REAL8) :: rel_error + real(REAL8), parameter :: TINY_DENOMINATOR = tiny(1.d0) + + logical :: verbose + + ! Set verbosity (default to false to keep logs clean) + verbose = .false. + if (present(verbose_opt)) verbose = verbose_opt + + local_mass = 0.d0 + + ! 1. Compute local domain sums directly (including area weighting) + do k=1,npz + do j=bd%js,bd%je + do i=bd%is,bd%ie + area_wt = gridstruct%area_64(i,j) + dp1 = (ple1(i,j,k+1) - ple1(i,j,k)) * area_wt + dp2 = (ple2(i,j,k+1) - ple2(i,j,k)) * area_wt + + do iq=1,nq + ! local_mass(iq) = Mass BEFORE advection + local_mass(iq) = local_mass(iq) + q1(i,j,k,iq) * dp1 + + ! local_mass(iq+nq) = Mass AFTER advection + local_mass(iq+nq) = local_mass(iq+nq) + q2(i,j,k,iq) * dp2 + enddo + enddo + enddo + enddo + + ! 2. ONE Single floating-point MPI sync + call mpp_sum(local_mass, 2*nq) + + ! 3. Calculate scaling factors and report adjustments + do iq=1,nq + if (local_mass(iq+nq) > TINY_DENOMINATOR) then + ! Calculate the scaling factor + scalings(iq) = REAL(local_mass(iq) / local_mass(iq+nq), KIND=kind(1.0)) + + ! Report exactly like the driver did (Master node only) + if (verbose) then + if (mpp_pe() == mpp_root_pe()) then + if (local_mass(iq) > 0.0_8) then + rel_error = (local_mass(iq+nq) - local_mass(iq)) / local_mass(iq) + + if (ABS(rel_error) >= epsilon(1.0)) then + write(6, 125) iq, rel_error + endif + endif + endif + endif + + else + scalings(iq) = 1.0 + end if + enddo + + 125 format('Mass Conservation Adjustment Required in offline_advection (Tracer ', I3, '): ', g21.14) + + end function calcScalingFactors + + function calcScalingFactors_BitReproducible(q1, q2, ple1, ple2, npz, nq, gridstruct, bd, domain, verbose_opt) result(scalings) + use mpp_mod, only: mpp_pe, mpp_root_pe + + integer, intent(in) :: npz, nq + type(fv_grid_bounds_type), intent(IN) :: bd + real, intent(in) :: q1(bd%is:bd%ie,bd%js:bd%je,npz,nq) + real, intent(in) :: q2(bd%is:bd%ie,bd%js:bd%je,npz,nq) + real, intent(in) :: ple1(bd%is:bd%ie,bd%js:bd%je,npz+1) + real, intent(in) :: ple2(bd%is:bd%ie,bd%js:bd%je,npz+1) + type(fv_grid_type), intent(IN) :: gridstruct type(domain2D), intent(INOUT) :: domain - type(fv_grid_type), intent(IN ) :: gridstruct - type(fv_flags_type), intent(INOUT) :: flagstruct - integer :: sflag - real :: scaling - - integer :: k - real(REAL8) :: qsum1(bd%is:bd%ie,bd%js:bd%je) - real(REAL8) :: qsum2(bd%is:bd%ie,bd%js:bd%je) - real(REAL8) :: globalSums(2) + logical, intent(in), optional :: verbose_opt + + real :: scalings(nq) + + integer :: i, j, k, iq + real(REAL8) :: qsum1(bd%is:bd%ie, bd%js:bd%je, nq) + real(REAL8) :: qsum2(bd%is:bd%ie, bd%js:bd%je, nq) + real(REAL8) :: global_mass1, global_mass2 + real(REAL8) :: dp1, dp2 + real(REAL8) :: rel_error real(REAL8), parameter :: TINY_DENOMINATOR = tiny(1.d0) - real(REAL8) :: scalingR8 - !------- - ! Compute partial sum on local array first to minimize communication. - ! This algorithm will not be strongly repdroducible under changes do domain - ! decomposition, but uses far less communication bandwidth (and memory BW) - ! then the preceding implementation. - !------- - qsum1(:,:) = 0.d0 - qsum2(:,:) = 0.d0 + + logical :: verbose + + ! Set verbosity (default to false) + verbose = .false. + if (present(verbose_opt)) verbose = verbose_opt + + qsum1 = 0.d0 + qsum2 = 0.d0 + + ! 1. Compute 2D vertical integrals for ALL tracers in ONE cache-friendly pass. + ! (No area weighting here, because g_sum_r8 does it internally) do k=1,npz - qsum1(:,:) = qsum1(:,:) + q1(:,:,k) * (ple1(:,:,k+1)-ple1(:,:,k)) - qsum2(:,:) = qsum2(:,:) + q2(:,:,k) * (ple2(:,:,k+1)-ple2(:,:,k)) + do j=bd%js,bd%je + do i=bd%is,bd%ie + dp1 = ple1(i,j,k+1) - ple1(i,j,k) + dp2 = ple2(i,j,k+1) - ple2(i,j,k) + + do iq=1,nq + qsum1(i,j,iq) = qsum1(i,j,iq) + q1(i,j,k,iq) * dp1 + qsum2(i,j,iq) = qsum2(i,j,iq) + q2(i,j,k,iq) * dp2 + enddo + enddo + enddo enddo - ! numerator - globalSums(1) = g_sum_r8(domain, qsum1, bd%is,bd%ie, bd%js,bd%je, 0, & - gridstruct%area_64(bd%is:bd%ie,bd%js:bd%je), 1) - ! denominator - globalSums(2) = g_sum_r8(domain, qsum2, bd%is,bd%ie, bd%js,bd%je, 0, & - gridstruct%area_64(bd%is:bd%ie,bd%js:bd%je), 1) - - if (globalSums(2) > TINY_DENOMINATOR) then - scalingR8 = globalSums(1) / globalSums(2) - !################################################################# - ! This line was added to ensure strong reproducibility of the code - !################################################################# - scaling = REAL(scalingR8, KIND=kind(1.00)) - else - scaling = 1.0 - end if + ! 2. Loop to call the Reproducible g_sum_r8 (2D slices) + do iq=1,nq + ! Mode 1 tells g_sum_r8 to divide by global_area internally + global_mass1 = g_sum_r8(domain, qsum1(:,:,iq), bd%is, bd%ie, bd%js, bd%je, 0, & + gridstruct%area_64(bd%is:bd%ie,bd%js:bd%je), 1) + + global_mass2 = g_sum_r8(domain, qsum2(:,:,iq), bd%is, bd%ie, bd%js, bd%je, 0, & + gridstruct%area_64(bd%is:bd%ie,bd%js:bd%je), 1) + + ! 3. Calculate scaling factor + if (global_mass2 > TINY_DENOMINATOR) then + scalings(iq) = REAL(global_mass1 / global_mass2, KIND=kind(1.0)) + + ! Report exactly like the driver did, if verbose is ON + if (verbose .and. mpp_pe() == mpp_root_pe()) then + if (global_mass1 > 0.0_8) then + rel_error = (global_mass2 - global_mass1) / global_mass1 + if (ABS(rel_error) >= epsilon(1.0)) then + write(6, 125) iq, rel_error + endif + endif + endif + else + scalings(iq) = 1.0 + end if + enddo + + 125 format('Mass Conservation Adjustment Required in offline_advection (Tracer ', I3, '): ', g21.14) - end function calcScalingFactor + end function calcScalingFactors_BitReproducible end module fv_tracer2d_mod diff --git a/model/nh_utils.F90 b/model/nh_utils.F90 index 3294f9860..1580845da 100644 --- a/model/nh_utils.F90 +++ b/model/nh_utils.F90 @@ -309,15 +309,15 @@ subroutine update_dz_d(ndif, damp, hord, is, ie, js, je, km, ng, npx, npy, area, !$OMP parallel do default(none) shared(is,ie,js,je,km,ws,zs,zh,rdt,dz_min) do j=js, je - do i=is,ie - ws(i,j) = ( zs(i,j) - zh(i,j,km+1) ) * rdt - enddo - do k=km, 1, -1 + do k=2, km+1 do i=is, ie ! Enforce monotonicity of height to prevent blowup - zh(i,j,k) = max( zh(i,j,k), zh(i,j,k+1) + dz_min ) + zh(i,j,k) = min( zh(i,j,k), zh(i,j,k-1) - dz_min ) enddo enddo + do i=is,ie + ws(i,j) = ( zs(i,j) - zh(i,j,km+1) ) * rdt + enddo enddo end subroutine update_dz_d @@ -347,14 +347,19 @@ subroutine Riem_Solver_c(ms, dt, is, ie, js, je, km, ng, & gama = 1./(1.-akap) rgrav = 1./grav - is1 = is - 1 - ie1 = ie + 1 - -!$OMP parallel do default(none) shared(js,je,is1,ie1,km,delp,pef,ptop,gz,rgrav,w3,pt, & +!$OMP parallel do default(none) shared(js,je,is,ie,km,delp,pef,ptop,gz,rgrav,w3,pt, & !$OMP a_imp,dt,gama,akap,ws,p_fac,scale_m,ms,hs,q_con,cappa) & -!$OMP private(cp2,gm2, dm, dz2, w2, pm2, pe2, pem, peg) +!$OMP private(is1, ie1, cp2, gm2, dm, dz2, w2, pm2, pe2, pem, peg) do 2000 j=js-1, je+1 + if ( (j < js) .or. (j > je) ) then + is1 = is + ie1 = ie + else + is1 = is - 1 + ie1 = ie + 1 + endif + do k=1,km do i=is1, ie1 dm(i,k) = delp(i,j,k) @@ -398,16 +403,23 @@ subroutine Riem_Solver_c(ms, dt, is, ie, js, je, km, ng, & enddo enddo - if ( a_imp < -0.01 ) then - call SIM3p0_solver(dt, is1, ie1, km, rdgas, gama, akap, pe2, dm, & - pem, w2, dz2, pt(is1:ie1,j,1:km), ws(is1,j), p_fac, scale_m) + call SIM3p0_solver(dt, is1, ie1, km, rdgas, gama, akap, & + pe2(is1:ie1,:), dm(is1:ie1,:), pem(is1:ie1,:), & + w2(is1:ie1,:), dz2(is1:ie1,:), & + pt(is1:ie1,j,1:km), ws(is1:ie1,j), p_fac, scale_m) + elseif ( a_imp <= 0.5 ) then - call RIM_2D(ms, dt, is1, ie1, km, rdgas, gama, gm2, pe2, & - dm, pm2, w2, dz2, pt(is1:ie1,j,1:km), ws(is1,j), .true.) + call RIM_2D(ms, dt, is1, ie1, km, rdgas, gama, gm2(is1:ie1,:), & + pe2(is1:ie1,:), dm(is1:ie1,:), pm2(is1:ie1,:), & + w2(is1:ie1,:), dz2(is1:ie1,:), & + pt(is1:ie1,j,1:km), ws(is1:ie1,j), .true.) + else - call SIM1_solver(dt, is1, ie1, km, rdgas, gama, gm2, cp2, akap, pe2, & - dm, pm2, pem, w2, dz2, pt(is1:ie1,j,1:km), ws(is1,j), p_fac) + call SIM1_solver(dt, is1, ie1, km, rdgas, gama, gm2(is1:ie1,:), & + cp2(is1:ie1,:), akap, pe2(is1:ie1,:), dm(is1:ie1,:), & + pm2(is1:ie1,:), pem(is1:ie1,:), w2(is1:ie1,:), dz2(is1:ie1,:), & + pt(is1:ie1,j,1:km), ws(is1:ie1,j), p_fac) endif do k=2,km+1 diff --git a/tools/fv_mp_mod.F90 b/tools/fv_mp_mod.F90 index 184f58154..e2de7e9fc 100644 --- a/tools/fv_mp_mod.F90 +++ b/tools/fv_mp_mod.F90 @@ -137,8 +137,9 @@ module fv_mp_mod integer :: isc, iec, jsc, jec public mp_start, mp_assign_gid, mp_barrier, mp_stop!, npes - public domain_decomp, mp_bcst, mp_reduce_max, mp_reduce_sum, mp_gather - public mp_reduce_min + public domain_decomp, mp_bcst, mp_gather + public mp_reduce_min, mp_reduce_max, mp_reduce_sum + public mp_ireduce_max, mp_wait_reduce_max public fill_corners, XDir, YDir public switch_current_domain, switch_current_Atm, broadcast_domains public is_master, setup_master @@ -201,6 +202,14 @@ module fv_mp_mod MODULE PROCEDURE mp_bcst_4d_i END INTERFACE + INTERFACE mp_ireduce_max + MODULE PROCEDURE mp_ireduce_max_r4_1d + END INTERFACE + + INTERFACE mp_wait_reduce_max + MODULE PROCEDURE mp_wait_reduce_max_r4_1d + END INTERFACE + !> The interface 'mp_reduce_min' contains routines that call SPMD_REDUCE. !! The routines compute the minima of values and place the !! absolute minimum value in a result. @@ -2396,6 +2405,45 @@ end subroutine mp_bcst_4d_i ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! !------------------------------------------------------------------------------- +!------------------------------------------------------------------------------- +! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv ! +! +! mp_ireduce_max_r4_1d :: START non-blocking SPMD REDUCE_MAX +! + subroutine mp_ireduce_max_r4_1d(mymax, npts, req) + integer, intent(IN) :: npts + real(kind=4), intent(IN) :: mymax(npts) + integer, intent(OUT) :: req ! MPI request handle + + req = MPI_REQUEST_NULL + ! Pass MPI_IN_PLACE as the send buffer, and mymax as the receive buffer + call MPI_IALLREDUCE( MPI_IN_PLACE, mymax, npts, MPI_REAL, MPI_MAX, & + commglobal, req, ierror ) + + end subroutine mp_ireduce_max_r4_1d +! +! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv ! +! +! mp_wait_reduce_max_r4_1d :: WAIT for non-blocking SPMD REDUCE_MAX +! + subroutine mp_wait_reduce_max_r4_1d(mymax, npts, req) + integer, intent(IN) :: npts + real(kind=4), intent(INOUT):: mymax(npts) + integer, intent(INOUT) :: req + + integer :: status(MPI_STATUS_SIZE) + + ! Wait for the IALLREDUCE to finish + call MPI_WAIT(req, status, ierror) + + end subroutine mp_wait_reduce_max_r4_1d +! +! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! +!------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv !