test_regridxg_const Subroutine

subroutine test_regridxg_const(rc)

Arguments

Type IntentOptional Attributes Name
integer, intent(out) :: rc

Source Code

    subroutine test_regridxg_const(rc)
        integer, intent(out)                :: rc
        integer                             :: localrc, i
        type(ESMF_XGrid)                    :: xgrid
        type(ESMF_Grid)                     :: sideA(2), sideB(1)
        type(ESMF_DistGrid)                 :: distgrid
        real(ESMF_KIND_R8)                  :: centroid(12,2), area(12)
        type(ESMF_XGridSpec)                :: sparseMatA2X(2), sparseMatX2B(1)

        type(ESMF_Grid)                     :: l_sideA(2), l_sideB(1)
        type(ESMF_DistGrid)                 :: l_sideAdg(2), l_sideBdg(1)
        real(ESMF_KIND_R8)                  :: l_centroid(12,2), l_area(12)
        type(ESMF_XGridSpec)                :: l_sparseMatA2X(2), l_sparseMatX2B(1)
        type(ESMF_Field)                    :: field, srcField(2), dstField(1)

        integer                             :: eleCount
        integer                             :: sideAGC, sideBGC, sideAMC, sideBMC
        integer                             :: elb, eub, ec, lpet, npet

        real(ESMF_KIND_R8), pointer         :: farrayPtr(:,:), xfarrayPtr(:)
        real(ESMF_KIND_R8)                  :: xgrid_area(12), B_area(2,2)
        integer                             :: xlb(1), xub(1)
        type(ESMF_RouteHandle)              :: rh_src2xgrid(2), rh_xgrid2dst(1)
        type(ESMF_VM)                       :: vm

        real(ESMF_KIND_R8)                  :: centroidA1X(2), centroidA1Y(2)
        real(ESMF_KIND_R8)                  :: centroidA2X(2), centroidA2Y(1)
        real(ESMF_KIND_R8)                  :: centroidBX(2), centroidBY(2)
        real(ESMF_KIND_R8), pointer         :: coordX(:), coordY(:)
        character(len=16)                   :: gridNameA(2), gridNameB(1)

        rc = ESMF_SUCCESS
        localrc = ESMF_SUCCESS

        call ESMF_VMGetCurrent(vm=vm, rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return

        call ESMF_VMGet(vm, petCount=npet, localPet=lpet, rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return

        if(npet /= 2) return

        gridNameA(1) = 'srcGrid 1'
        gridNameA(2) = 'srcGrid 2'
        sideA(1) = ESMF_GridCreateNoPeriDim(maxIndex=(/2,2/), &
            coordDep1=(/1/), &
            coordDep2=(/2/), &
            regDecomp=(/1,2/), &
            name=gridNameA(1), rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return
        sideA(2) = ESMF_GridCreateNoPeriDim(maxIndex=(/2,1/), &
            coordDep1=(/1/), &
            coordDep2=(/2/), &
            regDecomp=(/2,1/), &
            name=gridNameA(2), rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return
        do i = 1, 2
            call ESMF_GridAddCoord(sideA(i), staggerloc=ESMF_STAGGERLOC_CENTER, &
                rc=localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
        enddo

        ! SideA first grid
        centroidA1X=(/0.5, 1.5/)
        centroidA1Y=(/0.5, 1.5/)
        call ESMF_GridGetCoord(sideA(1), localDE=0, &
            staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=1, farrayPtr=coordX, &
            rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return
        coordX = centroidA1X
        call ESMF_GridGetCoord(sideA(1), localDE=0, &
            staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=2, farrayPtr=coordY, &
            rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return
        coordY = centroidA1Y(lpet+1)

        ! SideA second grid
        centroidA2X=(/0.5, 1.5/)
        centroidA2Y=(/2.5/)
        call ESMF_GridGetCoord(sideA(2), localDE=0, &
            staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=1, farrayPtr=coordX, &
            rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return
        coordX = centroidA2X(lpet+1)
        call ESMF_GridGetCoord(sideA(2), localDE=0, &
            staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=2, farrayPtr=coordY, &
            rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return
        if(lpet == 0) coordY = centroidA2Y

        gridNameB(1) = 'dstGrid 1'
        sideB(1) = ESMF_GridCreateNoPeriDim(maxIndex=(/2,2/), &
            coordDep1=(/1/), coordDep2=(/2/), &
            regDecomp=(/1,2/), &
            name=gridNameB(1), rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return
        do i = 1, 1
            call ESMF_GridAddCoord(sideB(i), &
                staggerloc=ESMF_STAGGERLOC_CENTER, rc=localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
        enddo

        ! SideB grid
        centroidBX=(/0.75, 1.75/)
        centroidBY=(/0.75, 2.25/)
        call ESMF_GridGetCoord(sideB(1), localDE=0, &
            staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=1, farrayPtr=coordX, &
            rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return
        coordX = centroidBX
        call ESMF_GridGetCoord(sideB(1), localDE=0, &
            staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=2, farrayPtr=coordY, &
            rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return
        coordY = centroidBY(lpet+1)

        if(lpet == 0) then
            allocate(sparseMatA2X(1)%factorIndexList(2,6), sparseMatA2X(1)%factorList(6))
        else
            allocate(sparseMatA2X(1)%factorIndexList(2,7:9), sparseMatA2X(1)%factorList(7:9))
        endif
        if(lpet == 0) then
            allocate(sparseMatA2X(2)%factorIndexList(2,1), sparseMatA2X(2)%factorList(1))
        else
            allocate(sparseMatA2X(2)%factorIndexList(2,2:3), sparseMatA2X(2)%factorList(2:3))
        endif
        if(lpet == 0) then
            allocate(sparseMatX2B(1)%factorIndexList(2,6), sparseMatX2B(1)%factorList(6))
        else
            allocate(sparseMatX2B(1)%factorIndexList(2,7:12), sparseMatX2B(1)%factorList(7:12))
        endif
    
        ! factorIndexList
        ! setting up mapping between A1 -> X
        if(lpet == 0) then
            sparseMatA2X(1)%factorIndexList(1,1)=1
            sparseMatA2X(1)%factorIndexList(1,2)=2
            sparseMatA2X(1)%factorIndexList(1,3)=2
            sparseMatA2X(1)%factorIndexList(1,4)=3
            sparseMatA2X(1)%factorIndexList(1,5)=4
            sparseMatA2X(1)%factorIndexList(1,6)=4
        else
            sparseMatA2X(1)%factorIndexList(1,7)=3
            sparseMatA2X(1)%factorIndexList(1,8)=4
            sparseMatA2X(1)%factorIndexList(1,9)=4
        endif
        if(lpet == 0) then
            sparseMatA2X(1)%factorIndexList(2,1)=1
            sparseMatA2X(1)%factorIndexList(2,2)=2
            sparseMatA2X(1)%factorIndexList(2,3)=3
            sparseMatA2X(1)%factorIndexList(2,4)=4
            sparseMatA2X(1)%factorIndexList(2,5)=5
            sparseMatA2X(1)%factorIndexList(2,6)=6
        else
            sparseMatA2X(1)%factorIndexList(2,7)=7
            sparseMatA2X(1)%factorIndexList(2,8)=8
            sparseMatA2X(1)%factorIndexList(2,9)=9
        endif
        ! setting up mapping between A2 -> X
        if(lpet == 0) then
            sparseMatA2X(2)%factorIndexList(1,1)=1
        else
            sparseMatA2X(2)%factorIndexList(1,2)=2
            sparseMatA2X(2)%factorIndexList(1,3)=2
        endif
        if(lpet == 0) then
            sparseMatA2X(2)%factorIndexList(2,1)=10
        else
            sparseMatA2X(2)%factorIndexList(2,2)=11
            sparseMatA2X(2)%factorIndexList(2,3)=12
        endif

        ! Note that the weights are dest area weighted
        ! factorList
        ! setting up mapping between A1 -> X
        if(lpet == 0) then
            sparseMatA2X(1)%factorList(1)=1
            sparseMatA2X(1)%factorList(2)=1
            sparseMatA2X(1)%factorList(3)=1
            sparseMatA2X(1)%factorList(4)=1
            sparseMatA2X(1)%factorList(5)=1
            sparseMatA2X(1)%factorList(6)=1
        else
            sparseMatA2X(1)%factorList(7)=1
            sparseMatA2X(1)%factorList(8)=1
            sparseMatA2X(1)%factorList(9)=1
        endif
        ! setting up mapping between A2 -> X
        if(lpet == 0) then
            sparseMatA2X(2)%factorList(1)=1
        else
            sparseMatA2X(2)%factorList(2)=1
            sparseMatA2X(2)%factorList(3)=1
        endif
    
        ! factorIndexList
        ! setting up mapping between X -> B
        if(lpet == 0) then
            sparseMatX2B(1)%factorIndexList(1,1)=1
            sparseMatX2B(1)%factorIndexList(1,2)=2
            sparseMatX2B(1)%factorIndexList(1,3)=3
            sparseMatX2B(1)%factorIndexList(1,4)=4
            sparseMatX2B(1)%factorIndexList(1,5)=5
            sparseMatX2B(1)%factorIndexList(1,6)=6
        else
            sparseMatX2B(1)%factorIndexList(1,7)=7
            sparseMatX2B(1)%factorIndexList(1,8)=8
            sparseMatX2B(1)%factorIndexList(1,9)=9
            sparseMatX2B(1)%factorIndexList(1,10)=10
            sparseMatX2B(1)%factorIndexList(1,11)=11
            sparseMatX2B(1)%factorIndexList(1,12)=12
        endif
        if(lpet == 0) then
            sparseMatX2B(1)%factorIndexList(2,1)=1
            sparseMatX2B(1)%factorIndexList(2,2)=1
            sparseMatX2B(1)%factorIndexList(2,3)=2
            sparseMatX2B(1)%factorIndexList(2,4)=1
            sparseMatX2B(1)%factorIndexList(2,5)=1
            sparseMatX2B(1)%factorIndexList(2,6)=2
        else
            sparseMatX2B(1)%factorIndexList(2,7)=3
            sparseMatX2B(1)%factorIndexList(2,8)=3
            sparseMatX2B(1)%factorIndexList(2,9)=4
            sparseMatX2B(1)%factorIndexList(2,10)=3
            sparseMatX2B(1)%factorIndexList(2,11)=3
            sparseMatX2B(1)%factorIndexList(2,12)=4
        endif

        ! factorList
        ! setting up mapping between X -> B
        if(lpet == 0) then
            sparseMatX2B(1)%factorList(1)=4./9
            sparseMatX2B(1)%factorList(2)=2./9
            sparseMatX2B(1)%factorList(3)=2./3
            sparseMatX2B(1)%factorList(4)=2./9
            sparseMatX2B(1)%factorList(5)=1./9
            sparseMatX2B(1)%factorList(6)=1./3
        else
            sparseMatX2B(1)%factorList(7)=2./9
            sparseMatX2B(1)%factorList(8)=1./9
            sparseMatX2B(1)%factorList(9)=1./3
            sparseMatX2B(1)%factorList(10)=4./9
            sparseMatX2B(1)%factorList(11)=2./9
            sparseMatX2B(1)%factorList(12)=2./3
        endif

        ! set up destination areas to adjust weighted flux
        xgrid_area(1) = 1.
        xgrid_area(2) = 0.5
        xgrid_area(3) = 0.5
        xgrid_area(4) = 0.5
        xgrid_area(5) = 0.25
        xgrid_area(6) = 0.25
        xgrid_area(7) = 0.5
        xgrid_area(8) = 0.25
        xgrid_area(9) = 0.25
        xgrid_area(10) = 1.
        xgrid_area(11) = 0.5
        xgrid_area(12) = 0.5

        B_area(1,1) = 9./4
        B_area(2,1) = 3./4
        B_area(1,2) = 9./4
        B_area(2,2) = 3./4

        ! Finally ready to do an flux exchange from A side to B side
        xgrid = ESMF_XGridCreateFromSparseMat(sideAGrid=sideA, sideBGrid=sideB, &
            area=xgrid_area, centroid=centroid, &
            sparseMatA2X=sparseMatA2X, sparseMatX2B=sparseMatX2B, rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return

        call ESMF_XGridGet(xgrid, &
            sideAGridCount=sideAGC, sideBGridCount=sideBGC, &
            sideAMeshCount=sideAMC, sideBMeshCount=sideBMC, &
            sideAGrid=l_sideA, sideBGrid=l_sideB, area=l_area, &
            centroid=l_centroid, distgridA=l_sideAdg, &
            distgridM = distgrid, sparseMatA2X=l_sparseMatA2X, &
            sparseMatX2B=l_sparseMatX2B, &
            rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return

        call ESMF_XGridGet(xgrid, localDe=0, elementCount=eleCount, &
            exclusiveCount=ec, exclusiveLBound=elb, exclusiveUBound=eub, &
            rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return

        !print *, eleCount, elb, eub

        call ESMF_DistGridPrint(distgrid, rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return

        do i = 1, 2
            call ESMF_DistGridPrint(l_sideAdg(i), rc=localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
        enddo

        call ESMF_XGridGet(xgrid, xgridSide=ESMF_XGRIDSIDE_A, gridIndex=1, &
            distgrid=distgrid, rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return
        call ESMF_XGridGet(xgrid, xgridSide=ESMF_XGRIDSIDE_A, gridIndex=2, &
            distgrid=distgrid, rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return
        call ESMF_XGridGet(xgrid, xgridSide=ESMF_XGRIDSIDE_B, gridIndex=1, &
            distgrid=distgrid, rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return

        field = ESMF_FieldCreate(xgrid, typekind=ESMF_TYPEKIND_R8, rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return
        call ESMF_FieldGet(field, farrayPtr=xfarrayPtr, &
            exclusiveLBound=xlb, exclusiveUBound=xub, rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return

        xfarrayPtr = 0.0
        call ESMF_XGridGet(xgrid, elementCount=eleCount, rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return

        call ESMF_FieldPrint(field, rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return

        ! setup and initialize src and dst Fields
        do i = 1, 2
            srcField(i) = ESMF_FieldCreate(sideA(i), typekind=ESMF_TYPEKIND_R8, rc=localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
            call ESMF_FieldGet(srcField(i), farrayPtr=farrayPtr, rc=localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
            farrayPtr = 1.
        enddo
        do i = 1, 1
            dstField(i) = ESMF_FieldCreate(sideB(i), typekind=ESMF_TYPEKIND_R8, rc=localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
            call ESMF_FieldGet(dstField(i), farrayPtr=farrayPtr, rc=localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
            farrayPtr = 0.0
        enddo

        ! use field on Xgrid to do smm from src to dst transformation
        ! dst = W'*W*src
        !print *, '- before SMM from A -> X'
        !print *, xlb, xub, xfarrayPtr

        ! from A -> X
        do i = 1, 2
            call ESMF_FieldRegridStore(xgrid, srcField(i), field, routehandle=rh_src2xgrid(i), &
                rc = localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
            call ESMF_FieldRegrid(srcField(i), field, routehandle=rh_src2xgrid(i), &
                zeroregion=ESMF_REGION_SELECT, &
                checkflag=.TRUE.,&
                rc = localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
            !print *, '- SMM from A -> X'
            !print *, xfarrayPtr
        enddo

        ! xfarrayPtr should be all 1. at this point
        ! To get the surface integral of flux on XGrid, adjust by dst area

        !do i = xlb(1), xub(1)
        !    xfarrayPtr(i) = xfarrayPtr(i) * xgrid_area(i) 
        !enddo

        !print *, '- after SMM from A -> X'
        !print *, xfarrayPtr ! should be xgrid_area

        !print *, '- B before SMM from X -> B'
        !print *, farrayPtr ! should be 0.

        ! from X -> B
        do i = 1, 1
            call ESMF_FieldRegridStore(xgrid, field, dstField(i), routehandle=rh_xgrid2dst(i), &
                rc = localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
            call ESMF_FieldRegrid(field, dstField(i), routehandle=rh_xgrid2dst(i), &
                zeroregion=ESMF_REGION_SELECT, &
                checkflag=.TRUE.,&
                rc = localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
        enddo

        !print *, '- B after SMM from X -> B'
        !print *, farrayPtr ! should be 1/B_area

        call ESMF_FieldDestroy(field, rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return

        call ESMF_XGridDestroy(xgrid, rc=localrc)
        if (ESMF_LogFoundError(localrc, &
            ESMF_ERR_PASSTHRU, &
            ESMF_CONTEXT, rcToReturn=rc)) return

        do i = 1, 2
            call ESMF_FieldDestroy(srcField(i), rc = localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
            call ESMF_GridDestroy(sideA(i), rc = localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
        enddo

        do i = 1, 1
            call ESMF_FieldDestroy(dstField(i), rc = localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
            call ESMF_GridDestroy(sideB(i), rc = localrc)
            if (ESMF_LogFoundError(localrc, &
                ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
        enddo

        deallocate(sparseMatA2X(1)%factorIndexList, sparseMatA2X(1)%factorList)
        deallocate(sparseMatA2X(2)%factorIndexList, sparseMatA2X(2)%factorList)
        deallocate(sparseMatX2B(1)%factorIndexList, sparseMatX2B(1)%factorList)

    end subroutine test_regridxg_const