utCFIO_Array.F90 Source File


This file depends on

sourcefile~~utcfio_array.f90~~EfferentGraph sourcefile~utcfio_array.f90 utCFIO_Array.F90 sourcefile~base_base.f90 Base_Base.F90 sourcefile~utcfio_array.f90->sourcefile~base_base.f90 sourcefile~constants.f90 Constants.F90 sourcefile~utcfio_array.f90->sourcefile~constants.f90 sourcefile~esmf_cfiomod.f90 ESMF_CFIOMod.F90 sourcefile~utcfio_array.f90->sourcefile~esmf_cfiomod.f90 sourcefile~mapl_cfio.f90 MAPL_CFIO.F90 sourcefile~utcfio_array.f90->sourcefile~mapl_cfio.f90 sourcefile~mapl_comms.f90 MAPL_Comms.F90 sourcefile~utcfio_array.f90->sourcefile~mapl_comms.f90

Source Code

!
! Simple unit test for CFIO Read/Write of Arrays
!

#include "MAPL_Generic.h"

   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