Program utCFIO
use ESMF
use MAPL_BaseMod
use MAPL_CommsMod
use MAPL_Constants
use ESMF_CfioMod
use MAPL_CfioMod
implicit NONE
type(ESMF_Grid) :: grid2d, grid3d
type (ESMF_VM) :: VM
type(ESMF_DELayout) :: layout
integer :: nymd, nhms
type(ESMF_Time) :: time2, time3, times
type(ESMF_Time) :: time2a, time2b, time2c
type(ESMF_Time) :: time3a, time3b, time3c
type(ESMF_CFIO) :: cfio
integer, parameter :: IM_WORLD = 288, JM_WORLD = 181, KM_WORLD = 72 ! global
integer :: i, j, k, im, jm, km, n ! local
character(len=*), parameter :: &
dirname = '/share/dasilva/fvInput/fvchem/b/', &
f2d_Filename = dirname // 'bian_biogenic_ch4_20010112.2x2.5.nc', &
f3d_Filename = dirname // 'bian_OHCH4_2x25x72.019498.nc', &
src_Filename = dirname // 'gocart.du_src.sfc.1971.hdf'
real, pointer, dimension(:,:) :: emcoterp, emconvoc
real, pointer, dimension(:,:,:) :: oh, ch4, xa, xb, xc, xm, xe
real, pointer :: lons(:), lats(:)
integer :: status, rc, dims(3)
logical :: IamRoot
character(len=*), parameter :: Iam = 'utCFIO'
! -----
call test_main()
CONTAINS
subroutine test_main()
! Initialize framework
! --------------------
call ESMF_Initialize (vm=vm, rc=status)
_VERIFY(status)
IamRoot = MAPL_am_I_root()
! Get the global vm
! -----------------
call ESMF_VMGetGlobal(vm, rc=status)
_VERIFY(status)
! Create a grid
! -------------
call MyGridCreate_ ( vm, grid2d, grid3d, rc=status )
_VERIFY(status)
! Set the time as the one on the hardwired file name
! --------------------------------------------------
call ESMF_CalendarSetDefault ( ESMF_CALKIND_GREGORIAN, rc=status )
_VERIFY(STATUS)
call ESMF_TimeSet(Time2, yy=2001, mm=1, dd=1, h=0, m=0, s=0, rc=status )
_VERIFY(STATUS)
call ESMF_TimeSet(Time3, yy=2001, mm=1, dd=14, h=12, m=0, s=0, rc=status )
_VERIFY(STATUS)
call ESMF_TimeSet(Time2a, yy=2001, mm=3, dd=1, h=12, m=0, s=0, rc=status )
_VERIFY(STATUS)
call ESMF_TimeSet(Times, yy=1971, mm=6, dd=5, h=0, m=0, s=0, rc=status )
_VERIFY(STATUS)
! Get Local dimensions
! --------------------
call MAPL_GridGet ( GRID3D, &
localCellCountPerDim=DIMS,RC=STATUS)
_VERIFY(STATUS)
! Allocate arrays
! ---------------
im = dims(1); jm = dims(2); km = dims(3)
allocate( emcoterp(im,jm), emconvoc(im,jm), &
oh(im,jm,km), ch4(im,jm,km), &
stat=STATUS )
_VERIFY(STATUS)
! Read 2D arrays from file
! ------------------------
! if ( IamRoot ) print *, 'Reading ' // fFilename
call ESMF_ioRead ( 'du_src', src_Filename, Times, grid3d, emcoterp, &
verbose=.true., force_regrid=.true., &
time_is_cyclic = .true., rc=STATUS )
_VERIFY(STATUS)
call ESMF_ioRead ( 'emcoterp', f2d_Filename, Time2, grid3d, emcoterp, &
verbose=.true., force_regrid=.true., &
time_is_cyclic = .true., rc=STATUS )
_VERIFY(STATUS)
call ESMF_ioRead ( 'emconvoc', f2d_Filename, Time2, grid3d, emconvoc, &
verbose=.true., force_regrid=.true., &
time_is_cyclic = .true., rc=STATUS )
_VERIFY(STATUS)
call ESMF_ioRead ( 'OH', f3d_Filename, Time3, grid3d, oh, &
verbose=.true., force_regrid=.true., rc=STATUS )
_VERIFY(STATUS)
call ESMF_ioRead ( 'CH4', f3d_Filename, Time3, grid3d, ch4, &
verbose=.true., force_regrid=.true., rc=STATUS )
_VERIFY(STATUS)
if ( IamRoot ) print *, 'Re-reading w/ time interpolation...'
call ESMF_ioRead ( 'emconvoc', f2d_Filename, Time2a, grid3d, emconvoc, &
verbose=.true., force_regrid=.true., &
time_is_cyclic = .true., time_interp=.true., &
rc=STATUS )
_VERIFY(STATUS)
deallocate(oh,ch4,emconvoc)
allocate( xa(im,jm,km), xb(im,jm,km), xc(im,jm,km), &
xm(im,jm,km), xe(im,jm,km), &
stat=STATUS )
_VERIFY(STATUS)
print *
! Testing time interpolation
! --------------------------
call ESMF_TimeSet(Time3a, yy=2001, mm=3, dd=16, h=12, m=0, s=0, rc=status )
call ESMF_TimeSet(Time3b, yy=2001, mm=3, dd=31, h=18, m=0, s=0, rc=status )
call ESMF_TimeSet(Time3c, yy=2001, mm=4, dd=16, h= 0, m=0, s=0, rc=status )
call ESMF_ioRead ( 'CH4', f3d_Filename, Time3a, grid3d, xa, &
verbose=.true., force_regrid=.true., &
time_interp = .true., &
rc=STATUS )
_VERIFY(STATUS)
call ESMF_ioRead ( 'CH4', f3d_Filename, Time3b, grid3d, xb, &
verbose=.true., force_regrid=.true., rc=STATUS, &
time_interp = .true. )
_VERIFY(STATUS)
call ESMF_ioRead ( 'CH4', f3d_Filename, Time3c, grid3d, xc, &
verbose=.true., force_regrid=.true., rc=STATUS, &
time_interp = .true. )
_VERIFY(STATUS)
xm = ( xa + xc ) / 2.
xe = abs ( xb - xm )
n = im * jm
i = im/2
if ( IamRoot ) then
print *
print *, 'Time interpolation stats: '
print *
do j = 1, jm, jm/4
print *, 'Lat = ', lats(j)
print *, 'CH4 at time a: ', (xa(i,j,k), k=1,km,km/4)
print *, 'CH4 at time b: ', (xb(i,j,k), k=1,km,km/4)
print *, 'CH4 m e a n : ', (xm(i,j,k), k=1,km,km/4)
print *, 'CH4 at time c: ', (xc(i,j,k), k=1,km,km/4)
print *, 'CH4 e r r o r: ', (xe(i,j,k), k=1,km,km/4)
print *
end do
print *, 'CH4 RMS error: ', sqrt(sum(xe*xe)/n)
endif
! All done
! --------
call ESMF_Finalize ( rc=status )
_VERIFY(STATUS)
end subroutine test_main
!........................................................................
subroutine MyGridCreate_ ( vm, grid2d, grid3d, rc)
type (ESMF_VM), intent(IN ) :: VM
type(ESMF_GRID), intent(INOUT) :: grid2d, grid3d
integer, optional, intent(OUT) :: rc
! Local vars
integer :: status
character(len=ESMF_MAXSTR), parameter :: IAm='MyGridCreate'
integer :: LM
integer :: L
integer :: NX, NY
integer, allocatable :: IMXY(:), JMXY(:)
character(len=ESMF_MAXSTR) :: gridname
real(ESMF_KIND_R8) :: minCoord(3)
real(ESMF_KIND_R8) :: deltaX, deltaY, deltaZ
real :: LON0, LAT0
real :: pi, d2r
! grid create
lm = KM_WORLD ! no. vertical layers
nx = 2
ny = 2
pi = 4.0 * atan ( 1.0 )
d2r = pi / 180.
LON0 = -180 * d2r
LAT0 = -90.0 * d2r
! Get the IMXY vector
! -------------------
allocate( imxy(0:nx-1) )
call MAPL_GET_LOCAL_DIMS ( IM_WORLD, imxy, nx )
! Get the JMXY vector
! -------------------
allocate( jmxy(0:ny-1) )
call MAPL_GET_LOCAL_DIMS ( JM_WORLD, jmxy, ny )
deltaX = 2.0*pi/IM_WORLD
deltaY = pi/(JM_WORLD-1)
deltaZ = 1.0
! Define South-West Corner of First Grid-Box
! ------------------------------------------
minCoord(1) = LON0 - deltaX/2
minCoord(2) = LAT0 - deltaY/2
minCoord(3) = deltaZ/2.
layout = ESMF_DELayoutCreate(vm, deCountList=(/NX, NY/), rc=status)
_VERIFY(STATUS)
if ( MAPL_Am_I_Root() ) then
print *
print *, ' Testing CFIO Read Array'
print *, ' -----------------------'
print *
print *, 'Decomposition Information:'
print *, ' X AXIS: nx, imxy = ', nx, imxy
print *, ' Y AXIS: ny, jmxy = ', ny, jmxy
print *, ' Z AXIS: nz = ', lm, ' (not decomposed)'
end if
! 2D grid
! -------
grid2d = ESMF_GridCreateHorzLatLonUni( &
counts = (/IM_WORLD, JM_WORLD/), &
minGlobalCoordPerDim=minCoord(1:2), &
deltaPerDim=(/deltaX, deltaY /), &
horzStagger=ESMF_Grid_Horz_Stagger_A, &
periodic=(/ESMF_TRUE, ESMF_FALSE/), &
name='Grid2d', rc=status)
_VERIFY(STATUS)
call ESMF_GridDistribute(grid2d, &
deLayout=layout, &
countsPerDEDim1=imxy, &
countsPerDEDim2=jmxy, &
rc=status)
_VERIFY(STATUS)
! 3D Grid
! -------
grid3d = ESMF_GridCreateHorzLatLonUni( &
counts = (/IM_WORLD, JM_WORLD/), &
minGlobalCoordPerDim=minCoord(1:2), &
deltaPerDim=(/deltaX, deltaY /), &
horzStagger=ESMF_Grid_Horz_Stagger_A, &
periodic=(/ESMF_TRUE, ESMF_FALSE/), &
name='Grid3d', rc=status)
_VERIFY(STATUS)
call ESMF_GridAddVertHeight(grid3d, &
delta=(/(deltaZ, L=1,LM) /), &
rc=status)
_VERIFY(STATUS)
call ESMF_GridDistribute(grid3d, &
deLayout=layout, &
countsPerDEDim1=imxy, &
countsPerDEDim2=jmxy, &
rc=status)
_VERIFY(STATUS)
call MAPL_GridGetLatLons ( grid3d, lons, lats )
if ( MAPL_Am_I_Root() ) then
print *
print *, 'Grid points:'
print *, ' - Longitudes: '
write(*,100) lons(:)
print *, ' - Latitudes:'
write(*,100) lats(:)
print *
100 format(( ' ',8F8.2))
endif
deallocate(imxy)
deallocate(jmxy)
_RETURN(STATUS)
end subroutine MyGridCreate_
subroutine MAPL_GridGetLatLons ( grid, lons, lats )
implicit NONE
type(ESMF_Grid) :: grid
real, pointer :: lons(:), lats(:)
! ---
type(ESMF_Array) :: eARRAY(2)
real(KIND=8), pointer :: R8D2(:,:)
real, pointer :: lons2d(:,:), lats2d(:,:)
real, pointer :: LONSLocal(:,:), LATSlocal(:,:)
integer :: IM_WORLD, JM_WORLD, dims(3)
! ----
! Get world dimensions
! --------------------
call ESMF_GridGet ( grid, horzRelloc=ESMF_CELL_CENTER, &
vertRelLoc=ESMF_CELL_CENTER, &
globalCellCountPerDim=DIMS, RC=STATUS)
_VERIFY(STATUS)
IM_WORLD = dims(1)
JM_WORLD = dims(2)
! Allocate memory for output if necessary
! ---------------------------------------
if ( .not. associated(lons) ) then
allocate(lons(IM_WORLD), stat=STATUS)
else
if(size(LONS,1) /= IM_WORLD) STATUS = 1
end if
_VERIFY(status)
if ( .not. associated(lats) ) then
allocate(lats(JM_WORLD), stat=STATUS)
else
if(size(LATS,1) /= JM_WORLD) STATUS = 1
end if
_VERIFY(status)
! Retrieve the ESMF array with coordinates
! ----------------------------------------
call ESMF_GridGetCoord ( grid, horzRelLoc =ESMF_CELL_CENTER, &
centerCoord=eARRAY, RC=STATUS )
_VERIFY(STATUS)
! Local work space
! ----------------
allocate(LONS2d(IM_WORLD,JM_WORLD), LATS2d(IM_WORLD,JM_WORLD), &
STAT=status)
_VERIFY(status)
! Get the local longitudes and gather them into a global array
! ------------------------------------------------------------
call ESMF_ArrayGetData(EARRAY(1), R8D2, RC=STATUS)
_VERIFY(STATUS)
allocate(LONSLOCAL(size(R8D2,1),size(R8D2,2)), STAT=status)
_VERIFY(status)
LONSLOCAL = R8D2*(180/MAPL_PI)
call ArrayGather(LONSLOCAL, LONS2D, GRID, RC=STATUS)
! Get the local longitudes and gather them into a global array
! ------------------------------------------------------------
call ESMF_ArrayGetData(eARRAY(2), R8D2, RC=STATUS)
_VERIFY(STATUS)
allocate(LATSLOCAL(size(R8D2,1),size(R8D2,2)), STAT=status)
_VERIFY(status)
LATSlocal = R8D2*(180/MAPL_PI)
call ArrayGather(LATSLOCAL, LATS2D, GRID, RC=STATUS)
_VERIFY(STATUS)
! Return 1D arrays
! ----------------
LONS = LONS2D(:,1)
LATS = LATS2D(1,:)
DEALLOCATE(LONSLOCAL, LATSLOCAL, LONS2d, LATS2d )
end subroutine MAPL_GridGetLatLons
end Program utCFIO