[BUG/ISSUE] Regrid_A2A_Mod is potentially broken for nested-grid simulations

jimmielin commented 4 years ago

Background

Regrid_A2A_Mod is used in GEOS-Chem Nested-Grid, WRF-GC and (temporarily) HEMCO_CESM to perform regridding of emission inventory data to the model grid.

In both WRF-GC and HEMCO_CESM I've been seeing horizontal regridding artifacts at CPU boundaries like so:

This led me to investigate Regrid_A2A_Mod.

Describe the bug

The issue is that Regrid_A2A_Mod when used within HEMCO's hco_interp_mod.F90 calls MAP_A2A with ig = 0, iv = 0 like so:

          ! Do the regridding
          CALL MAP_A2A( NX,      NY, LonEdgeI,    LatEdgeI, ORIG_2D,  &
                        HcoState%NX, HcoState%NY, LonEdgeO, LatEdgeO, &
                        REGR_2D, 0, 0, HCO_MISSVAL )

This triggers a particular treatment in Regrid_A2A_Mod.F90's YMAP_R*R* routine:

     !===================================================================
     ! Final processing for poles
     !===================================================================
     if ( ig .eq. 0 .and. iv .eq. 0 ) then

        ! South pole
        sum = 0.e+0_fp
        nlon= 0.0d0
        do i=1,im
           if(abs(q2(i,1)-miss)>tiny_r8 ) then
              sum = sum + q2(i,1)
              nlon= nlon + 1.0d0
           endif
        enddo

        if ( nlon > 0.0d0 ) sum = sum / nlon
        do i=1,im
           q2(i,1) = sum
        enddo

        ! North pole:
        sum = 0.e+0_fp
        nlon= 0.0d0
        do i=1,im
           if( abs(q2(i,jn)-miss)>tiny_r8 ) then
              sum = sum + q2(i,jn)
              nlon= nlon + 1.0d0
           endif
        enddo

        if ( nlon > 0.0d0 ) sum = sum / DBLE( im )
        do i=1,im
           q2(i,jn) = sum
        enddo

     endif

This is too wordy. Let me summarize the code:

! Final processing for poles
if(ig .eq. 0 .and. iv .eq. 0) then
    q2(:, 1) = sum(q2(:, 1) / im) ! im = number of horizontal dims
    q2(i, jn) = sum(q2, :, jn) / im) ! 1 is southernmost, jn is northernmost of the result/target grid (q2, dimensions im x jn)
endif

This is appropriate treating for poles as the earth is round. But this is not appropriate if the target domain as specified by LatEdgeO

    ! Get output grid edges from HEMCO state
    LonEdgeO(:) = HcoState%Grid%XEDGE%Val(:,1)
    LatEdgeO(:) = HcoState%Grid%YSIN%Val(1,:)

Is not spanning the full [-1, 1], aka when we are running on a regional grid.

What happens? The northern and southern edges of the resulting array are "averaged" into one, creating the artifact you see above.

Why hasn't it been noticed before?

In GEOS-Chem Nested Grid, it only averages the northern and southern edges of the domain. Generally the boundaries are fed with boundary conditions, so even if the emissions are averaged at the top/bottommost slices, it shouldn't be a problem.

In WRF-GC (and HEMCO_CESM), each CPU handles its own sub-domain; so every CPU-to-CPU northernmost/southernmost boundary is "smeared" with this "pole" treatment. It is not noticeable if the y-dimension is spread out wide enough, but it so happens that HEMCO_CESM has 80x2 dimensions as the decomposition, exacerbating this behavior.

Impacts

Emissions over northernmost / southernmost boundaries of a nested-grid simulation are not spatially accurate.
Emission total budgets are the same because it is simply an averaging.

Proposed solution

This section of the code should be blocked off for a nested-grid simulation, MODEL_WRF and MODEL_CESM. I think it should not be done via processor constants. I think the south pole and north pole treatments could be respectively blocked like

! South pole treatment
if(lon2(1, 1) .eq. -1) then   ! southernmost IS the south pole
...
endif

! North pole treatment
if(lon1(1, 1) .eq. 1) then   ! northernmost IS the north pole
...
endif

I'm happy to write a patch and further discuss the results. Thank you!

jimmielin commented 4 years ago

A tentative fix is in PR #274

I think it works. Still need to validate with GEOS-Chem main (verify it doesn't break global runs and correctly affects results of nested runs) and also merge into HEMCO/src/Shared/GeosUtil.

Pretty plot:

msulprizio commented 4 years ago

Pull request https://github.com/geoschem/geos-chem/pull/274 has now been merged into dev/12.8.1.

geoschem / geos-chem