Index: LMDZ4/branches/LMDZ4-dev/libf/bibio/netcdf95.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/bibio/netcdf95.F90	(revision 1262)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/netcdf95.F90	(revision 1263)
@@ -25,4 +25,8 @@
   ! criticisms for some (not all) procedures.
 
+  ! "nf95_get_att" is more secure than "nf90_get_att" because it
+  ! checks that the "values" argument is long enough and removes the
+  ! null terminator, if any.
+
   ! This module replaces some of the official NetCDF procedures.
   ! This module also provides the procedures "handle_err" and "nf95_gw_var".
@@ -35,4 +39,5 @@
   use nf95_gw_var_m
   use nf95_put_att_m
+  use nf95_get_att_m
   use simple
   use handle_err_m
Index: LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_def_var_m.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_def_var_m.F90	(revision 1262)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_def_var_m.F90	(revision 1263)
@@ -1,4 +1,10 @@
 ! $Id$
 module nf95_def_var_m
+
+  ! The generic procedure name "nf90_def_var" applies to
+  ! "nf90_def_var_Scalar" but we cannot apply the generic procedure name
+  ! "nf95_def_var" to "nf95_def_var_scalar" because of the additional
+  ! optional argument.
+  ! "nf95_def_var_scalar" cannot be distinguished from "nf95_def_var_oneDim".
 
   implicit none
@@ -9,7 +15,34 @@
 
   private
-  public nf95_def_var
+  public nf95_def_var, nf95_def_var_scalar
 
 contains
+
+  subroutine nf95_def_var_scalar(ncid, name, xtype, varid, ncerr)
+
+    use netcdf, only: nf90_def_var
+    use handle_err_m, only: handle_err
+
+    integer,               intent( in) :: ncid
+    character (len = *),   intent( in) :: name
+    integer,               intent( in) :: xtype
+    integer,               intent(out) :: varid
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_def_var(ncid, name, xtype, varid)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_def_var_scalar " // name, ncerr_not_opt, ncid)
+    end if
+
+  end subroutine nf95_def_var_scalar
+
+  !***********************
 
   subroutine nf95_def_var_oneDim(ncid, name, xtype, dimids, varid, ncerr)
Index: LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_get_att_m.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_get_att_m.F90	(revision 1263)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_get_att_m.F90	(revision 1263)
@@ -0,0 +1,60 @@
+! $Id$
+module nf95_get_att_m
+
+  implicit none
+
+  interface nf95_get_att
+     module procedure nf95_get_att_text
+  end interface
+
+  private
+  public nf95_get_att
+
+contains
+
+  subroutine nf95_get_att_text(ncid, varid, name, values, ncerr)
+
+    use netcdf, only: nf90_get_att, nf90_inquire_attribute, nf90_noerr
+    use handle_err_m, only: handle_err
+
+    integer,                          intent( in) :: ncid, varid
+    character(len = *),               intent( in) :: name
+    character(len = *),               intent(out) :: values
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+    integer att_len
+
+    !-------------------
+
+    ! Check that the length of "values" is large enough:
+    ncerr_not_opt = nf90_inquire_attribute(ncid, varid, name, len=att_len)
+    call handle_err("nf95_get_att_text nf90_inquire_attribute " &
+         // trim(name), ncerr_not_opt, ncid, varid)
+    if (len(values) < att_len) then
+       print *, "nf95_get_att_text"
+       print *, "varid = ", varid
+       print *, "attribute name: ", name
+       print *, 'length of "values" is not large enough'
+       print *, "len(values) = ", len(values)
+       print *, "number of characters in attribute: ", att_len
+       stop 1
+    end if
+
+    values = "" ! useless in NetCDF version 3.6.2 or better
+    ncerr_not_opt = nf90_get_att(ncid, varid, name, values)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_get_att_text", ncerr_not_opt, ncid, varid)
+    end if
+
+    if (att_len >= 1 .and. ncerr_not_opt == nf90_noerr) then
+       ! Remove null terminator, if any:
+       if (iachar(values(att_len:att_len)) == 0) values(att_len:att_len) = " "
+    end if
+
+  end subroutine nf95_get_att_text
+
+end module nf95_get_att_m
Index: LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_put_var_m.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_put_var_m.F90	(revision 1262)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_put_var_m.F90	(revision 1263)
@@ -5,5 +5,6 @@
 
   interface nf95_put_var
-     module procedure nf95_put_var_1D_FourByteReal, &
+     module procedure nf95_put_var_FourByteReal, nf95_put_var_FourByteInt, &
+          nf95_put_var_1D_FourByteReal, nf95_put_var_1D_FourByteInt, &
           nf95_put_var_2D_FourByteReal, nf95_put_var_3D_FourByteReal, &
           nf95_put_var_4D_FourByteReal
@@ -19,4 +20,58 @@
 contains
 
+  subroutine nf95_put_var_FourByteReal(ncid, varid, values, start, ncerr)
+
+    use netcdf, only: nf90_put_var
+    use handle_err_m, only: handle_err
+
+    integer, intent( in) :: ncid, varid
+    real, intent( in) :: values
+    integer, dimension(:), optional, intent( in) :: start
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_put_var(ncid, varid, values, start)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_put_var_FourByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_put_var_FourByteReal
+
+  !***********************
+
+  subroutine nf95_put_var_FourByteInt(ncid, varid, values, start, ncerr)
+
+    use netcdf, only: nf90_put_var
+    use handle_err_m, only: handle_err
+
+    integer, intent( in) :: ncid, varid
+    integer, intent( in) :: values
+    integer, dimension(:), optional, intent( in) :: start
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_put_var(ncid, varid, values, start)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_put_var_FourByteInt", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_put_var_FourByteInt
+
+  !***********************
+
   subroutine nf95_put_var_1D_FourByteReal(ncid, varid, values, start, count, &
        stride, map, ncerr)
@@ -45,4 +100,33 @@
 
   end subroutine nf95_put_var_1D_FourByteReal
+
+  !***********************
+
+  subroutine nf95_put_var_1D_FourByteInt(ncid, varid, values, start, count, &
+       stride, map, ncerr)
+
+    use netcdf, only: nf90_put_var
+    use handle_err_m, only: handle_err
+
+    integer,                         intent(in) :: ncid, varid
+    integer, intent(in) :: values(:)
+    integer, dimension(:), optional, intent(in) :: start, count, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_put_var(ncid, varid, values, start, count, stride, &
+         map)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_put_var_1D_FourByteInt", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_put_var_1D_FourByteInt
 
   !***********************
Index: LMDZ4/branches/LMDZ4-dev/libf/bibio/regr1_step_av_m.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/bibio/regr1_step_av_m.F90	(revision 1262)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/regr1_step_av_m.F90	(revision 1263)
@@ -17,5 +17,6 @@
      ! The difference between the procedures is the rank of the first argument.
 
-     module procedure regr11_step_av, regr12_step_av, regr13_step_av
+     module procedure regr11_step_av, regr12_step_av, regr13_step_av, &
+          regr14_step_av
   end interface
 
@@ -203,3 +204,65 @@
   end function regr13_step_av
 
+  !********************************************
+
+  function regr14_step_av(vs, xs, xt) result(vt)
+
+    ! "vs" has rank 4.
+
+    use assert_eq_m, only: assert_eq
+    use assert_m, only: assert
+    use interpolation, only: locate
+
+    real, intent(in):: vs(:, :, :, :) ! values of steps on the source grid
+    ! (Step "is" is between "xs(is)" and "xs(is + 1)".)
+
+    real, intent(in):: xs(:)
+    ! (edges of steps on the source grid, in strictly increasing order)
+
+    real, intent(in):: xt(:)
+    ! (edges of cells of the target grid, in strictly increasing order)
+
+    real vt(size(xt) - 1, size(vs, 2), size(vs, 3), size(vs, 4))
+    ! (average values on the target grid)
+    ! (Cell "it" is between "xt(it)" and "xt(it + 1)".)
+
+    ! Variables local to the procedure:
+    integer is, it, ns, nt
+    real left_edge
+
+    !---------------------------------------------
+
+    ns = assert_eq(size(vs, 1), size(xs) - 1, "regr14_step_av ns")
+    nt = size(xt) - 1
+
+    ! Quick check on sort order:
+    call assert(xs(1) < xs(2), "regr14_step_av xs bad order")
+    call assert(xt(1) < xt(2), "regr14_step_av xt bad order")
+
+    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
+         "regr14_step_av extrapolation")
+
+    is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
+    do it = 1, nt
+       ! 1 <= is <= ns
+       ! xs(is) <= xt(it) < xs(is + 1)
+       ! Compute "vt(it, :, :, :)":
+       left_edge = xt(it)
+       vt(it, :, :, :) = 0.
+       do while (xs(is + 1) < xt(it + 1))
+          ! 1 <= is <= ns - 1
+          vt(it, :, :, :) = vt(it, :, :, :) + (xs(is + 1) - left_edge) &
+               * vs(is, :, :, :)
+          is = is + 1
+          left_edge = xs(is)
+       end do
+       ! 1 <= is <= ns
+       vt(it, :, :, :) = (vt(it, :, :, :) + (xt(it + 1) - left_edge) &
+            * vs(is, :, :, :)) / (xt(it + 1) - xt(it))
+       if (xs(is + 1) == xt(it + 1)) is = is + 1
+       ! 1 <= is <= ns .or. it == nt
+    end do
+
+  end function regr14_step_av
+
 end module regr1_step_av_m
Index: LMDZ4/branches/LMDZ4-dev/libf/bibio/regr3_lint_m.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/bibio/regr3_lint_m.F90	(revision 1262)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/regr3_lint_m.F90	(revision 1263)
@@ -11,5 +11,5 @@
      ! input array.
      ! The difference betwwen the procedures is the rank of the first argument.
-     module procedure regr33_lint
+     module procedure regr33_lint, regr34_lint
   end interface
 
@@ -57,3 +57,44 @@
   end function regr33_lint
 
+  !*********************************************************
+
+  function regr34_lint(vs, xs, xt) result(vt)
+
+    ! "vs" has rank 4.
+
+    use assert_eq_m, only: assert_eq
+    use interpolation, only: hunt
+
+    real, intent(in):: vs(:, :, :, :)
+    ! (values of the function at source points "xs")
+
+    real, intent(in):: xs(:)
+    ! (abscissas of points in source grid, in strictly monotonic order)
+
+    real, intent(in):: xt(:)
+    ! (abscissas of points in target grid)
+
+    real vt(size(vs, 1), size(vs, 2), size(xt), size(vs, 4))
+    ! (values of the function on the target grid)
+
+    ! Variables local to the procedure:
+    integer is, it, ns
+    integer is_b ! "is" bound between 1 and "ns - 1"
+
+    !--------------------------------------
+
+    ns = assert_eq(size(vs, 3), size(xs), "regr34_lint ns")
+
+    is = -1 ! go immediately to bisection on first call to "hunt"
+
+    do it = 1, size(xt)
+       call hunt(xs, xt(it), is)
+       is_b = min(max(is, 1), ns - 1)
+       vt(:, :, it, :) = ((xs(is_b+1) - xt(it)) * vs(:, :, is_b, :) &
+            + (xt(it) - xs(is_b)) * vs(:, :, is_b+1, :)) &
+            / (xs(is_b+1) - xs(is_b))
+    end do
+
+  end function regr34_lint
+
 end module regr3_lint_m
Index: LMDZ4/branches/LMDZ4-dev/libf/bibio/simple.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/bibio/simple.F90	(revision 1262)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/simple.F90	(revision 1263)
@@ -118,5 +118,5 @@
     ! unknown in the calling procedure, before the call.
     ! Here we use a better solution: a pointer argument array.
-    ! This procedure associates and defines '"dimids" if it is present.
+    ! This procedure associates and defines "dimids" if it is present.
 
     use netcdf, only: nf90_inquire_variable, nf90_max_var_dims
Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/create_etat0_limit.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/dyn3d/create_etat0_limit.F	(revision 1262)
+++ LMDZ4/branches/LMDZ4-dev/libf/dyn3d/create_etat0_limit.F	(revision 1263)
@@ -47,5 +47,6 @@
       END IF
 
-      CALL Init_Phys_lmdz(iim,jjp1,llm,1,(jjm-1)*iim+2)
+      CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
+      PRINT *,'---> klon=',klon
       call InitComgeomphy
 
Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/etat0_netcdf.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/dyn3d/etat0_netcdf.F	(revision 1262)
+++ LMDZ4/branches/LMDZ4-dev/libf/dyn3d/etat0_netcdf.F	(revision 1263)
@@ -15,4 +15,5 @@
       USE filtreg_mod
       use regr_lat_time_climoz_m, only: regr_lat_time_climoz
+      use conf_phys_m, only: conf_phys
 #endif
 !#endif of #ifdef CPP_EARTH
@@ -143,5 +144,11 @@
       REAL      :: solarlong0
       real :: seuil_inversion
-      logical  read_climoz ! read ozone climatology
+
+      integer  read_climoz ! read ozone climatology
+C     Allowed values are 0, 1 and 2
+C     0: do not read an ozone climatology
+C     1: read a single ozone climatology that will be used day and night
+C     2: read two ozone climatologies, the average day and night
+C     climatology and the daylight climatology
 
       !
@@ -179,5 +186,4 @@
      &                 iflag_coupl,iflag_clos,iflag_wake, read_climoz )
 
-
 ! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value)
       co2_ppm0 = co2_ppm
@@ -194,5 +200,5 @@
 
       CALL inifilr()
-      CALL phys_state_var_init()
+      CALL phys_state_var_init(read_climoz)
       !
       latfi(1) = ASIN(1.0)
@@ -419,5 +425,6 @@
       !
 
-      if (read_climoz) call regr_lat_time_climoz ! ozone climatology
+!     Ozone climatology:
+      if (read_climoz >= 1) call regr_lat_time_climoz(read_climoz)
 
       varname = 'tsol'
Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/create_etat0_limit.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/create_etat0_limit.F	(revision 1262)
+++ LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/create_etat0_limit.F	(revision 1263)
@@ -9,5 +9,4 @@
        USE mod_phys_lmdz_para
        USE mod_const_mpi
-       USE phys_state_var_mod      
        USE infotrac
 #ifdef CPP_IOIPSL
@@ -38,5 +37,4 @@
 #include "indicesol.h"
 #include  "control.h"
-#include "clesphys.h"
       REAL :: masque(iip1,jjp1)
 !      REAL :: pctsrf(iim*(jjm-1)+2, nbsrf)
@@ -62,7 +60,6 @@
      &                 for 1 process and 1 task')
       ENDIF
-      CALL phys_state_var_init
       call InitComgeomphy
-      
+
 #ifdef CPP_IOIPSL
       call ioconf_calendar('360d')
Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/etat0_netcdf.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/etat0_netcdf.F	(revision 1262)
+++ LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/etat0_netcdf.F	(revision 1263)
@@ -15,4 +15,5 @@
       USE filtreg_mod
       use regr_lat_time_climoz_m, only: regr_lat_time_climoz
+      use conf_phys_m, only: conf_phys
 #endif
 !#endif of #ifdef CPP_EARTH
@@ -132,4 +133,5 @@
       REAL                 :: bl95_b0, bl95_b1
       real                 :: fact_cldcon, facttemps,ratqsbas,ratqshaut
+      real                 :: tau_ratqs
       integer              :: iflag_cldcon
       integer              :: iflag_ratqs
@@ -142,5 +144,11 @@
       REAL      :: solarlong0
       real :: seuil_inversion
-      logical  read_climoz ! read ozone climatology
+
+      integer  read_climoz ! read ozone climatology
+C     Allowed values are 0, 1 and 2
+C     0: do not read an ozone climatology
+C     1: read a single ozone climatology that will be used day and night
+C     2: read two ozone climatologies, the average day and night
+C     climatology and the daylight climatology
 
       !
@@ -170,5 +178,5 @@
      &                 fact_cldcon, facttemps,ok_newmicro,iflag_radia,  &
      &                 iflag_cldcon,                                    &
-     &                 iflag_ratqs,ratqsbas,ratqshaut,                  &
+     &                 iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,        &
      &                 ok_ade, ok_aie, aerosol_couple,                  &
      &                 flag_aerosol, new_aod,                           &
@@ -187,9 +195,10 @@
       CALL inigeom()
 
-      CALL inifilr()
 ! Initialisation pour traceurs
       call infotrac_init
       ALLOCATE(q3d(iip1, jjp1, llm, nqtot))
-!     CALL phys_state_var_init()
+
+      CALL inifilr()
+      CALL phys_state_var_init(read_climoz)
       !
       latfi(1) = ASIN(1.0)
@@ -416,5 +425,6 @@
       !
 
-      if (read_climoz) call regr_lat_time_climoz ! ozone climatology
+!     Ozone climatology:
+      if (read_climoz >= 1) call regr_lat_time_climoz(read_climoz)
 
       varname = 'tsol'
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/conf_phys.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/conf_phys.F90	(revision 1262)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/conf_phys.F90	(revision 1263)
@@ -5,4 +5,9 @@
 !
 !
+module conf_phys_m
+
+   implicit none
+
+contains
 
   subroutine conf_phys(ok_journe, ok_mensuel, ok_instan, ok_hf, &
@@ -22,6 +27,4 @@
    USE surface_data
    USE carbon_cycle_mod, ONLY : carbon_cycle_tr, carbon_cycle_cpl
-
-   implicit none
 
  include "conema3.h"
@@ -152,6 +155,10 @@
   LOGICAL,SAVE      :: carbon_cycle_cpl_omp
 
-  logical, intent(out):: read_climoz ! read ozone climatology, OpenMP shared
-!
+  integer, intent(out):: read_climoz ! read ozone climatology, OpenMP shared
+  ! Allowed values are 0, 1 and 2
+  ! 0: do not read an ozone climatology
+  ! 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
 
 !$OMP MASTER 
@@ -1311,5 +1318,5 @@
   call getin('ecrit_LES', ecrit_LES_omp)
 !
-  read_climoz = .false. ! default value
+  read_climoz = 0 ! default value
   call getin('read_climoz', read_climoz)
 
@@ -1601,4 +1608,5 @@
   end subroutine conf_phys
 
+end module conf_phys_m
 !
 !#################################################################
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_output_mod.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_output_mod.F90	(revision 1262)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_output_mod.F90	(revision 1263)
@@ -390,4 +390,5 @@
   type(ctrl_out),save :: o_rhum         = ctrl_out((/ 2, 10, 10, 10, 10 /),'rhum')
   type(ctrl_out),save :: o_ozone        = ctrl_out((/ 2, 10, 10, 10, 10 /),'ozone')
+  type(ctrl_out),save :: o_ozone_light        = ctrl_out((/ 2, 10, 10, 10, 10 /),'ozone_daylight')
   type(ctrl_out),save :: o_upwd         = ctrl_out((/ 2, 10, 10, 10, 10 /),'upwd')
   type(ctrl_out),save :: o_dtphy        = ctrl_out((/ 2, 10, 10, 10, 1 /),'dtphy')
@@ -483,7 +484,7 @@
   
   SUBROUTINE phys_output_open(jjmp1,nlevSTD,clevSTD,nbteta, &
-                              ctetaSTD,dtime, ok_veget, &
-                              type_ocean, iflag_pbl,ok_mensuel,ok_journe, &
-                              ok_hf,ok_instan,ok_LES,ok_ade,ok_aie)   
+       ctetaSTD,dtime, ok_veget, &
+       type_ocean, iflag_pbl,ok_mensuel,ok_journe, &
+       ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, read_climoz)
 
   USE iophy 
@@ -506,4 +507,12 @@
   logical                               :: ok_mensuel, ok_journe, ok_hf, ok_instan
   logical                               :: ok_LES,ok_ade,ok_aie
+
+  integer, intent(in)::  read_climoz ! read ozone climatology
+  !     Allowed values are 0, 1 and 2
+  !     0: do not read an ozone climatology
+  !     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
+
   real                                  :: dtime
   integer                               :: idayref
@@ -979,4 +988,7 @@
  CALL histdef3d(iff,o_rhum%flag,o_rhum%name, "Relative humidity", "-")
  CALL histdef3d(iff,o_ozone%flag,o_ozone%name, "Ozone mole fraction", "-")
+ if (read_climoz == 2) &
+      CALL histdef3d(iff,o_ozone_light%flag,o_ozone_light%name, &
+      "Daylight ozone mole fraction", "-")
  CALL histdef3d(iff,o_dtphy%flag,o_dtphy%name, "Physics dT", "K/s")
  CALL histdef3d(iff,o_dqphy%flag,o_dqphy%name, "Physics dQ", "(kg/kg)/s")
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_output_write.h
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_output_write.h	(revision 1262)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_output_write.h	(revision 1263)
@@ -1004,5 +1004,11 @@
       IF (o_ozone%flag(iff)<=lev_files(iff)) THEN
          CALL histwrite_phy(nid_files(iff), o_ozone%name, itau_w,
-     $        wo * dobson_u * 1e3 / zmasse / rmo3 * rmd)
+     $        wo(:, :, 1) * dobson_u * 1e3 / zmasse / rmo3 * rmd)
+      ENDIF
+
+      IF (o_ozone_light%flag(iff)<=lev_files(iff) .and.
+     $     read_climoz == 2) THEN
+         CALL histwrite_phy(nid_files(iff), o_ozone_light%name, itau_w,
+     $        wo(:, :, 2) * dobson_u * 1e3 / zmasse / rmo3 * rmd)
       ENDIF
 
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_state_var_mod.F90	(revision 1262)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_state_var_mod.F90	(revision 1263)
@@ -203,6 +203,9 @@
 !$OMP THREADPRIVATE(albsol1,albsol2)
 
-      REAL, ALLOCATABLE, SAVE:: wo(:,:)
+      REAL, ALLOCATABLE, SAVE:: wo(:, :, :)
       ! column-density of ozone in a layer, in kilo-Dobsons
+      ! Third dimension has size 1 or 2.
+      ! "wo(:, :, 1)" is for the average day-night field, 
+      ! "wo(:, :, 2)" is for daylight time.
       !$OMP THREADPRIVATE(wo)
 
@@ -284,8 +287,17 @@
 
 !======================================================================
-SUBROUTINE phys_state_var_init
+SUBROUTINE phys_state_var_init(read_climoz)
 use dimphy
 use aero_mod
 IMPLICIT NONE
+
+integer, intent(in)::  read_climoz
+! read ozone climatology
+! Allowed values are 0, 1 and 2
+! 0: do not read an ozone climatology
+! 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
+
 #include "indicesol.h"
 #include "control.h"
@@ -365,5 +377,12 @@
       ALLOCATE(paire_ter(klon))
       ALLOCATE(albsol1(klon), albsol2(klon))
-      ALLOCATE(wo(klon,klev))
+
+      if (read_climoz <= 1) then
+         ALLOCATE(wo(klon,klev, 1))
+      else
+         ! read_climoz == 2
+         ALLOCATE(wo(klon,klev, 2))
+      end if
+      
       ALLOCATE(clwcon0(klon,klev),rnebcon0(klon,klev))
       ALLOCATE(heat(klon,klev), heat0(klon,klev)) 
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F	(revision 1262)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F	(revision 1263)
@@ -13,6 +13,5 @@
 
       USE ioipsl, only: histbeg, histvert, histdef, histend, histsync,
-     $     histwrite, ju2ymds, ymds2ju
-!JG     $     histwrite, ju2ymds, ymds2ju, ioget_year_len
+     $     histwrite, ju2ymds, ymds2ju, ioget_year_len
       USE comgeomphy
       USE phys_cal_mod
@@ -33,9 +32,11 @@
       USE phys_output_mod
       use open_climoz_m, only: open_climoz ! ozone climatology from a file
-      use regr_pr, only: regr_pr_av
+      use regr_pr_av_m, only: regr_pr_av
       use netcdf95, only: nf95_close
       use mod_phys_lmdz_mpi_data, only: is_mpi_root
       USE aero_mod
       use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
+      use conf_phys_m, only: conf_phys
+      use radlwsw_m, only: radlwsw
 
       IMPLICIT none
@@ -1093,16 +1094,24 @@
       integer iunit
 
-      logical, save::  read_climoz ! read ozone climatology from a file
-      integer, save:: ncid_climoz ! NetCDF file containing ozone climatology
+      integer, save::  read_climoz ! read ozone climatology
+C     Allowed values are 0, 1 and 2
+C     0: do not read an ozone climatology
+C     1: read a single ozone climatology that will be used day and night
+C     2: read two ozone climatologies, the average day and night
+C     climatology and the daylight climatology
+
+      integer, save:: ncid_climoz ! NetCDF file containing ozone climatologies
 
       real, pointer, save:: press_climoz(:)
-!     edges of pressure intervals for ozone climatology, in Pa, in strictly
+!     edges of pressure intervals for ozone climatologies, in Pa, in strictly
 !     ascending order
 
-      integer,save:: co3i = 0 ! time index in NetCDF file of current ozone field
+      integer, save:: co3i = 0
+!     time index in NetCDF file of current ozone fields
 c$OMP THREADPRIVATE(co3i) 
 
       integer ro3i
-!     required time index in NetCDF file for the ozone field, between 1 and 360
+!     required time index in NetCDF file for the ozone fields, between 1
+!     and 360
 
 #include "YOMCST.h"
@@ -1165,5 +1174,19 @@
       print*, 'Allocation des variables locales et sauvegardees'
       call phys_local_var_init
-      call phys_state_var_init
+c     appel a la lecture du run.def physique
+      call conf_phys(ok_journe, ok_mensuel,
+     .     ok_instan, ok_hf,
+     .     ok_LES,
+     .     solarlong0,seuil_inversion,
+     .     fact_cldcon, facttemps,ok_newmicro,iflag_radia,
+     .     iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,
+     .     ok_ade, ok_aie, aerosol_couple, 
+     .     flag_aerosol, new_aod,
+     .     bl95_b0, bl95_b1,
+     .     iflag_thermals,nsplit_thermals,tau_thermals,
+     .     iflag_thermals_ed,iflag_thermals_optflux,
+c     nv flags pour la convection et les poches froides
+     .     iflag_coupl,iflag_clos,iflag_wake, read_climoz)
+      call phys_state_var_init(read_climoz)
       print*, '================================================='
 
@@ -1181,5 +1204,5 @@
 	first=.false.
 
-      endif  ! fisrt
+      endif  ! first
 
        modname = 'physiq'
@@ -1209,5 +1232,4 @@
 c
        IF (debut) THEN
-C
 !rv
 cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation 
@@ -1220,6 +1242,4 @@
          rain_con(:)=0.
          snow_con(:)=0.
-         bl95_b0=0.
-         bl95_b1=0.
          topswai(:)=0.
          topswad(:)=0.
@@ -1247,20 +1267,4 @@
          IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
 c
-c appel a la lecture du run.def physique
-c
-         call conf_phys(ok_journe, ok_mensuel,
-     .                  ok_instan, ok_hf,
-     .                  ok_LES,
-     .                  solarlong0,seuil_inversion,
-     .                  fact_cldcon, facttemps,ok_newmicro,iflag_radia,
-     .         iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,
-     .                  ok_ade, ok_aie, aerosol_couple, 
-     .                  flag_aerosol, new_aod,
-     .                  bl95_b0, bl95_b1,
-     .                  iflag_thermals,nsplit_thermals,tau_thermals,
-     .                  iflag_thermals_ed,iflag_thermals_optflux,
-cnv flags pour la convection et les poches froides
-     .                   iflag_coupl,iflag_clos,iflag_wake, read_climoz)
-
       print*,'iflag_coupl,iflag_clos,iflag_wake',
      .   iflag_coupl,iflag_clos,iflag_wake
@@ -1453,5 +1457,5 @@
      &                        ctetaSTD,dtime,ok_veget,
      &                        type_ocean,iflag_pbl,ok_mensuel,ok_journe,
-     &                        ok_hf,ok_instan,ok_LES,ok_ade,ok_aie)
+     &              ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, read_climoz)
 c$OMP END MASTER
 c$OMP BARRIER
@@ -1544,5 +1548,5 @@
 
 C$omp single
-      if (read_climoz) then
+      if (read_climoz >= 1) then
          call open_climoz(ncid_climoz, press_climoz)
       END IF
@@ -1711,13 +1715,9 @@
 c Prescrire l'ozone et calculer l'albedo sur l'ocean.
 c
-      if (read_climoz) then
+      if (read_climoz >= 1) then
 C        Ozone from a file
 !        Update required ozone index:
-!JG : ioget_year_len n'existe pas dans les versions officiels d'ioipsl
-!JG         ro3i = int((days_elapsed + jh_cur - jh_1jan)
-!JG     $        / ioget_year_len(year_cur) * 360.) + 1
-         CALL abort_gcm(modname,
-     $     'JG : read_climoz temporary desactivated',1)
-
+         ro3i = int((days_elapsed + jh_cur - jh_1jan)
+     $        / ioget_year_len(year_cur) * 360.) + 1
          if (ro3i == 361) ro3i = 360
 C        (This should never occur, except perhaps because of roundup
@@ -1725,9 +1725,18 @@
          if (ro3i /= co3i) then
 C           Update ozone field:
-            call regr_pr_av(ncid_climoz, "tro3", julien=ro3i, 
-     &           press_in_edg=press_climoz, paprs=paprs, v3=wo)
+            if (read_climoz == 1) then
+               call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i,
+     $              press_in_edg=press_climoz, paprs=paprs, v3=wo)
+            else
+C              read_climoz == 2
+               call regr_pr_av(ncid_climoz,
+     $              (/"tro3         ", "tro3_daylight"/),
+     $              julien=ro3i, press_in_edg=press_climoz, paprs=paprs,
+     $              v3=wo)
+            end if
 !           Convert from mole fraction of ozone to column density of ozone in a
 !           cell, in kDU:
-            wo = wo * rmo3 / rmd * zmasse / dobson_u / 1e3
+            forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l)
+     $           * rmo3 / rmd * zmasse / dobson_u / 1e3
 C           (By regridding ozone values for LMDZ only once every 360th of
 C           year, we have already neglected the variation of pressure in one
@@ -1738,5 +1747,5 @@
       elseif (MOD(itap-1,lmt_pas) == 0) THEN
 C        Once per day, update ozone from Royer:
-         wo = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
+         wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
       ENDIF
 c
@@ -2836,5 +2845,5 @@
      $                          u,
      $                          v,
-     $                          wo,
+     $                          wo(:, :, 1),
      $                          q_seri,
      $                          zxtsol,
@@ -2919,5 +2928,5 @@
      e        (kdlon,kflev,dist, rmu0, fract, solaire,
      e        paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
-     e        wo,
+     e        wo(:, :, 1),
      e        cldfra, cldemi, cldtau,
      s        heat,heat0,cool,cool0,radsol,albpla,
@@ -3529,5 +3538,5 @@
 !         close(97)
 C$OMP MASTER
-         if (read_climoz) then
+         if (read_climoz >= 1) then
             if (is_mpi_root) then
                call nf95_close(ncid_climoz)
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/radlwsw.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/radlwsw.F90	(revision 1262)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/radlwsw.F90	(revision 1263)
@@ -1,2 +1,8 @@
+module radlwsw_m
+
+  IMPLICIT NONE
+
+contains
+
 SUBROUTINE radlwsw( &
    dist, rmu0, fract, &
@@ -24,6 +30,6 @@
 
   USE DIMPHY
-
-  IMPLICIT NONE
+  use assert_m, only: assert
+
   !======================================================================
   ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
@@ -103,6 +109,8 @@
   REAL,    INTENT(in)  :: t(KLON,KLEV), q(KLON,KLEV)
 
-  REAL, INTENT(in)::wo(KLON,KLEV)
+  REAL, INTENT(in):: wo(:, :, :) ! dimension(KLON,KLEV, 1 or 2)
   ! column-density of ozone in a layer, in kilo-Dobsons
+  ! "wo(:, :, 1)" is for the average day-night field, 
+  ! "wo(:, :, 2)" is for daylight time.
 
   LOGICAL, INTENT(in)  :: ok_ade, ok_aie                                 ! switches whether to use aerosol direct (indirect) effects or not
@@ -157,5 +165,9 @@
   REAL(KIND=8) PTAVE(kdlon,kflev)
   REAL(KIND=8) PWV(kdlon,kflev), PQS(kdlon,kflev)
-  real(kind=8) POZON(kdlon,kflev) ! mass fraction of ozone
+
+  real(kind=8) POZON(kdlon, kflev, size(wo, 3)) ! mass fraction of ozone
+  ! "POZON(:, :, 1)" is for the average day-night field, 
+  ! "POZON(:, :, 2)" is for daylight time.
+
   REAL(KIND=8) PAER(kdlon,kflev,5)
   REAL(KIND=8) PCLDLD(kdlon,kflev)
@@ -187,4 +199,5 @@
   real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
 
+  call assert(size(wo, 1) == klon, size(wo, 2) == klev, "radlwsw wo")
   ! initialisation
   tauaero(:,:,:,:)=0.
@@ -245,5 +258,5 @@
         PWV(i,k) = MAX (q(iof+i,k), 1.0e-12)
         PQS(i,k) = PWV(i,k)
-        POZON(i,k) = wo(iof+i, k) * RG * dobson_u * 1e3 &
+        POZON(i,k, :) = wo(iof+i, k, :) * RG * dobson_u * 1e3 &
              / (paprs(iof+i, k) - paprs(iof+i, k+1))
         PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
@@ -298,8 +311,9 @@
     IF (iflag_rrtm == 0) THEN
        ! Old radiation scheme, used for AR4 runs
+       ! average day-night ozone for longwave
        CALL LW_LMDAR4(&
             PPMB, PDP,&
             PPSOL,PDT0,PEMIS,&
-            PTL, PTAVE, PWV, POZON, PAER,&
+            PTL, PTAVE, PWV, POZON(:, :, 1), PAER,&
             PCLDLD,PCLDLU,&
             PVIEW,&
@@ -309,5 +323,5 @@
             ZFLUP, ZFLDN, ZFLUP0,ZFLDN0)
 
-
+       ! daylight ozone, if we have it, for short wave
        IF (.NOT. new_aod) THEN 
           ! use old version
@@ -315,5 +329,5 @@
                PPMB, PDP, &
                PPSOL, PALBD, PALBP,&
-               PTAVE, PWV, PQS, POZON, PAER,&
+               PTAVE, PWV, PQS, POZON(:, :, size(wo, 3)), PAER,&
                PCLDSW, PTAU, POMEGA, PCG,&
                zheat, zheat0,&
@@ -330,5 +344,5 @@
                PPMB, PDP,&
                PPSOL, PALBD, PALBP,&
-               PTAVE, PWV, PQS, POZON, PAER,&
+               PTAVE, PWV, PQS, POZON(:, :, size(wo, 3)), PAER,&
                PCLDSW, PTAU, POMEGA, PCG,&
                zheat, zheat0,&
@@ -438,6 +452,5 @@
  ENDDO ! j = 1, nb_gr
 
-
 END SUBROUTINE radlwsw
 
-
+end module radlwsw_m
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_lat_time_climoz_m.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_lat_time_climoz_m.F90	(revision 1262)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_lat_time_climoz_m.F90	(revision 1263)
@@ -11,5 +11,5 @@
 contains
 
-  subroutine regr_lat_time_climoz
+  subroutine regr_lat_time_climoz(read_climoz)
 
     ! "regr_lat_time_climoz" stands for "regrid latitude time
@@ -23,7 +23,5 @@
 
     ! If the input field has missing values, they must be signaled by
-    ! the "missing_value" attribute. At each latitude and each date, the
-    ! missing values are replaced by the lowest valid value above missing
-    ! values.
+    ! the "missing_value" attribute.
 
     ! We assume that the input field is a step function of latitude
@@ -34,37 +32,51 @@
     ! The values of "rlatu" are taken to be the centers of intervals.
 
-    ! Regridding in time is by linear interpolation.
-    ! Monthly values are processed to get daily values, on the basis
-    ! of a 360-day calendar.
+    ! 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
-    ! We use the first December value to interpolate values between January
-    ! 1st and mid-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 than we assume
+    ! 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.
 
-    ! We assume that in the input file:
-    ! -- latitude is in degrees;
-    ! -- pressure is in hPa (even though we do not use pressure values
-    ! here, we write the unit of pressure in the NetCDF header, and we
-    ! will use the assumption later, when we regrid in pressure).
-    ! -- latitude and pressure are strictly monotonic (as all NetCDF
-    ! coordinate variables should be);
-    ! -- time increases (even though we do not use values of the input
-    ! time coordinate);
-    ! -- missing values are vertically contiguous, at the bottom of
-    ! the vertical domain;
-    ! -- there is no missing value at the top level of the atmosphere.
-
     use regr1_step_av_m, only: regr1_step_av
     use regr3_lint_m, only: regr3_lint
-    use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, handle_err, &
-         nf95_put_var, nf95_gw_var
-    use netcdf, only: nf90_nowrite, nf90_get_att, nf90_noerr
+    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
+
+    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:
@@ -79,7 +91,7 @@
     ! (for "pi")
 
-    integer ncid_in, ncid_out ! NetCDF IDs for input and output files
     integer n_plev ! number of pressure levels in the input data
-    integer n_lat! number of latitudes in the input data
+    integer n_lat ! number of latitudes in the input data
+    integer n_month ! number of months in the input data
 
     real, pointer:: latitude(:)
@@ -91,20 +103,26 @@
 
     real, pointer:: plev(:)
-    ! pressure levels of input data, sorted in strictly ascending order
+    ! 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, pointer:: o3_in(:, :, :) ! (n_lat, n_plev, 12 or 0:13)
-    ! (ozone climatology from the input file)
-    ! ("o3_in(j, k, l)" is at latitude "latitude(j)" and pressure
-    ! level "plev(k)". "l" is between 1 and 12 or between 0 and 13)
+    real, 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(:, :, :) ! (jjm + 1, n_plev, 0:13)
-    ! (mean of "o3_in" over a latitude interval of LMDZ)
-    ! (First dimension is latitude interval.
-    ! The latitude interval for "o3_regr_lat(j,:, :)" contains "rlatu(j)".
+    real, allocatable:: o3_regr_lat(:, :, :, :)
+    ! (jjm + 1, 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 "jjm" then the interval is:
     ! [rlatv(j), rlatv(j-1)]
@@ -114,17 +132,28 @@
     ! [- pi / 2, rlatv(jjm)]
     ! respectively.
-    ! "o3_regr_lat(:, k, :)" is for pressure level "plev(k)".
-    ! Last dimension is month number.)
-
-    real, allocatable:: o3_out(:, :, :) ! (jjm + 1, n_plev, 360)
-    ! (regridded ozone climatology)
-    ! ("o3_out(j, k, l)" is at latitude "rlatu(j)", pressure
+    ! "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(:, :, :, :)
+    ! (jjm + 1, 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.)
-
-    integer j, k, l
+    ! 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 varid_in, varid_out, varid_plev, varid_time, varid, ncerr
+    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)/)
@@ -141,4 +170,5 @@
 
     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)
@@ -172,4 +202,13 @@
     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:
@@ -183,42 +222,81 @@
     deallocate(plev) ! pointer
 
-    call nf95_inq_varid(ncid_in, "tro3", varid_in)
-    call nf95_gw_var(ncid_in, varid_in, o3_in)
-    if (desc_lat) o3_in = o3_in(n_lat:1:-1, :, :)
-    if (desc_plev) o3_in = o3_in(:, n_plev:1:-1, :)
-
-    ncerr = nf90_get_att(ncid_in, varid_in, "missing_value", missing_value)
-    if (ncerr == nf90_noerr) then
-       do l = 1, size(o3_in, 3)
-          do j = 1, n_lat
-             ! Find missing values, starting from top of atmosphere
-             ! and going down.
-             ! We assume that the highest level contains no missing value.
-             k = 2
-             do
-                if (o3_in(j, k, l) == missing_value .or. k == n_plev) exit
-                k = k + 1
+    ! Get the  number of months:
+    call nf95_inq_dimid(ncid_in, "time", dimid)
+    call nf95_inquire_dimension(ncid_in, dimid, len=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
-             ! Replace missing values with the valid value at the
-             ! lowest level above missing values:
-             if (o3_in(j, k, l) == missing_value) &
-                  o3_in(j, k:n_plev, l) = o3_in(j, k-1, l)
+             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
-       end do
-    else
-       print *, "regr_lat_time_climoz: no missing value attribute"
-    end if
-    
+       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(jjm + 1, n_plev, 0:13))
-    allocate(o3_out(jjm + 1, n_plev, 360))
+    allocate(o3_regr_lat(jjm + 1, n_plev, 0:13, read_climoz))
+    allocate(o3_out(jjm + 1, 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 (size(o3_in, 3) == 12) then
-       print *, "Found 12 months in ozone climatology, assuming periodicity..."
-       o3_regr_lat(jjm+1:1:-1, :, 1:12) = regr1_step_av(o3_in, &
+    if (n_month == 12) then
+       print *, &
+            "Found 12 months in ozone climatologies, assuming periodicity..."
+       o3_regr_lat(jjm+1:1:-1, :, 1:12, :) = regr1_step_av(o3_in, &
             xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(jjm:1:-1), pi / 2/)))
        ! (invert order of indices in "o3_regr_lat" because "rlatu" is
@@ -227,15 +305,13 @@
        ! 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)
+       o3_regr_lat(:, :, 0, :) = o3_regr_lat(:, :, 12, :)
+       o3_regr_lat(:, :, 13, :) = o3_regr_lat(:, :, 1, :)
     else
-       print *, "Using 14 months in ozone climatology..."
-       o3_regr_lat(jjm+1:1:-1, :, :) = regr1_step_av(o3_in, &
+       print *, "Using 14 months in ozone climatologies..."
+       o3_regr_lat(jjm+1:1:-1, :, :, :) = regr1_step_av(o3_in, &
             xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(jjm:1:-1), pi / 2/)))
        ! (invert order of indices in "o3_regr_lat" because "rlatu" is
        ! in descending order)
     end if
-    
-    deallocate(o3_in) ! pointer
 
     ! Regrid in time by linear interpolation:
@@ -243,6 +319,8 @@
 
     ! Write to file:
-    call nf95_put_var(ncid_out, varid_out, o3_out(jjm+1:1:-1, :, :))
-    ! (The order of "rlatu" is inverted in the output file)
+    do m = 1, read_climoz
+       call nf95_put_var(ncid_out, varid_out(m), o3_out(jjm+1:1:-1, :, :, m))
+       ! (The order of "rlatu" is inverted in the output file)
+    end do
 
     call nf95_close(ncid_out)
@@ -263,5 +341,9 @@
 
     integer, intent(in):: ncid_in, n_plev
-    integer, intent(out):: ncid_out, varid_out, varid_plev, varid_time
+    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:
@@ -307,11 +389,19 @@
     call nf95_put_att(ncid_out, varid_rlatu, "standard_name", "latitude")
 
-    ! Define the primary variable:
+    ! Define the primary variables:
 
     call nf95_def_var(ncid_out, "tro3", nf90_float, &
-         (/dimid_rlatu, dimid_plev, dimid_time/), varid_out)
-    call nf95_put_att(ncid_out, varid_out, "long_name", "ozone mole fraction")
-    call nf95_put_att(ncid_out, varid_out, "standard_name", &
+         (/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:
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_pr.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_pr.F90	(revision 1262)
+++ 	(revision )
@@ -1,193 +1,0 @@
-! $Id$
-module regr_pr
-
-  ! Author: Lionel GUEZ
-
-  ! In both procedures of this module:
-  ! -- the root process reads a 2D latitude-pressure field from a
-  !    NetCDF file, at a given day.
-  ! -- the field is packed to the LMDZ horizontal "physics"
-  !    grid and scattered to all threads of all processes;
-  ! -- in all the threads of all the processes, the field is regridded in
-  !    pressure to the LMDZ vertical grid.
-  ! We assume that, in the input file, the field has 3 dimensions:
-  ! latitude, pressure, julian day.
-  ! We assume that latitudes are in ascending order in the input file.
-
-
-  implicit none
-
-contains
-
-  subroutine regr_pr_av(ncid, name, julien, press_in_edg, paprs, v3)
-
-    ! "regr_pr_av" stands for "regrid pressure averaging".
-    ! The target vertical LMDZ grid is the grid of layer boundaries.
-    ! Regridding in pressure is done by averaging a step function of pressure.
-
-    use dimphy, only: klon
-    use netcdf95, only: nf95_inq_varid, handle_err
-    use netcdf, only: nf90_get_var
-    use assert_m, only: assert
-    use regr1_step_av_m, only: regr1_step_av
-    use mod_phys_lmdz_mpi_data, only: is_mpi_root
-
-    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 variable
-    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)
-    ! (regridded field on the partial "physics" grid)
-    ! ("v3(i, k)" is at longitude "xlon(i)", latitude
-    ! "xlat(i)", in pressure interval "[paprs(i, l+1), paprs(i, l)]".)
-
-    ! Variables local to the procedure:
-
-    include "dimensions.h"
-    integer varid, ncerr ! for NetCDF
-
-    real  v1(iim, jjm + 1, size(press_in_edg) - 1)
-    ! (input field at day "julien", on the global "dynamics" horizontal grid)
-    ! (First dimension is for longitude.
-    ! The value is the same for all longitudes.
-    ! "v1(:, j, k)" is at latitude "rlatu(j)" and for
-    ! pressure interval "[press_in_edg(k), press_in_edg(k+1)]".)
-
-    real v2(klon, size(press_in_edg) - 1)
-    ! (field scattered to the partial "physics" horizontal grid)
-    ! ("v2(i, k)" is at longitude "xlon(i)", latitude "xlat(i)" and
-    ! for pressure interval "[press_in_edg(k), press_in_edg(k+1)]".)
-
-    integer i
-
-    !--------------------------------------------
-
-    call assert(shape(v3) == (/klon, llm/), "regr_pr_av v3")
-    call assert(shape(paprs) == (/klon, llm+1/), "regr_pr_av paprs")
-
-    !$omp master
-    if (is_mpi_root) then
-       call nf95_inq_varid(ncid, name, varid)
-
-       ! Get data at the right day from the input file:
-       ncerr = nf90_get_var(ncid, varid, v1(1, :, :), start=(/1, 1, julien/))
-       call handle_err("regr_pr_av nf90_get_var " // name, ncerr, ncid)
-       ! 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, jjm+1:1:-1, :)
-
-       ! Duplicate on all longitudes:
-       v1(2:, :, :) = spread(v1(1, :, :), dim=1, ncopies=iim-1)
-    end if
-    !$omp end master
-
-    call scatter2d(v1, v2)
-
-    ! Regrid in pressure at each horizontal position:
-    do i = 1, klon
-       v3(i, llm:1:-1) = regr1_step_av(v2(i, :), press_in_edg, &
-            paprs(i, llm+1:1:-1))
-       ! (invert order of indices because "paprs" is in descending order)
-    end do
-
-  end subroutine regr_pr_av
-
-  !***************************************************************
-
-  subroutine regr_pr_int(ncid, name, julien, plev, pplay, top_value, v3)
-
-    ! "regr_pr_int" stands for "regrid pressure interpolation".
-    ! The target vertical LMDZ grid is the grid of mid-layers.
-    ! Regridding is by linear interpolation.
-
-    use dimphy, only: klon
-    use netcdf95, only: nf95_inq_varid, handle_err
-    use netcdf, only: nf90_get_var
-    use assert_m, only: assert
-    use regr1_lint_m, only: regr1_lint
-    use mod_phys_lmdz_mpi_data, only: is_mpi_root
-
-    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 variable
-    integer, intent(in):: julien ! jour julien, 1 <= julien <= 360
-
-    real, intent(in):: plev(:)
-    ! (pressure level of input data, in Pa, in strictly ascending order)
-
-    real, intent(in):: pplay(:, :) ! (klon, llm)
-    ! (pression pour le mileu de chaque couche, en Pa)
-
-    real, intent(in):: top_value
-    ! (extra value of field at 0 pressure)
-
-    real, intent(out):: v3(:, :) ! (klon, llm)
-    ! (regridded field on the partial "physics" grid)
-    ! ("v3(i, k)" is at longitude "xlon(i)", latitude
-    ! "xlat(i)", middle of layer "k".)
-
-    ! Variables local to the procedure:
-
-    include "dimensions.h"
-    integer varid, ncerr ! for NetCDF
-
-    real  v1(iim, jjm + 1, 0:size(plev))
-    ! (input field at day "julien", on the global "dynamics" horizontal grid)
-    ! (First dimension is for longitude.
-    ! The value is the same for all longitudes.
-    ! "v1(:, j, k >=1)" is at latitude "rlatu(j)" and pressure "plev(k)".)
-
-    real v2(klon, 0:size(plev))
-    ! (field scattered to the partial "physics" horizontal grid)
-    ! "v2(i, k >= 1)" is at longitude "xlon(i)", latitude "xlat(i)"
-    ! and pressure "plev(k)".)
-
-    integer i
-
-    !--------------------------------------------
-
-    call assert(shape(v3) == (/klon, llm/), "regr_pr_int v3")
-    call assert(shape(pplay) == (/klon, llm/), "regr_pr_int pplay")
-
-    !$omp master
-    if (is_mpi_root) then
-       call nf95_inq_varid(ncid, name, varid)
-
-       ! Get data at the right day from the input file:
-       ncerr = nf90_get_var(ncid, varid, v1(1, :, 1:), start=(/1, 1, julien/))
-       call handle_err("regr_pr_int nf90_get_var " // name, ncerr, ncid)
-       ! Latitudes are in ascending order in the input file while
-       ! "rlatu" is in descending order so we need to invert order:
-       v1(1, :, 1:) = v1(1, jjm+1:1:-1, 1:)
-
-       ! Complete "v1" with the value at 0 pressure:
-       v1(1, :, 0) = top_value
-
-       ! Duplicate on all longitudes:
-       v1(2:, :, :) = spread(v1(1, :, :), dim=1, ncopies=iim-1)
-    end if
-    !$omp end master
-
-    call scatter2d(v1, v2)
-
-    ! Regrid in pressure at each horizontal position:
-    do i = 1, klon
-       v3(i, llm:1:-1) = regr1_lint(v2(i, :), (/0., plev/), pplay(i, llm:1:-1))
-       ! (invert order of indices because "pplay" is in descending order)
-    end do
-
-  end subroutine regr_pr_int
-
-end module regr_pr
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_pr_av_m.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_pr_av_m.F90	(revision 1263)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_pr_av_m.F90	(revision 1263)
@@ -0,0 +1,121 @@
+! $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, excepth that latitudes are in ascending order in the input
+    ! file.
+    ! We assume that the inputs fields have the same pressure coordinate.
+
+    ! The target vertical LMDZ grid is the grid of layer boundaries.
+    ! Regridding in pressure is done by averaging a step function of pressure.
+
+    ! 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_step_av_m, only: regr1_step_av
+    use mod_phys_lmdz_mpi_data, only: is_mpi_root
+
+    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:
+
+    include "dimensions.h"
+    integer varid, ncerr ! for NetCDF
+
+    real  v1(iim, jjm + 1, 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) == llm, "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, llm+1/), "regr_pr_av paprs")
+
+    !$omp master
+    if (is_mpi_root) then
+       do i = 1, n_var
+          call nf95_inq_varid(ncid, 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 " // 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, jjm+1:1:-1, :, :)
+
+       ! Duplicate on all longitudes:
+       v1(2:, :, :, :) = spread(v1(1, :, :, :), dim=1, ncopies=iim-1)
+    end if
+    !$omp end master
+
+    call scatter2d(v1, v2)
+
+    ! Regrid in pressure at each horizontal position:
+    do i = 1, klon
+       v3(i, llm:1:-1, :) = regr1_step_av(v2(i, :, :), press_in_edg, &
+            paprs(i, llm+1:1:-1))
+       ! (invert order of indices because "paprs" is in descending order)
+    end do
+
+  end subroutine regr_pr_av
+
+end module regr_pr_av_m
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_pr_int_m.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_pr_int_m.F90	(revision 1263)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_pr_int_m.F90	(revision 1263)
@@ -0,0 +1,106 @@
+! $Id$
+module regr_pr_int_m
+
+  ! Author: Lionel GUEZ
+
+  implicit none
+
+contains
+
+  subroutine regr_pr_int(ncid, name, julien, plev, pplay, top_value, v3)
+
+    ! "regr_pr_int" stands for "regrid pressure interpolation".
+    ! In this procedure:
+    ! -- the root process reads a 2D latitude-pressure field from a
+    !    NetCDF file, at a given day.
+    ! -- the field is packed to the LMDZ horizontal "physics"
+    !    grid and scattered to all threads of all processes;
+    ! -- in all the threads of all the processes, the field is regridded in
+    !    pressure to the LMDZ vertical grid.
+    ! We assume that, in the input file, the field has 3 dimensions:
+    ! latitude, pressure, julian day.
+    ! We assume that latitudes are in ascending order in the input file.
+    ! The target vertical LMDZ grid is the grid of mid-layers.
+    ! Regridding is by linear interpolation.
+
+    use dimphy, only: klon
+    use netcdf95, only: nf95_inq_varid, handle_err
+    use netcdf, only: nf90_get_var
+    use assert_m, only: assert
+    use regr1_lint_m, only: regr1_lint
+    use mod_phys_lmdz_mpi_data, only: is_mpi_root
+
+    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 variable
+    integer, intent(in):: julien ! jour julien, 1 <= julien <= 360
+
+    real, intent(in):: plev(:)
+    ! (pressure level of input data, in Pa, in strictly ascending order)
+
+    real, intent(in):: pplay(:, :) ! (klon, llm)
+    ! (pression pour le mileu de chaque couche, en Pa)
+
+    real, intent(in):: top_value
+    ! (extra value of field at 0 pressure)
+
+    real, intent(out):: v3(:, :) ! (klon, llm)
+    ! (regridded field on the partial "physics" grid)
+    ! ("v3(i, k)" is at longitude "xlon(i)", latitude
+    ! "xlat(i)", middle of layer "k".)
+
+    ! Variables local to the procedure:
+
+    include "dimensions.h"
+    integer varid, ncerr ! for NetCDF
+
+    real  v1(iim, jjm + 1, 0:size(plev))
+    ! (input field at day "julien", on the global "dynamics" horizontal grid)
+    ! (First dimension is for longitude.
+    ! The value is the same for all longitudes.
+    ! "v1(:, j, k >=1)" is at latitude "rlatu(j)" and pressure "plev(k)".)
+
+    real v2(klon, 0:size(plev))
+    ! (field scattered to the partial "physics" horizontal grid)
+    ! "v2(i, k >= 1)" is at longitude "xlon(i)", latitude "xlat(i)"
+    ! and pressure "plev(k)".)
+
+    integer i
+
+    !--------------------------------------------
+
+    call assert(shape(v3) == (/klon, llm/), "regr_pr_int v3")
+    call assert(shape(pplay) == (/klon, llm/), "regr_pr_int pplay")
+
+    !$omp master
+    if (is_mpi_root) then
+       call nf95_inq_varid(ncid, name, varid)
+
+       ! Get data at the right day from the input file:
+       ncerr = nf90_get_var(ncid, varid, v1(1, :, 1:), start=(/1, 1, julien/))
+       call handle_err("regr_pr_int nf90_get_var " // name, ncerr, ncid)
+       ! Latitudes are in ascending order in the input file while
+       ! "rlatu" is in descending order so we need to invert order:
+       v1(1, :, 1:) = v1(1, jjm+1:1:-1, 1:)
+
+       ! Complete "v1" with the value at 0 pressure:
+       v1(1, :, 0) = top_value
+
+       ! Duplicate on all longitudes:
+       v1(2:, :, :) = spread(v1(1, :, :), dim=1, ncopies=iim-1)
+    end if
+    !$omp end master
+
+    call scatter2d(v1, v2)
+
+    ! Regrid in pressure at each horizontal position:
+    do i = 1, klon
+       v3(i, llm:1:-1) = regr1_lint(v2(i, :), (/0., plev/), pplay(i, llm:1:-1))
+       ! (invert order of indices because "pplay" is in descending order)
+    end do
+
+  end subroutine regr_pr_int
+
+end module regr_pr_int_m
