Index: LMDZ4/branches/LMDZ4-dev/libf/bibio/assert_eq_m.f90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/bibio/assert_eq_m.f90	(revision 1154)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/assert_eq_m.f90	(revision 1154)
@@ -0,0 +1,70 @@
+! $Id$
+MODULE assert_eq_m
+
+  implicit none
+
+  INTERFACE assert_eq
+     MODULE PROCEDURE assert_eq2,assert_eq3,assert_eq4,assert_eqn
+  END INTERFACE
+
+  private assert_eq2,assert_eq3,assert_eq4,assert_eqn
+
+CONTAINS
+
+  FUNCTION assert_eq2(n1,n2,string)
+    CHARACTER(LEN=*), INTENT(IN) :: string
+    INTEGER, INTENT(IN) :: n1,n2
+    INTEGER  assert_eq2
+    if (n1 == n2) then
+       assert_eq2=n1
+    else
+       write (*,*) 'nrerror: an assert_eq failed with this tag: ', &
+            string
+       print *, 'program terminated by assert_eq2'
+       stop 1
+    end if
+  END FUNCTION assert_eq2
+  !BL
+  FUNCTION assert_eq3(n1,n2,n3,string)
+    CHARACTER(LEN=*), INTENT(IN) :: string
+    INTEGER, INTENT(IN) :: n1,n2,n3
+    INTEGER  assert_eq3
+    if (n1 == n2 .and. n2 == n3) then
+       assert_eq3=n1
+    else
+       write (*,*) 'nrerror: an assert_eq failed with this tag: ', &
+            string
+       print *, 'program terminated by assert_eq3'
+       stop 1
+    end if
+  END FUNCTION assert_eq3
+  !BL
+  FUNCTION assert_eq4(n1,n2,n3,n4,string)
+    CHARACTER(LEN=*), INTENT(IN) :: string
+    INTEGER, INTENT(IN) :: n1,n2,n3,n4
+    INTEGER  assert_eq4
+    if (n1 == n2 .and. n2 == n3 .and. n3 == n4) then
+       assert_eq4=n1
+    else
+       write (*,*) 'nrerror: an assert_eq failed with this tag: ', &
+            string
+       print *, 'program terminated by assert_eq4'
+       stop 1
+    end if
+  END FUNCTION assert_eq4
+  !BL
+  FUNCTION assert_eqn(nn,string)
+    CHARACTER(LEN=*), INTENT(IN) :: string
+    INTEGER, DIMENSION(:), INTENT(IN) :: nn
+    INTEGER  assert_eqn
+    if (all(nn(2:) == nn(1))) then
+       assert_eqn=nn(1)
+    else
+       write (*,*) 'nrerror: an assert_eq failed with this tag: ', &
+            string
+       print *, 'program terminated by assert_eqn'
+       stop 1
+    end if
+  END FUNCTION assert_eqn
+
+END MODULE assert_eq_m
Index: LMDZ4/branches/LMDZ4-dev/libf/bibio/assert_m.f90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/bibio/assert_m.f90	(revision 1154)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/assert_m.f90	(revision 1154)
@@ -0,0 +1,69 @@
+! $Id$
+MODULE assert_m
+
+  implicit none
+
+  INTERFACE assert
+     MODULE PROCEDURE assert1,assert2,assert3,assert4,assert_v
+  END INTERFACE
+
+  private assert1,assert2,assert3,assert4,assert_v
+
+CONTAINS
+
+  SUBROUTINE assert1(n1,string)
+    CHARACTER(LEN=*), INTENT(IN) :: string
+    LOGICAL, INTENT(IN) :: n1
+    if (.not. n1) then
+       write (*,*) 'nrerror: an assertion failed with this tag:', &
+            string
+       print *, 'program terminated by assert1'
+       stop 1
+    end if
+  END SUBROUTINE assert1
+  !BL
+  SUBROUTINE assert2(n1,n2,string)
+    CHARACTER(LEN=*), INTENT(IN) :: string
+    LOGICAL, INTENT(IN) :: n1,n2
+    if (.not. (n1 .and. n2)) then
+       write (*,*) 'nrerror: an assertion failed with this tag:', &
+            string
+       print *, 'program terminated by assert2'
+       stop 1
+    end if
+  END SUBROUTINE assert2
+  !BL
+  SUBROUTINE assert3(n1,n2,n3,string)
+    CHARACTER(LEN=*), INTENT(IN) :: string
+    LOGICAL, INTENT(IN) :: n1,n2,n3
+    if (.not. (n1 .and. n2 .and. n3)) then
+       write (*,*) 'nrerror: an assertion failed with this tag:', &
+            string
+       print *, 'program terminated by assert3'
+       stop 1
+    end if
+  END SUBROUTINE assert3
+  !BL
+  SUBROUTINE assert4(n1,n2,n3,n4,string)
+    CHARACTER(LEN=*), INTENT(IN) :: string
+    LOGICAL, INTENT(IN) :: n1,n2,n3,n4
+    if (.not. (n1 .and. n2 .and. n3 .and. n4)) then
+       write (*,*) 'nrerror: an assertion failed with this tag:', &
+            string
+       print *, 'program terminated by assert4'
+       stop 1
+    end if
+  END SUBROUTINE assert4
+  !BL
+  SUBROUTINE assert_v(n,string)
+    CHARACTER(LEN=*), INTENT(IN) :: string
+    LOGICAL, DIMENSION(:), INTENT(IN) :: n
+    if (.not. all(n)) then
+       write (*,*) 'nrerror: an assertion failed with this tag:', &
+            string
+       print *, 'program terminated by assert_v'
+       stop 1
+    end if
+  END SUBROUTINE assert_v
+
+END MODULE assert_m
Index: LMDZ4/branches/LMDZ4-dev/libf/bibio/handle_err_m.f90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/bibio/handle_err_m.f90	(revision 1154)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/handle_err_m.f90	(revision 1154)
@@ -0,0 +1,46 @@
+! $Id$
+module handle_err_m
+
+  implicit none
+
+contains
+
+  subroutine handle_err(message, ncerr, ncid, varid)
+
+    use netcdf, only: nf90_strerror, nf90_noerr, nf90_close
+
+    character(len=*), intent(in):: message
+    ! (should include name of calling procedure)
+
+    integer, intent(in):: ncerr
+
+    integer, intent(in), optional :: ncid
+    ! (Provide this argument if you want "handle_err" to try to close
+    ! the file.)
+
+    integer, intent(in), optional :: varid
+
+    ! Variable local to the procedure:
+    integer ncerr_close
+
+    !-------------------
+
+    if (ncerr /= nf90_noerr) then
+       print *, message, ":"
+       if (present(varid)) print *, "varid = ", varid
+       print *, trim(nf90_strerror(ncerr))
+       if (present(ncid)) then
+          ! Try to close, to leave the file in a consistent state:
+          ncerr_close = nf90_close(ncid)
+          ! (do not call "nf95_close", we do not want to recurse)
+          if (ncerr_close /= nf90_noerr) then
+             print *, "nf90_close:"
+             print *, trim(nf90_strerror(ncerr_close))
+          end if
+       end if
+       stop 1
+    end if
+
+  end subroutine handle_err
+
+end module handle_err_m
Index: LMDZ4/branches/LMDZ4-dev/libf/bibio/interpolation.f90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/bibio/interpolation.f90	(revision 1154)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/interpolation.f90	(revision 1154)
@@ -0,0 +1,138 @@
+! $Id$
+module interpolation
+
+  ! From Press et al., 1996, version 2.10a
+  ! B3 Interpolation and Extrapolation
+
+  IMPLICIT NONE 
+
+contains
+
+  pure FUNCTION locate(xx,x)
+
+    REAL, DIMENSION(:), INTENT(IN) :: xx
+    REAL, INTENT(IN) :: x
+    INTEGER  locate
+
+    ! Given an array xx(1:N), and given a value x, returns a value j,
+    ! between 0 and N, such that x is between xx(j) and xx(j + 1). xx
+    ! must be monotonic, either increasing or decreasing. j = 0 or j =
+    ! N is returned to indicate that x is out of range. This
+    ! procedure should not be called with a zero-sized array argument.
+    ! See notes.
+
+    INTEGER  n,jl,jm,ju
+    LOGICAL  ascnd
+
+    !----------------------------
+
+    n=size(xx)
+    ascnd = (xx(n) >= xx(1))
+    ! (True if ascending order of table, false otherwise.)
+    ! Initialize lower and upper limits:
+    jl=0
+    ju=n+1
+    do while (ju-jl > 1)
+       jm=(ju+jl)/2 ! Compute a midpoint,
+       if (ascnd .eqv. (x >= xx(jm))) then
+          jl=jm ! and replace either the lower limit
+       else
+          ju=jm ! or the upper limit, as appropriate.
+       end if
+    end do
+    ! {ju == jl + 1}
+
+    ! {(ascnd .and. xx(jl) <= x < xx(jl+1)) 
+    !  .neqv. 
+    !  (.not. ascnd .and. xx(jl+1) <= x < xx(jl))}
+
+    ! Then set the output, being careful with the endpoints:
+    if (x == xx(1)) then
+       locate=1
+    else if (x == xx(n)) then
+       locate=n-1
+    else
+       locate=jl
+    end if
+
+  END FUNCTION locate
+
+  !***************************
+
+  pure SUBROUTINE hunt(xx,x,jlo)
+
+    ! Given an array xx(1:N ), and given a value x, returns a value
+    ! jlo such that x is between xx(jlo) and xx(jlo+1). xx must be
+    ! monotonic, either increasing or decreasing. jlo = 0 or jlo = N is
+    ! returned to indicate that x is out of range. jlo on input is taken as
+    ! the initial guess for jlo on output.
+    ! Modified so that it uses the information "jlo = 0" on input.
+
+    INTEGER, INTENT(INOUT) :: jlo
+    REAL, INTENT(IN) :: x
+    REAL, DIMENSION(:), INTENT(IN) :: xx
+    INTEGER  n,inc,jhi,jm
+    LOGICAL  ascnd, hunt_up
+
+    !-----------------------------------------------------
+
+    n=size(xx)
+    ascnd = (xx(n) >= xx(1))
+    ! (True if ascending order of table, false otherwise.)
+    if (jlo < 0 .or. jlo > n) then
+       ! Input guess not useful. Go immediately to bisection.
+       jlo=0
+       jhi=n+1
+    else
+       inc=1 ! Set the hunting increment.
+       if (jlo == 0) then
+          hunt_up = .true.
+       else
+          hunt_up = x >= xx(jlo) .eqv. ascnd
+       end if
+       if (hunt_up) then ! Hunt up:
+          do
+             jhi=jlo+inc
+             if (jhi > n) then ! Done hunting, since off end of table.
+                jhi=n+1
+                exit
+             else
+                if (x < xx(jhi) .eqv. ascnd) exit
+                jlo=jhi ! Not done hunting,
+                inc=inc+inc ! so double the increment
+             end if
+          end do ! and try again.
+       else ! Hunt down:
+          jhi=jlo
+          do
+             jlo=jhi-inc
+             if (jlo < 1) then ! Done hunting, since off end of table.
+                jlo=0
+                exit
+             else
+                if (x >= xx(jlo) .eqv. ascnd) exit
+                jhi=jlo ! Not done hunting,
+                inc=inc+inc ! so double the increment
+             end if
+          end do ! and try again.
+       end if
+    end if ! Done hunting, value bracketed.
+
+    do ! Hunt is done, so begin the final bisection phase:
+       if (jhi-jlo <= 1) then
+          if (x == xx(n)) jlo=n-1
+          if (x == xx(1)) jlo=1
+          exit
+       else
+          jm=(jhi+jlo)/2
+          if (x >= xx(jm) .eqv. ascnd) then
+             jlo=jm
+          else
+             jhi=jm
+          end if
+       end if
+    end do
+
+  END SUBROUTINE hunt
+
+end module interpolation
Index: LMDZ4/branches/LMDZ4-dev/libf/bibio/netcdf95.f90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/bibio/netcdf95.f90	(revision 1154)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/netcdf95.f90	(revision 1154)
@@ -0,0 +1,40 @@
+! $Id$
+module netcdf95
+
+  ! Author: Lionel GUEZ
+
+  ! Three criticisms may be made about the Fortran 90 NetCDF interface:
+
+  ! -- NetCDF procedures are usually functions with side effects.
+  ! First, they have "intent(out)" arguments.
+  ! Furthermore, there is obviously data transfer inside the procedures.
+  ! Any data transfer inside a function is considered as a side effect.
+
+  ! -- The caller of a NetCDF procedure usually has to handle the error
+  ! status. NetCDF procedures would be much friendlier if they behaved
+  ! like the Fortran input/output statements. That is, the error status
+  ! should be an optional output argument.
+  ! If the caller does not request the error status and there is an
+  ! error then the NetCDF procedure should produce an error message
+  ! and stop the program.
+
+  ! -- Some procedures use array arguments with assumed size.
+  ! It would be better to use the pointer attribute.
+
+  ! This module produces a NetCDF interface that answers those three
+  ! criticisms for some (not all) procedures.
+
+  ! This module replaces some of the official NetCDF procedures.
+  ! This module also provides the procedures "handle_err" and "nf95_gw_var".
+
+  ! This module provides only a partial replacement for some generic
+  ! procedures such as "nf90_def_var".
+
+  use nf95_def_var_m
+  use nf95_put_var_m
+  use nf95_gw_var_m
+  use nf95_put_att_m
+  use simple
+  use handle_err_m
+
+end module netcdf95
Index: LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_def_var_m.f90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_def_var_m.f90	(revision 1154)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_def_var_m.f90	(revision 1154)
@@ -0,0 +1,69 @@
+! $Id$
+module nf95_def_var_m
+
+  implicit none
+
+  interface nf95_def_var
+    module procedure nf95_def_var_oneDim, nf95_def_var_ManyDims
+  end interface
+
+  private
+  public nf95_def_var
+
+contains
+
+  subroutine nf95_def_var_oneDim(ncid, name, xtype, dimids, 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( in) :: dimids
+    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, dimids, varid)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_def_var_oneDim " // name, ncerr_not_opt, ncid)
+    end if
+
+  end subroutine nf95_def_var_oneDim
+
+  !***********************
+
+  subroutine nf95_def_var_ManyDims(ncid, name, xtype, dimids, 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, dimension(:), intent( in) :: dimids
+    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, dimids, varid)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_def_var_ManyDims " // name, ncerr_not_opt, ncid)
+    end if
+
+  end subroutine nf95_def_var_ManyDims
+
+end module nf95_def_var_m
Index: LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_gw_var_m.f90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_gw_var_m.f90	(revision 1154)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_gw_var_m.f90	(revision 1154)
@@ -0,0 +1,338 @@
+! $Id$
+module nf95_gw_var_m
+
+  implicit none
+
+  interface nf95_gw_var
+     ! "nf95_gw_var" stands for "NetCDF 1995 get whole variable".
+     ! These procedures read a whole NetCDF variable (coordinate or
+     ! primary) into an array.
+     ! The difference between the procedures is the rank of the array
+     ! and the type of Fortran values.
+     ! The procedures do not check the type of the NetCDF variable.
+
+!!$     module procedure nf95_gw_var_real_1d, nf95_gw_var_real_2d, &
+!!$          nf95_gw_var_real_3d, nf95_gw_var_real_4d, nf95_gw_var_dble_1d, &
+!!$          nf95_gw_var_dble_3d, nf95_gw_var_int_1d, nf95_gw_var_int_3d
+     module procedure nf95_gw_var_real_1d, nf95_gw_var_real_2d, &
+          nf95_gw_var_real_3d, nf95_gw_var_real_4d, nf95_gw_var_int_1d, &
+          nf95_gw_var_int_3d
+  end interface
+
+  private
+  public nf95_gw_var
+
+contains
+
+  subroutine nf95_gw_var_real_1d(ncid, varid, values)
+
+    ! Real type, the array has rank 1.
+
+    use netcdf, only: NF90_GET_VAR
+    use simple, only: nf95_inquire_variable, nf95_inquire_dimension
+    use handle_err_m, only: handle_err
+
+    integer, intent(in):: ncid
+    integer, intent(in):: varid
+    real, pointer:: values(:)
+
+    ! Variables local to the procedure:
+    integer ierr, len
+    integer, pointer :: dimids(:)
+
+    !---------------------
+
+    call nf95_inquire_variable(ncid, varid, dimids=dimids)
+
+    if (size(dimids) /= 1) then
+       print *, "nf95_gw_var_real_1d: NetCDF variable is not of rank 1"
+       stop 1
+    end if
+
+    call nf95_inquire_dimension(ncid, dimids(1), len=len)
+    deallocate(dimids) ! pointer
+
+    allocate(values(len))
+    if (len /= 0) then
+       ierr = NF90_GET_VAR(ncid, varid, values)
+       call handle_err("NF90_GET_VAR", ierr, ncid, varid)
+    end if
+
+  end subroutine nf95_gw_var_real_1d
+
+  !************************************
+
+  subroutine nf95_gw_var_real_2d(ncid, varid, values)
+
+    ! Real type, the array has rank 2.
+
+    use netcdf, only: NF90_GET_VAR
+    use simple, only: nf95_inquire_variable, nf95_inquire_dimension
+    use handle_err_m, only: handle_err
+
+    integer, intent(in):: ncid
+    integer, intent(in):: varid
+    real, pointer:: values(:, :)
+
+    ! Variables local to the procedure:
+    integer ierr, len1, len2
+    integer, pointer :: dimids(:)
+
+    !---------------------
+
+    call nf95_inquire_variable(ncid, varid, dimids=dimids)
+
+    if (size(dimids) /= 2) then
+       print *, "nf95_gw_var_real_2d: NetCDF variable is not of rank 2"
+       stop 1
+    end if
+
+    call nf95_inquire_dimension(ncid, dimids(1), len=len1)
+    call nf95_inquire_dimension(ncid, dimids(2), len=len2)
+    deallocate(dimids) ! pointer
+
+    allocate(values(len1, len2))
+    if (len1 /= 0 .and. len2 /= 0) then
+       ierr = NF90_GET_VAR(ncid, varid, values)
+       call handle_err("NF90_GET_VAR", ierr, ncid, varid)
+    end if
+
+  end subroutine nf95_gw_var_real_2d
+
+  !************************************
+
+  subroutine nf95_gw_var_real_3d(ncid, varid, values)
+
+    ! Real type, the array has rank 3.
+
+    use netcdf, only: NF90_GET_VAR
+    use simple, only: nf95_inquire_variable, nf95_inquire_dimension
+    use handle_err_m, only: handle_err
+
+    integer, intent(in):: ncid
+    integer, intent(in):: varid
+    real, pointer:: values(:, :, :)
+
+    ! Variables local to the procedure:
+    integer ierr, len1, len2, len3
+    integer, pointer :: dimids(:)
+
+    !---------------------
+
+    call nf95_inquire_variable(ncid, varid, dimids=dimids)
+
+    if (size(dimids) /= 3) then
+       print *, "nf95_gw_var_real_3d: NetCDF variable is not of rank 3"
+       stop 1
+    end if
+
+    call nf95_inquire_dimension(ncid, dimids(1), len=len1)
+    call nf95_inquire_dimension(ncid, dimids(2), len=len2)
+    call nf95_inquire_dimension(ncid, dimids(3), len=len3)
+    deallocate(dimids) ! pointer
+
+    allocate(values(len1, len2, len3))
+    if (len1 * len2 * len3 /= 0) then
+       ierr = NF90_GET_VAR(ncid, varid, values)
+       call handle_err("NF90_GET_VAR", ierr, ncid, varid)
+    end if
+
+  end subroutine nf95_gw_var_real_3d
+
+  !************************************
+
+  subroutine nf95_gw_var_real_4d(ncid, varid, values)
+
+    ! Real type, the array has rank 4.
+
+    use netcdf, only: NF90_GET_VAR
+    use simple, only: nf95_inquire_variable, nf95_inquire_dimension
+    use handle_err_m, only: handle_err
+
+    integer, intent(in):: ncid
+    integer, intent(in):: varid
+    real, pointer:: values(:, :, :, :)
+
+    ! Variables local to the procedure:
+    integer ierr, len_dim(4), i
+    integer, pointer :: dimids(:)
+
+    !---------------------
+
+    call nf95_inquire_variable(ncid, varid, dimids=dimids)
+
+    if (size(dimids) /= 4) then
+       print *, "nf95_gw_var_real_4d: NetCDF variable is not of rank 4"
+       stop 1
+    end if
+
+    do i = 1, 4
+       call nf95_inquire_dimension(ncid, dimids(i), len=len_dim(i))
+    end do
+    deallocate(dimids) ! pointer
+
+    allocate(values(len_dim(1), len_dim(2), len_dim(3), len_dim(4)))
+    if (all(len_dim /= 0)) then
+       ierr = NF90_GET_VAR(ncid, varid, values)
+       call handle_err("NF90_GET_VAR", ierr, ncid, varid)
+    end if
+
+  end subroutine nf95_gw_var_real_4d
+
+  !************************************
+
+!!$  subroutine nf95_gw_var_dble_1d(ncid, varid, values)
+!!$
+!!$    ! Double precision, the array has rank 1.
+!!$
+!!$    use netcdf, only: NF90_GET_VAR
+!!$    use simple, only: nf95_inquire_variable, nf95_inquire_dimension
+!!$    use handle_err_m, only: handle_err
+!!$
+!!$    integer, intent(in):: ncid
+!!$    integer, intent(in):: varid
+!!$    double precision, pointer:: values(:)
+!!$
+!!$    ! Variables local to the procedure:
+!!$    integer ierr, len
+!!$    integer, pointer :: dimids(:)
+!!$
+!!$    !---------------------
+!!$
+!!$    call nf95_inquire_variable(ncid, varid, dimids=dimids)
+!!$
+!!$    if (size(dimids) /= 1) then
+!!$       print *, "nf95_gw_var_dble_1d: NetCDF variable is not of rank 1"
+!!$       stop 1
+!!$    end if
+!!$
+!!$    call nf95_inquire_dimension(ncid, dimids(1), len=len)
+!!$    deallocate(dimids) ! pointer
+!!$
+!!$    allocate(values(len))
+!!$    if (len /= 0) then
+!!$       ierr = NF90_GET_VAR(ncid, varid, values)
+!!$       call handle_err("NF90_GET_VAR", ierr, ncid, varid)
+!!$    end if
+!!$
+!!$  end subroutine nf95_gw_var_dble_1d
+!!$
+!!$  !************************************
+!!$
+!!$  subroutine nf95_gw_var_dble_3d(ncid, varid, values)
+!!$
+!!$    ! Double precision, the array has rank 3.
+!!$
+!!$    use netcdf, only: NF90_GET_VAR
+!!$    use simple, only: nf95_inquire_variable, nf95_inquire_dimension
+!!$    use handle_err_m, only: handle_err
+!!$
+!!$    integer, intent(in):: ncid
+!!$    integer, intent(in):: varid
+!!$    double precision, pointer:: values(:, :, :)
+!!$
+!!$    ! Variables local to the procedure:
+!!$    integer ierr, len1, len2, len3
+!!$    integer, pointer :: dimids(:)
+!!$
+!!$    !---------------------
+!!$
+!!$    call nf95_inquire_variable(ncid, varid, dimids=dimids)
+!!$
+!!$    if (size(dimids) /= 3) then
+!!$       print *, "nf95_gw_var_dble_3d: NetCDF variable is not of rank 3"
+!!$       stop 1
+!!$    end if
+!!$
+!!$    call nf95_inquire_dimension(ncid, dimids(1), len=len1)
+!!$    call nf95_inquire_dimension(ncid, dimids(2), len=len2)
+!!$    call nf95_inquire_dimension(ncid, dimids(3), len=len3)
+!!$    deallocate(dimids) ! pointer
+!!$
+!!$    allocate(values(len1, len2, len3))
+!!$    if (len1 * len2 * len3 /= 0) then
+!!$       ierr = NF90_GET_VAR(ncid, varid, values)
+!!$       call handle_err("NF90_GET_VAR", ierr, ncid, varid)
+!!$    end if
+!!$
+!!$  end subroutine nf95_gw_var_dble_3d
+
+  !************************************
+
+  subroutine nf95_gw_var_int_1d(ncid, varid, values)
+
+    ! Integer type, the array has rank 1.
+
+    use netcdf, only: NF90_GET_VAR
+    use simple, only: nf95_inquire_variable, nf95_inquire_dimension
+    use handle_err_m, only: handle_err
+
+    integer, intent(in):: ncid
+    integer, intent(in):: varid
+    integer, pointer:: values(:)
+
+    ! Variables local to the procedure:
+    integer ierr, len
+    integer, pointer :: dimids(:)
+
+    !---------------------
+
+    call nf95_inquire_variable(ncid, varid, dimids=dimids)
+
+    if (size(dimids) /= 1) then
+       print *, "nf95_gw_var_int_1d: NetCDF variable is not of rank 1"
+       stop 1
+    end if
+
+    call nf95_inquire_dimension(ncid, dimids(1), len=len)
+    deallocate(dimids) ! pointer
+
+    allocate(values(len))
+    if (len /= 0) then
+       ierr = NF90_GET_VAR(ncid, varid, values)
+       call handle_err("NF90_GET_VAR", ierr, ncid, varid)
+    end if
+
+  end subroutine nf95_gw_var_int_1d
+
+  !************************************
+
+  subroutine nf95_gw_var_int_3d(ncid, varid, values)
+
+    ! Integer type, the array has rank 3.
+
+    use netcdf, only: NF90_GET_VAR
+    use simple, only: nf95_inquire_variable, nf95_inquire_dimension
+    use handle_err_m, only: handle_err
+
+    integer, intent(in):: ncid
+    integer, intent(in):: varid
+    integer, pointer:: values(:, :, :)
+
+    ! Variables local to the procedure:
+    integer ierr, len1, len2, len3
+    integer, pointer :: dimids(:)
+
+    !---------------------
+
+    call nf95_inquire_variable(ncid, varid, dimids=dimids)
+
+    if (size(dimids) /= 3) then
+       print *, "nf95_gw_var_int_3d: NetCDF variable is not of rank 3"
+       stop 1
+    end if
+
+    call nf95_inquire_dimension(ncid, dimids(1), len=len1)
+    call nf95_inquire_dimension(ncid, dimids(2), len=len2)
+    call nf95_inquire_dimension(ncid, dimids(3), len=len3)
+    deallocate(dimids) ! pointer
+
+    allocate(values(len1, len2, len3))
+    if (len1 * len2 * len3 /= 0) then
+       ierr = NF90_GET_VAR(ncid, varid, values)
+       call handle_err("NF90_GET_VAR", ierr, ncid, varid)
+    end if
+
+  end subroutine nf95_gw_var_int_3d
+
+end module nf95_gw_var_m
Index: LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_put_att_m.f90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_put_att_m.f90	(revision 1154)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_put_att_m.f90	(revision 1154)
@@ -0,0 +1,67 @@
+! $Id$
+module nf95_put_att_m
+
+  implicit none
+
+  interface nf95_put_att
+     module procedure nf95_put_att_text, nf95_put_att_one_FourByteInt
+  end interface
+
+  private
+  public nf95_put_att
+
+contains
+
+  subroutine nf95_put_att_text(ncid, varid, name, values, ncerr)
+
+    use netcdf, only: nf90_put_att
+    use handle_err_m, only: handle_err
+
+    integer, intent(in) :: ncid, varid
+    character(len = *), intent(in) :: name
+    character(len = *), intent(in) :: values
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_put_att(ncid, varid, name, values)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_put_att_text", ncerr_not_opt, ncid, varid)
+    end if
+
+  end subroutine nf95_put_att_text
+
+  !************************************
+
+  subroutine nf95_put_att_one_FourByteInt(ncid, varid, name, values, ncerr)
+
+    use netcdf, only: nf90_put_att
+    use handle_err_m, only: handle_err
+    use typesizes, only: FourByteInt
+
+    integer, intent(in) :: ncid, varid
+    character(len = *), intent(in) :: name
+    integer(kind = FourByteInt), intent(in) :: values
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_put_att(ncid, varid, name, values)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_put_att_one_FourByteInt", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_put_att_one_FourByteInt
+
+end module nf95_put_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 1154)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_put_var_m.f90	(revision 1154)
@@ -0,0 +1,199 @@
+! $Id$
+module nf95_put_var_m
+
+  implicit none
+
+  interface nf95_put_var
+     module procedure nf95_put_var_1D_FourByteReal, &
+          nf95_put_var_2D_FourByteReal, nf95_put_var_3D_FourByteReal, &
+          nf95_put_var_4D_FourByteReal
+!!$     module procedure nf95_put_var_1D_FourByteReal, &
+!!$          nf95_put_var_2D_FourByteReal, nf95_put_var_3D_FourByteReal, &
+!!$          nf95_put_var_4D_FourByteReal, nf90_put_var_1D_EightByteReal, &
+!!$          nf90_put_var_3D_EightByteReal
+  end interface
+
+  private
+  public nf95_put_var
+
+contains
+
+  subroutine nf95_put_var_1D_FourByteReal(ncid, varid, values, start, count, &
+       stride, map, ncerr)
+
+    use typesizes, only: FourByteReal
+    use netcdf, only: nf90_put_var
+    use handle_err_m, only: handle_err
+
+    integer,                         intent(in) :: ncid, varid
+    real(kind = FourByteReal), 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_FourByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_put_var_1D_FourByteReal
+
+  !***********************
+
+  subroutine nf95_put_var_2D_FourByteReal(ncid, varid, values, start, count, &
+       stride, map, ncerr)
+
+    use typesizes, only: FourByteReal
+    use netcdf, only: nf90_put_var
+    use handle_err_m, only: handle_err
+
+    integer,                         intent( in) :: ncid, varid
+    real (kind = FourByteReal), 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_2D_FourByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_put_var_2D_FourByteReal
+
+  !***********************
+
+  subroutine nf95_put_var_3D_FourByteReal(ncid, varid, values, start, count, &
+       stride, map, ncerr)
+
+    use typesizes, only: FourByteReal
+    use netcdf, only: nf90_put_var
+    use handle_err_m, only: handle_err
+
+    integer,                         intent( in) :: ncid, varid
+    real (kind = FourByteReal), 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_3D_FourByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_put_var_3D_FourByteReal
+
+  !***********************
+
+  subroutine nf95_put_var_4D_FourByteReal(ncid, varid, values, start, count, &
+       stride, map, ncerr)
+
+    use typesizes, only: FourByteReal
+    use netcdf, only: nf90_put_var
+    use handle_err_m, only: handle_err
+
+    integer,                         intent( in) :: ncid, varid
+    real (kind = FourByteReal), 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_4D_FourByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_put_var_4D_FourByteReal
+
+  !***********************
+
+!!$  subroutine nf90_put_var_1D_EightByteReal(ncid, varid, values, start, count, &
+!!$       stride, map, ncerr)
+!!$
+!!$    use typesizes, only: eightByteReal
+!!$    use netcdf, only: nf90_put_var
+!!$    use handle_err_m, only: handle_err
+!!$
+!!$    integer,                         intent( in) :: ncid, varid
+!!$    real (kind = EightByteReal),     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_eightByteReal", ncerr_not_opt, ncid, &
+!!$            varid)
+!!$    end if
+!!$
+!!$  end subroutine nf90_put_var_1D_EightByteReal
+!!$
+!!$  !***********************
+!!$
+!!$  subroutine nf90_put_var_3D_EightByteReal(ncid, varid, values, start, count, &
+!!$       stride, map, ncerr)
+!!$
+!!$    use typesizes, only: eightByteReal
+!!$    use netcdf, only: nf90_put_var
+!!$    use handle_err_m, only: handle_err
+!!$
+!!$    integer,                         intent( in) :: ncid, varid
+!!$    real (kind = EightByteReal),     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_3D_eightByteReal", ncerr_not_opt, ncid, &
+!!$            varid)
+!!$    end if
+!!$
+!!$  end subroutine nf90_put_var_3D_EightByteReal
+
+end module nf95_put_var_m
Index: LMDZ4/branches/LMDZ4-dev/libf/bibio/regr1_lint_m.f90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/bibio/regr1_lint_m.f90	(revision 1154)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/regr1_lint_m.f90	(revision 1154)
@@ -0,0 +1,98 @@
+! $Id$
+module regr1_lint_m
+
+  ! Author: Lionel GUEZ
+
+  implicit none
+
+  interface regr1_lint
+     ! Each procedure regrids by linear interpolation.
+     ! The regridding operation is done on the first dimension of the
+     ! input array.
+     ! The difference betwwen the procedures is the rank of the first argument.
+     module procedure regr11_lint, regr12_lint
+  end interface
+
+  private
+  public regr1_lint
+
+contains
+
+  function regr11_lint(vs, xs, xt) result(vt)
+
+    ! "vs" has rank 1.
+
+    use assert_eq_m, only: assert_eq
+    use interpolation, only: hunt !!, polint
+
+    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(xt)) ! 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), size(xs), "regr11_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)
+!!       call polint(xs(is_b:is_b+1), vs(is_b:is_b+1), xt(it), vt(it))
+       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 regr11_lint
+
+  !*********************************************************
+
+  function regr12_lint(vs, xs, xt) result(vt)
+
+    ! "vs" has rank 2.
+
+    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(xt), size(vs, 2)) ! 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, 1), size(xs), "regr12_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 regr12_lint
+
+end module regr1_lint_m
Index: LMDZ4/branches/LMDZ4-dev/libf/bibio/regr1_step_av_m.f90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/bibio/regr1_step_av_m.f90	(revision 1154)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/regr1_step_av_m.f90	(revision 1154)
@@ -0,0 +1,205 @@
+! $Id$
+module regr1_step_av_m
+
+  ! Author: Lionel GUEZ
+
+  implicit none
+
+  interface regr1_step_av
+
+     ! Each procedure regrids a step function by averaging it.
+     ! The regridding operation is done on the first dimension of the
+     ! input array.
+     ! Source grid contains edges of steps.
+     ! Target grid contains positions of cell edges.
+     ! The target grid should be included in the source grid: no
+     ! extrapolation is allowed.
+     ! The difference between the procedures is the rank of the first argument.
+
+     module procedure regr11_step_av, regr12_step_av, regr13_step_av
+  end interface
+
+  private
+  public regr1_step_av
+
+contains
+
+  function regr11_step_av(vs, xs, xt) result(vt)
+
+    ! "vs" has rank 1.
+
+    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 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) ! 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), size(xs) - 1, "regr11_step_av ns")
+    nt = size(xt) - 1
+    ! Quick check on sort order:
+    call assert(xs(1) < xs(2), "regr11_step_av xs bad order")
+    call assert(xt(1) < xt(2), "regr11_step_av xt bad order")
+
+    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
+         "regr11_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 regr11_step_av
+
+  !********************************************
+
+  function regr12_step_av(vs, xs, xt) result(vt)
+
+    ! "vs" has rank 2.
+
+    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)) ! 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, "regr12_step_av ns")
+    nt = size(xt) - 1
+
+    ! Quick check on sort order:
+    call assert(xs(1) < xs(2), "regr12_step_av xs bad order")
+    call assert(xt(1) < xt(2), "regr12_step_av xt bad order")
+
+    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
+         "regr12_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 regr12_step_av
+
+  !********************************************
+
+  function regr13_step_av(vs, xs, xt) result(vt)
+
+    ! "vs" has rank 3.
+
+    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)) 
+    ! (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, "regr13_step_av ns")
+    nt = size(xt) - 1
+
+    ! Quick check on sort order:
+    call assert(xs(1) < xs(2), "regr13_step_av xs bad order")
+    call assert(xt(1) < xt(2), "regr13_step_av xt bad order")
+
+    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
+         "regr13_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 regr13_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 1154)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/regr3_lint_m.f90	(revision 1154)
@@ -0,0 +1,59 @@
+! $Id$
+module regr3_lint_m
+
+  ! Author: Lionel GUEZ
+
+  implicit none
+
+  interface regr3_lint
+     ! Each procedure regrids by linear interpolation.
+     ! The regridding operation is done on the third dimension of the
+     ! input array.
+     ! The difference betwwen the procedures is the rank of the first argument.
+     module procedure regr33_lint
+  end interface
+
+  private
+  public regr3_lint
+
+contains
+
+  function regr33_lint(vs, xs, xt) result(vt)
+
+    ! "vs" has rank 3.
+
+    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))
+    ! (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), "regr33_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 regr33_lint
+
+end module regr3_lint_m
Index: LMDZ4/branches/LMDZ4-dev/libf/bibio/simple.f90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/bibio/simple.f90	(revision 1154)
+++ LMDZ4/branches/LMDZ4-dev/libf/bibio/simple.f90	(revision 1154)
@@ -0,0 +1,312 @@
+! $Id$
+module simple
+
+  implicit none
+
+contains
+
+  subroutine nf95_open(path, mode, ncid, chunksize, ncerr)
+
+    use netcdf, only: nf90_open
+    use handle_err_m, only: handle_err
+
+    character(len=*), intent(in):: path
+    integer, intent(in):: mode
+    integer, intent(out):: ncid
+    integer, intent(inout), optional:: chunksize
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_open(path, mode, ncid, chunksize)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_open " // path, ncerr_not_opt)
+    end if
+
+  end subroutine nf95_open
+
+  !************************
+
+  subroutine nf95_inq_dimid(ncid, name, dimid, ncerr)
+
+    use netcdf, only: nf90_inq_dimid
+    use handle_err_m, only: handle_err
+
+    integer,             intent( in) :: ncid
+    character (len = *), intent( in) :: name
+    integer,             intent(out) :: dimid
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_inq_dimid(ncid, name, dimid)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_inq_dimid", ncerr_not_opt, ncid)
+    end if
+
+  end subroutine nf95_inq_dimid
+
+  !************************
+
+  subroutine nf95_inquire_dimension(ncid, dimid, name, len, ncerr)
+
+    use netcdf, only: nf90_inquire_dimension
+    use handle_err_m, only: handle_err
+
+    integer,                       intent( in) :: ncid, dimid
+    character (len = *), optional, intent(out) :: name
+    integer,             optional, intent(out) :: len
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_inquire_dimension(ncid, dimid, name, len)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_inquire_dimension", ncerr_not_opt, ncid)
+    end if
+
+  end subroutine nf95_inquire_dimension
+
+  !************************
+
+  subroutine nf95_inq_varid(ncid, name, varid, ncerr)
+
+    use netcdf, only: nf90_inq_varid
+    use handle_err_m, only: handle_err
+
+    integer,             intent(in) :: ncid
+    character (len = *), intent(in) :: name
+    integer,             intent(out) :: varid
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_inq_varid(ncid, name, varid)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_inq_varid, name = " // name, ncerr_not_opt, ncid)
+    end if
+
+  end subroutine nf95_inq_varid
+
+  !************************
+
+  subroutine nf95_inquire_variable(ncid, varid, name, xtype, ndims, dimids, &
+       nAtts, ncerr)
+
+    ! In "nf90_inquire_variable", "dimids" is an assumed-size array.
+    ! This is the classical case of an array the size of which is
+    ! 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.
+
+    use netcdf, only: nf90_inquire_variable, nf90_max_var_dims
+    use handle_err_m, only: handle_err
+
+    integer, intent(in):: ncid, varid
+    character(len = *), optional, intent(out):: name
+    integer, optional, intent(out) :: xtype, ndims
+    integer, dimension(:), optional, pointer :: dimids
+    integer, optional, intent(out) :: nAtts
+    integer, intent(out), optional :: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+    integer dimids_local(nf90_max_var_dims)
+    integer ndims_not_opt
+
+    !-------------------
+
+    if (present(dimids)) then
+       ncerr_not_opt = nf90_inquire_variable(ncid, varid, name, xtype, &
+            ndims_not_opt, dimids_local, nAtts)
+       allocate(dimids(ndims_not_opt)) ! also works if ndims_not_opt == 0
+       dimids = dimids_local(:ndims_not_opt)
+       if (present(ndims)) ndims = ndims_not_opt
+    else
+       ncerr_not_opt = nf90_inquire_variable(ncid, varid, name, xtype, ndims, &
+            nAtts=nAtts)
+    end if
+
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_inquire_variable", ncerr_not_opt, ncid)
+    end if
+
+  end subroutine nf95_inquire_variable
+
+  !************************
+
+  subroutine nf95_create(path, cmode, ncid, initialsize, chunksize, ncerr)
+    
+    use netcdf, only: nf90_create
+    use handle_err_m, only: handle_err
+
+    character (len = *), intent(in   ) :: path
+    integer,             intent(in   ) :: cmode
+    integer,             intent(  out) :: ncid
+    integer, optional,   intent(in   ) :: initialsize
+    integer, optional,   intent(inout) :: chunksize
+    integer, intent(out), optional :: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_create(path, cmode, ncid, initialsize, chunksize)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_create " // path, ncerr_not_opt)
+    end if
+
+  end subroutine nf95_create
+
+  !************************
+
+  subroutine nf95_def_dim(ncid, name, len, dimid, ncerr)
+
+    use netcdf, only: nf90_def_dim
+    use handle_err_m, only: handle_err
+
+    integer,             intent( in) :: ncid
+    character (len = *), intent( in) :: name
+    integer,             intent( in) :: len
+    integer,             intent(out) :: dimid
+    integer, intent(out), optional :: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_def_dim(ncid, name, len, dimid)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_def_dim", ncerr_not_opt, ncid)
+    end if
+
+  end subroutine nf95_def_dim
+
+  !***********************
+
+  subroutine nf95_redef(ncid, ncerr)
+
+    use netcdf, only: nf90_redef
+    use handle_err_m, only: handle_err
+
+    integer, intent( in) :: ncid
+    integer, intent(out), optional :: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_redef(ncid)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_redef", ncerr_not_opt, ncid)
+    end if
+
+  end subroutine nf95_redef
+  
+  !***********************
+
+  subroutine nf95_enddef(ncid, h_minfree, v_align, v_minfree, r_align, ncerr)
+
+    use netcdf, only: nf90_enddef
+    use handle_err_m, only: handle_err
+
+    integer,           intent( in) :: ncid
+    integer, optional, intent( in) :: h_minfree, v_align, v_minfree, r_align
+    integer, intent(out), optional :: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_enddef(ncid, h_minfree, v_align, v_minfree, r_align)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_enddef", ncerr_not_opt, ncid)
+    end if
+
+  end subroutine nf95_enddef
+
+  !***********************
+
+  subroutine nf95_close(ncid, ncerr)
+
+    use netcdf, only: nf90_close
+    use handle_err_m, only: handle_err
+
+    integer, intent( in) :: ncid
+    integer, intent(out), optional :: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_close(ncid)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_close", ncerr_not_opt)
+    end if
+
+  end subroutine nf95_close
+
+  !***********************
+
+  subroutine nf95_copy_att(ncid_in, varid_in, name, ncid_out, varid_out, ncerr)
+
+    use netcdf, only: nf90_copy_att
+    use handle_err_m, only: handle_err
+
+    integer, intent( in):: ncid_in,  varid_in
+    character(len=*), intent( in):: name
+    integer, intent( in):: ncid_out, varid_out
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_copy_att(ncid_in, varid_in, name, ncid_out, varid_out)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_copy_att", ncerr_not_opt, ncid_out)
+    end if
+
+  end subroutine nf95_copy_att
+
+end module simple
Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/etat0_netcdf.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/dyn3d/etat0_netcdf.F	(revision 1153)
+++ LMDZ4/branches/LMDZ4-dev/libf/dyn3d/etat0_netcdf.F	(revision 1154)
@@ -14,10 +14,11 @@
       USE phys_state_var_mod
       USE filtreg_mod
+      use regr_lat_time_climoz_m, only: regr_lat_time_climoz
 #endif
 !#endif of #ifdef CPP_EARTH
+      use netcdf, only: nf90_open, NF90_NOWRITE, nf90_close
       !
       IMPLICIT NONE
       !
-#include "netcdf.inc"
 #include "dimensions.h"
 #include "paramet.h"
@@ -49,8 +50,9 @@
       REAL :: vvent(iip1, jjm, llm)
       REAL :: t3d(iip1, jjp1, llm), tpot(iip1, jjp1, llm)
-      REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: q3d
       REAL :: qsat(iip1, jjp1, llm)
+      REAL,ALLOCATABLE :: q3d(:, :, :,:)
       REAL :: tsol(klon), qsol(klon), sn(klon)
-      REAL :: tsolsrf(klon,nbsrf), qsolsrf(klon,nbsrf),snsrf(klon,nbsrf) 
+!!      REAL :: tsolsrf(klon,nbsrf)
+      real qsolsrf(klon,nbsrf),snsrf(klon,nbsrf) 
       REAL :: albe(klon,nbsrf), evap(klon,nbsrf)
       REAL :: alblw(klon,nbsrf)
@@ -72,5 +74,5 @@
       !
 
-      CHARACTER*80 :: varname
+      CHARACTER(len=80) :: varname
       !
       INTEGER :: i,j, ig, l, ji,ii1,ii2
@@ -102,9 +104,9 @@
       REAL :: w(ip1jmp1,llm)
       REAL ::phystep
-      REAL :: rugsrel(iip1*jjp1)
+CC      REAL :: rugsrel(iip1*jjp1)
       REAL :: fder(klon)
-      real zrel(iip1*jjp1),chmin,chmax
-
-      CHARACTER*80 :: visu_file
+!!      real zrel(iip1*jjp1),chmin,chmax
+
+!!      CHARACTER(len=80) :: visu_file
       INTEGER :: visuid
 
@@ -140,4 +142,5 @@
       REAL      :: solarlong0
       real :: seuil_inversion
+      logical  read_climoz ! read ozone climatology
 
       !
@@ -163,6 +166,5 @@
 !      CALL defrun_new(99,.TRUE.,clesphy0)
       CALL conf_gcm( 99, .TRUE. , clesphy0 )
-      call conf_phys(ok_journe, ok_mensuel, ok_instan,                  &
-     &                 ok_hf, ok_LES,                                   &
+      call conf_phys(  ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, &
      &                 solarlong0,seuil_inversion,                      &
      &                 fact_cldcon, facttemps,ok_newmicro,iflag_radia,  &
@@ -174,17 +176,15 @@
      &                 iflag_thermals,nsplit_thermals,tau_thermals,     &
      &                 iflag_thermals_ed,iflag_thermals_optflux,        &
-     &                 iflag_coupl,iflag_clos,iflag_wake )
+     &                 iflag_coupl,iflag_clos,iflag_wake, read_climoz )
+
       dtvr   = daysec/FLOAT(day_step)
       print*,'dtvr',dtvr
 
-
-
       CALL iniconst()
       CALL inigeom()
 
 ! Initialisation pour traceurs
-      CALL infotrac_init
-      ALLOCATE(q3d(iip1,jjp1,llm,nqtot))
-
+      call infotrac_init
+      ALLOCATE(q3d(iip1, jjp1, llm, nqtot))
 
       CALL inifilr()
@@ -246,5 +246,5 @@
 
       write(*,*)'Essai de lecture masque ocean'
-      iret = nf_open("o2a.nc", NF_NOWRITE, nid_o2a)
+      iret = nf90_open("o2a.nc", NF90_NOWRITE, nid_o2a)
       if (iret .ne. 0) then
         write(*,*)'ATTENTION!! pas de fichier o2a.nc trouve'
@@ -265,5 +265,5 @@
       else
         couple = .true.
-        iret = nf_close(nid_o2a)
+        iret = nf90_close(nid_o2a)
         call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp
      $    , nid_o2a)
@@ -402,5 +402,5 @@
      .                           maxval(qsat(:,:,:))
       !
-      WRITE(*,*) 'QSAT :', qsat(10,20,:)
+CC      WRITE(*,*) 'QSAT :', qsat(10,20,:)
       !
       varname = 'q'
@@ -413,4 +413,7 @@
       q3d(:,:,:,1) = qd(:,:,:)
       !
+
+      if (read_climoz) call regr_lat_time_climoz ! ozone climatology
+
       varname = 'tsol'
       ! This line needs to be replaced by a call to restget to get the values in the restart file
@@ -477,10 +480,10 @@
      .     jjm, rlonu, rlatv , interbar )
 c
-      rugsrel(:) = 0.0
-      IF(ok_orodr)  THEN
-        DO i = 1, iip1* jjp1
-         rugsrel(i) = MAX( 1.e-05, zstd(i)* zsig(i) /2. )
-        ENDDO
-      ENDIF
+cc      rugsrel(:) = 0.0
+cc      IF(ok_orodr)  THEN
+cc        DO i = 1, iip1* jjp1
+cc         rugsrel(i) = MAX( 1.e-05, zstd(i)* zsig(i) /2. )
+cc        ENDDO
+cc      ENDIF
 
 
@@ -716,4 +719,5 @@
       q_ancien = 0.
       agesno = 0.
+c
       frugs(1:klon,is_oce) = rugmer(1:klon)
       frugs(1:klon,is_ter) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
@@ -752,22 +756,20 @@
 
 C     Sortie Visu pour les champs dynamiques
-      if (1.eq.0 ) then
-      print*,'sortie visu'
-      time_step = 1.
-      t_ops = 2.
-      t_wrt = 2.
-      itau = 2.
-      visu_file='Etat0_visu.nc'
-      CALL initdynav(visu_file,dayref,anneeref,time_step,
-     .              t_ops, t_wrt, visuid)
-      CALL writedynav(visuid, itau,vvent ,
-     .                uvent,tpot,pk,phi,q3d,masse,psol,phis)
-      else
+cc      if (1.eq.0 ) then
+cc      print*,'sortie visu'
+cc      time_step = 1.
+cc      t_ops = 2.
+cc      t_wrt = 2.
+cc      itau = 2.
+cc      visu_file='Etat0_visu.nc'
+cc      CALL initdynav(visu_file,dayref,anneeref,time_step,
+cc     .              t_ops, t_wrt, visuid)
+cc      CALL writedynav(visuid, itau,vvent ,
+cc     .                uvent,tpot,pk,phi,q3d,masse,psol,phis)
+cc      else
          print*,'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
-      endif
+cc      endif
       print*,'entree histclo'
       CALL histclo
-
-      DEALLOCATE(q3d)
 
 #endif 
@@ -776,3 +778,2 @@
       !
       END SUBROUTINE etat0_netcdf
-
Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/temps.h
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/dyn3d/temps.h	(revision 1153)
+++ LMDZ4/branches/LMDZ4-dev/libf/dyn3d/temps.h	(revision 1154)
@@ -15,6 +15,6 @@
 
       INTEGER   itaufin
-      INTEGER*4 itau_dyn, itau_phy
-      INTEGER*4 day_ini, day_end, annee_ref, day_ref
+      INTEGER(kind=4) itau_dyn, itau_phy
+      INTEGER(kind=4) day_ini, day_end, annee_ref, day_ref
       REAL      dt
 
Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/etat0_netcdf.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/etat0_netcdf.F	(revision 1153)
+++ LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/etat0_netcdf.F	(revision 1154)
@@ -9,15 +9,16 @@
       USE ioipsl
       USE dimphy
+      USE infotrac
       USE fonte_neige_mod
       USE pbl_surface_mod
       USE phys_state_var_mod
       USE filtreg_mod
-      USE infotrac
+      use regr_lat_time_climoz_m, only: regr_lat_time_climoz
 #endif
 !#endif of #ifdef CPP_EARTH
+      use netcdf, only: nf90_open, NF90_NOWRITE, nf90_close
       !
       IMPLICIT NONE
       !
-#include "netcdf.inc"
 #include "dimensions.h"
 #include "paramet.h"
@@ -43,6 +44,6 @@
       ! local variables:
       REAL :: latfi(klon), lonfi(klon)
-      REAL :: orog(iip1,jjp1), rugo(iip1,jjp1),
-     . psol(iip1, jjp1), phis(iip1, jjp1)
+      REAL :: orog(iip1,jjp1), rugo(iip1,jjp1)
+      REAL :: psol(iip1, jjp1), phis(iip1, jjp1)
       REAL :: p3d(iip1, jjp1, llm+1)
       REAL :: uvent(iip1, jjp1, llm)
@@ -52,5 +53,6 @@
       REAL,ALLOCATABLE :: q3d(:, :, :,:)
       REAL :: tsol(klon), qsol(klon), sn(klon)
-      REAL :: tsolsrf(klon,nbsrf), qsolsrf(klon,nbsrf),snsrf(klon,nbsrf) 
+!!      REAL :: tsolsrf(klon,nbsrf)
+      real qsolsrf(klon,nbsrf),snsrf(klon,nbsrf) 
       REAL :: albe(klon,nbsrf), evap(klon,nbsrf)
       REAL :: alblw(klon,nbsrf)
@@ -72,8 +74,7 @@
       !
 
-      CHARACTER*80 :: varname
+      CHARACTER(len=80) :: varname
       !
       INTEGER :: i,j, ig, l, ji,ii1,ii2
-      INTEGER :: nq
       REAL :: xpi
       !
@@ -103,9 +104,9 @@
       REAL :: w(ip1jmp1,llm)
       REAL ::phystep
-      REAL :: rugsrel(iip1*jjp1)
+CC      REAL :: rugsrel(iip1*jjp1)
       REAL :: fder(klon)
-      real zrel(iip1*jjp1),chmin,chmax
-
-      CHARACTER*80 :: visu_file
+!!      real zrel(iip1*jjp1),chmin,chmax
+
+!!      CHARACTER(len=80) :: visu_file
       INTEGER :: visuid
 
@@ -141,4 +142,5 @@
       REAL      :: solarlong0
       real :: seuil_inversion
+      logical  read_climoz ! read ozone climatology
 
       !
@@ -174,5 +176,5 @@
      &                 iflag_thermals,nsplit_thermals,tau_thermals,     &
      &                 iflag_thermals_ed,iflag_thermals_optflux,        &
-     &                 iflag_coupl,iflag_clos,iflag_wake )
+     &                 iflag_coupl,iflag_clos,iflag_wake, read_climoz )
 
       dtvr   = daysec/FLOAT(day_step)
@@ -181,10 +183,10 @@
       CALL iniconst()
       CALL inigeom()
-      !
+
       CALL inifilr()
-C init pour traceurs
+! Initialisation pour traceurs
       call infotrac_init
-      ALLOCATE(q3d(iip1, jjp1, llm,nqtot))
-!      CALL phys_state_var_init()
+      ALLOCATE(q3d(iip1, jjp1, llm, nqtot))
+!     CALL phys_state_var_init()
       !
       latfi(1) = ASIN(1.0)
@@ -243,5 +245,5 @@
 
       write(*,*)'Essai de lecture masque ocean'
-      iret = nf_open("o2a.nc", NF_NOWRITE, nid_o2a)
+      iret = nf90_open("o2a.nc", NF90_NOWRITE, nid_o2a)
       if (iret .ne. 0) then
         write(*,*)'ATTENTION!! pas de fichier o2a.nc trouve'
@@ -262,5 +264,5 @@
       else
         couple = .true.
-        iret = nf_close(nid_o2a)
+        iret = nf90_close(nid_o2a)
         call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp
      $    , nid_o2a)
@@ -399,5 +401,5 @@
      .                           maxval(qsat(:,:,:))
       !
-      WRITE(*,*) 'QSAT :', qsat(10,20,:)
+CC      WRITE(*,*) 'QSAT :', qsat(10,20,:)
       !
       varname = 'q'
@@ -410,4 +412,7 @@
       q3d(:,:,:,1) = qd(:,:,:)
       !
+
+      if (read_climoz) call regr_lat_time_climoz ! ozone climatology
+
       varname = 'tsol'
       ! This line needs to be replaced by a call to restget to get the values in the restart file
@@ -474,10 +479,10 @@
      .     jjm, rlonu, rlatv , interbar )
 c
-      rugsrel(:) = 0.0
-      IF(ok_orodr)  THEN
-        DO i = 1, iip1* jjp1
-         rugsrel(i) = MAX( 1.e-05, zstd(i)* zsig(i) /2. )
-        ENDDO
-      ENDIF
+cc      rugsrel(:) = 0.0
+cc      IF(ok_orodr)  THEN
+cc        DO i = 1, iip1* jjp1
+cc         rugsrel(i) = MAX( 1.e-05, zstd(i)* zsig(i) /2. )
+cc        ENDDO
+cc      ENDIF
 
 
@@ -750,18 +755,18 @@
 
 C     Sortie Visu pour les champs dynamiques
-      if (1.eq.0 ) then
-      print*,'sortie visu'
-      time_step = 1.
-      t_ops = 2.
-      t_wrt = 2.
-      itau = 2.
-      visu_file='Etat0_visu.nc'
-      CALL initdynav(visu_file,dayref,anneeref,time_step,
-     .              t_ops, t_wrt, visuid)
-      CALL writedynav(visuid, itau,vvent ,
-     .                uvent,tpot,pk,phi,q3d,masse,psol,phis)
-      else
+cc      if (1.eq.0 ) then
+cc      print*,'sortie visu'
+cc      time_step = 1.
+cc      t_ops = 2.
+cc      t_wrt = 2.
+cc      itau = 2.
+cc      visu_file='Etat0_visu.nc'
+cc      CALL initdynav(visu_file,dayref,anneeref,time_step,
+cc     .              t_ops, t_wrt, visuid)
+cc      CALL writedynav(visuid, itau,vvent ,
+cc     .                uvent,tpot,pk,phi,q3d,masse,psol,phis)
+cc      else
          print*,'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
-      endif
+cc      endif
       print*,'entree histclo'
       CALL histclo
@@ -772,3 +777,2 @@
       !
       END SUBROUTINE etat0_netcdf
-
Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/temps.h
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/temps.h	(revision 1153)
+++ LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/temps.h	(revision 1154)
@@ -15,6 +15,6 @@
 
       INTEGER   itaufin
-      INTEGER*4 itau_dyn, itau_phy
-      INTEGER*4 day_ini, day_end, annee_ref, day_ref
+      INTEGER(kind=4) itau_dyn, itau_phy
+      INTEGER(kind=4) day_ini, day_end, annee_ref, day_ref
       REAL      dt
 !$OMP THREADPRIVATE(/temps/)
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/clesphys.h
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/clesphys.h	(revision 1153)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/clesphys.h	(revision 1154)
@@ -13,6 +13,6 @@
        INTEGER nbapp_rad, iflag_con
        REAL co2_ppm, solaire
-       REAL*8 RCO2, RCH4, RN2O, RCFC11, RCFC12  
-       REAL*8 CH4_ppb, N2O_ppb, CFC11_ppt, CFC12_ppt
+       REAL(kind=8) RCO2, RCH4, RN2O, RCFC11, RCFC12  
+       REAL(kind=8) CH4_ppb, N2O_ppb, CFC11_ppt, CFC12_ppt
 
 !OM ---> correction du bilan d'eau global
@@ -39,5 +39,5 @@
 !IM lev_histmth : niveau sorties mensuelles
        INTEGER lev_histhf, lev_histday, lev_histmth
-       CHARACTER*4 type_run
+       CHARACTER(len=4) type_run
 ! aer_type: pour utiliser un fichier constant dans readaerosol 
        CHARACTER*8 :: aer_type 
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/conf_phys.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/conf_phys.F90	(revision 1153)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/conf_phys.F90	(revision 1154)
@@ -7,15 +7,15 @@
 
   subroutine 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, &
- &                     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, &
- &                     iflag_coupl,iflag_clos,iflag_wake )
+                       ok_LES,&
+                       solarlong0,seuil_inversion, &
+                       fact_cldcon, facttemps,ok_newmicro,iflag_radia,&
+                       iflag_cldcon, &
+                       iflag_ratqs,ratqsbas,ratqshaut, &
+  		       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, &
+                       iflag_coupl,iflag_clos,iflag_wake, read_climoz)
 
    use IOIPSL
@@ -140,4 +140,6 @@
   LOGICAL,SAVE :: ok_strato_omp
   LOGICAL,SAVE :: ok_hines_omp
+
+  logical, intent(out):: read_climoz ! read ozone climatology
 !
 
@@ -1221,4 +1223,6 @@
   call getin('ecrit_LES', ecrit_LES_omp)
 !
+  read_climoz = .false. ! default value
+  call getin('read_climoz', read_climoz)
 
 !$OMP END MASTER
@@ -1467,10 +1471,11 @@
 
   write(numout,*)' lonmin lonmax latmin latmax bilKP_ins =',&
- & lonmin_ins, lonmax_ins, latmin_ins, latmax_ins
+   lonmin_ins, lonmax_ins, latmin_ins, latmax_ins
   write(numout,*)' ecrit_ hf, day, mth, reg, tra, ISCCP, LES',&
- & ecrit_hf, ecrit_day, ecrit_mth, ecrit_reg, ecrit_tra, ecrit_ISCCP, ecrit_LES
+   ecrit_hf, ecrit_day, ecrit_mth, ecrit_reg, ecrit_tra, ecrit_ISCCP, ecrit_LES
 
   write(numout,*) 'ok_strato = ', ok_strato
   write(numout,*) 'ok_hines = ',  ok_hines
+  write(numout,*) 'read_climoz = ', read_climoz
   
 !$OMP END MASTER
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/indicesol.h
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/indicesol.h	(revision 1153)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/indicesol.h	(revision 1154)
@@ -24,5 +24,5 @@
       PARAMETER (epsfra=1.0E-05)
 !
-      CHARACTER *3 clnsurf(nbsrf)
+      CHARACTER(len=3) clnsurf(nbsrf)
       DATA clnsurf/'ter', 'lic', 'oce', 'sic'/
       SAVE clnsurf
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/open_climoz_m.f90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/open_climoz_m.f90	(revision 1154)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/open_climoz_m.f90	(revision 1154)
@@ -0,0 +1,73 @@
+! $Id$
+module open_climoz_m
+
+  implicit none
+
+contains
+
+  subroutine open_climoz(ncid, press_in_edg)
+
+    ! This procedure should be called once per "gcm" run, by a single
+    ! thread of each MPI process.
+    ! The root process opens "climoz_LMDZ.nc", reads the pressure
+    ! levels and broadcasts them to the other processes.
+
+    ! We assume that, in "climoz_LMDZ.nc", the pressure levels are in hPa
+    ! and in strictly ascending order.
+
+    use netcdf95, only: nf95_open, nf95_close, nf95_gw_var, nf95_inq_varid
+    use netcdf, only: nf90_nowrite
+
+    use mod_phys_lmdz_mpi_data, only: is_mpi_root
+    use mod_phys_lmdz_mpi_transfert, only: bcast_mpi ! broadcast
+
+    integer, intent(out):: ncid ! of "climoz_LMDZ.nc"
+
+    real, pointer:: press_in_edg(:)
+    ! (edges of pressure intervals for ozone climatology, in Pa, in strictly
+    ! ascending order)
+
+    ! 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 k
+
+    !---------------------------------------
+
+    print *, "Call sequence information: open_climoz"
+
+    if (is_mpi_root) then
+       call nf95_open("climoz_LMDZ.nc", nf90_nowrite, ncid)
+
+       call nf95_inq_varid(ncid, "plev", varid)
+       call nf95_gw_var(ncid, varid, plev)
+       ! Convert from hPa to Pa because "paprs" and "pplay" are in Pa:
+       plev = plev * 100.
+       n_plev = size(plev)
+    end if
+
+    call bcast_mpi(n_plev)
+    if (.not. is_mpi_root) allocate(plev(n_plev))
+    call bcast_mpi(plev)
+    
+    ! Compute edges of pressure intervals:
+    allocate(press_in_edg(n_plev + 1))
+    if (is_mpi_root) then
+       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))
+       press_in_edg(n_plev + 1) = huge(0.)
+       ! (infinity, but any value guaranteed to be greater than the
+       ! surface pressure would do)
+    end if
+    call bcast_mpi(press_in_edg)
+    deallocate(plev) ! pointer
+
+  end subroutine open_climoz
+
+end module open_climoz_m
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F	(revision 1153)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F	(revision 1154)
@@ -28,4 +28,8 @@
       USE fonte_neige_mod, ONLY  : fonte_neige_get_vars
       USE phys_output_mod
+      use open_climoz_m, only: open_climoz ! ozone climatology
+      use regr_pr, only: regr_pr_av
+      use netcdf95, only: nf95_close
+      use mod_phys_lmdz_mpi_data, only: is_mpi_root
 
       IMPLICIT none
@@ -710,4 +714,7 @@
       SAVE lmt_pas                ! frequence de mise a jour
 c$OMP THREADPRIVATE(lmt_pas) 
+      real zmasse(klon, llm) 
+C     (column-density of mass of air in a cell, in kg m-2)
+      real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
 
 cIM sorties
@@ -1093,4 +1100,12 @@
       LOGICAL,SAVE :: first=.true.
 c$OMP THREADPRIVATE(first)
+
+      logical, save::  read_climoz ! read ozone climatology
+      integer, save:: ncid_climoz ! NetCDF file containing ozone climatology
+
+      real, pointer, save:: press_climoz(:)
+!     edges of pressure intervals for ozone climatology, in Pa, in strictly
+!     ascending order
+
 #include "YOMCST.h"
 #include "YOETHF.h"
@@ -1139,4 +1154,5 @@
 
       torsfc=0.
+      forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
 
       if (first) then 
@@ -1234,5 +1250,5 @@
      .                  iflag_thermals_ed,iflag_thermals_optflux,
 cnv flags pour la convection et les poches froides
-     .                   iflag_coupl,iflag_clos,iflag_wake)
+     .                   iflag_coupl,iflag_clos,iflag_wake, read_climoz)
 
       print*,'iflag_coupl,iflag_clos,iflag_wake',
@@ -1555,4 +1571,9 @@
       call iniradia(klon,klev,paprs(1,1:klev+1))
 
+C$omp single
+      if (read_climoz) then
+         call open_climoz(ncid_climoz, press_climoz)
+      END IF
+C$omp end single
       ENDIF
 !
@@ -1723,6 +1744,19 @@
 c
       IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
-         if(prt_level.ge.1) WRITE(lunout,*)' PHYS cond  julien ',julien
-         CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
+C        Once per day, update ozone:
+         if (read_climoz) then
+C           Ozone climatology from a NetCDF file
+            call regr_pr_av(ncid_climoz, "tro3", julien, press_climoz,
+     $           paprs, wo)
+!           Convert from mole fraction of ozone to column density of ozone in a
+!           cell, in kDU:
+            wo = wo * 48. / 29. * zmasse / dobson_u / 1e3
+C           (By regridding ozone values for LMDZ only once per day, we
+C           have already neglected the variation of pressure in one
+C           day. So do not recompute "wo" at each time step even if
+C           "zmasse" changes a little.)
+         else
+            CALL ozonecm(real(julien), rlat, paprs, wo)
+         end if
       ENDIF
 c
@@ -3496,4 +3530,12 @@
 !         write(97) u_seri,v_seri,t_seri,q_seri
 !         close(97)
+C$OMP single
+         if (read_climoz) then
+            if (is_mpi_root) then
+               call nf95_close(ncid_climoz)
+            end if
+            deallocate(press_climoz) ! pointer
+         end if
+C$OMP end single nowait
       ENDIF
       
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 1154)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_lat_time_climoz_m.f90	(revision 1154)
@@ -0,0 +1,321 @@
+! $Id$
+module regr_lat_time_climoz_m
+
+  ! Author: Lionel GUEZ
+
+  implicit none
+
+  private
+  public regr_lat_time_climoz
+
+contains
+
+  subroutine regr_lat_time_climoz
+
+    ! "regr_lat_time_climoz" stands for "regrid latitude time
+    ! climatology ozone".
+
+    ! This procedure reads a climatology of ozone from a NetCDF file,
+    ! regrids it in latitude and time, and writes the regridded field
+    ! to a new NetCDF file.
+
+    ! The input field depends on time, pressure level and latitude.
+    ! We assume that the input field is a step function of latitude
+    ! and that the input latitude coordinate gives the centers of steps.
+    ! Regridding in latitude is made by averaging, with a cosine of
+    ! latitude factor.
+    ! The target LMDZ latitude grid is the "scalar" grid: "rlatu".
+    ! The values of "rlatu" are taken to be the centers of intervals.
+
+    ! Regridding in time is by linear interpolation.
+    ! Monthly values are processed to get daily values, on the basis
+    ! of a 360-day calendar.
+    ! The input file may contain either values for 12 months or values
+    ! for 14 months.
+    ! If there are 14 months then we assume that we have (in that order):
+    ! December, January, February, ..., November, December, January
+    ! We use the first December value to interpolate values between January
+    ! 1st and mid-January.
+    ! We use the last January value to interpolate values between
+    ! mid-December and end of December.
+    ! If there are only 12 months in the input file than we assume
+    ! periodicity for interpolation at the beginning and at the end of the
+    ! year.
+
+    ! We assume that in the input file:
+    ! -- the latitude is in degrees and strictly monotonic (as all
+    ! NetCDF coordinate variables should be);
+    ! -- time increases (even though we do not use values of the input
+    ! time coordinate);
+    ! -- pressure is in hPa and in strictly ascending order (even
+    ! though we do not use pressure values here, we write the unit of
+    ! pressure in the NetCDF header, and we will use the assumptions later,
+    ! when we regrid in pressure).
+
+    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_var
+
+    ! Variables local to the procedure:
+
+    include "dimensions.h"
+    ! (for "jjm")
+    include "paramet.h"
+    ! (for the other included files)
+    include "comgeom2.h"
+    ! (for "rlatv")
+    include "comconst.h"
+    ! (for "pi")
+
+    integer ncid_in, ncid_out ! NetCDF IDs for input and output files
+    integer n_plev ! number of pressure levels in the input data
+    integer n_lat! number of latitudes in the input data
+
+    real, pointer:: latitude(:)
+    ! (of input data, converted to rad, sorted in strictly ascending order)
+
+    real, allocatable:: lat_in_edg(:)
+    ! (edges of latitude intervals for input data, in rad, in strictly
+    ! ascending order)
+
+    real, pointer:: plev(:) ! pressure level of input data
+    logical desc_lat ! latitude 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, l, month)" is at latitude "latitude(j)" and pressure
+    ! level "plev(l)". "month" is between 1 and 12 or between 0 and 13)
+
+    real, allocatable:: o3_regr_lat(:, :, :) ! (jjm + 1, n_plev, 0:13)
+    ! (mean of "o3_in" over a latitude interval of LMDZ)
+    ! (First dimension is latitude interval.
+    ! The latitude interval for "o3_regr_lat(j,:, :)" contains "rlatu(j)".
+    ! If "j" is between 2 and "jjm" then the interval is:
+    ! [rlatv(j), rlatv(j-1)]
+    ! If "j" is 1 or "jjm + 1" then the interval is:
+    ! [rlatv(1), pi / 2]
+    ! or:
+    ! [- pi / 2, rlatv(jjm)]
+    ! respectively.
+    ! "o3_regr_lat(:, l, :)" is for pressure level "plev(l)".
+    ! Last dimension is month number.)
+
+    real, allocatable:: o3_out(:, :, :) ! (jjm + 1, n_plev, 360)
+    ! (regridded ozone climatology)
+    ! ("o3_out(j, l, day)" is at latitude "rlatu(j)", pressure
+    ! level "plev(l)" and date "January 1st 0h" + "tmidday(day)", in a
+    ! 360-day calendar.)
+
+    integer j
+
+    integer varid_in, varid_out, varid_plev, varid_time
+    integer varid
+    ! (for NetCDF)
+
+    real, parameter:: tmidmonth(0:13) = (/(-15. + 30. * j, j = 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) = (/(j + 0.5, j = 0, 359)/)
+    ! (time to middle of day, in days since January 1st 0h, in a
+    ! 360-day calendar)
+
+    !---------------------------------
+
+    print *, "Call sequence information: regr_lat_time_climoz"
+
+    call nf95_open("climoz.nc", nf90_nowrite, ncid_in)
+
+    ! Get coordinates from the input file:
+
+    call nf95_inq_varid(ncid_in, "latitude", varid)
+    call nf95_gw_var(ncid_in, varid, latitude)
+    ! Convert from degrees to rad, because we will take the sine of latitude:
+    latitude = latitude / 180. * pi
+    n_lat = size(latitude)
+    ! We need to supply the latitudes to "regr1_step_av" in
+    ! ascending order, so invert order if necessary:
+    desc_lat = latitude(1) > latitude(n_lat)
+    if (desc_lat) latitude = latitude(n_lat:1:-1)
+
+    ! Compute edges of latitude intervals:
+    allocate(lat_in_edg(n_lat + 1))
+    lat_in_edg(1) = - pi / 2
+    forall (j = 2:n_lat) lat_in_edg(j) = (latitude(j - 1) + latitude(j)) / 2
+    lat_in_edg(n_lat + 1) = pi / 2
+    deallocate(latitude) ! pointer
+
+    call nf95_inq_varid(ncid_in, "plev", varid)
+    call nf95_gw_var(ncid_in, varid, plev)
+    n_plev = size(plev)
+    ! (We only need the pressure coordinate to copy it to the output file.)
+
+    ! Create the output file and get the variable IDs:
+    call prepare_out(ncid_in, n_plev, ncid_out, varid_out, varid_plev, &
+         varid_time)
+
+    ! Write remaining coordinate variables:
+    call nf95_put_var(ncid_out, varid_plev, plev)
+    call nf95_put_var(ncid_out, varid_time, tmidday)
+
+    deallocate(plev) ! pointer
+
+    call nf95_inq_varid(ncid_in, "tro3", varid_in)
+    call nf95_gw_var(ncid_in, varid_in, o3_in)
+    if (desc_lat) o3_in = o3_in(n_lat:1:-1, :, :)
+
+    call nf95_close(ncid_in)
+
+    allocate(o3_regr_lat(jjm + 1, n_plev, 0:13))
+    allocate(o3_out(jjm + 1, n_plev, 360))
+
+    ! Regrid in latitude:
+    ! We average with respect to sine of latitude, which is
+    ! equivalent to weighting by cosine of latitude:
+    if (size(o3_in, 3) == 12) then
+       print *, "Found 12 months in ozone climatology, assuming periodicity..."
+       o3_regr_lat(jjm+1:1:-1, :, 1:12) = regr1_step_av(o3_in, &
+            xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(jjm:1:-1), pi / 2/)))
+       ! (invert order of indices in "o3_regr_lat" because "rlatu" is
+       ! in descending order)
+
+       ! Duplicate January and December values, in preparation of time
+       ! interpolation:
+       o3_regr_lat(:, :, 0) = o3_regr_lat(:, :, 12)
+       o3_regr_lat(:, :, 13) = o3_regr_lat(:, :, 1)
+    else
+       print *, "Using 14 months in ozone climatology..."
+       o3_regr_lat(jjm+1:1:-1, :, :) = regr1_step_av(o3_in, &
+            xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(jjm:1:-1), pi / 2/)))
+       ! (invert order of indices in "o3_regr_lat" because "rlatu" is
+       ! in descending order)
+    end if
+    
+    deallocate(o3_in) ! pointer
+
+    ! Regrid in time by linear interpolation:
+    o3_out = regr3_lint(o3_regr_lat, tmidmonth, tmidday)
+
+    ! Write to file:
+    call nf95_put_var(ncid_out, varid_out, o3_out(jjm+1:1:-1, :, :))
+    ! (The order of "rlatu" is inverted in the output file)
+
+    call nf95_close(ncid_out)
+
+  end subroutine regr_lat_time_climoz
+
+  !********************************************
+
+  subroutine prepare_out(ncid_in, n_plev, ncid_out, varid_out, varid_plev, &
+       varid_time)
+
+    ! This subroutine creates the NetCDF output file, defines
+    ! dimensions and variables, and writes one of the coordinate variables.
+
+    use netcdf95, only: nf95_create, nf95_def_dim, nf95_def_var, &
+         nf95_put_att, nf95_enddef, nf95_copy_att, nf95_put_var
+    use netcdf, only: nf90_clobber, nf90_float, nf90_copy_att, nf90_global
+
+    integer, intent(in):: ncid_in, n_plev
+    integer, intent(out):: ncid_out, varid_out, varid_plev, varid_time
+
+    ! Variables local to the procedure:
+
+    include "dimensions.h"
+    ! (for "jjm")
+    include "paramet.h"
+    ! (for the other included files)
+    include "comgeom2.h"
+    ! (for "rlatu")
+    include "comconst.h"
+    ! (for "pi")
+
+    integer ncerr
+    integer dimid_rlatu, dimid_plev, dimid_time
+    integer varid_rlatu
+
+    !---------------------------
+
+    print *, "Call sequence information: prepare_out"
+
+    call nf95_create("climoz_LMDZ.nc", nf90_clobber, ncid_out)
+
+    ! Dimensions:
+    call nf95_def_dim(ncid_out, "time", 360, dimid_time)
+    call nf95_def_dim(ncid_out, "plev", n_plev, dimid_plev)
+    call nf95_def_dim(ncid_out, "rlatu", jjm + 1, dimid_rlatu)
+
+    ! Define coordinate variables:
+
+    call nf95_def_var(ncid_out, "time", nf90_float, dimid_time, varid_time)
+    call nf95_put_att(ncid_out, varid_time, "units", "days since 2000-1-1")
+    call nf95_put_att(ncid_out, varid_time, "calendar", "360_day")
+    call nf95_put_att(ncid_out, varid_time, "standard_name", "time")
+
+    call nf95_def_var(ncid_out, "plev", nf90_float, dimid_plev, varid_plev)
+    call nf95_put_att(ncid_out, varid_plev, "units", "millibar")
+    call nf95_put_att(ncid_out, varid_plev, "standard_name", "air_pressure")
+    call nf95_put_att(ncid_out, varid_plev, "long_name", "air pressure")
+
+    call nf95_def_var(ncid_out, "rlatu", nf90_float, dimid_rlatu, varid_rlatu)
+    call nf95_put_att(ncid_out, varid_rlatu, "units", "degrees_north")
+    call nf95_put_att(ncid_out, varid_rlatu, "standard_name", "latitude")
+
+    ! Define the primary variable:
+
+    call nf95_def_var(ncid_out, "tro3", nf90_float, &
+         (/dimid_rlatu, dimid_plev, dimid_time/), varid_out)
+    call nf95_put_att(ncid_out, varid_out, "long_name", "ozone mole fraction")
+    call nf95_put_att(ncid_out, varid_out, "standard_name", &
+         "mole_fraction_of_ozone_in_air")
+
+    ! Global attributes:
+
+    call nf95_put_att(ncid_out, nf90_global, "comment", "Regridded for LMDZ")
+
+    ! 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, "source", ncid_out, nf90_global, &
+         ncerr)
+    call handle_err_copy_att("source")
+
+    call nf95_enddef(ncid_out)
+
+    ! Write one of the coordinate variables:
+    call nf95_put_var(ncid_out, varid_rlatu, rlatu(jjm+1:1:-1) / pi * 180.)
+    ! (convert from rad to degrees and sort in ascending order)
+
+  contains
+
+    subroutine handle_err_copy_att(att_name)
+
+      use netcdf, only: nf90_noerr, nf90_strerror
+
+      character(len=*), intent(in):: att_name
+
+      !----------------------------------------
+
+      if (ncerr /= nf90_noerr) then
+         print *, "regr_lat_time_climoz_m prepare_out nf90_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: LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_pr.f90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_pr.f90	(revision 1154)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_pr.f90	(revision 1154)
@@ -0,0 +1,191 @@
+! $Id$
+module regr_pr
+
+  ! 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
