! $Id$

SUBROUTINE readaerosol_interp(id_aero, itap, pdtphys, r_day, first, pplay, paprs, t_seri, mass_out, pi_mass_out, load_src)

  ! This routine will return the mass concentration at actual day(mass_out) and
  ! the pre-industrial values(pi_mass_out) for aerosol corresponding to "id_aero".
  ! The mass concentrations for all aerosols are saved in this routine but each
  ! CALL to this routine only treats the aerosol "id_aero".

  ! 1) Read in data for the whole year, ONLY at first time step
  ! 2) Interpolate to the actual day, ONLY at new day
  ! 3) Interpolate to the model vertical grid (target grid), ONLY at new day
  ! 4) Test for negative mass values

  USE ioipsl
  USE dimphy, ONLY: klev, klon
  USE lmdz_phys_para, ONLY: mpi_rank
  USE readaerosol_mod
  USE aero_mod, ONLY: naero_spc, name_aero
  USE lmdz_write_field_phy
  USE phys_cal_mod
  USE lmdz_pres2lev
  USE lmdz_print_control, ONLY: lunout
  USE lmdz_abort_physic, ONLY: abort_physic

  IMPLICIT NONE

  INCLUDE "YOMCST.h"
  INCLUDE "chem.h"
  INCLUDE "clesphys.h"

  ! Input:
  !****************************************************************************************
  INTEGER, INTENT(IN) :: id_aero! Identity number for the aerosol to treat
  INTEGER, INTENT(IN) :: itap   ! Physic step count
  REAL, INTENT(IN) :: pdtphys! Physic day step
  REAL, INTENT(IN) :: r_day  ! Day of integration
  LOGICAL, INTENT(IN) :: first  ! First model timestep
  REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay  ! pression at model mid-layers
  REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs  ! pression between model layers
  REAL, DIMENSION(klon, klev), INTENT(IN) :: t_seri ! air temperature

  ! Output:
  !****************************************************************************************
  REAL, INTENT(OUT) :: mass_out(klon, klev)    ! Mass of aerosol (monthly mean data,from file) [ug AIBCM/m3]
  REAL, INTENT(OUT) :: pi_mass_out(klon, klev) ! Mass of preindustrial aerosol (monthly mean data,from file) [ug AIBCM/m3]
  REAL, INTENT(OUT) :: load_src(klon) ! Load of aerosol (monthly mean data,from file) [kg/m3]

  ! Local Variables:
  !****************************************************************************************
  INTEGER :: i, k, ierr
  INTEGER :: iday, iyr, lmt_pas
  !  INTEGER                         :: im, day1, day2, im2
  INTEGER :: im, im2
  REAL :: day1, day2
  INTEGER :: pi_klev_src ! Only for testing purpose
  INTEGER, SAVE :: klev_src    ! Number of vertical levles in source field
  !$OMP THREADPRIVATE(klev_src)

  REAL :: zrho      ! Air density [kg/m3]
  REAL :: volm      ! Volyme de melange [kg/kg]
  REAL, DIMENSION(klon) :: psurf_day, pi_psurf_day
  REAL, DIMENSION(klon) :: pi_load_src  ! Mass load at source grid
  REAL, DIMENSION(klon) :: load_tgt, load_tgt_test
  REAL, DIMENSION(klon, klev) :: delp ! pressure difference in each model layer

  REAL, ALLOCATABLE, DIMENSION(:, :) :: pplay_src ! pression mid-layer at source levels
  REAL, ALLOCATABLE, DIMENSION(:, :) :: tmp1, tmp2  ! Temporary variables
  REAL, ALLOCATABLE, DIMENSION(:, :, :, :), SAVE :: var_year    ! VAR in right dimension for the total year
  REAL, ALLOCATABLE, DIMENSION(:, :, :, :), SAVE :: pi_var_year ! pre-industrial VAR, -"-
  !$OMP THREADPRIVATE(var_year,pi_var_year)
  REAL, ALLOCATABLE, DIMENSION(:, :, :), SAVE :: var_day     ! VAR interpolated to the actual day and model grid
  REAL, ALLOCATABLE, DIMENSION(:, :, :), SAVE :: pi_var_day  ! pre-industrial VAR, -"-
  !$OMP THREADPRIVATE(var_day,pi_var_day)
  REAL, ALLOCATABLE, DIMENSION(:, :, :), SAVE :: psurf_year, pi_psurf_year ! surface pressure for the total year
  !$OMP THREADPRIVATE(psurf_year, pi_psurf_year)
  REAL, ALLOCATABLE, DIMENSION(:, :, :), SAVE :: load_year, pi_load_year   ! load in the column for the total year
  !$OMP THREADPRIVATE(load_year, pi_load_year)

  REAL, DIMENSION(:, :, :), POINTER :: pt_tmp      ! Pointer allocated in readaerosol
  REAL, POINTER, DIMENSION(:), SAVE :: pt_ap, pt_b ! Pointer for describing the vertical levels 
  !$OMP THREADPRIVATE(pt_ap, pt_b)
  INTEGER, SAVE :: nbr_tsteps ! number of time steps in file read
  REAL, DIMENSION(14), SAVE :: month_len, month_start, month_mid
  !$OMP THREADPRIVATE(nbr_tsteps, month_len, month_start, month_mid)
  REAL :: jDay

  LOGICAL :: lnewday      ! Indicates if first time step at a new day
  LOGICAL :: OLDNEWDAY
  LOGICAL, SAVE :: vert_interp  ! Indicates if vertical interpolation will be done
  LOGICAL, SAVE :: debug = .FALSE.! Debugging in this subroutine
  !$OMP THREADPRIVATE(vert_interp, debug)
  CHARACTER(len = 8) :: type
  CHARACTER(len = 8) :: filename


  !****************************************************************************************
  ! Initialization

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

  ! Calculation to find if it is a new day

  IF(mpi_rank == 0 .AND. debug)THEN
    PRINT*, 'CONTROL PANEL REGARDING TIME STEPING'
  ENDIF

  ! Use phys_cal_mod
  iday = day_cur
  iyr = year_cur
  im = mth_cur

  !  iday = INT(r_day)
  !  iyr  = iday/360
  !  iday = iday-iyr*360         ! day of the actual year
  !  iyr  = iyr + annee_ref      ! year of the run
  !  im   = iday/30 +1           ! the actual month
  CALL ymds2ju(iyr, im, iday, 0., jDay)
  !   CALL ymds2ju(iyr, im, iday-(im-1)*30, 0., jDay)

  IF(MOD(itap - 1, NINT(86400. / pdtphys)) == 0)THEN
    lnewday = .TRUE.
  ELSE
    lnewday = .FALSE.
  ENDIF

  IF(mpi_rank == 0 .AND. debug)THEN
    ! 0.02 is about 0.5/24, namly less than half an hour
    OLDNEWDAY = (r_day - REAL(iday) < 0.02)
    ! Once per day, update aerosol fields
    lmt_pas = NINT(86400. / pdtphys)
    PRINT*, 'r_day-REAL(iday) =', r_day - REAL(iday)
    PRINT*, 'itap =', itap
    PRINT*, 'pdtphys =', pdtphys
    PRINT*, 'lmt_pas =', lmt_pas
    PRINT*, 'iday =', iday
    PRINT*, 'r_day =', r_day
    PRINT*, 'day_cur =', day_cur
    PRINT*, 'mth_cur =', mth_cur
    PRINT*, 'year_cur =', year_cur
    PRINT*, 'NINT(86400./pdtphys) =', NINT(86400. / pdtphys)
    PRINT*, 'MOD(0,1) =', MOD(0, 1)
    PRINT*, 'lnewday =', lnewday
    PRINT*, 'OLDNEWDAY =', OLDNEWDAY
  ENDIF

  IF (.NOT. ALLOCATED(var_day)) THEN
    ALLOCATE(var_day(klon, klev, naero_spc), stat = ierr)
    IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 1', 1)
    ALLOCATE(pi_var_day(klon, klev, naero_spc), stat = ierr)
    IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 2', 1)

    ALLOCATE(psurf_year(klon, 12, naero_spc), pi_psurf_year(klon, 12, naero_spc), stat = ierr)
    IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 3', 1)

    ALLOCATE(load_year(klon, 12, naero_spc), pi_load_year(klon, 12, naero_spc), stat = ierr)
    IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 4', 1)

    lnewday = .TRUE.

    NULLIFY(pt_ap)
    NULLIFY(pt_b)
  ENDIF

  !****************************************************************************************
  ! 1) Read in data : corresponding to the actual year and preindustrial data.
  !    Only for the first day of the year.

  !****************************************************************************************
  IF ((first .OR. iday==0) .AND. lnewday) THEN
    NULLIFY(pt_tmp)

    ! Reading values corresponding to the closest year taking into count the choice of aer_type.
    ! For aer_type=scenario interpolation between 2 data sets is done in readaerosol.
    ! If aer_type=mix1, mix2 or mix3, the run type and file name depends on the aerosol.
    IF (aer_type=='preind' .OR. aer_type=='actuel' .OR. aer_type=='annuel' .OR. aer_type=='scenario') THEN
      ! Standard case
      filename = 'aerosols'
      type = aer_type
    ELSE IF (aer_type == 'mix1') THEN
      ! Special case using a mix of decenal sulfate file and annual aerosols(all aerosols except sulfate)
      IF (name_aero(id_aero) == 'SO4') THEN
        filename = 'so4.run '
        type = 'scenario'
      ELSE
        filename = 'aerosols'
        type = 'annuel'
      ENDIF
    ELSE IF (aer_type == 'mix2') THEN
      ! Special case using a mix of decenal sulfate file and natrual aerosols
      IF (name_aero(id_aero) == 'SO4') THEN
        filename = 'so4.run '
        type = 'scenario'
      ELSE
        filename = 'aerosols'
        type = 'preind'
      ENDIF
    ELSE IF (aer_type == 'mix3') THEN
      ! Special case using a mix of annual sulfate file and natrual aerosols
      IF (name_aero(id_aero) == 'SO4') THEN
        filename = 'aerosols'
        type = 'annuel'
      ELSE
        filename = 'aerosols'
        type = 'preind'
      ENDIF
    ELSE
      CALL abort_physic('readaerosol_interp', 'this aer_type not supported', 1)
    ENDIF

    CALL readaerosol(name_aero(id_aero), type, filename, iyr, klev_src, pt_ap, pt_b, pt_tmp, &
            psurf_year(:, :, id_aero), load_year(:, :, id_aero))
    IF (.NOT. ALLOCATED(var_year)) THEN
      ALLOCATE(var_year(klon, klev_src, 12, naero_spc), stat = ierr)
      IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 5', 1)
    ENDIF
    var_year(:, :, :, id_aero) = pt_tmp(:, :, :)

    ! Reading values corresponding to the preindustrial concentrations.
    type = 'preind'
    CALL readaerosol(name_aero(id_aero), type, filename, iyr, pi_klev_src, pt_ap, pt_b, pt_tmp, &
            pi_psurf_year(:, :, id_aero), pi_load_year(:, :, id_aero))

    ! klev_src must be the same in both files.
    ! Also supposing pt_ap and pt_b to be the same in the 2 files without testing.
    IF (pi_klev_src /= klev_src) THEN
      WRITE(lunout, *) 'Error! All forcing files for the same aerosol must have the same vertical dimension'
      WRITE(lunout, *) 'Aerosol : ', name_aero(id_aero)
      CALL abort_physic('readaerosol_interp', 'Differnt vertical axes in aerosol forcing files', 1)
    ENDIF

    IF (.NOT. ALLOCATED(pi_var_year)) THEN
      ALLOCATE(pi_var_year(klon, klev_src, 12, naero_spc), stat = ierr)
      IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 6', 1)
    ENDIF
    pi_var_year(:, :, :, id_aero) = pt_tmp(:, :, :)

    IF (debug) THEN
      CALL writefield_phy('var_year_jan', var_year(:, :, 1, id_aero), klev_src)
      CALL writefield_phy('var_year_dec', var_year(:, :, 12, id_aero), klev_src)
      CALL writefield_phy('psurf_src', psurf_year(:, :, id_aero), 1)
      CALL writefield_phy('pi_psurf_src', pi_psurf_year(:, :, id_aero), 1)
      CALL writefield_phy('load_year_src', load_year(:, :, id_aero), 1)
      CALL writefield_phy('pi_load_year_src', pi_load_year(:, :, id_aero), 1)
    ENDIF

    ! Pointer no more useful, deallocate.
    DEALLOCATE(pt_tmp)

    ! Test if vertical interpolation will be needed.
    IF (psurf_year(1, 1, id_aero)==not_valid .OR. pi_psurf_year(1, 1, id_aero)==not_valid) THEN
      ! Pressure=not_valid indicates old file format, see module readaerosol
      vert_interp = .FALSE.

      ! If old file format, both psurf_year and pi_psurf_year must be not_valid
      IF (psurf_year(1, 1, id_aero) /= pi_psurf_year(1, 1, id_aero)) THEN
        WRITE(lunout, *) 'Warning! All forcing files for the same aerosol must have the same structure'
        CALL abort_physic('readaerosol_interp', 'The aerosol files have not the same format', 1)
      ENDIF

      IF (klev /= klev_src) THEN
        WRITE(lunout, *) 'Old format of aerosol file do not allowed vertical interpolation'
        CALL abort_physic('readaerosol_interp', 'Old aerosol file not possible', 1)
      ENDIF

    ELSE
      vert_interp = .TRUE.
    ENDIF

    !    Calendar initialisation

    DO i = 2, 13
      month_len(i) = REAL(ioget_mon_len(year_cur, i - 1))
      CALL ymds2ju(year_cur, i - 1, 1, 0.0, month_start(i))
    ENDDO
    month_len(1) = REAL(ioget_mon_len(year_cur - 1, 12))
    CALL ymds2ju(year_cur - 1, 12, 1, 0.0, month_start(1))
    month_len(14) = REAL(ioget_mon_len(year_cur + 1, 1))
    CALL ymds2ju(year_cur + 1, 1, 1, 0.0, month_start(14))
    month_mid(:) = month_start (:) + month_len(:) / 2.

    IF (debug) THEN
      WRITE(lunout, *)' month_len = ', month_len
      WRITE(lunout, *)' month_mid = ', month_mid
    endif

  ENDIF  ! IF ( (first .OR. iday==0) .AND. lnewday ) THEN 

  !****************************************************************************************
  ! - 2) Interpolate to the actual day.
  ! - 3) Interpolate to the model vertical grid.

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

  IF (lnewday) THEN ! only if new day
    !****************************************************************************************
    ! 2) Interpolate to the actual day

    !****************************************************************************************
    ! Find which months and days to use for time interpolation
    nbr_tsteps = 12
    IF (nbr_tsteps == 12) THEN
      IF (jDay < month_mid(im + 1)) THEN
        im2 = im - 1
        day2 = month_mid(im2 + 1)
        day1 = month_mid(im + 1)
        IF (im2 <= 0) THEN
          ! the month is january, thus the month before december
          im2 = 12
        ENDIF
      ELSE
        ! the second half of the month
        im2 = im + 1
        day1 = month_mid(im + 1)
        day2 = month_mid(im2 + 1)
        IF (im2 > 12) THEN
          ! the month is december, the following thus january
          im2 = 1
        ENDIF
      ENDIF
    ELSE IF (nbr_tsteps == 14) THEN
      im = im + 1
      IF (jDay < month_mid(im)) THEN
        ! in the first half of the month use month before and actual month
        im2 = im - 1
        day2 = month_mid(im2)
        day1 = month_mid(im)
      ELSE
        ! the second half of the month
        im2 = im + 1
        day1 = month_mid(im)
        day2 = month_mid(im2)
      ENDIF
    ELSE
      CALL abort_physic('readaerosol_interp', 'number of months undefined', 1)
    ENDIF
    IF (debug) THEN
      WRITE(lunout, *)' jDay, day1, day2, im, im2 = ', jDay, day1, day2, im, im2
    endif


    ! Time interpolation, still on vertical source grid
    ALLOCATE(tmp1(klon, klev_src), tmp2(klon, klev_src), stat = ierr)
    IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 7', 1)

    ALLOCATE(pplay_src(klon, klev_src), stat = ierr)
    IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 8', 1)

    DO k = 1, klev_src
      DO i = 1, klon
        tmp1(i, k) = &
                var_year(i, k, im2, id_aero) - (jDay - day2) / (day1 - day2) * &
                        (var_year(i, k, im2, id_aero) - var_year(i, k, im, id_aero))

        tmp2(i, k) = &
                pi_var_year(i, k, im2, id_aero) - (jDay - day2) / (day1 - day2) * &
                        (pi_var_year(i, k, im2, id_aero) - pi_var_year(i, k, im, id_aero))
      ENDDO
    ENDDO

    ! Time interpolation for pressure at surface, still on vertical source grid
    DO i = 1, klon
      psurf_day(i) = &
              psurf_year(i, im2, id_aero) - (jDay - day2) / (day1 - day2) * &
                      (psurf_year(i, im2, id_aero) - psurf_year(i, im, id_aero))

      pi_psurf_day(i) = &
              pi_psurf_year(i, im2, id_aero) - (jDay - day2) / (day1 - day2) * &
                      (pi_psurf_year(i, im2, id_aero) - pi_psurf_year(i, im, id_aero))
    ENDDO

    ! Time interpolation for the load, still on vertical source grid
    DO i = 1, klon
      load_src(i) = &
              load_year(i, im2, id_aero) - (jDay - day2) / (day1 - day2) * &
                      (load_year(i, im2, id_aero) - load_year(i, im, id_aero))

      pi_load_src(i) = &
              pi_load_year(i, im2, id_aero) - (jDay - day2) / (day1 - day2) * &
                      (pi_load_year(i, im2, id_aero) - pi_load_year(i, im, id_aero))
    ENDDO

    !****************************************************************************************
    ! 3) Interpolate to the model vertical grid (target grid)

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

    IF (vert_interp) THEN

      ! - Interpolate variable tmp1 (on source grid) to var_day (on target grid)
      !********************************************************************************
      ! a) calculate pression at vertical levels for the source grid using the
      !    hybrid-sigma coordinates ap and b and the surface pressure, variables from file.
      DO k = 1, klev_src
        DO i = 1, klon
          pplay_src(i, k) = pt_ap(k) + pt_b(k) * psurf_day(i)
        ENDDO
      ENDDO

      IF (debug) THEN
        CALL writefield_phy('psurf_day_src', psurf_day(:), 1)
        CALL writefield_phy('pplay_src', pplay_src(:, :), klev_src)
        CALL writefield_phy('pplay', pplay(:, :), klev)
        CALL writefield_phy('day_src', tmp1, klev_src)
        CALL writefield_phy('pi_day_src', tmp2, klev_src)
      ENDIF

      ! b) vertical interpolation on pressure leveles
      CALL pres2lev(tmp1(:, :), var_day(:, :, id_aero), klev_src, klev, pplay_src, pplay, &
              1, klon, .FALSE.)

      IF (debug) CALL writefield_phy('day_tgt', var_day(:, :, id_aero), klev)

      ! c) adjust to conserve total aerosol mass load in the vertical pillar
      !    Calculate the load in the actual pillar and compare with the load
      !    read from aerosol file.

      ! Find the pressure difference in each model layer
      DO k = 1, klev
        DO i = 1, klon
          delp(i, k) = paprs(i, k) - paprs (i, k + 1)
        ENDDO
      ENDDO

      ! Find the mass load in the actual pillar, on target grid
      load_tgt(:) = 0.
      DO k = 1, klev
        DO i = 1, klon
          zrho = pplay(i, k) / t_seri(i, k) / RD       ! [kg/m3]
          volm = var_day(i, k, id_aero) * 1.E-9 / zrho ! [kg/kg]
          load_tgt(i) = load_tgt(i) + volm * delp(i, k) / RG
        ENDDO
      ENDDO

      ! Adjust, uniform
      DO k = 1, klev
        DO i = 1, klon
          var_day(i, k, id_aero) = var_day(i, k, id_aero) * load_src(i) / max(1.e-30, load_tgt(i))
        ENDDO
      ENDDO

      IF (debug) THEN
        load_tgt_test(:) = 0.
        DO k = 1, klev
          DO i = 1, klon
            zrho = pplay(i, k) / t_seri(i, k) / RD       ! [kg/m3]
            volm = var_day(i, k, id_aero) * 1.E-9 / zrho ! [kg/kg]
            load_tgt_test(i) = load_tgt_test(i) + volm * delp(i, k) / RG
          ENDDO
        ENDDO

        CALL writefield_phy('day_tgt2', var_day(:, :, id_aero), klev)
        CALL writefield_phy('load_tgt', load_tgt(:), 1)
        CALL writefield_phy('load_tgt_test', load_tgt_test(:), 1)
        CALL writefield_phy('load_src', load_src(:), 1)
      ENDIF

      ! - Interpolate variable tmp2 (source grid) to pi_var_day (target grid)
      !********************************************************************************
      ! a) calculate pression at vertical levels at source grid
      DO k = 1, klev_src
        DO i = 1, klon
          pplay_src(i, k) = pt_ap(k) + pt_b(k) * pi_psurf_day(i)
        ENDDO
      ENDDO

      IF (debug) THEN
        CALL writefield_phy('pi_psurf_day_src', pi_psurf_day(:), 1)
        CALL writefield_phy('pi_pplay_src', pplay_src(:, :), klev_src)
      ENDIF

      ! b) vertical interpolation on pressure leveles
      CALL pres2lev(tmp2(:, :), pi_var_day(:, :, id_aero), klev_src, klev, pplay_src, pplay, &
              1, klon, .FALSE.)

      IF (debug) CALL writefield_phy('pi_day_tgt', pi_var_day(:, :, id_aero), klev)

      ! c) adjust to conserve total aerosol mass load in the vertical pillar
      !    Calculate the load in the actual pillar and compare with the load
      !    read from aerosol file.

      ! Find the load in the actual pillar, on target grid
      load_tgt(:) = 0.
      DO k = 1, klev
        DO i = 1, klon
          zrho = pplay(i, k) / t_seri(i, k) / RD          ! [kg/m3]
          volm = pi_var_day(i, k, id_aero) * 1.E-9 / zrho ! [kg/kg]
          load_tgt(i) = load_tgt(i) + volm * delp(i, k) / RG
        ENDDO
      ENDDO

      DO k = 1, klev
        DO i = 1, klon
          pi_var_day(i, k, id_aero) = pi_var_day(i, k, id_aero) * pi_load_src(i) / max(1.e-30, load_tgt(i))
        ENDDO
      ENDDO

      IF (debug) THEN
        load_tgt_test(:) = 0.
        DO k = 1, klev
          DO i = 1, klon
            zrho = pplay(i, k) / t_seri(i, k) / RD          ! [kg/m3]
            volm = pi_var_day(i, k, id_aero) * 1.E-9 / zrho ! [kg/kg]
            load_tgt_test(i) = load_tgt_test(i) + volm * delp(i, k) / RG
          ENDDO
        ENDDO
        CALL writefield_phy('pi_day_tgt2', pi_var_day(:, :, id_aero), klev)
        CALL writefield_phy('pi_load_tgt', load_tgt(:), 1)
        CALL writefield_phy('pi_load_tgt_test', load_tgt_test(:), 1)
        CALL writefield_phy('pi_load_src', pi_load_src(:), 1)
      ENDIF

    ELSE   ! No vertical interpolation done

      var_day(:, :, id_aero) = tmp1(:, :)
      pi_var_day(:, :, id_aero) = tmp2(:, :)

    ENDIF ! vert_interp


    ! Deallocation
    DEALLOCATE(tmp1, tmp2, pplay_src, stat = ierr)

    !****************************************************************************************
    ! 4) Test for negative mass values

    !****************************************************************************************
    IF (MINVAL(var_day(:, :, id_aero)) < 0.) THEN
      DO k = 1, klev
        DO i = 1, klon
          ! Test for var_day
          IF (var_day(i, k, id_aero) < 0.) THEN
            IF (jDay - day2 < 0.) WRITE(lunout, *) 'jDay-day2=', jDay - day2
            IF (var_year(i, k, im2, id_aero) - var_year(i, k, im, id_aero) < 0.) THEN
              WRITE(lunout, *) trim(name_aero(id_aero)), '(i,k,im2)-', &
                      trim(name_aero(id_aero)), '(i,k,im)=', &
                      var_year(i, k, im2, id_aero) - var_year(i, k, im, id_aero)
            ENDIF
            WRITE(lunout, *) 'stop for aerosol : ', name_aero(id_aero)
            WRITE(lunout, *) 'day1, day2, jDay = ', day1, day2, jDay
            CALL abort_physic('readaerosol_interp', 'Error in interpolation 1', 1)
          ENDIF
        ENDDO
      ENDDO
    ENDIF

    IF (MINVAL(pi_var_day(:, :, id_aero)) < 0.) THEN
      DO k = 1, klev
        DO i = 1, klon
          ! Test for pi_var_day
          IF (pi_var_day(i, k, id_aero) < 0.) THEN
            IF (jDay - day2 < 0.) WRITE(lunout, *) 'jDay-day2=', jDay - day2
            IF (pi_var_year(i, k, im2, id_aero) - pi_var_year(i, k, im, id_aero) < 0.) THEN
              WRITE(lunout, *) trim(name_aero(id_aero)), '(i,k,im2)-', &
                      trim(name_aero(id_aero)), '(i,k,im)=', &
                      pi_var_year(i, k, im2, id_aero) - pi_var_year(i, k, im, id_aero)
            ENDIF

            WRITE(lunout, *) 'stop for aerosol : ', name_aero(id_aero)
            CALL abort_physic('readaerosol_interp', 'Error in interpolation 2', 1)
          ENDIF
        ENDDO
      ENDDO
    ENDIF

  ENDIF ! lnewday

  !****************************************************************************************
  ! Copy output from saved variables

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

  mass_out(:, :) = var_day(:, :, id_aero)
  pi_mass_out(:, :) = pi_var_day(:, :, id_aero)

END SUBROUTINE readaerosol_interp