From 36ca473b324e9a734e949df14e2c068f5416912d Mon Sep 17 00:00:00 2001 From: Scott Rabenhorst Date: Mon, 18 Dec 2023 11:11:41 -0500 Subject: [PATCH 1/2] Update fvdycore regression fix with model/fv_tracer2d.F90 --- model/fv_tracer2d.F90 | 90 +++++++++++++++++++++++++------------------ 1 file changed, 52 insertions(+), 38 deletions(-) diff --git a/model/fv_tracer2d.F90 b/model/fv_tracer2d.F90 index 85f8c8b83..f7970bf29 100644 --- a/model/fv_tracer2d.F90 +++ b/model/fv_tracer2d.F90 @@ -838,7 +838,7 @@ subroutine tracer_2d_nested(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, np end subroutine tracer_2d_nested -subroutine offline_tracer_advection(q, ple0, ple1, mfx, mfy, cx, cy, & +subroutine offline_tracer_advection(q, pleB, pleA, mfx, mfy, cx, cy, & gridstruct, flagstruct, bd, domain, & ak, bk, ptop, npx, npy, npz, & nq, dt) @@ -856,8 +856,8 @@ subroutine offline_tracer_advection(q, ple0, ple1, mfx, mfy, cx, cy, & type(domain2D), intent(INOUT) :: domain real, intent(IN ) :: dt - real, intent(IN ) ::ple0(bd%is:bd%ie,bd%js:bd%je,npz+1) ! DELP before dyn_core - real, intent(INOUT) ::ple1(bd%is:bd%ie,bd%js:bd%je,npz+1) ! DELP after dyn_core + real, intent(IN ) ::pleB(bd%is:bd%ie,bd%js:bd%je,npz+1) ! DELP before dyn_core + real, intent(INOUT) ::pleA(bd%is:bd%ie,bd%js:bd%je,npz+1) ! DELP after dyn_core real, intent(IN ) :: cx(bd%is:bd%ie,bd%js:bd%je,npz) ! Courant Number X-Dir real, intent(IN ) :: cy(bd%is:bd%ie,bd%js:bd%je,npz) ! Courant Number Y-Dir real, intent(IN ) :: mfx(bd%is:bd%ie,bd%js:bd%je,npz) ! Mass Flux X-Dir @@ -876,7 +876,7 @@ subroutine offline_tracer_advection(q, ple0, ple1, mfx, mfy, cx, cy, & 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 :: q1(bd%is :bd%ie ,bd%js :bd%je , npz )! 3D Tracers + real :: q1(bd%is :bd%ie , npz )! 2D Tracers real :: q3(bd%isd:bd%ied,bd%jsd:bd%jed, npz,nq)! 3D Tracers ! Local Buffer Arrarys real :: wbuffer(bd%js:bd%je,npz) @@ -886,7 +886,9 @@ subroutine offline_tracer_advection(q, ple0, ple1, mfx, mfy, cx, cy, & ! Local Remap Arrays real pe1(bd%is:bd%ie,npz+1) real pe2(bd%is:bd%ie,npz+1) - real dp2(bd%is:bd%ie,bd%js:bd%je,npz) + real dp2(bd%is:bd%ie,npz) + real ple1(bd%is:bd%ie,bd%js:bd%je,npz+1) + real ple2(bd%is:bd%ie,bd%js:bd%je,npz+1) integer kord_tracers(nq) ! Local indices @@ -935,7 +937,7 @@ subroutine offline_tracer_advection(q, ple0, ple1, mfx, mfy, cx, cy, & 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,:) = ple0(:,:,2:npz+1) - ple0(:,:,1:npz) + 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,:,:) call start_group_halo_update(i_pack, q3, domain) @@ -950,6 +952,28 @@ subroutine offline_tracer_advection(q, ple0, ple1, mfx, mfy, cx, cy, & flagstruct%nord_tr, flagstruct%trdm2, flagstruct%lim_fac, dpA=dpA) endif + + ! pressures mapping from (dpA is new delp after tracer_2d) + ple1(:,:,1) = ptop + do k=2,npz+1 + do j=js,je + do i=is,ie + ple1(i,j,k) = ple1(i,j,k-1) + dpA(i,j,k-1) + enddo + enddo + enddo + + ! pressures mapping to + ple2(:,:,1) = ptop + ple2(:,:,npz+1) = pleA(:,:,npz+1) + do k=2,npz + do j=js,je + do i=is,ie + ple2(i,j,k) = ak(k) + bk(k)*ple2(i,j,npz+1) + enddo + enddo + enddo + !------------------------------------------------------------------ ! Re-Map constituents !------------------------------------------------------------------ @@ -958,54 +982,44 @@ subroutine offline_tracer_advection(q, ple0, ple1, mfx, mfy, cx, cy, & kord_tracers(iq) = flagstruct%kord_tr enddo do j=js,je - ! pressures mapping from (dpA is new delp after tracer_2d) - pe1(:,1) = ptop - do k=2,npz+1 - pe1(:,k) = pe1(:,k-1) + dpA(:,j,k-1) - enddo - ! pressures mapping to - pe2(:,1) = ptop - pe2(:,npz+1) = pe1(:,npz+1) - do k=2,npz - pe2(: ,k) = ak(k) + bk(k)*pe1(:,npz+1) + do k=1,npz+1 + pe1(:,k) = ple1(:,j,k) + pe2(:,k) = ple2(:,j,k) enddo do k=1,npz - dp2(:,j,k) = pe2(:,k+1) - pe2(:,k) + dp2(:,k) = pe2(:,k+1) - pe2(:,k) enddo - call mapn_tracer(nq, npz, pe1, pe2, q3, dp2(:,j,:), kord_tracers, j, & + call mapn_tracer(nq, npz, pe1, pe2, q3, dp2, kord_tracers, j, & is, ie, isd, ied, jsd, jed, 0., flagstruct%fill) enddo elseif ( nq > 0 ) then do iq=1,nq do j=js,je - ! pressures mapping from (dpA is new delp after tracer_2d) - pe1(:,1) = ptop - do k=2,npz+1 - pe1(:,k) = pe1(:,k-1) + dpA(:,j,k-1) - enddo - ! pressures mapping to - pe2(:,1) = ptop - pe2(:,npz+1) = pe1(:,npz+1) - do k=2,npz - pe2(: ,k) = ak(k) + bk(k)*pe1(:,npz+1) + do k=1,npz+1 + pe1(:,k) = ple1(:,j,k) + pe2(:,k) = ple2(:,j,k) enddo do k=1,npz - dp2(:,j,k) = pe2(:,k+1) - pe2(:,k) + dp2(:,k) = pe2(:,k+1) - pe2(:,k) enddo call map1_q2(npz, pe1, q3(isd,jsd,1,iq), & - npz, pe2, q1(:,j,:), dp2(:,j,:), & + npz, pe2, q1, dp2, & is, ie, 0, flagstruct%kord_tr, j,& isd, ied, jsd, jed, 0.) - if (flagstruct%fill) call fillz(ie-is+1, npz, 1, q1(:,j,:), dp2(:,j,:)) - q3(is:ie,j,1:npz,iq) = q1(:,j,:) + if (flagstruct%fill) call fillz(ie-is+1, npz, 1, q1, dp2) + do k=1,npz + do i=is,ie + q3(i,j,k,iq) = q1(i,k) + enddo + enddo enddo enddo end if - ! Rescale tracers based on ple1 at destination timestep + ! Rescale tracers based on pleA at destination timestep !------------------------------------------------------ do iq=1,nq - scalingFactor = calcScalingFactor(q3(is:ie,js:je,1:npz,iq), dp2, ple1, npx, npy, npz, gridstruct, bd) + scalingFactor = calcScalingFactor(q3(is:ie,js:je,1:npz,iq), ple2, pleA, npx, npy, npz, gridstruct, bd) ! Return tracers !--------------- q(is:ie,js:je,1:npz,iq) = q3(is:ie,js:je,1:npz,iq) * scalingFactor @@ -1015,14 +1029,14 @@ end subroutine offline_tracer_advection !------------------------------------------------------------------------------------ - function calcScalingFactor(q1, dp2, ple1, npx, npy, npz, gridstruct, bd) result(scaling) + function calcScalingFactor(q1, ple1, ple2, npx, npy, npz, gridstruct, bd) result(scaling) use mpp_mod, only: mpp_sum integer, intent(in) :: npx integer, intent(in) :: npy integer, intent(in) :: npz real, intent(in) :: q1(:,:,:) - real, intent(in) :: dp2(:,:,:) real, intent(in) :: ple1(:,:,:) + real, intent(in) :: ple2(:,:,:) type(fv_grid_type), intent(IN ) :: gridstruct type(fv_grid_bounds_type), intent(IN ) :: bd real :: scaling @@ -1039,9 +1053,9 @@ function calcScalingFactor(q1, dp2, ple1, npx, npy, npz, gridstruct, bd) result( !------- do k = 1, npz ! numerator - partialSums(1,k) = sum(q1(:,:,k)*dp2(:,:,k)*gridstruct%area(bd%is:bd%ie,bd%js:bd%je)) + partialSums(1,k) = sum(q1(:,:,k)*(ple1(:,:,k+1)-ple1(:,:,k))*gridstruct%area(bd%is:bd%ie,bd%js:bd%je)) ! denominator - partialSums(2,k) = sum(q1(:,:,k)*(ple1(:,:,k+1)-ple1(:,:,k))*gridstruct%area(bd%is:bd%ie,bd%js:bd%je)) + partialSums(2,k) = sum(q1(:,:,k)*(ple2(:,:,k+1)-ple2(:,:,k))*gridstruct%area(bd%is:bd%ie,bd%js:bd%je)) end do globalSums(1) = sum(partialSums(1,:)) From bec1af399ce51c3f091b95c64c5b1ca59cbdf3d7 Mon Sep 17 00:00:00 2001 From: Scott Rabenhorst Date: Mon, 18 Dec 2023 12:50:28 -0500 Subject: [PATCH 2/2] Update geos/develop with fixes to offline tracer mass conservation and regression --- model/fv_grid_utils.F90 | 63 ++++++++++++++++++++--- model/fv_tracer2d.F90 | 109 +++++++++++++++++++++++----------------- 2 files changed, 119 insertions(+), 53 deletions(-) diff --git a/model/fv_grid_utils.F90 b/model/fv_grid_utils.F90 index be853b98d..c509717c1 100644 --- a/model/fv_grid_utils.F90 +++ b/model/fv_grid_utils.F90 @@ -84,7 +84,7 @@ module fv_grid_utils_mod use mpp_parameter_mod, only: CORNER, SCALAR_PAIR use fv_arrays_mod, only: fv_atmos_type, fv_grid_type, fv_grid_bounds_type, & - R_GRID + R_GRID, FVPRC, REAL4, REAL8 use fv_eta_mod, only: set_eta use fv_mp_mod, only: ng, is_master use fv_mp_mod, only: mp_reduce_sum, mp_reduce_min, mp_reduce_max @@ -108,12 +108,15 @@ module fv_grid_utils_mod real, parameter:: ptop_min=1.d-8 + logical, SAVE :: g_sum_initialized = .false. + real(kind=R_GRID), SAVE :: global_area + public f_p public ptop_min, big_number !CLEANUP: OK to keep since they are constants? public cos_angle public latlon2xyz, gnomonic_grids, gnomonic_grids_local, & global_mx, unit_vect_latlon, & - cubed_to_latlon, c2l_ord2, g_sum, global_qsum, great_circle_dist, & + cubed_to_latlon, c2l_ord2, g_sum, g_sum_r8, global_qsum, great_circle_dist, & v_prod, get_unit_vect2, project_sphere_v public mid_pt_sphere, mid_pt_cart, vect_cross, grid_utils_init, grid_utils_end, & spherical_angle, cell_center2, get_area, inner_prod, fill_ghost, direct_transform, & @@ -2900,13 +2903,11 @@ real function g_sum(domain, p, ifirst, ilast, jfirst, jlast, ngc, area, mode, re integer, intent(IN) :: jfirst, jlast, ngc integer, intent(IN) :: mode ! if ==1 divided by area logical, intent(in), optional :: reproduce - real, intent(IN) :: p(ifirst:ilast,jfirst:jlast) ! field to be summed + real(kind=REAL4), intent(IN) :: p(ifirst:ilast,jfirst:jlast) ! field to be summed real(kind=R_GRID), intent(IN) :: area(ifirst-ngc:ilast+ngc,jfirst-ngc:jlast+ngc) type(domain2d), intent(IN) :: domain integer :: i,j - real gsum - logical, SAVE :: g_sum_initialized = .false. - real(kind=R_GRID), SAVE :: global_area + real(kind=REAL4) gsum real :: tmp(ifirst:ilast,jfirst:jlast) integer :: err @@ -2946,6 +2947,56 @@ real function g_sum(domain, p, ifirst, ilast, jfirst, jlast, ngc, area, mode, re endif end function g_sum + !>@brief The function 'g_sum_r8' is the fast version of 'globalsum'. + real(kind=REAL8) function g_sum_r8(domain, p, ifirst, ilast, jfirst, jlast, ngc, area, mode, reproduce) + integer, intent(IN) :: ifirst, ilast + integer, intent(IN) :: jfirst, jlast, ngc + integer, intent(IN) :: mode ! if ==1 divided by area + logical, intent(in), optional :: reproduce + real(kind=REAL8), intent(IN) :: p(ifirst:ilast,jfirst:jlast) ! field to be summed + real(kind=R_GRID), intent(IN) :: area(ifirst-ngc:ilast+ngc,jfirst-ngc:jlast+ngc) + type(domain2d), intent(IN) :: domain + integer :: i,j + real(kind=REAL8) :: gsum + real(kind=REAL8) :: tmp(ifirst:ilast,jfirst:jlast) + integer :: err + + if ( .not. g_sum_initialized ) then + global_area = mpp_global_sum(domain, area, flags=BITWISE_EFP_SUM) + if ( is_master() ) write(*,*) 'Global Area=',global_area + g_sum_initialized = .true. + end if + +!------------------------- +! FMS global sum algorithm: +!------------------------- + if ( present(reproduce) ) then + if (reproduce) then + gsum = mpp_global_sum(domain, p(:,:)*area(ifirst:ilast,jfirst:jlast), & + flags=BITWISE_EFP_SUM) + else + gsum = mpp_global_sum(domain, p(:,:)*area(ifirst:ilast,jfirst:jlast)) + endif + else +!------------------------- +! Quick local sum algorithm +!------------------------- + gsum = 0. + do j=jfirst,jlast + do i=ifirst,ilast + gsum = gsum + p(i,j)*area(i,j) + enddo + enddo + call mp_reduce_sum(gsum) + endif + + if ( mode==1 ) then + g_sum_r8 = gsum / global_area + else + g_sum_r8 = gsum + endif + end function g_sum_r8 + !>@brief The function 'global_qsum' computes the quick global sum without area weighting. real function global_qsum(p, ifirst, ilast, jfirst, jlast) integer, intent(IN) :: ifirst, ilast diff --git a/model/fv_tracer2d.F90 b/model/fv_tracer2d.F90 index f7970bf29..0703effbe 100644 --- a/model/fv_tracer2d.F90 +++ b/model/fv_tracer2d.F90 @@ -67,11 +67,12 @@ module fv_tracer2d_mod 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 - use mpp_domains_mod, only: mpp_update_domains, CGRID_NE, domain2d, mpp_get_boundary + use mpp_domains_mod, only: mpp_update_domains, CGRID_NE, domain2d, mpp_get_boundary, mpp_global_sum, BITWISE_EFP_SUM use fv_timing_mod, only: timing_on, timing_off use boundary_mod, only: nested_grid_BC_apply_intT - use fv_arrays_mod, only: fv_grid_type, fv_flags_type, fv_nest_type, fv_atmos_type, fv_grid_bounds_type + use fv_arrays_mod, only: fv_grid_type, fv_flags_type, fv_nest_type, fv_atmos_type, fv_grid_bounds_type, REAL8 use mpp_mod, only: mpp_error, FATAL, mpp_broadcast, mpp_send, mpp_recv, mpp_sum, mpp_max + use fv_grid_utils_mod, only: g_sum_r8 implicit none private @@ -888,8 +889,10 @@ subroutine offline_tracer_advection(q, pleB, pleA, mfx, mfy, cx, cy, & real pe2(bd%is:bd%ie,npz+1) real dp2(bd%is:bd%ie,npz) real ple1(bd%is:bd%ie,bd%js:bd%je,npz+1) - real ple2(bd%is:bd%ie,bd%js:bd%je,npz+1) integer kord_tracers(nq) +! Local mass conservation + real(REAL8) :: qsum(bd%is:bd%ie,bd%js:bd%je) + real(REAL8) :: g_massB, g_massA ! Local indices integer :: i,j,k,n,iq @@ -963,17 +966,6 @@ subroutine offline_tracer_advection(q, pleB, pleA, mfx, mfy, cx, cy, & enddo enddo - ! pressures mapping to - ple2(:,:,1) = ptop - ple2(:,:,npz+1) = pleA(:,:,npz+1) - do k=2,npz - do j=js,je - do i=is,ie - ple2(i,j,k) = ak(k) + bk(k)*ple2(i,j,npz+1) - enddo - enddo - enddo - !------------------------------------------------------------------ ! Re-Map constituents !------------------------------------------------------------------ @@ -984,7 +976,7 @@ subroutine offline_tracer_advection(q, pleB, pleA, mfx, mfy, cx, cy, & do j=js,je do k=1,npz+1 pe1(:,k) = ple1(:,j,k) - pe2(:,k) = ple2(:,j,k) + pe2(:,k) = pleA(:,j,k) enddo do k=1,npz dp2(:,k) = pe2(:,k+1) - pe2(:,k) @@ -997,7 +989,7 @@ subroutine offline_tracer_advection(q, pleB, pleA, mfx, mfy, cx, cy, & do j=js,je do k=1,npz+1 pe1(:,k) = ple1(:,j,k) - pe2(:,k) = ple2(:,j,k) + pe2(:,k) = pleA(:,j,k) enddo do k=1,npz dp2(:,k) = pe2(:,k+1) - pe2(:,k) @@ -1016,61 +1008,84 @@ subroutine offline_tracer_advection(q, pleB, pleA, mfx, mfy, cx, cy, & enddo end if - ! Rescale tracers based on pleA at destination timestep - !------------------------------------------------------ - do iq=1,nq - scalingFactor = calcScalingFactor(q3(is:ie,js:je,1:npz,iq), ple2, pleA, npx, npy, npz, gridstruct, bd) - ! Return tracers - !--------------- - q(is:ie,js:je,1:npz,iq) = q3(is:ie,js:je,1:npz,iq) * scalingFactor - enddo +! Compute Global Mass +! ------------------- + qsum(:,:) = 0.d0 + do k=1,npz + qsum(:,:) = qsum(:,:) + (pleB(:,:,k+1)-pleB(:,:,k)) + enddo + g_massB = g_sum_r8(domain, qsum, is,ie, js,je, 0, & + gridstruct%area_64(is:ie,js:je), 1, .true.) + qsum(:,:) = 0.d0 + do k=1,npz + qsum(:,:) = qsum(:,:) + (pleA(:,:,k+1)-pleA(:,:,k)) + enddo + g_massA = g_sum_r8(domain, qsum, is,ie, js,je, 0, & + gridstruct%area_64(is:ie,js:je), 1, .true.) + + ! Rescale tracers based on pleB and pleA for mass conservation + !------------------------------------------------------------- + do iq=1,nq + scalingFactor = calcScalingFactor(q(is:ie,js:je,1:npz,iq), q3(is:ie,js:je,1:npz,iq), & + pleB, pleA, g_massB, g_massA, npz, domain, gridstruct, bd) + ! Return tracers + !--------------- + q(is:ie,js:je,1:npz,iq) = q3(is:ie,js:je,1:npz,iq) * scalingFactor + enddo end subroutine offline_tracer_advection !------------------------------------------------------------------------------------ - function calcScalingFactor(q1, ple1, ple2, npx, npy, npz, gridstruct, bd) result(scaling) - use mpp_mod, only: mpp_sum - integer, intent(in) :: npx - integer, intent(in) :: npy + function calcScalingFactor(q1, q2, ple1, ple2, g_mass1, g_mass2, npz, domain, gridstruct, bd) result(scaling) integer, intent(in) :: npz - real, intent(in) :: q1(:,:,:) - real, intent(in) :: ple1(:,:,:) - real, intent(in) :: ple2(:,:,:) - type(fv_grid_type), intent(IN ) :: gridstruct 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) + 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) + real(REAL8), intent(in) :: g_mass1, g_mass2 + type(domain2D), intent(INOUT) :: domain + type(fv_grid_type), intent(IN ) :: gridstruct real :: scaling integer :: k - real :: partialSums(2,npz), globalSums(2) - real, parameter :: TINY_DENOMINATOR = tiny(1.0) - + 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) + 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. !------- - do k = 1, npz - ! numerator - partialSums(1,k) = sum(q1(:,:,k)*(ple1(:,:,k+1)-ple1(:,:,k))*gridstruct%area(bd%is:bd%ie,bd%js:bd%je)) - ! denominator - partialSums(2,k) = sum(q1(:,:,k)*(ple2(:,:,k+1)-ple2(:,:,k))*gridstruct%area(bd%is:bd%ie,bd%js:bd%je)) - end do + qsum1(:,:) = 0.d0 + qsum2(:,:) = 0.d0 + do k=1,npz + qsum1(:,:) = qsum1(:,:) + q1(:,:,k) * (ple1(:,:,k+1)-ple1(:,:,k)) + qsum2(:,:) = qsum2(:,:) + q2(:,:,k) * (ple2(:,:,k+1)-ple2(:,:,k)) + enddo - globalSums(1) = sum(partialSums(1,:)) - globalSums(2) = sum(partialSums(2,:)) + ! 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, .true.) + ! 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, .true.) - call mpp_sum(globalSums, 2) + if (g_mass1 > 0.0) globalSums(1) = globalSums(1)/g_mass1 + if (g_mass2 > 0.0) globalSums(2) = globalSums(2)/g_mass2 if (globalSums(2) > TINY_DENOMINATOR) then - scaling = globalSums(1) / globalSums(2) + scalingR8 = globalSums(1) / globalSums(2) !################################################################# ! This line was added to ensure strong reproducibility of the code !################################################################# - scaling = REAL(scaling, KIND=kind(1.00)) + scaling = REAL(scalingR8, KIND=kind(1.00)) else - scaling = 1.d0 + scaling = 1.0 end if end function calcScalingFactor