MODULE regr_horiz_time_climoz_m

  USE interpolation, ONLY: locate
  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, grid_type, unstructured
  USE nrtype, ONLY: pi
  USE netcdf, ONLY: nf90_clobber, nf90_float, nf90_open, &
          nf90_nowrite, nf90_noerr, nf90_get_att, nf90_global
  USE netcdf95, ONLY: nf95_def_dim, nf95_inq_dimid, nf95_inquire_dimension, &
          nf95_def_var, nf95_inq_varid, nf95_inquire_variable, &
          nf95_open, nf95_create, nf95_get_att, nf95_gw_var, nf95_get_var, &
          nf95_close, nf95_enddef, nf95_put_att, nf95_put_var, nf95_copy_att
  USE print_control_mod, ONLY: lunout
  USE dimphy
  IMPLICIT NONE
  PRIVATE
  PUBLIC :: regr_horiz_time_climoz
  REAL, PARAMETER :: deg2rad = pi / 180.
  CHARACTER(LEN = 13), PARAMETER :: vars_in(2) = ['tro3         ', 'tro3_daylight']

  INTEGER :: nlat_ou, nlon_ou
  REAL, ALLOCATABLE :: latitude_glo(:)
  !$OMP THREADPRIVATE(latitude_glo)
  INTEGER, ALLOCATABLE :: ind_cell_glo_glo(:)
  !$OMP THREADPRIVATE(ind_cell_glo_glo)

CONTAINS

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

  SUBROUTINE regr_horiz_time_climoz(read_climoz, interpt)

    !-------------------------------------------------------------------------------
    ! Purpose: Regrid horizontally and in time zonal or 3D ozone climatologies.
    !   * Read ozone climatology from netcdf file
    !   * Regrid it horizontaly to LMDZ grid (quasi-conservative method)
    !   * If interpt=T, interpolate linearly in time (one record each day)
    !     If interpt=F, keep original time sampling  (14 months).
    !   * Save it to a new netcdf file.
    !-------------------------------------------------------------------------------
    ! Remarks:
    !   * Up to 2 variables treated: "tro3" and "tro3_daylight" (if read_climoz=2)
    !   * Input fields coordinates: (longitudes, latitudes, pressure_levels, time)
    !   * Output grid cells centers coordinates given by [rlonv,] rlatu.
    !   * Output grid cells edges   coordinates given by [rlonu,] rlatv.
    !   * Input file [longitudes and] latitudes given in degrees.
    !   * Input file pressure levels are given in Pa or hPa.
    !   * All coordinates variables are stricly monotonic.
    !   * Monthly fields are interpolated linearly in time to get daily values.
    !   * Fields are known at the middle of the months, so interpolation requires an
    !     additional record both for 1st half of january and 2nd half of december:
    !     - For a 14-records "climoz.nc": records 1 and 14.
    !     - For 12-records files:
    !       record 12 of "climoz_m.nc" if available, or record 1  of "climoz.nc".
    !       record 1  of "climoz_p.nc" if available, or record 12 of "climoz.nc".
    !   * Calendar is taken into account to get one record each day (not 360 always).
    !   * Missing values are filled in from sky to ground by copying lowest valid one.
    !     Attribute "missing_value" or "_FillValue" must be present in input file.
    !-------------------------------------------------------------------------------
    USE assert_m, ONLY: assert
    USE cal_tools_m, ONLY: year_len, mid_month
    !!  USE control_mod,        ONLY: anneeref
    USE time_phylmdz_mod, ONLY: annee_ref
    USE ioipsl, ONLY: ioget_year_len, ioget_calendar
    USE regr_conserv_m, ONLY: regr_conserv
    USE regr_lint_m, ONLY: regr_lint
    USE regular_lonlat_mod, ONLY: boundslon_reg, boundslat_reg, south, west, east
    USE slopes_m, ONLY: slopes
    USE lmdz_xios
    USE mod_phys_lmdz_para, ONLY: is_mpi_root, is_master, is_omp_master, gather, gather_mpi, bcast_mpi, klon_mpi
    USE geometry_mod, ONLY: latitude_deg, ind_cell_glo
    USE mod_grid_phy_lmdz, ONLY: klon_glo

    !-------------------------------------------------------------------------------
    ! Arguments:
    INTEGER, INTENT(IN) :: read_climoz ! read ozone climatology, 1 or 2
    !                         1: read a single ozone climatology used day and night
    !                         2: same + read also a daylight climatology
    LOGICAL, INTENT(IN) :: interpt     ! TRUE  => daily interpolation
    ! FALSE => no interpolation (14 months)
    !-------------------------------------------------------------------------------
    ! Local variables:

    !--- Input files variables
    INTEGER :: nlon_in                       ! Number of longitudes
    INTEGER :: nlat_in                       ! Number of latitudes
    INTEGER :: nlev_in                       ! Number of pressure levels
    INTEGER :: nmth_in                       ! Number of months
    REAL, ALLOCATABLE :: lon_in(:)           ! Longitudes   (ascending order, rad)
    REAL, ALLOCATABLE :: lat_in(:)           ! Latitudes    (ascending order, rad)
    REAL, ALLOCATABLE :: lev_in(:)           ! Pressure levels (ascen. order, hPa)
    REAL, ALLOCATABLE :: lon_in_edge(:)      ! Longitude intervals edges
    !              (ascending order,  / )
    REAL, ALLOCATABLE :: sinlat_in_edge(:)   ! Sinus of latitude intervals edges
    !              (ascending order,  / )
    LOGICAL :: ldec_lon, ldec_lat, ldec_lev  ! Decreasing order in input file
    CHARACTER(LEN = 20) :: cal_in              ! Calendar
    REAL, ALLOCATABLE :: o3_in3(:, :, :, :, :)   ! Ozone climatologies
    REAL, ALLOCATABLE :: o3_in3bis(:, :, :, :, :)   ! Ozone climatologies
    REAL, ALLOCATABLE :: o3_in2  (:, :, :, :)   ! Ozone climatologies
    REAL, ALLOCATABLE :: o3_in2bis(:, :, :, :, :)   ! Ozone climatologies
    ! last index: 1 for the day-night average, 2 for the daylight field.
    REAL :: NaN

    !--- Partially or totally regridded variables      (:,:,nlev_in,:,read_climoz)
    REAL, ALLOCATABLE :: o3_regr_lon   (:, :, :, :, :) ! (nlon_ou,nlat_in,:,0:13   ,:)
    REAL, ALLOCATABLE :: o3_regr_lonlat(:, :, :, :, :) ! (nlon_ou,nlat_ou,:,0:13   ,:)
    REAL, ALLOCATABLE :: o3_out3       (:, :, :, :, :) ! (nlon_ou,nlat_ou,:,ntim_ou,:)
    REAL, ALLOCATABLE :: o3_out3_glo   (:, :, :, :) !   (nbp_lat,:,ntim_ou,:)
    REAL, ALLOCATABLE :: o3_regr_lat     (:, :, :, :) !         (nlat_in,:,0:13   ,:)
    REAL, ALLOCATABLE :: o3_out2         (:, :, :, :) !         (nlat_ou,:,ntim_ou,:)
    REAL, ALLOCATABLE :: o3_out2_glo     (:, :, :, :) !         (nbp_lat,:,ntim_ou,:)
    REAL, ALLOCATABLE :: o3_out          (:, :, :, :) !         (nbp_lat,:,ntim_ou,:)
    ! Dimension number  | Interval                | Contains  | For variables:
    !   1 (longitude)   | [rlonu(i-1), rlonu(i)]  | rlonv(i)  | all
    !   2 (latitude)    | [rlatv(j), rlatv(j-1)]  | rlatu(j)  | all but o3_regr_lon
    !   3 (press level) |                         |   lev(k)  | all
    ! Note that rlatv(0)=pi/2 and rlatv(nlat_ou)=-pi/2.
    ! Dimension 4 is: month number                             (all vars but o3_out)
    !                 days elapsed since Jan. 1st 0h at mid-day (o3_out only)
    REAL, ALLOCATABLE :: v1(:)

    !--- For NetCDF:
    INTEGER :: fID_in_m, fID_in, levID_ou, dimid, vID_in(read_climoz), ntim_ou
    INTEGER :: fID_in_p, fID_ou, timID_ou, varid, vID_ou(read_climoz), ndims, ncerr
    INTEGER, ALLOCATABLE :: dIDs(:)
    CHARACTER(LEN = 20) :: cal_ou     !--- Calendar; no time inter => same as input
    CHARACTER(LEN = 80) :: press_unit !--- Pressure unit
    REAL :: tmidmonth(0:13)      !--- Elapsed days since Jan-1 0h at mid-months
    ! Additional records 0, 13 for interpolation
    REAL, ALLOCATABLE :: tmidday(:) !--- Output times (mid-days since Jan 1st 0h)
    LOGICAL :: lprev, lnext         !--- Flags: previous/next files are present
    LOGICAL :: l3D, l2D             !--- Flag:  input fields are 3D or zonal
    INTEGER :: ii, i, j, k, l, m, dln, ib, ie, iv, dx1, dx2
    INTEGER, ALLOCATABLE :: sta(:), cnt(:)
    CHARACTER(LEN = 80) :: sub, dim_nam, msg
    REAL :: null_array(0)
    LOGICAL, SAVE :: first = .TRUE.
    !$OMP THREADPRIVATE(first)
    REAL, ALLOCATABLE :: test_o3_in(:, :)
    REAL, ALLOCATABLE :: test_o3_out(:)

    IF (grid_type==unstructured) THEN
      IF (first) THEN
        IF (is_master) THEN
          ALLOCATE(latitude_glo(klon_glo))
          ALLOCATE(ind_cell_glo_glo(klon_glo))
        ELSE
          ALLOCATE(latitude_glo(0))
          ALLOCATE(ind_cell_glo_glo(0))
        ENDIF
        CALL gather(latitude_deg, latitude_glo)
        CALL gather(ind_cell_glo, ind_cell_glo_glo)
      ENDIF
    ENDIF

    IF (is_omp_master) THEN
      nlat_ou = nbp_lat
      nlon_ou = nbp_lon

      !-------------------------------------------------------------------------------
      IF (is_mpi_root) THEN
        sub = "regr_horiz_time_climoz"
        WRITE(lunout, *)"Call sequence information: " // TRIM(sub)
        CALL assert(read_climoz == 1 .OR. read_climoz == 2, "regr_lat_time_climoz")

        CALL  nf95_open("climoz.nc", nf90_nowrite, fID_in)
        lprev = nf90_open("climoz_m.nc", nf90_nowrite, fID_in_m)==nf90_noerr
        lnext = nf90_open("climoz_p.nc", nf90_nowrite, fID_in_p)==nf90_noerr

        !--- Get coordinates from the input file. Converts lon/lat in radians.
        !    Few inversions because "regr_conserv" and gcm need ascending vectors.
        CALL nf95_inq_varid(fID_in, vars_in(1), varid)
        CALL nf95_inquire_variable(fID_in, varid, dimids = dIDs, ndims = ndims)
        l3D = ndims==4; l2D = ndims==3
        IF(l3D) WRITE(lunout, *)"Input files contain full 3D ozone fields."
        IF(l2D) WRITE(lunout, *)"Input files contain zonal 2D ozone fields."
        DO i = 1, ndims
          CALL nf95_inquire_dimension(fID_in, dIDs(i), name = dim_nam, nclen = dln)
          CALL nf95_inq_varid(fID_in, dim_nam, varid)
          ii = i; IF(l2D) ii = i + 1                              !--- ndims==3:NO LONGITUDE
          SELECT CASE(ii)
          CASE(1)                                         !--- LONGITUDE
            CALL nf95_gw_var(fID_in, varid, lon_in)
            ldec_lon = lon_in(1)>lon_in(dln); IF(ldec_lon) lon_in = lon_in(dln:1:-1)
            nlon_in = dln; lon_in = lon_in * deg2rad
          CASE(2)                                         !--- LATITUDE
            CALL nf95_gw_var(fID_in, varid, lat_in)
            ldec_lat = lat_in(1)>lat_in(dln); IF(ldec_lat) lat_in = lat_in(dln:1:-1)
            nlat_in = dln; lat_in = lat_in * deg2rad
          CASE(3)                                         !--- PRESSURE LEVELS
            CALL nf95_gw_var(fID_in, varid, lev_in)
            ldec_lev = lev_in(1)>lev_in(dln); IF(ldec_lev) lev_in = lev_in(dln:1:-1)
            nlev_in = dln
            CALL nf95_get_att(fID_in, varid, "units", press_unit)
            k = LEN_TRIM(press_unit)
            DO WHILE(ICHAR(press_unit(k:k))==0)
              press_unit(k:k) = ' '; k = LEN_TRIM(press_unit) !--- REMOVE NULL END CHAR
            END DO
            IF(press_unit ==  "Pa") THEN
              lev_in = lev_in / 100.                        !--- CONVERT TO hPa
            ELSE IF(press_unit /= "hPa") THEN
              CALL abort_physic(sub, "the only recognized units are Pa and hPa.", 1)
            END IF
          CASE(4)                                         !--- TIME
            CALL nf95_inquire_dimension(fID_in, dIDs(i), nclen = nmth_in)
            cal_in = 'gregorian'
            IF(nf90_get_att(fID_in, varid, 'calendar', cal_in)/=nf90_noerr)        &
                    WRITE(lunout, *)'WARNING: missing "calendar" attribute for "' // &
                            TRIM(dim_nam) // '" in "climoz.nc". Choosing default: "gregorian".'
            k = LEN_TRIM(cal_in)
            DO WHILE(ICHAR(cal_in(k:k))==0)
              cal_in(k:k) = ' '; k = LEN_TRIM(cal_in)         !--- REMOVE NULL END CHAR
            END DO
          END SELECT
        END DO

        !--- Prepare quantities for time interpolation
        tmidmonth = mid_month(annee_ref, cal_in)
        IF(interpt) THEN
          ntim_ou = ioget_year_len(annee_ref)
          ALLOCATE(tmidday(ntim_ou))
          tmidday = [(REAL(k) - 0.5, k = 1, ntim_ou)]
          CALL ioget_calendar(cal_ou)
        ELSE
          ntim_ou = 14
          cal_ou = cal_in
        END IF
      ENDIF

      IF (grid_type==unstructured) THEN
        CALL bcast_mpi(nlon_in)
        CALL bcast_mpi(nlat_in)
        CALL bcast_mpi(nlev_in)
        CALL bcast_mpi(l3d)
        CALL bcast_mpi(tmidmonth)
        IF(interpt) CALL bcast_mpi(tmidday)
        CALL bcast_mpi(ntim_ou)

        IF (is_mpi_root) THEN
          CALL xios_set_domain_attr("domain_climoz", nj_glo = nlat_in, nj = nlat_in, jbegin = 0, latvalue_1d = lat_in / deg2rad)
          IF (l3D) THEN
            CALL xios_set_domain_attr("domain_climoz", ni_glo = nlon_in, ni = nlon_in, ibegin = 0, lonvalue_1d = lon_in / deg2rad)
          ELSE
            CALL xios_set_domain_attr("domain_climoz", ni_glo = 8, ni = 8, ibegin = 0, lonvalue_1d = (/ 0., 45., 90., 135., 180., 225., 270., 315. /))
          ENDIF
        ELSE
          CALL xios_set_domain_attr("domain_climoz", nj_glo = nlat_in, nj = 0, jbegin = 0, latvalue_1d = null_array)
          IF (l3D) THEN
            CALL xios_set_domain_attr("domain_climoz", ni_glo = nlon_in, ni = 0, ibegin = 0, lonvalue_1d = null_array)
          ELSE
            CALL xios_set_domain_attr("domain_climoz", ni_glo = 8, ni = 0, ibegin = 0, lonvalue_1d = null_array)
          ENDIF
        ENDIF
        CALL  xios_set_axis_attr("axis_climoz", n_glo = nlev_in)
        CALL  xios_set_axis_attr("time_axis_climoz", n_glo = ntim_ou)
        CALL  xios_set_axis_attr("time_axis_climoz", n_glo = ntim_ou)
        CALL  xios_set_axis_attr("tr_climoz", n_glo = read_climoz)
        CALL  xios_set_field_attr("tro3_out", enabled = .TRUE.)
        CALL  xios_set_field_attr("tro3_out", enabled = .TRUE.)

        IF (first) THEN
          first = .FALSE.
          RETURN
        ENDIF
      ENDIF

      IF (is_mpi_root) THEN
        !--- Longitudes management:
        !    * Need to shift data if the origin of input file longitudes /= -pi
        !    * Need to add some margin in longitude to ensure input interval contains
        !      all the output intervals => at least one longitudes slice has to be
        !      duplicated, possibly more for undersampling.
        IF(l3D) THEN
          IF (grid_type==unstructured) THEN
            dx2 = 0
          ELSE
            !--- Compute input edges longitudes vector (no end point yet)
            ALLOCATE(v1(nlon_in + 1))
            v1(1) = (lon_in(nlon_in) + lon_in(1)) / 2. - pi
            FORALL(i = 2:nlon_in) v1(i) = (lon_in(i - 1) + lon_in(i)) / 2.
            v1(nlon_in + 1) = v1(1) + 2. * pi
            DEALLOCATE(lon_in)

            !--- Shift input longitudes vector until it contains first output point boundslon_reg(1,west)
            v1 = v1 + 2 * pi * REAL(FLOOR((boundslon_reg(1, west) - v1(1)) / (2. * pi)))

            !--- Ensure first input longitudes interval contains first output point boundslon_reg(1,west)
            dx1 = locate(v1, boundslon_reg(1, west)) - 1
            v1 = CSHIFT(v1, SHIFT = dx1, DIM = 1)
            v1(nlon_in - dx1 + 2:) = v1(nlon_in - dx1 + 2:) + 2. * pi

            !--- Extend input longitudes vector until last interval contains boundslon_reg(nlat_ou,east)
            dx2 = 0; DO WHILE(v1(1 + dx2) + 2. * pi<=boundslon_reg(nlon_ou, east)); dx2 = dx2 + 1;
            END DO

            !--- Final edges longitudes vector (with margin and end point)
            ALLOCATE(lon_in_edge(nlon_in + dx2 + 1)); lon_in_edge = [v1, v1(2:1 + dx2) + 2. * pi]
            DEALLOCATE(v1)
          ENDIF
        END IF

        !--- Compute sinus of intervals edges latitudes:
        ALLOCATE(sinlat_in_edge(nlat_in + 1))
        sinlat_in_edge(1) = -1. ; sinlat_in_edge(nlat_in + 1) = 1.
        FORALL(j = 2:nlat_in) sinlat_in_edge(j) = SIN((lat_in(j - 1) + lat_in(j)) / 2.)
        DEALLOCATE(lat_in)



        !--- Check for contiguous years:
        ib = 0; ie = 13
        IF(nmth_in == 14) THEN; lprev = .FALSE.; lnext = .FALSE.
        WRITE(lunout, *)'Using 14 months ozone climatology "climoz.nc"...'
        ELSE
          IF(lprev) WRITE(lunout, *)'Using "climoz_m.nc" last record (previous year).'
          IF(.NOT.lprev) WRITE(lunout, *)"No previous year file ; assuming periodicity."
          IF(lnext) WRITE(lunout, *)'Using "climoz_p.nc" first record (next year).'
          IF(.NOT.lnext) WRITE(lunout, *)"No next year file ; assuming periodicity."
          IF(.NOT.lprev) ib = 1
          IF(.NOT.lnext) ie = 12
        END IF
        ALLOCATE(sta(ndims), cnt(ndims)); sta(:) = 1
        IF(l3D) cnt = [nlon_in, nlat_in, nlev_in, 1]
        IF(l2D) cnt = [        nlat_in, nlev_in, 1]
        IF(l3D) ALLOCATE(o3_in3(nlon_in + dx2, nlat_in, nlev_in, ib:ie, read_climoz))
        IF(l2D) ALLOCATE(o3_in2(nlat_in, nlev_in, ib:ie, read_climoz))

        !--- Read full current file and one record each available contiguous file
        DO iv = 1, read_climoz
          CALL nf95_inq_varid(fID_in, vars_in(1), vID_in(iv))
          IF(l3D) call nf95_get_var(fID_in, vID_in(iv), o3_in3(1:nlon_in, :, :, 1:12, iv))
          IF(l2D) call nf95_get_var(fID_in, vID_in(iv), o3_in2(:, :, 1:12, iv))
          IF(lprev) THEN; sta(ndims) = 12
          CALL nf95_inq_varid(fID_in_m, vars_in(1), vID_in(iv))
          IF(l3D) call nf95_get_var(fID_in_m, vID_in(iv), o3_in3(1:nlon_in, :, :, 0, iv), sta, cnt)
          IF(l2d) call nf95_get_var(fID_in_m, vID_in(iv), o3_in2(:, :, 0, iv), sta, cnt)
          END IF
          IF(lnext) THEN; sta(ndims) = 1
          CALL nf95_inq_varid(fID_in_p, vars_in(1), vID_in(iv))
          IF(l3D) call nf95_get_var(fID_in_p, vID_in(iv), o3_in3(1:nlon_in, :, :, 13, iv), sta, cnt)
          IF(l2D) call nf95_get_var(fID_in_p, vID_in(iv), o3_in2(:, :, 13, iv), sta, cnt)
          END IF
        END DO
        IF(lprev.OR.lnext) DEALLOCATE(sta, cnt)
        IF(lprev) CALL nf95_close(fID_in_m)
        IF(lnext) CALL nf95_close(fID_in_p)

        !--- Revert decreasing coordinates vector
        IF(l3D) THEN
          IF(ldec_lon) o3_in3(1:nlon_in, :, :, :, :) = o3_in3(nlon_in:1:-1, :, :, :, :)
          IF(ldec_lat) o3_in3 = o3_in3(:, nlat_in:1:-1, :, :, :)
          IF(ldec_lev) o3_in3 = o3_in3(:, :, nlev_in:1:-1, :, :)

          IF (grid_type /= unstructured) THEN
            !--- Shift values for longitude and duplicate some longitudes slices
            o3_in3(1:nlon_in, :, :, :, :) = CSHIFT(o3_in3(1:nlon_in, :, :, :, :), SHIFT = dx1, DIM = 1)
            o3_in3(nlon_in + 1:nlon_in + dx2, :, :, :, :) = o3_in3(1:dx2, :, :, :, :)
          ENDIF
        ELSE
          IF(ldec_lat) o3_in2 = o3_in2(nlat_in:1:-1, :, :, :)
          IF(ldec_lev) o3_in2 = o3_in2(:, nlev_in:1:-1, :, :)
        END IF

        !--- Deal with missing values
        DO m = 1, read_climoz
          WRITE(msg, '(a,i0)')"regr_lat_time_climoz: field Nr.", m
          IF(nf90_get_att(fID_in, vID_in(m), "missing_value", NaN)/= nf90_noerr) THEN
            IF(nf90_get_att(fID_in, vID_in(m), "_FillValue", NaN)/= nf90_noerr) THEN
              WRITE(lunout, *)TRIM(msg) // ": no missing value attribute found."; CYCLE
            END IF
          END IF
          WRITE(lunout, *)TRIM(msg) // ": missing value attribute found."
          WRITE(lunout, *)"Trying to fill in NaNs ; a full field would be better."

          !--- Check top layer contains no NaNs & search NaNs from top to ground
          msg = TRIM(sub) // ": NaNs in top layer !"
          IF(l3D) THEN
            IF(ANY(o3_in3(:, :, 1, :, m)==NaN)) CALL abort_physic(sub, msg, 1)
            DO k = 2, nlev_in
              WHERE(o3_in3(:, :, k, :, m)==NaN) o3_in3(:, :, k, :, m) = o3_in3(:, :, k - 1, :, m)
            END DO
          ELSE
            IF(ANY(o3_in2(:, 1, :, m)==NaN)) THEN
              WRITE(lunout, *)msg
              !--- Fill in latitudes where all values are missing
              DO l = 1, nmth_in
                !--- Next to south pole
                j = 1;       DO WHILE(o3_in2(j, 1, l, m)==NaN); j = j + 1;
                END DO
                IF(j>1) &
                        o3_in2(:j - 1, :, l, m) = SPREAD(o3_in2(j, :, l, m), DIM = 1, ncopies = j - 1)
                !--- Next to north pole
                j = nlat_in; DO WHILE(o3_in2(j, 1, l, m)==NaN); j = j + 1;
                END DO
                IF(j<nlat_in) &
                        o3_in2(j + 1:, :, l, m) = SPREAD(o3_in2(j, :, l, m), DIM = 1, ncopies = nlat_in - j)
              END DO
            END IF

            !--- Fill in high latitudes missing values
            !--- Highest level been filled-in, so has always valid values.
            DO k = 2, nlev_in
              WHERE(o3_in2(:, k, :, m)==NaN) o3_in2(:, k, :, m) = o3_in2(:, k - 1, :, m)
            END DO
          END IF
        END DO

      ENDIF

      !=============================================================================
      IF(l3D) THEN                                                   !=== 3D FIELDS
        !=============================================================================
        IF (grid_type==unstructured) THEN
          nlat_ou = klon_mpi

          IF (is_mpi_root) THEN
            ALLOCATE(o3_in3bis(nlon_in, nlat_in, nlev_in, 0:13, read_climoz))
            o3_in3bis(:, :, :, ib:ie, :) = o3_in3(1:nlon_in, :, :, ib:ie, :)
          ELSE
            ALLOCATE(o3_in3bis(0, 0, 0, 0, read_climoz))
          ENDIF
          ALLOCATE(o3_regr_lonlat(1, nlat_ou, nlev_in, 0:13, read_climoz))

          CALL xios_send_field("tro3_in", o3_in3bis(:, :, :, :, :))
          CALL xios_recv_field("tro3_out", o3_regr_lonlat(1, :, :, :, :))
        ELSE

          !--- Regrid in longitude
          ALLOCATE(o3_regr_lon(nlon_ou, nlat_in, nlev_in, ie - ib + 1, read_climoz))
          CALL regr_conserv(1, o3_in3, xs = lon_in_edge, &
                  xt = [boundslon_reg(1, west), boundslon_reg(:, east)], &
                  vt = o3_regr_lon, slope = slopes(1, o3_in3, lon_in_edge))
          DEALLOCATE(o3_in3)

          !--- Regrid in latitude: averaging with respect to SIN(lat) is
          !                        equivalent to weighting by COS(lat)
          !--- (inverted indices in "o3_regr_lonlat" because "rlatu" is decreasing)
          ALLOCATE(o3_regr_lonlat(nlon_ou, nlat_ou, nlev_in, 0:13, read_climoz))
          CALL regr_conserv(2, o3_regr_lon, xs = sinlat_in_edge, &
                  xt = [- 1., SIN(boundslat_reg(nlat_ou - 1:1:-1, south)), 1.], &
                  vt = o3_regr_lonlat(:, nlat_ou:1:- 1, :, ib:ie, :), &
                  slope = slopes(2, o3_regr_lon, sinlat_in_edge))
          DEALLOCATE(o3_regr_lon)

        ENDIF

        !--- Duplicate previous/next record(s) if they are not available
        IF(.NOT.lprev) o3_regr_lonlat(:, :, :, 0, :) = o3_regr_lonlat(:, :, :, 12, :)
        IF(.NOT.lnext) o3_regr_lonlat(:, :, :, 13, :) = o3_regr_lonlat(:, :, :, 1, :)

        !--- Regrid in time by linear interpolation:
        ALLOCATE(o3_out3(nlon_ou, nlat_ou, nlev_in, ntim_ou, read_climoz))
        IF(interpt) CALL regr_lint(4, o3_regr_lonlat, tmidmonth, tmidday, o3_out3)
        IF(.NOT.interpt) o3_out3 = o3_regr_lonlat
        DEALLOCATE(o3_regr_lonlat)

        nlat_ou = nbp_lat
        IF (grid_type==unstructured) THEN
          CALL xios_send_field('o3_out', o3_out3)
          ndims = 3
          ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
          CALL gather_mpi(o3_out3(1, :, :, :, :), o3_out3_glo)
        ENDIF

        !--- Create the output file and get the variable IDs:
        CALL prepare_out(fID_in, nlev_in, ntim_ou, fID_ou, levID_ou, timID_ou, vID_ou, &
                ndims, cal_ou)

        IF (is_mpi_root) THEN
          !--- Write remaining coordinate variables:
          CALL nf95_put_var(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
          IF(interpt) CALL nf95_put_var(fID_ou, timID_ou, tmidday)
          IF(.NOT.interpt) CALL nf95_put_var(fID_ou, timID_ou, tmidmonth)

          !--- Write to file (the order of "rlatu" is inverted in the output file):
          IF (grid_type==unstructured) THEN

            ALLOCATE(o3_out(nlat_ou, nlev_in, ntim_ou, read_climoz))
            DO i = 1, klon_glo
              o3_out(ind_cell_glo_glo(i), :, :, :) = o3_out3_glo(i, :, :, :)
            ENDDO

            DO m = 1, read_climoz
              CALL nf95_put_var(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1, :, :, m))
            END DO

          ELSE
            DO m = 1, read_climoz
              CALL nf95_put_var(fID_ou, vID_ou(m), o3_out3(:, nlat_ou:1:-1, :, :, m))
            END DO
          ENDIF
          CALL nf95_close(fID_ou)

        ENDIF


        !=============================================================================
      ELSE                                                         !=== ZONAL FIELDS
        !=============================================================================

        IF (grid_type==unstructured) THEN
          nlat_ou = klon_mpi

          IF (is_mpi_root) THEN
            ALLOCATE(o3_in2bis(8, nlat_in, nlev_in, 0:13, read_climoz))
            o3_in2bis(:, :, :, ib:ie, :) = SPREAD(o3_in2, 1, 8)
          ELSE
            ALLOCATE(o3_in2bis(0, 0, 0, 0, read_climoz))
          ENDIF
          ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
          CALL xios_send_field("tro3_in", o3_in2bis(:, :, :, :, :))
          CALL xios_recv_field("tro3_out", o3_regr_lat(:, :, :, :))
          IF(.NOT.lprev) o3_regr_lat(:, :, 0, :) = o3_regr_lat(:, :, 12, :)
          IF(.NOT.lnext) o3_regr_lat(:, :, 13, :) = o3_regr_lat(:, :, 1, :)

        ELSE
          !--- Regrid in latitude: averaging with respect to SIN(lat) is
          !                        equivalent to weighting by COS(lat)
          !--- (inverted indices in "o3_regr_lat" because "rlatu" is decreasing)
          ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
          CALL regr_conserv(1, o3_in2, xs = sinlat_in_edge, &
                  xt = [- 1., SIN(boundslat_reg(nlat_ou - 1:1:-1, south)), 1.], &
                  vt = o3_regr_lat(nlat_ou:1:- 1, :, ib:ie, :), &
                  slope = slopes(1, o3_in2, sinlat_in_edge))
          DEALLOCATE(o3_in2)

          !--- Duplicate previous/next record(s) if they are not available
          IF(.NOT.lprev) o3_regr_lat(:, :, 0, :) = o3_regr_lat(:, :, 12, :)
          IF(.NOT.lnext) o3_regr_lat(:, :, 13, :) = o3_regr_lat(:, :, 1, :)

        ENDIF

        !--- Regrid in time by linear interpolation:
        ALLOCATE(o3_out2(nlat_ou, nlev_in, ntim_ou, read_climoz))
        IF(interpt) CALL regr_lint(3, o3_regr_lat, tmidmonth, tmidday, o3_out2)
        IF(.NOT.interpt) o3_out2 = o3_regr_lat
        DEALLOCATE(o3_regr_lat)

        nlat_ou = nbp_lat

        IF (grid_type==unstructured) THEN
          ndims = 3
          ALLOCATE(o3_out2_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
          CALL gather_mpi(o3_out2, o3_out2_glo)
        ENDIF

        !--- Create the output file and get the variable IDs:
        CALL prepare_out(fID_in, nlev_in, ntim_ou, fID_ou, levID_ou, timID_ou, vID_ou, &
                ndims, cal_ou)

        IF (is_mpi_root) THEN

          !--- Write remaining coordinate variables:
          CALL nf95_put_var(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
          IF(interpt) CALL nf95_put_var(fID_ou, timID_ou, tmidday)
          IF(.NOT.interpt) CALL nf95_put_var(fID_ou, timID_ou, tmidmonth)

          IF (grid_type==unstructured) THEN

            ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
            DO i = 1, klon_glo
              o3_out(ind_cell_glo_glo(i), :, :, :) = o3_out2_glo(i, :, :, :)
            ENDDO

            DO m = 1, read_climoz
              CALL nf95_put_var(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1, :, :, m))
            END DO
          ELSE
            !--- Write to file (the order of "rlatu" is inverted in the output file):
            DO m = 1, read_climoz
              CALL nf95_put_var(fID_ou, vID_ou(m), o3_out2(nlat_ou:1:-1, :, :, m))
            END DO
          ENDIF

          CALL nf95_close(fID_ou)

        ENDIF

        !=============================================================================
      END IF
      !=============================================================================

      IF (is_mpi_root) CALL nf95_close(fID_in)

    ENDIF ! is_omp_master

    first = .FALSE.
  END SUBROUTINE regr_horiz_time_climoz

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


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

  SUBROUTINE prepare_out(fID_in, nlev_in, ntim_ou, fID_ou, vlevID, vtimID, &
          vID_ou, ndims, cal_ou)
    !-------------------------------------------------------------------------------
    ! Purpose:  This subroutine creates the NetCDF output file, defines
    !     dimensions and variables, and writes some of the coordinate variables.
    !-------------------------------------------------------------------------------
    USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
    USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
    USE mod_phys_lmdz_para, ONLY: is_mpi_root
    USE mod_grid_phy_lmdz, ONLY: klon_glo

    !-------------------------------------------------------------------------------
    ! Arguments:
    INTEGER, INTENT(IN) :: fID_in, nlev_in, ntim_ou
    INTEGER, INTENT(OUT) :: fID_ou, vlevID, vtimID
    INTEGER, INTENT(OUT) :: vID_ou(:)      ! dim(1/2) 1: O3day&night 2: O3daylight
    INTEGER, INTENT(IN) :: ndims          ! fields rank (3 or 4)
    CHARACTER(LEN = *), INTENT(IN) :: cal_ou ! calendar
    !-------------------------------------------------------------------------------
    ! Local variables:
    INTEGER :: dlonID, dlatID, dlevID, dtimID, dIDs(4)
    INTEGER :: vlonID, vlatID, ncerr, is
    REAL, ALLOCATABLE :: latitude_glo_(:)
    CHARACTER(LEN = 80) :: sub
    INTEGER :: i


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

    IF (is_mpi_root) THEN
      sub = "prepare_out"
      WRITE(lunout, *)"CALL sequence information: " // TRIM(sub)
      CALL nf95_create("climoz_LMDZ.nc", nf90_clobber, fID_ou)

      !--- Dimensions:
      IF(ndims==4) &
              CALL nf95_def_dim(fID_ou, "rlonv", nlon_ou, dlonID)
      CALL nf95_def_dim(fID_ou, "rlatu", nlat_ou, dlatID)
      CALL nf95_def_dim(fID_ou, "plev", nlev_in, dlevID)
      CALL nf95_def_dim(fID_ou, "time", ntim_ou, dtimID)

      !--- Define coordinate variables:
      IF(ndims==4) &
              CALL nf95_def_var(fID_ou, "rlonv", nf90_float, dlonID, vlonID)
      CALL nf95_def_var(fID_ou, "rlatu", nf90_float, dlatID, vlatID)
      CALL nf95_def_var(fID_ou, "plev", nf90_float, dlevID, vlevID)
      CALL nf95_def_var(fID_ou, "time", nf90_float, dtimID, vtimID)
      IF(ndims==4) &
              CALL nf95_put_att(fID_ou, vlonID, "units", "degrees_east")
      CALL nf95_put_att(fID_ou, vlatID, "units", "degrees_north")
      CALL nf95_put_att(fID_ou, vlevID, "units", "millibar")
      CALL nf95_put_att(fID_ou, vtimID, "units", "days since 2000-1-1")
      IF(ndims==4) &
              CALL nf95_put_att(fID_ou, vlonID, "standard_name", "longitude")
      CALL nf95_put_att(fID_ou, vlatID, "standard_name", "latitude")
      CALL nf95_put_att(fID_ou, vlevID, "standard_name", "air_pressure")
      CALL nf95_put_att(fID_ou, vtimID, "standard_name", "time")
      CALL nf95_put_att(fID_ou, vlevID, "long_name", "air pressure")
      CALL nf95_put_att(fID_ou, vtimID, "calendar", cal_ou)

      !--- Define the main variables:
      IF(ndims==3) dIDs(1:3) = [ dlatID, dlevID, dtimID]
      IF(ndims==4) dIDs = [dlonID, dlatID, dlevID, dtimID]
      CALL nf95_def_var(fID_ou, vars_in(1), nf90_float, dIDs(1:ndims), vID_ou(1))
      CALL nf95_put_att(fID_ou, vID_ou(1), "long_name", "ozone mole fraction")
      CALL nf95_put_att(fID_ou, vID_ou(1), "standard_name", "mole_fraction_of_ozone&
              _in_air")
      IF(SIZE(vID_ou) == 2) THEN
        CALL nf95_def_var(fID_ou, vars_in(2), nf90_float, dIDs(1:ndims), vID_ou(2))
        CALL nf95_put_att(fID_ou, vID_ou(2), "long_name", "ozone mole fraction in da&
                ylight")
      END IF

      !--- Global attributes:
      ! The following commands, copying attributes, may fail. That is OK.
      ! It should just mean that the attribute is not defined in the input file.
      CALL nf95_copy_att(fID_in, nf90_global, "Conventions", fID_ou, nf90_global, ncerr)
      CALL handle_err_copy_att("Conventions")
      CALL nf95_copy_att(fID_in, nf90_global, "title", fID_ou, nf90_global, ncerr)
      CALL handle_err_copy_att("title")
      CALL nf95_copy_att(fID_in, nf90_global, "institution", fID_ou, nf90_global, ncerr)
      CALL handle_err_copy_att("institution")
      CALL nf95_copy_att(fID_in, nf90_global, "source", fID_ou, nf90_global, ncerr)
      CALL handle_err_copy_att("source")
      CALL nf95_put_att (fID_ou, nf90_global, "comment", "Regridded for LMDZ")
      CALL nf95_enddef(fID_ou)

      IF (grid_type==unstructured) THEN
        ALLOCATE(latitude_glo_(klon_glo))
        DO i = 1, klon_glo
          latitude_glo_(ind_cell_glo_glo(i)) = latitude_glo(i)
        ENDDO
        CALL nf95_put_var(fID_ou, vlatID, latitude_glo_)
      ELSE
        !--- Write one of the coordinate variables:
        IF(ndims==4) CALL nf95_put_var(fID_ou, vlonID, lon_reg / deg2rad)
        CALL nf95_put_var(fID_ou, vlatID, lat_reg(nlat_ou:1:-1) / deg2rad)
        !    (convert from rad to degrees and sort in ascending order)
      ENDIF
    ENDIF

  CONTAINS

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

    SUBROUTINE handle_err_copy_att(att_name)

      !-------------------------------------------------------------------------------
      USE netcdf, ONLY: nf90_noerr, nf90_strerror
      !-------------------------------------------------------------------------------
      ! Arguments:
      CHARACTER(LEN = *), INTENT(IN) :: att_name
      !-------------------------------------------------------------------------------
      IF(ncerr /= nf90_noerr) &
              WRITE(lunout, *)TRIM(sub) // " prepare_out nf95_copy_att " // TRIM(att_name) // &
                      " -- " // TRIM(nf90_strerror(ncerr))

    END SUBROUTINE handle_err_copy_att

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

  END SUBROUTINE prepare_out

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

END MODULE regr_horiz_time_climoz_m

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