! $Id: ESMF_FieldRegridMaskEx.F90,v 1.24 2011/07/02 05:53:55 oehmke Exp $ ! ! Earth System Modeling Framework ! Copyright 2002-2011, University Corporation for Atmospheric Research, ! Massachusetts Institute of Technology, Geophysical Fluid Dynamics ! Laboratory, University of Michigan, National Centers for Environmental ! Prediction, Los Alamos National Laboratory, Argonne National Laboratory, ! NASA Goddard Space Flight Center. ! Licensed under the University of Illinois-NCSA License. ! !============================================================================== ! program ESMF_FieldRegridEx !============================================================================== !ESMF_MULTI_PROC_EXAMPLE String used by test script to count examples. !============================================================================== !------------------------------------------------------------------------------ #include "ESMF_Macros.inc" ! !USES: use ESMF use ESMF_TestMod ! test methods use ESMF_RegridMod use ESMF_FieldMod use ESMF_GridUtilMod use ESMF_FieldGetMod implicit none !------------------------------------------------------------------------------ ! The following line turns the CVS identifier string into a printable variable. character(*), parameter :: version = & '$Id: ESMF_FieldRegridMaskEx.F90,v 1.24 2011/07/02 05:53:55 oehmke Exp $' !------------------------------------------------------------------------------ ! individual test result code integer :: localrc, rc, petCount,localPet ! individual test failure message character(ESMF_MAXSTR) :: name logical :: correct type(ESMF_Grid) :: gridSrc type(ESMF_Grid) :: gridDst type(ESMF_Field) :: srcField type(ESMF_Field) :: dstField type(ESMF_Array) :: dstArray type(ESMF_Array) :: srcArray type(ESMF_RouteHandle) :: routeHandle type(ESMF_ArraySpec) :: arrayspec type(ESMF_VM) :: vm real(ESMF_KIND_R8), pointer :: farrayPtrXC(:,:) real(ESMF_KIND_R8), pointer :: farrayPtrYC(:,:) real(ESMF_KIND_R8), pointer :: farrayPtr(:,:) integer(ESMF_KIND_I4), pointer :: maskSrc(:,:), maskDst(:,:) integer :: clbnd(2),cubnd(2) integer :: fclbnd(2),fcubnd(2) integer :: i1,i2 integer :: lDE, localDECount integer src_nx, src_ny, dst_nx, dst_ny integer num_arrays real(ESMF_KIND_R8) :: src_dx, src_dy real(ESMF_KIND_R8) :: dst_dx, dst_dy real(ESMF_KIND_R8) :: ctheta, stheta real(ESMF_KIND_R8) :: theta, d2rad, xtmp, x, y integer(ESMF_KIND_I4), pointer :: indices(:,:) real(ESMF_KIND_R8), pointer :: weights(:) integer :: spherical_grid ! result code integer :: finalrc finalrc = ESMF_SUCCESS call ESMF_Initialize(vm=vm, defaultlogfilename="FieldRegridMaskEx.Log", & logkindflag=ESMF_LOGKIND_MULTI, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) call ESMF_VMGet(vm, localPet=localPet, petCount=petCount, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) !----------------------------------------------------------------------------- !NEX_Ex write(name, *) "Test GridToMesh" ! init success flag correct=.true. rc=ESMF_SUCCESS ! Establish the resolution of the grids src_nx = 100; src_ny = 100; dst_nx = 75; dst_ny = 50; ! Source mesh covers [0,1]x[0,2] src_dx = 1. / (REAL(src_nx)+1.) src_dy = 1. / (REAL(src_ny)+1.) dst_dx = 0.5 / (REAL(dst_nx)+1.) dst_dy = 0.5 / (REAL(dst_ny)+1.) ! if petCount >1, setup petMap gridSrc=ESMF_GridCreateNoPeriDim(minIndex=(/1,1/),maxIndex=(/src_nx,src_ny/),regDecomp=(/petCount,1/), & indexflag=ESMF_INDEX_GLOBAL, & rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) gridDst=ESMF_GridCreateNoPeriDim(minIndex=(/1,1/),maxIndex=(/dst_nx,dst_ny/),regDecomp=(/1,petCount/), & indexflag=ESMF_INDEX_GLOBAL, & rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) ! Create source/destination fields call ESMF_ArraySpecSet(arrayspec, 2, ESMF_TYPEKIND_R8, rc=rc) if (rc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) srcField = ESMF_FieldCreate(gridSrc, arrayspec, & staggerloc=ESMF_STAGGERLOC_CENTER, name="source", rc=localrc) if (rc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) dstField = ESMF_FieldCreate(gridDst, arrayspec, & staggerloc=ESMF_STAGGERLOC_CENTER, name="dest", rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) ! Allocate coordinates call ESMF_GridAddCoord(gridSrc, staggerloc=ESMF_STAGGERLOC_CENTER, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) call ESMF_GridAddCoord(gridDst, staggerloc=ESMF_STAGGERLOC_CENTER, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) ! Allocate items for masking call ESMF_GridAddItem(gridSrc, staggerloc=ESMF_STAGGERLOC_CENTER, & itemflag=ESMF_GRIDITEM_MASK, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) call ESMF_GridAddItem(gridDst, staggerloc=ESMF_STAGGERLOC_CENTER, & itemflag=ESMF_GRIDITEM_MASK, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) ! Get number of local DEs call ESMF_GridGet(gridSrc, localDECount=localDECount, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) ! Get arrays ! dstArray call ESMF_FieldGet(dstField, array=dstArray, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) ! srcArray call ESMF_FieldGet(srcField, array=srcArray, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) ! Get memory and set coords for src do lDE=0,localDECount-1 !! get coord 1 call ESMF_GridGetCoord(gridSrc, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=1, & computationalLBound=clbnd, computationalUBound=cubnd, farrayPtr=farrayPtrXC, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) call ESMF_GridGetCoord(gridSrc, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=2, & computationalLBound=clbnd, computationalUBound=cubnd, farrayPtr=farrayPtrYC, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) call ESMF_GridGetItem(gridSrc, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, & itemflag=ESMF_GRIDITEM_MASK, farrayPtr=maskSrc, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) call ESMF_FieldGet(srcField, lDE, farrayPtr, computationalLBound=fclbnd, & computationalUBound=fcubnd, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) if (clbnd(1) .ne. fclbnd(1)) print *, 'Error clbnd != fclbnd' if (clbnd(2) .ne. fclbnd(2)) print *, 'Error clbnd != fclbnd' if (cubnd(1) .ne. fcubnd(1)) print *, 'Error cubnd != fcubnd' if (cubnd(2) .ne. fcubnd(2)) print *, 'Error cubnd != fcubnd' !! set coords, interpolated function do i1=clbnd(1),cubnd(1) do i2=clbnd(2),cubnd(2) farrayPtrXC(i1,i2) = REAL((i1-1)*src_dx) farrayPtrYC(i1,i2) = REAL((i2-1)*src_dx) x = farrayPtrXC(i1, i2) y = farrayPtrYC(i1,i2) ! Set src mask as a circle of radius 0.5 around origin if (sqrt(x*x+y*y) < 0.5) then maskSrc(i1,i2)=1 else maskSrc(i1,i2)=0 endif ! Function farrayPtr(i1, i2) = sin(x*10*3.145)+cos(y*4*3.145) enddo enddo enddo ! lDE ! Get number of local DEs call ESMF_GridGet(gridDst, localDECount=localDECount, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) ! Get memory and set coords for dst do lDE=0,localDECount-1 !! get coord 1 call ESMF_GridGetCoord(gridDst, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=1, & computationalLBound=clbnd, computationalUBound=cubnd, farrayPtr=farrayPtrXC, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) call ESMF_GridGetCoord(gridDst, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=2, & computationalLBound=clbnd, computationalUBound=cubnd, farrayPtr=farrayPtrYC, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) call ESMF_GridGetItem(gridDst, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, & itemflag=ESMF_GRIDITEM_MASK, farrayPtr=maskDst, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) call ESMF_FieldGet(dstField, lDE, farrayPtr, computationalLBound=fclbnd, & computationalUBound=fcubnd, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) if (clbnd(1) .ne. fclbnd(1)) print *, 'Error clbnd != fclbnd' if (clbnd(2) .ne. fclbnd(2)) print *, 'Error clbnd != fclbnd' if (cubnd(1) .ne. fcubnd(1)) print *, 'Error cubnd != fcubnd' if (cubnd(2) .ne. fcubnd(2)) print *, 'Error cubnd != fcubnd' !! set coords d2rad = 0.01745329251994329547 theta = 45. ctheta = cos(theta*d2rad) stheta = sin(theta*d2rad) do i1=clbnd(1),cubnd(1) do i2=clbnd(2),cubnd(2) x = REAL((i1-1)*dst_dx) y = REAL((i2-1)*dst_dy) farrayPtrXC(i1,i2) = x-0.25 farrayPtrYC(i1,i2) = y-0.03*cos(y*3.145/0.5)*cos(x*2*3.145/0.5)-0.25 !! Set dst mask as anything .25 from y-axis if (abs(x) < 0.25) then maskDst(i1,i2)=1 else maskDst(i1,i2)=0 endif !! Now apply the transformation xtmp = farrayPtrXC(i1,i2) farrayPtrXC(i1,i2) = ctheta*farrayPtrXC(i1,i2)-stheta*farrayPtrYC(i1,i2)+0.5 farrayPtrYC(i1,i2) = stheta*xtmp+ctheta*farrayPtrYC(i1,i2)+0.5 farrayPtr(i1,i2) = 0. ! set destination field to zero enddo enddo ! Set field values enddo ! lDE !BOE ! !\subsubsection{Precompute a regridding operation using masks} ! As before, to create the sparse matrix regrid operator we call the ! {\tt ESMF\_FieldRegridStore()} routine. ! However, in this case we apply masking to the regrid operation. ! The mask value for each index location in the Grids may be set using ! the {\tt ESMF\_GridAddItem()} call (see Section~\ref{sec:usage:items} ! and Section~\ref{sec:usage:items:accessing}). Mask values may be set independantly ! for the source and destination Grids. If no mask values have been set in a Grid, then it is ! assumed no masking should be used for that Grid. The {\tt srcMaskValues} ! parameter allows the user to set the list of values which indicate ! that a source location should be masked out. The {\tt dstMaskValues} ! parameter allows the user to set the list of values which indicate ! that a destination location should be masked out. The absence of ! one of these parameters indicates that no masking should be used ! for that Field (e.g no {\tt srcMaskValue} parameter indicates that source ! masking shouldn't occur). The {\tt unmappedaction} flag may be ! used with or without masking and indicates what should occur ! if destination points can not be mapped to a source cell. ! Here the {\tt ESMF\_UNMAPPEDACTION\_IGNORE} value indicates that unmapped ! destination points are to be ignored and no sparse matrix entries should be ! generated for them. !EOE !BOC call ESMF_FieldRegridStore(srcField=srcField, srcMaskValues=(/1/), & dstField=dstField, dstMaskValues=(/1/), & unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, & routeHandle=routeHandle, & indices=indices, weights=weights, & regridmethod=ESMF_REGRIDMETHOD_BILINEAR, & rc=localrc) !EOC if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) !BOE ! ! The {\tt ESMF\_FieldRegrid} and {\tt ESMF\_FieldRegridRelease} calls ! may then be applied as in the previous example. !EOE call ESMF_FieldRegrid(srcField, dstField, routeHandle, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) call ESMF_FieldRegridRelease(routeHandle, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) #if 0 ! Write results to a mesh num_arrays = 1 spherical_grid = 0 ! Uncomment these calls to see some actual regrid results call ESMF_MeshIO(vm, GridSrc, ESMF_STAGGERLOC_CENTER, & "srcmesh", srcArray, & spherical=spherical_grid, rc=localrc) call ESMF_MeshIO(vm, Griddst, ESMF_STAGGERLOC_CENTER, & "dstmesh", dstArray, & spherical=spherical_grid, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) ! Uncomment print statement to print the weights if (associated(indices)) then do i1 = 1, size(indices,1) !print *, indices(i1,1), indices(i1,2) , ':', weights(i1) enddo endif #endif call ESMF_FieldDestroy(srcField, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) call ESMF_FieldDestroy(dstField, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) call ESMF_GridDestroy(gridSrc, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) call ESMF_GridDestroy(gridDst, rc=localrc) if (localrc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT) ! NOTE: The below should probably eventually be its own file, but it ! Seemed weird to do that now without any actual code in it. !BOE ! !\subsubsection{Regrid troubleshooting guide} ! ! The below is a list of problems users commonly encounter with regridding and potential solutions. ! This is by no means an exhaustive list, so if none of these problems fit your case, or if the solutions ! don't fix your problem, please feel free to email esmf support (esmf\_support@list.woc.noaa.gov). ! ! \bigskip ! ! {\bf Problem:} Regridding is too slow. ! ! \medskip ! ! {\bf Possible Cause:} The {\tt ESMF\_FieldRegridStore()} method is called more than is necessary. \newline ! The {\tt ESMF\_FieldRegridStore()} operation is a complex one and can be ! relatively slow for some cases (large Grids, 3D grids, etc.) ! ! \smallskip ! ! {\bf Solution:} Reduce the number of {\tt ESMF\_FieldRegridStore()} calls to the minimum necessary. The ! routeHandle generated by the {\tt ESMF\_FieldRegridStore()} call depends on only four factors: the ! stagger locations that the input Fields are created on, the coordinates in the Grids the input Fields ! are built on at those stagger locations, the padding of the input Fields ! (specified by the {\tt totalWidth} arguments in {\tt FieldCreate}) and the size of the tensor ! dimensions in the input Fields (specified by the {\tt ungridded} arguments in {\tt FieldCreate}). ! For any pair of Fields which share these attributes with the Fields used in the ! {\tt ESMF\_FieldRegridStore} call the same routeHandle can be used. Note, that the data in the ! Fields does NOT matter, the same routeHandle can be used no matter how the data in the Fields changes. ! ! \smallskip ! ! In particular: ! \begin{itemize} ! ! \item If Grid coordinates do not change during a run, then the {\tt ESMF\_FieldRegridStore()} call can be ! done once between a pair of Fields at the beginning and the resulting routeHandle used for each ! timestep during the run. ! ! \item If a pair of Fields was created with exactly the same arguments to {\tt ESMF\_FieldCreate()} as the ! pair of Fields used during an {\tt ESMF\_FieldRegridStore()} call, then the resulting routeHandle can ! also be used between that pair of Fields. ! \end{itemize} ! ! \bigskip ! ! {\bf Problem:} Distortions in destination Field at periodic boundary. ! ! \medskip ! ! {\bf Possible Cause:} The Grid overlaps itself. With a periodic Grid, the regrid system expects ! the first point to not be a repeat of the last point. In other words, ! regrid constructs its own connection and overlap between the first and last points of the ! periodic dimension and so the Grid doesn't need to contain these. If the Grid does, then this ! can cause problems. ! ! \smallskip ! ! {\bf Solution:} Define the Grid so that it doesn't contain the overlap point. This typically means simply making ! the Grid one point smaller in the periodic dimension. If a Field ! constructuted on the Grid needs to contain these overlap points then the user can use the ! {\tt totalWidth} arguments to include this extra padding in the Field. Note, however, ! that the regrid won't update these extra points, so the user will have to do a copy to fill the points ! in the overlap region in the Field. ! !EOE 10 continue call ESMF_Finalize(rc=rc) if (rc/=ESMF_SUCCESS) finalrc = ESMF_FAILURE if (finalrc==ESMF_SUCCESS) then print *, "PASS: ESMF_FieldRegridEx.F90" else print *, "FAIL: ESMF_FieldRegridEx.F90" endif end program ESMF_FieldRegridEx