! $Id$
module regr_lat_time_climoz_m

  ! Author: Lionel GUEZ

  implicit none

  private
  public regr_lat_time_climoz

contains

  subroutine regr_lat_time_climoz

    ! "regr_lat_time_climoz" stands for "regrid latitude time
    ! climatology ozone".

    ! This procedure reads a climatology of ozone from a NetCDF file,
    ! regrids it in latitude and time, and writes the regridded field
    ! to a new NetCDF file.

    ! The input field depends on time, pressure level and latitude.

    ! If the input field has missing values, they must be signaled by
    ! the "missing_value" attribute. At each latitude and each date, the
    ! missing values are replaced by the lowest valid value above missing
    ! values.

    ! We assume that the input field is a step function of latitude
    ! and that the input latitude coordinate gives the centers of steps.
    ! Regridding in latitude is made by averaging, with a cosine of
    ! latitude factor.
    ! The target LMDZ latitude grid is the "scalar" grid: "rlatu".
    ! The values of "rlatu" are taken to be the centers of intervals.

    ! Regridding in time is by linear interpolation.
    ! Monthly values are processed to get daily values, on the basis
    ! of a 360-day calendar.
    ! The input file may contain either values for 12 months or values
    ! for 14 months.
    ! If there are 14 months then we assume that we have (in that order):
    ! December, January, February, ..., November, December, January
    ! We use the first December value to interpolate values between January
    ! 1st and mid-January.
    ! We use the last January value to interpolate values between
    ! mid-December and end of December.
    ! If there are only 12 months in the input file than we assume
    ! periodicity for interpolation at the beginning and at the end of the
    ! year.

    ! We assume that in the input file:
    ! -- latitude is in degrees;
    ! -- pressure is in hPa (even though we do not use pressure values
    ! here, we write the unit of pressure in the NetCDF header, and we
    ! will use the assumption later, when we regrid in pressure).
    ! -- latitude and pressure are strictly monotonic (as all NetCDF
    ! coordinate variables should be);
    ! -- time increases (even though we do not use values of the input
    ! time coordinate);
    ! -- missing values are vertically contiguous, at the bottom of
    ! the vertical domain;
    ! -- there is no missing value at the top level of the atmosphere.

    use regr1_step_av_m, only: regr1_step_av
    use regr3_lint_m, only: regr3_lint
    use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, handle_err, &
         nf95_put_var, nf95_gw_var
    use netcdf, only: nf90_nowrite, nf90_get_att, nf90_noerr

    ! Variables local to the procedure:

    include "dimensions.h"
    ! (for "jjm")
    include "paramet.h"
    ! (for the other included files)
    include "comgeom2.h"
    ! (for "rlatv")
    include "comconst.h"
    ! (for "pi")

    integer ncid_in, ncid_out ! NetCDF IDs for input and output files
    integer n_plev ! number of pressure levels in the input data
    integer n_lat! number of latitudes in the input data

    real, pointer:: latitude(:)
    ! (of input data, converted to rad, sorted in strictly ascending order)

    real, allocatable:: lat_in_edg(:)
    ! (edges of latitude intervals for input data, in rad, in strictly
    ! ascending order)

    real, pointer:: plev(:)
    ! pressure levels of input data, sorted in strictly ascending order

    logical desc_lat ! latitude in descending order in the input file
    logical desc_plev ! pressure levels in descending order in the input file

    real, pointer:: o3_in(:, :, :) ! (n_lat, n_plev, 12 or 0:13)
    ! (ozone climatology from the input file)
    ! ("o3_in(j, k, l)" is at latitude "latitude(j)" and pressure
    ! level "plev(k)". "l" is between 1 and 12 or between 0 and 13)

    real missing_value

    real, allocatable:: o3_regr_lat(:, :, :) ! (jjm + 1, n_plev, 0:13)
    ! (mean of "o3_in" over a latitude interval of LMDZ)
    ! (First dimension is latitude interval.
    ! The latitude interval for "o3_regr_lat(j,:, :)" contains "rlatu(j)".
    ! If "j" is between 2 and "jjm" then the interval is:
    ! [rlatv(j), rlatv(j-1)]
    ! If "j" is 1 or "jjm + 1" then the interval is:
    ! [rlatv(1), pi / 2]
    ! or:
    ! [- pi / 2, rlatv(jjm)]
    ! respectively.
    ! "o3_regr_lat(:, k, :)" is for pressure level "plev(k)".
    ! Last dimension is month number.)

    real, allocatable:: o3_out(:, :, :) ! (jjm + 1, n_plev, 360)
    ! (regridded ozone climatology)
    ! ("o3_out(j, k, l)" is at latitude "rlatu(j)", pressure
    ! level "plev(k)" and date "January 1st 0h" + "tmidday(l)", in a
    ! 360-day calendar.)

    integer j, k, l

    ! For NetCDF:
    integer varid_in, varid_out, varid_plev, varid_time, varid, ncerr

    real, parameter:: tmidmonth(0:13) = (/(-15. + 30. * l, l = 0, 13)/)
    ! (time to middle of month, in days since January 1st 0h, in a
    ! 360-day calendar)
    ! (We add values -15 and 375 so that, for example, day 3 of the year is
    ! interpolated between the December and the January value.)

    real, parameter:: tmidday(360) = (/(l + 0.5, l = 0, 359)/)
    ! (time to middle of day, in days since January 1st 0h, in a
    ! 360-day calendar)

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

    print *, "Call sequence information: regr_lat_time_climoz"

    call nf95_open("climoz.nc", nf90_nowrite, ncid_in)

    ! Get coordinates from the input file:

    call nf95_inq_varid(ncid_in, "latitude", varid)
    call nf95_gw_var(ncid_in, varid, latitude)
    ! Convert from degrees to rad, because we will take the sine of latitude:
    latitude = latitude / 180. * pi
    n_lat = size(latitude)
    ! We need to supply the latitudes to "regr1_step_av" in
    ! ascending order, so invert order if necessary:
    desc_lat = latitude(1) > latitude(n_lat)
    if (desc_lat) latitude = latitude(n_lat:1:-1)

    ! Compute edges of latitude intervals:
    allocate(lat_in_edg(n_lat + 1))
    lat_in_edg(1) = - pi / 2
    forall (j = 2:n_lat) lat_in_edg(j) = (latitude(j - 1) + latitude(j)) / 2
    lat_in_edg(n_lat + 1) = pi / 2
    deallocate(latitude) ! pointer

    call nf95_inq_varid(ncid_in, "plev", varid)
    call nf95_gw_var(ncid_in, varid, plev)
    n_plev = size(plev)
    ! We only need the pressure coordinate to copy it to the output file.
    ! The program "gcm" will assume that pressure levels are in
    ! ascending order in the regridded climatology so invert order if
    ! necessary:
    desc_plev = plev(1) > plev(n_plev)
    if (desc_plev) plev = plev(n_plev:1:-1)

    ! Create the output file and get the variable IDs:
    call prepare_out(ncid_in, n_plev, ncid_out, varid_out, varid_plev, &
         varid_time)

    ! Write remaining coordinate variables:
    call nf95_put_var(ncid_out, varid_plev, plev)
    call nf95_put_var(ncid_out, varid_time, tmidday)

    deallocate(plev) ! pointer

    call nf95_inq_varid(ncid_in, "tro3", varid_in)
    call nf95_gw_var(ncid_in, varid_in, o3_in)
    if (desc_lat) o3_in = o3_in(n_lat:1:-1, :, :)
    if (desc_plev) o3_in = o3_in(:, n_plev:1:-1, :)

    ncerr = nf90_get_att(ncid_in, varid_in, "missing_value", missing_value)
    if (ncerr == nf90_noerr) then
       do l = 1, size(o3_in, 3)
          do j = 1, n_lat
             ! Find missing values, starting from top of atmosphere
             ! and going down.
             ! We assume that the highest level contains no missing value.
             k = 2
             do
                if (o3_in(j, k, l) == missing_value .or. k == n_plev) exit
                k = k + 1
             end do
             ! Replace missing values with the valid value at the
             ! lowest level above missing values:
             if (o3_in(j, k, l) == missing_value) &
                  o3_in(j, k:n_plev, l) = o3_in(j, k-1, l)
          end do
       end do
    else
       print *, "regr_lat_time_climoz: no missing value attribute"
    end if
    
    call nf95_close(ncid_in)

    allocate(o3_regr_lat(jjm + 1, n_plev, 0:13))
    allocate(o3_out(jjm + 1, n_plev, 360))

    ! Regrid in latitude:
    ! We average with respect to sine of latitude, which is
    ! equivalent to weighting by cosine of latitude:
    if (size(o3_in, 3) == 12) then
       print *, "Found 12 months in ozone climatology, assuming periodicity..."
       o3_regr_lat(jjm+1:1:-1, :, 1:12) = regr1_step_av(o3_in, &
            xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(jjm:1:-1), pi / 2/)))
       ! (invert order of indices in "o3_regr_lat" because "rlatu" is
       ! in descending order)

       ! Duplicate January and December values, in preparation of time
       ! interpolation:
       o3_regr_lat(:, :, 0) = o3_regr_lat(:, :, 12)
       o3_regr_lat(:, :, 13) = o3_regr_lat(:, :, 1)
    else
       print *, "Using 14 months in ozone climatology..."
       o3_regr_lat(jjm+1:1:-1, :, :) = regr1_step_av(o3_in, &
            xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(jjm:1:-1), pi / 2/)))
       ! (invert order of indices in "o3_regr_lat" because "rlatu" is
       ! in descending order)
    end if
    
    deallocate(o3_in) ! pointer

    ! Regrid in time by linear interpolation:
    o3_out = regr3_lint(o3_regr_lat, tmidmonth, tmidday)

    ! Write to file:
    call nf95_put_var(ncid_out, varid_out, o3_out(jjm+1:1:-1, :, :))
    ! (The order of "rlatu" is inverted in the output file)

    call nf95_close(ncid_out)

  end subroutine regr_lat_time_climoz

  !********************************************

  subroutine prepare_out(ncid_in, n_plev, ncid_out, varid_out, varid_plev, &
       varid_time)

    ! This subroutine creates the NetCDF output file, defines
    ! dimensions and variables, and writes one of the coordinate variables.

    use netcdf95, only: nf95_create, nf95_def_dim, nf95_def_var, &
         nf95_put_att, nf95_enddef, nf95_copy_att, nf95_put_var
    use netcdf, only: nf90_clobber, nf90_float, nf90_global

    integer, intent(in):: ncid_in, n_plev
    integer, intent(out):: ncid_out, varid_out, varid_plev, varid_time

    ! Variables local to the procedure:

    include "dimensions.h"
    ! (for "jjm")
    include "paramet.h"
    ! (for the other included files)
    include "comgeom2.h"
    ! (for "rlatu")
    include "comconst.h"
    ! (for "pi")

    integer ncerr
    integer dimid_rlatu, dimid_plev, dimid_time
    integer varid_rlatu

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

    print *, "Call sequence information: prepare_out"

    call nf95_create("climoz_LMDZ.nc", nf90_clobber, ncid_out)

    ! Dimensions:
    call nf95_def_dim(ncid_out, "time", 360, dimid_time)
    call nf95_def_dim(ncid_out, "plev", n_plev, dimid_plev)
    call nf95_def_dim(ncid_out, "rlatu", jjm + 1, dimid_rlatu)

    ! Define coordinate variables:

    call nf95_def_var(ncid_out, "time", nf90_float, dimid_time, varid_time)
    call nf95_put_att(ncid_out, varid_time, "units", "days since 2000-1-1")
    call nf95_put_att(ncid_out, varid_time, "calendar", "360_day")
    call nf95_put_att(ncid_out, varid_time, "standard_name", "time")

    call nf95_def_var(ncid_out, "plev", nf90_float, dimid_plev, varid_plev)
    call nf95_put_att(ncid_out, varid_plev, "units", "millibar")
    call nf95_put_att(ncid_out, varid_plev, "standard_name", "air_pressure")
    call nf95_put_att(ncid_out, varid_plev, "long_name", "air pressure")

    call nf95_def_var(ncid_out, "rlatu", nf90_float, dimid_rlatu, varid_rlatu)
    call nf95_put_att(ncid_out, varid_rlatu, "units", "degrees_north")
    call nf95_put_att(ncid_out, varid_rlatu, "standard_name", "latitude")

    ! Define the primary variable:

    call nf95_def_var(ncid_out, "tro3", nf90_float, &
         (/dimid_rlatu, dimid_plev, dimid_time/), varid_out)
    call nf95_put_att(ncid_out, varid_out, "long_name", "ozone mole fraction")
    call nf95_put_att(ncid_out, varid_out, "standard_name", &
         "mole_fraction_of_ozone_in_air")

    ! 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(ncid_in, nf90_global, "Conventions", ncid_out, &
         nf90_global, ncerr)
    call handle_err_copy_att("Conventions")

    call nf95_copy_att(ncid_in, nf90_global, "title", ncid_out, nf90_global, &
         ncerr)
    call handle_err_copy_att("title")

    call nf95_copy_att(ncid_in, nf90_global, "institution", ncid_out, &
         nf90_global, ncerr)
    call handle_err_copy_att("institution")

    call nf95_copy_att(ncid_in, nf90_global, "source", ncid_out, nf90_global, &
         ncerr)
    call handle_err_copy_att("source")

    call nf95_put_att(ncid_out, nf90_global, "comment", "Regridded for LMDZ")

    call nf95_enddef(ncid_out)

    ! Write one of the coordinate variables:
    call nf95_put_var(ncid_out, varid_rlatu, rlatu(jjm+1:1:-1) / pi * 180.)
    ! (convert from rad to degrees and sort in ascending order)

  contains

    subroutine handle_err_copy_att(att_name)

      use netcdf, only: nf90_noerr, nf90_strerror

      character(len=*), intent(in):: att_name

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

      if (ncerr /= nf90_noerr) then
         print *, "regr_lat_time_climoz_m prepare_out nf95_copy_att " &
              // att_name // " -- " // trim(nf90_strerror(ncerr))
      end if

    end subroutine handle_err_copy_att

  end subroutine prepare_out

end module regr_lat_time_climoz_m