Index: /LMDZ4/trunk/arch/arch-SX8_BRODIE.fcm
===================================================================
--- /LMDZ4/trunk/arch/arch-SX8_BRODIE.fcm	(revision 1424)
+++ /LMDZ4/trunk/arch/arch-SX8_BRODIE.fcm	(revision 1425)
@@ -2,5 +2,5 @@
 %LINK                sxmpif90
 %AR                  sxar
-%MAKE                make
+%MAKE                sxgmake
 %FPP_FLAGS           -P -traditional
 %FPP_DEF             NC_DOUBLE BLAS SGEMV=DGEMV SGEMM=DGEMM FFT_MATHKEISAN
Index: /LMDZ4/trunk/arch/arch-g95.fcm
===================================================================
--- /LMDZ4/trunk/arch/arch-g95.fcm	(revision 1425)
+++ /LMDZ4/trunk/arch/arch-g95.fcm	(revision 1425)
@@ -0,0 +1,15 @@
+%COMPILER            g95
+%LINK                g95
+%AR                  ar
+%MAKE                make
+%FPP_FLAGS           -P -traditional
+%FPP_DEF             PLOUF
+%BASE_FFLAGS         -c
+%PROD_FFLAGS         -O3
+%DEV_FFLAGS          -O
+%DEBUG_FFLAGS        -fbounds-check -freal=nan -ftrace=full -g -O0 -Wall
+%MPI_FFLAGS
+%OMP_FFLAGS          
+%BASE_LD             
+%MPI_LD
+%OMP_LD              
Index: /LMDZ4/trunk/arch/arch-g95.path
===================================================================
--- /LMDZ4/trunk/arch/arch-g95.path	(revision 1425)
+++ /LMDZ4/trunk/arch/arch-g95.path	(revision 1425)
@@ -0,0 +1,6 @@
+NETCDF_LIBDIR=...
+NETCDF_INCDIR=...
+IOIPSL_INCDIR=...
+IOIPSL_LIBDIR=...
+ORCH_INCDIR=...
+ORCH_LIBDIR=...
Index: /LMDZ4/trunk/arch/arch-linux-32bit.fcm
===================================================================
--- /LMDZ4/trunk/arch/arch-linux-32bit.fcm	(revision 1424)
+++ /LMDZ4/trunk/arch/arch-linux-32bit.fcm	(revision 1425)
@@ -1,4 +1,4 @@
-%COMPILER            pgf90
-%LINK                pgf90
+%COMPILER            pgf95
+%LINK                pgf95
 %AR                  ar
 %MAKE                make
@@ -7,9 +7,9 @@
 %BASE_FFLAGS         
 %PROD_FFLAGS         -fast
-%DEV_FFLAGS          -g
-%DEBUG_FFLAGS        -g
+%DEV_FFLAGS          -g -O1
+%DEBUG_FFLAGS        -g -O0 -Kieee -Ktrap=fp -Mbounds
 %MPI_FFLAGS
 %OMP_FFLAGS          
-%BASE_LD             -Wl,-Bstatic -L/usr/lib/gcc-lib/i386-linux/2.95.2 
+%BASE_LD             
 %MPI_LD
 %OMP_LD              
Index: /LMDZ4/trunk/bld.cfg
===================================================================
--- /LMDZ4/trunk/bld.cfg	(revision 1424)
+++ /LMDZ4/trunk/bld.cfg	(revision 1425)
@@ -103,2 +103,3 @@
 # extension for module output
 bld::outfile_ext::mod .mod
+bld::tool::SHELL   /bin/bash
Index: /LMDZ4/trunk/libf/bibio/arth.F90
===================================================================
--- /LMDZ4/trunk/libf/bibio/arth.F90	(revision 1425)
+++ /LMDZ4/trunk/libf/bibio/arth.F90	(revision 1425)
@@ -0,0 +1,82 @@
+MODULE arth_m
+
+  IMPLICIT NONE
+
+  INTEGER, PARAMETER, private:: NPAR_ARTH=16, NPAR2_ARTH=8
+
+  INTERFACE arth
+     ! Returns an arithmetic progression, given a first term "first", an
+     ! increment and a number of terms "n".
+
+     MODULE PROCEDURE arth_r, arth_i
+     ! The difference between the procedures is the type of
+     ! arguments "first" and "increment" and of function result.
+  END INTERFACE
+
+  private arth_r, arth_i
+
+CONTAINS
+
+  pure FUNCTION arth_r(first,increment,n)
+
+    REAL, INTENT(IN) :: first,increment
+    INTEGER, INTENT(IN) :: n
+    REAL, DIMENSION(n) :: arth_r
+
+    ! Variables local to the procedure:
+
+    INTEGER :: k,k2
+    REAL :: temp
+
+    !---------------------------------------
+
+    if (n > 0) arth_r(1)=first
+    if (n <= NPAR_ARTH) then
+       do k=2,n
+          arth_r(k)=arth_r(k-1)+increment
+       end do
+    else
+       do k=2,NPAR2_ARTH
+          arth_r(k)=arth_r(k-1)+increment
+       end do
+       temp=increment*NPAR2_ARTH
+       k=NPAR2_ARTH
+       do
+          if (k >= n) exit
+          k2=k+k
+          arth_r(k+1:min(k2,n)) = temp + arth_r(1:min(k,n-k))
+          temp=temp+temp
+          k=k2
+       end do
+    end if
+  END FUNCTION arth_r
+
+  !*************************************
+
+  pure FUNCTION arth_i(first,increment,n)
+
+    INTEGER, INTENT(IN) :: first,increment,n
+    INTEGER, DIMENSION(n) :: arth_i
+    INTEGER :: k,k2,temp
+    if (n > 0) arth_i(1)=first
+    if (n <= NPAR_ARTH) then
+       do k=2,n
+          arth_i(k)=arth_i(k-1)+increment
+       end do
+    else
+       do k=2,NPAR2_ARTH
+          arth_i(k)=arth_i(k-1)+increment
+       end do
+       temp=increment*NPAR2_ARTH
+       k=NPAR2_ARTH
+       do
+          if (k >= n) exit
+          k2=k+k
+          arth_i(k+1:min(k2,n))=temp+arth_i(1:min(k,n-k))
+          temp=temp+temp
+          k=k2
+       end do
+    end if
+  END FUNCTION arth_i
+
+END MODULE arth_m
Index: /LMDZ4/trunk/libf/bibio/chfev.F
===================================================================
--- /LMDZ4/trunk/libf/bibio/chfev.F	(revision 1425)
+++ /LMDZ4/trunk/libf/bibio/chfev.F	(revision 1425)
@@ -0,0 +1,155 @@
+*DECK CHFEV
+      SUBROUTINE CHFEV (X1, X2, F1, F2, D1, D2, NE, XE, FE, NEXT, IERR)
+C***BEGIN PROLOGUE  CHFEV
+C***PURPOSE  Evaluate a cubic polynomial given in Hermite form at an
+C            array of points.  While designed for use by PCHFE, it may
+C            be useful directly as an evaluator for a piecewise cubic
+C            Hermite function in applications, such as graphing, where
+C            the interval is known in advance.
+C***LIBRARY   SLATEC (PCHIP)
+C***CATEGORY  E3
+C***TYPE      SINGLE PRECISION (CHFEV-S, DCHFEV-D)
+C***KEYWORDS  CUBIC HERMITE EVALUATION, CUBIC POLYNOMIAL EVALUATION,
+C             PCHIP
+C***AUTHOR  Fritsch, F. N., (LLNL)
+C             Lawrence Livermore National Laboratory
+C             P.O. Box 808  (L-316)
+C             Livermore, CA  94550
+C             FTS 532-4275, (510) 422-4275
+C***DESCRIPTION
+C
+C          CHFEV:  Cubic Hermite Function EValuator
+C
+C     Evaluates the cubic polynomial determined by function values
+C     F1,F2 and derivatives D1,D2 on interval (X1,X2) at the points
+C     XE(J), J=1(1)NE.
+C
+C ----------------------------------------------------------------------
+C
+C  Calling sequence:
+C
+C        INTEGER  NE, NEXT(2), IERR
+C        REAL  X1, X2, F1, F2, D1, D2, XE(NE), FE(NE)
+C
+C        CALL  CHFEV (X1,X2, F1,F2, D1,D2, NE, XE, FE, NEXT, IERR)
+C
+C   Parameters:
+C
+C     X1,X2 -- (input) endpoints of interval of definition of cubic.
+C           (Error return if  X1.EQ.X2 .)
+C
+C     F1,F2 -- (input) values of function at X1 and X2, respectively.
+C
+C     D1,D2 -- (input) values of derivative at X1 and X2, respectively.
+C
+C     NE -- (input) number of evaluation points.  (Error return if
+C           NE.LT.1 .)
+C
+C     XE -- (input) real array of points at which the function is to be
+C           evaluated.  If any of the XE are outside the interval
+C           [X1,X2], a warning error is returned in NEXT.
+C
+C     FE -- (output) real array of values of the cubic function defined
+C           by  X1,X2, F1,F2, D1,D2  at the points  XE.
+C
+C     NEXT -- (output) integer array indicating number of extrapolation
+C           points:
+C            NEXT(1) = number of evaluation points to left of interval.
+C            NEXT(2) = number of evaluation points to right of interval.
+C
+C     IERR -- (output) error flag.
+C           Normal return:
+C              IERR = 0  (no errors).
+C           "Recoverable" errors:
+C              IERR = -1  if NE.LT.1 .
+C              IERR = -2  if X1.EQ.X2 .
+C                (The FE-array has not been changed in either case.)
+C
+C***REFERENCES  (NONE)
+C***ROUTINES CALLED  XERMSG
+C***REVISION HISTORY  (YYMMDD)
+C   811019  DATE WRITTEN
+C   820803  Minor cosmetic changes for release 1.
+C   890411  Added SAVE statements (Vers. 3.2).
+C   890531  Changed all specific intrinsics to generic.  (WRB)
+C   890703  Corrected category record.  (WRB)
+C   890703  REVISION DATE from Version 3.2
+C   891214  Prologue converted to Version 4.0 format.  (BAB)
+C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
+C***END PROLOGUE  CHFEV
+C  Programming notes:
+C
+C     To produce a double precision version, simply:
+C        a. Change CHFEV to DCHFEV wherever it occurs,
+C        b. Change the real declaration to double precision, and
+C        c. Change the constant ZERO to double precision.
+C
+C  DECLARE ARGUMENTS.
+C
+      INTEGER  NE, NEXT(2), IERR
+      REAL  X1, X2, F1, F2, D1, D2, XE(*), FE(*)
+C
+C  DECLARE LOCAL VARIABLES.
+C
+      INTEGER  I
+      REAL  C2, C3, DEL1, DEL2, DELTA, H, X, XMI, XMA, ZERO
+      SAVE ZERO
+      DATA  ZERO /0./
+C
+C  VALIDITY-CHECK ARGUMENTS.
+C
+C***FIRST EXECUTABLE STATEMENT  CHFEV
+      IF (NE .LT. 1)  GO TO 5001
+      H = X2 - X1
+      IF (H .EQ. ZERO)  GO TO 5002
+C
+C  INITIALIZE.
+C
+      IERR = 0
+      NEXT(1) = 0
+      NEXT(2) = 0
+      XMI = MIN(ZERO, H)
+      XMA = MAX(ZERO, H)
+C
+C  COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1).
+C
+      DELTA = (F2 - F1)/H
+      DEL1 = (D1 - DELTA)/H
+      DEL2 = (D2 - DELTA)/H
+C                                           (DELTA IS NO LONGER NEEDED.)
+      C2 = -(DEL1+DEL1 + DEL2)
+      C3 = (DEL1 + DEL2)/H
+C                               (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.)
+C
+C  EVALUATION LOOP.
+C
+      DO 500  I = 1, NE
+         X = XE(I) - X1
+         FE(I) = F1 + X*(D1 + X*(C2 + X*C3))
+C          COUNT EXTRAPOLATION POINTS.
+         IF ( X.LT.XMI )  NEXT(1) = NEXT(1) + 1
+         IF ( X.GT.XMA )  NEXT(2) = NEXT(2) + 1
+C        (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.)
+  500 CONTINUE
+C
+C  NORMAL RETURN.
+C
+      RETURN
+C
+C  ERROR RETURNS.
+C
+ 5001 CONTINUE
+C     NE.LT.1 RETURN.
+      IERR = -1
+      CALL XERMSG ('SLATEC', 'CHFEV',
+     +   'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1)
+      RETURN
+C
+ 5002 CONTINUE
+C     X1.EQ.X2 RETURN.
+      IERR = -2
+      CALL XERMSG ('SLATEC', 'CHFEV', 'INTERVAL ENDPOINTS EQUAL', IERR,
+     +   1)
+      RETURN
+C------------- LAST LINE OF CHFEV FOLLOWS ------------------------------
+      END
Index: /LMDZ4/trunk/libf/bibio/fdump.F
===================================================================
--- /LMDZ4/trunk/libf/bibio/fdump.F	(revision 1425)
+++ /LMDZ4/trunk/libf/bibio/fdump.F	(revision 1425)
@@ -0,0 +1,31 @@
+*DECK FDUMP
+      SUBROUTINE FDUMP
+C***BEGIN PROLOGUE  FDUMP
+C***PURPOSE  Symbolic dump (should be locally written).
+C***LIBRARY   SLATEC (XERROR)
+C***CATEGORY  R3
+C***TYPE      ALL (FDUMP-A)
+C***KEYWORDS  ERROR, XERMSG
+C***AUTHOR  Jones, R. E., (SNLA)
+C***DESCRIPTION
+C
+C        ***Note*** Machine Dependent Routine
+C        FDUMP is intended to be replaced by a locally written
+C        version which produces a symbolic dump.  Failing this,
+C        it should be replaced by a version which prints the
+C        subprogram nesting list.  Note that this dump must be
+C        printed on each of up to five files, as indicated by the
+C        XGETUA routine.  See XSETUA and XGETUA for details.
+C
+C     Written by Ron Jones, with SLATEC Common Math Library Subcommittee
+C
+C***REFERENCES  (NONE)
+C***ROUTINES CALLED  (NONE)
+C***REVISION HISTORY  (YYMMDD)
+C   790801  DATE WRITTEN
+C   861211  REVISION DATE from Version 3.2
+C   891214  Prologue converted to Version 4.0 format.  (BAB)
+C***END PROLOGUE  FDUMP
+C***FIRST EXECUTABLE STATEMENT  FDUMP
+      RETURN
+      END
Index: /LMDZ4/trunk/libf/bibio/i1mach.F
===================================================================
--- /LMDZ4/trunk/libf/bibio/i1mach.F	(revision 1425)
+++ /LMDZ4/trunk/libf/bibio/i1mach.F	(revision 1425)
@@ -0,0 +1,126 @@
+*DECK I1MACH
+      INTEGER FUNCTION I1MACH (I)
+C***BEGIN PROLOGUE  I1MACH
+C***PURPOSE  Return integer machine dependent constants.
+C***LIBRARY   SLATEC
+C***CATEGORY  R1
+C***TYPE      INTEGER (I1MACH-I)
+C***KEYWORDS  MACHINE CONSTANTS
+C***AUTHOR  Fox, P. A., (Bell Labs)
+C           Hall, A. D., (Bell Labs)
+C           Schryer, N. L., (Bell Labs)
+C***DESCRIPTION
+C
+C   I1MACH can be used to obtain machine-dependent parameters for the
+C   local machine environment.  It is a function subprogram with one
+C   (input) argument and can be referenced as follows:
+C
+C        K = I1MACH(I)
+C
+C   where I=1,...,16.  The (output) value of K above is determined by
+C   the (input) value of I.  The results for various values of I are
+C   discussed below.
+C
+C   I/O unit numbers:
+C     I1MACH( 1) = the standard input unit.
+C     I1MACH( 2) = the standard output unit.
+C     I1MACH( 3) = the standard punch unit.
+C     I1MACH( 4) = the standard error message unit.
+C
+C   Words:
+C     I1MACH( 5) = the number of bits per integer storage unit.
+C     I1MACH( 6) = the number of characters per integer storage unit.
+C
+C   Integers:
+C     assume integers are represented in the S-digit, base-A form
+C
+C                sign ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) )
+C
+C                where 0 .LE. X(I) .LT. A for I=0,...,S-1.
+C     I1MACH( 7) = A, the base.
+C     I1MACH( 8) = S, the number of base-A digits.
+C     I1MACH( 9) = A**S - 1, the largest magnitude.
+C
+C   Floating-Point Numbers:
+C     Assume floating-point numbers are represented in the T-digit,
+C     base-B form
+C                sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
+C
+C                where 0 .LE. X(I) .LT. B for I=1,...,T,
+C                0 .LT. X(1), and EMIN .LE. E .LE. EMAX.
+C     I1MACH(10) = B, the base.
+C
+C   Single-Precision:
+C     I1MACH(11) = T, the number of base-B digits.
+C     I1MACH(12) = EMIN, the smallest exponent E.
+C     I1MACH(13) = EMAX, the largest exponent E.
+C
+C   Double-Precision:
+C     I1MACH(14) = T, the number of base-B digits.
+C     I1MACH(15) = EMIN, the smallest exponent E.
+C     I1MACH(16) = EMAX, the largest exponent E.
+C
+C   To alter this function for a particular environment, the desired
+C   set of DATA statements should be activated by removing the C from
+C   column 1.  Also, the values of I1MACH(1) - I1MACH(4) should be
+C   checked for consistency with the local operating system.
+C
+C***REFERENCES  P. A. Fox, A. D. Hall and N. L. Schryer, Framework for
+C                 a portable library, ACM Transactions on Mathematical
+C                 Software 4, 2 (June 1978), pp. 177-188.
+C***ROUTINES CALLED  (NONE)
+C***REVISION HISTORY  (YYMMDD)
+C   750101  DATE WRITTEN
+C   891012  Added VAX G-floating constants.  (WRB)
+C   891012  REVISION DATE from Version 3.2
+C   891214  Prologue converted to Version 4.0 format.  (BAB)
+C   900618  Added DEC RISC constants.  (WRB)
+C   900723  Added IBM RS 6000 constants.  (WRB)
+C   901009  Correct I1MACH(7) for IBM Mainframes. Should be 2 not 16.
+C           (RWC)
+C   910710  Added HP 730 constants.  (SMR)
+C   911114  Added Convex IEEE constants.  (WRB)
+C   920121  Added SUN -r8 compiler option constants.  (WRB)
+C   920229  Added Touchstone Delta i860 constants.  (WRB)
+C   920501  Reformatted the REFERENCES section.  (WRB)
+C   920625  Added Convex -p8 and -pd8 compiler option constants.
+C           (BKS, WRB)
+C   930201  Added DEC Alpha and SGI constants.  (RWC and WRB)
+C   930618  Corrected I1MACH(5) for Convex -p8 and -pd8 compiler
+C           options.  (DWL, RWC and WRB).
+C   100623  Use Fortran 95 intrinsic functions (Lionel GUEZ)
+C***END PROLOGUE  I1MACH
+C
+      INTEGER IMACH(16),OUTPUT
+      SAVE IMACH
+      EQUIVALENCE (IMACH(4),OUTPUT)
+C***FIRST EXECUTABLE STATEMENT  I1MACH
+      IMACH( 1) =         5
+      IMACH( 2) =         6
+      IMACH( 3) =         6
+      IMACH( 4) =         6
+      IMACH( 5) =        bit_size(0)
+      IMACH( 6) =         IMACH( 5) / 8
+      IMACH( 7) =         radix(0)
+      IMACH( 8) =        digits(0)
+      IMACH( 9) =     huge(0)
+      IMACH(10) =         radix(0.)
+      IMACH(11) =        digits(0.)
+      IMACH(12) =      minexponent(0.)
+      IMACH(13) =       maxexponent(0.)
+      IMACH(14) =        digits(0d0)
+      IMACH(15) =      minexponent(0d0)
+      IMACH(16) =       maxexponent(0d0)
+      IF (I .LT. 1  .OR.  I .GT. 16) GO TO 10
+C
+      I1MACH = IMACH(I)
+      RETURN
+C
+   10 CONTINUE
+      WRITE (UNIT = OUTPUT, FMT = 9000)
+ 9000 FORMAT ('1ERROR    1 IN I1MACH - I OUT OF BOUNDS')
+C
+C     CALL FDUMP
+C
+      STOP
+      END
Index: /LMDZ4/trunk/libf/bibio/j4save.F
===================================================================
--- /LMDZ4/trunk/libf/bibio/j4save.F	(revision 1425)
+++ /LMDZ4/trunk/libf/bibio/j4save.F	(revision 1425)
@@ -0,0 +1,65 @@
+*DECK J4SAVE
+      FUNCTION J4SAVE (IWHICH, IVALUE, ISET)
+C***BEGIN PROLOGUE  J4SAVE
+C***SUBSIDIARY
+C***PURPOSE  Save or recall global variables needed by error
+C            handling routines.
+C***LIBRARY   SLATEC (XERROR)
+C***TYPE      INTEGER (J4SAVE-I)
+C***KEYWORDS  ERROR MESSAGES, ERROR NUMBER, RECALL, SAVE, XERROR
+C***AUTHOR  Jones, R. E., (SNLA)
+C***DESCRIPTION
+C
+C     Abstract
+C        J4SAVE saves and recalls several global variables needed
+C        by the library error handling routines.
+C
+C     Description of Parameters
+C      --Input--
+C        IWHICH - Index of item desired.
+C                = 1 Refers to current error number.
+C                = 2 Refers to current error control flag.
+C                = 3 Refers to current unit number to which error
+C                    messages are to be sent.  (0 means use standard.)
+C                = 4 Refers to the maximum number of times any
+C                     message is to be printed (as set by XERMAX).
+C                = 5 Refers to the total number of units to which
+C                     each error message is to be written.
+C                = 6 Refers to the 2nd unit for error messages
+C                = 7 Refers to the 3rd unit for error messages
+C                = 8 Refers to the 4th unit for error messages
+C                = 9 Refers to the 5th unit for error messages
+C        IVALUE - The value to be set for the IWHICH-th parameter,
+C                 if ISET is .TRUE. .
+C        ISET   - If ISET=.TRUE., the IWHICH-th parameter will BE
+C                 given the value, IVALUE.  If ISET=.FALSE., the
+C                 IWHICH-th parameter will be unchanged, and IVALUE
+C                 is a dummy parameter.
+C      --Output--
+C        The (old) value of the IWHICH-th parameter will be returned
+C        in the function value, J4SAVE.
+C
+C***SEE ALSO  XERMSG
+C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+C                 Error-handling Package, SAND82-0800, Sandia
+C                 Laboratories, 1982.
+C***ROUTINES CALLED  (NONE)
+C***REVISION HISTORY  (YYMMDD)
+C   790801  DATE WRITTEN
+C   891214  Prologue converted to Version 4.0 format.  (BAB)
+C   900205  Minor modifications to prologue.  (WRB)
+C   900402  Added TYPE section.  (WRB)
+C   910411  Added KEYWORDS section.  (WRB)
+C   920501  Reformatted the REFERENCES section.  (WRB)
+C***END PROLOGUE  J4SAVE
+      LOGICAL ISET
+      INTEGER IPARAM(9)
+      SAVE IPARAM
+      DATA IPARAM(1),IPARAM(2),IPARAM(3),IPARAM(4)/0,2,0,10/
+      DATA IPARAM(5)/1/
+      DATA IPARAM(6),IPARAM(7),IPARAM(8),IPARAM(9)/0,0,0,0/
+C***FIRST EXECUTABLE STATEMENT  J4SAVE
+      J4SAVE = IPARAM(IWHICH)
+      IF (ISET) IPARAM(IWHICH) = IVALUE
+      RETURN
+      END
Index: /LMDZ4/trunk/libf/bibio/pchdf.F
===================================================================
--- /LMDZ4/trunk/libf/bibio/pchdf.F	(revision 1425)
+++ /LMDZ4/trunk/libf/bibio/pchdf.F	(revision 1425)
@@ -0,0 +1,106 @@
+*DECK PCHDF
+      REAL FUNCTION PCHDF (K, X, S, IERR)
+C***BEGIN PROLOGUE  PCHDF
+C***SUBSIDIARY
+C***PURPOSE  Computes divided differences for PCHCE and PCHSP
+C***LIBRARY   SLATEC (PCHIP)
+C***TYPE      SINGLE PRECISION (PCHDF-S, DPCHDF-D)
+C***AUTHOR  Fritsch, F. N., (LLNL)
+C***DESCRIPTION
+C
+C          PCHDF:   PCHIP Finite Difference Formula
+C
+C     Uses a divided difference formulation to compute a K-point approx-
+C     imation to the derivative at X(K) based on the data in X and S.
+C
+C     Called by  PCHCE  and  PCHSP  to compute 3- and 4-point boundary
+C     derivative approximations.
+C
+C ----------------------------------------------------------------------
+C
+C     On input:
+C        K      is the order of the desired derivative approximation.
+C               K must be at least 3 (error return if not).
+C        X      contains the K values of the independent variable.
+C               X need not be ordered, but the values **MUST** be
+C               distinct.  (Not checked here.)
+C        S      contains the associated slope values:
+C                  S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1.
+C               (Note that S need only be of length K-1.)
+C
+C     On return:
+C        S      will be destroyed.
+C        IERR   will be set to -1 if K.LT.2 .
+C        PCHDF  will be set to the desired derivative approximation if
+C               IERR=0 or to zero if IERR=-1.
+C
+C ----------------------------------------------------------------------
+C
+C***SEE ALSO  PCHCE, PCHSP
+C***REFERENCES  Carl de Boor, A Practical Guide to Splines, Springer-
+C                 Verlag, New York, 1978, pp. 10-16.
+C***ROUTINES CALLED  XERMSG
+C***REVISION HISTORY  (YYMMDD)
+C   820503  DATE WRITTEN
+C   820805  Converted to SLATEC library version.
+C   870813  Minor cosmetic changes.
+C   890411  Added SAVE statements (Vers. 3.2).
+C   890411  REVISION DATE from Version 3.2
+C   891214  Prologue converted to Version 4.0 format.  (BAB)
+C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
+C   900328  Added TYPE section.  (WRB)
+C   910408  Updated AUTHOR and DATE WRITTEN sections in prologue.  (WRB)
+C   920429  Revised format and order of references.  (WRB,FNF)
+C   930503  Improved purpose.  (FNF)
+C***END PROLOGUE  PCHDF
+C
+C**End
+C
+C  DECLARE ARGUMENTS.
+C
+      INTEGER  K, IERR
+      REAL  X(K), S(K)
+C
+C  DECLARE LOCAL VARIABLES.
+C
+      INTEGER  I, J
+      REAL  VALUE, ZERO
+      SAVE ZERO
+      DATA  ZERO /0./
+C
+C  CHECK FOR LEGAL VALUE OF K.
+C
+C***FIRST EXECUTABLE STATEMENT  PCHDF
+      IF (K .LT. 3)  GO TO 5001
+C
+C  COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL.
+C
+      DO 10  J = 2, K-1
+         DO 9  I = 1, K-J
+            S(I) = (S(I+1)-S(I))/(X(I+J)-X(I))
+    9    CONTINUE
+   10 CONTINUE
+C
+C  EVALUATE DERIVATIVE AT X(K).
+C
+      VALUE = S(1)
+      DO 20  I = 2, K-1
+         VALUE = S(I) + VALUE*(X(K)-X(I))
+   20 CONTINUE
+C
+C  NORMAL RETURN.
+C
+      IERR = 0
+      PCHDF = VALUE
+      RETURN
+C
+C  ERROR RETURN.
+C
+ 5001 CONTINUE
+C     K.LT.3 RETURN.
+      IERR = -1
+      CALL XERMSG ('SLATEC', 'PCHDF', 'K LESS THAN THREE', IERR, 1)
+      PCHDF = ZERO
+      RETURN
+C------------- LAST LINE OF PCHDF FOLLOWS ------------------------------
+      END
Index: /LMDZ4/trunk/libf/bibio/pchfe.F
===================================================================
--- /LMDZ4/trunk/libf/bibio/pchfe.F	(revision 1425)
+++ /LMDZ4/trunk/libf/bibio/pchfe.F	(revision 1425)
@@ -0,0 +1,308 @@
+*DECK PCHFE
+      SUBROUTINE PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
+C***BEGIN PROLOGUE  PCHFE
+C***PURPOSE  Evaluate a piecewise cubic Hermite function at an array of
+C            points.  May be used by itself for Hermite interpolation,
+C            or as an evaluator for PCHIM or PCHIC.
+C***LIBRARY   SLATEC (PCHIP)
+C***CATEGORY  E3
+C***TYPE      SINGLE PRECISION (PCHFE-S, DPCHFE-D)
+C***KEYWORDS  CUBIC HERMITE EVALUATION, HERMITE INTERPOLATION, PCHIP,
+C             PIECEWISE CUBIC EVALUATION
+C***AUTHOR  Fritsch, F. N., (LLNL)
+C             Lawrence Livermore National Laboratory
+C             P.O. Box 808  (L-316)
+C             Livermore, CA  94550
+C             FTS 532-4275, (510) 422-4275
+C***DESCRIPTION
+C
+C          PCHFE:  Piecewise Cubic Hermite Function Evaluator
+C
+C     Evaluates the cubic Hermite function defined by  N, X, F, D  at
+C     the points  XE(J), J=1(1)NE.
+C
+C     To provide compatibility with PCHIM and PCHIC, includes an
+C     increment between successive values of the F- and D-arrays.
+C
+C ----------------------------------------------------------------------
+C
+C  Calling sequence:
+C
+C        PARAMETER  (INCFD = ...)
+C        INTEGER  N, NE, IERR
+C        REAL  X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE)
+C        LOGICAL  SKIP
+C
+C        CALL  PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
+C
+C   Parameters:
+C
+C     N -- (input) number of data points.  (Error return if N.LT.2 .)
+C
+C     X -- (input) real array of independent variable values.  The
+C           elements of X must be strictly increasing:
+C                X(I-1) .LT. X(I),  I = 2(1)N.
+C           (Error return if not.)
+C
+C     F -- (input) real array of function values.  F(1+(I-1)*INCFD) is
+C           the value corresponding to X(I).
+C
+C     D -- (input) real array of derivative values.  D(1+(I-1)*INCFD) is
+C           the value corresponding to X(I).
+C
+C     INCFD -- (input) increment between successive values in F and D.
+C           (Error return if  INCFD.LT.1 .)
+C
+C     SKIP -- (input/output) logical variable which should be set to
+C           .TRUE. if the user wishes to skip checks for validity of
+C           preceding parameters, or to .FALSE. otherwise.
+C           This will save time in case these checks have already
+C           been performed (say, in PCHIM or PCHIC).
+C           SKIP will be set to .TRUE. on normal return.
+C
+C     NE -- (input) number of evaluation points.  (Error return if
+C           NE.LT.1 .)
+C
+C     XE -- (input) real array of points at which the function is to be
+C           evaluated.
+C
+C          NOTES:
+C           1. The evaluation will be most efficient if the elements
+C              of XE are increasing relative to X;
+C              that is,   XE(J) .GE. X(I)
+C              implies    XE(K) .GE. X(I),  all K.GE.J .
+C           2. If any of the XE are outside the interval [X(1),X(N)],
+C              values are extrapolated from the nearest extreme cubic,
+C              and a warning error is returned.
+C
+C     FE -- (output) real array of values of the cubic Hermite function
+C           defined by  N, X, F, D  at the points  XE.
+C
+C     IERR -- (output) error flag.
+C           Normal return:
+C              IERR = 0  (no errors).
+C           Warning error:
+C              IERR.GT.0  means that extrapolation was performed at
+C                 IERR points.
+C           "Recoverable" errors:
+C              IERR = -1  if N.LT.2 .
+C              IERR = -2  if INCFD.LT.1 .
+C              IERR = -3  if the X-array is not strictly increasing.
+C              IERR = -4  if NE.LT.1 .
+C             (The FE-array has not been changed in any of these cases.)
+C               NOTE:  The above errors are checked in the order listed,
+C                   and following arguments have **NOT** been validated.
+C
+C***REFERENCES  (NONE)
+C***ROUTINES CALLED  CHFEV, XERMSG
+C***REVISION HISTORY  (YYMMDD)
+C   811020  DATE WRITTEN
+C   820803  Minor cosmetic changes for release 1.
+C   870707  Minor cosmetic changes to prologue.
+C   890531  Changed all specific intrinsics to generic.  (WRB)
+C   890831  Modified array declarations.  (WRB)
+C   890831  REVISION DATE from Version 3.2
+C   891214  Prologue converted to Version 4.0 format.  (BAB)
+C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
+C***END PROLOGUE  PCHFE
+C  Programming notes:
+C
+C     1. To produce a double precision version, simply:
+C        a. Change PCHFE to DPCHFE, and CHFEV to DCHFEV, wherever they
+C           occur,
+C        b. Change the real declaration to double precision,
+C
+C     2. Most of the coding between the call to CHFEV and the end of
+C        the IR-loop could be eliminated if it were permissible to
+C        assume that XE is ordered relative to X.
+C
+C     3. CHFEV does not assume that X1 is less than X2.  thus, it would
+C        be possible to write a version of PCHFE that assumes a strict-
+C        ly decreasing X-array by simply running the IR-loop backwards
+C        (and reversing the order of appropriate tests).
+C
+C     4. The present code has a minor bug, which I have decided is not
+C        worth the effort that would be required to fix it.
+C        If XE contains points in [X(N-1),X(N)], followed by points .LT.
+C        X(N-1), followed by points .GT.X(N), the extrapolation points
+C        will be counted (at least) twice in the total returned in IERR.
+C
+C  DECLARE ARGUMENTS.
+C
+      INTEGER  N, INCFD, NE, IERR
+      REAL  X(*), F(INCFD,*), D(INCFD,*), XE(*), FE(*)
+      LOGICAL  SKIP
+C
+C  DECLARE LOCAL VARIABLES.
+C
+      INTEGER  I, IERC, IR, J, JFIRST, NEXT(2), NJ
+C
+C  VALIDITY-CHECK ARGUMENTS.
+C
+C***FIRST EXECUTABLE STATEMENT  PCHFE
+      IF (SKIP)  GO TO 5
+C
+      IF ( N.LT.2 )  GO TO 5001
+      IF ( INCFD.LT.1 )  GO TO 5002
+      DO 1  I = 2, N
+         IF ( X(I).LE.X(I-1) )  GO TO 5003
+    1 CONTINUE
+C
+C  FUNCTION DEFINITION IS OK, GO ON.
+C
+    5 CONTINUE
+      IF ( NE.LT.1 )  GO TO 5004
+      IERR = 0
+      SKIP = .TRUE.
+C
+C  LOOP OVER INTERVALS.        (   INTERVAL INDEX IS  IL = IR-1  . )
+C                              ( INTERVAL IS X(IL).LE.X.LT.X(IR) . )
+      JFIRST = 1
+      IR = 2
+   10 CONTINUE
+C
+C     SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS.
+C
+         IF (JFIRST .GT. NE)  GO TO 5000
+C
+C     LOCATE ALL POINTS IN INTERVAL.
+C
+         DO 20  J = JFIRST, NE
+            IF (XE(J) .GE. X(IR))  GO TO 30
+   20    CONTINUE
+         J = NE + 1
+         GO TO 40
+C
+C     HAVE LOCATED FIRST POINT BEYOND INTERVAL.
+C
+   30    CONTINUE
+         IF (IR .EQ. N)  J = NE + 1
+C
+   40    CONTINUE
+         NJ = J - JFIRST
+C
+C     SKIP EVALUATION IF NO POINTS IN INTERVAL.
+C
+         IF (NJ .EQ. 0)  GO TO 50
+C
+C     EVALUATE CUBIC AT XE(I),  I = JFIRST (1) J-1 .
+C
+C       ----------------------------------------------------------------
+        CALL CHFEV (X(IR-1),X(IR), F(1,IR-1),F(1,IR), D(1,IR-1),D(1,IR),
+     *              NJ, XE(JFIRST), FE(JFIRST), NEXT, IERC)
+C       ----------------------------------------------------------------
+         IF (IERC .LT. 0)  GO TO 5005
+C
+         IF (NEXT(2) .EQ. 0)  GO TO 42
+C        IF (NEXT(2) .GT. 0)  THEN
+C           IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE
+C           RIGHT OF X(IR).
+C
+            IF (IR .LT. N)  GO TO 41
+C           IF (IR .EQ. N)  THEN
+C              THESE ARE ACTUALLY EXTRAPOLATION POINTS.
+               IERR = IERR + NEXT(2)
+               GO TO 42
+   41       CONTINUE
+C           ELSE
+C              WE SHOULD NEVER HAVE GOTTEN HERE.
+               GO TO 5005
+C           ENDIF
+C        ENDIF
+   42    CONTINUE
+C
+         IF (NEXT(1) .EQ. 0)  GO TO 49
+C        IF (NEXT(1) .GT. 0)  THEN
+C           IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(1) TO THE
+C           LEFT OF X(IR-1).
+C
+            IF (IR .GT. 2)  GO TO 43
+C           IF (IR .EQ. 2)  THEN
+C              THESE ARE ACTUALLY EXTRAPOLATION POINTS.
+               IERR = IERR + NEXT(1)
+               GO TO 49
+   43       CONTINUE
+C           ELSE
+C              XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST
+C              EVALUATION INTERVAL.
+C
+C              FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1).
+               DO 44  I = JFIRST, J-1
+                  IF (XE(I) .LT. X(IR-1))  GO TO 45
+   44          CONTINUE
+C              NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR
+C                     IN CHFEV.
+               GO TO 5005
+C
+   45          CONTINUE
+C              RESET J.  (THIS WILL BE THE NEW JFIRST.)
+               J = I
+C
+C              NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY.
+               DO 46  I = 1, IR-1
+                  IF (XE(J) .LT. X(I)) GO TO 47
+   46          CONTINUE
+C              NB-- CAN NEVER DROP THROUGH HERE, SINCE XE(J).LT.X(IR-1).
+C
+   47          CONTINUE
+C              AT THIS POINT, EITHER  XE(J) .LT. X(1)
+C                 OR      X(I-1) .LE. XE(J) .LT. X(I) .
+C              RESET IR, RECOGNIZING THAT IT WILL BE INCREMENTED BEFORE
+C              CYCLING.
+               IR = MAX(1, I-1)
+C           ENDIF
+C        ENDIF
+   49    CONTINUE
+C
+         JFIRST = J
+C
+C     END OF IR-LOOP.
+C
+   50 CONTINUE
+      IR = IR + 1
+      IF (IR .LE. N)  GO TO 10
+C
+C  NORMAL RETURN.
+C
+ 5000 CONTINUE
+      RETURN
+C
+C  ERROR RETURNS.
+C
+ 5001 CONTINUE
+C     N.LT.2 RETURN.
+      IERR = -1
+      CALL XERMSG ('SLATEC', 'PCHFE',
+     +   'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
+      RETURN
+C
+ 5002 CONTINUE
+C     INCFD.LT.1 RETURN.
+      IERR = -2
+      CALL XERMSG ('SLATEC', 'PCHFE', 'INCREMENT LESS THAN ONE', IERR,
+     +   1)
+      RETURN
+C
+ 5003 CONTINUE
+C     X-ARRAY NOT STRICTLY INCREASING.
+      IERR = -3
+      CALL XERMSG ('SLATEC', 'PCHFE', 'X-ARRAY NOT STRICTLY INCREASING'
+     +   , IERR, 1)
+      RETURN
+C
+ 5004 CONTINUE
+C     NE.LT.1 RETURN.
+      IERR = -4
+      CALL XERMSG ('SLATEC', 'PCHFE',
+     +   'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1)
+      RETURN
+C
+ 5005 CONTINUE
+C     ERROR RETURN FROM CHFEV.
+C   *** THIS CASE SHOULD NEVER OCCUR ***
+      IERR = -5
+      CALL XERMSG ('SLATEC', 'PCHFE',
+     +   'ERROR RETURN FROM CHFEV -- FATAL', IERR, 2)
+      RETURN
+C------------- LAST LINE OF PCHFE FOLLOWS ------------------------------
+      END
Index: /LMDZ4/trunk/libf/bibio/pchfe_95.F90
===================================================================
--- /LMDZ4/trunk/libf/bibio/pchfe_95.F90	(revision 1425)
+++ /LMDZ4/trunk/libf/bibio/pchfe_95.F90	(revision 1425)
@@ -0,0 +1,76 @@
+module PCHFE_95_m
+
+  implicit none
+
+contains
+
+  SUBROUTINE PCHFE_95(X, F, D, SKIP, XE, FE, IERR)
+
+    ! PURPOSE  Evaluate a piecewise cubic Hermite function at an array of
+    !            points.  May be used by itself for Hermite interpolation,
+    !            or as an evaluator for PCHIM or PCHIC.
+    ! CATEGORY  E3
+    ! KEYWORDS  CUBIC HERMITE EVALUATION, HERMITE INTERPOLATION, PCHIP,
+    !             PIECEWISE CUBIC EVALUATION
+
+    !          PCHFE:  Piecewise Cubic Hermite Function Evaluator
+    ! Evaluates the cubic Hermite function defined by  X, F, D  at
+    ! the points  XE.
+
+    use assert_eq_m, only: assert_eq
+
+    REAL, intent(in):: X(:) ! real array of independent variable values
+    ! The elements of X must be strictly increasing.
+
+    REAL, intent(in):: F(:) ! real array of function values
+    ! F(I) is the value corresponding to X(I).
+
+    REAL, intent(in):: D(:) ! real array of derivative values
+    ! D(I) is the value corresponding to X(I).
+
+    LOGICAL, intent(inout):: SKIP 
+    ! request to skip checks for validity of "x"
+    ! If "skip" is false then "pchfe" will check that size(x) >= 2 and
+    ! "x" is in strictly ascending order.
+    ! Setting "skip" to true will save time in case these checks have
+    ! already been performed (say, in "PCHIM" or "PCHIC").
+    ! "SKIP" will be set to TRUE on normal return.
+
+    real, intent(in):: XE(:) ! points at which the function is to be evaluated
+    ! NOTES:
+    ! 1. The evaluation will be most efficient if the elements of XE
+    ! are increasing relative to X.
+    ! That is,   XE(J) .GE. X(I)
+    ! implies    XE(K) .GE. X(I),  all K.GE.J
+    ! 2. If any of the XE are outside the interval [X(1),X(N)], values
+    ! are extrapolated from the nearest extreme cubic, and a warning
+    ! error is returned.
+
+    real, intent(out):: FE(:) ! values of the cubic Hermite function
+    ! defined by X, F, D at the points XE
+
+    integer, intent(out):: IERR ! error flag
+    ! Normal return:
+    ! IERR = 0  no error
+    ! Warning error:
+    ! IERR > 0  means that extrapolation was performed at IERR points
+    ! "Recoverable" errors:
+    !              IERR = -1  if N < 2
+    !              IERR = -3  if the X-array is not strictly increasing
+    !              IERR = -4  if NE < 1
+    ! NOTE: The above errors are checked in the order listed, and
+    ! following arguments have **NOT** been validated.
+
+    ! Variables local to the procedure:
+
+    INTEGER  N, NE
+
+    !---------------------------------------
+
+    n = assert_eq(size(x), size(f), size(d), "PCHFE_95 n")
+    ne = assert_eq(size(xe), size(fe), "PCHFE_95 ne")
+    call PCHFE(N, X, F, D, 1, SKIP, NE, XE, FE, IERR)
+
+  end SUBROUTINE PCHFE_95
+
+end module PCHFE_95_m
Index: /LMDZ4/trunk/libf/bibio/pchsp.F
===================================================================
--- /LMDZ4/trunk/libf/bibio/pchsp.F	(revision 1425)
+++ /LMDZ4/trunk/libf/bibio/pchsp.F	(revision 1425)
@@ -0,0 +1,388 @@
+*DECK PCHSP
+      SUBROUTINE PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
+C***BEGIN PROLOGUE  PCHSP
+C***PURPOSE  Set derivatives needed to determine the Hermite represen-
+C            tation of the cubic spline interpolant to given data, with
+C            specified boundary conditions.
+C***LIBRARY   SLATEC (PCHIP)
+C***CATEGORY  E1A
+C***TYPE      SINGLE PRECISION (PCHSP-S, DPCHSP-D)
+C***KEYWORDS  CUBIC HERMITE INTERPOLATION, PCHIP,
+C             PIECEWISE CUBIC INTERPOLATION, SPLINE INTERPOLATION
+C***AUTHOR  Fritsch, F. N., (LLNL)
+C             Lawrence Livermore National Laboratory
+C             P.O. Box 808  (L-316)
+C             Livermore, CA  94550
+C             FTS 532-4275, (510) 422-4275
+C***DESCRIPTION
+C
+C          PCHSP:   Piecewise Cubic Hermite Spline
+C
+C     Computes the Hermite representation of the cubic spline inter-
+C     polant to the data given in X and F satisfying the boundary
+C     conditions specified by IC and VC.
+C
+C     To facilitate two-dimensional applications, includes an increment
+C     between successive values of the F- and D-arrays.
+C
+C     The resulting piecewise cubic Hermite function may be evaluated
+C     by PCHFE or PCHFD.
+C
+C     NOTE:  This is a modified version of C. de Boor's cubic spline
+C            routine CUBSPL.
+C
+C ----------------------------------------------------------------------
+C
+C  Calling sequence:
+C
+C        PARAMETER  (INCFD = ...)
+C        INTEGER  IC(2), N, NWK, IERR
+C        REAL  VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(NWK)
+C
+C        CALL  PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
+C
+C   Parameters:
+C
+C     IC -- (input) integer array of length 2 specifying desired
+C           boundary conditions:
+C           IC(1) = IBEG, desired condition at beginning of data.
+C           IC(2) = IEND, desired condition at end of data.
+C
+C           IBEG = 0  to set D(1) so that the third derivative is con-
+C              tinuous at X(2).  This is the "not a knot" condition
+C              provided by de Boor's cubic spline routine CUBSPL.
+C              < This is the default boundary condition. >
+C           IBEG = 1  if first derivative at X(1) is given in VC(1).
+C           IBEG = 2  if second derivative at X(1) is given in VC(1).
+C           IBEG = 3  to use the 3-point difference formula for D(1).
+C                     (Reverts to the default b.c. if N.LT.3 .)
+C           IBEG = 4  to use the 4-point difference formula for D(1).
+C                     (Reverts to the default b.c. if N.LT.4 .)
+C          NOTES:
+C           1. An error return is taken if IBEG is out of range.
+C           2. For the "natural" boundary condition, use IBEG=2 and
+C              VC(1)=0.
+C
+C           IEND may take on the same values as IBEG, but applied to
+C           derivative at X(N).  In case IEND = 1 or 2, the value is
+C           given in VC(2).
+C
+C          NOTES:
+C           1. An error return is taken if IEND is out of range.
+C           2. For the "natural" boundary condition, use IEND=2 and
+C              VC(2)=0.
+C
+C     VC -- (input) real array of length 2 specifying desired boundary
+C           values, as indicated above.
+C           VC(1) need be set only if IC(1) = 1 or 2 .
+C           VC(2) need be set only if IC(2) = 1 or 2 .
+C
+C     N -- (input) number of data points.  (Error return if N.LT.2 .)
+C
+C     X -- (input) real array of independent variable values.  The
+C           elements of X must be strictly increasing:
+C                X(I-1) .LT. X(I),  I = 2(1)N.
+C           (Error return if not.)
+C
+C     F -- (input) real array of dependent variable values to be inter-
+C           polated.  F(1+(I-1)*INCFD) is value corresponding to X(I).
+C
+C     D -- (output) real array of derivative values at the data points.
+C           These values will determine the cubic spline interpolant
+C           with the requested boundary conditions.
+C           The value corresponding to X(I) is stored in
+C                D(1+(I-1)*INCFD),  I=1(1)N.
+C           No other entries in D are changed.
+C
+C     INCFD -- (input) increment between successive values in F and D.
+C           This argument is provided primarily for 2-D applications.
+C           (Error return if  INCFD.LT.1 .)
+C
+C     WK -- (scratch) real array of working storage.
+C
+C     NWK -- (input) length of work array.
+C           (Error return if NWK.LT.2*N .)
+C
+C     IERR -- (output) error flag.
+C           Normal return:
+C              IERR = 0  (no errors).
+C           "Recoverable" errors:
+C              IERR = -1  if N.LT.2 .
+C              IERR = -2  if INCFD.LT.1 .
+C              IERR = -3  if the X-array is not strictly increasing.
+C              IERR = -4  if IBEG.LT.0 or IBEG.GT.4 .
+C              IERR = -5  if IEND.LT.0 of IEND.GT.4 .
+C              IERR = -6  if both of the above are true.
+C              IERR = -7  if NWK is too small.
+C               NOTE:  The above errors are checked in the order listed,
+C                   and following arguments have **NOT** been validated.
+C             (The D-array has not been changed in any of these cases.)
+C              IERR = -8  in case of trouble solving the linear system
+C                         for the interior derivative values.
+C             (The D-array may have been changed in this case.)
+C             (             Do **NOT** use it!                )
+C
+C***REFERENCES  Carl de Boor, A Practical Guide to Splines, Springer-
+C                 Verlag, New York, 1978, pp. 53-59.
+C***ROUTINES CALLED  PCHDF, XERMSG
+C***REVISION HISTORY  (YYMMDD)
+C   820503  DATE WRITTEN
+C   820804  Converted to SLATEC library version.
+C   870707  Minor cosmetic changes to prologue.
+C   890411  Added SAVE statements (Vers. 3.2).
+C   890703  Corrected category record.  (WRB)
+C   890831  Modified array declarations.  (WRB)
+C   890831  REVISION DATE from Version 3.2
+C   891214  Prologue converted to Version 4.0 format.  (BAB)
+C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
+C   920429  Revised format and order of references.  (WRB,FNF)
+C***END PROLOGUE  PCHSP
+C  Programming notes:
+C
+C     To produce a double precision version, simply:
+C        a. Change PCHSP to DPCHSP wherever it occurs,
+C        b. Change the real declarations to double precision, and
+C        c. Change the constants ZERO, HALF, ... to double precision.
+C
+C  DECLARE ARGUMENTS.
+C
+      INTEGER  IC(2), N, INCFD, NWK, IERR
+      REAL  VC(2), X(*), F(INCFD,*), D(INCFD,*), WK(2,*)
+C
+C  DECLARE LOCAL VARIABLES.
+C
+      INTEGER  IBEG, IEND, INDEX, J, NM1
+      REAL  G, HALF, ONE, STEMP(3), THREE, TWO, XTEMP(4), ZERO
+      SAVE ZERO, HALF, ONE, TWO, THREE
+      REAL  PCHDF
+C
+      DATA  ZERO /0./,  HALF /0.5/,  ONE /1./,  TWO /2./,  THREE /3./
+C
+C  VALIDITY-CHECK ARGUMENTS.
+C
+C***FIRST EXECUTABLE STATEMENT  PCHSP
+      IF ( N.LT.2 )  GO TO 5001
+      IF ( INCFD.LT.1 )  GO TO 5002
+      DO 1  J = 2, N
+         IF ( X(J).LE.X(J-1) )  GO TO 5003
+    1 CONTINUE
+C
+      IBEG = IC(1)
+      IEND = IC(2)
+      IERR = 0
+      IF ( (IBEG.LT.0).OR.(IBEG.GT.4) )  IERR = IERR - 1
+      IF ( (IEND.LT.0).OR.(IEND.GT.4) )  IERR = IERR - 2
+      IF ( IERR.LT.0 )  GO TO 5004
+C
+C  FUNCTION DEFINITION IS OK -- GO ON.
+C
+      IF ( NWK .LT. 2*N )  GO TO 5007
+C
+C  COMPUTE FIRST DIFFERENCES OF X SEQUENCE AND STORE IN WK(1,.). ALSO,
+C  COMPUTE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN WK(2,.).
+      DO 5  J=2,N
+         WK(1,J) = X(J) - X(J-1)
+         WK(2,J) = (F(1,J) - F(1,J-1))/WK(1,J)
+    5 CONTINUE
+C
+C  SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL.
+C
+      IF ( IBEG.GT.N )  IBEG = 0
+      IF ( IEND.GT.N )  IEND = 0
+C
+C  SET UP FOR BOUNDARY CONDITIONS.
+C
+      IF ( (IBEG.EQ.1).OR.(IBEG.EQ.2) )  THEN
+         D(1,1) = VC(1)
+      ELSE IF (IBEG .GT. 2)  THEN
+C        PICK UP FIRST IBEG POINTS, IN REVERSE ORDER.
+         DO 10  J = 1, IBEG
+            INDEX = IBEG-J+1
+C           INDEX RUNS FROM IBEG DOWN TO 1.
+            XTEMP(J) = X(INDEX)
+            IF (J .LT. IBEG)  STEMP(J) = WK(2,INDEX)
+   10    CONTINUE
+C                 --------------------------------
+         D(1,1) = PCHDF (IBEG, XTEMP, STEMP, IERR)
+C                 --------------------------------
+         IF (IERR .NE. 0)  GO TO 5009
+         IBEG = 1
+      ENDIF
+C
+      IF ( (IEND.EQ.1).OR.(IEND.EQ.2) )  THEN
+         D(1,N) = VC(2)
+      ELSE IF (IEND .GT. 2)  THEN
+C        PICK UP LAST IEND POINTS.
+         DO 15  J = 1, IEND
+            INDEX = N-IEND+J
+C           INDEX RUNS FROM N+1-IEND UP TO N.
+            XTEMP(J) = X(INDEX)
+            IF (J .LT. IEND)  STEMP(J) = WK(2,INDEX+1)
+   15    CONTINUE
+C                 --------------------------------
+         D(1,N) = PCHDF (IEND, XTEMP, STEMP, IERR)
+C                 --------------------------------
+         IF (IERR .NE. 0)  GO TO 5009
+         IEND = 1
+      ENDIF
+C
+C --------------------( BEGIN CODING FROM CUBSPL )--------------------
+C
+C  **** A TRIDIAGONAL LINEAR SYSTEM FOR THE UNKNOWN SLOPES S(J) OF
+C  F  AT X(J), J=1,...,N, IS GENERATED AND THEN SOLVED BY GAUSS ELIM-
+C  INATION, WITH S(J) ENDING UP IN D(1,J), ALL J.
+C     WK(1,.) AND WK(2,.) ARE USED FOR TEMPORARY STORAGE.
+C
+C  CONSTRUCT FIRST EQUATION FROM FIRST BOUNDARY CONDITION, OF THE FORM
+C             WK(2,1)*S(1) + WK(1,1)*S(2) = D(1,1)
+C
+      IF (IBEG .EQ. 0)  THEN
+         IF (N .EQ. 2)  THEN
+C           NO CONDITION AT LEFT END AND N = 2.
+            WK(2,1) = ONE
+            WK(1,1) = ONE
+            D(1,1) = TWO*WK(2,2)
+         ELSE
+C           NOT-A-KNOT CONDITION AT LEFT END AND N .GT. 2.
+            WK(2,1) = WK(1,3)
+            WK(1,1) = WK(1,2) + WK(1,3)
+            D(1,1) =((WK(1,2) + TWO*WK(1,1))*WK(2,2)*WK(1,3)
+     *                        + WK(1,2)**2*WK(2,3)) / WK(1,1)
+         ENDIF
+      ELSE IF (IBEG .EQ. 1)  THEN
+C        SLOPE PRESCRIBED AT LEFT END.
+         WK(2,1) = ONE
+         WK(1,1) = ZERO
+      ELSE
+C        SECOND DERIVATIVE PRESCRIBED AT LEFT END.
+         WK(2,1) = TWO
+         WK(1,1) = ONE
+         D(1,1) = THREE*WK(2,2) - HALF*WK(1,2)*D(1,1)
+      ENDIF
+C
+C  IF THERE ARE INTERIOR KNOTS, GENERATE THE CORRESPONDING EQUATIONS AND
+C  CARRY OUT THE FORWARD PASS OF GAUSS ELIMINATION, AFTER WHICH THE J-TH
+C  EQUATION READS    WK(2,J)*S(J) + WK(1,J)*S(J+1) = D(1,J).
+C
+      NM1 = N-1
+      IF (NM1 .GT. 1)  THEN
+         DO 20 J=2,NM1
+            IF (WK(2,J-1) .EQ. ZERO)  GO TO 5008
+            G = -WK(1,J+1)/WK(2,J-1)
+            D(1,J) = G*D(1,J-1)
+     *                  + THREE*(WK(1,J)*WK(2,J+1) + WK(1,J+1)*WK(2,J))
+            WK(2,J) = G*WK(1,J-1) + TWO*(WK(1,J) + WK(1,J+1))
+   20    CONTINUE
+      ENDIF
+C
+C  CONSTRUCT LAST EQUATION FROM SECOND BOUNDARY CONDITION, OF THE FORM
+C           (-G*WK(2,N-1))*S(N-1) + WK(2,N)*S(N) = D(1,N)
+C
+C     IF SLOPE IS PRESCRIBED AT RIGHT END, ONE CAN GO DIRECTLY TO BACK-
+C     SUBSTITUTION, SINCE ARRAYS HAPPEN TO BE SET UP JUST RIGHT FOR IT
+C     AT THIS POINT.
+      IF (IEND .EQ. 1)  GO TO 30
+C
+      IF (IEND .EQ. 0)  THEN
+         IF (N.EQ.2 .AND. IBEG.EQ.0)  THEN
+C           NOT-A-KNOT AT RIGHT ENDPOINT AND AT LEFT ENDPOINT AND N = 2.
+            D(1,2) = WK(2,2)
+            GO TO 30
+         ELSE IF ((N.EQ.2) .OR. (N.EQ.3 .AND. IBEG.EQ.0))  THEN
+C           EITHER (N=3 AND NOT-A-KNOT ALSO AT LEFT) OR (N=2 AND *NOT*
+C           NOT-A-KNOT AT LEFT END POINT).
+            D(1,N) = TWO*WK(2,N)
+            WK(2,N) = ONE
+            IF (WK(2,N-1) .EQ. ZERO)  GO TO 5008
+            G = -ONE/WK(2,N-1)
+         ELSE
+C           NOT-A-KNOT AND N .GE. 3, AND EITHER N.GT.3 OR  ALSO NOT-A-
+C           KNOT AT LEFT END POINT.
+            G = WK(1,N-1) + WK(1,N)
+C           DO NOT NEED TO CHECK FOLLOWING DENOMINATORS (X-DIFFERENCES).
+            D(1,N) = ((WK(1,N)+TWO*G)*WK(2,N)*WK(1,N-1)
+     *                  + WK(1,N)**2*(F(1,N-1)-F(1,N-2))/WK(1,N-1))/G
+            IF (WK(2,N-1) .EQ. ZERO)  GO TO 5008
+            G = -G/WK(2,N-1)
+            WK(2,N) = WK(1,N-1)
+         ENDIF
+      ELSE
+C        SECOND DERIVATIVE PRESCRIBED AT RIGHT ENDPOINT.
+         D(1,N) = THREE*WK(2,N) + HALF*WK(1,N)*D(1,N)
+         WK(2,N) = TWO
+         IF (WK(2,N-1) .EQ. ZERO)  GO TO 5008
+         G = -ONE/WK(2,N-1)
+      ENDIF
+C
+C  COMPLETE FORWARD PASS OF GAUSS ELIMINATION.
+C
+      WK(2,N) = G*WK(1,N-1) + WK(2,N)
+      IF (WK(2,N) .EQ. ZERO)   GO TO 5008
+      D(1,N) = (G*D(1,N-1) + D(1,N))/WK(2,N)
+C
+C  CARRY OUT BACK SUBSTITUTION
+C
+   30 CONTINUE
+      DO 40 J=NM1,1,-1
+         IF (WK(2,J) .EQ. ZERO)  GO TO 5008
+         D(1,J) = (D(1,J) - WK(1,J)*D(1,J+1))/WK(2,J)
+   40 CONTINUE
+C --------------------(  END  CODING FROM CUBSPL )--------------------
+C
+C  NORMAL RETURN.
+C
+      RETURN
+C
+C  ERROR RETURNS.
+C
+ 5001 CONTINUE
+C     N.LT.2 RETURN.
+      IERR = -1
+      CALL XERMSG ('SLATEC', 'PCHSP',
+     +   'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
+      RETURN
+C
+ 5002 CONTINUE
+C     INCFD.LT.1 RETURN.
+      IERR = -2
+      CALL XERMSG ('SLATEC', 'PCHSP', 'INCREMENT LESS THAN ONE', IERR,
+     +   1)
+      RETURN
+C
+ 5003 CONTINUE
+C     X-ARRAY NOT STRICTLY INCREASING.
+      IERR = -3
+      CALL XERMSG ('SLATEC', 'PCHSP', 'X-ARRAY NOT STRICTLY INCREASING'
+     +   , IERR, 1)
+      RETURN
+C
+ 5004 CONTINUE
+C     IC OUT OF RANGE RETURN.
+      IERR = IERR - 3
+      CALL XERMSG ('SLATEC', 'PCHSP', 'IC OUT OF RANGE', IERR, 1)
+      RETURN
+C
+ 5007 CONTINUE
+C     NWK TOO SMALL RETURN.
+      IERR = -7
+      CALL XERMSG ('SLATEC', 'PCHSP', 'WORK ARRAY TOO SMALL', IERR, 1)
+      RETURN
+C
+ 5008 CONTINUE
+C     SINGULAR SYSTEM.
+C   *** THEORETICALLY, THIS CAN ONLY OCCUR IF SUCCESSIVE X-VALUES   ***
+C   *** ARE EQUAL, WHICH SHOULD ALREADY HAVE BEEN CAUGHT (IERR=-3). ***
+      IERR = -8
+      CALL XERMSG ('SLATEC', 'PCHSP', 'SINGULAR LINEAR SYSTEM', IERR,
+     +   1)
+      RETURN
+C
+ 5009 CONTINUE
+C     ERROR RETURN FROM PCHDF.
+C   *** THIS CASE SHOULD NEVER OCCUR ***
+      IERR = -9
+      CALL XERMSG ('SLATEC', 'PCHSP', 'ERROR RETURN FROM PCHDF', IERR,
+     +   1)
+      RETURN
+C------------- LAST LINE OF PCHSP FOLLOWS ------------------------------
+      END
Index: /LMDZ4/trunk/libf/bibio/pchsp_95.F90
===================================================================
--- /LMDZ4/trunk/libf/bibio/pchsp_95.F90	(revision 1425)
+++ /LMDZ4/trunk/libf/bibio/pchsp_95.F90	(revision 1425)
@@ -0,0 +1,116 @@
+module pchsp_95_m
+
+  implicit none
+
+contains
+
+  function pchsp_95(x, f, ibeg, iend, vc_beg, vc_end)
+
+    ! PURPOSE: Set derivatives needed to determine the Hermite
+    ! representation of the cubic spline interpolant to given data,
+    ! with specified boundary conditions.
+
+    ! Part of the "pchip" package.
+
+    ! CATEGORY: E1A
+
+    ! KEYWORDS: cubic hermite interpolation, piecewise cubic
+    ! interpolation, spline interpolation
+
+    ! DESCRIPTION: "pchsp" stands for "Piecewise Cubic Hermite Spline"
+    ! Computes the Hermite representation of the cubic spline
+    ! interpolant to the data given in X and F satisfying the boundary
+    ! conditions specified by Ibeg, iend, vc_beg and VC_end.
+
+    ! The resulting piecewise cubic Hermite function may be evaluated
+    ! by "pchfe" or "pchfd".
+
+    ! NOTE: This is a modified version of C. de Boor's cubic spline
+    ! routine "cubspl".
+
+    ! REFERENCE: Carl de Boor, A Practical Guide to Splines, Springer,
+    ! 2001, pages 43-47
+
+    use assert_eq_m, only: assert_eq
+
+    real, intent(in):: x(:)
+    ! independent variable values
+    ! The elements of X must be strictly increasing:
+    !                X(I-1) < X(I),  I = 2...N.
+    !           (Error return if not.)
+    ! (error if size(x) < 2)
+
+    real, intent(in):: f(:)
+    !     dependent variable values to be interpolated
+    !  F(I) is value corresponding to X(I).
+
+    INTEGER, intent(in):: ibeg
+    !     desired boundary condition at beginning of data
+
+    !        IBEG = 0  to set pchsp_95(1) so that the third derivative is con-
+    !              tinuous at X(2).  This is the "not a knot" condition
+    !              provided by de Boor's cubic spline routine CUBSPL.
+    !              This is the default boundary condition.
+    !        IBEG = 1  if first derivative at X(1) is given in VC_BEG.
+    !        IBEG = 2  if second derivative at X(1) is given in VC_BEG.
+    !        IBEG = 3  to use the 3-point difference formula for pchsp_95(1).
+    !              (Reverts to the default boundary condition if size(x) < 3 .)
+    !        IBEG = 4  to use the 4-point difference formula for pchsp_95(1).
+    !              (Reverts to the default boundary condition if size(x) < 4 .)
+
+    !          NOTES:
+    !           1. An error return is taken if IBEG is out of range.
+    !           2. For the "natural" boundary condition, use IBEG=2 and
+    !              VC_BEG=0.
+
+    INTEGER, intent(in):: iend
+    !           IC(2) = IEND, desired condition at end of data.
+    !  IEND may take on the same values as IBEG, but applied to
+    !  derivative at X(N). In case IEND = 1 or 2, The value is given in VC_END.
+
+    !          NOTES:
+    !           1. An error return is taken if IEND is out of range.
+    !           2. For the "natural" boundary condition, use IEND=2 and
+    !              VC_END=0.
+
+    REAL, intent(in), optional:: vc_beg
+    ! desired boundary value, as indicated above.
+    !           VC_BEG need be set only if IBEG = 1 or 2 .
+
+    REAL, intent(in), optional:: vc_end
+    ! desired boundary value, as indicated above.
+    !           VC_END need be set only if Iend = 1 or 2 .
+
+    real pchsp_95(size(x))
+    ! derivative values at the data points
+    !           These values will determine the cubic spline interpolant
+    !           with the requested boundary conditions.
+    !           The value corresponding to X(I) is stored in
+    !                PCHSP_95(I),  I=1...N.
+
+    ! LOCAL VARIABLES:
+    real wk(2, size(x)) ! real array of working storage
+    INTEGER n ! number of data points
+    integer ierr, ic(2)
+    real vc(2)
+
+    !-------------------------------------------------------------------
+
+    n = assert_eq(size(x), size(f), "pchsp_95 n")
+    if ((ibeg == 1 .or. ibeg == 2) .and. .not. present(vc_beg)) then
+       print *, "vc_beg required for IBEG = 1 or 2"
+       stop 1
+    end if
+    if ((iend == 1 .or. iend == 2) .and. .not. present(vc_end)) then
+       print *, "vc_end required for IEND = 1 or 2"
+       stop 1
+    end if
+    ic = (/ibeg, iend/)
+    if (present(vc_beg)) vc(1) = vc_beg
+    if (present(vc_end)) vc(2) = vc_end
+    call PCHSP(IC, VC, N, X, F, pchsp_95, 1, WK, size(WK), IERR)
+    if (ierr /= 0) stop 1
+
+  END function pchsp_95
+
+end module pchsp_95_m
Index: /LMDZ4/trunk/libf/bibio/xercnt.F
===================================================================
--- /LMDZ4/trunk/libf/bibio/xercnt.F	(revision 1425)
+++ /LMDZ4/trunk/libf/bibio/xercnt.F	(revision 1425)
@@ -0,0 +1,60 @@
+*DECK XERCNT
+      SUBROUTINE XERCNT (LIBRAR, SUBROU, MESSG, NERR, LEVEL, KONTRL)
+C***BEGIN PROLOGUE  XERCNT
+C***SUBSIDIARY
+C***PURPOSE  Allow user control over handling of errors.
+C***LIBRARY   SLATEC (XERROR)
+C***CATEGORY  R3C
+C***TYPE      ALL (XERCNT-A)
+C***KEYWORDS  ERROR, XERROR
+C***AUTHOR  Jones, R. E., (SNLA)
+C***DESCRIPTION
+C
+C     Abstract
+C        Allows user control over handling of individual errors.
+C        Just after each message is recorded, but before it is
+C        processed any further (i.e., before it is printed or
+C        a decision to abort is made), a call is made to XERCNT.
+C        If the user has provided his own version of XERCNT, he
+C        can then override the value of KONTROL used in processing
+C        this message by redefining its value.
+C        KONTRL may be set to any value from -2 to 2.
+C        The meanings for KONTRL are the same as in XSETF, except
+C        that the value of KONTRL changes only for this message.
+C        If KONTRL is set to a value outside the range from -2 to 2,
+C        it will be moved back into that range.
+C
+C     Description of Parameters
+C
+C      --Input--
+C        LIBRAR - the library that the routine is in.
+C        SUBROU - the subroutine that XERMSG is being called from
+C        MESSG  - the first 20 characters of the error message.
+C        NERR   - same as in the call to XERMSG.
+C        LEVEL  - same as in the call to XERMSG.
+C        KONTRL - the current value of the control flag as set
+C                 by a call to XSETF.
+C
+C      --Output--
+C        KONTRL - the new value of KONTRL.  If KONTRL is not
+C                 defined, it will remain at its original value.
+C                 This changed value of control affects only
+C                 the current occurrence of the current message.
+C
+C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+C                 Error-handling Package, SAND82-0800, Sandia
+C                 Laboratories, 1982.
+C***ROUTINES CALLED  (NONE)
+C***REVISION HISTORY  (YYMMDD)
+C   790801  DATE WRITTEN
+C   861211  REVISION DATE from Version 3.2
+C   891214  Prologue converted to Version 4.0 format.  (BAB)
+C   900206  Routine changed from user-callable to subsidiary.  (WRB)
+C   900510  Changed calling sequence to include LIBRARY and SUBROUTINE
+C           names, changed routine name from XERCTL to XERCNT.  (RWC)
+C   920501  Reformatted the REFERENCES section.  (WRB)
+C***END PROLOGUE  XERCNT
+      CHARACTER*(*) LIBRAR, SUBROU, MESSG
+C***FIRST EXECUTABLE STATEMENT  XERCNT
+      RETURN
+      END
Index: /LMDZ4/trunk/libf/bibio/xerhlt.F
===================================================================
--- /LMDZ4/trunk/libf/bibio/xerhlt.F	(revision 1425)
+++ /LMDZ4/trunk/libf/bibio/xerhlt.F	(revision 1425)
@@ -0,0 +1,39 @@
+*DECK XERHLT
+      SUBROUTINE XERHLT (MESSG)
+C***BEGIN PROLOGUE  XERHLT
+C***SUBSIDIARY
+C***PURPOSE  Abort program execution and print error message.
+C***LIBRARY   SLATEC (XERROR)
+C***CATEGORY  R3C
+C***TYPE      ALL (XERHLT-A)
+C***KEYWORDS  ABORT PROGRAM EXECUTION, ERROR, XERROR
+C***AUTHOR  Jones, R. E., (SNLA)
+C***DESCRIPTION
+C
+C     Abstract
+C        ***Note*** machine dependent routine
+C        XERHLT aborts the execution of the program.
+C        The error message causing the abort is given in the calling
+C        sequence, in case one needs it for printing on a dayfile,
+C        for example.
+C
+C     Description of Parameters
+C        MESSG is as in XERMSG.
+C
+C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+C                 Error-handling Package, SAND82-0800, Sandia
+C                 Laboratories, 1982.
+C***ROUTINES CALLED  (NONE)
+C***REVISION HISTORY  (YYMMDD)
+C   790801  DATE WRITTEN
+C   861211  REVISION DATE from Version 3.2
+C   891214  Prologue converted to Version 4.0 format.  (BAB)
+C   900206  Routine changed from user-callable to subsidiary.  (WRB)
+C   900510  Changed calling sequence to delete length of character
+C           and changed routine name from XERABT to XERHLT.  (RWC)
+C   920501  Reformatted the REFERENCES section.  (WRB)
+C***END PROLOGUE  XERHLT
+      CHARACTER*(*) MESSG
+C***FIRST EXECUTABLE STATEMENT  XERHLT
+      STOP
+      END
Index: /LMDZ4/trunk/libf/bibio/xermsg.F
===================================================================
--- /LMDZ4/trunk/libf/bibio/xermsg.F	(revision 1425)
+++ /LMDZ4/trunk/libf/bibio/xermsg.F	(revision 1425)
@@ -0,0 +1,364 @@
+*DECK XERMSG
+      SUBROUTINE XERMSG (LIBRAR, SUBROU, MESSG, NERR, LEVEL)
+C***BEGIN PROLOGUE  XERMSG
+C***PURPOSE  Process error messages for SLATEC and other libraries.
+C***LIBRARY   SLATEC (XERROR)
+C***CATEGORY  R3C
+C***TYPE      ALL (XERMSG-A)
+C***KEYWORDS  ERROR MESSAGE, XERROR
+C***AUTHOR  Fong, Kirby, (NMFECC at LLNL)
+C***DESCRIPTION
+C
+C   XERMSG processes a diagnostic message in a manner determined by the
+C   value of LEVEL and the current value of the library error control
+C   flag, KONTRL.  See subroutine XSETF for details.
+C
+C    LIBRAR   A character constant (or character variable) with the name
+C             of the library.  This will be 'SLATEC' for the SLATEC
+C             Common Math Library.  The error handling package is
+C             general enough to be used by many libraries
+C             simultaneously, so it is desirable for the routine that
+C             detects and reports an error to identify the library name
+C             as well as the routine name.
+C
+C    SUBROU   A character constant (or character variable) with the name
+C             of the routine that detected the error.  Usually it is the
+C             name of the routine that is calling XERMSG.  There are
+C             some instances where a user callable library routine calls
+C             lower level subsidiary routines where the error is
+C             detected.  In such cases it may be more informative to
+C             supply the name of the routine the user called rather than
+C             the name of the subsidiary routine that detected the
+C             error.
+C
+C    MESSG    A character constant (or character variable) with the text
+C             of the error or warning message.  In the example below,
+C             the message is a character constant that contains a
+C             generic message.
+C
+C                   CALL XERMSG ('SLATEC', 'MMPY',
+C                  *'THE ORDER OF THE MATRIX EXCEEDS THE ROW DIMENSION',
+C                  *3, 1)
+C
+C             It is possible (and is sometimes desirable) to generate a
+C             specific message--e.g., one that contains actual numeric
+C             values.  Specific numeric values can be converted into
+C             character strings using formatted WRITE statements into
+C             character variables.  This is called standard Fortran
+C             internal file I/O and is exemplified in the first three
+C             lines of the following example.  You can also catenate
+C             substrings of characters to construct the error message.
+C             Here is an example showing the use of both writing to
+C             an internal file and catenating character strings.
+C
+C                   CHARACTER*5 CHARN, CHARL
+C                   WRITE (CHARN,10) N
+C                   WRITE (CHARL,10) LDA
+C                10 FORMAT(I5)
+C                   CALL XERMSG ('SLATEC', 'MMPY', 'THE ORDER'//CHARN//
+C                  *   ' OF THE MATRIX EXCEEDS ITS ROW DIMENSION OF'//
+C                  *   CHARL, 3, 1)
+C
+C             There are two subtleties worth mentioning.  One is that
+C             the // for character catenation is used to construct the
+C             error message so that no single character constant is
+C             continued to the next line.  This avoids confusion as to
+C             whether there are trailing blanks at the end of the line.
+C             The second is that by catenating the parts of the message
+C             as an actual argument rather than encoding the entire
+C             message into one large character variable, we avoid
+C             having to know how long the message will be in order to
+C             declare an adequate length for that large character
+C             variable.  XERMSG calls XERPRN to print the message using
+C             multiple lines if necessary.  If the message is very long,
+C             XERPRN will break it into pieces of 72 characters (as
+C             requested by XERMSG) for printing on multiple lines.
+C             Also, XERMSG asks XERPRN to prefix each line with ' *  '
+C             so that the total line length could be 76 characters.
+C             Note also that XERPRN scans the error message backwards
+C             to ignore trailing blanks.  Another feature is that
+C             the substring '$$' is treated as a new line sentinel
+C             by XERPRN.  If you want to construct a multiline
+C             message without having to count out multiples of 72
+C             characters, just use '$$' as a separator.  '$$'
+C             obviously must occur within 72 characters of the
+C             start of each line to have its intended effect since
+C             XERPRN is asked to wrap around at 72 characters in
+C             addition to looking for '$$'.
+C
+C    NERR     An integer value that is chosen by the library routine's
+C             author.  It must be in the range -99 to 999 (three
+C             printable digits).  Each distinct error should have its
+C             own error number.  These error numbers should be described
+C             in the machine readable documentation for the routine.
+C             The error numbers need be unique only within each routine,
+C             so it is reasonable for each routine to start enumerating
+C             errors from 1 and proceeding to the next integer.
+C
+C    LEVEL    An integer value in the range 0 to 2 that indicates the
+C             level (severity) of the error.  Their meanings are
+C
+C            -1  A warning message.  This is used if it is not clear
+C                that there really is an error, but the user's attention
+C                may be needed.  An attempt is made to only print this
+C                message once.
+C
+C             0  A warning message.  This is used if it is not clear
+C                that there really is an error, but the user's attention
+C                may be needed.
+C
+C             1  A recoverable error.  This is used even if the error is
+C                so serious that the routine cannot return any useful
+C                answer.  If the user has told the error package to
+C                return after recoverable errors, then XERMSG will
+C                return to the Library routine which can then return to
+C                the user's routine.  The user may also permit the error
+C                package to terminate the program upon encountering a
+C                recoverable error.
+C
+C             2  A fatal error.  XERMSG will not return to its caller
+C                after it receives a fatal error.  This level should
+C                hardly ever be used; it is much better to allow the
+C                user a chance to recover.  An example of one of the few
+C                cases in which it is permissible to declare a level 2
+C                error is a reverse communication Library routine that
+C                is likely to be called repeatedly until it integrates
+C                across some interval.  If there is a serious error in
+C                the input such that another step cannot be taken and
+C                the Library routine is called again without the input
+C                error having been corrected by the caller, the Library
+C                routine will probably be called forever with improper
+C                input.  In this case, it is reasonable to declare the
+C                error to be fatal.
+C
+C    Each of the arguments to XERMSG is input; none will be modified by
+C    XERMSG.  A routine may make multiple calls to XERMSG with warning
+C    level messages; however, after a call to XERMSG with a recoverable
+C    error, the routine should return to the user.  Do not try to call
+C    XERMSG with a second recoverable error after the first recoverable
+C    error because the error package saves the error number.  The user
+C    can retrieve this error number by calling another entry point in
+C    the error handling package and then clear the error number when
+C    recovering from the error.  Calling XERMSG in succession causes the
+C    old error number to be overwritten by the latest error number.
+C    This is considered harmless for error numbers associated with
+C    warning messages but must not be done for error numbers of serious
+C    errors.  After a call to XERMSG with a recoverable error, the user
+C    must be given a chance to call NUMXER or XERCLR to retrieve or
+C    clear the error number.
+C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+C                 Error-handling Package, SAND82-0800, Sandia
+C                 Laboratories, 1982.
+C***ROUTINES CALLED  FDUMP, J4SAVE, XERCNT, XERHLT, XERPRN, XERSVE
+C***REVISION HISTORY  (YYMMDD)
+C   880101  DATE WRITTEN
+C   880621  REVISED AS DIRECTED AT SLATEC CML MEETING OF FEBRUARY 1988.
+C           THERE ARE TWO BASIC CHANGES.
+C           1.  A NEW ROUTINE, XERPRN, IS USED INSTEAD OF XERPRT TO
+C               PRINT MESSAGES.  THIS ROUTINE WILL BREAK LONG MESSAGES
+C               INTO PIECES FOR PRINTING ON MULTIPLE LINES.  '$$' IS
+C               ACCEPTED AS A NEW LINE SENTINEL.  A PREFIX CAN BE
+C               ADDED TO EACH LINE TO BE PRINTED.  XERMSG USES EITHER
+C               ' ***' OR ' *  ' AND LONG MESSAGES ARE BROKEN EVERY
+C               72 CHARACTERS (AT MOST) SO THAT THE MAXIMUM LINE
+C               LENGTH OUTPUT CAN NOW BE AS GREAT AS 76.
+C           2.  THE TEXT OF ALL MESSAGES IS NOW IN UPPER CASE SINCE THE
+C               FORTRAN STANDARD DOCUMENT DOES NOT ADMIT THE EXISTENCE
+C               OF LOWER CASE.
+C   880708  REVISED AFTER THE SLATEC CML MEETING OF JUNE 29 AND 30.
+C           THE PRINCIPAL CHANGES ARE
+C           1.  CLARIFY COMMENTS IN THE PROLOGUES
+C           2.  RENAME XRPRNT TO XERPRN
+C           3.  REWORK HANDLING OF '$$' IN XERPRN TO HANDLE BLANK LINES
+C               SIMILAR TO THE WAY FORMAT STATEMENTS HANDLE THE /
+C               CHARACTER FOR NEW RECORDS.
+C   890706  REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO
+C           CLEAN UP THE CODING.
+C   890721  REVISED TO USE NEW FEATURE IN XERPRN TO COUNT CHARACTERS IN
+C           PREFIX.
+C   891013  REVISED TO CORRECT COMMENTS.
+C   891214  Prologue converted to Version 4.0 format.  (WRB)
+C   900510  Changed test on NERR to be -9999999 < NERR < 99999999, but
+C           NERR .ne. 0, and on LEVEL to be -2 < LEVEL < 3.  Added
+C           LEVEL=-1 logic, changed calls to XERSAV to XERSVE, and
+C           XERCTL to XERCNT.  (RWC)
+C   920501  Reformatted the REFERENCES section.  (WRB)
+C***END PROLOGUE  XERMSG
+      CHARACTER*(*) LIBRAR, SUBROU, MESSG
+      CHARACTER*8 XLIBR, XSUBR
+      CHARACTER*72  TEMP
+      CHARACTER*20  LFIRST
+C***FIRST EXECUTABLE STATEMENT  XERMSG
+      LKNTRL = J4SAVE (2, 0, .FALSE.)
+      MAXMES = J4SAVE (4, 0, .FALSE.)
+C
+C       LKNTRL IS A LOCAL COPY OF THE CONTROL FLAG KONTRL.
+C       MAXMES IS THE MAXIMUM NUMBER OF TIMES ANY PARTICULAR MESSAGE
+C          SHOULD BE PRINTED.
+C
+C       WE PRINT A FATAL ERROR MESSAGE AND TERMINATE FOR AN ERROR IN
+C          CALLING XERMSG.  THE ERROR NUMBER SHOULD BE POSITIVE,
+C          AND THE LEVEL SHOULD BE BETWEEN 0 AND 2.
+C
+      IF (NERR.LT.-9999999 .OR. NERR.GT.99999999 .OR. NERR.EQ.0 .OR.
+     *   LEVEL.LT.-1 .OR. LEVEL.GT.2) THEN
+         CALL XERPRN (' ***', -1, 'FATAL ERROR IN...$$ ' //
+     *      'XERMSG -- INVALID ERROR NUMBER OR LEVEL$$ '//
+     *      'JOB ABORT DUE TO FATAL ERROR.', 72)
+         CALL XERSVE (' ', ' ', ' ', 0, 0, 0, KDUMMY)
+         CALL XERHLT (' ***XERMSG -- INVALID INPUT')
+         RETURN
+      ENDIF
+C
+C       RECORD THE MESSAGE.
+C
+      I = J4SAVE (1, NERR, .TRUE.)
+      CALL XERSVE (LIBRAR, SUBROU, MESSG, 1, NERR, LEVEL, KOUNT)
+C
+C       HANDLE PRINT-ONCE WARNING MESSAGES.
+C
+      IF (LEVEL.EQ.-1 .AND. KOUNT.GT.1) RETURN
+C
+C       ALLOW TEMPORARY USER OVERRIDE OF THE CONTROL FLAG.
+C
+      XLIBR  = LIBRAR
+      XSUBR  = SUBROU
+      LFIRST = MESSG
+      LERR   = NERR
+      LLEVEL = LEVEL
+      CALL XERCNT (XLIBR, XSUBR, LFIRST, LERR, LLEVEL, LKNTRL)
+C
+      LKNTRL = MAX(-2, MIN(2,LKNTRL))
+      MKNTRL = ABS(LKNTRL)
+C
+C       SKIP PRINTING IF THE CONTROL FLAG VALUE AS RESET IN XERCNT IS
+C       ZERO AND THE ERROR IS NOT FATAL.
+C
+      IF (LEVEL.LT.2 .AND. LKNTRL.EQ.0) GO TO 30
+      IF (LEVEL.EQ.0 .AND. KOUNT.GT.MAXMES) GO TO 30
+      IF (LEVEL.EQ.1 .AND. KOUNT.GT.MAXMES .AND. MKNTRL.EQ.1) GO TO 30
+      IF (LEVEL.EQ.2 .AND. KOUNT.GT.MAX(1,MAXMES)) GO TO 30
+C
+C       ANNOUNCE THE NAMES OF THE LIBRARY AND SUBROUTINE BY BUILDING A
+C       MESSAGE IN CHARACTER VARIABLE TEMP (NOT EXCEEDING 66 CHARACTERS)
+C       AND SENDING IT OUT VIA XERPRN.  PRINT ONLY IF CONTROL FLAG
+C       IS NOT ZERO.
+C
+      IF (LKNTRL .NE. 0) THEN
+         TEMP(1:21) = 'MESSAGE FROM ROUTINE '
+         I = MIN(LEN(SUBROU), 16)
+         TEMP(22:21+I) = SUBROU(1:I)
+         TEMP(22+I:33+I) = ' IN LIBRARY '
+         LTEMP = 33 + I
+         I = MIN(LEN(LIBRAR), 16)
+         TEMP(LTEMP+1:LTEMP+I) = LIBRAR (1:I)
+         TEMP(LTEMP+I+1:LTEMP+I+1) = '.'
+         LTEMP = LTEMP + I + 1
+         CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72)
+      ENDIF
+C
+C       IF LKNTRL IS POSITIVE, PRINT AN INTRODUCTORY LINE BEFORE
+C       PRINTING THE MESSAGE.  THE INTRODUCTORY LINE TELLS THE CHOICE
+C       FROM EACH OF THE FOLLOWING THREE OPTIONS.
+C       1.  LEVEL OF THE MESSAGE
+C              'INFORMATIVE MESSAGE'
+C              'POTENTIALLY RECOVERABLE ERROR'
+C              'FATAL ERROR'
+C       2.  WHETHER CONTROL FLAG WILL ALLOW PROGRAM TO CONTINUE
+C              'PROG CONTINUES'
+C              'PROG ABORTED'
+C       3.  WHETHER OR NOT A TRACEBACK WAS REQUESTED.  (THE TRACEBACK
+C           MAY NOT BE IMPLEMENTED AT SOME SITES, SO THIS ONLY TELLS
+C           WHAT WAS REQUESTED, NOT WHAT WAS DELIVERED.)
+C              'TRACEBACK REQUESTED'
+C              'TRACEBACK NOT REQUESTED'
+C       NOTICE THAT THE LINE INCLUDING FOUR PREFIX CHARACTERS WILL NOT
+C       EXCEED 74 CHARACTERS.
+C       WE SKIP THE NEXT BLOCK IF THE INTRODUCTORY LINE IS NOT NEEDED.
+C
+      IF (LKNTRL .GT. 0) THEN
+C
+C       THE FIRST PART OF THE MESSAGE TELLS ABOUT THE LEVEL.
+C
+         IF (LEVEL .LE. 0) THEN
+            TEMP(1:20) = 'INFORMATIVE MESSAGE,'
+            LTEMP = 20
+         ELSEIF (LEVEL .EQ. 1) THEN
+            TEMP(1:30) = 'POTENTIALLY RECOVERABLE ERROR,'
+            LTEMP = 30
+         ELSE
+            TEMP(1:12) = 'FATAL ERROR,'
+            LTEMP = 12
+         ENDIF
+C
+C       THEN WHETHER THE PROGRAM WILL CONTINUE.
+C
+         IF ((MKNTRL.EQ.2 .AND. LEVEL.GE.1) .OR.
+     *       (MKNTRL.EQ.1 .AND. LEVEL.EQ.2)) THEN
+            TEMP(LTEMP+1:LTEMP+14) = ' PROG ABORTED,'
+            LTEMP = LTEMP + 14
+         ELSE
+            TEMP(LTEMP+1:LTEMP+16) = ' PROG CONTINUES,'
+            LTEMP = LTEMP + 16
+         ENDIF
+C
+C       FINALLY TELL WHETHER THERE SHOULD BE A TRACEBACK.
+C
+         IF (LKNTRL .GT. 0) THEN
+            TEMP(LTEMP+1:LTEMP+20) = ' TRACEBACK REQUESTED'
+            LTEMP = LTEMP + 20
+         ELSE
+            TEMP(LTEMP+1:LTEMP+24) = ' TRACEBACK NOT REQUESTED'
+            LTEMP = LTEMP + 24
+         ENDIF
+         CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72)
+      ENDIF
+C
+C       NOW SEND OUT THE MESSAGE.
+C
+      CALL XERPRN (' *  ', -1, MESSG, 72)
+C
+C       IF LKNTRL IS POSITIVE, WRITE THE ERROR NUMBER AND REQUEST A
+C          TRACEBACK.
+C
+      IF (LKNTRL .GT. 0) THEN
+         WRITE (TEMP, '(''ERROR NUMBER = '', I8)') NERR
+         DO 10 I=16,22
+            IF (TEMP(I:I) .NE. ' ') GO TO 20
+   10    CONTINUE
+C
+   20    CALL XERPRN (' *  ', -1, TEMP(1:15) // TEMP(I:23), 72)
+         CALL FDUMP
+      ENDIF
+C
+C       IF LKNTRL IS NOT ZERO, PRINT A BLANK LINE AND AN END OF MESSAGE.
+C
+      IF (LKNTRL .NE. 0) THEN
+         CALL XERPRN (' *  ', -1, ' ', 72)
+         CALL XERPRN (' ***', -1, 'END OF MESSAGE', 72)
+         CALL XERPRN ('    ',  0, ' ', 72)
+      ENDIF
+C
+C       IF THE ERROR IS NOT FATAL OR THE ERROR IS RECOVERABLE AND THE
+C       CONTROL FLAG IS SET FOR RECOVERY, THEN RETURN.
+C
+   30 IF (LEVEL.LE.0 .OR. (LEVEL.EQ.1 .AND. MKNTRL.LE.1)) RETURN
+C
+C       THE PROGRAM WILL BE STOPPED DUE TO AN UNRECOVERED ERROR OR A
+C       FATAL ERROR.  PRINT THE REASON FOR THE ABORT AND THE ERROR
+C       SUMMARY IF THE CONTROL FLAG AND THE MAXIMUM ERROR COUNT PERMIT.
+C
+      IF (LKNTRL.GT.0 .AND. KOUNT.LT.MAX(1,MAXMES)) THEN
+         IF (LEVEL .EQ. 1) THEN
+            CALL XERPRN
+     *         (' ***', -1, 'JOB ABORT DUE TO UNRECOVERED ERROR.', 72)
+         ELSE
+            CALL XERPRN(' ***', -1, 'JOB ABORT DUE TO FATAL ERROR.', 72)
+         ENDIF
+         CALL XERSVE (' ', ' ', ' ', -1, 0, 0, KDUMMY)
+         CALL XERHLT (' ')
+      ELSE
+         CALL XERHLT (MESSG)
+      ENDIF
+      RETURN
+      END
Index: /LMDZ4/trunk/libf/bibio/xerprn.F
===================================================================
--- /LMDZ4/trunk/libf/bibio/xerprn.F	(revision 1425)
+++ /LMDZ4/trunk/libf/bibio/xerprn.F	(revision 1425)
@@ -0,0 +1,228 @@
+*DECK XERPRN
+      SUBROUTINE XERPRN (PREFIX, NPREF, MESSG, NWRAP)
+C***BEGIN PROLOGUE  XERPRN
+C***SUBSIDIARY
+C***PURPOSE  Print error messages processed by XERMSG.
+C***LIBRARY   SLATEC (XERROR)
+C***CATEGORY  R3C
+C***TYPE      ALL (XERPRN-A)
+C***KEYWORDS  ERROR MESSAGES, PRINTING, XERROR
+C***AUTHOR  Fong, Kirby, (NMFECC at LLNL)
+C***DESCRIPTION
+C
+C This routine sends one or more lines to each of the (up to five)
+C logical units to which error messages are to be sent.  This routine
+C is called several times by XERMSG, sometimes with a single line to
+C print and sometimes with a (potentially very long) message that may
+C wrap around into multiple lines.
+C
+C PREFIX  Input argument of type CHARACTER.  This argument contains
+C         characters to be put at the beginning of each line before
+C         the body of the message.  No more than 16 characters of
+C         PREFIX will be used.
+C
+C NPREF   Input argument of type INTEGER.  This argument is the number
+C         of characters to use from PREFIX.  If it is negative, the
+C         intrinsic function LEN is used to determine its length.  If
+C         it is zero, PREFIX is not used.  If it exceeds 16 or if
+C         LEN(PREFIX) exceeds 16, only the first 16 characters will be
+C         used.  If NPREF is positive and the length of PREFIX is less
+C         than NPREF, a copy of PREFIX extended with blanks to length
+C         NPREF will be used.
+C
+C MESSG   Input argument of type CHARACTER.  This is the text of a
+C         message to be printed.  If it is a long message, it will be
+C         broken into pieces for printing on multiple lines.  Each line
+C         will start with the appropriate prefix and be followed by a
+C         piece of the message.  NWRAP is the number of characters per
+C         piece; that is, after each NWRAP characters, we break and
+C         start a new line.  In addition the characters '$$' embedded
+C         in MESSG are a sentinel for a new line.  The counting of
+C         characters up to NWRAP starts over for each new line.  The
+C         value of NWRAP typically used by XERMSG is 72 since many
+C         older error messages in the SLATEC Library are laid out to
+C         rely on wrap-around every 72 characters.
+C
+C NWRAP   Input argument of type INTEGER.  This gives the maximum size
+C         piece into which to break MESSG for printing on multiple
+C         lines.  An embedded '$$' ends a line, and the count restarts
+C         at the following character.  If a line break does not occur
+C         on a blank (it would split a word) that word is moved to the
+C         next line.  Values of NWRAP less than 16 will be treated as
+C         16.  Values of NWRAP greater than 132 will be treated as 132.
+C         The actual line length will be NPREF + NWRAP after NPREF has
+C         been adjusted to fall between 0 and 16 and NWRAP has been
+C         adjusted to fall between 16 and 132.
+C
+C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+C                 Error-handling Package, SAND82-0800, Sandia
+C                 Laboratories, 1982.
+C***ROUTINES CALLED  I1MACH, XGETUA
+C***REVISION HISTORY  (YYMMDD)
+C   880621  DATE WRITTEN
+C   880708  REVISED AFTER THE SLATEC CML SUBCOMMITTEE MEETING OF
+C           JUNE 29 AND 30 TO CHANGE THE NAME TO XERPRN AND TO REWORK
+C           THE HANDLING OF THE NEW LINE SENTINEL TO BEHAVE LIKE THE
+C           SLASH CHARACTER IN FORMAT STATEMENTS.
+C   890706  REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO
+C           STREAMLINE THE CODING AND FIX A BUG THAT CAUSED EXTRA BLANK
+C           LINES TO BE PRINTED.
+C   890721  REVISED TO ADD A NEW FEATURE.  A NEGATIVE VALUE OF NPREF
+C           CAUSES LEN(PREFIX) TO BE USED AS THE LENGTH.
+C   891013  REVISED TO CORRECT ERROR IN CALCULATING PREFIX LENGTH.
+C   891214  Prologue converted to Version 4.0 format.  (WRB)
+C   900510  Added code to break messages between words.  (RWC)
+C   920501  Reformatted the REFERENCES section.  (WRB)
+C***END PROLOGUE  XERPRN
+      CHARACTER*(*) PREFIX, MESSG
+      INTEGER NPREF, NWRAP
+      CHARACTER*148 CBUFF
+      INTEGER IU(5), NUNIT
+      CHARACTER*2 NEWLIN
+      PARAMETER (NEWLIN = '$$')
+C***FIRST EXECUTABLE STATEMENT  XERPRN
+      CALL XGETUA(IU,NUNIT)
+C
+C       A ZERO VALUE FOR A LOGICAL UNIT NUMBER MEANS TO USE THE STANDARD
+C       ERROR MESSAGE UNIT INSTEAD.  I1MACH(4) RETRIEVES THE STANDARD
+C       ERROR MESSAGE UNIT.
+C
+      N = I1MACH(4)
+      DO 10 I=1,NUNIT
+         IF (IU(I) .EQ. 0) IU(I) = N
+   10 CONTINUE
+C
+C       LPREF IS THE LENGTH OF THE PREFIX.  THE PREFIX IS PLACED AT THE
+C       BEGINNING OF CBUFF, THE CHARACTER BUFFER, AND KEPT THERE DURING
+C       THE REST OF THIS ROUTINE.
+C
+      IF ( NPREF .LT. 0 ) THEN
+         LPREF = LEN(PREFIX)
+      ELSE
+         LPREF = NPREF
+      ENDIF
+      LPREF = MIN(16, LPREF)
+      IF (LPREF .NE. 0) CBUFF(1:LPREF) = PREFIX
+C
+C       LWRAP IS THE MAXIMUM NUMBER OF CHARACTERS WE WANT TO TAKE AT ONE
+C       TIME FROM MESSG TO PRINT ON ONE LINE.
+C
+      LWRAP = MAX(16, MIN(132, NWRAP))
+C
+C       SET LENMSG TO THE LENGTH OF MESSG, IGNORE ANY TRAILING BLANKS.
+C
+      LENMSG = LEN(MESSG)
+      N = LENMSG
+      DO 20 I=1,N
+         IF (MESSG(LENMSG:LENMSG) .NE. ' ') GO TO 30
+         LENMSG = LENMSG - 1
+   20 CONTINUE
+   30 CONTINUE
+C
+C       IF THE MESSAGE IS ALL BLANKS, THEN PRINT ONE BLANK LINE.
+C
+      IF (LENMSG .EQ. 0) THEN
+         CBUFF(LPREF+1:LPREF+1) = ' '
+         DO 40 I=1,NUNIT
+            WRITE(IU(I), '(A)') CBUFF(1:LPREF+1)
+   40    CONTINUE
+         RETURN
+      ENDIF
+C
+C       SET NEXTC TO THE POSITION IN MESSG WHERE THE NEXT SUBSTRING
+C       STARTS.  FROM THIS POSITION WE SCAN FOR THE NEW LINE SENTINEL.
+C       WHEN NEXTC EXCEEDS LENMSG, THERE IS NO MORE TO PRINT.
+C       WE LOOP BACK TO LABEL 50 UNTIL ALL PIECES HAVE BEEN PRINTED.
+C
+C       WE LOOK FOR THE NEXT OCCURRENCE OF THE NEW LINE SENTINEL.  THE
+C       INDEX INTRINSIC FUNCTION RETURNS ZERO IF THERE IS NO OCCURRENCE
+C       OR IF THE LENGTH OF THE FIRST ARGUMENT IS LESS THAN THE LENGTH
+C       OF THE SECOND ARGUMENT.
+C
+C       THERE ARE SEVERAL CASES WHICH SHOULD BE CHECKED FOR IN THE
+C       FOLLOWING ORDER.  WE ARE ATTEMPTING TO SET LPIECE TO THE NUMBER
+C       OF CHARACTERS THAT SHOULD BE TAKEN FROM MESSG STARTING AT
+C       POSITION NEXTC.
+C
+C       LPIECE .EQ. 0   THE NEW LINE SENTINEL DOES NOT OCCUR IN THE
+C                       REMAINDER OF THE CHARACTER STRING.  LPIECE
+C                       SHOULD BE SET TO LWRAP OR LENMSG+1-NEXTC,
+C                       WHICHEVER IS LESS.
+C
+C       LPIECE .EQ. 1   THE NEW LINE SENTINEL STARTS AT MESSG(NEXTC:
+C                       NEXTC).  LPIECE IS EFFECTIVELY ZERO, AND WE
+C                       PRINT NOTHING TO AVOID PRODUCING UNNECESSARY
+C                       BLANK LINES.  THIS TAKES CARE OF THE SITUATION
+C                       WHERE THE LIBRARY ROUTINE HAS A MESSAGE OF
+C                       EXACTLY 72 CHARACTERS FOLLOWED BY A NEW LINE
+C                       SENTINEL FOLLOWED BY MORE CHARACTERS.  NEXTC
+C                       SHOULD BE INCREMENTED BY 2.
+C
+C       LPIECE .GT. LWRAP+1  REDUCE LPIECE TO LWRAP.
+C
+C       ELSE            THIS LAST CASE MEANS 2 .LE. LPIECE .LE. LWRAP+1
+C                       RESET LPIECE = LPIECE-1.  NOTE THAT THIS
+C                       PROPERLY HANDLES THE END CASE WHERE LPIECE .EQ.
+C                       LWRAP+1.  THAT IS, THE SENTINEL FALLS EXACTLY
+C                       AT THE END OF A LINE.
+C
+      NEXTC = 1
+   50 LPIECE = INDEX(MESSG(NEXTC:LENMSG), NEWLIN)
+      IF (LPIECE .EQ. 0) THEN
+C
+C       THERE WAS NO NEW LINE SENTINEL FOUND.
+C
+         IDELTA = 0
+         LPIECE = MIN(LWRAP, LENMSG+1-NEXTC)
+         IF (LPIECE .LT. LENMSG+1-NEXTC) THEN
+            DO 52 I=LPIECE+1,2,-1
+               IF (MESSG(NEXTC+I-1:NEXTC+I-1) .EQ. ' ') THEN
+                  LPIECE = I-1
+                  IDELTA = 1
+                  GOTO 54
+               ENDIF
+   52       CONTINUE
+         ENDIF
+   54    CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
+         NEXTC = NEXTC + LPIECE + IDELTA
+      ELSEIF (LPIECE .EQ. 1) THEN
+C
+C       WE HAVE A NEW LINE SENTINEL AT MESSG(NEXTC:NEXTC+1).
+C       DON'T PRINT A BLANK LINE.
+C
+         NEXTC = NEXTC + 2
+         GO TO 50
+      ELSEIF (LPIECE .GT. LWRAP+1) THEN
+C
+C       LPIECE SHOULD BE SET DOWN TO LWRAP.
+C
+         IDELTA = 0
+         LPIECE = LWRAP
+         DO 56 I=LPIECE+1,2,-1
+            IF (MESSG(NEXTC+I-1:NEXTC+I-1) .EQ. ' ') THEN
+               LPIECE = I-1
+               IDELTA = 1
+               GOTO 58
+            ENDIF
+   56    CONTINUE
+   58    CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
+         NEXTC = NEXTC + LPIECE + IDELTA
+      ELSE
+C
+C       IF WE ARRIVE HERE, IT MEANS 2 .LE. LPIECE .LE. LWRAP+1.
+C       WE SHOULD DECREMENT LPIECE BY ONE.
+C
+         LPIECE = LPIECE - 1
+         CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
+         NEXTC  = NEXTC + LPIECE + 2
+      ENDIF
+C
+C       PRINT
+C
+      DO 60 I=1,NUNIT
+         WRITE(IU(I), '(A)') CBUFF(1:LPREF+LPIECE)
+   60 CONTINUE
+C
+      IF (NEXTC .LE. LENMSG) GO TO 50
+      RETURN
+      END
Index: /LMDZ4/trunk/libf/bibio/xersve.F
===================================================================
--- /LMDZ4/trunk/libf/bibio/xersve.F	(revision 1425)
+++ /LMDZ4/trunk/libf/bibio/xersve.F	(revision 1425)
@@ -0,0 +1,155 @@
+*DECK XERSVE
+      SUBROUTINE XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL,
+     +   ICOUNT)
+C***BEGIN PROLOGUE  XERSVE
+C***SUBSIDIARY
+C***PURPOSE  Record that an error has occurred.
+C***LIBRARY   SLATEC (XERROR)
+C***CATEGORY  R3
+C***TYPE      ALL (XERSVE-A)
+C***KEYWORDS  ERROR, XERROR
+C***AUTHOR  Jones, R. E., (SNLA)
+C***DESCRIPTION
+C
+C *Usage:
+C
+C        INTEGER  KFLAG, NERR, LEVEL, ICOUNT
+C        CHARACTER * (len) LIBRAR, SUBROU, MESSG
+C
+C        CALL XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL, ICOUNT)
+C
+C *Arguments:
+C
+C        LIBRAR :IN    is the library that the message is from.
+C        SUBROU :IN    is the subroutine that the message is from.
+C        MESSG  :IN    is the message to be saved.
+C        KFLAG  :IN    indicates the action to be performed.
+C                      when KFLAG > 0, the message in MESSG is saved.
+C                      when KFLAG=0 the tables will be dumped and
+C                      cleared.
+C                      when KFLAG < 0, the tables will be dumped and
+C                      not cleared.
+C        NERR   :IN    is the error number.
+C        LEVEL  :IN    is the error severity.
+C        ICOUNT :OUT   the number of times this message has been seen,
+C                      or zero if the table has overflowed and does not
+C                      contain this message specifically.  When KFLAG=0,
+C                      ICOUNT will not be altered.
+C
+C *Description:
+C
+C   Record that this error occurred and possibly dump and clear the
+C   tables.
+C
+C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+C                 Error-handling Package, SAND82-0800, Sandia
+C                 Laboratories, 1982.
+C***ROUTINES CALLED  I1MACH, XGETUA
+C***REVISION HISTORY  (YYMMDD)
+C   800319  DATE WRITTEN
+C   861211  REVISION DATE from Version 3.2
+C   891214  Prologue converted to Version 4.0 format.  (BAB)
+C   900413  Routine modified to remove reference to KFLAG.  (WRB)
+C   900510  Changed to add LIBRARY NAME and SUBROUTINE to calling
+C           sequence, use IF-THEN-ELSE, make number of saved entries
+C           easily changeable, changed routine name from XERSAV to
+C           XERSVE.  (RWC)
+C   910626  Added LIBTAB and SUBTAB to SAVE statement.  (BKS)
+C   920501  Reformatted the REFERENCES section.  (WRB)
+C***END PROLOGUE  XERSVE
+      PARAMETER (LENTAB=10)
+      INTEGER LUN(5)
+      CHARACTER*(*) LIBRAR, SUBROU, MESSG
+      CHARACTER*8  LIBTAB(LENTAB), SUBTAB(LENTAB), LIB, SUB
+      CHARACTER*20 MESTAB(LENTAB), MES
+      DIMENSION NERTAB(LENTAB), LEVTAB(LENTAB), KOUNT(LENTAB)
+      SAVE LIBTAB, SUBTAB, MESTAB, NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG
+      DATA KOUNTX/0/, NMSG/0/
+C***FIRST EXECUTABLE STATEMENT  XERSVE
+C
+      IF (KFLAG.LE.0) THEN
+C
+C        Dump the table.
+C
+         IF (NMSG.EQ.0) RETURN
+C
+C        Print to each unit.
+C
+         CALL XGETUA (LUN, NUNIT)
+         DO 20 KUNIT = 1,NUNIT
+            IUNIT = LUN(KUNIT)
+            IF (IUNIT.EQ.0) IUNIT = I1MACH(4)
+C
+C           Print the table header.
+C
+            WRITE (IUNIT,9000)
+C
+C           Print body of table.
+C
+            DO 10 I = 1,NMSG
+               WRITE (IUNIT,9010) LIBTAB(I), SUBTAB(I), MESTAB(I),
+     *            NERTAB(I),LEVTAB(I),KOUNT(I)
+   10       CONTINUE
+C
+C           Print number of other errors.
+C
+            IF (KOUNTX.NE.0) WRITE (IUNIT,9020) KOUNTX
+            WRITE (IUNIT,9030)
+   20    CONTINUE
+C
+C        Clear the error tables.
+C
+         IF (KFLAG.EQ.0) THEN
+            NMSG = 0
+            KOUNTX = 0
+         ENDIF
+      ELSE
+C
+C        PROCESS A MESSAGE...
+C        SEARCH FOR THIS MESSG, OR ELSE AN EMPTY SLOT FOR THIS MESSG,
+C        OR ELSE DETERMINE THAT THE ERROR TABLE IS FULL.
+C
+         LIB = LIBRAR
+         SUB = SUBROU
+         MES = MESSG
+         DO 30 I = 1,NMSG
+            IF (LIB.EQ.LIBTAB(I) .AND. SUB.EQ.SUBTAB(I) .AND.
+     *         MES.EQ.MESTAB(I) .AND. NERR.EQ.NERTAB(I) .AND.
+     *         LEVEL.EQ.LEVTAB(I)) THEN
+                  KOUNT(I) = KOUNT(I) + 1
+                  ICOUNT = KOUNT(I)
+                  RETURN
+            ENDIF
+   30    CONTINUE
+C
+         IF (NMSG.LT.LENTAB) THEN
+C
+C           Empty slot found for new message.
+C
+            NMSG = NMSG + 1
+            LIBTAB(I) = LIB
+            SUBTAB(I) = SUB
+            MESTAB(I) = MES
+            NERTAB(I) = NERR
+            LEVTAB(I) = LEVEL
+            KOUNT (I) = 1
+            ICOUNT    = 1
+         ELSE
+C
+C           Table is full.
+C
+            KOUNTX = KOUNTX+1
+            ICOUNT = 0
+         ENDIF
+      ENDIF
+      RETURN
+C
+C     Formats.
+C
+ 9000 FORMAT ('0          ERROR MESSAGE SUMMARY' /
+     +   ' LIBRARY    SUBROUTINE MESSAGE START             NERR',
+     +   '     LEVEL     COUNT')
+ 9010 FORMAT (1X,A,3X,A,3X,A,3I10)
+ 9020 FORMAT ('0OTHER ERRORS NOT INDIVIDUALLY TABULATED = ', I10)
+ 9030 FORMAT (1X)
+      END
Index: /LMDZ4/trunk/libf/bibio/xgetua.F
===================================================================
--- /LMDZ4/trunk/libf/bibio/xgetua.F	(revision 1425)
+++ /LMDZ4/trunk/libf/bibio/xgetua.F	(revision 1425)
@@ -0,0 +1,51 @@
+*DECK XGETUA
+      SUBROUTINE XGETUA (IUNITA, N)
+C***BEGIN PROLOGUE  XGETUA
+C***PURPOSE  Return unit number(s) to which error messages are being
+C            sent.
+C***LIBRARY   SLATEC (XERROR)
+C***CATEGORY  R3C
+C***TYPE      ALL (XGETUA-A)
+C***KEYWORDS  ERROR, XERROR
+C***AUTHOR  Jones, R. E., (SNLA)
+C***DESCRIPTION
+C
+C     Abstract
+C        XGETUA may be called to determine the unit number or numbers
+C        to which error messages are being sent.
+C        These unit numbers may have been set by a call to XSETUN,
+C        or a call to XSETUA, or may be a default value.
+C
+C     Description of Parameters
+C      --Output--
+C        IUNIT - an array of one to five unit numbers, depending
+C                on the value of N.  A value of zero refers to the
+C                default unit, as defined by the I1MACH machine
+C                constant routine.  Only IUNIT(1),...,IUNIT(N) are
+C                defined by XGETUA.  The values of IUNIT(N+1),...,
+C                IUNIT(5) are not defined (for N .LT. 5) or altered
+C                in any way by XGETUA.
+C        N     - the number of units to which copies of the
+C                error messages are being sent.  N will be in the
+C                range from 1 to 5.
+C
+C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+C                 Error-handling Package, SAND82-0800, Sandia
+C                 Laboratories, 1982.
+C***ROUTINES CALLED  J4SAVE
+C***REVISION HISTORY  (YYMMDD)
+C   790801  DATE WRITTEN
+C   861211  REVISION DATE from Version 3.2
+C   891214  Prologue converted to Version 4.0 format.  (BAB)
+C   920501  Reformatted the REFERENCES section.  (WRB)
+C***END PROLOGUE  XGETUA
+      DIMENSION IUNITA(5)
+C***FIRST EXECUTABLE STATEMENT  XGETUA
+      N = J4SAVE(5,0,.FALSE.)
+      DO 30 I=1,N
+         INDEX = I+4
+         IF (I.EQ.1) INDEX = 3
+         IUNITA(I) = J4SAVE(INDEX,0,.FALSE.)
+   30 CONTINUE
+      RETURN
+      END
Index: /LMDZ4/trunk/libf/dyn3d/abort_gcm.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3d/abort_gcm.F	(revision 1424)
+++ /LMDZ4/trunk/libf/dyn3d/abort_gcm.F	(revision 1425)
@@ -26,6 +26,5 @@
       character(len=*) message
 
-!      write(lunout,*) 'in abort_gcm'
-      write(6,*) 'in abort_gcm'
+      write(lunout,*) 'in abort_gcm'
 #ifdef CPP_IOIPSL
       call histclo
@@ -37,18 +36,12 @@
 c     call histclo(4)
 c     call histclo(5)
-c     write(lunout,*) 'Stopping in ', modname
-c     write(lunout,*) 'Reason = ',message
-c     if (ierr .eq. 0) then
-c       write(lunout,*) 'Everything is cool'
-c     else
-c       write(lunout,*) 'Houston, we have a problem ', ierr
-c     endif
-      write(6,*) 'Stopping in ', modname
-      write(6,*) 'Reason = ',message
+      write(lunout,*) 'Stopping in ', modname
+      write(lunout,*) 'Reason = ',message
       if (ierr .eq. 0) then
-        write(6,*) 'Everything is cool'
+        write(lunout,*) 'Everything is cool'
+        stop
       else
-        write(6,*) 'Houston, we have a problem ', ierr
+        write(lunout,*) 'Houston, we have a problem ', ierr
+        stop 1
       endif
-      STOP
       END
Index: /LMDZ4/trunk/libf/dyn3d/ce0l.F90
===================================================================
--- /LMDZ4/trunk/libf/dyn3d/ce0l.F90	(revision 1424)
+++ /LMDZ4/trunk/libf/dyn3d/ce0l.F90	(revision 1425)
@@ -93,5 +93,4 @@
 #endif
 ! of #ifndef CPP_EARTH #else
-  STOP
 
 END PROGRAM ce0l
Index: /LMDZ4/trunk/libf/dyn3d/etat0_netcdf.F90
===================================================================
--- /LMDZ4/trunk/libf/dyn3d/etat0_netcdf.F90	(revision 1424)
+++ /LMDZ4/trunk/libf/dyn3d/etat0_netcdf.F90	(revision 1425)
@@ -106,5 +106,5 @@
   INTEGER :: iflag_thermals_ed, iflag_thermals_optflux
   REAL    :: tau_thermals, solarlong0,  seuil_inversion
-  INTEGER :: read_climoz ! read ozone climatology 
+  INTEGER :: read_climoz ! read ozone climatology
 !  Allowed values are 0, 1 and 2
 !     0: do not read an ozone climatology
@@ -140,5 +140,5 @@
                    iflag_thermals_ed,iflag_thermals_optflux,            &
                    iflag_coupl,iflag_clos,iflag_wake, read_climoz,      &
-                   alp_offset )
+                   alp_offset)
 
 ! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value)
@@ -291,19 +291,19 @@
   q3d(:,:,:,:) = 0.0 ; q3d(:,:,:,1) = qd(:,:,:)
 
-!     Parameterization of ozone chemistry:
-!     Look for ozone tracer:
-      i = 1
-      DO
-         found = tname(i)=="O3" .OR. tname(i)=="o3"
-         if (found .or. i == nqtot) exit
-         i = i + 1
-      end do
-      if (found) then
-         call regr_lat_time_coefoz
-         call press_coefoz
-         call regr_pr_o3(p3d, q3d(:, :, :, i))
-!     Convert from mole fraction to mass fraction:
-         q3d(:, :, :, i) = q3d(:, :, :, i)  * 48. / 29.
-      end if
+! Parameterization of ozone chemistry:
+! Look for ozone tracer:
+  i = 1
+  DO
+    found = tname(i)=="O3" .OR. tname(i)=="o3"
+    if (found .or. i == nqtot) exit
+    i = i + 1
+  end do
+  if (found) then
+    call regr_lat_time_coefoz
+    call press_coefoz
+    call regr_pr_o3(p3d, q3d(:, :, :, i))
+!   Convert from mole fraction to mass fraction:
+    q3d(:, :, :, i) = q3d(:, :, :, i)  * 48. / 29.
+  end if
 
 !--- OZONE CLIMATOLOGY
Index: /LMDZ4/trunk/libf/dyn3d/limit_netcdf.F90
===================================================================
--- /LMDZ4/trunk/libf/dyn3d/limit_netcdf.F90	(revision 1424)
+++ /LMDZ4/trunk/libf/dyn3d/limit_netcdf.F90	(revision 1425)
@@ -267,5 +267,4 @@
 #endif
 ! of #ifdef CPP_EARTH
-  STOP
 
 
@@ -281,5 +280,5 @@
 SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask, lCPL)
 !
-!-------------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
 ! Comments:
 !   There are two assumptions concerning the NetCDF files, that are satisfied
@@ -287,11 +286,15 @@
 !     1) The last dimension of the variables used is the time record.
 !     2) Dimensional variables have the same names as corresponding dimensions.
-!-------------------------------------------------------------------------------
-  USE netcdf, ONLY : NF90_OPEN,    NF90_INQ_VARID, NF90_INQUIRE_VARIABLE,      &
-                     NF90_CLOSE,   NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION,     &
-                     NF90_GET_VAR, NF90_GET_ATT
+!-----------------------------------------------------------------------------
+  USE netcdf, ONLY: NF90_OPEN, NF90_INQ_VARID, NF90_INQUIRE_VARIABLE, &
+       NF90_CLOSE, NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, NF90_GET_VAR, &
+       NF90_GET_ATT
   USE dimphy, ONLY : klon
   USE phys_state_var_mod, ONLY : pctsrf
   USE control_mod
+  use pchsp_95_m, only: pchsp_95
+  use pchfe_95_m, only: pchfe_95
+  use arth_m, only: arth
+
   IMPLICIT NONE
 #include "dimensions.h"
@@ -300,5 +303,5 @@
 #include "indicesol.h"
 #include "iniprint.h"
-!-------------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
 ! Arguments:
   CHARACTER(LEN=*),  INTENT(IN)     :: fnam     ! NetCDF file name
@@ -306,12 +309,12 @@
   LOGICAL,           INTENT(IN)     :: ibar     ! interp on pressure levels
   INTEGER,           INTENT(IN)     :: ndays    ! current year number of days
-  REAL,    POINTER,  DIMENSION(:,:) :: champo   ! output field = f(t)
-  LOGICAL, OPTIONAL, INTENT(IN)     :: flag     ! extrapol. (SST)  old ice (SIC)
-  REAL,    OPTIONAL, DIMENSION(iim,jjp1), INTENT(IN) :: mask
+  REAL,    POINTER,  DIMENSION(:, :) :: champo   ! output field = f(t)
+  LOGICAL, OPTIONAL, INTENT(IN)     :: flag     ! extrapol. (SST) old ice (SIC)
+  REAL,    OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask
   LOGICAL, OPTIONAL, INTENT(IN)     :: lCPL     ! Coupled model flag (for ICE)
-!-------------------------------------------------------------------------------
+!------------------------------------------------------------------------------
 ! Local variables:
 !--- NetCDF
-  INTEGER :: ierr, ncid, varid                  ! NetCDF identifiers
+  INTEGER :: ncid, varid                  ! NetCDF identifiers
   CHARACTER(LEN=30)               :: dnam       ! dimension name
   CHARACTER(LEN=80)               :: varname    ! NetCDF variable name
@@ -323,10 +326,9 @@
 !--- fields
   INTEGER :: imdep, jmdep, lmdep                ! dimensions of 'champ'
-  REAL, ALLOCATABLE, DIMENSION(:,:) :: champ      ! wanted field on initial grid
+  REAL, ALLOCATABLE, DIMENSION(:, :) :: champ   ! wanted field on initial grid
   REAL, ALLOCATABLE, DIMENSION(:) :: yder, timeyear
-  REAL,              DIMENSION(iim,jjp1) :: champint   ! interpolated field
-  REAL, ALLOCATABLE, DIMENSION(:,:,:)    :: champtime
-  REAL, ALLOCATABLE, DIMENSION(:,:,:)    :: champan
-  REAL                              :: time, by
+  REAL,              DIMENSION(iim, jjp1) :: champint   ! interpolated field
+  REAL, ALLOCATABLE, DIMENSION(:, :, :)    :: champtime
+  REAL, ALLOCATABLE, DIMENSION(:, :, :)    :: champan
 !--- input files
   CHARACTER(LEN=20)                 :: cal_in   ! calendar
@@ -334,10 +336,13 @@
 !--- misc
   INTEGER :: i, j, k, l                         ! loop counters
-  REAL, ALLOCATABLE, DIMENSION(:,:) :: work     ! used for extrapolation
+  REAL, ALLOCATABLE, DIMENSION(:, :) :: work     ! used for extrapolation
   CHARACTER(LEN=25)                 :: title    ! for messages
   LOGICAL                           :: extrp    ! flag for extrapolation
   REAL                              :: chmin, chmax
-!-------------------------------------------------------------------------------
-!---Variables depending on keyword 'mode' --------------------------------------
+  INTEGER ierr
+  integer n_extrap ! number of extrapolated points
+  logical skip
+!------------------------------------------------------------------------------
+!---Variables depending on keyword 'mode' -------------------------------------
   NULLIFY(champo)
   SELECT CASE(mode)
@@ -352,70 +357,75 @@
   END IF
 
-!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ------------------------------
-  ierr=NF90_OPEN(fnam,NF90_NOWRITE,ncid);                  CALL ncerr(ierr,fnam)
-  ierr=NF90_INQ_VARID(ncid,varname,varid);                 CALL ncerr(ierr,fnam)
-  ierr=NF90_INQUIRE_VARIABLE(ncid,varid,dimids=dids);      CALL ncerr(ierr,fnam)
+!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
+  ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid);             CALL ncerr(ierr, fnam)
+  ierr=NF90_INQ_VARID(ncid, varname, varid);            CALL ncerr(ierr, fnam)
+  ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam)
 
 !--- Longitude
-  ierr=NF90_INQUIRE_DIMENSION(ncid,dids(1),name=dnam,len=imdep)
-  CALL ncerr(ierr,fnam); ALLOCATE(dlon_ini(imdep),dlon(imdep))
-  ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
-  ierr=NF90_GET_VAR(ncid,varid,dlon_ini);                  CALL ncerr(ierr,fnam)
-  WRITE(lunout,*) 'variable ', dnam, 'dimension ', imdep
+  ierr=NF90_INQUIRE_DIMENSION(ncid, dids(1), name=dnam, len=imdep)
+  CALL ncerr(ierr, fnam); ALLOCATE(dlon_ini(imdep), dlon(imdep))
+  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
+  ierr=NF90_GET_VAR(ncid, varid, dlon_ini);              CALL ncerr(ierr, fnam)
+  WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep
 
 !--- Latitude
-  ierr=NF90_INQUIRE_DIMENSION(ncid,dids(2),name=dnam,len=jmdep)
-  CALL ncerr(ierr,fnam); ALLOCATE(dlat_ini(jmdep),dlat(jmdep))
-  ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
-  ierr=NF90_GET_VAR(ncid,varid,dlat_ini);                  CALL ncerr(ierr,fnam)
-  WRITE(lunout,*) 'variable ', dnam, 'dimension ', jmdep
+  ierr=NF90_INQUIRE_DIMENSION(ncid, dids(2), name=dnam, len=jmdep)
+  CALL ncerr(ierr, fnam); ALLOCATE(dlat_ini(jmdep), dlat(jmdep))
+  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
+  ierr=NF90_GET_VAR(ncid, varid, dlat_ini);              CALL ncerr(ierr, fnam)
+  WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep
 
 !--- Time (variable is not needed - it is rebuilt - but calendar is)
-  ierr=NF90_INQUIRE_DIMENSION(ncid,dids(3),name=dnam,len=lmdep)
-  CALL ncerr(ierr,fnam); ALLOCATE(timeyear(lmdep))
-  ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
+  ierr=NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep)
+  CALL ncerr(ierr, fnam); ALLOCATE(timeyear(lmdep))
+  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
   cal_in=' '
-  ierr=NF90_GET_ATT(ncid,varid,'calendar',cal_in)
+  ierr=NF90_GET_ATT(ncid, varid, 'calendar', cal_in)
   IF(ierr/=NF90_NOERR) THEN
     SELECT CASE(mode)
-      CASE('RUG','ALB'); cal_in='360d'
-      CASE('SIC','SST'); cal_in='gregorian'
+      CASE('RUG', 'ALB'); cal_in='360d'
+      CASE('SIC', 'SST'); cal_in='gregorian'
     END SELECT
-    WRITE(lunout,*)'ATTENTION: variable ''time'' sans attribut ''calendrier'' d&
-     &ans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
+    WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &
+         // 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
   END IF
-  WRITE(lunout,*) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', cal_in
-
-!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION ----------------------
+  WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &
+       cal_in
+
+!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION --------------------
   !--- Determining input file number of days, depending on calendar
-  ndays_in=year_len(anneeref,cal_in)
+  ndays_in=year_len(anneeref, cal_in)
 
 !--- Time vector reconstruction (time vector from file is not trusted)
 !--- If input records are not monthly, time sampling has to be constant !
-  timeyear=mid_months(anneeref,cal_in,lmdep)
-  IF(lmdep/=12) WRITE(lunout,'(a,i3,a)')'Note: les fichiers de '//TRIM(mode)   &
-    //' ne comportent pas 12, mais ',lmdep,' enregistrements.'
-
-!--- GETTING THE FIELD AND INTERPOLATING IT ------------------------------------
-  ALLOCATE(champ(imdep,jmdep),champtime(iim,jjp1,lmdep))
-  IF(extrp) ALLOCATE(work(imdep,jmdep))
-
-  WRITE(lunout,*)
-  WRITE(lunout,'(a,i3,a)')'LECTURE ET INTERPOLATION HORIZ. DE ',lmdep,' CHAMPS.'
-  ierr=NF90_INQ_VARID(ncid,varname,varid);                 CALL ncerr(ierr,fnam)
-  DO l=1,lmdep
-    ierr=NF90_GET_VAR(ncid,varid,champ,(/1,1,l/),(/imdep,jmdep,1/))
-    CALL ncerr(ierr,fnam)
-    CALL conf_dat2d(title,imdep,jmdep,dlon_ini,dlat_ini,dlon,dlat,champ,ibar)
-
-    IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
+  timeyear=mid_months(anneeref, cal_in, lmdep)
+  IF (lmdep /= 12) WRITE(lunout, '(a, i3, a)') 'Note : les fichiers de ' &
+       // TRIM(mode) // ' ne comportent pas 12, mais ', lmdep, &
+       ' enregistrements.'
+
+!--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
+  ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep))
+  IF(extrp) ALLOCATE(work(imdep, jmdep))
+
+  WRITE(lunout, *)
+  WRITE(lunout, '(a, i3, a)')'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, &
+       ' CHAMPS.'
+  ierr=NF90_INQ_VARID(ncid, varname, varid);             CALL ncerr(ierr, fnam)
+  DO l=1, lmdep
+    ierr=NF90_GET_VAR(ncid, varid, champ, (/1, 1, l/), (/imdep, jmdep, 1/))
+    CALL ncerr(ierr, fnam)
+    CALL conf_dat2d(title, imdep, jmdep, dlon_ini, dlat_ini, dlon, dlat, &
+         champ, ibar)
+
+    IF (extrp) CALL extrapol(champ, imdep, jmdep, 999999., .TRUE., .TRUE., 2, &
+         work)
 
     IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN
       IF(l==1) THEN
-        WRITE(lunout,*)                                                        &
+        WRITE(lunout, *)                                                      &
   '-------------------------------------------------------------------------'
-        WRITE(lunout,*)                                                        &
-  '$$$ Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$'
-        WRITE(lunout,*)                                                        &
+        WRITE(lunout, *)                                                     &
+  'Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$'
+        WRITE(lunout, *)                                                      &
   '-------------------------------------------------------------------------'
       END IF
@@ -433,59 +443,64 @@
         CASE('SIC');       CALL sea_ice (imdep, jmdep, dlon, dlat, champ,    &
                                     iim, jjp1, rlonv, rlatu, champint)
-        CASE('SST','ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ,    &
+        CASE('SST', 'ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ,    &
                                     iim, jjp1, rlonv, rlatu, champint)
       END SELECT
     END IF
-    champtime(:,:,l)=champint
+    champtime(:, :, l)=champint
   END DO
-  ierr=NF90_CLOSE(ncid);                                   CALL ncerr(ierr,fnam)
-
-  DEALLOCATE(dlon_ini,dlat_ini,dlon,dlat,champ)
+  ierr=NF90_CLOSE(ncid);                                 CALL ncerr(ierr, fnam)
+
+  DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
   IF(extrp) DEALLOCATE(work)
 
-!--- TIME INTERPOLATION --------------------------------------------------------
-  WRITE(lunout,*)
-  WRITE(lunout,*)'INTERPOLATION TEMPORELLE.'
-  WRITE(lunout,"(2x,' Vecteur temps en entree: ',10f6.1)") timeyear
-  WRITE(lunout,"(2x,' Vecteur temps en sortie de 0 a ',i3)") ndays
-  ALLOCATE(yder(lmdep),champan(iip1,jjp1,ndays))
-  DO j=1,jjp1
-    DO i=1,iim
-      CALL spline(timeyear,champtime(i,j,:),lmdep,1.e30,1.e30,yder)
-      DO k=1,ndays
-        time=FLOAT((k-1)*ndays_in)/FLOAT(ndays)
-        CALL splint(timeyear,champtime(i,j,:),yder,lmdep,time,by)
-        champan(i,j,k) = by
-      END DO
+!--- TIME INTERPOLATION ------------------------------------------------------
+  WRITE(lunout, *)
+  WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
+  WRITE(lunout, "(2x, ' Vecteur temps en entree: ', 10f6.1)") timeyear
+  WRITE(lunout, "(2x, ' Vecteur temps en sortie de 0 a ', i3)") ndays
+  ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays))
+  skip = .false.
+  n_extrap = 0
+  DO j=1, jjp1
+    DO i=1, iim
+      yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, &
+           vc_beg=0., vc_end=0.)
+      CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, &
+           arth(0., real(ndays_in) / ndays, ndays), champan(i, j, :), ierr)
+      if (ierr < 0) stop 1
+      n_extrap = n_extrap + ierr
     END DO
   END DO
-  champan(iip1,:,:)=champan(1,:,:)
-  DEALLOCATE(yder,champtime,timeyear)
+  if (n_extrap /= 0) then
+     print *, "get_2Dfield pchfe_95: n_extrap = ", n_extrap
+  end if
+  champan(iip1, :, :)=champan(1, :, :)
+  DEALLOCATE(yder, champtime, timeyear)
 
 !--- Checking the result
-  DO j=1,jjp1
-    CALL minmax(iip1,champan(1,j,10),chmin,chmax)
-    WRITE(lunout,*)' '//TRIM(title)//' au temps 10 ',chmin,chmax,j
+  DO j=1, jjp1
+    CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
+    WRITE(lunout, *)' '//TRIM(title)//' au temps 10 ', chmin, chmax, j
   END DO
 
-!--- SPECIAL FILTER FOR SST: SST>271.38 ----------------------------------------
+!--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
   IF(mode=='SST') THEN
-    WRITE(lunout,*) 'Filtrage de la SST: SST >= 271.38'
+    WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'
     WHERE(champan<271.38) champan=271.38
   END IF
 
-!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 ---------------------------------------
+!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
   IF(mode=='SIC') THEN
-    WRITE(lunout,*) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
-    IF(.NOT.lCPL) champan(:,:,:)=champan(:,:,:)/100.
-    champan(iip1,:,:)=champan(1,:,:)
+    WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
+    IF(.NOT.lCPL) champan(:, :, :)=champan(:, :, :)/100.
+    champan(iip1, :, :)=champan(1, :, :)
     WHERE(champan>1.0) champan=1.0
     WHERE(champan<0.0) champan=0.0
   END IF
 
-!--- DYNAMICAL TO PHYSICAL GRID ------------------------------------------------
-  ALLOCATE(champo(klon,ndays))
-  DO k=1,ndays
-    CALL gr_dyn_fi(1,iip1,jjp1,klon,champan(1,1,k),champo(1,k))
+!--- DYNAMICAL TO PHYSICAL GRID ----------------------------------------------
+  ALLOCATE(champo(klon, ndays))
+  DO k=1, ndays
+    CALL gr_dyn_fi(1, iip1, jjp1, klon, champan(1, 1, k), champo(1, k))
   END DO
   DEALLOCATE(champan)
Index: DZ4/trunk/libf/dyn3d/spline.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3d/spline.F	(revision 1424)
+++ 	(revision )
@@ -1,79 +1,0 @@
-!
-! $Header$
-!
-      subroutine spline(x,y,n,yp1,ypn,y2)
-     
-c
-     
-c     Routine to set up the interpolating function for a cubic spline
-     
-c     interpolation (see "Numerical Recipes" for details).
-     
-c
-	  implicit real (a-h,o-z)
-	  implicit integer (i-n)
-     
-      parameter(nllm=4096)
-     
-      dimension x(n),y(n),y2(n),u(nllm)
-     
-c
-c	write(6,*)(x(i),i=1,n)
-c	write(6,*)(y(i),i=1,n)
-     
-      if(yp1.gt.0.99E30) then
-c the lower boundary condition is set
-       y2(1)=0.
-c either to be "natural"
-       u(1)=0.
-     
-      else
-c or else to have a specified first
-       y2(1)=-0.5
-c derivative
-       u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
-     
-      end if
-     
-      do 11 i=2,n-1
-c decomposition loop of the tridiagonal
-       sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
-c algorithm. Y2 and U are used
-       p=sig*y2(i-1)+2.
-c for temporary storage of the decompo-
-       y2(i)=(sig-1.)/p
-c sed factors
-       u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))
-     
-     . /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
-     
- 11   continue
-     
-      if(ypn.gt.0.99E30) then
-c the upper boundary condition is set
-       qn=0.
-c either to be "natural"
-       un=0.
-     
-      else
-c or else to have a specified first
-       qn=0.5
-c derivative
-       un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
-     
-      end if
-     
-      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
-     
-      do 12 k=n-1,1,-1
-c this is the backsubstitution loop of
-       y2(k)=y2(k)*y2(k+1)+u(k)
-c the tridiagonal algorithm
- 12   continue
-     
-c
-     
-      return
-     
-      end
-     
Index: DZ4/trunk/libf/dyn3d/splint.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3d/splint.F	(revision 1424)
+++ 	(revision )
@@ -1,56 +1,0 @@
-!
-! $Header$
-!
-     
-      SUBROUTINE splint(xa,ya,y2a,n,x,y)
-     
-c
-c     Routine to compute a cubic-spline interpolated value Y given the
-c     value of X, the arrays XA, YA and the 2nd derivative array Y2A
-c     computed by SUBROUTINE SPLINE. See "Numerical Recipes" for details
-c
-     
-      IMPLICIT REAL (a-h,o-z)
-      IMPLICIT INTEGER (i-n)
-      DIMENSION xa(n),ya(n),y2a(n)
-     
-      kl0=1
-     
-      khi=n
-c means of bisection
- 1    IF(khi-kl0.gt.1) THEN
-     
-       k=(khi+kl0)/2
-     
-       IF(xa(k).gt.x) THEN
-     
-        khi=k
-     
-       ELSE
-     
-        kl0=k
-     
-       END IF
-     
-       GO TO 1
-     
-      END IF
-c KL0 and KHI now bracket the X
-      h=xa(khi)-xa(kl0)
-     
-      IF(h.eq.0.0) STOP
-      a=(xa(khi)-x)/h
-c evaluation of cubic spline polynomial
-      b=(x-xa(kl0))/h
-     
-      y=a*ya(kl0)+b*ya(khi)+((a**3-a)*y2a(kl0)+(b**3-b)*y2a(khi))*(h**2)
-     
-     ./6.
-     
-c
-     
-      RETURN
-     
-      END
-     
-
Index: /LMDZ4/trunk/libf/dyn3d/startvar.F90
===================================================================
--- /LMDZ4/trunk/libf/dyn3d/startvar.F90	(revision 1424)
+++ /LMDZ4/trunk/libf/dyn3d/startvar.F90	(revision 1425)
@@ -643,8 +643,10 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2,lon_in2,&
-                          lat_in2, pls_in, var3d, ibar)
-!
-!-------------------------------------------------------------------------------
+SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2, &
+     lon_in2, lat_in2, pls_in, var3d, ibar)
+
+  use pchsp_95_m, only: pchsp_95
+  use pchfe_95_m, only: pchfe_95
+
 ! Arguments:
   CHARACTER(LEN=*),             INTENT(IN)    :: varname
@@ -655,73 +657,71 @@
   REAL, DIMENSION(iml),         INTENT(IN)    :: lon_in2
   REAL, DIMENSION(jml2),        INTENT(IN)    :: lat_in2
-  REAL, DIMENSION(iml,jml,lml), INTENT(IN)    :: pls_in
-  REAL, DIMENSION(iml,jml,lml), INTENT(OUT)   :: var3d
+  REAL, DIMENSION(iml, jml, lml), INTENT(IN)    :: pls_in
+  REAL, DIMENSION(iml, jml, lml), INTENT(OUT)   :: var3d
   LOGICAL,                      INTENT(IN)    :: ibar
-!-------------------------------------------------------------------------------
+!----------------------------------------------------------------------------
 ! Local variables:
 #include "iniprint.h"
-  LOGICAL               :: check=.TRUE.
-  REAL                  :: bx, by, chmin, chmax
-  INTEGER               :: ii, ij, il
-  REAL,    DIMENSION(:,:,:), ALLOCATABLE :: var_tmp3d
+  LOGICAL:: check=.TRUE., skip
+  REAL                  chmin, chmax
+  INTEGER ii, ij, il, ierr
+  integer n_extrap ! number of extrapolated points
+  REAL, DIMENSION(iml, jml, llm_dyn):: var_tmp3d
   REAL,    DIMENSION(:),     ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
-  REAL,    DIMENSION(:),     ALLOCATABLE :: lev_dyn, ax, ay, yder
-  INTEGER, DIMENSION(:),     ALLOCATABLE :: lind
-!-------------------------------------------------------------------------------
-  IF(check) WRITE(lunout,*)'Going into flinget to extract the 3D  field.'
-  IF(check) WRITE(lunout,*)fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn
-  IF(check) WRITE(lunout,*)'Allocating space for interpolation',iml,jml,llm_dyn
-
-  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn,jml_dyn,llm_dyn))
-  CALL flinget(fid_dyn,varname,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d)
+  REAL, DIMENSION(llm_dyn):: lev_dyn, ax, ay, yder
+
+!---------------------------------------------------------------------------
+  IF(check) WRITE(lunout, *)'Going into flinget to extract the 3D  field.'
+  IF(check) WRITE(lunout, *) fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, &
+       ttm_dyn
+  IF(check) WRITE(lunout, *) 'Allocating space for interpolation', iml, jml, &
+       llm_dyn
+
+  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
+  CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, 1, 1, &
+       var_ana3d)
 
 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
-  ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
-  lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
-  lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
+  ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn))
+  lon_ini(:)=lon_dyn(:, 1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
+  lat_ini(:)=lat_dyn(1, :); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
 
 !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
-  ALLOCATE(lon_rad(iml_dyn),lat_rad(jml_dyn),lev_dyn(llm_dyn))
-  CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini,       &
+  ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn))
+  CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini,      &
                    levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
-  DEALLOCATE(lon_ini,lat_ini)
+  DEALLOCATE(lon_ini, lat_ini)
 
 !--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
-  ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
-  DO il=1,llm_dyn
-    CALL interp_startvar(varname, ibar, il==1,                                 &
-      iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana3d(:,:,il), iml, jml, jml2,   &
-      lon_in,  lat_in,  lon_in2, lat_in2, var_tmp3d(:,:,il))
+  DO il=1, llm_dyn
+    CALL interp_startvar(varname, ibar, il==1, iml_dyn, jml_dyn, lon_rad, &
+         lat_rad, var_ana3d(:, :, il), iml, jml, jml2, lon_in, lat_in, &
+         lon_in2, lat_in2, var_tmp3d(:, :, il))
   END DO
-  DEALLOCATE(lon_rad,lat_rad)
-
-  ALLOCATE(lind(llm_dyn))
-  DO il=1,llm_dyn
-    lind(il) = llm_dyn-il+1
+  DEALLOCATE(lon_rad, lat_rad)
+
+!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
+  ax = lev_dyn(llm_dyn:1:-1) 
+  skip = .false.
+  n_extrap = 0
+  DO ij=1, jml
+    DO ii=1, iml-1
+      ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
+      yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.)
+      CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), &
+           var3d(ii, ij, lml:1:-1), ierr)
+      if (ierr < 0) stop 1
+      n_extrap = n_extrap + ierr
+    END DO
   END DO
-
-!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
-  ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn))
-  DO ij=1,jml
-    DO ii=1,iml-1
-      ax(:)=lev_dyn(lind(:)) 
-      ay(:)=var_tmp3d(ii,ij,lind(:))
-      CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder)
-      DO il=1,lml
-        bx=pls_in(ii,ij,il)
-        CALL SPLINT(ax, ay, yder, llm_dyn, bx, by)
-        var3d(ii,ij,il) = by
-      END DO
-    END DO
-    var3d(iml,ij,:) = var3d(1,ij,:) 
+  if (n_extrap /= 0) then
+     print *, "start_inter_3d pchfe_95: n_extrap = ", n_extrap
+  end if
+  var3d(iml, :, :) = var3d(1, :, :) 
+
+  DO il=1, lml
+    CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax)
+    WRITE(lunout, *)' '//TRIM(varname)//'  min max l ', il, chmin, chmax
   END DO
-  DEALLOCATE(var_tmp3d,lev_dyn,ax,ay,yder,lind)
-
-  DO il=1,lml
-    CALL minmax(iml*jml,var3d(1,1,il),chmin,chmax)
-    WRITE(lunout,*)' '//TRIM(varname)//'  min max l ',il,chmin,chmax
-  END DO
-
-  RETURN
 
 END SUBROUTINE start_inter_3d
Index: /LMDZ4/trunk/libf/dyn3dpar/abort_gcm.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/abort_gcm.F	(revision 1424)
+++ /LMDZ4/trunk/libf/dyn3dpar/abort_gcm.F	(revision 1425)
@@ -23,7 +23,7 @@
 C         ierr    = severity of situation ( = 0 normal )
 
-      character (len=*) :: modname
+      character(len=*) modname
       integer ierr
-      character (len=*) :: message
+      character(len=*) message
 
       write(lunout,*) 'in abort_gcm'
@@ -45,7 +45,8 @@
       if (ierr .eq. 0) then
         write(lunout,*) 'Everything is cool'
+        stop
       else
         write(lunout,*) 'Houston, we have a problem ', ierr
-      STOP
+        stop 1
       endif
       END
Index: /LMDZ4/trunk/libf/dyn3dpar/ce0l.F90
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/ce0l.F90	(revision 1424)
+++ /LMDZ4/trunk/libf/dyn3dpar/ce0l.F90	(revision 1425)
@@ -5,5 +5,4 @@
 !
 PROGRAM ce0l
-!
 !-------------------------------------------------------------------------------
 ! Purpose: Calls etat0, creates initial states and limit_netcdf
@@ -104,5 +103,4 @@
 #endif
 ! of #ifndef CPP_EARTH #else
-  STOP
 
 END PROGRAM ce0l
Index: /LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F90
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F90	(revision 1424)
+++ /LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F90	(revision 1425)
@@ -11,4 +11,5 @@
 ! Note: This routine is designed to work for Earth
 !-------------------------------------------------------------------------------
+  USE control_mod
 #ifdef CPP_EARTH
   USE startvar
@@ -28,5 +29,4 @@
   USE netcdf, ONLY : NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR
 #endif
-  USE control_mod
   IMPLICIT NONE
 !-------------------------------------------------------------------------------
@@ -114,5 +114,5 @@
 !-------------------------------------------------------------------------------
   REAL    :: alp_offset
-  LOGICAL found
+  logical found
 
 !--- Constants
@@ -505,4 +505,9 @@
   wake_cstar(:) = 0.
   wake_fip(:) = 0.
+  wake_pe = 0.
+  fm_therm = 0.
+  entr_therm = 0.
+  detr_therm = 0.
+
   CALL fonte_neige_init(run_off_lic_0)
   CALL pbl_surface_init( qsol, fder, snsrf, qsolsrf, evap, frugs, agesno, tsoil )
Index: /LMDZ4/trunk/libf/dyn3dpar/limit_netcdf.F90
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/limit_netcdf.F90	(revision 1424)
+++ /LMDZ4/trunk/libf/dyn3dpar/limit_netcdf.F90	(revision 1425)
@@ -20,4 +20,5 @@
 !  *    12/2009: D. Cugnet   (f77->f90, calendars, files from coupled runs)
 !-------------------------------------------------------------------------------
+  USE control_mod
 #ifdef CPP_EARTH
   USE dimphy
@@ -27,8 +28,7 @@
                    NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT,     &
                    NF90_NOERR,   NF90_NOWRITE, NF90_DOUBLE,  NF90_GLOBAL,      &
-		   NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED
+		   NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED, NF90_FLOAT
   USE inter_barxy_m, only: inter_barxy
 #endif
-  USE control_mod
   IMPLICIT NONE
 !-------------------------------------------------------------------------------
@@ -267,5 +267,4 @@
 #endif
 ! of #ifdef CPP_EARTH
-  STOP
 
 
@@ -281,5 +280,5 @@
 SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask, lCPL)
 !
-!-------------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
 ! Comments:
 !   There are two assumptions concerning the NetCDF files, that are satisfied
@@ -287,11 +286,15 @@
 !     1) The last dimension of the variables used is the time record.
 !     2) Dimensional variables have the same names as corresponding dimensions.
-!-------------------------------------------------------------------------------
-  USE netcdf, ONLY : NF90_OPEN,    NF90_INQ_VARID, NF90_INQUIRE_VARIABLE,      &
-                     NF90_CLOSE,   NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION,     &
-                     NF90_GET_VAR, NF90_GET_ATT
+!-----------------------------------------------------------------------------
+  USE netcdf, ONLY: NF90_OPEN, NF90_INQ_VARID, NF90_INQUIRE_VARIABLE, &
+       NF90_CLOSE, NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, NF90_GET_VAR, &
+       NF90_GET_ATT
   USE dimphy, ONLY : klon
   USE phys_state_var_mod, ONLY : pctsrf
   USE control_mod
+  use pchsp_95_m, only: pchsp_95
+  use pchfe_95_m, only: pchfe_95
+  use arth_m, only: arth
+
   IMPLICIT NONE
 #include "dimensions.h"
@@ -300,5 +303,5 @@
 #include "indicesol.h"
 #include "iniprint.h"
-!-------------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
 ! Arguments:
   CHARACTER(LEN=*),  INTENT(IN)     :: fnam     ! NetCDF file name
@@ -306,12 +309,12 @@
   LOGICAL,           INTENT(IN)     :: ibar     ! interp on pressure levels
   INTEGER,           INTENT(IN)     :: ndays    ! current year number of days
-  REAL,    POINTER,  DIMENSION(:,:) :: champo   ! output field = f(t)
-  LOGICAL, OPTIONAL, INTENT(IN)     :: flag     ! extrapol. (SST)  old ice (SIC)
-  REAL,    OPTIONAL, DIMENSION(iim,jjp1), INTENT(IN) :: mask
+  REAL,    POINTER,  DIMENSION(:, :) :: champo   ! output field = f(t)
+  LOGICAL, OPTIONAL, INTENT(IN)     :: flag     ! extrapol. (SST) old ice (SIC)
+  REAL,    OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask
   LOGICAL, OPTIONAL, INTENT(IN)     :: lCPL     ! Coupled model flag (for ICE)
-!-------------------------------------------------------------------------------
+!------------------------------------------------------------------------------
 ! Local variables:
 !--- NetCDF
-  INTEGER :: ierr, ncid, varid                  ! NetCDF identifiers
+  INTEGER :: ncid, varid                  ! NetCDF identifiers
   CHARACTER(LEN=30)               :: dnam       ! dimension name
   CHARACTER(LEN=80)               :: varname    ! NetCDF variable name
@@ -323,10 +326,9 @@
 !--- fields
   INTEGER :: imdep, jmdep, lmdep                ! dimensions of 'champ'
-  REAL, ALLOCATABLE, DIMENSION(:,:) :: champ      ! wanted field on initial grid
+  REAL, ALLOCATABLE, DIMENSION(:, :) :: champ   ! wanted field on initial grid
   REAL, ALLOCATABLE, DIMENSION(:) :: yder, timeyear
-  REAL,              DIMENSION(iim,jjp1) :: champint   ! interpolated field
-  REAL, ALLOCATABLE, DIMENSION(:,:,:)    :: champtime
-  REAL, ALLOCATABLE, DIMENSION(:,:,:)    :: champan
-  REAL                              :: time, by
+  REAL,              DIMENSION(iim, jjp1) :: champint   ! interpolated field
+  REAL, ALLOCATABLE, DIMENSION(:, :, :)    :: champtime
+  REAL, ALLOCATABLE, DIMENSION(:, :, :)    :: champan
 !--- input files
   CHARACTER(LEN=20)                 :: cal_in   ! calendar
@@ -334,10 +336,13 @@
 !--- misc
   INTEGER :: i, j, k, l                         ! loop counters
-  REAL, ALLOCATABLE, DIMENSION(:,:) :: work     ! used for extrapolation
+  REAL, ALLOCATABLE, DIMENSION(:, :) :: work     ! used for extrapolation
   CHARACTER(LEN=25)                 :: title    ! for messages
   LOGICAL                           :: extrp    ! flag for extrapolation
   REAL                              :: chmin, chmax
-!-------------------------------------------------------------------------------
-!---Variables depending on keyword 'mode' --------------------------------------
+  INTEGER ierr
+  integer n_extrap ! number of extrapolated points
+  logical skip
+!------------------------------------------------------------------------------
+!---Variables depending on keyword 'mode' -------------------------------------
   NULLIFY(champo)
   SELECT CASE(mode)
@@ -347,76 +352,80 @@
     CASE('ALB'); varname='ALBEDO'; title='Albedo'
   END SELECT
-  extrp=.FALSE.
-  IF ( PRESENT(flag) ) THEN
-    IF ( flag .AND. mode=='SST' ) extrp=.TRUE.  
-  END IF  
-
-
-!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ------------------------------
-  ierr=NF90_OPEN(fnam,NF90_NOWRITE,ncid);                  CALL ncerr(ierr,fnam)
-  ierr=NF90_INQ_VARID(ncid,varname,varid);                 CALL ncerr(ierr,fnam)
-  ierr=NF90_INQUIRE_VARIABLE(ncid,varid,dimids=dids);      CALL ncerr(ierr,fnam)
+  extrp=.FALSE. 
+  IF ( PRESENT(flag) ) THEN 
+    IF ( flag .AND. mode=='SST' ) extrp=.TRUE. 
+  END IF
+
+!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
+  ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid);             CALL ncerr(ierr, fnam)
+  ierr=NF90_INQ_VARID(ncid, varname, varid);            CALL ncerr(ierr, fnam)
+  ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam)
 
 !--- Longitude
-  ierr=NF90_INQUIRE_DIMENSION(ncid,dids(1),name=dnam,len=imdep)
-  CALL ncerr(ierr,fnam); ALLOCATE(dlon_ini(imdep),dlon(imdep))
-  ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
-  ierr=NF90_GET_VAR(ncid,varid,dlon_ini);                  CALL ncerr(ierr,fnam)
-  WRITE(lunout,*) 'variable ', dnam, 'dimension ', imdep
+  ierr=NF90_INQUIRE_DIMENSION(ncid, dids(1), name=dnam, len=imdep)
+  CALL ncerr(ierr, fnam); ALLOCATE(dlon_ini(imdep), dlon(imdep))
+  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
+  ierr=NF90_GET_VAR(ncid, varid, dlon_ini);              CALL ncerr(ierr, fnam)
+  WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep
 
 !--- Latitude
-  ierr=NF90_INQUIRE_DIMENSION(ncid,dids(2),name=dnam,len=jmdep)
-  CALL ncerr(ierr,fnam); ALLOCATE(dlat_ini(jmdep),dlat(jmdep))
-  ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
-  ierr=NF90_GET_VAR(ncid,varid,dlat_ini);                  CALL ncerr(ierr,fnam)
-  WRITE(lunout,*) 'variable ', dnam, 'dimension ', jmdep
+  ierr=NF90_INQUIRE_DIMENSION(ncid, dids(2), name=dnam, len=jmdep)
+  CALL ncerr(ierr, fnam); ALLOCATE(dlat_ini(jmdep), dlat(jmdep))
+  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
+  ierr=NF90_GET_VAR(ncid, varid, dlat_ini);              CALL ncerr(ierr, fnam)
+  WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep
 
 !--- Time (variable is not needed - it is rebuilt - but calendar is)
-  ierr=NF90_INQUIRE_DIMENSION(ncid,dids(3),name=dnam,len=lmdep)
-  CALL ncerr(ierr,fnam); ALLOCATE(timeyear(lmdep))
-  ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
+  ierr=NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep)
+  CALL ncerr(ierr, fnam); ALLOCATE(timeyear(lmdep))
+  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
   cal_in=' '
-  ierr=NF90_GET_ATT(ncid,varid,'calendar',cal_in)
+  ierr=NF90_GET_ATT(ncid, varid, 'calendar', cal_in)
   IF(ierr/=NF90_NOERR) THEN
     SELECT CASE(mode)
-      CASE('RUG','ALB'); cal_in='360d'
-      CASE('SIC','SST'); cal_in='gregorian'
+      CASE('RUG', 'ALB'); cal_in='360d'
+      CASE('SIC', 'SST'); cal_in='gregorian'
     END SELECT
-    WRITE(lunout,*)'ATTENTION: variable ''time'' sans attribut ''calendrier'' d&
-     &ans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
+    WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &
+         // 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
   END IF
-  WRITE(lunout,*) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', cal_in
-
-!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION ----------------------
+  WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &
+       cal_in
+
+!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION --------------------
   !--- Determining input file number of days, depending on calendar
-  ndays_in=year_len(anneeref,cal_in)
+  ndays_in=year_len(anneeref, cal_in)
 
 !--- Time vector reconstruction (time vector from file is not trusted)
 !--- If input records are not monthly, time sampling has to be constant !
-  timeyear=mid_months(anneeref,cal_in,lmdep)
-  IF(lmdep/=12) WRITE(lunout,'(a,i3,a)')'Note: les fichiers de '//TRIM(mode)   &
-    //' ne comportent pas 12, mais ',lmdep,' enregistrements.'
-
-!--- GETTING THE FIELD AND INTERPOLATING IT ------------------------------------
-  ALLOCATE(champ(imdep,jmdep),champtime(iim,jjp1,lmdep))
-  IF(extrp) ALLOCATE(work(imdep,jmdep))
-
-  WRITE(lunout,*)
-  WRITE(lunout,'(a,i3,a)')'LECTURE ET INTERPOLATION HORIZ. DE ',lmdep,' CHAMPS.'
-  ierr=NF90_INQ_VARID(ncid,varname,varid);                 CALL ncerr(ierr,fnam)
-  DO l=1,lmdep
-    ierr=NF90_GET_VAR(ncid,varid,champ,(/1,1,l/),(/imdep,jmdep,1/))
-    CALL ncerr(ierr,fnam)
-    CALL conf_dat2d(title,imdep,jmdep,dlon_ini,dlat_ini,dlon,dlat,champ,ibar)
-
-    IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
+  timeyear=mid_months(anneeref, cal_in, lmdep)
+  IF (lmdep /= 12) WRITE(lunout, '(a, i3, a)') 'Note : les fichiers de ' &
+       // TRIM(mode) // ' ne comportent pas 12, mais ', lmdep, &
+       ' enregistrements.'
+
+!--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
+  ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep))
+  IF(extrp) ALLOCATE(work(imdep, jmdep))
+
+  WRITE(lunout, *)
+  WRITE(lunout, '(a, i3, a)')'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, &
+       ' CHAMPS.'
+  ierr=NF90_INQ_VARID(ncid, varname, varid);             CALL ncerr(ierr, fnam)
+  DO l=1, lmdep
+    ierr=NF90_GET_VAR(ncid, varid, champ, (/1, 1, l/), (/imdep, jmdep, 1/))
+    CALL ncerr(ierr, fnam)
+    CALL conf_dat2d(title, imdep, jmdep, dlon_ini, dlat_ini, dlon, dlat, &
+         champ, ibar)
+
+    IF (extrp) CALL extrapol(champ, imdep, jmdep, 999999., .TRUE., .TRUE., 2, &
+         work)
 
     IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN
       IF(l==1) THEN
-        WRITE(lunout,*)                                                        &
+        WRITE(lunout, *)                                                      &
   '-------------------------------------------------------------------------'
-        WRITE(lunout,*)                                                        &
-  '$$$ Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$'
-        WRITE(lunout,*)                                                        &
+        WRITE(lunout, *)                                                     &
+  'Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$'
+        WRITE(lunout, *)                                                      &
   '-------------------------------------------------------------------------'
       END IF
@@ -434,59 +443,64 @@
         CASE('SIC');       CALL sea_ice (imdep, jmdep, dlon, dlat, champ,    &
                                     iim, jjp1, rlonv, rlatu, champint)
-        CASE('SST','ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ,    &
+        CASE('SST', 'ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ,    &
                                     iim, jjp1, rlonv, rlatu, champint)
       END SELECT
     END IF
-    champtime(:,:,l)=champint
+    champtime(:, :, l)=champint
   END DO
-  ierr=NF90_CLOSE(ncid);                                   CALL ncerr(ierr,fnam)
-
-  DEALLOCATE(dlon_ini,dlat_ini,dlon,dlat,champ)
+  ierr=NF90_CLOSE(ncid);                                 CALL ncerr(ierr, fnam)
+
+  DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
   IF(extrp) DEALLOCATE(work)
 
-!--- TIME INTERPOLATION --------------------------------------------------------
-  WRITE(lunout,*)
-  WRITE(lunout,*)'INTERPOLATION TEMPORELLE.'
-  WRITE(lunout,"(2x,' Vecteur temps en entree: ',10f6.1)") timeyear
-  WRITE(lunout,"(2x,' Vecteur temps en sortie de 0 a ',i3)") ndays
-  ALLOCATE(yder(lmdep),champan(iip1,jjp1,ndays))
-  DO j=1,jjp1
-    DO i=1,iim
-      CALL spline(timeyear,champtime(i,j,:),lmdep,1.e30,1.e30,yder)
-      DO k=1,ndays
-        time=FLOAT((k-1)*ndays_in)/FLOAT(ndays)
-        CALL splint(timeyear,champtime(i,j,:),yder,lmdep,time,by)
-        champan(i,j,k) = by
-      END DO
+!--- TIME INTERPOLATION ------------------------------------------------------
+  WRITE(lunout, *)
+  WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
+  WRITE(lunout, "(2x, ' Vecteur temps en entree: ', 10f6.1)") timeyear
+  WRITE(lunout, "(2x, ' Vecteur temps en sortie de 0 a ', i3)") ndays
+  ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays))
+  skip = .false.
+  n_extrap = 0
+  DO j=1, jjp1
+    DO i=1, iim
+      yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, &
+           vc_beg=0., vc_end=0.)
+      CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, &
+           arth(0., real(ndays_in) / ndays, ndays), champan(i, j, :), ierr)
+      if (ierr < 0) stop 1
+      n_extrap = n_extrap + ierr
     END DO
   END DO
-  champan(iip1,:,:)=champan(1,:,:)
-  DEALLOCATE(yder,champtime,timeyear)
+  if (n_extrap /= 0) then
+     print *, "get_2Dfield pchfe_95: n_extrap = ", n_extrap
+  end if
+  champan(iip1, :, :)=champan(1, :, :)
+  DEALLOCATE(yder, champtime, timeyear)
 
 !--- Checking the result
-  DO j=1,jjp1
-    CALL minmax(iip1,champan(1,j,10),chmin,chmax)
-    WRITE(lunout,*)' '//TRIM(title)//' au temps 10 ',chmin,chmax,j
+  DO j=1, jjp1
+    CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
+    WRITE(lunout, *)' '//TRIM(title)//' au temps 10 ', chmin, chmax, j
   END DO
 
-!--- SPECIAL FILTER FOR SST: SST>271.38 ----------------------------------------
+!--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
   IF(mode=='SST') THEN
-    WRITE(lunout,*) 'Filtrage de la SST: SST >= 271.38'
+    WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'
     WHERE(champan<271.38) champan=271.38
   END IF
 
-!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 ---------------------------------------
+!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
   IF(mode=='SIC') THEN
-    WRITE(lunout,*) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
-    IF(.NOT.lCPL) champan(:,:,:)=champan(:,:,:)/100.
-    champan(iip1,:,:)=champan(1,:,:)
+    WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
+    IF(.NOT.lCPL) champan(:, :, :)=champan(:, :, :)/100.
+    champan(iip1, :, :)=champan(1, :, :)
     WHERE(champan>1.0) champan=1.0
     WHERE(champan<0.0) champan=0.0
   END IF
 
-!--- DYNAMICAL TO PHYSICAL GRID ------------------------------------------------
-  ALLOCATE(champo(klon,ndays))
-  DO k=1,ndays
-    CALL gr_dyn_fi(1,iip1,jjp1,klon,champan(1,1,k),champo(1,k))
+!--- DYNAMICAL TO PHYSICAL GRID ----------------------------------------------
+  ALLOCATE(champo(klon, ndays))
+  DO k=1, ndays
+    CALL gr_dyn_fi(1, iip1, jjp1, klon, champan(1, 1, k), champo(1, k))
   END DO
   DEALLOCATE(champan)
Index: DZ4/trunk/libf/dyn3dpar/spline.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/spline.F	(revision 1424)
+++ 	(revision )
@@ -1,79 +1,0 @@
-!
-! $Header$
-!
-      subroutine spline(x,y,n,yp1,ypn,y2)
-     
-c
-     
-c     Routine to set up the interpolating function for a cubic spline
-     
-c     interpolation (see "Numerical Recipes" for details).
-     
-c
-	  implicit real (a-h,o-z)
-	  implicit integer (i-n)
-     
-      parameter(nllm=4096)
-     
-      dimension x(n),y(n),y2(n),u(nllm)
-     
-c
-c	write(6,*)(x(i),i=1,n)
-c	write(6,*)(y(i),i=1,n)
-     
-      if(yp1.gt.0.99E30) then
-c the lower boundary condition is set
-       y2(1)=0.
-c either to be "natural"
-       u(1)=0.
-     
-      else
-c or else to have a specified first
-       y2(1)=-0.5
-c derivative
-       u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
-     
-      end if
-     
-      do 11 i=2,n-1
-c decomposition loop of the tridiagonal
-       sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
-c algorithm. Y2 and U are used
-       p=sig*y2(i-1)+2.
-c for temporary storage of the decompo-
-       y2(i)=(sig-1.)/p
-c sed factors
-       u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))
-     
-     . /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
-     
- 11   continue
-     
-      if(ypn.gt.0.99E30) then
-c the upper boundary condition is set
-       qn=0.
-c either to be "natural"
-       un=0.
-     
-      else
-c or else to have a specified first
-       qn=0.5
-c derivative
-       un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
-     
-      end if
-     
-      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
-     
-      do 12 k=n-1,1,-1
-c this is the backsubstitution loop of
-       y2(k)=y2(k)*y2(k+1)+u(k)
-c the tridiagonal algorithm
- 12   continue
-     
-c
-     
-      return
-     
-      end
-     
Index: DZ4/trunk/libf/dyn3dpar/splint.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/splint.F	(revision 1424)
+++ 	(revision )
@@ -1,56 +1,0 @@
-!
-! $Header$
-!
-     
-      SUBROUTINE splint(xa,ya,y2a,n,x,y)
-     
-c
-c     Routine to compute a cubic-spline interpolated value Y given the
-c     value of X, the arrays XA, YA and the 2nd derivative array Y2A
-c     computed by SUBROUTINE SPLINE. See "Numerical Recipes" for details
-c
-     
-      IMPLICIT REAL (a-h,o-z)
-      IMPLICIT INTEGER (i-n)
-      DIMENSION xa(n),ya(n),y2a(n)
-     
-      kl0=1
-     
-      khi=n
-c means of bisection
- 1    IF(khi-kl0.gt.1) THEN
-     
-       k=(khi+kl0)/2
-     
-       IF(xa(k).gt.x) THEN
-     
-        khi=k
-     
-       ELSE
-     
-        kl0=k
-     
-       END IF
-     
-       GO TO 1
-     
-      END IF
-c KL0 and KHI now bracket the X
-      h=xa(khi)-xa(kl0)
-     
-      IF(h.eq.0.0) STOP
-      a=(xa(khi)-x)/h
-c evaluation of cubic spline polynomial
-      b=(x-xa(kl0))/h
-     
-      y=a*ya(kl0)+b*ya(khi)+((a**3-a)*y2a(kl0)+(b**3-b)*y2a(khi))*(h**2)
-     
-     ./6.
-     
-c
-     
-      RETURN
-     
-      END
-     
-
Index: /LMDZ4/trunk/libf/dyn3dpar/startvar.F90
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/startvar.F90	(revision 1424)
+++ /LMDZ4/trunk/libf/dyn3dpar/startvar.F90	(revision 1425)
@@ -643,8 +643,10 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2,lon_in2,&
-                          lat_in2, pls_in, var3d, ibar)
-!
-!-------------------------------------------------------------------------------
+SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2, &
+     lon_in2, lat_in2, pls_in, var3d, ibar)
+
+  use pchsp_95_m, only: pchsp_95
+  use pchfe_95_m, only: pchfe_95
+
 ! Arguments:
   CHARACTER(LEN=*),             INTENT(IN)    :: varname
@@ -655,73 +657,71 @@
   REAL, DIMENSION(iml),         INTENT(IN)    :: lon_in2
   REAL, DIMENSION(jml2),        INTENT(IN)    :: lat_in2
-  REAL, DIMENSION(iml,jml,lml), INTENT(IN)    :: pls_in
-  REAL, DIMENSION(iml,jml,lml), INTENT(OUT)   :: var3d
+  REAL, DIMENSION(iml, jml, lml), INTENT(IN)    :: pls_in
+  REAL, DIMENSION(iml, jml, lml), INTENT(OUT)   :: var3d
   LOGICAL,                      INTENT(IN)    :: ibar
-!-------------------------------------------------------------------------------
+!----------------------------------------------------------------------------
 ! Local variables:
 #include "iniprint.h"
-  LOGICAL               :: check=.TRUE.
-  REAL                  :: bx, by, chmin, chmax
-  INTEGER               :: ii, ij, il
-  REAL,    DIMENSION(:,:,:), ALLOCATABLE :: var_tmp3d
+  LOGICAL:: check=.TRUE., skip
+  REAL                  chmin, chmax
+  INTEGER ii, ij, il, ierr
+  integer n_extrap ! number of extrapolated points
+  REAL, DIMENSION(iml, jml, llm_dyn):: var_tmp3d
   REAL,    DIMENSION(:),     ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
-  REAL,    DIMENSION(:),     ALLOCATABLE :: lev_dyn, ax, ay, yder
-  INTEGER, DIMENSION(:),     ALLOCATABLE :: lind
-!-------------------------------------------------------------------------------
-  IF(check) WRITE(lunout,*)'Going into flinget to extract the 3D  field.'
-  IF(check) WRITE(lunout,*)fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn
-  IF(check) WRITE(lunout,*)'Allocating space for interpolation',iml,jml,llm_dyn
-
-  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn,jml_dyn,llm_dyn))
-  CALL flinget(fid_dyn,varname,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d)
+  REAL, DIMENSION(llm_dyn):: lev_dyn, ax, ay, yder
+
+!---------------------------------------------------------------------------
+  IF(check) WRITE(lunout, *)'Going into flinget to extract the 3D  field.'
+  IF(check) WRITE(lunout, *) fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, &
+       ttm_dyn
+  IF(check) WRITE(lunout, *) 'Allocating space for interpolation', iml, jml, &
+       llm_dyn
+
+  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
+  CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, 1, 1, &
+       var_ana3d)
 
 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
-  ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
-  lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
-  lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
+  ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn))
+  lon_ini(:)=lon_dyn(:, 1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
+  lat_ini(:)=lat_dyn(1, :); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
 
 !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
-  ALLOCATE(lon_rad(iml_dyn),lat_rad(jml_dyn),lev_dyn(llm_dyn))
-  CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini,       &
+  ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn))
+  CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini,      &
                    levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
-  DEALLOCATE(lon_ini,lat_ini)
+  DEALLOCATE(lon_ini, lat_ini)
 
 !--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
-  ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
-  DO il=1,llm_dyn
-    CALL interp_startvar(varname, ibar, il==1,                                 &
-      iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana3d(:,:,il), iml, jml, jml2,   &
-      lon_in,  lat_in,  lon_in2, lat_in2, var_tmp3d(:,:,il))
+  DO il=1, llm_dyn
+    CALL interp_startvar(varname, ibar, il==1, iml_dyn, jml_dyn, lon_rad, &
+         lat_rad, var_ana3d(:, :, il), iml, jml, jml2, lon_in, lat_in, &
+         lon_in2, lat_in2, var_tmp3d(:, :, il))
   END DO
-  DEALLOCATE(lon_rad,lat_rad)
-
-  ALLOCATE(lind(llm_dyn))
-  DO il=1,llm_dyn
-    lind(il) = llm_dyn-il+1
+  DEALLOCATE(lon_rad, lat_rad)
+
+!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
+  ax = lev_dyn(llm_dyn:1:-1) 
+  skip = .false.
+  n_extrap = 0
+  DO ij=1, jml
+    DO ii=1, iml-1
+      ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
+      yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.)
+      CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), &
+           var3d(ii, ij, lml:1:-1), ierr)
+      if (ierr < 0) stop 1
+      n_extrap = n_extrap + ierr
+    END DO
   END DO
-
-!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
-  ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn))
-  DO ij=1,jml
-    DO ii=1,iml-1
-      ax(:)=lev_dyn(lind(:)) 
-      ay(:)=var_tmp3d(ii,ij,lind(:))
-      CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder)
-      DO il=1,lml
-        bx=pls_in(ii,ij,il)
-        CALL SPLINT(ax, ay, yder, llm_dyn, bx, by)
-        var3d(ii,ij,il) = by
-      END DO
-    END DO
-    var3d(iml,ij,:) = var3d(1,ij,:) 
+  if (n_extrap /= 0) then
+     print *, "start_inter_3d pchfe_95: n_extrap = ", n_extrap
+  end if
+  var3d(iml, :, :) = var3d(1, :, :) 
+
+  DO il=1, lml
+    CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax)
+    WRITE(lunout, *)' '//TRIM(varname)//'  min max l ', il, chmin, chmax
   END DO
-  DEALLOCATE(var_tmp3d,lev_dyn,ax,ay,yder,lind)
-
-  DO il=1,lml
-    CALL minmax(iml*jml,var3d(1,1,il),chmin,chmax)
-    WRITE(lunout,*)' '//TRIM(varname)//'  min max l ',il,chmin,chmax
-  END DO
-
-  RETURN
 
 END SUBROUTINE start_inter_3d
Index: /LMDZ4/trunk/libf/phylmd/physiq.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/physiq.F	(revision 1424)
+++ /LMDZ4/trunk/libf/phylmd/physiq.F	(revision 1425)
@@ -1127,4 +1127,5 @@
 
       integer, save::  read_climoz ! read ozone climatology
+C     (let it keep the default OpenMP shared attribute)
 C     Allowed values are 0, 1 and 2
 C     0: do not read an ozone climatology
@@ -1134,6 +1135,8 @@
 
       integer, save:: ncid_climoz ! NetCDF file containing ozone climatologies
+C     (let it keep the default OpenMP shared attribute)
 
       real, pointer, save:: press_climoz(:)
+C     (let it keep the default OpenMP shared attribute)
 !     edges of pressure intervals for ozone climatologies, in Pa, in strictly
 !     ascending order
