DownBit.F90 Source File


This file depends on

sourcefile~~downbit.f90~~EfferentGraph sourcefile~downbit.f90 DownBit.F90 sourcefile~mapl_exceptionhandling.f90 MAPL_ExceptionHandling.F90 sourcefile~downbit.f90->sourcefile~mapl_exceptionhandling.f90 sourcefile~mapl_errorhandling.f90 MAPL_ErrorHandling.F90 sourcefile~mapl_exceptionhandling.f90->sourcefile~mapl_errorhandling.f90 sourcefile~mapl_throw.f90 MAPL_Throw.F90 sourcefile~mapl_exceptionhandling.f90->sourcefile~mapl_throw.f90 sourcefile~mapl_errorhandling.f90->sourcefile~mapl_throw.f90

Files dependent on this one

sourcefile~~downbit.f90~~AfferentGraph sourcefile~downbit.f90 DownBit.F90 sourcefile~griddedio.f90 GriddedIO.F90 sourcefile~griddedio.f90->sourcefile~downbit.f90 sourcefile~mapl_epochswathmod.f90 MAPL_EpochSwathMod.F90 sourcefile~mapl_epochswathmod.f90->sourcefile~downbit.f90 sourcefile~mapl_historygridcomp.f90 MAPL_HistoryGridComp.F90 sourcefile~mapl_historygridcomp.f90->sourcefile~downbit.f90 sourcefile~mapl_historygridcomp.f90->sourcefile~mapl_epochswathmod.f90 sourcefile~mapl_generic.f90 MAPL_Generic.F90 sourcefile~mapl_historygridcomp.f90->sourcefile~mapl_generic.f90 sourcefile~mapl_historycollection.f90 MAPL_HistoryCollection.F90 sourcefile~mapl_historygridcomp.f90->sourcefile~mapl_historycollection.f90 sourcefile~maplshared.f90 MaplShared.F90 sourcefile~maplshared.f90->sourcefile~downbit.f90 sourcefile~extdata_iobundlemod.f90 ExtData_IOBundleMod.F90 sourcefile~extdata_iobundlemod.f90->sourcefile~griddedio.f90 sourcefile~extdata_iobundlemod.f90~2 ExtData_IOBundleMod.F90 sourcefile~extdata_iobundlemod.f90~2->sourcefile~griddedio.f90 sourcefile~extdatadrivergridcomp.f90 ExtDataDriverGridComp.F90 sourcefile~extdatadrivergridcomp.f90->sourcefile~mapl_historygridcomp.f90 sourcefile~extdataroot_gridcomp.f90 ExtDataRoot_GridComp.F90 sourcefile~extdataroot_gridcomp.f90->sourcefile~maplshared.f90 sourcefile~fieldbundleread.f90 FieldBundleRead.F90 sourcefile~fieldbundleread.f90->sourcefile~griddedio.f90 sourcefile~fieldbundlewrite.f90 FieldBundleWrite.F90 sourcefile~fieldbundlewrite.f90->sourcefile~griddedio.f90 sourcefile~mapl_capgridcomp.f90 MAPL_CapGridComp.F90 sourcefile~mapl_capgridcomp.f90->sourcefile~mapl_historygridcomp.f90 sourcefile~mapl_capgridcomp.f90->sourcefile~mapl_generic.f90 sourcefile~mapl_generic.f90->sourcefile~maplshared.f90 sourcefile~mapl_historycollection.f90->sourcefile~griddedio.f90 sourcefile~mapl_historycollection.f90->sourcefile~mapl_epochswathmod.f90 sourcefile~mapl_historycollection.f90->sourcefile~mapl_generic.f90 sourcefile~capdriver.f90 CapDriver.F90 sourcefile~capdriver.f90->sourcefile~extdataroot_gridcomp.f90 sourcefile~comp_testing_driver.f90 Comp_Testing_Driver.F90 sourcefile~comp_testing_driver.f90->sourcefile~mapl_capgridcomp.f90 sourcefile~comp_testing_driver.f90->sourcefile~mapl_generic.f90 sourcefile~extdata_iobundlevectormod.f90 ExtData_IOBundleVectorMod.F90 sourcefile~extdata_iobundlevectormod.f90->sourcefile~extdata_iobundlemod.f90 sourcefile~extdata_iobundlevectormod.f90~2 ExtData_IOBundleVectorMod.F90 sourcefile~extdata_iobundlevectormod.f90~2->sourcefile~extdata_iobundlemod.f90~2 sourcefile~extdatadriver.f90 ExtDataDriver.F90 sourcefile~extdatadriver.f90->sourcefile~extdatadrivergridcomp.f90 sourcefile~extdatadriver.f90->sourcefile~extdataroot_gridcomp.f90 sourcefile~extdatadrivermod.f90 ExtDataDriverMod.F90 sourcefile~extdatadrivermod.f90->sourcefile~extdatadrivergridcomp.f90 sourcefile~extdatadrivermod.f90->sourcefile~extdataroot_gridcomp.f90 sourcefile~extdatagridcompmod.f90 ExtDataGridCompMod.F90 sourcefile~extdatagridcompmod.f90->sourcefile~extdata_iobundlemod.f90 sourcefile~extdatagridcompmod.f90->sourcefile~mapl_generic.f90 sourcefile~extdatagridcompng.f90 ExtDataGridCompNG.F90 sourcefile~extdatagridcompng.f90->sourcefile~extdata_iobundlemod.f90~2 sourcefile~extdatagridcompng.f90->sourcefile~mapl_generic.f90 sourcefile~mapl.f90 MAPL.F90 sourcefile~mapl.f90->sourcefile~fieldbundleread.f90 sourcefile~mapl.f90->sourcefile~fieldbundlewrite.f90 sourcefile~mapl.f90->sourcefile~mapl_generic.f90 sourcefile~mapl_bundleio_test.f90 mapl_bundleio_test.F90 sourcefile~mapl_bundleio_test.f90->sourcefile~fieldbundleread.f90 sourcefile~mapl_bundleio_test.f90->sourcefile~fieldbundlewrite.f90 sourcefile~mapl_cap.f90 MAPL_Cap.F90 sourcefile~mapl_cap.f90->sourcefile~mapl_capgridcomp.f90 sourcefile~mapl_geosatmaskmod.f90 MAPL_GeosatMaskMod.F90 sourcefile~mapl_geosatmaskmod.f90->sourcefile~mapl_generic.f90 sourcefile~mapl_orbgridcompmod.f90 MAPL_OrbGridCompMod.F90 sourcefile~mapl_orbgridcompmod.f90->sourcefile~mapl_generic.f90 sourcefile~mapl_stationsamplermod.f90 MAPL_StationSamplerMod.F90 sourcefile~mapl_stationsamplermod.f90->sourcefile~mapl_generic.f90 sourcefile~mapl_trajectorymod.f90 MAPL_TrajectoryMod.F90 sourcefile~mapl_trajectorymod.f90->sourcefile~mapl_generic.f90 sourcefile~regrid_util.f90 Regrid_Util.F90 sourcefile~regrid_util.f90->sourcefile~fieldbundleread.f90 sourcefile~regrid_util.f90->sourcefile~fieldbundlewrite.f90

Source Code

#include "MAPL_Generic.h"
module MAPL_DownbitMod
   use, intrinsic :: iso_c_binding, only: c_f_pointer, c_loc, c_ptr
   use mpi
   use MAPL_ExceptionHandling

   implicit none
   private

   public :: DownBit

   interface DownBit
      module procedure DownBit1D
      module procedure DownBit2D
      module procedure DownBit3D
   end interface DownBit

contains

!--------------------------------------------------------------------------
!>
! The routine `DownBit3D` returns a lower precision version of the input array
! `x` which retains `nbits_to_keep` of precision. 
! See routine `ESMF_CFIODownBit2D` or additional details. This version for
! rank 3 arrays, calls `ESMF_CFIODownBit2D()` for each vertical level.
!
!### History
!- 06Dec2006  da Silva  Initial version.
!
   subroutine DownBit3D ( x, xr, nbits_to_keep, undef, flops, mpi_comm, rc )

     implicit NONE

!
! !INPUT PARAMETERS:
!
     real, intent(in)    ::  x(:,:,:)       !! input array
     integer, intent(in) :: nbits_to_keep   !! number of bits per word to retain
                                            !! - no action if nbits_to_keep<1
     real, OPTIONAL, intent(in) :: undef    !! missing value
     logical, OPTIONAL, intent(in) :: flops !! if true, uses slower float point
                                            !!  based algorithm
     integer, optional, intent(in) :: mpi_comm
!
! !OUTPUT PARAMETERS:
!
     real, intent(out)   :: xr(:,:,:)       !! precision reduced array; can
                                            !! share storage with input array
                                            !! if it has same kind
     integer, optional, intent(out)  :: rc  !! error code
                                            !!  = 0 - all is well
                                            !! /= 0 - something went wrong
!
!------------------------------------------------------------------------------

   integer :: k

   do k = lbound(x,3), ubound(x,3)
      call DownBit2D ( x(:,:,k), xr(:,:,k), nbits_to_keep, &
                                 undef=undef, flops=flops, mpi_comm=mpi_comm, rc=rc )
   end do

   end subroutine DownBit3D

!---------------------------------------------------------------------------
!>
! This routine returns a lower precision version of the input array
! `x` which retains `nbits_to_keep` of precision. Two algorithms are
! implemented: 1) a fast one writen in C which downgrades precision
! by shifting `xbits = 24 - nbits_to_keep` bits of the mantissa, and 2) a slower
! float point based algorithm which is the same algorithm as GRIB
! with fixed number of bits packing. Notice that as in GRIB the scaling
! factor is forced to be a power of 2 rather than a generic float.
! Using this power of 2 binary scaling has the advantage of improving
! the GZIP compression rates.
!
! This routine returns an array of the same type and kind as the input array,
! so no data compression has taken place. The goal here is to reduce the
! entropy in the input array, thereby improving compression rates
! by the lossless algorithms implemented internally by HDF-4/5 when writing
! these data to a file. In fact, these GZIP'ed and pre-conditioned files
! have sizes comparable to the equivalent GRIB file, while being a bonafide
! self-describing HDF/NetCDF file.
!
! @todo
! Perhaps implement GRIB decimal scaling (variable number of bits).
!@endtodo
!
!#### History
!- 06Dec2006  da Silva  Initial version.
! 
   subroutine DownBit2D ( x, xr, nbits_to_keep, undef, flops, mpi_comm, rc )

     implicit NONE

!
! !INPUT PARAMETERS:
!
     real, intent(in)    ::  x(:,:)         !! input array
     integer, intent(in) :: nbits_to_keep   !! number of bits per word to retain
     real, OPTIONAL, intent(in) :: undef    !! missing value
     logical, OPTIONAL, intent(in) :: flops !! if true, uses slower float point
                                            !!  based algorithm
     integer, optional, intent(in) :: mpi_comm
!
! !OUTPUT PARAMETERS:
!
     real, intent(out)   :: xr(:,:)          !! precision reduced array; can
                                             !!  share storage with input array
                                             !!  if it has same kind
     integer, optional, intent(out)  :: rc   !! error code
                                             !!  = 0 - all is well
                                             !! /= 0 - something went wrong
!
!------------------------------------------------------------------------------
    integer   :: E, xbits, has_undef, passed_minmax
    real    :: scale, xmin, xmax, tol, undef_
    logical   :: shave_mantissa
    integer, external :: MAPL_ShaveMantissa32
    real :: min_value, max_value
    integer :: useable_mpi_comm,status

    rc = 0

    if (present(mpi_comm)) then
       useable_mpi_comm = mpi_comm
    else
       useable_mpi_comm = MPI_COMM_NULL
    end if
!   Defaults for optinal arguments
!   ------------------------------
    if ( present(undef) ) then
         undef_ = undef
         has_undef = 1
    else
         undef_ = 1.0
         undef_ = huge(undef_)   ! why not?
         has_undef = 0
    endif
    if ( present(flops) ) then
         shave_mantissa = .not. flops
    else
         shave_mantissa = .true.
    endif

!   Fast, bit shifting in C
!   -----------------------
    if ( shave_mantissa ) then

       xr = x   ! compiled r8 this will convert to r4.
       xbits = 24 - nbits_to_keep
       call compute_min_max(xr,min_value,max_value,undef_,useable_mpi_comm,_RC)
       if (useable_mpi_comm/=MPI_COMM_NULL) passed_minmax = 1
       rc = MAPL_ShaveMantissa32 ( xr, xr, size(x), xbits, has_undef, undef_, size(x), passed_minmax, min_value, max_value )
       return

!   Slow, flops in FORTRAN (GRIB inspired)
!   --------------------------------------
    else

       if ( nbits_to_keep < 1 ) then
          xr = x
          rc = 1
          return
       end if

       tol = 0.0001 * undef_
       xmin = minval(x,mask=(abs(undef_-x)>tol))
       xr = x - xmin     ! As in GRIB, force non-negative values
       xmax = maxval(xr,mask=(abs(undef_-x)>tol)) ! max of positive

       if ( xmax <= 0.0 ) then
            xr = x
            rc = 0
            return  ! this means field is constant
       end if

       E = nint(log(xmax)/log(2.)) - nbits_to_keep ! GRIB binary scale factor
       scale = 2.**E                       ! GRIB requires power of 2

       if ( present(undef) ) then
          where ( abs(x - undef_) > tol )
             xr  = xmin + nint(xr/scale) * scale
          endwhere
       else
          xr  = xmin + nint(xr/scale) * scale
       end if

    end if

   end subroutine DownBit2D

   subroutine DownBit1D ( x, xr, nbits_to_keep, undef, flops, mpi_comm, rc )
     implicit NONE
!
! !INPUT PARAMETERS:
!
     real,target, intent(in)    ::  x(:)         ! input array
     integer, intent(in) :: nbits_to_keep           ! number of bits per word to retain
     real, OPTIONAL, intent(in) :: undef    ! missing value
     logical, OPTIONAL, intent(in) :: flops ! if true, uses slower float point
                                            !  based algorithm
     integer, optional, intent(in) :: mpi_comm
!
! !OUTPUT PARAMETERS:
!
     real, target, intent(inout)   :: xr(:)   ! precision reduced array; can
     integer, optional, intent(out) :: rc


     real, pointer :: x_tmp(:,:)
     real, pointer :: xr_tmp(:,:)
     type(c_ptr) :: x_ptr, xr_ptr

     x_ptr = c_loc(x(1))
     call c_f_pointer(x_ptr,  x_tmp,[size(x),1])
     xr_ptr = c_loc(xr(1))
     call c_f_pointer(xr_ptr, xr_tmp,[size(x),1])

     call Downbit2d(x_tmp(:,:), xr_tmp(:,:), nbits_to_keep, undef=undef, flops=flops, mpi_comm=mpi_comm, rc=rc)

   end subroutine Downbit1D

   subroutine compute_min_max(array,min_value,max_value,undef_value,mpi_comm,rc)
      real, intent(in)  :: array(:,:)
      real, intent(out) :: min_value
      real, intent(out) :: max_value
      real, intent(in ) :: undef_value
      integer, intent(in ) :: mpi_comm
      integer, optional, intent(out) :: rc

      real :: local_min_value, local_max_value
      logical, allocatable :: undef_mask(:,:)
      integer :: status

      allocate(undef_mask(size(array,1),size(array,2)))
      undef_mask = .false.
      where(array /= undef_value) undef_mask = .true.

      local_min_value = minval(array,mask=undef_mask)
      local_max_value = maxval(array,mask=undef_mask)
      if (mpi_comm /= MPI_COMM_NULL) then
         call MPI_AllReduce(local_min_value,min_value,1,MPI_FLOAT,MPI_MIN,mpi_comm,status)
         _VERIFY(status)
         call MPI_AllReduce(local_max_value,max_value,1,MPI_FLOAT,MPI_MAX,mpi_comm,status)
         _VERIFY(status)
      else
         min_value = local_min_value
         max_value = local_max_value
      end if
      _RETURN(_SUCCESS)
   end subroutine 

end module MAPL_DownbitMod