Index: LMDZ5/trunk/libf/phylmd/clesphys.h
===================================================================
--- LMDZ5/trunk/libf/phylmd/clesphys.h	(revision 2786)
+++ LMDZ5/trunk/libf/phylmd/clesphys.h	(revision 2788)
@@ -87,4 +87,6 @@
        LOGICAL :: ok_qch4
        LOGICAL :: ok_conserv_q
+       LOGICAL :: adjust_tropopause
+       LOGICAL :: ok_daily_climoz
 
        COMMON/clesphys/                                                 &
@@ -130,5 +132,6 @@
      &     , iflag_rrtm, ok_strato,ok_hines, ok_qch4                    &
      &     , iflag_ice_thermo, ok_gwd_rando, NSW, iflag_albedo          &
-     &     , ok_chlorophyll,ok_conserv_q, ok_all_xml
+     &     , ok_chlorophyll,ok_conserv_q, adjust_tropopause             &
+     &     , ok_daily_climoz, ok_all_xml
      
        save /clesphys/
Index: LMDZ5/trunk/libf/phylmd/conf_phys_m.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/conf_phys_m.F90	(revision 2786)
+++ LMDZ5/trunk/libf/phylmd/conf_phys_m.F90	(revision 2788)
@@ -222,4 +222,6 @@
     LOGICAL,SAVE  :: carbon_cycle_tr_omp
     LOGICAL,SAVE  :: carbon_cycle_cpl_omp
+    LOGICAL,SAVE  :: adjust_tropopause_omp
+    LOGICAL,SAVE  :: ok_daily_climoz_omp
 
     INTEGER, INTENT(OUT):: read_climoz ! read ozone climatology, OpenMP shared
@@ -2035,4 +2037,24 @@
     !Config Help = ... 
     !
+    !
+    adjust_tropopause = .FALSE.
+    CALL getin('adjust_tropopause', adjust_tropopause_omp)
+    !
+    !Config Key  = adjust_tropopause
+    !Config Desc = Adjust the ozone field from the climoz file by stretching its
+    !              tropopause so that it matches the one of LMDZ.
+    !Config Def  = .FALSE.
+    !Config Help = Ensure tropospheric ozone column conservation.
+    !
+    !
+    ok_daily_climoz = .TRUE.
+    CALL getin('ok_daily_climoz', ok_daily_climoz_omp)
+    !
+    !Config Key  = ok_daily_climoz
+    !Config Desc = Interpolate in time the ozone forcings within ce0l.
+    !              .TRUE. if backward compatibility is needed.
+    !Config Def  = .TRUE.
+    !Config Help = .FALSE. ensure much fewer (no calendar dependency)
+    !  and lighter monthly climoz files, inetrpolated in time at gcm run time.
     !
     ecrit_LES_omp = 1./8.
@@ -2291,4 +2313,6 @@
     callstats = callstats_omp
     ecrit_LES = ecrit_LES_omp
+    adjust_tropopause = adjust_tropopause_omp
+    ok_daily_climoz = ok_daily_climoz_omp
     carbon_cycle_tr = carbon_cycle_tr_omp
     carbon_cycle_cpl = carbon_cycle_cpl_omp
@@ -2485,5 +2509,5 @@
     write(lunout,*)' aerosol_couple = ', aerosol_couple
     write(lunout,*)' flag_aerosol = ', flag_aerosol
-    write(lunout,*)' flag_aerosol_strat = ', flag_aerosol_strat
+    write(lunout,*)' flag_aerosol_strat= ', flag_aerosol_strat
     write(lunout,*)' new_aod = ', new_aod
     write(lunout,*)' aer_type = ',aer_type
@@ -2568,4 +2592,6 @@
     write(lunout,*) 'SSO gkwake=',gkwake
     write(lunout,*) 'SSO gklift=',gklift
+    write(lunout,*) 'adjust_tropopause = ', adjust_tropopause
+    write(lunout,*) 'ok_daily_climoz = ',ok_daily_climoz
     write(lunout,*) 'read_climoz = ', read_climoz
     write(lunout,*) 'carbon_cycle_tr = ', carbon_cycle_tr
Index: LMDZ5/trunk/libf/phylmd/limit_read_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/limit_read_mod.F90	(revision 2786)
+++ LMDZ5/trunk/libf/phylmd/limit_read_mod.F90	(revision 2788)
@@ -223,5 +223,5 @@
           ierr=NF90_INQUIRE(nid, UnlimitedDimID=ndimid)
           ierr=NF90_INQUIRE_DIMENSION(nid, ndimid, len=nn)
-          WRITE(abort_message,'(a,2(i0,a))')'limit.nc records number (',nn,') does no'//&
+          WRITE(abort_message,'(a,2(i3,a))')'limit.nc records number (',nn,') does no'//&
             't match year length (',year_len,')'
           IF(nn/=year_len) CALL abort_physic(modname,abort_message,1)
Index: LMDZ5/trunk/libf/phylmd/open_climoz_m.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/open_climoz_m.F90	(revision 2786)
+++ LMDZ5/trunk/libf/phylmd/open_climoz_m.F90	(revision 2788)
@@ -6,10 +6,10 @@
 contains
 
-  subroutine open_climoz(ncid, press_in_edg)
+  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 broadcasts them to the other processes.
+    ! levels and the times and broadcasts them to the other processes.
 
     ! We assume that, in "climoz_LMDZ.nc", the pressure levels are in hPa
@@ -21,24 +21,29 @@
     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
 
-    real, pointer:: press_in_edg(:)
-    ! edges of pressure intervals for ozone climatology, in Pa, in strictly
-    ! ascending order, 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:
 
-    real, pointer:: plev(:)
-    ! (pressure levels for ozone climatology, converted to Pa, in strictly
-    ! ascending order)
-
     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
 
     !---------------------------------------
 
-    print *, "Call sequence information: open_climoz"
+    modname="open_climoz"
+    print *, "Call sequence information: "//TRIM(modname)
 
     if (is_mpi_root) then
@@ -46,14 +51,27 @@
 
        call nf95_inq_varid(ncid, "plev", varid)
-       call nf95_gw_var(ncid, varid, plev)
+       call nf95_gw_var(ncid, varid, press_in_cen)
        ! Convert from hPa to Pa because "paprs" and "pplay" are in Pa:
-       plev = plev * 100.
-       n_plev = size(plev)
+       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(plev(n_plev))
-    call bcast_mpi(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))
@@ -61,5 +79,5 @@
        press_in_edg(1) = 0.
        ! We choose edges halfway in logarithm:
-       forall (k = 2:n_plev) press_in_edg(k) = sqrt(plev(k - 1) * plev(k))
+       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
@@ -67,5 +85,5 @@
     end if
     call bcast_mpi(press_in_edg)
-    deallocate(plev) ! pointer
+!    deallocate(press_in_cen) ! pointer
 
   end subroutine open_climoz
Index: LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 2786)
+++ LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 2788)
@@ -487,4 +487,6 @@
 !$OMP THREADPRIVATE(budg_sed_part)
 #endif
+      REAL, ALLOCATABLE, SAVE, DIMENSION(:) :: pr_tropopause
+!$OMP THREADPRIVATE(pr_tropopause)
 
 CONTAINS
@@ -764,4 +766,5 @@
       ALLOCATE (vsed_aer(klon,klev))
 #endif
+      ALLOCATE (pr_tropopause(klon))
 
 END SUBROUTINE phys_local_var_init
@@ -1018,4 +1021,5 @@
       DEALLOCATE (budg_sed_part)
 #endif
+      DEALLOCATE (pr_tropopause)
 
 END SUBROUTINE phys_local_var_end
Index: LMDZ5/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/physiq_mod.F90	(revision 2786)
+++ LMDZ5/trunk/libf/phylmd/physiq_mod.F90	(revision 2788)
@@ -35,4 +35,5 @@
     USE change_srf_frac_mod
     USE surface_data,     ONLY : type_ocean, ok_veget, ok_snow
+    USE tropopause_m,     ONLY: dyn_tropopause
 #ifdef CPP_Dust
     USE phytracr_spl_mod, ONLY: phytracr_spl
@@ -191,5 +192,6 @@
        beta_prec,  &
        rneb,  &
-       zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic
+       zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic, &
+       pr_tropopause
        !
     USE phys_state_var_mod ! Variables sauvegardees de la physique
@@ -205,5 +207,5 @@
     USE phys_output_ctrlout_mod
     use open_climoz_m, only: open_climoz ! ozone climatology from a file
-    use regr_pr_av_m, only: regr_pr_av
+    use regr_pr_time_av_m, only: regr_pr_time_av
     use netcdf95, only: nf95_close
     !IM for NMC files
@@ -1006,5 +1008,7 @@
     !$OMP THREADPRIVATE(first)
 
-    integer, save::  read_climoz ! read ozone climatology
+    ! VARIABLES RELATED TO OZONE CLIMATOLOGIES ; all are OpenMP shared
+    ! Note that pressure vectors are in Pa and in stricly ascending order
+    INTEGER,SAVE :: read_climoz                ! Read ozone climatology
     !     (let it keep the default OpenMP shared attribute)
     !     Allowed values are 0, 1 and 2
@@ -1013,16 +1017,11 @@
     !     2: read two ozone climatologies, the average day and night
     !     climatology and the daylight climatology
-
-    integer, save:: ncid_climoz ! NetCDF file containing ozone climatologies
-    !     (let it keep the default OpenMP shared attribute)
-
-    real, pointer, save:: press_climoz(:)
-    ! (let it keep the default OpenMP shared attribute)
-    ! edges of pressure intervals for ozone climatologies, in Pa, in strictly
-    ! ascending order
-
-    integer ro3i
-    !     required time index in NetCDF file for the ozone fields, between 1
-    !     and 360
+    INTEGER,SAVE :: ncid_climoz                ! NetCDF file identifier
+    REAL, POINTER, SAVE :: press_cen_climoz(:) ! Pressure levels
+    REAL, POINTER, SAVE :: press_edg_climoz(:) ! Edges of pressure intervals
+    REAL, POINTER, SAVE :: time_climoz(:)      ! Time vector
+    CHARACTER(LEN=13), PARAMETER :: vars_climoz(2) &
+                                  = ["tro3         ","tro3_daylight"]
+    ! vars_climoz(1:read_climoz): variables names in climoz file.
 
     INTEGER ierr
@@ -1671,7 +1670,6 @@
 
        !$omp single
-       IF (read_climoz >= 1) THEN
-          CALL open_climoz(ncid_climoz, press_climoz)
-       ENDIF
+       IF (read_climoz >= 1) CALL open_climoz(ncid_climoz, press_cen_climoz,   &
+                                         press_edg_climoz, time_climoz)
        !$omp end single
        !
@@ -1978,24 +1976,22 @@
           wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
        ELSE
-          ro3i = int((days_elapsed + jh_cur - jh_1jan) / year_len * 360.) + 1
-
-          IF (ro3i == 361) ro3i = 360
-          ! (This should never occur, except perhaps because of roundup
-          ! error. See documentation.)
-
-          IF (read_climoz == 1) THEN
-             CALL regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
-                  press_in_edg=press_climoz, paprs=paprs, v3=wo)
+          IF(adjust_tropopause) THEN
+            WRITE(*,*)' >> Adjusting O3 field to match gcm tropopause.'
+            CALL dyn_tropopause(t_seri, ztsol, paprs, pplay, presnivs, rot,  &
+                 pr_tropopause)
+            CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),    &
+                 days_elapsed+jh_cur-jh_1jan, press_edg_climoz, paprs, wo,   &
+                 time_climoz,  latitude_deg,  press_cen_climoz, pr_tropopause)
           ELSE
-             ! read_climoz == 2
-             CALL regr_pr_av(ncid_climoz, (/"tro3         ", &
-                  "tro3_daylight"/), julien=ro3i, press_in_edg=press_climoz, &
-                  paprs=paprs, v3=wo)
-          ENDIF
+            WRITE(*,*)' >> Interpolating O3 field directly on gcm levels.'
+            CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),    &
+                 days_elapsed+jh_cur-jh_1jan, press_edg_climoz, paprs, wo,   &
+                 time_climoz)
+          END IF
           ! Convert from mole fraction of ozone to column density of ozone in a
           ! cell, in kDU:
           FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
                * zmasse / dobson_u / 1e3
-          ! (By regridding ozone values for LMDZ only once per day, we
+          ! (By regridding ozone values for LMDZ only once a day, we
           ! have already neglected the variation of pressure in one
           ! day. So do not recompute "wo" at each time step even if
@@ -4465,5 +4461,6 @@
              CALL nf95_close(ncid_climoz)
           ENDIF
-          DEALLOCATE(press_climoz) ! pointer
+          DEALLOCATE(press_edg_climoz) ! pointer
+          DEALLOCATE(press_cen_climoz) ! pointer
        ENDIF
        !$OMP END MASTER
Index: LMDZ5/trunk/libf/phylmd/regr_horiz_time_climoz_m.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/regr_horiz_time_climoz_m.F90	(revision 2788)
+++ LMDZ5/trunk/libf/phylmd/regr_horiz_time_climoz_m.F90	(revision 2788)
@@ -0,0 +1,510 @@
+MODULE regr_horiz_time_climoz_m
+
+  USE interpolation,     ONLY: locate
+  USE mod_grid_phy_lmdz, ONLY: nlon_ou => nbp_lon, nlat_ou => nbp_lat
+  USE nrtype,            ONLY: pi
+  USE netcdf,   ONLY: NF90_CLOBBER, NF90_FLOAT,     NF90_GET_VAR, 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,  HANDLE_ERR,  &
+          NF95_CLOSE, NF95_ENDDEF,  NF95_PUT_ATT,   NF95_PUT_VAR, NF95_COPY_ATT
+  USE print_control_mod, ONLY: lunout
+  IMPLICIT NONE
+  PRIVATE
+  PUBLIC :: regr_horiz_time_climoz
+  REAL, PARAMETER :: deg2rad=pi/180.
+  CHARACTER(LEN=13), PARAMETER :: vars_in(2)=['tro3         ','tro3_daylight']
+
+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 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
+!-------------------------------------------------------------------------------
+! 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, POINTER     :: lon_in(:)           ! Longitudes   (ascending order, rad)
+  REAL, POINTER     :: lat_in(:)           ! Latitudes    (ascending order, rad)
+  REAL, POINTER     :: 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_in2  (:,:,:,:)   ! 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_regr_lat     (:,:,:,:) !         (nlat_in,:,0:13   ,:)
+  REAL, ALLOCATABLE :: o3_out2         (:,:,:,:) !         (nlat_ou,:,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, POINTER :: 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
+!-------------------------------------------------------------------------------
+  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
+
+  !--- 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
+    !--- Compute input edges longitudes vector (no end point yet)
+    ALLOCATE(v1(nlon_in))
+    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.
+    DEALLOCATE(lon_in)
+
+    !--- Shift input longitudes vector until it contains first output point boundslon_reg(1,east)
+    v1=v1+2*pi*REAL(FLOOR((boundslon_reg(1,east)-v1(1))/(2.*pi)))
+
+    !--- Ensure first input longitudes interval contains first output point boundslon_reg(1,east)
+    dx1=locate(v1,boundslon_reg(1,east))-1
+    v1=CSHIFT(v1,SHIFT=dx1,DIM=1); v1(nlon_in-dx1+1:)=v1(nlon_in-dx1+1:)+2.*pi
+
+    !--- Extend input longitudes vector until last interval contains boundslon_reg(nlat_ou,west)
+    dx2=0; DO WHILE(v1(1+dx2)+2.*pi<boundslon_reg(nlat_ou,west)); 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(1:1+dx2)+2.*pi]
+    DEALLOCATE(v1)
+  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)
+
+  !--- Prepare quantities for time interpolation
+  tmidmonth=mid_month(anneeref, cal_in)
+  IF(interpt) THEN
+    ntim_ou=ioget_year_len(anneeref)
+    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
+
+  !--- 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)
+
+  !--- 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)
+
+  !--- 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
+    msg=TRIM(sub)//" NF90_GET_VAR "//TRIM(vars_in(iv))
+    CALL NF95_INQ_VARID(fID_in, vars_in(1), vID_in(iv))
+    IF(l3D) ncerr=NF90_GET_VAR(fID_in, vID_in(iv), o3_in3(1:nlon_in,:,:,1:12,iv))
+    IF(l2D) ncerr=NF90_GET_VAR(fID_in, vID_in(iv), o3_in2(          :,:,1:12,iv))
+    CALL handle_err(TRIM(msg), ncerr, fID_in)
+    IF(lprev) THEN; sta(ndims)=12
+      CALL NF95_INQ_VARID(fID_in_m, vars_in(1), vID_in(iv))
+      IF(l3D) ncerr=NF90_GET_VAR(fID_in_m,vID_in(iv),o3_in3(1:nlon_in,:,:, 0,iv),sta,cnt)
+      IF(l2d) ncerr=NF90_GET_VAR(fID_in_m,vID_in(iv),o3_in2(          :,:, 0,iv),sta,cnt)
+      CALL handle_err(TRIM(msg)//" previous", ncerr, fID_in_m)
+    END IF
+    IF(lnext) THEN; sta(ndims)=1
+      CALL NF95_INQ_VARID(fID_in_p, vars_in(1), vID_in(iv))
+      IF(l3D) ncerr=NF90_GET_VAR(fID_in_p,vID_in(iv),o3_in3(1:nlon_in,:,:,13,iv),sta,cnt)
+      IF(l2D) ncerr=NF90_GET_VAR(fID_in_p,vID_in(iv),o3_in2(          :,:,13,iv),sta,cnt)
+      CALL handle_err(TRIM(msg)//" next", ncerr, fID_in_p)
+    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,:,:)
+    !--- 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,:,:,:,:)
+  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
+  CALL NF95_CLOSE(fID_in)
+
+  !=============================================================================
+  IF(l3D) THEN                                                   !=== 3D FIELDS
+  !=============================================================================
+    !--- 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,east),boundslon_reg(:,west)],    &
+                        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)
+
+    !--- 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)
+
+    !--- 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_out3(:,nlat_ou:1:-1,:,:,m))
+    END DO
+
+  !=============================================================================
+  ELSE                                                         !=== ZONAL FIELDS
+  !=============================================================================
+    !--- 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,:)
+
+    !--- 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)
+
+    !--- 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
+
+  !=============================================================================
+  END IF
+  !=============================================================================
+
+  CALL NF95_CLOSE(fID_ou)
+
+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
+!-------------------------------------------------------------------------------
+! 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
+  CHARACTER(LEN=80) :: sub
+!-------------------------------------------------------------------------------
+  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)
+
+  !--- 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)
+
+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
+!
+!-------------------------------------------------------------------------------
Index: LMDZ5/trunk/libf/phylmd/regr_lat_time_climoz_m.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/regr_lat_time_climoz_m.F90	(revision 2786)
+++ 	(revision )
@@ -1,450 +1,0 @@
-! $Id$
-module regr_lat_time_climoz_m
-
-  ! Author: Lionel GUEZ
-
-  implicit none
-
-  private
-  public regr_lat_time_climoz
-
-contains
-
-  subroutine regr_lat_time_climoz(read_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.
-
-    ! 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.
-
-    ! We assume that in the input file:
-
-    ! -- Latitude is in degrees.
-
-    ! -- Latitude and pressure are strictly monotonic (as all NetCDF
-    ! coordinate variables should be).
-
-    ! -- The time coordinate is in ascending order (even though we do
-    ! not use its values).
-    ! 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
-
-    ! -- Missing values are contiguous, at the bottom of
-    ! the vertical domain and at the latitudinal boundaries.
-
-    ! If values are all missing at a given latitude and date, then we
-    ! replace those missing values by values at the closest latitude,
-    ! equatorward, with valid values.
-    ! Then, at each latitude and each date, the missing values are replaced
-    ! by the lowest valid value above missing values.
-
-    ! Regridding in time is by linear interpolation.
-    ! Monthly values are processed to get daily values, on the basis
-    ! of a 360-day calendar.
-    ! If there are 14 months, 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 then we assume
-    ! periodicity for interpolation at the beginning and at the end of the
-    ! year.
-
-    use mod_grid_phy_lmdz, ONLY : nbp_lat
-    use regr1_conserv_m, only: regr1_conserv
-    use regr3_lint_m, only: regr3_lint
-    use netcdf95, only: handle_err, nf95_close, nf95_get_att, nf95_gw_var, &
-         nf95_inq_dimid, nf95_inq_varid, nf95_inquire_dimension, nf95_open, &
-         nf95_put_var
-    use netcdf, only: nf90_get_att, nf90_get_var, nf90_noerr, nf90_nowrite
-    use assert_m, only: assert
-    use regular_lonlat_mod, only : boundslat_reg, south
-    use nrtype, only: pi
-    use slopes_m, only: slopes
-
-    integer, intent(in):: read_climoz ! read ozone climatology
-    ! Allowed values are 1 and 2
-    ! 1: read a single ozone climatology that will be used day and night
-    ! 2: read two ozone climatologies, the average day and night
-    ! climatology and the daylight climatology
-
-    ! Variables local to the procedure:
-
-    integer n_plev ! number of pressure levels in the input data
-    integer n_lat ! number of latitudes in the input data
-    integer n_month ! number of months in the input data
-
-    real, pointer:: latitude(:)
-    ! (of input data, converted to rad, sorted in strictly ascending order)
-
-    real, allocatable:: sin_lat_in_edg(:)
-    ! (sine of 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, converted to hPa
-
-    logical desc_lat ! latitude in descending order in the input file
-    logical desc_plev ! pressure levels in descending order in the input file
-
-    real, allocatable:: o3_in(:, :, :, :)
-    ! (n_lat, n_plev, n_month, read_climoz)
-    ! ozone climatologies from the input file
-    ! "o3_in(j, k, :, :)" is at latitude "latitude(j)" and pressure
-    ! level "plev(k)".
-    ! Third dimension is month index, first value may be December or January.
-    ! "o3_in(:, :, :, 1)" is for the day- night average, "o3_in(:, :, :, 2)"
-    ! is for daylight.
-
-    real missing_value
-
-    real, allocatable:: o3_regr_lat(:, :, :, :)
-    ! (nbp_lat, n_plev, 0:13, read_climoz)
-    ! 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 "nbp_lat - 1" then the interval is:
-    ! [rlatv(j), rlatv(j-1)]
-    ! If "j" is 1 or "nbp_lat" then the interval is:
-    ! [rlatv(1), pi / 2]
-    ! or:
-    ! [- pi / 2, rlatv(nbp_lat - 1)]
-    ! respectively.
-    ! "o3_regr_lat(:, k, :, :)" is for pressure level "plev(k)".
-    ! Third dimension is month number, 1 for January.
-    ! "o3_regr_lat(:, :, :, 1)" is average day and night,
-    ! "o3_regr_lat(:, :, :, 2)" is for daylight.
-
-    real, allocatable:: o3_out(:, :, :, :)
-    ! (nbp_lat, n_plev, 360, read_climoz)
-    ! 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.
-    ! "o3_out(:, :, :, 1)" is average day and night,
-    ! "o3_out(:, :, :, 2)" is for daylight.
-
-    integer j, k, l,m
-
-    ! For NetCDF:
-    integer ncid_in, ncid_out ! IDs for input and output files
-    integer varid_plev, varid_time, varid, ncerr, dimid
-    character(len=80) press_unit ! pressure unit
-
-    integer varid_in(read_climoz), varid_out(read_climoz)
-    ! index 1 is for average ozone day and night, index 2 is for
-    ! daylight ozone.
-
-    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 assert(read_climoz == 1 .or. read_climoz == 2, "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_conserv" 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(sin_lat_in_edg(n_lat + 1))
-    sin_lat_in_edg(1) = - 1.
-    forall (j = 2:n_lat) sin_lat_in_edg(j) = sin((latitude(j - 1) &
-         + latitude(j)) / 2.)
-    sin_lat_in_edg(n_lat + 1) = 1.
-    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)
-    call nf95_get_att(ncid_in, varid, "units", press_unit)
-    if (press_unit == "Pa") then
-       ! Convert to hPa:
-       plev = plev / 100.
-    elseif (press_unit /= "hPa") then
-       print *, "regr_lat_time_climoz: the only recognized units are Pa " &
-            // "and hPa."
-       stop 1
-    end if
-
-    ! 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
-
-    ! Get the  number of months:
-    call nf95_inq_dimid(ncid_in, "time", dimid)
-    call nf95_inquire_dimension(ncid_in, dimid, nclen=n_month)
-
-    allocate(o3_in(n_lat, n_plev, n_month, read_climoz))
-
-    call nf95_inq_varid(ncid_in, "tro3", varid_in(1))
-    ncerr = nf90_get_var(ncid_in, varid_in(1), o3_in(:, :, :, 1))
-    call handle_err("regr_lat_time_climoz nf90_get_var tro3", ncerr, ncid_in)
-
-    if (read_climoz == 2) then
-       call nf95_inq_varid(ncid_in, "tro3_daylight", varid_in(2))
-       ncerr = nf90_get_var(ncid_in, varid_in(2), o3_in(:, :, :, 2))
-       call handle_err("regr_lat_time_climoz nf90_get_var tro3_daylight", &
-            ncerr, ncid_in, varid_in(2))
-    end if
-
-    if (desc_lat) o3_in = o3_in(n_lat:1:-1, :, :, :)
-    if (desc_plev) o3_in = o3_in(:, n_plev:1:-1, :, :)
-
-    do m = 1, read_climoz
-       ncerr = nf90_get_att(ncid_in, varid_in(m), "missing_value", &
-            missing_value)
-       if (ncerr == nf90_noerr) then
-          do l = 1, n_month
-             ! Take care of latitudes where values are all missing:
-
-             ! Next to the south pole:
-             j = 1
-             do while (o3_in(j, 1, l, m) == missing_value)
-                j = j + 1
-             end do
-             if (j > 1) o3_in(:j-1, :, l, m) = &
-                  spread(o3_in(j, :, l, m), dim=1, ncopies=j-1)
-             
-             ! Next to the north pole:
-             j = n_lat
-             do while (o3_in(j, 1, l, m) == missing_value)
-                j = j - 1
-             end do
-             if (j < n_lat) o3_in(j+1:, :, l, m) = &
-                  spread(o3_in(j, :, l, m), dim=1, ncopies=n_lat-j)
-
-             ! Take care of missing values at high pressure:
-             do j = 1, n_lat
-                ! Find missing values, starting from top of atmosphere
-                ! and going down.
-                ! We have already taken care of latitudes full of
-                ! missing values so the highest level has a valid value.
-                k = 2
-                do while  (o3_in(j, k, l, m) /= missing_value .and. k < n_plev)
-                   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, m) == missing_value) &
-                     o3_in(j, k:n_plev, l, m) = o3_in(j, k-1, l, m)
-             end do
-          end do
-       else
-          print *, "regr_lat_time_climoz: field ", m, &
-               ", no missing value attribute"
-       end if
-    end do
-
-    call nf95_close(ncid_in)
-
-    allocate(o3_regr_lat(nbp_lat, n_plev, 0:13, read_climoz))
-    allocate(o3_out(nbp_lat, n_plev, 360, read_climoz))
-
-    ! Regrid in latitude:
-    ! We average with respect to sine of latitude, which is
-    ! equivalent to weighting by cosine of latitude:
-    if (n_month == 12) then
-       print *, &
-            "Found 12 months in ozone climatologies, assuming periodicity..."
-       call regr1_conserv(o3_in, xs = sin_lat_in_edg, &
-            xt = (/- 1., sin(boundslat_reg(nbp_lat - 1:1:- 1, south)), 1./), &
-            vt = o3_regr_lat(nbp_lat:1:- 1, :, 1:12, :), &
-            slope = slopes(o3_in, sin_lat_in_edg))
-       ! (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 climatologies..."
-       call regr1_conserv(o3_in, xs = sin_lat_in_edg, &
-            xt = (/- 1., sin(boundslat_reg(nbp_lat - 1:1:- 1, south)), 1./), &
-            vt = o3_regr_lat(nbp_lat:1:- 1, :, :, :), &
-            slope = slopes(o3_in, sin_lat_in_edg))
-       ! (invert order of indices in "o3_regr_lat" because "rlatu" is
-       ! in descending order)
-    end if
-
-    ! Regrid in time by linear interpolation:
-    o3_out = regr3_lint(o3_regr_lat, tmidmonth, tmidday)
-
-    ! Write to file:
-    do m = 1, read_climoz
-       call nf95_put_var(ncid_out, varid_out(m), o3_out(nbp_lat:1:-1, :, :, m))
-       ! (The order of "rlatu" is inverted in the output file)
-    end do
-
-    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 mod_grid_phy_lmdz, ONLY : nbp_lat
-    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
-    use nrtype, only: pi
-    use regular_lonlat_mod, only : lat_reg
-
-    integer, intent(in):: ncid_in, n_plev
-    integer, intent(out):: ncid_out, varid_plev, varid_time
-
-    integer, intent(out):: varid_out(:) ! dim(1 or 2)
-    ! "varid_out(1)" is for average ozone day and night,
-    ! "varid_out(2)" is for daylight ozone.
-
-    ! Variables local to the procedure:
-
-    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", nbp_lat, 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 variables:
-
-    call nf95_def_var(ncid_out, "tro3", nf90_float, &
-         (/dimid_rlatu, dimid_plev, dimid_time/), varid_out(1))
-    call nf95_put_att(ncid_out, varid_out(1), "long_name", &
-         "ozone mole fraction")
-    call nf95_put_att(ncid_out, varid_out(1), "standard_name", &
-         "mole_fraction_of_ozone_in_air")
-
-    if (size(varid_out) == 2) then
-       call nf95_def_var(ncid_out, "tro3_daylight", nf90_float, &
-            (/dimid_rlatu, dimid_plev, dimid_time/), varid_out(2))
-       call nf95_put_att(ncid_out, varid_out(2), "long_name", &
-            "ozone mole fraction in daylight")
-    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(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, lat_reg(nbp_lat: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
Index: LMDZ5/trunk/libf/phylmd/regr_lat_time_coefoz_m.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/regr_lat_time_coefoz_m.F90	(revision 2786)
+++ LMDZ5/trunk/libf/phylmd/regr_lat_time_coefoz_m.F90	(revision 2788)
@@ -41,6 +41,6 @@
 
     use mod_grid_phy_lmdz, ONLY : nbp_lat
-    use regr1_conserv_m, only: regr1_conserv
-    use regr3_lint_m, only: regr3_lint
+    use regr_conserv_m, only: regr_conserv
+    use regr_lint_m, only: regr_lint
     use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, handle_err, &
          nf95_put_var, nf95_gw_var
@@ -84,5 +84,5 @@
     ! Last dimension is month number.)
 
-    real, allocatable:: o3_par_out(:, :, :) ! (jjm + 1, n_plev, 360)
+    real, allocatable:: o3_par_out(:, :, :) ! (jjm + 1, n_plev, ndays)
     ! (regridded ozone parameter)
     ! ("o3_par_out(j, l, day)" is at latitude "rlatu(j)", pressure
@@ -162,5 +162,5 @@
     latitude = latitude / 180. * pi
     n_lat = size(latitude)
-    ! We need to supply the latitudes to "regr1_conserv" in
+    ! We need to supply the latitudes to "regr_conserv" in
     ! ascending order, so invert order if necessary:
     desc_lat = latitude(1) > latitude(n_lat)
@@ -209,5 +209,5 @@
        ! We average with respect to sine of latitude, which is
        ! equivalent to weighting by cosine of latitude:
-       call regr1_conserv(o3_par_in, xs = sin(lat_in_edg), &
+       call regr_conserv(1, o3_par_in, xs = sin(lat_in_edg), &
             xt = (/-1., sin((/boundslat_reg(nbp_lat-1:1:-1,south)/)), 1./), &
             vt = v_regr_lat(nbp_lat:1:-1, :, 1:12))
@@ -221,5 +221,5 @@
 
        ! Regrid in time by linear interpolation:
-       o3_par_out = regr3_lint(v_regr_lat, tmidmonth, tmidday)
+       call regr_lint(3, v_regr_lat, tmidmonth, tmidday, o3_par_out)
 
        ! Write to file:
Index: LMDZ5/trunk/libf/phylmd/regr_pr_av_m.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/regr_pr_av_m.F90	(revision 2786)
+++ 	(revision )
@@ -1,124 +1,0 @@
-! $Id$
-module regr_pr_av_m
-
-  ! Author: Lionel GUEZ
-
-  implicit none
-
-contains
-
-  subroutine regr_pr_av(ncid, name, julien, press_in_edg, paprs, v3)
-
-    ! "regr_pr_av" stands for "regrid pressure averaging".
-    ! In this procedure:
-    ! -- the root process reads 2D latitude-pressure fields from a
-    !    NetCDF file, at a given day.
-    ! -- the fields are packed to the LMDZ horizontal "physics"
-    !    grid and scattered to all threads of all processes;
-    ! -- in all the threads of all the processes, the fields are regridded in
-    !    pressure to the LMDZ vertical grid.
-    ! We assume that, in the input file, the fields have 3 dimensions:
-    ! latitude, pressure, julian day.
-    ! We assume that the input fields are already on the "rlatu"
-    ! latitudes, except that latitudes are in ascending order in the input
-    ! file.
-    ! We assume that all the inputs fields have the same coordinates.
-
-    ! The target vertical LMDZ grid is the grid of layer boundaries.
-    ! Regridding in pressure is conservative, second order.
-
-    ! All the fields are regridded as a single multi-dimensional array
-    ! so it saves CPU time to call this procedure once for several NetCDF
-    ! variables rather than several times, each time for a single
-    ! NetCDF variable.
-
-    use dimphy, only: klon
-    use netcdf95, only: nf95_inq_varid, handle_err
-    use netcdf, only: nf90_get_var
-    use assert_m, only: assert
-    use assert_eq_m, only: assert_eq
-    use regr1_conserv_m, only: regr1_conserv
-    use slopes_m, only: slopes
-    use mod_phys_lmdz_mpi_data, only: is_mpi_root
-    use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat, nbp_lev
-    use mod_phys_lmdz_transfert_para, only: scatter2d
-    ! (pack to the LMDZ horizontal "physics" grid and scatter)
-
-    integer, intent(in):: ncid ! NetCDF ID of the file
-    character(len=*), intent(in):: name(:) ! of the NetCDF variables
-    integer, intent(in):: julien ! jour julien, 1 <= julien <= 360
-
-    real, intent(in):: press_in_edg(:)
-    ! edges of pressure intervals for input data, in Pa, in strictly
-    ! ascending order
-
-    real, intent(in):: paprs(:, :) ! (klon, llm + 1)
-    ! (pression pour chaque inter-couche, en Pa)
-
-    real, intent(out):: v3(:, :, :) ! (klon, llm, size(name))
-    ! regridded fields on the partial "physics" grid
-    ! "v3(i, k, l)" is at longitude "xlon(i)", latitude
-    ! "xlat(i)", in pressure interval "[paprs(i, k+1), paprs(i, k)]",
-    ! for NetCDF variable "name(l)".
-
-    ! Variables local to the procedure:
-
-    integer varid, ncerr ! for NetCDF
-
-    real  v1(nbp_lon, nbp_lat, size(press_in_edg) - 1, size(name))
-    ! input fields at day "julien", on the global "dynamics" horizontal grid
-    ! First dimension is for longitude.
-    ! The values are the same for all longitudes.
-    ! "v1(:, j, k, l)" is at latitude "rlatu(j)", for
-    ! pressure interval "[press_in_edg(k), press_in_edg(k+1)]" and
-    ! NetCDF variable "name(l)".
-
-    real v2(klon, size(press_in_edg) - 1, size(name))
-    ! fields scattered to the partial "physics" horizontal grid
-    ! "v2(i, k, l)" is at longitude "xlon(i)", latitude "xlat(i)",
-    ! for pressure interval "[press_in_edg(k), press_in_edg(k+1)]" and
-    ! NetCDF variable "name(l)".
-
-    integer i, n_var
-
-    !--------------------------------------------
-
-    call assert(size(v3, 1) == klon, size(v3, 2) == nbp_lev, &
-         "regr_pr_av v3 klon")
-    n_var = assert_eq(size(name), size(v3, 3), "regr_pr_av v3 n_var")
-    call assert(shape(paprs) == (/klon, nbp_lev+1/), "regr_pr_av paprs")
-
-    !$omp master
-    if (is_mpi_root) then
-       do i = 1, n_var
-          call nf95_inq_varid(ncid, trim(name(i)), varid)
-          
-          ! Get data at the right day from the input file:
-          ncerr = nf90_get_var(ncid, varid, v1(1, :, :, i), &
-               start=(/1, 1, julien/))
-          call handle_err("regr_pr_av nf90_get_var " // trim(name(i)), ncerr, &
-               ncid)
-       end do
-       
-       ! Latitudes are in ascending order in the input file while
-       ! "rlatu" is in descending order so we need to invert order:
-       v1(1, :, :, :) = v1(1, nbp_lat:1:-1, :, :)
-
-       ! Duplicate on all longitudes:
-       v1(2:, :, :, :) = spread(v1(1, :, :, :), dim=1, ncopies=nbp_lon-1)
-    end if
-    !$omp end master
-
-    call scatter2d(v1, v2)
-
-    ! Regrid in pressure at each horizontal position:
-    do i = 1, klon
-       call regr1_conserv(v2(i, :, :), press_in_edg, &
-            paprs(i, nbp_lev + 1:1:-1), v3(i, nbp_lev:1:-1, :), &
-            slopes(v2(i, :, :), press_in_edg))
-       ! (invert order of indices because "paprs" is in descending order)
-    end do
-
-  end subroutine regr_pr_av
-
-end module regr_pr_av_m
Index: LMDZ5/trunk/libf/phylmd/regr_pr_comb_coefoz_m.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/regr_pr_comb_coefoz_m.F90	(revision 2786)
+++ LMDZ5/trunk/libf/phylmd/regr_pr_comb_coefoz_m.F90	(revision 2788)
@@ -76,5 +76,5 @@
     use dimphy, only: klon
     use mod_phys_lmdz_mpi_data, only: is_mpi_root
-    use regr_pr_av_m, only: regr_pr_av
+    use regr_pr_time_av_m, only: regr_pr_time_av
     use regr_pr_int_m, only: regr_pr_int
     use press_coefoz_m, only: press_in_edg, plev
@@ -128,6 +128,6 @@
     !$omp end master
 
-    call regr_pr_av(ncid, (/"a2       ", "a4       ", "a6       ", &
-         "P_net_Mob", "r_Mob    ", "temp_Mob ", "R_Het    "/), julien, &
+    call regr_pr_time_av(ncid, (/"a2       ", "a4       ", "a6       ", &
+         "P_net_Mob", "r_Mob    ", "temp_Mob ", "R_Het    "/), REAL(julien-1), &
          press_in_edg, paprs, coefoz)
     a2 = coefoz(:, :, 1)
Index: LMDZ5/trunk/libf/phylmd/regr_pr_int_m.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/regr_pr_int_m.F90	(revision 2786)
+++ LMDZ5/trunk/libf/phylmd/regr_pr_int_m.F90	(revision 2788)
@@ -28,5 +28,5 @@
     use netcdf, only: nf90_get_var
     use assert_m, only: assert
-    use regr1_lint_m, only: regr1_lint
+    use regr_lint_m, only: regr_lint
     use mod_phys_lmdz_mpi_data, only: is_mpi_root
     use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat, nbp_lev
@@ -97,5 +97,6 @@
     ! Regrid in pressure at each horizontal position:
     do i = 1, klon
-       v3(i, nbp_lev:1:-1) = regr1_lint(v2(i, :), (/0., plev/), pplay(i, nbp_lev:1:-1))
+       call regr_lint(1, v2(i, :), (/0., plev/), pplay(i, nbp_lev:1:-1), &
+                         v3(i, nbp_lev:1:-1))
        ! (invert order of indices because "pplay" is in descending order)
     end do
Index: LMDZ5/trunk/libf/phylmd/regr_pr_o3_m.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/regr_pr_o3_m.F90	(revision 2786)
+++ LMDZ5/trunk/libf/phylmd/regr_pr_o3_m.F90	(revision 2788)
@@ -28,5 +28,5 @@
     use netcdf, only:  nf90_nowrite, nf90_get_var
     use assert_m, only: assert
-    use regr1_conserv_m, only: regr1_conserv
+    use regr_conserv_m, only: regr_conserv
     use press_coefoz_m, only: press_in_edg
     use time_phylmdz_mod, only: day_ref
@@ -75,5 +75,5 @@
     ! Poles:
     do j = 1, nbp_lat, nbp_lat-1
-       call regr1_conserv(r_mob(j, :), press_in_edg, &
+       call regr_conserv(1, r_mob(j, :), press_in_edg, &
             p3d(1, j, nbp_lev + 1:1:-1), o3_mob_regr(1, j, nbp_lev:1:-1))
        ! (invert order of indices because "p3d" is in descending order)
@@ -83,5 +83,5 @@
     do j = 2, nbp_lat-1
        do i = 1, nbp_lon
-          call regr1_conserv(r_mob(j, :), press_in_edg, &
+          call regr_conserv(1, r_mob(j, :), press_in_edg, &
                p3d(i, j, nbp_lev + 1:1:-1), o3_mob_regr(i, j, nbp_lev:1:-1))
           ! (invert order of indices because "p3d" is in descending order)
Index: LMDZ5/trunk/libf/phylmd/regr_pr_time_av_m.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/regr_pr_time_av_m.F90	(revision 2788)
+++ LMDZ5/trunk/libf/phylmd/regr_pr_time_av_m.F90	(revision 2788)
@@ -0,0 +1,390 @@
+! $Id$
+MODULE regr_pr_time_av_m
+
+  USE write_field_phy
+  USE mod_phys_lmdz_transfert_para, ONLY: bcast
+  IMPLICIT NONE
+
+!-------------------------------------------------------------------------------
+! Purpose: Regrid pressure with an averaging method. Operations done:
+!  * on the root process: read a NetCDF 3D or 4D field at a given day.
+!  * pack the fields to the LMDZ "horizontal "physics" grid and scatter
+!    to all threads of all processes;
+!  * in all the threads of all the processes, regrid the fields in pressure
+!    to the LMDZ vertical grid.
+!  * input files fields are stretched in the viscinity (+/- 5 kms) of the mean
+!    tropopause (geometrical mean of LMDZ and input fields tropopauses) to force
+!    the resulting field tropopause height to the one of LMDZ. To switch this
+!    feature on, the following arguments must be present:
+!     - "lat_in":   input file latitudes.
+!     - "pcen_in":  input file cells center pressure levels.
+!     - "ptrop_ou": target grid (LMDZ) tropopause pressure.
+!    Note that the ozone quantity conservation is not ensured.
+!-------------------------------------------------------------------------------
+! Initial routine: regr_pr_av_m module (L. Guez).
+! Present version: David Cugnet ; corresponding additions:
+!    - time interpolation
+!    - 3D compliant
+!    - vertical stretching of the field to allow tropopause and ground pressure
+!      matching between input field and current lmdz field.
+!-------------------------------------------------------------------------------
+! Remarks:
+!  * 3D fields are zonal means, with dimensions (latitude, pressure, julian day)
+!  * 4D fields have the dimensions:  (longitude, latitude, pressure, julian day)
+!  * All the fields are already on the horizontal grid (rlatu) or (rlatu,rlonv),
+!    except that the latitudes are in ascending order in the input file.
+!  * We assume that all the input fields have the same coordinates.
+!  * The target vertical LMDZ grid is the grid of layer boundaries.
+!  * Regridding in pressure is conservative, second order.
+!  * All the fields are regridded as a single multi-dimensional array, so it
+!    saves CPU time to call this procedure once for several NetCDF variables
+!    rather than several times, each time for a single NetCDF variable.
+!  * The input file pressure at tropopause can be (in decreasing priority):
+!    1) read from the file if "tropopause_air_pressure" field is available.
+!    2) computed from the input file ozone field using:
+!       - o3 concentration at tropopause if "tro3_at_tropopause" is available.
+!       - a default value (100ppbv) if not.
+!  * If available, the field "ps" (input file pressure at surface) is used.
+!    If not, the current LMDZ ground pressure is taken instead.
+!  * The O3 threshold for tropopause is defined using 3 coefficients:
+!    o3t(ppbV)=co1+co2*SIN(PI*(month-2)/6)*TANH(lat_deg/c)
+!     => co1 and co2 are in ppmV, and co3 in degrees.
+!     => co3 allow a smooth transition between north and south hemispheres.
+!-------------------------------------------------------------------------------
+! * Fields with suffix "m"/"p" are at the closest records earlier/later than the
+!   middle of the julian day "julien", on the global "dynamics" horizontal grid.
+! * Fields(i,j,k,l) are at longitude-latitude "rlonv(i)-rlatu(j)", pressure
+!   interval "pint_in(k:k+1)]" and variable "nam(l)".
+! * In the 2D file case, the values are the same for all longitudes.
+! * The following fields are on the global "dynamics" grid, as read from files:
+  REAL,    SAVE, ALLOCATABLE :: v1m(:,:,:,:)       !--- Previous ozone fields
+  REAL,    SAVE, ALLOCATABLE :: v1p(:,:,:,:)       !--- Next ozone fields
+  REAL,    SAVE, ALLOCATABLE :: psm(:,:), psp(:,:) !--- Surface pressure
+  REAL,    SAVE, ALLOCATABLE :: ptm(:,:), ptp(:,:) !--- Tropopause pressure
+  REAL,    SAVE, ALLOCATABLE :: otm(:,:), otp(:,:) !--- Tropopause o3 mix. ratio
+  INTEGER, SAVE :: irec                            !--- Current time index
+!      * for daily   input files: current julian day number
+!      * for monthly input files: julien is in [time_in(irec),time_in(irec+1)]
+  LOGICAL, SAVE :: linterp                         !--- Interpolation in time
+  LOGICAL, SAVE :: lPrSurf                         !--- Surface pressure flag
+  LOGICAL, SAVE :: lPrTrop                         !--- Tropopause pressure flag
+  LOGICAL, SAVE :: lO3Trop                         !--- Tropopause ozone flag
+  LOGICAL, SAVE :: lfirst=.TRUE.                   !--- First call flag
+  REAL, PARAMETER :: co3(3)=[91.,28.,20.]          !--- Coeffs for o3 threshold
+  REAL, PARAMETER :: prtrop=5.E+3  !--- Value lower than the tropopause pressure
+  REAL, PARAMETER :: delta =5.     !--- Dist. around tropopause for strain (kms)
+
+CONTAINS
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE regr_pr_time_av(fID, nam, julien, pint_in, pcen_ou, v3, time_in,    &
+                                     lat_in, pcen_in, ptrop_ou)
+!
+!-------------------------------------------------------------------------------
+  USE dimphy,         ONLY: klon
+  USE netcdf95,       ONLY: NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, handle_err
+  USE netcdf,         ONLY: NF90_INQ_VARID, NF90_GET_VAR, NF90_NOERR
+  USE assert_m,       ONLY: assert
+  USE assert_eq_m,    ONLY: assert_eq
+  USE comvert_mod,    ONLY: scaleheight
+  USE interpolation,  ONLY: locate
+  USE regr_conserv_m, ONLY: regr_conserv
+  USE regr_lint_m,    ONLY: regr_lint
+  USE slopes_m,       ONLY: slopes
+  USE mod_phys_lmdz_mpi_data,       ONLY: is_mpi_root
+  USE mod_grid_phy_lmdz,            ONLY: nbp_lon, nbp_lat, nbp_lev
+  USE mod_phys_lmdz_transfert_para, ONLY: scatter2d, scatter
+  USE phys_cal_mod,                 ONLY: calend, year_len, days_elapsed, jH_cur
+!-------------------------------------------------------------------------------
+! Arguments:
+  INTEGER, INTENT(IN)  :: fID           !--- NetCDF file ID
+  CHARACTER(LEN=13), INTENT(IN) :: nam(:)     !--- NetCDF variables names
+  REAL,    INTENT(IN)  :: julien        !--- Days since Jan 1st
+  REAL,    INTENT(IN)  :: pint_in(:)    !--- Interfaces Pressure, Pa, ascending
+  REAL,    INTENT(IN)  :: pcen_ou(:,:)  !--- Mid-layers LMDZ Pres, Pa (klon,llm+1)
+  REAL,    INTENT(OUT) :: v3(:,:,:)     !--- Regridded field (klon,llm,SIZE(nam))
+  REAL,    INTENT(IN), OPTIONAL :: time_in(:) !--- Records times, in days
+                                              !    since Jan 1 of current year
+  REAL,    INTENT(IN), OPTIONAL :: lat_in(:)  !--- Input/output latitudes vector
+  REAL,    INTENT(IN), OPTIONAL :: pcen_in(:) !--- Mid-layers pressure
+  REAL,    INTENT(IN), OPTIONAL :: ptrop_ou(:)!--- Output tropopause pressure (klon)
+!-------------------------------------------------------------------------------
+! Local variables:
+  include "clesphys.h"
+  include "YOMCST.h"
+  CHARACTER(LEN=80)  :: sub
+  CHARACTER(LEN=320) :: msg
+  INTEGER :: vID, ncerr, n_var, nlev_in, ndim, i, ibot, itop, itrp, itrp0
+  LOGICAL :: ltrop                            !--- Need to adjust tropopause
+  REAL    :: y_frac                           !--- Elapsed year fraction
+  REAL    :: alpha, beta, al                  !--- For strectching/interpolation
+  REAL    :: SigT_in, SigT_ou, SigT_mn        !--- Tropopause: in/out/mean
+  REAL    :: SigA_bot, SigA_top               !--- Strained domain bounds
+  REAL    :: Sig_in (SIZE(pint_in))           !--- Input field sigma levels
+  REAL    :: phi    (SIZE(pint_in))           !--- Stretching exponent anomaly
+  REAl    :: Pint_st(SIZE(pint_in))           !--- Stretched pressure levels
+  REAL    :: v1(nbp_lon,nbp_lat,SIZE(pint_in)-1,SIZE(nam))
+  REAL    :: v2(klon,           SIZE(pint_in)-1,SIZE(nam))
+! v1: Field read/interpol at time "julien" on the global "dynamics" horiz. grid.
+! v2: Field scattered to the partial "physics" horizontal grid.
+!     In the 2D file case, the values are the same for all longitudes.
+!     "v2(i,    k, l)" is at longitude-latitude  "xlon(i)-xlat(i)".
+! Both are for pressure interval "press_in_edg(k:k+1)]" and variable "nam(l)"
+  REAL, DIMENSION(nbp_lon, nbp_lat) ::  ps1, pt1, ot1
+  REAL, DIMENSION(klon)             ::  ps2, pt2, ot2
+!-------------------------------------------------------------------------------
+  sub="regr_pr_time_av"
+  CALL assert(SIZE(v3,1)==klon,                 TRIM(sub)//" v3 klon")
+  CALL assert(SIZE(v3,2)==nbp_lev,              TRIM(sub)//" v3 nbp_lev")
+  n_var = assert_eq(SIZE(nam), SIZE(v3,3),      TRIM(sub)//" v3 n_var")
+  CALL assert(SHAPE(pcen_ou)==[klon, nbp_lev+1],TRIM(sub)//" pcen_ou")
+  IF(PRESENT(lat_in))   CALL assert(SIZE(lat_in  )==klon,TRIM(sub)//" lat_in klon")
+  IF(PRESENT(ptrop_ou)) CALL assert(SIZE(ptrop_ou)==klon,TRIM(sub)//" ptrop_ou klon")
+  ltrop=PRESENT(lat_in).AND.PRESENT(pcen_in).AND.PRESENT(ptrop_ou)
+  nlev_in=SIZE(pint_in)-1
+
+  !$OMP MASTER
+  IF(is_mpi_root) THEN
+
+    !=== CHECK WHICH FIELDS ARE AVAILABLE IN THE INPUT FILE
+    IF(lfirst) THEN
+      lPrSurf=NF90_INQ_VARID(fID,"ps"                     ,vID)==NF90_NOERR
+      lPrTrop=NF90_INQ_VARID(fID,"tropopause_air_pressure",vID)==NF90_NOERR
+      lO3Trop=NF90_INQ_VARID(fID,"tro3_at_tropopause"     ,vID)==NF90_NOERR
+      linterp=PRESENT(time_in)
+      IF(linterp) linterp=SIZE(time_in,1)==14
+      IF(linterp) THEN
+        ALLOCATE(v1m(nbp_lon,nbp_lat,nlev_in,n_var))
+        ALLOCATE(v1p(nbp_lon,nbp_lat,nlev_in,n_var))
+        ALLOCATE(psm(nbp_lon,nbp_lat),psp(nbp_lon,nbp_lat))
+        ALLOCATE(ptm(nbp_lon,nbp_lat),ptp(nbp_lon,nbp_lat))
+        IF(lO3Trop) ALLOCATE(otm(nbp_lon,nbp_lat),otp(nbp_lon,nbp_lat))
+      END IF
+      IF(PRESENT(pcen_in)) itrp0=locate(pcen_in,prtrop) !--- Above tropopause: 50hPa
+      IF(lPrSurf) WRITE(*,*)' >> Using GROUND PRESSURE from input O3 forcing file.'
+      IF(linterp) WRITE(*,*)' >> Monthly O3 files => ONLINE TIME INTERPOLATION.'
+      WRITE(*,*)' >> o3 forcing file tropopause location uses:'
+      IF(lPrTrop) THEN;      WRITE(*,*)'    PRESSURE AT TROPOPAUSE from file.'
+      ELSE IF(lO3Trop) THEN; WRITE(*,*)'    O3 CONCENTRATION AT TROPOPAUSE from file.'
+      ELSE;                  WRITE(*,*)'    PARAMETRIZATION of O3 concentration at tropopause.'
+      END IF
+    END IF
+
+    !=== UPDATE (ALWAYS FOR DAILY FILES, EACH MONTH FOR MONTHLY FILES)
+    IF(.NOT.linterp.OR.(linterp.AND.(lfirst.OR.julien>=time_in(irec+1)))) &
+      CALL update_fields()
+
+    !=== TIME INTERPOLATION FOR MONTHLY INPUT FILES
+    IF(linterp) THEN
+      WRITE(*,'(3(a,f7.3))')'  >> Interpolating O3 at julian day ',julien,' fr&
+      &om forcing fields at times ',time_in(irec),' and ', time_in(irec+1)
+      al=(time_in(irec+1)-julien)/(time_in(irec+1)-time_in(irec))
+      v1=al*v1m+(1.-al)*v1p
+      IF(lPrSurf) ps1=al*psm+(1.-al)*psp
+      IF(lPrTrop) pt1=al*ptm+(1.-al)*ptp
+      IF(lO3Trop) ot1=al*otm+(1.-al)*otp
+    END IF
+  END IF
+  !$OMP END MASTER
+  !$OMP BARRIER
+  CALL scatter2d(v1,v2)
+  IF(PRESENT(pcen_in)) CALL bcast(itrp0)
+  !--- "ps" not in input file => assume it is equal to current LMDZ value
+  IF(lPrSurf) THEN; CALL scatter2d(ps1,ps2); ELSE; ps2=pcen_ou(:,1); END IF
+  IF(lPrTrop) CALL scatter2d(pt1,pt2)
+  IF(lO3Trop) CALL scatter2d(ot1,ot2)
+
+  !--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing
+  IF(.NOT.ltrop) THEN
+    DO i=1,klon
+      CALL regr_conserv(1,v2(i,:,:) , Pint_in , Pcen_ou(i,nbp_lev+1:1:-1),    &
+                          v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),pint_in))
+    END DO
+  ELSE
+    y_frac=(REAL(days_elapsed)+jH_cur)/year_len
+    DO i=1,klon
+      SigT_in = get_SigTrop(i,itrp)        !--- input  (file) tropopause
+      SigT_ou = ptrop_ou(i)/pcen_ou(i,1)   !--- output (lmdz) tropopause
+      SigT_mn = SQRT(SigT_in*SigT_ou)      !--- mean tropopause>strained domain
+
+      !--- DEFINE THE VERTICAL LAYER IN WHICH THE PROFILE IS STRECHED
+      beta=EXP(delta/scaleheight); Sig_in(:) = [pint_in(1:nlev_in)/ps2(i),1.]
+      SigA_bot = SigT_mn*beta ; ibot=locate(Sig_in(:),SigA_bot)
+      SigA_top = SigT_mn/beta ; itop=locate(Sig_in(:),SigA_top)+1
+
+      !--- HAT FUNCTION phi (/=0 in [SigA_bot,SigA_top] only)
+      phi(:)=0.
+      phi(itop:itrp)=(Sig_in(itop:itrp)-SigA_top)/(SigT_in-SigA_top)
+      phi(itrp:ibot)=(SigA_bot-Sig_in(itrp:ibot))/(SigA_bot-SigT_in)
+
+      !--- STRAINED INPUT logP PROFILE ; alpha: MAX. STRETCH. EXPONENT INCREMENT
+      alpha = LOG(SigT_ou/SigT_in)/LOG(SigT_in)
+      Pint_st(:) = pcen_ou(i,1) * Sig_in(:)**(1+alpha*phi(:))
+
+      !--- REGRID INPUT PROFILE ON STRAINED VERTICAL LEVELS
+      CALL regr_conserv(1,v2(i,:,:) , Pint_st,  Pcen_ou(i,nbp_lev+1:1:-1),     &
+                          v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),Pint_st))
+    END DO
+  END IF
+
+
+CONTAINS
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE update_fields()
+!
+!-------------------------------------------------------------------------------
+  IF(.NOT.linterp) THEN                 !=== DAILY FILES: NO TIME INTERPOLATION
+    WRITE(*,*)' >> Updating Ozone forcing field: read from file.'
+    irec=INT(julien)+1                  !--- irec is just the julian day number
+    CALL get_3Dfields(v1)               !--- Read ozone field(s)
+    IF(ltrop) THEN                      !--- Additional files for fields strain
+      IF(lPrSurf) CALL get_2Dfield(ps1,"ps")
+      IF(lPrTrop) CALL get_2Dfield(pt1,"tropopause_air_pressure")
+      IF(lO3Trop) CALL get_2Dfield(ot1,"tro3_at_tropopause")
+    END IF
+  ELSE                                  !=== MONTHLY FILES: GET 2 NEAREST RECS
+    WRITE(*,*)' >> Refreshing adjacent Ozone forcing fields.'
+    IF(lfirst) THEN                     !=== READ EARLIEST TIME FIELDS
+      irec=locate(time_in,julien)       !--- Need to locate surrounding times
+      CALL get_3Dfields(v1m)            !--- Read ozone field(s)
+      IF(ltrop) THEN                    !--- Additional files for fields strain
+        IF(lPrSurf) CALL get_2Dfield(psm,"ps")
+        IF(lPrTrop) CALL get_2Dfield(ptm,"tropopause_air_pressure")
+        IF(lO3Trop) CALL get_2Dfield(otm,"tro3_at_tropopause")
+      END IF
+    ELSE                                !=== SHIFT FIELDS
+      v1m=v1p                           !--- Ozone fields
+      IF(ltrop) THEN                    !--- Additional files for fields strain
+        IF(lPrSurf) psm=psp             !--- Surface pressure
+        IF(lPrTrop) ptm=ptp             !--- Tropopause pressure
+        IF(lO3Trop) otm=otp             !--- Tropopause ozone
+      END IF
+    END IF
+    irec=irec+1
+    CALL get_3Dfields(v1p)              !--- Read ozone field(s)
+    IF(ltrop) THEN                      !--- Additional files for fields strain
+      IF(lPrSurf) CALL get_2Dfield(psp,"ps")
+      IF(lPrTrop) CALL get_2Dfield(ptp,"tropopause_air_pressure")
+      IF(lO3Trop) CALL get_2Dfield(otp,"tro3_at_tropopause")
+    END IF
+    IF(lfirst) irec=irec-1
+  END IF
+  lfirst=.FALSE.
+
+END SUBROUTINE update_fields
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE get_2Dfield(v,var)
+!
+!-------------------------------------------------------------------------------
+! Purpose: Shortcut to get the "irec"th record of the surface field named "var"
+!          from the input file.
+! Remark: In case the field is zonal, it is duplicated along first dimension.
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL,             INTENT(INOUT) :: v(:,:)
+  CHARACTER(LEN=*), INTENT(IN)    :: var
+!-------------------------------------------------------------------------------
+  CALL NF95_INQ_VARID(fID, TRIM(var), vID)
+  CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=ndim)
+  IF(ndim==2) ncerr=NF90_GET_VAR(fID,vID,v(1,:), start=[  1,irec])
+  IF(ndim==3) ncerr=NF90_GET_VAR(fID,vID,v(:,:), start=[1,1,irec])
+  CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(var),ncerr,fID)
+
+  !--- Flip latitudes: ascending in input file, descending in "rlatu".
+  IF(ndim==3) THEN
+    v(1,:) = v(1,nbp_lat:1:-1)
+    v(2:,:)= SPREAD(v(1,:),DIM=1,ncopies=nbp_lon-1)  !--- Duplication
+  ELSE
+    v(:,:) = v(:,nbp_lat:1:-1)
+  END IF
+
+END SUBROUTINE get_2Dfield
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE get_3Dfields(v)
+!
+!-------------------------------------------------------------------------------
+! Purpose: Shortcut to get the "irec"th record of the 3D fields named "nam"
+!          from the input file. Fields are stacked on fourth dimension.
+! Remark: In case the field is zonal, it is duplicated along first dimension.
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(INOUT) :: v(:,:,:,:)
+!-------------------------------------------------------------------------------
+  DO i=1,SIZE(nam)
+    CALL NF95_INQ_VARID(fID, TRIM(nam(i)), vID)
+    CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=ndim)
+    IF(ndim==3) ncerr=NF90_GET_VAR(fID,vID,v(1,:,:,i), start=[  1,1,irec])
+    IF(ndim==4) ncerr=NF90_GET_VAR(fID,vID,v(:,:,:,i), start=[1,1,1,irec])
+    CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(nam(i)),ncerr,fID)
+  END DO
+
+  !--- Flip latitudes: ascending in input file, descending in "rlatu".
+  IF(ndim==3) THEN
+    v(1,:,:,:) = v(1,nbp_lat:1:-1,:,:)
+    v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nbp_lon-1)  !--- Duplication
+  ELSE
+    v(:,:,:,:) = v(:,nbp_lat:1:-1,:,:)
+  END IF
+
+END SUBROUTINE get_3Dfields
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+FUNCTION get_SigTrop(ih,it)
+!
+!-------------------------------------------------------------------------------
+! Arguments:
+  INTEGER, INTENT(IN)  :: ih
+  INTEGER, INTENT(OUT) :: it
+  REAL                 :: get_Sigtrop
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: ii                                !--- Idx of layer containing o3t
+  REAL    :: o3t                               !--- Ozone concent. at tropopause
+  REAL    :: prt                               !--- Air pressure   at tropopause
+  REAL    :: al                                !--- Interpolation coefficient
+  REAL    :: y_frac                            !--- Elapsed year fraction
+  REAL    :: coef                              !--- Coef: North/South transition
+!-------------------------------------------------------------------------------
+  !--- DETERMINE PRESSURE AT TROPOPAUSE AND DIVIDE IT BY GROUND PRESSURE
+  IF(lPrTrop) THEN                             !--- PrTrop KNOWN FROM FILE
+    get_SigTrop=pt2(ih)/ps2(ih)
+    it=locate(pint_in(:),pt2(ih))
+  ELSE                                         !--- O3 THRESHOLD
+    coef = TANH(lat_in(i)/co3(3))              !--- Latitude dependant ponderat.
+    coef = SIN (2.*RPI*(y_frac-1./6.)) * coef  !--- Time-dependant ponderation
+    o3t  = 1.E-9 * (co3(1)+co3(2)*coef)        !--- Value from parametrization
+    IF(lO3Trop) o3t=ot2(ih)                    !--- Value from file
+    !--- Starts from 50hPa and go down.
+    it=itrp0; DO WHILE(v2(ih,it+1,1)>=o3t); it=it+1; END DO
+    al=(v2(ih,it,1)-o3t)/(v2(ih,it,1)-v2(ih,it+1,1))
+    get_SigTrop = ( pcen_in(it)**(1.-al) * pcen_in(it+1)**al )/ps2(ih)
+  END IF
+
+END FUNCTION get_SigTrop
+!
+!-------------------------------------------------------------------------------
+
+
+END SUBROUTINE regr_pr_time_av
+!
+!-------------------------------------------------------------------------------
+
+END MODULE regr_pr_time_av_m
Index: LMDZ5/trunk/libf/phylmd/tropopause_m.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/tropopause_m.F90	(revision 2788)
+++ LMDZ5/trunk/libf/phylmd/tropopause_m.F90	(revision 2788)
@@ -0,0 +1,104 @@
+MODULE tropopause_m
+
+!  USE phys_local_var_mod, ONLY: ptrop => pr_tropopause
+  PRIVATE
+  PUBLIC :: dyn_tropopause
+
+CONTAINS
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE dyn_tropopause(t,ts,paprs,pplay,presnivs,rot,ptrop,tpot0,pvor0,pmin0,pmax0)
+!
+!-------------------------------------------------------------------------------
+  USE assert_m,      ONLY: assert
+  USE assert_eq_m,   ONLY: assert_eq
+  USE comvert_mod,   ONLY: preff
+  USE dimphy,        ONLY: klon, klev
+  USE geometry_mod,  ONLY: latitude_deg, longitude_deg
+  USE interpolation, ONLY: locate
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)  ::      t(:,:) !--- Cells-centers temperature
+  REAL, INTENT(IN)  ::     ts(:)   !--- Surface       temperature
+  REAL, INTENT(IN)  ::  paprs(:,:) !--- Cells-edges   pressure
+  REAL, INTENT(IN)  ::  pplay(:,:) !--- Cells-centers pressure
+  REAL, INTENT(IN)  :: presnivs(:) !--- Cells-centers nominal pressure
+  REAL, INTENT(IN)  ::    rot(:,:) !--- Cells-centers relative vorticity
+  REAL, INTENT(OUT) :: ptrop(klon) !--- Tropopause pressure
+  REAL, INTENT(IN), OPTIONAL :: tpot0, pvor0, pmin0, pmax0
+!-------------------------------------------------------------------------------
+! Local variables:
+  include "YOMCST.h"
+  CHARACTER(LEN=80)  :: sub
+  INTEGER :: it, i, nsrf, k, nh, kmin, kmax
+  REAL    :: p1, p2, tpo0, pvo0, pmn0, pmx0, al
+  REAL, DIMENSION(klon,klev)   :: t_edg, t_pot
+  REAL, DIMENSION(klon,klev-1) :: pvort
+!-------------------------------------------------------------------------------
+  sub='dyn_tropopause'
+  CALL assert(SIZE(t ,1)==klon, TRIM(sub)//" t klon")
+  CALL assert(SIZE(t ,2)==klev, TRIM(sub)//" t klev")
+  CALL assert(SIZE(ts,1)==klon, TRIM(sub)//" ts klon")
+  CALL assert(SIZE(presnivs,1)==klev, TRIM(sub)//" presnivs klev")
+  CALL assert(SHAPE(paprs)==[klon,klev+1],TRIM(sub)//" paprs shape")
+  CALL assert(SHAPE(pplay)==[klon,klev  ],TRIM(sub)//" pplay shape")
+  CALL assert(SHAPE(rot)  ==[klon,klev  ],TRIM(sub)//" rot shape")
+
+  !--- DEFAULT THRESHOLDS
+  tpo0=380.;   IF(PRESENT(tpot0)) tpo0=tpot0  !--- In kelvins
+  pvo0=2.0;    IF(PRESENT(pvor0)) pvo0=pvor0  !--- In PVU
+  pmn0= 8000.; IF(PRESENT(pmin0)) pmn0=pmin0  !--- In pascals
+  pmx0=50000.; IF(PRESENT(pmax0)) pmx0=pmax0  !--- In pascals
+  kmin=klev-locate(presnivs(klev:1:-1),pmx0)+1
+  kmax=klev-locate(presnivs(klev:1:-1),pmn0)+1
+
+  !--- POTENTIAL TEMPERATURE AT CELLS CENTERS
+  DO k= 1, klev
+    DO i = 1, klon
+      t_pot(i,k) = t(i,k)*(preff/pplay(i,k))**RKAPPA
+    END DO
+  END DO
+
+  !--- TEMPERATURE AT CELLS INTERFACES (except in top layer)
+  t_edg(:,1) = ts(:)
+  DO k= 2, klev
+    DO i = 1, klon
+      al = LOG(pplay(i,k-1)/paprs(i,k))/LOG(pplay(i,k-1)/pplay(i,k))
+      t_edg(i,k) = t(i,k) + al*(t(i,k-1)-t(i,k))
+    END DO
+  END DO
+
+  !--- ERTEL POTENTIAL VORTICITY AT CELLS CENTERS (except in top layer)
+  DO k= 1, klev-1
+    DO i = 1, klon
+      IF(paprs(i,k)==paprs(i,k+1)) THEN; pvort(i,k)=HUGE(1.); CYCLE; END IF
+      pvort(i,k) = -1.E6 * RG * 2.*ROMEGA*ABS(SIN(latitude_deg(i)*RPI/180.))   &
+                         * ( t_edg(i,k+1)*(preff/paprs(i,k+1) )**RKAPPA        &
+                           - t_edg(i,k  )*(preff/paprs(i,k  ) )**RKAPPA)       &
+                         / ( paprs(i,k+1)  -     paprs(i,k  ) )
+    END DO
+  END DO
+
+  !--- LOCATE TROPOPAUSE
+  ptrop(:)=0.
+  DO i = 1, klon
+    !--- Dynamical tropopause
+    it=kmax; DO WHILE(pvort(i,it)>pvo0.AND.it>=kmin); it=it-1; END DO
+    IF(pvort(i,it+1)/=pvort(i,it).AND.ALL([kmax,kmin-1]/=it) &
+      .AND.ALL([pvort(i,it),pvort(i,it+1)]/=HUGE(1.))) THEN
+      al = (pvort(i,it+1)-pvo0)/(pvort(i,it+1)-pvort(i,it))
+      ptrop(i) = MAX(ptrop(i),pplay(i,it)**al * pplay(i,it+1)**(1.-al))
+    END IF
+    !--- Potential temperature iso-surface
+    it = kmin-1+locate(t_pot(i,kmin:kmax),tpo0)
+    IF(t_pot(i,it+1)/=t_pot(i,it).AND.ALL([kmin-1,kmax]/=it)) THEN
+      al = (t_pot(i,it+1)-tpo0)/(t_pot(i,it+1)-t_pot(i,it))
+      ptrop(i) = MAX(ptrop(i),pplay(i,it)**al * pplay(i,it+1)**(1.-al))
+    END IF
+  END DO
+
+END SUBROUTINE dyn_tropopause
+
+END MODULE tropopause_m
