Index: LMDZ6/trunk/libf/misc/cbrt.F
===================================================================
--- LMDZ6/trunk/libf/misc/cbrt.F	(revision 5214)
+++ 	(revision )
@@ -1,11 +1,0 @@
-      FUNCTION cbrt(x)
-      ! Return the cubic root for positive or negative x
-      IMPLICIT NONE
-
-      REAL x,cbrt
-
-      cbrt=sign(1.,x)*(abs(x)**(1./3.))
-
-      RETURN
-      END
-
Index: LMDZ6/trunk/libf/misc/cbrt.f90
===================================================================
--- LMDZ6/trunk/libf/misc/cbrt.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/cbrt.f90	(revision 5246)
@@ -0,0 +1,11 @@
+FUNCTION cbrt(x)
+  ! ! Return the cubic root for positive or negative x
+  IMPLICIT NONE
+
+  REAL :: x,cbrt
+
+  cbrt=sign(1.,x)*(abs(x)**(1./3.))
+
+  RETURN
+END FUNCTION cbrt
+
Index: LMDZ6/trunk/libf/misc/chfev.F
===================================================================
--- LMDZ6/trunk/libf/misc/chfev.F	(revision 5214)
+++ 	(revision )
@@ -1,155 +1,0 @@
-*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: LMDZ6/trunk/libf/misc/chfev.f90
===================================================================
--- LMDZ6/trunk/libf/misc/chfev.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/chfev.f90	(revision 5246)
@@ -0,0 +1,155 @@
+!DECK CHFEV
+SUBROUTINE CHFEV (X1, X2, F1, F2, D1, D2, NE, XE, FE, NEXT, IERR)
+  !***BEGIN PROLOGUE  CHFEV
+  !***PURPOSE  Evaluate a cubic polynomial given in Hermite form at an
+         ! array of points.  While designed for use by PCHFE, it may
+         ! be useful directly as an evaluator for a piecewise cubic
+         ! Hermite function in applications, such as graphing, where
+         ! the interval is known in advance.
+  !***LIBRARY   SLATEC (PCHIP)
+  !***CATEGORY  E3
+  !***TYPE      SINGLE PRECISION (CHFEV-S, DCHFEV-D)
+  !***KEYWORDS  CUBIC HERMITE EVALUATION, CUBIC POLYNOMIAL EVALUATION,
+         !  PCHIP
+  !***AUTHOR  Fritsch, F. N., (LLNL)
+         !  Lawrence Livermore National Laboratory
+         !  P.O. Box 808  (L-316)
+         !  Livermore, CA  94550
+         !  FTS 532-4275, (510) 422-4275
+  !***DESCRIPTION
+  !
+  !      CHFEV:  Cubic Hermite Function EValuator
+  !
+  ! Evaluates the cubic polynomial determined by function values
+  ! F1,F2 and derivatives D1,D2 on interval (X1,X2) at the points
+  ! XE(J), J=1(1)NE.
+  !
+  ! ----------------------------------------------------------------------
+  !
+  !  Calling sequence:
+  !
+  !    INTEGER  NE, NEXT(2), IERR
+  !    REAL  X1, X2, F1, F2, D1, D2, XE(NE), FE(NE)
+  !
+  !    CALL  CHFEV (X1,X2, F1,F2, D1,D2, NE, XE, FE, NEXT, IERR)
+  !
+  !   Parameters:
+  !
+  ! X1,X2 -- (input) endpoints of interval of definition of cubic.
+  !       (Error return if  X1.EQ.X2 .)
+  !
+  ! F1,F2 -- (input) values of function at X1 and X2, respectively.
+  !
+  ! D1,D2 -- (input) values of derivative at X1 and X2, respectively.
+  !
+  ! NE -- (input) number of evaluation points.  (Error return if
+  !       NE.LT.1 .)
+  !
+  ! XE -- (input) real array of points at which the function is to be
+  !       evaluated.  If any of the XE are outside the interval
+  !       [X1,X2], a warning error is returned in NEXT.
+  !
+  ! FE -- (output) real array of values of the cubic function defined
+  !       by  X1,X2, F1,F2, D1,D2  at the points  XE.
+  !
+  ! NEXT -- (output) integer array indicating number of extrapolation
+  !       points:
+  !        NEXT(1) = number of evaluation points to left of interval.
+  !        NEXT(2) = number of evaluation points to right of interval.
+  !
+  ! IERR -- (output) error flag.
+  !       Normal return:
+  !          IERR = 0  (no errors).
+  !       "Recoverable" errors:
+  !          IERR = -1  if NE.LT.1 .
+  !          IERR = -2  if X1.EQ.X2 .
+  !            (The FE-array has not been changed in either case.)
+  !
+  !***REFERENCES  (NONE)
+  !***ROUTINES CALLED  XERMSG
+  !***REVISION HISTORY  (YYMMDD)
+  !   811019  DATE WRITTEN
+  !   820803  Minor cosmetic changes for release 1.
+  !   890411  Added SAVE statements (Vers. 3.2).
+  !   890531  Changed all specific intrinsics to generic.  (WRB)
+  !   890703  Corrected category record.  (WRB)
+  !   890703  REVISION DATE from Version 3.2
+  !   891214  Prologue converted to Version 4.0 format.  (BAB)
+  !   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
+  !***END PROLOGUE  CHFEV
+  !  Programming notes:
+  !
+  ! To produce a double precision version, simply:
+  !    a. Change CHFEV to DCHFEV wherever it occurs,
+  !    b. Change the real declaration to double precision, and
+  !    c. Change the constant ZERO to double precision.
+  !
+  !  DECLARE ARGUMENTS.
+  !
+  INTEGER :: NE, NEXT(2), IERR
+  REAL :: X1, X2, F1, F2, D1, D2, XE(*), FE(*)
+  !
+  !  DECLARE LOCAL VARIABLES.
+  !
+  INTEGER :: I
+  REAL :: C2, C3, DEL1, DEL2, DELTA, H, X, XMI, XMA, ZERO
+  SAVE ZERO
+  DATA  ZERO /0./
+  !
+  !  VALIDITY-CHECK ARGUMENTS.
+  !
+  !***FIRST EXECUTABLE STATEMENT  CHFEV
+  IF (NE .LT. 1)  GO TO 5001
+  H = X2 - X1
+  IF (H .EQ. ZERO)  GO TO 5002
+  !
+  !  INITIALIZE.
+  !
+  IERR = 0
+  NEXT(1) = 0
+  NEXT(2) = 0
+  XMI = MIN(ZERO, H)
+  XMA = MAX(ZERO, H)
+  !
+  !  COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1).
+  !
+  DELTA = (F2 - F1)/H
+  DEL1 = (D1 - DELTA)/H
+  DEL2 = (D2 - DELTA)/H
+                                        ! (DELTA IS NO LONGER NEEDED.)
+  C2 = -(DEL1+DEL1 + DEL2)
+  C3 = (DEL1 + DEL2)/H
+                            ! (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.)
+  !
+  !  EVALUATION LOOP.
+  !
+  DO  I = 1, NE
+     X = XE(I) - X1
+     FE(I) = F1 + X*(D1 + X*(C2 + X*C3))
+       ! COUNT EXTRAPOLATION POINTS.
+     IF ( X.LT.XMI )  NEXT(1) = NEXT(1) + 1
+     IF ( X.GT.XMA )  NEXT(2) = NEXT(2) + 1
+     ! (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.)
+  END DO
+  !
+  !  NORMAL RETURN.
+  !
+  RETURN
+  !
+  !  ERROR RETURNS.
+  !
+ 5001   CONTINUE
+  ! NE.LT.1 RETURN.
+  IERR = -1
+  CALL XERMSG ('SLATEC', 'CHFEV', &
+        'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1)
+  RETURN
+  !
+ 5002   CONTINUE
+  ! X1.EQ.X2 RETURN.
+  IERR = -2
+  CALL XERMSG ('SLATEC', 'CHFEV', 'INTERVAL ENDPOINTS EQUAL', IERR, &
+        1)
+  RETURN
+  !------------- LAST LINE OF CHFEV FOLLOWS ------------------------------
+END SUBROUTINE CHFEV
Index: LMDZ6/trunk/libf/misc/cray.F
===================================================================
--- LMDZ6/trunk/libf/misc/cray.F	(revision 5214)
+++ 	(revision )
@@ -1,42 +1,0 @@
-!
-! $Header$
-!
-#ifdef CRAY
-      SUBROUTINE riencray
-      END
-#else
-      subroutine scopy(n,sx,incx,sy,incy)
-c
-      IMPLICIT NONE
-c
-      integer n,incx,incy,ix,iy,i
-      real sx((n-1)*incx+1),sy((n-1)*incy+1)
-c
-      iy=1
-      ix=1
-      do 10 i=1,n
-         sy(iy)=sx(ix)
-         ix=ix+incx
-         iy=iy+incy
-10    continue
-c
-      return
-      end
-
-      function ssum(n,sx,incx)
-c
-      IMPLICIT NONE
-c
-      integer n,incx,i,ix
-      real ssum,sx((n-1)*incx+1)
-c
-      ssum=0.
-      ix=1
-      do 10 i=1,n
-         ssum=ssum+sx(ix)
-         ix=ix+incx
-10    continue
-c
-      return
-      end
-#endif
Index: LMDZ6/trunk/libf/misc/cray.F90
===================================================================
--- LMDZ6/trunk/libf/misc/cray.F90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/cray.F90	(revision 5246)
@@ -0,0 +1,42 @@
+!
+! $Header$
+!
+#ifdef CRAY
+SUBROUTINE riencray
+END SUBROUTINE riencray
+#else
+subroutine scopy(n,sx,incx,sy,incy)
+  !
+  IMPLICIT NONE
+  !
+  integer :: n,incx,incy,ix,iy,i
+  real :: sx((n-1)*incx+1),sy((n-1)*incy+1)
+  !
+  iy=1
+  ix=1
+  do i=1,n
+     sy(iy)=sx(ix)
+     ix=ix+incx
+     iy=iy+incy
+  end do
+  !
+  return
+end subroutine scopy
+
+function ssum(n,sx,incx)
+  !
+  IMPLICIT NONE
+  !
+  integer :: n,incx,i,ix
+  real :: ssum,sx((n-1)*incx+1)
+  !
+  ssum=0.
+  ix=1
+  do i=1,n
+     ssum=ssum+sx(ix)
+     ix=ix+incx
+  end do
+  !
+  return
+end function ssum
+#endif
Index: LMDZ6/trunk/libf/misc/fdump.F
===================================================================
--- LMDZ6/trunk/libf/misc/fdump.F	(revision 5214)
+++ 	(revision )
@@ -1,31 +1,0 @@
-*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: LMDZ6/trunk/libf/misc/fdump.f90
===================================================================
--- LMDZ6/trunk/libf/misc/fdump.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/fdump.f90	(revision 5246)
@@ -0,0 +1,31 @@
+!DECK FDUMP
+SUBROUTINE FDUMP
+  !***BEGIN PROLOGUE  FDUMP
+  !***PURPOSE  Symbolic dump (should be locally written).
+  !***LIBRARY   SLATEC (XERROR)
+  !***CATEGORY  R3
+  !***TYPE      ALL (FDUMP-A)
+  !***KEYWORDS  ERROR, XERMSG
+  !***AUTHOR  Jones, R. E., (SNLA)
+  !***DESCRIPTION
+  !
+  !    ***Note*** Machine Dependent Routine
+  !    FDUMP is intended to be replaced by a locally written
+  !    version which produces a symbolic dump.  Failing this,
+  !    it should be replaced by a version which prints the
+  !    subprogram nesting list.  Note that this dump must be
+  !    printed on each of up to five files, as indicated by the
+  !    XGETUA routine.  See XSETUA and XGETUA for details.
+  !
+  ! Written by Ron Jones, with SLATEC Common Math Library Subcommittee
+  !
+  !***REFERENCES  (NONE)
+  !***ROUTINES CALLED  (NONE)
+  !***REVISION HISTORY  (YYMMDD)
+  !   790801  DATE WRITTEN
+  !   861211  REVISION DATE from Version 3.2
+  !   891214  Prologue converted to Version 4.0 format.  (BAB)
+  !***END PROLOGUE  FDUMP
+  !***FIRST EXECUTABLE STATEMENT  FDUMP
+  RETURN
+END SUBROUTINE FDUMP
Index: LMDZ6/trunk/libf/misc/formcoord.F
===================================================================
--- LMDZ6/trunk/libf/misc/formcoord.F	(revision 5214)
+++ 	(revision )
@@ -1,54 +1,0 @@
-!
-! $Header$
-!
-      subroutine formcoord(unit,n,x,a,rev,text)
-      implicit none
-      integer n,unit,ndec
-      logical rev
-      real x(n),a
-      character*4 text
-
-      integer i,id,i1,i2,in
-      real dx,dxmin
-
-      if(rev) then
-         id=-1
-         i1=n
-         i2=n-1
-         in=1
-         write(unit,3000) text(1:1)
-      else
-         id=1
-         i1=1
-         i2=2
-         in=n
-      endif
-
-      if (n.lt.2) then
-         ndec=1
-         write(unit,1000) text,n,x(1)*a
-      else
-         dxmin=abs(x(2)-x(1))
-         do i=2,n-1
-            dx=abs(x(i+1)-x(i))
-            if (dx.lt.dxmin) dxmin=dx
-         enddo
-
-         ndec=-log10(dxmin)+2
-         if(mod(n,6).eq.1) then
-            write(unit,1000) text,n,x(i1)*a
-            write(unit,2000) (x(i)*a,i=i2,in,id)
-         else
-            write(unit,1000) text,n
-            write(unit,2000) (x(i)*a,i=i1,in,id)
-         endif
-      endif
-
-1000  format(a4,2x,i4,' LEVELS',43x,f12.2)
-2000  format(6f12.2)
-c1000  format(a4,2x,i4,' LEVELS',43x,f12.<ndec>)
-c2000  format(6f12.<ndec>)
-3000  format('FORMAT ',a1,'REV')
-      return
-
-      end
Index: LMDZ6/trunk/libf/misc/formcoord.f90
===================================================================
--- LMDZ6/trunk/libf/misc/formcoord.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/formcoord.f90	(revision 5246)
@@ -0,0 +1,54 @@
+!
+! $Header$
+!
+subroutine formcoord(unit,n,x,a,rev,text)
+  implicit none
+  integer :: n,unit,ndec
+  logical :: rev
+  real :: x(n),a
+  character(len=4) :: text
+
+  integer :: i,id,i1,i2,in
+  real :: dx,dxmin
+
+  if(rev) then
+     id=-1
+     i1=n
+     i2=n-1
+     in=1
+     write(unit,3000) text(1:1)
+  else
+     id=1
+     i1=1
+     i2=2
+     in=n
+  endif
+
+  if (n.lt.2) then
+     ndec=1
+     write(unit,1000) text,n,x(1)*a
+  else
+     dxmin=abs(x(2)-x(1))
+     do i=2,n-1
+        dx=abs(x(i+1)-x(i))
+        if (dx.lt.dxmin) dxmin=dx
+     enddo
+
+     ndec=-log10(dxmin)+2
+     if(mod(n,6).eq.1) then
+        write(unit,1000) text,n,x(i1)*a
+        write(unit,2000) (x(i)*a,i=i2,in,id)
+     else
+        write(unit,1000) text,n
+        write(unit,2000) (x(i)*a,i=i1,in,id)
+     endif
+  endif
+
+1000   format(a4,2x,i4,' LEVELS',43x,f12.2)
+2000   format(6f12.2)
+  !1000  format(a4,2x,i4,' LEVELS',43x,f12.<ndec>)
+  !2000  format(6f12.<ndec>)
+3000   format('FORMAT ',a1,'REV')
+  return
+
+end subroutine formcoord
Index: LMDZ6/trunk/libf/misc/i1mach.F
===================================================================
--- LMDZ6/trunk/libf/misc/i1mach.F	(revision 5214)
+++ 	(revision )
@@ -1,128 +1,0 @@
-*DECK I1MACH
-      INTEGER FUNCTION I1MACH (I)
-      IMPLICIT NONE
-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)
-      INTEGER I
-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: LMDZ6/trunk/libf/misc/i1mach.f90
===================================================================
--- LMDZ6/trunk/libf/misc/i1mach.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/i1mach.f90	(revision 5246)
@@ -0,0 +1,128 @@
+!DECK I1MACH
+INTEGER FUNCTION I1MACH (I)
+  IMPLICIT NONE
+  !***BEGIN PROLOGUE  I1MACH
+  !***PURPOSE  Return integer machine dependent constants.
+  !***LIBRARY   SLATEC
+  !***CATEGORY  R1
+  !***TYPE      INTEGER (I1MACH-I)
+  !***KEYWORDS  MACHINE CONSTANTS
+  !***AUTHOR  Fox, P. A., (Bell Labs)
+        ! Hall, A. D., (Bell Labs)
+        ! Schryer, N. L., (Bell Labs)
+  !***DESCRIPTION
+  !
+  !   I1MACH can be used to obtain machine-dependent parameters for the
+  !   local machine environment.  It is a function subprogram with one
+  !   (input) argument and can be referenced as follows:
+  !
+  !    K = I1MACH(I)
+  !
+  !   where I=1,...,16.  The (output) value of K above is determined by
+  !   the (input) value of I.  The results for various values of I are
+  !   discussed below.
+  !
+  !   I/O unit numbers:
+  ! I1MACH( 1) = the standard input unit.
+  ! I1MACH( 2) = the standard output unit.
+  ! I1MACH( 3) = the standard punch unit.
+  ! I1MACH( 4) = the standard error message unit.
+  !
+  !   Words:
+  ! I1MACH( 5) = the number of bits per integer storage unit.
+  ! I1MACH( 6) = the number of characters per integer storage unit.
+  !
+  !   Integers:
+  ! assume integers are represented in the S-digit, base-A form
+  !
+  !            sign ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) )
+  !
+  !            where 0 .LE. X(I) .LT. A for I=0,...,S-1.
+  ! I1MACH( 7) = A, the base.
+  ! I1MACH( 8) = S, the number of base-A digits.
+  ! I1MACH( 9) = A**S - 1, the largest magnitude.
+  !
+  !   Floating-Point Numbers:
+  ! Assume floating-point numbers are represented in the T-digit,
+  ! base-B form
+  !            sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
+  !
+  !            where 0 .LE. X(I) .LT. B for I=1,...,T,
+  !            0 .LT. X(1), and EMIN .LE. E .LE. EMAX.
+  ! I1MACH(10) = B, the base.
+  !
+  !   Single-Precision:
+  ! I1MACH(11) = T, the number of base-B digits.
+  ! I1MACH(12) = EMIN, the smallest exponent E.
+  ! I1MACH(13) = EMAX, the largest exponent E.
+  !
+  !   Double-Precision:
+  ! I1MACH(14) = T, the number of base-B digits.
+  ! I1MACH(15) = EMIN, the smallest exponent E.
+  ! I1MACH(16) = EMAX, the largest exponent E.
+  !
+  !   To alter this function for a particular environment, the desired
+  !   set of DATA statements should be activated by removing the C from
+  !   column 1.  Also, the values of I1MACH(1) - I1MACH(4) should be
+  !   checked for consistency with the local operating system.
+  !
+  !***REFERENCES  P. A. Fox, A. D. Hall and N. L. Schryer, Framework for
+  !             a portable library, ACM Transactions on Mathematical
+  !             Software 4, 2 (June 1978), pp. 177-188.
+  !***ROUTINES CALLED  (NONE)
+  !***REVISION HISTORY  (YYMMDD)
+  !   750101  DATE WRITTEN
+  !   891012  Added VAX G-floating constants.  (WRB)
+  !   891012  REVISION DATE from Version 3.2
+  !   891214  Prologue converted to Version 4.0 format.  (BAB)
+  !   900618  Added DEC RISC constants.  (WRB)
+  !   900723  Added IBM RS 6000 constants.  (WRB)
+  !   901009  Correct I1MACH(7) for IBM Mainframes. Should be 2 not 16.
+  !       (RWC)
+  !   910710  Added HP 730 constants.  (SMR)
+  !   911114  Added Convex IEEE constants.  (WRB)
+  !   920121  Added SUN -r8 compiler option constants.  (WRB)
+  !   920229  Added Touchstone Delta i860 constants.  (WRB)
+  !   920501  Reformatted the REFERENCES section.  (WRB)
+  !   920625  Added Convex -p8 and -pd8 compiler option constants.
+  !       (BKS, WRB)
+  !   930201  Added DEC Alpha and SGI constants.  (RWC and WRB)
+  !   930618  Corrected I1MACH(5) for Convex -p8 and -pd8 compiler
+  !       options.  (DWL, RWC and WRB).
+  !   100623  Use Fortran 95 intrinsic functions (Lionel GUEZ)
+  !***END PROLOGUE  I1MACH
+  !
+  INTEGER :: IMACH(16),OUTPUT
+  SAVE IMACH
+  EQUIVALENCE (IMACH(4),OUTPUT)
+  INTEGER :: I
+  !***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
+  !
+  I1MACH = IMACH(I)
+  RETURN
+  !
+   10   CONTINUE
+  WRITE (UNIT = OUTPUT, FMT = 9000)
+ 9000   FORMAT ('1ERROR    1 IN I1MACH - I OUT OF BOUNDS')
+  !
+  ! CALL FDUMP
+  !
+  STOP
+END FUNCTION I1MACH
Index: LMDZ6/trunk/libf/misc/ismax.F
===================================================================
--- LMDZ6/trunk/libf/misc/ismax.F	(revision 5214)
+++ 	(revision )
@@ -1,24 +1,0 @@
-!
-! $Header$
-!
-      function ismax(n,sx,incx)
-c
-      IMPLICIT NONE
-c
-      INTEGER n,i,incx,ismax,ix
-      real sx((n-1)*incx+1),sxmax
-c
-      ix=1
-      ismax=1
-      sxmax=sx(1)
-      do 10 i=1,n-1
-       ix=ix+incx
-       if(sx(ix).gt.sxmax) then
-         sxmax=sx(ix)
-         ismax=i+1
-       endif
-10    continue
-c
-      return
-      end
-
Index: LMDZ6/trunk/libf/misc/ismax.f90
===================================================================
--- LMDZ6/trunk/libf/misc/ismax.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/ismax.f90	(revision 5246)
@@ -0,0 +1,24 @@
+!
+! $Header$
+!
+function ismax(n,sx,incx)
+  !
+  IMPLICIT NONE
+  !
+  INTEGER :: n,i,incx,ismax,ix
+  real :: sx((n-1)*incx+1),sxmax
+  !
+  ix=1
+  ismax=1
+  sxmax=sx(1)
+  do i=1,n-1
+   ix=ix+incx
+   if(sx(ix).gt.sxmax) then
+     sxmax=sx(ix)
+     ismax=i+1
+   endif
+  end do
+  !
+  return
+end function ismax
+
Index: LMDZ6/trunk/libf/misc/ismin.F
===================================================================
--- LMDZ6/trunk/libf/misc/ismin.F	(revision 5214)
+++ 	(revision )
@@ -1,24 +1,0 @@
-!
-! $Header$
-!
-      FUNCTION ismin(n,sx,incx)
-c
-      IMPLICIT NONE
-c
-      integer n,i,incx,ismin,ix
-      real sx((n-1)*incx+1),sxmin
-c
-      ix=1
-      ismin=1
-      sxmin=sx(1)
-      DO i=1,n-1
-         ix=ix+incx
-         if(sx(ix).lt.sxmin) then
-             sxmin=sx(ix)
-             ismin=i+1
-         endif
-      ENDDO
-c
-      return
-      end
-C
Index: LMDZ6/trunk/libf/misc/ismin.f90
===================================================================
--- LMDZ6/trunk/libf/misc/ismin.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/ismin.f90	(revision 5246)
@@ -0,0 +1,24 @@
+!
+! $Header$
+!
+FUNCTION ismin(n,sx,incx)
+  !
+  IMPLICIT NONE
+  !
+  integer :: n,i,incx,ismin,ix
+  real :: sx((n-1)*incx+1),sxmin
+  !
+  ix=1
+  ismin=1
+  sxmin=sx(1)
+  DO i=1,n-1
+     ix=ix+incx
+     if(sx(ix).lt.sxmin) then
+         sxmin=sx(ix)
+         ismin=i+1
+     endif
+  ENDDO
+  !
+  return
+end function ismin
+!
Index: LMDZ6/trunk/libf/misc/j4save.F
===================================================================
--- LMDZ6/trunk/libf/misc/j4save.F	(revision 5214)
+++ 	(revision )
@@ -1,67 +1,0 @@
-*DECK J4SAVE
-      FUNCTION J4SAVE (IWHICH, IVALUE, ISET)
-      IMPLICIT NONE
-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/
-      INTEGER J4SAVE,IWHICH,IVALUE
-C***FIRST EXECUTABLE STATEMENT  J4SAVE
-      J4SAVE = IPARAM(IWHICH)
-      IF (ISET) IPARAM(IWHICH) = IVALUE
-      RETURN
-      END
Index: LMDZ6/trunk/libf/misc/j4save.f90
===================================================================
--- LMDZ6/trunk/libf/misc/j4save.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/j4save.f90	(revision 5246)
@@ -0,0 +1,67 @@
+!DECK J4SAVE
+FUNCTION J4SAVE (IWHICH, IVALUE, ISET)
+  IMPLICIT NONE
+  !***BEGIN PROLOGUE  J4SAVE
+  !***SUBSIDIARY
+  !***PURPOSE  Save or recall global variables needed by error
+         ! handling routines.
+  !***LIBRARY   SLATEC (XERROR)
+  !***TYPE      INTEGER (J4SAVE-I)
+  !***KEYWORDS  ERROR MESSAGES, ERROR NUMBER, RECALL, SAVE, XERROR
+  !***AUTHOR  Jones, R. E., (SNLA)
+  !***DESCRIPTION
+  !
+  ! Abstract
+  !    J4SAVE saves and recalls several global variables needed
+  !    by the library error handling routines.
+  !
+  ! Description of Parameters
+  !  --Input--
+  !    IWHICH - Index of item desired.
+  !            = 1 Refers to current error number.
+  !            = 2 Refers to current error control flag.
+  !            = 3 Refers to current unit number to which error
+  !                messages are to be sent.  (0 means use standard.)
+  !            = 4 Refers to the maximum number of times any
+  !                 message is to be printed (as set by XERMAX).
+  !            = 5 Refers to the total number of units to which
+  !                 each error message is to be written.
+  !            = 6 Refers to the 2nd unit for error messages
+  !            = 7 Refers to the 3rd unit for error messages
+  !            = 8 Refers to the 4th unit for error messages
+  !            = 9 Refers to the 5th unit for error messages
+  !    IVALUE - The value to be set for the IWHICH-th parameter,
+  !             if ISET is .TRUE. .
+  !    ISET   - If ISET=.TRUE., the IWHICH-th parameter will BE
+  !             given the value, IVALUE.  If ISET=.FALSE., the
+  !             IWHICH-th parameter will be unchanged, and IVALUE
+  !             is a dummy parameter.
+  !  --Output--
+  !    The (old) value of the IWHICH-th parameter will be returned
+  !    in the function value, J4SAVE.
+  !
+  !***SEE ALSO  XERMSG
+  !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+  !             Error-handling Package, SAND82-0800, Sandia
+  !             Laboratories, 1982.
+  !***ROUTINES CALLED  (NONE)
+  !***REVISION HISTORY  (YYMMDD)
+  !   790801  DATE WRITTEN
+  !   891214  Prologue converted to Version 4.0 format.  (BAB)
+  !   900205  Minor modifications to prologue.  (WRB)
+  !   900402  Added TYPE section.  (WRB)
+  !   910411  Added KEYWORDS section.  (WRB)
+  !   920501  Reformatted the REFERENCES section.  (WRB)
+  !***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/
+  INTEGER :: J4SAVE,IWHICH,IVALUE
+  !***FIRST EXECUTABLE STATEMENT  J4SAVE
+  J4SAVE = IPARAM(IWHICH)
+  IF (ISET) IPARAM(IWHICH) = IVALUE
+  RETURN
+END FUNCTION J4SAVE
Index: LMDZ6/trunk/libf/misc/juldate.F
===================================================================
--- LMDZ6/trunk/libf/misc/juldate.F	(revision 5214)
+++ 	(revision )
@@ -1,39 +1,0 @@
-!
-! $Id$
-!
-	subroutine juldate(ian,imoi,ijou,oh,om,os,tjd,tjdsec)
-c	Sous-routine de changement de date:
-c	gregorien>>>date julienne
-c	En entree:an,mois,jour,heure,min.,sec.
-c	En sortie:tjd
-        IMPLICIT NONE
-        INTEGER,INTENT(IN) :: ian,imoi,ijou,oh,om,os
-        REAL,INTENT(OUT) :: tjd,tjdsec
-        
-        REAL frac,year,rmon,cf,a,b
-        INTEGER ojou
-        
-	frac=((os/60.+om)/60.+oh)/24.
-	ojou=dble(ijou)+frac
-	    year=dble(ian)
-	    rmon=dble(imoi)
-	if (imoi .le. 2) then
-	    year=year-1.
-	    rmon=rmon+12.
-	endif
-	cf=year+(rmon/100.)+(ojou/10000.)
-	if (cf .ge. 1582.1015) then
-	    a=int(year/100)
-	    b=2-a+int(a/4)
-	else
-	    b=0
-	endif
-	tjd=int(365.25*year)+int(30.6001*(rmon+1))+int(ojou)
-     +   +1720994.5+b
-        tjdsec=(ojou-int(ojou))+(tjd-int(tjd))
-        tjd=int(tjd)+int(tjdsec)
-	tjdsec=tjdsec-int(tjdsec)
-	return
-	end
-
-
Index: LMDZ6/trunk/libf/misc/juldate.f90
===================================================================
--- LMDZ6/trunk/libf/misc/juldate.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/juldate.f90	(revision 5246)
@@ -0,0 +1,39 @@
+!
+! $Id$
+!
+  subroutine juldate(ian,imoi,ijou,oh,om,os,tjd,tjdsec)
+  ! Sous-routine de changement de date:
+  ! gregorien>>>date julienne
+  ! En entree:an,mois,jour,heure,min.,sec.
+  ! En sortie:tjd
+    IMPLICIT NONE
+    INTEGER,INTENT(IN) :: ian,imoi,ijou,oh,om,os
+    REAL,INTENT(OUT) :: tjd,tjdsec
+
+    REAL :: frac,year,rmon,cf,a,b
+    INTEGER :: ojou
+
+    frac=((os/60.+om)/60.+oh)/24.
+    ojou=dble(ijou)+frac
+        year=dble(ian)
+        rmon=dble(imoi)
+    if (imoi .le. 2) then
+        year=year-1.
+        rmon=rmon+12.
+    endif
+    cf=year+(rmon/100.)+(ojou/10000.)
+    if (cf .ge. 1582.1015) then
+        a=int(year/100)
+        b=2-a+int(a/4)
+    else
+        b=0
+    endif
+    tjd=int(365.25*year)+int(30.6001*(rmon+1))+int(ojou) &
+          +1720994.5+b
+    tjdsec=(ojou-int(ojou))+(tjd-int(tjd))
+    tjd=int(tjd)+int(tjdsec)
+    tjdsec=tjdsec-int(tjdsec)
+    return
+end subroutine juldate
+
+
Index: LMDZ6/trunk/libf/misc/minmax.F
===================================================================
--- LMDZ6/trunk/libf/misc/minmax.F	(revision 5214)
+++ 	(revision )
@@ -1,23 +1,0 @@
-!
-! $Header$
-!
-       SUBROUTINE minmax(imax, xi, zmin, zmax )
-c
-c      P. Le Van
-
-       INTEGER imax
-       REAL    xi(imax)
-       REAL    zmin,zmax
-       INTEGER i
-
-       zmin = xi(1)
-       zmax = xi(1)
-
-       DO i = 2, imax
-         zmin = MIN( zmin,xi(i) )
-         zmax = MAX( zmax,xi(i) )
-       ENDDO
-
-       RETURN
-       END
-
Index: LMDZ6/trunk/libf/misc/minmax.f90
===================================================================
--- LMDZ6/trunk/libf/misc/minmax.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/minmax.f90	(revision 5246)
@@ -0,0 +1,23 @@
+!
+! $Header$
+!
+ SUBROUTINE minmax(imax, xi, zmin, zmax )
+  !
+  !  P. Le Van
+
+   INTEGER :: imax
+   REAL :: xi(imax)
+   REAL :: zmin,zmax
+   INTEGER :: i
+
+   zmin = xi(1)
+   zmax = xi(1)
+
+   DO i = 2, imax
+     zmin = MIN( zmin,xi(i) )
+     zmax = MAX( zmax,xi(i) )
+   ENDDO
+
+   RETURN
+END SUBROUTINE minmax
+
Index: LMDZ6/trunk/libf/misc/minmax2.F
===================================================================
--- LMDZ6/trunk/libf/misc/minmax2.F	(revision 5214)
+++ 	(revision )
@@ -1,20 +1,0 @@
-!
-! $Header$
-!
-       SUBROUTINE minmax2(imax, jmax, lmax, xi, zmin, zmax )
-c
-       INTEGER lmax,jmax,imax
-       REAL xi(imax*jmax*lmax) 
-       REAL zmin,zmax
-       INTEGER i
-    
-       zmin = xi(1)
-       zmax = xi(1)
-
-       DO i = 2, imax*jmax*lmax
-         zmin = MIN( zmin,xi(i) )
-         zmax = MAX( zmax,xi(i) )
-       ENDDO
-
-       RETURN
-       END
Index: LMDZ6/trunk/libf/misc/minmax2.f90
===================================================================
--- LMDZ6/trunk/libf/misc/minmax2.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/minmax2.f90	(revision 5246)
@@ -0,0 +1,20 @@
+!
+! $Header$
+!
+ SUBROUTINE minmax2(imax, jmax, lmax, xi, zmin, zmax )
+  !
+   INTEGER :: lmax,jmax,imax
+   REAL :: xi(imax*jmax*lmax)
+   REAL :: zmin,zmax
+   INTEGER :: i
+
+   zmin = xi(1)
+   zmax = xi(1)
+
+   DO i = 2, imax*jmax*lmax
+     zmin = MIN( zmin,xi(i) )
+     zmax = MAX( zmax,xi(i) )
+   ENDDO
+
+   RETURN
+END SUBROUTINE minmax2
Index: LMDZ6/trunk/libf/misc/pchdf.F
===================================================================
--- LMDZ6/trunk/libf/misc/pchdf.F	(revision 5214)
+++ 	(revision )
@@ -1,106 +1,0 @@
-*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: LMDZ6/trunk/libf/misc/pchdf.f90
===================================================================
--- LMDZ6/trunk/libf/misc/pchdf.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/pchdf.f90	(revision 5246)
@@ -0,0 +1,106 @@
+!DECK PCHDF
+REAL FUNCTION PCHDF (K, X, S, IERR)
+  !***BEGIN PROLOGUE  PCHDF
+  !***SUBSIDIARY
+  !***PURPOSE  Computes divided differences for PCHCE and PCHSP
+  !***LIBRARY   SLATEC (PCHIP)
+  !***TYPE      SINGLE PRECISION (PCHDF-S, DPCHDF-D)
+  !***AUTHOR  Fritsch, F. N., (LLNL)
+  !***DESCRIPTION
+  !
+  !      PCHDF:   PCHIP Finite Difference Formula
+  !
+  ! Uses a divided difference formulation to compute a K-point approx-
+  ! imation to the derivative at X(K) based on the data in X and S.
+  !
+  ! Called by  PCHCE  and  PCHSP  to compute 3- and 4-point boundary
+  ! derivative approximations.
+  !
+  ! ----------------------------------------------------------------------
+  !
+  ! On input:
+  !    K      is the order of the desired derivative approximation.
+  !           K must be at least 3 (error return if not).
+  !    X      contains the K values of the independent variable.
+  !           X need not be ordered, but the values **MUST** be
+  !           distinct.  (Not checked here.)
+  !    S      contains the associated slope values:
+  !              S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1.
+  !           (Note that S need only be of length K-1.)
+  !
+  ! On return:
+  !    S      will be destroyed.
+  !    IERR   will be set to -1 if K.LT.2 .
+  !    PCHDF  will be set to the desired derivative approximation if
+  !           IERR=0 or to zero if IERR=-1.
+  !
+  ! ----------------------------------------------------------------------
+  !
+  !***SEE ALSO  PCHCE, PCHSP
+  !***REFERENCES  Carl de Boor, A Practical Guide to Splines, Springer-
+  !             Verlag, New York, 1978, pp. 10-16.
+  !***ROUTINES CALLED  XERMSG
+  !***REVISION HISTORY  (YYMMDD)
+  !   820503  DATE WRITTEN
+  !   820805  Converted to SLATEC library version.
+  !   870813  Minor cosmetic changes.
+  !   890411  Added SAVE statements (Vers. 3.2).
+  !   890411  REVISION DATE from Version 3.2
+  !   891214  Prologue converted to Version 4.0 format.  (BAB)
+  !   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
+  !   900328  Added TYPE section.  (WRB)
+  !   910408  Updated AUTHOR and DATE WRITTEN sections in prologue.  (WRB)
+  !   920429  Revised format and order of references.  (WRB,FNF)
+  !   930503  Improved purpose.  (FNF)
+  !***END PROLOGUE  PCHDF
+  !
+  !**End
+  !
+  !  DECLARE ARGUMENTS.
+  !
+  INTEGER :: K, IERR
+  REAL :: X(K), S(K)
+  !
+  !  DECLARE LOCAL VARIABLES.
+  !
+  INTEGER :: I, J
+  REAL :: VALUE, ZERO
+  SAVE ZERO
+  DATA  ZERO /0./
+  !
+  !  CHECK FOR LEGAL VALUE OF K.
+  !
+  !***FIRST EXECUTABLE STATEMENT  PCHDF
+  IF (K .LT. 3)  GO TO 5001
+  !
+  !  COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL.
+  !
+  DO  J = 2, K-1
+     DO  I = 1, K-J
+        S(I) = (S(I+1)-S(I))/(X(I+J)-X(I))
+     END DO
+  END DO
+  !
+  !  EVALUATE DERIVATIVE AT X(K).
+  !
+  VALUE = S(1)
+  DO  I = 2, K-1
+     VALUE = S(I) + VALUE*(X(K)-X(I))
+  END DO
+  !
+  !  NORMAL RETURN.
+  !
+  IERR = 0
+  PCHDF = VALUE
+  RETURN
+  !
+  !  ERROR RETURN.
+  !
+ 5001   CONTINUE
+  ! K.LT.3 RETURN.
+  IERR = -1
+  CALL XERMSG ('SLATEC', 'PCHDF', 'K LESS THAN THREE', IERR, 1)
+  PCHDF = ZERO
+  RETURN
+  !------------- LAST LINE OF PCHDF FOLLOWS ------------------------------
+END FUNCTION PCHDF
Index: LMDZ6/trunk/libf/misc/pchfe.F
===================================================================
--- LMDZ6/trunk/libf/misc/pchfe.F	(revision 5214)
+++ 	(revision )
@@ -1,308 +1,0 @@
-*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: LMDZ6/trunk/libf/misc/pchfe.f90
===================================================================
--- LMDZ6/trunk/libf/misc/pchfe.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/pchfe.f90	(revision 5246)
@@ -0,0 +1,308 @@
+!DECK PCHFE
+SUBROUTINE PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
+  !***BEGIN PROLOGUE  PCHFE
+  !***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.
+  !***LIBRARY   SLATEC (PCHIP)
+  !***CATEGORY  E3
+  !***TYPE      SINGLE PRECISION (PCHFE-S, DPCHFE-D)
+  !***KEYWORDS  CUBIC HERMITE EVALUATION, HERMITE INTERPOLATION, PCHIP,
+         !  PIECEWISE CUBIC EVALUATION
+  !***AUTHOR  Fritsch, F. N., (LLNL)
+         !  Lawrence Livermore National Laboratory
+         !  P.O. Box 808  (L-316)
+         !  Livermore, CA  94550
+         !  FTS 532-4275, (510) 422-4275
+  !***DESCRIPTION
+  !
+  !      PCHFE:  Piecewise Cubic Hermite Function Evaluator
+  !
+  ! Evaluates the cubic Hermite function defined by  N, X, F, D  at
+  ! the points  XE(J), J=1(1)NE.
+  !
+  ! To provide compatibility with PCHIM and PCHIC, includes an
+  ! increment between successive values of the F- and D-arrays.
+  !
+  ! ----------------------------------------------------------------------
+  !
+  !  Calling sequence:
+  !
+  !    PARAMETER  (INCFD = ...)
+  !    INTEGER  N, NE, IERR
+  !    REAL  X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE)
+  !    LOGICAL  SKIP
+  !
+  !    CALL  PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
+  !
+  !   Parameters:
+  !
+  ! N -- (input) number of data points.  (Error return if N.LT.2 .)
+  !
+  ! X -- (input) real array of independent variable values.  The
+  !       elements of X must be strictly increasing:
+  !            X(I-1) .LT. X(I),  I = 2(1)N.
+  !       (Error return if not.)
+  !
+  ! F -- (input) real array of function values.  F(1+(I-1)*INCFD) is
+  !       the value corresponding to X(I).
+  !
+  ! D -- (input) real array of derivative values.  D(1+(I-1)*INCFD) is
+  !       the value corresponding to X(I).
+  !
+  ! INCFD -- (input) increment between successive values in F and D.
+  !       (Error return if  INCFD.LT.1 .)
+  !
+  ! SKIP -- (input/output) logical variable which should be set to
+  !       .TRUE. if the user wishes to skip checks for validity of
+  !       preceding parameters, or to .FALSE. otherwise.
+  !       This 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.
+  !
+  ! NE -- (input) number of evaluation points.  (Error return if
+  !       NE.LT.1 .)
+  !
+  ! XE -- (input) real array of 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.
+  !
+  ! FE -- (output) real array of values of the cubic Hermite function
+  !       defined by  N, X, F, D  at the points  XE.
+  !
+  ! IERR -- (output) error flag.
+  !       Normal return:
+  !          IERR = 0  (no errors).
+  !       Warning error:
+  !          IERR.GT.0  means that extrapolation was performed at
+  !             IERR points.
+  !       "Recoverable" errors:
+  !          IERR = -1  if N.LT.2 .
+  !          IERR = -2  if INCFD.LT.1 .
+  !          IERR = -3  if the X-array is not strictly increasing.
+  !          IERR = -4  if NE.LT.1 .
+  !         (The FE-array has not been changed in any of these cases.)
+  !           NOTE:  The above errors are checked in the order listed,
+  !               and following arguments have **NOT** been validated.
+  !
+  !***REFERENCES  (NONE)
+  !***ROUTINES CALLED  CHFEV, XERMSG
+  !***REVISION HISTORY  (YYMMDD)
+  !   811020  DATE WRITTEN
+  !   820803  Minor cosmetic changes for release 1.
+  !   870707  Minor cosmetic changes to prologue.
+  !   890531  Changed all specific intrinsics to generic.  (WRB)
+  !   890831  Modified array declarations.  (WRB)
+  !   890831  REVISION DATE from Version 3.2
+  !   891214  Prologue converted to Version 4.0 format.  (BAB)
+  !   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
+  !***END PROLOGUE  PCHFE
+  !  Programming notes:
+  !
+  ! 1. To produce a double precision version, simply:
+  !    a. Change PCHFE to DPCHFE, and CHFEV to DCHFEV, wherever they
+  !       occur,
+  !    b. Change the real declaration to double precision,
+  !
+  ! 2. Most of the coding between the call to CHFEV and the end of
+  !    the IR-loop could be eliminated if it were permissible to
+  !    assume that XE is ordered relative to X.
+  !
+  ! 3. CHFEV does not assume that X1 is less than X2.  thus, it would
+  !    be possible to write a version of PCHFE that assumes a strict-
+  !    ly decreasing X-array by simply running the IR-loop backwards
+  !    (and reversing the order of appropriate tests).
+  !
+  ! 4. The present code has a minor bug, which I have decided is not
+  !    worth the effort that would be required to fix it.
+  !    If XE contains points in [X(N-1),X(N)], followed by points .LT.
+  !    X(N-1), followed by points .GT.X(N), the extrapolation points
+  !    will be counted (at least) twice in the total returned in IERR.
+  !
+  !  DECLARE ARGUMENTS.
+  !
+  INTEGER :: N, INCFD, NE, IERR
+  REAL :: X(*), F(INCFD,*), D(INCFD,*), XE(*), FE(*)
+  LOGICAL :: SKIP
+  !
+  !  DECLARE LOCAL VARIABLES.
+  !
+  INTEGER :: I, IERC, IR, J, JFIRST, NEXT(2), NJ
+  !
+  !  VALIDITY-CHECK ARGUMENTS.
+  !
+  !***FIRST EXECUTABLE STATEMENT  PCHFE
+  IF (SKIP)  GO TO 5
+  !
+  IF ( N.LT.2 )  GO TO 5001
+  IF ( INCFD.LT.1 )  GO TO 5002
+  DO  I = 2, N
+     IF ( X(I).LE.X(I-1) )  GO TO 5003
+  END DO
+  !
+  !  FUNCTION DEFINITION IS OK, GO ON.
+  !
+    5   CONTINUE
+  IF ( NE.LT.1 )  GO TO 5004
+  IERR = 0
+  SKIP = .TRUE.
+  !
+  !  LOOP OVER INTERVALS.        (   INTERVAL INDEX IS  IL = IR-1  . )
+  !                              ( INTERVAL IS X(IL).LE.X.LT.X(IR) . )
+  JFIRST = 1
+  IR = 2
+   10   CONTINUE
+  !
+  ! SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS.
+  !
+     IF (JFIRST .GT. NE)  GO TO 5000
+  !
+  ! LOCATE ALL POINTS IN INTERVAL.
+  !
+     DO  J = JFIRST, NE
+        IF (XE(J) .GE. X(IR))  GO TO 30
+     END DO
+     J = NE + 1
+     GO TO 40
+  !
+  ! HAVE LOCATED FIRST POINT BEYOND INTERVAL.
+  !
+   30   CONTINUE
+     IF (IR .EQ. N)  J = NE + 1
+  !
+   40   CONTINUE
+     NJ = J - JFIRST
+  !
+  ! SKIP EVALUATION IF NO POINTS IN INTERVAL.
+  !
+     IF (NJ .EQ. 0)  GO TO 50
+  !
+  ! EVALUATE CUBIC AT XE(I),  I = JFIRST (1) J-1 .
+  !
+  !   ----------------------------------------------------------------
+    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)
+    ! ----------------------------------------------------------------
+     IF (IERC .LT. 0)  GO TO 5005
+  !
+     IF (NEXT(2) .EQ. 0)  GO TO 42
+     ! IF (NEXT(2) .GT. 0)  THEN
+     !    IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE
+     !    RIGHT OF X(IR).
+  !
+        IF (IR .LT. N)  GO TO 41
+        ! IF (IR .EQ. N)  THEN
+        !    THESE ARE ACTUALLY EXTRAPOLATION POINTS.
+           IERR = IERR + NEXT(2)
+           GO TO 42
+   41    CONTINUE
+        ! ELSE
+        !    WE SHOULD NEVER HAVE GOTTEN HERE.
+           GO TO 5005
+        ! ENDIF
+     ! ENDIF
+   42   CONTINUE
+  !
+     IF (NEXT(1) .EQ. 0)  GO TO 49
+     ! IF (NEXT(1) .GT. 0)  THEN
+     !    IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(1) TO THE
+     !    LEFT OF X(IR-1).
+  !
+        IF (IR .GT. 2)  GO TO 43
+        ! IF (IR .EQ. 2)  THEN
+        !    THESE ARE ACTUALLY EXTRAPOLATION POINTS.
+           IERR = IERR + NEXT(1)
+           GO TO 49
+   43    CONTINUE
+        ! ELSE
+        !    XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST
+        !    EVALUATION INTERVAL.
+  !
+  !          FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1).
+           DO  I = JFIRST, J-1
+              IF (XE(I) .LT. X(IR-1))  GO TO 45
+           END DO
+           ! NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR
+           !        IN CHFEV.
+           GO TO 5005
+  !
+   45       CONTINUE
+           ! RESET J.  (THIS WILL BE THE NEW JFIRST.)
+           J = I
+  !
+  !          NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY.
+           DO  I = 1, IR-1
+              IF (XE(J) .LT. X(I)) GO TO 47
+           END DO
+           ! NB-- CAN NEVER DROP THROUGH HERE, SINCE XE(J).LT.X(IR-1).
+  !
+   47       CONTINUE
+           ! AT THIS POINT, EITHER  XE(J) .LT. X(1)
+           !    OR      X(I-1) .LE. XE(J) .LT. X(I) .
+           ! RESET IR, RECOGNIZING THAT IT WILL BE INCREMENTED BEFORE
+           ! CYCLING.
+           IR = MAX(1, I-1)
+        ! ENDIF
+     ! ENDIF
+   49   CONTINUE
+  !
+     JFIRST = J
+  !
+  ! END OF IR-LOOP.
+  !
+   50   CONTINUE
+  IR = IR + 1
+  IF (IR .LE. N)  GO TO 10
+  !
+  !  NORMAL RETURN.
+  !
+ 5000   CONTINUE
+  RETURN
+  !
+  !  ERROR RETURNS.
+  !
+ 5001   CONTINUE
+  ! N.LT.2 RETURN.
+  IERR = -1
+  CALL XERMSG ('SLATEC', 'PCHFE', &
+        'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
+  RETURN
+  !
+ 5002   CONTINUE
+  ! INCFD.LT.1 RETURN.
+  IERR = -2
+  CALL XERMSG ('SLATEC', 'PCHFE', 'INCREMENT LESS THAN ONE', IERR, &
+        1)
+  RETURN
+  !
+ 5003   CONTINUE
+  ! X-ARRAY NOT STRICTLY INCREASING.
+  IERR = -3
+  CALL XERMSG ('SLATEC', 'PCHFE', 'X-ARRAY NOT STRICTLY INCREASING' &
+        , IERR, 1)
+  RETURN
+  !
+ 5004   CONTINUE
+  ! NE.LT.1 RETURN.
+  IERR = -4
+  CALL XERMSG ('SLATEC', 'PCHFE', &
+        'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1)
+  RETURN
+  !
+ 5005   CONTINUE
+  ! ERROR RETURN FROM CHFEV.
+  !   *** THIS CASE SHOULD NEVER OCCUR ***
+  IERR = -5
+  CALL XERMSG ('SLATEC', 'PCHFE', &
+        'ERROR RETURN FROM CHFEV -- FATAL', IERR, 2)
+  RETURN
+  !------------- LAST LINE OF PCHFE FOLLOWS ------------------------------
+END SUBROUTINE PCHFE
Index: LMDZ6/trunk/libf/misc/pchsp.F
===================================================================
--- LMDZ6/trunk/libf/misc/pchsp.F	(revision 5214)
+++ 	(revision )
@@ -1,388 +1,0 @@
-*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: LMDZ6/trunk/libf/misc/pchsp.f90
===================================================================
--- LMDZ6/trunk/libf/misc/pchsp.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/pchsp.f90	(revision 5246)
@@ -0,0 +1,388 @@
+!DECK PCHSP
+SUBROUTINE PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
+  !***BEGIN PROLOGUE  PCHSP
+  !***PURPOSE  Set derivatives needed to determine the Hermite represen-
+         ! tation of the cubic spline interpolant to given data, with
+         ! specified boundary conditions.
+  !***LIBRARY   SLATEC (PCHIP)
+  !***CATEGORY  E1A
+  !***TYPE      SINGLE PRECISION (PCHSP-S, DPCHSP-D)
+  !***KEYWORDS  CUBIC HERMITE INTERPOLATION, PCHIP,
+         !  PIECEWISE CUBIC INTERPOLATION, SPLINE INTERPOLATION
+  !***AUTHOR  Fritsch, F. N., (LLNL)
+         !  Lawrence Livermore National Laboratory
+         !  P.O. Box 808  (L-316)
+         !  Livermore, CA  94550
+         !  FTS 532-4275, (510) 422-4275
+  !***DESCRIPTION
+  !
+  !      PCHSP:   Piecewise Cubic Hermite Spline
+  !
+  ! Computes the Hermite representation of the cubic spline inter-
+  ! polant to the data given in X and F satisfying the boundary
+  ! conditions specified by IC and VC.
+  !
+  ! To facilitate two-dimensional applications, includes an increment
+  ! between successive values of the F- and D-arrays.
+  !
+  ! 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.
+  !
+  ! ----------------------------------------------------------------------
+  !
+  !  Calling sequence:
+  !
+  !    PARAMETER  (INCFD = ...)
+  !    INTEGER  IC(2), N, NWK, IERR
+  !    REAL  VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(NWK)
+  !
+  !    CALL  PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
+  !
+  !   Parameters:
+  !
+  ! IC -- (input) integer array of length 2 specifying desired
+  !       boundary conditions:
+  !       IC(1) = IBEG, desired condition at beginning of data.
+  !       IC(2) = IEND, desired condition at end of data.
+  !
+  !       IBEG = 0  to set D(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(1).
+  !       IBEG = 2  if second derivative at X(1) is given in VC(1).
+  !       IBEG = 3  to use the 3-point difference formula for D(1).
+  !                 (Reverts to the default b.c. if N.LT.3 .)
+  !       IBEG = 4  to use the 4-point difference formula for D(1).
+  !                 (Reverts to the default b.c. if N.LT.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(1)=0.
+  !
+  !       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(2).
+  !
+  !      NOTES:
+  !       1. An error return is taken if IEND is out of range.
+  !       2. For the "natural" boundary condition, use IEND=2 and
+  !          VC(2)=0.
+  !
+  ! VC -- (input) real array of length 2 specifying desired boundary
+  !       values, as indicated above.
+  !       VC(1) need be set only if IC(1) = 1 or 2 .
+  !       VC(2) need be set only if IC(2) = 1 or 2 .
+  !
+  ! N -- (input) number of data points.  (Error return if N.LT.2 .)
+  !
+  ! X -- (input) real array of independent variable values.  The
+  !       elements of X must be strictly increasing:
+  !            X(I-1) .LT. X(I),  I = 2(1)N.
+  !       (Error return if not.)
+  !
+  ! F -- (input) real array of dependent variable values to be inter-
+  !       polated.  F(1+(I-1)*INCFD) is value corresponding to X(I).
+  !
+  ! D -- (output) real array of 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
+  !            D(1+(I-1)*INCFD),  I=1(1)N.
+  !       No other entries in D are changed.
+  !
+  ! INCFD -- (input) increment between successive values in F and D.
+  !       This argument is provided primarily for 2-D applications.
+  !       (Error return if  INCFD.LT.1 .)
+  !
+  ! WK -- (scratch) real array of working storage.
+  !
+  ! NWK -- (input) length of work array.
+  !       (Error return if NWK.LT.2*N .)
+  !
+  ! IERR -- (output) error flag.
+  !       Normal return:
+  !          IERR = 0  (no errors).
+  !       "Recoverable" errors:
+  !          IERR = -1  if N.LT.2 .
+  !          IERR = -2  if INCFD.LT.1 .
+  !          IERR = -3  if the X-array is not strictly increasing.
+  !          IERR = -4  if IBEG.LT.0 or IBEG.GT.4 .
+  !          IERR = -5  if IEND.LT.0 of IEND.GT.4 .
+  !          IERR = -6  if both of the above are true.
+  !          IERR = -7  if NWK is too small.
+  !           NOTE:  The above errors are checked in the order listed,
+  !               and following arguments have **NOT** been validated.
+  !         (The D-array has not been changed in any of these cases.)
+  !          IERR = -8  in case of trouble solving the linear system
+  !                     for the interior derivative values.
+  !         (The D-array may have been changed in this case.)
+  !         (             Do **NOT** use it!                )
+  !
+  !***REFERENCES  Carl de Boor, A Practical Guide to Splines, Springer-
+  !             Verlag, New York, 1978, pp. 53-59.
+  !***ROUTINES CALLED  PCHDF, XERMSG
+  !***REVISION HISTORY  (YYMMDD)
+  !   820503  DATE WRITTEN
+  !   820804  Converted to SLATEC library version.
+  !   870707  Minor cosmetic changes to prologue.
+  !   890411  Added SAVE statements (Vers. 3.2).
+  !   890703  Corrected category record.  (WRB)
+  !   890831  Modified array declarations.  (WRB)
+  !   890831  REVISION DATE from Version 3.2
+  !   891214  Prologue converted to Version 4.0 format.  (BAB)
+  !   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
+  !   920429  Revised format and order of references.  (WRB,FNF)
+  !***END PROLOGUE  PCHSP
+  !  Programming notes:
+  !
+  ! To produce a double precision version, simply:
+  !    a. Change PCHSP to DPCHSP wherever it occurs,
+  !    b. Change the real declarations to double precision, and
+  !    c. Change the constants ZERO, HALF, ... to double precision.
+  !
+  !  DECLARE ARGUMENTS.
+  !
+  INTEGER :: IC(2), N, INCFD, NWK, IERR
+  REAL :: VC(2), X(*), F(INCFD,*), D(INCFD,*), WK(2,*)
+  !
+  !  DECLARE LOCAL VARIABLES.
+  !
+  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
+  !
+  DATA  ZERO /0./,  HALF /0.5/,  ONE /1./,  TWO /2./,  THREE /3./
+  !
+  !  VALIDITY-CHECK ARGUMENTS.
+  !
+  !***FIRST EXECUTABLE STATEMENT  PCHSP
+  IF ( N.LT.2 )  GO TO 5001
+  IF ( INCFD.LT.1 )  GO TO 5002
+  DO  J = 2, N
+     IF ( X(J).LE.X(J-1) )  GO TO 5003
+  END DO
+  !
+  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
+  !
+  !  FUNCTION DEFINITION IS OK -- GO ON.
+  !
+  IF ( NWK .LT. 2*N )  GO TO 5007
+  !
+  !  COMPUTE FIRST DIFFERENCES OF X SEQUENCE AND STORE IN WK(1,.). ALSO,
+  !  COMPUTE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN WK(2,.).
+  DO  J=2,N
+     WK(1,J) = X(J) - X(J-1)
+     WK(2,J) = (F(1,J) - F(1,J-1))/WK(1,J)
+  END DO
+  !
+  !  SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL.
+  !
+  IF ( IBEG.GT.N )  IBEG = 0
+  IF ( IEND.GT.N )  IEND = 0
+  !
+  !  SET UP FOR BOUNDARY CONDITIONS.
+  !
+  IF ( (IBEG.EQ.1).OR.(IBEG.EQ.2) )  THEN
+     D(1,1) = VC(1)
+  ELSE IF (IBEG .GT. 2)  THEN
+     ! PICK UP FIRST IBEG POINTS, IN REVERSE ORDER.
+     DO  J = 1, IBEG
+        INDEX = IBEG-J+1
+        ! INDEX RUNS FROM IBEG DOWN TO 1.
+        XTEMP(J) = X(INDEX)
+        IF (J .LT. IBEG)  STEMP(J) = WK(2,INDEX)
+     END DO
+              ! --------------------------------
+     D(1,1) = PCHDF (IBEG, XTEMP, STEMP, IERR)
+              ! --------------------------------
+     IF (IERR .NE. 0)  GO TO 5009
+     IBEG = 1
+  ENDIF
+  !
+  IF ( (IEND.EQ.1).OR.(IEND.EQ.2) )  THEN
+     D(1,N) = VC(2)
+  ELSE IF (IEND .GT. 2)  THEN
+     ! PICK UP LAST IEND POINTS.
+     DO  J = 1, IEND
+        INDEX = N-IEND+J
+        ! INDEX RUNS FROM N+1-IEND UP TO N.
+        XTEMP(J) = X(INDEX)
+        IF (J .LT. IEND)  STEMP(J) = WK(2,INDEX+1)
+     END DO
+              ! --------------------------------
+     D(1,N) = PCHDF (IEND, XTEMP, STEMP, IERR)
+              ! --------------------------------
+     IF (IERR .NE. 0)  GO TO 5009
+     IEND = 1
+  ENDIF
+  !
+  ! --------------------( BEGIN CODING FROM CUBSPL )--------------------
+  !
+  !  **** A TRIDIAGONAL LINEAR SYSTEM FOR THE UNKNOWN SLOPES S(J) OF
+  !  F  AT X(J), J=1,...,N, IS GENERATED AND THEN SOLVED BY GAUSS ELIM-
+  !  INATION, WITH S(J) ENDING UP IN D(1,J), ALL J.
+  ! WK(1,.) AND WK(2,.) ARE USED FOR TEMPORARY STORAGE.
+  !
+  !  CONSTRUCT FIRST EQUATION FROM FIRST BOUNDARY CONDITION, OF THE FORM
+  !         WK(2,1)*S(1) + WK(1,1)*S(2) = D(1,1)
+  !
+  IF (IBEG .EQ. 0)  THEN
+     IF (N .EQ. 2)  THEN
+        ! NO CONDITION AT LEFT END AND N = 2.
+        WK(2,1) = ONE
+        WK(1,1) = ONE
+        D(1,1) = TWO*WK(2,2)
+     ELSE
+        ! 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
+     ! SLOPE PRESCRIBED AT LEFT END.
+     WK(2,1) = ONE
+     WK(1,1) = ZERO
+  ELSE
+     ! 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
+  !
+  !  IF THERE ARE INTERIOR KNOTS, GENERATE THE CORRESPONDING EQUATIONS AND
+  !  CARRY OUT THE FORWARD PASS OF GAUSS ELIMINATION, AFTER WHICH THE J-TH
+  !  EQUATION READS    WK(2,J)*S(J) + WK(1,J)*S(J+1) = D(1,J).
+  !
+  NM1 = N-1
+  IF (NM1 .GT. 1)  THEN
+     DO 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))
+     END DO
+  ENDIF
+  !
+  !  CONSTRUCT LAST EQUATION FROM SECOND BOUNDARY CONDITION, OF THE FORM
+  !       (-G*WK(2,N-1))*S(N-1) + WK(2,N)*S(N) = D(1,N)
+  !
+  ! IF SLOPE IS PRESCRIBED AT RIGHT END, ONE CAN GO DIRECTLY TO BACK-
+  ! SUBSTITUTION, SINCE ARRAYS HAPPEN TO BE SET UP JUST RIGHT FOR IT
+  ! AT THIS POINT.
+  IF (IEND .EQ. 1)  GO TO 30
+  !
+  IF (IEND .EQ. 0)  THEN
+     IF (N.EQ.2 .AND. IBEG.EQ.0)  THEN
+        ! 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
+        ! EITHER (N=3 AND NOT-A-KNOT ALSO AT LEFT) OR (N=2 AND *NOT*
+        ! 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
+        ! NOT-A-KNOT AND N .GE. 3, AND EITHER N.GT.3 OR  ALSO NOT-A-
+        ! KNOT AT LEFT END POINT.
+        G = WK(1,N-1) + WK(1,N)
+        ! 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
+     ! 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
+  !
+  !  COMPLETE FORWARD PASS OF GAUSS ELIMINATION.
+  !
+  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)
+  !
+  !  CARRY OUT BACK SUBSTITUTION
+  !
+   30   CONTINUE
+  DO 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)
+  END DO
+  ! --------------------(  END  CODING FROM CUBSPL )--------------------
+  !
+  !  NORMAL RETURN.
+  !
+  RETURN
+  !
+  !  ERROR RETURNS.
+  !
+ 5001   CONTINUE
+  ! N.LT.2 RETURN.
+  IERR = -1
+  CALL XERMSG ('SLATEC', 'PCHSP', &
+        'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
+  RETURN
+  !
+ 5002   CONTINUE
+  ! INCFD.LT.1 RETURN.
+  IERR = -2
+  CALL XERMSG ('SLATEC', 'PCHSP', 'INCREMENT LESS THAN ONE', IERR, &
+        1)
+  RETURN
+  !
+ 5003   CONTINUE
+  ! X-ARRAY NOT STRICTLY INCREASING.
+  IERR = -3
+  CALL XERMSG ('SLATEC', 'PCHSP', 'X-ARRAY NOT STRICTLY INCREASING' &
+        , IERR, 1)
+  RETURN
+  !
+ 5004   CONTINUE
+  ! IC OUT OF RANGE RETURN.
+  IERR = IERR - 3
+  CALL XERMSG ('SLATEC', 'PCHSP', 'IC OUT OF RANGE', IERR, 1)
+  RETURN
+  !
+ 5007   CONTINUE
+  ! NWK TOO SMALL RETURN.
+  IERR = -7
+  CALL XERMSG ('SLATEC', 'PCHSP', 'WORK ARRAY TOO SMALL', IERR, 1)
+  RETURN
+  !
+ 5008   CONTINUE
+  ! SINGULAR SYSTEM.
+  !   *** THEORETICALLY, THIS CAN ONLY OCCUR IF SUCCESSIVE X-VALUES   ***
+  !   *** ARE EQUAL, WHICH SHOULD ALREADY HAVE BEEN CAUGHT (IERR=-3). ***
+  IERR = -8
+  CALL XERMSG ('SLATEC', 'PCHSP', 'SINGULAR LINEAR SYSTEM', IERR, &
+        1)
+  RETURN
+  !
+ 5009   CONTINUE
+  ! ERROR RETURN FROM PCHDF.
+  !   *** THIS CASE SHOULD NEVER OCCUR ***
+  IERR = -9
+  CALL XERMSG ('SLATEC', 'PCHSP', 'ERROR RETURN FROM PCHDF', IERR, &
+        1)
+  RETURN
+  !------------- LAST LINE OF PCHSP FOLLOWS ------------------------------
+END SUBROUTINE PCHSP
Index: LMDZ6/trunk/libf/misc/q_sat.F
===================================================================
--- LMDZ6/trunk/libf/misc/q_sat.F	(revision 5214)
+++ 	(revision )
@@ -1,72 +1,0 @@
-!
-! $Header$
-!
-!
-!
-
-      subroutine q_sat(np,temp,pres,qsat)
-!
-      IMPLICIT none
-!======================================================================
-! Autheur(s): Z.X. Li (LMD/CNRS)
-!  reecriture vectorisee par F. Hourdin.
-! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
-!======================================================================
-! Arguments:
-! kelvin---input-R: temperature en Kelvin
-! millibar--input-R: pression en mb
-!
-! q_sat----output-R: vapeur d'eau saturante en kg/kg
-!======================================================================
-!
-      integer np
-      REAL temp(np),pres(np),qsat(np)
-!
-      REAL r2es
-      PARAMETER (r2es=611.14 *18.0153/28.9644)
-!
-      REAL r3les, r3ies, r3es
-      PARAMETER (R3LES=17.269)
-      PARAMETER (R3IES=21.875)
-!
-      REAL r4les, r4ies, r4es
-      PARAMETER (R4LES=35.86)
-      PARAMETER (R4IES=7.66)
-!
-      REAL rtt
-      PARAMETER (rtt=273.16)
-!
-      REAL retv
-      PARAMETER (retv=28.9644/18.0153 - 1.0)
-
-      real zqsat
-      integer ip
-!
-!     ------------------------------------------------------------------
-!
-!
-
-      do ip=1,np
-
-!      write(*,*)'kelvin,millibar=',kelvin,millibar
-!       write(*,*)'temp,pres=',temp(ip),pres(ip)
-!
-         IF (temp(ip) .LE. rtt) THEN
-            r3es = r3ies
-            r4es = r4ies
-         ELSE
-            r3es = r3les
-            r4es = r4les
-         ENDIF
-!
-         zqsat=r2es/pres(ip)*EXP(r3es*(temp(ip)-rtt)/(temp(ip)-r4es))
-         zqsat=MIN(0.5,ZQSAT)
-         zqsat=zqsat/(1.-retv *zqsat)
-!
-         qsat(ip)= zqsat
-!      write(*,*)'qsat=',qsat(ip)
-
-      enddo
-!
-      RETURN
-      END
Index: LMDZ6/trunk/libf/misc/q_sat.f90
===================================================================
--- LMDZ6/trunk/libf/misc/q_sat.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/q_sat.f90	(revision 5246)
@@ -0,0 +1,72 @@
+!
+! $Header$
+!
+!
+!
+
+subroutine q_sat(np,temp,pres,qsat)
+  !
+  IMPLICIT none
+  !======================================================================
+  ! Autheur(s): Z.X. Li (LMD/CNRS)
+  !  reecriture vectorisee par F. Hourdin.
+  ! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
+  !======================================================================
+  ! Arguments:
+  ! kelvin---input-R: temperature en Kelvin
+  ! millibar--input-R: pression en mb
+  !
+  ! q_sat----output-R: vapeur d'eau saturante en kg/kg
+  !======================================================================
+  !
+  integer :: np
+  REAL :: temp(np),pres(np),qsat(np)
+  !
+  REAL :: r2es
+  PARAMETER (r2es=611.14 *18.0153/28.9644)
+  !
+  REAL :: r3les, r3ies, r3es
+  PARAMETER (R3LES=17.269)
+  PARAMETER (R3IES=21.875)
+  !
+  REAL :: r4les, r4ies, r4es
+  PARAMETER (R4LES=35.86)
+  PARAMETER (R4IES=7.66)
+  !
+  REAL :: rtt
+  PARAMETER (rtt=273.16)
+  !
+  REAL :: retv
+  PARAMETER (retv=28.9644/18.0153 - 1.0)
+
+  real :: zqsat
+  integer :: ip
+  !
+  ! ------------------------------------------------------------------
+  !
+  !
+
+  do ip=1,np
+
+   ! write(*,*)'kelvin,millibar=',kelvin,millibar
+   !  write(*,*)'temp,pres=',temp(ip),pres(ip)
+  !
+     IF (temp(ip) .LE. rtt) THEN
+        r3es = r3ies
+        r4es = r4ies
+     ELSE
+        r3es = r3les
+        r4es = r4les
+     ENDIF
+  !
+     zqsat=r2es/pres(ip)*EXP(r3es*(temp(ip)-rtt)/(temp(ip)-r4es))
+     zqsat=MIN(0.5,ZQSAT)
+     zqsat=zqsat/(1.-retv *zqsat)
+  !
+     qsat(ip)= zqsat
+   ! write(*,*)'qsat=',qsat(ip)
+
+  enddo
+  !
+  RETURN
+END SUBROUTINE q_sat
Index: LMDZ6/trunk/libf/misc/ran1.F
===================================================================
--- LMDZ6/trunk/libf/misc/ran1.F	(revision 5214)
+++ 	(revision )
@@ -1,38 +1,0 @@
-!
-! $Id$
-!
-      FUNCTION RAN1(IDUM)
-      IMPLICIT NONE
-      REAL RAN1
-      REAL,SAVE :: R(97)
-      REAL,PARAMETER :: RM1=3.8580247E-6,RM2=7.4373773E-6
-      INTEGER,SAVE :: IFF=0
-      integer,save :: ix1,ix2,ix3
-      INTEGER,PARAMETER :: M1=259200,IA1=7141,IC1=54773
-      INTEGER,PARAMETER :: M2=134456,IA2=8121,IC2=28411
-      INTEGER,PARAMETER :: M3=243000,IA3=4561,IC3=51349
-      INTEGER :: IDUM,J
-
-      IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
-        IFF=1
-        IX1=MOD(IC1-IDUM,M1)
-        IX1=MOD(IA1*IX1+IC1,M1)
-        IX2=MOD(IX1,M2)
-        IX1=MOD(IA1*IX1+IC1,M1)
-        IX3=MOD(IX1,M3)
-        DO 11 J=1,97
-          IX1=MOD(IA1*IX1+IC1,M1)
-          IX2=MOD(IA2*IX2+IC2,M2)
-          R(J)=(REAL(IX1)+REAL(IX2)*RM2)*RM1
-11      CONTINUE
-        IDUM=1
-      ENDIF
-      IX1=MOD(IA1*IX1+IC1,M1)
-      IX2=MOD(IA2*IX2+IC2,M2)
-      IX3=MOD(IA3*IX3+IC3,M3)
-      J=1+(97*IX3)/M3
-      IF(J.GT.97.OR.J.LT.1) stop 1
-      RAN1=R(J)
-      R(J)=(REAL(IX1)+REAL(IX2)*RM2)*RM1
-      RETURN
-      END
Index: LMDZ6/trunk/libf/misc/ran1.f90
===================================================================
--- LMDZ6/trunk/libf/misc/ran1.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/ran1.f90	(revision 5246)
@@ -0,0 +1,38 @@
+!
+! $Id$
+!
+FUNCTION RAN1(IDUM)
+  IMPLICIT NONE
+  REAL :: RAN1
+  REAL,SAVE :: R(97)
+  REAL,PARAMETER :: RM1=3.8580247E-6,RM2=7.4373773E-6
+  INTEGER,SAVE :: IFF=0
+  integer,save :: ix1,ix2,ix3
+  INTEGER,PARAMETER :: M1=259200,IA1=7141,IC1=54773
+  INTEGER,PARAMETER :: M2=134456,IA2=8121,IC2=28411
+  INTEGER,PARAMETER :: M3=243000,IA3=4561,IC3=51349
+  INTEGER :: IDUM,J
+
+  IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
+    IFF=1
+    IX1=MOD(IC1-IDUM,M1)
+    IX1=MOD(IA1*IX1+IC1,M1)
+    IX2=MOD(IX1,M2)
+    IX1=MOD(IA1*IX1+IC1,M1)
+    IX3=MOD(IX1,M3)
+    DO J=1,97
+      IX1=MOD(IA1*IX1+IC1,M1)
+      IX2=MOD(IA2*IX2+IC2,M2)
+      R(J)=(REAL(IX1)+REAL(IX2)*RM2)*RM1
+    END DO
+    IDUM=1
+  ENDIF
+  IX1=MOD(IA1*IX1+IC1,M1)
+  IX2=MOD(IA2*IX2+IC2,M2)
+  IX3=MOD(IA3*IX3+IC3,M3)
+  J=1+(97*IX3)/M3
+  IF(J.GT.97.OR.J.LT.1) stop 1
+  RAN1=R(J)
+  R(J)=(REAL(IX1)+REAL(IX2)*RM2)*RM1
+  RETURN
+END FUNCTION RAN1
Index: LMDZ6/trunk/libf/misc/sort.F
===================================================================
--- LMDZ6/trunk/libf/misc/sort.F	(revision 5214)
+++ 	(revision )
@@ -1,37 +1,0 @@
-!
-! $Header$
-!
-C
-C
-      SUBROUTINE sort(n,d)
-c
-c     P.Le Van
-c      
-c...  cette routine met le tableau d  dans l'ordre croissant  ....
-cc   ( pour avoir l'ordre decroissant,il suffit de remplacer l'instruc
-c      tion  situee + bas  IF(d(j).LE.p)  THEN     par
-c                           IF(d(j).GE.p)  THEN
-c
-
-      INTEGER n
-      REAL d(n) , p
-      INTEGER i,j,k
-
-      DO i=1,n-1
-        k=i
-        p=d(i)
-        DO j=i+1,n
-         IF(d(j).LE.p) THEN
-           k=j
-           p=d(j)
-         ENDIF
-        ENDDO
-
-       IF(k.ne.i) THEN
-         d(k)=d(i)
-         d(i)=p
-       ENDIF
-      ENDDO
-
-       RETURN
-       END
Index: LMDZ6/trunk/libf/misc/sort.f90
===================================================================
--- LMDZ6/trunk/libf/misc/sort.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/sort.f90	(revision 5246)
@@ -0,0 +1,37 @@
+!
+! $Header$
+!
+!
+!
+SUBROUTINE sort(n,d)
+  !
+  ! P.Le Van
+  !
+  !...  cette routine met le tableau d  dans l'ordre croissant  ....
+  !c   ( pour avoir l'ordre decroissant,il suffit de remplacer l'instruc
+  !  tion  situee + bas  IF(d(j).LE.p)  THEN     par
+  !                       IF(d(j).GE.p)  THEN
+  !
+
+  INTEGER :: n
+  REAL :: d(n) , p
+  INTEGER :: i,j,k
+
+  DO i=1,n-1
+    k=i
+    p=d(i)
+    DO j=i+1,n
+     IF(d(j).LE.p) THEN
+       k=j
+       p=d(j)
+     ENDIF
+    ENDDO
+
+   IF(k.ne.i) THEN
+     d(k)=d(i)
+     d(i)=p
+   ENDIF
+  ENDDO
+
+   RETURN
+END SUBROUTINE sort
Index: LMDZ6/trunk/libf/misc/xercnt.F
===================================================================
--- LMDZ6/trunk/libf/misc/xercnt.F	(revision 5214)
+++ 	(revision )
@@ -1,62 +1,0 @@
-*DECK XERCNT
-      SUBROUTINE XERCNT (LIBRAR, SUBROU, MESSG, NERR, LEVEL, KONTRL)
-      IMPLICIT NONE
-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
-      INTEGER NERR, LEVEL, KONTRL
-C***FIRST EXECUTABLE STATEMENT  XERCNT
-      RETURN
-      END
Index: LMDZ6/trunk/libf/misc/xercnt.f90
===================================================================
--- LMDZ6/trunk/libf/misc/xercnt.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/xercnt.f90	(revision 5246)
@@ -0,0 +1,62 @@
+!DECK XERCNT
+SUBROUTINE XERCNT (LIBRAR, SUBROU, MESSG, NERR, LEVEL, KONTRL)
+  IMPLICIT NONE
+  !***BEGIN PROLOGUE  XERCNT
+  !***SUBSIDIARY
+  !***PURPOSE  Allow user control over handling of errors.
+  !***LIBRARY   SLATEC (XERROR)
+  !***CATEGORY  R3C
+  !***TYPE      ALL (XERCNT-A)
+  !***KEYWORDS  ERROR, XERROR
+  !***AUTHOR  Jones, R. E., (SNLA)
+  !***DESCRIPTION
+  !
+  ! Abstract
+  !    Allows user control over handling of individual errors.
+  !    Just after each message is recorded, but before it is
+  !    processed any further (i.e., before it is printed or
+  !    a decision to abort is made), a call is made to XERCNT.
+  !    If the user has provided his own version of XERCNT, he
+  !    can then override the value of KONTROL used in processing
+  !    this message by redefining its value.
+  !    KONTRL may be set to any value from -2 to 2.
+  !    The meanings for KONTRL are the same as in XSETF, except
+  !    that the value of KONTRL changes only for this message.
+  !    If KONTRL is set to a value outside the range from -2 to 2,
+  !    it will be moved back into that range.
+  !
+  ! Description of Parameters
+  !
+  !  --Input--
+  !    LIBRAR - the library that the routine is in.
+  !    SUBROU - the subroutine that XERMSG is being called from
+  !    MESSG  - the first 20 characters of the error message.
+  !    NERR   - same as in the call to XERMSG.
+  !    LEVEL  - same as in the call to XERMSG.
+  !    KONTRL - the current value of the control flag as set
+  !             by a call to XSETF.
+  !
+  !  --Output--
+  !    KONTRL - the new value of KONTRL.  If KONTRL is not
+  !             defined, it will remain at its original value.
+  !             This changed value of control affects only
+  !             the current occurrence of the current message.
+  !
+  !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+  !             Error-handling Package, SAND82-0800, Sandia
+  !             Laboratories, 1982.
+  !***ROUTINES CALLED  (NONE)
+  !***REVISION HISTORY  (YYMMDD)
+  !   790801  DATE WRITTEN
+  !   861211  REVISION DATE from Version 3.2
+  !   891214  Prologue converted to Version 4.0 format.  (BAB)
+  !   900206  Routine changed from user-callable to subsidiary.  (WRB)
+  !   900510  Changed calling sequence to include LIBRARY and SUBROUTINE
+  !       names, changed routine name from XERCTL to XERCNT.  (RWC)
+  !   920501  Reformatted the REFERENCES section.  (WRB)
+  !***END PROLOGUE  XERCNT
+  CHARACTER(len=*) :: LIBRAR, SUBROU, MESSG
+  INTEGER :: NERR, LEVEL, KONTRL
+  !***FIRST EXECUTABLE STATEMENT  XERCNT
+  RETURN
+END SUBROUTINE XERCNT
Index: LMDZ6/trunk/libf/misc/xerhlt.F
===================================================================
--- LMDZ6/trunk/libf/misc/xerhlt.F	(revision 5214)
+++ 	(revision )
@@ -1,39 +1,0 @@
-*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: LMDZ6/trunk/libf/misc/xerhlt.f90
===================================================================
--- LMDZ6/trunk/libf/misc/xerhlt.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/xerhlt.f90	(revision 5246)
@@ -0,0 +1,39 @@
+!DECK XERHLT
+SUBROUTINE XERHLT (MESSG)
+  !***BEGIN PROLOGUE  XERHLT
+  !***SUBSIDIARY
+  !***PURPOSE  Abort program execution and print error message.
+  !***LIBRARY   SLATEC (XERROR)
+  !***CATEGORY  R3C
+  !***TYPE      ALL (XERHLT-A)
+  !***KEYWORDS  ABORT PROGRAM EXECUTION, ERROR, XERROR
+  !***AUTHOR  Jones, R. E., (SNLA)
+  !***DESCRIPTION
+  !
+  ! Abstract
+  !    ***Note*** machine dependent routine
+  !    XERHLT aborts the execution of the program.
+  !    The error message causing the abort is given in the calling
+  !    sequence, in case one needs it for printing on a dayfile,
+  !    for example.
+  !
+  ! Description of Parameters
+  !    MESSG is as in XERMSG.
+  !
+  !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+  !             Error-handling Package, SAND82-0800, Sandia
+  !             Laboratories, 1982.
+  !***ROUTINES CALLED  (NONE)
+  !***REVISION HISTORY  (YYMMDD)
+  !   790801  DATE WRITTEN
+  !   861211  REVISION DATE from Version 3.2
+  !   891214  Prologue converted to Version 4.0 format.  (BAB)
+  !   900206  Routine changed from user-callable to subsidiary.  (WRB)
+  !   900510  Changed calling sequence to delete length of character
+  !       and changed routine name from XERABT to XERHLT.  (RWC)
+  !   920501  Reformatted the REFERENCES section.  (WRB)
+  !***END PROLOGUE  XERHLT
+  CHARACTER(len=*) :: MESSG
+  !***FIRST EXECUTABLE STATEMENT  XERHLT
+  STOP
+END SUBROUTINE XERHLT
Index: LMDZ6/trunk/libf/misc/xermsg.F
===================================================================
--- LMDZ6/trunk/libf/misc/xermsg.F	(revision 5214)
+++ 	(revision )
@@ -1,368 +1,0 @@
-*DECK XERMSG
-      SUBROUTINE XERMSG (LIBRAR, SUBROU, MESSG, NERR, LEVEL)
-      IMPLICIT NONE
-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
-      INTEGER NERR, LEVEL, LKNTRL
-      INTEGER J4SAVE, MAXMES, KDUMMY, I, KOUNT, LERR, LLEVEL
-      INTEGER MKNTRL, LTEMP
-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: LMDZ6/trunk/libf/misc/xermsg.f90
===================================================================
--- LMDZ6/trunk/libf/misc/xermsg.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/xermsg.f90	(revision 5246)
@@ -0,0 +1,368 @@
+!DECK XERMSG
+SUBROUTINE XERMSG (LIBRAR, SUBROU, MESSG, NERR, LEVEL)
+  IMPLICIT NONE
+  !***BEGIN PROLOGUE  XERMSG
+  !***PURPOSE  Process error messages for SLATEC and other libraries.
+  !***LIBRARY   SLATEC (XERROR)
+  !***CATEGORY  R3C
+  !***TYPE      ALL (XERMSG-A)
+  !***KEYWORDS  ERROR MESSAGE, XERROR
+  !***AUTHOR  Fong, Kirby, (NMFECC at LLNL)
+  !***DESCRIPTION
+  !
+  !   XERMSG processes a diagnostic message in a manner determined by the
+  !   value of LEVEL and the current value of the library error control
+  !   flag, KONTRL.  See subroutine XSETF for details.
+  !
+  !    LIBRAR   A character constant (or character variable) with the name
+  !         of the library.  This will be 'SLATEC' for the SLATEC
+  !         Common Math Library.  The error handling package is
+  !         general enough to be used by many libraries
+  !         simultaneously, so it is desirable for the routine that
+  !         detects and reports an error to identify the library name
+  !         as well as the routine name.
+  !
+  !    SUBROU   A character constant (or character variable) with the name
+  !             of the routine that detected the error.  Usually it is the
+  !         name of the routine that is calling XERMSG.  There are
+  !         some instances where a user callable library routine calls
+  !         lower level subsidiary routines where the error is
+  !         detected.  In such cases it may be more informative to
+  !         supply the name of the routine the user called rather than
+  !         the name of the subsidiary routine that detected the
+  !         error.
+  !
+  !    MESSG    A character constant (or character variable) with the text
+  !         of the error or warning message.  In the example below,
+  !         the message is a character constant that contains a
+  !         generic message.
+  !
+  !               CALL XERMSG ('SLATEC', 'MMPY',
+  !              *'THE ORDER OF THE MATRIX EXCEEDS THE ROW DIMENSION',
+  !              *3, 1)
+  !
+  !         It is possible (and is sometimes desirable) to generate a
+  !         specific message--e.g., one that contains actual numeric
+  !         values.  Specific numeric values can be converted into
+  !         character strings using formatted WRITE statements into
+  !         character variables.  This is called standard Fortran
+  !         internal file I/O and is exemplified in the first three
+  !         lines of the following example.  You can also catenate
+  !         substrings of characters to construct the error message.
+  !         Here is an example showing the use of both writing to
+  !         an internal file and catenating character strings.
+  !
+  !               CHARACTER*5 CHARN, CHARL
+  !               WRITE (CHARN,10) N
+  !               WRITE (CHARL,10) LDA
+  !            10 FORMAT(I5)
+  !               CALL XERMSG ('SLATEC', 'MMPY', 'THE ORDER'//CHARN//
+  !              *   ' OF THE MATRIX EXCEEDS ITS ROW DIMENSION OF'//
+  !              *   CHARL, 3, 1)
+  !
+  !         There are two subtleties worth mentioning.  One is that
+  !         the // for character catenation is used to construct the
+  !         error message so that no single character constant is
+  !         continued to the next line.  This avoids confusion as to
+  !         whether there are trailing blanks at the end of the line.
+  !         The second is that by catenating the parts of the message
+  !             as an actual argument rather than encoding the entire
+  !         message into one large character variable, we avoid
+  !         having to know how long the message will be in order to
+  !         declare an adequate length for that large character
+  !         variable.  XERMSG calls XERPRN to print the message using
+  !         multiple lines if necessary.  If the message is very long,
+  !         XERPRN will break it into pieces of 72 characters (as
+  !         requested by XERMSG) for printing on multiple lines.
+  !         Also, XERMSG asks XERPRN to prefix each line with ' *  '
+  !         so that the total line length could be 76 characters.
+  !         Note also that XERPRN scans the error message backwards
+  !         to ignore trailing blanks.  Another feature is that
+  !         the substring '$$' is treated as a new line sentinel
+  !         by XERPRN.  If you want to construct a multiline
+  !         message without having to count out multiples of 72
+  !         characters, just use '$$' as a separator.  '$$'
+  !         obviously must occur within 72 characters of the
+  !         start of each line to have its intended effect since
+  !         XERPRN is asked to wrap around at 72 characters in
+  !         addition to looking for '$$'.
+  !
+  !    NERR     An integer value that is chosen by the library routine's
+  !         author.  It must be in the range -99 to 999 (three
+  !         printable digits).  Each distinct error should have its
+  !         own error number.  These error numbers should be described
+  !         in the machine readable documentation for the routine.
+  !         The error numbers need be unique only within each routine,
+  !         so it is reasonable for each routine to start enumerating
+  !         errors from 1 and proceeding to the next integer.
+  !
+  !    LEVEL    An integer value in the range 0 to 2 that indicates the
+  !         level (severity) of the error.  Their meanings are
+  !
+  !        -1  A warning message.  This is used if it is not clear
+  !            that there really is an error, but the user's attention
+  !            may be needed.  An attempt is made to only print this
+  !            message once.
+  !
+  !         0  A warning message.  This is used if it is not clear
+  !            that there really is an error, but the user's attention
+  !            may be needed.
+  !
+  !         1  A recoverable error.  This is used even if the error is
+  !            so serious that the routine cannot return any useful
+  !            answer.  If the user has told the error package to
+  !            return after recoverable errors, then XERMSG will
+  !            return to the Library routine which can then return to
+  !            the user's routine.  The user may also permit the error
+  !            package to terminate the program upon encountering a
+  !            recoverable error.
+  !
+  !         2  A fatal error.  XERMSG will not return to its caller
+  !            after it receives a fatal error.  This level should
+  !            hardly ever be used; it is much better to allow the
+  !            user a chance to recover.  An example of one of the few
+  !            cases in which it is permissible to declare a level 2
+  !            error is a reverse communication Library routine that
+  !            is likely to be called repeatedly until it integrates
+  !            across some interval.  If there is a serious error in
+  !            the input such that another step cannot be taken and
+  !            the Library routine is called again without the input
+  !            error having been corrected by the caller, the Library
+  !            routine will probably be called forever with improper
+  !            input.  In this case, it is reasonable to declare the
+  !            error to be fatal.
+  !
+  !    Each of the arguments to XERMSG is input; none will be modified by
+  !    XERMSG.  A routine may make multiple calls to XERMSG with warning
+  !    level messages; however, after a call to XERMSG with a recoverable
+  !    error, the routine should return to the user.  Do not try to call
+  !    XERMSG with a second recoverable error after the first recoverable
+  !    error because the error package saves the error number.  The user
+  !    can retrieve this error number by calling another entry point in
+  !    the error handling package and then clear the error number when
+  !    recovering from the error.  Calling XERMSG in succession causes the
+  !    old error number to be overwritten by the latest error number.
+  !    This is considered harmless for error numbers associated with
+  !    warning messages but must not be done for error numbers of serious
+  !    errors.  After a call to XERMSG with a recoverable error, the user
+  !    must be given a chance to call NUMXER or XERCLR to retrieve or
+  !    clear the error number.
+  !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+  !             Error-handling Package, SAND82-0800, Sandia
+  !             Laboratories, 1982.
+  !***ROUTINES CALLED  FDUMP, J4SAVE, XERCNT, XERHLT, XERPRN, XERSVE
+  !***REVISION HISTORY  (YYMMDD)
+  !   880101  DATE WRITTEN
+  !   880621  REVISED AS DIRECTED AT SLATEC CML MEETING OF FEBRUARY 1988.
+  !       THERE ARE TWO BASIC CHANGES.
+  !       1.  A NEW ROUTINE, XERPRN, IS USED INSTEAD OF XERPRT TO
+  !           PRINT MESSAGES.  THIS ROUTINE WILL BREAK LONG MESSAGES
+  !           INTO PIECES FOR PRINTING ON MULTIPLE LINES.  '$$' IS
+  !           ACCEPTED AS A NEW LINE SENTINEL.  A PREFIX CAN BE
+  !           ADDED TO EACH LINE TO BE PRINTED.  XERMSG USES EITHER
+  !           ' ***' OR ' *  ' AND LONG MESSAGES ARE BROKEN EVERY
+  !           72 CHARACTERS (AT MOST) SO THAT THE MAXIMUM LINE
+  !           LENGTH OUTPUT CAN NOW BE AS GREAT AS 76.
+  !       2.  THE TEXT OF ALL MESSAGES IS NOW IN UPPER CASE SINCE THE
+  !           FORTRAN STANDARD DOCUMENT DOES NOT ADMIT THE EXISTENCE
+  !           OF LOWER CASE.
+  !   880708  REVISED AFTER THE SLATEC CML MEETING OF JUNE 29 AND 30.
+  !       THE PRINCIPAL CHANGES ARE
+  !       1.  CLARIFY COMMENTS IN THE PROLOGUES
+  !       2.  RENAME XRPRNT TO XERPRN
+  !       3.  REWORK HANDLING OF '$$' IN XERPRN TO HANDLE BLANK LINES
+  !           SIMILAR TO THE WAY FORMAT STATEMENTS HANDLE THE /
+  !           CHARACTER FOR NEW RECORDS.
+  !   890706  REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO
+  !       CLEAN UP THE CODING.
+  !   890721  REVISED TO USE NEW FEATURE IN XERPRN TO COUNT CHARACTERS IN
+  !       PREFIX.
+  !   891013  REVISED TO CORRECT COMMENTS.
+  !   891214  Prologue converted to Version 4.0 format.  (WRB)
+  !   900510  Changed test on NERR to be -9999999 < NERR < 99999999, but
+  !       NERR .ne. 0, and on LEVEL to be -2 < LEVEL < 3.  Added
+  !       LEVEL=-1 logic, changed calls to XERSAV to XERSVE, and
+  !       XERCTL to XERCNT.  (RWC)
+  !   920501  Reformatted the REFERENCES section.  (WRB)
+  !***END PROLOGUE  XERMSG
+  CHARACTER(len=*) :: LIBRAR, SUBROU, MESSG
+  CHARACTER(len=8) :: XLIBR, XSUBR
+  CHARACTER(len=72) :: TEMP
+  CHARACTER(len=20) :: LFIRST
+  INTEGER :: NERR, LEVEL, LKNTRL
+  INTEGER :: J4SAVE, MAXMES, KDUMMY, I, KOUNT, LERR, LLEVEL
+  INTEGER :: MKNTRL, LTEMP
+  !***FIRST EXECUTABLE STATEMENT  XERMSG
+  LKNTRL = J4SAVE (2, 0, .FALSE.)
+  MAXMES = J4SAVE (4, 0, .FALSE.)
+  !
+  !   LKNTRL IS A LOCAL COPY OF THE CONTROL FLAG KONTRL.
+  !   MAXMES IS THE MAXIMUM NUMBER OF TIMES ANY PARTICULAR MESSAGE
+  !      SHOULD BE PRINTED.
+  !
+  !   WE PRINT A FATAL ERROR MESSAGE AND TERMINATE FOR AN ERROR IN
+  !      CALLING XERMSG.  THE ERROR NUMBER SHOULD BE POSITIVE,
+  !      AND THE LEVEL SHOULD BE BETWEEN 0 AND 2.
+  !
+  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
+  !
+  !   RECORD THE MESSAGE.
+  !
+  I = J4SAVE (1, NERR, .TRUE.)
+  CALL XERSVE (LIBRAR, SUBROU, MESSG, 1, NERR, LEVEL, KOUNT)
+  !
+  !   HANDLE PRINT-ONCE WARNING MESSAGES.
+  !
+  IF (LEVEL.EQ.-1 .AND. KOUNT.GT.1) RETURN
+  !
+  !   ALLOW TEMPORARY USER OVERRIDE OF THE CONTROL FLAG.
+  !
+  XLIBR  = LIBRAR
+  XSUBR  = SUBROU
+  LFIRST = MESSG
+  LERR   = NERR
+  LLEVEL = LEVEL
+  CALL XERCNT (XLIBR, XSUBR, LFIRST, LERR, LLEVEL, LKNTRL)
+  !
+  LKNTRL = MAX(-2, MIN(2,LKNTRL))
+  MKNTRL = ABS(LKNTRL)
+  !
+  !   SKIP PRINTING IF THE CONTROL FLAG VALUE AS RESET IN XERCNT IS
+  !   ZERO AND THE ERROR IS NOT FATAL.
+  !
+  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
+  !
+  !   ANNOUNCE THE NAMES OF THE LIBRARY AND SUBROUTINE BY BUILDING A
+  !   MESSAGE IN CHARACTER VARIABLE TEMP (NOT EXCEEDING 66 CHARACTERS)
+  !   AND SENDING IT OUT VIA XERPRN.  PRINT ONLY IF CONTROL FLAG
+  !   IS NOT ZERO.
+  !
+  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
+  !
+  !   IF LKNTRL IS POSITIVE, PRINT AN INTRODUCTORY LINE BEFORE
+  !   PRINTING THE MESSAGE.  THE INTRODUCTORY LINE TELLS THE CHOICE
+  !   FROM EACH OF THE FOLLOWING THREE OPTIONS.
+  !   1.  LEVEL OF THE MESSAGE
+  !          'INFORMATIVE MESSAGE'
+  !          'POTENTIALLY RECOVERABLE ERROR'
+  !          'FATAL ERROR'
+  !   2.  WHETHER CONTROL FLAG WILL ALLOW PROGRAM TO CONTINUE
+  !          'PROG CONTINUES'
+  !          'PROG ABORTED'
+  !   3.  WHETHER OR NOT A TRACEBACK WAS REQUESTED.  (THE TRACEBACK
+  !       MAY NOT BE IMPLEMENTED AT SOME SITES, SO THIS ONLY TELLS
+  !       WHAT WAS REQUESTED, NOT WHAT WAS DELIVERED.)
+  !          'TRACEBACK REQUESTED'
+  !          'TRACEBACK NOT REQUESTED'
+  !   NOTICE THAT THE LINE INCLUDING FOUR PREFIX CHARACTERS WILL NOT
+  !   EXCEED 74 CHARACTERS.
+  !   WE SKIP THE NEXT BLOCK IF THE INTRODUCTORY LINE IS NOT NEEDED.
+  !
+  IF (LKNTRL .GT. 0) THEN
+  !
+  !   THE FIRST PART OF THE MESSAGE TELLS ABOUT THE LEVEL.
+  !
+     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
+  !
+  !   THEN WHETHER THE PROGRAM WILL CONTINUE.
+  !
+     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
+  !
+  !   FINALLY TELL WHETHER THERE SHOULD BE A TRACEBACK.
+  !
+     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
+  !
+  !   NOW SEND OUT THE MESSAGE.
+  !
+  CALL XERPRN (' *  ', -1, MESSG, 72)
+  !
+  !   IF LKNTRL IS POSITIVE, WRITE THE ERROR NUMBER AND REQUEST A
+  !      TRACEBACK.
+  !
+  IF (LKNTRL .GT. 0) THEN
+     WRITE (TEMP, '(''ERROR NUMBER = '', I8)') NERR
+     DO I=16,22
+        IF (TEMP(I:I) .NE. ' ') GO TO 20
+     END DO
+  !
+   20   CALL XERPRN (' *  ', -1, TEMP(1:15) // TEMP(I:23), 72)
+     CALL FDUMP
+  ENDIF
+  !
+  !   IF LKNTRL IS NOT ZERO, PRINT A BLANK LINE AND AN END OF MESSAGE.
+  !
+  IF (LKNTRL .NE. 0) THEN
+     CALL XERPRN (' *  ', -1, ' ', 72)
+     CALL XERPRN (' ***', -1, 'END OF MESSAGE', 72)
+     CALL XERPRN ('    ',  0, ' ', 72)
+  ENDIF
+  !
+  !   IF THE ERROR IS NOT FATAL OR THE ERROR IS RECOVERABLE AND THE
+  !   CONTROL FLAG IS SET FOR RECOVERY, THEN RETURN.
+  !
+   30   IF (LEVEL.LE.0 .OR. (LEVEL.EQ.1 .AND. MKNTRL.LE.1)) RETURN
+  !
+  !   THE PROGRAM WILL BE STOPPED DUE TO AN UNRECOVERED ERROR OR A
+  !   FATAL ERROR.  PRINT THE REASON FOR THE ABORT AND THE ERROR
+  !   SUMMARY IF THE CONTROL FLAG AND THE MAXIMUM ERROR COUNT PERMIT.
+  !
+  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 SUBROUTINE XERMSG
Index: LMDZ6/trunk/libf/misc/xerprn.F
===================================================================
--- LMDZ6/trunk/libf/misc/xerprn.F	(revision 5214)
+++ 	(revision )
@@ -1,231 +1,0 @@
-*DECK XERPRN
-      SUBROUTINE XERPRN (PREFIX, NPREF, MESSG, NWRAP)
-      IMPLICIT NONE
-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 = '$$')
-      INTEGER N, I1MACH, I, LPREF, LWRAP, LENMSG, NEXTC
-      INTEGER LPIECE, IDELTA
-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: LMDZ6/trunk/libf/misc/xerprn.f90
===================================================================
--- LMDZ6/trunk/libf/misc/xerprn.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/xerprn.f90	(revision 5246)
@@ -0,0 +1,231 @@
+!DECK XERPRN
+SUBROUTINE XERPRN (PREFIX, NPREF, MESSG, NWRAP)
+  IMPLICIT NONE
+  !***BEGIN PROLOGUE  XERPRN
+  !***SUBSIDIARY
+  !***PURPOSE  Print error messages processed by XERMSG.
+  !***LIBRARY   SLATEC (XERROR)
+  !***CATEGORY  R3C
+  !***TYPE      ALL (XERPRN-A)
+  !***KEYWORDS  ERROR MESSAGES, PRINTING, XERROR
+  !***AUTHOR  Fong, Kirby, (NMFECC at LLNL)
+  !***DESCRIPTION
+  !
+  ! This routine sends one or more lines to each of the (up to five)
+  ! logical units to which error messages are to be sent.  This routine
+  ! is called several times by XERMSG, sometimes with a single line to
+  ! print and sometimes with a (potentially very long) message that may
+  ! wrap around into multiple lines.
+  !
+  ! PREFIX  Input argument of type CHARACTER.  This argument contains
+  !     characters to be put at the beginning of each line before
+  !     the body of the message.  No more than 16 characters of
+  !     PREFIX will be used.
+  !
+  ! NPREF   Input argument of type INTEGER.  This argument is the number
+  !     of characters to use from PREFIX.  If it is negative, the
+  !     intrinsic function LEN is used to determine its length.  If
+  !         it is zero, PREFIX is not used.  If it exceeds 16 or if
+  !     LEN(PREFIX) exceeds 16, only the first 16 characters will be
+  !     used.  If NPREF is positive and the length of PREFIX is less
+  !     than NPREF, a copy of PREFIX extended with blanks to length
+  !     NPREF will be used.
+  !
+  ! MESSG   Input argument of type CHARACTER.  This is the text of a
+  !     message to be printed.  If it is a long message, it will be
+  !     broken into pieces for printing on multiple lines.  Each line
+  !     will start with the appropriate prefix and be followed by a
+  !     piece of the message.  NWRAP is the number of characters per
+  !     piece; that is, after each NWRAP characters, we break and
+  !     start a new line.  In addition the characters '$$' embedded
+  !     in MESSG are a sentinel for a new line.  The counting of
+  !     characters up to NWRAP starts over for each new line.  The
+  !     value of NWRAP typically used by XERMSG is 72 since many
+  !     older error messages in the SLATEC Library are laid out to
+  !     rely on wrap-around every 72 characters.
+  !
+  ! NWRAP   Input argument of type INTEGER.  This gives the maximum size
+  !     piece into which to break MESSG for printing on multiple
+  !     lines.  An embedded '$$' ends a line, and the count restarts
+  !     at the following character.  If a line break does not occur
+  !     on a blank (it would split a word) that word is moved to the
+  !     next line.  Values of NWRAP less than 16 will be treated as
+  !     16.  Values of NWRAP greater than 132 will be treated as 132.
+  !     The actual line length will be NPREF + NWRAP after NPREF has
+  !     been adjusted to fall between 0 and 16 and NWRAP has been
+  !     adjusted to fall between 16 and 132.
+  !
+  !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+  !             Error-handling Package, SAND82-0800, Sandia
+  !             Laboratories, 1982.
+  !***ROUTINES CALLED  I1MACH, XGETUA
+  !***REVISION HISTORY  (YYMMDD)
+  !   880621  DATE WRITTEN
+  !   880708  REVISED AFTER THE SLATEC CML SUBCOMMITTEE MEETING OF
+  !       JUNE 29 AND 30 TO CHANGE THE NAME TO XERPRN AND TO REWORK
+  !       THE HANDLING OF THE NEW LINE SENTINEL TO BEHAVE LIKE THE
+  !       SLASH CHARACTER IN FORMAT STATEMENTS.
+  !   890706  REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO
+  !       STREAMLINE THE CODING AND FIX A BUG THAT CAUSED EXTRA BLANK
+  !       LINES TO BE PRINTED.
+  !   890721  REVISED TO ADD A NEW FEATURE.  A NEGATIVE VALUE OF NPREF
+  !       CAUSES LEN(PREFIX) TO BE USED AS THE LENGTH.
+  !   891013  REVISED TO CORRECT ERROR IN CALCULATING PREFIX LENGTH.
+  !   891214  Prologue converted to Version 4.0 format.  (WRB)
+  !   900510  Added code to break messages between words.  (RWC)
+  !   920501  Reformatted the REFERENCES section.  (WRB)
+  !***END PROLOGUE  XERPRN
+  CHARACTER(len=*) :: PREFIX, MESSG
+  INTEGER :: NPREF, NWRAP
+  CHARACTER(len=148) :: CBUFF
+  INTEGER :: IU(5), NUNIT
+  CHARACTER(len=2) :: NEWLIN
+  PARAMETER (NEWLIN = '$$')
+  INTEGER :: N, I1MACH, I, LPREF, LWRAP, LENMSG, NEXTC
+  INTEGER :: LPIECE, IDELTA
+  !***FIRST EXECUTABLE STATEMENT  XERPRN
+  CALL XGETUA(IU,NUNIT)
+  !
+  !   A ZERO VALUE FOR A LOGICAL UNIT NUMBER MEANS TO USE THE STANDARD
+  !   ERROR MESSAGE UNIT INSTEAD.  I1MACH(4) RETRIEVES THE STANDARD
+  !   ERROR MESSAGE UNIT.
+  !
+  N = I1MACH(4)
+  DO I=1,NUNIT
+     IF (IU(I) .EQ. 0) IU(I) = N
+  END DO
+  !
+  !   LPREF IS THE LENGTH OF THE PREFIX.  THE PREFIX IS PLACED AT THE
+  !   BEGINNING OF CBUFF, THE CHARACTER BUFFER, AND KEPT THERE DURING
+  !   THE REST OF THIS ROUTINE.
+  !
+  IF ( NPREF .LT. 0 ) THEN
+     LPREF = LEN(PREFIX)
+  ELSE
+     LPREF = NPREF
+  ENDIF
+  LPREF = MIN(16, LPREF)
+  IF (LPREF .NE. 0) CBUFF(1:LPREF) = PREFIX
+  !
+  !   LWRAP IS THE MAXIMUM NUMBER OF CHARACTERS WE WANT TO TAKE AT ONE
+  !   TIME FROM MESSG TO PRINT ON ONE LINE.
+  !
+  LWRAP = MAX(16, MIN(132, NWRAP))
+  !
+  !   SET LENMSG TO THE LENGTH OF MESSG, IGNORE ANY TRAILING BLANKS.
+  !
+  LENMSG = LEN(MESSG)
+  N = LENMSG
+  DO I=1,N
+     IF (MESSG(LENMSG:LENMSG) .NE. ' ') GO TO 30
+     LENMSG = LENMSG - 1
+  END DO
+   30   CONTINUE
+  !
+  !   IF THE MESSAGE IS ALL BLANKS, THEN PRINT ONE BLANK LINE.
+  !
+  IF (LENMSG .EQ. 0) THEN
+     CBUFF(LPREF+1:LPREF+1) = ' '
+     DO I=1,NUNIT
+        WRITE(IU(I), '(A)') CBUFF(1:LPREF+1)
+     END DO
+     RETURN
+  ENDIF
+  !
+  !   SET NEXTC TO THE POSITION IN MESSG WHERE THE NEXT SUBSTRING
+  !   STARTS.  FROM THIS POSITION WE SCAN FOR THE NEW LINE SENTINEL.
+  !   WHEN NEXTC EXCEEDS LENMSG, THERE IS NO MORE TO PRINT.
+  !   WE LOOP BACK TO LABEL 50 UNTIL ALL PIECES HAVE BEEN PRINTED.
+  !
+  !   WE LOOK FOR THE NEXT OCCURRENCE OF THE NEW LINE SENTINEL.  THE
+  !   INDEX INTRINSIC FUNCTION RETURNS ZERO IF THERE IS NO OCCURRENCE
+  !   OR IF THE LENGTH OF THE FIRST ARGUMENT IS LESS THAN THE LENGTH
+  !   OF THE SECOND ARGUMENT.
+  !
+  !   THERE ARE SEVERAL CASES WHICH SHOULD BE CHECKED FOR IN THE
+  !   FOLLOWING ORDER.  WE ARE ATTEMPTING TO SET LPIECE TO THE NUMBER
+  !   OF CHARACTERS THAT SHOULD BE TAKEN FROM MESSG STARTING AT
+  !   POSITION NEXTC.
+  !
+  !   LPIECE .EQ. 0   THE NEW LINE SENTINEL DOES NOT OCCUR IN THE
+  !                   REMAINDER OF THE CHARACTER STRING.  LPIECE
+  !                   SHOULD BE SET TO LWRAP OR LENMSG+1-NEXTC,
+  !                   WHICHEVER IS LESS.
+  !
+  !   LPIECE .EQ. 1   THE NEW LINE SENTINEL STARTS AT MESSG(NEXTC:
+  !                   NEXTC).  LPIECE IS EFFECTIVELY ZERO, AND WE
+  !                   PRINT NOTHING TO AVOID PRODUCING UNNECESSARY
+  !                   BLANK LINES.  THIS TAKES CARE OF THE SITUATION
+  !                   WHERE THE LIBRARY ROUTINE HAS A MESSAGE OF
+  !                   EXACTLY 72 CHARACTERS FOLLOWED BY A NEW LINE
+  !                   SENTINEL FOLLOWED BY MORE CHARACTERS.  NEXTC
+  !                   SHOULD BE INCREMENTED BY 2.
+  !
+  !   LPIECE .GT. LWRAP+1  REDUCE LPIECE TO LWRAP.
+  !
+  !   ELSE            THIS LAST CASE MEANS 2 .LE. LPIECE .LE. LWRAP+1
+  !                   RESET LPIECE = LPIECE-1.  NOTE THAT THIS
+  !                   PROPERLY HANDLES THE END CASE WHERE LPIECE .EQ.
+  !                   LWRAP+1.  THAT IS, THE SENTINEL FALLS EXACTLY
+  !                   AT THE END OF A LINE.
+  !
+  NEXTC = 1
+   50   LPIECE = INDEX(MESSG(NEXTC:LENMSG), NEWLIN)
+  IF (LPIECE .EQ. 0) THEN
+  !
+  !   THERE WAS NO NEW LINE SENTINEL FOUND.
+  !
+     IDELTA = 0
+     LPIECE = MIN(LWRAP, LENMSG+1-NEXTC)
+     IF (LPIECE .LT. LENMSG+1-NEXTC) THEN
+        DO I=LPIECE+1,2,-1
+           IF (MESSG(NEXTC+I-1:NEXTC+I-1) .EQ. ' ') THEN
+              LPIECE = I-1
+              IDELTA = 1
+              GOTO 54
+           ENDIF
+        END DO
+     ENDIF
+   54   CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
+     NEXTC = NEXTC + LPIECE + IDELTA
+  ELSEIF (LPIECE .EQ. 1) THEN
+  !
+  !   WE HAVE A NEW LINE SENTINEL AT MESSG(NEXTC:NEXTC+1).
+  !   DON'T PRINT A BLANK LINE.
+  !
+     NEXTC = NEXTC + 2
+     GO TO 50
+  ELSEIF (LPIECE .GT. LWRAP+1) THEN
+  !
+  !   LPIECE SHOULD BE SET DOWN TO LWRAP.
+  !
+     IDELTA = 0
+     LPIECE = LWRAP
+     DO I=LPIECE+1,2,-1
+        IF (MESSG(NEXTC+I-1:NEXTC+I-1) .EQ. ' ') THEN
+           LPIECE = I-1
+           IDELTA = 1
+           GOTO 58
+        ENDIF
+     END DO
+   58   CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
+     NEXTC = NEXTC + LPIECE + IDELTA
+  ELSE
+  !
+  !   IF WE ARRIVE HERE, IT MEANS 2 .LE. LPIECE .LE. LWRAP+1.
+  !   WE SHOULD DECREMENT LPIECE BY ONE.
+  !
+     LPIECE = LPIECE - 1
+     CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
+     NEXTC  = NEXTC + LPIECE + 2
+  ENDIF
+  !
+  !   PRINT
+  !
+  DO I=1,NUNIT
+     WRITE(IU(I), '(A)') CBUFF(1:LPREF+LPIECE)
+  END DO
+  !
+  IF (NEXTC .LE. LENMSG) GO TO 50
+  RETURN
+END SUBROUTINE XERPRN
Index: LMDZ6/trunk/libf/misc/xersve.F
===================================================================
--- LMDZ6/trunk/libf/misc/xersve.F	(revision 5214)
+++ 	(revision )
@@ -1,159 +1,0 @@
-*DECK XERSVE
-      SUBROUTINE XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL,
-     +   ICOUNT)
-      IMPLICIT NONE
-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
-      INTEGER,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/
-      INTEGER NERR,LEVEL,KONTRL
-      INTEGER NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG
-      INTEGER KFLAG, ICOUNT, NUNIT, KUNIT, IUNIT, I1MACH, I
-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: LMDZ6/trunk/libf/misc/xersve.f90
===================================================================
--- LMDZ6/trunk/libf/misc/xersve.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/xersve.f90	(revision 5246)
@@ -0,0 +1,159 @@
+!DECK XERSVE
+SUBROUTINE XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL, &
+        ICOUNT)
+  IMPLICIT NONE
+  !***BEGIN PROLOGUE  XERSVE
+  !***SUBSIDIARY
+  !***PURPOSE  Record that an error has occurred.
+  !***LIBRARY   SLATEC (XERROR)
+  !***CATEGORY  R3
+  !***TYPE      ALL (XERSVE-A)
+  !***KEYWORDS  ERROR, XERROR
+  !***AUTHOR  Jones, R. E., (SNLA)
+  !***DESCRIPTION
+  !
+  ! *Usage:
+  !
+  !    INTEGER  KFLAG, NERR, LEVEL, ICOUNT
+  !    CHARACTER * (len) LIBRAR, SUBROU, MESSG
+  !
+  !    CALL XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL, ICOUNT)
+  !
+  ! *Arguments:
+  !
+  !    LIBRAR :IN    is the library that the message is from.
+  !    SUBROU :IN    is the subroutine that the message is from.
+  !    MESSG  :IN    is the message to be saved.
+  !    KFLAG  :IN    indicates the action to be performed.
+  !                  when KFLAG > 0, the message in MESSG is saved.
+  !                  when KFLAG=0 the tables will be dumped and
+  !                  cleared.
+  !                  when KFLAG < 0, the tables will be dumped and
+  !                  not cleared.
+  !    NERR   :IN    is the error number.
+  !    LEVEL  :IN    is the error severity.
+  !    ICOUNT :OUT   the number of times this message has been seen,
+  !                  or zero if the table has overflowed and does not
+  !                  contain this message specifically.  When KFLAG=0,
+  !                  ICOUNT will not be altered.
+  !
+  ! *Description:
+  !
+  !   Record that this error occurred and possibly dump and clear the
+  !   tables.
+  !
+  !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+  !             Error-handling Package, SAND82-0800, Sandia
+  !             Laboratories, 1982.
+  !***ROUTINES CALLED  I1MACH, XGETUA
+  !***REVISION HISTORY  (YYMMDD)
+  !   800319  DATE WRITTEN
+  !   861211  REVISION DATE from Version 3.2
+  !   891214  Prologue converted to Version 4.0 format.  (BAB)
+  !   900413  Routine modified to remove reference to KFLAG.  (WRB)
+  !   900510  Changed to add LIBRARY NAME and SUBROUTINE to calling
+  !       sequence, use IF-THEN-ELSE, make number of saved entries
+  !       easily changeable, changed routine name from XERSAV to
+  !       XERSVE.  (RWC)
+  !   910626  Added LIBTAB and SUBTAB to SAVE statement.  (BKS)
+  !   920501  Reformatted the REFERENCES section.  (WRB)
+  !***END PROLOGUE  XERSVE
+  INTEGER,PARAMETER :: LENTAB=10
+  INTEGER :: LUN(5)
+  CHARACTER(len=*) :: LIBRAR, SUBROU, MESSG
+  CHARACTER(len=8) :: LIBTAB(LENTAB), SUBTAB(LENTAB), LIB, SUB
+  CHARACTER(len=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/
+  INTEGER :: NERR,LEVEL,KONTRL
+  INTEGER :: NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG
+  INTEGER :: KFLAG, ICOUNT, NUNIT, KUNIT, IUNIT, I1MACH, I
+  !***FIRST EXECUTABLE STATEMENT  XERSVE
+  !
+  IF (KFLAG.LE.0) THEN
+  !
+  !    Dump the table.
+  !
+     IF (NMSG.EQ.0) RETURN
+  !
+  !    Print to each unit.
+  !
+     CALL XGETUA (LUN, NUNIT)
+     DO KUNIT = 1,NUNIT
+        IUNIT = LUN(KUNIT)
+        IF (IUNIT.EQ.0) IUNIT = I1MACH(4)
+  !
+  !       Print the table header.
+  !
+        WRITE (IUNIT,9000)
+  !
+  !       Print body of table.
+  !
+        DO I = 1,NMSG
+           WRITE (IUNIT,9010) LIBTAB(I), SUBTAB(I), MESTAB(I), &
+                 NERTAB(I),LEVTAB(I),KOUNT(I)
+        END DO
+  !
+  !       Print number of other errors.
+  !
+        IF (KOUNTX.NE.0) WRITE (IUNIT,9020) KOUNTX
+        WRITE (IUNIT,9030)
+     END DO
+  !
+  !    Clear the error tables.
+  !
+     IF (KFLAG.EQ.0) THEN
+        NMSG = 0
+        KOUNTX = 0
+     ENDIF
+  ELSE
+  !
+  !    PROCESS A MESSAGE...
+  !    SEARCH FOR THIS MESSG, OR ELSE AN EMPTY SLOT FOR THIS MESSG,
+  !    OR ELSE DETERMINE THAT THE ERROR TABLE IS FULL.
+  !
+     LIB = LIBRAR
+     SUB = SUBROU
+     MES = MESSG
+     DO 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
+     END DO
+  !
+     IF (NMSG.LT.LENTAB) THEN
+  !
+  !       Empty slot found for new message.
+  !
+        NMSG = NMSG + 1
+        LIBTAB(I) = LIB
+        SUBTAB(I) = SUB
+        MESTAB(I) = MES
+        NERTAB(I) = NERR
+        LEVTAB(I) = LEVEL
+        KOUNT (I) = 1
+        ICOUNT    = 1
+     ELSE
+  !
+  !       Table is full.
+  !
+        KOUNTX = KOUNTX+1
+        ICOUNT = 0
+     ENDIF
+  ENDIF
+  RETURN
+  !
+  ! Formats.
+  !
+ 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 SUBROUTINE XERSVE
Index: LMDZ6/trunk/libf/misc/xgetua.F
===================================================================
--- LMDZ6/trunk/libf/misc/xgetua.F	(revision 5214)
+++ 	(revision )
@@ -1,53 +1,0 @@
-*DECK XGETUA
-      SUBROUTINE XGETUA (IUNITA, N)
-      IMPLICIT NONE
-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)
-      INTEGER IUNITA, N, J4SAVE, INDEX, I
-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: LMDZ6/trunk/libf/misc/xgetua.f90
===================================================================
--- LMDZ6/trunk/libf/misc/xgetua.f90	(revision 5246)
+++ LMDZ6/trunk/libf/misc/xgetua.f90	(revision 5246)
@@ -0,0 +1,53 @@
+!DECK XGETUA
+SUBROUTINE XGETUA (IUNITA, N)
+  IMPLICIT NONE
+  !***BEGIN PROLOGUE  XGETUA
+  !***PURPOSE  Return unit number(s) to which error messages are being
+         ! sent.
+  !***LIBRARY   SLATEC (XERROR)
+  !***CATEGORY  R3C
+  !***TYPE      ALL (XGETUA-A)
+  !***KEYWORDS  ERROR, XERROR
+  !***AUTHOR  Jones, R. E., (SNLA)
+  !***DESCRIPTION
+  !
+  ! Abstract
+  !    XGETUA may be called to determine the unit number or numbers
+  !    to which error messages are being sent.
+  !    These unit numbers may have been set by a call to XSETUN,
+  !    or a call to XSETUA, or may be a default value.
+  !
+  ! Description of Parameters
+  !  --Output--
+  !    IUNIT - an array of one to five unit numbers, depending
+  !            on the value of N.  A value of zero refers to the
+  !            default unit, as defined by the I1MACH machine
+  !            constant routine.  Only IUNIT(1),...,IUNIT(N) are
+  !            defined by XGETUA.  The values of IUNIT(N+1),...,
+  !            IUNIT(5) are not defined (for N .LT. 5) or altered
+  !            in any way by XGETUA.
+  !    N     - the number of units to which copies of the
+  !            error messages are being sent.  N will be in the
+  !            range from 1 to 5.
+  !
+  !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+  !             Error-handling Package, SAND82-0800, Sandia
+  !             Laboratories, 1982.
+  !***ROUTINES CALLED  J4SAVE
+  !***REVISION HISTORY  (YYMMDD)
+  !   790801  DATE WRITTEN
+  !   861211  REVISION DATE from Version 3.2
+  !   891214  Prologue converted to Version 4.0 format.  (BAB)
+  !   920501  Reformatted the REFERENCES section.  (WRB)
+  !***END PROLOGUE  XGETUA
+  DIMENSION IUNITA(5)
+  INTEGER :: IUNITA, N, J4SAVE, INDEX, I
+  !***FIRST EXECUTABLE STATEMENT  XGETUA
+  N = J4SAVE(5,0,.FALSE.)
+  DO I=1,N
+     INDEX = I+4
+     IF (I.EQ.1) INDEX = 3
+     IUNITA(I) = J4SAVE(INDEX,0,.FALSE.)
+  END DO
+  RETURN
+END SUBROUTINE XGETUA
