source: LMDZ6/trunk/libf/phylmd/ecrad/ifs/easy_netcdf_read_mpi.F90 @ 5450

! easy_netcdf_read_mpi.f90 - Read netcdf file on one task and share with other tasks
!
! (C) Copyright 2017- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.
!
! Author:  Robin Hogan
! Email:   r.j.hogan@ecmwf.int
!

module easy_netcdf_read_mpi

  use easy_netcdf,   only : netcdf_file_raw => netcdf_file
  use parkind1,      only : jpim, jprb
  use radiation_io,  only : nulout, nulerr, my_abort => radiation_abort

  implicit none

  ! MPI tag for radiation and physics communication
  integer(kind=jpim), parameter :: mtagrad = 2800

  !---------------------------------------------------------------------
  ! An object of this type provides convenient read or write access to
  ! a NetCDF file
  type netcdf_file
    type(netcdf_file_raw) :: file
    logical :: is_master_task = .true.
  contains
    procedure :: open => open_netcdf_file
    procedure :: close => close_netcdf_file
    procedure :: get_real_scalar
    procedure :: get_real_vector
    procedure :: get_real_matrix
    procedure :: get_real_array3
    generic   :: get => get_real_scalar, get_real_vector, &
         &              get_real_matrix, get_real_array3
    procedure :: get_global_attribute

    procedure :: set_verbose
    procedure :: transpose_matrices
    procedure :: exists
  end type netcdf_file

contains

  ! --- GENERIC SUBROUTINES ---

  !---------------------------------------------------------------------
  ! Open a NetCDF file with name "file_name", optionally specifying the
  ! verbosity level (0-5)
  subroutine open_netcdf_file(this, file_name, iverbose)

    USE MPL_MODULE, ONLY : MPL_BROADCAST, MPL_MYRANK

    class(netcdf_file)            :: this
    character(len=*), intent(in)  :: file_name
    integer, intent(in), optional :: iverbose

    integer                       :: istatus

    ! Store verbosity level in object
    if (present(iverbose)) then
      this%file%iverbose = iverbose
    else
      ! By default announce files being opened and closed, but not
      ! variables read/written
      this%file%iverbose = 2
    end if

    ! By default we don't transpose 2D arrays on read
    this%file%do_transpose_2d = .false.

    if (MPL_MYRANK() == 1) then
      this%is_master_task = .true.
      call this%file%open(file_name, iverbose)
    else
      this%is_master_task = .false.
    end if

  end subroutine open_netcdf_file


  !---------------------------------------------------------------------
  ! Close the NetCDF file
  subroutine close_netcdf_file(this)
    class(netcdf_file) :: this
    integer            :: istatus

    if (this%is_master_task) then
      call this%file%close()
    end if

  end subroutine close_netcdf_file


  !---------------------------------------------------------------------
  ! Set the verbosity level from 0 to 5, where the codes have the
  ! following meaning: 0=errors only, 1=warning, 2=info, 3=progress,
  ! 4=detailed, 5=debug
  subroutine set_verbose(this, ival)
    class(netcdf_file) :: this
    integer, optional  :: ival

    if (present(ival)) then
      this%file%iverbose = ival
    else
      this%file%iverbose = 2
    end if

  end subroutine set_verbose



  !---------------------------------------------------------------------
  ! Specify whether 2D arrays should be transposed on read
  subroutine transpose_matrices(this, do_transpose)
    class(netcdf_file) :: this
    logical, optional  :: do_transpose

    if (present(do_transpose)) then
      this%file%do_transpose_2d = do_transpose
    else
      this%file%do_transpose_2d = .true.
    end if

  end subroutine transpose_matrices



  ! --- READING SUBROUTINES ---

  !---------------------------------------------------------------------
  ! Return true if the variable is present, false otherwise
  function exists(this, var_name) result(is_present)

    USE MPL_MODULE, ONLY : MPL_BROADCAST, MPL_NPROC

    class(netcdf_file)           :: this
    character(len=*), intent(in) :: var_name

    logical :: is_present

    if (this%is_master_task) then
      is_present = this%file%exists(var_name)
    end if

    if (MPL_NPROC() > 1) then
      CALL MPL_BROADCAST(is_present, mtagrad, 1, &
           &  CDSTRING='EASY_NETCDF_READ_MPI:EXISTS')
    end if

  end function exists


  !---------------------------------------------------------------------
  ! The method "get" will read either a scalar, vector or matrix
  ! depending on the rank of the output argument. This version reads a
  ! scalar.
  subroutine get_real_scalar(this, var_name, scalar)

    USE MPL_MODULE, ONLY : MPL_BROADCAST, MPL_NPROC

    class(netcdf_file)           :: this
    character(len=*), intent(in) :: var_name
    real(jprb), intent(out)      :: scalar

    if (this%is_master_task) then
      call this%file%get(var_name, scalar)
    end if

    if (MPL_NPROC() > 1) then
      CALL MPL_BROADCAST(scalar, mtagrad, 1, &
           &  CDSTRING='EASY_NETCDF_READ_MPI:GET_REAL_SCALAR')
    end if

  end subroutine get_real_scalar


  !---------------------------------------------------------------------
  ! Read a 1D array into "vector", which must be allocatable and will
  ! be reallocated if necessary
  subroutine get_real_vector(this, var_name, vector)

    USE MPL_MODULE, ONLY : MPL_BROADCAST, MPL_NPROC

    class(netcdf_file)           :: this
    character(len=*), intent(in) :: var_name
    real(jprb), allocatable, intent(out) :: vector(:)

    integer                      :: n  ! Length of vector

    n = 0

    if (this%is_master_task) then
      call this%file%get(var_name, vector)
      n = size(vector)
    end if

    if (MPL_NPROC() > 1) then
      CALL MPL_BROADCAST(n, mtagrad, 1, &
           &  CDSTRING='EASY_NETCDF_READ_MPI:GET_REAL_VECTOR:SIZE')

      if (.not. this%is_master_task) then
        allocate(vector(n))
      end if

      CALL MPL_BROADCAST(vector, mtagrad, 1, &
           &  CDSTRING='EASY_NETCDF_READ_MPI:GET_REAL_VECTOR')
    end if

  end subroutine get_real_vector


  !---------------------------------------------------------------------
  ! Read 2D array into "matrix", which must be allocatable and will be
  ! reallocated if necessary.  Whether to transpose is specifed by the
  ! final optional argument, but can also be specified by the
  ! do_transpose_2d class data member.
  subroutine get_real_matrix(this, var_name, matrix, do_transp)

    USE MPL_MODULE, ONLY : MPL_BROADCAST, MPL_NPROC

    class(netcdf_file)           :: this
    character(len=*), intent(in) :: var_name
    real(jprb), allocatable, intent(out) :: matrix(:,:)
    logical, optional, intent(in):: do_transp ! Transpose data?

    integer                      :: n(2)

    n = 0

    if (this%is_master_task) then
      call this%file%get(var_name, matrix, do_transp)
      n = shape(matrix)
    end if

    if (MPL_NPROC() > 1) then
      CALL MPL_BROADCAST(n, mtagrad, 1, &
           &  CDSTRING='EASY_NETCDF_READ_MPI:GET_REAL_MATRIX:SIZE')

      if (.not. this%is_master_task) then
        allocate(matrix(n(1),n(2)))
      end if

      CALL MPL_BROADCAST(matrix, mtagrad, 1, &
           &  CDSTRING='EASY_NETCDF_READ_MPI:GET_REAL_MATRIX')
    end if

  end subroutine get_real_matrix


  !---------------------------------------------------------------------
  ! Read 3D array into "var", which must be allocatable and will be
  ! reallocated if necessary.  Whether to pemute is specifed by the
  ! final optional argument
  subroutine get_real_array3(this, var_name, var, ipermute)

    USE MPL_MODULE, ONLY : MPL_BROADCAST, MPL_NPROC

    class(netcdf_file)                   :: this
    character(len=*), intent(in)         :: var_name
    real(jprb), allocatable, intent(out) :: var(:,:,:)
    integer, optional, intent(in)        :: ipermute(3)

    integer                              :: n(3)

    n = 0

    if (this%is_master_task) then
      call this%file%get(var_name, var, ipermute)
      n = shape(var)
    end if

    if (MPL_NPROC() > 1) then
      CALL MPL_BROADCAST(n, mtagrad, 1, &
           &  CDSTRING='EASY_NETCDF_READ_MPI:GET_REAL_ARRAY3:SIZE')

      if (.not. this%is_master_task) then
        allocate(var(n(1),n(2),n(3)))
      end if

      CALL MPL_BROADCAST(var, mtagrad, 1, &
           &  CDSTRING='EASY_NETCDF_READ_MPI:GET_REAL_ARRAY3')
    end if

  end subroutine get_real_array3


  !---------------------------------------------------------------------
  ! Get a global attribute as a character string
  subroutine get_global_attribute(this, attr_name, attr_str)

    USE MPL_MODULE, ONLY : MPL_BROADCAST, MPL_NPROC

    class(netcdf_file) :: this

    character(len=*), intent(in)    :: attr_name
    character(len=*), intent(inout) :: attr_str

    if (this%is_master_task) then
      call this%file%get_global_attribute(attr_name, attr_str)
    end if

    if (MPL_NPROC() > 1) then
      CALL MPL_BROADCAST(attr_str, mtagrad, 1, &
           &  CDSTRING='EASY_NETCDF_READ_MPI:GET_GLOBAL_ATTRIBUTE')
    end if

  end subroutine get_global_attribute

end module easy_netcdf_read_mpi