! $Id$
module open_climoz_m

  implicit none

contains

  subroutine open_climoz(ncid, press_in_cen, press_in_edg, time_in)

    ! This procedure should be called once per "gcm" run, by a single
    ! thread of each MPI process.
    ! The root MPI process opens "climoz_LMDZ.nc", reads the pressure
    ! levels and the times and broadcasts them to the other processes.

    ! We assume that, in "climoz_LMDZ.nc", the pressure levels are in hPa
    ! and in strictly ascending order.

    use netcdf95, only: nf95_open, nf95_close, nf95_gw_var, nf95_inq_varid
    use netcdf, only: nf90_nowrite

    use mod_phys_lmdz_mpi_data, only: is_mpi_root
    use mod_phys_lmdz_mpi_transfert, only: bcast_mpi ! broadcast
    USE phys_cal_mod, ONLY: calend, year_len, year_cur

    ! OpenMP shares arguments:
    integer, intent(out):: ncid ! of "climoz_LMDZ.nc", OpenMP shared

    ! pressure levels for ozone climatology, in Pa, strictly ascending order
    real, pointer:: press_in_cen(:) !--- at cells centers
    real, pointer:: press_in_edg(:) !--- at the interfaces (pressure intervals)

    real, pointer:: time_in(:)
    ! times of ozone climatology records, in days since Jan. 1st 0h

    ! Variables local to the procedure:

    integer varid ! for NetCDF
    integer n_plev ! number of pressure levels in the input data
    integer n_time ! number of time records    in the input field
    integer k
    CHARACTER(LEN=80)  :: modname
    CHARACTER(LEN=320) :: msg

    !---------------------------------------

    modname="open_climoz"
    print *, "Call sequence information: "//TRIM(modname)

    if (is_mpi_root) then
       call nf95_open("climoz_LMDZ.nc", nf90_nowrite, ncid)

       call nf95_inq_varid(ncid, "plev", varid)
       call nf95_gw_var(ncid, varid, press_in_cen)
       ! Convert from hPa to Pa because "paprs" and "pplay" are in Pa:
       press_in_cen = press_in_cen * 100.
       n_plev = size(press_in_cen)
       CALL NF95_INQ_VARID(ncID, "time", varID)
       CALL NF95_GW_VAR(ncid, varid, time_in)
       n_time = SIZE(time_in)
       IF(ALL([year_len,14]/=n_time)) THEN             !--- Check records number
         WRITE(msg,'(a,3(i4,a))')TRIM(modname)//': expected input file could be&
         & monthly or daily (14 or ',year_len,' records for year ',year_cur,' a&
         &nd calendar "'//TRIM(calend)//'"), but ',n_time,' records were found'
         CALL abort_physic(modname, msg, 1)
       END IF

    end if

    call bcast_mpi(n_plev)
    if (.not. is_mpi_root) allocate(press_in_cen(n_plev))
    call bcast_mpi(press_in_cen)
    call bcast_mpi(n_time)
    if (.not. is_mpi_root) allocate(time_in(n_time))
    call bcast_mpi(time_in)

    ! Compute edges of pressure intervals:
    allocate(press_in_edg(n_plev + 1))
    if (is_mpi_root) then
       press_in_edg(1) = 0.
       ! We choose edges halfway in logarithm:
       forall (k = 2:n_plev) press_in_edg(k) = sqrt(press_in_cen(k - 1) * press_in_cen(k))
       press_in_edg(n_plev + 1) = huge(0.)
       ! (infinity, but any value guaranteed to be greater than the
       ! surface pressure would do)
    end if
    call bcast_mpi(press_in_edg)
!    deallocate(press_in_cen) ! pointer

  end subroutine open_climoz

end module open_climoz_m